绘制概率 P(Xn=s) 如何随时间变化？答案

【问题标题】：Plot how the probability P(Xn=s) changes as a function of time?绘制概率 P(Xn=s) 如何随时间变化？
【发布时间】：2020-02-08 22:38:11
【问题描述】：

整个问题：我们有一个带有 5 个状态的马尔可夫链模型：s, t, m, f, r

TPM 如下：

      P <- matrix(c(.84,.03,.01,.03,.03,
                    .11,.80,.15,.19,.09,
                    .01,.04,.70,.02,.05,
                    .04,.10,.07,.75,.00,
                    .00,.03,.07,.01,.83),
                  nrow=5

                   )

通过矩阵乘法，极限分布变为：

(0.1478365,0.4149259,0.09555939,0.2163813,0.1252968)

我正在尝试绘制 P(Xn = s) 如何随时间变化。

给定初始分布是 P(X0 = i) = 1/5 即：

                  s   t   m   f   r  
           α = ( 1/5 1/5 1/5 1/5 1/5 )

我需要针对 n = 0、1、2、3、4、5（x 轴）绘制 P(Xn = s)（在 y 轴上）。

【问题讨论】：

标签： r markov-chains markov

【解决方案1】：

expm 包具有矩阵幂函数 %^%。我首先检查了第 10 个状态，但它还没有达到稳定，所以我选择了第 20 个状态

每个连续的状态都在进行中a %*% P, a %*% (P%^%2), ...., a %*% (P%^%20)

library(expm)

#Attaching package: ‘expm’

#The following object is masked from ‘package:Matrix’:

 #   expm

所以我想我已经有了它，无论如何。

evolve <- sapply(1:20, function(n){ a%*% (P %^% n)})  # 5 x 20 matrix
png(); matplot( 1:20, t(evolve) );dev.off()  # matplot needs data in rows

【讨论】：

【解决方案2】：

这是一个基础 R 解决方案，它可以进行矩阵幂运算并查看进化进度

n <- 20
res <- do.call(rbind,Reduce(`%*%`,c(list(a),replicate(n,P,simplify = FALSE)),accumulate = T))

然后，类似于@42- 的绘图方法，您可以使用 w.r.t 来查看演变。 1:n

matplot(seq(nrow(r)),r)

【讨论】：