R中手动模拟马尔可夫链（二）

Question

考虑具有状态空间的马尔可夫链S = {1, 2}，转移矩阵

和初始分布α = (1/2, 1/2)。

1. 模拟马尔可夫链的5个步骤（即模拟X₀、X₁、. . . , X₅).重复模拟 100 次。

我的解决方案：

states <- c(1, 2)
alpha <- c(1, 1)/2
mat <- matrix(c(1/2, 1/2, 0, 1), nrow = 2, ncol = 2) 

nextX <- function(X, pMat)
{
    probVec <- vector()

    if(X == states[1])
    {
        probVec <- pMat[1,]
    }
    if(X==states[2])
    {
        probVec <- pMat[2,]
    }

    return(sample(states, 1, replace=TRUE, prob=probVec))
}

steps <- function(alpha1, mat1, n1)
{
    X0 <- sample(states, 1, replace=TRUE, prob=alpha1)

    if(n1 <=0)
    {
        return (X0)
    }
    else
    {
        vec <- vector(mode="numeric", length=n1)

        for (i in 1:n1) 
        {
            X <- nextX(X0, mat1)
            vec[i] <- X
        }

        return (vec)
    }
}

# steps(alpha1=alpha, mat1=mat, n1=5)

simulate <- function(alpha1, mat1, n1)
{
    for (i in 1:n1) 
    {
        vec <- steps(alpha1, mat1, 5)
        print(vec)
    }
}

simulate(alpha, mat, 100)

输出

> simulate(alpha, mat, 100)
[1] 1 2 2 2 2
[1] 2 1 2 2 2
[1] 1 1 1 1 1
[1] 1 1 1 1 1
[1] 2 2 2 2 2
[1] 2 1 1 2 2
[1] 1 1 1 1 1
[1] 2 2 1 2 2
[1] 1 1 1 1 1
[1] 2 2 1 1 2
[1] 2 2 1 2 2
[1] 1 1 1 1 1
[1] 1 1 1 1 1
[1] 1 1 2 2 2
[1] 1 1 1 1 1
[1] 2 1 2 2 2
[1] 1 1 1 1 1
[1] 1 1 1 1 1
[1] 2 1 2 1 1
[1] 1 1 1 1 1
[1] 1 1 1 1 1
[1] 1 1 1 1 1
[1] 1 1 1 1 1
[1] 1 2 1 2 2
[1] 2 2 2 1 2
[1] 1 1 1 1 1
[1] 1 1 1 1 1
[1] 1 1 1 1 1
[1] 2 1 2 2 1
[1] 1 2 2 2 2
[1] 1 1 2 2 2
[1] 1 2 2 1 2
[1] 1 1 1 1 1
[1] 2 2 1 2 2
[1] 2 2 2 1 1
[1] 1 1 1 1 1
[1] 1 1 1 1 1
[1] 1 1 1 1 1
[1] 1 1 1 2 2
[1] 1 2 1 1 2
[1] 2 2 1 1 1
[1] 1 1 1 1 1
[1] 2 2 2 2 1
[1] 1 1 1 1 1
[1] 1 1 1 1 1
[1] 1 2 2 2 2
[1] 2 1 1 2 2
[1] 1 1 1 1 1
[1] 1 2 1 2 1
[1] 1 1 1 1 1
[1] 1 1 1 1 1
[1] 2 1 2 1 2
[1] 1 1 1 1 1
[1] 1 1 1 1 1
[1] 1 1 1 1 1
[1] 1 1 1 1 1
[1] 2 1 2 1 2
[1] 1 1 1 1 1
[1] 1 1 1 1 1
[1] 2 2 1 1 2
[1] 1 1 1 1 1
[1] 1 1 1 1 1
[1] 2 2 1 1 2
[1] 1 1 1 1 1
[1] 1 2 1 1 2
[1] 1 1 1 1 1
[1] 2 1 1 2 1
[1] 1 1 1 1 1
[1] 2 1 2 2 2
[1] 1 1 1 1 1
[1] 1 1 1 1 1
[1] 1 2 2 2 2
[1] 1 1 1 1 1
[1] 2 2 2 1 2
[1] 2 2 2 1 1
[1] 1 1 2 2 2
[1] 1 1 1 1 1
[1] 2 2 2 1 2
[1] 1 1 1 1 1
[1] 2 2 1 2 2
[1] 2 2 1 2 1
[1] 1 1 1 1 1
[1] 2 2 2 2 2
[1] 1 1 1 1 1
[1] 1 2 1 2 2
[1] 1 1 1 1 1
[1] 2 1 1 2 1
[1] 2 2 2 2 1
[1] 2 2 2 2 2
[1] 1 1 1 1 1
[1] 2 2 2 1 1
[1] 2 2 2 2 2
[1] 1 1 1 1 1
[1] 1 1 1 1 1
[1] 2 1 2 2 1
[1] 2 2 1 1 1
[1] 1 1 1 1 1
[1] 2 2 1 2 2
[1] 2 1 2 2 2
[1] 1 1 1 1 1

如您所见，我在每次迭代中得到相同的输出。

如何修复我的代码？

Answer 1

有两个问题：

转置矩阵

如果你检查你输入的矩阵，它就是你想要的转置：

> mat
     [,1] [,2]
[1,]  0.5    0
[2,]  0.5    1

所以，改变它。

状态不是链式的

在step函数中，返回的状态不用于启动后续状态。相反，X0 只是不断地被传入：

for (i in 1:n1) 
{
    X <- nextX(X0, mat1)
    vec[i] <- X
}

老实说，你根本不需要X0。只需将 step 函数中的所有 X0s 更改为 X 即可。