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

Question

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)

Answer 1

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

Answer 2

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