R симуляция цепи Маркова - PullRequest
Новая
       60

R симуляция цепи Маркова

2 голосов
/ 10 марта 2020 
data1=data.frame("group"=c(1,2,3,4,5),
"t11"=c(0.01,0.32,0.25,0.37,0.11),
"t12"=c(0.48,0.45,0.61,0.29,0.23),
    "t13"=c(0.51,0.23,0.14,0.3,0.67),
    "t22"=c(0.13,0.91,0.41,0.69,0.42),
    "t23"=c(0.87,0.09,0.59,0.31,0.58))
set.seed(1)
    data2=data.frame("student"=c(1:20),
    "group"=c(sample(1:5,rep=T,20)))

Я стараюсь оценить, как ученики переходят через классы. Это примерные данные, где t11 = остаться в классе, t12 = продвинуться вперед, а t13 соответствует окончанию. И так далее. Это data1.

Я видел несколько очень сложных пакетов для имитации результатов для этого типа матрицы вероятностей - мне интересно, есть ли более простой способ симулировать эти 10-кратные шаги, используя data2 в качестве тела студента и перевести их с помощью data1?

1 Ответ

2 голосов
/ 10 марта 2020

Вот базовое решение R.

  • Сначала вы можете определить свою функцию markov для матрицы перехода, т. Е. 
markov <- function(x, n) {
  m <- matrix(0,nrow = 3,ncol = 3)
  m[lower.tri(m,diag = TRUE)] <- c(unlist(x),1)
  r<-(u<-Reduce(`%*%`,replicate(n,m,simplify = FALSE)))[lower.tri(u,diag = TRUE)][-5]
}
  • Затем вы можете добавить получившиеся вероятности к data1, получив data1_10step, и объединить data1_10step с data2
data1_10step <- data1
data1_10step[-1]<-t(apply(data1[-1], 1, markov,10))
data2out <- merge(data2,data1_10step)

так, чтобы 

> data2out
   group student          t11          t12       t13          t22 t23
1      1       1 1.000000e-20 5.514340e-09 1.0000000 1.378585e-09   1
2      1      10 1.000000e-20 5.514340e-09 1.0000000 1.378585e-09   1
3      1       3 1.000000e-20 5.514340e-09 1.0000000 1.378585e-09   1
4      1      18 1.000000e-20 5.514340e-09 1.0000000 1.378585e-09   1
5      1      15 1.000000e-20 5.514340e-09 1.0000000 1.378585e-09   1
6      1      19 1.000000e-20 5.514340e-09 1.0000000 1.378585e-09   1
7      2       4 1.125900e-05 2.970037e-01 0.7029850 3.894161e-01   1
8      2      13 1.125900e-05 2.970037e-01 0.7029850 3.894161e-01   1
9      2      14 1.125900e-05 2.970037e-01 0.7029850 3.894161e-01   1
10     2       7 1.125900e-05 2.970037e-01 0.7029850 3.894161e-01   1
11     3       9 9.536743e-07 5.081030e-04 0.9994909 1.342266e-04   1
12     3       6 9.536743e-07 5.081030e-04 0.9994909 1.342266e-04   1
13     3       8 9.536743e-07 5.081030e-04 0.9994909 1.342266e-04   1
14     4       2 4.808584e-05 2.212506e-02 0.9778269 2.446194e-02   1
15     5      12 1.000000e-10 1.227639e-04 0.9998772 1.708020e-04   1
16     5       5 1.000000e-10 1.227639e-04 0.9998772 1.708020e-04   1
17     5      16 1.000000e-10 1.227639e-04 0.9998772 1.708020e-04   1
18     5      17 1.000000e-10 1.227639e-04 0.9998772 1.708020e-04   1
19     5      20 1.000000e-10 1.227639e-04 0.9998772 1.708020e-04   1
20     5      11 1.000000e-10 1.227639e-04 0.9998772 1.708020e-04   1

РЕДАКТИРОВАТЬ

Если вы отслеживаете процесс markov, вы можете векторизовать markov, то есть 

markov <- Vectorize(function(x, n) {
  m <- matrix(0,nrow = 3,ncol = 3)
  m[lower.tri(m,diag = TRUE)] <- c(unlist(x),1)
  r<-(u<-Reduce(`%*%`,replicate(n,m,simplify = FALSE)))[lower.tri(u,diag = TRUE)][-5]
})

, а затем вы можете отслеживать n от 1 до 10 с помощью 

markov(x,seq(10))
...