本次作業使用豐原站的觀測記錄，分成 train set 跟 test set，train set是豐原站每個月的前20天所有資料，test set 則是從豐原站剩下的資料中取樣出來。
train.csv: 每個月前20天每個小時的氣象資料(共24小时，每小時有18種測資)。共12個月。
test.csv: 從剩下的資料當中取樣出連續的10小時為一筆，前九小時的所有觀測數據當作 feature，第十小時的PM2.5當作 answer。一共取出 240 筆不重複的 test data，請根據 feauure 預測這240筆的PM2.5。
数据预处理
import numpy as np
import pandas as pd
# 导入训练数据
df_train = pd.read_csv('train.csv')
df_train.head(10)
可以看出在数据集合中有很多其他的观测变量，这里我们的任务主要是预测 PM2.5，所以将无关的记录剔除，同时前三列的数据也是与任务无关的。
df_train = df_train[df_train['observation'] == 'PM2.5']
df_train = df_train.iloc[:,3:]
df_train
最终得到的训练数据如下：
接着看一看测试数据：
# 导入测试数据
df_test = pd.read_csv('./test.csv')
df_test.head(10)
我们的测试数据是连续9个小时的PM2.5观测值，作为预测任务的 features，第 10 个小时的值为需要预测的 label。
所以需要对训练数据进一步处理，选择连续的9小时数据作为feature，第十个小时的观测数据作为预测，这样我们需要选择连续的10列数据，总共 24 列，那么要 24-10+1 ，即 15 次。
train_X = []
train_y = []
for i in range(24-10+1):
x = df_train.iloc[:, i:i+9]
x.columns=np.array(range(9))
y = df_train.iloc[:, i+9]
y.columns=np.array(range(1))
train_X.append(x)
train_y.append(y)
train_X=pd.concat(train_y) pd.concat(train_X)
train_y=pd.concat(train_y)
构建模型
对于时序预测问题，可以将前九个小时的值看作是特征，建立多元线性回归模型。
class LinearRegression:
def __init__(self):
self.coef_ = None
self.intercept_ = None
self._theta = None
def fit_normal(self, X_train, y_train):
assert X_train.shape[0] == y_train.shape[0], \
"the size of X_train must be equal to the size of y_train"
X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
self._theta = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y_train)
self.intercept_ = self._theta[0]
self.coef_ = self._theta[1:]
return self
def fit_gd(self, X_train, y_train, eta=0.01, n_iters=1e4):
assert X_train.shape[0] == y_train.shape[0], \
"the size of X_train must be equal to the size of y_train"
def J(theta, X_b, y):
try:
return np.sum((y - X_b.dot(theta)) ** 2) / len(y)
except:
return float('inf')
def dJ(theta, X_b, y):
return X_b.T.dot(X_b.dot(theta) - y) * 2. / len(y)
def gradient_descent(X_b, y, initial_theta, eta, n_iters=1e4, epsilon=1e-8):
theta = initial_theta
cur_iter = 0
while cur_iter < n_iters:
gradient = dJ(theta, X_b, y)
last_theta = theta
theta = theta - eta * gradient
if (abs(J(theta, X_b, y) - J(last_theta, X_b, y)) < epsilon):
break
cur_iter += 1
return theta
X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
initial_theta = np.zeros(X_b.shape[1])
self._theta = gradient_descent(X_b, y_train, initial_theta, eta, n_iters)
self.intercept_ = self._theta[0]
self.coef_ = self._theta[1:]
return self
def predict(self, X_predict):
assert self.intercept_ is not None and self.coef_ is not None, \
"must fit before predict!"
assert X_predict.shape[1] == len(self.coef_), \
"the feature number of X_predict must be equal to X_train"
X_b = np.hstack([np.ones((len(X_predict), 1)), X_predict])
return X_b.dot(self._theta)
def score(self, X_test, y_test):
y_predict = self.predict(X_test)
return r2_score(y_test, y_predict)
def __repr__(self):
return "LR()"
train_X = train_X.astype('float')
train_y = train_y.astype('float')
def r2_score(y_true, y_predict):
MSE = np.sum((y_true - y_predict) ** 2) / len(y_true)
return 1 - MSE / np.var(y_true)
LR = LinearRegression().fit_gd(train_X, train_y)
LR.score(train_X, train_y)
result = LR.predict(test_x)
