最小二乘法拟合直线c++代码

最近公司的一个项目需要计算TVDI（Temperature Vegetation Dryness Index ，温度植被干旱指数） ，TVDI的计算公式如下（具体原理自行百度）：



其中，

为任意像元的地表温度；

为某一NDVI对应的最小地表温度，对应的是湿边

；

为某一NDVI对应的最大地表温度，对应的是干边

；a，b为湿边的拟合方程系数，c，d为干边的拟合方程系数。

在拟合干边和湿边的过程中，需要利用最小二乘方法来对有效的NDVI和Lst数据来进行线性拟合。因此，本文记录在工作中用C++实现的最小二乘拟合直线，关键是理解最小二乘拟合直线的基本原理，实现起来比较简单。具体的最小二乘原理再此不做过多的阐述，网上有大量的介绍资料，这里只给出形如

的线性回归计算a，b系数以及r^2的最终计算公式，相关代码如下：

[cpp] view
plain copy

/************************************************************************* 

 最小二乘法拟合直线，y = a*x + b; n组数据; r-相关系数[-1,1],fabs(r)->1,说明x,y之间线性关系好，fabs(r)->0，x,y之间无线性关系，拟合无意义 

 a = (n*C - B*D) / (n*A - B*B) 

 b = (A*D - B*C) / (n*A - B*B) 

 r = E / F 

 其中： 

 A = sum(Xi * Xi) 

 B = sum(Xi) 

 C = sum(Xi * Yi) 

 D = sum(Yi) 

 E = sum((Xi - Xmean)*(Yi - Ymean)) 

 F = sqrt(sum((Xi - Xmean)*(Xi - Xmean))) * sqrt(sum((Yi - Ymean)*(Yi - Ymean))) 

 

**************************************************************************/  

void LineFitLeastSquares(float *data_x, float *data_y, int data_n, vector<float> &vResult)  

{  

    float A = 0.0;  

    float B = 0.0;  

    float C = 0.0;  

    float D = 0.0;  

    float E = 0.0;  

    float F = 0.0;  

  

    for (int i=0; i<data_n; i++)  

    {  

        A += data_x[i] * data_x[i];  

        B += data_x[i];  

        C += data_x[i] * data_y[i];  

        D += data_y[i];  

    }  

  

    // 计算斜率a和截距b  

    float a, b, temp = 0;  

    if( temp = (data_n*A - B*B) )// 判断分母不为0  

    {  

        a = (data_n*C - B*D) / temp;  

        b = (A*D - B*C) / temp;  

    }  

    else  

    {  

        a = 1;  

        b = 0;  

    }  

  

    // 计算相关系数r  

    float Xmean, Ymean;  

    Xmean = B / data_n;  

    Ymean = D / data_n;  

  

    float tempSumXX = 0.0, tempSumYY = 0.0;  

    for (int i=0; i<data_n; i++)  

    {  

        tempSumXX += (data_x[i] - Xmean) * (data_x[i] - Xmean);  

        tempSumYY += (data_y[i] - Ymean) * (data_y[i] - Ymean);  

        E += (data_x[i] - Xmean) * (data_y[i] - Ymean);  

    }  

    F = sqrt(tempSumXX) * sqrt(tempSumYY);  

  

    float r;  

    r = E / F;  

  

    vResult.push_back(a);  

    vResult.push_back(b);  

    vResult.push_back(r*r);  

}  

为了验证该算法的有效性，给出如下测试数据，数据来源为某论文的实验数据：

[cpp] view
plain copy

float pY[25] = { 10.98, 11.13, 12.51, 8.40, 9.27,  

                      8.73, 6.36, 8.50, 7.82, 9.14,  

                      8.24, 12.19, 11.88, 9.57, 10.94,  

                      9.58, 10.09, 8.11, 6.83, 8.88,  

                      7.68, 8.47, 8.86, 10.38, 11.08 };  

  

float pX[25] = { 35.3, 29.7, 30.8, 58.8, 61.4,  

                      71.3, 74.4, 76.6, 70.7, 57.5,  

                      46.4, 28.9, 28.1, 39.1, 46.8,  

                      48.5, 59.3, 70.0, 70.0, 74.5,  

                      72.1, 58.1, 44.6, 33.4, 28.6 };  

该数据在Excel的拟合结果为

，其中

。

转载地址   http://blog.csdn.net/pl20140910/article/details/51926886