aaaa.rar_曲线拟合

#include "matrixf.h" #include <complex> #include <iostream> #include<math.h> main() { int M,N,i,j; double s=0.0,sy=0.0,sy2=0.0,sup=0.0,sdow1=0.0,sdow2=0.0,sdow=0.0,r; do { printf("请输入拟和次数：");//输入拟和次数 scanf("%d",&N); printf("有几对数据："); scanf("%d",&M); if(M<N+1) {printf("数据对数不够请再次输入数据:\n");} }while(M<N+1); N=N+1; MatrixTwo X,XT,MX,Y,MMX,b,Y2;//申请空间 MatrixOne T(M); X.resize(M); for(i=0;i<M;i++) { X[i].resize(N); } Y.resize(M); for(i=0;i<M;i++) { Y[i].resize(1); } XT.resize(N); for(i=0;i<N;i++) { XT[i].resize(M); } MX.resize(N); for(i=0;i<N;i++) { MX[i].resize(N); } MMX.resize(N); for(i=0;i<N;i++) { MMX[i].resize(M); } b.resize(N); for(i=0;i<N;i++) { b[i].resize(1); } Y2.resize(M); for(i=0;i<M;i++) { Y2[i].resize(1); } for(i = 0; i < M; i++)//输入原始的X数据 { printf("please input X[%d][1]=",i); scanf("%lf",&X[i][1]); } /*for(i = 0; i < M; i++) { for (j = 0; j < N; j++) { printf("%6.2lf ", X[i][j]); } printf("\n"); } printf("\n");*/ for(i=0;i<M;i++)//扩展为X矩阵 { for(j=0;j<N;j++) { if(j==1) { } else { X[i][j]=pow(X[i][1],j); } } } for(i = 0; i < M; i++) { for (j = 0; j < N; j++) { printf("%6.2lf ", X[i][j]); } printf("\n"); } for(i = 0; i < M; i++)//输入Y矩阵 { printf("please input Y[%d][0]=",i); scanf("%lf",&Y[i][0]); } for(i = 0; i < M; i++) { printf("%6.2lf ", Y[i][0]); printf("\n"); } MatrixTrans(X,XT,M,N);//X矩阵的转置 /*printf("转置后的矩阵为\n"); for(i = 0; i < N; i++) { for (j = 0; j < M; j++) { printf("%6.2lf ", XT[i][j]); } printf("\n"); }*/ MatrixMulti(XT,X,MX,N,M,N);//X矩阵转置左乘X矩阵 //printf("相乘后的矩阵为：\n"); /*for(i = 0; i < N; i++) { for (j = 0; j < N; j++) { printf("%6.2lf ", MX[i][j]); } printf("\n"); }*/ AntiMatrix(MX,T,N,N);//求逆 printf("矩阵求逆为：\n"); for(i = 0; i < N; i++) { for (j = 0; j < N; j++) { printf("%6.2lf ", MX[i][j]); } printf("\n"); } MatrixMulti(MX,XT,MMX,N,N,M); //printf("再次相乘后矩阵为：\n"); /*for(i = 0; i < N; i++) { for (j = 0; j < M; j++) { printf("%6.2lf ", MMX[i][j]); } printf("\n"); }*/ MatrixMulti(MMX,Y,b,N,M,1);//得出b的估计值 printf("b的最小二乘估计值为：\n"); for(i = 0; i < N; i++) { for (j = 0; j < 1; j++) { printf("b%d=%lf ",i,b[i][j]); } printf("\n"); } MatrixMulti(X,b,Y2,M,N,1);//得出Y的测量值 printf("Y的测量值为：\n"); for(i = 0; i < M; i++) { for (j = 0; j < 1; j++) { printf("Y[%d]=%lf ",i,Y2[i][j]); } printf("\n"); } for(i=0;i<M;i++)//误差的平方和 { for(j=0;j<1;j++) { s+=pow((Y[i][j]-Y2[i][j]),2); } } printf("各点误差的平方和为：s=%lf\n",s); for(i=0;i<M;i++)//相关系数求解 { for(j=0;j<1;j++) { sy+=Y[i][j]; sy2+=Y2[i][j]; } } sy=sy/M; sy2=sy2/M; for(i=0;i<M;i++) { for(j=0;j<1;j++) { sup+=(Y[i][j]-sy)*(Y2[i][j]-sy2); sdow1+=pow((Y[i][j]-sy),2); sdow2+=pow((Y2[i][j]-sy2),2); } } sdow=sqrt((sdow1*sdow2)); r=sup/sdow; printf("相关系数为%lf\n",r); return 0; }