fmain.cpp

非线形最小二乘程序,希望对大家有用哈,我现在急需APRIOR算法的C++程序,有的请给我发啊一份,

#include <iostream.h>
#include <math.h>
#include "Matrix.h"

double fun(double t, double m, double a, double& dm, double& da)
{
	double	fenzi = 282 ;

	double	fenmu_1 = 1 ;
	double	fenmu_2 = m*exp(0-a*t) ;
	double	fenmu = 1+fenmu_2 ;
	dm = 0-(1/(fenmu*fenmu))*282*exp(0-a*t) ;
	da = 1/(fenmu*fenmu)* 282*m*t*exp(0-a*t) ;
//	cout<<fenzi/fenmu<<endl;

	return	fenzi / fenmu ;
}

#define	M  2//方程的维数（系数矩阵的维数）
#define	F  9
void main()
{
	double w = 0.5 ; //当ｄ>100时的乘数ｗ值
	//自变量t, 函数值y,　参数b，ｂ0，残差ｅ，残差平方和Ｑ，偏导数ｄＢ，系数矩阵Ａ,校正值Ｘ
	double t[F] = { 11, 21, 31, 41, 51, 61, 66, 69, 72 } ;
	double y[F] = { 20.5, 47.1, 83.5, 171, 226.7, 259.7, 271.3, 275.3, 280.3 } ;
	double b[M] = { 80, 0.2 };
	double b0[M] = { 80, 0.2 };
	double b1[M] ;
	//double y[9] = { 14.75057909,42.00247589,100.6422118,179.8121533,239.1396291,266.913152,273.3240223,275.8071105,277.5909959 } ;
	//double b[2] = { 64.4764473660417, 0.1154 };
	//double b[2] = { 64.47645, 0.1154 };
	
//	double t[F] = { 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31 } ;
//	double y[F] = { 0.67,0.85,1.28,1.75,2.27,2.75,3.69,4.71,6.36,7.73,9.91,12.75,16.55,20.1,27.35,32.55,37.55,44.75,53.38,71.61,83.89,97.46,112.73,135.12,153.6,160.32,167.05,174.9,177.87,180.19,180.79} ;
//	double b[M] = { 393, 0.22965 };
//	double b0[M] = { 393, 0.22965 };//记录ｂ的初值
	double eps ; //允许误差
	int	maxp ; //最大迭代次数
	double e[F];
	double Q = 0 ;
	double dB[M];
	CMatrix A(M,M+1) ;
	CVector X(M) ;
	CVector	H(M+1) ; //过渡数组，标准化Ａ用
	bool bflag =false ;

		
		double d = 0.001 ;
		double c =10 ;
		int a = -1 ;
		d = d*pow(c,a);
		CMatrix Atest ;

	for ( int nrow=0; nrow<M; nrow++ ) 
	{
		dB[nrow] = 0;
		for ( int ncol=0; ncol<M+1; ncol++ ) 
		{
			A[nrow][ncol]=0 ;
		}
	}

	for ( int i=0; i<F; i++) 
	{
		//残差
		e[i] = y[i] - fun( t[i], b[0], b[1], dB[0], dB[1] ) ;
				
		//残差平方和
		Q = Q + e[i] * e[i] ;
	}

	double Q1 = Q ;
	double Q2 = Q ;
	double Q3 ;

	do
	{
		if (d<100) 
		{
			for ( int nrow=0; nrow<M; nrow++ ) 
			{
				b1[nrow] = b[nrow] ;
				for ( int ncol=0; ncol<M+1; ncol++ ) 
				{
					A[nrow][ncol]=0 ;
				}
			}

			for ( int i=0; i<F; i++) 
			{
				fun( t[i], b[0], b[1], dB[0], dB[1] ) ;
				//计算Ａij，Ａiy
				for ( int drow=0; drow<M; drow++ ) 
				{
					A[drow][M] = A[drow][M] + dB[drow] * e[i];
					
					for ( int dcol=0; dcol<M; dcol++ ) 
					{
						A[drow][dcol] = A[drow][dcol] + dB[drow] * dB[dcol];
					}
				}
			}



			//标准化
			A.StandardNonLinear();

			for (int ib=0;ib<M;ib++) 
			{
				//cout<<A[0][0]<<endl;
				A[ib][M] = A[ib][M]/sqrt(A[ib][ib])/sqrt(Q1) ;
				//cout<<A[ib][M]<<endl;
			}

			//增加参数d
			for ( int iadd1 = 0; iadd1 < M; iadd1++ ) 
			{
				for ( int iadd2 = 0; iadd2 < M; iadd2++ ) 
				{
					A[iadd1][iadd2]=A[iadd1][iadd2]+d ;
				}
			}

			//解校正值方程
			X = A.ColMax() ;
			/////
			
			//校正参数值
			for(int itest=0;itest<M;itest++)
			{
				b[itest] = b[itest] + X[itest] ;
			//	cout<<"b"<<itest+1<<"="<<b[itest]<<endl;
			}

			////
			Q3 = 0 ;
			Q3 = fabs(X[0]/b[0]) ;
			//计算新Q值
			Q2 = 0 ;
			for ( int i2=0; i2<F; i2++) 
			{
				//残差
				e[i2] = y[i2] - fun( t[i2], b[0], b[1], dB[0], dB[1] ) ;
						
				//残差平方和
				Q2 = Q2 + e[i2] * e[i2] ;
			}
			if (Q2<Q1) 
			{
			//	A=Atest;
				Q1=Q2 ;
				bflag = true ;
			}
			else
			{
				bflag = false ;
				a++;
				d=pow(c,a);
				for (int ib1=0;ib1<M;ib1++) 
				{
					b[ib1]=b1[ib1];
				}
			}

			/////////////////////////////////////////////////////////////////////////
		}
		else
		{
			for(int i0=0;i0<M;i0++)
			{
				b[i0] = b1[i0] + w*X[i0] ;
			//	cout<<"b"<<itest+1<<"="<<b[itest]<<endl;
			}
			a = -1 ;
			d=0.001 ;
			d = d*pow(c,a);
		}

	}
	while( (Q3>0.001)&&(!bflag) );
		//校正参数值
		for(int iout=0;iout<2;iout++)
		{
			cout<<"b"<<iout+1<<"="<<b[iout]<<endl;
		}
	
	
 }

//void main()
//{
//	double dA[2][3] = {1,2,4,6,3,8};//增广矩阵
//	CMatrix A(2,3) ;
//	for ( int i = 0; i < 2; i++ ) 
//	{
//		for ( int j = 0; j < 3; j++ ) 
//		{
//			A[i][j] = dA[i][j] ;
//		}
//	}
//	A.ColMax();
// }