最小二乘法拟合直线

为了求取两组数据之间的关联时，经常会采用曲线拟合的方式去求得两组数据之间的关系表达模型，曲线拟合中最基本和最常用的是直线拟合。

理论部分引用自：http://blog.sina.com.cn/s/blog_648868460100hevs.html

设x和y之间的函数关系为：

                   y＝a+bx

式中有两个待定参数，a代表截距，b代表斜率。对于等精度测量所得到的N组数据（x_i，y_i），i＝1，2……，N，x_i值被认为是准确的，所有的误差只联系着y_i。下面利用最小二乘法把观测数据拟合为直线。     

用最小二乘法估计参数时，要求观测值y_i的偏差的加权平方和为最小。对于等精度观测值的直线拟合来说，可使下式的值最小：

上式分别对a、b求偏导得：
    

整理后得到方程组

    解上述方程组便可求得直线参数a和b的最佳估计值。

相关系数r:

最小二乘法处理数据除给出a、b外，常常还给出相关系数r,  r定义为

 

C++实现代码如下：

#ifndef LSM_H#define LSM_H#include <iostream>#include <vector>using namespace std;class LSM{public:    LSM();    void initLSM(vector<int> src, vector<int> dest);       //获取样本点数    double geta();                       //获取斜率    double getb();                       //获取截距    double getr();                       //获取相关系数    ~LSM();protected:    void readCalcData(vector<int> src, vector<int> dest);        //读取数据并计算相关系数    void getTiTba(vector<int> data, double avg,double& sumTitba,double& sumTiba2);    double calAvg(vector<int> ptr);                //计算数组的平均值    int calSum(vector<int> ptr);        //计算整数和    int calXYsum(vector<int> ptrx,vector<int> ptry);    int calX2sum(vector<int> ptr);private:    int PointNum;      //样本点的总数    int sumX;       //所有样本点X的累加和    int sumY;       //所有样本点Y的累加和    int sumXY;      //所有样本点XY的累加和    int sumX2;      //X^2的和    double xAvg;    //x的平均值    double yAvg;    //y的平均值    double sumXiXba;    //计算相关系数用，（xi-x）或者(yi-y)    double sumXiXba2;   //计算相关系数用，（xi-x）^2或者(yi-y)^2    double sumYiYba;    double sumYiYba2;    double a;          //斜率    double b;          //截距    double r;          //相关系数};#endif // LSM_H

</pre><pre code_snippet_id="1639905" snippet_file_name="blog_20160409_1_2240911" name="code" class="cpp">

#include "lsm.h"#include <math.h>LSM::LSM(){    PointNum =0;    sumX = 0;    sumX2 = 0;    sumXY = 0;    sumY = 0;    xAvg = 0.0;    yAvg = 0.0;    sumXiXba = 0.0;    sumXiXba2 = 0.0;    sumYiYba = 0.0;    sumYiYba2 =0.0;    a =0.0;    b=0.0;    r=0.0;}//获取样本数据,读取指定文件的数据//必须采用此函数初始化LSM类void LSM::initLSM(vector<int> src, vector<int> dest){    readCalcData(src,dest);}double LSM::calAvg(vector<int> ptr){    double avg=0;    int Num = ptr.size();    for(int i=0;i<Num;i++)        avg += ptr[i];    avg = avg/Num;    return avg;}int LSM::calSum(vector<int> ptr){    int sum = 0;    int Num  = ptr.size();    for(int i=0;i<Num;i++)    {        sum += ptr[i];    }    return sum;}int LSM::calXYsum(vector<int> ptrx, vector<int> ptry){    int sum=0;    int Num = ptrx.size();    for(int i=0;i<Num;i++)    {        sum += ptrx[i]*ptry[i];    }    return sum;}int LSM::calX2sum(vector<int> ptr){    int sum=0;    int Num = ptr.size();    for(int i=0;i<Num;i++)    {        sum += ptr[i]*ptr[i];    }    return sum;}double LSM::geta(){    if(PointNum>0)        a = ((double)(PointNum*sumXY-sumY*sumX))/((double)(PointNum*sumX2-sumX*sumX));    return a;}double LSM::getb(){    if(PointNum>0)        b = (sumX2*sumY-sumX*sumXY)/((double)(PointNum*sumX2-sumX*sumX));    return b;}double LSM::getr(){    if(PointNum>0)        r = (sumXiXba*sumYiYba)/(sqrt(sumXiXba2)*sqrt(sumYiYba2));    return r;}LSM::~LSM(){}void LSM::readCalcData(vector<int> src,vector<int>dest){    PointNum = src.size();    if(PointNum != dest.size())    {        cout<<"uncompared vectors!\n"<<endl;        exit(-1);    }    sumX = calSum(src);   //计算x的和    sumY = calSum(dest);    sumX2 = calX2sum(src);    sumXY = calXYsum(src,dest);    xAvg = calAvg(src);    yAvg = calAvg(dest);    getTiTba(src,xAvg,sumXiXba,sumXiXba2);    getTiTba(dest,yAvg,sumYiYba,sumYiYba2);}void LSM::getTiTba(vector<int> data, double avg,double& sumTitba,double& sumTiba2){    sumTitba =0.0;    sumTiba2 = 0.0;    for(int i=0;i<data.size();i++)    {        sumTitba += data[i]-avg;        sumTiba2 += (data[i]-avg)*(data[i]-avg);    }}