Markov Chains: Dice Problem I'm not sure how to start

Question

The Problem: You start with five dice. Roll all the dice and put aside those dice that come up 6. Then roll the remaining dice, putting aside those dice that come up 6. And so on. Let $X_n$ be the number of dice that are sixes after $n$ rolls.

a) Describe the transition matrix $P$ for this Markov Chain. I'm just trying to figure out how to construct this transition matrix..

At first I thought I can just have something like (I just put the first row for simplicity since I realized this was wrong)

$$ \begin{bmatrix} 1 & 1/4 &1/4 & 1/4 &1/4 \\ \end{bmatrix} $$

I realize this is wrong since the sum exceeds 1. I also realize that I was trying to construct a 5x5 matrix but I would just have the probability that the other dice come up six after rolling one individual dice.

Can someone give me some advice on how to start this? Thank you in advance.

Edit: The only thing I've thought of is creating a massive matrix that has a row for the probability of a success with dice 1, 2, 3, 4, 5, 1&2, 1&3, 1&4 etc. etc. I.e. each row represents a different possibility. But this seems a bit too large..

Edit2: I think I realize that I can have it be a 6x6 matrix. So the first row will be conditioned on not getting any sixes the first roll.

Answer 1

The Markov chain determines $X_n$ which takes values in $S = \{0, 1, 2, 3, 4, 5\}$. The $6 \times 6$ transition matrix $P$ will determine the probability mass over $S$ for $X_{n}$ when $X_{n-1}$ takes a certain value. That is, you need the matrix $P$ that defines for all $x \in S$ and all $i \in S$, $$\Pr(X_n = i \mid X_{n-1} = x)\,. $$

Note that since $X_n$ counts the number of 6s till time $n$, it cannot decrease. So in the transition matrix, the probabilities are zero for whenever $i$ is less than $x$. So let's start filling the transition matrix, starting from the bottom row. When $x = 5$, that is $X_{n-1} = 5$, then the next state must be 5 because all dice are out of play, so all the mass is on state $5$. $$ P = \left[ \begin{array}{cccccc} . & . & . &. &.&.\\ . & . & . &. &.&.\\ . & . & . &. &.&.\\ . & . & . &. &.&.\\ . & . & . &. &.&.\\ 0 & 0 & 0 &0 &0&1\\ \end{array} \right]\,. $$

Now, lets fill the 5th row (that is $X_{n-1} = 4$. Since the number of 6s so far is four, that means, that the number of dice still in play is 1. When I roll do the experiment this $n$th time, I am only rolling the one die. So $\Pr(X_n = 5 \mid X_{n-1} = 4) = \Pr(\text{rolling a 6 for a single die}) = 1/6$. The rest of the probability goes to the state $4.$ So,

$$ P = \left[ \begin{array}{cccccc} . & . & . &. &.&.\\ . & . & . &. &.&.\\ . & . & . &. &.&.\\ . & . & . &. &.&.\\ 0 & 0 & 0 & 0 & \frac{5}{6}& \frac{1}{6}\\ 0 & 0 & 0 &0 &0&1\\ \end{array} \right]\,. $$

I'll fill one more row for you. for when $X_{n-1} = 3$, then there are two dice to roll. For $X_n = 5$, both dice should be 6, which happens with probability $1/36$. For $X_n = 4$, exactly one die has to be a 6, which happens with probability $2 \times \frac{1}{6} \times \frac{5}{6} = \frac{10}{36}$. The rest of the probability goes to the state $X_n = 3$.

$$ P = \left[ \begin{array}{cccccc} . & . & . &. &.&.\\ . & . & . &. &.&.\\ . & . & . &. &.&.\\ 0 & 0 & 0 & \frac{25}{36} & \frac{10}{36}& \frac{1}{36}\\ 0 & 0 & 0 & 0 & \frac{5}{6}& \frac{1}{6}\\ 0 & 0 & 0 &0 &0&1\\ \end{array} \right]\,. $$

Since this is self-study, I will let you fill in the rest of the rows.

EDIT Seeing Chaconne's comment, here I check the final matrix $P$ obtained by looking at $P^t$ at $t = 1000$. I know that eventually all the mass should be at the state $X_n = 5$ irrespective of the starting value. The output matches this result.

> Pt <- P
> for(t in 2:1000)
+ {
+   Pt <- Pt%*%P
+ }
> Pt
     [,1]          [,2]          [,3]          [,4]         [,5] [,6]
[1,]    0 9.418588e-317 2.859314e-237 4.340182e-158 3.294003e-79    1
[2,]    0 1.883718e-317 1.143726e-237 2.604109e-158 2.635202e-79    1
[3,]    0  0.000000e+00 2.859314e-238 1.302054e-158 1.976402e-79    1
[4,]    0  0.000000e+00  0.000000e+00 4.340182e-159 1.317601e-79    1
[5,]    0  0.000000e+00  0.000000e+00  0.000000e+00 6.588005e-80    1
[6,]    0  0.000000e+00  0.000000e+00  0.000000e+00 0.000000e+00    1

Answer 2

If you follow through the examples in previous answer, you may notice that the probability of going from $n$ dice to $m$ dice in one roll is the number of ways you can choose $m$ of $n$ dice, multiplied by $p^mq^{n-m}$, where $p=5/6$ and $q=1/6$. This suggests looking at coefficients of $x^m$ in the polynomial $(p+q\cdot x)^m$, which may allow you to more easily set up a Markov transition probability matrix.

For further explanation of this idea, see Probability Generating Functions at Wikipedia; in particular, see the second item in the Examples section, which treats the probability generating function of a binomial random variable.

In Python you can use SymPy, and in Julia can use SymPy or the Polynomials package, to compute powers of polynomials. In the Julia program below the four lines that compute v, x, vv, P compute the Markov matrix for a given number of dice. As an example, here's the program's output when it's run with a parameter of 4 on the command line:

 P = Markov transition probabilities
 0.3164  0.4219  0.2109  0.0469  0.0039
 0.0000  0.4219  0.4219  0.1406  0.0156
 0.0000  0.0000  0.5625  0.3750  0.0625
 0.0000  0.0000  0.0000  0.7500  0.2500
 0.0000  0.0000  0.0000  0.0000  1.0000

 N = inv(I-Q) = expected durations in states
 1.4629  1.0675  1.7347  3.4768
 0.0000  1.7297  1.6680  3.4749
 0.0000  0.0000  2.2857  3.4286
 0.0000  0.0000  0.0000  4.0000
Expected length from S1 to Sz is 7.7418   nd=4

As you can see in the output, the program computes the Markov matrix and also the expected number of times the Markov process visits each node before falling into the absorption case of zero dice, when starting in any given state. (See Theorem 1 in Appendix D of M. M. Conroy's A Collection of Dice Problems with solutions and useful appendices regarding this computation.) The expected number of steps to get from 4 dice to 0 dice is the sum of the first row of N, or 7.7418.

#!/usr/bin/env julia
# -*- mode: julia; -*- 
# ~ jiw ~  25 Aug 2022 ~  ~

# Re: https://stats.stackexchange.com/questions/327141/
# markov-chains-dice-problem ...  "start with five dice. Roll all the
# dice and put aside those dice that come up 6. Then roll the
# remaining dice, putting aside those dice that come up 6" ...

# Method: Use theorem 1 in Appendix D of M.M.Conroy,
# <https://www.madandmoonly.com/doctormatt/mathematics/dice1.pdf>

module expPathLen
using Formatting, Polynomials, LinearAlgebra

function printMat(M, label)
    println(label)
    for i in 1:size(M,1)
        println(join(map(x->fmt("7.4f", x), M[i,:])," "))
    end
end

function pathLen(nd)
    # Build Markov transition matrix for given number of dice.  See eg
    # <https://en.wikipedia.org/wiki/Probability-generating_function#Examples>
    # where second example is generating function of a binomial RV
    v = Polynomial([(nd-1)/nd,1/nd])   # single-die probabilities
    x = Polynomial([ 0, 1])            # polynomial f(x)=x
    vv =  [ ((x^j)*(v^(nd-j))).coeffs for j in 0:nd ]
    P = reduce(vcat,transpose.(vv)) # Reshape vector of vectors → matrix

    # Find expected path length from first state (nd dice) to last (0 dice)
    Q = P[1:nd, 1:nd]         # Get reg. part of transition matrix
    N = inv(I-Q)              # Solve for expected path lengths
    printMat(P, "\n P = Markov transition probabilities")
    printMat(N, "\n N = inv(I-Q) = expected durations in states")
    # Find sum of row 1 (expected encounter counts, starting at state 1)
    expectL = sum(N[1,:]) # expectL = time to fall into absorbing state
    printfmtln("Expected length from S1 to Sz is {:5.4f}   nd={}",
                       expectL, nd)
end
ndice = length(ARGS)>=1 ? parse(Int64, ARGS[1]) : 5 # Get param
pathLen(ndice)
end # module

^{(I specified <!-- language: python --> before this code, as <!-- language: julia --> apparently had no effect.)}