在计算P ^ n时,matrixpower()和markov()之间有什么区别? [英] What is the difference between matrixpower() and markov() when it comes to computing P^n?
问题描述
考虑具有状态空间S = {1, 2, 3, 4}
和转移矩阵的马尔可夫链
Consider a Markov chain with state space S = {1, 2, 3, 4}
and transition matrix
P = 0.1 0.2 0.4 0.3
0.4 0.0 0.4 0.2
0.3 0.3 0.0 0.4
0.2 0.1 0.4 0.3
然后,看看下面的源代码:
And, take a look at the following source code:
# markov function
markov <- function(init,mat,n,labels)
{
if (missing(labels))
{
labels <- 1:length(init)
}
simlist <- numeric(n+1)
states <- 1:length(init)
simlist[1] <- sample(states,1,prob=init)
for (i in 2:(n+1))
{
simlist[i] <- sample(states, 1, prob = mat[simlist[i-1],])
}
labels[simlist]
}
# matrixpower function
matrixpower <- function(mat,k)
{
if (k == 0) return (diag(dim(mat)[1]))
if (k == 1) return(mat)
if (k > 1) return( mat %*% matrixpower(mat, k-1))
}
tmat = matrix(c(0.1, 0.2, 0.4, 0.3,
0.4, 0.0, 0.4, 0.2,
0.3, 0.3, 0.0, 0.4,
0.2, 0.1, 0.4, 0.3), nrow=4, ncol=4, byrow=TRUE)
p10 = matrixpower(mat = tmat, k=10)
rowMeans(p10)
nn <- 10
alpha <- c(0.25, 0.25, 0.25, 0.25)
set.seed(1)
steps <- markov(init=alpha, mat=tmat, n=nn)
table(steps)/(nn + 1)
输出
> rowMeans(p10)
[1] 0.25 0.25 0.25 0.25
>
.
.
.
> table(steps)/(nn + 1)
steps
1 2 3 4
0.09090909 0.18181818 0.18181818 0.54545455
> ?rowMeans
为什么结果如此不同?
在计算 P n 时,使用matrixpower()
和markov()
有什么区别?
What is the difference between using matrixpower()
and markov()
when it come to compute Pn?
推荐答案
当前,您正在比较完全不同的事物.首先,我将不专注于计算P n ,而是关注A * P n ,其中A是初始分布.在这种情况下,matrixpower
可以完成工作:
Currently you are comparing completely different things. First, I'll focus not on computing Pn, but rather A*Pn, where A is the initial distribution. In that case matrixpower
does the job:
p10 <- matrixpower(mat = tmat, k = 10)
alpha <- c(0.25, 0.25, 0.25, 0.25)
alpha %*% p10
# [,1] [,2] [,3] [,4]
# [1,] 0.2376945 0.1644685 0.2857105 0.3121265
分别是经过10个步骤(使用A进行的初始绘制之后)分别处于状态1、2、3、4的真实概率.
those are the true probabilities of being in states 1, 2, 3, 4, respectively, after 10 steps (after the initial draw made using A).
同时,markov(init = alpha, mat = tmat, n = nn)
只是长度nn + 1
的单个实现,并且此实现的最后一个数字与A * P n 有关.因此,为了尝试获得与理论值相似的数字,我们需要使用nn <- 10
进行许多实现,如
Meanwhile, markov(init = alpha, mat = tmat, n = nn)
is only a single realization of length nn + 1
and only the last number of this realization is relevant for A*Pn. So, as to try to get similar numbers to the theoretical ones, we need many realizations with nn <- 10
, as in
table(replicate(markov(init = alpha, mat = tmat, n = nn)[nn + 1], n = 10000)) / 10000
#
# 1 2 3 4
# 0.2346 0.1663 0.2814 0.3177
我在其中模拟10000个实现,并仅获取每个实现的最后状态.
where I simulate 10000 realizations and take only the last state of each realization.
这篇关于在计算P ^ n时,matrixpower()和markov()之间有什么区别?的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!