用Python动态模拟马尔可夫链：从随机游走到平稳分布的视觉之旅

马尔可夫链作为描述状态转移的数学模型，在金融预测、自然语言处理、生物信息学等领域有着广泛应用。但对于初学者而言，那些抽象的概率转移公式和收敛性证明往往让人望而生畏。本文将通过Python代码构建一个可视化实验室，让读者能够 亲眼目睹 状态概率如何逐步演化并最终稳定——这种直观感受比任何数学推导都更能帮助理解马尔可夫链的本质。

我们将以天气预测为例：假设某地天气只有"晴天"和"雨天"两种状态，构建一个简单的齐次马尔可夫链模型。通过 numpy 进行矩阵运算，用 matplotlib 动态展示概率分布变化，最终解释特征值分解与收敛性的关系。所有代码均可直接复制到Jupyter Notebook中运行。

1. 环境准备与基础概念

在开始模拟之前，我们需要明确几个核心概念：

• 马尔可夫性质 ：系统下一时刻的状态仅取决于当前状态，与历史路径无关。就像今天的天气只受昨天影响，而与上周天气无关。
• 状态转移矩阵 ：描述各状态间转换概率的方阵，其每行元素之和必须为1。例如晴天转雨天概率为0.3，则保持晴天概率就是0.7。
• 平稳分布 ：经过足够多步转移后，系统状态概率分布趋于稳定的分布，后续转移不再改变这个分布。

先导入必要的库并设置可视化参数：

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
plt.style.use('seaborn')

# 设置中文显示
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

2. 构建天气预测模型

假设我们有以下天气转移规则：

• 晴天后有70%概率仍为晴天，30%转为雨天
• 雨天后有60%概率仍为雨天，40%转为晴天

对应的状态转移矩阵为：

当前状态 \ 下一状态	晴天	雨天
晴天	0.7	0.3
雨天	0.4	0.6

用NumPy实现这个矩阵：

P = np.array([[0.7, 0.3],  # 晴天行
              [0.4, 0.6]]) # 雨天行

验证每行和是否为1：

assert np.allclose(P.sum(axis=1), np.ones(2))

3. 模拟状态演化过程

我们从初始分布开始（假设第一天有80%概率是晴天），观察20天内的变化：

def simulate_markov_chain(transition_matrix, initial_dist, steps):
    """模拟马尔可夫链状态演化"""
    distributions = [initial_dist]
    for _ in range(steps):
        distributions.append(distributions[-1] @ transition_matrix)
    return np.array(distributions)

initial_dist = np.array([0.8, 0.2])  # 初始分布
steps = 20
distributions = simulate_markov_chain(P, initial_dist, steps)

可视化结果：

def plot_distribution_evolution(distributions):
    fig, ax = plt.subplots(figsize=(10, 5))
    ax.plot(distributions[:, 0], label='晴天概率', marker='o')
    ax.plot(distributions[:, 1], label='雨天概率', marker='s')
    ax.set_xlabel('时间步')
    ax.set_ylabel('概率')
    ax.set_title('天气状态概率演化')
    ax.legend()
    ax.grid(True)
    plt.show()

plot_distribution_evolution(distributions)

运行这段代码，你会看到两条曲线逐渐趋近于某个固定值——这就是平稳分布的视觉表现。

4. 特征值分解与收敛性分析

为什么系统会收敛？数学上这与转移矩阵的特征值密切相关。让我们分解这个矩阵：

eigenvalues, eigenvectors = np.linalg.eig(P.T)  # 注意要转置
print("特征值:", eigenvalues)
print("特征向量:\n", eigenvectors)

输出示例：

特征值: [1.   0.3]
特征向量:
 [[ 0.8       -0.70710678]
 [ 0.6        0.70710678]]

关键观察点：

1. 存在特征值1（这是所有随机矩阵的必然性质）
2. 其他特征值的绝对值小于1
3. 对应特征值1的特征向量(0.8, 0.6)经过归一化后就是平稳分布

特征值小于1的项会随着矩阵乘幂而衰减，最终只有特征值1对应的成分保留下来

我们可以通过计算验证平稳分布：

stationary_dist = eigenvectors[:, 0] / eigenvectors[:, 0].sum()
print("理论平稳分布:", stationary_dist)

5. 动态可视化收敛过程

为了让理解更直观，我们创建一个动画展示概率分布如何"寻找"平稳分布：

def create_animation(transition_matrix, initial_dist, steps):
    fig, ax = plt.subplots(figsize=(10, 5))
    ax.set_xlim(0, steps)
    ax.set_ylim(0, 1)
    ax.set_xlabel('时间步')
    ax.set_ylabel('概率')
    ax.set_title('马尔可夫链收敛过程动态演示')
    
    line1, = ax.plot([], [], label='晴天概率', lw=2)
    line2, = ax.plot([], [], label='雨天概率', lw=2)
    ax.legend()
    ax.grid(True)
    
    def init():
        line1.set_data([], [])
        line2.set_data([], [])
        return line1, line2
    
    def animate(i):
        current_dist = initial_dist @ np.linalg.matrix_power(transition_matrix, i)
        line1.set_data(range(i+1), [initial_dist[0]] + [d[0] for d in distributions[1:i+1]])
        line2.set_data(range(i+1), [initial_dist[1]] + [d[1] for d in distributions[1:i+1]])
        return line1, line2
    
    anim = FuncAnimation(fig, animate, frames=steps, init_func=init, blit=True, interval=500)
    plt.close()
    return anim

anim = create_animation(P, initial_dist, steps)
from IPython.display import HTML
HTML(anim.to_jshtml())

这个动画清晰地展示了无论从什么初始分布出发，系统最终都会稳定在相同的概率分布上。

6. 不同初始分布的比较实验

为了验证马尔可夫链的收敛性与初始状态无关，我们可以测试几个不同的初始条件：

initial_conditions = [
    [0.8, 0.2],  # 大概率晴天
    [0.5, 0.5],  # 均等概率
    [0.2, 0.8],  # 大概率雨天
    [0.99, 0.01] # 几乎确定是晴天
]

plt.figure(figsize=(10, 5))
for init_dist in initial_conditions:
    dist = simulate_markov_chain(P, np.array(init_dist), 10)
    plt.plot(dist[:, 0], label=f'初始分布={init_dist}')
plt.axhline(y=0.57, color='r', linestyle='--', label='平稳分布(晴天)')
plt.legend()
plt.title('不同初始条件下的收敛情况')
plt.xlabel('时间步')
plt.ylabel('晴天概率')
plt.grid(True)
plt.show()

所有曲线最终都收敛到相同的水平线（平稳分布），这验证了马尔可夫链的一个重要性质——收敛结果与初始状态无关。

7. MCMC的实用启示

虽然我们的天气模型很简单，但它揭示了MCMC（马尔可夫链蒙特卡洛）方法的核心思想：

1. 设计转移矩阵 ：构造一个满足特定条件的马尔可夫链
2. 确保收敛性 ：转移矩阵必须能使链收敛到平稳分布
3. 采样目标分布 ：平稳分布就是我们想要采样的目标分布

在实际应用中（如贝叶斯统计），我们经常：

# 伪代码展示MCMC典型流程
current_state = initial_state
samples = []
for _ in range(num_samples):
    # 根据转移规则提出新状态
    proposed_state = propose_new_state(current_state)  
    
    # 计算接受概率
    acceptance_prob = calculate_acceptance(current_state, proposed_state)
    
    # 决定是否接受新状态
    if random() < acceptance_prob:
        current_state = proposed_state
    
    samples.append(current_state)

关键技巧在于设计合适的转移规则，使得最终的平稳分布就是我们想要的目标分布

通过这个简单的天气模型实验，我深刻体会到可视化工具对于理解抽象数学概念的重要性。当看到那些概率曲线最终稳定在理论计算值附近时，马尔可夫链的收敛性从数学公式变成了直观感受。在实际项目中，我通常会先用这样的小模型验证理解，再应用到复杂问题中。