本次作業(yè)使用豐原站的觀測記錄，分成 train set 跟 test set，train set是豐原站每個月的前20天所有資料，test set 則是從豐原站剩下的資料中取樣出來。

train.csv: 每個月前20天每個小時的氣象資料(共24小時，每小時有18種測資)。共12個月。

test.csv: 從剩下的資料當(dāng)中取樣出連續(xù)的10小時為一筆，前九小時的所有觀測數(shù)據(jù)當(dāng)作 feature，第十小時的PM2.5當(dāng)作 answer。一共取出 240 筆不重複的 test data，請根據(jù) feauure 預(yù)測這240筆的PM2.5。

數(shù)據(jù)預(yù)處理

import numpy as np
import pandas as pd

# 導(dǎo)入訓(xùn)練數(shù)據(jù)
df_train = pd.read_csv('train.csv')
df_train.head(10)

可以看出在數(shù)據(jù)集合中有很多其他的觀測變量，這里我們的任務(wù)主要是預(yù)測 PM2.5，所以將無關(guān)的記錄剔除，同時前三列的數(shù)據(jù)也是與任務(wù)無關(guān)的。

df_train = df_train[df_train['observation'] == 'PM2.5']
df_train = df_train.iloc[:,3:]
df_train

最終得到的訓(xùn)練數(shù)據(jù)如下：

接著看一看測試數(shù)據(jù)：

# 導(dǎo)入測試數(shù)據(jù)
df_test = pd.read_csv('./test.csv')
df_test.head(10)

我們的測試數(shù)據(jù)是連續(xù)9個小時的PM2.5觀測值，作為預(yù)測任務(wù)的 features，第 10 個小時的值為需要預(yù)測的 label。

所以需要對訓(xùn)練數(shù)據(jù)進一步處理，選擇連續(xù)的9小時數(shù)據(jù)作為feature，第十個小時的觀測數(shù)據(jù)作為預(yù)測，這樣我們需要選擇連續(xù)的10列數(shù)據(jù)，總共 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) 

構(gòu)建模型

對于時序預(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)