最近開始學(xué)習(xí)機(jī)器學(xué)習(xí)，聽同學(xué)介紹，選擇了臺大李宏毅老師的視頻，https://www.bilibili.com/video/BV1JE411g7XF
老師布置了很多作業(yè)，這里記錄下自己寫作業(yè)的過程和一些參考資料
本篇是作業(yè)一缓呛，根據(jù)豐源站歷史的空氣數(shù)據(jù)，預(yù)測下一小時(shí)的pm2.5的值杭隙。
本次作業(yè)說明視頻參見 https://www.bilibili.com/video/BV1gE411F7td
import numpy as np
import pandas as pd
本次作業(yè)要求只用numpy,pandas手寫線性回歸哟绊，所以只導(dǎo)入pandas和numpy兩個(gè)包，
學(xué)習(xí)這兩個(gè)包推薦一個(gè)gitbook痰憎， 利用python進(jìn)行數(shù)據(jù)分析
另外這本書有本實(shí)體版的票髓，是另外一個(gè)人翻譯的，翻譯的非常差铣耘，非常不建議購買洽沟，看這個(gè)gitbook版的就好。
本次作業(yè)代碼參考 https://zhuanlan.zhihu.com/p/115335196
上面知乎專欄的大佬寫的比較簡單蜗细，我添加了更多的細(xì)節(jié)裆操，便于新人理解。
data = pd.read_csv('data/train.csv', encoding='big5')
首先使用pandas讀取訓(xùn)練文件鳄乏，并做一些預(yù)處理跷车，李宏毅老師是灣灣人棘利，
用的是繁體橱野，所以要指定下編碼為big5，否則里面的中文為亂碼善玫。
看一下文件的結(jié)構(gòu)水援，有27列，第一列為日期茅郎，第三列為檢測項(xiàng)蜗元，后面的是0-23點(diǎn)每一點(diǎn)的數(shù)值。
然后每天有18行系冗，即18個(gè)檢測項(xiàng)奕扣，其中第10行是pm2.5。
一共12月*20天*18項(xiàng)=4320行
data.head(20)
然后前三項(xiàng)是沒用的掌敬，我們用切片去掉惯豆。shape就變?yōu)榱?320*24
數(shù)據(jù)中還有很多是NR，我們替換為0奔害。
接下來就是numpy的處理范圍了楷兽，使用to_numpy轉(zhuǎn)換為矩陣。
data = data.iloc[:,3:]
data.replace('NR',0,inplace=True)
row_data = data.to_numpy()
row_data
參考上圖华临，我們需要把數(shù)據(jù)shape做一個(gè)轉(zhuǎn)換芯杀，把240天內(nèi)相同的測量項(xiàng)放在同一行。
也就是從(12月*20天*18項(xiàng)),24小時(shí)變?yōu)?18項(xiàng),(24小時(shí)*20天),12月
下面的兩個(gè)循環(huán)就是在做轉(zhuǎn)換，每個(gè)月的數(shù)據(jù)寫在以月份為鍵名的字典里的揭厚。
month_data={}
for month in range(12):
sample = np.empty([18, 480])
for day in range(20):
sample[:,day*24:(day+1)*24] = row_data[18*(month*20+day):18*(month*20+(day+1)),:]
month_data[month] = sample
然后我們以10小時(shí)為一個(gè)data却特，把前9個(gè)小時(shí)的data作訓(xùn)練，第10個(gè)小時(shí)的值為target筛圆。
即0-9點(diǎn)的data來預(yù)測10點(diǎn)的pm2.5核偿，第1-10點(diǎn)的和data來預(yù)測11點(diǎn)的pm2.5
最后一筆是471-479來預(yù)測480的pm2.5，所以一共有471筆data
shape從變?yōu)?img class="math-inline" src="https://math.jianshu.com/math?formula=%5Ccolor%7Bred%7D%7B(12%E6%9C%88*471%E7%AC%94)%2C(18%E9%A1%B9*9%E5%B0%8F%E6%97%B6)%7D" alt="\color{red}{(12月*471筆),(18項(xiàng)*9小時(shí))}" mathimg="1">
y就是第10個(gè)小時(shí)的pm2.5的值顽染，shape為
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)
y[month * 471 + day * 24 + hour, 0] = month_data[month][9, day * 24 + hour + 9]
x
接下來做Normalize漾岳，教學(xué)視頻里說明了為什么要做Normalize，不知道的請復(fù)習(xí)視頻
按照上面的公式粉寞，需要求平均值和標(biāo)準(zhǔn)差尼荆，然后x處理對應(yīng)位置的值
mean_x = np.mean(x, axis = 0)
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]
x
接下來就是最重要的部分了，我們開始做梯度下降
我們有9小時(shí)*18項(xiàng)的變量唧垦，所以需要9*18個(gè)w值捅儒，還有一個(gè)bias:b
注意下面代碼的第三行，這里是整個(gè)矩陣前拼接了一列1值振亮，這一列在點(diǎn)乘后就是b
所以這就是為什么dim有18*9+1項(xiàng)
dim = 18 * 9 + 1
w = np.zeros([dim, 1])
x = np.concatenate((np.ones([12 * 471, 1]), x), axis = 1).astype(float)
learning_rate = 100
iter_time = 1000
adagrad = np.zeros([dim, 1])
eps = 0.0000000001
然后定義loss函數(shù)巧还，這里知乎作者用的是均方根誤差來做loss函數(shù)
貼一段百度百科的解釋：
均方根誤差是預(yù)測值與真實(shí)值偏差的平方與觀測次數(shù)n比值的平方根，
在實(shí)際測量中坊秸，觀測次數(shù)n總是有限的麸祷，真值只能用最可信賴（最佳）值來代替。
標(biāo)準(zhǔn)誤差對一組測量中的特大或特小誤差反映非常敏感褒搔，
所以阶牍，標(biāo)準(zhǔn)誤差能夠很好地反映出測量的精密度。
這正是標(biāo)準(zhǔn)誤差在工程測量中廣泛被采用的原因星瘾。
因此走孽，標(biāo)準(zhǔn)差是用來衡量一組數(shù)自身的離散程度，
而均方根誤差是用來衡量觀測值同真值之間的偏差琳状，
它們的研究對象和研究目的不同磕瓷，但是計(jì)算過程類似。
因?yàn)槌?shù)項(xiàng)的存在念逞，所以 dimension (dim) 需要多加一欄困食；
eps 項(xiàng)是避免 adagrad 的分母為 0 而加的極小數(shù)值。
rmse公式：
然后adagrid的公式：
這里對gt的計(jì)算我沒有看懂肮柜，根據(jù)公式應(yīng)該是loss對w的偏導(dǎo)陷舅，以后搞懂了再來補(bǔ)說明吧
adagrid這一行就是上面的累計(jì)求和
然后對w做更新，這樣就可以實(shí)現(xiàn)梯度下降了
最后保存計(jì)算出的w矩陣
for t in range(iter_time):
loss = np.sqrt(np.sum(np.power(np.dot(x, w) - y, 2))/471/12) #rmse
if(t%100==0):
print(str(t) + ":" + str(loss))
gradient = 2 * np.dot(x.transpose(), np.dot(x, w) - y) #dim*1
adagrad += gradient ** 2
w = w - learning_rate * gradient / np.sqrt(adagrad + eps)
np.save('weight.npy', w)
這樣我們的w和b就計(jì)算完了审洞，接下來拿測試集做測試
還是和train.csv一樣先做預(yù)處理和Normalize莱睁，然后用w和測試集做點(diǎn)乘就可以得到預(yù)測值了
記得在拼接1列1作為b
testdata = pd.read_csv('./data/test.csv', header = None, encoding = 'big5')
test_data = testdata.iloc[:, 2:].copy()
test_data.replace('NR', 0, inplace=True)
test_data = test_data.to_numpy()
test_x = np.empty([240, 18*9], dtype = float)
for i in range(240):
test_x[i, :] = test_data[18 * i: 18* (i + 1), :].reshape(1, -1)
for i in range(len(test_x)):
for j in range(len(test_x[0])):
if std_x[j] != 0:
test_x[i][j] = (test_x[i][j] - mean_x[j]) / std_x[j]
test_x = np.concatenate((np.ones([240, 1]), test_x), axis = 1).astype(float)
w = np.load('weight.npy')
ans_y = np.dot(test_x, w)
ans_y[:10]
ans_y就是預(yù)測值了待讳，然后按照作業(yè)要求寫成csv文件。
這里轉(zhuǎn)回pandas處理仰剿，因?yàn)閜andas的to_csv方法很好用
res = pd.DataFrame(ans_y)
res['id'] = ['id_' + str(x) for x in res.index]
res.rename(columns={0:'value'}, inplace=True)
res=res[['id', 'value']] # 調(diào)換順序
res.to_csv('submit.csv', index=False)
最后创淡，可以調(diào)整 learning rate、iter_time (iteration 次數(shù))南吮、
取用 features 的多少(取幾個(gè)小時(shí)琳彩，取哪些特征)
甚至是不同的 model 來超越 baseline。
