如何建立马尔可夫链转移概率矩阵答案

【问题标题】：How to build a Markov's chain transition probability matrix如何建立马尔可夫链转移概率矩阵
【发布时间】：2023-03-14 19:28:01
【问题描述】：

我正在自学 R，但在尝试使用 markovchain 包在 Rstudio 中构建转移概率矩阵时遇到了一些麻烦。首先，我尝试计算 DNA 序列的转换概率。

ATTCAACACATCCAGCCACATGCTCCGAGAGGAGGCAGAGGGCCCCCGGAATGATGCTTACCGAGATTCTTGTTTTTATCCTCGTGGTTGTTTAAAAACGAGTTGAAACTGACGGCATGTCGGACTATAAGCTACTTACTCACCATAGACGTGACCATAGGCCCTAAAACGTTACCGAGATATTCACTTCTAATAACAGTTGTCGGCAGAGCCAAAAGGCCGGGTGATAATACTTTAAAAAGGGAGTTGATTGTTGTATCTAATCCTAGAATGTCAAGAGCGACCATAACAAGATAATTCGGCAGAGCCAGAAAGCGTTCAAGGACTAGAACCATACCGAGACGCAAACGTTCAGGTCGAACTCTAATACCGATTAGT

但是如何以这样的顺序计算转移概率矩阵，我正在考虑使用 R 索引但我真的不知道如何计算这些转移概率。

有没有办法在 R 中做到这一点？我猜矩阵中这些概率的输出应该是这样的：

     A    T    C    G
  A 0.60 0.10 0.10 0.20
  T 0.10 0.50 0.30 0.10
  C 0.05 0.20 0.70 0.05
  G 0.40 0.05 0.05 0.50

【问题讨论】：

标签： r markov-chains

【解决方案1】：

您可以使用markovchain 包来获得帮助。首先，你的数据

seq <- "ATTCAACACATCCAGCCACATGCTCCGAGAGGAGGCAGAGGGCCCCCGGAATGATGCTTACCGAGATTCTTGTTTTTATCCTCGTGGTTGTTTAAAAACGAGTTGAAACTGACGGCATGTCGGACTATAAGCTACTTACTCACCATAGACGTGACCATAGGCCCTAAAACGTTACCGAGATATTCACTTCTAATAACAGTTGTCGGCAGAGCCAAAAGGCCGGGTGATAATACTTTAAAAAGGGAGTTGATTGTTGTATCTAATCCTAGAATGTCAAGAGCGACCATAACAAGATAATTCGGCAGAGCCAGAAAGCGTTCAAGGACTAGAACCATACCGAGACGCAAACGTTCAGGTCGAACTCTAATACCGATTAGT"

然后使用包

library(markovchain)
base_sequence <- strsplit(seq, "")[[1]]
mcX <- markovchainFit(base_sequence)$estimate
mcX

#           A         C         G         T
# A 0.3000000 0.2250000 0.2583333 0.2166667
# C 0.2857143 0.2619048 0.2380952 0.2142857
# G 0.3764706 0.1882353 0.2117647 0.2235294
# T 0.3068182 0.2159091 0.1818182 0.2954545

【讨论】：

【解决方案2】：

创造DNA

DNA <- "ATTCAACACATCCAGCCACATGCTCCGAGAGGAGGCAGAGGGCCCCCGGAATGATGCTTACCGAGATTCTTGTTTTTATCCTCGTGGTTGTTTAAAAACGAGTTGAAACTGACGGCATGTCGGACTATAAGCTACTTACTCACCATAGACGTGACCATAGGCCCTAAAACGTTACCGAGATATTCACTTCTAATAACAGTTGTCGGCAGAGCCAAAAGGCCGGGTGATAATACTTTAAAAAGGGAGTTGATTGTTGTATCTAATCCTAGAATGTCAAGAGCGACCATAACAAGATAATTCGGCAGAGCCAGAAAGCGTTCAAGGACTAGAACCATACCGAGACGCAAACGTTCAGGTCGAACTCTAATACCGATTAGT"

逐个字符分割

DNA_list <- unlist(strsplit(DNA, split = ""))

检索唯一元素

DNA_unique <- unique(DNA_list)

创建一个空矩阵

matrix <- matrix(0, ncol = length(DNA_unique), nrow=length(DNA_unique))

填充它：到elt i和元素i + 1并在矩阵的相应单元格中添加一个。

for (i in 1:(length(DNA_list) - 1)){
  index_of_i <- DNA_unique == DNA_list[i]
  index_of_i_plus_1 <- DNA_unique == DNA_list[i + 1]
  matrix[index_of_i, index_of_i_plus_1] = matrix[index_of_i, index_of_i_plus_1] + 1
}

标准化

matrix <- matrix / rowSums(matrix)

> matrix
          [,1]      [,2]      [,3]      [,4]
[1,] 0.3000000 0.2166667 0.2250000 0.2583333
[2,] 0.3068182 0.2954545 0.2159091 0.1818182
[3,] 0.2857143 0.2142857 0.2619048 0.2380952
[4,] 0.3764706 0.2235294 0.1882353 0.2117647

注意：如果您要计算的 DNA 非常大，可能有一种方法可以更快地执行它。但是这里似乎已经足够快了。

【讨论】：