企业官网建设流程全解析

用Python动画和天气模型玩转齐次马尔可夫链

想象你正在规划一周的户外活动，但当地的天气像孩子的脸说变就变——今天的晴天可能明天就转为雨天。这种"明日天气只取决于今日状态"的特性，正是齐次马尔可夫链的精髓所在。不同于传统概率论教材里令人望而生畏的矩阵运算，我们将用Python代码构建一个可视化天气预测引擎，通过动态图表展现概率如何随时间"流动"并最终达到稳定状态。你会发现，理解这个MCMC（马尔可夫链蒙特卡洛）的基石概念，就像观察天气预报图一样直观有趣。

1. 从天气模型理解马尔可夫性质

我们用一个简单的二元天气系统作为认知锚点：假设天气只有"晴"和"雨"两种状态。根据历史数据统计发现：

• 如果今天是晴天，明天有70%概率继续晴天，30%概率转为雨天
• 如果今天是雨天，明天有60%概率继续雨天，40%概率转晴

这种"未来状态仅依赖当前状态"的特性就是马尔可夫性质。用数学语言表达即：

P(明日天气|今日天气, 昨日天气,...) = P(明日天气|今日天气)

用Python字典可以直观表示这个转移规则：

weather_transition = { '晴': {'晴': 0.7, '雨': 0.3}, '雨': {'晴': 0.4, '雨': 0.6} }

齐次性体现在转移概率不随时间变化——无论是周一转到周二，还是周五转到周六，晴转雨的概率恒定为30%。这种时间无关性让模型大幅简化，也是后续平稳分布存在的前提条件。

2. 构建概率转移矩阵与状态模拟

将天气模型转化为数学矩阵，就是概率论中核心的状态转移矩阵。对于我们的二元天气系统：

这个矩阵的每个元素Pᵢⱼ表示从状态i转移到状态j的概率。用NumPy实现并模拟10天的天气变化：

import numpy as np transition_matrix = np.array([[0.7, 0.3], [0.4, 0.6]]) def simulate_weather(days, initial_state): states = ['晴', '雨'] current_state = initial_state sequence = [current_state] for _ in range(days): prob = transition_matrix[states.index(current_state)] current_state = np.random.choice(states, p=prob) sequence.append(current_state) return sequence # 从晴天开始模拟10天 print(simulate_weather(10, '晴'))

运行结果可能类似：['晴', '晴', '雨', '晴', '晴', '雨', '雨', '晴', '晴', '晴', '雨']。每次运行结果都会不同，这正是随机过程的魅力所在。

3. 可视化状态分布演变过程

单纯看状态序列难以把握整体规律，我们改用概率向量表示系统状态。假设初始时有90%概率是晴天，10%概率是雨天，初始状态向量为π₀=[0.9, 0.1]。通过矩阵乘法计算后续状态：

import matplotlib.pyplot as plt def plot_prob_evolution(initial_prob, steps): prob = initial_prob history = [prob] for _ in range(steps): prob = prob @ transition_matrix # 矩阵乘法 history.append(prob) plt.figure(figsize=(10,6)) plt.plot(np.array(history)[:,0], label='晴') plt.plot(np.array(history)[:,1], label='雨') plt.xlabel('时间步') plt.ylabel('概率') plt.legend() plt.show() plot_prob_evolution(np.array([0.9, 0.1]), 20)

观察图表会发现，无论初始概率如何，经过约10步迭代后，系统都会收敛到平稳分布π=[0.57, 0.43]。这意味着长期来看，57%的时间是晴天，43%是雨天。计算平稳分布的数学方法是求解特征向量：

eigenvalues, eigenvectors = np.linalg.eig(transition_matrix.T) stationary = eigenvectors[:, np.isclose(eigenvalues, 1)][:,0] stationary = stationary / stationary.sum() # 归一化 print(stationary) # 输出 [0.57142857 0.42857143]

4. 进阶应用与MCMC联系

这个简单天气模型揭示了MCMC方法的核心思想——通过构造具有特定平稳分布的马尔可夫链来采样复杂分布。当我们需要从某个复杂分布p(x)采样时，可以：

1. 设计转移矩阵Q，使其平稳分布等于p(x)
2. 从任意初始状态开始运行马尔可夫链
3. 经过足够长的"烧入期"(burn-in)后，收集的状态序列即为p(x)的样本

def mcmc_sampling(transition_matrix, num_samples, burn_in): samples = [] current_state = 0 # 假设从状态0开始 for i in range(burn_in + num_samples): current_state = np.random.choice( len(transition_matrix), p=transition_matrix[current_state] ) if i >= burn_in: samples.append(current_state) return samples

实际应用中，转移矩阵的设计需要满足细致平衡条件（detailed balance）：πᵢPᵢⱼ = πⱼPⱼᵢ。这保证了链的可逆性和平稳分布的存在性。虽然我们的天气模型是离散状态的，但相同原理也适用于连续状态空间，这正是Metropolis-Hastings等MCMC算法的基础。