Some time ago I asked about the average waiting time to get two successive 6 when throwing a dice.  To find an answer to this, we could start by simulating the problem in a programming language or in EMT.

>function sim (n=100000) ...
$  sum=0; k=1000;
$  loop 1 to n;
$     v=intrandom(1,k,6);
$     sum=sum+
$         min(nonzeros(v[1:k-1]==6 && v[2:k]==6))+1;
$  end;
$  return sum/n;
$endfunction
>sim
 41.86826

This is not effective since we generate 1000 random numbers for each run to make sure we get two 6. It should be improved by using a C or Python function for one run. Note that loops in EMT or such a language are not efficient.

>function tinyc sim (n) ...
$  ARG_INT(n);
$  double sum=0.0;
$  int i,j;
$  for (i=0; i<n; i++)
$  {
$     j=1;
$     int k1=rand()%6;
$     while (1)
$     {
$         int k2=rand()%6;
$         j++;
$         if (k1==5 && k2==5) { sum+=j; break; }
$         k1=k2;
$     }
$  }
$  new_real(sum/n);
$endfunction
>sim(1000000)
 42.003386

The answer looks like 42, the answer to everything.

How to compute this? The best way is to formulate the problem in terms of random transitions between states. We have three states here.

1. We had two successive 6.
2. Not 1, but we have just thrown a 6.
3. Not 1 or 2.

Then we can compute the probabilities to get from one state to another with a random throw. We put this in a matrix.

\(P = \begin{pmatrix} 1 & 1/6 & 0 \\ 0 & 0 & 1/6 \\ 0 & 5/6 & 5/6 \end{pmatrix}\)

How do we read this matrix? The columns are the states we are in, and the rows the state we get into. E.g., we have 1/6 chance to get into state 1 from state 2, and 5/6 chance to fall into state 3 from state 2. If we are in state 1, we stay there. This is a probability matrix. The sums of the columns are 1.

To simulate the random changes between states, we can multiply the matrix P with a start vector p=(p1,p2,p3) over and over again. The values p1, p2, p3 are the distribution in the states at time 0. For our dice throws the reasonable start value is (0,0,1). Using this, we can determine the expected fraction of experiments which reach state 1 after 42 dice throws.

>P=[1,1/6,0;0,0,1/6;0,5/6,5/6]; fraction P
         1       1/6         0 
         0         0       1/6 
         0       5/6       5/6 
>matrixpower(P,42).[0;0;1]
      0.636651 
     0.0530119 
      0.310337

For a simulation, we have to modify our C function a little bit. We do the throwing 42 times and count how often two successive 6 appears. Such an experiment would reach state 1.

>function tinyc sim (n) ...
$  ARG_INT(n);
$  double sum=0.0;
$  int i,j;
$  for (i=0; i<n; i++)
$  {
$     int k1=rand()%6;
$     for (j=1; j<42; j++)
$     {
$         int k2=rand()%6;
$         if (k1==5 && k2==5) { sum+=1; break; }
$         k1=k2;
$     }
$  }
$  new_real(sum/n);
$endfunction
>sim(1000000)
 0.636486

Now, we can answer the question of the average waiting time for two successive throws of 6. For this, we define a matrix w(i,j), which is the average waiting time to get from state j to state i.  We clearly have w(i,i)=0, and for i different from j

\(w_{i,j} = p_{i,j} + \sum_{k \ne i} p_{k,j} (w_{i,k}+1) = 1 + \sum_k p_{k,j} w_{i,k}\)

For to go from j to i, we can either go directly with 1 step or over k with an average waiting time of w(i,k)+1. We are only interested in the case i=1. Waiting times from i=1 to other states are infinite. So we get the equations

\(w_{1,2} = 1/6 + 5/6 \, (w_{1,3}+1), \quad w_{1,3} = 5/6 \, (w_{1,3}+1) + 1/6 \, (w_{1,2}+1)\)

With the solution

\(w_{1,2}=36, \quad w_{1,3}=42\)

Using matrices, this can be written as follows, if we take the w(i,j) on the left hand side.

>M=id(3)-P'
             0             0             0 
     -0.166667             1     -0.833333 
             0     -0.166667      0.166667 
>k=2:3; M[k,k]\ones(2,1)
            36 
            42

In fact, we do not need the first column and row of the matrix P, if state 1 is our goal. So, for k successive throws of 6, we use the following code.

>function waiting (k) ...
$  P=setdiag(zeros(k-1,k)_5/6,1,1/6);
$  return (id(k)-P')\ones(k,1);
$endfunction
>waiting(2)
            36 
            42 
>waiting(3)
           216 
           252 
           258

A variation of the simulation in C confirms the average waiting time of 258.

>function tinyc sim (n) ...
$ARG_INT(n);
$double sum=0.0;
$int i,j;
$for (i=0; i<n; i++)
${
$   j=2;
$   int k1=rand()%6,k2=rand()%6;
$   while (1)
$   {
$       int k3=rand()%6;
$       j++;
$       if (k1==5 && k2==5 && k3==5)
$       {
$           sum+=j; break;
$       }
$       k1=k2; k2=k3;
$   }
$}
$new_real(sum/n);
$endfunction
>sim(1000000)
 257.734547