Markov chains for collector problems

Question

Consider the following problem, from Tijms's Understanding Probability:

The Bubble Company offers a picture of one of 25 popstars in a pack of chewing gum. John and Peter each buy one pack every week. They pool the pictures of the popstars. Assuming equal chances of getting any of the 25 pictures with one purchase, denote by the random variable $N$ the number of weeks until John and Peter have collected two complete sets of 25 pictures. Calculate the expected value of N and calculate the probability $P (N > n)$ for $n = 50$, 75, 100, 125, and 150.

This problem should be solvable with an absorbing Markov chain.

I define the system to be in the state $(i,j)$ when $j$ popstars have appeared two or more times so far, and $i-j$ popstars have appeared exactly one time. Clearly, $1\le i\le 25$ and $0\le j\le i$.

The state $(25,25)$ is absorbing, since it corresponds to a completed collection.

The transition probabilities for the non-absorbing states are: $$ \begin{gathered} p_{(i,j),(i+1,j)} = \frac{25-i}{25} \\ p_{(i,j),(i,j+1)} = \frac{i-j}{25} \text{ for } j+1\le i\\ p_{(i,j),(i,j)} = \frac{j}{25} \text{ for } (i,j)\ne(25,25). \end{gathered} $$

To solve a problem with some matrix computations, we need a way to index the states, and a possible choice I guess is: $$ \text{index}(i,j) = i + 25j - \frac12 (j-1)(j-2). $$ Then, to answer e.g. the question on the probability $P(N>50)$ I should compute the 99-th power (because there are two draws for every week, and we start with state $(1,0)$ which corresponds with a single draw), but this seems quite complicated for a textbook problem. Also, the matrix is quite large so computing all these matrix products seems expensive.

Am I missing something?

Answer 1

The "Double Dixie-Cup" paper that I mentioned in a comment proves that the expected number of packs needed to get $m$ full sets of $n$ pictures is $$n\log{n}+(m-1)\log\log{n} +O(n)\text{ as } n\to\infty.$$

That is the first full set "costs" $n\log{n}$, but each additional set costs only $\log\log{n}$.

This gives $25\log{25}+\log\log{25} \approx 81.641$, and since they buy two packs a week, about $41$ weeks.

Since $25$ is a long way from $\infty$, I did a simulation to see what actually happens. In my model, they make two piles, and the first time they get a picture, they put it in pile A, the second time in pile B. The third and later times are ignored.

from random import choice
from collections import defaultdict

tests = 100000
population = range(25)

weeks = defaultdict(int)
for _ in range(tests):
    A = 25*[0]
    B = 25*[0]
    week = 0
    while sum(B)<25:
        for i in range(2):
            s = choice(population)
            if A[s] == 0:
                A[s] = 1
            elif B[s] == 0:
                B[s] = 1
        week +=1
    weeks[week] += 1

expect = sum(k*weeks[k] for k in weeks)/tests
print(tests, 'trials')
print('Mean number of weeks', expect)
for n in range(50, 175, 25):
    p = sum(weeks[k] for k in weeks if k > n)/tests
    print('Prob of more than %d weeks: %s'%(n,p))

This produced the output

100000 trials
Mean number of weeks 71.35634
Prob of more than 50 weeks: 0.92788
Prob of more than 75 weeks: 0.3365
Prob of more than 100 weeks: 0.06485
Prob of more than 125 weeks: 0.01089
Prob of more than 150 weeks: 0.0018 

which is a waiting time about $12.5\%$ shorter than the estimate for large $n$.

Though I haven't read the book, I would guess the author has the theoretical result in mind, but if I'm not mistaken, that actually gives a very poor answer. Perhaps I do him an injustice, and he means you to do the simulation. Although the Markov chain method is certainly correct, as you say, it does seem like a lot of calculation. Still, I may do it, if only to convince myself that I haven't made an error in the simulation.

EDIT

I ran the Markov chain calculations today, and they agree very closely with the simulation.

import numpy as np

def idx(n,m):
    return n*(n+1)//2+m
states = idx(25,25)+1
P = np.zeros((states, states))
for n in range(26):
    for m in range(n+1):
        s = idx(n,m)       # start state
        sa = idx(n+1, m) # add to pile A
        sb = idx(n, m+1) # add to pile B
        if n==m==25:
            P[s,s] = 1
        elif n ==25:
            P[s,sb] = (25-m)/25
            P[s,s] = m/25
        elif n == m:
            P[s, sa] = (25-n)/25
            P[s,s] = n/25
        else:
            P[s,sa] = (25-n)/25  # card not in pile A
            P[s,sb] = (n-m)/25   # card in pile A, not in pile B
            P[s,s] = m/25           # card in both piles
P=P@P #two cards per week
N=np.linalg.inv(np.eye(states-1)-P[:-1, :-1])
expect= ([email protected](states-1))[0]
print('Mean number of weeks', expect)        
P25 = np.linalg.matrix_power(P, 25)
start = np.zeros(states)
start[0] = 1

for n in range(50, 175, 25):
    m = n//25
    Q = np.linalg.matrix_power(P25,m)
    done  = (start@Q)[states-1]
    print('Prob of more than %d weeks: %s'%(n,1-done))

produced the output

Mean number of weeks 71.43383331673435
Prob of more than 50 weeks: 0.9288061487085777
Prob of more than 75 weeks: 0.33953755904361504
Prob of more than 100 weeks: 0.06481946982508657
Prob of more than 125 weeks: 0.010513629113593548
Prob of more than 150 weeks: 0.0016196726615316237