線性迴歸預測PM2.5----臺大李宏毅機器學習做業1(HW1)

1、做業說明

　　給定訓練集train.csv，要求根據前9個小時的空氣監測狀況預測第10個小時的PM2.5含量。

訓練集介紹：

　　(1)、CSV文件，包含臺灣豐原地區240天的氣象觀測資料(取每月前20天的數據作訓練集，12月X20天=240天，每個月後10天數據用於測試，對學生不可見);

　　(2)、天天的監測時間點爲0時，1時......到23時，共24個時間節點;

　　(3)、天天的檢測指標包括CO、NO、PM2.五、PM10等氣體濃度，是否降雨、颳風等氣象信息，共計18項；

　　(4)、數據集地址：https://pan.baidu.com/s/1o2Yx42dZBJZFZqCa5y3WzQ，提取碼：qgtm。 

 用excel打開，繁體字會出現亂碼：

用notepad++能夠正常打開：

2、思路分析及代碼實現

　　前注：下文中提到的「數據幀」並不是指pandas庫中的數據結構DataFrame,而是指一個二維的數據包。

2.1 數據預處理

　　訓練集中數據排列形式符合人類觀察數據的習慣，但並不能直接拿來餵給模型進行訓練，所以須要對數據進行預處理。

　　瀏覽數據可知，數據中存在必定量的空數據NR，且多存在於RAINFALL一項。對於空數據，常規的處理方法無非就是刪除法和補全法兩種。查閱資料後發現，RAINFALL表示當天對應時間點是否降雨，有降雨值爲1，無降雨值爲NR，相似於布爾變量。所以能夠採用補全法處理空數據：將空數據NR所有補爲0便可。

　　根據做業要求可知，須要用到連續9個時間點的氣象觀測數據，來預測第10個時間點的PM2.5含量。針對每一天來講，其包含的信息維度爲(18,24)(18項指標，24個時間節點)。能夠將0到8時的數據截取出來，造成一個維度爲(18,9)的數據幀，做爲訓練數據，將9時的PM2.5含量取出來，做爲該訓練數據對應的label；同理可取1到9時的數據做爲訓練用的數據幀，10時的PM2.5含量做爲label......以此分割，可將天天的信息分割爲15個shape爲(18,9)的數據幀和與之對應的15個label。

　　訓練集中共包含240天的數據，所以共可得到240X15=3600個數據幀和與之對應的3600個label。

# 數據預處理
def dataProcess(df):
    x_list, y_list = [], []
    # df替換指定元素，將空數據填充爲0
    df = df.replace(['NR'], [0.0])
    # astype() 轉換array中元素數據類型
    array = np.array(df).astype(float)
    # 將數據集拆分爲多個數據幀
    for i in range(0, 4320, 18):
        for j in range(24-9):
            mat = array[i:i+18, j:j+9]
            label = array[i+9, j+9] # 第10行是PM2.5
            x_list.append(mat)
            y_list.append(label)
    x = np.array(x_list)
    y = np.array(y_list)
    
    return x, y, array

2.2 模型創建

　　若是對相關領域比較熟悉的話，能夠根據PM2.5與PM十、SO、NO的濃度關係選擇合適的模型。

　　若是對數據比較敏感的話，能夠從數據中發現規律並以此爲依據創建模型。

　　不過筆者對氣象領域並不熟悉，對數據也不夠敏感，只能採用最簡單、最low的線性迴歸模型。不過既然是做業嘛，就應該容許學生隨意發揮，不見得就存在標準答案。

2.2.1 迴歸模型

　　採用最普通的線性迴歸模型，並無用上訓練集中全部的數據，只用到了每一個數據幀樣本中的9個PM2.5含量值：

爲對應數據幀中第i個PM2.5含量，爲其對應的權重值，爲偏置項，爲該數據幀樣本的預測結果。                                   

2.2.2 損失函數

　　用預測值與label之間的平均歐式距離來衡量預測的準確程度，並充當損失函數（這裏的損失指的是平均損失；乘1/2是爲了在後續求梯度過程當中保證梯度項係數爲1，方便計算）：

 

爲第n個label，爲第n個數據幀的預測結果，爲參加訓練的數據幀樣本個數。

　　爲了防止過擬合，加入正則項：

 

爲正則項，爲正則項係數。

2.2.3 梯度更新

　　梯度計算：需明確此時的目標是使Loss最小，而可優化的參數爲權重w和偏置值b，所以須要求Loss在w上的偏微分和Loss在b上的偏微分。

　　計算出梯度後，經過梯度降低法實現參數更新。

爲權重w更新時的學習率，爲偏置b更新時的學習率。

2.2.3 學習率更新

　　爲了在不影響模型效果的前提下提升學習速度，能夠對學習率進行實時更新：即讓學習率的值在學習初期較大，以後逐漸減少。這裏採用比較經典的adagrad算法來更新學習率。

爲更新後的學習率，爲更新前的學習率。爲在此以前全部梯度平方和的二次根。

# 更新參數，訓練模型
def train(x_train, y_train, epoch):
    bias = 0 # 偏置值初始化
    weights = np.ones(9) # 權重初始化
    learning_rate = 1 # 初始學習率
    reg_rate = 0.001 # 正則項係數
    bg2_sum = 0 # 用於存放偏置值的梯度平方和
    wg2_sum = np.zeros(9) # 用於存放權重的梯度平方和

    for i in range(epoch):
        b_g = 0
        w_g = np.zeros(9)
        # 在全部數據上計算Loss_label的梯度
        for j in range(3200):
            b_g += (y_train[j] - weights.dot(x_train[j, 9, :]) - bias) * (-1)
            for k in range(9):
                w_g[k] += (y_train[j] - weights.dot(x_train[j, 9, :]) - bias) * (-x_train[j, 9, k]) 
        # 求平均    
        b_g /= 3200
        w_g /= 3200
        #  加上Loss_regularization在w上的梯度
        for m in range(9):
            w_g[m] += reg_rate * weights[m]
        
        # adagrad
        bg2_sum += b_g**2
        wg2_sum += w_g**2
        # 更新權重和偏置
        bias -= learning_rate/bg2_sum**0.5 * b_g
        weights -= learning_rate/wg2_sum**0.5 * w_g

    return weights, bias

 3、代碼分享與結果分析

3.1 源代碼

 1 import pandas as pd
 2 import numpy as np
 3 
 4 # 數據預處理
 5 def dataProcess(df):
 6     x_list, y_list = [], []
 7     # df替換指定元素，將空數據填充爲0
 8     df = df.replace(['NR'], [0.0])
 9     # astype() 轉換array中元素數據類型
10     array = np.array(df).astype(float)
11     # 將數據集拆分爲多個數據幀
12     for i in range(0, 4320, 18):
13         for j in range(24-9):
14             mat = array[i:i+18, j:j+9]
15             label = array[i+9, j+9] # 第10行是PM2.5
16             x_list.append(mat)
17             y_list.append(label)
18     x = np.array(x_list)
19     y = np.array(y_list)
20 
21     '''
22     # 將每行數據都scale到0到1的範圍內，有利於梯度降低，但經嘗試發現效果並很差
23     for i in range(18):
24         if(np.max(x[:, i, :]) != 0):
25             x[: , i, :] /= np.max(x[:, i, :])
26     '''
27     return x, y, array
28 
29 # 更新參數，訓練模型
30 def train(x_train, y_train, epoch):
31     bias = 0 # 偏置值初始化
32     weights = np.ones(9) # 權重初始化
33     learning_rate = 1 # 初始學習率
34     reg_rate = 0.001 # 正則項係數
35     bg2_sum = 0 # 用於存放偏置值的梯度平方和
36     wg2_sum = np.zeros(9) # 用於存放權重的梯度平方和
37 
38     for i in range(epoch):
39         b_g = 0
40         w_g = np.zeros(9)
41         # 在全部數據上計算Loss_label的梯度
42         for j in range(3200):
43             b_g += (y_train[j] - weights.dot(x_train[j, 9, :]) - bias) * (-1)
44             for k in range(9):
45                 w_g[k] += (y_train[j] - weights.dot(x_train[j, 9, :]) - bias) * (-x_train[j, 9, k]) 
46         # 求平均    
47         b_g /= 3200
48         w_g /= 3200
49         #  加上Loss_regularization在w上的梯度
50         for m in range(9):
51             w_g[m] += reg_rate * weights[m]
52         
53         # adagrad
54         bg2_sum += b_g**2
55         wg2_sum += w_g**2
56         # 更新權重和偏置
57         bias -= learning_rate/bg2_sum**0.5 * b_g
58         weights -= learning_rate/wg2_sum**0.5 * w_g
59 
60         # 每訓練200輪，輸出一次在訓練集上的損失
61         if i%200 == 0:
62             loss = 0
63             for j in range(3200):
64                 loss += (y_train[j] - weights.dot(x_train[j, 9, :]) - bias)**2
65             print('after {} epochs, the loss on train data is:'.format(i), loss/3200)
66 
67     return weights, bias
68 
69 # 驗證模型效果
70 def validate(x_val, y_val, weights, bias):
71     loss = 0
72     for i in range(400):
73         loss += (y_val[i] - weights.dot(x_val[i, 9, :]) - bias)**2
74     return loss / 400
75 
76 def main():
77     # 從csv中讀取有用的信息
78     # 因爲你們獲取數據集的渠道不一樣，因此數據集的編碼格式可能不一樣
79     # 若讀取失敗，可在參數欄中加入encoding = 'gb18030'
80     df = pd.read_csv('train.csv', usecols=range(3,27))
81     x, y, _ = dataProcess(df)
82     #劃分訓練集與驗證集
83     x_train, y_train = x[0:3200], y[0:3200]
84     x_val, y_val = x[3200:3600], y[3200:3600]
85     epoch = 2000 # 訓練輪數
86     # 開始訓練
87     w, b = train(x_train, y_train, epoch)
88     # 在驗證集上看效果
89     loss = validate(x_val, y_val, w, b)
90     print('The loss on val data is:', loss)
91 
92 if __name__ == '__main__':
93     main()

View Code

3.1 結果展現

　　能夠看出，模型在驗證集上的損失爲40左右，即預測值與label之間的平均差別在6到7之間，因而可知，模型的總體效果仍是比較差的。

3.3 模型改進的方向

　　(1)在從csv文件中提取數據幀和label時，本文以天爲單位，天天分割出15個數據幀和15個label。事實上，時間是連續的，能夠將每個月的20天首尾鏈接，再從其中分割數據幀和label，可以使數據幀樣本數量大大提高，可能會使模型效果更優。

　　(2)在構建模型時，應充分考慮PM2.5與其餘大氣成分之間的關係，構建更合理的模型。

　　(3)分割訓練集和驗證集時，應該按照比例隨機抽取數據幀做爲訓練集和驗證集，而不是像本文那樣簡單地把前3200個數據樣本做爲訓練集，後400個做爲驗證集。

 

參考資料：

　　李宏毅老師機器學習課程視頻：https://www.bilibili.com/video/av10590361

　　李宏毅老師機器學習課程講義資料：http://speech.ee.ntu.edu.tw/~tlkagk/courses_ML17_2.html

　　邱錫鵬老師《神經網絡與深度學習》： https://nndl.github.io/