small_secend.c

这是最小二乘法

/*--以下代码从王一波上位机中摘取出来,一种算法之一,以下为他的源代码*/

void CDlgStdCurve::StraightParabolaFitting( int   m,  double x[],  double y[], 	 double a[])
{
     /* 
	 //正规方程组
a0----b
a1----k
stdr---R
stdr2---R2(stdr*stdr)

	 m*a0+xsum*a1+x2sum*a2=ysum
	 xsum*a0+x2sum*a1+x3sum*a2=xysum
	 x2sum*a0+x3sum*a1+x4sum*a2=x2ysum

     A1=m;     B1=xsum;  C1=x2sum;
     A2=xsum;  B2=x2sum; C2=x3sum;
     A3=x2sum; B3=x3sum; C3=x4sum;
	 y1=ysum;  y2=xysum; y3=x2ysum;  
     
	 //求解

	 tp1=(A2*y3-A3*y2)*(A1*B2-A2*B1)
	 -(A1*y2-A2*y1)*(A2*B3-A3*B2);
	 tp2=(A2*C3-A3*C2)*(A1*B2-A2*B1)
	 -(A1*C2-A2*C1)*(A2*B3-A3*B2);	 
	 a2=tp1/tp2;

	 tp1=A1*B2-A2*B1;
	 a1=((A1*y2-A2*y1)-(A1*C2-A2*C1)*a2)/tp1;

	 a0=(y1-(B1*a1+C1*a2))/A1;

	 tp=(m*x2sum-xsum*xsum)*(m*y2sum-ysum*ysum);
	 stdr = (m*xysum-xsum*ysum)/(sqrt(tp));
	*/
	
	double xsum,ysum,x2sum,x3sum,x4sum,y2sum,xysum,x2ysum,tp1,tp2,stdr,a0,a1,a2,
		   A1,A2,A3,B1,B2,B3,C1,C2,C3,y1,y2,y3;
	xsum=ysum=x2sum=x3sum=x4sum=y2sum=xysum=x2ysum=0;

	m=4;
	x[0]=1;
	x[1]=2;
	x[2]=3;
	x[3]=4;

	y[0]=2;
	y[1]=6;
	y[2]=12;
	y[3]=20;

	for(int i=0;i<m;i++)
	{
		xsum +=x[i];
		ysum +=y[i];
        	x2sum+=x[i]*x[i];
        	x3sum+=x[i]*x[i]*x[i];
		x4sum+=x[i]*x[i]*x[i]*x[i];
		y2sum+=y[i]*y[i];
		xysum+=x[i]*y[i];
        	x2ysum+=x[i]*x[i]*y[i];
	}
	if (m>0) 
	{
		A1=m ;    B1=xsum;  C1=x2sum;
		A2=xsum;  B2=x2sum; C2=x3sum;
		A3=x2sum; B3=x3sum; C3=x4sum;
		y1=ysum;  y2=xysum; y3=x2ysum;

		tp1=(A2*y3-A3*y2)*(A1*B2-A2*B1)
			-(A1*y2-A2*y1)*(A2*B3-A3*B2);
        tp2=(A2*C3-A3*C2)*(A1*B2-A2*B1)
			-(A1*C2-A2*C1)*(A2*B3-A3*B2);

		if (tp2!=0) 
		{
			a2=tp1/tp2;
			tp1=A1*B2-A2*B1;
			if (tp1!=0) 
			{
				a1=((A1*y2-A2*y1)-(A1*C2-A2*C1)*a2)
					/tp1;
			}			
			a0=(y1-(B1*a1+C1*a2))/A1;
		}		
	
	}
	tp1=(m*x2sum-xsum*xsum)*(m*y2sum-ysum*ysum);
	if (tp1>0) 
	{
        stdr = (m*xysum-xsum*ysum)/(sqrt(tp1));
		stdr =stdr*stdr;
	}
	
	a[0]=a0;
	a[1]=a1;
	a[2]=a2;
	
	CString s;
	s.Format("a0=%.4f a1=%.4f a2=%.4f r^2=%.4f",a0,a1,a2,stdr);
	ShowMS(s);
	//AfxMessageBox(s);
}

/*依据王一波算法转换而成的C算法*/
//----算法缺点:
//	1、R2值与EXCEL基本一致，略大，因为EXCEL只取四位小数，而C编译后是以浮点数计算。
//	2、A、B值与EXCEL计算出来以后，略大或略小，不基本一致，小则差几，大则差别上百。

void smallsuanfa(uint m,float x[6],float y[6])
{
	uint i;
	float tp1,tp2;
	float stdr,stdr2;
	float xsum,ysum,x2sum,x3sum,x4sum;
	float y2sum,xysum,x2ysum;
	float a0,a1,a2;
	float A1,A2,A3,B1,B2,B3,C1,C2,C3;
	float y1,y2,y3;

	xsum=0;
	ysum=0;
	x2sum=0;
	x3sum=0;
	x4sum=0;
	y2sum=0;
	xysum=0;
	x2ysum=0;
//----------------------------------
	for(i=0;i<m;i++)
	{
		xsum +=x[i];
		ysum +=y[i];
        x2sum+=x[i]*x[i];
        x3sum+=x[i]*x[i]*x[i];
		x4sum+=x[i]*x[i]*x[i]*x[i];
		y2sum+=y[i]*y[i];
		xysum+=x[i]*y[i];
        x2ysum+=x[i]*x[i]*y[i];
	}
//----------------------------------
	if (m>0) 
	{
		A1=m ;    B1=xsum;  C1=x2sum;
		A2=xsum;  B2=x2sum; C2=x3sum;
		A3=x2sum; B3=x3sum; C3=x4sum;
		y1=ysum;  y2=xysum; y3=x2ysum;

		tp1=(A2*y3-A3*y2)*(A1*B2-A2*B1)-(A1*y2-A2*y1)*(A2*B3-A3*B2);
        tp2=(A2*C3-A3*C2)*(A1*B2-A2*B1)-(A1*C2-A2*C1)*(A2*B3-A3*B2);

		if (tp2!=0) 
		{
			a2=tp1/tp2;
			tp1=A1*B2-A2*B1;
			if (tp1!=0) 
			{
				a1=((A1*y2-A2*y1)-(A1*C2-A2*C1)*a2)/tp1;
			}			
			a0=(y1-(B1*a1+C1*a2))/A1;
		}		
	
	}
//------------------------------------
	tp1=(m*x2sum-xsum*xsum)*(m*y2sum-ysum*ysum);
	if (tp1>0) 
	{
        stdr = (m*xysum-xsum*ysum)/(sqrt(tp1));
		stdr2 =stdr*stdr;
	}
//-------------------------------------		

}