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

Question

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?

Answer 1

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

Answer 2

This is not a full answer, but it may be helpful.

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