Рассчитать вероятность наблюдения последовательности, используя пакет markovchain

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

Выглядит достаточно близко!