给定一个转移概率矩阵，如何找到马尔可夫过程的稳态分布

Question

有两个与此问题相关的线程在Stack Overflow上：

• 如何在给定转移概率矩阵的情况下获得马尔可夫链的稳态分布阐述了什么是转移概率矩阵，并演示了通过使用该矩阵的幂来达到稳态分布；
• 如何找到矩阵何时收敛的循环使用R循环来确定矩阵幂何时收敛。

以上方法直接，但非常昂贵。如果我们有一个顺序为n的转移矩阵，则每次迭代都要计算矩阵-矩阵乘法，其成本为O（n ^ 3）。

有更有效的方法吗？我想到的一件事是使用特征值分解。已知马尔可夫矩阵：

• 可以在复数域中对角化：A = E * D * E^{-1}；
• 具有实特征值为1和其他（复数）特征值的长度小于1。

稳态分布是与特征值1相关联的特征向量，即第一个特征向量。

好吧，理论很好，但我无法使其工作。将第一个链接问题中的矩阵P取出：

P <- structure(c(0, 0.1, 0, 0, 0, 0, 0, 0.1, 0.2, 0, 0, 0, 0, 0, 0.2, 
0.3, 0, 0, 0.5, 0.4, 0.3, 0.5, 0.4, 0, 0, 0, 0, 0, 0.6, 0.4, 
0.5, 0.4, 0.3, 0.2, 0, 0.6), .Dim = c(6L, 6L))

如果我执行以下操作：

Re(eigen(P)$vectors[, 1])
# [1] 0.4082483 0.4082483 0.4082483 0.4082483 0.4082483 0.4082483

根据之前的问题，现在发生了什么？静态分布是：

# [1] 0.002590673 0.025906737 0.116580322 0.310880848 0.272020713 0.272020708

- Zheyuan Li

Answer 1

使用特征分解，我们需要处理 t(P)。

转移概率矩阵的定义在概率/统计学和线性代数中有所不同。在统计学中，P 的所有行之和为1，而在线性代数中，P 的所有列之和为1。因此，我们需要使用 eigen(t(P)) 而不是 eigen(P)：

e <- Re(eigen(t(P))$vectors[, 1])
e / sum(e)
# [1] 0.002590673 0.025906737 0.116580322 0.310880848 0.272020713 0.272020708

我们可以看到，我们只使用了第一个特征向量，即最大特征值的特征向量。因此，没有必要使用eigen来计算所有特征值/特征向量。可以使用幂方法来找到最大特征值的特征向量。让我们在R中实现它：

stydis1 <- function (A) {
  n <- dim(A)[1L]
  ## checking
  if (any(.rowSums(A, n, n) != 1)) 
    stop (" 'A' is not a Markov matrix")
  ## implement power method
  e <- runif (n)
  oldnorm <- sqrt(c(crossprod(e)))
  repeat {
    e <- crossprod(A, e)
    newnorm <- sqrt(c(crossprod(e)))
    if (abs(newnorm / oldnorm - 1) < 1e-8) break
    e <- e / newnorm
    oldnorm <- newnorm
    }
  ## rescale `e` so that it sums up to 1
  c(e / sum(e))
  }

stydis1 (P)
# [1] 0.002590673 0.025906737 0.116580322 0.310880848 0.272020713 0.272020708

结果是正确的。

---

实际上，我们不必利用特征分解。我们可以调整您第二个链接问题中使用的方法。在那里，我们采用了矩阵幂，正如您所评论的那样，这是昂贵的; 但为什么不将其重新转换为矩阵-向量乘法呢？

stydis2 <- function (A) {
  n <- dim(A)[1L]
  ## checking
  if (any(.rowSums(A, n, n) != 1)) 
    stop (" 'A' is not a Markov matrix")
  ## direct computation
  b <- A[1, ]
  oldnorm <- sqrt(c(crossprod(b)))
  repeat {
    b <- crossprod(A, b)
    newnorm <- sqrt(c(crossprod(b)))
    if (abs(newnorm / oldnorm - 1) < 1e-8) break
    oldnorm <- newnorm
    }
  ## return stationary distribution
  c(b)
  }

stydis2 (P)
# [1] 0.002590673 0.025906737 0.116580322 0.310880848 0.272020713 0.272020708

我们从任意的初始分布开始，比如说A[1, ]，然后迭代地应用转移矩阵直到分布收敛。同样地，结果是正确的。

- Zheyuan Li

---

如果我没错的话，这就是计算离散时间马尔可夫链稳态分布的方法。对于连续时间马尔可夫链，我们可以使用LU分解等方法。 - 989

Answer 2

您的向量y = Re(eigen(P)$vectors[, 1])不是一个分布（因为它们不相加），解决方案是P'y = y，而不是x'P = x。 您链接的问答中提供的解决方案大致上可以解决后者：

x = c(0.00259067357512953, 0.0259067357512953, 0.116580310880829, 
0.310880829015544, 0.272020725388601, 0.272020725388601)
all(abs(x %*% P - x) < 1e-10) # TRUE

通过转置P，您可以使用特征值方法：

x2 = Re(eigen(t(P))$vectors[, 1])
x2 <- x2/sum(x2) 
(function(x) all(abs(x %*% P - x) < 1e-10))(
  x2
) # TRUE

在这种情况下，它寻找的是一个不同的静止向量。

- Frank

---

感谢您发布答案。几个小时前我被问到这个问题，现在我想将其作为一个问题和答案发布以与大家分享。但是我决定留出30分钟的时间间隔，看看是否有其他人获得了奖金。我认为您忘记了对特征向量进行归一化处理。特征向量的长度/范数为1，但不是总和为1。 - Zheyuan Li

谢谢。我希望我早知道你计划自问自答，但这并没有花费太多我的时间，而且你的回答看起来很有用；我已经将其加星标以便以后查阅。感谢你的问答。 - Frank

谢谢 :) 但我认为你的回答更实用，值得被采纳，当然，这还要看你的决定。我只是讲解了基本的编程和数学定义，而你的回答则更详细地探讨了两个领域。 - Frank

Answer 3

根据静态概率向量的定义，它是转移概率矩阵的一个左特征向量，其特征值为1。我们可以通过计算矩阵的特征分解来找到这种对象，识别单位特征值，然后为每个单位特征值计算静态概率向量。以下是一个在R中执行此操作的函数。

stationary <- function(P) {
  
  #Get matrix information
  K     <- nrow(P)
  NAMES <- rownames(P)
  
  #Compute the eigendecomposition
  EIGEN <- eigen(P)
  VALS  <- EIGEN$values
  RVECS <- EIGEN$vectors
  LVECS <- solve(VECS)
  
  #Find the unit eigenvalue(s)
  RES <- zapsmall(Mod(VALS - as.complex(rep(1, K))))
  IND <- which(RES == 0)
  N   <- length(IND)
  
  #Find the stationary vector(s)
  OUT <- matrix(0, nrow = N, ncol = K)
  rownames(OUT) <- sprintf('Stationary[%s]', 1:N)
  colnames(OUT) <- NAMES
  for (i in 1:length(IND)) { 
    SSS     <- Re(eigen(t(P))$vectors[, IND[i]])
    OUT[i,] <- SSS/sum(SSS) }
  
  #Give the output
  OUT }

(注意: 使用 eigen 计算的特征分解存在一些数值误差，因此不存在一个特征值完全等于一。为此，我们使用 zapsmall 函数将模数偏差从一中减去以识别单位特征向量。只要没有真正的特征值小于一但又非常接近于一被“zapped”到一，这将给我们正确的答案。)

将此函数应用于您的转移概率矩阵可以正确地识别出在这种情况下唯一的稳态概率向量。计算中存在一些数值误差，但在大多数情况下，这应该是可管理的。

#Compute the stationary probability vector
S <- stationary(P)

#Show this vector and confirm stationarity
S
                     [,1]       [,2]      [,3]      [,4]      [,5]      [,6]
Stationary[1] 0.002590674 0.02590674 0.1165803 0.3108808 0.2720207 0.2720207

S %*% P
                     [,1]       [,2]      [,3]      [,4]      [,5]      [,6]
Stationary[1] 0.002590674 0.02590674 0.1165803 0.3108808 0.2720207 0.2720207

#Show error in computation
c(S %*% P - S)
[1]  4.336809e-17  2.775558e-17  1.110223e-16 -2.775558e-16  1.665335e-16 -5.551115e-17

- Ben