Imagine you roll $d$ sided dice until the sum of the dice is gt;d$. Let $N$ be the number of dice rolls you need. i.e. if we denote the random variables be $X_1,X_2,\ldots \in \{1,\ldots,d\}$ then $ N = \min \{n : X_1 + \ldots +X_n > d \}$ # Question: What is $\mathbb{E}[N]$? It turns out that the answer is $\mathbb{E}[N]=\left( 1+\frac{1}{d} \right)^d$ which means that in the limit it takes $e$ dice rolls on average: $\lim_{d\to \infty}\mathbb{E}[N] = e = 2.718...$ # Combinatorics Solution We use the formula $\mathbb{E}[N] = \sum_{k=0}^n \mathbb{P}(N > k )$ Note that the event $N > k$ happens if and only if $X_1 + \ldots + X_k \leq d$. So $\mathbb{P}(N > k ) = \mathbb{P}(X_1+\ldots +X_k \leq d)$. To evaluate this probability, we must simply count the configurations of $X_1,\ldots,X_k$ such that $X_1+\ldots+X_k \leq d$. We claim that the number of such configurations is: $ \left| \left\{ (X_1,X_2,\ldots,X_k) | X_1+\ldots+X_k \leq d \right\} \right| = \binom{d}{k} $ This is because there is a bjiection between $k$ element subset of $\{1,\ldots,d\}$ and this set by mapping any $k$ element subset $A_1 < A_2 <\ldots < A_k$ by the correspondence: $ A_k := \sum_{i=1}^k X_i$or conversely: $X_i = A_{i+1}-A_i$ Hence the probability we are interested in is: $\begin{align} \mathbb{P}\left( X_1 + \ldots + X_k \leq d\right) &= \frac{\left| \left\{ (X_1,X_2,\ldots,X_k) | X_1+\ldots+X_k \leq d \right\} \right|}{d^k}\\ &= \frac{1}{d^k}\binom{d}{k} \\&= \frac{1(1-\frac{1}{d})(1-\frac{2}{d})\cdots(1-\frac{k-1}{d})}{k!}\end{align}$ (Notice that in the limit $d \to \infty$, that this $\to \frac{1}{k!}$.) Hence $\mathbb{E}[N]$ is the sum of these probabilities. To simplify thing, we can actually see from the Binomial Theorem that $\mathbb{E}[N] = \sum_{k=0}^d \frac{1}{d^k}\binom{d}{k} = \left(1+\frac{1}{d}\right)^d$ # Functional Equation Solution Consider the following function for any $x\in\{0,1,2,\ldots,d\}$: $f_d(x) = \mathbb{E}[\text{Number rolls needed to exceed a sum of $x$}]$ The original problem is to compute $f_d(d)$. Note that when we roll the dice one time, the sum remaining my decrease by $1$ or $2$ or $3$ etc each with probability $1/d$. We hence see that $f_d(x) = 1 + \frac{1}{d} \left( f_d(x-1) + f_d( x-2 ) + \ldots + f_d(x-d) \right)$ (Note that here we use the convention that $f_d(0)=1$ and $f_d(x)=0$ for all $x<0$ .) If we apply the same idea to evaluate $f_d(x-1)$ we get $f_d(x-1) = 1 + \frac{1}{d} \left( f_d(x-2) + f_d( x-3 ) + \ldots + f_d(x-d) + 0 \right)$ (Note the last terms is 0 because $f_d(x-d-1)$ is always $0$ for $x \in \{0,1,\ldots,d\}$) So subtracting these two equations we remain with: $f_d(x) - f_d(x-1) = \frac{1}{d}f_d(x-1)$This can be rearranged to get $f_d(x)=\frac{d}{d-1}f_d(x-1)$, which along with $f_d(0)=1$, gives a type of "discrete'' exponential growth solution for the function $f_d(x) = \left( \frac{d+1}{d} \right)^x$, so $f_d(d)=\left( \frac{d+1}{d} \right)^d$ ## Differential Equation Limit To get a nice limit you can rescale the $f$: For $q \in (0,1)$, consider $g_d(q) = f_d(qd)$. This is the expected number of rolls to exceed $x=qd$ (i.e. $q$ is the fraction we need to exceed) We can move $\frac{1}{d}$ to the other side, we have: $\frac{g_d(x) - g_d(x-\frac{1}{d})}{\frac{1}{d}} = g_d(x-\frac{1}{d})$ In the limit we will see that $\frac{d}{dx}g (x) = g(x)$ which again gives the solution $g(x) = e^x$. # Python Code for simulations This uses the above to formulas and also runs a Monte Carlo simulation to confirm its working as intended! ```python import numpy as np d = 100 # number of sides on the dice ### Using combinatorial formula max_terms = min(12,d) #maximum numbe of terms to compute. You will be off by 1/max_terms! in the limit j = np.arange(max_terms+1) #the sequence 0,1,2,3...d one_minus_j_over_d = 1 - j/d #the sequence 1, 1-1/d, 1-2/d, ... k_factorial = np.cumprod(j+1) #the sequence 1!,2!,3!,... #the sequence of probabilities according the above combinatorial formula Prob_geq_k = np.cumprod(one_minus_j_over_d)/k_factorial #P(N>1),P(N>2) etc E_N_combinat = 1+np.sum(Prob_geq_k) #the expected value (must add 1 for P(N>0)) print(f"{E_N_combinat=}") ### Using the functional equation solution E_N_function = (1.0+1.0/d)**(d) print(f"{E_N_function=}") #### Monte Carlo Trial N_Monte_Carlo = 2**22 max_dice_rolls_per_trial = max_terms #give up after this many rolls! Again will cause us to be off by 1/max_terms! dice = np.random.randint(1,d+1,size=(max_dice_rolls_per_trial,N_Monte_Carlo)) dice_sums = np.cumsum(dice,axis=0) dice_sums_leq_d = (dice_sums <= d) N_rolls_needed = 1+np.sum(dice_sums_leq_d,axis=0) E_N_MonteCarlo = np.mean(N_rolls_needed) print(f"{E_N_MonteCarlo=}") ``` Sample output when $d=100$: ``` E_N_combinat=2.70481382941684 E_N_function=2.7048138294215285 E_N_MonteCarlo=2.705012321472168 ``` # Setup for an Explainer Video 1. Rolling physical dice: How many rolls to exceed the sum, approximately 3? 2. Is it pi? No! 3. Time lapse of calculating the probabilities (slow way!) 4. One weird trick to figure it out: The Bellman Equation also used in RL 5. Set up with a function and the solution