用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%概率保持状态180%概率转移到状态20%概率转移到状态3。矩阵的每一行必须满足概率归一化条件总和为1。注意实际应用中转移矩阵通常通过统计历史数据获得本文为演示使用人工构造的矩阵。2. 不可约性验证与模拟不可约性(Irreducibility)是指马尔可夫链中任意两个状态都可以相互到达。要验证这一点我们可以计算转移矩阵的高次幂def is_irreducible(P, max_power10): 检查马尔可夫链是否不可约 cumulative P.copy() for _ in range(2, max_power1): 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(fS{i1}, fS{j1}, weightP[i,j]) pos nx.circular_layout(G) nx.draw(G, pos, with_labelsTrue, node_size800) edge_labels nx.get_edge_attributes(G, weight) nx.draw_networkx_edge_labels(G, pos, edge_labelsedge_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, piNone): 检查马尔可夫链是否可逆 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), pP[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_countsTrue) # 忽略前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, initialNone): 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]})这个简单模型展示了如何将马尔可夫链应用于实际问题。在实际项目中你可能需要基于历史数据估计转移概率增加更多天气状态如雪、雾等考虑季节性因素对转移矩阵的影响结合其他气象数据进行更精确的预测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(maxsizeNone) 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)在处理实际问题时我经常遇到状态空间定义不当的问题。一个常见错误是将本应区分的状态合并导致模型精度下降。例如在用户行为分析中将浏览商品和查看详情合并为查看行为可能会丢失重要信息。