5.1.3 马尔可夫链

什么是马尔可夫链？

未来只依赖于现在，而不依赖于过去。

马尔可夫链是状态空间中从一个状态到另一个状态的转换的随机过程，用来描述一系列状态之间的转移概率。在转移过程中，下一个状态$S_{t+1}$只和当前状态$S_t$有关，和再之前的状态$S_{t-1}$无关。

$$P(x_{t+1}|x_t,x_{t-1},x_{t-2},...,x_t) = P(x_{t+1}|x_t)$$

马尔可夫链：状态转移、稳态分布和细致平稳条件

在马尔可夫链中，有两个核心概念：状态转移与稳态分布。

---

我们下面用一个情绪变化的例子来理解这两个概念。

假定人的情绪有三种状态：高兴、一般和悲伤。在每个时间点，人的情绪可以在这三种状态之间转换，且此刻的情绪只和上一个时刻的情绪状态有关，与再之前的情绪无关。

我们可以用状态转移矩阵$Q$来描述这种转换。

在马尔可夫链中，随着状态转移过程，状态分布会发生什么样的变化呢？

假设在$t_0$时刻的初始情绪概率分布为$p_0$，状态转移矩阵为$Q$。我们可以借助代码一步一步计算，来跟踪情绪的演化。

首先，我们定义一个情绪状态的初始分布$p_0=[0.5,0.3,0.2]$

import numpy as np
# 定义初始状态p0
e = np.array([0.5, 0.3, 0.2])

然后，我们定义一个$3\times 3$的转移矩阵$Q$，第$i$行代表从第$i$种情绪状态转移到三种情绪状态的概率

import numpy as np
# 定义状态转移矩阵Q
Q = np.array([0.9, 0.075, 0.025, 0.15, 0.8, 0.05, 0.25, 0.25, 0.5])
Q = Q.reshape(3,3) # 得到3*3矩阵 

接着，从初始状态开始，以转移矩阵进行40次状态转移：

# 状态转移
for i in range(40): 
    e = e.dot(Q)
    print(f'interation i={i}, {e[0]:.3f} {e[1]:.3f} {e[2]:.3f}\n')

我们可以得到如下结果：

interation i=0, 0.545 0.328 0.128
interation i=1, 0.572 0.335 0.094
interation i=2, 0.588 0.334 0.078
interation i=3, 0.599 0.331 0.070
interation i=4, 0.606 0.327 0.067
...
interation i=36, 0.625 0.313 0.063
interation i=37, 0.625 0.313 0.063
interation i=38, 0.625 0.313 0.063
interation i=39, 0.625 0.313 0.063

你发现了吗？经过多次状态转移，情绪状态的概率分布最后会保持不变：$p=[0.625,0.313,0.063]$。

那如果我们换一个初始条件呢？

# 定义初始状态
e = np.array([0.1, 0.1, 0.8])
# 状态转移
for i in range(40): 
    e = e.dot(Q)
    print(f'interation i={i}, {e[0]:.3f} {e[1]:.3f} {e[2]:.3f}\n')

输出：

interation i=0, 0.305 0.288 0.408
interation i=1, 0.420 0.355 0.226
interation i=2, 0.487 0.372 0.141
interation i=3, 0.530 0.369 0.101
interation i=4, 0.557 0.360 0.082
...
interation i=36, 0.625 0.313 0.063
interation i=37, 0.625 0.313 0.063
interation i=38, 0.625 0.313 0.063
interation i=39, 0.625 0.313 0.063

改变了初始条件后，经过状态转移过程，不仅分布依然稳定下来，而且得到了一模一样的分布！

---

由此，我们引出马尔可夫链的两个性质：

1. 对于初始的状态分布，如果固定了状态转移矩阵（离散状态）或状态转移函数（连续状态），那么多次迭代之后会收敛到一个最终的稳定的分布。

2. 无论初始状态分布如何，都会收敛到那个稳态分布。

也就是说，如果我们想要得到一个目标分布，只需要找到对应的状态转移矩阵，然后随机初始状态，通过不断迭代就能够达到收敛。

因此，我们最后介绍马尔可夫链的细致平稳条件：

1. 满足以下条件的状态分布$p$就是状态转移矩阵（或函数）$Q$的细致平稳分布：

$$p(x_1)*Q(x_2|x_1)=p(x_2)*Q(x_1|x_2)$$

马尔可夫链在MCMC中的作用

在MCMC中，我们旨在构造一个马尔可夫链，其稳态分布就是我们想要采样的目标后验分布。

具体来说，沿着马尔可夫链采样：

1. 每一步只需要当前状态，不需要追踪过去整个状态历史。

2. 无论目标分布是否已知，只要找到正确的转移矩阵（或函数），理论上就能收敛到目标分布。

3. 细致平稳条件确保了我们设计的转移矩阵最终能够收敛到目标分布上。

4. 可以在复杂、高维状态空间中缓慢但稳定地探索，避免维度灾难。

下面我们还是用前面提到的情绪变化的例子，示例如何用代码进行马尔可夫蒙特卡洛采样。

定义初始情绪状态：

# 定义一个初始状态
s = 0 # 3种情绪状态, s=0,1,2

设置烧录（burn-in）和采样的样本量（烧录的目的是保证取样时分布已经稳定，所以一般会取比较大的数）：

burnin = 1000000 # 被丢弃的burn-in samples
nSample = 1000000
samples = []

进行状态转移：

for i in range(burnin + nSample):
    u = np.random.rand()
    match s:
        case 0:
            if u <= 0.9:
                s = 0
            elif (0.9<u) & (u <= 0.9+0.075):
                s = 1
            elif 0.9+0.075 < u:
                s = 2
        case 1:
            if u <=0.15:
                s = 0
            elif (0.15<u) & (u <= 0.15+0.8):
                s = 1
            elif 0.15+0.8 < u:
                s = 2
        case 2:
            if u <=0.25:
                s = 0
            elif (0.25<u) & (u <= 0.25+0.5):
                s = 1
            elif 0.25+0.5 < u:
                s = 2
    samples.append(s)

得到采样分布：

bins = [ -0.5, 0.5, 1.5, 2.5 ]
a, _, _ = plt.hist(samples[burnin:], bins=bins)
plt.xlabel('State')
plt.ylabel('Count')
print(a/nSample)

采样得到的分布接近马尔可夫链的稳态分布[0.625, 0.313, 0.062]：

[0.617709 0.339236 0.043055]