当前位置:   article > 正文

最小二乘法原理

最小二乘法原理

转自:作者:金良(golden1314521@gmail.com) csdn博客:http://blog.csdn.net/u012176591

1.线性代数模型

首先给出最小二乘解的矩阵形式的公式:


推导过程:


条件:

矩阵必须是列满秩矩阵,否则的逆就不会存在。

若A为m×n的矩阵,b为m×1的矩阵,则Ax=b表达了一个线性方程组,它的normal equation的形式为ATAx=ATb。
当Ax=b有解时(即矩阵[A|b]的秩与A的秩相同),Ax=b与ATAx=ATb的解集是一样。
而当Ax=b无解时,ATAx=ATb仍然有解,其解集即最小二乘解(least squares solution),即使得(Ax-b)T(Ax-b)的值最小的解,可以理解为使方程组Ax=b近似成立且误差最小的解。

Python语言写的一个例子:

[python]  view plain copy 在CODE上查看代码片 派生到我的代码片
  1. #encoding=UTF-8  
  2. ''''' 
  3. Created on 2014年6月30日 
  4.  
  5. @author: jin 
  6. '''  
  7. from numpy import *  
  8. import matplotlib.pyplot as plt  
  9. from random import *  
  10.   
  11. def loadData():  
  12.     x = arange(-1,1,0.02)  
  13.     y = ((x*x-1)**3+1)*(cos(x*2)+0.6*sin(x*1.3))  
  14.     #生成的曲线上的各个点偏移一下,并放入到xa,ya中去  
  15.     xr=[];yr=[];i = 0  
  16.     for xx in x:  
  17.         yy=y[i]  
  18.         d=float(randint(80,120))/100  
  19.         i+=1  
  20.         xr.append(xx*d)  
  21.         yr.append(yy*d)    
  22.     return x,y,xr,yr  
  23. def XY(x,y,order):  
  24.     X=[]  
  25.     for i in range(order+1):  
  26.         X.append(x**i)  
  27.     X=mat(X).T  
  28.     Y=array(y).reshape((len(y),1))  
  29.     return X,Y  
  30. def figPlot(x1,y1,x2,y2):  
  31.     plt.plot(x1,y1,color='g',linestyle='-',marker='')  
  32.     plt.plot(x2,y2,color='m',linestyle='',marker='.')  
  33.     plt.show()  
  34. def Main():      
  35.     x,y,xr,yr = loadData()  
  36.     X,Y = XY(x,y,9)  
  37.     XT=X.transpose()#X的转置  
  38.     B=dot(dot(linalg.inv(dot(XT,X)),XT),Y)#套用最小二乘法公式  
  39.     myY=dot(X,B)  
  40.     figPlot(x,myY,xr,yr)  
  41. Main()  

程序截图:



MATLAB写的例子:

[plain]  view plain copy 在CODE上查看代码片 派生到我的代码片
  1. clear  
  2. clc  
  3. Y=[33815    33981   34004   34165   34212   34327   34344   34458   34498   34476   34483   34488   34513   34497   34511   34520   34507   34509   34521   34513   34515   34517   34519   34519   34521   34521   34523   34525   34525   34527]  
  4. T=[1    2   3   4   5   6   7   8   9   10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30]  
  5. % 线性化处理  
  6. for t = 1:30,   
  7.    x(t)=exp(-t);  
  8.    y(t)=1/Y(t);  
  9. end  
  10. % 计算,并输出回归系数B  
  11. c=zeros(30,1)+1;  
  12. X=[c,x'];  
  13. B=inv(X'*X)*X'*y'  
  14. for i=1:30,  
  15. % 计算回归拟合值      
  16.     z(i)=B(1,1)+B(2,1)*x(i);  
  17. end  
  18. Y2=[]  
  19. for j=1:30,  
  20.     Y2(j)=1/(B(1,1)+B(2,1)*exp(-j));  
  21. end  
  22. plot(T,Y2)  
  23. hold on  
  24. plot(T,Y,'r.')  

截图:


2、一般线性模型

一般线性模型,即普通最小二乘法( Ordinary  Least Square,OLS):所选择的回归模型应该使所有观察值的残差平方和达到最小。
以简单线性模型 y =  b1t +b0 作为例子。

回归模型:


最优化目标函数:




则目标函数可以简化成如下形式:


对简化后的目标函数进行求解,得到表达式:

下面是C++的实现例子:
[cpp]  view plain copy 在CODE上查看代码片 派生到我的代码片
  1. #include<iostream>  
  2. #include<fstream>  
  3. #include<vector>  
  4. #include<cstdlib>  
  5. #include<ctime>  
  6. using namespace std;  
  7.   
  8. class LeastSquare{  
  9.     double b0, b1;  
  10.     public:  
  11.         LeastSquare(const vector<double>& x, const vector<double>& y)  
  12.         {   //下面是最小二乘法的核心过程  
  13.             double xi_xi=0, xi=0, xi_yi=0, yi=0;  
  14.             for(int i=0; i<x.size(); ++i)  
  15.             {  
  16.                 xi_xi += x[i]*x[i];  
  17.                 xi += x[i];  
  18.                 xi_yi += x[i]*y[i];  
  19.                 yi += y[i];  
  20.             }  
  21.             b1= (xi_yi*x.size() - xi*yi) / (xi_xi*x.size() - xi*xi);    
  22.             b0 = (xi_xi*yi - xi*xi_yi) / (xi_xi*x.size() - xi*xi);         
  23.         }  
  24.       
  25.         double getY(const double x) const  
  26.         {  
  27.             return b0+b1*x;  
  28.         }  
  29.         void print() const  
  30.         {  
  31.             cout<<"y = "<<b0<<"+"<<b1<<"x  "<<"\n";  
  32.         }  
  33.     };  
  34.   
  35. int main()  
  36. {  
  37.         srand((unsigned int)(time(NULL)));//  
  38.         vector<double> x,y;  
  39.         double xi = 0,yi,xin,xout;  
  40.         for (int i = 0;i <10; i++) {       
  41.             yi = 2*xi +1;//原模型  
  42.             yi += 0.05*rand()/RAND_MAX*yi;//添加噪声  
  43.             y.push_back(yi);  
  44.             x.push_back(xi);  
  45.   
  46.             xi += 5.0;  
  47.         }  
  48.         LeastSquare lsObj(x, y);//用样本数据实例化对象  
  49.         lsObj.print();          //输出最小二乘法得到的模型  
  50.         cout<<"Input x:\n";  
  51.         while(cin>>xin)  
  52.         {  
  53.             xout = lsObj.getY(xin);//利用得到的模型计算因变量  
  54.             cout<<"y = "<<xout<<endl;  
  55.             cout<<"Input x:\n";  
  56.         }  
  57. }  

执行效果截图:

3.最小二乘法和梯度下降算法

相同点:
1 、本质相同:两种方法都是在给定已知数据(因变量 & 自变量)的前提下对因变量算出出一个一般性的估值函数。然后对给定的新的自变量用估值函数对其因变量进行估算。
2、 目标相同:都是在已知数据的框架内,使得估算值与实际值的差的平方和尽量更小(事实上未必一定要使用平方),估算值与实际值的差的平方和的公式为:


不同点:
1、实现方法和结果不同:
最小二乘法是直接对error求导找出全局最小,是非迭代法。
而梯度下降法是一种迭代法,有一个学习的过程,先由给定参数计算一个error,然后向该error下降最快的方向调整参数值,在若干次迭代之后找到局部最小。
梯度下降法的缺点是到最小点的时候收敛速度变慢,并且对初始点的选择极为敏感,其改进大多是在这两方面下功夫。
声明:本文内容由网友自发贡献,不代表【wpsshop博客】立场,版权归原作者所有,本站不承担相应法律责任。如您发现有侵权的内容,请联系我们。转载请注明出处:https://www.wpsshop.cn/w/IT小白/article/detail/117430
推荐阅读
相关标签
  

闽ICP备14008679号