最小二乘法拟合直线

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

理论部分引用自：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);
    }
}