李巨集毅機器學習HM1（梯度下降）
問題引入
作業所給的資料是某地的觀測記錄，每個月取前20天的資料，觀測資料共有18個指標，每小時記錄這18個指標的值，共記錄12個月。
現在從剩下的資料中取出連續的9小時的觀測資料，請預測第10個小時的PM2.5指標的值。
資料處理
先將csv檔案內容讀入進來，首先需要注意的是RAINFALL指標還有NR，得把它替換成0。
1 def read_csv(): 2 df = pd.read_csv(CSV_file_path) 3 df = df.iloc[:, 3:] # 前面3列的內容沒用 4 df.replace('NR', 0, inplace=True) # 將降雨量中的NR替換為0 5 return df
接下來生成測試集，因為最後預測資料的特徵是連續9小時的觀測資料，所以我們每次取連續的10小時觀測資料，前9小時的資料作為特徵，最後一小時的PM2.5指標觀測值作為標籤值。
為方便生成資料，先將資料作如下處理：
一共生成12行這樣的資料，注意每個月不能連起來，因為每個月只取了前20天，並不連續。
接下來生成資料即可，程式碼如下：
1 def get_train_data(df): 2 data = df.to_numpy() 3 month_data = {}4 5 for month in range(12): 6 sample = np.empty([18, 480]) # 一共18個觀測指標，24*20=480 7 for day in range(20): # 將每個月的20天資料連線起來 8 sample[:, day*24:(day+1)*24] = data[month*20*18+day*18:month*20*18+(day+1)*18, :] 9 month_data[month] = sample 10 11 x_set = np.empty((12*471, 18*9)) #每10小時可取出一組資料，共471組，一共12個月，所以總的資料量為12*471，屬性值一共有18*9 12 y_set = np.empty((12*471, 1)) 13 14 for month in range(12): 15 for day in range(20): 16 for hour in range(24): 17 if day == 19 and hour > 14: 18 continue 19 x_set[month * 471 + day * 24 + hour, :] = month_data[month][:, day * 24 + hour: day * 24 + hour + 9].reshape(1, -1) # 將資料 重組成一行 20 y_set[month * 471 + day * 24 + hour, 0] = month_data[month][9, day * 24 + hour + 9] # 只取第9行的PM2.5觀測值 21 return x_set, y_set
為了能獲得更好的效果，在梯度下降之前先進行標準化處理。
這裡使用的公式為z-score 標準化：
$x_{i}^{r} = \frac{x_{i}^{r}-m_{i}}{\sigma_i}$
其中$x_{i}^{r} $表示第$r$組資料中的第$i$個值，$m_{i}$表示所有樣本中第$i$個值的均值，$\sigma_i$表示所有樣本中第$i$個值的標準差。
1 def normalization(x_set): 2 x_mean = np.mean(x_set, axis=0) 3 x_std = np.std(x_set, axis=0) 4 for row in range(len(x_set)): 5 for col in range(len(x_set[0])): 6 if x_std[col] != 0: # 標準差為0表示資料基本無波動 7 x_set[row][col] = (x_set[row][col] - x_mean[col])/x_std[col] 8 return x_set, x_mean, x_std
在得到訓練資料後，我們再拆分一部分資料出來作為驗證集。
1 def split_data(x_set, y_set): 2 x_train_set = x_set[: math.floor(len(x_set) * 0.8), :] 3 y_train_set = y_set[: math.floor(len(y_set) * 0.8), :] 4 x_validation = x_set[math.floor(len(x_set) * 0.8):, :] 5 y_validation = y_set[math.floor(len(y_set) * 0.8):, :] 6 return x_train_set, y_train_set, x_validation, y_validation
梯度下降
我們首先假設線性迴歸的函式為：
$H(x)=w_{0}+w_{1}x_{1}+w_{2}x_{2}+···+w_{n}x_{n}$
損失函式使用均方誤差，李巨集毅老師的示例作業裡使用了均方根誤差，但是程式碼中算梯度的式子我沒看懂，還請懂的人可以告訴我下，所以這裡我就用了均方誤差。
$L(w)=\frac{1}{2m}\sum_{i=1}^{m}(H(x^{(i)})-y^{(i)})^{2}$
現在對引數$w_{j}$求偏導，即
$\frac{\partial L(w)}{w_{j}}=\frac{1}{m}\sum_{i=1}^{m}(H(x^{(i)})-y^{(i)})x_{j}^{(i)}$
這式子是可以轉換成矩陣運算的，具體的請看下面的程式碼。
在得到梯度之後，需要不斷更新$w$引數的值，在這裡使用李巨集毅老師所講的Adagrad方法，公式為：
$w^{t+1} = w^{t}-\frac{\eta^{t}}{\sigma^{t}}g^{t}$
其中$\eta^{t}=\frac{\eta}{\sqrt{t+1}}$，這可以使得我們一開始以儘量快得到速度靠近目標，然後逐漸減小學習率，$g^{t}=\frac{\partial L(w)}{w}$，\sigma^{t}=\sqrt{\frac{1}{t+1}\sum_{i=0}^{t}(g^{i})^2}。
經過化簡計算，即為：
$w^{t+1} = w^{t}-\frac{\eta}{\sqrt{\sum_{i=0}^{t}(g^{i})^2}}g^{t}$
1 def training(x_train_set, y_train_set): 2 dim = 18 * 9 + 1 # w引數的維度，+1是可以把b也當成一個w 3 w = np.ones([dim, 1]) 4 b = np.ones([len(x_train_set), 1]) 5 x = np.concatenate((b, x_train_set), axis=1).astype(float) # 將b初始化為1，載入樣本屬性值的最前面 6 7 learning_rate = 100 # 學習率 8 iter_time = 20000 # 迭代次數 9 adagrad = np.zeros([dim, 1]) 10 eps = 0.0000000001 # 新的學習率是learning_rate/sqrt(sum_of_pre_grads**2),而adagrad=sum_of_grads**2,所以處在分母上而迭代時adagrad可能為0，所以加上一個極小數，使其不除0 11 12 for i in range(iter_time): 13 loss = np.sum(np.power(np.dot(x, w) - y_train_set, 2)) / len(x_train_set)/2 # 均方誤差 14 if i % 100 == 0: # 每迭代100次輸出loss值 15 print(str(i)+':'+str(loss)) 16 gradient = np.dot(x.T, np.dot(x, w) - y_train_set)/len(x_train_set) # 計算梯度 17 adagrad += gradient ** 2 # 累加adagrad值 18 w = w - learning_rate * gradient / np.sqrt(adagrad+eps) # 更新引數 19 return w
驗證預測
將訓練得到的$w$引數在之前拆分得到的驗證集上進行計算。
1 def validation(w, x_validation, y_validation): 2 b = np.ones([len(x_validation), 1]) 3 x = np.concatenate((b, x_validation), axis=1).astype(float) 4 loss = np.sum(np.power(np.dot(x, w) - y_validation, 2)) / len(x_validation)/2 5 print('驗證的loss為' + str(loss))
最後進行預測，讀入test.csv檔案，也和訓練資料一樣，先處理一下資料，然後直接np.dot(x,w)即可。
1 def predict(w, x_mean, x_std): 2 df = pd.read_csv(test_file_path, header=None) 3 df = df.iloc[:, 2:] 4 df.replace('NR', 0, inplace=True) 5 data = df.to_numpy() 6 pre_data = np.empty((240, 18*9)) 7 for i in range(240): 8 pre_data[i, :] = data[18*i:18*(i+1), :].reshape(1, -1) 9 for row in range(len(pre_data)): # 需要標準化，而且均值和標準差需要使用之前的 10 for col in range(len(pre_data[0])): 11 if x_std[col] != 0: 12 pre_data[row][col] = (pre_data[row][col] - x_mean[col]) / x_std[col] 13 b = np.ones([len(pre_data), 1]) 14 pre_data = np.concatenate((b, pre_data), axis=1).astype(float) 15 result = np.dot(pre_data, w) 16 17 file = open('result.csv', 'w') 18 for i in range(240): 19 file.write('id_' + str(i) + ',' + str(result[i][0])) 20 file.write('\n') 21 file.close()
預測的結果為：