Search code examples
rmarkov-chainsmarkov

Finding conditional and joint probabilities from a simulation

Consider the Markov chain with state space S = {1, 2} and transition matrix

enter image description here

and initial distribution α = (1/2, 1/2).

Suppose, the source code for simulation is the following:

alpha <- c(1, 1) / 2
mat <- matrix(c(1 / 2, 0, 1 / 2, 1), nrow = 2, ncol = 2) 

chainSim <- function(alpha, mat, n) 
{
  out <- numeric(n)
  out[1] <- sample(1:2, 1, prob = alpha)
  for(i in 2:n)
    out[i] <- sample(1:2, 1, prob = mat[out[i - 1], ])
  out
}

Suppose the following is the result of a 5-step Markov Chain simulation repeated 10 times:

> sim
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]    2    1    1    2    2    2    1    1    1     2
[2,]    2    1    2    2    2    2    2    1    1     2
[3,]    2    1    2    2    2    2    2    1    2     2
[4,]    2    2    2    2    2    2    2    1    2     2
[5,]    2    2    2    2    2    2    2    2    2     2
[6,]    2    2    2    2    2    2    2    2    2     2

What would be the values of the following?

  1. P(X1 = 1, X3 = 1)

  2. P(X5 = 2 | X0 = 1, X2 = 1)

  3. E(X2)

I tried them as follows:

  1. mean(sim[4, ] == 1 && sim[2, ]== 1)
  2. ?
  3. c(1,2) * mean(sim[2, ])

What would be (2)? Am I correct with the rest?

Kindly explain your response.

Solution

  • You are almost correct about 1: there is a difference whether you use && or &, see 

    ?`&&`

    It should be

    mean(sim[1 + 1, ] == 1 & sim[1 + 3, ] == 1)

    Then 2 is given by

    mean(sim[1 + 5, sim[1 + 0, ] == 1 & sim[1 + 2, ] == 1] == 2)

    where you may get NaN in the case if the conditional event {X0 = 1, X2 = 1} doesn't appear in your simulation.

    Lastly, point 3 is 

    mean(sim[1 + 2, ])

    since a natural estimator of the expected value is simply the sample average.