C语言版的线性回归分析函数

C语言版的线性回归分析函数

分类：C/C++2007-08-0323:3913840人阅读评论(31)收藏举报

语言c数学计算delphisystem

前几天，清理出一些十年以前DOS下的程序与代码•看来目前也没什么用了，想打个包刻在光碟

上，却发现有些代码现在可能还能起作用，其中就有计算一元回归和多元回归的代码，一看代码文件时间，

居然是1993年的，于是稍作整理，存放在这，分析虽不十分完整，但一般应用是没问题的，最起码■可提

供给那些刚学C的学生们参考。

先看看一元线性回归函数代码：

//求线性回归方程：Y=a+bx

//dada[rows*2]数组：X,Y；rows:数据行数；a,b：返回回归系数

//SquarePoor[4]：返回方差分析指标：回归平方和，剩余平方和，回归平方差，

剌余平方差

//返回值：0求解成功，-1错误

intLinearRegression(doub1e*data,introws,double/a,double沏do

uble*SquarePoor)

{

intni;

double*p,Lxx=0.0,Lxy=0.0,xa=0.0,ya=0.0;

if(data==0Ia==011b==0IIrows<1)

return-1;

for(p=data,n=0;m<rows;m++)

xa+=*p++;

ya+=电++:

}

xa/=rows;//X平均值

ya/=rows;//Y平均值

for(p=data,in=0;m<rows;m++,p+=2)

{

Lxx+=((*p-xa)*(*p-xa));//Lxx=Sum((X-

Xa)平方)

Lxy+=((*p-xa)*(*(p+1)-ya));//Lxy=Sum((X-

Xa)(Y-Ya))

)

米b=Lxy/Lxx;//b=Lxy/Lxx

=ya-*b*xa;//a=Ya-bxXa

if(SquarePoor==0)

return0;

//方差分析

SquarePoorfO]=SquarePoor[1]=0.0;

for(p=data,m=0;m<rows;m++,p++)

(

Lxy=*a+xb**p++;

SquarePoor[0]+=((Lxy-ya)*(Lxy-ya));//U(回归平方和)

SquarePoor[l]+=((配-Lxy)*(*p-Lxy));//Q(剩余平方和)

)

SquarePoor[2]=SquarePoor[0];//回归方差

SquarePoor[3]=SquarePoor[1]/(rows-2);//剩余方差

return0;

)

为了理解代码，把几个与代码有关的公式写在下面（回归理论和公式推导就免了，网上搜索到欠是，

下面的公式图E也是网上拽的，有曲公式图形网上没找到或者不公适，可参见后面多元回扫中的公式）：

1、回归方程式：

2'回均系数：

其中：

3'回归平方和：

4'轲余平方和：

实例计算：

doubledatal[121[2]={

//XY

{187.1,25.4},

{179.5,22.8},

(157.0,20.6),

{197.0,21.8},

{239.4,32.4},

{217.8,24.4},

{227.1,29.3),

(233.4,27.9(,

{242.0,27.8),

{251.9,34.2),

{230.0,29.2),

{271.8,30.0}

voidDisplay(double*dat,double*Answer,double米SquarePoor,introw

s,intcols)

(

doublev,*p;

inti,j;

printfC回归方程式：Y=Answer[0]);

for(i=1;i<cols;i++)

printf("+%.51f*X%d",Answer[i],i);

printfCn);

printf("回归显著性检验：")；

printf("回归平方和：％12.41f回归方差：%12.41f

",SquarePoor[OLSquarePoor[2]);

printfC剩余平方和：％12.41f剩余方差：%12.41f

",SquarePoor[1],SquarePoor[3]);

printfC离差平方和：%12.41f标准误差：%12.41f

",SquarePoorLOJ+SquarePoorL1J,sqrt(SquarePoor[3J));

printfC'F检验：％12.41f相关系数：％12.41f

",SquarePoor[2]/SquarePoor[3],

sqrt(SquarePoor[0]/(SquarePoor[0]+SquarePoor[l])));

printfC剩余分析：");

printfC观察值估计值剩余值剩余平方");

for(i=0,p=dat;i<rows;i++,p++)

(

v=Answer[0];

for(j=1;j<cols;j++,p++)

v+=*p*Answer[j];

printfC%12.21f%12.21f%12.21f%12.21f

",*p,v,*p-V,(xp-v)*(迎-v));

)

systemC'pause");

)

intmain()

doubleAnswer[2],SquarePoor[4];

if(LinearRegression((double*)data1,12,&Answer[0],&Answer[1],

SquarePoor)==0)

Disp1ay((doub1e*)data1,Answer,SquarePoor,12,2);

return0;

1

运行结果：

3.41297*0.10814MX1

日至性检验;

方差

141.6641141.6641

54.6059差5.4606

¥和：196.2700就2.3368

检25.94300.8496

观察值估计值剩余值剩余平方

25.4023.651.753.08

22.8022.82-0.020.00

20.6020.390.210.04

21.8024.72-2.928.51

32.4029.303.109.60

24.4026.97-2.576.59

29.3027.971.331.76

27.9028.65-0.750.57

27.8029.58-1.783.18

34.2030.653.5512.58

29.2028.290.910.84

30.0032.81-2.817.87

上面的函数和例子程序X仅计算了回归方程式，还计算了显著性检验指标，例如F检脸指标，我们可

以在统计F分布表上查到F0.01(1,10)=10.04(注：括号里的1,10分别为回归平方和和轲余平方和所拥有

的自由度)，小于计算的F检验值25.94，可以认为该回归例子高度显著。

如果使用图形界面，可以根据原始数据和计算结果绘制各种图表，如散点图、趋势图、控制图等。很

多非线性方程可以借助数学计算，转化为直线方程进行回归分析。

同一元线性回归相比，多元线性回归分析代码可就复杂多了，必须求解线性方程，因此本代码中包含

一个可独立使用的线性方程求解函数：

voidFreeData(double米*dat,double*d,intcount)

(

inti,j;

free(d);

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

free(dat[i])；

free(dat);

)

//解线性方程。data[count米(count+1)]矩阵数组；count:方程元数；

//Answer[count]­求解数组。返回：0求解成功1-1无解或者无穷解

intLinearEquations(doub1e*data,intcount,double*Answer)

(

intj,m,n;

doubletmp,糕dat,*d=data;

dat=(double**)malloc(count*sizeof(double*));

for(m=0;m<count;m++,d+=(count+1))

(

dat[m]=(doub1e*)ma11oc((count+1)*sizeof(double));

memcpy(dat[n],d,(count+1)*sizeof(double));

d=(double*)malloc((count+1)*sizeof(double));

for(m=0;m<count-1;m++)

(

//如果主对由线元素为0，行交换

for(n=m-1;n<count&&dat[m][m]==0.0;n++)

(

if(dat[n][m]!=0.0)

{

memcpy(d,dat[m],(count+1)*sizeof(double));

memcpy(dat[in],dat[n],(count+1)米sizeof(double));

memcpy(dat[n],d,(count+1)*sizeof(double));

)

)

//行交换后，主对有线元素仍然为0，无解，返回-1

if(dat[m][m]==0.0)

(

FreeData(dat,d,count);

return-1;

)

//消元

for(n=m-1;n<count;n++)

tmp=dat[n][m]/dat[m][m];

for(j=m;j<=count;j++)

dat[n][j]-=tmp*dat[m][j];

)

)

for(j=0;j<count;j++)

d[j]=0.0;

//求得count-1的元

Answer[count-1]=dat[count-1Hcount]/dat[count-1][count

-1]；

//逐行代入求各元

for(m=count-2;m>=0;m­)

{

for(j=count-1;j>m;j-)

d[m]+=Answer[j]*dat[m][j];

Answer[m]=(dat[m][count]-d[m])/dat[m][m];

)

FreeData(dat,d,count);

return0;

//求多元回归方程：Y=BO+B1X1+B2X2+...BnXn

//data[rows*cols]二维数组；XIi,X2i,...Xni,Yi(i=0torows-1)

//rows:数据行数；cols数据列数；Answer[cols]:返回回归系数数组

(B0,BL..Bn)

//SquarePoor[4」：返回方差分析指标：回归平方和，剩余平方和，回归平方差，

剩余平方差

//返回值：0求解成功，-1错误

intMu11ip1eRegression(doub1e*data,introws,intcols,double*Answ

er,double*SquarePoor)

(

intm,n,i,count=cols-1;

doubleMat,迎，a,b;

if(data==0IAnswer==0IIrows<2IIcols<2)

return-1;

dat=(double*)nalloc(cols*(cols+1)*sizeof(double));

dat[0]=(double)rows;

for(n=0;n<count;n++)//n=0tocols

-2

(

a=b=0.0;

for(p=data+n,m=0;m<rows;m++,p+=cols)

a+=*p：

b+=(的)**p)；

(

dat[n+1]；a;//dat[0,n+1]=

Sum(Xn)

dat[(n+1)*(cols+1)]=a;//dat[n+L0]=

Sum(Xn)

dat[(n+1)*(cols+1)+n+1]=b;//dat[n+l,n-1]

=Sum(Xn*Xn)

for(i=n1;i<count;i++)//i=n+1toco

Is-2

(

for(a:0.0,p=data,m=0;m<rows;m++,p+=cols)

a+：(p[n]*p[i]);

dat[(n1)*(cols+1)+i+1]=a;//dat[n+l,i+1]

=Sum(Xn*Xi)

dat[(i1)*(cols+1)+n+1]=a;//dat[i+1,n+1]

=Sum(Xn*Xi)

)

)

for(b=0.0,m=0,p=data+n;m<rows;m++,p+=cols)

b+=*p;

dat[cols]=b;//dat[0,cols]

=Sum(Y)

for(n=0;n<count;n++)

(

for(a=0.0,p=data,m=0;m<rows;m++,p+=cols)

a+=(p[n]*p[count]);

dat[(n+1)*(cols+1)+cols]=a;//dat[n+l,cols

]=Sum(Xn*Y)

)

n=LinearEquations(dat,cols,Answer)://计算方程式

//方差分析

if(n==0&&SquarePoor)

{

b=b/rows;//b=Y的平均

值

SquarePoor[0]=SquarePoor[1]=0.0;

p=data;

for(m=0;m<rows;m++,p++)

(

for(i=1,a=Answer[0];i<cols;i++,P++)

a+=(*p*Answer[i])；//a=Ym的估计

值

SquarePoor[0J+=((a-b)*(a-b));//U(回归平方和)

SquarePoor[1]+=((*p-a)*(xp-a));//Q(剩余平方

和)(*p=Ym)

)

SquarePoorL2]=SquarePoorLOJ/count;//回归方差

if(rows-cols>0.0)

SquarePoor[3]=SquarePoor[l]/(rows-cols);//剩余方差

else

SquarePoor[3]=0.0;

)

free(dat);

returnn;

)

为了理解代码，同样贴几个主要公式在下面，其中回归平方和和剩余平方和

公式和一元回归一样：

1、回归方程式：，

2、回归系数方程组：

Z)IA-J)u-l)U-1)

%卬…[£&}・■£3i

.*-1

3'F检验：

4、相关系数：，其中，Syy是离差平方和（回归平方和与剩余平方和之和）。

该公式其实就是U/（U+Q）的平方根（没找到这个公式的图）。

5、回归方差：U/m，m为回归方程式中自变量的个数（没找到图）。

6'剩余方差：Q/（n-m-1），n为观察数据的样本数，m同上（没找至U

图）。

7、标准误差：也叫标准误，就是剩余方差的平方根（没找到图）。

下面是一个多元回归的例子：

doubledata[15][5]={

//XIX2X3X4Y

{316,1536,874,981,3894},

(385,1771,777,1386,4628},

{299,1565,678,1672,4569},

(326,1970,785,1864,5340},

(441,1890,785,2143,5449),

{460,2050,709,2176,5599),

{470,1873,673,1769,5010),

{504,1955,793,2207,5694},

{348,2016,968,2251,5792),

{400,2199,944,2390,6126},

{496,1328,749,2287,5025},

{497,1920,952,2388,5924),

{533,1400,1452,2093,5657),

{506,1612,1587,2083,6019),

{458,1613,1485,2390,6141},

voidDisplay(double*dat,double*Answer,double米SquarePoor,introw

s,intcols)

(

doublev,*p;

inti,j;

printf("回归方程式：Y=%,51f",Answer[0]);

for(i=1;i<cols;i++)

printf("+%.51f*X%d",Answer[i],i);

printfCn);

printf("回归显著性检验:"

printf("回归平方和：%12.41f回归方差:%12.41f

,SquarePoor[0],SquarePoor[2]);

printfC剩余平方和：％12.41f剩余方差：%12.41f

,SquarePoorf1],SquarePoor[3〕);

printfC离差平方和：％12.41f标准误差：%12.41f

,SquarePoor[0]+SquarePoor[l],sqrt(SquarePoor[3]));

printfC'F检验：％12.41f相关系数：％12.41f

,SquarePoor[2]/SquarePoor[3],

sqrt(SquarePoor[0J/(SquarePoor[0]+SquarePoor"])));

printfC剩余分析：”)；

printfC观察值估计值剩余值剩余平方");

for(i=0,p=dat;i<rows;i++,p++)

{

v=Answer[0];

for(j=1;j<cols;j++,p++)

v+=*Answer[j];

printf("%12.21f%12.21f%12.21f%12.21f

,电v,知-v,(xp-v)*(而-v));

systemCpause"):

intmain()

doub1eAnswer[5],SquarePoor[4];

if(Mu11ip1eRegression((doub1e*)data,15,5,Answer,SquarePoor)

==0)

Display((double*)data,Answer,SquarePoor,15,5);

return0;

}

运行结果见下图，同上面,查F分布表，F检验远远大于F0.005(4,10)的

7.34,可以说是极度回归显著。

CBSainpleMultipleRegressionMultiRegreTes1.exe

］回归方程式：

Y=448.64758+0.54412*X1+1.01215*X2+0.99077*X3+0.981??*X4

［回归显著性检验：

心归；：方和：5889270.0128回归方差：1L472317.5032

牌U余军方和：42501.7206R力麦：4250.1721

1离差"方和：5931771.7333标用圭误差：65.1933

F检验：346.4136相关系数：0.9964

剩余分析:

观察值估计值剩余值剩余平方

3894.004004.29-110.2912164.51

4628.004581.2046.802189.87

4569.004508.6160.393647.29

5340.005227.73112.2712604.62

5449.005483.25-34.251172.74

5599.005612.63-13.63185.71

5010.005003.676.3340.06

5694.005654.0739.931594.09

5792.005847.51-55.513081.78