短期气候预测实习6

xiaoxiao啊

基于jupyter notebook
数据获取方式如下
链接：https://pan.baidu.com/s/17TSNCN-XXV5zrMXhIG9iTA
提取码：0un5


import pandas as pd
import numpy as np
from numpy.linalg import inv
import statsmodels.api as sm
import matplotlib.pyplot as plt


plt.rcParams['font.sans-serif'] = ['Times New Roman']
plt.rcParams['axes.unicode_minus'] = False


df = pd.read_csv('shixi.txt', sep='\s+', names=['y', 'x1', 'x2', 'x3', 'x4', 'x5'])

df1 = df.loc[1952:2001]
y = df1.y
X = df1.iloc[:, 1:].to_numpy()
Xd = X - X.mean(0)

df2 = df.loc[2002:2008]
y2 = df2.y
X2 = df2.iloc[:, 1:].to_numpy()
Xd2 = X2 - X2.mean(0)

b = np.dot(inv(np.dot(Xd.T, Xd)), np.dot(Xd.T, y))

def pred(var):
    var_pred = b[0] * var[:, 0] + b[1] * var[:, 1] + b[2] * var[:, 2] + b[3] * var[:, 3] + b[4] * var[:, 4]
    return var_pred
y_pred = pred(Xd)
y2_pred = pred(Xd2)


fig = plt.figure(figsize=(6, 3), facecolor='w')
ax = fig.add_subplot()
ax.plot(range(1952, 2002), y, label='y', c='r', lw=1)
ax.plot(range(1952, 2002), y_pred, label='y_pred', c='b', lw=1)
ax.set_xlabel('year')
ax.set_ylabel('RP[%]')
ax.set_title('Precip Ano Percent (1952-2001)', loc='left')
ax.legend()
ax.set_xlim(1950, 2003)
ax.set_ylim(-40, 60)
for i in ['bottom', 'top', 'left', 'right']:
    ax.spines【i】.set_linewidth(0.5)     # 记得改成英文的中括号，气象家园无法显示英文中括号加i
ax.tick_params(which='major', width=0.5, length=5)


# 使用库来做多元回归
x = sm.add_constant(Xd)  # 生成自变量
model = sm.OLS(y, x)  # 生成模型
result = model.fit()  # 模型拟合
print(result.summary())
print(result.summary())