Марковское моделирование цепей, расчет предельного распределения

Question

У меня есть цепь Маркова с состояниями S = ​​{1,2,3,4} и матрицей вероятности

Р = (. 180, 0,274, 0,426, 0,120) (0,171, 0,368, 0,274, 0,188) (0,161, 0,339, 0,375, 0,125) (0,079, 0,355, 0,384, 0,182)

Первый, второй, третий, четвертый ряд соответственно.

При оценке различных степеней P предельное распределение составляет (.155, .342, .351, .155)

Вот мой подход для реализации этого в R с использованием симуляции:

f<-function(Nsim)
{

x<-numeric(Nsim)
x[1]=1 #the seed

ones<-numeric(1)
twos<-numeric(1)
thres<-numeric(1)
fours<-numeric(1)

for(i in 2:Nsim)
{
  if(x[i-1]==1)
    x[i]=sample(1:4,1,prob=c(.180,.274,.426,.120))
  if(x[i-1]==2)
    x[i]=sample(1:4,1,prob=c(.171,.368,.274,.188))
  if(x[i-1]==3)
    x[i]=sample(1:4,1,prob=c(.161,.339,.375,.125))
  if(x[i-1]==4)
    x[i]=sample(1:4,1,prob=c(.079,.355,.384,.182))

}
x

for(i in 1:Nsim)
{
  if(x[i]==1)
    ones<-ones+1
  if(x[i]==2)
    twos<-twos+1
  if(x[i]==3)
    thres<-thres+1
  else
    fours<-fours+1
}

prop1<-1/ones
prop2<-2/twos
prop3<-3/thres
prop4<-4/fours

list<-c(prop1,prop2,prop3,prop4)
return(list)
}

Код не помечает ошибок, к счастью :), но он не возвращает то, что ожидается с (.155,.342,.351,.155).

Например, f(1000) возвращает [1] 0.006993007 0.006172840 0.008620690 0.006134969

Может кто-нибудь сказать мне, что я делаю не так?

Answer 1

В вашем коде есть две ошибки:

  for(i in 1:Nsim)
  {
    if(x[i]==1)
      ones<-ones+1
    else if(x[i]==2) # this 'else' was missing
      twos<-twos+1
    else if(x[i]==3) # this 'else' was missing
      thres<-thres+1
    else
      fours<-fours+1
  }

  prop1<- ones/Nsim # not 1/ones
  prop2<- twos/Nsim # not 2/twos
  prop3<- thres/Nsim # not 3/thres
  prop4<- fours/Nsim # not 4/fours

Answer 2

Ваша функция правильно сохраняет единственную реализацию цепочки Маркова длиной от Nsim до x, но тогда prop1, ..., prop4 на самом деле не являются пропорциями единиц, ..., четверок; похоже, они больше связаны с ожидаемой ценностью во всей этой цепочке. Вы также переоцениваете число четверок, но ответ @ StéphaneLaurent также имеет отношение к этому.

Затем, после исправления, ваш подход с очень большим Nsim работает, потому что, начиная, скажем, с шага 30, мы уже близки к стационарному распределению, и хотя начальные 30 значений "шумят", они становятся незначительными с большой Nsim.

Альтернативным подходом было бы сосредоточиться на P ^k для некоторых больших и фиксированных k, что должно быть менее эффективным, но, вероятно, более интуитивным. В частности, в этом случае мы моделируем много (для закона большого числа для работы) реализации относительно длинной (как для чего-то близкого к предельному распределению, которое будет работать) Маркова цепей. Также симуляция может быть написана гораздо более компактно. В частности, рассмотрим обобщение моего другого ответа :

chainSim <- function(alpha, mat, n) {
  out <- numeric(n)
  out[1] <- sample(1:ncol(mat), 1, prob = alpha)
  for(i in 2:n)
    out[i] <- sample(1:ncol(mat), 1, prob = mat[out[i - 1], ])
  out
}

Теперь давайте смоделируем 30000 цепочек длиной 30, снова начиная с состояния 1, как в вашем случае. Это дает (см. Также здесь )

set.seed(1)
k <- 30
n <- 30000
table(replicate(chainSim(c(1, 0, 0, 0), M, k), n = n)[k, ]) / n
#         1         2         3         4 
# 0.1557333 0.3442333 0.3490333 0.1510000 

, где

M
#       [,1]  [,2]  [,3]  [,4]
# [1,] 0.180 0.274 0.426 0.120
# [2,] 0.171 0.368 0.274 0.188
# [3,] 0.161 0.339 0.375 0.125
# [4,] 0.079 0.355 0.384 0.182

с

M <- structure(c(0.18, 0.171, 0.161, 0.079, 0.274, 0.368, 0.339, 0.355, 
                 0.426, 0.274, 0.375, 0.384, 0.12, 0.188, 0.125, 0.182), .Dim = c(4L, 4L))

Таким образом, мы аппроксимируем стационарное распределение, используя n наблюдения состояния на k-м шаге.