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++實現

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