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
    If I understand correctly, we stop rolling the dice the moment that the sequence of all rolls contains a sequence of the form $1^a 2^b 3^c 4^d 5^e 6$, where $a$ through $e$ are positive integers, and the question is the expected number of dice we roll before stopping. Is that right?2012-03-18
  • 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
    Any help getting mathjax to understand tabulars? It's not math mode, so how do I signal non-math TeX?2012-03-18
  • 0
    Try "array" instead. (Why should it not be math?)2012-03-18
  • 0
    I wanted the top row and the first column to give headers about what is in those rows/columns. It makes it easier to interpret than just putting a matrix.2012-03-18
  • 1
    Just replacing the "tabular" by "array" in your original version gives you a table with lines between rows/columns like you wanted.2012-03-18
  • 0
    I'm not sure that the transition matrix is right. When you have some numbers in the sequence, rolling a 1 breaks the sequence but does not send you back to the start, just back to 1.2012-03-18
  • 0
    Ah yes, you're right about that. If you roll a 1 at any stage, you go back to the 1 state. I'll modify the probabilities.2012-03-18
  • 0
    I have an answer but cannot get my simulation to agree with it. If I overlook the return to one in both simulation and analytics then I get good agreement.2012-03-18
  • 0
    @EMS I'm afraid I don't see how to get from the transition matrix to the expected number of rolls. Can you explain this?2012-03-19
  • 0
    The idea is that by making the transition matrix entry in the bottom left equal to 1.0, we make the chain as a whole recurrent. In the long run, the chain only sits in the 'accept' state for 1 "time unit" before flipping back to the start. Thus, in the long run, the probability of seeing the chain in that accept state should be equal to the reciprocal of the expected number of rolls to get there. So then you need to use either the system of equations method or the eigenvalue/eigenvector method with the transition matrix to get that long-term probability.2012-03-19
  • 0
    Ahh, also a friend just pointed out that I had an error. The last row should really be $[5/6, 1/6, 0, 0, 0, 0, 0]$, because you don't automatically go back to start if the next play of the game began with a roll of a $1$. That should fix issues with the eigenvector of my transition matrix. I've edited the matrix to reflect that.2012-03-19
  • 0
    I've added some Python code that verifies this with eigenvectors and Monte Carlo, and added a bit of discussion. Hopefully it clears things up now.2012-03-19
  • 0
    Excellent! 22500 is in good agreement 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.

  • 6
    Would the downvoter care to comment on why?2012-03-19
  • 0
    A neat approach. I'm curious as to why your intermediate answers differ from EMS's? $6,30,150,750,3750,22500$ vs $6,36,180,900,4500,22500$2012-03-20
  • 0
    I noticed that EMS's intermediates seem to be the differences of mine. I have not tried to analyze EMS's answer yet.2012-03-20
  • 0
    I fixed the typo in my algebra and I get 6,36,180,900,4500,22500 as the intermediate times. I'll post my solution once I have a bit more time. @robjohn How did you create your diagram?2012-03-21
  • 0
    @jcoe: I use [Intaglio](http://purgatorydesign.com/Intaglio/) for the Mac. There is no Windows port so far.2012-03-21
  • 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$