/*------------------------------------------------------------------------------*
 * File Name: Stat_Utils.c														*
 * Creation: GJL 4/1/2003														*
 * Purpose: Source c file containing statistics LLOC functions.					*
 * Copyright (c) OriginLab Corp. 2003-2007										*
 * All Rights Reserved															*
 * 																				*
 * Modification Log:															*
 *------------------------------------------------------------------------------*/

////////////////////////////////////////////////////////////////////////////////////
// Included header files
//////////////////////////////////////////////////////////////////////////////////
//
// System includes
#include <Origin.h>				// Most Origin related classes
//#include <OC_const.h>			// OC Constants
//#include <Common.h>				// Basic types
//#include <Data.h>				// Vectors and matrices
//#include <Tree.h>				// Tree class
//#include <Tree_Utils.h>			// Tree utilities
//#include <NAG\nag_types.h>		// NAG structures
//#include <NAG\OCN_g01.h>		// NAG Basic Statistics functions
//#include <NAG\OCN_g02.h>		// NAG nag_regsn_mult_linear function
//#include <GetNBox.h>			// GetNBox
#include "Stat_Utils.h"			// SU constants, non-localized strings, and function prototypes

////////////////////////////////////////////////////////////////////////////////////

// Include definitions of all localized strings. $ causes CodeBuilder to look for 
// correct version of Local.h depending on Windows Regional settings.
#include "$Local.h"

////////////////////////////////////////////////////////////////////////////////////
// Function definitions
//////////////////////////////////////////////////////////////////////////////////
//
/**
	Function to compute summary statistics for a vector passing results back in a tree.
*/
bool stat_descriptive(const vector& vData, TreeNode& trWhole, int nOffset, const vector* pvWeight) // nOffset = 0, pvWeight = NULL
{
	if( trWhole )
	{
		TreeNode tr = trWhole.Settings;		
		TreeNode trOther = trWhole.OtherSettings;

		// *** Use NAG to compute summary statistics as needed ***
		if( tr.N.Enable			||
			tr.Sum.Enable		||
			tr.Mean.Enable		||
			tr.SD.Enable		||
			tr.Kurtosis.Enable	||
			tr.Skewness.Enable	||
			tr.SE.Enable		||
			tr.Variance.Enable	||
			tr.CoefVar.Enable	||
			tr.CIL.Enable		||
			tr.CIU.Enable		||
			tr.Min.Enable		||
			tr.Max.Enable		||
			tr.Range.Enable )
		{
			SumStats ss;
			if( stat_summary(ss, vData, pvWeight, trOther.ConfLevel.dVal) )
			{
				if( tr.N.Enable )
					tr.N.dVal = (double) ss.N;
				if( tr.Sum.Enable )
					tr.Sum.dVal = ss.dSum;
				if( tr.Mean.Enable )
					tr.Mean.dVal = ss.dMean;
				if( tr.SD.Enable )
					tr.SD.dVal = ss.dSD;
				if( tr.Kurtosis.Enable )
					tr.Kurtosis.dVal = ss.dKurtosis;
				if( tr.Skewness.Enable )
					tr.Skewness.dVal = ss.dSkewness;
				if( tr.SE.Enable )
					tr.SE.dVal = ss.dSE;
				if( tr.Variance.Enable )
					tr.Variance.dVal = ss.dVariance;
				if( tr.CoefVar.Enable )
					tr.CoefVar.dVal = ss.dCoefVar;
				if( tr.CIL.Enable )
					tr.CIL.dVal = ss.dCIL;
				if( tr.CIU.Enable )
					tr.CIU.dVal = ss.dCIU;
				if( tr.Min.Enable )
					tr.Min.dVal = ss.dMin;
				if( tr.Max.Enable )
					tr.Max.dVal = ss.dMax;
				if( tr.Range.Enable )
					tr.Range.dVal = ss.dRange;
				tr.Missing.dVal = (double) ss.nMissing;
			}
			else
			{
				if( tr.N.Enable )
					tr.N.dVal = (double) vData.GetSize();
				if( tr.Sum.Enable )
					tr.Sum.dVal = NANUM;
				if( tr.Mean.Enable )
					tr.Mean.dVal = NANUM;
				if( tr.SD.Enable )
					tr.SD.dVal = NANUM;
				if( tr.Kurtosis.Enable )
					tr.Kurtosis.dVal = NANUM;
				if( tr.Skewness.Enable )
					tr.Skewness.dVal = NANUM;
				if( tr.SE.Enable )
					tr.SE.dVal = NANUM;
				if( tr.Variance.Enable )
					tr.Variance.dVal = NANUM;
				if( tr.CoefVar.Enable )
					tr.CoefVar.dVal = NANUM;
				if( tr.CIL.Enable )
					tr.CIL.dVal = NANUM;
				if( tr.CIU.Enable )
					tr.CIU.dVal = NANUM;
				if( tr.Min.Enable )
					tr.Min.dVal = NANUM;
				if( tr.Max.Enable )
					tr.Max.dVal = NANUM;
				if( tr.Range.Enable )
					tr.Range.dVal = NANUM;
				vector vTemp;
				vTemp = vData;
				vTemp.Trim();
				tr.Missing.dVal = (double) (vData.GetSize() - vTemp.GetSize());
			}
		}

		// *** Use vectorbase GetMinMax to get iMin and iMax as needed *** 
		if( tr.iMin.Enable || tr.iMax.Enable )
		{
			double dMin, dMax;
			uint iMin, iMax;
			if( vData.GetMinMax(dMin, dMax, &iMin, &iMax) )
			{
				if( tr.iMin.Enable )
					tr.iMin.dVal = c_index_to_labtalk_index(iMin, nOffset);
				if( tr.iMax.Enable )
					tr.iMax.dVal = c_index_to_labtalk_index(iMax, nOffset);
			}
			else
			{
				if( tr.iMin.Enable )
					tr.iMin.dVal = NANUM;
				if( tr.iMax.Enable )
					tr.iMax.dVal = NANUM;
			}
		}

		// *** Use Percentiles function to compute quartiles and percentiles as needed ***
		int nPercents, ii;
		vector vPercents;

		ii = 0;
		vPercents.SetSize(ii);

		if( tr.Q25.Enable || tr.IQR.Enable )
		{
			vPercents.SetSize(++ii);
			vPercents[ii - 1] = 25.0;
		}

		if(	tr.Median.Enable )
		{
			vPercents.SetSize(++ii);
			vPercents[ii - 1] = 50.0;
		}

		if( tr.Q75.Enable || tr.IQR.Enable )
		{
			vPercents.SetSize(++ii);
			vPercents[ii - 1] = 75.0;
		}

		if(	trOther.Percentile.Enable )
		{
			foreach( TreeNode trN in trOther.Percents.Children )
				vPercents.Add( trN.dVal );
		}

		nPercents = vPercents.GetSize();
		if( nPercents )
		{
			vector vPercentiles;
			vPercentiles.SetSize(nPercents);
			vPercentiles = NANUM;
			if( stat_percentiles(vPercentiles, vData, vPercents, trOther.Interpolate.nVal) )
			{
				ii = 0;
				if( tr.Q25.Enable )
					tr.Q25.dVal = vPercentiles[ii];
				if( tr.Q25.Enable || tr.IQR.Enable )
					ii++;
				if( tr.Median.Enable )
				{
					tr.Median.dVal = vPercentiles[ii];
					ii++;
				}
				if( tr.Q75.Enable )
					tr.Q75.dVal = vPercentiles[ii];
				if( tr.Q75.Enable || tr.IQR.Enable )
					ii++;
				if( tr.IQR.Enable )
					tr.IQR.dVal = tr.Q75.dVal - tr.Q25.dVal;
				if( trOther.Percentile.Enable )
				{
					TreeNode trSubNode;
					foreach( TreeNode trN in trOther.Percents.Children )
					{
						trSubNode = tr.GetNode(trN.tagName);
						trSubNode.dVal = vPercentiles[ii++];
					}
				}
			}
			else
			{
				if( tr.Q25.Enable )
					tr.Q25.dVal = NANUM;
				if( tr.Median.Enable )
					tr.Median.dVal = NANUM;
				if( tr.Q75.Enable )
					tr.Q75.dVal = NANUM;
				if( tr.IQR.Enable )
					tr.IQR.dVal = NANUM;
				if( trOther.Percentile.Enable )
				{
					TreeNode trSubNode;
					foreach( TreeNode trN in trOther.Percents.Children )
					{
						trSubNode = tr.GetNode(trN.tagName);
						trSubNode.dVal = NANUM;
					}
				}
			}
		}
	}

	return true;
}

/**
	Function to compute summary statistics for a vector passing results back in a SumStats structure.
*/
bool stat_summary(SumStats& ss, const vector& vData, const vector* pvWeight, double dConfLevel) // pvWeight = NULL, dConfLevel = SU_DEFAULT_CONF_LEVEL
{
	bool bRet;
	int ii, iRet;
	uint nPts;
	double dValue, dVary, dN, dDOF, dSigLevel, dHalfWidth;
	NagError neErr;

	nPts = vData.GetSize(); // Number of data points in vector
	
	if( nPts < SU_MIN_NPTS ) // NAG basic stats requires minimum number of points (2)
	{
		ss.N = nPts;
		ss.nMissing = 0;
		ss.dWeightSum = NANUM;

		if( nPts ) // If at least 1 point
		{
			dValue = vData[0];
			if( NANUM == dValue )
			{
				ss.nMissing = 1;
				dVary = NANUM;
			}
			else
				dVary = 0;
			//if( pvWeight )
				//ss.dWeightSum = *pvWeight; // GJLFix #4509 Can not dereference elements of vector using pointer to vector
		}
		else
		{
			dValue = NANUM;
			dVary = NANUM;
		}

		ss.dSum = dValue;
		ss.dMean = dValue;
		ss.dSD = dVary;
		ss.dKurtosis = NANUM;
		ss.dSkewness = NANUM;
		ss.dSE = dVary;
		ss.dVariance = dVary;
		ss.dCoefVar = dVary;
		ss.dCIL = dValue;
		ss.dCIU = dValue;
		ss.dMin = dValue;
		ss.dMax = dValue;
		ss.dRange = dVary;
		bRet = true;
	}
	else
	{
		// Use NAG to compute core statistics
		if( pvWeight ) // If weighting vector is specified
			iRet = nag_summary_stats_1var(nPts, vData, *pvWeight, &ss.N, &ss.dMean, &ss.dSD,
				&ss.dSkewness, &ss.dKurtosis, &ss.dMin, &ss.dMax, &ss.dWeightSum);
		else // Else no weighting vector 
			iRet = nag_summary_stats_1var(nPts, vData, NULL, &ss.N, &ss.dMean, &ss.dSD,
				&ss.dSkewness, &ss.dKurtosis, &ss.dMin, &ss.dMax, &ss.dWeightSum);
	
		if( !iRet ) // If Nag function was successful...
		{
			dN = (double) ss.N;
			ss.dSum = ss.dMean * dN; // Reverse compute sum of data values
		
			ss.dVariance = ss.dSD * ss.dSD; // Compute Variance
		
			if( ss.dMean )
				ss.dCoefVar = ss.dSD / ss.dMean; // If mean exists and is not 0 Compute Coefficient of Variation
			else
				ss.dCoefVar = NANUM;
	
			ss.dRange = ss.dMax - ss.dMin; // Compute Range
		
			dDOF = dN - 1.0;
			dSigLevel = 1.0 - dConfLevel;
			ss.dSE = ss.dSD / sqrt(dN);
			dHalfWidth = nag_deviates_students_t(Nag_TwoTailSignif, dSigLevel, dDOF, &neErr);
			if( NE_NOERROR == neErr.code ) // Compute Upper and Lower Confidence Intervals
			{
				dHalfWidth *= ss.dSE;
				ss.dCIL = ss.dMean - dHalfWidth;
				ss.dCIU = ss.dMean + dHalfWidth;		
			}
			else
			{
				ss.dCIL = NANUM;
				ss.dCIU = NANUM;	
			}
	
			ss.nMissing = nPts - ss.N; // Count number of missing values = total number of points - number of valid points
	
			bRet = true;
		}
		else
			bRet = false;
	}

	return bRet;
}

/**
	Function to compute percentiles for a vector.
*/
bool stat_percentiles(vector& vPercentiles, const vector& vData, const vector& vPercents, bool bInterpolate) // bInterpolate = false
{
	bool bRet;
	uint ii, nPercents, nPts, iIndex;
	vector vTemp;
	double dNpts, dIndex, dFractionalPart;

	nPercents = vPercents.GetSize();
	vPercentiles.SetSize(nPercents); // Size Percentiles vector
	vPercentiles = NANUM; // Initialize Percentiles to NANUM

	if( nPercents ) // If there is at least one Percent then continue
	{
		vTemp = vData; // Copy input data
		vTemp.Trim(); // Remove missing values from vector
		nPts = vTemp.GetSize(); // Get number of non-missing values
		dNpts = (double) nPts; // Cast to double for computation
		if( nPts > 0 ) // If at least one data point
		{
			vTemp.Sort(); // Sort temporary vector in ascending order

			for( ii = 0; ii < nPercents; ii++ )
			{
				if( vPercents[ii] < 0. || vPercents[ii] > 100. ) 
					vPercentiles[ii] = NANUM; // If Percent is less than 0 or greater than 100 Percentile is NANUM
				else
				{
					// *** This algorithm for computing Percentiles is based on the ***
					// *** LabTalk Percentile function in MoStatVc.cpp              ***
					dIndex = dNpts * vPercents[ii] / 100.0; // Absolute index

					iIndex = (uint) dIndex; // Integral part of index
					dFractionalPart = dIndex - (double) iIndex; // Fractional part of index

					if( iIndex == 0 ) // If iIndex is 0 define Percentile = vTemp[0]
						vPercentiles[ii] = vTemp[0];
					else if( iIndex == nPts ) // Else if iIndex is nPts define Percentile = vTemp[nPts - 1]
						vPercentiles[ii] = vTemp[nPts - 1];
					else // Else will need to compute Percentile 
					{
						if( bInterpolate ) // If interpolation is specified
						{
							if( dFractionalPart == 0 ) // If fractional part is 0 (exact hit)
							{
								if( vPercents[ii] == 100 )
									vPercentiles[ii] = vTemp[iIndex - 1]; // If Percent is 100 then no interpolation
								else
									vPercentiles[ii] = ( vTemp[iIndex - 1] + vTemp[iIndex] ) / 2.0; // Else interpolate
							}
							else // Else take next larger item
								vPercentiles[ii] = vTemp[iIndex];
						}
						else
						{
							if( dFractionalPart == 0 ) // If fractional part is 0 (exact hit)
								vPercentiles[ii] = vTemp[iIndex - 1];
							else // Else take next larger item
								vPercentiles[ii] = vTemp[iIndex];
						}
					}
				}
			}
			bRet = true; // Successful computation of Percentiles
		}
		else // Else no data so return error 
			bRet = false;
	}
	else // Else no Percents so return error
		bRet = false;

	return bRet; // Return true on success and false on failure
}

/**
		Function to perform multiple linear regression using the NAG function nag_regsn_mult_linear.
	Parameters:
		nPts=Input number of rows (or data points) in the input matrix
		nTdx=Input number of columns (or independent variables) in the input matrix
		mX=Input matrix containing data points of the independent variables
		vY=Input vector containing data points of dependent variable
		vWT=Input vector containing weights of data points, NULL if weighting not used
		dRss=Output residual sum of squares for the regression 
		dDf=Output degrees of freedom associated with dRss
		vB=Output least-squares estimates of the parameters of the regression model
		vSE=Output standard errors of the parameter estimates given in vB
		vCOV=Output variance-covariance matrix of estimated parameters in vB
		vRES=Output weighted residuals
		vH=Output diagonal elements of H
		mQ=Output results of the QR decomposition
		bSvd=Output flag set TRUE if singular value decomposition was performed, otherwise FALSE
		iRank=Output rank of the independent variables
		vP=Output details of QR decomposition if SVD is used
		vCOMAR=Output information which is needed by nag_regsn_mult_linear_newyvar_if bSvd is TRUE 
		dTol=Input tolerance used to decide rank of independent variables, default is MR_DEFAULT_TOLERANCE 
		bThroughZero=Input flag forcing regression line through origin, default is FALSE
	Return:
		Returns the results of the NAG nag_regsn_mult_linear function.
*/
int stat_multiple_linear_regression( int nPts, int nTdx, matrix& mX, vector& vY, vector& vWT,
	double& dRss, double& dDf, vector& vB, vector& vSE, vector& vCOV, vector& vRES,
	vector& vH, matrix& mQ, BOOL& bSvd, int& iRank, vector& vP,	vector& vCOMAR,
	double dTol, BOOL bThroughZero)
{
	// Declare local variable holding error code
	int iErr;

	// *** Declare, size, prepare, and check all NAG Linear Regression arguments (as needed) ***
	// Force regression line through zero or not...
	int iP;
	Nag_IncludeMean iMean;						// NAG variable type from header file \OriginC\system\NAG\nag_types.h
	if( bThroughZero )							// If input flag says to force through zero...
	{
		iMean = Nag_MeanZero;					// Set NAG constant to force regression line through zero
		iP = nTdx;								// Set number P of independent variables in the model not including the mean (or intercept)		
	}
	else
	{
		iMean = Nag_MeanInclude;				// Set NAG constant to NOT force regression line through zero
		iP = nTdx + 1;							// Set number P of independent variables in the model including the mean
	}

	// NAG nag_regsn_mult_linear function requires nPts > 1
	if( nPts < SU_MIN_NPTS )					// If less than 2 points return error
		return ERROR_TO_FEW_PTS;

	// The number of independent variables equals the second dimension of the matrix
	int nM = nTdx;
	
	// Indicates which of the potential independent variables are to be included in the model...use all so set to 1
	vector<int> vSX;
	vSX.SetSize( nM );
	vSX = 1;
	
	// The least-squares estimates of the parameters of the regression model
	vB.SetSize( iP );
	
	// The standard errors of the iP parameter estimates given in B
	vSE.SetSize( iP );
	
	// The variance-covariance matrix of estimated parameters in B
	vCOV.SetSize( iP * ( iP + 1 ) / 2);
	
	// The (weighted) residuals
	vRES.SetSize( nPts );
	
	// The diagonal elements of H
	vH.SetSize( nPts );
	
	// The second dimension of the array Q
	int nTdq;
	nTdq = iP + 1;
	
	// The results of the QR decomposition
	mQ.SetSize( nPts, nTdq );
	
	// Details of the QR decomposition and SVD if used
	vP.SetSize( 2 * iP + iP * iP );
	
	// Information which is needed by nag_regsn_mult_linear_newyvar_if bSvd is TRUE
	vCOMAR.SetSize( 5 * ( iP - 1 ) + iP * iP );
	
	// *** Call NAG function nag_regsn_mult_linear to perform linear regression ***
	// From <NAG\OCN_g02.h>: g02dac nag_regsn_mult_linear: Perform multiple linear regression
	iErr = nag_regsn_mult_linear( iMean, nPts, mX, nTdx, nM, vSX, iP, vY, vWT, &dRss, &dDf,
		vB, vSE, vCOV, vRES, vH, mQ, nTdq, &bSvd, &iRank, vP, dTol, vCOMAR );

	return iErr;
}

/*
		Perform a linear regression 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
		trLR=Input settings & Output results of the linear regression
			trLR.GUI.Fit.ThroughZero.nVal				// Force fit to pass thru zero
			trLR.GUI.Fit.FixSlope.nVal					// force slope to be fixed
			trLR.GUI.Fit.FixSlopeAt.dVal				// fixed value of slope (Note: cannot fix slope and thru zero at same time, should return error from LLOC)
			trLR.GUI.Fit.ErrBarWeight.nVal				// Use error column for wt - use (1/err^2) as wt factor
			trLR.GUI.Fit.UseReducedChiSq.nVal			// Scale parameter errors with reduced chisqr

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

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

			trLR.Calculation.Parameters.P1.Name.strVal	// Parameter Name...A
			trLR.Calculation.Parameters.P1.Value.dVal	// Parameter Value
			trLR.Calculation.Parameters.P1.Error.dVal	// Parameter Error
			trLR.Calculation.Parameters.P1.Vary.nVal	// Parameter Vary
			trLR.Calculation.Parameters.P1.tValue.dVal	// Parameter t-Value
			trLR.Calculation.Parameters.P1.LCI.dVal		// Parameter LCI
			trLR.Calculation.Parameters.P1.UCI.dVal		// Parameter UCI
			trLR.Calculation.Parameters.P1.Prob.dVal	// Parameter Prob > |t|
	
			trLR.Calculation.Parameters.P2.Name.strVal	// Parameter Name...B
			trLR.Calculation.Parameters.P2.Value.dVal
			trLR.Calculation.Parameters.P2.Error.dVal
			trLR.Calculation.Parameters.P2.Vary.nVal
			trLR.Calculation.Parameters.P2.tValue.dVal
			trLR.Calculation.Parameters.P2.LCI.dVal
			trLR.Calculation.Parameters.P2.UCI.dVal
			trLR.Calculation.Parameters.P2.Prob.dVal
	
			trLR.Calculation.Stat.Rvalue.dVal			// R-value
			trLR.Calculation.Stat.RSqCOD.dVal			// RSq-COD value
			trLR.Calculation.Stat.AdjRSq.dVal			// Adjusted RSq
			trLR.Calculation.Stat.RMSESD.dVal			// RMSE-Standard Dev.
			trLR.Calculation.Stat.N.nVal				// No. of points in fit
	
			trLR.Calculation.ANOVA.Model.DOF.nVal		// Model DOF 
			trLR.Calculation.ANOVA.Model.SSq.dVal		// Model SS
			trLR.Calculation.ANOVA.Model.MeanSq.dVal	// Model Mean Sq
			trLR.Calculation.ANOVA.Error.DOF.nVal		// Error DOF
			trLR.Calculation.ANOVA.Error.SSq.dVal		// Error SS
			trLR.Calculation.ANOVA.Error.MeanSq.dVal	// Error Mean Sq
			trLR.Calculation.ANOVA.Total.DOF.dVal		// Total degrees of freedom
			trLR.Calculation.ANOVA.Total.SSq.dVal		// Total SS
			trLR.Calculation.ANOVA.Fvalue.dVal			// F-value
			trLR.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_X_RANGE:			All the elements of X are equal
		ERROR_SETTINGS:			'ThroughZero' and 'FixSlope' can not be both enabled
*/
int stat_linear_fit(const curvebase& crvData, matrix& mResultCurves, TreeNode& trLR) 
{
	int nRet = ERROR_NO_ERROR;

	// Check the input arguments
	if( !crvData.IsValid() )
		return ERROR_INVALID_CURVE;

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

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

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

	if( trFit.ThroughZero.nVal && trFit.FixSlope.nVal )
		return ERROR_SETTINGS;

	double dIntercept, dSlope;	// The intercept and slope of the fitted line
	double da_Err, db_Err;		// The standard error of the intercept and slope
	double rsq, rss, dof;		// The R-square, residual sum of squares, and the degrees of freedom
	int mdof;					// The degrees of freedom of the Model

	double dSumWeight;			// The sum of weights
	double dSxx, dSxy, dSyy;
	double dxAve, dyAve;		// The average value of x and y

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

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

	if( trFit.FixSlope.nVal == 1 ) // Fix slope
	{
		trLR.Calculation.Parameters.P1.Vary.nVal = 1;
		trLR.Calculation.Parameters.P2.Vary.nVal = 0;

		if( 0 != stat_linear_regrssion_sum(vX, vY, vWeight, dSumWeight, dxAve, dyAve, dSxx, dSxy, dSyy) )
			return ERROR_X_RANGE;

		// Slope
		dSlope = trFit.FixSlopeAt.dVal;
		db_Err = 0;

		// Intercept		
		dIntercept = dyAve - dSlope * dxAve;

		// Residual sum of squares
		rss = dSyy - 2 * dSlope * dSxy + dSlope^2 * dSxx;

		// Degrees of freedom
		dof = nPts - 1;
		mdof = 0;

		// Standard error of intercept
		da_Err = sqrt(rss / dof / dSumWeight);

		// R-square
		rsq = dSxy^2 / dSxx / dSyy;
	}
	else // Slope is not fixed
	{
		mdof = 1;
		Nag_SumSquare mean; // Whether or not force the line through zero
		trLR.Calculation.Parameters.P2.Vary.nVal = 1;

		if( trFit.ThroughZero.nVal == 0 ) // Without constant 
		{
			trLR.Calculation.Parameters.P1.Vary.nVal = 1;
			mean = Nag_AboutMean;
			dof = nPts - 2;
		}
		else // With constant
		{
			trLR.Calculation.Parameters.P1.Vary.nVal = 0;
			mean = Nag_AboutZero;
			dof = nPts - 1;
		}
	
		// Linear regression by NAG function
		double df;
		nRet = nag_simple_linear_regression(mean, nPts, vX, vY, vWeight, &dIntercept,
			&dSlope, &da_Err, &db_Err, &rsq, &rss,&df);
		if( nRet != 0 )
			return -1 * nRet;

		// Get other statistics sums
		if( 0 != stat_linear_regrssion_sum(vX, vY, vWeight, dSumWeight, dxAve, dyAve, dSxx, dSxy, dSyy) )
			return ERROR_X_RANGE;

		// DOF_RELATED_DIFFERENCE
		if( trFit.ErrBarWeight.nVal )
		{
			da_Err = da_Err * sqrt(df / dof);
			db_Err = db_Err * sqrt(df / dof);
			rsq = dSxy^2 / dSxx / dSyy;
		}
	}
	
	// Set the dimension of the result matrix
	int FitPoints = trResultCurves.Points.nVal;
	mResultCurves.SetSize(FitPoints, FITLR_NUM_COLS);
	
	// Fit line
	double Xmin, Xmax, Xinc;
	if( trLR.Calculation.Control.UseDataXRange.nVal > 0 )
	{
		vX.GetMinMax(Xmin, Xmax);
	}
	else
	{
		Xmin = trLR.Calculation.Control.X1.dVal;
		Xmax = trLR.Calculation.Control.X2.dVal;
	}
	Xinc = (Xmax - Xmin) / (FitPoints - 1);
	
	vector vFitX(FitPoints), vFitY(FitPoints);
	vFitX[0] = Xmin;
	for( int ii = 1; ii < FitPoints; ii++)
	{
		vFitX[ii] = vFitX[ii - 1] + Xinc;
	}
	vFitY = vFitX * dSlope + dIntercept;
	
	mResultCurves.SetColumn(vFitX, FITLR_X);
	mResultCurves.SetColumn(vFitY, FITLR_Y);
	
	// Standard deviation
	double SD = sqrt(rss / dof);
	
	// t(1-alpha/2,dof) ------- alpha: significance level
	double tvalue = tTable(1 - (1.0 - trResultCurves.Confidence.dVal) / 2, dof);
	
	// Confidence bands
	if( trResultCurves.ConfBands.nVal )
	{
		vector vConf(FitPoints), vSerr(FitPoints);
		
		if( trFit.ThroughZero.nVal )
		{
			vSerr = fabs(vFitX) * db_Err;
		}
		else
		{
			vSerr = sqrt(1.0 / dSumWeight + (vFitX - dxAve)^2 / dSxx) * SD;
		}

		vConf = vFitY - tvalue * vSerr;
		mResultCurves.SetColumn(vConf, FITLR_UCB);
		vConf = vFitY + tvalue * vSerr;	
		mResultCurves.SetColumn(vConf, FITLR_UCB);
	}

	// Prediction bands
	if( trResultCurves.PredBands.nVal )
	{
		vector vPre(FitPoints), vSerr(FitPoints);
		vSerr = sqrt(1.0 + 1.0 / dSumWeight + (vFitX - dxAve)^2 / dSxx) * SD;
		vPre=vFitY - tvalue * vSerr;
		mResultCurves.SetColumn(vPre, FITLR_LPB);
		vPre=vFitY + tvalue * vSerr;	
		mResultCurves.SetColumn(vPre, FITLR_UPB);
	}

	// Simple statistics results
	if( trFit.ErrBarWeight.nVal && !trFit.UseReducedChiSq.nVal ) 
	{
		double dtemp = sqrt(rss / dof);
		da_Err /= dtemp;
		db_Err /= dtemp;
	}
	
	string str = LR_COL_NAMES;
	// Intercept
	trLR.Calculation.Parameters.P1.Name.strVal = str.GetToken(0, SU_DELIM);
	trLR.Calculation.Parameters.P1.Value.dVal = dIntercept;
	trLR.Calculation.Parameters.P1.Error.dVal = da_Err; 
	trLR.Calculation.Parameters.P1.tValue.dVal = dIntercept / da_Err;
	trLR.Calculation.Parameters.P1.LCI.dVal = dIntercept - da_Err * tvalue;
	trLR.Calculation.Parameters.P1.UCI.dVal = dIntercept + da_Err * tvalue;
	trLR.Calculation.Parameters.P1.Prob.dVal = 2 * (1 - invt(fabs(trLR.Calculation.Parameters.P1.tValue.dVal), dof));
	
	// Slope
	trLR.Calculation.Parameters.P2.Name.strVal = str.GetToken(1, SU_DELIM);
	trLR.Calculation.Parameters.P2.Value.dVal = dSlope;
	trLR.Calculation.Parameters.P2.Error.dVal = db_Err;
	trLR.Calculation.Parameters.P2.tValue.dVal = dSlope / db_Err;
	trLR.Calculation.Parameters.P2.LCI.dVal = dSlope - db_Err * tvalue;
	trLR.Calculation.Parameters.P2.UCI.dVal = dSlope + db_Err * tvalue;
	trLR.Calculation.Parameters.P2.Prob.dVal = 2 * (1 - invt(fabs(trLR.Calculation.Parameters.P2.tValue.dVal), dof));
	
	// Stat - Correlation coefficient
	trLR.Calculation.Stat.Rvalue.dVal = sqrt(rsq);
	trLR.Calculation.Stat.RSqCOD.dVal = rsq;
	trLR.Calculation.Stat.AdjRSq.dVal = 1 - (1 - rsq) * (dof + 1) / dof;
	trLR.Calculation.Stat.RMSESD.dVal = SD;
	trLR.Calculation.Stat.N.nVal = nPts;
	
	// ANOVA
	trLR.Calculation.ANOVA.Error.DOF.nVal = dof;
	trLR.Calculation.ANOVA.Error.SSq.dVal = rss;
	trLR.Calculation.ANOVA.Error.MeanSq.dVal = rss /dof;
	trLR.Calculation.ANOVA.Total.DOF.nVal = dof + mdof;
	trLR.Calculation.ANOVA.Total.SSq.dVal = dSyy;
	trLR.Calculation.ANOVA.Model.DOF.nVal = mdof;
	trLR.Calculation.ANOVA.Model.SSq.dVal = trLR.Calculation.ANOVA.Total.SSq.dVal - trLR.Calculation.ANOVA.Error.SSq.dVal;
	trLR.Calculation.ANOVA.Model.MeanSq.dVal = trLR.Calculation.ANOVA.Model.SSq.dVal;
	trLR.Calculation.ANOVA.Fvalue.dVal = (dSyy - rss) / rss * dof;
	trLR.Calculation.ANOVA.Pvalue.dVal = 1 - invf(trLR.Calculation.ANOVA.Fvalue.dVal, 1, dof);

	return nRet;
}

int stat_linear_regrssion_sum(const vectorbase& vX, const vectorbase& vY, const vectorbase& vW,
	double& dSumWeight, double& dxAve, double& dyAve, double& dSxx, double& dSxy, double& dSyy)
{
	vector vtemp;	

	vW.Sum(dSumWeight);

	// Average value of y
	vtemp = vY * vW;
	vtemp.Sum(dyAve);
	dyAve /= dSumWeight;

	// Syy
	vtemp *= vY;
	vtemp.Sum(dSyy);
	dSyy -= dyAve^2 * dSumWeight;

	// Average value of x
	vtemp = vX * vW;
	vtemp.Sum(dxAve);
	dxAve /= dSumWeight;

	// Sxx
	vtemp *= vX;
	vtemp.Sum(dSxx);
	dSxx -= dxAve^2 * dSumWeight;
	if( dSxx < 1e-15 )
		return -1;

	// Sxy
	vtemp = vW * vY * vX;
	vtemp.Sum(dSxy);
	dSxy -= dxAve * dSumWeight * dyAve;

	return 0;
}

/*
		Perform a multiple linear regression using the NAG function nag_regsn_mult_linear.
	Parameters:
		vDepData=Input, vector holding input dependent data
		mIndepData=Input, matrix holding input independent data
			trMR: Input settings and Output results, tree having the nodes:
			trMR.GUI.Fit.ThroughZero.nVal					// Force fit to pass thru zero
	
			trMR.Calculation.Parameters.P1.Name.strVal		// Parameter Name...Y-Intercept
			trMR.Calculation.Parameters.P1.Value.dVal		// Parameter Value
			trMR.Calculation.Parameters.P1.Error.dVal		// Parameter Error
			trMR.Calculation.Parameters.P1.tValue.dVal		// Parameter t-Value
			trMR.Calculation.Parameters.P1.Prob.dVal		// Parameter Prob > |t|
	
			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;
}
