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

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

李哲源.. 5

那么,要使用Eigen分解,我们需要使用t(P).

转移概率矩阵的定义在概率/统计和线性代数之间不同.统计中所有行的P总和为1,而在线性代数中,所有列的P总和为1.因此eigen(P),我们需要eigen(t(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, ],迭代地应用转换矩阵,直到分布收敛.再次,结果是正确的.

那么,要使用Eigen分解,我们需要使用t(P).

转移概率矩阵的定义在概率/统计和线性代数之间不同.统计中所有行的P总和为1,而在线性代数中,所有列的P总和为1.因此eigen(P),我们需要eigen(t(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, ],迭代地应用转换矩阵,直到分布收敛.再次,结果是正确的.