445
線性回歸之——最小二乘法
一、引言
這段時間學習《機器學習》，學到第5章的“Logistic回歸”，感覺相當吃力。追本溯源，從“Logistic回歸”到“線性回歸”，再到“最小二乘法”。最終定格到了《高等數學》（第六版·下冊）第九章第10節“最小二乘法”，這才了解到最小二乘法背後的數學原理出自哪裏。
“最小二乘法”是最優化問題中建立經驗公式的一種實現方法。了解它的原理，對於了解“Logistic回歸”和“支持向量機的學習”都很有裨益。
二、背景知識
“最小二乘法”出現的曆史背景是很有意思的。（以下文字摘錄維基百科）
1801年，意大利天文學家朱賽普·皮亞齊發現了第一顆小行星穀神星。經過40天的跟蹤觀測後，由於穀神星運行至太陽背後，使得皮亞齊失去了穀神星的位置。隨後全世界的科學家利用皮亞齊的觀測數據開始尋找穀神星，但是根據大多數人計算的結果來尋找穀神星都沒有結果。時年24歲的高斯也計算了穀神星的軌道。奧地利天文學家海因裏希·奧爾伯斯根據高斯計算出來的軌道重新發現了穀神星。
高斯使用的最小二乘法的方法發表於1809年他的著作《天體運動論》中，而法國科學家勒讓德於1806年獨立發現“最小二乘法”，但因不為時人所知而默默無聞。兩人曾為誰最早創立最小二乘法原理發生爭執。
1829年，高斯提供了最小二乘法的優化效果強於其他方法的證明，見高斯-馬爾可夫定理。
三、知識運用
“最小二乘法”的核心就是保證所有數據偏差的平方和最小。（“平方”的在古時侯的稱謂為“二乘”）
假設我們收集到一些戰艦的長度與寬度數據
|
1 |
2 |
3 |
4 |
5 |
6 |
7 |
8 |
9 |
10 |
|
|
長度（m） |
208 |
152 |
113 |
227 |
137 |
238 |
178 |
104 |
191 |
130 |
|
寬度(m) |
21.6 |
15.5 |
10.4 |
31.0 |
13.0 |
32.4 |
19.0 |
10.4 |
19.0 |
11.8 |
根據這些數據我們用python畫出散點圖：
畫散點圖的代碼如下：
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
|
# -*- coding: utf-8 -*import numpy as np
import os
import matplotlib.pyplot as plt
def drawScatterDiagram(fileName):
#改變工作路徑到數據文件存放的地方
os.chdir("d:/workspace_ml")
xcord=[];ycord=[]
fr=open(fileName)
for line in fr.readlines():
lineArr=line.strip().split()
xcord.append(float(lineArr[1]));ycord.append(float(lineArr[2]))
plt.scatter(xcord,ycord,s=30,c='red',marker='s')
plt.show()
|
假如我們取前兩個點（238,32.4）（152, 15.5）就可以得到兩個方程
152*a+b=15.5
328*a+b=32.4
那樣的話，我們可以得到這樣的擬合圖：
好了，新的問題來了，這樣的a,b是不是最優解呢？用專業的說法就是：a,b是不是模型的最優化參數？在回答這個問題之前，我們先解決另外一個問題：a,b滿足什麼條件才是最好的？答案是：保證所有數據偏差的平方和最小。至於原理，我們會在後麵講，先來看看怎麼利用這個工具來計算最好的a和b。
假設所有數據的平方和為M，則
我們現在要做的就是求使得M最小的a和b。請注意這個方程中，我們已知yi和xi
那其實這個方程就是一個以（a,b）為自變量，M為因變量的二元函數。
回想一下高數中怎麼對一元函數就極值。我們用的是導數這個工具。那麼在二元函數中，
我們依然用導數。隻不過這裏的導數有了新的名字“偏導數”。偏導數就是把兩個變量中的一個視為常數來求導。
通過對M來求偏導數，我們得到一個方程組
=0
=0
這兩個方程中xi和yi都是知道的。
很容易就求得a和b了。由於采用的是維基百科的數據，我這裏就直接用答案來畫出擬合圖像：
代碼如下：
|
1
|
# -*- coding: utf-8 -*importnumpy as npimportosimportmatplotlib.pyplot as pltdefdrawScatterDiagram(fileName):#改變工作路徑到數據文件存放的地方os.chdir("d:/workspace_ml")xcord=[];ycord=[]fr=open(fileName)forline infr.readlines():lineArr=line.strip().split()xcord.append(float(lineArr[1]));ycord.append(float(lineArr[2]))plt.scatter(xcord,ycord,s=30,c='red',marker='s')#a=0.1965;b=-14.486a=0.1612;b=-8.6394x=np.arange(90.0,250.0,0.1)y=a*x+bplt.plot(x,y)plt.show() |
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
|
# -*- coding: utf-8 -*import numpy as np
import os
import matplotlib.pyplot as plt
def drawScatterDiagram(fileName):
#改變工作路徑到數據文件存放的地方
os.chdir("d:/workspace_ml")
xcord=[];ycord=[]
fr=open(fileName)
for line in fr.readlines():
lineArr=line.strip().split()
xcord.append(float(lineArr[1]));ycord.append(float(lineArr[2]))
plt.scatter(xcord,ycord,s=30,c='red',marker='s')
#a=0.1965;b=-14.486
a=0.1612;b=-8.6394
x=np.arange(90.0,250.0,0.1)
y=a*x+b
plt.plot(x,y)
plt.show()
|
四、原理探究
數據擬合中，為什麼要讓模型的預測數據與實際數據之差的平方而不是絕對值和最小來優化模型參數？
這個問題已經有人回答了，見鏈接（https://blog.sciencenet.cn/blog-430956-621997.html）
個人感覺這個解釋是非常有意思的。特別是裏麵的假設：所有偏離f(x)的點都是有噪音的。
一個點偏離越遠說明噪音越大，這個點出現的概率也越小。那麼偏離程度x與出現概率f(x)滿足什麼關係呢？——正態分布。
已知N個點（用D來表示），求直線（用h來表示）出現的概率就可以表示為：P(h|D)
根據貝葉斯定理：P(h|D)=P(D|h)*P(h)/P(D)即P(h|D)∝P(D|h)*P(h) (∝表示“正比於”)
這就是一個生成模型了——由直線h生成點集D。
我們再作一個假設：h生成D中的每一個點都是獨立的（如果了解貝葉斯文本分類的話，這裏就很好理解了），那麼P(D|h)=p(d1|h)*p(d2|h)…
結合前麵正態分布，我們可以寫出這樣的式子：p(di|h)∝ exp(-(ΔYi)^2)
那麼P(D|h)∝EXP[-(ΔY1)^2]* EXP[-(ΔY2)^2] * EXP[-(ΔY3)^2] * ..
又因為：ea*eb*ec=ea+b+c
所以p(D|h)∝ EXP{-[(ΔY1)^2 +(ΔY2)^2 + (ΔY3)^2 + ..]}
我們知道f(x)=ex的分布圖像為：
因為e的指數永遠小於0，所以，想要p(D|h)最大，就必須使[(ΔY1)^2 + (ΔY2)^2 + (ΔY3)^2 + ..]無限接近於0，即：最大化p(D|h)就是要最小化[(ΔY1)^2 + (ΔY2)^2 + (ΔY3)^2 + ..]
到此，最小二乘法的原理得到了詮釋。
五、拓展延伸
上麵講的都是二維的情況，也就是隻有一個自變量。但現實世界中影響最後結果的都是多種因素的疊加，即自變量會有多個的情況。
對於一般N元線性函數，用《線性代數》中的逆矩陣來求解就OK了；由於暫時沒有找到合適的例子，就作為一個引子，留在這裏了。
當然自然界更多的是多項式擬合，而非簡單的線性，那就是更高級的內容了。
參考文獻
《高等數學》（第六版）（高等教育出版社）
《線性代數》（北京大學出版社）
互動百科：最小二乘法
維基百科：最小二乘法
科學網：最小二乘法？為神馬不是差的絕對值
最後更新：2017-07-21 21:32:59