背景与目标

每天变化多端的天气让我们的生活从满了未知的可能，极大的影响了我们的生活
天气预测在许多领域中发挥着至关重要的作用，从农业到交通运输再到灾害管理，准确的天气预报能够帮助各行业提前应对可能的天气变化，减少损失，优化资源分配。随着数据科学和机器学习技术的发展，基于历史天气数据的预测模型已经成为现代天气预报的重要工具。通过利用大量的历史气象数据，构建和训练机器学习模型，可以更精确地预测未来的天气情况。

本项目的目标是构建和评估多种机器学习模型，以预测某些地点是否会在第二天降雨（RainTomorrow）。我们使用了一系列气象特征数据，例如最高气温、最低气温、降雨量、风速、湿度、气压等，通过不同的机器学习算法进行训练和测试。具体目标包括：

  1. 数据预处理与标准化：清洗和规范化气象数据，以确保模型训练的有效性。
  2. 模型构建与训练：使用多种机器学习算法，包括决策树、随机森林、逻辑回归、K近邻、支持向量机和神经网络，构建预测模型。
  3. 模型评估与比较：通过多种评估指标（如准确率、精确率、召回率、F1值和AUC值）评估各模型的性能，找出最优的预测模型。
  4. 结果可视化：使用图表展示各地点的模型性能，以便直观地比较不同模型在不同地点的预测效果。
  5. 进一步优化与改进：根据评估结果，对模型进行优化，以提高预测的准确性和可靠性。
    通过实现上述目标，我们希望能够构建出一个高效、可靠的天气预测系统，为相关领域提供有价值的参考。

数据说明

列名 含义
Date 日期
Location 地点
MinTemp 最低气温（摄氏度）
MaxTemp 最高气温（摄氏度）
Rainfall 降雨量（毫米）
Evaporation 蒸发量（毫米）
Sunshine 日照时数（小时）
WindGustDir 最大阵风的风向
WindGustSpeed 最大阵风的风速（公里/小时）
WindDir9am 上午9点的风向
WindDir3pm 下午3点的风向
WindSpeed9am 上午9点的风速（公里/小时）
WindSpeed3pm 下午3点的风速（公里/小时）
Humidity9am 上午9点的湿度（百分比）
Humidity3pm 下午3点的湿度（百分比）
Pressure9am 上午9点的气压（百帕）
Pressure3pm 下午3点的气压（百帕）
Cloud9am 上午9点的云量（八分制）
Cloud3pm 下午3点的云量（八分制）
Temp9am 上午9点的温度（摄氏度）
Temp3pm 下午3点的温度（摄氏度）
RainToday 是否在今天下雨（是/否）
RISK_MM 明天降雨的风险量（毫米）
RainTomorrow 是否在明天下雨（是/否）

预处理

首先我们对数据进行读取以及预处理

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.tsa.stattools import adfuller, acf, pacf
from statsmodels.tsa.arima.model import ARIMA
from sklearn.metrics import mean_squared_error
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
# 读取数据
df = pd.read_csv('weatherAUS.csv')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']  # 使用黑体
plt.rcParams['axes.unicode_minus'] = False  # 正常显示负号

# 按日期排序
# df = df.sort_values(by='Date')
df['Date'] = pd.to_datetime(df['Date'])
# 设置日期为索引
df.set_index('Date', inplace=True)

显而易见的，我要预测一个数据的最高温度的数据，以及地点，因为每个地点对应的时间是唯一的 ，如果不分类地点的话
就会出现一个时间有很多数据最高温度的 情况
之前没分类 就走错了很多弯路 导致步骤对了反而没有成功的数据

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
# 只保留需要的列
df = df[['Location', 'MaxTemp']]

# 去除缺失值
df = df.dropna()

# 选择一个地点进行分析，比如 "Albury"
location = 'Albury'
location_data = df[df['Location'] == location]['MaxTemp']

# 确保时间序列数据具有频率
location_data = location_data.asfreq('D')

# 使用插值法填补缺失值
location_data = location_data.interpolate(method='linear')

# 删除依然存在的缺失值
location_data = location_data.dropna()

1
2
3
4
5
# 绘制时间序列图
plt.figure(figsize=(10, 6))
plt.plot(location_data)
plt.title(f'{location} 的最高温度时间序列')
plt.show()

png

显示最高温度的时间序列，我们发现时间序列图中对于这个时间序列跨度很大，，序列是按照以年份来周期循环
很平稳且有规律进行周期在一定的范围内稳定的波动，并不是在某个常量值附近波动，根据时序图检验平稳性的原理来判断，该序列属于非平稳的随机序列

1
2
3
4
5
6
7
8
# ADF平稳性检验
adf_result = adfuller(location_data)
print('ADF 检验结果:', adf_result)

# 分解时间序列，指定周期，比如一年365天
decomposition = sm.tsa.seasonal_decompose(location_data, model='additive', period=365)
fig = decomposition.plot()
plt.show()
1
ADF 检验结果: (-3.3666814947988617, 0.012148512496123455, 12, 3116, {'1%': -3.432450351481275, '5%': -2.862468004787272, '10%': -2.567263999190068}, 15487.167404620237)

png

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
# ADF平稳性检验
adf_result = adfuller(location_data)
print('ADF 检验结果:', adf_result)
print(f'ADF统计量: {adf_result[0]}')
print(f'p值: {adf_result[1]}')
print(f'使用的滞后数: {adf_result[2]}')
print(f'样本量: {adf_result[3]}')
print('临界值:')
for key, value in adf_result[4].items():
    print(f'    {key}: {value}')

if adf_result[1] < 0.05:
    print("时间序列是平稳的")
else:
    print("时间序列不是平稳的")

# 绘制自相关函数（ACF）图
plt.figure(figsize=(10, 6))
plot_acf(location_data, lags=50)
plt.title(f'{location} 的自相关函数（ACF）')
plt.show()

1
2
3
4
5
6
7
8
9
10
ADF 检验结果: (-3.3666814947988617, 0.012148512496123455, 12, 3116, {'1%': -3.432450351481275, '5%': -2.862468004787272, '10%': -2.567263999190068}, 15487.167404620237)
ADF统计量: -3.3666814947988617
p值: 0.012148512496123455
使用的滞后数: 12
样本量: 3116
临界值:
    1%: -3.432450351481275
    5%: -2.862468004787272
    10%: -2.567263999190068
时间序列是平稳的

png

从自相关图看出，自相关系数在X轴 的上方 ，并不是在x轴上下波动，是单调趋势序列的典型特征，具有持久的正自相关性，且数值基本无变化，意味着时间序列中存在某种周期性或趋势。

这个序列是有规律的以年为周期进行波动

根据ADF检验结果，p值为0.0121，小于0.05，因此时间序列是平稳的。即使ACF图显示自相关系数从1.0慢慢降低到0.6，这也可能表示该时间序列中存在一些长期依赖性或季节性。

尽管ADF检验表明时间序列是平稳的，但是ACF图的结果可以提供更多信息。ACF图显示自相关系数从1.0逐渐降低到0.6，表明存在一些长期依赖性，这可能是由于季节性或其他长期趋势引起的。

解释与处理建议
长期依赖性：ACF图显示的长期依赖性表明时间序列中存在长期趋势或周期性成分。即使时间序列整体上是平稳的，但这些成分可能需要进一步的处理或建模。

季节性：如果时间序列中存在季节性成分，可以考虑使用季节性差分来去除季节性。例如，对于月度数据，可以进行12期差分。

模型选择：在进行时间序列建模时，可以考虑使用带有季节性成分的模型，如SARIMA模型。

所以我们要进行差分来处理数据，由于数据是一年循环波动，那么属于季节周期性
所以使用季节差分，以一年365为参数进行降维

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
# 进行季节性差分
seasonal_diff = location_data.diff(365).dropna()

# 绘制差分后的时间序列图
plt.figure(figsize=(10, 6))
plt.plot(seasonal_diff)
plt.title(f'{location} 的最高温度季节性差分时间序列')
plt.show()

# 绘制差分后的自相关函数（ACF）图
plt.figure(figsize=(18, 10))
plot_acf(seasonal_diff, lags=50)
plt.title(f'{location} 的最高温度季节性差分自相关函数（ACF）')
plt.show()

# 绘制差分后的偏自相关函数（PACF）图
plt.figure(figsize=(10, 6))
plot_pacf(seasonal_diff, lags=50)
plt.title(f'{location} 的最高温度季节性差分偏自相关函数（PACF）')
plt.show()

png

png

png 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16



# 季节性差分后的ADF检验
result = adfuller(seasonal_diff.dropna())
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
print('Critical Values:')
for key, value in result[4].items():
    print('\t%s: %.3f' % (key, value))

# 结果解读
if result[1] > 0.05:
    print("季节性差分后的数据仍然是非平稳的")
else:
    print("季节性差分后的数据是平稳的")
1
2
3
4
5
6
7
ADF Statistic: -11.476252
p-value: 0.000000
Critical Values:
	1%: -3.433
	5%: -2.863
	10%: -2.567
季节性差分后的数据是平稳的
1
2
3
4
5
6
7
from statsmodels.stats.diagnostic import acorr_ljungbox
# 去除NaN值
max_temp_diff = seasonal_diff.dropna()


# 白噪声检验
print("差分序列的白噪声检验结果：\n" , acorr_ljungbox(max_temp_diff, lags=1))
1
2
3
差分序列的白噪声检验结果：
        lb_stat      lb_pvalue
1  1056.822391  8.006870e-232

模型预测

找出最佳参数

1
2
3
4
5
6
7
8
9
10
# 选择ARIMA模型参数
# max_temp_diff = location_data.dropna()



max_temp_diff  = seasonal_diff.dropna()


order_select = sm.tsa.arma_order_select_ic(max_temp_diff, ic=['bic'], trend='n', max_ar=5, max_ma=5)
print('BIC选择的最佳参数:', order_select.bic_min_order)
BIC选择的最佳参数: (1, 4)

BIC选择的最佳参数: (1, 4)

1
2
3
4
# 拟合ARIMA模型
model = ARIMA(df['MaxTemp'], order=(order_select.bic_min_order[0], 0, order_select.bic_min_order[1]))
model_fit = model.fit()
print(model_fit.summary())
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
                               SARIMAX Results                                
==============================================================================
Dep. Variable:                MaxTemp   No. Observations:               141871
Model:                 ARIMA(1, 0, 4)   Log Likelihood             -366137.686
Date:                Tue, 02 Jul 2024   AIC                         732289.372
Time:                        04:42:15   BIC                         732358.411
Sample:                             0   HQIC                        732310.000
                             - 141871                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         23.2269      0.377     61.576      0.000      22.488      23.966
ar.L1          0.9950      0.000   3094.406      0.000       0.994       0.996
ma.L1         -0.3711      0.002   -178.084      0.000      -0.375      -0.367
ma.L2         -0.3006      0.002   -127.512      0.000      -0.305      -0.296
ma.L3         -0.1166      0.002    -48.575      0.000      -0.121      -0.112
ma.L4         -0.0158      0.002     -7.138      0.000      -0.020      -0.011
sigma2        10.2129      0.028    369.185      0.000      10.159      10.267
===================================================================================
Ljung-Box (L1) (Q):                   0.03   Jarque-Bera (JB):             30952.95
Prob(Q):                              0.87   Prob(JB):                         0.00
Heteroskedasticity (H):               1.04   Skew:                            -0.15
Prob(H) (two-sided):                  0.00   Kurtosis:                         5.27
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
# 拆分训练集和测试集
train_data = seasonal_diff[:-365]
test_data = seasonal_diff[-365:]

# 训练ARIMA模型
model = ARIMA(train_data, order=(2, 0, 5))
model_fit = model.fit()

# 预测
forecast_result = model_fit.get_forecast(steps=len(test_data))
forecast = forecast_result.predicted_mean
conf_int = forecast_result.conf_int()

# 打印预测结果
print(f'预测未来30天的最高温度: {forecast}')
print(f'预测的置信区间: {conf_int}')

# 绘制预测结果
plt.figure(figsize=(10, 6))
plt.plot(train_data.index, train_data, label='训练数据')
plt.plot(test_data.index, test_data, label='实际值')
plt.plot(test_data.index, forecast, label='预测值')
plt.fill_between(test_data.index, conf_int.iloc[:, 0], conf_int.iloc[:, 1], color='pink', alpha=0.3)
plt.legend()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
预测未来30天的最高温度: 2016-06-26   -1.574331
2016-06-27   -0.481513
2016-06-28    0.067825
2016-06-29    0.390562
2016-06-30    0.363980
                ...   
2017-06-21    0.144490
2017-06-22    0.144482
2017-06-23    0.144474
2017-06-24    0.144466
2017-06-25    0.144458
Freq: D, Name: predicted_mean, Length: 365, dtype: float64
预测的置信区间:             lower MaxTemp  upper MaxTemp
2016-06-26      -9.334773       6.186111
2016-06-27      -9.752366       8.789340
2016-06-28      -9.567745       9.703394
2016-06-29      -9.345149      10.126273
2016-06-30      -9.393319      10.121280
...                   ...            ...
2017-06-21      -9.795905      10.084886
2017-06-22      -9.795913      10.084877
2017-06-23      -9.795922      10.084869
2017-06-24      -9.795930      10.084861
2017-06-25      -9.795938      10.084853

[365 rows x 2 columns]

png

1
2
3
4
5
6
7
8
9
10
11
# 只显示最后90天的数据
last_360_days = seasonal_diff[-700:]

# 绘制预测结果
plt.figure(figsize=(10, 6))
plt.plot(last_360_days.index, last_360_days, label='实际值', color='blue')
plt.plot(test_data.index, forecast, label='预测值', color='red')
plt.fill_between(test_data.index, conf_int.iloc[:, 0], conf_int.iloc[:, 1], color='pink', alpha=0.3)
plt.title(f'{location} 的最高温度预测 (最后360天)')
plt.legend()
plt.show()

png

返回原本的值

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
# 提取原始的最高温度数据
original_data = df[df['Location'] == location]['MaxTemp']

# 确保时间序列数据具有频率
original_data = original_data.asfreq('D')

# 使用插值法填补缺失值
original_data = original_data.interpolate(method='linear')

# 删除依然存在的缺失值
original_data = original_data.dropna()

# 获取季节性差分的基础值
seasonal_base = original_data.shift(365).dropna()

# 将差分后的预测值还原为原始的时间序列数据
restored_forecast = forecast + seasonal_base[-len(forecast):]

# 绘制原始数据和还原后的预测值的对比图
plt.figure(figsize=(12, 8))
plt.plot(original_data.index, original_data, label='原始数据', color='blue')
plt.plot(restored_forecast.index, restored_forecast, label='还原后的预测值', color='red')
plt.fill_between(test_data.index, conf_int.iloc[:, 0] + seasonal_base[-len(forecast):],
                 conf_int.iloc[:, 1] + seasonal_base[-len(forecast):], color='pink', alpha=0.3)
plt.title(f'{location} 的最高温度原始数据与还原后的预测值对比')
plt.legend()
plt.show()

# 打印还原后的预测值
print(f'还原后的预测值: {restored_forecast}')
print(f'预测的置信区间: {conf_int}')

png

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
还原后的预测值: 2016-06-26    11.625669
2016-06-27    11.918487
2016-06-28    15.067825
2016-06-29     8.990562
2016-06-30    10.963980
                ...    
2017-06-21    11.744490
2017-06-22    13.844482
2017-06-23    13.344474
2017-06-24     8.344466
2017-06-25    10.644458
Freq: D, Length: 365, dtype: float64
预测的置信区间:             lower MaxTemp  upper MaxTemp
2016-06-26      -9.334773       6.186111
2016-06-27      -9.752366       8.789340
2016-06-28      -9.567745       9.703394
2016-06-29      -9.345149      10.126273
2016-06-30      -9.393319      10.121280
...                   ...            ...
2017-06-21      -9.795905      10.084886
2017-06-22      -9.795913      10.084877
2017-06-23      -9.795922      10.084869
2017-06-24      -9.795930      10.084861
2017-06-25      -9.795938      10.084853

[365 rows x 2 columns]

发现测试结果和 实际结果相差不大

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
# 提取原始的最高温度数据
original_data = df[df['Location'] == location]['MaxTemp']

# 确保时间序列数据具有频率
original_data = original_data.asfreq('D')

# 使用插值法填补缺失值
original_data = original_data.interpolate(method='linear')

# 删除依然存在的缺失值
original_data = original_data.dropna()

# 获取季节性差分的基础值
seasonal_base = original_data.shift(365).dropna()

# 将差分后的预测值还原为原始的时间序列数据
restored_forecast = forecast + seasonal_base[-len(forecast):]

# 设置时间范围，显示测试集的最后一段时间
last_days = 365 # 可以根据需要调整
time_range = original_data.index[-last_days:]

# 绘制原始数据和还原后的预测值的详细对比图
plt.figure(figsize=(14, 10))
plt.plot(original_data[time_range], label='原始数据', color='blue', marker='o', markersize=3)
plt.plot(restored_forecast[time_range], label='还原后的预测值', color='red', linestyle='--', marker='x', markersize=3)
plt.fill_between(time_range, conf_int.iloc[:, 0][time_range] + seasonal_base[-len(forecast):][time_range],
                 conf_int.iloc[:, 1][time_range] + seasonal_base[-len(forecast):][time_range], color='pink', alpha=0.3)
plt.title(f'{location} 的最高温度原始数据与还原后的预测值对比 (详细)')
plt.xlabel('日期')
plt.ylabel('最高温度')
plt.grid(True)
plt.legend()
plt.show()

# 打印还原后的预测值
print(f'还原后的预测值: {restored_forecast}')
print(f'预测的置信区间: {conf_int}')

png

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
还原后的预测值: 2016-06-26    11.625669
2016-06-27    11.918487
2016-06-28    15.067825
2016-06-29     8.990562
2016-06-30    10.963980
                ...    
2017-06-21    11.744490
2017-06-22    13.844482
2017-06-23    13.344474
2017-06-24     8.344466
2017-06-25    10.644458
Freq: D, Length: 365, dtype: float64
预测的置信区间:             lower MaxTemp  upper MaxTemp
2016-06-26      -9.334773       6.186111
2016-06-27      -9.752366       8.789340
2016-06-28      -9.567745       9.703394
2016-06-29      -9.345149      10.126273
2016-06-30      -9.393319      10.121280
...                   ...            ...
2017-06-21      -9.795905      10.084886
2017-06-22      -9.795913      10.084877
2017-06-23      -9.795922      10.084869
2017-06-24      -9.795930      10.084861
2017-06-25      -9.795938      10.084853

[365 rows x 2 columns]

查看更详细的图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
# 提取原始的最高温度数据
original_data = df[df['Location'] == location]['MaxTemp']

# 确保时间序列数据具有频率
original_data = original_data.asfreq('D')

# 使用插值法填补缺失值
original_data = original_data.interpolate(method='linear')

# 删除依然存在的缺失值
original_data = original_data.dropna()

# 获取季节性差分的基础值
seasonal_base = original_data.shift(365).dropna()

# 将差分后的预测值还原为原始的时间序列数据
restored_forecast = forecast + seasonal_base[-len(forecast):]

# 设置时间范围，显示测试集的最后一段时间
last_days = 30 # 可以根据需要调整
time_range = original_data.index[-last_days:]

# 绘制原始数据和还原后的预测值的详细对比图
plt.figure(figsize=(14, 10))
plt.plot(original_data[time_range], label='原始数据', color='blue', marker='o', markersize=3)
plt.plot(restored_forecast[time_range], label='还原后的预测值', color='red', linestyle='--', marker='x', markersize=3)
plt.fill_between(time_range, conf_int.iloc[:, 0][time_range] + seasonal_base[-len(forecast):][time_range],
                 conf_int.iloc[:, 1][time_range] + seasonal_base[-len(forecast):][time_range], color='pink', alpha=0.3)
plt.title(f'{location} 的最高温度原始数据与还原后的预测值对比 (详细)')
plt.xlabel('日期')
plt.ylabel('最高温度')
plt.grid(True)
plt.legend()
plt.show()

# 打印还原后的预测值
print(f'还原后的预测值: {restored_forecast}')
print(f'预测的置信区间: {conf_int}')

png

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
还原后的预测值: 2016-06-26    11.625669
2016-06-27    11.918487
2016-06-28    15.067825
2016-06-29     8.990562
2016-06-30    10.963980
                ...    
2017-06-21    11.744490
2017-06-22    13.844482
2017-06-23    13.344474
2017-06-24     8.344466
2017-06-25    10.644458
Freq: D, Length: 365, dtype: float64
预测的置信区间:             lower MaxTemp  upper MaxTemp
2016-06-26      -9.334773       6.186111
2016-06-27      -9.752366       8.789340
2016-06-28      -9.567745       9.703394
2016-06-29      -9.345149      10.126273
2016-06-30      -9.393319      10.121280
...                   ...            ...
2017-06-21      -9.795905      10.084886
2017-06-22      -9.795913      10.084877
2017-06-23      -9.795922      10.084869
2017-06-24      -9.795930      10.084861
2017-06-25      -9.795938      10.084853

[365 rows x 2 columns]

2所有地区

全部模型预测&详细信息

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.tsa.stattools import adfuller, acf, pacf
from statsmodels.tsa.arima.model import ARIMA
from sklearn.metrics import mean_squared_error
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.stats.diagnostic import acorr_ljungbox

# 读取数据
df = pd.read_csv('weatherAUS.csv')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']  # 使用黑体
plt.rcParams['axes.unicode_minus'] = False  # 正常显示负号

# 按日期排序
df['Date'] = pd.to_datetime(df['Date'])
df.set_index('Date', inplace=True)

# 去除缺失值
df = df.dropna(subset=['MaxTemp'])

# 获取所有地点的列表
locations = df['Location'].unique()

# 创建一个字典来存储预测结果
forecasts = {}

# 遍历每个地点
for location in locations:
    print(f"正在处理地点: {location}")
    
    # 获取特定地点的数据
    location_data = df[df['Location'] == location]['MaxTemp']
    
    # 确保时间序列数据具有频率
    location_data = location_data.asfreq('D')
    
    # 使用插值法填补缺失值
    location_data = location_data.interpolate(method='linear')
    
    # 删除依然存在的缺失值
    location_data = location_data.dropna()
    
    # 绘制时间序列图
    plt.figure(figsize=(14, 6))
    plt.plot(location_data)
    plt.title(f'{location} 的最高温度时间序列')
    plt.xlabel('日期')
    plt.ylabel('最高温度')
    plt.grid(True)
    plt.show()
    
    # ADF平稳性检验
    adf_result = adfuller(location_data)
    print('ADF 检验结果:', adf_result)
    print(f'ADF统计量: {adf_result[0]}')
    print(f'p值: {adf_result[1]}')
    print(f'使用的滞后数: {adf_result[2]}')
    print(f'样本量: {adf_result[3]}')
    print('临界值:')
    for key, value in adf_result[4].items():
        print(f'    {key}: {value}')
    
    if adf_result[1] < 0.05:
        print("时间序列是平稳的")
    else:
        print("时间序列不是平稳的")
    
    # 分解时间序列，指定周期，比如一年365天
    decomposition = sm.tsa.seasonal_decompose(location_data, model='additive', period=365)
    fig = decomposition.plot()
    plt.show()
    
    # 绘制自相关函数（ACF）图
    plt.figure(figsize=(14, 6))
    plot_acf(location_data, lags=50)
    plt.title(f'{location} 的自相关函数（ACF）')
    plt.xlabel('滞后数')
    plt.ylabel('自相关')
    plt.grid(True)
    plt.show()
    
    # 进行季节性差分
    seasonal_diff = location_data.diff(365).dropna()
    
    # 绘制差分后的时间序列图
    plt.figure(figsize=(14, 6))
    plt.plot(seasonal_diff)
    plt.title(f'{location} 的最高温度季节性差分时间序列')
    plt.xlabel('日期')
    plt.ylabel('差分后的最高温度')
    plt.grid(True)
    plt.show()
    
    # 绘制差分后的自相关函数（ACF）图
    plt.figure(figsize=(14, 6))
    plot_acf(seasonal_diff, lags=50)
    plt.title(f'{location} 的最高温度季节性差分自相关函数（ACF）')
    plt.xlabel('滞后数')
    plt.ylabel('自相关')
    plt.grid(True)
    plt.show()
    
    # 绘制差分后的偏自相关函数（PACF）图
    plt.figure(figsize=(14, 6))
    plot_pacf(seasonal_diff, lags=50)
    plt.title(f'{location} 的最高温度季节性差分偏自相关函数（PACF）')
    plt.xlabel('滞后数')
    plt.ylabel('偏自相关')
    plt.grid(True)
    plt.show()
    
    # 季节性差分后的ADF检验
    result = adfuller(seasonal_diff.dropna())
    print('ADF Statistic: %f' % result[0])
    print('p-value: %f' % result[1])
    print('Critical Values:')
    for key, value in result[4].items():
        print('\t%s: %.3f' % (key, value))
    
    # 结果解读
    if result[1] > 0.05:
        print("季节性差分后的数据仍然是非平稳的")
    else:
        print("季节性差分后的数据是平稳的")
    
    # 白噪声检验
    print("差分序列的白噪声检验结果：\n", acorr_ljungbox(seasonal_diff.dropna(), lags=1))
    
    # 选择ARIMA模型参数
    order_select = sm.tsa.arma_order_select_ic(seasonal_diff.dropna(), ic=['bic'], trend='n', max_ar=5, max_ma=5)
    print('BIC选择的最佳参数:', order_select.bic_min_order)
    
    # 拆分训练集和测试集
    train_data = seasonal_diff[:-365]
    test_data = seasonal_diff[-365:]
    
    # 训练ARIMA模型
    model = ARIMA(train_data, order=(order_select.bic_min_order[0], 0, order_select.bic_min_order[1]))
    model_fit = model.fit()
    
    # 预测
    forecast_result = model_fit.get_forecast(steps=len(test_data))
    forecast = forecast_result.predicted_mean
    conf_int = forecast_result.conf_int()
    
    # 获取季节性差分的基础值
    seasonal_base = location_data.shift(365).dropna()
    
    # 将差分后的预测值还原为原始的时间序列数据
    restored_forecast = forecast + seasonal_base[-len(forecast):]
    
    # 存储预测结果
    forecasts[location] = {
        'forecast': restored_forecast,
        'conf_int': conf_int,
        'train_data': train_data,
        'test_data': test_data,
        'original_data': location_data
    }
    
    # 绘制原始数据和还原后的预测值的对比图（所有数据）
    plt.figure(figsize=(14, 10))
    plt.plot(location_data.index, location_data, label='原始数据', color='blue', marker='o', markersize=3)
    plt.plot(restored_forecast.index, restored_forecast, label='还原后的预测值', color='red', linestyle='--', marker='x', markersize=3)
    plt.fill_between(test_data.index, conf_int.iloc[:, 0] + seasonal_base[-len(forecast):],
                     conf_int.iloc[:, 1] + seasonal_base[-len(forecast):], color='pink', alpha=0.3)
    plt.title(f'{location} 的最高温度原始数据与还原后的预测值对比（所有数据）')
    plt.xlabel('日期')
    plt.ylabel('最高温度')
    plt.grid(True)
    plt.legend()
    plt.show()
    
    # 绘制原始数据和还原后的预测值的对比图（最后700天的数据）
    plt.figure(figsize=(14, 10))
    plt.plot(location_data.index[-700:], location_data[-700:], label='原始数据', color='blue', marker='o', markersize=3)
    plt.plot(restored_forecast.index[-365:], restored_forecast, label='还原后的预测值', color='red', linestyle='--', marker='x', markersize=3)
    plt.fill_between(test_data.index[-365:], conf_int.iloc[:, 0] + seasonal_base[-365:],
                     conf_int.iloc[:, 1] + seasonal_base[-365:], color='pink', alpha=0.3)
    plt.title(f'{location} 的最高温度原始数据与还原后的预测值对比（最后700天）')
    plt.xlabel('日期')
    plt.ylabel('最高温度')
    plt.grid(True)
    plt.legend()
    plt.show()
    
    plt.figure(figsize=(14, 10))
    plt.plot(location_data.index[-30:], location_data[-30:], label='原始数据', color='blue', marker='o', markersize=3)
    plt.plot(restored_forecast.index[-30:], restored_forecast[-30:], label='还原后的预测值', color='red', linestyle='--', marker='x', markersize=3)
    plt.fill_between(test_data.index[-30:], conf_int.iloc[:, 0].iloc[-30:] + seasonal_base.iloc[-30:],
                     conf_int.iloc[:, 1].iloc[-30:] + seasonal_base.iloc[-30:], color='pink', alpha=0.3)
    plt.title(f'{location} 的最高温度原始数据与还原后的预测值对比（最后30天）')
    plt.xlabel('日期')
    plt.ylabel('最高温度')
    plt.grid(True)
    plt.legend()
    plt.show()
    
    
    
    # 打印还原后的预测值
    print(f'还原后的预测值: {restored_forecast}')
    print(f'预测的置信区间: {conf_int}')

# 打印所有地点的预测结果
for location, result in forecasts.items():
    print(f"地点: {location}")
    print(f"还原后的预测值: {result['forecast']}")
    print(f"预测的置信区间: {result['conf_int']}")
    print("\n")

全部地区图片太多不舍得流量钱和COS存储详细请看github

​
​