马尔可夫链模拟，计算极限分布

Question

我有一个状态为 S={1,2,3,4} 和概率矩阵

的马尔可夫链

P=(.180,.274,.426,.120) (.171,.368,.274,.188) (.161,.339,.375,.125) (.079,.355,.384,.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)
}

代码标记没有错误，幸运的是 :)，但它没有 return 预期的是 (.155,.342,.351,.155)。

例如，f(1000) returns [1] 0.006993007 0.006172840 0.008620690 0.006134969

有人可以告诉我我做错了什么吗？

Answer 1

您的函数正确存储了长度为 Nsim 到 x 的单个马尔可夫链实现，但是 prop1、...、prop4 并不是真正的比例个，...，四个；它们似乎与整个链条中的预期价值更相关。您还高估了四的数量，但@StéphaneLaurent 的回答也涉及到这一点。

然后，一旦固定，你的方法非常大 Nsim 工作因为从，比方说，第 30 步我们已经接近平稳分布，而最初的 30 个值是 "noisy", 它们变得可以忽略不计 Nsim.

另一种方法是针对一些较大且固定的 k 关注 P^k，这应该效率较低，但可能更直观。特别是，在那种情况下，我们模拟 many（为了使大数定律起作用）实现 相对较长的（对于接近限制分布起作用）马尔可夫链。此外，模拟可以写得更紧凑。特别是，考虑我的 :

的概括

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))

通过这种方式，我们在 k-th 步骤中使用 n 对状态的观察来近似平稳分布。

Answer 2

您的代码中有两个错误：

  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

马尔可夫链模拟，计算极限分布

Markov Chain simulation, calculating limit distribution

simulation

statistics

r

markov-chains