多元最小二乘法

多元最小二乘法[C语言版]2007年02月13日 星期二 13:19/*======关键词:最小二乘法,逆矩阵,矩阵,伴随矩阵,矩阵乘法,c语言,数组==*/
/*======春节在厦大过年了,将最小二乘法用C语言实现,也算一点收获.:)==*/
/*=====由于采用递归,这个算法实际测试不能高于9维,否则,近似死机=====*/
/*====改进的算法请参考Blog内其它内容:单位矩阵法和Gaues-Jordon法===*/
/*======Copyright@时光牧者( https://www.360docs.net/doc/5713889756.html,/junus ).==========*/
/*======Email: firt123@https://www.360docs.net/doc/5713889756.html, ==============================*/
/*=====================main.c=============================*/

#include <io.h>
#include <..\Array\OLS.c>
FILE *f;

main()
{
int i,j;

if ((f=fopen("array.txt","rt"))==NULL)
{printf("Open file error!Please check the file.\n");
}
fscanf(f,"%d %d\n",&ncount,&mcount);/*得到矩阵维数*/
mcount++;/*列数加1*/
for(i=1;i<=ncount;i++){/*读X矩阵数据n行*m列*/
Aarray[i][1]=1;/*第一列置1*/
for(j=2;j<=mcount;j++)
{
fscanf(f,"%lf ",&Aarray[i][j]);
fscanf(f,"\n");
}
}
fscanf(f,"\n");
for(i=1;i<=ncount;i++){/*读Y矩阵数据n列*/
fscanf(f,"%lf\n",&Barray[i]);
}
printf("All data read over.\n");

OLS();
}

/*=======================================================*/

/*========================OLS.c=============================*/

#include <..\Array\getvalue.c>


OLS()/*多元最小二乘法*/
{
int i,j,k;
double vtemp=0;

for(i=1;i<=ncount;i++){/*X矩阵数据n*m */
for(j=1;j<=mcount;j++)
{
printf("%5.8f ",Aarray[i][j]);
}
printf("\n");
}
for(i=1;i<=mcount;i++){/*计算X'X矩阵数据m*m */
for(j=1;j<=mcount;j++)
{
array[i][j]=0;
for(k=1;k<=ncount;k++){
array[i][j]=array[i][j]+Aarray[k][i]*Aarray[k][j];
}
}
}

printf("X'X over\n");
for(i=1;i<=mcount;i++){/*显示X'X矩阵数据m*m */
for(j=1;j<=mcount;j++)
{
printf("%5.8f ",array[i][j]);
}
printf("\n");
}

printf("The value is: %10.15f\n",getvalue());
vtemp=1/getvalue();
printf("The 1/value is: %10.15f\n",vtemp);
for(i=1;i<=mcount;i++){/*生成逆矩阵数据*/
for(j=1;j<=mcount;j++)
{
Carray[i][j]=getavalue(j,i)*getov(i+j)*vtemp;
/*printf("%5.2f ",Carray[i][j]);*/
}
}
printf("Iarray over\n");
for(i=1;i<=mcount;i++){/*显示逆矩阵数据m*m */
for(j=1;j<=mcount;j++)
{
printf("%10.15f ",Carray[i][j]);
}
printf("\n");
}

for(i=1;i<=mcount;i++){/*计算AX'矩阵数据M行*N列 */
for(j=1;j<=ncount;j++)
{
for(k=1;k<=mcount;k++){
Darray[i][j]=Darray[i][j]+Carray[i][k]*Aarray[j][k];
}
}
}
printf("The result is:");
for(i=1;i<=mcount;i++){/*计算DY矩阵数据M行*/
for(k=1;k&

lt;=ncount;k++){
Rarray[i]=Rarray[i]+Darray[i][k] * Barray[k];
}
printf("%10.15f ",Rarray[i]);
}
printf("\n");
}

/*=======================================================*/

/*========================getvalue.c=============================*/
#include <..\array\head.h>
#include <..\array\getord
er.c>

double value=0.0;

double getvalue()/*得到矩阵值*/
{ value=0.0;
creatorder(mcount);
return value;
};

creatorder(int k)/*生成全排列序列*/
{
int i,j,m;
double val;/*一排列元素的值*/
if(k>=1)/*未生成完毕*/
{
for(i=mcount;i>=1;i--)
{
if(check(i,k))/*若该序数未分配,则分配*/
{
order[k]=i;
creatorder(k-1);/*分配下一个*/
}
}
}
else/*一个序列生成完毕*/
{
val=getorder(order,mcount);
for(m=1;m<=mcount;m++){
val=val*array[m][order[m]];
}
value=value+val;
}
};
double getavalue(int X,int Y)/*得到伴随矩阵值*/
{ value=0.0;
creataorder(mcount,X,Y);
return value;
};

creataorder(int k,int X,int Y)/*生成全排列序列*/
{
int i,j,m;
double val;/*一排列元素的值*/
if(k>=1)/*未生成完毕*/
{
for(i=mcount;i>=1;i--)
{
if(check(i,k))/*若该序数未分配,则分配*/
{
order[k]=i;
creataorder(k-1,X,Y);/*分配下一个*/
}
}
}
else/*一个序列生成完毕*/
{
val=getaorder(order,mcount,X);
for(m=1;m<=mcount;m++){
if(m!=X){
if(order[m]!=Y)
val=val*array[m][order[m]];
else
val=0.0;
}
}

/*
for(m=1;m<=mcount;m++){
printf("%d ",order[m]);
}
printf("=> ");
for(m=1;m<=mcount;m++){
if(m!=X){
if(order[m]!=Y)
printf("%3.2f ",array[m][order[m]]);
else
printf("invaluable");
}
}
printf("\n");*/
value=value+val;
}
};

int check(int m,int k)/*检查该序数是否已分配*/
{
int i;
for(i=mcount;i>k;i--)
if(m==order[i])
return 0;
return 1;
}

/*=======================================================*/

/*========================getorder.c=============================*/

#include <stdio.h>

int getorder(int Order[],int COUNT)/*得到序数*/
{
int counter=0;
int i=1,j=1;
for(i=1;i<=COUNT;i++)
{
for(j=1;j<=i;j++)
{

if(Order[i]<Order[j])
{
counter++;
}
}
}
return getov(counter);
};

int getaorder(int Order[],int COUNT,int X)/*得到伴随矩阵的序数,X为除去的序数*/
{
int counter=0;
int i=1,j=1;
for(i=1;i<=COUNT;i++)
{
for(j=1;j<=i;j++)
{
if(i!=X && j!=X)
if(Order[i]<Order[j])
counter++;
}
}
return getov(counter);
};

int getov(int c){
if((c%2)==0) return 1;/*return(counter);得到序数的奇偶,奇则为-1否则为1*/
return -1;
}

/*=======================================================*/

/*========================head.h=============================*/

double Aarray[25][25];
double Barray[25];
double Rarray[25];
double Carray[25][25];
double

Darray[25][25];
double array[25][25];
int ncount,mcount;/*维数*/
int order[25];