最小二乘法直线拟合

点击打开链接

1.原理 
在现实中经常遇到这样的问题，一个函数并不是以某个数学表达式的形式给出，而是以一些自变量与因变量的对应表给出，老师讲课的时候举的个例子是犯罪人的身高和留下的脚印长，可以测出一些人的数据然后得到一张表，它反应的是一个函数，回归的意思就是将它还原成数学表达式，这个式子也称为经验表达式，之所以叫经验就是说它不完全是实际中的那样准确，是有一定偏差的，只是偏差很小罢了。

最小二乘法 
    设经验 
方程是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) 
这就是我们要进行的算法。 
2.C++实现 
/* 
* ===================================================================================== 
* 
*       Filename:  nihe.cpp 
* 
*    Description:  A least square method for fitting a curve 
* 
*        Version:  1.0 
*        Created:  03/21/2009 12:32:56 PM 
*       Revision:  none 
*       Compiler:  gcc 
* 
*         Author:  Futuredaemon (BUPT), gnuhpc@gmail.com 
*        Company:  BUPT_UNITED 
* 
* ===================================================================================== 
*/ 
#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; 
} 
3.OpenCV结构实现 
#include "cv.h" 
#include <iostream> 
using namespace std; 
int main(int argc, char *argv[]) 
{ 
  int i=0; 
  int j=0; 
  int num; 
  double A,B,C,D; 
  double k,b,tmp=0; 
  cout <<"Input how many numbers you want to calculate:"; 
  cin >>num; 
  CvMat *mat1=cvCreateMat(1,num,CV_64FC1); 
  CvMat *mat2=cvCreateMat(1,num,CV_64FC1); 
  CvMat *mattmp=cvCreateMat(1,num,CV_64FC1); 
  for (j=0;j<mat1->cols;j++) 
    { 
      cout << "data X"<<j<<"="; 
      cin>>CV_MAT_ELEM(*mat1,double,0,j); 
      cout << "data Y"<<j<<"="; 
      cin>>CV_MAT_ELEM(*mat2,double,0,j); 
    } 
  for (j=0;j<mat1->cols;j++) 
    { 
      cout<<"X="<<CV_MAT_ELEM(*mat1,double,0,j) 
          <<",Y="<<CV_MAT_ELEM(*mat2,double,0,j)<<endl; 
    } 
  cvMul(mat1,mat1,mattmp,1); 
  A = cvSum(mattmp).val[0]; 
  B = cvSum(mat1).val[0]; 
  cvMul(mat1,mat2,mattmp,1); 
  C = cvSum(mattmp).val[0]; 
  D = cvSum(mat2).val[0]; 
  tmp = A*mat1->cols-B*B; 
  k = (C*mat1->cols-B*D)/tmp; 
  b = (A*D-C*B)/tmp; 
  cout << "k=" << k <<endl; 
  cout << "b=" << b <<endl; 
  cvReleaseMat(&mat1); 
  cvReleaseMat(&mat2); 
  return 0; 
}