用Python从零实现一个天气预测模型：马尔可夫链实战（附完整代码）

用Python从零实现一个天气预测模型：马尔可夫链实战（附完整代码）

天气预报看似神秘，实则背后隐藏着优雅的数学原理。想象一下，如果仅凭过去几天的天气数据，就能预测明天的天气状况，这不仅是数据科学的魅力所在，更是马尔可夫链在实际应用中的完美体现。本文将带你用Python从零构建一个天气预测模型，无需复杂的高等数学，只需基础的编程知识和一颗探索的心。

我们将使用Jupyter Notebook作为开发环境，从模拟天气数据开始，逐步构建转移概率矩阵，最终实现未来多天的天气预测。这个项目特别适合希望将概率论知识转化为实际应用的开发者，或是想通过实践理解随机过程的数据科学爱好者。你会发现，那些课本上抽象的"无记忆性"概念，在代码中变得直观而生动。

1. 环境准备与数据模拟

在开始构建模型前，我们需要准备好Python开发环境和模拟天气数据。这个基础步骤将确保后续所有分析都能顺利进行。

首先安装必要的Python库：

pip install numpy pandas matplotlib seaborn

对于天气数据，我们可以模拟一个包含三种天气状态的简单数据集：晴天(Sunny)、阴天(Cloudy)和雨天(Rainy)。以下是生成模拟数据的代码：

import numpy as np import pandas as pd # 设置随机种子保证结果可复现 np.random.seed(42) # 定义天气状态 weather_states = ['Sunny', 'Cloudy', 'Rainy'] # 生成30天的模拟天气数据 def generate_weather_data(days=30): weather_sequence = [] current_weather = np.random.choice(weather_states) # 定义简单的转移规则 transition_rules = { 'Sunny': {'Sunny': 0.7, 'Cloudy': 0.2, 'Rainy': 0.1}, 'Cloudy': {'Sunny': 0.3, 'Cloudy': 0.4, 'Rainy': 0.3}, 'Rainy': {'Sunny': 0.2, 'Cloudy': 0.5, 'Rainy': 0.3} } for _ in range(days): weather_sequence.append(current_weather) probabilities = list(transition_rules[current_weather].values()) current_weather = np.random.choice(weather_states, p=probabilities) return weather_sequence weather_data = generate_weather_data() print(f"模拟天气数据: {weather_data}")

这段代码生成了一个30天的天气序列，其中每天的天气状态取决于前一天的状态，这正是马尔可夫链的核心特性。我们使用了预设的转移概率来模拟真实世界天气变化的随机性。

提示：在实际应用中，你可以替换这部分代码，从真实天气API获取历史数据，如OpenWeatherMap或WeatherAPI提供的接口。

2. 构建转移概率矩阵

转移概率矩阵是马尔可夫链的核心，它量化了系统从一种状态转移到另一种状态的可能性。对于我们的天气模型，这个矩阵将捕捉晴天、阴天和雨天之间的转换规律。

计算转移概率矩阵的步骤如下：

1. 统计所有相邻天气状态对的出现次数
2. 将计数转换为概率分布
3. 构建规范的转移概率矩阵

以下是实现代码：

from collections import defaultdict def calculate_transition_matrix(sequence): # 初始化计数字典 transition_counts = defaultdict(lambda: defaultdict(int)) # 统计状态转移次数 for i in range(len(sequence)-1): current_state = sequence[i] next_state = sequence[i+1] transition_counts[current_state][next_state] += 1 # 转换为概率矩阵 transition_matrix = {} for current_state, transitions in transition_counts.items(): total = sum(transitions.values()) transition_matrix[current_state] = { next_state: count/total for next_state, count in transitions.items() } # 确保所有状态都有定义 for state in weather_states: if state not in transition_matrix: transition_matrix[state] = {s: 0 for s in weather_states} else: for s in weather_states: if s not in transition_matrix[state]: transition_matrix[state][s] = 0 return transition_matrix # 计算并打印转移概率矩阵 transition_matrix = calculate_transition_matrix(weather_data) print("转移概率矩阵:") pd.DataFrame(transition_matrix).T

这个矩阵的每一行代表当前天气状态，每一列代表下一天可能的天气状态，矩阵中的数值表示转移概率。例如，矩阵中"Sunny"行"Cloudy"列的值表示从晴天转为阴天的概率。

注意：当数据量较小时，某些转移可能没有出现在样本中，导致概率为0。在实际应用中，可以考虑使用拉普拉斯平滑等技术处理零概率问题。

3. 实现天气预测功能

有了转移概率矩阵，我们就可以预测未来天气了。马尔可夫链预测的基本原理是：当前状态的概率分布乘以转移矩阵，得到下一状态的概率分布。

实现多步预测的代码如下：

def predict_weather(current_state, transition_matrix, steps=1): # 初始化当前状态概率分布 current_dist = {state: 0 for state in weather_states} current_dist[current_state] = 1 predictions = [] for _ in range(steps): # 计算下一步的概率分布 next_dist = {state: 0 for state in weather_states} for current in weather_states: for next_state in weather_states: next_dist[next_state] += current_dist[current] * transition_matrix[current][next_state] # 更新当前分布 current_dist = next_dist predictions.append(current_dist) return predictions # 示例：假设今天是晴天，预测未来3天的天气概率 current_weather = 'Sunny' future_predictions = predict_weather(current_weather, transition_matrix, steps=3) for i, dist in enumerate(future_predictions, 1): print(f"第{i}天预测:") for state, prob in dist.items(): print(f" {state}: {prob:.2%}")

这段代码展示了如何从当前状态出发，通过连续应用转移矩阵，预测未来多天的天气概率分布。预测结果以概率形式呈现，反映了天气系统的不确定性。

为了更直观地理解预测结果，我们可以用matplotlib进行可视化：

import matplotlib.pyplot as plt def plot_predictions(predictions): days = len(predictions) fig, ax = plt.subplots(figsize=(10, 6)) # 准备绘图数据 states = weather_states day_numbers = range(1, days+1) prob_data = {state: [pred[state] for pred in predictions] for state in states} # 绘制堆叠柱状图 bottom = np.zeros(days) for state, probs in prob_data.items(): ax.bar(day_numbers, probs, label=state, bottom=bottom) bottom += probs ax.set_xlabel('预测天数') ax.set_ylabel('概率') ax.set_title('未来天气概率预测') ax.legend(title='天气状态') ax.set_xticks(day_numbers) ax.set_ylim(0, 1) plt.grid(True, axis='y', alpha=0.3) plt.show() plot_predictions(future_predictions)

可视化图表清晰地展示了未来几天各种天气状态的概率变化，帮助我们理解天气系统的演变趋势。

4. 模型评估与优化

构建模型只是第一步，我们需要评估其预测性能并探索优化方法。对于马尔可夫链模型，常用的评估方法包括：

• 对数似然评估：衡量模型对观测序列的拟合程度
• 预测准确率：比较预测状态与实际观测状态
• 稳定性分析：检查模型是否收敛到稳态分布

以下是实现这些评估方法的代码示例：

def calculate_log_likelihood(sequence, transition_matrix): log_likelihood = 0 for i in range(len(sequence)-1): current = sequence[i] next_state = sequence[i+1] prob = transition_matrix[current][next_state] log_likelihood += np.log(prob) if prob > 0 else -np.inf return log_likelihood # 计算模拟数据的对数似然 ll = calculate_log_likelihood(weather_data, transition_matrix) print(f"模型对数似然: {ll:.2f}") def check_steady_state(transition_matrix, tolerance=1e-6): # 将转移矩阵转换为numpy数组 states = weather_states P = np.array([[transition_matrix[i][j] for j in states] for i in states]) # 计算特征值和特征向量 eigenvalues, eigenvectors = np.linalg.eig(P.T) # 寻找接近1的特征值 steady_index = np.where(np.abs(eigenvalues - 1) < tolerance)[0] if len(steady_index) == 0: return None # 提取稳态分布 steady_vector = eigenvectors[:, steady_index[0]].real steady_vector = steady_vector / steady_vector.sum() return {state: steady_vector[i] for i, state in enumerate(states)} steady_state = check_steady_state(transition_matrix) if steady_state: print("\n稳态分布:") for state, prob in steady_state.items(): print(f" {state}: {prob:.2%}") else: print("未找到稳态分布")

模型优化可以从以下几个方面考虑：

1. 数据量扩展：使用更长时间跨度的历史数据提高统计可靠性
2. 高阶马尔可夫链：考虑前几天的天气对今天的影响
3. 季节性调整：为不同季节建立不同的转移矩阵
4. 外部因素整合：引入温度、气压等额外变量

高阶马尔可夫链的实现示例：

def build_nth_order_matrix(sequence, order=2): from itertools import product # 生成所有可能的状态组合 state_combinations = list(product(weather_states, repeat=order)) # 初始化计数 transition_counts = {comb: defaultdict(int) for comb in state_combinations} context_counts = {comb: 0 for comb in state_combinations} # 统计转移次数 for i in range(len(sequence)-order): context = tuple(sequence[i:i+order]) next_state = sequence[i+order] transition_counts[context][next_state] += 1 context_counts[context] += 1 # 计算概率 transition_matrix = {} for context in state_combinations: total = context_counts[context] if total > 0: transition_matrix[context] = { next_state: count/total for next_state, count in transition_counts[context].items() } else: transition_matrix[context] = {state: 1/len(weather_states) for state in weather_states} return transition_matrix # 构建二阶马尔可夫链转移矩阵 order_2_matrix = build_nth_order_matrix(weather_data, order=2) print("二阶转移矩阵示例:") print(order_2_matrix[('Sunny', 'Sunny')])

5. 高级应用与可视化

为了让模型更加实用，我们可以添加一些高级功能和可视化展示。状态转移图是理解马尔可夫链的绝佳工具，它能直观展示不同状态之间的转换关系。

使用networkx和matplotlib绘制状态转移图：

import networkx as nx def plot_transition_graph(transition_matrix): G = nx.DiGraph() # 添加节点和边 for source in transition_matrix: for target in transition_matrix[source]: prob = transition_matrix[source][target] if prob > 0: G.add_edge(source, target, weight=prob, label=f"{prob:.2f}") # 设置图形布局 pos = nx.spring_layout(G) # 绘制图形 plt.figure(figsize=(10, 8)) nx.draw_networkx_nodes(G, pos, node_size=2000, node_color='lightblue', alpha=0.9) nx.draw_networkx_labels(G, pos, font_size=12, font_weight='bold') # 绘制带权重的边 edges = G.edges(data=True) curved_edges = [edge for edge in edges if edge[0] == edge[1]] straight_edges = [edge for edge in edges if edge[0] != edge[1]] nx.draw_networkx_edges(G, pos, edgelist=straight_edges, width=1.5, edge_color='gray', arrowsize=20) # 绘制自环边 for edge in curved_edges: nx.draw_networkx_edges(G, pos, edgelist=[edge], width=1.5, edge_color='gray', arrowsize=20, connectionstyle=f"arc3,rad={0.3}") # 添加边标签 edge_labels = nx.get_edge_attributes(G, 'label') nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels, font_size=10, label_pos=0.3) plt.title("天气状态转移图", fontsize=15) plt.axis('off') plt.tight_layout() plt.show() plot_transition_graph(transition_matrix)

另一个有用的可视化是天气状态的时间序列分析，展示天气变化的实际模式和模型预测的对比：

def plot_weather_timeline(sequence, predictions=None): state_to_num = {state: i for i, state in enumerate(weather_states)} numeric_seq = [state_to_num[s] for s in sequence] plt.figure(figsize=(12, 6)) plt.plot(numeric_seq, 'o-', label='实际天气', markersize=8) if predictions: pred_states = [max(pred.items(), key=lambda x: x[1])[0] for pred in predictions] pred_numeric = [state_to_num[s] for s in pred_states] pred_days = range(len(sequence), len(sequence)+len(predictions)) plt.plot(pred_days, pred_numeric, 's--', label='预测天气', markersize=8) plt.yticks(range(len(weather_states)), weather_states) plt.xlabel('天数') plt.ylabel('天气状态') plt.title('天气时间序列与预测') plt.grid(True, alpha=0.3) plt.legend() plt.show() # 假设我们有新的5天观测数据 new_data = generate_weather_data(5) plot_weather_timeline(weather_data + new_data, future_predictions)

对于更复杂的分析，我们可以实现多步预测的蒙特卡洛模拟，展示预测的不确定性：

def monte_carlo_simulation(current_state, transition_matrix, steps=7, simulations=100): all_paths = [] for _ in range(simulations): path = [current_state] current = current_state for _ in range(steps): next_state = np.random.choice( weather_states, p=[transition_matrix[current][s] for s in weather_states] ) path.append(next_state) current = next_state all_paths.append(path) # 统计每天的状态分布 daily_dist = [] for day in range(steps+1): day_states = [path[day] for path in all_paths] dist = {state: day_states.count(state)/simulations for state in weather_states} daily_dist.append(dist) return daily_dist # 运行蒙特卡洛模拟 mc_results = monte_carlo_simulation('Rainy', transition_matrix) # 可视化结果 plot_predictions(mc_results)

这些可视化工具不仅增强了模型的可解释性，也为进一步优化提供了直观依据。在实际项目中，这样的分析可以帮助我们识别模型的局限性和改进方向。