Modeling a probabilistic game with changing probabilities

Question

I am trying to model a simple game (a more complex problem I am working on reduces to it).

1. You roll an $n$-sided dice labeled 1 to $n$. (This can be a normal 6-sided dice if you'd like.)
2. If the dice shows 1, you draw a card from a deck of $m$ cards labeled 1 to $m$. (This can be a normal 52-card deck if you'd like.) You record which card you drew, replace it in the deck, and shuffle.
3. This process repeats until you draw a card you have previously seen.

The question: what is the distribution of the expected length of this game, in terms of the number of dice rolls?

---

I am able to model the individual sub-games without an issue:

Game 1, the dice game, is geometrically distributed with $p=1/n$. On average, we expect to see a 1 every $n$ throws.

Game 2, the card game, is a bit more complex, but we can model it like this:

• On our $(i+1)$-th draw, we have already seen $i$ cards, all of which are unique (if not, the game would already be over).
• Therefore, the probability the given card $i+1$ matches one of the previous cards is $\frac{i}{m}$ for a deck of size $m$.
• Factoring in the probability of not having already won, the distribution of winning the game on draw $i$ is

$$ P(i) = \frac{m!/(m-i)!}{m^i} \cdot \frac{i}{m} $$

(The numerator of the left fraction is just the number of permutations of $i$ items from $m$ items.)

I can't find a closed form for the expected value of this PDF, but it matches simulations given a value of $m$. For example, consider the standard case of $n=6$ and $m=52$.

• We expect a 1 on the dice every 6 rolls.
• We expect to draw a duplicate card after 9.7 draws, on average. (Solved numerically in Mathematica and confirmed in simulation.)
• So, we expect the game to last $\approx 58$ rolls on average, since we can multiply expectations of independent events.

But what is the actual distribution of game lengths?

Here's some simulation data. I play the game 100,000 times:

The simulation code looks like this:

def play_game():
    cards_seen = set()
    for i in it.count(1):
        dice = random.randint(1, 6)
        if dice != 1:
            continue

        card = random.randint(1, 52)
        if card in cards_seen:
            return i, 1 + len(cards_seen)

        cards_seen.add(card)

Statistics from the simulation:

• Mean number of turns: 58.19
• ...median: 54
• Mean number of cards drawn: 9.69
• ...median: 9
• Implied number of turns/per card (i.e., dice rolls 1): 6.005

This looks very familiar (i.e., vaguely negative-binomial), but the changing probability of pulling a matching card makes such an analysis tricky.

I know the combination of distributions can be computed via convolution, but I'm not clear how I would apply that here, since these two distributions are not really independent — whether we play the card subgame is directly dependent on the results of the dice subgame.

I believe this game can be modeled as a Markov chain, but I'm not sure if that helps with the analysis.

• We have dice states $D_k$ and card states $C_k$.
• We start at state $D_1$ and move to state $C_1$ with probability $1/6$, or remain at state $D_1$ with probability $5/6$.
• At some state $C_k$, we end the game with probability $P(k)$ from above, and move to state $D_{k+1}$ with probability $1-P(k)$.

Now the question reduces to finding the distribution of the number of steps in this chain.

Answer 1

Lets start with deck game only.

Call $T$ the random time (# of deck extraction) where we see the first duplicate.
We can use this formula for expectation: $$E[T]=\sum_{k=0}^{m}P(T>k)$$ and observe that $P(T>k)=\frac{m!}{(m-k)!}\frac{1}{m^k}$, as it requires that in k extractions, we get k different cards.

Adding the die-part, will just delay the inevitable by a factor of $n$ on average.

$$$$

For the distribution part lets introduce $S$ for the random time of the overall game and $X_k$ for the numbers of 1s in the first k rolls (it is a Binomial).

So we have that: $$ P(S=k)=\sum_{p=0}^{\min\{k-1, m\}} P(S=k \ \bigcap \ X_{k-1}=p) $$

A bit of intuition. If S=k, it means that the k-th iteration gets a 1 with a card already seen, while all previous extractions were unique.
But how many extractions? Lets cycle through them leveraging $X_k-1$.

Finally we get the following (unsimplified) expression: $$ P(S=k \ \bigcap \ X_{k-1}=p) = \binom{k-1}{p}\ \left(\frac{1}{n}\right)^p\ \left(1-\frac{1}{n}\right)^ {(k-1)-p}\ \frac{1}{n}\ \ \frac{1}{m^p}\ \ \binom{m}{p} \ p\ !\ \ \frac{p}{m} $$

$$$$ I leave some R code that I have used to view results.

prob_k_p <- function(k, p, n=6, m=52){
  choose(k-1, p) * (n^(-p)) * (((n-1)/n)^(k-1-p)) / n *
  m^(-p) * choose(m, p) * factorial(p) * p/m
}

prob_k <- function(k, n=6, m=52){
  temp.1 <- sapply(0:min(k-1,m), prob_k_p, k=k, n=n, m=m)
  return(sum(temp.1))
}


n<- 6; m <- 52

n * sum(factorial(m)/factorial(m-(0:m))*m^(-(0:m)))
plot(2:200, sapply(2:200, prob_k, n=n, m=m))

Answer 2

You can model the general case as a Markov chain whose state at any time is the number of distinct cards that have been observed up to that time or (conveniently) $\ m+1\ $ if a duplicate card has been observed. This chain has $\ m+2\ $ states, $\ 0,1,\dots,m+1\ $, with transition matrix $\ P\ $ given by $$ P(S_{t+1}=x\,|\,S_t=y)=\cases{\frac{n-1}{n}&if $\ 0\le x=y\le m$\\ \frac{m-y}{mn}&if $\ 0\le x=y+1\le m$\\ \frac{y}{mn}&if $\ x=m+1\ \text{ and }\ y\le m$\\ 1&if $\ x=y=m+1\ $\\ 0&otherwise.} $$ If $\ e_k\ $ is the expected number of further dice rolls until the game terminates given that the current state is $\ k\ $, then $\ e_k\ $ must satisfy the first order difference equation $$ e_k=1+\left(\frac{n-1}{n}\right)e_k+\left(\frac{m-k}{mn}\right)e_{k+1} $$ with terminal condition $\ e_{m+1}=0\ $. This has solution $$ e_t=n\sum_{k=t}^m\frac{(m-t)!}{(m-k)!\,m^{k-t}}\ , $$ and, in particular, $$ e_0=n\sum_{k=0}^m~\frac{m!}{(m-k)!\,m^k}\ . $$ For $\ n=6,m=52\ $ this gives $$ e_0\approx58.31\ , $$ correct to two decimal places.

Correction

My original expression for $\ e_t\ $ was missing a factor $\ \frac{m^t(m-t)!}{m!}\ $, which I have now included. Since the value of this factor is $\ 1\ $at $\ t=0\ $, the correction doesn't affect the value given for $\ e_0\ $.

Answer 3

This example maybe gives some insight how to calculate it. Using the same Markov chain as lonza leggiera, we can calculate the Jordan form of the transition matrix (the non-absorbing part). Unfortunately there are unit-blocks, when you exponentiate you get binomial coefficients.

$$ \left(\begin{array}{rrrrrr} 25 & 5 & 0 & 0 & 0 & 0 \\ 0 & 25 & 4 & 0 & 0 & 0 \\ 0 & 0 & 25 & 3 & 0 & 0 \\ 0 & 0 & 0 & 25 & 2 & 0 \\ 0 & 0 & 0 & 0 & 25 & 1 \\ 0 & 0 & 0 & 0 & 0 & 25 \end{array}\right) = \\ \left(\begin{array}{rrrrrr} 120 & 0 & 0 & 0 & 0 & 0 \\ 0 & 24 & 0 & 0 & 0 & 0 \\ 0 & 0 & 6 & 0 & 0 & 0 \\ 0 & 0 & 0 & 2 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{array}\right) \left(\begin{array}{rrrrrr} 25 & 1 & 0 & 0 & 0 & 0 \\ 0 & 25 & 1 & 0 & 0 & 0 \\ 0 & 0 & 25 & 1 & 0 & 0 \\ 0 & 0 & 0 & 25 & 1 & 0 \\ 0 & 0 & 0 & 0 & 25 & 1 \\ 0 & 0 & 0 & 0 & 0 & 25 \end{array}\right) \left(\begin{array}{rrrrrr} \frac{1}{120} & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{1}{24} & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{1}{6} & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{2} & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{array}\right) $$

Here's the code (It has some remnants of my previous method where I made an absorbing state for each amount of cards seen before being absorbed because I thought that might be helpful):

def f(n, m):
    #state = (distinct cards seen, seen same), seen sames are absorbing
    states = sorted([(k,0) for k in range(m+1)], key=lambda t:t[1]) + [(m+1,1)]
    #old way
    #sorted([(k,w) for k in range(m+1) for w in range(2)], key=lambda t:t[1])
    sToI = {s:i for i,s in enumerate(states)}
    mat = matrix(QQ, len(states))
    for i,s in enumerate(states):
        if s[1]==1:
            mat.add_to_entry(i, i, 1)
        else:
            mat.add_to_entry(i, i, (n-1)/n)
            if s[0]<m: mat.add_to_entry(i, sToI[(s[0]+1, 0)], (m-s[0])/(n*m))
            mat.add_to_entry(i, sToI[(m+1, 1)], s[0]/(n*m))
    print(states)
    show(1/(m*n), m*n*mat)
    D, S = (m*n*mat[:m+1, :m+1]).jordan_form(transformation=True)
    show(S * D * S^(-1))
    show(S, D, S^(-1))
    fundMat = (matrix.identity(m+1)-mat[:m+1, :m+1])^(-1)
    R = mat[:m+1, m+1:]
    #show(fundMat, R)
    pMat = fundMat*R
    #show(pMat)
    #return tuple(pMat[0])

    
    
print (f(6, 5))

I added the computer calculation of the probabilities (and cleaned up the code a bit):

def cumuF(n, m, uptoK=100):
    mat = matrix(QQ, m+2)
    for j in range(m+1):
        mat.add_to_entry(j, j, (n-1)/n)
        if j<m: mat.add_to_entry(j, j+1, (m-j)/(n*m))
        mat.add_to_entry(j, m+1, j/(n*m))
    mat.add_to_entry(m+1, m+1, 1)
    ret = []
    A = mat
    for k in range(uptoK):
        ret.append(A[0][-1])
        A *= mat
    return ret

    
F = cumuF(6, 52, 200)
ps = [F[i]-F[i-1] for i in range(1, len(F))]

show(bar_chart(ps))

Answer 4

Let's say one step in the process consists of either (1) rolling a number greater than $1$, with probability $1-1/n$; or (2) rolling a $1$ followed by drawing card $i$, for $1 \le i \le m$, with probability $1/(mn)$. Define $T$ to be the number of the step on which we first draw a duplicate card, and define $p_k$ to be the probability that we have not drawn any duplicates by the $k$th step; equivalently, $p_k = P(T > k)$. By a well-known theorem, $$E(T) = \sum_{k=0}^{\infty} P(T>k) = \sum_{k=0}^{\infty} p_k$$

Now let $f(x)$ be the exponential generating function of $p_k$, i.e. $$f(x) = \sum_{k=0}^{\infty} \frac{p_k}{k!} x^k$$ In a sequence of $k$ steps with no duplicates, we can have any number of occurrences of rolling a number greater than $1$, and any card drawn either zero or one times (combined with rolling a $1$). So $$f(x) = e^{(1-1/n)x}\left(1 + \frac{x}{mn} \right)^m$$ Taking advantage of the indentity $$\int_0^{\infty} x^k e^{-x} \; dx = k!$$ we have $$\begin{align} E(T) &= \sum_{k=0}^{\infty} p_k \\ &= \int_0^{\infty} e^{-x} f(x) \; dx \\ &= \int_0^{\infty} e^{-x} \cdot e^{(1-1/n)x}\left(1 + \frac{x}{mn} \right)^m \; dx \\ &= \int_0^{\infty} e^{-x/n}\left(1 + \frac{x}{mn} \right)^m \; dx \end{align}$$ In the case $n=6$, $m=52$, Mathematica evaluates the integral as $58.3102$.