线性最小二乘拟合算法实现-附C++源码

最近由于工作需要要做一套线性方程的拟合算法，在网上找了很多，不是用gsl就是用matlab做的，都不符合要求。实在找不到合适的于是就决定自己从头写，于是自己写了一个，用了一段时间决得还行，就决定贴出来和大家分享一下，一方面可以顺便骗点积分，另一方面也算是做个备忘吧。

源代码可以在百度云盘中下载：https://pan.baidu.com/s/1nvgjf41 密码：8kqs

如果觉得积分用不完也可以在http://download.csdn.net/download/lxwyw10902/10178153下载。

1、常用的一些函数的拟合算法，包括多项式、指数、对数和幂函数等
2、可以设置截距、权重
3、不依赖第三方库，跨平台运行

一、实例模型

         有一种传感器（假设就是温度传感器），传感器工作的时候会输出一个电压信号u（mV），我们需要根据这个电压信号预测出传感器感知到的温度t（℃）。假设通过大量的实验数据我们知道温度和电压信号呈如下的函数关系：

 

         每个传感器由于自身的特性导致各自的a、b、c值都不一样，因此在传感器出厂时我们需要对这个传感器进行标定，标定过程中，我们先在传感器的量程上均匀地挑选几个点，然后分别用被测传感器和标准仪器测量各个点的值，预测方程式中的a、b、c值。假设最终测得的这些离散的点为  ，我们的目标就是找到一条曲线，使得标准仪器测量的结果和通过换算后的被测传感器的差距的总和最小。这里我们令y=t，x=u，

我们令

我们的目标就是求出a,b,c的值使得g(a,b,c)取得最小值。但是上面的式子中含有绝对值，使得式子变成一个非线性的函数而无法求解，可以令

推广到更一般的形式，假设

最终目标就是计算

 

二、计算原理

         通过上面的转换，一个实际的问题就变成了一个求极值的问题。如下所示，这一类的函数是可以通过求每个分量的偏导数来计算极小值的，并且极小值就是最小值。证明过程就不写出来了，如果感兴趣可以自己去查（其实我也不知道怎么证明，哈哈）

对各个分量求偏导数

展开上面的等式，得到了一个m+1维的一次方程组，然后就需要求解这个一次方程组了

写成增广矩阵的形式，如下所示

观察这个矩阵发现左边的矩阵为一个对角矩阵，得到这个矩阵后就可以根据高斯消元的方法求解，这里就不详细介绍了，最终化简的结果就是

 

三、代码实现

有了上面的推算就可以编写计算多项式线性回归的算法代码了，C++编写的算法很方便，捣鼓了一阵子，就用C++写了一套算法，这些代码我分为三个部分，

1、高斯消元的部分，对应的文件是augmentedmatrix.h和augmentedmatrix.cpp；

2、计算多项式的部分，对应的文件是polynomialfit.h和polynomialfit.cpp；

3、基于多项式拟合的算法计算指数、对数和幂函数的拟合值，对应的文件是linearfit.h和linearfit.cpp。

除此之外还添加了

1、 权重的算法

2、 部分函数还可以设置截距

3、 拟合优度的计算方式

但是测试过程中，发现除了不设置截距的多项式拟合的优度和Excel计算的优度能够对应上之外，其他函数拟合的优度都无法和Excel对应上，也不知道Excel是如何计算拟合优度的。另外Excel也不能做权重的拟合，所以加了权重的拟合方式也没有测试，但是加权值的效果确实做到了，优度也可以体现拟合结果的好坏。

 

四、测试效果

intmain(int ,char**)

{

    PolyfitTest();

//   PowerTest();

//   ExponentTest();

//   LogarithmTest();

 

    return0;

}

 

// 多项式拟合测试

voidPolyfitTest()

{

    intdegree=2;

    intcount=5;

    doublesx[]={3.35,100.35,258.23,526.24,935.6};

    doublesy[]={0.5225,10200.8225,67000.1929,277000.0176,697000.56};

    doublew[]={1,0.02,0.01,0.005,0.001};

    Polynomialpf;

    pf.setAttribute(degree,false,1.0);

    if(pf.setSample(sx,sy,count,false,w)&&

            pf.process()){

        pf.print();

        std::cout<<"\ngoodness: "<<pf.getGoodNess()<<std::endl;

    }else

        std::cout<<"failed";

}

 

运行结果：

Polynomial Fit :

8.47470e+0119.82936e+008 1.22904e+006 6.91397e+011

9.82936e+0081.22904e+006 1823.77 8.16207e+008

1.22904e+0061823.77 5 1.0512e+006

results:0.58924    208.84    -10774.8

goodness: 0.998259

下面贴出Excel的计算结果

其他阶次的多项式也试过了，Excel最多好像只能拟合到6阶，这套代码设置最高到32阶，感觉也够了。

 

就写到这里，也不知道有没有讲清楚，由于时间仓促、水平有限，代码里面当然也难免会有不足之处，欢迎大家批评指正。