aa.cpp

直线拟合的几种算法

		}
		e[mm-2]=fg[0];
		it=it-1;
		}
		else
		{ if (ks==mm)
		{ kk=kk+1;
		fg[1]=e[mm-2]; e[mm-2]=0.0;
		for (ll=kk; ll<=mm-1; ll++)
		{ i=mm+kk-ll-1;
		fg[0]=s[i-1];
		sss(fg,cs);
		s[i-1]=fg[0];
		if (i!=kk)
		{ fg[1]=-cs[1]*e[i-2];
		e[i-2]=cs[0]*e[i-2];
		}
		if ((cs[0]!=1.0)||(cs[1]!=0.0))
			for (j=1; j<=n; j++)
			{ ix=(j-1)*n+i-1;
			iy=(j-1)*n+mm-1;
			d=cs[0]*v[ix]+cs[1]*v[iy];
			v[iy]=-cs[1]*v[ix]+cs[0]*v[iy];
			v[ix]=d;
			}
		}
		}
		else
		{ kk=ks+1;
		fg[1]=e[kk-2];
		e[kk-2]=0.0;
		for (i=kk; i<=mm; i++)
		{ fg[0]=s[i-1];
		sss(fg,cs);
		s[i-1]=fg[0];
		fg[1]=-cs[1]*e[i-1];
		e[i-1]=cs[0]*e[i-1];
		if ((cs[0]!=1.0)||(cs[1]!=0.0))
			for (j=1; j<=m; j++)
			{ ix=(j-1)*m+i-1;
			iy=(j-1)*m+kk-2;
			d=cs[0]*u[ix]+cs[1]*u[iy];
			u[iy]=-cs[1]*u[ix]+cs[0]*u[iy];
			u[ix]=d;
			}
		}
		}
		}
          }
      }
	  free(s); free(e); free(w);
	  return(1);
}

int dngin(int m, int n, double eps1, double eps2, double *x, int ka)
{
 	//函数介绍：非线性方程组最小二乘解的广义逆法；
 	//m:非线性方程组中方程个数；n：非线性方程组中未知数个数
 	//eps1,控制最小二乘解的精度要求；esp2:用于奇异值分解中的控制精度要求
 	//x[]:一维数组，长度n，存放初始近似解，返回时存放最小二乘解a，当m=n时，即为非线性方程组的解；
 	//ka=max(m,n)+1
     int i,j,k,l,kk,jt;
     double y[15],b[15];
 	double alpha,z,h2,y1,y2,y3,y0,h1;
     double *p,*d,*pp,*dx,*u,*v,*w;
     p=(double *)malloc(m*n*sizeof(double));
     d=(double *)malloc(m*sizeof(double));
     pp=(double *)malloc(n*m*sizeof(double));
     dx=(double *)malloc(n*sizeof(double));
     u=(double *)malloc(m*m*sizeof(double));
     v=(double *)malloc(n*n*sizeof(double));
     w=(double *)malloc(ka*sizeof(double));
     l=60; alpha=1.0;
     while (l>0)
 	{ dnginf(m,n,x,d);
 	dngins(m,n,x,p);
 	jt=::agmiv(p,m,n,d,dx,pp,eps2,u,v,ka);
 	if (jt<0)
 	{ 
 		free(p); free(d); free(pp); free(w);
 		free(dx); free(u); free(v); return(jt);
 	}
 	j=0; jt=1; h2=0.0;
 	while (jt==1)
 	{ 
 		jt=0;
 	if (j<=2) 
 		z=alpha+0.01*j;
 	else z=h2;
 	for (i=0; i<=n-1; i++) 
 		w[i]=x[i]-z*dx[i];
 	dnginf(m,n,w,d);
 	y1=0.0;
 	for (i=0; i<=m-1; i++) 
 		y1=y1+d[i]*d[i];
 	for (i=0; i<=n-1; i++)
 		w[i]=x[i]-(z+0.00001)*dx[i];
 	dnginf(m,n,w,d);
 	y2=0.0;
 	for (i=0; i<=m-1; i++) 
 		y2=y2+d[i]*d[i];
 	y0=(y2-y1)/0.00001;
 	if (fabs(y0)>1.0e-10)
 	{ h1=y0; h2=z;
 	if (j==0) { y[0]=h1; b[0]=h2;}
 	else
 	{ 
 		y[j]=h1; 
 	kk=0; 
 	k=0;
 	while ((kk==0)&&(k<=j-1))
 	{ 
 		y3=h2-b[k];
 	if (fabs(y3)+1.0==1.0) 
 		kk=1;
 	else 
 		h2=(h1-y[k])/y3;
 	k=k+1;
 	}
 	b[j]=h2;
 	if (kk!=0) 
 		b[j]=1.0e+35;
 	h2=0.0;
 	for (k=j-1; k>=0; k--)
 		h2=-y[k]/(b[k+1]+h2);
 		h2=h2+b[0];
 	}
 	j=j+1;
 	if 
 		(j<=7) jt=1;
 	else z=h2;
 	}
 	}
 	alpha=z; y1=0.0; y2=0.0;
 	for (i=0; i<=n-1; i++)
 	{ dx[i]=-alpha*dx[i];
 	x[i]=x[i]+dx[i];
 	y1=y1+fabs(dx[i]);
 	y2=y2+fabs(x[i]);
 	}
 	if (y1<eps1*y2)
 	{ free(p); free(pp); free(d); free(w);
 	free(dx); free(u); free(v); return(1);
 	}
 	l=l-1;
 	}
     free(p); free(pp); free(d); free(dx);
     free(u); free(v); free(w); return(0);
}

void dnginf(int m, int n, double x[], double d[])
{
 	int i;
 	int num=numbers;/*PointNum*/;
	double ddelt;
	for(i=0; i<num; i++)
	{
		ddelt = (x[0]*xxx[i]+x[1]-yyy[i]);
		d[i] = ddelt * ddelt;
	}
}

void dngins(int m, int n, double x[], double *p)
{
 	int i;
 	int num=numbers;/*PointNum*/;
	double ddelt;
	for(i=0; i<num; i++)
	{
		ddelt = (x[0]*xxx[i]+x[1]-yyy[i]);
		p[2*i+0] = 2 * ddelt * xxx[i];
		p[2*i+1] = 2 * ddelt;
	}
}

int dngin2(int m, int n, double eps1, double eps2, double *x, int ka)
{
 	//函数介绍：非线性方程组最小二乘解的广义逆法；
 	//m:非线性方程组中方程个数；n：非线性方程组中未知数个数
 	//eps1,控制最小二乘解的精度要求；esp2:用于奇异值分解中的控制精度要求
 	//x[]:一维数组，长度n，存放初始近似解，返回时存放最小二乘解a，当m=n时，即为非线性方程组的解；
 	//ka=max(m,n)+1
     int i,j,k,l,kk,jt;
     double y[15],b[15];
 	double alpha,z,h2,y1,y2,y3,y0,h1;
     double *p,*d,*pp,*dx,*u,*v,*w;
     p=(double *)malloc(m*n*sizeof(double));
     d=(double *)malloc(m*sizeof(double));
     pp=(double *)malloc(n*m*sizeof(double));
     dx=(double *)malloc(n*sizeof(double));
     u=(double *)malloc(m*m*sizeof(double));
     v=(double *)malloc(n*n*sizeof(double));
     w=(double *)malloc(ka*sizeof(double));
     l=60; alpha=1.0;
     while (l>0)
 	{ dnginf2(m,n,x,d);
 	dngins2(m,n,x,p);
 	jt=::agmiv(p,m,n,d,dx,pp,eps2,u,v,ka);
 	if (jt<0)
 	{ 
 		free(p); free(d); free(pp); free(w);
 		free(dx); free(u); free(v); return(jt);
 	}
 	j=0; jt=1; h2=0.0;
 	while (jt==1)
 	{ 
 		jt=0;
 	if (j<=2) 
 		z=alpha+0.01*j;
 	else z=h2;
 	for (i=0; i<=n-1; i++) 
 		w[i]=x[i]-z*dx[i];
 	dnginf2(m,n,w,d);
 	y1=0.0;
 	for (i=0; i<=m-1; i++) 
 		y1=y1+d[i]*d[i];
 	for (i=0; i<=n-1; i++)
 		w[i]=x[i]-(z+0.00001)*dx[i];
 	dnginf2(m,n,w,d);
 	y2=0.0;
 	for (i=0; i<=m-1; i++) 
 		y2=y2+d[i]*d[i];
 	y0=(y2-y1)/0.00001;
 	if (fabs(y0)>1.0e-10)
 	{ h1=y0; h2=z;
 	if (j==0) { y[0]=h1; b[0]=h2;}
 	else
 	{ 
 		y[j]=h1; 
 	kk=0; 
 	k=0;
 	while ((kk==0)&&(k<=j-1))
 	{ 
 		y3=h2-b[k];
 	if (fabs(y3)+1.0==1.0) 
 		kk=1;
 	else 
 		h2=(h1-y[k])/y3;
 	k=k+1;
 	}
 	b[j]=h2;
 	if (kk!=0) 
 		b[j]=1.0e+35;
 	h2=0.0;
 	for (k=j-1; k>=0; k--)
 		h2=-y[k]/(b[k+1]+h2);
 		h2=h2+b[0];
 	}
 	j=j+1;
 	if 
 		(j<=7) jt=1;
 	else z=h2;
 	}
 	}
 	alpha=z; y1=0.0; y2=0.0;
 	for (i=0; i<=n-1; i++)
 	{ dx[i]=-alpha*dx[i];
 	x[i]=x[i]+dx[i];
 	y1=y1+fabs(dx[i]);
 	y2=y2+fabs(x[i]);
 	}
 	if (y1<eps1*y2)
 	{ free(p); free(pp); free(d); free(w);
 	free(dx); free(u); free(v); return(1);
 	}
 	l=l-1;
 	}
     free(p); free(pp); free(d); free(dx);
     free(u); free(v); free(w); return(0);
}

void dnginf2(int m, int n, double x[], double d[])
{
 	//函数介绍：列入等号左边的行向量
 	int i;
 	int num=numbers;/*PointNum*/;
	for(i=0; i<num; i++)
	{
		d[i] = 1000*(x[0]*xxx[i]+x[1]-yyy[i]);
	}
}

void dngins2(int m, int n, double x[], double *p)
{
 	int i,j;
 	int num=numbers;/*PointNum*/;
	for(i=0;i<m;i++)
		for(j=0;j<n;j++)
			p[n*i+j]=0;
	for(i=0; i<num; i++)
	{
		p[2*i+0] = 1000*xxx[i];
		p[2*i+1] = 1000;
	}
}

#define PI 3.1415926535897932384626433832795
#define CAAngle(x) fabs(atan(x) * 180.0 / PI)

int main(int argc, char* argv[])
{
	xxx = new double[numbers];
	yyy = new double[numbers];
	FILE* fp = fopen("aa.txt","r");
	if(fp)
	{
		for(int i = 0; i < numbers; i ++)
		{
			fscanf(fp,"%lf %lf",xxx + i, yyy +i);
		}
		fclose(fp);
	}

	//////////////////////////////////////////////////////////////////////////
	// 斜率平均
	{
		double avgx = 0.0,avgy = 0.0;
		for(int i = 0; i < numbers; i ++)
		{
			avgx += xxx[i];
			avgy += yyy[i];
		}
		avgx /= numbers;
		avgy /= numbers;
		
		double kavg = 0.0;
		double angavg = 0.0;
		for(i = 0; i < numbers; i ++)
		{
			kavg = fabs((avgy - yyy[i]) / (avgx - xxx[i]));
			// 		angavg += atan(kavg) * 180 / 3.1415926;
		}
		
		kavg /= numbers;

		printf("\n斜率平均: %lf\n",CAAngle(kavg));
	}
	//////////////////////////////////////////////////////////////////////////
	// 线性最小二乘
	{
		double sumX = 0.0, sumY = 0.0, sumXX = 0.0, sumXY = 0.0;
		for(int i=0; i<numbers; i++)
		{
			sumX += xxx[i];
			sumY += yyy[i];
			sumXX += xxx[i] * xxx[i];
			sumXY += yyy[i] * yyy[i];
		}
		
		double	qa=(sumY*sumX-sumXY*numbers)/(sumX*sumX-sumXX*numbers);
		double qb=(sumY-sumX*qa)/numbers;
		printf("\n线性最小二乘: %lf\n",CAAngle(qa));
	}
	/////////////////////////////////////////////////////////////////////////
	// 非线性最小二乘(平方)
	{
		int num = numbers;
		int m,n,ka;
		double eps1,eps2;
		
		m=num;
		n=2;
		ka=m>n ? m:n;
		ka=ka+1;
		eps1=0.00001;eps2=0.00001;
		double x[2];
		x[0]=-30.0; x[1]=0.0;
		int i=dngin(m,n,eps1,eps2,x,ka);
		
		double qa = x[0];
		double qb = x[1];
		printf("\n非线性最小二乘(平方): %lf\n",CAAngle(qa));
	}
	/////////////////////////////////////////////////////////////////////////
	// 非线性最小二乘(相等)
	{
		int num = numbers;
		int m,n,ka;
		double eps1,eps2;
		
		m=num;
		n=2;
		ka=m>n ? m:n;
		ka=ka+1;
		eps1=0.00001;eps2=0.00001;
		double x[2];
		x[0]=-30.0; x[1]=0.0;
		int i=dngin2(m,n,eps1,eps2,x,ka);
		
		double qa = x[0];
		double qb = x[1];
		printf("\n非线性最小二乘(相等): %lf\n",CAAngle(qa));
	}
	delete [] xxx;
	delete [] yyy;

	getchar();
	return 0;
}