Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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 <iostream>
#include <random>

int main(int argc, char* argv[])
{
  int seed = 101;
  if ( argc>1 )
    seed = atoi(argv[1]);

  std::uniform_int_distribution<int> distribution(1,6);
  std::mt19937 engine(seed);
  auto generator = std::bind(distribution,engine);

  int rollForSequenceSum = 0;
  int multiples = 1000;
  for ( int i=0; i<multiples; ++i )
  {
    int rollCount = 0;
    int nextInSequence = 1;

    while ( nextInSequence <= 6 )
    {
      ++rollCount;
      int random = generator();
      if ( random == nextInSequence )
      {
        ++nextInSequence;
      } 
      else if ( random == (nextInSequence-1) )
      {
        //Do nothing
      } 
      else if ( random == 1 )
      {
        nextInSequence = 2;
      }
      else
      {
        nextInSequence = 1;
      }
    }
    rollForSequenceSum += rollCount;
  }
  double mean = (double) rollForSequenceSum / (double) multiples;
  std::cout << mean << std::endl;
}
share|improve this question
    
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? –  Tanner Swett Mar 18 '12 at 16:39
    
yes that's correct. –  j coe Mar 18 '12 at 17:01

3 Answers 3

up vote 10 down vote accepted

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.

share|improve this answer
    
Any help getting mathjax to understand tabulars? It's not math mode, so how do I signal non-math TeX? –  EMS Mar 18 '12 at 16:27
    
Try "array" instead. (Why should it not be math?) –  TMM Mar 18 '12 at 16:29
    
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. –  EMS Mar 18 '12 at 16:30
1  
Just replacing the "tabular" by "array" in your original version gives you a table with lines between rows/columns like you wanted. –  TMM Mar 18 '12 at 16:32
    
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. –  j coe Mar 18 '12 at 16:36

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.

share|improve this answer
6  
Would the downvoter care to comment on why? –  robjohn Mar 19 '12 at 19:22
    
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$ –  j coe Mar 20 '12 at 12:02
    
I noticed that EMS's intermediates seem to be the differences of mine. I have not tried to analyze EMS's answer yet. –  robjohn Mar 20 '12 at 15:16
    
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? –  j coe Mar 21 '12 at 10:03
    
@jcoe: I use Intaglio for the Mac. There is no Windows port so far. –  robjohn Mar 21 '12 at 12:16

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$

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.