在 R 和序列搜索中模拟马尔可夫链

Question

所以我正在用 R 模拟马尔可夫链，其中状态为晴天 (S)、多云 (C) 和下雨 (R)，我想计算出晴天之后是晴天的概率连续两个阴天。

这是我目前的情况：

    P = matrix(c(0.7, 0.3, 0.2, 0.2, 0.5, 0.6, 0.1, 0.2, 0.2), 3)
    print(P)
    x = c("S", "C", "R") 
    n = 10000
    states = character(n+100)
    states[1] = "C"

    for (i in 2:(n+100)){
    if (states[i-1] == "S") {cond.prob = P[1,]}
    else if (states[i-1] == "C") {cond.prob = P[2,]}
    else {cond.prob = P[3,]}
    states[i]=sample(x, 1, prob = cond.prob )
    }

    print(states[1:100])
    states = states[-(1:100)]  
    head(states)
    tail(states)
    states[1:200]

最后我得到了一系列状态。我希望将此序列分成三个状态组（对于链中的三天），然后计算等于 SCC 的这三个集合状态的数量。

我运行对我将如何做这件事一无所知，任何帮助将不胜感激！！

Answer 1

假设您想要滑动 window（即 SCC 可能出现在位置 1-3 或 2-4 等），将状态折叠为字符串和正则表达式搜索应该执行技巧：

collapsed <- paste(states, collapse="")
length(gregexpr("SCC", collapsed)[[1]])

另一方面，如果您不想滑动 window（即 SCC 必须位于 1-3、4-6 或 7-9 等位置），那么您可以使用 tapply:

分割序列

indexer <- rep(1:(ceiling(length(states)/3)), each=3, length=length(states))
chopped <- tapply(states, indexer, paste0, collapse="")
sum(chopped == "SCC")

Answer 2

Eric 提供了正确答案，但只是为了完整起见：您可以使用链的均衡分布获得您正在寻找的概率：

# the equilibrium distribution
e <- eigen(t(P))$vectors[,1]
e.norm <- e/sum(e)
# probability of SCC
e.norm[1] * P[1,2] * P[2,2]

这在计算上更便宜，并且会给您更准确的概率估计，因为您的模拟将偏向链的初始状态。

在 R 和序列搜索中模拟马尔可夫链

Simulating a Markov Chain in R and Sequence Search

r

probability

markov-chains