|  | 
 
 
 楼主|
发表于 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;
 }
 
 | 
 |