|
|
楼主 |
发表于 2005-1-19 23:09:03
|
显示全部楼层
//多元线性回归分析
//x-长度为m*n的数组，每列存放m个自变量的一组观测值
//y-长度为n的数组，存放n个观测点y的值
//m-自变量的个数
//n-观测数据的组数
//a-长度为m+1的数组，返回回归系数a[0],a[1],a[2]...a[m]
//dt-长度为6的数组，
// dt[0]-返回偏差平方和 q=sum[(y-a x-b)^2]
// dt[1]-返回平均标准偏差 s=sqrt[q/n]
// dt[2]-返回复相关系数 r
// dt[3]-返回回归平方和 u
int chchol(double a[],int n,int m,double d[]);
void jcsqt(double x[],double y[],int m,int n,double a[],double dt[],double v[])
{
int i,j,k,l,mm;
double q,e,u,p,yy,s,r,pp,*b;
b=(double *)malloc((m+1)*(m+1)*sizeof(double));
mm=m+1;
b[mm*mm-1]=n;
for (j=0; j<=m-1; j++)
{
p=0.0;
for (i=0; i<=n-1; i++)
{
p=p+x[j*n+i];
}
b[m*mm+j]=p;
b[j*mm+m]=p;
}
for (i=0; i<=m-1; i++)
{
for (j=i; j<=m-1; j++)
{
p=0.0;
for (k=0; k<=n-1; k++)
{
p=p+x[i*n+k]*x[j*n+k];
}
b[j*mm+i]=p;
b[i*mm+j]=p;
}
}
a[m]=0.0;
for (i=0; i<=n-1; i++)
{
a[m]=a[m]+y;
}
for (i=0; i<=m-1; i++)
{
a=0.0;
for (j=0; j<=n-1; j++)
{
a=a+x[i*n+j]*y[j];
}
}
chchol(b,mm,1,a);
yy=0.0;
for (i=0; i<=n-1; i++)
{
yy=yy+y/n;
}
q=0.0; e=0.0; u=0.0;
for (i=0; i<=n-1; i++)
{
p=a[m];
for (j=0; j<=m-1; j++)
{
p=p+a[j]*x[j*n+i];
}
q=q+(y-p)*(y-p);
e=e+(y-yy)*(y-yy);
u=u+(yy-p)*(yy-p);
}
s=sqrt(q/n);
r=sqrt(1.0-q/e);
for (j=0; j<=m-1; j++)
{
p=0.0;
for (i=0; i<=n-1; i++)
{
pp=a[m];
for (k=0; k<=m-1; k++)
{
if (k!=j)
{
pp=pp+a[k]*x[k*n+i];
}
}
p=p+(y-pp)*(y-pp);
}
v[j]=sqrt(1.0-q/p);
}
dt[0]=q;
dt[1]=s;
dt[2]=r;
dt[3]=u;
free(b);
return;
}
int chchol(double a[],int n,int m,double d[])
{
int i,j,k,u,v;
if ((a[0]+1.0==1.0)||(a[0]<0.0))
{
printf("fail\n");
return(-2);
}
a[0]=sqrt(a[0]);
for (j=1; j<=n-1; j++)
{
a[j]=a[j]/a[0];
}
for (i=1; i<=n-1; i++)
{
u=i*n+i;
for (j=1; j<=i; j++)
{
v=(j-1)*n+i;
a=a-a[v]*a[v];
}
if ((a+1.0==1.0)||(a<0.0))
{
printf("fail\n");
return(-2);
}
a=sqrt(a);
if (i!=(n-1))
{
for (j=i+1; j<=n-1; j++)
{
v=i*n+j;
for (k=1; k<=i; k++)
{
a[v]=a[v]-a[(k-1)*n+i]*a[(k-1)*n+j];
}
a[v]=a[v]/a;
}
}
}
for (j=0; j<=m-1; j++)
{
d[j]=d[j]/a[0];
for (i=1; i<=n-1; i++)
{
u=i*n+i; v=i*m+j;
for (k=1; k<=i; k++)
{
d[v]=d[v]-a[(k-1)*n+i]*d[(k-1)*m+j];
}
d[v]=d[v]/a;
}
}
for (j=0; j<=m-1; j++)
{
u=(n-1)*m+j;
d=d/a[n*n-1];
for (k=n-1; k>=1; k--)
{
u=(k-1)*m+j;
for (i=k; i<=n-1; i++)
{
v=(k-1)*n+i;
d=d-a[v]*d[i*m+j];
}
v=(k-1)*n+k-1;
d=d/a[v];
}
}
return(2);
} |
|
发表于 2005-1-19 23:08:47
>选自<<徐世良数值计算程序集(C)>></P>