最小二乘法的曲线拟合

最小二乘法的曲线拟合

BOOL CDataDistillView::LeastDoubleMultiplication(LONG *pX, LONG *pY, 
								long M, long N, double *result, double &warp)
{
	//data		--  X,Y两轴的坐标
	//M			--  结果变量组数
	//N			--  采样数目
	//result	--  结果参数
	//warp		--	偏差平方和

	register long i, j, k;
	double Z, D1, D2, C, P, G, Q;
	double *B, *T, *S;

	B = new double[N];
	T = new double[N];
	S = new double[N];

	if(M > N)
	{
		M = N;
	}

	for(i = 0; i < M; i++)
	{
		result[i] = 0;
	}

	Z = 0;
	B[0] = 1;
	D1 = N;
	P = 0;
	C = 0;

	for(i = 0; i < N; i++)
	{
		P = P + pX[i] - Z;
		C = C + pY[i];
	}

	C = C / D1;
	P = P / D1;
	result[0] = C * B[0];

	if(M > 1)
	{
		T[1] = 1;
		T[0] = -P;
		D2 = 0;
		C = 0;
		G = 0;

		for(i = 0; i < N; i++)
		{
			Q = pX[i] - Z - P;
			D2 = D2 + Q * Q;
			C = pY[i] * Q + C;
			G = (pX[i] - Z) * Q * Q + G;
		}

		C = C / D2;
		P = G / D2;
		Q = D2 / D1;
		D1 = D2;
		result[1] = C * T[1];
		result[0] = C * T[0] + result[0];
	}

	for(j = 2; j < M; 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;
		C = 0;
		G = 0;

		for(i = 0; i < N; i++)
		{
			Q=S[j];
			
			for(k = j - 1; k >= 0; k--)
			{
				Q = Q * (pX[i] - Z) + S[k];
			}

			D2 = D2 + Q * Q;
			C = pY[i] * Q + C;
			G = (pX[i] - Z) * Q * Q + G;
		}

		C = C / D2;
		P = G / D2;
		Q = D2 / D1;
		D1 = D2;
		result[j] = C * S[j];
		T[j] = S[j];

		for(k = j - 1;k >= 0;k--)
		{
			result[k] = C * S[k] + result[k];
			B[k] = T[k];
			T[k] = S[k];
		}
	}

	warp = 0;
	double dt;
 	for(i = 0; i < N; i++)   
	{
		Q = result[M - 1];   
		
		for(k = M-1; k>=1; k--)   
		{
			Q = Q * (pX[i] - Z) + result[k - 1];   
		}

		dt = Q - pY[i];   
		
		warp = warp + dt * dt;   
	}

	delete [] B;
	delete [] T;
	delete [] S;
	
	return TRUE;
}