# Expected of number of discrete uniform variables whose sum is bigger than k (from characteristic function of discrete Irwin–Hall distribution?)

## 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
}


• 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. – whuber Jan 21 at 19:06
• Thank you for supplying the results of your simulation! – whuber Jan 23 at 15:02

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.$$

$$\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$$).

• 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? – ManuHaq Jan 22 at 13:38
• (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. – whuber Jan 22 at 14:31
• 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!) – ManuHaq Jan 22 at 16:41
• 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. – whuber Jan 22 at 17:03
• 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$". – ManuHaq Jan 22 at 17:40

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$$.

• 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?) – ManuHaq Jan 22 at 13:45
• BTW, I hade made a mistake in the code (now corrected); the sum of values should be more than 10 (so 11 upwards.) – ManuHaq Jan 22 at 14:02
• 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$. – jbowman Jan 22 at 14:34