李宏毅作业1 PM2.5预测

Homework 1: Linear Regression

本次目標：由前 9 個小時的 18 個 features (包含 PM2.5)預測的 10 個小時的 PM2.5。

Load ‘train.csv’

train.csv 的資料為 12 個月中，每個月取 20 天，每天 24 小時的資料(每小時資料有 18 個 features)。

import sys
import pandas as pd
import numpy as np

data = pd.read_csv('./data/train.csv', encoding = 'big5')

data.iloc[0:4,0:12] #看看数据

	日期	測站	測項	0	1	2	3	4	5	6	7	8
0	2014/1/1	豐原	AMB_TEMP	14	14	14	13	12	12	12	12	15
1	2014/1/1	豐原	CH4	1.8	1.8	1.8	1.8	1.8	1.8	1.8	1.8	1.8
2	2014/1/1	豐原	CO	0.51	0.41	0.39	0.37	0.35	0.3	0.37	0.47	0.78
3	2014/1/1	豐原	NMHC	0.2	0.15	0.13	0.12	0.11	0.06	0.1	0.13	0.26

Preprocessing

將 ‘RAINFALL’ 欄位全部補 0。数据RAINFALL中有一些值是NR，把他变成0

data=data.iloc[:,3:] #去掉数据的前面不包含数据的三列

data.head(3)

	0	1	2	3	4	5	6	7	8	9	...	14	15	16	17	18	19	20	21	22	23
0	14	14	14	13	12	12	12	12	15	17	...	22	22	21	19	17	16	15	15	15	15
1	1.8	1.8	1.8	1.8	1.8	1.8	1.8	1.8	1.8	1.8	...	1.8	1.8	1.8	1.8	1.8	1.8	1.8	1.8	1.8	1.8
2	0.51	0.41	0.39	0.37	0.35	0.3	0.37	0.47	0.78	0.74	...	0.37	0.37	0.47	0.69	0.56	0.45	0.38	0.35	0.36	0.32

3 rows × 24 columns

data[data=='NR']=0 #把值为NR的置为0
raw_data=data.to_numpy()

raw_data[0:3,0:6]

array([['14', '14', '14', '13', '12', '12'],
       ['1.8', '1.8', '1.8', '1.8', '1.8', '1.8'],
       ['0.51', '0.41', '0.39', '0.37', '0.35', '0.3']], dtype=object)

提取特征

原始的数据是(4320, 24),现在把它分为12个月。每个月18(features) * 480 (hours) .只有480=20天*24个小时

raw_data.shape

(4320, 24)

month_data={} #声明字典
for month in range(12):
    sample=np.empty([18,480]) #暂存一个月的数据
    for day in range(20): #一个月只要20天
        sample[:,day*24:(day+1)*24]=raw_data[18 * (20 * month + day) : 18 * (20 * month + day + 1), :]
    month_data[month]=sample

每個月會有 480hrs，每 9 小時形成一個 data，每個月會有 471 個 data，故總資料數為 471 * 12 筆，而每筆 data 有 9 * 18 的 features (一小時 18 個 features * 9 小時)。

對應的 target 則有 471 * 12 個(第 10 個小時的 PM2.5)

x = np.empty([12 * 471, 18 * 9], dtype = float)
y = np.empty([12 * 471, 1], dtype = float)
for month in range(12):
    for day in range(20):
        for hour in range(24):
            if day == 19 and hour > 14:
                continue
            x[month * 471 + day * 24 + hour, :] = month_data[month][:,day * 24 + hour : day * 24 + hour + 9].reshape(1, -1) #vector dim:18*9 (9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9)
            y[month * 471 + day * 24 + hour, 0] = month_data[month][9, day * 24 + hour + 9] #value
print(x.shape)
print(y.shape)

(5652, 162)
(5652, 1)

mean_x = np.mean(x, axis = 0) #18 * 9 
std_x = np.std(x, axis = 0) #18 * 9 
for i in range(len(x)): #12 * 471
    for j in range(len(x[0])): #18 * 9 
        if std_x[j] != 0:
            x[i][j] = (x[i][j] - mean_x[j]) / std_x[j]

Normalize

np.set_printoptions(suppress=True) # 使输出不是科学计数法的形式，而是浮点型的数据输出。
print(mean_x[0:5])
np.mean(x,axis=0)[0:5]

[22.53651805 22.54324133 22.54978769 22.55368011 22.55598018]





array([ 0., -0.,  0., -0., -0.])

划分数据集 “train_set” and “validation_set”

這部分是針對作業中 report 的第二題、第三題做的簡單示範，以生成比較中用來訓練的 train_set 和不會被放入訓練、只是用來驗證的 validation_set。

import math
x_train_set=x[:math.floor(len(x)*0.8),:] #train set是前80%的数据
y_train_set=y[:math.floor(len(x)*0.8),:]
x_validation=x[math.floor(len(x)*0.8):,:]
y_validation=y[math.floor(len(x)*0.8):,:]

训练

dim=18*9+1
w=np.zeros([dim,1])
x=np.concatenate( (np.ones([12*471,1]) , x), axis=1).astype(float) #增加b
learning_rate=100
iter_time=1000
adagrad=np.zeros([dim,1])
eps = 0.0000000001
loss_all=[]
for i in range(iter_time):
    loss= np.sqrt(np.sum((np.dot(x,w)-y)**2)/471/12)
    if(i%50==0):
        loss_all.append(loss)
    gradient=2*np.dot(x.T, np.dot(x,w)-y)
    adagrad+=gradient**2
    w=w-learning_rate*gradient/np.sqrt(adagrad+eps)
np.save('weight.npy', w)

import matplotlib.pyplot as plt
plt.plot(loss_all)

np.mean(np.dot(x,w)-y)

-3.3389977108613344e-15

dim=18*9+1
w=np.zeros([dim,1])
x=np.concatenate( (np.ones([len(x_train_set),1]) , x_train_set), axis=1).astype(float) #增加b
learning_rate=100
iter_time=1000
adagrad=np.zeros([dim,1])
eps = 0.0000000001
loss_all=[]
y=y_train_set
for i in range(iter_time):
    loss= np.sqrt(np.sum((np.dot(x,w)-y)**2)/471/12)
    if(i%50==0):
        loss_all.append(loss)
    gradient=2*np.dot(x.T, np.dot(x,w)-y)
    adagrad+=gradient**2
    w=w-learning_rate*gradient/np.sqrt(adagrad+eps)
np.save('weight.npy', w)

np.mean(np.dot(x,w)-y)
np.sqrt(np.sum((np.dot(x,w)-y)**2)/471/12)

5.342550298780461

x_validation=np.concatenate( (np.ones([len(x_validation),1]) , x_validation), axis=1).astype(float) #增加b
np.mean(np.dot(x_validation,w)-y_validation)
np.sqrt(np.sum((np.dot(x_validation,w)-y_validation)**2)/471/12)

2.6447202742907496