/*
代码作者：不详
代码整理者：设计天下 MySDN网站 算法天下工作室
功能：多元线性回归分析
*/
#include "math.h"
#include "stdio.h"
#include "stdlib.h"
typedef struct _rmatrix {
int row; /* 行数 */
int col; /* 列数 */
double *data; /* 数据区 */
} RM, *RMP; /* RM: 实矩阵类型，RMP: 实矩阵类型指针*/
typedef struct _cnumber
{
double rpart;
double ipart;
} CNUM, *CNUMP; /* CNUM:复数类型, CNUMP 复数类型指针*/
typedef struct _cmatrix {
int row; /* 行数 */
int col; /* 列数 */
CNUMP data; /* 数据区 */
} CM, *CMP; /* CM: 复矩阵类型，CMP: 复矩阵类型指针*/
typedef struct _cmatrix2 {
int row; /* 行数 */
int col; /* 列数 */
double *rdata; /* 实部数据区 */
double *idata; /* 虚部数据区 */
} CM2, *CMP2; /* CM2: 复矩阵类型，CMP2: 复矩阵类型指针*/
typedef struct _tmatrix {
int row; /* 阶数 */
double *tdata; /* 数据区: T型阵的元素 t0, t1, ... tn */
double *ttdata; /* 数据区: 0, T型阵的元素 tt1, ... ttn */
} TM, *TMP; /* TM: 托伯利兹矩阵类型，TMP: 托伯利兹矩阵类型指针*/
typedef struct _trimatrix {
int row; /* 阶数 */
double *bdata; /* 数据区: 对称三角阵的主对角线元素 b0, b1, ... bn */
double *cdata; /* 数据区: 对称三角阵的主对角线元素 c0, c1, ... cn-1 */
} TRIM, *TRIMP; /* TRIM: 对称三角阵阵类型，TRIMP: 对称三角阵阵类型指针*/
typedef struct _tridiagonal {
int row; /* 阶数 */
double *data; /* 数据区: 三对角线元素 */
} TDM, *TDMP; /* TDM: 三对角线矩阵类型, TDMP: 三对角线矩阵类型指针 */
int equations_square_root_lis(RMP ap, RMP bp)
{
int i,j,m,n,k,u,v;
double *a,*d;
n=ap->row;
m=bp->col;
a=ap->data;
d=bp->data;
if (a[0]<0.0000001)
{
printf("fail/n");
return -1;
}
a[0]=sqrt(a[0]);
for (j=1; j<n; j++)
{
a[j]/=a[0];
}
for (i=1; i<n; i++)
{
u=i*n+i;
for (j=1; j<=i; j++)
{
v=(j-1)*n+i;
a[u]-=a[v]*a[v];
}
if (a[u]<0.0000001)
{
printf("fail/n");
return -2;
}
a[u]=sqrt(a[u]);
if (i!=(n-1))
{
for (j=i+1; j<n; j++)
{
v=i*n+j;
for (k=1; k<=i; k++)
{
a[v]-=a[(k-1)*n+i]*a[(k-1)*n+j];
}
a[v]/=a[u];
}
}
}
for (j=0; j<m; j++)
{
d[j]/=a[0];
for (i=1; i<n; i++)
{
u=i*n+i;
v=i*m+j;
for (k=1; k<=i; k++)
{
d[v]-=a[(k-1)*n+i]*d[(k-1)*m+j];
}
d[v]/=a[u];
}
}
for (j=0; j<m; j++)
{
u=(n-1)*m+j;
d[u]/=a[n*n-1];
for (k=n-1; k>0; k--)
{
u=(k-1)*m+j;
for (i=k; i<n; i++)
{
v=(k-1)*n+i;
d[u]-=a[v]*d[i*m+j];
}
v=(k-1)*n+k-1;
d[u]/=a[v];
}
}
return 0;
}
void mulliregression_lis(double x[],double y[],int m,int n,double a[],double dt[],double v[])
{
int i,j,k,m1;
double q,e,u,p,ay,s,r,pp,*b;
RM ma, mb;
m1=m+1;
b=malloc(m1*m1*sizeof(double));
b[m1*m1-1]=n;
for (j=0; j<m; j++)
{
for (i=0,p=0.0; i<n; i++)
{
p+=x[j*n+i];
}
b[m*m1+j]=b[j*m1+m]=p;
}
for (i=0; i<m; i++)
{
for (j=i; j<m; j++)
{
for (k=0,p=0.0; k<=n-1; k++)
{
p+=x[i*n+k]*x[j*n+k];
}
b[j*m1+i]=b[i*m1+j]=p;
}
}
for (i=0,a[m]=0.0; i<n; i++)
{
a[m]+=y[i];
}
for (i=0; i<m; i++)
{
for (j=0,a[i]=0.0; j<n; j++)
{
a[i]+=x[i*n+j]*y[j];
}
}
ma.row=m1; ma.col=1; ma.data=a;
mb.row=mb.col=m1; mb.data=b;
equations_square_root_lis(&mb, &ma);
for (i=0,ay=0.0; i<n; i++)
{
ay+=y[i];
}
ay/=n;
for (i=0,q=e=u=0.0; i<n; i++)
{
for (j=0,p=a[m]; j<m; j++)
{
p+=a[j]*x[j*n+i];
}
q+=(y[i]-p)*(y[i]-p);
e+=(y[i]-ay)*(y[i]-ay);
u+=(ay-p)*(ay-p);
}
s=sqrt(q/n);
r=sqrt(1.0-q/e);
for (j=0; j<m; j++)
{
for (i=0,p=0.0; i<n; i++)
{
for (k=0,pp=a[m]; k<=m-1; k++)
{
if (k!=j)
{
pp+=a[k]*x[k*n+i];
}
}
p+=(y[i]-pp)*(y[i]-pp);
}
v[j]=sqrt(1.0-q/p);
}
dt[0]=q; /*偏差平方和*/
dt[1]=s; /*平均标准差*/
dt[2]=r; /*复相关系数*/
dt[3]=u; /*回归平方和*/
free(b); /*释放动态分配的内存*/
}
int main()
{
double a[4],v[3],dt[4];
double x[3][4]={ {2.2,2.3,2.2,2.0},{1.9,1.8,2.0,2.1},{3.1,3.1,3.0,2.9}};
double y[4]={8.1,8.0,8.1,8.0};
int i, m=3, n=4;
mulliregression_lis((double*)x,y,m,n,a,dt,v);
printf("/n");
for (i=0; i<=3; i++)
{
printf("a[%d]=%e ",i,a[i]);
}
printf("/n");
printf("q=%e s=%e r=%e",dt[0],dt[1],dt[2]);
printf("/n");
for (i=0; i<=2; i++)
{
printf("v[%d]=%e ",i,v[i]);
}
printf("/n");
printf("u=%e",dt[3]);
return 0;
}