3
$\begingroup$

The problem

Imagine we keep uniformly drawing $n$ integers $X_i$ from {0, 1, ..., 9} so that their sum is more than $10$. For instance, one draw would be {1, 0, 2, 5, 3}, hence $n=5$, and repeat this procedure over and over. What would the expected value of $n$ be (analytically calculated)? By simulation (over 10 million trials), it is 3.063577.

Apparently, if that were uniformly drawn from $[0,10]$, this expected value would be equal to the Euler number $e$.

[I very recently asked it on Mathematics StackExchange and am waiting for answers (and will post them here), but I thought maybe it may be more suitable for CV. Sorry if that makes it off-topic.]

What I have been trying to do

Consulting the paper Polynomial coefficients and distribution of the sum of discrete uniform variables by Caiado & Rathie (2007), I suspect (given $Y = \sum_{i=1}^{n}X_i$) the characteristic function of the distribution of $Y$ is something of the following form (Equation 2.3 in the paper)

$$ \Phi_Y(t) = \left( \sum_{p=0}^{k} \frac{e^{i.t.p}}{k+1} \right)^n , \forall t \in \mathbb{R}, i=\sqrt{-1} $$

If I understand it correctly, $k$ should be equal to 9 (right?)

If I have been correct, I tried calculating the inverse Fourier transform of $\Phi_Y(t)$ for $k=9$ and calculate its expected value but it is getting too complicated—and I suspect I'm very wrong here.

I went through multiple "similar" questions on SE (e.g., +, +, +, +, +, and +) but I'm too confused to infere something useful from them.

Is there an answer to my question?

Many thanks in advance!


Some R code for numerical estimation

This is the R code I used to calculate it numerically (edited since the first submission):

N <- 1e+7
s.list <- n.list <- rep(NA, N)

for(i in 1:N){
  s <- 0
  n <- 0
  seed <- i
  while(s < 11){
    s <- s + (sample(10, 1, replace = TRUE) - 1)
    n <- n + 1
  }
  if(!(i %% 10000)) print(paste("At iteration", round(i/1000,1), "K, s is", s, "and n is", n))
  s.list[i] <- s
  n.list[i] <- n
}

answer <- mean(n.list)
$\endgroup$
  • 3
    $\begingroup$ The thread at stats.stackexchange.com/questions/145621 concerns the closely related problem of finding the distribution of the sum. The methods described there can be applied to finding the distribution of the number of rolls. $\endgroup$ – whuber Jan 21 at 19:06
  • 1
    $\begingroup$ Thank you for supplying the results of your simulation! $\endgroup$ – whuber Jan 23 at 15:02
2
$\begingroup$

The probability generating function (pgf) of $X_i$ is

$$p_{b}(x) = (1 + x^1 + x^2 + \cdots + x^{b-1})/b = b^{-1} \frac{1-x^b}{1-x}$$

with $b=10.$ (For related examples and explanations see, inter alia, 1, 2, 3, and 4; or consult Wikipedia.)

Because the $X_i$ are independent, the pgf of the sum of $n$ of them is

$$p_b(x)^n = b^{-n} \frac{(1-x^b)^n}{(1-x)^n}.$$

Let's generalize the question by letting $N$ be the smallest $n$ for which $X_1+X_2+\cdots+X_n \ge a$ for a given natural number $a.$ (We will focus on $a=b$ and $a=b+1$ for simplicity.) We can find the distribution of $N$ by focusing on sums that are less than $a.$ This permits us to work with the pgf modulo $x^a$ (see, e.g., 5, 6, or 7).

The Binomial Theorem, applied separately to the numerator and denominator of $p_b(x)$ modulo $x^a,$ states

$$\eqalign{ p_b(x)^n &= b^{-n} (1-x^b)^n (1-x)^{-n} \\ &= b^{-n}\left(\sum_{i=0}^{\min(n,(a-1)/b)} \binom{n}{i} (-1)^i x^{bi}\right)\left(\sum_{k=0}^{a-1} (-1)^k \binom{-n}{k} x^k\right) \mod x^a.}$$

The chance that the sum is still less than $a$ (that is, the survival function $S_a$ of the random variable $N$) after $n$ rolls is the sum of these coefficients, easily evaluated by plugging in $x=1.$ I will derive explicit solutions for two cases:

$$S_b(n) = b^{-n} \sum_{k=0}^{b-1} (-1)^k \binom{-n}{k} = b^{-n} \sum_{k=0}^{b-1} \binom{n-1+k}{k} = b^{-n} \binom{n-1+b}{b-1}$$

and

$$\eqalign{p_b(x)^{n} &= (1-n x^b)\, b^{-n} \sum_{k=0}^{b} (-1)^k \binom{-n}{k}x^k \mod x^{b+1}\\ & = b^{-n} \left(-n x^b + \sum_{k=0}^{b} \binom{n-1+k}{k}x^k\right) \mod x^{b+1}}$$

whence

$$S_{b+1}(n) = b^{-n} \left(-n + \binom{n+b}{b}\right).$$

We have thereby obtained a formula for the entire distribution of the number of rolls needed to equal or exceed $a.$

Since the expectation of a discrete non-negative random variable is the sum of the values of its survival function, the answers are

$$\sum_{n=0}^\infty S_b(n) =\sum_{n=0}^\infty b^{-n}\binom{n-1+b}{b-1} = \left(1-\frac{1}{b}\right)^{-b}.$$

$$\sum_{n=0}^\infty S_{b+1}(n) =\sum_{n=0}^\infty b^{-n}\left(-n + \binom{n+b}{b}\right) = \frac{b}{b-1}\left(\left(1-\frac{1}{b}\right)^{-b}-\frac{1}{b-1}\right).$$

With $b=10$ the answers are

$$(1 - 1/10)^{-10} = \frac{2^{10}\, 5^{10}}{3^{20}} = 2.867\,971\,990\,792\,441\,313\,322\, \ldots$$

and

$$ \frac{10}{9}\left(\left(\frac{9}{10}\right)^{-10} - \frac{1}{9}\right) = \frac{96125795110}{31381059609} \approx 3.063\,178\,755\,201\,478\,002\,456\,\ldots$$

Obviously as $b\to\infty,$ $(1-1/b)^{-b}\to e.$ The difference between this expectation and $e = \exp(1)$ is positive and behaves asymptotically as $O(1/b).$ Thus, asymptotically the expectations are both $e + O(1/b),$ as suggested in the question (by comparison to the continuous uniform distribution of $X_i$).

$\endgroup$
  • $\begingroup$ Thanks a lot @whuber--I have some questions regarding your answer. Firstly, comparing to this derivation, shouldn't the pgf have an additional $x$ in the numerator? Also, I'm lost after the second line of math--ain't we interested in sums more than b? $\endgroup$ – ManuHaq Jan 22 at 13:38
  • 1
    $\begingroup$ (1) The pgf is correct, because the values of $X_i$ range from $0$ through $b-1.$ (2) When you know the chance that a sum is less than $b,$ you can easily compute the chance it equals or exceeds $b.$ Another proof is in the pudding: the answer is exact, simple to compute, and agrees with the simulations. I have applied further simplifications (compared to my original post) to demonstrate this. $\endgroup$ – whuber Jan 22 at 14:31
  • 1
    $\begingroup$ Ah, I see (1). My confusion about (2) was that I couldn't spot $1-P(s>10)$ in the derivation. Otherwise, the pudding is super delicious (and thanks for the added references, I'm reading them now!) $\endgroup$ – ManuHaq Jan 22 at 16:41
  • $\begingroup$ There is a slight difference in interpretation: the edits to the question clarify you are asking for the expected number of rolls to exceed $b$ rather than equal or exceed $b$ (my interpretation). This is a complication, but not a serious one: you have to work modulo $x^{b+1},$ for which the pdf is congruent to $(1-x^b)/(1-x).$ Raising it to the $n^\text{th}$ power modulo $x^{b+1}$ gives $(1-x^b)/(1-x)^n=1/(1-x)^n - x^b/(1-x)^n.$ That second term "corrects" the first term (by shifting everything by $b$ and subtracting). You can still write down a closed formula for the solution. $\endgroup$ – whuber Jan 22 at 17:03
  • $\begingroup$ Hmmm. I should still admit that I have difficulty pushing that correction forward in derivation. Moreover, in the main answer (after the last edit 42mins ago), I cannot see how you get from "The chance that the sum is still less than $b$ ..." to "We have thereby ... needed to equal or exceed $b$". $\endgroup$ – ManuHaq Jan 22 at 17:40
4
$\begingroup$

You can calculate this recursively, using a (backwards) dynamic programming-like iteration.

Let us construct a state variable $s \in \{0, \dots, 10\}$. $s$ represents the current value of the sum, except for state $10$, which represents any sum $>=10$. You start in state $s=0$.

If you are in state $10$, obviously you are done, and the expected number of remaining transitions to reach state $10$ is zero.

If you are in state $9$, you have a probability $0.9$ of transitioning to state $10$ and $0.1$ of remaining in state $9$. Using well-known results, this means your expected time-to-$10$ is $1.111\dots$ draws.

If you are in state $8$, using a similar logic as above, your expected time to transition to some higher state in $1.111\dots$ draws. When you do, you will wind up in state $9$ with probability $1/9$ and state $10$ with probability $8/9$. Taking the appropriate weighted average of the time-to-$10$ of states $9$ and $10$ gives an expected time-to-$10$ of state $8$ equal to $1.111\dots + 1.111\dots/9 = 1.234568$.

If you are in state $7$, using the same logic as above, ...

In each state, you have an expected time-to-transition to some higher state of $1.111$ draws, and when you do, you just take the transition probability-weighted average of the time-to-$10$ of the higher states to derive the time-to-$10$ of the current state.

This recursion can be easily implemented; I've chosen to do so in C++ using Rcpp, as the indexes in C/C++ start at $0$ but in R start at $1$, requiring some minor mental agility to correct for given that your state space starts with $0$:

> source_code <- '
+ #include <Rcpp.h>
+ using namespace Rcpp;
+ 
+ // [[Rcpp::export]]
+ NumericVector e_draws() 
+ {
+   NumericVector s(10);
+   s[9] = 10.0/9.0;
+   for (int i = 8; i >= 0; --i) {
+     s[i] = 10.0/9.0;
+     for (int j = i+1; j <= 9; ++j) {
+       s[i] += s[j]/9;
+     }
+   }
+   return s;
+ }
+ '
> 
> sourceCpp(code=source_code)
> 
> e_draws()
 [1] 2.867972 2.581175 2.323057 2.090752 1.881676 1.693509 1.524158 1.371742 1.234568 1.111111

... and the answer is the first element of the vector, $2.867972$.

$\endgroup$
  • $\begingroup$ Thanks for the interesting approach to it! But how does it not converge to the result @whuber reaches above? (And isn't there a closed-form, analytical solution to it?) $\endgroup$ – ManuHaq Jan 22 at 13:45
  • $\begingroup$ BTW, I hade made a mistake in the code (now corrected); the sum of values should be more than 10 (so 11 upwards.) $\endgroup$ – ManuHaq Jan 22 at 14:02
  • 1
    $\begingroup$ It does converge to the same result; you may be looking at @whuber's last calculation, which is for $b=10^7$. Partway through his answer is the number $2.867 971 990 792 44$, which rounds off to $2.867972$. $\endgroup$ – jbowman Jan 22 at 14:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.