2
$\begingroup$

I am trying to find out the number of times a die needs to be rolled before observing a 4 followed by a 6. I would like to model this problem using a discrete time Markov chain with 3 states:

  • State 1: Start
  • State 2: Rolling a 4
  • State 3: Rolling a 6 Immediately after Rolling a 4 (Absorbing State)

I recognize that :

  • State 1 has a 5/6 probability of going to State 1, and State 1 has a 1/6 probability of going to State 2
  • State 2 has a 4/6 probability of going to State 1, State 2 has a 4/6 probability of going to State 2 and State 2 has a 1/6 probability of going to State 3
  • State 3 has a probability of 1 of going to State 3

Using the R programming language. I wrote a small simulation to demonstrate which state the Markov chain is likely to be in after each step:

    library(expm)
library(ggplot2)
library(reshape2)

 mat1.data <- c(5/6, 1/6, 0, 4/6, 1/6, 1/6, 0,0,1)
 mat <- matrix(mat1.data,nrow=3,ncol=3,byrow=TRUE)

V1<- c(1,0,0)

 results <- list()
for (i in 1:100)

{

iteration = i

my_vec_i =  V1%*% (mat %^% i)

vec_1_i = my_vec_i[1]
vec_2_i = my_vec_i[2]
vec_3_i = my_vec_i[3]

results_tmp = data.frame(iteration,vec_1_i, vec_2_i, vec_3_i)

results[[i]] <- results_tmp

}

results_df <- do.call(rbind.data.frame, results)

colnames(results_df )[2] <- 'State_1'
colnames(results_df )[3] <- 'State_2'
colnames(results_df )[4] <- 'State_3'

mdf <- reshape2::melt(results_df, id.var = "iteration")
ggplot(mdf, aes(x = iteration, y = value, colour = variable)) + 
    geom_point() + 
    geom_line()

I plotted the graph below:

enter image description here

My Question: Its clear to see that as the number of iterations ("steps") go on, we are more and more likely to a roll a 4 followed by a 6 – but is there some way we can determine the "average number of times we need to roll a die before observing a 4 followed by a 6"?

I tried to set up a Markov chain to try and find "number of rolls" needed before observing a 4 followed by a 6, but I can't seem to figure this out. Can someone please show me how to mathematically determine the average number of rolls needed before observing a 4 followed by a 6?

$\endgroup$

2 Answers 2

4
$\begingroup$

Actually, this is (perhaps surprisingly) fairly straightforward.

The probability of rolling a 4 is $1/6$, and each roll is independent. The number of rolls is therefore distributed according to the Geometric distribution with probability parameter $1/6$; the mean of the Geometric is 6, so, on average, it will take us 6 rolls to see a 4.

The probability of rolling a 6 on the roll immediately following a 4 is also $1/6$, and each roll is independent. Therefore, on average, it will take us 6 rolls to see a 6 after a 4. Note that these rolls only occur after a 4 has been observed, which, as we have found above, only occurs on $1/6$ of the rolls.

The final number: $6\cdot6 = 36$.

A little R code provides some extra confidence in our result:

n_rolls <- rep(0,100000)
for (i in seq_along(n_rolls)) {
  r1 <- ceiling(6*runif(1))
  n_rolls[i] <- n_rolls[i] + 1
  repeat {
    r2 <- ceiling(6*runif(1)) 
    n_rolls[i] <- n_rolls[i] + 1
    if (r1 == 4 & r2 == 6) break
    r1 <- r2
  }
}
t.test(n_rolls, mu=36)

with result:

One Sample t-test

data:  n_rolls
t = -1.1521, df = 99999, p-value = 0.2493
alternative hypothesis: true mean is not equal to 36
95 percent confidence interval:
 35.66185 36.08779
sample estimates:
mean of x 
 35.87482 

If we want to do this with a Markov chain, we need to construct an appropriate state space, transition matrix, and starting state probability vector. The state space can take on 36 values, e.g., $\{(1,1), (1,2), \dots, (5,6), (6,6)\}$, with the first value of each pair being the "first" of the two dice rolls and the second value being the "second" of the two dice rolls. The transition matrix is pretty straightforward; state $(1,1)$ transitions to state $(1,1)$ with probability $1/6$, to state $(1,2)$ with probability $1/6$, etc.; we can only transition to states whose first element corresponds to the second element of the current state (e.g., state $(1,3)$ can only transition to states $(3,1), (3,2), \dots, (3,6)$.) We assume we start out in any of the 36 states with equal probability, and the starting number of dice rolls equals $2$ to reflect the first two dice rolls that determine our starting state.

# Construct transition matrix
A <- matrix(0,36,36)
for (i in 1:6) {
  for (j in 1:6) {
    k1 <- 6*(i-1) + j
    k2 <- 6*(j-1)
    A[k1, k2+(1:6)] <- 1/6
  }
}

# Run chain until reach state 24 (corresponding to (4,6))
n_rolls <- rep(2,100000)
for (i in seq_along(n_rolls)) {
  s1 <- ceiling(36*runif(1)) 
  while (s1 != 24) {
    s1 <- sample(1:36, 1, prob=A[s1,])
    n_rolls[i] <- n_rolls[i] + 1
  }
}
t.test(n_rolls, mu=36)

with result:

    One Sample t-test

data:  n_rolls
t = 0.35505, df = 99999, p-value = 0.7226
alternative hypothesis: true mean is not equal to 36
95 percent confidence interval:
 35.82407 36.25377
sample estimates:
mean of x 
 36.03892 
$\endgroup$
2
  • $\begingroup$ @ jbowman: Thank you so much! Can this same thing be shown using a Markov Chain? $\endgroup$
    – stats_noob
    Commented Jun 10, 2022 at 20:16
  • $\begingroup$ Do you mean a code example or a mathematical proof? $\endgroup$
    – jbowman
    Commented Jun 10, 2022 at 20:18
3
$\begingroup$

You can also do an exact calculation. Let $K(x,y)$ be the transition probability for going from point $x$ to point $y$ in some Markov chain. For points $a,b$, let $\tau_{a,b}$ be the expected time to go from $a$ to $b$, with $\tau_{a,a}=0$ by convention. We then have the following recurrence relationships for $a \neq b$:

$\tau_{a,b} = 1 + \sum_{c} K(a,c) \tau_{c,b}$ by looking at what happens after one step.

If your Markov chain has $n$ states, this gives you a collection of $\frac{n(n-1)}{2}$ equations in the same number of unknowns. You can solve using the usual linear algebra approach.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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