HomeWork1 PM2.5 Prediction

課程資料：

從印象筆記移過來修改的，懶得編輯太細了，CSDN又不讓直接插圖，CSDN的編輯器老出出現各種錯誤，太麻煩，反正程式碼能跑。

做題思路：

資料處理

id_x是時間點，每一個具體的x，是一個時間點。

怎麼抽取訓練集？1-9訓練，第十個驗證，2-10訓練，第11個驗證？這樣每單次訓練都要驗證一次？

不對，這個方法是抽訓練集，訓練集包括input和output，前9個id_x是輸入，第10個id_x是output。10個數據加起來，才是一個訓練資料，而不是單條觀測記錄。

第一個簡陋版本（有理解偏差）：

for迴圈，從時刻0到時刻n-9迭代：

    抽0-8作為input，抽9為output

    把抽到的資料想好一個格式存起來。
 

PS：9比1的輸入輸出，也是自己可選的，不是固定的，考慮到給的test集是9:1預測，所以按9:1來設計。

難點：除了天數，還有每天的小時，不是按天預測的，有點複雜了，按小時能不能跨天？可能模型不要這麼複雜？（錯了，要跨）

跨天和不跨天的差距：

模型偏科，就是說從來不預測1點、2點、到9點的資料，只預測10點到24點的資料，這就會產生——用不針對1-9點的模型來預測1-9點，肯定有問題

不跨天的另一個問題就是大幅縮減了資料量——資料集是12個月，每個月20天，每天24小時，不跨天的話12*20*（24-9）==3600，而跨天是12*（20*24-9）==5652

另外，降雨的資料中有“NR”，這個要自己處理轉換一下，用bool？不行，非NR是浮點數，NR是不降雨，NR標0就行了！！！

第二版本：

k=9，k是打算用來預測的前邊的input，因為test集提供了9個，所以為了最大化自己的預測，肯定不會去用5個預測！k以test的最大數為準！！

大迴圈：

    按日期的不同，來做分割

    中迴圈：

        把所有不同的指標都讀一遍，分別加標記以區分

        小for迴圈，i從時刻0到時刻n-9迭代（停止指標，i+9沒資料了）：

            抽i到i+（k-1）作為input，抽i+k為output

            把抽到的資料想好一個格式存起來。

第三版本，參考了作業補充ppt

需要跨天的，他也沒演示跨月怎麼處理，我估計每個月都算一個獨立的資料集？？？可能他簡化了，不用？用！！，比如你可以用隨機森林什麼的，把多個月拼接起來

跨天的話還面臨另一個建模問題，時段問題（其實不跨天也有這問題，所以這不是問題）

特徵抽取虛擬碼：

最外層迴圈-月份i迭代：

    每個月i內，每天之間要拼接起來。形成小時序列

    第二層迴圈，按小時i迭代抽取特徵：

        從j到j+8取出來，作為input儲存（train_x.append）這8個組合起來才是1個input，不是8個input，每一個都是向量，所以intpu是矩陣？

        j+9作為output，儲存（train_y.append）

        直到j+10沒資料停止

    直到月份取光結束

整個train_x，每筆加個bias

實操：參考程式碼片段，有分月份，只不過分到第二步去做了！！！

模型建立

建模

0.真正做預測，我肯定考慮要跨天了，但是操作也挺麻煩的，而且他的資料是斷的，每個月只有20天（剩下10天可能當測試集了），每個月至少還要做個截斷，操作繁瑣。

1.特徵選擇，我肯定先用純PM2.5了

2.然後可以考慮PM2.5+PM10.

3.然後我會考慮風，因為北京就和風有關

4.可以考慮一下降雨因素

5.然後考慮所有的都用？要不要正則化？

6.每天的不同時段明顯不一樣，也就是說，半夜那部分的資料，肯定影響你對白天那部分的預測？

既然是機器學習，可能就都考慮進去了，資料足夠多，模型足夠大，應該所有部分都能正確預測？

PS：應該西風或者西北風好點？

這都不絕對，因為你不知道臺灣的汙染源分佈？

也許受大陸影響，西北風會更髒？

東邊的海風也許更乾淨？

也許因為中部有山，所以東風不好用？

所以這個建模，跟實際地形、外部環境，等等，都有關係。你不瞭解詳細資訊，不會建模建的很好。

訓練

之前把資料都處理好了，現在是分N批（月份），每批一堆input和一堆output。開始寫模型訓練。

資料處理轉存完，應該可以直接訓練了。把建模部分想好的模型的公式列出來，然後用這個公式去訓練，然後用訓練集輸入一遍，然後把test輸入，輸出，做準確率對比。

第一版：

只算所有PM2.5，

隨意初始化weights（多少層？）

模型：y = bias+w1*x1+w2*x2+...+w9*x9（深層？一層？）

訓練過程：梯度下降：

η=0.01（先隨便吧）

epoch = 1000

for e in epoch:

    還是要算出Loss才能進行下邊的update的（主要問題就是這個gradient在python怎麼表達）

    

       #theta = gradientDescent(X, y, theta, alpha, iterations);#octave程式碼

        #theta = theta - alpha/m * X' *(X*theta - y)#這是octave的，X'好像是逆之類的，而且這裡X可能是矩陣，所以應該是通過公式轉換了

    批量update：

        wij=wij-η*gradient#每個wij應該單獨有一個gradient吧？

        ......

第二版：參考答案

(Pseudo code)

1. 宣告weight vector、初始learning rate、# of iteration

2. for i_th iteration :

3. y’ = train_x 和 weight vector 的 內積

4. L = y’ - train_y

5. gra = 2*np.dot( (train_x)’ , L )

6. weight vector -= learning rate * gra

模型公式示意和步驟示意（不過這裡的y'沒有+b了呢。。。grad+2*x乘以L）

測試

給了參考

# use close form to check whether ur gradient descent is good 
# however, this cannot be used in hw1.sh 
# w = np.matmul(np.matmul(inv(np.matmul(x.transpose(),x)),x.transpose()),y)

但是這個也不是分類問題，我應該怎麼比結果？取整？不是一個整數就算不準？不是吧？！

但是要的是準確了就該是個比例，而不是其他的

看一下訓練過程？那個Loss的計算方法應該可以。

下邊是照著參考PPT和程式碼片段做的

import csv
import numpy as np
from numpy.linalg import inv
#import sys
import math
#import random


#############################################
#parse data
#initial
data = []
for i in range(18):
        data.append([])
n_row = 0
#print("data:",data)
#print("type(data):",type(data))

text = open('train.csv','r',encoding='big5')
row = csv.reader(text,delimiter=",")

print("type(text):",type(text))
print("type(row):",type(row))


for r in row:
        #0 cols
        if n_row != 0:
                for i in range(3,27):
                        if r[i] != "NR":
                                data[(n_row-1)%18].append(float(r[i]))
                        else:
                                data[(n_row-1)%18].append(float(0))
        n_row = n_row + 1
text.close()
print("data[0]:",data[0][:12])
print("data[1]:",data[1][:12])


#divide data
x = []
y = []
const_month = 12
src_data_in_month = int(len(data[0]) / const_month)     #480=12*20
data_per_month = src_data_in_month - 9          #471
#per month
for i in range(12):
        for j in range(data_per_month):
                x.append([])
                for t in range(18):
                        for s in range(9):
                                x[471*i+j].append(data[t][480*i+j+s])
                y.append(data[9][480*i+j+9])    #s=9
print("after append,x[0]:",x[0][:10])
print("len(x):",len(x))

x=np.array(x)
y=np.array(y)
print("len(x):",len(x))

x=np.concatenate((np.ones((x.shape[0],1)),x),axis=1)
print("len(x):",len(x))
w2 = np.matmul(np.matmul(inv(np.matmul(x.transpose(),x)),x.transpose()),y)

################################################
#train
w = np.zeros(len(x[0]))
print("len(x[0]):",len(x[0]))
l_rate=10
repeat=10000
#print("w[0]:",w[0][:10])
#print("w[1]:",w[1][:10])

x_t = x.transpose()
s_gra = np.zeros(len(x[0]))

for i in range(repeat):
        hypo = np.dot(x,w)
        loss = hypo - y
        cost = np.sum(loss**2)/len(x)
        cost_a = math.sqrt(cost)
        gra = np.dot(x_t,loss)
        s_gra += gra**2
        ada = np.sqrt(s_gra)
        w = w - l_rate * gra/ada
        print('iteration:%d | Cost: %f ' %(i,cost_a))

np.save('model.npy',w)

w2 = np.load('model.npy')


#######################################################
#test
test_x = []
n_row = 0
text = open('test.csv',"r")
row = csv.reader(text, delimiter=",")

for r in row:
        if n_row % 18 == 0:
                test_x.append([])
                for i in range(2,11):
                        test_x[n_row//18].append(float(r[i]))
        else:
                for i in range(2,11):
                        if r[i] != "NR":
                                test_x[n_row//18].append(float(r[i]))
                        else:
                                test_x[n_row//18].append(0)
        n_row = n_row + 1
text.close()
test_x = np.array(test_x)
test_x = np.concatenate((np.ones((test_x.shape[0],1)),test_x),axis=1)


#######################################################3
#test
ans = []
for i in range(len(test_x)):
        ans.append(["id_"+str(i)])
        a = np.dot(w,test_x[i])
        ans[i].append(a)
#############
'''
        a = np.dot(w,test_x[i])
        a2 = np.dot(w2,test_x[i])
        loss = a - a2
        cost = np.sum(loss**2)/len(test_x[i])
        cost_a = math.sqrt(cost)
        print("test Cost:%f"%(cost_a))
'''
########################
filename = "predict.csv"
text = open(filename,"w+")
s = csv.writer(text,delimiter=',',lineterminator='\n')
s.writerow(["id","value"])
for i in range(len(ans)):
        s.writerow(ans[i])
text.close()

順便加入了自己寫的和“標準答案”對比的程式碼，因為這是他們臺大的作業，kaggle已經過期，沒法直接拿去跑，只能用這個答案對比。為什麼能比？這個模型他給的就是線性的，凸函式，可以用close form（翻譯是直接公式法？）來找到答案。

import csv
import math

ans = []
n_row = 0
text = open('ans.csv',"r")
row = csv.reader(text,delimiter=",")
print(text)
print(row)
#print("len(row):",len(row))
for r in row:
        #print("r[0] is ",r[0])
        #print("r[1] is ",r[1])
        if r[1] != 'value':
                ans.append(round(float(r[1])))
print(ans)

predict = []
n_row = 0
text = open('predict.csv',"r")
row = csv.reader(text,delimiter=",")
for r in row:
        #print("r[0] is ",r[0])
        #print("r[1] is ",r[1])
        if r[1] != 'value':
                predict.append(round(float(r[1])))
print(predict)
print(type(predict))


sum = 0
for i in range(len(ans)):
        sum += abs(predict[i] - ans[i])
err = sum / len(ans)
#err = abs(predict - ans).sum()/len(ans)
print(err)

方法二：array操作，最後求和平均

#######################################
#implementation 2

predict_ndarray = np.array(predict)
ans_ndarray = np.array(ans)
err_ndarray = predict_ndarray - ans_ndarray
err_ndarray = np.abs(err_ndarray)
sum2 = np.ndarray.sum(err_ndarray)
print(sum2)
err2 = sum2 / len(ans)
print(err2)