7
$\begingroup$

Similar to: "What is the expected number of dice one needs to roll to get 1,2,3,4,5,6 in order?" but we allow repeats so 1,1,2,2,3,4,4,4,4,5,5,6 would count.

My answer (or simulation) is flawed as I cannot get reasonable agreement.

My C++11 simulation code is below. I fear that the 'bug' is more likely to be in my algebra

#include  #include   int main(int argc, char* argv[]) {   int seed = 101;   if ( argc>1 )     seed = atoi(argv[1]);    std::uniform_int_distribution distribution(1,6);   std::mt19937 engine(seed);   auto generator = std::bind(distribution,engine);    int rollForSequenceSum = 0;   int multiples = 1000;   for ( int i=0; i
  • 0
    yes that's correct.2012-03-18

3 Answers 3

11

This is not too bad if you set it up as a discrete Markov chain. Let the different possible states of the game be "S" (for start), and then 1, 2, 3, 4, 5, "6/A" (for accept). Start means you failed to get the monotone run of numbers and had to start over from scratch at the beginning of the game. Accept means you got to a 6 only after a monotone sequence of the other numbers.

You want to know the long-term probability of being in state 6/A, and from this you can easily figure out what the expected number of die rolls would be (I'll leave that to you because it's satisfying.)

But to get you started, you can write down a transition matrix for this game:

$\begin{array}{| c || c | c | c | c | c | c | c |} \hline & S & 1 & 2 & 3 & 4 & 5 & 6/A \\ \hline S & 5/6 & 1/6 & 0 & 0 & 0 & 0 & 0 \\ \hline 1 & 4/6 & 1/6 & 1/6 & 0 & 0 & 0 & 0 \\ \hline 2 & 1/2 & 1/6 & 1/6 & 1/6 & 0 & 0 & 0 \\ \hline 3 & 1/2 & 1/6 & 0 & 1/6 & 1/6 & 0 & 0 \\ \hline 4 & 1/2 & 1/6 & 0 & 0 & 1/6 & 1/6 & 0 \\ \hline 5 & 1/2 & 1/6 & 0 & 0 & 0 & 1/6 & 1/6 \\ \hline 6/A & 5/6 & 1/6 & 0 & 0 & 0 & 0 & 0 \\ \hline \end{array}$

Now it's a matter of getting the right eigenvector for this matrix, to tell you the long-term probabilities, and then figuring out how to turn that long term probability into an expected number of rolls. Note that most numerical eigenvector routines only compute "right eigenvectors", which occur on the right side of the matrix in question. For a typical stochastic matrix, you want to compute the left eigenvector corresponding to the eigenvalue of 1. So to get this numerically, you first need to take the transpose of the transition matrix and then feed that to your linear algebra package.

Added

After fixing the typos I had made in the above transition matrix (which is correct now), I used the following Python code to check what the eigenvector method yielded as a solution, and I translated your Monte Carlo simulation into Python to check numerically.

The result I got for the long-term stationary distribution of the game is $[ 0.79168, 0.16667, 0.03333, 0.00667, 0.00133, 0.000267, 0.0000444],$

and so the reciprocals of these elements gives the expected number of rolls to reach each state: $[1.26, 6.0, 30, 150, 750, 3750, 22500]$ (there could be some roundoff error depending on how good NumPy's linalg package is; I'm not the one to ask about that). Note that the second entry, 6 rolls to get the first 1, fits in with intuition for a uniform 6-sided die, so that's comforting.

So according to this, the "analytical" solution is that it should take about 22500 rolls to reach a monotone sequence ending with a 6. My Monte Carlo simulations agree with this, as I obtained an expected number of rolls of 22609.618 in a simulation of 10,000 flips (somewhat slow in Python but shouldn't be an issue in C++).

Here is the code for my solution:

import numpy as np  # Set up the transition matrix. tmat = np.asarray([ [5.0/6.0, 1.0/6.0, 0, 0, 0, 0, 0],                      [4.0/6.0, 1.0/6.0, 1.0/6.0, 0, 0, 0, 0],                     [3.0/6.0, 1.0/6.0, 1.0/6.0, 1.0/6.0, 0, 0, 0],                     [1.0/2.0, 1.0/6.0, 0, 1.0/6.0, 1.0/6.0, 0, 0],                     [1.0/2.0, 1.0/6.0, 0, 0, 1.0/6.0, 1.0/6.0, 0],                     [1.0/2.0, 1.0/6.0, 0, 0, 0, 1.0/6.0, 1.0/6.0],                     [5.0/6.0, 1.0/6.0, 0, 0, 0, 0, 0]]).astype('float64')   # Compute long-run probability and expected rolls # with eigenvalue method. w,v = np.linalg.eig(np.transpose(tmat))  # Get stationary distribution as eigenvector # corresponding to eigenvalue 1. # Take real part to get rid of Python's # complex notation, and divide by vector sum # just in case round off error makes it not # quite sum up to 1.0  stationary_probs = np.real(v[:,0]) print stationary_probs/sum(stationary_probs)  # Print the element-wise reciprocals which are # the expected number of rolls to each state. print 1/(stationary_probs/sum(stationary_probs))  # Monte Carlo simulation to compare results.  total_rolls = 0; num_repeats = 10000; iter_rolls = np.zeros((1,num_repeats)) for ii in range(num_repeats):      cur_roll = np.random.randint(1,7)     num_rolls = 1      if cur_roll == 1:         next_val = 2     else:         next_val = 1      while( next_val <= 6 ):         prev_roll = np.copy(cur_roll)         cur_roll = np.random.randint(1,7)         num_rolls =  num_rolls + 1          if(cur_roll == next_val):             next_val = next_val + 1         elif(cur_roll == prev_roll):             pass         elif(cur_roll == 1):             next_val = 2         else:             next_val = 1      iter_rolls[0,ii] = num_rolls     total_rolls = total_rolls + num_rolls     print "Finished iteration %d of %d"%(ii,num_repeats)  print "Monte Carlo mean: %f"%(np.mean(iter_rolls[0,:])) 

As an aside, this linked book chapter does a better job explaining what's going on than almost any other reference I was able to Google. In particular, the section where it defines the "mean passage time" and proves that these are the reciprocals of the stationary distribution entries is exactly the background information relevant for problems like these.

Wikipedia's treatment of this kind of common discrete Markov chain issue is paltry and disappointing. If you want to "pay back" this answer, consider editing the Wikipedia entry on Stochastic Matrix to include this example. The only example currently there is one with a trivial absorbing state which is totally unhelpful for someone learning, especially someone trying to learn how to model problems with Markov chains rather than to just solve pedagogical Markov chains given to you by teachers.

  • 0
    Excellent! 225$0$0 is i$n$ good agreeme$n$t with my simulation. I'll read through the notes you attached and consider updating the wikipedia page.2012-03-19
10

Consider the following diagram. The top row represents the expected number of rolls until the boxed number is rolled in sequence. The second row is the expected number of rolls until another number is rolled. $\boxed{0}$ represents the state before a $1$ is first rolled. The numbers by the arrows, are probabilities.

$\hskip{3mm}$flow diagram

Since there is a $\frac16$ chance of rolling a $1$, it takes, on average, $6$ rolls to get from $\boxed{0}$ to $\boxed{1}$.

Once a number is rolled, it will take, on average, $\frac65$ rolls to roll something else.

To get from $\boxed{1}$ to $\boxed{2}$, there is a $\frac15$ chance it will take $\frac65$ rolls. There is a $\frac45\cdot\frac15$ chance it will take $6+2\cdot\frac65$ rolls, and a $\left(\frac45\right)^2\cdot\frac15$ chance it will take $2\cdot6+3\cdot\frac65$, etc. Thus, it takes, on average, $ \begin{align} \sum_{k=0}^\infty\frac15\left(\frac45\right)^k\left(6k+\frac65(k+1)\right) &=4\cdot6+5\cdot\frac65\\ &=30\tag{1} \end{align} $ rolls to get from $\boxed{1}$ to $\boxed{2}$.

Getting from $\boxed{2}$ to $\boxed{3}$ is the common case. If a roll misses its mark, there is a $\frac14$ chance that a $1$ was rolled and we go back to $\boxed{1}$, saving an average of $6$ rolls; otherwise, we go back to $\boxed{0}$. Thus, missing the mark means an extra $ \frac14\cdot(36-6)+\frac34\cdot36=36-\frac32\tag{2} $ rolls, on average, to get back to $\boxed{2}$. Using the same reasoning as in $(1)$, to get from $\boxed{2}$ to $\boxed{3}$, takes $ \begin{align} \sum_{k=0}^\infty\frac15\left(\frac45\right)^k\left(\left(36-\frac32\right)k+\frac65(k+1)\right) &=4\cdot\left(36-\frac32\right)+5\cdot\frac65\\ &=4\cdot36-6+6\\ &=4\cdot36\tag{3} \end{align} $ Therefore, getting to $\boxed{3}$ takes $5$ times as long as getting to $\boxed{2}$.

The identical thing happens for $\boxed{3}$, $\boxed{4}$, $\boxed{5}$, and $\boxed{6}$. This leaves us with an average of $22500$ rolls until we get to $\boxed{6}$.


Matrix Representation for the Diagram

Alternatively, consider the following State Transition Matrix $ M_7=\frac{1}{6} \begin{bmatrix} 5 & 1 & 0 & 0 & 0 & 0 & 0 \\ 4 & 1 & 1 & 0 & 0 & 0 & 0 \\ 3 & 1 & 1 & 1 & 0 & 0 & 0 \\ 3 & 1 & 0 & 1 & 1 & 0 & 0 \\ 3 & 1 & 0 & 0 & 1 & 1 & 0 \\ 3 & 1 & 0 & 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 6 \end{bmatrix}\tag{4} $ and the initial State Vector $ S=\begin{bmatrix}1&0&0&0&0&0&0\end{bmatrix}\tag{5} $ The state vector $SM_7^n$ represents the expected distribution of states $\boxed{0}-\boxed{6}$ after $n$ rolls.

To compute the expected duration to roll a $6$, in sequence, let $M_6$ be $M_7$ with row $7$ set to $0$. Then $SM_6^n$ represents the same states as $SM_7^n$ except that when state $\boxed{6}$ is reached, it is vacated the next roll. Thus, element $7$ of $SM_6^n$ is the probability that a $6$ is first rolled, in sequence, on roll $n$. Therefore, element $7$ of $ \sum_{n=1}^\infty nSM_6^n\tag{6} $ is the expected duration until a $6$ is rolled, in sequence. Indeed, Mathematica says that $ SM_6(I-M_6)^{-2}=\begin{bmatrix}400687782 & 84353250 & 16871550 & 3374490 & 674934 & 134994 & \color{red}{22500}\end{bmatrix} $ We can perform the same computation with $M_k$, $M_7$ with all but the first $k$ rows set to $0$: $ \begin{align} SM_5(I-M_5)^{-2}&=\begin{bmatrix}16016832 & 3371550 & 674490 & 134934 & 26994 & \color{red}{4500} & 0\end{bmatrix}\\ SM_4(I-M_4)^{-2}&=\begin{bmatrix}639222 & 134490 & 26934 & 5394 & \color{red}{900} & 0 & 0\end{bmatrix}\\ SM_3(I-M_3)^{-2}&=\begin{bmatrix}25416 & 5334 & 1074 & \color{red}{180} & 0 & 0 & 0\end{bmatrix}\\ SM_2(I-M_2)^{-2}&=\begin{bmatrix}1014 & 210 & \color{red}{36} & 0 & 0 & 0 & 0\end{bmatrix}\\ SM_1(I-M_1)^{-2}&=\begin{bmatrix}30 & \color{red}{6} & 0 & 0 & 0 & 0 & 0\end{bmatrix}\\ SM_0(I-M_0)^{-2}&=\begin{bmatrix}\color{red}{0} & 0 & 0 & 0 & 0 & 0 & 0\end{bmatrix} \end{align} $ all of which agree with the diagram above.

  • 0
    Running my simulation with seeds 1-9 I get: 22422.4, 23750.7, 22047.1, 21315.7, 22979, 23112.5, 23550.9, 22797.9, 22965.82012-03-21
1

I'm posting my algebra and reasoning below. Thanks to robjohn and EMS for independently reaching an answer which agrees with my simulation and corrected algebra.

Let $N_i$ be the expected number of throws taken to have got $i$ in the sequence.

$N_6$ is the expected number of throws for the game to end.

$N_1 = \frac{1}{6}.1 + \frac{5}{6}(1+N_1)$

as there is a 1 in 6 chance of getting to $N_1$ in a single throw and a $\frac{5}{6}$ chance of wasting a throw and needing to throw the dice another $N_1$ times.

$N_2 = N_1 + \frac{1}{6}.1 + \frac{1}{6}(1+N_2-N_1) + \frac{4}{6}(1+N_2)$

as one must be at 1 already which takes $N_1$ throws and one has a 1 in 6 chance of rolling a 2, a 1 in 6 chance of rolling a 1 thus wasting a throw but only having to start from 1, and a 4 in 6 chance of rolling anything other than 1 or 2 in which case a throw is wasted and one has to start again from scratch.

$N_3 = N_2 +\frac{1}{6}.1 + \frac{1}{6}(1+N_3-N_1) + \frac{1}{6}(1+N_3-N_2) + \frac{3}{6}(1+N_3)$

as one must be at 2 already which takes $N_2$ throws and one has a 1 in 6 chance of rolling a 3, a 1 in 6 chance of rolling a 1 thus wasting a throw but only having to start from 1, a 1 in 6 chance of rolling a 2 thus wasting a throw but only having to start from 2, and a 3 in 6 chance of rolling anything other than 1, 2 or 3 in which case a throw is wasted and one has to start again from scratch.

Similarly

$N_4 = N_3+\frac{1}{6}.1 + \frac{1}{6}(1+N_4-N_1) + \frac{1}{6}(1+N_4-N_3) + \frac{3}{6}(1+N_4),$

$N_5 = N_4 +\frac{1}{6}.1 + \frac{1}{6}(1+N_5-N_1) + \frac{1}{6}(1+N_5-N_4) + \frac{3}{6}(1+N_5),$

$N_6 = N_5 +\frac{1}{6}.1 + \frac{1}{6}(1+N_6-N_1) + \frac{1}{6}(1+N_6-N_5) + \frac{3}{6}(1+N_6).$

These are linear simultaneous equations which one can solve to get:

$N_1 = 6$

$N_2 = 5.N_1 + 6$

$N_3 = 5.N_2$

$N_4 = 5.N_3$

$N_5 = 5.N_4$

$N_6 = 5.N_5$

or

$N_1 = 6$

$N_2 = 36$

$N_3 = 180$

$N_4 = 900$

$N_5 = 4500$

$N_6 = 22500$