stat_utils.c

图像处理的压缩算法

	
			trMR.Calculation.Parameters.P2.Name.strVal		// Parameter Name...<Dataset Name>
			trMR.Calculation.Parameters.P2.Value.dVal
			trMR.Calculation.Parameters.P2.Error.dVal
			trMR.Calculation.Parameters.P2.tValue.dVal
			trMR.Calculation.Parameters.P2.Prob.dVal
	
			etc. depending on number of dependent variables
	
			trMR.Calculation.Stat.Rvalue.dVal				// R-value
			trMR.Calculation.Stat.RSqCOD.dVal				// RSq-COD value
			trMR.Calculation.Stat.AdjRSq.dVal				// Adjusted RSq
			trMR.Calculation.Stat.RMSESD.dVal				// RMSE-Standard Dev.
			trMR.Calculation.Stat.N.nVal					// No. of points in fit
	
			trMR.Calculation.ANOVA.Model.DOF.nVal			// Model DOF 
			trMR.Calculation.ANOVA.Model.SSq.dVal			// Model SS
			trMR.Calculation.ANOVA.Model.MeanSq.dVal		// Model Mean Sq
			trMR.Calculation.ANOVA.Error.DOF.nVal			// Error DOF
			trMR.Calculation.ANOVA.Error.SSq.dVal			// Error SS
			trMR.Calculation.ANOVA.Error.MeanSq.dVal		// Error Mean Sq
			trMR.Calculation.ANOVA.Total.DOF.nVal			// Total degrees of freedom
			trMR.Calculation.ANOVA.Total.SSq.dVal			// Total SS
			trMR.Calculation.ANOVA.Fvalue.dVal				// F-value
			trMR.Calculation.ANOVA.Pvalue.dVal				// P-value
		vWeightData=Input, Optional vector containing input weighting data if used
		mCov=Output, the variance-covariance matrix
	Return:
		Returns ERROR_NO_ERROR on successful exit and an error code on failure.
		NAG error codes:		Negative 1 * NAG error code
		ERROR_INVALID_TREENODE:	Invalid TreeNode
		ERROR_TO_FEW_PTS:		To few data points
		ERROR_UNEQUAL_N:		Unequal N where equal N is required
*/
int stat_multiple_linear_regression(const vector& vDepData, const matrix& mIndepData, 
	TreeNode& trMR, const vector* vWeightData, matrix* mCov)
{
	// Check the input arguments
	if( !trMR.IsValid() )
		return ERROR_INVALID_TREENODE;

	// The number of observations
	int nPts = vDepData.GetSize();

	// NAG nag_regsn_mult_linear function requires nPts > 1
	if( nPts < SU_MIN_NPTS ) 
		return ERROR_TO_FEW_PTS;

	if( mIndepData.GetNumRows() != nPts )
		return ERROR_UNEQUAL_N;

	// The number of independent variables
	Nag_IncludeMean iMean;
	int nPstart; // First parameter that needs to be estimated
	int iTdx = mIndepData.GetNumCols();					
	int iP;
	if( trMR.GUI.Fit.ThroughZero.nVal == 1 )
	{
		iP = iTdx;
		iMean = Nag_MeanZero;
		nPstart = 2;
	}
	else
	{
		iP = iTdx + 1;
		iMean = Nag_MeanInclude;
		nPstart = 1;
	}
	int iTdq = iP + 1;

	// Declare NAG Linear regression arguments
	vector<int> vSx(iTdx);
	vSx = 1;
	double rss, df;					// Resisual sum of squares and the degrees of freedom with rss
	vector vParameter(iP);			// Estimates of the parameters of the regression model
	vector vPErr(iP);				// Standard errors of the parameters
	vector vCov(iP * iTdq / 2); 	// Upper triangular part of the variance-covariance matrix
	vector vRes(nPts);				// Residuals 
	vector vLev(nPts);				// Leverages
	matrix mQ(nPts,iTdq);			// Results of the QR decomposition
	BOOL bSVD;						
	int rank;						// Rank of the independent variables
	vector vP(2 * iP + iP * iP);	// Details of the QR decomposition and SVD
	vector vCOMAR(5 * iTdx + iP * iP);

	int nRet = nag_regsn_mult_linear(iMean, nPts, mIndepData, iTdx, iTdx, vSx, iP, vDepData, *vWeightData,
		&rss, &df, vParameter, vPErr, vCov, vRes, vLev, mQ, iTdq, &bSVD, &rank, vP, SU_DEFAULT_TOLERANCE, vCOMAR);

	if( nRet != NE_NOERROR )
		return -1 * nRet;

	// Calculate the total sum of squares
	double yAve, TSS; // Average value of the dependent variable and the total sum of squares
	vector vtemp(nPts);
	if(vWeightData != NULL)
	{
		double dsumW;
		vtemp = vDepData * (*vWeightData);
		vtemp.Sum(yAve);
		vWeightData->Sum(dsumW);
		yAve /= dsumW;
		vtemp = (vDepData - yAve)^2 * (*vWeightData);
	}
	else
	{
		vtemp = vDepData;
		vtemp.Sum(yAve);
		yAve /= nPts;
		vtemp = (vDepData - yAve)^2;
	}
	vtemp.Sum(TSS);

	// Output the results to trMR
	df = nPts - iTdx - 1; // reset the degree of freedom associated with rss to NumObservation-NumXVar-1

	// Regression parameters - must insure parent node exists
	if( !trMR.Calculation.Parameters.IsValid() )
		trMR.Calculation.AddNode(SU_PARAMS_NODE);

	TreeNode trPara;
	string strTrNode;
	if( trMR.GUI.Fit.ThroughZero.nVal == 1 )
	{
		trPara = trMR.Calculation.Parameters.GetNode(SU_P1);
		if( !trPara.IsValid() )
			trPara = trMR.Calculation.Parameters.AddNode(SU_P1);
		trPara.Value.dVal = 0;
		trPara.Error.dVal = NANUM;
		trPara.tValue.dVal = NANUM;
		trPara.Prob.dVal = NANUM;
	}

	for( int ii = 0; ii < iP; ii++ )
	{
		strTrNode.Format(SU_PN,ii + nPstart);
		trPara = trMR.Calculation.Parameters.GetNode(strTrNode);
		if( !trPara.IsValid() )
			trPara = trMR.Calculation.Parameters.AddNode(strTrNode);
		trPara.Value.dVal = vParameter[ii];
		trPara.Error.dVal = vPErr[ii];
		trPara.tValue.dVal = vParameter[ii] / vPErr[ii];
		trPara.Prob.dVal = 2 * (1 - invt(fabs(trPara.tValue.dVal),df));
	}

	// Stat - Correlation coefficient
	trMR.Calculation.Stat.Rvalue.dVal = sqrt(1 - rss / TSS);				// R value
	trMR.Calculation.Stat.RSqCOD.dVal = 1 - rss / TSS;						// RSq-COD value
	trMR.Calculation.Stat.AdjRSq.dVal = 1 - (rss / TSS * (nPts - 1) / df);	// Adjusted RSq
	trMR.Calculation.Stat.RMSESD.dVal = sqrt(rss / df);						// RMSE-Standard Dev.
	trMR.Calculation.Stat.N.nVal = nPts;									// N

	// ANOVA
	trMR.Calculation.ANOVA.Error.DOF.nVal = df;								// Error DOF
	trMR.Calculation.ANOVA.Error.SSq.dVal = rss;							// Error SS
	trMR.Calculation.ANOVA.Error.MeanSq.dVal = rss/df;						// Error Mean Sq
	trMR.Calculation.ANOVA.Total.DOF.nVal = nPts-1;							// Total degrees of freedom
	trMR.Calculation.ANOVA.Total.SSq.dVal = TSS;							// Total SS
	trMR.Calculation.ANOVA.Model.DOF.nVal = iTdx;							// Model DOF 
	trMR.Calculation.ANOVA.Model.SSq.dVal = TSS - rss;						// Model SS
	trMR.Calculation.ANOVA.Model.MeanSq.dVal = (TSS - rss) / iTdx;			// Model Mean Sq
	trMR.Calculation.ANOVA.Fvalue.dVal = trMR.Calculation.ANOVA.Model.MeanSq.dVal / trMR.Calculation.ANOVA.Error.MeanSq.dVal; // F
	trMR.Calculation.ANOVA.Pvalue.dVal = 1 - invf(trMR.Calculation.ANOVA.Fvalue.dVal, iTdx, df);					// P

	// Fill the variance-covariance matrix
	if( mCov != NULL )
	{
		matrix mTemp(iP,iP);
		mCov->SetSize(iP,iP);
		for(int jj=0; jj < iP; jj++)
			for(int kk = jj; kk < iP; kk++)
			{
				int nn = kk * (kk + 1) / 2 + jj;
				mTemp[jj][kk] = vCov[nn];
				mTemp[kk][jj] = vCov[nn];
			}
		*mCov = mTemp;
	}

	return 0;
}

/*
		Perform a polynomial fit on the input curve.
	Parameters:
		crvData=Input, the X and Y coordinates of the sample data
		mResultCurves=Output, the result curves where
			Col(0) is the X coordinates of the fit line
			Col(1) is the Y coordinates of the fit line
			Col(2) is the lower confidence band
			Col(3) is the upper confidence band
			Col(4) is the lower prediction band
			Col(5) is the upper prediction band
		trPR=Input settings & Output results of the linear regression
			trPR.GUI.Fit.PolynomialOrder.nVal			// The order of the polynomial
			trPR.GUI.Fit.ThroughZero.nVal				// Force fit to pass thru zero
			trPR.GUI.Fit.ErrBarWeight.nVal				// Use error column for wt - use (1/err^2) as wt factor
			trPR.GUI.Fit.UseReducedChiSq.nVal			// Scale parameter errors with reduced chisqr

			trPR.GUI.ResultCurves.Points.nVal			// Number of points in fit curve
			trPR.GUI.ResultCurves.ConfBands.nVal		// Create confidence bands - if not set, then matrix will have empty columns
			trPR.GUI.ResultCurves.PredBands.nVal		// Create prediction bands - if not set, then matrix will have empty columns
			trPR.GUI.ResultCurves.Confidence.dVal		// Confidence value to be used 

			trPR.Calculation.Control.UseDataXRange.nVal	// Option = 1 to use data range, = 0 to use X1, X2
			trPR.Calculation.Control.X1.dVal			// Default X minimum for fit curve
			trPR.Calculation.Control.X2.dVal			// Default X maximum for fit curve

			trPR.Calculation.Parameters.P1.Name.strVal	// Parameter Name...A
			trPR.Calculation.Parameters.P1.Value.dVal	// Parameter Value
			trPR.Calculation.Parameters.P1.Error.dVal	// Parameter Error
			trPR.Calculation.Parameters.P1.Vary.nVal	// Parameter Vary
			trPR.Calculation.Parameters.P1.tValue.dVal	// Parameter t-Value
			trPR.Calculation.Parameters.P1.LCI.dVal		// Parameter LCI
			trPR.Calculation.Parameters.P1.UCI.dVal		// Parameter UCI
			trPR.Calculation.Parameters.P1.Prob.dVal	// Parameter Prob > |t|

			trPR.Calculation.Parameters.P2.Name.strVal	// Parameter Name...B1
			trPR.Calculation.Parameters.P2.Value.dVal
			trPR.Calculation.Parameters.P2.Error.dVal
			trPR.Calculation.Parameters.P2.Vary.nVal
			trPR.Calculation.Parameters.P2.tValue.dVal
			trPR.Calculation.Parameters.P2.LCI.dVal
			trPR.Calculation.Parameters.P2.UCI.dVal
			trPR.Calculation.Parameters.P2.Prob.dVal

			etc. per order of Polynomial

			trPR.Calculation.Stat.Rvalue.dVal			// R-value
			trPR.Calculation.Stat.RSqCOD.dVal			// RSq-COD value
			trPR.Calculation.Stat.AdjRSq.dVal			// Adjusted RSq
			trPR.Calculation.Stat.RMSESD.dVal			// RMSE-Standard Dev.
			trPR.Calculation.Stat.N.nVal				// No. of points in fit
	
			trPR.Calculation.ANOVA.Model.DOF.nVal		// Model DOF 
			trPR.Calculation.ANOVA.Model.SSq.dVal		// Model SS
			trPR.Calculation.ANOVA.Model.MeanSq.dVal	// Model Mean Sq
			trPR.Calculation.ANOVA.Error.DOF.nVal		// Error DOF
			trPR.Calculation.ANOVA.Error.SSq.dVal		// Error SS
			trPR.Calculation.ANOVA.Error.MeanSq.dVal	// Error Mean Sq
			trPR.Calculation.ANOVA.Total.DOF.dVal		// Total degrees of freedom
			trPR.Calculation.ANOVA.Total.SSq.dVal		// Total SS
			trPR.Calculation.ANOVA.Fvalue.dVal			// F-value
			trPR.Calculation.ANOVA.Pvalue.dVal			// P-value
	Return:
		Returns ERROR_NO_ERROR on successful exit and an error code on failure.
		NAG error codes:		Negative 1 * NAG error code
		ERROR_INVALID_CURVE:	Invalid input curve
		ERROR_INVALID_TREENODE:	Invalid TreeNode
		ERROR_TO_FEW_PTS:		Not enough points in the curve
		ERROR_SETTINGS:			Polynomial order invalid
*/
int stat_polynomial_fit(const curvebase& crvData, matrix& mResultCurves, TreeNode& trPR)
{
	// Check the input arguments
	if( !crvData.IsValid() )
		return ERROR_INVALID_CURVE;

	if( !trPR.IsValid() )
		return ERROR_INVALID_TREENODE;

	TreeNode trFit = trPR.GUI.Fit;
	if( !trFit.IsValid() )
		return ERROR_INVALID_TREENODE;

	TreeNode trResultCurves = trPR.GUI.ResultCurves;
	if( !trResultCurves.IsValid() )
		return ERROR_INVALID_TREENODE;

	vector vDepData, vX, vWeight, vxWeight;
	if( !crvData.CopyData(vX, vDepData, vWeight, vxWeight) )
		return ERROR_INVALID_CURVE;

	// Get the number of points of the curve
	int nPts = vDepData.GetSize();
	if( nPts < SU_MIN_NPTS ) 
		return ERROR_TO_FEW_PTS;	

	// The weights
	vector* pvWeight = NULL;
	if( trFit.ErrBarWeight.nVal )
	{
		vWeight = 1 / vWeight^2;
		pvWeight = &vWeight;
	}

	// Get the associated independent variable
	matrix mIndepData;
	int nOrder = trFit.PolynomialOrder.nVal; // The order of the polynomial
	if( nOrder < 1 ) 
		return ERROR_SETTINGS;
	mIndepData.SetSize(nPts, nOrder);
	int nOrder1 = nOrder + 1;

	if( !crvData.HasX() )
		vX.Data(0, nPts - 1);

	vector vtemp(nPts);
	vtemp = 1;
	for( int ii = 0; ii < nOrder; ii++)
	{
		vtemp *= vX;
		mIndepData.SetColumn(vtemp, ii);
	}

	// The variance-covariance matrix
	matrix mCov;

	// Use the multiple regression function for polynomial fitting
	int nRet = stat_multiple_linear_regression(vDepData, mIndepData, trPR, pvWeight, &mCov);
	if( nRet )
		return nRet;

	// Results
	// Standard deviation
	double SD = trPR.Calculation.Stat.RMSESD.dVal;

	// t(1 - alpha / 2, dof) ------- alpha: significance level
	double tvalue = tTable(1.0 - (1.0 - trResultCurves.Confidence.dVal) / 2.0, trPR.Calculation.ANOVA.Error.DOF.nVal);

	// Use reduced chi^2?
	double dfactor;
	if( trFit.UseReducedChiSq.nVal )
		dfactor = SD;
	else
		dfactor = 1.0;

	// Summary of the fitting
	string strTrNode, strParamName;
	strParamName = PR_P0_NAME;
	TreeNode trPara;
	vector vPara(nOrder1);
	for( ii = 0; ii < nOrder1; ii++ )
	{
		strTrNode.Format(SU_PN, ii + 1);
		trPara = trPR.Calculation.Parameters.GetNode(strTrNode);
		trPara.Name.strVal = strParamName;
		vPara[ii] = trPara.Value.dVal;
		trPara.Error.dVal = trPara.Error.dVal / dfactor;
		trPara.LCI.dVal = vPara[ii] - trPara.Error.dVal * tvalue;
		trPara.UCI.dVal = vPara[ii] + trPara.Error.dVal * tvalue;
		strParamName.Format(PR_PN_NAME, ii + 1);
	}

	// Set the dimension of the result matrix
	int FitPoints = trResultCurves.Points.nVal;
	mResultCurves.SetSize(FitPoints, PR_N_OUTPUT_COLS);

	// Fitted line
	double Xmin, Xmax, Xinc;
	if( trPR.Calculation.Control.UseDataXRange.nVal )
		vX.GetMinMax(Xmin, Xmax);
	else
	{
		Xmin = trPR.Calculation.Control.X1.dVal;
		Xmax = trPR.Calculation.Control.X2.dVal;
	}
	Xinc = (Xmax - Xmin) / (FitPoints - 1);

	vector vFitX(FitPoints), vFitY(FitPoints);
	matrix mVar(FitPoints, nOrder1); // Independent variables include the constant column: transpose([1 1 1...1])

	vFitX[0] = Xmin;
	for( ii = 1; ii < FitPoints; ii++ )
		vFitX[ii] = vFitX[ii - 1] + Xinc;

	vFitY = trPR.Calculation.Parameters.P1.Value.dVal;
	vtemp.SetSize(FitPoints);
	vtemp = 1;
	int jj = 0;
	if( trFit.ThroughZero.nVal == 0 )
		mVar.SetColumn(vtemp, jj++);

	for(ii = 1; ii < nOrder1; ii++)
	{
		double dtemp = vPara[ii];
		vtemp *= vFitX;
		vFitY += vtemp * dtemp;
		mVar.SetColumn(vtemp, jj++);
	}

	mResultCurves.SetColumn(vFitX, PR_FITX_COL);
	mResultCurves.SetColumn(vFitY, PR_FITY_COL);

	// The standard error of each estimated value
	vector vSerr(FitPoints);
	vSerr = 0;
	int nCovSize = mCov.GetNumRows();
	for( ii = 0; ii < FitPoints; ii++)
	{
		for( jj = 0; jj < nCovSize; jj++)
		{
			for(int kk = 0; kk < nCovSize; kk++)
			{
				vSerr[ii] += mVar[ii][jj] * mVar[ii][kk] * mCov[jj][kk];
			}
		}
	}

	// Confidence bands
	if( trResultCurves.ConfBands.nVal )
	{
		vector vConf(FitPoints); 
		vConf = 0;

		vtemp = sqrt(vSerr) * tvalue;

		vConf = vFitY - vtemp;
		mResultCurves.SetColumn(vConf, PR_LCB_COL);

		vConf = vFitY + vtemp;	
		mResultCurves.SetColumn(vConf, PR_UCB_COL);
	}

	// Prediction bands
	if( trResultCurves.PredBands.nVal )
	{
		vector vPred(FitPoints); 
		vPred = 0;

		vtemp = sqrt(vSerr + SD^2) * tvalue;

		vPred = vFitY - vtemp;
		mResultCurves.SetColumn(vPred, PR_LPB_COL);

		vPred = vFitY + vtemp;	
		mResultCurves.SetColumn(vPred, PR_UPB_COL);
	}

	return ERROR_NO_ERROR;
}