一、說明

給定訓(xùn)練集train.csv，要求根據(jù)前9個小時的空氣監(jiān)測情況預(yù)測第10個小時的PM2.5含量。

訓(xùn)練集介紹：

(1)、CSV文件，包含臺灣豐原地區(qū)240天的氣象觀測資料(取每個月前20天的數(shù)據(jù)做訓(xùn)練集，12月X20天=240天，每月后10天數(shù)據(jù)用于測試，對學(xué)生不可見);

(2)、每天的監(jiān)測時間點為0時，1時......到23時，共24個時間節(jié)點;

(3)、每天的檢測指標(biāo)包括CO、NO、PM2.5、PM10等氣體濃度，是否降雨、刮風(fēng)等氣象信息，共計18項；

二、數(shù)據(jù)處理

根據(jù)要求，要用前9個小時的數(shù)據(jù)，來預(yù)測第10個小時的PM2.5值。一筆訓(xùn)練數(shù)據(jù)如下圖所示：

數(shù)據(jù)中存在一定量的空數(shù)據(jù)NR，且多存在于RAINFALL一項。對于空數(shù)據(jù)，常規(guī)的處理方法無非就是刪除法和補(bǔ)全法兩種。

RAINFALL表示當(dāng)天對應(yīng)時間點是否降雨，因此可以采用補(bǔ)全法處理空數(shù)據(jù)：將空數(shù)據(jù)NR全部補(bǔ)為0即可。

#將NR替換成 0

data = data.replace(['NR'], [0.0])

我們先將數(shù)據(jù)進(jìn)行預(yù)處理，得到 每筆訓(xùn)練數(shù)據(jù) 和對應(yīng)的 結(jié)果label。

1. 由于每個月前20天的數(shù)據(jù)檢測是連續(xù)24小時進(jìn)行的，所以為了得到多筆數(shù)據(jù)，先將每個月20天數(shù)據(jù) “連起來”，如下圖示：

每個月的數(shù)據(jù)就是18行480(24*20)列,一共12個月的數(shù)據(jù)。

#將每個月20天數(shù)據(jù)連成一大行

month_data =[]for month in range(12):#每個月的數(shù)據(jù)

sub_data = np.empty([18, 20*24])for day in range(20):#每一天的數(shù)據(jù)

sub_data[:, day*24:(day+1)*24] = data[(month*18*20+day*18):(month*18*20+(day+1)*18), :]

month_data.append(sub_data)

2. 對于連續(xù)的10個小時，可以取一筆 9小時訓(xùn)練數(shù)據(jù)(輸入) 和 第10小時對應(yīng)PM2.5值(結(jié)果)。

每個月20天，20*24=480小時， 480-9=471，每個月可以取471筆數(shù)據(jù)。

#將每個月中20天，相鄰9個小時生成一筆數(shù)據(jù)，第10個小時的pm2.5值，生成一個label

for i in range(12):

sub_data=month_data[i]for j in range(20*24-9):#相鄰9小時的數(shù)據(jù)

x_list.append(sub_data[:, j:j+9])#第10小時的 pm2.5

y_list.append(sub_data[9, j+9])

完整數(shù)據(jù)處理代碼：

defdata_process(data):

x_list , y_list=[], []#將NR替換成 0

data = data.replace(['NR'], [0.0])#astype() 轉(zhuǎn)換為float

data =np.array(data).astype(float)#將每個月20天數(shù)據(jù)連成一大行

month_data =[]for month in range(12):#每個月的數(shù)據(jù)

sub_data = np.empty([18, 20*24])for day in range(20):#每一天的數(shù)據(jù)

sub_data[:, day*24:(day+1)*24] = data[(month*18*20+day*18):(month*18*20+(day+1)*18), :]

month_data.append(sub_data)#將每個月中20天，相鄰9個小時生成一筆數(shù)據(jù)，第10個小時的pm2.5值，生成一個label

for i in range(12):

sub_data=month_data[i]for j in range(20*24-9):#相鄰9小時的數(shù)據(jù)

x_list.append(sub_data[:, j:j+9])#第10小時的 pm2.5

y_list.append(sub_data[9, j+9])

x=np.array(x_list)

y=np.array(y_list)return x, y, month_data

View Code

二、模型建立

如題所說，使用的是最簡單的線性回歸模型，作為課程作業(yè)沒有那么難,但也從中學(xué)到不少東西。

2.1 線性回歸模型

如果把b作為w0，加到權(quán)值向量前面，可以得到向量運算的形式，如下：

h(X) = WTX? ? ? ? #W為權(quán)值， X為輸入。

2.2?損失函數(shù)

用預(yù)測值與label之間的平均歐式距離來衡量預(yù)測的準(zhǔn)確程度，并充當(dāng)損失函數(shù)。

這里的損失指的是平均損失；乘1/2是為了在后續(xù)求梯度過程中保證梯度項系數(shù)為1，方便計算。

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

完整的損失函數(shù)：

2.3 梯度下降

對參數(shù) w 和 b 求偏導(dǎo)：

對參數(shù)進(jìn)行更新：

ηw、ηb為學(xué)習(xí)率。

2.4 學(xué)習(xí)率更新

為了在不影響模型效果的前提下提高學(xué)習(xí)速度，可以對學(xué)習(xí)率進(jìn)行實時更新：即讓學(xué)習(xí)率的值在學(xué)習(xí)初期較大，之后逐漸減小。

這里采用比較經(jīng)典的adagrad算法來更新學(xué)習(xí)率：

? ? ?根號下為梯度的累加值。

2.5 矩陣加速計算推導(dǎo)

因為python中使用矩陣的計算速度非常快，遠(yuǎn)遠(yuǎn)快于循環(huán)計算，所以這里我們推導(dǎo)一下利用矩陣計算梯度值的寫法。

前文提到過，如果把b作為w0，加到權(quán)值向量前面可以得到如下形式：

h(X) = WTX? ? ? ? #??W為權(quán)值， X為輸入。? W = [b, w0, w1, ...]

對于損失函數(shù)的轉(zhuǎn)化改造：

對于

1. 將平方求和改造成向量模的平方：

假設(shè)我們有個向量?

所以有

將平方項展開：? ? ? ??

對W進(jìn)行求偏導(dǎo)：?與上面展開對應(yīng)

所以我們要求的梯度就是：? ?

*這里將分子的2去掉是因為與的損失函數(shù)分母抵消。

#計算梯度 W = X轉(zhuǎn)置.(XW-Y)

w_1 = np.dot(X.transpose(), X.dot(W)-y_train)

三、訓(xùn)練模型

3.0 數(shù)據(jù)轉(zhuǎn)化

將訓(xùn)練數(shù)據(jù)分成兩部分(8：2)，一部分用來訓(xùn)練，一部分用來驗證效果。

#8:2 cross validation

x_train = x[:(int)(x.shape[0]*0.8)]

y_train= y[:(int)(x.shape[0]*0.8)]

x_val= x[(int)(x.shape[0]*0.8+0.5):]

y_val= y[(int)(y.shape[0]*0.8+0.5):]

由于參數(shù)太多，也可以取其中的幾類進(jìn)行訓(xùn)練，比如下文中將選取NO、NO2、NOx、O3、PM10、PM2.5作為輸入。

其中的一筆數(shù)據(jù)如下：

首先將每筆數(shù)據(jù)的輸入轉(zhuǎn)化成一行，并在前面加上 1， 對應(yīng)于bias項。

#定義參數(shù) b，w b作為w0

W = np.ones(1+9*6)#將訓(xùn)練數(shù)據(jù)轉(zhuǎn)化成 每一筆數(shù)據(jù)一行，并且前面添加 1，作為b的權(quán)值 [[1, ...], [1, ...],...,[1, ...]]

X = np.empty([n, W.size-1])for i inrange(n):

X[i]= x_train[i][4:10].reshape(1, -1)#添加 1

X = np.concatenate((np.ones([n, 1]), X), axis=1)

3.1 訓(xùn)練函數(shù)

完整的訓(xùn)練函數(shù)代碼如下，具體請看注釋：

deftrain(x_train, y_train, times):#定義參數(shù) b，w b作為w0

W = np.ones(1+9*6)#多少筆數(shù)據(jù)

n =y_train.size#學(xué)習(xí)率

learning_rate = 100

#正則項大小

reg_rate = 0.011

#將訓(xùn)練數(shù)據(jù)轉(zhuǎn)化成 每一筆數(shù)據(jù)一行，并且前面添加 1，作為b的權(quán)值 [[1, ...], [1, ...],...,[1, ...]]

X = np.empty([n, W.size-1])for i inrange(n):

X[i]= x_train[i][4:10].reshape(1, -1)#添加 1

X = np.concatenate((np.ones([n, 1]), X), axis=1)

# 累加正則項

adagrad=0#正則項的選擇矩陣， 去掉bias部分

reg_mat=np.concatenate((np.array([0]), np.ones([9*6,])), axis=0)for t inrange(times):#計算梯度 W = X轉(zhuǎn)置.(XW-Y)

w_1 = np.dot(X.transpose(), X.dot(W)-y_train)#加正則項

w_1 += reg_rate * W *reg_mat#正則項參數(shù)更新

adagrad += sum(w_1**2)**0.5

#梯度下降

W -= learning_rate/adagrad *w_1#每200次迭代輸出一次

if t%200==0:

loss=0for j inrange(n):

loss+= (y_train[j]-X[j].dot(W))**2

print(t)print('times', loss/n)return W

3.2 驗證

defvalidate(x_val, y_val, w):

n=y_val.size#轉(zhuǎn)化成一行，并加一列 1

X = np.empty([n, w.size - 1])for i inrange(n):

X[i]= x_val[i][4:10].reshape(1, -1)

X= np.concatenate((np.ones([n, 1]), X), axis=1)

loss=0#計算loss

for j inrange(n):

loss+= (y_val[j] - X[j].dot(W)) ** 2

return loss/n

四、結(jié)果分析

運行輸出的結(jié)果看，loss還是挺大的，還有改進(jìn)的空間。

改進(jìn)思路：

1. 分割訓(xùn)練集和驗證集時，應(yīng)該按照比例隨機(jī)抽取數(shù)據(jù)幀作為訓(xùn)練集和驗證集，選取loss最小的模型。

2. 充分考慮其他參數(shù)對空氣PM2.5的影響，加入更加復(fù)雜的高次項。

五、預(yù)測結(jié)果

對test集的數(shù)據(jù)進(jìn)行結(jié)果預(yù)測

## 計算預(yù)測值 ##

Y =X_test.dot(W)#預(yù)測值寫入

data_test =np.array(data_test)

data_test= np.concatenate((data_test, np.zeros([n, 1])), axis=1)for j in range(0, n, 18):

data_test[j+9, 11] = int(Y[int(j/18)]+0.5)

為了方便查看，將數(shù)據(jù)寫回源文件格式。

#保存結(jié)果

data_test =pd.DataFrame(data_test)

data_test.to_csv('test_res.csv')

第一筆數(shù)據(jù)的預(yù)測值：

六、程序代碼

**在項目根目錄存放‘train.csv’、'test.csv'

**每次訓(xùn)練后會保存參數(shù)，下次訓(xùn)練時請事先刪除根目錄文件‘weight_2.npy’

importpandas as pdimportnumpy as npdefdata_process(data):

x_list , y_list=[], []#將NR替換成 0

data = data.replace(['NR'], [0.0])#astype() 轉(zhuǎn)換為float

data =np.array(data).astype(float)#將每個月20天數(shù)據(jù)連成一大行

month_data =[]for month in range(12):#每個月的數(shù)據(jù)

sub_data = np.empty([18, 20*24])for day in range(20):#每一天的數(shù)據(jù)

sub_data[:, day*24:(day+1)*24] = data[(month*18*20+day*18):(month*18*20+(day+1)*18), :]

month_data.append(sub_data)#將每個月中20天，相鄰9個小時生成一筆數(shù)據(jù)，第10個小時的pm2.5值，生成一個label

for i in range(12):

sub_data=month_data[i]for j in range(20*24-9):#相鄰9小時的數(shù)據(jù)

x_list.append(sub_data[:, j:j+9])#第10小時的 pm2.5

y_list.append(sub_data[9, j+9])

x=np.array(x_list)

y=np.array(y_list)returnx, y, month_datadeftrain(x_train, y_train, times):#定義參數(shù) b，w b作為w0

W = np.ones(1+9*6)#多少筆數(shù)據(jù)

n =y_train.size#學(xué)習(xí)率

learning_rate = 100

#正則項大小

reg_rate = 0.011

#將訓(xùn)練數(shù)據(jù)轉(zhuǎn)化成 每一筆數(shù)據(jù)一行，并且前面添加 1，作為b的權(quán)值 [[1, ...], [1, ...],...,[1, ...]]

X = np.empty([n, W.size-1])for i inrange(n):

X[i]= x_train[i][4:10].reshape(1, -1)#添加 1

X = np.concatenate((np.ones([n, 1]), X), axis=1)#data_X = pd.DataFrame(X)

#data_X.to_csv('data.csv')

adagrad=0#正則項的選擇矩陣， 去掉bias部分

reg_mat=np.concatenate((np.array([0]), np.ones([9*6,])), axis=0)for t inrange(times):#計算梯度 W = X轉(zhuǎn)置.(XW-Y)

w_1 = np.dot(X.transpose(), X.dot(W)-y_train)#加正則項

w_1 += reg_rate * W *reg_mat#正則項參數(shù)更新

adagrad += sum(w_1**2)**0.5

#梯度下降

W -= learning_rate/adagrad *w_1#每200次迭代輸出一次

if t%200==0:

loss=0for j inrange(n):

loss+= (y_train[j]-X[j].dot(W))**2

print('After', t,'times loss=', loss/n)returnWdefvalidate(x_val, y_val, w):

n=y_val.size#轉(zhuǎn)化成一行，并加一列 1

X = np.empty([n, w.size - 1])for i inrange(n):

X[i]= x_val[i][4:10].reshape(1, -1)

X= np.concatenate((np.ones([n, 1]), X), axis=1)

loss=0#計算loss

for j inrange(n):

loss+= (y_val[j] - X[j].dot(W)) ** 2

return loss/nif __name__ == '__main__':

data= pd.read_csv('./train.csv', encoding='big5')#去掉前三列

data = data.iloc[:, 3:]

[x, y, month_data]=data_process(data)#8:2 cross validation

x_train = x[:(int)(x.shape[0]*0.8)]

y_train= y[:(int)(x.shape[0]*0.8)]

x_val= x[(int)(x.shape[0]*0.8+0.5):]

y_val= y[(int)(y.shape[0]*0.8+0.5):]try:

W= np.load('weight_2.npy')except:#迭代次數(shù)

times = 10000W=train(x_train, y_train, times)

np.save('weight_2.npy', W)## 計算在val上的loss ##

loss =validate(x_val, y_val, W)print('validate loss=', loss)## 在test上進(jìn)行驗證 ##

#header=None 無表頭讀入

data_test = pd.read_csv('./test.csv', header=None, encoding='big5')#去掉前兩列

test = data_test.iloc[:, 2:]

test= test.replace(['NR'], [0.0])#處理數(shù)據(jù)

test =np.array(test).astype(float)

[n, m]=test.shape#讀出參數(shù)值

X_test = np.empty([int(n/18), 9*6])for i in range(0, n, 18):

X_test[int(i/18), :] = test[i+4:i+10, :].reshape(1, -1)

[n_test, m_test]=X_test.shape#加一列 1

X_test = np.concatenate((np.ones([n_test, 1]), X_test), axis=1)## 計算預(yù)測值 ##

Y =X_test.dot(W)#預(yù)測值寫入

data_test =np.array(data_test)

data_test= np.concatenate((data_test, np.zeros([n, 1])), axis=1)for j in range(0, n, 18):

data_test[j+9, 11] = int(Y[int(j/18)]+0.5)#保存結(jié)果

data_test =pd.DataFrame(data_test)

data_test.to_csv('test_res.csv')

View Code

感謝閱讀，如有錯誤歡迎留言指正。

ps：本文實現(xiàn)參照以下兩篇博客：

總結(jié)

以上是生活随笔為你收集整理的python线性回归预测pm2.5_线性回归--PM2.5预测--李宏毅机器学习的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。