# Expected number of rolls until a number appears $k$ times

Roll a fair die, what is the expected number of rolls until a number appear $$k$$ times? Not necessarily consecutive.

Let $$N$$ be the number of rolls until a number appear $$k$$ times. For $$k=2$$, we know that the largest possible value for $$N$$ is $$7$$. Hence we have \begin{align} &P(N=1)=0\\ &P(N=2)=1/6\\ &P(N=3)=5/6\cdot 2/6\\ &P(N=4)=5/6\cdot 4/6\cdot 3/6\\ &P(N=5)=5/6\cdot 4/6\cdot 3/6 \cdot 4/6\\ &P(N=6)=5/6\cdot 4/6\cdot 3/6 \cdot 2/6\cdot 5/6\\ &P(N=7)=5/6\cdot 4/6\cdot 3/6 \cdot 2/6\cdot 1/6 \end{align} However I do not know how to generalise it for any $$k$$. Could someone please help?

• Is this the birthday problem? en.wikipedia.org/wiki/Birthday_problem Jul 6, 2020 at 21:51
• @eric_kernfeld it is the birthday problem when $k=2$ Jul 6, 2020 at 21:56
• @eric_kernfeld Sure. In the simple case $k=2$. I only ask for the expectation of $N$. But following the birthday problem, I will have to calculate/approximate the entire distribution of $N$. I was hoping for a simpler method. This becomes more important for general $k$. Jul 6, 2020 at 22:09

It's sometimes useful to recast a problem in terms that yield better search engine results. Here is an alternative formulation of your problem:

We throw balls at random into $$n=6$$ urns, with equal probability. How many balls do we expect to throw until one urn contains $$k$$ balls?

And there is actually a closed form solution to this question at Balls are placed into 3 urns. Expected time until some urn has 100 balls. Namely,

$$n\int_0^\infty\bigg(\frac{\Gamma(k,a)}{\Gamma(k)}\bigg)^n\,da = \frac{n}{(k-1)!^n}\int_0^\infty\Gamma(k,a)^n\,da.$$

The calculation at the link works just as well for an $$n=6$$-sided die as for $$d=3$$ urns, in particular using the relevant property of the upper incomplete gamma function.

You can evaluate this improper integral numerically (like this for $$k=2$$ at WolframAlpha) or use it for subsequent analyses as-is. The numerical evaluation (by WolframAlpha as above) is reassuringly close to simulation results for $$n=6$$ and $$1\leq k\leq 10$$:

 k  Numerical  Simulation
----------------------
1   1          1
2   3.77469    3.77777
3   7.29554    7.29863
4  11.2138    11.21731
5  15.3858    15.37895
6  19.7374    19.75814
7  24.2245    24.23791
8  28.8185    28.79771
9  33.4995    33.48532
10  38.2533    38.21238


Simulation R code:

n_sides <- 6
kk_max <- 10
expectation_sim <- structure(rep(0,kk_max),.Names=1:kk_max)
n_sims <- 1e5
pb <- winProgressBar(max=kk_max)
for ( kk in 1:kk_max ) {
setWinProgressBar(pb,kk,paste(kk,"of",kk_max))
for ( ii in 1:n_sims ) {
state <- rep(0,n_sides)
counter <- 0
while ( all(state<kk) ) {
roll <- sample(1:n_sides,1)
state[roll] <- state[roll]+1
counter <- counter+1
}
expectation_sim[kk] <- expectation_sim[kk]+counter
}
}
close(pb)
expectation_sim <- expectation_sim/n_sims
expectation_sim


We can model your problem as an absorbing Markov Chain. The possible states are $$n$$-tuples of numbers between $$0$$ and $$k$$,

$$\mathcal{S} := \{0,\dots,k\}^n,$$

each state marking how often each number between $$1$$ and $$n$$ has already come up. (Of course, $$n=6$$.) The transient states are those where all entries are strictly smaller than $$k$$,

$$\mathcal{S}_t := \{0,\dots,k-1\}^n\subset\mathcal{S},$$

and the absorbing states are those with at least one entry equal to $$k$$,

$$\mathcal{S_a} := \{s\in\mathcal{S}\,|\,\exists i\colon s_i=k\}=\mathcal{S}\setminus\mathcal{S}_t.$$

We start in the state $$(\underbrace{0,\dots,0}_{n \text{ times}})$$.

In principle, it's easy to set up the transition matrix $$P$$, but it's painful. There are $$(k+1)^n$$ states, which already for $$n=6$$ and $$k=2$$ is $$3^6=729$$. There are various orderings possible on $$\mathcal{S}$$, but none that appear to make the transition matrix $$P$$ very easy to work with abstractly. However, it should not be too hard to set $$P$$ up for a concrete (small) choice of $$n$$ and $$k$$. (I'll admit that I didn't manage to make my calculations match simulations. It's rather late here.)

However, once we do have $$P$$, we can use a standard result on the expected number of steps to reach an absorbing state. Namely, we can reorder the states with the absorbing ones at the end and express $$P$$ in block diagonal form,

$$P = \begin{pmatrix} Q & R \\ 0 & I_{|\mathcal{S}_a|} \end{pmatrix},$$

where $$Q$$ corresponds to transition probabilities between transient states only, $$R$$ to transition probabilities from transient to absorbing states, and $$I_{|\mathcal{S}_a|}$$ is an identity matrix (of size equal to the number of absorbing states $$|\mathcal{S}_a|$$).

Now, let $$N:=(I_{|\mathcal{S}_t|}-Q)^{-1}$$, and multiply $$N$$ by a vector of ones, $$N1$$. The $$i$$-th entry of this vector gives the expected number of steps until we reach an absorbing state when starting from the $$i$$-th state. So we can just read off the entry in this vector that corresponds to our starting state.

So, no formula, and unfortunately, I didn't get my little program to give me results that matched a quick simulation. However, you may be able to write your own program, or looking at the Markov chain literature may be helpful. (Note that $$\mathcal{S}$$ is a kind of $$n$$-dimensional discrete cube, which may also be helpful in searching.)