“马尔科夫链”的版本间的差异
小 (→时间可逆性与细致平衡条件) |
|||
| 第269行: | 第269行: | ||
称<math>\pi_0</math>为'''平稳分布'''。 | 称<math>\pi_0</math>为'''平稳分布'''。 | ||
| − | ==== | + | ====时间可逆的马尔科夫链==== |
| + | |||
前面的讨论都沿着时间正向在讨论,下面考虑时间反向的过程。 | 前面的讨论都沿着时间正向在讨论,下面考虑时间反向的过程。 | ||
| 第287行: | 第288行: | ||
此时,反向过程是齐次的。 | 此时,反向过程是齐次的。 | ||
| − | + | =====时间可逆性的定义(细致平衡条件)及可逆分布===== | |
进一步,考虑一种特殊情形,若<math>q_{ij}=p_{ij}</math>,即对于任意的状态<math>i,j</math>,若存在<math>S</math>上的概率分布<math>\pi</math>满足 | 进一步,考虑一种特殊情形,若<math>q_{ij}=p_{ij}</math>,即对于任意的状态<math>i,j</math>,若存在<math>S</math>上的概率分布<math>\pi</math>满足 | ||
<math>\pi_i p_{ij}=\pi_jp_{ji}</math>, | <math>\pi_i p_{ij}=\pi_jp_{ji}</math>, | ||
| − | 此时,称马尔科夫链满足细致平衡条件,也称该马尔科夫链是'''时间可逆的''',称概率分布<math>\pi</math>为<math>P</math> | + | 此时,称马尔科夫链满足细致平衡条件,也称该马尔科夫链是'''时间可逆的''',称概率分布<math>\pi</math>为<math>P</math>的一个'''可逆分布'''。 |
注:<math>\pi_i p_{ij}</math>可以理解为状态<math>i</math>向状态<math>j</math>的输送,同理,<math>\pi_jp_{ji}</math>可以理解为状态<math>j</math>向状态<math>i</math>的输送。细致平衡条件可以理解为任意的状态<math>i</math>与状态<math>j</math>之间的相互输送是相同的。 | 注:<math>\pi_i p_{ij}</math>可以理解为状态<math>i</math>向状态<math>j</math>的输送,同理,<math>\pi_jp_{ji}</math>可以理解为状态<math>j</math>向状态<math>i</math>的输送。细致平衡条件可以理解为任意的状态<math>i</math>与状态<math>j</math>之间的相互输送是相同的。 | ||
| − | + | =====可逆分布与平稳分布===== | |
进一步,由细致平衡条件,则可逆分布是平稳分布,即 | 进一步,由细致平衡条件,则可逆分布是平稳分布,即 | ||
2025年5月12日 (一) 09:32的版本
首先,介绍一个经典的骰子游戏——Yahtzee。在游戏中共有5个六面均匀的骰子,游戏参与者至多有三次投掷机会。为了简化,此处假设目标是使得所有骰子的结果相同(称为获得Yahtzee)。具体规则是在第一次和第二次投掷后,可以仅重新投掷部分数量骰子,而保持其他骰子不动。例如,第一次投掷结果为11234,由于有两个1,可以保持对应的两个骰子不动,重新投掷234对应的三个骰子。[1]
实际上,是否获得Yahtzee是看投掷结果中相同的骰子个数是否等于5。每次投掷中结果相同的骰子个数都是一个随机变量,其可能的取值均为0,1,2,3,4,5。只不过,后面一次投掷结果相同的骰子个数小于前面投掷结果相同的骰子个数发生的概率为0。在整个过程中,投掷结果相同的骰子个数在0-5六个数字间进行转移,所以称集合[math]\displaystyle{ \{0,1,2,3,4,5\} }[/math]为状态空间。此外,当投掷两次骰子后,分析第三次投掷结果中相同的骰子个数的发生概率,可以发现这个概率与第一次投掷结果中相同的骰子个数的发生概率无关。即在给定当前状态下,对于未来状态的预测不再依赖于过去状态。这个特性被称为马尔科夫特性,这一类特殊的随机过程被称为马尔科夫过程。
马尔科夫链是一种应用十分广泛的马尔科夫过程,能够用于描述天气变化和基因调控等现象。本词条主要参考北京大学李东风老师的应用随机过程讲义[2]以及清华大学张颢老师的书籍《随机过程及其应用》[3]对相关基础概念进行梳理,为后续相关应用做铺垫。
离散(或连续)时间马尔科夫链的定义
上述Yahtzee例子中,时间和状态空间均是离散的且满足马尔科夫性质,称这类随机过程为离散时间马尔科夫链(简称马尔科夫链,Markov chain),严格的数学定义如下:
称状态空间[math]\displaystyle{ S=\{1,2,\ldots\} }[/math]上的随机过程[math]\displaystyle{ \{X_n,\ n=0,1,2,\ldots\} }[/math]是马尔科夫链[4]或具有马尔科夫性质,若满足
[math]\displaystyle{ \forall n,\ \ P(X_n=j|X_{n-1}=i,X_{n-2}=x_{n-2},\ldots,X_0=x_0)=P(X_n=j|X_{n-1}=i), }[/math]
其中,[math]\displaystyle{ X_n=j\in S }[/math]表示过程在时刻[math]\displaystyle{ n }[/math]处于状态[math]\displaystyle{ j }[/math],称[math]\displaystyle{ P(X_n=j|X_{n-1}=i) }[/math]为马尔科夫链的一步转移概率,并引入记号[math]\displaystyle{ P_{ij}(m,n)=P(X_n=j|X_m=i) }[/math]。
此外,当时间是连续的但状态空间仍是离散的且满足马尔科夫性质,称这类随机过程为连续时间马尔科夫链,具体定义为
称状态空间[math]\displaystyle{ S=\{1,2,\ldots\} }[/math]上的随机过程[math]\displaystyle{ \{X_t,\ t\geq 0\} }[/math]是连续时间马尔科夫链[5] ,若满足 [math]\displaystyle{ \forall n,\ 0=t_0\lt t_1\lt \ldots\lt t_{n-1}\lt t_{n},\ \ P(X_{t_n}=j|X_{t_{n-1}}=i,X_{t_{n-2}}=x_{t_{n-2}},\ldots,X_{t_{0}}=x_{t_{0}})=P(X_{t_n}=j|X_{t_{n-1}}=i) }[/math]。
日常生活中的许多现象都可以利用马尔科夫链来描述,例如天气的变化。简单考虑两种天气状态,下雨和不下雨。那么,对应的马尔科夫假设就是明天是否下雨只取决于今天的天气,而与之前的天气无关。
离散时间马尔科夫链
基本概念
多步转移概率的遗忘性
上述定义中强调了一步转移概率对于历史的遗忘性,那么对于多步转移概率是否也满足?答案是肯定的。
考虑离散时间马尔科夫链,对任意的[math]\displaystyle{ n\geq0 }[/math]和[math]\displaystyle{ m\geq1 }[/math],有
[math]\displaystyle{ \begin{aligned} P(X_{n+m}=j|X_n=i,X_{n-1}=x_{n-1},\ldots,X_0=x_0)=P(X_{n+m}=j|X_n=i), \end{aligned} }[/math]
即当[math]\displaystyle{ X_n }[/math]给定后,[math]\displaystyle{ X_{n+m} }[/math]和[math]\displaystyle{ X_{n-1},\ldots,X_0 }[/math]条件独立。
注:可以对[math]\displaystyle{ m }[/math]进行数学归纳得证上述性质。
多步转移概率的计算—Chapman-Kolmogorov方程
由于上述性质,即给定现在,过去与未来独立,对于多步转移概率也是成立的。那么,任意步转移概率的计算就可以根据这个性质进行分解计算。
具体来说,对于[math]\displaystyle{ n\geq2 }[/math]步转移概率[math]\displaystyle{ P(X_n=j|X_0=i) }[/math],可以先给定某一中间时刻[math]\displaystyle{ 0\lt m\lt n }[/math]的状态[math]\displaystyle{ k }[/math],此时由条件独立性,从[math]\displaystyle{ i }[/math]先经[math]\displaystyle{ m }[/math]步至[math]\displaystyle{ k }[/math],再经[math]\displaystyle{ n-m }[/math]步至[math]\displaystyle{ j }[/math]的概率就可以被分解为两个转移概率的乘积。接下来,只需要遍历状态[math]\displaystyle{ k }[/math]的所有可能,即
上述过程可写为如下表达式
[math]\displaystyle{ P(X_n=j|X_0=i)=\sum_{k\in S}P(X_n=j|X_m=k)P(X_m=k|X_0=i), }[/math]
称上式为Chapman-Kolmogorov方程。
时齐性与非时齐性
由上述Chapman-Kolmogorov方程可知,多步转移概率其实是一步转移概率的乘积(和)。那么,自然地会想到,在计算乘积的过程中,一步转移概率是否会随时间变化?
在一个经典概率模型——波利亚坛子模型(Polya's urn scheme)中,一步转移概率就与时间有关。模型假设坛子中装有[math]\displaystyle{ b\gt 0 }[/math]个黑球和[math]\displaystyle{ w\gt 0 }[/math]个白球。每次随机地从坛子中有放回地抽出一个球,并且再放入一个同颜色的球。记[math]\displaystyle{ X_n }[/math]为第[math]\displaystyle{ n }[/math]次抽取后坛子中黑球的个数,则一步转移概率
[math]\displaystyle{ P(X_n=m+1|X_{n-1}=m)=\frac{m}{b+w+n-1},\ \ P(X_n=m|X_{n-1}=m)=\frac{b+w+n-1-m}{b+w+n-1}. }[/math]
显然,上述概率是随着[math]\displaystyle{ n }[/math]的变化而变化的。
一般的,当一步转移概率[math]\displaystyle{ P(X_n=j|X_{n-1}=i) }[/math]只与状态[math]\displaystyle{ i,j }[/math]有关,而与时刻[math]\displaystyle{ n }[/math]无关(平稳性假设),称此马尔科夫链为时齐马尔科夫链(homogeneous);若[math]\displaystyle{ P(X_n=j|X_{n-1}=i) }[/math]与状态[math]\displaystyle{ i,j }[/math]和时刻[math]\displaystyle{ n }[/math]均有关,则称此马尔科夫链为非时齐的。
下面主要讨论离散时间时齐马尔科夫链。
时齐马尔科夫链的转移过程——初始状态确定
多步转移概率的计算与一步转移概率矩阵
由时齐性假设,记一步转移概率 [math]\displaystyle{ P(X_n=j|X_{n-1}=i)=p_{ij} }[/math]。
此时,转移过程如下
对应的Chapman-kolmogorov方程的具体表达式如下,
[math]\displaystyle{ \begin{aligned} P_{ij}(0,n) &=\sum_{k_1\in S}\ldots\sum_{k_{n-1}\in S}p_{ik_1}p_{k_1k_2}\ldots p_{k_{n-1}j}. \end{aligned} }[/math]
不难看出,上述方程与起始时刻无关,仅与时间间隔有关,即如果将时间窗任意平移[math]\displaystyle{ l }[/math]个时刻,有
[math]\displaystyle{ P_{ij}(l,n+l)=\sum_{k_1\in S}\ldots\sum_{k_{n-1}\in S}p_{ik_1}p_{k_1k_2}\ldots p_{k_{n-1}j}=P_{ij}(0,n) }[/math]。
所以,一般的,记[math]\displaystyle{ n-m }[/math]步转移概率 [math]\displaystyle{ P_{ij}(m,n)=P_{ij}(n-m) }[/math]。
观察上述方程容易发现,求和项是对处在中间的下角标进行遍历,这恰恰对应着矩阵乘法。 将上述例子中一步转移概率[math]\displaystyle{ p_{ij}(i,j\in S) }[/math]排列为一个矩阵
[math]\displaystyle{ P=(p_{ij})_{i,j\in S}=\begin{bmatrix} p_{11}& p_{12}&\ldots\\ p_{21}& p_{22}&\ldots\\ \vdots&\vdots&\ddots \end{bmatrix}, }[/math]
那么,Chapman-Kolmogorov方程变形为,
[math]\displaystyle{ P(n)=P\times P(n-1)=\ldots=P^n, }[/math]
其中,[math]\displaystyle{ P(n) }[/math]称为[math]\displaystyle{ n }[/math]步转移概率矩阵,[math]\displaystyle{ P }[/math]称为一步转移概率矩阵(transition probability matrix, TPM),下面简称为转移概率矩阵。
通过观察可知,矩阵[math]\displaystyle{ P }[/math]满足如下性质
[math]\displaystyle{ \begin{aligned} p_{ij}&\geq0,\ \forall\ i,j\in S,\\ \sum_{j\in S}p_{ij}&=1,\ \forall\ i\in S。 \end{aligned} }[/math]
一般的,称满足上述性质的矩阵为随机矩阵(stochastic matrix)。
基于转移行为的状态分类
可达、互通与可约性
下面举一个简单的例子具体观察状态间的转移行为,下图中(a)和(b)分别为马尔科夫链的转移概率矩阵和图表示。
由Chapman-Kolmogorov方程,可得100步转移概率矩阵和200步转移概率矩阵为
[math]\displaystyle{ \begin{aligned} P(100)=\left(\begin{matrix}0&0.9&0.1\\0.8&0.2&0\\ 0&0&1 \end{matrix}\right)^{100}=\left(\begin{matrix}0.0041&0.0050&0.9909\\0.0044&0.0052&0.9904\\ 0&0&1 \end{matrix}\right),\ \ \ \ P(200)=\left(\begin{matrix}0&0.9&0.1\\0.8&0.2&0\\ 0&0&1 \end{matrix}\right)^{200}=\left(\begin{matrix}0&0&1\\0&0&1\\ 0&0&1 \end{matrix}\right)。\\ \end{aligned} }[/math]
通过观察马尔科夫链的图表示可以发现,状态1和状态2之间存在互相到达的路径,但状态3与状态1和状态2之间均不能互相到达。此外,容易看出,随着转移步数的增加,对应的转移概率矩阵的第一列和第二列均趋近0,即状态1和状态2的转移行为逐渐一致,但与状态3的转移行为存在明显差异。那么,可以利用这种互相可达的特性来将相同转移行为的状态划分为一类,状态1和状态2是互通的,可被归为一类,而状态3只可达自身,所以状态3为另一个类。下面给出可达和互通的严格定义。
首先,需要明确可达的定义。图中的路径的权重是大于0的概率值,所以定义如下:
[math]\displaystyle{ \exists n\geq 0,\ \ P_{ij}(n)\gt 0 }[/math],则称状态[math]\displaystyle{ i }[/math]可达状态[math]\displaystyle{ j }[/math],记为[math]\displaystyle{ i\to j }[/math]。易知,可达是传递关系,即 [math]\displaystyle{ i\to j,\ j\to k\Longrightarrow i\to k }[/math]。
进一步,若[math]\displaystyle{ i\to j }[/math]且[math]\displaystyle{ j\to i }[/math],则称状态[math]\displaystyle{ i }[/math]和[math]\displaystyle{ j }[/math]互通,记为[math]\displaystyle{ i\leftrightarrow j }[/math]。易知,互通是一种等价关系。那么,可以将互通的状态归为一类,且任一状态仅属于一个类。特别的,若马尔科夫链中只存在一个类,则称它是不可约的(irreducible);否则,称为是可约的(reducible)。
若对例1中的状态3的转移分布稍作修改,
[math]\displaystyle{ \begin{aligned} P'=\left(\begin{matrix}0&0.9&0.1\\0.8&0.2&0\\ 0&0.1&0.9 \end{matrix}\right), \end{aligned} }[/math]
那么,3个状态间彼此互通,即均属于同一个类,则其对应的马尔科夫链是不可约的。
周期性
例2. 下面考察一个简单的两状态马尔科夫链,其转移概率矩阵[math]\displaystyle{ P=\begin{bmatrix}0&1\\1&0\end{bmatrix} }[/math]。
计算[math]\displaystyle{ n }[/math]步转移概率矩阵,得
[math]\displaystyle{ P(n)=\begin{bmatrix}0&1\\1&0\end{bmatrix}^n=\begin{cases} \begin{bmatrix}0&1\\1&0\end{bmatrix},&n=2k+1,\ k\in\mathcal{Z},\\ \begin{bmatrix}1&0\\0&1\end{bmatrix},& n=2k,\ k\in\mathcal{Z}。 \end{cases} }[/math]
不难发现这个马尔科夫链的状态转移存在着周期性——初始状态不论是状态1还是状态2,每转移2步,都能够以正概率返回自身状态,即周期为2。进一步,如果有任一状态能够以一定周期、正概率返回自身,转移矩阵则不存在极限,即状态的周期性成为了一个转移矩阵收敛的判断条件。下面介绍状态的周期性的严格定义:
记使得转移概率[math]\displaystyle{ P_{ii}(n)\gt 0 }[/math]的转移步数[math]\displaystyle{ n }[/math]所组成的集合为[math]\displaystyle{ A_i=\{n:n\geq1,P_{ii}(n)\gt 0\} }[/math]。若[math]\displaystyle{ A_i }[/math]非空,则称其最大公约数为状态[math]\displaystyle{ i }[/math]的周期,记为[math]\displaystyle{ d_i }[/math]。若[math]\displaystyle{ d_i\gt 1 }[/math],称[math]\displaystyle{ i }[/math]是周期的;若[math]\displaystyle{ d_i=1 }[/math],称[math]\displaystyle{ i }[/math]是非周期的。
根据[math]\displaystyle{ n }[/math]步转移概率矩阵可知,[math]\displaystyle{ P_{ii}(n) }[/math]序列为[math]\displaystyle{ 0,1,0,1,\ldots }[/math],[math]\displaystyle{ P_{ij}(n) }[/math]序列为[math]\displaystyle{ 1,0,1,0\ldots }[/math],四个转移概率均不存在极限,即转移矩阵没有极限。由此可见,周期性对于转移矩阵是否具有极限的判断起着重要的作用。
回顾上节例1中的三状态马尔科夫链,由于[math]\displaystyle{ P_{33}=1\gt P_{22}=0.2\gt 0 }[/math],那么[math]\displaystyle{ A_{2}=A_{3}=\{1,2,\ldots\} }[/math],所以它们都是非周期的。对于状态1,其与状态2互通且状态2一步可达自身,那么[math]\displaystyle{ A_{1}=\{2,3,\ldots\} }[/math],所以状态1也是非周期的。可以看出,互通的状态1和状态2的周期性是相同的。实际上,可以通过严格证明得出,互通状态的周期是相同的,即周期性也是一个类性质。那么,互通的状态只需要判断其中任一状态的周期性即可,这大大减少了判断的次数。
上述分类方法是基于局部性质进行两两判断分类。此外,还可以从整体转移机制的角度对状态进行划分,详细介绍参见成块性(Lumpability)。
时齐马尔科夫链的演化——初始状态随机
演化过程
上节中关于转移过程的内容,实际上是固定了初始时刻的状态。一般的,初始时刻的状态也存在着一定的随机性,即可以使用一个分布列来刻画。
记[math]\displaystyle{ X_n(n\geq0) }[/math]的分布列
[math]\displaystyle{ \pi_n=(\pi_n(1),\pi_n(2),\ldots) }[/math],
其中,[math]\displaystyle{ \pi_n(i)=P(X_n=i),\ i\in S }[/math]。
根据全概率公式,[math]\displaystyle{ X_{n+1} }[/math]的分布列
[math]\displaystyle{ \pi_{n+1}=\pi_{n}P }[/math]。
进一步,递推可得,
[math]\displaystyle{ \forall\ n\geq0,\ \ \pi_n=\pi_0P^n }[/math]。
注:[math]\displaystyle{ \pi }[/math]是一个行向量,上述方程为左手方程(left hand equation)。
例1中的马尔科夫链如果在0时刻三个状态等概率地出现,那么在1时刻和2时刻的状态分布分别为
[math]\displaystyle{ \pi_1=(1/3,1/3,1/3)\left(\begin{matrix}0&0.9&0.1\\0.8&0.2&0\\ 0&0&1 \end{matrix}\right)=(0.2667,0.3667,0.3667),\ \ \ \ \ }[/math] [math]\displaystyle{ \pi_2=(1/3,1/3,1/3)\left(\begin{matrix}0&0.9&0.1\\0.8&0.2&0\\ 0&0&1 \end{matrix}\right)^2=(0.2933,0.3133,0.3933) }[/math]。
若将左手方程的左右两边同时进行转置,得
[math]\displaystyle{ \pi_{n}^{\rm T}=(P^{\rm T})^n\pi_{0}^{\rm T} }[/math]。
这是一个一阶线性齐次常微分方程。由递推关系,自然地会联想到,马尔科夫链的长期行为本质上是一个特征值问题。
平稳分布与转移概率矩阵的谱
下面考察[math]\displaystyle{ P }[/math]的特征值,回顾之前例1和例2中的转移矩阵
[math]\displaystyle{ \begin{aligned} \left(\begin{matrix}0&0.9&0.1\\0.8&0.2&0\\ 0&0&1 \end{matrix}\right),\ \ \ \ \left(\begin{matrix}0&1\\1&0\end{matrix}\right),\\ \end{aligned} }[/math]
它们的特征值分别为[math]\displaystyle{ -0.7544,0.9544,1 }[/math]以及[math]\displaystyle{ -1,1 }[/math]。
通过观察不难发现两个共同的性质:1)1都是两个转移矩阵的特征值;2)特征值的模长均不超过1,这两条性质对于一般的随机矩阵仍是成立的。下面,简要分析特征值1。
由于[math]\displaystyle{ P }[/math]的行和为[math]\displaystyle{ 1 }[/math],那么,[math]\displaystyle{ P\mathbf{1}=\mathbf{1} }[/math],即[math]\displaystyle{ 1 }[/math]是[math]\displaystyle{ P }[/math]的特征值,且对应的特征向量为[math]\displaystyle{ \mathbf{1} }[/math]。
另一方面,由于
[math]\displaystyle{ det(P-\lambda I)=det((P-\lambda I)^{\rm T})=det(P^{\rm T}-\lambda I) }[/math]。
所以,[math]\displaystyle{ 1 }[/math]也是[math]\displaystyle{ P^{\rm T} }[/math]的特征值。
注:即方程[math]\displaystyle{ \pi=\pi P }[/math]总有非零解,这与马尔科夫链是否有极限没有关系。
若[math]\displaystyle{ \pi_0 }[/math]是方程[math]\displaystyle{ \pi=\pi P }[/math]的解,那么由递推关系,有
[math]\displaystyle{ \pi_n=\pi_0P^n=\pi_0 }[/math]。
即若[math]\displaystyle{ P(X_0=i)=\pi_0(i) }[/math],有
[math]\displaystyle{ \forall\ n,\ P(X_n=i)=\pi_0(i) }[/math],
称[math]\displaystyle{ \pi_0 }[/math]为平稳分布。
时间可逆的马尔科夫链
前面的讨论都沿着时间正向在讨论,下面考虑时间反向的过程。
[math]\displaystyle{ P(X_n=j|X_{n+1}=i,X_{n+2},\ldots)= \frac{P(X_n=j|X_{n+1}=i)P(X_{n+2},\ldots|X_{n+1}=i,X_n=j)}{P(X_{n+2},\ldots|X_{n+1}=i)} =P(X_n=j|X_{n+1}=i) =\frac{P(X_n=j)P(X_{n+1}=i|X_n=j)}{P(X_{n+1}=i)} }[/math],
即反向过程仍满足马尔科夫性质。这其实是马尔科夫性质的等价定义——给定现在的条件下,过去与未来条件独立。并且,一般来讲,上述一步转移概率是与时间[math]\displaystyle{ n }[/math]有关的,即反向过程是非齐次的。
下面考察有限状态空间上不可约、非周期的马尔科夫链,其唯一的平稳分布为[math]\displaystyle{ \pi }[/math]。若从平稳分布[math]\displaystyle{ \pi }[/math]出发,有
[math]\displaystyle{ P(X_n=j|X_{n+1}=i,X_{n+2},\ldots) =\frac{P(X_n=j)P(X_{n+1}=i|X_n=j)}{P(X_{n+1}=i)}=\frac{\pi_jp_{ji}}{\pi_i}\triangleq q_{ij} }[/math],
此时,反向过程是齐次的。
时间可逆性的定义(细致平衡条件)及可逆分布
进一步,考虑一种特殊情形,若[math]\displaystyle{ q_{ij}=p_{ij} }[/math],即对于任意的状态[math]\displaystyle{ i,j }[/math],若存在[math]\displaystyle{ S }[/math]上的概率分布[math]\displaystyle{ \pi }[/math]满足
[math]\displaystyle{ \pi_i p_{ij}=\pi_jp_{ji} }[/math],
此时,称马尔科夫链满足细致平衡条件,也称该马尔科夫链是时间可逆的,称概率分布[math]\displaystyle{ \pi }[/math]为[math]\displaystyle{ P }[/math]的一个可逆分布。
注:[math]\displaystyle{ \pi_i p_{ij} }[/math]可以理解为状态[math]\displaystyle{ i }[/math]向状态[math]\displaystyle{ j }[/math]的输送,同理,[math]\displaystyle{ \pi_jp_{ji} }[/math]可以理解为状态[math]\displaystyle{ j }[/math]向状态[math]\displaystyle{ i }[/math]的输送。细致平衡条件可以理解为任意的状态[math]\displaystyle{ i }[/math]与状态[math]\displaystyle{ j }[/math]之间的相互输送是相同的。
可逆分布与平稳分布
进一步,由细致平衡条件,则可逆分布是平稳分布,即
[math]\displaystyle{ \pi_i=\sum_k\pi_kP_{ki}\iff \pi=\pi P }[/math]。
证:[math]\displaystyle{ \pi_i=\pi_i\sum_kP_{ik}=\sum_k\pi_iP_{ik}=\sum_k\pi_kP_{ki} }[/math]。
这为求解平稳分布提供了一种新的方法。考虑如下概率转移矩阵
[math]\displaystyle{ P=\left(\begin{matrix} 0&\frac{1}{2}&0&0&0&\frac{1}{2}\\ \frac{1}{3}&0&\frac{1}{3}&0&0&\frac{1}{3}\\ 0&\frac{1}{2}&0&\frac{1}{2}&0&0\\ 0&0&\frac{1}{3}&0&\frac{1}{3}&\frac{1}{3}\\ 0&0&0&\frac{1}{2}&0&\frac{1}{2}\\ \frac{1}{4}&\frac{1}{4}&0&\frac{1}{4}&\frac{1}{4}&0\\ \end{matrix}\right) }[/math]。
利用细致平衡条件,得
[math]\displaystyle{ \begin{aligned} \frac{\pi_1}{2}&=\frac{\pi_2}{3},\ \ \frac{\pi_2}{3}=\frac{\pi_3}{2},\ \ \frac{\pi_3}{2}=\frac{\pi_4}{3},\ \ \frac{\pi_4}{3}=\frac{\pi_5}{2},\\ \frac{\pi_1}{2}&=\frac{\pi_6}{4},\ \ \frac{\pi_2}{3}=\frac{\pi_6}{4},\ \ \frac{\pi_4}{3}=\frac{\pi_6}{4},\ \ \frac{\pi_5}{2}=\frac{\pi_6}{4}。 \end{aligned} }[/math]
进一步,由求和为1,得可逆分布(也是平稳分布)
[math]\displaystyle{ \pi=\left(\frac{1}{8},\frac{3}{16},\frac{1}{8},\frac{3}{16},\frac{1}{8},\frac{1}{4}\right) }[/math]。
但上述结论反之则不然,考虑如下概率转移矩阵
[math]\displaystyle{ P=\left(\begin{matrix}\frac{1}{2}&\frac{1}{2}&0\\0&\frac{1}{2}&\frac{1}{2}\\ \frac{1}{2}&0&\frac{1}{2}\end{matrix}\right) }[/math]。
易知,其存在平稳分布[math]\displaystyle{ \pi=\left(\frac{1}{3},\frac{1}{3},\frac{1}{3}\right) }[/math],但
[math]\displaystyle{ \pi_1p_{12}-\pi_2p_{21}=1/6\neq 0 }[/math],
即该链不是可逆的。
时间可逆性与动力学可逆性
时间可逆性描述的是在达到平稳后,时间反转前后的转移机制的一致性,即当前状态给定后,不论时间方向是朝向未来还是过去,它以同一个转移分布进行转移。与时间可逆性不同,动力学可逆性是指转移概率矩阵可逆且其逆矩阵仍是一个转移概率矩阵,它描述的则是转移机制的一一对应性。进一步,近似动力学可逆性可以作为一种因果度量,详细内容参见词条基于可逆性的因果涌现理论。
离散时间马尔科夫链的粗粒化和度量
离散时间马尔科夫链的粗粒化
一般的,在状态分类的基础上,对分类后状态空间(宏观尺度上)重新定义马尔可夫转移矩阵的过程被称为马尔科夫概率转移矩阵的粗粒化(coarse-graining),详细介绍参见词条马尔科夫链的粗粒化。下面,简要介绍两种粗粒化方法——成块性和基于SVD分解的粗粒化方法。
马尔科夫链的成块性
可约是基于局部连通性进行两两判断从而对状态空间进行分类,但有时这个划分粒度仍很粗糙,即连通的状态间存在其他性质或结构可以继续划分。
例如,考虑如下马尔科夫链,
[math]\displaystyle{ P=\begin{bmatrix} 0.3&0.3&0.3&0.1\\ 0.2&0.4&0.2&0.2\\ 0&0.2&0.3&0.5\\ 0.2&0&0.2&0.6 \end{bmatrix}, }[/math]
易知,这个马尔科夫链是不可约的。此外,无论从状态1还是2出发,下一时刻转移至状态1或状态2的概率均为0.6,转移至其他两个状态3和4的概率均为0.4。另一方面,状态3和状态4也是如此。那么,基于这样的转移机制,可以将状态1和2划分为1组,状态3和状态4划分为一组,形成两个宏观状态,得到一个新的两状态的马尔科夫转移矩阵
[math]\displaystyle{ P=\begin{bmatrix} 0.6&0.4\\ 0.2&0.8\\ \end{bmatrix}, }[/math]
即同属于一个宏观状态的原始微观状态转移至宏观状态的机制是相同的,这种性质被称为成块性(Lumpability)。定义[math]\displaystyle{ A = \{A_i \} }[/math]为宏观态,成块性的数学表达式为:
[math]\displaystyle{ \forall A_i, A_j, \forall m, n \in A_i, \sum_{k \in A_j} p_{mk} = \sum_{k \in A_j} p_{nk} }[/math]
详细介绍参见词条马尔科夫链的成块性。
马尔科夫链的SVD粗粒化
从转移矩阵谱的角度出发,张江团队提出了基于奇异值分解(SVD)的粗粒化方法[6]。
第一步,对状态空间进行划分。首先,考虑[math]\displaystyle{ N }[/math]阶转移矩阵[math]\displaystyle{ P }[/math]的截断奇异值分解
[math]\displaystyle{ U\Sigma_{r\times r} V^{\rm T} }[/math],
其中,[math]\displaystyle{ r }[/math]为粗粒化后的状态个数。之后,利用矩阵[math]\displaystyle{ V }[/math]对转移矩阵[math]\displaystyle{ P }[/math]做投影操作,得
[math]\displaystyle{ \tilde{P}=PV }[/math]。接下来,利用Kmeans等聚类算法将[math]\displaystyle{ \tilde{P} }[/math]的行向量聚为[math]\displaystyle{ r }[/math]类,得到指示矩阵[math]\displaystyle{ \Phi }[/math],其元素表示如下
[math]\displaystyle{ \Phi_{ij}=\begin{cases} 1,&j=a(i),\\ 0,&j\neq a(i), \end{cases} }[/math]
其中,[math]\displaystyle{ a(i)\in\{1,2,\ldots,r\} }[/math]为[math]\displaystyle{ \tilde{P} }[/math]的第[math]\displaystyle{ i }[/math]行行向量所属类别。
第二步,定义新的转移概率矩阵。假设有限维转移矩阵[math]\displaystyle{ P }[/math]对应的马尔科夫链是不可约且非周期的,则存在唯一平稳分布[math]\displaystyle{ \pi }[/math]。考虑达到平稳后的两步联合概率
[math]\displaystyle{ P(X_n=j,X_{n-1}=i)=P(X_{n-1}=i)P(X_n=j|X_{n-1}=i)=\pi_ip_{ij},\ \ i,j=1,2,\ldots,N }[/math],
称其组成的矩阵[math]\displaystyle{ F }[/math]为平稳流矩阵。接下来,基于平稳流矩阵和指示矩阵,定义简化平稳流矩阵
[math]\displaystyle{ F'=\Phi^{\rm T}F\Phi }[/math]。
为了满足行和为1的约束,对矩阵[math]\displaystyle{ F' }[/math]进行行归一化,得到粗粒化后转移矩阵[math]\displaystyle{ P' }[/math],其元素为
[math]\displaystyle{ P'_{ij}=\frac{F'_{ij}}{\sum_kF'_{ik}} }[/math]。
详细介绍参见词条基于可逆性的因果涌现理论。
离散时间马尔科夫链的度量
转移概率矩阵的有效性息
马尔科夫链可以视为是通信领域中一种特殊的信道——离散无记忆信道。信道的输入和输出分别记为[math]\displaystyle{ X }[/math]和[math]\displaystyle{ Y }[/math],当前输入[math]\displaystyle{ X }[/math]给定的条件下,对应的输出[math]\displaystyle{ Y }[/math]与之前的输入和输出独立。
为了描述信道能够传输的最大信息量,信息论领域提出了信道容量的概念,具体定义为
[math]\displaystyle{ C=\max_{\pi_X}I(X;Y)=\max_{\pi_X} H(Y)-H(Y|X) }[/math]。
基于上述信道容量的概念,Erik Hoel[7]提出了离散马尔科夫动力学的有效信息度量
[math]\displaystyle{ C=\max_{U_X}I(X;Y)=\max_{U_X} H(Y)-H(Y|X), }[/math]
其中,[math]\displaystyle{ U_X }[/math]为均匀分布。
通过粗粒化操作,可以得到宏观尺度下的转移矩阵[math]\displaystyle{ P_M }[/math],并计算有效信息[math]\displaystyle{ {\rm EI}(P_M) }[/math]。进一步,可以通过宏观与微观尺度下的动力学有效信息差值来度量马尔科夫动力学的因果效应,并称这个差值
[math]\displaystyle{ {\rm CE}={\rm EI}(P_M)−{\rm EI}(P_m), }[/math]
为有效信息增益。如果存在某个宏观尺度,使得有效信息增益大于0,则认为该系统发生了因果涌现。详细内容参见词条有效信息。
转移概率矩阵的动力学可逆性
易知,当且仅当[math]\displaystyle{ P }[/math]为置换矩阵,有效信息的值最大。并且,置换矩阵的逆矩阵也是一个置换矩阵,同时也是一个转移概率矩阵,即[math]\displaystyle{ P }[/math]是可逆矩阵且[math]\displaystyle{ P^{−1} }[/math]也是一个转移概率矩阵,则称这类转移概率矩阵是严格动力学可逆的。但是,只有置换矩阵满足严格动力学可逆。所以,需要定义一个可以度量接近动力学可逆程度的指标。
通过观察可以发现置换矩阵具有两方面特性:一方面,它是满秩矩阵;另一方面,它每行或列均为独热向量,具有稀疏性。这两方面都可以通过转移概率矩阵[math]\displaystyle{ P }[/math]的奇异值进行度量。考虑一个[math]\displaystyle{ N }[/math]阶转移概率矩阵[math]\displaystyle{ P }[/math],设其奇异值为[math]\displaystyle{ \sigma_1\geq\sigma_2\geq\ldots\ge\sigma_N\geq0 }[/math],则[math]\displaystyle{ P }[/math]的秩有如下表示
[math]\displaystyle{ r=\sum_{i=1}^N\sigma_i^0, }[/math]
并且[math]\displaystyle{ P }[/math]的稀疏性可以通过矩阵的Frobenius范数度量,
[math]\displaystyle{ \|P\|_F^2=\sum_{i=1}^N\sigma_i^2 }[/math]。
基于上述特性,张江团队定义转移概率矩阵[math]\displaystyle{ P }[/math]的[math]\displaystyle{ \alpha\in(0,2) }[/math]阶近似动力学可逆性
[math]\displaystyle{ \Gamma_\alpha=\sum_{i=1}^N\sigma_i^\alpha }[/math]
作为动力学可逆性的近似[6]。在应用中,常采用[math]\displaystyle{ \Gamma_{\alpha=1} }[/math]即转移概率矩阵的核范数来度量接近动力学可逆的程度。
基于[math]\displaystyle{ \Gamma_{\alpha} }[/math],文中还给出了清晰因果涌现和模糊因果涌现程度的计算方法。值得注意的是,它们均与粗粒化方法无关。详细内容参见词条基于可逆性的因果涌现理论。
连续时间马尔科夫链
连续时间马尔科夫链和离散时间马尔科夫链同样具有马尔科夫性质。它们之间的区别在于,离散时间马尔科夫链可以用不变的转移概率矩阵[math]\displaystyle{ P }[/math]描述,但在连续时间马尔科夫链中,转移概率矩阵[math]\displaystyle{ P(t) }[/math]会依赖于状态之间的时间间隔[math]\displaystyle{ t ,t\in \R }[/math]。[8]
处理实际问题时,[math]\displaystyle{ P(t) }[/math]并不容易得到,所以我们常用容易得到的不依赖于[math]\displaystyle{ t }[/math]的转移速率矩阵[math]\displaystyle{ Q }[/math]来刻画连续时间马尔科夫链的性质,对转移速率矩阵[math]\displaystyle{ Q }[/math]更严格的定义会在后面给出。这部分主要目的是将离散时间马尔科夫链中的概念扩展到连续时间马尔科夫链中,更详细的内容见连续时间马尔科夫链。
回顾一下之前的定义,我们设[math]\displaystyle{ S }[/math]是状态空间,[math]\displaystyle{ \{X(t)\}=\{X(t)|t\geq 0\} }[/math]是以[math]\displaystyle{ S }[/math]为状态空间的连续时间随机过程。如果对于任意正整数[math]\displaystyle{ n,t_0\lt t_1\lt \dots\lt t_{n+1} }[/math]和[math]\displaystyle{ i,j,i_0,i_1,\dots,i_{n-1}\in S }[/math],有
[math]\displaystyle{ \mathbb{P}(X(t_{n+1})=j|X(t_n)=i,X(t_{n-1})=i_{n-1},\dots,X(t_0)=i_0)=\mathbb{P}(X(t_{n+1})=j|X(t_n)=i) }[/math]
则称[math]\displaystyle{ \{X(t)\} }[/math]是连续时间离散状态的马尔科夫链,简称为连续时间马尔科夫链,后面用英文缩写CTMC指代。其中的“链”指状态空间[math]\displaystyle{ S }[/math]是离散的。
基本性质
历史遗忘性
与离散情形中一样,CTMC具有历史遗忘性。即对于[math]\displaystyle{ t\gt 0 }[/math],已知[math]\displaystyle{ X(t)=i }[/math]的条件下,将来随机过程[math]\displaystyle{ \{X(u)|u\gt t\} }[/math]与过去随机过程[math]\displaystyle{ \{X(v)|0\leq v\lt t\} }[/math]独立。
时齐性
和离散时间马尔科夫链相同,我们定义:
[math]\displaystyle{ \mathbb{P}(X(t+s)=j|X(s)=i)=\mathbb{P}(X(t)=j|X(0)=i),s,t\geq 0 }[/math]
的CTMC具有时齐性。后面我们讨论的都是具有时齐性的CTMC。
转移概率与转移概率矩阵
由时齐性,我们可知转移概率与起始时间[math]\displaystyle{ s }[/math]无关:
[math]\displaystyle{ p_{ij}(t)=\mathbb{P}(X(t+s)=j|X(s)=i) }[/math]
并且当[math]\displaystyle{ t=0 }[/math]时,转移概率[math]\displaystyle{ p_{ij}(0) }[/math]是[math]\displaystyle{ \delta }[/math]函数:
[math]\displaystyle{ p_{ij}(0)=\delta_{ij}=\begin{cases} 1 &,j=i\in S\\ 0 &,i\neq j \end{cases} }[/math]
我们定义CTMC[math]\displaystyle{ \{X(t)\} }[/math]的转移概率矩阵为:
[math]\displaystyle{ P(t)=\left(p_{ij}(t)\right)_{i,j\in I} }[/math]
Chapman-Kolmogorov方程
与离散情况中一样,CTMC中也成立Chapman-Kolmogorov方程:
[math]\displaystyle{ p_{ij}(t+s)=\sum_{k\in I}p_{ik}(t)p_{kj}(s),\forall t,s\geq 0 }[/math]
或者写成矩阵形式有:
[math]\displaystyle{ P(t+s)=P(t)P(s)=P(s)P(t) }[/math]
决定CTMC的参数
[math]\displaystyle{ X(t) }[/math]的概率分布[math]\displaystyle{ \mathbf{p}(t) }[/math]由转移概率矩阵[math]\displaystyle{ P(t) }[/math]和[math]\displaystyle{ X(0) }[/math]的概率分布[math]\displaystyle{ \mathbf{p}(0)=(p_1(0),p_2(0),\dots) }[/math]唯一决定:
[math]\displaystyle{ \mathbf{p}(t)=\mathbf{p}(0)P(t) }[/math],
其中[math]\displaystyle{ \mathbf{p}(t)=(p_1(t),p_2(t),\dots),t\geq 0,p_i(t)=\mathbb{P}(X(t)=i),i\in S }[/math]。
转移速率矩阵
离散时间马尔科夫链的一步转移概率矩阵[math]\displaystyle{ P }[/math]可以唯一决定[math]\displaystyle{ n }[/math]步转移概率矩阵:[math]\displaystyle{ P(n)=P^n }[/math]。但是对于CTMC,并不存在一个最小的正数[math]\displaystyle{ \delta }[/math],使得对于任意正数[math]\displaystyle{ t }[/math],存在正整数[math]\displaystyle{ n }[/math],有[math]\displaystyle{ P(t)=P(\delta)^n }[/math]成立。但是对于任意[math]\displaystyle{ \epsilon\gt 0 }[/math],可以证明[math]\displaystyle{ \{P(s)|0\lt s\leq\epsilon\} }[/math]可以决定所有的[math]\displaystyle{ P(t) }[/math]。这就提示我们,[math]\displaystyle{ P(t) }[/math]可能被[math]\displaystyle{ P'(0)=\lim_{s\to 0^+}(P(s)-P(0))/s }[/math]唯一决定。这个极限是否存在呢?
在任何有限之间内只有有限次转移的CTMC,称作规则的CTMC。为了简化讨论,以后无特殊声明时,讨论的都是规则的CTMC。CTMC的轨迹是阶梯函数,并且在跳跃点是右连续的,更重要的是,它的概率转移矩阵[math]\displaystyle{ P(t) }[/math]在[math]\displaystyle{ t=0 }[/math]时,存在右导数:
[math]\displaystyle{ \lim_{t\to 0^+}=\frac{p_{ij}(t)-p_{ij}(0)}{t}=q_{ij} }[/math]
由于[math]\displaystyle{ p_{ij}(t) }[/math]是CTMC从[math]\displaystyle{ i }[/math]出发,经过[math]\displaystyle{ t }[/math]时间后处于[math]\displaystyle{ j }[/math]的概率,所以相应的,我们称[math]\displaystyle{ q_{ij} }[/math]是某一时刻从[math]\displaystyle{ i }[/math]出发,下一步向[math]\displaystyle{ j }[/math]转移的速率或强度。称[math]\displaystyle{ Q=(q_{ij})_{i,j\in S} }[/math]为CTMC的转移速率矩阵或者转移强度矩阵.转移速率矩阵具有下面三个重要的性质:
[math]\displaystyle{ (1)q_{ij}\geq 0,i\neq j }[/math]
[math]\displaystyle{ (2)q_{ii}\leq 0 }[/math]
[math]\displaystyle{ (3)q_{ii}=-\sum_{i\neq j}q_{ij} }[/math]
Kolmogorov方程
CTMC的转移速率矩阵[math]\displaystyle{ Q }[/math]由转移概率矩阵[math]\displaystyle{ P(t) }[/math]唯一决定:[math]\displaystyle{ Q=P'(0) }[/math]。自然要问:在什么条件下转移速率矩阵[math]\displaystyle{ Q }[/math]可以决定转移概率矩阵[math]\displaystyle{ P(t) }[/math]?下面就通过Kolmogorov方程回答这个问题。
回忆Chapman-Kolmogorov方程:
[math]\displaystyle{ p_{ij}(t+s)=\sum_{k\in S}p_{ik}(t)p_{kj}(s),i,j\in S }[/math]
我们将右边称为前,把左边称为后,在上式两边对后面的[math]\displaystyle{ t=0 }[/math]求导数,形式上得到
[math]\displaystyle{ p'_{ij}(s)=\sum_{k\in S}q_{ik}p{kj}(s) }[/math]
写成矩阵形式得到[math]\displaystyle{ \pmb{P'(s)=QP(s)} }[/math],我们上面这两个式子称作CTMC的向后方程。同样可以对Chapman-Kolmogorov方程前面的[math]\displaystyle{ s=0 }[/math]求到,形式上得到:
[math]\displaystyle{ p'_{ij}(t)=\sum_{k\in S}p_{ik}(t)q_{kj} }[/math]
写成矩阵形式得到[math]\displaystyle{ \pmb{P'(t)=P(t)Q} }[/math],这两个式子称作CTMC的向前方程。向前方程和向后方程统称为Kolmogorov方程。
如果CTMC[math]\displaystyle{ \{X(t)\} }[/math]的状态空间[math]\displaystyle{ S }[/math]是有限集合,则称[math]\displaystyle{ \{X(t)\} }[/math]是有限状态CTMC。设状态空间[math]\displaystyle{ S }[/math]有[math]\displaystyle{ n }[/math]个状态,取:
[math]\displaystyle{ q=\max_{i,j\in S}q_{ij}=\max_{i\in S}q_{ii} }[/math]
则矩阵[math]\displaystyle{ Qt }[/math]元素的绝对值都小于[math]\displaystyle{ qt }[/math],[math]\displaystyle{ (Qt)^2 }[/math]元素的绝对值都小于[math]\displaystyle{ n(qt)^2 }[/math],…,[math]\displaystyle{ (Qt)^k }[/math]元素的绝对值都小于[math]\displaystyle{ n^{k-1}(qt)^k }[/math],于是下面级数中每个元素都收敛
[math]\displaystyle{ e^{Qt}\equiv\sum^{\infty}_{k=0}\frac{1}{k!}(Qt)^k }[/math]
再关于t逐项求导得到:
[math]\displaystyle{ (e^{Qt})'=\left(\sum^{\infty}_{k=0}\frac{1}{k!}(Qt)^k\right)'=\sum^{\infty}_{k=1}\frac{1}{(k-1)!}(Qt)^{k-1}Q=e^{Qt}Q }[/math]
利用这个等式以及向前方程(也可以用向后方程证明),我们可以证明有限状态CTMC的转移概率矩阵[math]\displaystyle{ P(t) }[/math]由转移速率矩阵[math]\displaystyle{ Q }[/math]唯一决定,并且有:
[math]\displaystyle{ \pmb{P(t)=\sum^{\infty}_{k=0}\frac{1}{k!}(Qt)^k=e^{Qt},t\geq 0} }[/math]
具体例子
我们可以用CTMC来描述超市收银台的顾客数。假设超市中只有一个收银台,顾客随机到达收银台并排队结账。顾客到达人数程服从泊松过程,平均到达速率为[math]\displaystyle{ \lambda }[/math];而服务时间服从指数分布,平均服务速率为[math]\displaystyle{ \mu }[/math],状态空间为收银台前的顾客数目,此时我们有转移速率矩阵[math]\displaystyle{ Q }[/math]:
[math]\displaystyle{ Q = \begin{bmatrix} -\lambda & \lambda & 0 & 0 & \cdots \\ \mu & -(\lambda + \mu) & \lambda & 0 & \cdots \\ 0 & \mu & -(\lambda + \mu) & \lambda & \cdots \\ 0 & 0 & \mu & -(\lambda + \mu) & \cdots \\ \vdots & \vdots & \vdots & \vdots & \ddots \\ \end{bmatrix} }[/math]
很容易它满足转移速率矩阵的三条性质。其中对角线元素表示离开状态n的总速率为[math]\displaystyle{ -(\lambda+\mu) }[/math],超对角线表示顾客到达的速率,次对角线表示服务完成的速率。矩阵中每行仅有3个非零元素,体现“局部状态转移”特性。在下一节中,我们还会介绍如何求出它的平稳分布。
下面我们给一个维度更少的例子来说明转移速率矩阵在CTMC中相对于转移概率矩阵的优势:
[math]\displaystyle{ Q= \begin{pmatrix} -3 & 2 & 1 \\ 3 & -4 & 1 \\ 0 & 0 & 0 \end{pmatrix} }[/math]
固定时间间隔为1,可以通过数值计算的方式计算矩阵指数,从[math]\displaystyle{ Q }[/math]得到[math]\displaystyle{ P(1) }[/math]:
[math]\displaystyle{ P(1) = e^{Q \cdot 1} \approx \begin{pmatrix} 0.213 & 0.354 & 0.433 \\ 0.354 & 0.223 & 0.423 \\ 0 & 0 & 1 \end{pmatrix} }[/math]
可以看到用[math]\displaystyle{ Q }[/math]来描述连续时间马尔科夫链很直接,而试图用[math]\displaystyle{ P(t) }[/math]描述连续时间马尔科夫链的话,对每一个t都要计算一次矩阵指数。[math]\displaystyle{ P(1) }[/math]第3行[math]\displaystyle{ [0, 0, 1] }[/math],表示一旦进入第三个状态,将永远停留,这也可以直接从[math]\displaystyle{ Q }[/math]中直接读出来。[math]\displaystyle{ Q }[/math]还可以很方便地描述[math]\displaystyle{ P(t) }[/math]中无法直接描述的比例成块性,详见马尔科夫链的成块性。
转移过程和演化性质
可达、互通和可约性
和离散时间马尔科夫链的一样,为了研究CTMC的极限分布以及平稳分布存在的条件,需要引入可达概念。设[math]\displaystyle{ S }[/math]是[math]\displaystyle{ \{X(t)\} }[/math]的状态空间,如果存在[math]\displaystyle{ t\gt 0 }[/math]使得[math]\displaystyle{ p_{ij}(t)\gt 0 }[/math],则称[math]\displaystyle{ i }[/math]可达[math]\displaystyle{ j }[/math],记作[math]\displaystyle{ i\to j }[/math]。如果[math]\displaystyle{ i\to j }[/math]且[math]\displaystyle{ j\to i }[/math],则称[math]\displaystyle{ i,j }[/math]互通,记作[math]\displaystyle{ i\leftrightarrow j }[/math]。如果[math]\displaystyle{ S }[/math]的所有状态互通,则称[math]\displaystyle{ \{X(t)\} }[/math]互通。上面这个例子中,所有状态之间均存在转移路径,因此所有状态互通,该CTMC是不可约的。
对于任何[math]\displaystyle{ h\gt 0 }[/math],称[math]\displaystyle{ \{X(nh)\}=\{X(nh)|n=0,1,2,\dots\} }[/math]为[math]\displaystyle{ \{X(t)\} }[/math]的一个离散骨架或h骨架。容易看出,[math]\displaystyle{ \{X(t)\} }[/math]的h骨架是有1步转移概率[math]\displaystyle{ p_{ij}(h)=\mathbb{P}(X(h)=j|X(0)=i) }[/math]的离散时间马尔科夫链,其状态空间和[math]\displaystyle{ \{X(t)\} }[/math]的状态空间相同。因为总有[math]\displaystyle{ p_{ii}(nh)\gt 0 }[/math],所以对于h骨架来讲,从[math]\displaystyle{ i }[/math]出发可以在任意时刻回到[math]\displaystyle{ i }[/math],说明它的所有状态都是非周期的,该性质为以后的讨论带来许多方便。我们可以用[math]\displaystyle{ \{X(t)\} }[/math]的任意一个h骨架来描述可达性:[math]\displaystyle{ i\to j }[/math]的充分必要条件是对存在一个h骨架[math]\displaystyle{ i \to j }[/math]。
平稳过程
类似离散时间马尔科夫链,可以引入CTMC的平稳分布的定义。设[math]\displaystyle{ P(t) }[/math]是[math]\displaystyle{ \{X(t)\} }[/math]的转移概率矩阵,如果概率分布列[math]\displaystyle{ \mathbf{p}=(p_0,p_1,\dots) }[/math]满足[math]\displaystyle{ \mathbf{p}=\mathbf{p}P(t) }[/math],则称[math]\displaystyle{ \mathbf(p) }[/math]是[math]\displaystyle{ \{X(t)\} }[/math]的平稳分布,或平稳不变分布。
更进一步,我们可以定义平稳过程。如果对于任何[math]\displaystyle{ n\geq1,s\gt 0,0\leq t_0\lt t_1\lt \dots\lt t_n }[/math]随机向量[math]\displaystyle{ (X(t_0+s),X(t_1+s),\dots,X(t_n+s)) }[/math]和[math]\displaystyle{ (X(t_0),X(t_1),\dots,X(t_n)) }[/math]同分布,则称[math]\displaystyle{ \{X(t)\} }[/math]是平稳过程。显然,当初始分布为平稳分布[math]\displaystyle{ \mathbf{p} }[/math]时,[math]\displaystyle{ \{X(t)\} }[/math]为平稳过程。
一般来讲,在实际问题中得到转移速率矩阵[math]\displaystyle{ Q }[/math]要容易得多,但是得到转移概率矩阵[math]\displaystyle{ P(t) }[/math]相对困难一些。所以要得到CTMC的平稳分布,最好能从转移速率矩阵[math]\displaystyle{ Q }[/math]入手。平稳分布的定义[math]\displaystyle{ \mathbf{p}P(t)=\mathbf{p} }[/math]提示我们,形式上关于t求导应当有[math]\displaystyle{ \mathbf{p}Q=\mathbf{p}P'(t)|_{t=0}=\mathbf{p}'=0 }[/math]。事实上,如果互通的CTMC[math]\displaystyle{ \{X(t)\} }[/math]有转移速率矩阵[math]\displaystyle{ Q }[/math],则[math]\displaystyle{ \{X(t)\} }[/math]有唯一的平稳分布的充分必要条件是下面方程组有唯一非负解,这时[math]\displaystyle{ \mathbf{p} }[/math]是[math]\displaystyle{ \{X(t)\} }[/math]的唯一平稳分布。:
[math]\displaystyle{ \begin{cases} \mathbf{p}Q=0,\\ \sum_{j\in S}p_j=1 \end{cases} }[/math]
平稳可逆过程
在离散时间马尔科夫链[math]\displaystyle{ \{X_n\} }[/math]中,如果概率分布[math]\displaystyle{ \mathbf{\pi}=(\pi_0,\pi_1,\dots) }[/math]满足细致平衡条件[math]\displaystyle{ \pi_ip_{ij}=\pi_jp_{ji} }[/math],则称[math]\displaystyle{ \mathbf{\pi} }[/math]是[math]\displaystyle{ \{X_n\} }[/math]的平稳可逆分布,以平稳可逆分布[math]\displaystyle{ \mathbf{\pi} }[/math]为初始分布的[math]\displaystyle{ \{X_(n)\} }[/math]称为平稳可逆的。CTCM中也有类似结论。
设[math]\displaystyle{ \{X(t)\} }[/math]是平稳过程,如果对于任何[math]\displaystyle{ 0\leq t_1\lt t_2\lt \dots\lt t_n\lt T }[/math],随机向量[math]\displaystyle{ (X(T-t_1),X(T-t_2),\dots,X(T-t_n)) }[/math]和[math]\displaystyle{ (X(t_1),X(t_2),\dots,X(t_n)) }[/math]同分布,则称[math]\displaystyle{ \{X(t)\} }[/math]是平稳可逆过程。如果[math]\displaystyle{ \{X(t)\} }[/math]是平稳可逆过程,则可以把时间方向倒过来观测这个过程,并且倒过来观测到的这个过程也是CTMC,转移概率和原来的过程一样。
设[math]\displaystyle{ \{X(t)\} }[/math]的转移概率矩阵是[math]\displaystyle{ P(t) }[/math],[math]\displaystyle{ \mathbf{p}=(p_0,p_1,\dots) }[/math]是概率分布,称
[math]\displaystyle{ p_ip_{ij}(t)=p_jp_{ji}(t)\ i,j\in S,t\geq 0 }[/math]
为细致平衡条件,若[math]\displaystyle{ \mathbf{p} }[/math]满足细致平衡条件,则称[math]\displaystyle{ \mathbf{p} }[/math]为[math]\displaystyle{ \{X(t)\} }[/math]的平稳可逆分布或可逆分布。当[math]\displaystyle{ \mathbf{p} }[/math]为[math]\displaystyle{ \{X(t)\} }[/math]的平稳可逆分布时,[math]\displaystyle{ \mathbf{p} }[/math]也是[math]\displaystyle{ \{X(t)\} }[/math]的平稳分布,并且[math]\displaystyle{ \{X(t)\} }[/math]以[math]\displaystyle{ \mathbf{p} }[/math]为初始分布时,[math]\displaystyle{ \{X(t)\} }[/math]是平稳可逆过程。
同样的,由于实际问题中转移速率矩阵[math]\displaystyle{ Q }[/math]比转移概率矩阵[math]\displaystyle{ P(t) }[/math]更好得到,所以我们也想从转移速率矩阵[math]\displaystyle{ Q }[/math]入手得到平稳可逆分布。对上面给出的细致平衡条件的[math]\displaystyle{ t }[/math]求导数,并利用[math]\displaystyle{ q_{ij}=p'_{ij}(0) }[/math]得到
[math]\displaystyle{ p_iq_{ij}=p_jq_{ji},\ i,j\in S }[/math]
显然若细致平衡条件有解,这个方程就有解;那这个方程有解时,细致平衡条件是否有解呢?事实上,若[math]\displaystyle{ \{X(t)\} }[/math]互通,有转移速率矩阵[math]\displaystyle{ Q }[/math]和平稳分布,则细致平衡条件和这个方程有相同的的非负解,满足这个方程的概率分布[math]\displaystyle{ \mathbf{p} }[/math]是[math]\displaystyle{ \{X(t)\} }[/math]的平稳可逆分布。所以我们也称[math]\displaystyle{ p_iq_{ij}=p_jq_{ji} }[/math]为细致平衡条件。
具体例子
对于一个转移速率矩阵如下的CTMC
[math]\displaystyle{ Q = \begin{pmatrix} -5 & 2 & 3 \\ 1 & -4 & 3 \\ 4 & 1 & -5 \end{pmatrix} }[/math]
通过解平稳过程一节中给出的线性方程组,可以得到CTMC的唯一平稳分布[math]\displaystyle{ \mathbf{p}=(\frac{5}{18},\frac{4}{18},\frac{1}{2}) }[/math]。平稳分布是不是可逆的呢?以状态0和1为例,有[math]\displaystyle{ \mathbf{p}_0q_{01}=\frac{5}{9}\neq\frac{2}{9}=\mathbf{p}_1q_{10} }[/math],不满足细致平衡条件,所以这个CTMC是不可逆的。
而对于另一个转移速率矩阵:
[math]\displaystyle{ Q = \begin{pmatrix} -2 & 1 & 1 \\ 1 & -2 & 1 \\ 1 & 1 & -2 \end{pmatrix} }[/math]
用同样的方法算出它的平稳分布为均匀分布[math]\displaystyle{ \mathbf{p}=(\frac{1}{3},\frac{1}{3},\frac{1}{3}) }[/math],并且它满足[math]\displaystyle{ \mathbf{p}_iq_{ij}=\pi_jq_{ji} }[/math],所以[math]\displaystyle{ \mathbf{p} }[/math]是该转移速率矩阵的平稳可逆分布。
我们再回到最开始给出的超市排队的例子,它的转移速率矩阵为:
[math]\displaystyle{ Q = \begin{bmatrix} -\lambda & \lambda & 0 & 0 & \cdots \\ \mu & -(\lambda + \mu) & \lambda & 0 & \cdots \\ 0 & \mu & -(\lambda + \mu) & \lambda & \cdots \\ 0 & 0 & \mu & -(\lambda + \mu) & \cdots \\ \vdots & \vdots & \vdots & \vdots & \ddots \\ \end{bmatrix} }[/math]
读者可以自行求解平稳分布满足的方程组,会发现只有当[math]\displaystyle{ \lambda\lt \mu }[/math]时,平稳分布的方程组才有唯一非负解,可以解得[math]\displaystyle{ \mathbf{p}_n=(1-\rho)\rho^n,\rho=\lambda/\mu }[/math],当[math]\displaystyle{ \rho\geq 1 }[/math]时,队伍会无限增长,不存在稳态分布。这个平稳分布是否是可逆的呢?由于这个CTMC的状态是局部转移的,我们可以对细致平衡条件进行局部验证:[math]\displaystyle{ \mathbf{p}_n\lambda=\mathbf{p}_{n+1}\mu }[/math],带入[math]\displaystyle{ \mathbf{p}_n=(1-\rho)\rho^n }[/math]可知等式成立,所以它的平稳分布[math]\displaystyle{ \mathbf{p} }[/math]是可逆的。
连续时间马尔科夫链的粗粒化
在连续时间马尔科夫链中粗粒化的本质是构造映射[math]\displaystyle{ \Phi:S \rightarrow S' }[/math],将原始状态空间[math]\displaystyle{ S }[/math]划分为互不相交的区块[math]\displaystyle{ \{B_1,B_2,\dots,B_m\} }[/math],并定义新的转移速率矩阵[math]\displaystyle{ Q' }[/math]描述区块间的转移动力学。理想的粗粒化需要满足三个要求:1. 维度压缩([math]\displaystyle{ |S'|\lt |S| }[/math]);2. 信息保留(关键统计量与动力学特征不变);3. 结构一致性(粗粒度过程仍然为CTMC)。具体实现中,不同方法侧重不同目标,形成多样化的技术路线。
强成块性
当CTMC满足强成块性时,即存在状态划分是的任意区块[math]\displaystyle{ B_k }[/math]内的状态对外的集体转移速率一致。数学上,划分[math]\displaystyle{ \{B_k\} }[/math]需满足:
[math]\displaystyle{ \sum_{j\in B_l}Q^{ij}=\sum_{j\in B_l}Q^{i'j},\forall i,i'\in B_k,\forall l\neq k }[/math]
此时粗粒化后的转移速率矩阵[math]\displaystyle{ Q' }[/math]的元素为:
[math]\displaystyle{ Q'_{kl}=\sum_{j\in B_l}Q^{ij},i\in B_k }[/math]
强成块性确保了区块内状态的不可区分性,粗粒化过程严格保持瞬态与稳态分布[9]。
时间尺度分解
对于有快慢分离特性的系统,快变量可很快达到准稳态,慢变量主导宏观动力学。设快变量状态集[math]\displaystyle{ S_f }[/math]和慢变量状态集[math]\displaystyle{ S_s }[/math],粗粒化后的转移速率矩阵[math]\displaystyle{ Q' }[/math]的元素计算为:
[math]\displaystyle{ Q'_{kl}=\frac{\sum_{i\in B_k}\sum_{j\in B_l}\pi_i^{B_k}Q_{ij}}{\sum_{I\in B_k}\pi_i^{B_k}} }[/math]
其中[math]\displaystyle{ \pi_i^{B_k} }[/math]为区块[math]\displaystyle{ B_k }[/math]内状态的准稳态分布,该方法在生物化学网络和通信协议中广泛应用[10]。
参考文献
- ↑ http://www.stat.yale.edu/~pollard/Courses/251.spring2013/Handouts/Chang-MarkovChains.pdf
- ↑ 李东风. 应用随机过程, 2024, 163-215. https://www.math.pku.edu.cn/teachers/lidf/course/stochproc/stochprocnotes/html/_book/StocProc.pdf
- ↑ 陆大金,张颢. 随机过程及其应用, 清华出版社,2012, 176-255.
- ↑ Brémaud P.Markov chains: Gibbs fields, Monte Carlo simulation, and queues.Springer Science & Business Media, 1999, 53-156.
- ↑ Kendall D G, Batchelor G K, Bingham N H, Hayman W K, Hyland J M E, Lorentz G G, Moffatt H K, Parry W, Razborov A A, Robinson C A and Whittle P. Andrei Nikolaevich Kolmogorov (1903–1987). Bulletin of the London Mathematical Society, 1990, 22(1): 31-100.
- ↑ 6.0 6.1 Zhang J, Tao R Y, Leong H K, Yang M Z, Yuan B. Dynamical reversibility and a new theory of causal emergence based on SVD. npj Complexity, 2025, 2(3): 1-11.
- ↑ Hoel E P, Albantakis L, Tononi G. Quantifying causal emergence shows that macro can beat micro. Proceedings of the National Academy of Sciences, 2013, 110(49): 19790-19795.
- ↑ 何书元. 随机过程, 北京大学出版社,2008, 204-253.
- ↑ Luca Cardelli and Radu Grosu and Kim G. Larsen and Mirco Tribastone and Max Tschaikowski and Andrea Vandin. Lumpability for Uncertain Continuous-Time Markov Chains. QEST, 2021, 391-409.
- ↑ Daria Stepanova, Helen M. Byrne, Philip K. Maini, and Tomás Alarcón. A Method to Coarse-Grain Multiagent Stochastic Systems with Regions of Multistability. Multiscale Model. Simul. 2022, 20(1):404-432.