14、时间可逆马尔可夫链与马尔可夫链蒙特卡罗方法解析

时间可逆马尔可夫链与马尔可夫链蒙特卡罗方法解析

时间可逆马尔可夫链

考虑一个平稳遍历的马尔可夫链，它已经运行了很长时间，具有转移概率 $P_{ij}$ 和平稳概率 $\pi_i$。若从某个时间开始逆向追踪状态序列，即从时间 $n$ 开始，考虑状态序列 $X_n, X_{n - 1}, X_{n - 2}, \cdots$，可以发现这个逆向序列本身也是一个马尔可夫链，其转移概率 $Q_{ij}$ 定义为：
[Q_{ij} = \frac{\pi_j P_{ji}}{\pi_i}]

要证明逆向过程确实是马尔可夫链，假设当前时间为 $m + 1$。由于 $X_0, X_1, X_2, \cdots$ 是马尔可夫链，给定当前状态 $X_{m + 1}$ 时，未来状态 $X_{m + 2}, X_{m + 3}, \cdots$ 的条件分布与过去状态 $X_m$ 无关。因为独立性是对称关系，所以给定 $X_{m + 1}$ 时，$X_m$ 与 $X_{m + 2}, X_{m + 3}, \cdots$ 也无关，这就验证了逆向过程是马尔可夫链。

若对于所有的 $i$ 和 $j$ 都有 $Q_{ij} = P_{ij}$，则称该马尔可夫链是时间可逆的。时间可逆的条件 $Q_{ij} = P_{ij}$ 也可以表示为：
[\pi_i P_{ij} = \pi_j P_{ji}, \quad \text{对于所有的 } i, j]
这个条件意味着，对于所有状态 $i$ 和 $j$，过程从 $i$ 转移到 $j$ 的速率（即 $\pi_i P_{ij}$）等于从 $j$ 转移到 $i$ 的速率（即 $\pi_j P_{ji}$）。

如果能找到非负且和为 1 的数满足上述方程，那么该马尔可夫链是时间可逆的，且这些数就是极限概率。

下面通过几个例子来进一步说明：
- 随机游走示例 ：考虑一个状态为 $0, 1, \cdots, M$ 的随机游走，该马尔可夫链只能在相邻状态之间转移。由于从 $i$ 到 $i + 1$ 的转移次数与从 $i + 1$ 到 $i$ 的转移次数始终相差不超过 1，所以从 $i$ 到 $i + 1$ 的转移速率等于从 $i + 1$ 到 $i$ 的转移速率，该过程是时间可逆的。通过使每个状态 $i = 0, 1, \cdots, M - 1$ 从 $i$ 到 $i + 1$ 的转移速率与从 $i + 1$ 到 $i$ 的转移速率相等，可以得到极限概率。
- 埃伦费斯特瓮模型 ：假设有 $M$ 个分子分布在两个瓮中，每次随机选择一个分子并将其从一个瓮移到另一个瓮。瓮 I 中的分子数是上述随机游走马尔可夫链的一个特殊情况。利用相关公式可以得到该模型的极限概率，这些概率是二项式概率，这意味着从长远来看，每个球的位置是独立的，并且每个球在两个瓮中的可能性相等。
- 图上的粒子移动示例 ：考虑一个任意的连通图，每条弧 $(i, j)$ 都有一个关联的数 $w_{ij}$。粒子在图的节点之间移动，从节点 $i$ 移动到节点 $j$ 的概率为：
[P_{ij} = \frac{w_{ij}}{\sum_{k} w_{ik}}]
通过时间可逆性方程可以推导出该过程是时间可逆的，并得到其极限概率。

对于一个任意的马尔可夫链，若尝试求解方程 $\sum_{k} x_k P_{kj} = x_j$，通常可能无解。时间可逆的一个必要条件是：
[P_{ik} P_{kj} P_{ji} = P_{ij} P_{jk} P_{ki}, \quad \text{对于所有的 } i, j, k]
这意味着从状态 $i$ 出发，路径 $i \to k \to j \to i$ 的概率与反向路径 $i \to j \to k \to i$ 的概率相同。

下面是一个关于时间可逆马尔可夫链的判断流程图：

graph TD;
    A[给定马尔可夫链] --> B[检查是否满足PikPkjPji = PijPjkPki];
    B -- 是 --> C[可能时间可逆];
    B -- 否 --> D[不是时间可逆];
    C --> E[进一步验证其他条件];
    E -- 通过 --> F[时间可逆];
    E -- 未通过 --> D;

时间可逆性定理

一个遍历的马尔可夫链，当 $P_{ji} = 0$ 时 $P_{ij} = 0$，该马尔可夫链时间可逆的充要条件是：从状态 $i$ 出发，任何回到 $i$ 的路径的概率与反向路径的概率相同。

元素排序示例

假设有 $n$ 个元素，编号为 $1$ 到 $n$，要将它们排列成一个有序列表。每个时间单位会请求检索其中一个元素，元素 $i$ 被请求的概率为 $P_i$。请求后元素会被放回，但不一定在原来的位置。这里考虑两种规则：
- 靠近前端规则 ：被请求的元素会向列表前端移动一位。
- 移到前端规则 ：被请求的元素总是被移到列表前端。

可以将这个过程建模为一个有 $n!$ 个状态的马尔可夫链。通过定理 4.2 可以证明靠近前端规则下的马尔可夫链是时间可逆的，并且可以证明在靠近前端规则下，被请求元素的平均位置小于在移到前端规则下的平均位置。

逆向链的应用

即使过程不是时间可逆的，逆向链的概念也很有用。有如下命题：
考虑一个不可约的马尔可夫链，其转移概率为 $P_{ij}$。如果能找到正数 $\pi_i$（$i \geq 0$），它们的和为 1，以及一个转移概率矩阵 $Q = [Q_{ij}]$，使得：
[\pi_i P_{ij} = \pi_j Q_{ji}]
那么 $Q_{ij}$ 就是逆向链的转移概率，$\pi_i$ 是原始链和逆向链的平稳概率。

通过思考逆向链的性质，有时可以猜测逆向链的形式，然后利用上述方程得到平稳概率和逆向链的转移概率。

例如，一个房间需要一个灯泡照明，当灯泡损坏时，会在第二天开始时更换一个新的。设 $X_n$ 表示第 $n$ 天开始时灯泡的使用天数。假设每个灯泡在第 $i$ 天损坏的概率为 $p_i$，可以发现该马尔可夫链的逆向链有特定的转移概率形式。通过验证相关方程，可以确定逆向链的转移概率和平稳概率。

下面是逆向链应用的步骤列表：
1. 确定原始马尔可夫链的转移概率 $P_{ij}$。
2. 猜测逆向链的转移概率形式 $Q_{ij}$。
3. 寻找满足 $\pi_i P_{ij} = \pi_j Q_{ji}$ 的正数 $\pi_i$。
4. 验证 $\pi_i$ 满足和为 1 的条件。
5. 检查 $\pi_i$ 和 $Q_{ij}$ 是否满足其他相关条件，以确定它们是逆向链的平稳概率和转移概率。

时间可逆马尔可夫链与马尔可夫链蒙特卡罗方法解析

马尔可夫链蒙特卡罗方法

在计算离散随机向量相关的某些值时，常常会遇到困难。设 $X$ 是一个离散随机向量，其可能取值集合为 ${x_j, j \geq 1}$，概率质量函数为 $P(X = x_j)$，我们感兴趣的是计算：
[\theta = \sum_{j = 1}^{\infty} h(x_j) P(X = x_j)]
对于一些指定的函数 $h$。通常的蒙特卡罗模拟方法是生成具有概率质量函数 $P(X = x_j)$ 的独立同分布随机向量序列 $X_1, X_2, \cdots, X_n$，然后利用大数定律，用 $\sum_{i = 1}^{n} h(X_i) / n$ 来估计 $\theta$。

但在实际中，生成具有指定概率质量函数的随机向量往往很困难，尤其是当 $X$ 是由相关随机变量组成的向量时。而且，概率质量函数有时以 $P(X = x_j) = Cb_j$ 的形式给出，其中 $b_j$ 已知，但 $C$ 需要计算，而在很多应用中，对 $b_j$ 求和来确定 $C$ 是不可行的。

幸运的是，有一种替代方法，即生成一个向量值马尔可夫链 $X_1, X_2, \cdots$，其平稳概率为 $P(X = x_j)$。如果能做到这一点，那么根据相关命题，$\sum_{i = 1}^{n} h(X_i) / n$ 仍然可以作为 $\theta$ 的估计值。

下面介绍如何生成具有任意平稳概率（可能仅指定到一个乘法常数）的马尔可夫链。设 $b(j)$（$j = 1, \cdots$）是正数，其和 $B = \sum_{j = 1}^{\infty} b(j)$ 有限。哈斯汀斯 - 梅特罗波利斯算法可以用来生成一个时间可逆的马尔可夫链，其平稳概率为：
[\pi(j) = \frac{b(j)}{B}]

算法步骤如下：
1. 选择一个任意的不可约马尔可夫转移概率矩阵 $Q$，其中 $q(i, j)$ 表示 $Q$ 的第 $i$ 行第 $j$ 列元素。
2. 当 $X_n = i$ 时，生成一个随机变量 $Y$，使得 $P(Y = j) = q(i, j)$。
3. 如果 $Y = j$，则以概率 $a(i, j)$ 将 $X_{n + 1}$ 设置为 $j$，以概率 $1 - a(i, j)$ 将其设置为 $i$，其中：
[a(i, j) = \min \left(\frac{\pi(j) q(j, i)}{\pi(i) q(i, j)}, 1\right)]

在这些条件下，状态序列构成一个马尔可夫链，其转移概率 $\tilde{P} {ij}$ 为：
[\tilde{P} {ij} = \begin{cases} q(i, j) a(i, j), & \text{如果 } j \neq i \ q(i, i) + \sum_{k \neq i} q(i, k)(1 - a(i, k)), & \text{如果 } j = i \end{cases}]

下面是哈斯汀斯 - 梅特罗波利斯算法的流程图：

graph TD;
    A[当前状态Xn = i] --> B[生成随机变量Y，P(Y = j) = q(i, j)];
    B --> C{Y = j?};
    C -- 是 --> D[计算a(i, j)];
    D --> E{以概率a(i, j)决定};
    E -- 是 --> F[Xn+1 = j];
    E -- 否 --> G[Xn+1 = i];
    C -- 否 --> G;

具体示例

• 均匀分布排列示例 ：假设要生成集合 $S$ 中均匀分布的元素，$S$ 是满足 $\sum_{j = 1}^{n} jx_j > a$ 的所有排列 $(x_1, \cdots, x_n)$ 的集合。可以利用哈斯汀斯 - 梅特罗波利斯算法，首先定义 $S$ 中元素的“相邻”概念，构建一个图，然后定义转移概率。具体来说，候选下一个状态从当前状态 $s$ 的邻居中随机选择，根据邻居的数量和一些条件决定是否转移到候选状态。如果 $n$ 较小，可以直接随机选择邻居；如果 $n$ 较大，可以扩展状态空间到所有 $n!$ 个排列。
• 吉布斯采样器 ：设 $X = (X_1, \cdots, X_n)$ 是一个随机向量，其概率质量函数 $p(x)$ 可能仅指定到一个乘法常数，我们想生成 $X$ 在 $X \in \Omega$ 条件下的条件分布的随机向量。吉布斯采样器假设对于任意 $i$（$i = 1, \cdots, n$）和值 $x_j$（$j \neq i$），可以生成具有概率质量函数 $P(X_i = x | X_j = x_j, j \neq i)$ 的随机变量 $X$。

其操作步骤如下：
1. 考虑一个状态为 $x$ 的马尔可夫链。
2. 随机选择一个坐标 $i$（$i$ 等可能地为 $1, \cdots, n$）。
3. 生成具有概率质量函数 $P(X_i = x | X_j = x_j, j \neq i)$ 的随机变量 $X$。
4. 若 $X = x$，则将状态 $y = (x_1, \cdots, x_{i - 1}, X, x_{i + 1}, \cdots, x_n)$ 作为转移的候选状态。
5. 以概率 $a(x, y) = \min \left(\frac{f(y) q(y, x)}{f(x) q(x, y)}, 1\right)$ 接受 $y$ 作为新状态，否则保持当前状态 $x$ 不变，其中 $f(x) = \frac{p(x)}{P(X \in \Omega)}$。

下面是吉布斯采样器的步骤表格：
|步骤|操作|
|----|----|
|1|选择当前状态 $x$|
|2|随机选坐标 $i$|
|3|生成 $X$ 满足 $P(X_i = x | X_j = x_j, j \neq i)$|
|4|得到候选状态 $y$|
|5|根据概率 $a(x, y)$ 决定是否更新状态|

吉布斯采样器示例

• 圆内点分布示例 ：假设要生成圆内 $n$ 个均匀分布的点，条件是任意两点之间的距离不小于 $d$。可以使用吉布斯采样器，从满足条件的 $n$ 个点开始，随机选择一个点，生成一个新的圆内点，如果该点与其他 $n - 1$ 个点的距离都不小于 $d$，则替换原来的点，否则保持不变。经过大量迭代后，这 $n$ 个点的分布近似满足要求。
• 指数分布随机向量示例 ：设 $X_i$（$i = 1, \cdots, n$）是独立的指数分布随机变量，速率为 $\lambda_i$，$S = \sum_{i = 1}^{n} X_i$。要生成 $X = (X_1, \cdots, X_n)$ 在 $S > c$ 条件下的随机向量。可以从满足 $\sum_{i = 1}^{n} x_i > c$ 且 $x_i > 0$ 的初始向量 $x$ 开始，随机选择一个坐标 $I$，生成满足条件的随机变量替换 $x_I$，如果新向量仍然满足 $S > c$，则更新状态，否则保持不变。

马尔可夫链蒙特卡罗方法为解决复杂的概率计算问题提供了有效的途径，通过构建合适的马尔可夫链，利用其平稳概率来进行估计，在实际应用中具有重要的价值。