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

【问题标题】：Simulating a Markov Chain in R and Sequence Search在 R 和序列搜索中模拟马尔可夫链
【发布时间】：2017-11-23 01:03:18
【问题描述】：

因此，我正在使用 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 的状态的数量。

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

【问题讨论】：

标签： r probability markov-chains

【解决方案1】：

假设您想要一个滑动窗口（即 SCC 可能出现在位置 1-3 或 2-4 等），将状态折叠到字符串和正则表达式搜索应该可以解决问题：

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

另一方面，如果您不想要滑动窗口（即 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")

【讨论】：

非常感谢！非常感谢！

【解决方案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]

这会降低计算成本，并为您提供更准确的概率估计，因为您的模拟将偏向链的初始状态。

【讨论】：

谢谢！！这也有帮助！