R：构建二阶转换矩阵和得分序列

Question

其他问题
有another question询问如何建立二阶转移矩阵，但答案似乎没有产生二阶转移矩阵。

二阶转移矩阵和序列评分
让我们使用这个数据集：

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

构建二阶转换矩阵的最佳方法是什么，这样我就可以轻松地对遇到的新序列进行评分，如 here 所讨论的那样。例如，这样我就可以计算观察到AAABCAD 的概率。

---

对朱利叶斯·维诺拉的反应

set.seed(1)
mat <-data.frame(replicate(100,sample(c("AAA", "BBB", "CCC","DDD", "ABC", 'ABD'), size = 5, replace=TRUE)))

aux <- apply(mat, 2, function(col) rbind(paste0(head(col, -2), head(col[-1], -1)), col[-1:-2]))
aux <- data.frame(t(matrix(aux, nrow = 2)))
names(aux) <- c("From", "To")
head(aux, 3)
TM <- table(aux)
TM <- TM / rowSums(TM)


x <- as.character(unlist(mat[1,]))
transitions <- cbind(paste0(head(x, -2), head(x[-1], -1)), x[-1:-2])

prAA <- 1 / (4 * 4)
prAA * prod(TM[transitions])

当我运行此代码时，它给了我0 的概率，但是我计算概率的序列也用于构建转换矩阵（即 df 的第一行，此处为mat）。我想这不应该发生，因为序列是用来构建转换矩阵的，所以没有一个转换可以是零，对吧？

此外，当我将垫子创建更改为这一行时：

mat <-data.frame(replicate(10,sample(c("AAA", "BBB", "CCC","DDD", "ABC", 'ABD'), size = 5, replace=TRUE)))

会报错Error in [.default (TM, transitions) : subscript out of bounds

Answer 1

让我们从矩阵格式的数据开始：

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

为了估计二阶转移矩阵，我们提取以下观察到的转移：

aux <- apply(dat, 2, function(col) rbind(paste0(head(col, -2), head(col[-1], -1)), col[-1:-2]))
aux <- data.frame(t(matrix(aux, nrow = 2)))
names(aux) <- c("From", "To")
head(aux, 3)
#   From To
# 1   DD  D
# 2   DD  B
# 3   DB  A

然后可以估计转移矩阵

TM <- table(aux)
(TM <- TM / rowSums(TM)) # As expected, everything around 0.25
#     To
# From         A         B         C         D
#   AA 0.2459016 0.2950820 0.2049180 0.2540984
#   AB 0.2222222 0.3037037 0.1925926 0.2814815
#   AC 0.3162393 0.1794872 0.1709402 0.3333333
#   AD 0.3211679 0.2189781 0.1824818 0.2773723
#   BA 0.2066116 0.2066116 0.2727273 0.3140496
#   BB 0.2517483 0.2587413 0.2167832 0.2727273
#   BC 0.2647059 0.2745098 0.2254902 0.2352941
#   BD 0.3007519 0.2180451 0.2105263 0.2706767
#   CA 0.2500000 0.2931034 0.2068966 0.2500000
#   CB 0.2178218 0.3168317 0.2178218 0.2475248
#   CC 0.2584270 0.2247191 0.2359551 0.2808989
#   CD 0.3083333 0.2583333 0.2500000 0.1833333
#   DA 0.2402597 0.2727273 0.2272727 0.2597403
#   DB 0.2689076 0.2605042 0.2016807 0.2689076
#   DC 0.2416667 0.2750000 0.2166667 0.2666667
#   DD 0.2442748 0.2213740 0.2671756 0.2671756

---

在您的示例中，我们有由给出的序列和转换

x <- c("A", "A", "A", "B", "C", "A", "D")
(transitions <- cbind(paste0(head(x, -2), head(x[-1], -1)), x[-1:-2]))
#      [,1] [,2]
# [1,] "AA" "A" 
# [2,] "AA" "B" 
# [3,] "AB" "C" 
# [4,] "BC" "A" 
# [5,] "CA" "D" 

类似于我的其他答案，

prAA <- 1 / (4 * 4)
prAA * prod(TM[transitions])
# [1] 6.223154e-05

是观察到x 的概率，其中prAA 是观察到序列的前两个元素AA 的概率（由用户指定）。

---

泛化：n 阶马尔可夫链。

n <- 3
aux <- apply(dat, 2, function(col) {
  from <- head(apply(embed(col, n)[, n:1], 1, paste, collapse = ""), -1)
  to <- col[-1:-n]
  rbind(from, to)
})
aux <- data.frame(t(matrix(aux, nrow = 2)))
names(aux) <- c("From", "To")
TM <- table(aux)
TM <- TM / rowSums(TM)
head(TM)
#      To
# From          A         B         C         D
#   AAA 0.3541667 0.2083333 0.2083333 0.2291667
#   AAB 0.3103448 0.3103448 0.1724138 0.2068966
#   AAC 0.2142857 0.2857143 0.2857143 0.2142857
#   AAD 0.1463415 0.3902439 0.2439024 0.2195122
#   ABA 0.1200000 0.4800000 0.2000000 0.2000000
#   ABB 0.2424242 0.2727273 0.1515152 0.3333333    

x <- c("A", "A", "A", "B", "C", "A", "D")
(transitions <- cbind(head(apply(embed(x, n)[, n:1], 1, paste, collapse = ""), -1), x[-1:-n]))
#      [,1]  [,2]
# [1,] "AAA" "B" 
# [2,] "AAB" "C" 
# [3,] "ABC" "A" 
# [4,] "BCA" "D" 
prAAA <- 1 / 4^n
prAAA * prod(TM[transitions])
# [1] 3.048129e-05