【算法研究与实现】最小二乘法直线拟合
作者：gnuhpc
出处：http://www.cnblogs.com/gnuhpc/
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;
}
