使用markovchain包计算观察序列的概率

Question

让我们使用来自this question的数据集：

dat<-data.frame(replicate(20,sample(c("A", "B", "C","D"), size = 100, replace=TRUE)))

那么我们就可以建立转移矩阵和马尔可夫链了：

# Build transition matrix
trans.matrix <- function(X, prob=T)
{
  tt <- table( c(X[,-ncol(X)]), c(X[,-1]) )
  if(prob) tt <- tt / rowSums(tt)
  tt
}
trans.mat <- trans.matrix(as.matrix(dat))
attributes(trans.mat)$class <- 'matrix'

# Build markovchain
library(markovchain)
chain <- new('markovchain', transitionMatrix = trans.mat)

如果我现在遇到一个新序列，比如说AAABCAD，那么我可以在给定这个马尔可夫链的情况下计算观察到这个序列的概率吗？

Answer 1

我在markovchain 中看不到一个功能，但它也可以很容易地手动完成。但是有一个警告：转换矩阵不提供观察第一个A 的概率，这需要您提供。让它为 0.25，如果所有四个状态的可能性相同（在您的示例中是这样）。

然后可以得到观察链中的转换

cbind(head(obs, -1), obs[-1])
#      [,1] [,2]
# [1,] "A"  "A" 
# [2,] "A"  "A" 
# [3,] "A"  "B" 
# [4,] "B"  "C" 
# [5,] "C"  "A" 
# [6,] "A"  "D" 

每个转换的概率是

trans.mat[cbind(head(obs, -1), obs[-1])]
# [1] 0.2268722 0.2268722 0.2268722 0.2926316 0.2791165 0.2665198

最终答案是 0.25 *（上述向量的乘积）：

0.25 * prod(trans.mat[cbind(head(obs, -1), obs[-1])])
# [1] 6.355069e-05

---

为了比较，我们可以通过生成许多长度为 7 的链来估计这个概率：

dat <- replicate(2000000, paste(sample(c("A", "B", "C", "D"), size = 7, replace = TRUE), collapse = ""))
mean(dat == "AAABCAD")
# [1] 6.55e-05

看起来够近了！