别再死记硬背公式了！用Python模拟一个天气预测的马尔可夫链（附完整代码）

用Python实战马尔可夫链：从天气预测到商业决策

天气预报总是让人又爱又恨——明明说今天会下雨，结果阳光明媚；或者预测晴天，却突然倾盆大雨。但如果我们告诉你，只需几十行Python代码，就能自己建立一个简单的天气预测模型，你会不会觉得概率论突然变得有趣起来？这就是马尔可夫链的魅力所在。

1. 马尔可夫链的核心思想

想象你正在观察一个城市的天气变化。如果今天的天气只取决于昨天的天气，而与更早的天气无关，这就是典型的马尔可夫性质。这种"无记忆性"（或者说"无后效性"）是马尔可夫链的核心特征。

关键概念解析：

• 状态：系统在某一时刻的表现（如晴天、雨天）
• 转移概率：从一个状态转移到另一个状态的概率
• 状态空间：所有可能状态的集合

用数学语言来说，如果{Xₙ}是一个随机过程，满足： P(Xₙ₊₁ = j | Xₙ = i, Xₙ₋₁ = iₙ₋₁, ..., X₀ = i₀) = P(Xₙ₊₁ = j | Xₙ = i)

这就构成了一个马尔可夫链。这种简洁的性质让它成为建模各种随机系统的强大工具。

2. 构建天气预测模型

让我们用Python实现一个简单的两状态天气模型。假设天气只有"晴"和"雨"两种状态，我们可以定义转移概率矩阵：

import numpy as np # 定义转移概率矩阵 # 行表示当前状态，列表示下一状态 # 顺序：[晴，雨] weather_matrix = np.array([ [0.8, 0.2], # 今天是晴天，明天晴天的概率0.8，雨天0.2 [0.3, 0.7] # 今天是雨天，明天晴天的概率0.3，雨天0.7 ]) # 初始状态（假设第一天是晴天） current_state = 0 # 0表示晴，1表示雨

这个矩阵的每一行之和必须为1，因为系统必然会转移到某个状态。我们可以用这个模型来模拟未来几天的天气变化：

def simulate_weather(days, initial_state): states = [initial_state] for _ in range(days-1): next_state = np.random.choice( [0, 1], p=weather_matrix[states[-1]] ) states.append(next_state) return states # 模拟10天天气 weather_sequence = simulate_weather(10, current_state) print("天气序列（0=晴，1=雨）：", weather_sequence)

运行结果可能像这样：天气序列（0=晴，1=雨）： [0, 0, 1, 1, 0, 0, 0, 1, 0, 0]

3. 计算多步转移概率

马尔可夫链的一个强大特性是可以通过矩阵乘法计算n步转移概率。根据Chapman-Kolmogorov方程，n步转移矩阵等于一步转移矩阵的n次幂。

def n_step_transition(matrix, steps): return np.linalg.matrix_power(matrix, steps) # 计算3天后的转移概率 three_day_matrix = n_step_transition(weather_matrix, 3) print("3天转移概率矩阵：\n", three_day_matrix)

输出示例：

3天转移概率矩阵： [[0.662 0.338] [0.507 0.493]]

这意味着如果今天是晴天，3天后晴天的概率约为66.2%；如果今天是雨天，3天后晴天的概率约为50.7%。

4. 可视化天气预测结果

为了让预测更直观，我们可以用Matplotlib进行可视化：

import matplotlib.pyplot as plt def plot_weather_simulation(states): days = len(states) plt.figure(figsize=(10, 4)) plt.plot(range(days), states, 'o-') plt.yticks([0, 1], ['晴', '雨']) plt.xlabel('天数') plt.ylabel('天气状态') plt.title('30天天气模拟') plt.grid(True) plt.show() # 模拟并可视化30天天气 long_simulation = simulate_weather(30, current_state) plot_weather_simulation(long_simulation)

5. 马尔可夫链的稳态分布

有趣的是，许多马尔可夫链会收敛到一个稳态分布，即经过足够长的时间后，系统处于各状态的概率趋于稳定。对于我们的天气模型，可以通过求解特征向量找到这个稳态：

def find_steady_state(transition_matrix): eigenvalues, eigenvectors = np.linalg.eig(transition_matrix.T) steady_state = eigenvectors[:, np.isclose(eigenvalues, 1)] steady_state = steady_state[:, 0] / steady_state[:, 0].sum() return steady_state.real steady_dist = find_steady_state(weather_matrix) print("稳态分布：", steady_dist)

输出示例：稳态分布： [0.6 0.4]

这意味着长期来看，这个城市有60%的概率是晴天，40%的概率是雨天，无论初始天气如何。

6. 马尔可夫链的商业应用实例

马尔可夫链远不止用于天气预测。以下是一些实际应用场景：

客户行为分析：

• 将客户划分为不同状态（新客户、活跃客户、流失客户等）
• 分析客户状态间的转移概率
• 预测客户留存率和生命周期价值

# 客户状态转移矩阵示例 customer_matrix = np.array([ [0.7, 0.2, 0.1], # 新客户 → 新/活跃/流失 [0.0, 0.8, 0.2], # 活跃 → 活跃/流失 [0.1, 0.0, 0.9] # 流失 → 流失/新(重新激活) ])

库存管理：

• 预测产品需求状态（低需求、中需求、高需求）
• 优化库存策略以减少缺货或过剩

金融市场：

• 建模股票市场的状态（牛市、熊市、横盘）
• 分析信用评级的迁移概率

7. 高级应用：隐马尔可夫模型

当状态不可直接观察时，可以使用隐马尔可夫模型(HMM)。比如语音识别中，声音信号是可见的，但对应的单词或音素是隐藏状态。

from hmmlearn import hmm # 创建一个简单的HMM模型 model = hmm.GaussianHMM(n_components=2, covariance_type="diag") model.startprob_ = np.array([0.6, 0.4]) # 初始状态概率 model.transmat_ = np.array([[0.7, 0.3], [0.4, 0.6]]) # 转移矩阵 model.means_ = np.array([[0.0], [3.0]]) # 各状态的均值 model.covars_ = np.array([[0.5], [1.0]]) # 各状态的方差

8. 模型评估与优化

建立马尔可夫模型后，我们需要评估其性能：

1. 验证马尔可夫性：检查当前状态是否真的只依赖前一状态
2. 参数估计：使用最大似然估计从数据中学习转移概率
3. 模型比较：使用AIC或BIC准则选择最佳状态数

def estimate_transition_matrix(sequence, n_states): matrix = np.zeros((n_states, n_states)) for (i, j) in zip(sequence[:-1], sequence[1:]): matrix[i][j] += 1 return matrix / matrix.sum(axis=1, keepdims=True) # 从模拟数据估计转移矩阵 estimated_matrix = estimate_transition_matrix(long_simulation, 2) print("估计的转移矩阵：\n", estimated_matrix)

9. 实际应用中的注意事项

虽然马尔可夫链强大，但也有局限：

• 马尔可夫假设：现实系统可能具有更长记忆
• 状态定义：不恰当的状态划分会导致模型失效
• 数据需求：需要足够数据来可靠估计转移概率

改进方法：

• 使用高阶马尔可夫链（考虑更多历史状态）
• 结合其他模型（如马尔可夫决策过程）
• 引入半马尔可夫模型（考虑状态持续时间）

10. 完整代码示例

以下是整合了所有功能的完整天气预测系统：

import numpy as np import matplotlib.pyplot as plt from hmmlearn import hmm class WeatherPredictor: def __init__(self, sunny_to_sunny=0.8, rain_to_sunny=0.3): self.transition_matrix = np.array([ [sunny_to_sunny, 1 - sunny_to_sunny], [rain_to_sunny, 1 - rain_to_sunny] ]) def simulate(self, days, initial_state): states = [initial_state] for _ in range(days-1): next_state = np.random.choice( [0, 1], p=self.transition_matrix[states[-1]] ) states.append(next_state) return states def n_step_probability(self, steps): return np.linalg.matrix_power(self.transition_matrix, steps) def steady_state(self): eigenvalues, eigenvectors = np.linalg.eig(self.transition_matrix.T) steady = eigenvectors[:, np.isclose(eigenvalues, 1)] return (steady[:, 0] / steady[:, 0].sum()).real def plot_simulation(self, states): plt.figure(figsize=(10, 4)) plt.plot(range(len(states)), states, 'o-') plt.yticks([0, 1], ['晴', '雨']) plt.xlabel('天数') plt.ylabel('天气状态') plt.title(f'{len(states)}天天气模拟') plt.grid(True) plt.show() # 使用示例 predictor = WeatherPredictor() simulation = predictor.simulate(30, 0) predictor.plot_simulation(simulation) print("5天转移概率：\n", predictor.n_step_probability(5)) print("稳态分布：", predictor.steady_state())

这个简单的框架可以扩展到更复杂的应用场景，只需调整状态空间和转移矩阵的定义。