一. 線性回歸是什么？
線性回歸就是線性的回歸。線性是形容詞，回歸是本質。
我對於視覺記憶比較深刻，所以我們先上圖。
這張圖就是一個線性回歸的實例，紅色的點是實際的值，藍色為估計的線性方程
我們回歸的目的就是研究橫坐標和縱坐標的關系，當然我們首先考慮這個關系是不是線性的，換句話說這些點關系可不可以用多項式表示
w, b 分別是直線的斜率和截據，也是線性回歸最終需要獲取的結果。
這張圖是線性回歸最簡單的形式，一維，只有一個自變量，一個特征（Feature）
但是現實生活中，並不是所有的東西都只有一個特征，可能是好幾個特征決定一個結果
例如，成績總分是由所有學科的分數相加，各個學科就是不同的特征，總分就是最終想要的結果，並不能用單個成績來預測總分
線性回歸的公式是：
用成績來說，語文x1，數學x2，英語x3三門學科成績為輸入，總分y為輸出
那么可以得到這樣的模型 
（w1, w2, w3 均為1）
二. 線性回歸模型評估
評估函數的建立
模型建立完成后，我們是不知道它是不是真的優秀。
想要知道模型是否優秀，就需要對模型進行評估度量。
評估是什么意思呢，就是預測值y_preddict和我們真實數據y的差距。通過這個值的大小來判斷模型的好壞。
機器學習代碼中經常看到的Loss損失值，就是我們的評估度量模型的函數，輸入預測值和真實值，輸出損失
在統計學中，有很多度量的方法，但是統計學幾乎忘沒了 T_T 要慢慢地多掌握些統計內容。
目前我接觸到線性回歸使用最多的是平方差之和
還是先通過直觀的案例認識一下什么是平方差之和
平方差之和就是每個真實點到預測直線之間距離的平方之和，每個紅點到藍線的距離的平方 累加
平方差之和越大，那么真實值距離預測直線越遠，那么這個模型就不好
所以我們希望這個平方差之和是越小越好的，這個思想就是最小二乘法
使用公式可以表示為
基礎的損失模型建立后，可以加入正則化部分（regularization）P(w)
P(w) 通常是w的第一范式或第二范式
第一范式為所有參數和
第二范式為所有參數平方和
為什么要加入正則化部分呢，主要是為了讓預測曲線更加的平滑，讓更多的參數接近於0
三. 線性回歸模型求解
建立好模型，下一步我們就要找到最好的參數w*,b*，使我們的損失函數Loss最小
這些參數的初始值是我們設定的，目前我看到的要么是設0要么是設1.
顯然，0，1都不可能是我們真實數據的最佳參數
我們需要通過Loss(w,b)函數求使 Loss值最小時對應的 w 和 b 值
那么使用什么方法呢？
假如是高中的數學題，那么我們下意識就想要對w,b求導，然后令求導式為零，就得到w,b
但是，在真實數據中，我們不可能得到求導式，讓它為0，那么該怎么操作呢？
GradientDescent！就是他，梯度下降！
梯度在一元線性回歸中可簡化為斜率
先給出一張圖，假設它是損失值隨參數w變化的曲線，我們來模擬一下梯度下降的過程。
在點A處，斜率是小於零的，可知有更小值在A右側。於是更新參數值，使它向更小值靠近。
參數更新按照下述公式，梯度小於零，乘一個負數就是使w值增大，即向右側移動。
#更新函數 
為Loss對w 在A點的偏微分
點B處，斜率是大於零，可知有更大值在B的左側。於是更新參數值，按照更新函數，w減小，即向左移動
那按照上面這樣的想法，是不是總能找到最小值呢？
我們來看看C,D兩點，采用梯度下降只能找到極小值的位置，卡在極小值到達不了最小值
所以，使用梯度下降時，如果初始點沒選好或是學習率設置太大，是根本找不到最佳解的
如果學習率過高，已A為例，有可能會跳過最低點，飛到很遠
高中數學中令求導項為0，在計算機中可以通過大量的進行獲得
更新關鍵就是獲得Loss對w的偏微分。
於是參數更新函數為，這里的參數2可以加也可以不加
確定好參數更新函數后，接下來做的就是大量的循環，暴力求解啦
四：實際案例
一元線性回歸
單純的X，Y線性關系，畫出的散點圖就是上面做案例的圖
回歸效果：
下面是w,b參數在訓練中的變化情況
代碼部分：
首先，將數據加載進來
#1.加載數據 data_source = "data/fire_theft.xls" #excel表形式 book = xlrd.open_workbook(data_source,encoding_override="utf-8") #通過索引獲取內容 sheet = book.sheet_by_index(0) #print(sheet) #讀取每一行，將每一行內容提取作為list,再將所有list作為np.array存儲 data = np.asarray([sheet.row_values(i) for i in range(1,sheet.nrows)]) #print(data) x = [] y = [] for i in range(len(data)): [x_in,y_in] = data[i] x.append(x_in) y.append(y_in)
將數據分割成訓練集，和測試集
#2.數據分割 #選擇%數據進行訓練，其中70%Training Set 30%Test Set train = int(len(x)*0.7) test = int(len(x)*0.3) x_train = x[0:train] y_train = y[0:train] #[train+valid:-1]無法讀取最后一個數 a:b 讀取到b-1位置停止，若a: 沒有指定尾坐標，直接取到最后一個 x_test = x[train+valid:] y_test = y[train+valid:]
進行訓練，這里使用的是SGD(Stochstic Gradient Descent), 不是獲取全部數據后更新，而是獲取一個數據就更新一次，這樣計算時間更快
#3.1 使用Gradient Descent 求解 y = w*x + b #記錄每次迭代的w,b值 w_history = [] b_history = [] lr_w = 0.0 lr_b = 0.0 for i in range(iteration): w_grad = 0.0 b_grad = 0.0 Loss = 0.0 for j in range(len(x_train)): #損失函數 y_pred = w*x_train[j] +b Loss = Loss + np.square(y_pred-y_train) #梯度下降 w(i+1) = w(i) - Loss偏微分 w_grad = w_grad - 2.0*(y_train[j] - y_pred)*1.0 b_grad = b_grad - 2.0*(y_train[j] - y_pred)*x_train[j] #adagram 更新 lr_w = lr_w + w_grad**2 #w = w - lr/np.sqrt(lr_w) * w_grad lr_b = lr_b + b_grad**2 w = w - lr * w_grad b = b - lr * b_grad w_history.append(w) b_history.append(b) #3.2.展示回歸數據 y_end=[] for i in range(len(x_train)): y_end.append(x_train[i]*w + b) #print('w: %f , b: %f' %(w,b))
plt.plot(x_train,y_train,'ro') plt.plot(x_train,y_end) plt.title("Self GD w:%f b:%f" %(w,b)) plt.show()
plt.plot(w_history,label = 'w',color ='r') plt.plot(b_history,label = 'b',color = 'b') plt.title("the run chart of w and b") plt.legend() #加入左上角顯示框 plt.show() #3.3 測試數據 y_test_end=[] y_test_loss = 0.0 for i in range(len(x_test)): y_test_end.append(x_test[i]*w + b) y_test_loss += np.square(y_test[i] - y_test_end[i]) #print('w: %f , b: %f' %(w,b)) plt.plot(x_test,y_test,'ro') plt.plot(x_test,y_test_end) plt.title("Self GD w:%f b:%f,loss:%f" %(w,b,y_test_loss/len(x_test))) plt.show()
多元線性回歸
使用的是李宏毅老師的PM2.5預測
將不同的檢測指標作為特征 xi ，預測值y 為PM2.5值
先看線性回歸對PM2.5預測的效果
紅色點：實際PM2.5值
紅色曲線：真實PM2.5變化曲線
藍色曲線：模型預測曲線
能看出，藍色曲線能模擬基本的趨勢，第二張圖就能說明，大致趨勢是相同的
測試集的損失值為39，效果並不是很好
代碼：
import numpy as np import pandas as pd import matplotlib.pyplot as plt #------PM2.5線性回歸預測-------------- #1.數據讀取、處理 file = open('data/weather_train.csv',encoding = 'gb18030') #'gbk' codec can't decode byte 0xac in position 9: illegal multibyte sequence file = pd.read_csv(file, usecols=range(3, 27)) # usecols #變為函數 def get_data(file): f = file.replace(['NR'], [0.0]) data = f.values.astype(float) # 分析數據，提取數據 # 需要提取240天 每天24 - 9 = 15 個數據集 進行時間片的預測 # 首先取數據模塊，18個數據為一組 x_list = [] y_list = [] for i in range(0, data.shape[0]-18, 18): # 取行 18 行 for j in range(24 - 9): # 循環取值 取15次 外圍控制輸入為3開始 temp_x = data[i:i + 18, j:j + 9] temp_y = data[i + 9, j + 9] x_list.append(temp_x) y_list.append(temp_y) x = np.array(x_list) y = np.array(y_list) return x,y #2.模型設計 不用tensorflow #使用SGD def train_model(x,y,learning_rate=1,epoch=1000): #設置初始參數 weight = np.ones(9) bias = 0.0 #正則參數 regu_rate = 0.001 w_sum = np.zeros(9) #所有訓練綜合，因此放在最外面 b_sum = 0.0for ep in range(epoch): for i in range(len(x)): w_in = (y[i] - weight.dot(x[i, 9, :]) - bias) * (-x[i, 9, :]) b_in = (y[i] - weight.dot(x[i, 9, :]) - bias) * (-1) w_sum += w_in ** 2 b_sum += b_in ** 2 # 進行函數更新 weight += - learning_rate / w_sum ** 0.5 * (w_in + regu_rate*np.sum(weight)) #加入regu_rate 正則化函數 bias -= learning_rate / b_sum ** 0.5 * b_in w_out.append(weight[0]) b_out.append(bias) # 輸出損失值 loss = 0.0 for m in range(len(x)): loss += (y[m] - weight.dot(x[m,9,:]) - bias)**2 print("Epoch:%d loss:%f"%(ep,loss/len(x))) x,y = get_data(file) x_train = x[0:3200] y_train = y[0:3200] x_valid = x[3200:] y_valid = y[3200:] weight,bias=train_model(x_train,y_train,epoch=1000) #測試集訓練 y_pred = [] loss_pred = 0.0 for i in range(len(x_valid)): y_pred.append( weight.dot(x_valid[i,9,:]) + bias) loss_pred+=(y_valid[i] - y_pred[i])**2 print("the cost of test:",loss_pred/len(x_valid)) plt.plot(range(len(y_valid)),y_valid,'r') plt.plot(range(len(y_valid)),y_pred,'b') plt.show()
第一次博客園發博成就 get