用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)

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