線性迴歸 李巨集毅機器學習HW1
阿新 • • 發佈：2018-11-23
本文是李巨集毅老師機器學習的第一次大作業，參考網上程式碼，寫了一下自己的思路。
李巨集毅 HM1:
要求：本次作業使用豐原站的觀測記錄，分成train set跟test set，train set是豐原站每個月的 前20天所有資料。test set則是從豐原站剩下的資料中取樣出來。 train.csv：每個月前20天的完整資料。 test_X.csv：從剩下的10天資料中取樣出連續的10小時為一筆，前九小時的所有觀測 資料當作feature，第十小時的PM2.5當作answer。一共取出240筆不重複的test data，請 根據feauure預測這240筆的PM2.5。
1review
首先回顧下線性迴歸的基本知識：
https://blog.csdn.net/ZHOUJIAN_TANK/article/details/83750105
關鍵的點：
step1:set a model ，很顯然採用線性模型
step2:goodness of function 均方誤差
step3:pick the best function Gradient descent
本次作業處理思路：
利用訓練集對引數進行學習；對test data進行預測。
關鍵的地方（1）對train data 進行資料處理 ，資料如何儲存和讀取
（python 基礎）
（2）向量的操作方法
下面補充一點向量的知識：
式中，N為訓練集的樣本數，為實際的值，x為輸入值(可以為向量），w與b為訓練引數（parament）。
2、分析：
訓練集中每個月20天，每天24個小時。根據前9個小時觀測的所有資料（共18中汙染物）， 也就是18*9=162個feature輸入來預測第10個小時的P,2.5值。由於資料的不連續性， 在連續取九個小時的過程中，每個月第20天的後9個小時是沒有對應的實際對應的測量值的， 所以實際上每個月有 20*24-9=471筆連續9小時，也就是每個月有471筆訓練資料， 一年就有 12*471=5652筆訓練資料;每筆訓練資料包含9個小時的18種汙染物18*9=162， 也就是每筆資料有162個feature，需要162個引數外加一個bais需要163各引數。
function set 可設定為下：
上式中，b為bias,w1,pm為PM的第一個小時觀測值，w2,pm為PM的第二個小時觀測值，
w3,pm為PM的第三個小時觀測值。。。。。w9,PM為PM的第9個小時觀測值weight值，。。。。。。。。。。。。以此類推，wi,cat 即為種類為cat的汙染物第i小時的觀測值的權重。。。。。。。。 18*9+1=163個paraments。
轉換成向量的形式，即為：
其中 phi為訓練資料集，5652筆訓練資料，163個引數；第一個引數為bias，其餘的為feature（汙染物）。
theta為訓練的引數向量，
3、程式碼分析
3.1對資料的處理
train資料的原始儲存形式為：
處理原始的訓練資料，方便對資料的使用。上圖為一天的所有訓練資料，為了方便操作，可以把18種汙染物的所有資料，按照種類放在一個數組中data18， 每種汙染物共202412=5670個數據也就是：原始表格轉換為：
在資料轉換為Data185670 基礎上，再取的訓練集與實際輸出值。訓練集為每個月所有可能連續9小時，共471種可能，每個訓練裡有 189=162個數據，怎麼樣取得162個數據？按照順序 [第一種汙染物連續9小時，第二種汙染物連續9小時，第三種，。。。。。。。。第18種連續9小時]。得道轉換後的訓練集資料 X5652*162
程式碼如下：
################完成原始資料按按照每種汙染物的類別進行儲存##############
data = [] #data存什麼？按類別儲存 feature 共18種特徵
for i in range(18):
data.append([])
n_row = 0
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): #取所有的小時，注意資料在xls表格中的格式
if r[i] != "NR": #對 rainfall中的資料進行處理
data[(n_row - 1) % 18].append(float(r[i]))
else:
data[(n_row - 1) % 18].append(float(0)) #若不下雨則 rain 那一行加0
n_row = n_row + 1
text.close()
以上程式碼實現原始資料到 新的表格資料的儲存轉換。
# divide data
x = []
y = []
const_month = 12
a=len(data[0]) #每項指標所有的資料
src_data_in_month = int(len(data[0]) / const_month) # 5760/12
data_per_month = src_data_in_month - 9 # 471 每種特徵資料每月儲存的個數
# per month
a=0
for i in range(12):#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]) #每月的後9小時不參與預測
y.append(data[9][480 * i + j + 9]) # s=9 ；y裡儲存的是 PM2.5,每月第一天的前九個資料不參與預測
print("after append,x[0]:", x[0])
print("len(x):", len(x))
x = np.array(x)
y = np.array(y)
print("len(x):", len(x))
以上程式碼完成X資料的儲存，如下表：
Y資料的儲存為：
x = np.concatenate((np.ones((x.shape[0], 1)), x), axis=1) #add bias 在原來x的基礎上在第一列之前新增一列每個元素值為1
第一列元素作為Base。
至此得到完成資料的處理工作。
3.2訓練資料
w = np.zeros(len(x[0])) #引數的初始值 設定為 0，共163個引數
print("len(x[0]):", len(x[0]))
l_rate = 10 #learning-rate
repeat = 10000 #迭代次數
x_t = x.transpose()
s_gra = np.zeros(len(x[0])) #grad
for i in range(repeat): #迭代次數
hypo = np.dot(x, w) #（5652,1），這一步得到是所有訓練集組成的一個預測值組成的 矩陣
loss = hypo - y #（5652,1） 預測數值-真實數值，
cost = np.sum(loss ** 2) / len(x) #loss 均方誤差和
cost_a = math.sqrt(cost) #這句話和上面三句話用來求LOSS function均方根，代表的都是資料擬合的好壞
#梯度下降
gra = np.dot(x_t, loss) #這裡不懂，應該是矩陣的導數
s_gra += gra ** 2 #sum_grad
ada = np.sqrt(s_gra)
w = w - l_rate * gra / ada #完成引數的一次迭代更新，具體參考上面的公式
print('iteration:%d | Cost: %f ' % (i, cost_a))
#save model
np.save('model.npy', w)
#read model
w2 = np.load('model.npy')
#完成迭代更新，得到訓練的引數w。
3.3測試資料
# test
test_x = []
n_row = 0
text = open('test.csv', "r")
row = csv.reader(text, delimiter=",")
for r in row: #讀取測試集的資料儲存到test_x[]中
if n_row % 18 == 0: #按指每次的第一行，也就是18的倍數行
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)
以上實現讀取測試集原始資料，完成資料的轉換。
ans = []
yy=[]
for i in range(len(test_x)):
ans.append(["id_" + str(i)])
a = np.dot(w, test_x[i])
ans[i].append(a) #得到預測值
yy.append(a)
#以上得到最終的預測值。
**
程式碼附錄：
**
#_author_:'alex zhou'
# Time: 2018/10/18 0018
import csv #讀寫CSV資料
import numpy as np
import math
import matplotlib.pyplot as plt
# parse data
# initial
data = [] #data存什麼？按類別儲存 feature 共18種特徵
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": #對 rainfall中的資料進行處理
data[(n_row - 1) % 18].append(float(r[i]))
else:
data[(n_row - 1) % 18].append(float(0)) #若不下雨則 rain 那一行加0
n_row = n_row + 1
text.close()
# print("data[0]:", data[0][:12])
# print("data[1]:", data[1][:12])
# with open('data.csv','w') as myfile:
# mywriter=csv.writer(myfile)
# mywriter.writerows(data)
# divide data
x = []
y = []
const_month = 12
a=len(data[0]) #每項指標所有的資料
src_data_in_month = int(len(data[0]) / const_month) # 5760/12
data_per_month = src_data_in_month - 9 # 471 每種特徵資料每月儲存的個數
# per month
a=0
for i in range(12):#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 ；y裡儲存的是 PM2.5,第一天的前九個資料不參與預測
print("after append,x[0]:", x[0])
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) #add bias 在原來x的基礎上在第一列之前新增一列每個元素值為1 #拼接 x.shape[0]為5652 shape函式是numpy.core.fromnumeric中的函式，它的功能是檢視矩陣或者陣列的維數。
# with open('123.csv','w') as myfile:
# mywriter=csv.writer(myfile)
# mywriter.writerows(x)
print(len(x[0])) #為163
w2 = np.matmul(np.matmul(inv(np.matmul(x.transpose(), x)), x.transpose()), y) #matual矩陣相乘
################################################
# train
w = np.zeros(len(x[0]))
print("len(x[0]):", len(x[0]))
l_rate = 10 #learning-rate
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) #（5652,1）
loss = hypo - y #（5652,1） 預測數值-真實數值
cost = np.sum(loss ** 2) / len(x) #loss
cost_a = math.sqrt(cost) #這句話和上面三句話用來求LOSS function均方根，代表的都是資料擬合的好壞
#梯度下降
gra = np.dot(x_t, loss) #這裡不懂，應該是矩陣的導數
s_gra += gra ** 2 #sum_grad
ada = np.sqrt(s_gra)
w = w - l_rate * gra / ada
print('iteration:%d | Cost: %f ' % (i, cost_a))
#save model
np.save('model.npy', w)
#read model
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: #讀取測試集的資料儲存到test_x[]中
if n_row % 18 == 0: #按指每次的第一行，也就是18的倍數行
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 = []
yy=[]
for i in range(len(test_x)):
ans.append(["id_" + str(i)])
a = np.dot(w, test_x[i])
ans[i].append(a) #得到預測值
yy.append(a)
#############
## 這一部分是把得出的 pm2.5 的數值寫入 predict.csv檔案
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()
#plot
row_n = 0
y = []
text = open('ans.csv', 'r', encoding = 'big5')
row = csv.reader(text, delimiter = ',') # 此處 delimiter = ',' 有沒有無所謂，這句話用來當一行中用 ','時分割句子的
for r in row:
if row_n != 0:
y.append(r[1])
row_n = row_n + 1
text.close()
plt.figure(figsize = (13, 7))
plt.plot(np.arange(0, 240, 1), yy, 'r', label = 'prediction pm2.5')
plt.plot(np.arange(0, 240, 1), y, 'b', label = 'ans pm2.5')
plt.legend()
plt.show()