最小二乘法拟合直线公式推导及vc实现[转]

 

最小二乘拟合直线方程：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)
这就是我们要进行的算法。

 

-----------------------------------------------------------------------------------

1.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
#include
#include

using namespace std;

int main(int argc, char *argv[])
{
int num = 0;

cout << " Input how many numbers you want to calculate:";
cin >> num;

valarray data_x(num);
valarray 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="< cout <<"b="<

return 0;
}

 

2.OpenCV结构实现

#include "cv.h"
#include

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;jcols;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;jcols;j++)
{

cout<<"X="< <<",Y="< }

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 < cout << "b=" << b <

cvReleaseMat(&mat1);
cvReleaseMat(&mat2);

return 0;
}