最小二乘法及C語言實現

      我們以最簡單的一元線性模型來解釋最小二乘法。什么是一元線性模型呢？ 監督學習中，如果預測的變量是離散的，我們稱其為分類（如決策樹，支持向量機等），如果預測的變量是連續的，我們稱其為回歸。回歸分析中，如果只包括一個自變量和一個因變量，且二者的關系可用一條直線近似表示，這種回歸分析稱為一元線性回歸分析。如果回歸分析中包括兩個或兩個以上的自變量，且因變量和自變量之間是線性關系，則稱為多元線性回歸分析。對於二維空間線性是一條直線；對於三維空間線性是一個平面，對於多維空間線性是一個超平面...

   對於一元線性回歸模型, 假設從總體中獲取了n組觀察值（X1，Y1），（X2，Y2）， …，（Xn，Yn）。對於平面中的這n個點，可以使用無數條曲線來擬合。要求樣本回歸函數盡可能好地擬合這組值。綜合起來看，這條直線處於樣本數據的中心位置最合理。 選擇最佳擬合曲線的標准可以確定為：使總的擬合誤差（即總殘差）達到最小。有以下三個標准可以選擇：

        （1）用“殘差和最小”確定直線位置是一個途徑。但很快發現計算“殘差和”存在相互抵消的問題。
        （2）用“殘差絕對值和最小”確定直線位置也是一個途徑。但絕對值的計算比較麻煩。
        （3）最小二乘法的原則是以“殘差平方和最小”確定直線位置。用最小二乘法除了計算比較方便外，得到的估計量還具有優良特性。這種方法對異常值非常敏感。

　 最常用的是普通最小二乘法（ Ordinary  Least Square，OLS）：所選擇的回歸模型應該使所有觀察值的殘差平方和達到最小。（Q為殘差平方和）- 即采用平方損失函數。

 　樣本回歸模型：

                                     其中ei為樣本（Xi, Yi）的誤差

   平方損失函數：

                      

   則通過Q最小確定這條直線，即確定，以為變量，把它們看作是Q的函數，就變成了一個求極值的問題，可以通過求導數得到。求Q對兩個待估參數的偏導數：

                       

    根據數學知識我們知道，函數的極值點為偏導為0的點。

    解得：

                   

 

這就是最小二乘法的解法，就是求得平方損失函數的極值點。

    LeastSquare(const vector<double>& x, const vector<double>& y)  
    {  
        double t1=0, t2=0, t3=0, t4=0;  
        for(int i=0; i<x.size(); ++i)  
        {  
            t1 += x[i]*x[i];  
            t2 += x[i];  
            t3 += x[i]*y[i];  
            t4 += y[i];  
        }  
        a = (t3*x.size() - t2*t4) / (t1*x.size() - t2*t2);  // 求得β1   
        b = (t1*t4 - t2*t3) / (t1*x.size() - t2*t2);        // 求得β2  
    }  

最小二乘法 
    設經驗 
方程是y=F(x)，方程中含有一些待定系數an，給出真實值{(xi,yi)|i=1,2,...n},將這些x,y值代入方程然后作 
差，可以描述誤差：yi-F(xi)，為了考慮整體的誤差，可以取平方和，之所以要平方是考慮到誤差可正可負直接相加可以相互抵消，所以記誤差為：

e=∑(yi-F(xi))^2

    它是一個多元函數，有an共n個未知量，現在要求的是最小值。所以必然滿足對各變量的偏導等於0，於是得到n個方程：

de/da1=0 
de/da2=0 
... 
de/dan=0

n個方程確定n個未知量為常量是理論上可以解出來的。用這種誤差分析的方法進行回歸方程的方法就是最小二乘法。

線性回歸 
如果經驗方程是線性的，形如y=ax+b，就是線性回歸。按上面的分析，誤差函數為：

e=∑(yi-axi-b)^2

各偏導為：

de/da=2∑(yi-axi-b)xi=0 
de/db=-2∑(yi-axi-b)=0

於是得到關於a,b的線性方程組：

(∑xi^2)a+(∑xi)b=∑yixi 
(∑xi)a+nb=∑yi

設A=∑xi^2,B=∑xi,C=∑yixi,D=∑yi，則方程化為：

Aa+Bb=C 
Ba+nb=D

解出a,b得：

a=(Cn-BD)/(An-BB) 
b=(AD-CB)/(An-BB) 

#include  <stdlib.h> 
#include  <iostream> 
#include  <valarray> 
using namespace std; 
int main(int argc, char *argv[]) 
{ 
    int num = 0; 
    cout << " Input how many numbers you want to calculate:"; 
    cin >> num; 
    valarray<double> data_x(num); 
    valarray<double> data_y(num); 
    while( num ) 
    { 
        cout << "Input the "<< num <<" of x:"; 
        cin >> data_x[num-1]; 
        cout << "Input the "<< num <<" of y:"; 
        cin >> data_y[num-1]; 
        num--; 
    } 
    double A =0.0; 
    double B =0.0; 
    double C =0.0; 
    double D =0.0; 
    A = (data_x*data_x).sum(); 
    B = data_x.sum(); 
    C = (data_x*data_y).sum(); 
    D = data_y.sum(); 
    double k,b,tmp =0; 
    if(tmp=(A*data_x.size()-B*B)) 
    { 
        k = (C*data_x.size()-B*D)/tmp; 
        b = (A*D-C*B)/tmp; 
    } 
    else 
    { 
        k=1; 
        b=0; 
    } 
    cout <<"k="<<k<<endl; 
    cout <<"b="<<b<<endl; 
    return 0; 
}