R中马尔可夫链的手动仿真(二)
Manual simulation of Markov Chain in R (2)
Consider the Markov chain with state space S = {1, 2}, transition matrix
and initial distribution α = (1/2, 1/2).
- Simulate 5 steps of the Markov chain (that is, simulate X0, X1, . . . , X5). Repeat the simulation 100 times.
我的解决方案:
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
如您所见,我在每次迭代中得到相同的输出。
如何修复我的代码?
有两个问题:
转置矩阵
如果你检查输入的矩阵,它是你想要的转置:
> 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
函数中的所有 X0
更改为 X
即可。
Consider the Markov chain with state space S = {1, 2}, transition matrix
and initial distribution α = (1/2, 1/2).
- Simulate 5 steps of the Markov chain (that is, simulate X0, X1, . . . , X5). Repeat the simulation 100 times.
我的解决方案:
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
如您所见,我在每次迭代中得到相同的输出。
如何修复我的代码?
有两个问题:
转置矩阵
如果你检查输入的矩阵,它是你想要的转置:
> 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
函数中的所有 X0
更改为 X
即可。