在给定转移概率矩阵的情况下寻找马尔可夫过程的平稳分布

Question

Stack Overflow 上有两个与此问题相关的话题：

描述什么是转移概率矩阵，并演示如何通过取该矩阵的幂达到平稳分布；
使用 R 循环来确定矩阵幂何时收敛。

以上很简单，但是非常昂贵。如果我们有一个 n 阶的转移矩阵，那么在每次迭代中，我们都会以 O(n ^ 3).

的成本计算矩阵-矩阵乘法

有没有更有效的方法来做到这一点？我想到的一件事是使用 Eigen 分解。已知马尔可夫矩阵：

在复数域中可对角化：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

Answer 1

您的向量 y = Re(eigen(P)$vectors[, 1]) 不是分布（因为它加起来不等于 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

不过，它在这个实例中找到了一个不同的固定向量。

Answer 2

嗯，要使用 Eigen 分解，我们需要使用 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 计算所有特征值/向量。 power method 可以用来求一个最大特征值的特征向量。让我们在 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

结果是正确的。

事实上，我们不必利用特征分解。我们可以调整您的第二个链接问题中使用的方法。在那里，我们采用了矩阵电源，正如您所评论的那样昂贵；但为什么不 re-cast 它变成 matrix-vector 乘法呢？

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, ]，并迭代应用转移矩阵，直到分布收敛。同样，结果是正确的。

Answer 3

根据平稳概率向量的定义，是一个left-eigenvector的转移概率矩阵单位特征值 .我们可以通过计算矩阵的特征分解、识别单位特征值，然后计算每个单位特征值的平稳概率向量来找到此类对象。这是 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 计算的特征分解存在一些数值误差，因此不存在完全等于 1 的特征值。因此我们zapsmall 与 1 的模偏差以识别单位特征向量。只要没有小于 1 的真实特征值，这将给我们正确的答案，但非常接近 1，它也得到“变成了一个。）

在这种情况下，将此函数应用于您的转移概率矩阵可以正确识别唯一的固定概率向量。计算中存在少量数值误差，但在大多数情况下应该是可以控制的。

#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

在给定转移概率矩阵的情况下寻找马尔可夫过程的平稳分布

Finding stationary distribution of a markov process given a transition probability matrix

statistics

r

probability

matrix

markov-chains