加权最小二乘回归方法的程序实现范例

      加权最小二乘回归方法-程序范例

在一般的线性回归求相关系数时候，大都选择最小二乘回归分析方法来拟合。它的原理是：假设拟合方程为：，将实测值Yi与预测值（Yj=a0+a1X）的离差（Yi-Yj）的平方和最小作为“优化判据”，那么对a0、a1求偏导数，令这两个偏导数等于零：

a0 = （∑Yi) / n - a1（∑Xi) / n，

a1 = [n∑(Xi Yi) –(∑Xi∑Yi)] / [n∑(Xi^2) – (∑Xi)^2 )]

 

或者用最小二乘公式推理如下：

∑（X-X_平）(Y-Y_平）=∑（XY-X_平Y-XY_平+X_平Y_平）

=∑XY-X_平∑Y-Y_平∑X+nX_平Y_平=∑XY-nX_平Y_平-nX_平Y_平+nX_平Y_平=∑XY-nX_平Y_平；

 

∑（X -X_平）²=∑（X²-2XX_平+X_平²)=∑X²-2X_平∑X+nX_平²=∑X²-2nX_平²+nX_平²

=∑X²-nX_平²；

(1) K=∑（X-X_平）(Y-Y_平）/∑（X -X_平）²⁼(∑XY-nX_平Y_平) /(∑X²-nX_平²)

(2) b=Y_平-kX_平

 

1）原理解决后，就是代码了，参考一个气候模型中回归函数的代码，主要目的是解决下一个问题：加权回归。

float regression(float  z[],float  temp[],intndata)

{

 float  sz, st;

 float   sz_media, st_media, sz_acum,st_minus_st_media, sz_minus_sz_media;

 int i;

 

float beta0=0.0;

 floatbeta1=0.0;

 

 sz=0.0;

 st=0.0;

 

 for(i=0;i<ndata;i++)

   {

     sz += z[i];

     st += temp[i];

   }

 

 /*sz_media =sz/ndata;*/

 sz_media=sz;

 st_media=st;

 //st_media =st/ndata;

 sz_acum = 0.0;

//printf("\n sz = %f,  st = %f,sz_media = %f, st_media = %f", sz, st, sz_media, st_media );

 for(i=0;i<ndata;i++)

  {

     sz_minus_sz_media = z[i]- sz_media;

     sz_acum += (sz_minus_sz_media *sz_minus_sz_media) ;

     st_minus_st_media = temp[i]- st_media;

     beta1 += (sz_minus_sz_media *st_minus_st_media);

   }

 

 beta1 /= ( sz_acum);

 

 beta0 = st_media - beta1 * sz_media;

 

//printf("\n beta0 = %f,  beta1 =%f", *beta0, *beta1);

 

 return beta1;

}

 

2）在需要考虑到不同数据对线性回归影响大小，比如分析浓度与检测峰值的回归关系的时候，浓度越大的数据样本对建模回归结果影响较大；在求解温度递减率的时候，温差大的数据站点的影响力度要大。所以简单最小二乘回归要考虑单个样本的权重大小。如何在程序中加入权重w数组。仔细对比下面代码和上面的，就可以了解如何处理。PS：权重大小用小数表示，权重总和是1。

float regression(floatz[],float temp[],float weight[],int ndata)

{

    float  sz, st;floatsw;float beta1=0.0;floatbeta0=0.0;

    float   sz_media, st_media, sz_acum,st_minus_st_media, sz_minus_sz_media;

    int i;

 

 

    sz=0.0;

    st=0.0;

 

    for(i=0;i<ndata;i++)

    {

        sz += weight[i]*z[i];

        st += weight[i]*temp[i];

    }

 

    //sz_media =sz/ndata;

    sz_media = sz;

    //st_media =st/ndata;

    st_media = st;

    sz_acum = 0.0;

    // printf("\nsz = %f,  st = %f, sz_media = %f,st_media = %f", sz, st, sz_media, st_media );

    for(i=0;i<ndata;i++)

    {

        sz_minus_sz_media = z[i]- sz_media;

        sz_acum += weight[i]*(sz_minus_sz_media* sz_minus_sz_media) ;

        st_minus_st_media = temp[i]- st_media;

        beta1 += weight[i]*(sz_minus_sz_media *st_minus_st_media);

    }

 

    beta1 /= sz_acum;

 

    beta0 = st_media - beta1 * sz_media;

 

    // printf("\nbeta0 = %f,  beta1 = %f", *beta0,*beta1);

 

    returnbeta1;

 

 

}