|
楼主 |
发表于 2020-6-13 01:59:24
|
显示全部楼层
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define N 200
#define P 4
int MatrixInvSymmetric(double *a,int n)//求逆
{
int i,j,k,m;
double w,g,*b;
b = (double*)malloc(n*sizeof(double));
for (k=0; k<n; k++)
{
w=a[0];
if (fabs(w)<1.0e-8)
{
free(b);
return(-2);
}
m=n-k-1;
for (i=1; i<n; i++)
{
g=a[i*n];
b[i]=g/w;
if (i<=m)
b[i]=-b[i];
for (j=1; j<=i; j++)
a[(i-1)*n+j-1]=a[i*n+j]+g*b[j];
}
a[n*n-1]=1.0/w;
for (i=1; i<n; i++)
a[(n-1)*n+i-1]=b[i];
}
for (i=0; i<n-1; i++)
for (j=i+1; j<n; j++)
a[i*n+j]=a[j*n+i];
free(b);
return(2);
}
void readBinaryFile(char *filenames, double *arr, int n)
{
//数据的读取
FILE *fp;
if((fp=fopen(filenames,"rb"))==NULL)
{
printf("cannot open file!\n");
exit(0);
}
fread(arr,sizeof(double),n,fp);
fclose(fp);
}
void EstimateBeta(double *β,double *Y,double *X)
{
double M[N*N]= {0},Q[P*1]= {0};
int i,j,k;
for(i=0; i<N; i++)
{
for(j=0; j<N; j++)
{
for(k=0; k<N; k++)
{
M[N*i+j]+=X[i*N+k]*X[j*N+k];
}
}
}
MatrixInvSymmetric(M,N);
for(i=0; i<P; i++)
{
for(j=0; j<N; j++)
{
Q[i]+=X[j*N+i]*Y[j];
}
}
for(i=0; i<P; i++)
{
for(j=0; j<P; j++)
{
β[i]+=M[i*P+j]*Q[j];
}
}
}
int main()
{
double Y[N], X[N*P],β[P]= {0};
readBinaryFile("yb.dat",Y,N);
readBinaryFile("xb.dat",X,N*P);
EstimateBeta(β,Y,X);
return 0;
}
|
|