Likan Zhan

马科夫模型和马科夫链

战立侃 · 2019-08-13

本帖子需要的 R 包

suppressMessages(library(markovchain))
suppressMessages(library(igraph)) # getSlots("markovchain")
suppressMessages(library(expm))   #  %^%
suppressMessages(require(matlab)) # canonicForm
library(ggplot2)

解释生物系统发生的变化并不一定需要知道其内在机制。实际上,变化可被看作是在不同状态(states)间发生的一系列转化(transitions)。为了描述从一种状态转化为另一状态的机率(chance),每种转化都可以被赋予一个概率(probability)。所有可能状态及其转化概率综合到一起就形成了一个具有马科夫特征(Markov property)的随机模型(stochastic model):不同的转化概率仅仅依赖于当前事件状态;即当当前事件已知时,未来事件与过去事件是无关的。

具有马科夫特征的最简单模型是马科夫链(Markov Chain)。例如,在生物学领域,一个细胞会存在三种状态:生长(Growth,G),分裂(Mitosis,M), 和发育停滞(Arrest,A)。在任何一个给定时间点,某个特定细胞的状态可以用一个随机变量 X 来表示。该随机变量只有三个可能值:G、M、或 A,这三个值出现的概率分别是 \(p_G\)\(p_M\)、和 \(p_A\)。这个系统形成的马科夫链是这样一个序列:(\(\textbf{X}_0, \textbf{X}_1, \textbf{X}_2, ...\)),其中 \(X_i\) 指在给定时间点i系统分别处于不同状态时的所有概率组成的一个数组。另外这些概率从 \(\textbf{X}_i\)\(\textbf{X}_{i+1}\) 的变化仅依赖于系统在 \(X_i\) 时不同状态的观测值(G、M、或A)。这个马科夫链的实际实现(Realization)或观测(Observation)是一系列观测值组成的集合,如 G、M、G、G、 …

马科夫链的最简单形式包含两种状态:G 和 M,它是时间上离散 (discrete-time,时间上的步骤是相同的) 和同质 (time-homogeneous, 转化比率是恒定的) 的。在一个给定的时间步骤 (time step) 上,假定位于 G 状态的细胞转化为 M 状态的概率为 p\(_{GM}\) = 0.2,则在该步骤上细胞保持 G 状态的概率为 \(p_{GG} = 1-0.2 = 0.8\)。另外我们假定 \(p_{MG} = 1\), 即细胞只要在M状态,那么它在下一步就会变回 G 状态 (图1)。这些概率定义了一个 2 \(\times\) 2 的转化概率矩阵 (Transition probability matrix)、T。位于该矩阵第 i 行第 j 列的元素 \(T_{ij}\) 描述的就是从状态 i 转化为 状态 j 的机率。因为矩阵中的值是概率,所以每一行的和应该是 1。

如果一个细胞的起始状态是 G,那么其起始状态概率(initial state probability)可以用下面的数组表示:\(\textbf{X}_0 = [p_G = 1, p_m = 0]\)。该细胞将来不同状态的概率可以通过把当前状态概率数组乘以细胞的转化概率矩阵来实现:第二组状态是 \(\textbf{X}_1 = \textbf{X}_0\textbf{T}=[0.8, 0.2]\),第三组状态是 \(\textbf{X}_2 = \textbf{X}_1\textbf{T} = \textbf{X}_0\textbf{T}^2 =[0.84, 0.16]\)。每一次转化都会乘上一个转化矩阵 ,那么 \(n\) 次转化就可以通过乘以 转化矩阵 的 \(n\) 次方来实现,即 \(\textbf{T}^n\)。矩阵中的元素 \((\textbf{T}^n)_{ij}\) 即初始状态 i 经过 n 次转化后变成状态 j 的概率。

X0 <- c(1, 0)
TI2 <- matrix(c(0.8, 0.2, 1, 0), nrow = 2, byrow = TRUE)
(X1 <- X0 %*% TI2)
##      [,1] [,2]
## [1,]  0.8  0.2
(X2 <- X1 %*% TI2)
##      [,1] [,2]
## [1,] 0.84 0.16

为了研究多次转化后的结果,我们可以观察当转化次数趋近于无穷时 \(n\rightarrow \infty\),转化概率矩阵\(\textbf{T}^n\)的变化情况。在我们的双状态系统中, \(\textbf{T}^n\) 的收敛(converge)是非常迅速的(当然其收敛速度取决于概率转化矩阵 \(\textbf{T}\) 本身):在精确到小数点后两位的情况下,\(\textbf{T}^4\) 就是一个非常准确的近似值了。这个转化矩阵中的第 i 行表示的是初始状态 i 经过无限多个步骤后终结于每一种状态的概率。这些描述每一种状态的极限的终结概率合称为极限分布(limiting distribution)。前述例子始终收敛为 \([0.83, 0.17]\)。这说明无论细胞的初始状态是什么,在经过长时间变化后,该细胞会有\(17\%\)的可能性处于\(M\)状态。这个例子中,极限分布也定义了该马科夫链的的稳定状态行为,即稳定分布(stationary distribution)。一个马科夫链一旦进入了稳定分布状态,它就不再受新的转化影响了。

round(TI2 %^% 4, 2)
##      [,1] [,2]
## [1,] 0.83 0.17
## [2,] 0.83 0.17
round(TI2 %^% 100, 2)
##      [,1] [,2]
## [1,] 0.83 0.17
## [2,] 0.83 0.17

一个马科夫链可以根据其可化简性和周期性分为不同的类型。 一个非周期不可化简(aperiodic and irreducible)马科夫链的极限分布和稳定分布是唯一而且相同的。另外,此类马科夫链的极限分布也是独立于其初始状态的,其转化概率矩阵\(\textbf{T}^n\)的每一行都相同。在一个不可化简链(Irreducible Chain)中,所有的状态彼此之间都是可`沟通的’(communicate),即任何一种状态都可以通过零次或更多次的转化。与此相对,一个非周期性链(Aperiodic chain)从一种状态转化为另一种状态的次数则是不固定的:例如,他从 G 到 M 再到 G 的次数可以是任何一个数字。周期性链有一个稳定分布但没有极限分布:如果其起始状态是稳定分布,则链将始终处于稳定分布状态,如果其起始状态不是稳定分布,则链将永远无法到达稳定分布。我们的例子中,如果 \(p_{GM} = p_{MG} = 1\), 则该链就是周期性的:如果链的起始状态是 G,则每一次转化链的状态都会发生反转。此时该链的稳定分布存在,\([0.5, 0.5]\), 但是其转化函数矩阵\(\textbf{T}^n\)不收敛。一个可化简链的极限分布依赖于其起始状态。我们的例子中,如果 \(p_{GG} = p_{MM} = 1\),那么当起始状态是G或M时,该链的稳定分布分别为\([0, 1]\)\([1, 0]\)

# 非周期不可化简的(irreducible and aperiodic)
(XI2 <- matrix(c(1, 0, 0, 1), nrow = 2, byrow = TRUE))
##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1
round(XI2 %*% (TI2 %^% 20), 2)
##      [,1] [,2]
## [1,] 0.83 0.17
## [2,] 0.83 0.17
round(XI2 %*% (TI2 %^% 21), 2)
##      [,1] [,2]
## [1,] 0.83 0.17
## [2,] 0.83 0.17
# 周期性的(periodic)
(TP2 <- matrix(c(0, 1, 1, 0), nrow = 2, byrow = TRUE))
##      [,1] [,2]
## [1,]    0    1
## [2,]    1    0
XI2 %*% (TP2 %^% 20)
##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1
XI2 %*% (TP2 %^% 21)
##      [,1] [,2]
## [1,]    0    1
## [2,]    1    0
# 可化简的(reducible)
(TR2 <- matrix(c(1, 0, 0, 1), nrow = 2, byrow = TRUE))
##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1
XI2 %*% (TR2 %^% 20)
##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1
XI2 %*% (TR2 %^% 21)
##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1

我们在模型中添加状态 A,细胞从 M 状态进入 A 状态的概率是 \(p_{MA} = 0.2\),其他情况下细胞将回到 G 状态。这种发育停滞状态是暂时的,因为在较低概率情况下,细胞能够从 A 状态回到 G 状态,即 \(p_{AG}= 0.1\)。因为这个链也是非周期和不可还价的,\(\textbf{T}^n\) 也会聚合,只是过程会比前面的例子要慢一些。该例极限分布是 \([0.625, 0.125, 0.25]\)。经过很多次转化后,细胞处于 M 状态的概率 \(12.5\%\)

# 三种状态、非周期不可化简的(irreducible and aperiodic)
(XI3 <- matrix(c(1, 0, 0, 0, 1, 0, 0, 0, 1), 
    nrow = 3, byrow = TRUE))
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1
(TI3 <- matrix(c(0.8, 0.2, 0, 0.8, 0, 0.2, 0.1, 0, 0.9), 
    nrow = 3, byrow = TRUE))
##      [,1] [,2] [,3]
## [1,]  0.8  0.2  0.0
## [2,]  0.8  0.0  0.2
## [3,]  0.1  0.0  0.9
round(XI3 %*% (TI3 %^% 50), 3)
##       [,1]  [,2] [,3]
## [1,] 0.625 0.125 0.25
## [2,] 0.625 0.125 0.25
## [3,] 0.625 0.125 0.25

我们可以通过让停滞状态A保持恒定,而让该马科夫链变成一个可化简的链。此时从状态 A 离开的概率均为 0,即 \(p_{AG} = p_{AM} = 0\),链一旦进入状态 A 便无法离开。该链中的 A 状态变成了一个吸入状态(absorbing state),该链也变成了一个吸入性、可化简、和非周期性的链。此类马科夫链具有如下特征:该马科夫链的进化时间越长,其进入吸入状态的概率就越高,而进入 G 或 M 状态的概率就越低。在我们的例子中,无论其起始状态是什么,在经过 63 次转化后 \(p_G\)\(p_M\) 均小于 0.1。 而在 153 次转化后,在两两个小数点的精确水平上,该马科夫链的极限分布就变成了 \([0, 0, 1]\)

# 三种状态、可化简的
(TR3 <- matrix(c(0.8, 0.2, 0, 0.8, 0, 0.2, 0, 0, 1), 
    nrow = 3, byrow = TRUE))
##      [,1] [,2] [,3]
## [1,]  0.8  0.2  0.0
## [2,]  0.8  0.0  0.2
## [3,]  0.0  0.0  1.0
round(XI3 %*% (TR3 %^% 63), 6)
##          [,1]     [,2]     [,3]
## [1,] 0.094597 0.019592 0.885811
## [2,] 0.078367 0.016230 0.905403
## [3,] 0.000000 0.000000 1.000000
round(XI3 %*% (TR3 %^% 153), 2)
##      [,1] [,2] [,3]
## [1,]    0    0    1
## [2,]    0    0    1
## [3,]    0    0    1

细胞进入发育停滞期的比率受其起始状态影响吗?细胞从G或M状态进入发育停滞期的时间有差异吗?所有这些问题都可以通过其转化概率矩阵的基础矩阵(fundamental matrix)来解决:\(\textbf{Fn} = (\textbf{I} - \textbf{Q})^{-1}\),其中 \(\textbf{Q}\) 是概率转化矩阵中不包含吸收状态的子矩阵,\(\textbf{I}\) 是一个单位矩阵(identity matrix)。基础矩阵每一行的和都给出了该状态进入吸收状态的期望次数。在我们的例子中,G 进入吸入状态的次数大约为 30 步。这 30 步中,其中有 25 步处于 G 状态,5 步处于 M 状态。如果起始状态是 M,那么链进入吸收状态的次数则缩减为 25 步,其中 20 步处于 G 状态, 5 步处于 M 状态。

MTR3 <- as(TR3, "markovchain")
States <- c("G", "M", "A")
MTR3@states <- States
attr(MTR3@transitionMatrix, "dimnames") <- list(States, States)
ExtMatr <- extractMatrices(MTR3)
# names(ExtMatr)
FN <- ExtMatr[["N"]]
FN
##    G M
## G 25 5
## M 20 5
FN %*% c(1, 1)
##   [,1]
## G   30
## M   25

从 G 进入吸入状态比从 M 进入吸入状态要多 5 步,原因在于:(a) 系统在开始时处于状态 G 的平均长度是 5 步;(b)一旦链进入 M 状态,其转化概率就与始于 M 状态的情况一样了,因为马科夫链的将来状态仅与现在状态相关。 举个例子,假如一个有偏硬币掷出正面的概率是 \(p_H = 0.2\)。如果把正面 H 作为吸入状态, 那么该链进入吸入状态的时间(absorption time)即为我们第一次掷得正面的平均次数。该平均次数为相应几何分布的平均值,即 \(1/p_H = 5\)。另外,这个硬币掷的反面的概率是 \(Q = p_T = 0.8\),其基础矩阵是 \(Fn = (1 - Q)^{-1} = 5\)

为了进一步考察其变化情况,我们随机产生了2500个起始状态为 G 的马科夫链,其中 \(p_{MA} = 0.2\)。随着时间的延长,A 状态出现的比率稳定增加。模拟数据中,模型进入停滞状态的平均步数是 30.9 步,与用转化概率矩阵\(\textbf{Fn}\)计算的期望值 30 很接近。接下来假定转化概率矩阵相同,但链的起始状态为 M。因为 \(p_{MA} = 0.2\), 所以该马科夫链会有 \(20\%\) 的可能性在第二步就进入停滞状态,另外 \(80\%\) 的可能性回到 \(G\) 状态。高。而同样的转化概率矩阵告诉我们,这 \(80\%\) 的链将与起始点为 \(G\) 的链一样发生进化。所以,如果我们把起始状态为 \(M\) 的马科夫链的步数增加 5, 则这两种不同起始状态马科夫链进入停滞状态的概率几乎完全相同。事实上,如果把链中第一次出现 \(G\) 的位置假定为第一次模拟,那么模拟马科夫链实际上与起始于 \(G\) 的马科夫链具有完全相同的统计特征。

FM <- function(T0){
LL <- sapply(1:2500, 
  function(i) grep("A", rmarkovchain(60, MTR3, t0 = T0, include.t0 = TRUE))) # markovchainSequence
LL <- LL[lengths(LL) > 0L]
#LL[lengths(LL) == 0L] <- 60
LLS <- sapply(LL, min)
mean(LLS)   
}
FM("G")
## [1] 22.57484
FM("M")
## [1] 18.30045
FM("G") - FM("M")
## [1] 3.792273
Simulate <- function(start) {
  LL <- lapply(1:2500, function(i) rmarkovchain(60, MTR3, t0 = start))
  LL <- do.call(rbind, LL)
  LL <- gsub("A", 1, LL)
  LL <- gsub("G|M", 0, LL)
  mode(LL) <- "numeric"
  Probs <- colMeans(LL) # colSums(LL) / 2500
  SS <- data.frame(
      Prob  = Probs, Steps = 1:60, 
      Start = start, stringsAsFactors = FALSE)
  return(SS)
}

DS1 <- Simulate("G")
mean(DS1 $ Prob * DS1 $ Steps)
## [1] 21.73269
DS2 <- Simulate("M")
mean(DS2 $ Prob * DS2 $ Steps)
## [1] 23.40962
DS2 $ Steps <- DS2 $ Steps + 5
mean(DS2 $ Prob * DS2 $ Steps)
## [1] 26.66659
DS  <- rbind(DS1, DS2)

如果我们把吸收概率(absorption probability)提升为 \(p_{MA} = 0.4\),那么基础矩阵中的值将减半,马科夫链进入停滞状态的平均时间也将减半:起始状态为 G 的链需要15步,而起始状态为 M 的链则需要 10 步。与前面的例子一样,起始状态为 M 的链在第二步将有 \(60\%\) 的可能性回到 G 状态,而这些状态的进一步变化将于那些始于 G 状态的变化完全一样。

MTR3@transitionMatrix[2, ] <- c(0.6, 0, 0.4)
ExtMatr <- extractMatrices(MTR3)
# names(ExtMatr)
FN <- ExtMatr[["N"]]
FN
##      G   M
## G 12.5 2.5
## M  7.5 2.5
FN %*% c(1, 1)
##   [,1]
## G   15
## M   10

马科夫链实际上可以进一步拓展到将来情景也受过去短暂几个步骤影响的情景中去。例如,第三条警戒线的位置可以同时受到第一条警戒线和第二条警戒线位置的影响。这种马科夫链需要记住过去两种状态来预测将来第三种状态,被称作二阶马科夫链(second-order markov chain)。

参考文献:

Grewal, J. K., Krzywinski, M., & Altman, N. (2019). Points of significance: Markov models—Markov chains. Nature Methods, 16(8), 661. doi:10.1038/s41592-019-0463-2