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; }