最小二乘曲线拟和.cpp

简单的最小二乘曲线拟合方法

#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;


}