用Python和NumPy手把手模拟马尔可夫链：从不可约到平稳分布的完整代码实战

用Python和NumPy手把手模拟马尔可夫链：从不可约到平稳分布的完整代码实战

马尔可夫链作为概率论与统计学中的重要工具，在自然语言处理、金融建模、生物信息学等领域有着广泛应用。但对于许多初学者来说，那些抽象的理论概念往往让人望而生畏。本文将用Python代码带你直观理解马尔可夫链的核心性质，通过NumPy实现从状态转移矩阵到平稳分布的全过程模拟。

1. 环境准备与基础概念

在开始编码前，我们需要确保环境配置正确。建议使用Python 3.8+版本，并安装以下依赖库：

pip install numpy matplotlib sympy

马尔可夫链的核心是状态转移矩阵，它定义了系统从一个状态转移到另一个状态的概率。用数学表示就是：

import numpy as np # 示例转移矩阵 P = np.array([ [0.2, 0.8, 0.0], [0.3, 0.5, 0.2], [0.0, 0.6, 0.4] ])

这个3×3矩阵表示：当系统处于状态1时，有20%概率保持状态1，80%概率转移到状态2，0%概率转移到状态3。矩阵的每一行必须满足概率归一化条件（总和为1）。

注意：实际应用中，转移矩阵通常通过统计历史数据获得，本文为演示使用人工构造的矩阵。

2. 不可约性验证与模拟

不可约性(Irreducibility)是指马尔可夫链中任意两个状态都可以相互到达。要验证这一点，我们可以计算转移矩阵的高次幂：

def is_irreducible(P, max_power=10): """检查马尔可夫链是否不可约""" cumulative = P.copy() for _ in range(2, max_power+1): cumulative += np.linalg.matrix_power(P, _) return np.all(cumulative > 0) # 测试不可约矩阵 P_irreducible = np.array([[0.5, 0.5], [0.3, 0.7]]) print(f"矩阵是否不可约: {is_irreducible(P_irreducible)}") # 测试可约矩阵 P_reducible = np.array([[1.0, 0.0], [0.2, 0.8]]) print(f"矩阵是否不可约: {is_irreducible(P_reducible)}")

运行结果会显示第一个矩阵不可约，而第二个可约。我们可以进一步可视化状态转移：

import matplotlib.pyplot as plt import networkx as nx def plot_markov_chain(P): G = nx.DiGraph() n = P.shape[0] for i in range(n): for j in range(n): if P[i,j] > 0: G.add_edge(f"S{i+1}", f"S{j+1}", weight=P[i,j]) pos = nx.circular_layout(G) nx.draw(G, pos, with_labels=True, node_size=800) edge_labels = nx.get_edge_attributes(G, 'weight') nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels) plt.show() plot_markov_chain(P_irreducible) plot_markov_chain(P_reducible)

3. 周期性分析与平稳分布计算

周期性(Periodicity)描述了状态返回自身的规律性。非周期链更容易收敛到平稳分布。我们可以用以下方法检测周期性：

def get_period(P, state): """计算指定状态的周期""" from math import gcd from functools import reduce n = P.shape[0] reachable_times = set() current = np.zeros(n) current[state] = 1 for t in range(1, 100): current = current @ P if current[state] > 1e-6: reachable_times.add(t) return reduce(gcd, reachable_times) # 测试非周期链 P_aperiodic = np.array([[0.5, 0.5], [0.9, 0.1]]) print(f"状态0的周期: {get_period(P_aperiodic, 0)}") # 测试周期链 P_periodic = np.array([[0, 1], [1, 0]]) print(f"状态0的周期: {get_period(P_periodic, 0)}")

平稳分布是马尔可夫链长期行为的重要特征，可以通过求解特征向量得到：

def stationary_distribution(P): """计算平稳分布""" eigvals, eigvecs = np.linalg.eig(P.T) stationary = eigvecs[:, np.isclose(eigvals, 1)] stationary = stationary[:,0].real return stationary / stationary.sum() # 计算并验证平稳分布 P_example = np.array([[0.9, 0.1], [0.4, 0.6]]) pi = stationary_distribution(P_example) print(f"平稳分布: {pi}") # 验证是否满足πP=π print(f"验证结果: {np.allclose(pi @ P_example, pi)}")

4. 可逆性验证与收敛模拟

可逆性(Reversibility)意味着链在时间反转后统计特性不变。我们可以通过细致平衡条件验证：

def is_reversible(P, pi=None): """检查马尔可夫链是否可逆""" if pi is None: pi = stationary_distribution(P) n = P.shape[0] for i in range(n): for j in range(n): if not np.isclose(pi[i]*P[i,j], pi[j]*P[j,i]): return False return True P_reversible = np.array([[0.6, 0.4], [0.2, 0.8]]) print(f"是否可逆: {is_reversible(P_reversible)}") P_non_reversible = np.array([[0.7, 0.3], [0.4, 0.6]]) print(f"是否可逆: {is_reversible(P_non_reversible)}")

最后，我们模拟马尔可夫链的收敛过程：

def simulate_markov_chain(P, initial_state, steps): """模拟马尔可夫链演化""" states = [initial_state] for _ in range(steps): new_state = np.random.choice(len(P), p=P[states[-1]]) states.append(new_state) return states # 模拟并可视化 initial = 0 states = simulate_markov_chain(P_reversible, initial, 1000) plt.figure(figsize=(12,4)) plt.plot(states[:100], 'o-') plt.xlabel('时间步') plt.ylabel('状态') plt.title('马尔可夫链前100步模拟') plt.show() # 计算状态占比 unique, counts = np.unique(states[500:], return_counts=True) # 忽略前500步 empirical = counts / counts.sum() print(f"经验分布: {empirical}") print(f"理论平稳分布: {stationary_distribution(P_reversible)}")

5. 综合案例：天气预测模型

让我们构建一个简单的天气预测模型，假设天气只有"晴"、"阴"、"雨"三种状态：

# 定义天气转移矩阵 weather_P = np.array([ [0.6, 0.3, 0.1], # 晴→晴,晴→阴,晴→雨 [0.4, 0.4, 0.2], # 阴→晴,阴→阴,阴→雨 [0.2, 0.3, 0.5] # 雨→晴,雨→阴,雨→雨 ]) # 分析模型特性 print(f"是否不可约: {is_irreducible(weather_P)}") print(f"状态0的周期: {get_period(weather_P, 0)}") weather_pi = stationary_distribution(weather_P) print(f"平稳分布: {weather_pi}") print(f"是否可逆: {is_reversible(weather_P, weather_pi)}") # 长期预测 def predict_weather(days, initial=None): if initial is None: initial = np.random.choice(3) states = simulate_markov_chain(weather_P, initial, days) return states[-1] # 预测未来7天后的天气 weather_map = {0: "晴", 1: "阴", 2: "雨"} current_weather = 0 # 今天晴天 for day in range(1, 8): current_weather = predict_weather(1, current_weather) print(f"第{day}天天气: {weather_map[current_weather]}")

这个简单模型展示了如何将马尔可夫链应用于实际问题。在实际项目中，你可能需要：

1. 基于历史数据估计转移概率
2. 增加更多天气状态（如雪、雾等）
3. 考虑季节性因素对转移矩阵的影响
4. 结合其他气象数据进行更精确的预测

6. 性能优化与实用技巧

当处理大型状态空间时，我们需要考虑计算效率。以下是几个实用技巧：

稀疏矩阵优化：对于大多数转移概率为零的情况

from scipy.sparse import csr_matrix # 创建稀疏转移矩阵 large_P = np.zeros((1000,1000)) large_P[0,1] = 0.5 large_P[0,0] = 0.5 # ...其他非零元素... sparse_P = csr_matrix(large_P) # 稀疏矩阵乘法 result = sparse_P.dot(sparse_P)

并行计算：使用多核加速矩阵运算

from multiprocessing import Pool def parallel_power(args): P, power = args return np.linalg.matrix_power(P, power) with Pool() as p: results = p.map(parallel_power, [(P,2), (P,3), (P,4)])

缓存计算结果：对于需要多次访问的中间结果

from functools import lru_cache @lru_cache(maxsize=None) def matrix_power_cached(P_tuple, power): return np.linalg.matrix_power(np.array(P_tuple), power) # 使用时先将numpy数组转为元组 P_tuple = tuple(map(tuple, P)) P_10 = matrix_power_cached(P_tuple, 10)

在处理实际问题时，我经常遇到状态空间定义不当的问题。一个常见错误是将本应区分的状态合并，导致模型精度下降。例如在用户行为分析中，将"浏览商品"和"查看详情"合并为"查看行为"可能会丢失重要信息。