R: Построить матрицу переходов второго порядка и оценочные последовательности

Question

Другие вопросы
Существует еще один вопрос , в котором задается вопрос о том, как построить матрицу переходов второго порядка, однако, похоже, ответ не дает матрицы переходов второго порядка.

Матрица переходов второго порядка и оценка последовательности
Давайте использовать этот набор данных:

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

Какой наилучший способ построить матрицу переходов второго порядка так, чтобыЯ легко могу набрать новую последовательность, с которой сталкиваюсь, как обсуждено здесь .Например, такой, что я могу рассчитать вероятность наблюдения 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