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

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

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

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

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

Answer 1

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

然后，一旦固定，您具有非常大Nsim的方法就可以工作，因为从步骤 30 开始，我们已经接近平稳分布，虽然最初的 30 个值是"嘈杂的"，但它们在大Nsim下变得可以忽略不计。

另一种方法是关注一些大而固定的 k 的 P^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))

通过这种方式，我们使用第 k 步中状态的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