李宏毅机器学习作业1-使用线性回归预测PM2.5

题目

  给定台湾丰原地区一年（12个月）的PM2.5及各项污染指标的数据，抽取每个月前二十天数据组成训练集train.csv，取每月后十天的数据组成测试集test.csv，其中每10小时的数据为一组，前九小时作为输入，要求根据输入预测第十小时的输出。
  由于编码不同，我怒改了一波表头并改了文件编码，这不影响结果。

图1 train.csv中的数据

正文

数据预处理
  观察train.csv文件发现，对应RAINFALL一项有很多缺失值，用NR表示。分析得出，RAINFALL应该是一个$[01]$之间的浮点数，表示下雨还是不下，那么NR很有可能代表的是0，这里先将NR都表示为0.0.

数据存储
  由于题目要求是根据前九个小时预测第十个小时的PM2.5值，因此在一天的24小时中，都可以以$9+1$为滑动窗口采集数据，则一天可采集15组数据，总共可采集$12\times 20\times (24-9)=3600$组数据；
  数据集是否可以扩大呢？想象一下，前面说到每个月采取前20天数据作为已知，而这前20天又都是连续的，因此可以以20天为1组，将每20天的小时数据首尾相连，这样通过滑动窗口获得的数据数量为$12\times (20\times 24-9)=5652$组。

回归模型
  采用多元线性回归模型，其目标为：$\arg {\min }_{w,b}\frac{1}{\text{ n }}{\sum }_{i=1}^n{\left({\hat{y}}^i-{\mathbf{x}}^{iT}\cdot \mathbf{w}-b\right)}^2$
其中$\mathbf{x}$，$\mathbf{w}$均为向量，且均属于$R^{9\times 1}$.

梯度下降法求解
  经调参发现，学习率为0.00001，迭代次数为1000次时结果尚佳，训练集平均损失为：4.473245434194629，测试集平均损失为：4.97302969793434。
  另外，本模型采用adagrad进行了优化。

线性代数方法求标准解（基于上面的模型，只考虑PM2.5）
  计算得出训练集平均损失为：4.443108788904938

代码

# -*- coding:utf-8 -*-
import numpy as np
import pandas as pd


def data_processing(df):
    df = df.replace('NR', 0.0)
    temp = np.array(df).astype(float)

    # 每隔18取一次PM2.5，分成12月（组），每组480小时
    temp = temp[9:len(temp):18]
    new = []
    for i in range(0, len(temp), 20):
        t = []
        for j in range(i, i+20):
            t.extend(temp[j].tolist())
        new.append(t)

    # 下面取训练集
    x_data, y_data = [], []
    for i in range(len(new)):
        for j in range(9, len(new[0])):
            x_data.append(new[i][j-9:j])
            y_data.append(new[i][j])

    return x_data, y_data


def linear_algbra(x_data, y_data):
    y = np.array(y_data).T
    x = np.array(x_data)
    x = np.insert(x, 0, len(x)*[1], axis=1)
    # print(x)
    temp = np.matrix(x.T.dot(x))
    weight = temp.I.dot(x.T).dot(y)
    return weight, list(x)


def count_loss(standard, x_data, y_data):
    sum = 0.0
    for i in range(len(y_data)):
        x = np.array(x_data[i])
        st = np.array(standard)
        sum += np.abs(y_data[i] - x.dot(st.T))
    sum /= len(y_data)
    print('average loss: ', float(sum))


def gradient_descent(x_data, y_data):
    # 初始化
    lr = 1
    epoch = 1000
    weight = np.zeros(9)
    bias = 0
    b_sum, w_sum = 0.0, 0.0
    # 梯度下降
    for e in range(epoch):
        tw = np.zeros(9)
        tb = 0.0
        for i in range(len(y_data)):
            tb -= 2.0*(y_data[i]-bias-weight.T.dot(x_data[i]))
            tw -= 2.0*(y_data[i]-bias-weight.T.dot(x_data[i]))*np.array(x_data[i])
        # 下面不加的话会上溢
        tw /= len(y_data)
        tb /= len(y_data)
        # anagrad
        b_sum += tb**2
        w_sum += tw.dot(tw.T)
        bias -= lr/np.sqrt(b_sum)*tb
        weight -= lr/np.sqrt(w_sum)*tw

    return np.insert(weight, 0, [bias])


def test(standard):
    ans = pd.read_csv('ans.csv', usecols=[1])
    y_test = ans.astype(float)
    df = pd.read_csv('test.csv', usecols=[1])
    res = []
    df = df.replace('NR', 0.0)
    x_test = df.astype(float)
    x_test = x_test[8::18]
    t = np.copy(x_test)
    x_test = np.insert(t, 0, len(t)*[1], axis=1)
    sum = 0.0
    for i in range(len(x_test)):
        temp = standard.T.dot(x_test[i])
        sum += abs(np.array(y_test)[i] - temp)
        res.append(temp)
    sum /= len(x_test)
    return res, sum


if __name__ == '__main__':
    df = pd.read_csv('train.csv', usecols=np.arange(3, 27).tolist())  # usecols控制读取的列，这里不读前三列
    x_data, y_data = data_processing(df)
    standard, x1 = linear_algbra(x_data, y_data)  # 计算标准的weight
    gd = gradient_descent(x_data, y_data)         # 使用梯度下降计算出的权重
    print('linear algobra', end=' ')
    count_loss(standard, x1, y_data)              # 4.443108788904938
    print('gradient descent', end=' ')
    count_loss(gd, x1, y_data)                    # 5.057686349538455
    res, los = test(gd)
    print('loss in test.csv: ', los)

结果