# Expected number of dice rolls require to make a sum greater than or equal to K?

A 6 sided die is rolled iteratively. What is the expected number of rolls required to make a sum greater than or equal to K?

Before Edit

P(Sum>=1 in exactly 1 roll)=1
P(Sum>=2 in exactly 1 roll)=5/6
P(Sum>=2 in exactly 2 rolls)=1/6
P(Sum>=3 in exactly 1 roll)=5/6
P(Sum>=3 in exactly 2 rolls)=2/6
P(Sum>=3 in exactly 3 rolls)=1/36
P(Sum>=4 in exactly 1 roll)=3/6
P(Sum>=4 in exactly 2 rolls)=3/6
P(Sum>=4 in exactly 3 rolls)=2/36
P(Sum>=4 in exactly 4 rolls)=1/216


After Edit

P(Sum>=1 in atleast 1 roll)=1
P(Sum>=2 in atleast 1 roll)=5/6
P(Sum>=2 in atleast 2 rolls)=1
P(Sum>=3 in atleast 1 roll)=4/6
P(Sum>=3 in atleast 2 rolls)=35/36
P(Sum>=3 in atleast 3 rolls)=1
P(Sum>=4 in atleast 1 roll)=3/6
P(Sum>=4 in atleast 2 rolls)=33/36
P(Sum>=4 in atleast 3 rolls)=212/216
P(Sum>=4 in atleast 4 rolls)=1


I am not sure this is correct first of all and but I think this probability is related to the expected number of rolls?

But I don't know how to proceed further. Am I proceeding in the right direction?

• How did you get $P(S\geq 2 \text{ in 2 rolls})$? – Glen_b -Reinstate Monica Apr 13 '15 at 13:43
• @Glen_b You have to get a number less than 2 in the first roll which is 1. So probability of getting 1 is 1/6 and the second roll can be any number. if you get a number greater or equal to 2 in the first roll, then you wont go for a second roll. – Usual Suspect Apr 13 '15 at 13:50
• Ah, I see what's going on. You don't describe that as "P(S\geq 2 in 2 rolls)"; that expression implies the number of rolls is fixed. What you want is either "P(exactly 2 rolls required to get $S\geq 2$)" or "P(at least 2 rolls required to get $S\geq 2$)". – Glen_b -Reinstate Monica Apr 13 '15 at 14:08
• @Glen_b Yeah thats the confusion. P(exactly 2 rolls required to get S>2) i guess. All i ultimately want to calculate is Expected number of rolls to reach a sum greater than K? – Usual Suspect Apr 13 '15 at 14:12
• @Glen_b should I use atleast or exactly for this purpose? And how to calculate expected number of rolls for larger sum like 10000? – Usual Suspect Apr 14 '15 at 1:41

This is so far only some ideas for another, more exact, approach, based on the same observation that my first answer. With time I will extend this ...

First, some notation. Let $K$ be some given, positive (large) integer. We want the distribution of $N$, which is the minimum number of throws of an ordinary dice to get sum at least $K$. So, first we define $X_i$ as the outcome of dice throw $i$, and $X^{(n)}=X_1+\dots+X_n$. If we can find the distribution of $X^{(n)}$ for all $n$ then we can find the distribution of $N$ by using $$P(N \ge n)= P(X_1+\dots+X_n \le K),$$ and we are done.

Now, the possible values for $X_1+\dots+X_n$ are $n,n+1,n+2,\dots,6n$, and for $k$ in that range, to find the probability $P(X_1+\dots+X_n=k)$, we need to find the total number of ways to write $k$ as a sum of exactly $n$ integers, all in the range $1,2,\dots,6$. But that is called an restricted integer composition, a problem well studied in combinatorics. Some related questions on math SE is found by https://math.stackexchange.com/search?q=integer+compositions

So searching and studying that combinatorics literature we can get quiet precise results. I will follow up on that, but later ...

There is a simple closed formula in terms of the roots of a degree-6 polynomial.

It's actually a little easier to consider a general fair die with $$d\ge 2$$ faces labeled with the numbers $$1,2,\ldots, d.$$

Let $$e_k$$ be the expected number of rolls needed to equal or exceed $$k.$$ For $$k\le 0,$$ $$e_k=0.$$ Otherwise the expectation is one more than the expectation of the number of rolls to reach the immediately preceding value, which would be among $$k-d, k-d+1, \ldots, k-1,$$ whence

$$e_k = 1 + \frac{1}{d}\left(e_{k-d} + e_{k-d+1} + \cdots + e_{k-1}\right).\tag{1}$$

This linear recurrence relation has a solution in the form

$$e_k = \frac{2k}{d+1} + \sum_{i=1}^d a_i \lambda_i^k\tag{2}$$

where the $$\lambda_i$$ are the $$d$$ complex roots of the polynomial

$$T^d - \frac{1}{d}(T^{d-1} + T^{d-2} + \cdots + T + 1).\tag{3}$$

The constants $$a_i$$ are found by applying the solution $$(2)$$ to the values $$k=-(d-1), -(d-2), \ldots, -1, 0$$ where $$e_k=0$$ in every case. This gives a set of $$d$$ linear equations in the $$d$$ constants and it has a unique solution. That the solution works can be demonstrated by verifying the recurrence $$(1)$$ using the fact that every root satisfies $$(3):$$

\eqalign{ 1 + \frac{1}{d}\sum_{j=1}^{d} e_{k-j} &= 1 + \frac{1}{d}\sum_{j=1}^{d} \left(\frac{2(k-j)}{d+1} + \sum_{i=1}^d a_i \lambda_i^{k-j}\right) \\ &= \frac{2k}{d+1} + \sum_{i=1}^d a_i \lambda_i^{k-d}\left[\frac{1}{d}(1 + \lambda_i + \cdots + \lambda_i^{d-1})\right] \\ &= \frac{2k}{d+1} + \sum_{i=1}^d a_i \lambda_i^{k-d}\lambda_i^d \\ &= \frac{2k}{d+1} + \sum_{i=1}^d a_i \lambda_i^k = e_k. }

This closed form solution gives us good ways to approximate the answer as well as to evaluate it accurately. (For small to modest values of $$k,$$ direct application of the recurrence is an effective computational technique.) For example, with $$d=6$$ we can readily compute

$$e_{1\,000\,000} = 285714.761905\ldots$$

For approximations, there will be a unique largest root $$\lambda_{+}=1$$ so eventually (for sufficiently large $$k$$) the term $$\lambda_{+}^k$$ will dominate the $$d$$ terms in $$(2).$$ The error will decrease exponentially according to the second smallest norm of the roots. Continuing the example with $$k=6,$$ the coefficient of $$\lambda_{+}$$ is $$a_{+}=0.4761905$$ and the next-smallest norm is $$0.7302500.$$ (Incidentally, the other $$a_i$$ tend to be very close to $$1$$ in size.) Thus we may approximate the previous value as

$$e_{1\,000\,000} \approx \frac{2\times 10^6}{6+1} + 0.4761905 = 285714.761905\ldots$$

with an error on the order of $$0.7302500^{10^6} \approx 10^{-314\,368}.$$

To demonstrate how practical this solution is, here is R code that returns a function to evaluate $$e_k$$ for any $$k$$ (within the scope of double precision floating point calculations) and not overly large $$d$$ (it will bog down once $$d\gg 100$$):

die <- function(d, mult=1, cnst=1, start=rep(0,d)) {
# Create the companion matrix (its eigenvalues are the lambdas).
X <- matrix(c(0,1,rep(0,d-1)),d,d+1)
X[, d] <- mult/d
lambda <- eigen(X[, 1:d], symmetric=FALSE, only.values=TRUE)$values # Find the coefficients that agree with the starting values. u <- 2*cnst/(d+1) a <- solve(t(outer(lambda, 1:d, ^)), start - u*((1-d):0)) # This function assumes the starting values are all real numbers. f <- Vectorize(function(i) Re(sum(a * lambda ^ (i+d))) + u*i) list(f=f, lambda=lambda, a=a, multiplier=mult, offset=cnst) }  As an example of its use, here it calculates the expectations for $$k=1,2,\ldots, 16:$$ round(die(6)$f(1:10), 3)


1.000 1.167 1.361 1.588 1.853 2.161 2.522 2.775 3.043 3.324 3.613 3.906 4.197 4.476 4.760 5.046

The object it returns includes the roots $$\lambda_i$$ and their multipliers $$a_i$$ for further analysis. The first component of the multipliers array is the useful coefficient $$a_{+}.$$

(If you're curious what the other parameters of die are for, execute die(2, 2, 0, c(1,0))$f(1:10) and see whether you recognize the output ;-). This generalization assisted in developing and testing the function.) • +1. The function die gives an error for me: object 'phi' not found. – COOLSerdash Sep 13 at 11:38 • @COOL Thanks for checking. A last-minute change of variable name (from phi to a) to match the text was the culprit. I have fixed (and checked) it. – whuber Sep 13 at 12:32 there is no way to get exact expected number of rolls in general, but for a K. Let N be the event of of expected rolling to get sum=>K. for K=1, E(N)=1 for K=2,$E(N)=(\frac{5}{6}+2*1)/(\frac{5}{6}+1)=\frac{17}{11}$and so on. It will be going difficult to get E(N) for large K. for example,for K=20 you'll need to expect from (4 rolls,20 rolls) Central Limit Theorem will be more benefitiary with some % confidence. as we know occurrence is uniformly distributed, for large values of K. $$K(Sum)~follows~N(3.5N,\frac{35N}{12})$$(Normal Distribution) Now you need "N" to get Sum at least K.... we convert it in standard normal distribution.$$\frac{K-3.5N}{\sqrt{\frac{35N}{12}}}=Z_\alpha$$ where$\alpha=1-confidence$% You can get Z values from "Standard Normal Tables" or from here for example$Z_{0.01}=2.31,Z_{0.001}=2.98$You know K,Z(at any error) ........ then you can get N=E(N) at some confidence % by solving equation. • How you calculated those probabilities? How did you arrive at that E(N) equation? – Usual Suspect Apr 13 '15 at 13:58 • @UsualSuspect P(Sum>=2 in 1 roll)=5/6( you know) P(Sum>=2 in 2 rolls)=1(because you must get sum at least 2 from 2 rollings) and for E(N).........its just a Expected mean – Hemant Rupani Apr 13 '15 at 14:00 • Sorry i dint mention. Its not atleast, exactly 2 rolls. I understood the E(N) equation now. – Usual Suspect Apr 13 '15 at 14:20 • @UsualSuspect ohh! by the way if you need E(N) for any particular K, then I can make it :). – Hemant Rupani Apr 13 '15 at 14:48 • i need for k=20 and k=10000. Its better if you explain me rather than straightly giving answers. – Usual Suspect Apr 13 '15 at 14:57 I will give one method to find an approximate solution. First, let $$X_i$$ be the random variable, "result of throw $$i$$ with the dice" and let $$N$$ be the number of throws necessary to reach a sum at least $$k$$. Then we have that $$P(N \ge n) = P(X_1+X_2+\dots+X_n \le k)$$ so to find the distribution of $$N$$ we need to find the convolutions of the distributions of the $$X_i$$ for $$i=1,2,\dots,n$$, for all $$n$$. Those convolutions can be found numerically, but for large $$n$$ it might be much work, so we try instead to approximate the cumulative distribution function for the convolutions, using saddlepoint methods. For another example of saddlepoint methods, see my answer to Generic sum of Gamma random variables We will use the Lugannini-Rice approximation for the discrete case, and follows R Butler: "Saddlepoint Approximations with Applications", page 18 (second continuity correction). First, we need the moment generating function of the $$X_i$$, which is $$M(T) = E e^{tX_i}= \frac16 (e^t+e^{2t}+e^{3t}+e^{4t}+e^{5t}+e^{6t})$$ Then the cumulant generating function for the sum of $$n$$ independent dice becomes $$K_n(t)=n \cdot log(\frac16\sum_{i=1}^6 e^{it})$$ and we also need the first few derivatives of $$K$$, but we will find those symbolically using R. The code is the following:  DD <- function(expr, name, order = 1) { if(order < 1) stop("'order' must be >= 1") if(order == 1) D(expr, name) else DD(D(expr, name), name, order - 1) } make_cumgenfun <- function() { fun0 <- function(n, t) n*log(mean(exp((1:6)*t))) fun1 <- function(n, t) {} fun2 <- function(n, t) {} fun3 <- function(n, t) {} d1 <- DD(expression(n*log((1/6)*(exp(t)+exp(2*t)+exp(3*t)+exp(4*t)+exp(5*t)+exp(6*t)))), "t", 1) d2 <- DD(expression(n*log((1/6)*(exp(t)+exp(2*t)+exp(3*t)+exp(4*t)+exp(5*t)+exp(6*t)))), "t", 2) d3 <- DD(expression(n*log((1/6)*(exp(t)+exp(2*t)+exp(3*t)+exp(4*t)+exp(5*t)+exp(6*t)))), "t", 3) body(fun1) <- d1 body(fun2) <- d2 body(fun3) <- d3 return(list(fun0, fun1, fun2, fun3)) }  Next, we must solve the saddlepoint equation. That is done by the following code: funlist <- make_cumgenfun() # To solve the saddlepoint equation for n, k: solve_speq <- function(n, k) {# note that n+1 <= k <= 6n is needed Kd <- function(t) funlist[[2]](n, t) k <- k-0.5 uniroot(function(s) Kd(s)-k, lower=-100, upper=1, extendInt="upX")$root
}


Note that the above code is not very robust, for values of $$k$$ far in either tail of the distribution it will not work. Then some code for actually calculating the tail probability function, approximately, by the Luganini-Rice approximation, following Butler, page 18, (second continuity correction):

# Function for returning the tail probability:

#

Ghelp  <-  function(n, k) {
stilde  <-  solve_speq(n, k)
K  <-  function(t) funlist[[1]](n, t)
Kd <-  function(t) funlist[[2]](n, t)
Kdd <- function(t) funlist[[3]](n, t)
Kddd <- function(t) funlist[[4]](n, t)
w2tilde  <-  sign(stilde)*sqrt(2*(stilde*(k-0.5)-K(stilde)))
u2tilde  <-  2*sinh(stilde/2)*sqrt(Kdd(stilde))
mu  <-  Kd(0)
result  <- if (abs(mu-(k-0.5)) <= 0.001) 0.5-Kddd(0)/(6*sqrt(2*pi)*Kdd(0)^(3/2))  else
1-pnorm(w2tilde)-dnorm(w2tilde)*(1/w2tilde - 1/u2tilde)
return(result)
}
G  <- function(n, k) {
fun  <- function(k) Ghelp(n, k)
Vectorize(fun)(k)
}


Then let us try to use this to calculate a table of the distribution, based on the formula $$P(N \ge n) = P(X_1+X_2+\dots+X_n \le k) \\ = 1-P(X_1+\dots+X_n \ge k+1) \\ = 1-G(n,k+1)$$ where $$G$$ is the function fron the R code above.

Now, let us answer the original question with $$K=20$$. Then the minimum number of rolls is 4 and the maximum number of rolls is 20. The probability that 20 rolls is needed is very small, and can be calculated exactly from the binomial formula, I leave that to the reader. (the approximation above will not work for $$n=20$$).

So the probability that $$N \ge 19$$ is approximated by

> 1-G(20, 21)
[1] 2.220446e-16


The probability that $$N\ge 10$$ is approximated by:

> 1-G(10, 21)
[1] 0.002880649


And so on. Using all this, you can get an approximation for the expectation yourself. This should be much better than the approximations based on the central limit theorem.