機器學習-最小二乘法

1、背景知識

「最小二乘法」出現的歷史背景是頗有意思的。（如下文字摘錄維基百科）

1801年，意大利天文學家朱賽普·皮亞齊發現了第一顆小行星穀神星。通過40天的跟蹤觀測後，因爲穀神星運行至太陽背後，使得皮亞齊失去了穀神星的位置。隨後全世界的科學家利用皮亞齊的觀測數據開始尋找穀神星，可是根據大多數人計算的結果來尋找穀神星都沒有結果。時年24歲的高斯也計算了穀神星的軌道。奧地利天文學家海因裏希·奧爾伯斯根據高斯計算出來的軌道從新發現了穀神星。

高斯使用的最小二乘法的方法發表於1809年他的著做《天體運動論》中，而法國科學家勒讓德於1806年獨立發現「最小二乘法」，但因不爲時人所知而默默無聞。兩人曾爲誰最先創立最小二乘法原理髮生爭執。

1829年，高斯提供了最小二乘法的優化效果強於其餘方法的證實，見高斯-馬爾可夫定理。

2、知識運用

「最小二乘法」的核心就是保證全部數據誤差的平方和最小。（「平方」的在古時侯的稱謂爲「二乘」）

假設咱們收集到一些戰艦的長度與寬度數據

	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畫出散點圖：

畫散點圖的代碼以下：

# -*- 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=0.197,b=-14.48

那樣的話，咱們能夠獲得這樣的擬合圖：

好了，新的問題來了，這樣的a,b是否是最優解呢？用專業的說法就是：a,b是否是模型的最優化參數？在回答這個問題以前，咱們先解決另一個問題：a,b知足什麼條件纔是最好的？答案是：保證全部數據誤差的平方和最小。至於原理，咱們會在後面講，先來看看怎麼利用這個工具來計算最好的a和b。

假設全部數據的平方和爲M，則

咱們如今要作的就是求使得M最小的a和b。請注意這個方程中，咱們已知yⁱ和xⁱ

那其實這個方程就是一個以（a,b）爲自變量，M爲因變量的二元函數。

回想一下高數中怎麼對一元函數就極值。咱們用的是導數這個工具。那麼在二元函數中，

咱們依然用導數。只不過這裏的導數有了新的名字「偏導數」。偏導數就是把兩個變量中的一個視爲常數來求導。

經過對M來求偏導數，咱們獲得一個方程組

=0

=0

這兩個方程中x_i和y_i都是知道的。

很容易就求得a和b了。因爲採用的是維基百科的數據，我這裏就直接用答案來畫出擬合圖像：

# -*- 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()

3、原理探究

數據擬閤中，爲何要讓模型的預測數據與實際數據之差的平方而不是絕對值和最小來優化模型參數？

這個問題已經有人回答了，見連接（http://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] * ..

又由於：e^a_{*e^b*e^c=e^a+b+c}

因此p(D|h)∝ EXP{-[(ΔY1)^2 +(ΔY2)^2 + (ΔY3)^2 + ..]}

咱們知道f(x)=e^x的分佈圖像爲：

 由於e的指數永遠小於0，因此，想要p(D|h)最大，就必須使[(ΔY1)^2 + (ΔY2)^2 + (ΔY3)^2 + ..]無限接近於0，即：最大化p(D|h)就是要最小化[(ΔY1)^2 + (ΔY2)^2 + (ΔY3)^2 + ..]

到此，最小二乘法的原理獲得了詮釋。

4、最小二乘法與梯度降低法

   最小二乘法跟梯度降低法都是經過求導來求損失函數的最小值，那它們有什麼區別呢。

   相同
　　1.本質相同：兩種方法都是在給定已知數據（independent & dependent variables）的前提下對dependent variables算出出一個通常性的估值函數。而後對給定新數據的dependent variables進行估算。
　　2.目標相同：都是在已知數據的框架內，使得估算值與實際值的總平方差儘可能更小（事實上未必必定要使用平方），估算值與實際值的總平方差的公式爲：

                             

   其中爲第i組數據的independent variable，爲第i組數據的dependent variable，爲係數向量。
   不一樣
　　1.實現方法和結果不一樣：最小二乘法是直接對求導找出全局最小，是非迭代法。而梯度降低法是一種迭代法，先給定一個，而後向降低最快的方向調整，在若干次迭代以後找到局部最小。梯度降低法的缺點是到最小點的時候收斂速度變慢，而且對初始點的選擇極爲敏感，其改進大可能是在這兩方面下功夫。

5、拓展延伸

上面講的都是二維的狀況，也就是隻有一個自變量。但現實世界中影響最後結果的都是多種因素的疊加，即自變量會有多個的狀況。

對於通常N元線性函數，用《線性代數》中的逆矩陣來求解就OK了；因爲暫時沒有找到合適的例子，就做爲一個引子，留在這裏了。

固然天然界更多的是多項式擬合，而非簡單的線性，那就是更高級的內容了。