leastsquare1.cpp

最小二乘法源代码；可直接用于多项式拟和。

// TT.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "math.h"
#include "stdlib.h"
/*******************************
//求多项市值
*******************************/
double summiax(double *a,int n,double x)
{
	double ret=0;
	for(int i=0;i<=n;i++)//横向元素数目 n+1
	{
		ret+=a[i]*pow(x,i);
	}
	return ret;
}
//指针数值 
#define pa(j,i) (*((double*)a+(j)*(n+1)+(i)))
/*******************************
//累计幂
*******************************/
double summi(double *a,int m,int k)
{
	double re=0;
	for(int i=0;i<m;i++)
	{
		re+=pow(*(a+i),k);
	}
	return re;
}
/*******************************
输入测量数据生成可计算的矩阵
输入数据组数m，x^(n-1)次多项式，m>=n,生成n*(n+1)矩阵
a0+a1*x+a2*x^2+.....ai*x^i+...a(n-1)*x^(n-1)
x,y,n,m,a;
当结果矩阵a使用完成后需要手动释放
*******************************/
double* creatematrix(double *x,double*y,int n,int m)
{
	double *a=new double[n*n+n];
	double *b=new double[2*n+1];
	int i,j;
	for(i=0;i<=2*n;i++)
		*(b+i)=summi(x,m,i);
	for(j=0;j<n;j++)
	{
		for(i=0;i<n;i++)
		{
			pa(j,i)=*(b+i+j);
		}
	}
	for(j=0;j<n;j++)
	{
		pa(j,n)=0;
		for(int i=0;i<m;i++)
		{
			pa(j,n)+=pow(*(x+i),j)**(y+i);
		}
	}
	delete b;
	return a;
}
/*******************************
专门处理矩阵计算,
计算后，矩阵结果将放置在最后一列 a[0,...n-1][n];
********************************/
void dealmatrix(double *a,int n)
{
	double aaa=pa(0,0);
	double b=pa(0,1);
	double c=b/aaa;
	//c=pa(0,1)/pa(0,0);
	//pa(0,1)-=pa(1,1);
	
	int i,j,k;
	//消去参数矩阵左下角、
	for(k=1;k<n;k++)
	{
		for(j=k;j<n;j++)
		{
			double delt =pa(j,k-1)/pa(k-1,k-1);
			pa(j,k-1)=0;//this matix item must be zero
			for(i=k;i<=n;i++)//横向元素数目 n+1
			{
				pa(j,i)-=delt*pa(k-1,i);
			}
		}
	}
	//消去参数矩阵右上角、
	for(k=n-2;k>=0;k--)
	{
		for(j=k;j>=0;j--)
		{
			double delt=pa(j,k+1)/pa(k+1,k+1);
			//a(k+1,j)=0;
			for(i=k;i<n+1;i++)
			{
				pa(j,i)-=delt*pa(k+1,i);
			}
		}
	}
	for( i=0;i<n;i++)
	{
		pa(i,n)/=pa(i,i);
		pa(i,i)=1;
	}
}

/**********************************
最小二乘法曲线拟和
double *x
double *y
int m 测试数据项
int n n-1次多项式
double *outa 返回结果
double *outA 绝对误差
double *outR 相关系数
**********************************/
void CurvefitLeastSquare(double *x,double *y,int m,int n,double *outa,double *outA,double *outR)
{
	double *a=creatematrix(x,y,n,m);
	dealmatrix(a,n);
	for(int i=0;i<n;i++)
	{
		*(outa+i)=pa(i,n);
	}

	//double yn;
	double *yy=new double[m];
	for(i=0;i<m;i++)
	{
		 *(yy+i)=0;
		for(int j=0;j<n;j++)
		{
			*(yy+i)+=pa(j,n)*pow(x[i],j);
		}
	}
	printf("拟和方程:\nY=");
	for( i=0;i<n-1;i++)
	{
		printf("%f*x^%d+",outa[i],i);
	}
	printf("%f*x^%d;",outa[i],i);


	for(i=0;i<m;i++)
	{
		printf("\n真:%f\t计算:%f\t偏差:%f",y[i],yy[i],(yy[i]-y[i])/y[i]);
	}
	if(outA!=NULL||outR!=NULL)
	{
		/*
		//计算绝对误差
		*/
	}
	delete a;
}
/*******************************
多项式求解
a 多项式参数
n 生成n-1项式
xa,ya,xb,yb//多项式单调区间
y
--------
求对应的x
//多项式模拟数据
*******************************/
int main(int argc, char* argv[])
{
	double x[6]={0,2.5,5,10,20,40};
	double y[6]={264,4511,20492,79271,243996,570705};
	double range=100;//测量浓度范围
	int m=6;//测量项数
	double deltb=0;//测量误差范围
	//simulate4paramlogisic(x,y,5,range,deltb);
	int n=2;
	double y0[]={264,4511,20492,79271,243996,570705};
	double x0[]={264,4511,20492,79271,243996,570705};

	double *outa=new double[100];
	CurvefitLeastSquare(x+1,y+1,m-1,n,outa,NULL,NULL);

	for(int i=1;i<m;i++)
	{
	//	y[i]=parseMultiFun(outa,n,x[1],y[1],x[5],y[5],y[i],y+i);
	}
//	Curvefit4ParamLogistic(x,y,m,outa,NULL,NULL);
//	Curvefit3ParamLogistic(x,y,m,outa,NULL,NULL);
//	printf("\n返回参数:d\ta\tb\tc\n%f\t%f\t%f\t%f\t",outa[0],outa[1],outa[2],outa[3]);
	delete outa;
	getchar();
	return 0;
}