在R中手动模拟马尔可夫链

Question

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

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

Simulate 5 steps of the Markov chain (that is, simulate X₀, X₁, . . . , X₅). Repeat the simulation 100 times. Use the results of your simulations to solve the following problems.

Estimate P(X₁ = 1|X₀ = 1). Compare your result with the exact probability.

我的解决方案：

# returns Xn 
func2 <- function(alpha1, mat1, n1) 
{
  xn <- alpha1 %*% matrixpower(mat1, n1+1)

  return (xn)
}

alpha <- c(0.5, 0.5)
mat <- matrix(c(0.5, 0.5, 0, 1), nrow=2, ncol=2)
n <- 10


for (variable in 1:100) 
{
   print(func2(alpha, mat, n))
}

如果我运行此代码一次或 100 次（如问题陈述中所述）有什么区别？

我怎样才能从这里找到条件概率？

Answer 1

让

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

是初始分布和转移矩阵。您的 func2 只找到不需要的第 n 步分布，并且不模拟任何东西。相反，我们可以使用

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
}

其中 out[1] 仅使用初始分布生成，然后对于后续项，我们使用转移矩阵。

然后我们有

set.seed(1)
# Doing once
chainSim(alpha, mat, 1 + 5)
# [1] 2 2 2 2 2 2

因此链从 2 开始并由于指定的转移概率而卡在那里。

我们做 100 次

# Doing 100 times
sim <- replicate(chainSim(alpha, mat, 1 + 5), n = 100)
rowMeans(sim - 1)
# [1] 0.52 0.78 0.87 0.94 0.99 1.00

最后一行显示了我们以状态 2 而不是状态 1 结束的频率。这给出了一个（出于许多）原因，为什么 100 次重复提供更多信息：我们在状态 2 中卡住了，只做了一次模拟，在重复 100 次的同时，我们探索了更多可能的路径。

那么条件概率可以用

求出

mean(sim[2, sim[1, ] == 1] == 1)
# [1] 0.4583333

而真实概率为 0.5（由转移矩阵的左上角条目给出）。

在R中手动模拟马尔可夫链

Manual simulation of Markov Chain in R

statistics

r

probability

markov-chains