【最小二乘曲线拟合】

//least square curve fitting.cpp
#include <iostream>
#include <cmath>
using namespace std;
//x[n]	Store the X coordinate of a given data point	
//y[n]	Store the Y coordinate of a given data point	
//n		Given the number of data points	
//a[m]	Returns the coefficients of the m-1 degree fitting polynomial
//m	The Number of terms of the fitting polynomial. Requirementsm<=min(n,20)	
/*	
dt[3]	Return the sum of the squared error, the sum of the absolute value 
of the error and the maximum value of the absolute value of the error
*/
void pirl(double x[], double y[], int n, double a[], int m, double dt[])
{
	int i,j,k;
	double p,c,g,q,d1,d2,s[20],t[20],b[20]; 
	for(i=0;i<=m-1;i++)a[i]=0.0; 
	if(m>n)m=n; 
	if(m>20)m=20;
	b[0]=1.0;
	d1=1.0*n;
	p=0.0;c=0.0; 
	for (i=0;i<=n-1;i++)
		{
			p=p+x[i]; c=c+y[i];
		}
	c=c/d1;
	p=p/d1; 
	a[0]=c*b[0];
	if(m>1)	
	{
		t[1]=1.0;t[0]=-p; 
		d2=0.0;c=0.0;g=0.0; 
		for (i=0; i<=n-1;i++)
			{
				q=x[i]-p;
				d2=d2+q*q; 
				c=c+y[i]*q; 
				g=g+x[i]*q*q;
			}
		c=c/d2;p=g/d2;q=d2/d1; 
		d1=d2;
		a[1]=c*t[1];
		a[0]=c*t[0]+a[0];
	}
	for (j=2; j<=m-1; j++)
		{
			s[j]=t[j-1];	
			s[j-1]=-p*t[j-1]+t[j-2];	
			if (j>=3)
				for (k=j-2; k>=1; k--)
				s[k]=-p*t[k]+t[k-1]-q*b[k]; 
			s[0]=-p*t[0]-q*b[0]; d2=0.0; 
			c=0.0;g=0.0; 
			for (i=0; i<=n-1; i++)
				{
					q=s[j];
					for (k=j-1;k>=0;k--) 	q=q*x[i]+s[k]; 
					d2=d2+q*q; c=c+y[i]*q; 
					g=g+x[i]*q*q;
				}

			c=c/d2;p=g/d2;q=d2/d1; 
			d1=d2;
			a[j]=c*s[j];
			t[j]=s[j]; 
			for (k=j-1; k>=0; k--)
			{
				a[k]=c*s[k]+a[k]; 
				b[k]=t[k];
				t[k]=s[k];
			}
		}
		dt[0]=0.0; dt[1]=0.0; 
		dt[2]=0.0; 
		for (i=0; i<=n-1; i++)
		{
			q=a[m-1];
			for(k=m-2; k>=0; k--) 
			q=a[k]+q*x[i]; 
			p=q-y[i];
			if (fabs(p)>dt[2]) dt[2]=fabs(p);
			dt[0]=dt[0]+p*p;
			dt[1]=dt[1]+fabs(p);
		}
	return ;
}
int main()
{
	int i;
	double x[20],y[20],a[6],dt[3]; 
	for (i=0;i<=19;i++)
		{
		x[i]=0.1*i;
		y[i]=x[i]-exp(-x[i]);
		}
	pirl(x,y,20,a,6,dt);
	cout<<"拟合多项式系数:"<<endl; 
	for(i=0;i<=5;i++)
		cout <<"a (" <<i <<") =" <<a[i]<<endl; 
	cout<<"误差平方和 =" <<dt[0] <<endl; 
	cout<<"误差绝对值和="<<dt[1] <<endl;
	cout<<"误差绝对值最大值="<<dt[2] <<endl; 
	return 0;
}

计算结果：