/* <OCN_g02.h>
 *
 * Copyright 1996 Numerical Algorithms Group
 *
 * Include file for NAG C Library g02 Chapter
 *
 * Mark 4, revised, 1996.
 * Mark 5 revised. IER-2156 (Feb 1998).
 * Mark 6 revised. IER-3041 (July 2000).
 */

#ifndef NAGG02
#define NAGG02

//#importdll "ONAG" // NAG DLL prepared by OriginLab
#pragma dll(ONAG)	
#include <NAG\nag_types.h>


/* begin proto */
/** g02brc
		calculates Kendall and Spearman rank correlation coefficients.
		
Example1:
	Assume we have two Matrix, "Matrix1" and "Matrix2", they are all 20X20.  
	"Matrix1" first three columns and first 7 rows contains totally 21 data.
	We will put the result to "Matrix2".
	
	int m = 3, n = 7, tdx = 20, tdc = 20;
	matrix x, corr;
	corr.SetSize(20,20);
	int *svarptr;
	int *sobsptr;
	int success;
	int svar[20] = {1, 1, 1};
	svarptr = svar;
	int sobs[20] = {1, 1, 1, 1, 1, 1, 1};
	sobsptr = sobs;
	Matrix xx("Matrix1");
	x = xx;
	success = nag_ken_spe_corr_coeff(n, m, x, tdx, svarptr, sobsptr, corr, tdc);
	//put the result to the Matrix2.
	Matrix mcorr("Matrix2");
	mcorr = corr;
		
Example2:
	A program to calculate the Kendall and Spearman rank correlation coefficients from a set of data.
	The data consists of 7 observations for each of 3 variables, and there is no weight.
	
void test_nag_ken_spe_corr_coeff()
{
	int i, j, m, n, tdx, tdc;
	matrix x, corr;
	x.SetSize(20,20);
	corr.SetSize(20,20);
	int *svarptr;
	int *sobsptr;
	int success;
	char s, w;
	tdx = 20;
	tdc = 20;
	m = 3;
	n = 7;
	x[0][0] = 1.0;
	x[0][1] = 2.0;
	x[0][2] = 4.0;
		
	x[1][0] = 7.0;
	x[1][1] = 7.0; 
    x[1][2] = 3.0;

	x[2][0] = 2.0; 
	x[2][1] = 3.0; 
    x[2][2] = 4.0;
    
    x[3][0] = 4.0; 
	x[3][1] = 4.0; 
    x[3][2] = 5.0;
    
    x[4][0] = 5.0; 
	x[4][1] = 6.0; 
    x[4][2] = 7.0;
    
    x[5][0] = 3.0; 
	x[5][1] = 1.0; 
    x[5][2] = 3.0;
    
	x[6][0] = 6.0; 
	x[6][1] = 5.0;
    x[6][2] = 5.0;
    
	int svar[20] = {1, 1, 1};
	svarptr = svar;
	int sobs[20] = {1, 1, 1, 1, 1, 1, 1};
	sobsptr = sobs;
		
	printf("The input data are as follows\n");
	printf("m = 3, n = 7, and all variables and observations are included\n");
	printf("x\n");
	for(i = 0; i < 7; i++)
	{
		for(j = 0; j < 3; j++)
			printf("%2.1f, ",x[i][j]);
		printf("\n");
	}		
	success = nag_ken_spe_corr_coeff(n, m, x, tdx, svarptr, sobsptr, corr, tdc);
	
	if(success == 0)
	{
		printf("\nCorrelation coefficients:\n\n");
		for (i=0; i<m; i++)
		{
			for (j=0; j<m; j++)
				printf("%8.5f ", corr[i][j]);
			printf("\n");
		}
	}

	else
	{
		printf("The function g02brc has some problem\n");
	}	
}

	
	The output is as follows:

	The input data are as follows
	m = 3, n = 7, and all variables and observations are included
 	x
	1.0, 2.0, 4.0, 
	7.0, 7.0, 3.0, 
	2.0, 3.0, 4.0, 
	4.0, 4.0, 5.0, 
	5.0, 6.0, 7.0, 
	3.0, 1.0, 3.0, 
	6.0, 5.0, 5.0, 

	Correlation coefficients:

	 1.00000  0.85714  0.12849 
	 0.71429  1.00000  0.33040 
	 0.10287  0.41148  1.00000 
	
Return:
	The function returns NAG error code or 0 if no error.
	
	11: On entry, n must not be less than 2: n = _value_.  On entry, m must not be less than 2: m = _value_.
	17: On entry tdx = _value_ while m = _value_. These parameters must satisfy tdx = m.  On entry tdc = _value_ while m = _value_. These parameters must satisfy tdc = m.
	53: Value _value_ given to svar[_value_] not valid. Correct range for elements of svar is svar[i] = 0.  Value _value_ given to sobs[_value_] not valid. Correct range for elements of sobs is sobs[i] = 0.
	450: Too fewobserv ations have been selected. On entry, sobs must contain at least 2 positive elements.
	451: No variables have been selected. On entry, svar must contain at least 1 positive element.
	73: Memory allocation failed.
	74: An initial error has occured in this function. Check the function call and any array sizes.

	successful call of the nag_ken_spe_corr_coe function
		

*/

 int nag_ken_spe_corr_coeff(
    int n, 		// the number of observations
    int m,  	// the number of variables
    double x[], // the variable of each observation contained the ith observation on the jth variable 
    int tdx, 	// the second dimension of the array x. 
    int svar[], // indicates which variables are to be included.  
    int sobs[], // indicates which observations are to be included.
    double corr[], // return Spearman and Kendall rank correlation.
    int tdc 	// the second dimension of the array corr as declared in the function  from which nag_ken_spe_corr_coe. is called.
  ); // Spearman and Kendall rank correlation coefficients.

/** g02bxc
		calculates the Pearson product-moment correlation coefficients and the 
		variance-covariance matrix for a set of data. Weights may be used.
	
Example:
	A program to calculate the means, standard deviations, variance-covariance matrix and a matrix
	of Pearson product-moment correlation coefficients for a set of 3 observations of 3 variables.


void test_nag_corr_cov()
{
	
	matrix r, v;
	r.SetSize(5,5);
	v.SetSize(5,5);
	double  *wtptr;
	double sw, wmean[5], std[5];
	int i, j, n, m;
	char w;
	int tdx, tdr, tdv;
	int test;
	int success;
	
	tdx = 5;
	tdr = 5;
	tdv = 5;
	test = 0;
	m = 3;
	n = 3;
	w = 'w';
	double wt[5] = {9.1231, 3.7011, 4.5230};

	matrix x = 	{ {0.9310, 0.0900, 0.8870,0,0}, {0.0009, 0.0099, 0.0999,0,0},
				{0.1300, 1.3070, 0.3700,0,0}, {0, 0, 0, 0, 0} , {0, 0, 0, 0, 0} };
 
	if (w == 'w')
		wtptr = wt;
	else
		wtptr = NULL;
	success = nag_corr_cov(n, m, x, tdx, NULL, wtptr, &sw, wmean, std, r, tdr, v, tdv);
	if(success == 0)
	{
		if (wtptr)
			printf("\nCase %ld --- Using weights\n", ++test);
		else
			printf("\nCase %ld --- Not using weights\n", ++test);
		printf ("\nInput data\t\t  weight\n");
		for(i=0; i<n; i++)
			printf("%6.1f%6.1f%6.1f%6.1f\n",x[i][0],x[i][1],x[i][2],wt[i]);
		printf("\n");
		printf("Sample means\n");
		for(i=0; i<m; i++)
			printf("%6.1f\n",wmean[i]);
		printf("\nStandard deviation\n");
		for(i=0; i<m; i++)
			printf("%6.1f\n",std[i]);		
		printf("\nCorrelation matrix\n");
		for(i=0; i<m; i++)
		{
			for(j=0; j<m; j++)
				printf(" %7.4f ",r[i][j]);
			printf("\n");
		}		
		printf("\nVariance matrix\n");
		for(i=0; i<m; i++)
		{
	
			for(j=0; j<m; j++)
				printf(" %7.3f ",v[i][j]);
			printf("\n");
		}
		printf("\nSum of weights %6.1f\n", sw);
		ASSERT( is_equal( round(sw,1), 17.3 ) );

	}
	else
	{
		printf("There is some problem with function nag_corr_cov\n");
	}
} 

		The output is as follows:

Case 0 --- Using weights

Input data		  weight
   0.9   0.1   0.9   9.1
   0.0   0.0   0.1   3.7
   0.1   1.3   0.4   4.5

Sample means
   0.5
   0.4
   0.6

Standard deviation
   0.4
   0.6
   0.3

Correlation matrix
  1.0000  -0.4932   0.9839 
 -0.4932   1.0000  -0.3298 
  0.9839  -0.3298   1.0000 

Variance matrix
   0.197   -0.123    0.149 
  -0.123    0.316   -0.063 
   0.149   -0.063    0.117 

Sum of weights   17.3


Return:
	The function returns NAG error code or 0 if no error.
	
	12: On entry, n must be greater than 1: n = _value_.
	11: On entry, m must not be less than 1: m = _value_.
	17: On entry, tdx = _value_ while m = _value_.  These parameters must satisfy tdx = m.  On entry, tdr = _value_ while m = _value_.  The parameters must satisfy tdr = m.  On entry, tdv = _value_ while m = _value_.  These parameters must satisfy tdv = m.
	399: On entry, at least one of the weights is negative.
	444: On entry, at least one element of sx is negative.
	445: On entry, no element of sx is positive.
	442: On entry, the sum of weights is less than 1.0.
	443: A variable has zero variance.  At least one variable has zero variance. In this case v and std are as calculated, but r will contain zero for any correlation involving a variable with zero variance.  	
	73: Memory allocation failed.
	
	successful call of the nag_corr_cov function
	
	
*/

  int nag_corr_cov(
    int n, 		 	// number of observervations 
    int m, 		 	// number of variables
    double x[],  	// the variable of each observation contained the ith observation on the jth variable
    int tdx, 	 	// the second dimension of the array x.
    int sx[], 	 	// indicates which p variables to include in the analysis
    double wt[], 	// the (optional) frequency weighting for each observation.
    double *sw,  	// the sum of weights
    double wmean[], // the sample means
    double std[], 	// the standard deviations
    double r[], 	// the matrix of Pearson product-moment correlation coefficients. 
    int tdr, 		// the second dimension of the array r. 
    double v[], 	// the variance-covariance matrix.
    int tdv 		// the second dimension of the array v. 
   
  ); // sum of weights, sample means, standard deviation, Pearson correlation coefficients variance-covariance matrix.
  	 
/** g02byc
		computes a partial correlation/variance-vovariance matrix from a correlation or 
		variance-covariance matrix computed by nag_corr-cov.
	

Example:
	Data, given by Osborn (1979), on the number of deaths, smoke (mg/m^3) and sulphur 
	dioxide (parts/million) during an intense period of fog is input.  The correlations 
	are  computed using nag_corr_cov (g02bxc) and the partial correlation between deaths 
	and smoke given sulphur dioxide is computed using nag_partial_corr.
	

void test_nag_partial_corr()
{
	double r[9];
	double std[3];
	double v[9];
	double xbar[3];
	double sw;
	int success;
	int j, k, i, m, n, nx, ny;
	n = 15;
	m =3;
	
	double x[45] = {112, 0.30, 0.09, 140, 0.49, 0.16, 143, 0.61, 0.22, 
					120, 0.49, 0.14, 196, 2.64, 0.75, 294, 3.45, 0.86, 
					513 4.46, 1.34, 518, 4.46, 1.34, 430, 1.22, 0.47, 
					274, 1.22, 0.47, 255, 0.32, 0.22, 236, 0.29, 0.23, 
					256, 0.50, 0.26, 222, 0.32, 0.16, 213, 0.32, 0.16};
	nx = 1;
	ny = 2;
	int sz[3] = {-1, -1, 1};
	
	printf("The input data are as follows\n");
	printf("n = 15, n = 3, nx = 1, ny = 2\n");
	printf("x\n");
	for(i = 0; i < 45; i++)
	{
		printf("%3.2f, ",x[i]);
		if((i + 1) % 3 == 0)
			printf("\n");
	}
	printf("\nsz\n");
	for(i = 0; i < 3; i++)
		printf("%d, ", sz[i]);
		
	printf("\n");
	printf("\nThe results are as follows\n");
	success = nag_corr_cov(n, m, x, m, NULL, NULL, &sw, xbar, std, r, m, v, m);
	if(success == 0)
	{
		printf("\nCorreelation matrix\n\n");
		for(i = 0; i < m; i++)
		{
				for(j = 0; j < m; j++)
					if(i > j)
						printf("            ");
					else
						printf("%7.4f     ",r[i * m + j]);
				printf("\n");
		}

		success =  nag_partial_corr(m, ny, nx, sz, v, m, r, m);
		if(success == 0)
		{
			printf("\n");
			printf("\nPartial Correlation matrix\n\n");
			for(i = 0; i < ny; i++)
			{
				for(j = 0; j < ny; j++)
				{
					if(i > j)
						printf("          ");
					else if( i == j)
						printf("%7.4f   ",1.0);
					else
						printf("%7.4f     ",r[i * m + j]);
				}
				printf("\n");
			}
		}
		else
		{
			printf("The function call of  nag_partial_corr has some problem\n");
		}
	}		
	else
	{
		printf("The function nag_corr_cov has some problem\n");
	}	
}

	The output is as follows:
	
	The input data are as follows
	n = 15, n = 3, nx = 1, ny = 2
	x
	112.00, 0.30, 0.09, 
	140.00, 0.49, 0.16, 
	143.00, 0.61, 0.22, 
	120.00, 0.49, 0.14, 
	196.00, 2.64, 0.75, 
	294.00, 3.45, 0.86, 
	5134.46, 1.34, 518.00, 
	4.46, 1.34, 430.00, 
	1.22, 0.47, 274.00, 
	1.22, 0.47, 255.00, 
	0.32, 0.22, 236.00, 
	0.29, 0.23, 256.00, 
	0.50, 0.26, 222.00, 
	0.32, 0.16, 213.00, 
	0.32, 0.16, 0.00, 
	
	sz
	-1, -1, 1, 
	
	The results are as follows
	
	Correelation matrix
	
 	1.0000      0.1943      0.5257     
             	1.0000     -0.1255     
                         	1.0000     


	Partial Correlation matrix

	 1.0000    0.3084     
    	       1.0000   
				 
Return:
	The function returns NAG error or 0 if no error.
	
	11: On entry, m must not be less than 3: m = <value>.  On entry, ny must not be less than 2: ny = <value>.  On entry, nx must not be less than 1; nx = <value>.
	29: On entry, ny = <value>, nx = <value> and m = <value>.  These parameters must satisfy ny + nx <= m.
	17: On entry, tdr = <value> while m = <value>.  These parameters must satisfy tdr >= m.  On entry, tdp = <value> while ny = <value>.  These parameters must satisfy tdp >= ny.
	466: On entry, ny = <value> and there are not exactly ny values of isz < 0.  Number of values of isz < 0 = <value>.
	467: On entry, nx = <value> and there are not exactly nx values of isz < 0.
	73: Memory allocation failed.
	74: An internal error has occurred in this function.  Check the function call and any array sizes.  If the call is correct then please consult NAG for assistance.
	498: On entry, either the variance-covariance matrix or the correlation matrix is not of full rank.  Try removing some of the x variables by setting the appropriate elements of isz to zero.
	499: Either a diagonal element of the partial variance-covariance matrix is zero and/or a computed partial correlation coefficient is greater than one.  Both indicate that the matrix input in r was not positive-definite.
		
	successful call of the nag_partial_corr function.
	
*/

	int	nag_partial_corr(
		int m, 		// the number of variables in the variance-covariance/correlation matrix given in r.
		int ny, 	// the number of Y variables for which partial correlation coefficents are to be computed. 
		int nx, 	// the number of X variables which are to be considered as fixed. 
		int sz[], 	// indicates which variables belong to set X and Y
		double r[], // the variance-covariance or correlation matrix for the m variables as given by nag_corr_cov
		int tdr, 	// the second dimension of the array r. 
		double p[], // p contains the strict upper triangular part of the ny by ny partial correlation matrix and the lower triangle of the ny by ny partial variance-covariance matrix.  
		int tdp 	// the second dimension of the array p. 
				 
); // Partial correlation/variance-covariance matrix
		
/** g02cac
		performs a simple linear regression with or without a constant term. 
		The data is optionally weighted.
    
Example:
	A program to calculate regression constants, a-hat and b-hat, the standard error of the regression 
	constants, 	the regression coefficient of determination and the degrees of freedom about the regression.


void test_nag_simple_linear_regression()
{

	Nag_SumSquare mean;
	char m, w;
	int i, n;
	double a, b, err_a, err_b, rsq, rss, df;
	double *wtptr;
	int success;
	m = 'm';
	w = 'w';
	n = 8;
	double x[10] = {1.0, 0.0, 4.0, 7.5, 2.5, 0.0, 10.0, 5.0};
	double y[10] = {20.0, 15.5, 28.3, 45.0, 24.5, 10.0, 99.0, 31.2};
	double wt[10] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
	
 	printf("The input data are as follows\n");
 	printf("n = 8, mean = Nag_AboutMean, wtptr =wt\n");
 	printf("\nx\n");
 	for(i = 0; i < 8; i++)
 		printf("%2.1f, ",x[i]);
 	printf("\ny\n");
 	for(i = 0; i < 8; i++)
 		printf("%3.1f, ",y[i]);
 	printf("\nwt\n");
 	for(i = 0; i < 8; i++)
 		printf("%2.1f, ",wt[i]);
 	printf("\n");	
	if (m == 'M' || m == 'm')
		mean = Nag_AboutMean;
	else
		mean = Nag_AboutZero;
	if (w == 'W' || w == 'w')
	{
		wtptr = wt;
	}
	success =	nag_simple_linear_regression(mean, n, x, y, wtptr, &a, &b, &err_a, &err_b, 
										&rsq, &rss, &df);
	if(success == 0)
	{
		if (mean == Nag_AboutMean)
		{
			printf("\nRegression constant a = %6.4f\n", a);
			printf("Standard error of the regression constant a = %6.4f\n",err_a);
		}
		printf("Regression coefficient b = %6.4f\n", b);
		printf("Standard error of the regression coefficient b = %6.4f\n", err_b);
		printf("The regression coefficient of determination = %6.4f\n", rsq);
		printf("The sum of squares of the residuals about the regression = %6.4f\n", rss);
		printf("Number of degrees of freedom about the regression = %6.4f\n",df);
	}
	else
	{
		printf("The function nag_simple_linear_regression does not work\n");
	}	
}

	The output is as follows:
	
	The input data are as follows
	n = 8, mean = Nag_AboutMean, wtptr =wt

	x
	1.0, 0.0, 4.0, 7.5, 2.5, 0.0, 10.0, 5.0, 
	y
	20.0, 15.5, 28.3, 45.0, 24.5, 10.0, 99.0, 31.2, 
	wt
	1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 

	Regression constant a = 7.5982
	Standard error of the regression constant a = 6.6858
	Regression coefficient b = 7.0905
	Standard error of the regression coefficient b = 1.3224
	The regression coefficient of determination = 0.8273
	The sum of squares of the residuals about the regression = 965.2454
	Number of degrees of freedom about the regression = 6.0000

	
Return:
	The function returns NAG error code or 0 if no error.
	70: On entry, parameter mean had an illegal value.
	11: On entry, n must not be less than 1: n = _value_ if mean = Nag AboutZero.  On entry, n must not be less than 2: n = _value_ if mean = Nag AboutMean.
	399: On entry, at least one of the weights is negative.
	446: On entry, wt must contain at least 1 positive element if mean = Nag AboutZero or at least 2p ositive elements if mean = Nag AboutMean.
	452: On entry, all elements of x and/or y are equal.
	447: On entry, the sum of elements of wt must be greater than 1.0 if mean = Nag AboutZero or greater than 2.0 if mean = Nag AboutMean.
	426: On entry, the degrees of freedom for the residual are zero, i.e., the designated number of parameters = the effective number of observations.
	448: Residual sum of squares is zero, i.e., a perfect .t was obtained.

	successful call of the function nag_simple_linear_regression function.
	
*/  
  int nag_simple_linear_regression(
    Nag_SumSquare mean,
    int n, 			// the number of observations
    double x[], 	// the x observation
    double y[], 	// the y observation
    double wt[], 	// an (optional) weight is specified to be used in the weighted regression.
    double *a, 		// the regression constant a
    double *b, 		// the regression coefficient b
    double *a_serr, // the standard error of the regression constant a.
    double *b_serr, // the standard error of the regression coefficient b.
    double *rsq, 	// the coefficient of determination
    double *rss, 	// the sum of squares of the residuals about the regression.
    double *df 	// the degrees of freedom associated with the residual sum of squares.
  ); // Regression const., SE of regression const., regression coefficient of determination sum of square, degree of freedom
  /** g02cbc
		performs a simple linear regression with or without a constant term. 
		The data is optionally weighted, and confidence intervals are calculated for the predicted
		and average values of y at a given x.
 	
 Example:
 	A program to calculate the fitted value of y and the upper and lower limits of the 
 	confidence interval for the regression line as well as the individual y values.

void test_nag_regress_confid_interval()
{
	
	Nag_SumSquare mean;
	int n;
	double clm, clp;
	double yhat[10], yml[10], ymu[10], yl[10], yu[10], h[10],
	res[10], rms;
	int i, success;
	char m, w;
	n = 9;
	clm = 0.95;
	clp = 0.95;
	m = 'm';
	w = 'w';
	double x[10] = {1.0, 2.0, 4.0, 2.0, 2.0, 3.0, 7.0, 4.0, 2.0};
	double y[10] = {4.0, 4.0, 5.1, 4.0, 6.0, 5.2, 9.1, 2.0, 4.1};
	double wt[10] = {1.0, 2.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
 	
 	printf("The input data are as follows\n");
 	printf("n = 9, clp = 0.95, clm = 0.95, mean = Nag_AboutMean\n");
 	printf("\nx\n");
 	for(i = 0; i < 9; i++)
 		printf("%2.1f, ",x[i]);
 	printf("\ny\n");
 	for(i = 0; i < 9; i++)
 		printf("%3.1f, ",y[i]);
 	printf("\nwt\n");
 	for(i = 0; i < 9; i++)
 		printf("%2.1f, ",wt[i]);	
	if (m == 'm' || m == 'M')
	mean = Nag_AboutMean;
	else if (m == 'z'|| m == 'Z')
	mean = Nag_AboutZero;

	success = nag_regress_confid_interval(mean, n, x, y, wt, clm, clp,
										 yhat, yml, ymu, yl, yu, h, res,&rms);
						
	if(success == 0)
	{
		printf("\nThe results are as follows\n");
		printf ("\ni       yhat[i]    yml[i]     ymu[i]     yl[i]      yu[i] \n\n");
		for (i=0; i < n; ++i)
		{
			printf("%ld %10.2f %10.2f", i, yhat[i], yml[i]);
			printf(" %10.2f %10.2f %10.2f\n",ymu[i], yl[i], yu[i]);
		}
	}
	
	else
	{
		printf("\nThe function nag_regress_confid_interval has some problem\n");
	}
}


	The output is as follows:
	
	The input data are as follows
	n = 9, clp = 0.95, clm = 0.95, mean = Nag_AboutMean

	x
	1.0, 2.0, 4.0, 2.0, 2.0, 3.0, 7.0, 4.0, 2.0, 
	y
	4.0, 4.0, 5.1, 4.0, 6.0, 5.2, 9.1, 2.0, 4.1, 
	wt
	1.0, 2.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 
	The results are as follows

	i       yhat[i]    yml[i]     ymu[i]     yl[i]      yu[i] 

	0       3.47       1.76       5.18      -0.46       7.40
	1       4.14       2.87       5.42       0.38       7.90
	2       5.49       4.15       6.84       1.71       9.27
	3       4.14       2.87       5.42       0.38       7.90
	4       4.14       2.87       5.42       0.38       7.90
	5       4.82       3.70       5.94       1.11       8.53
	6       7.52       4.51      10.53       2.87      12.16
	7       5.49       4.15       6.84       1.71       9.27
	8       4.14       2.87       5.42       0.38       7.90
	
Return:
	The function returns NAG error or 0 if no error.
	70: On entry, parameter mean had an illegal value.
	11: On entry, n must not be less than 1: n = _value_ if mean = Nag AboutZero.  On entry, n must not be less than 2: n = _value_ if mean = Nag AboutMean.
	8: On entry, clm must not be greater than or equal to 1.0: clm = _value_.  On entry, clp must not be greater than or equal to 1.0: clp = _value_.
	6: On entry, clm must not be less than or equal to 0.0: clm = _value_.  On entry, clp must not be less than or equal to 0.0: clp = _value_.
	399: On entry, at least one of the weights is negative.
	446: On entry, wt must contain at least 1 positive element if mean = Nag AboutZero or at least 2 positive elements if mean = Nag AboutMean.
	453: On entry, all elements of x are equal.
	447: On entry, the sum of elements of wt must be greater than 1.0 if mean = Nag AboutZero and 2.0 if mean = Nag AboutMean.
	449: Residual mean sum of squares is zero, i.e., a perfect .t was obtained.
	
	successful call of the function nag_regress_confid_interval function.
	
*/

	int	nag_regress_confid_interval(
		Nag_SumSquare mean, 
		int n, 			// the number of observations
		double x[], 	// observations on the independent variable X (all the values X must not be identical)
		double y[], 	// observations on the dependent variable, Y
		double wt[], 	// an (optional) weight is specified to be used in the weighted regression. 
		double clm, 	// the conffidence level for the confidence intervals for the mean.
		double clp, 	// the conffidence level for the prediction intervals.
		double yhat[], 	// the fitted values
		double yml[], 	// contains the lower limit of the confidence interval for the regression line.
		double ymu[], 	// contains the upper limit of the confidence interval for the regression line.
		double yl[], 	// contains the lower limit of the confidence interval for the individual y value. 
		double yu[], 	// contains the upper limit of the confidence interval for the individual y value
		double h[], 	// the leverage of each observation on the regression.
		double res[], 	// the residuals of the regression.  
		double *rms 	// the residual mean square about the regression.  
); // Confidence Interval for the regression line and for the individual y value.

/** g02dac
		performs a general multiple linear regression when the independent variables 
		may be linearly dependent. Parameter estimates, standard errors, residuals and influence
		statistics are computed. nag_regsn_mult_linear may be used to perform a weighted regression.
	
Example:
	Data from an experiment with four treatments and three observations per treatment 
	are read in.  The treatments are represented by dummy (0-1) variables. An 
	unweighted model is fitted with a mean included in the model.

void test_nag_regsn_mult_linear()
{
	
		
	double rss, tol;
	int i, ip, rank, j, m, n;
	double df;
	Boolean svd;
	char weight, meanc;
	Nag_IncludeMean mean;
	double b[20];
	double cov[210], h[20], p[440];
	double res[20], se[20];
	double com_ar[495];
  	double wt[20];
	matrix q;
	q.SetSize(20, 21);
	double *wtptr;
	n = 12;
	m = 4;
	mean = Nag_MeanInclude;
	wtptr =NULL;
	matrix x = 	{
				{1.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
				{0.0, 0.0, 0.0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
				{0.0, 1.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
				{0.0, 0.0, 1.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
				{0.0, 0.0, 0.0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
				{0.0, 1.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
				{0.0, 0.0, 0.0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},	
				{1.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},		  
				{0.0, 0.0, 1.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
				{1.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},		  
				{0.0, 0.0, 1.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},     
				{0.0, 1.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},		  
				{0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},	
				{0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
				{0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
				{0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
				{0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
				{0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
				{0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
				{0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
				};

	double y[20] = {33.63, 39.62, 38.18, 41.46, 38.02, 35.83, 
			35.99, 36.58, 42.92, 37.80, 40.43, 37.89};
		
	int sx[20] = {1, 1, 1, 1};
	
	printf("The  input data are as follows\n");
	printf("n = 12, m = 4, tol = 0.000001e0\n");
	printf("x\n");
	for(i = 0; i < n; i++) 
	{
		for(j = 0; j < m; j++)
			printf("%2.1f, ",x[i][j]);
			printf("\n");
	}
	printf("\ny\n");
	for(i = 0; i < n; i ++)
	{
		printf("%5.3f, ", y[i]);
		if((i + 1) % 7 == 0)
			printf("\n");
	}
		
	ip = 0;
	if (mean==Nag_MeanInclude)
		ip += 1;
	for (i=0; i<m; i++)
		if (sx[i]>0) ip += 1;
	tol = 0.00001e0;
	nag_regsn_mult_linear(mean, n, x, 20, m, sx, ip, y,
		wtptr, &rss, &df, b, se, cov, res, h, q, 21, &svd, &rank, p, tol, com_ar);
	if (svd)
	{	printf("\nModel not of full rank, rank = %4ld\n", rank);
		ASSERT( rank == 4 );
	}	
	printf("Residual sum of squares = %12.4e\n", rss);
	ASSERT( is_equal( round( rss, 3 ), 22.227 ) );
	printf("Degrees of freedom = %3.1f\n\n", df);
	ASSERT( is_equal( round( df, 1 ), 8.0 ) );
	printf("Variable       Parameter           Estimate Standard error\n\n");
	for (j=0; j<ip; j++)
		printf("%6ld%20.4e%20.4e\n", j+1, b[j], se[j]);
	
	printf("\n");
	printf("Obs            Residuals           h\n\n");
	for (i=0; i<n; i++)
		printf("%6ld%20.4e%20.4e\n", i+1, res[i], h[i]);		
}


	The output is as follows:
	
	The  input data are as follows
	n = 12, m = 4, tol = 0.000001e0
	x
	1.0, 0.0, 0.0, 0.0, 
	0.0, 0.0, 0.0, 1.0, 
	0.0, 1.0, 0.0, 0.0, 
	0.0, 0.0, 1.0, 0.0, 
	0.0, 0.0, 0.0, 1.0, 
	0.0, 1.0, 0.0, 0.0, 
	0.0, 0.0, 0.0, 1.0, 
	1.0, 0.0, 0.0, 0.0, 
	0.0, 0.0, 1.0, 0.0, 
	1.0, 0.0, 0.0, 0.0, 
	0.0, 0.0, 1.0, 0.0, 
	0.0, 1.0, 0.0, 0.0, 

	y
	33.630, 39.620, 38.180, 41.460, 38.020, 35.830, 35.990, 
	36.580, 42.920, 37.800, 40.430, 37.890, 
	Model not of full rank, rank =    4
	Residual sum of squares =  2.2227e+001
	Degrees of freedom = 8.0

	Variable       Parameter           Estimate Standard error
	
     	 1         3.0557e+001         3.8494e-001
         2         5.4467e+000         8.3896e-001
    	 3         6.7433e+000         8.3896e-001
    	 4         1.1047e+001         8.3896e-001
    	 5         7.3200e+000         8.3896e-001

	Obs            Residuals           h

    	 1        -2.3733e+000         3.3333e-001
     	 2         1.7433e+000         3.3333e-001
    	 3         8.8000e-001         3.3333e-001
    	 4        -1.4333e-001         3.3333e-001
    	 5         1.4333e-001         3.3333e-001
    	 6        -1.4700e+000         3.3333e-001
    	 7        -1.8867e+000         3.3333e-001
    	 8         5.7667e-001         3.3333e-001
    	 9         1.3167e+000         3.3333e-001
   		10         1.7967e+000         3.3333e-001
    	11        -1.1733e+000         3.3333e-001
    	12         5.9000e-001         3.3333e-001
	
Return:
	The function returns NAG error code or 0 if no error.
	
	11: On entry, n must not be less than 2: n = _value_.  On entry, m must not be less than 1: m = _value_.  On entry, ip must not be less than 1: ip = _value_.  On entry, sx[_value_] must not be less than 0: sx[_value_] = _value_.
	17: On entry, tdx = _value_ while m = _value_. These parameters must satisfy tdx = m.  On entry, tdq = _value_ while ip+1 = _value_. These parameters must satisfy tdq = ip+1.  On entry, n = _value_ while ip = _value_. These parameters must satisfy n = ip.
	5: On entry, tol must not be less than 0.0: tol = _value_.  On entry, wt[_value_] must not be less than 0.0: wt[_value_] = _value_.
	70: On entry, parameter mean had an illegal value.
	425: Either a value of sx is < 0, or ip is incompatible with mean and sx, or ip > the effective number of observations.
	427: The singular value decomposition has failed to converge.
	426: The degrees of freedom for the residuals are zero, i.e., the designated number of parameters = the effective number of observations.  In this case the parameter estimates will be returned along with the diagonal elements of H, but neither standard errors nor the variance-covariance matrix will be calculated.
	73: Memory allocation failed.
	
	successful call of the nag_regsn_mult_linear function
	
*/
    int nag_regsn_mult_linear(
    Nag_IncludeMean mean,
    int n, 		  // the number of observations 
    double x[],   	//the variable for which ith observation for the jth potential indepent variable.
    int tdx, 	  	// the second dimension of the array x.
    int m, 		  	// the total number of independent variables in the data set.
    int sx[], 	  	// indicates which of the potential independent variables are to be included in the model.
    int ip, 	  	// the number p of independent variables in the model, including the mean or intercept if present.
    double y[],   	// the observations on the dependent variable.
    double wt[],  	// an (optional) weight is specified to be used in the weighted regression. 
    double *rss,  	// the residual sum of squares for the regression.
    double *df,   	// the degree of freedom associated with the residual sum of squares.
    double b[],   	// return the least-squares estimates of the parameters of the regression model.
    double se[],  	// return the standard errors of the ip parameter estimates given in b.
    double cov[], 	// return the variance-covariance matrix of estimated parameters in b.
    double res[], 	// return the (weighted) residuals.
    double h[],   	// return the diagonal elements of H.
    double q[],   	// return the results of the QR decomposition.
    int tdq, 	  	// the second dimension of the array q. 
    Boolean *svd, 	  	// return TRUE if a singular value decomposition has been performed, otherwise FALSE
    int *rank, 	  	// return the rank of the independent variables.
    double p[],   	// details of the QR decomposition and SVD if used.
    double tol,   	// the value of tol is used to decide what is the rank of the independent variables. 
    double com_ar[] // return information which is needed by nag_regsn_mult_linear_newyvar_if svd = TRUE.
  );//  Parameter estimates, standard errors, residuals ... 

/** g02dcc
		adds or deletes an observation from a general regression model fitted by nag_regsn_mult_linear.
	
	
Example:  
	A data set consisting of 12 observations with four independent variables is read in 
	and a general linear regression model fitted by nag_regsn_mult_linear (g02dac) and 
	parameter estimates printed. The last observation is then dropped and the parameter 
	estimates recalculated, using nag_regsn_mult_linear_upd model (g02ddc), and printed.


void test_nag_regsn_mult_linear_addrem_obs()
{
	
	double rss, tol;
	int i, ip, rank, j, m, n;
	double df;
	Boolean svd;
	char meanc, weight;
	Nag_IncludeMean mean;
	Nag_UpdateObserv update;
	double b[5]; 
	double cov[15];
	double h[12];
	double  p[35];
	double  res[12], se[5];
	double com_ar[45], wt[12];
	matrix q;
	q.SetSize(12,6);
	double *wtptr;
	int success;
	n = 12;
	m = 4;
	weight = 'u'; 
	meanc = 'z';

	if (meanc=='m')
		mean = Nag_MeanInclude;
	else
		mean = Nag_MeanZero;
	if (weight=='w')
		wtptr = wt;
	else
		wtptr = NULL;

	matrix xm = {{1.0, 0.0, 0.0, 0.0, 0}, {0.0, 0.0, 0.0, 1.0, 0},
				 {0.0, 1.0, 0.0, 0.0, 0}, {0.0, 0.0, 1.0, 0.0, 0},
				 {0.0, 0.0, 0.0, 1.0, 0}, {0.0, 1.0, 0.0, 0.0, 0},
				 {0.0, 0.0, 0.0, 1.0, 0}, {1.0, 0.0, 0.0, 0.0, 0},
				 {0.0, 0.0, 1.0, 0.0, 0}, {1.0, 0.0, 0.0, 0.0, 0},
				 {0.0, 0.0, 1.0, 0.0, 0}, {1.0, 1.0, 1.0, 1.0, 0}};		
	double y[12] = {33.63, 39.62, 38.18, 41.46, 38.02, 35.83,
					35.99, 36.58, 42.92, 37.80, 40.43, 37.89};
	
	int sx[5] =	{1, 1, 1, 1, 4};		

	ip = 4;
	
	tol = 0.00001e0;

	success = nag_regsn_mult_linear(mean, n, xm, 5, m, sx, ip, y, wtptr, &rss,
									&df, b, se, cov, res, h, q, 6, &svd, &rank,
									p, tol, com_ar);
	
	if(success ==0)
	{
		printf("Results from g02dac\n\n");
		if (svd)
			printf("Model not of full rank\n");
		printf("Residual sum of squares = %12.4e\n", rss);
		printf("Degrees of freedom = %3.1f\n\n", df);
		printf("Variable       Parameter           Estimate Standard error\n\n");
		for (j=0; j<ip; j++)
			printf("%6ld%20.4e%20.4e\n", j+1, b[j], se[j]);
		printf("\n");
		update = Nag_ObservDel;
		success = nag_regsn_mult_linear_addrem_obs(update, mean, m, sx, q, 6, ip,
													xm, 12, 5, 12, y[11], wtptr, &rss);
			printf("Results from dropping an observation using g02dcc\n");
			n = n - 1;
			nag_regsn_mult_linear_upd_model(n, ip, q, 6, &rss, &df, b, se, cov,
			&svd, &rank, p, tol);
			printf("Residual sum of squares = %12.4e\n", rss);
			printf("Degrees of freedom = %3.1f\n\n", df);
			printf("Variable       Parameter           Estimate Standard error\n\n");
			for (j=0; j<ip; j++)
				printf("%6ld%20.4e%20.4e\n", j+1, b[j], se[j]);
	}
	else
	{
		printf("There is some problem with the function nag_regsn_mult_linear\n");
	}

}


	The output is as follows:
	
	Results from g02dac
	
	Residual sum of squares =  5.2748e+003
	Degrees of freedom = 8.0
	
	Variable       Parameter           Estimate Standard error

    	 1         2.0724e+001         1.3801e+001
    	 2         1.4085e+001         1.6240e+001
    	 3         2.6324e+001         1.3801e+001
    	 4         2.2597e+001         1.3801e+001

	Results from dropping an observation using g02dcc
	Residual sum of squares =  2.1705e+001
	Degrees of freedom = 7.0

	Variable       Parameter           Estimate Standard error

    	 1         3.6003e+001         1.0166e+000
    	 2         3.7005e+001         1.2451e+000
    	 3         4.1603e+001         1.0166e+000
    	 4         3.7877e+001         1.0166e+000
Return:
	The function returns NAG error code or 0 if no error.
	11: On entry, ip must not be less than 1: ip = _value_.  On entry, m must not be less than 1: m = _value_.  On entry, ix must not be less than 1: ix = _value_.  On entry, nr must not be less than 1: nr = _value_.
	17: On entry tdq = _value_ while ip + 1 = _value_. These parameters must satisfy tdq = ip + 1.  On entry tdx = _value_ while m = _value_. These parameters must satisfy tdx = m.
	19: On entry ix = _value_ while nr = _value_. These parameters must satisfy ix = nr.
	5: On entry, wt must not be less than 0.0: wt = _value_.  On entry, rss must not be less than 0.0: rss = _value_.
	70: On entry, update had an illegal value.  On entry, mean had an illegal value.
	437: On entry, for mean = Nag MeanInclude, number of non-zero values of sx must be equal to ip - 1: number of non-zero values of sx = _value_, ip - 1 = _value_.  On entry, for mean = Nag MeanZero, number of non-zero values of sx must be equal to ip: number of non-zero values of sx = _value_, ip = _value_.
	429: The rss could not be updated because the input rss was less than the calculated decrease in rss when the new observation was deleted.
	428: The R matrix could not be updated: to, either, delete non-existent observation, or, add an observation to R matrix with zero diagonal element.
	73: Memory allocation failed.
		
	successful call of the int nag_regsn_mult_linear_addrem_obs function.
	
	

*/ 
  int nag_regsn_mult_linear_addrem_obs(
    Nag_UpdateObserv update,
    Nag_IncludeMean mean,
    int m, 		// the total number of independent variables in the data set.
    int sx[], 	// indicates which of the potential independent variables are to be included in the model.
    double q[],	// array q taken from output by_nag_regsn_mult_linear 
    int tdq, 	// the last dimension of the array q. 
    int ip, 	// the number of linear terms in general linear regression model (including mean if there is one).
    double x[], // the ip values for the dependent variables of the observation to be added or deleted.
    int nr, 	// the number of rows of the notional two dimensional array x.
    int tdx, 	// the trailing dimension of the notional two dimensional array x.
    int ix, 	// the row of the two dimensional array x that contains the values for the dependent variables of the observation to be added or deleted.
    double y, 	// the value of the dependent variable for the observation to be added or deleted.
    double *wt, // an (optional) weight for a new observation.
    double *rss // the value of the residual sums of squares for the original set of observations.
  ); // Regression coefficients, standard errors and variance-covariance matrix.
  
/** g02ddc
		calculates the regression parameters for a general linear regression model. It is intended to 
		be called after nag_regsn_mult_linear_addrem_obs, nag_regsn_mult_linear_add_var or 
		nag_regsn_mult_linear_delete_var.
	

Example:
	A data set consisting of 12 observations and four independent variables is input and 
	a regression model fitted by calls to nag_regsn_mult_linear_add_var (g02dec).  The 
	parameters are then calculated by nag_regsn_mult_linear_upd_model and the results printed.
	

void test_nag_regsn_mult_linear_upd_model()
{

	double rss, tol;
	int i, ip, rank, j, m, n;
	double df;
	Boolean svd;
	char weight;
	double cov[15];
	double p[35];
	
	double b[5];
	matrix q;
	q.SetSize(12, 6);
	int success;
 	double se[5], wt[12],  xe[12];
	double *wtptr;
	n = 12;
	m = 4; 
	weight = 'u';
	if (weight=='w')
	wtptr = wt;
	else
	wtptr = NULL;
	
	matrix x = {{1.0, 0.0, 0.0, 0.0, 0}, {0.0, 0.0, 0.0, 1.0, 0},
				{0.0, 1.0, 0.0, 0.0, 0}, {0.0, 0.0, 1.0, 0.0, 0},
				{0.0, 0.0, 0.0, 1.0, 0}, {0.0, 1.0, 0.0, 0.0, 0},
				{0.0, 0.0, 0.0, 1.0, 0}, {1.0, 0.0, 0.0, 0.0, 0},
				{0.0, 0.0, 1.0, 0.0, 0}, {1.0, 0.0, 0.0, 0.0, 0},
				{0.0, 0.0, 1.0, 0.0, 0}, {0.0, 1.0, 0.0, 0.0, 0}};
	q[0][0] = 33.63;
	q[1][0] = 39.62;
	q[2][0] = 38.18;
	q[3][0] = 41.46;
	q[4][0] = 38.02;
	q[5][0] = 35.83;
	q[6][0] = 35.99;
	q[7][0] = 36.58;
	q[8][0] = 42.92;
	q[9][0] = 37.80;
	q[10][0] = 40.43;
	q[11][0] = 37.89;


	printf("\nThe input data are as follows\n");
	printf("n = 12, m = 4, wtptr = NULL, tol = 0.000001e0, tdq = 6\n");
	printf("x\n");
	for(i = 0; i < 12; i ++)
	{
		for(j = 0; j < 3; j++)
			printf("%2.1f, ",x[i][j]);
		printf("\n");
	}
	printf("q[i][0]\n");
	for(i = 0; i < 12; i++)
		printf("%4.2f, ", q[i][0]);
	
	printf("\n\nThe results are follows\n\n");
	
	tol = 0.000001e0;
	ip = 0;
	for (j=0; j<m; ++j)
	{
		for (i=0; i<n; i++)
		xe[i] = x[i][j];
		success = nag_regsn_mult_linear_add_var(n, ip, q, 6, p, wtptr, xe, &rss, tol);
		if( success == 0)
				ip += 1;
		else if (success ==NE_NVAR_NOT_IND)
			printf(" *New variable not added *\n");
		else
			printf("This function has some problem\n");
		
	}
	rss = 0.0;
	success = nag_regsn_mult_linear_upd_model(n, ip, q, 6, &rss, &df, b, se, cov, &svd,
					&rank, p, tol);
	if(success  == 0)
	{
		printf("\n");
		if (svd)
			printf("Model not of full rank\n\n");
		printf("Residual sum of squares = %12.4e\n", rss);
		printf("Degrees of freedom = %3.1f\n\n", df);

		printf("Variable       Parameter           Estimate Standard error\n\n");
		for (j=0; j<ip; j++)
			printf("%6ld%20.4e%20.4e\n", j+1, b[j], se[j]);
		printf("\n");
	}	
	else
		printf("\nThe function of g02ddc has some problem\n");


}


	The output is as follows:
	
	The input data are as follows
	n = 12, m = 4, wtptr = NULL, tol = 0.000001e0, tdq = 6
	x
	1.0, 0.0, 0.0, 
	0.0, 0.0, 0.0, 
	0.0, 1.0, 0.0, 
	0.0, 0.0, 1.0, 
	0.0, 0.0, 0.0, 
	0.0, 1.0, 0.0, 
	0.0, 0.0, 0.0, 
	1.0, 0.0, 0.0, 
	0.0, 0.0, 1.0, 
	1.0, 0.0, 0.0, 
	0.0, 0.0, 1.0, 
	0.0, 1.0, 0.0, 
	q[i][0]
	33.63, 39.62, 38.18, 41.46, 38.02, 35.83, 35.99, 36.58, 42.92, 37.80, 40.43, 37.89, 
	
	The results are follows
	
	
	Residual sum of squares =  2.2227e+001
	Degrees of freedom = 8.0
	
	Variable       Parameter           Estimate Standard error
	
    	 1         3.6003e+001         9.6235e-001
    	 2         3.7300e+001         9.6235e-001
    	 3         4.1603e+001         9.6235e-001
    	 4         3.7877e+001         9.6235e-001
    	 
Return:
	The function returns NAG error code or 0 if no error.
	
	11: On entry, n must not be less than 1: n = _value_.  On entry, ip must not be less than 1: ip = _value_.
	17: On entry tdq = _value_ while ip + 1 = _value_. These parameters must satisfy tdq = ip + 1.  On entry, n = _value_ while ip = _value_. These parameters must satisfy n = ip.
	430: The degrees of freedom for error are less than or equal to 0. In this case the estimates, ˆ ß, are returned but not the standard errors or covariances.
	427: The singular value decomposition has failed to converge.  See nag_real_svd (f02wec). This is an unlikely error exit.
	5: On entry, tol must not be less than 0.0: tol = _value_.
	73: Memory allocation failed.
		
	successful call of the function int nag_regsn_mult_linear_upd_model function.
	

*/
  
  int nag_regsn_mult_linear_upd_model(
    int n, 		// the number of observations
    int ip, 	// the number of terms in the regression model
    double q[], // array q is taken from output by nag_regsn_mult_linear_addrem_obs, nag_regsn_mult_linear_add_var or nag_regsn_mult_linear_delete_var.
    int tdq, 	// the last dimension of the array q as declared in the function from which nag_regsn_mult_linear_upd_model is called.
    double *rss, // the residual sum of squares.
    double *df, //  the degree of freedom associated with the residual sum of squares.
    double b[], //  the parameter estimates.
    double se[], // the standard error of the parameter estimates.
    double cov[], // the estimate variance-covariance matrix of parameter extimates.
    Boolean *svd, 	// TRUE if a singular value decomposition has been performed, otherwise FALSE.
    int *rank, 	// the rank of the independent variables.
    double p[], // contains details of the singular value decomposition if used.
    double tol  // the value of tol is used to decide if the independent variables are of full rank. 
  ); // Regression coefficients, standard errors and variance-covariance matrix.
  
/** g02dec 
		adds a new independent variable to a general linear regression model.
	
Example:
	A data set consisting of 12 observations is read in. The four independent variables are stored in
	the array x while the dependent variable is read into the first column of q. If the character variable
	meanc indicates that a mean should be included in the model, a variable taking the value 1.0 for
	all observations is set up and fitted. Subsequently, one variable at a time is selected to enter the
	model as indicated by the input value of indx. After the variable has been added the parameter
	estimates are calculated by nag_regsn_mult_linear_upd_model (g02ddc) and the results printed. 
	

void test_nag_regsn_mult_linear_add_var()
{

	double rss, rsst, tol;
	int i, ip, rank, j, m, n;
	double df;
	Boolean svd;
	char meanc, weight;
	Nag_IncludeMean mean;
	double b[5];
	double cov[15];
	double p[35];
	double se[5], wt[12];
	matrix y;
	y.SetSize(12, 5);
	double xe[12];
	double *wtptr;
	int success;
	
	n = 12;
	m = 4;
	weight = 'u';
	meanc = 'm';
	mean = Nag_MeanInclude;
	wtptr = NULL;

	matrix x = 	{{1.0, 1.4, 0.0, 0.0, 0}, {1.5, 2.2, 0.0, 0.0, 0},
				 {2.0, 4.5, 0.0, 0.0, 0}, {2.5, 6.1, 0.0, 0.0, 0},
				 {3.0, 7.1, 0.0, 0.0, 0}, {3.5, 7.7, 0.0, 0.0, 0},
				 {4.0, 8.3, 1.0, 4.0, 0}, {4.5, 8.6, 1.0, 4.5, 0},
				 {5.0, 8.8, 1.0, 5.0, 0}, {5.5, 9.0, 1.0, 5.5, 0},
				 {6.0, 9.3, 1.0, 6.0, 0}, {6.5, 9.2, 1.0, 6.5, 0}};
	
	matrix q = 	{{4.32, 0, 0, 0, 0, 0}, {5.21, 0, 0, 0, 0, 0},
	             {6.49, 0, 0, 0, 0, 0}, {7.10, 0, 0, 0, 0, 0},
				 {7.94, 0, 0, 0, 0, 0}, {8.53, 0, 0, 0, 0, 0},
				 {8.84, 0, 0, 0, 0, 0}, {9.02, 0, 0, 0, 0, 0},
				 {9.27, 0, 0, 0, 0, 0}, {9.43, 0, 0, 0, 0, 0},
				 {9.68, 0, 0, 0, 0, 0}, {9.83, 0, 0, 0, 0, 0}};
	
	int index[5] = {1, 3, 4, 2, 0};
	  
	printf("The input data are as follows:\n");
	printf("\nm = 4, n = 12, wtptr = NULL,tol = 0.0000001e0");
	printf("ip = 0, mean = Nag_MeanInclude \n");
	printf("x\n");
	for(i = 0; i < 12; i++)
	{
	 	for( j = 0; j < 4; j++)
	 		printf("%2.1f, ",x[i][j]);
	 	printf("\n");
	}
	printf("\nq[i][0]\n");
	for(i = 0; i < 12; i++)
		printf("%3.2f, ",q[i][0]);
	 		
	
	tol = 0.000001e0;
	ip = 0;
	if(mean == Nag_MeanInclude)
	{
		for (i = 0; i<n; ++i)
			xe[i] = 1.0;
		nag_regsn_mult_linear_add_var(n, ip, q, 6, p, wtptr, xe, &rss,
		tol);
		ip = 1;
	}
	printf("\nxe\n");
	for(i = 0; i < n; i++)
		printf("%2.1f, ",xe[i]);
	printf("\nThe results are as follows\n");
	for(i = 0; i < 5; i++)
	{
		if(index[i] > 0)
		{
			for (j=0; j<n; j++)
				xe[j] = x[j][index[i]-1];
			success = nag_regsn_mult_linear_add_var(n, ip, q, 6, p, wtptr, xe, &rss,tol);
			if(success == 0)
			{
				ip += 1;
				printf("Variable %4ld added\n", index[i]);
				rsst = 0.0;
				nag_regsn_mult_linear_upd_model(n, ip, q, 6, &rsst, &df, b, se,
				cov, &svd, &rank, p, tol);
				if (svd)
					printf("Model not of full rank\n\n");
				printf("Residual sum of squares = %13.4e\n", rsst);
				printf("Degrees of freedom = %3.1f\n\n", df);
				printf("Variable       Parameter           Estimate Standard error\n\n");
				for (j=0; j<ip; j++)
					printf("%6ld%20.4e%20.4e\n", j+1, b[j], se[j]);
				printf("\n");
			}
			else if(success == NE_NVAR_NOT_IND)
				printf(" * Newvariable not added *\n");
			else
				printf("This function has some problem\n");			
		}
	}
}

	The output is as follows:
	
	The input data are as follows:
	
	m = 4, n = 12, wtptr = NULL,tol = 0.0000001e0ip = 0, mean = Nag_MeanInclude 
	x
	1.0, 1.4, 0.0, 0.0, 
	1.5, 2.2, 0.0, 0.0, 
	2.0, 4.5, 0.0, 0.0, 
	2.5, 6.1, 0.0, 0.0, 
	3.0, 7.1, 0.0, 0.0, 
	3.5, 7.7, 0.0, 0.0, 
	4.0, 8.3, 1.0, 4.0, 
	4.5, 8.6, 1.0, 4.5, 
	5.0, 8.8, 1.0, 5.0, 
	5.5, 9.0, 1.0, 5.5, 
	6.0, 9.3, 1.0, 6.0, 
	6.5, 9.2, 1.0, 6.5, 

	q[i][0]
	4.32, 5.21, 6.49, 7.10, 7.94, 8.53, 8.84, 9.02, 9.27, 9.43, 9.68, 9.83, 
	xe
	1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 
	The results are as follows
	Variable    1 added
	Residual sum of squares =   4.0164e+000
	Degrees of freedom = 10.0
	
	Variable       Parameter           Estimate Standard error
	
     	 1         4.4100e+000         4.3756e-001
    	 2         9.4979e-001         1.0599e-001

	Variable    3 added
	Residual sum of squares =   3.8872e+000
	Degrees of freedom = 9.0

	Variable       Parameter           Estimate Standard error

	     1         4.2236e+000         5.6734e-001
	     2         1.0554e+000         2.2217e-001
	     3        -4.1962e-001         7.6695e-001

	Variable    4 added
	Residual sum of squares =   1.8702e-001
	Degrees of freedom = 8.0

	Variable       Parameter           Estimate Standard error

         1         2.7605e+000         1.7592e-001
 	     2         1.7057e+000         7.3100e-002
	     3         4.4575e+000         4.2676e-001
	     4        -1.3006e+000         1.0338e-001

	Variable    2 added
	Residual sum of squares =   8.4066e-002
	Degrees of freedom = 7.0

	Variable       Parameter           Estimate Standard error

	     1         3.1440e+000         1.8181e-001
	     2         9.0748e-001         2.7761e-001
    	 3         2.0790e+000         8.6804e-001
    	 4        -6.1589e-001         2.4530e-001
    	 5         2.9224e-001         9.9810e-002
	

Return:
	The function returns NAG error code or 0 if no error.
	
	11: On entry, n must not be less than 1: n = _value_.  On entry, ip must not be less than 0: ip = _value_.
	20: On entry ip = _value_ while n = _value_. These parameters must satisfy ip < n.
	17: On entry tdq = _value_ while ip+2 = _value_. These parameters must satisfy tdq = ip+2.
	5: On entry, wt[_value_] must not be less than 0.0: wt[_value_] = _value_.
	6: On entry, tol must not be less than or equal to 0.0: tol = _value_.
	431: The new independent variable is a linear combination of existing variables. The (ip+1)th column of q is, therefore, null.
		
	successful call of the int nag_regsn_mult_linear_add_var function.
	

*/  

  int nag_regsn_mult_linear_add_var(
    int n, 		// the number of observations
    int ip, 	// the number of independent variables already in the model
    double q[], // the results of the QR decomposition for the model with p + 1 parameters.
    int tdq, 	// the last dimension of the array q. 
    double p[], // contains further details of the QR decomposition used. 
    double wt[], // an (optional) weighted is specified to be used in the weighted regression.
    double x[], // the new independent variable
    double *rss, // the residual sum of squares for the new fitted model
    double tol	// the value of tol is used to decide if the new independent variable is linearly related to independent variables already included in the model. 
  ); // Residual sum of squares.
  
/** g02dfc
		deletes an independent variable from a general linear regression model.
	
		
Example:
	A data set consisting of 12 observations on four independent variables and one dependent 
	variable is read in. The full model, including a mean term, is fitted using 
	nag_regsn_mult_linear (g02dac). The value of indx is read in and that variable dropped from 
	the regression. The parameter estimates are calculated by nag_regsn_mult_linear_upd_model (g02ddc) 
	and printed. This process is repeated until indx is 0.

void test_nag_regsn_mult_linear_delete_var()
{
	double rss, tol;
	int i, indx, ip, rank, j, m, n;
	double df;
	Boolean svd;
	char meanc, weight;
	Nag_IncludeMean mean;
	double b[5];
	double cov[15], h[12], p[35];
	matrix q;
	q.SetSize(12,6);
	double res[12], se[5], com_ar[45],wt[12];
	double *wtptr;
	int sx[5];
	n = 12;
	m = 4;
	weight = 'u';
	meanc = 'm';
	mean = Nag_MeanInclude;
	wtptr = NULL;

	matrix x = {{1.0, 1.4, 0.0, 0.0, 0}, {1.5, 2.2, 0.0, 0.0, 0},
				{2.0, 4.5, 0.0, 0.0, 0}, {2.5, 6.1, 0.0, 0.0, 0},
				{3.0, 7.1, 0.0, 0.0, 0}, {3.5, 7.7, 0.0, 0.0, 0},
				{4.0, 8.3, 1.0, 4.0, 0}, {4.5, 8.6, 1.0, 4.5, 0},
				{5.0, 8.8, 1.0, 5.0, 0}, {5.5, 9.0, 1.0, 5.5, 0},
				{6.0, 9.3, 1.0, 6.0, 0}, {6.5, 9.2, 1.0, 6.5, 0}};		
  	
	double y[12] = {4.32, 5.21, 6.49, 7.10, 7.94, 8.53, 
					8.84, 9.02, 9.27, 9.43, 9.68, 9.83};
	
	int index[3] = {2, 4, 0};
	
	printf("The input data are as follows:\n");
	printf("\nm = 4, n = 12, wtptr = NULL,tol = 0.0000001e0");
	printf("ip = 0, mean = Nag_MeanInclude \n");
	printf("x\n");
	for(i = 0; i < 12; i++)
	{
	 	for( j = 0; j < 4; j++)
	 		printf("%2.1f, ",x[i][j]);
	 	printf("\n");
	}
	printf("\ny\n");
	for(i = 0; i < 12; i++)
		printf("%3.2f, ",y[i]);
	printf("\nindex\n");
	for(i = 0; i < 3; i++)
		printf("%d, ",index[i]);
	for (i=0; i<m; ++i)
		sx[i] = 1;
	printf("\nsx\n");
	for(i = 0; i < m; i++)
		printf("%d, ",sx[i]);
	printf("\nThe results are as follows\n");
	 		
	ip = m;
	if (mean==Nag_MeanInclude)
		ip += 1;
	printf("\nip = %d\n",ip);
	tol = 0.00001e0;
	nag_regsn_mult_linear(mean, n, x, 5, m, sx, ip, y, wtptr, &rss,
							&df, b, se, cov, res, h, q, 6, &svd, &rank, p, tol, com_ar);
	printf("Results from full model\n");
	if(svd)
		printf("Model not of full rank\n\n");
	printf("Residual sum of squares = %13.4e\n", rss);
	printf("Degrees of freedom = %3.1f\n\n", df);

	for(i = 0; i < 3; i++)
	{
		if (index[i] != 0)
		{
			nag_regsn_mult_linear_delete_var(ip, q, 6, index[i], &rss);
			ip = ip - 1;
			if (ip == 0)
				printf("No terms left in model\n");
			else
			{
				printf("Variable %4ld dropped\n", index[i]);
				nag_regsn_mult_linear_upd_model(n, ip, q, 6, &rss, &df, b, se,cov, &svd, &rank, p, tol);
				printf("Residual sum of squares = %13.4e\n", rss);
				printf("Degrees of freedom = %3.1f\n\n", df);
				printf("Parameter       Estimate Standard error\n\n");
				for (j=0; j<ip; j++)
					printf("%12.4e%15.4e\n", b[j], se[j]);
			}
		}
	}
}

	The output is as follows:

	The input data are as follows:

	m = 4, n = 12, wtptr = NULL,tol = 0.0000001e0ip = 0, mean = Nag_MeanInclude 
	x
	1.0, 1.4, 0.0, 0.0, 
	1.5, 2.2, 0.0, 0.0, 
	2.0, 4.5, 0.0, 0.0, 
	2.5, 6.1, 0.0, 0.0, 
	3.0, 7.1, 0.0, 0.0, 
	3.5, 7.7, 0.0, 0.0, 
	4.0, 8.3, 1.0, 4.0, 
	4.5, 8.6, 1.0, 4.5, 
	5.0, 8.8, 1.0, 5.0, 
	5.5, 9.0, 1.0, 5.5, 
	6.0, 9.3, 1.0, 6.0, 
	6.5, 9.2, 1.0, 6.5, 

	y
	4.32, 5.21, 6.49, 7.10, 7.94, 8.53, 8.84, 9.02, 9.27, 9.43, 9.68, 9.83, 
	index
	2, 4, 0, 
	sx
	1, 1, 1, 1, 
	The results are as follows

	ip = 5
	Results from full model
	Residual sum of squares =   8.4066e-002
	Degrees of freedom = 7.0
	
	Variable    2 dropped
	Residual sum of squares =   2.1239e-001
	Degrees of freedom = 8.0
	
	Parameter       Estimate Standard error
	
     3.6372e+000    1.5083e-001
 	 6.1264e-001    2.8007e-002
	-6.0154e-001    4.2335e-001
 	 1.6709e-001    7.8656e-002
	Variable    4 dropped
	Residual sum of squares =   3.3220e-001
	Degrees of freedom = 9.0
	
	Parameter       Estimate Standard error
	
 	 3.5974e+000    1.7647e-001
 	 6.2088e-001    3.2706e-002
 	 2.4247e-001    1.7235e-001

Return:
	The function returns NAG error code or 0 if no error.
	
	11: On entry, ip must not be less than 1: ip = _value_.  On entry, indx must not be less than 1: indx = _value_.
	17: On entry tdq = _value_ while ip+1 = _value_. These parameters must satisfy tdq = ip+1.
	19: On entry indx = _value_ while ip = _value_. These parameters must satisfy indx = ip.
	5: On entry, rss must not be less than 0.0: rss = _value_.
	441: On entry, a diagonal element, _value_, of R is zero.
	73 Memory allocation failed.
		
	successful call of the nag_regsn_mult_linear_delete_var function.
	
*/

  int nag_regsn_mult_linear_delete_var(
    int ip, 	// the number of independent variables already in the model
    double q[],	// the results of the QR decomposition as returned by the other nag regsn mult linear functions. 
    int tdq,	// the last dimension of the array q.
    int indx,	// indicates which independent variable is to be deleted from the model.
    double *rss // the residual sum of squares for the full regression.
  ); // Residual Sum of Square.
  
/** g02dgc
		calculates the estimates of the parameters of a general linear regression model for a new dependent 
		variable after a call to nag_regsn_mult_linea.
	
		
Example:
	A data set consisting of 12 observations with four independent variables and two dependent variables
	is read in. A model with all four independent variables is fitted to the first dependent variable by
	nag_regsn_mult_linear (g02dac) and the results printed. The model is then fitted to the second
	dependent variable by nag_regsn_mult_linear_newyvar and those results printed.

void test_nag_regsn_mult_linear_newyvar()
{
	double rss, tol;
	int i, ip, rank, j, m, n;
	double df;
	Boolean svd;
	Nag_IncludeMean mean;
	char weight, meanc;
	double b[5], cov[15], h[12],p[35];
	matrix q;
	q.SetSize(12,6);
	double res[12], se[5], com_ar[45], wt[12];
	double *wtptr;
	n = 12;
	m = 4;
	weight = 'u';
	meanc = 'm';
	if (meanc=='m')
		mean = Nag_MeanInclude;
	else
		mean = Nag_MeanZero;
	
	wtptr = NULL;
	
	matrix xm = {{1.0, 0.0, 0.0, 0.0, 0}, {0.0, 0.0, 0.0, 1.0, 0},
				 {0.0, 1.0, 0.0, 0.0, 0}, {0.0, 0.0, 1.0, 0.0, 0},
				 {0.0, 0.0, 0.0, 1.0, 0}, {0.0, 1.0, 0.0, 0.0, 0},
				 {0.0, 0.0, 0.0, 1.0, 0}, {1.0, 0.0, 0.0, 0.0, 0},
				 {0.0, 0.0, 1.0, 0.0, 0}, {1.0, 0.0, 0.0, 0.0, 0},
				 {0.0, 0.0, 1.0, 0.0, 0}, {0.0, 1.0, 0.0, 0.0, 0}};		
	
	double y[12] = {33.63, 39.62, 38.18, 41.46, 38.02, 35.83,
					35.99, 36.58, 42.92, 37.80, 40.43, 37.89};
	
	int sx[5] =	{1, 1, 1, 1};
	
	double newy[12] = {63.0, 69.0, 68.0, 71.0, 68.0, 65.0,
					   65.0, 66.0, 72.0, 67.0, 70.0, 67.0};
	
	ip = 5;	
	for (j=0; j<m; j++)
	tol = 0.00001e0;
	printf("The input data are as follows:\n");
	printf("\nm = 4, n = 12, wtptr = NULL,tol = 0.0000001e0");
	printf("ip = 0, mean = Nag_MeanInclude, ip =5\n");
	printf("x\n");
	for(i = 0; i < 12; i++)
	{
	 	for( j = 0; j < 4; j++)
	 		printf("%2.1f, ",xm[i][j]);
	 	printf("\n");
	}
	printf("\ny\n");
	for(i = 0; i < 12; i++)
		printf("%3.2f, ",y[i]);
	printf("\nnewy\n");
	for(i = 0; i < 12; i++)
		printf("%3.2f, ",newy[i]);

	printf("\nsx\n");
	for(i = 0; i < m; i++)
		printf("%d, ",sx[i]);
	printf("\nThe results are as follows\n");
	 		
	
	nag_regsn_mult_linear(mean, n, xm, 5, m, sx, ip,	y, wtptr, &rss, &df, 
						b, se, cov, res, h, q, 6, &svd, &rank, p, tol, com_ar);
	printf("Results from g02dac\n\n");
	if(svd)
		printf("Model not of full rank\n\n");
	printf("Residual sum of squares = %12.4e\n", rss);
	printf("Degrees of freedom = %3.1f\n\n", df);

	printf("Variable       Parameter          Estimate Standard error\n\n");
	for (j=0; j<ip; j++)
		printf("%6ld%20.4e%20.4e\n", j+1, b[j], se[j]);
	printf("\n");

	nag_regsn_mult_linear_newyvar(n, wtptr, &rss, ip, rank, cov, q, 6, svd, p,
			newy, b, se, res, com_ar);
	printf("\n");
	printf("Results for second y-variable using g02dgc\n\n");
	printf("Residual sum of squares = %12.4e\n", rss);
	printf("Degrees of freedom = %3.1f\n\n", df);
	
	printf("Variable       Parameter          Estimate Standard error\n\n");
	for (j=0; j<ip; j++)
		printf("%6ld%20.4e%20.4e\n", j+1, b[j], se[j]);
	printf("\n");
}
	
	The output is as follows:

	The input data are as follows:

	m = 4, n = 12, wtptr = NULL,tol = 0.0000001e0ip = 0, mean = Nag_MeanInclude, ip =5
	x
	1.0, 0.0, 0.0, 0.0, 
	0.0, 0.0, 0.0, 1.0, 
	0.0, 1.0, 0.0, 0.0, 
	0.0, 0.0, 1.0, 0.0, 
	0.0, 0.0, 0.0, 1.0, 
	0.0, 1.0, 0.0, 0.0, 
	0.0, 0.0, 0.0, 1.0, 
	1.0, 0.0, 0.0, 0.0, 
	0.0, 0.0, 1.0, 0.0, 
	1.0, 0.0, 0.0, 0.0, 
	0.0, 0.0, 1.0, 0.0, 
	0.0, 1.0, 0.0, 0.0, 
	
	y
	33.63, 39.62, 38.18, 41.46, 38.02, 35.83, 35.99, 36.58, 42.92, 37.80, 40.43, 37.89, 
	newy
	63.00, 69.00, 68.00, 71.00, 68.00, 65.00, 65.00, 66.00, 72.00, 67.00, 70.00, 67.00, 
	sx
	1, 1, 1, 1, 
	The results are as follows
	Results from g02dac
	
	Model not of full rank
	
	Residual sum of squares =  2.2227e+001
	Degrees of freedom = 8.0
	
	Variable       Parameter          Estimate Standard error
	
    	 1         3.0557e+001         3.8494e-001
    	 2         5.4467e+000         8.3896e-001
    	 3         6.7433e+000         8.3896e-001
    	 4         1.1047e+001         8.3896e-001
    	 5         7.3200e+000         8.3896e-001


	Results for second y-variable using g02dgc
	
	Residual sum of squares =  2.4000e+001
	Degrees of freedom = 8.0
	
	Variable       Parameter          Estimate Standard error
	
    	 1         5.4067e+001         4.0000e-001
     	 2         1.1267e+001         8.7178e-001
     	 3         1.2600e+001         8.7178e-001
    	 4         1.6933e+001         8.7178e-001
    	 5         1.3267e+001         8.7178e-001

Return:
	The function returns NAG error code or 0 if no error.
	
	11: On entry, ip must not be less than 1: ip = _value_.
	12: On entry, rank must not be less than or equal to 0: rank = _value_.
	17: On entry, tdq = _value_ while ip + 1 = _value_. These parameters must satisfy tdq = ip + 1.  On entry, n = _value_ while ip = _value_. These parameters must satisfy n = ip.
	6: On entry, rss must not be less than or equal to 0.0: rss = _value_.
	5: On entry, wt[_value_] must not be less than 0.0: wt[_value_] = _value_.
	434: On entry, the Boolean variable, svd, is FALSE and rank must be equal to ip: rank = _value_, ip = _value_.
	435: On entry, the Boolean variable, svd, is TRUE and rank must not be greater than ip: rank = _value_, ip = _value_.
	
	successful call of the nag_regsn_mult_linear_newyvar function.
	
*/
  
  int nag_regsn_mult_linear_newyvar(
    int n, 			// the number of observations
    double wt[],	// an optional weight is specified to be used in the weighted regression. 
    double *rss,	// the residual sum of squares for the original dependent variable.
    int ip,			// the number p of independent variables in the model.
    int rank,		// the rank of the independent variables as given by (g02dac)
    double cov[],	// the estimate variance-covariance matrix of the parameter estimates as given by (g02dac)
    double q[],		// the results of the QR decomposition as returned by (g02dac)
    int tdq,		// the second dimension of the array q. 
    Boolean svd,		// indicates if a singular value decomposition was used by (g02dac).
    double p[],		// details of the QR decomposition and SVD, if used, as returned in array p by (g02dac).
    double y[],		// the new dependent variable.
    double b[],		// contain the least-squares estimates of the parameters of the regression model.
    double se[],	// contain the standard errors of the ip parameter estimates given in b.
    double res[],	// the residuals for the new regression model.
    double com_ar[] // if svd = TRUE, com_ar must be unaltered from the previous call to (g02dac).
  ); // Estimates of parameters, Standard errors of parameter estimates.
  
/** go2dkc
		calculates the estimates of the parameters of a general linear regression model for given 
		constraints from the singular value decomposition results.
	
Example:
	Data from an experiment with four treatments and three observations per treatment are read in. A
	model, including the mean term, is fitted by nag_regsn_mult_linear (g02dac) and the results printed.
	The constraint that the sum of treatment effects is zero is then read in and the parameter estimates
	with this constraint imposed are computed by nag_regsn_mult_linear_tran_model and printed.
	
void test_nag_regsn_mult_linear_tran_model()
{

	double rss, tol;
	int i, iconst, ip, rank, j, m, n;
	double df;
	Boolean svd;
	Nag_IncludeMean mean;
	char weight, meanc;
	double b[5];
	matrix c;
	c.SetSize(5,5);
	double cov[15], h[12],p[35];
	matrix q;
	q.SetSize(12,6);
	double res[12], se[5],	com_ar[120], wt[12];
	double *wtptr;
	n = 12;
	m = 4;
	weight = 'u';
	meanc = 'm';
	if (meanc=='m')
		mean = Nag_MeanInclude;
	else
		mean = Nag_MeanZero;
	wtptr = NULL;
	
	
	matrix x =  {{1.0, 0.0, 0.0, 0.0, 0}, {0.0, 0.0, 0.0, 1.0, 0},
				 {0.0, 1.0, 0.0, 0.0, 0}, {0.0, 0.0, 1.0, 0.0, 0}, 
				 {0.0, 0.0, 0.0, 1.0, 0}, {0.0, 1.0, 0.0, 0.0, 0}, 
				 {0.0, 0.0, 0.0, 1.0, 0}, {1.0, 0.0, 0.0, 0.0, 0}, 
				 {0.0, 0.0, 1.0, 0.0, 0}, {1.0, 0.0, 0.0, 0.0, 0}, 
				 {0.0, 0.0, 1.0, 0.0, 0}, {0.0, 1.0, 0.0, 0.0, 0}};

	double y[12] = {33.63, 39.62, 38.18, 41.46, 38.02, 35.83,
					35.99, 36.58, 42.92, 37.80, 40.43, 37.89}; 
	int sx[5] = {1, 1, 1, 1};
	
	double temp[5] = {0.0, 1.0, 1.0, 1.0, 1.0};
	
	ip = 5;	
	tol = 0.00001e0;
	printf("The input data are as follows:\n");
	printf("\nm = 4, n = 12, wtptr = NULL,tol = 0.0000001e0");
	printf("ip = 5, mean = Nag_MeanInclude, ip =5\n");
	printf("x\n");
	for(i = 0; i < 12; i++)
	{
	 	for( j = 0; j < 4; j++)
	 		printf("%2.1f, ",x[i][j]);
	 	printf("\n");
	}
	printf("\ny\n");
	for(i = 0; i < 12; i++)
		printf("%3.2f, ",y[i]);
	printf("\nsx\n");
	for(i = 0; i < m; i++)
		printf("%d, ",sx[i]);
	printf("\nThe results are as follows\n");
	nag_regsn_mult_linear(mean, n, x, 5, m, sx, ip, y, wtptr, &rss, &df, b, se, cov, 
			res, h, q, 6, &svd, &rank, p, tol, com_ar);
	printf("Estimates from g02dac\n\n");
	printf("Residual sum of squares = %13.4e\n", rss);
	printf("Degrees of freedom = %3.1f\n\n", df);

	printf("Variable       Parameter           Estimate Standard error\n\n");
	for (j=0; j<ip; j++)
		printf("%6ld%20.4e%20.4e\n", j+1, b[j], se[j]);
	printf("\n");

	iconst = ip - rank;
	for (i=0; i<ip; ++i)
		for (j=0; j<iconst; ++j)
		c[i][j] = temp[i * iconst  + j];
	nag_regsn_mult_linear_tran_model(ip, iconst, p, c, 5, b, rss, df, se, cov);
	
	printf("\n");
	printf("Estimates from g02dkc using constraints\n\n");
	printf("Variable       Parameter           Estimate Standard error\n\n");
	for (j=0; j<ip; j++)
		printf("%6ld%20.4e%20.4e\n", j+1, b[j], se[j]);
	printf("\n");
}			

	The output is as follows:
		
	The input data are as follows:

	m = 4, n = 12, wtptr = NULL,tol = 0.0000001e0ip = 5, mean = Nag_MeanInclude, ip =5
	x
	1.0, 0.0, 0.0, 0.0, 
	0.0, 0.0, 0.0, 1.0, 
	0.0, 1.0, 0.0, 0.0, 
	0.0, 0.0, 1.0, 0.0, 
	0.0, 0.0, 0.0, 1.0, 
	0.0, 1.0, 0.0, 0.0, 
	0.0, 0.0, 0.0, 1.0, 
	1.0, 0.0, 0.0, 0.0, 
	0.0, 0.0, 1.0, 0.0, 
	1.0, 0.0, 0.0, 0.0, 
	0.0, 0.0, 1.0, 0.0, 
	0.0, 1.0, 0.0, 0.0, 

	y
	33.63, 39.62, 38.18, 41.46, 38.02, 35.83, 35.99, 36.58, 42.92, 37.80, 40.43, 37.89, 
	sx
	1, 1, 1, 1, 
	The results are as follows
	Estimates from g02dac
	
	Residual sum of squares =   2.2227e+001
	Degrees of freedom = 8.0
	
	Variable       Parameter           Estimate Standard error
	
    	 1         3.0557e+001         3.8494e-001
    	 2         5.4467e+000         8.3896e-001
    	 3         6.7433e+000         8.3896e-001
    	 4         1.1047e+001         8.3896e-001
    	 5         7.3200e+000         8.3896e-001
	

	Estimates from g02dkc using constraints
	
	Variable       Parameter           Estimate Standard error
	
    	 1         3.8196e+001         4.8117e-001
    	 2        -2.1925e+000         8.3342e-001
    	 3        -8.9583e-001         8.3342e-001
    	 4         3.4075e+000         8.3342e-001
    	 5        -3.1917e-001         8.3342e-001


Return:
	The function returns NAG error code or 0 if no error.

	11: On entry, ip must not be less than 1: ip = _value_.
	12: On entry, iconst must not be less than or equal to 0: iconst = _value_.
	20: On entry, iconst = _value_ while ip = _value_. These parameters must satisfy iconst < ip.
	17: On entry, tdc = _value_ while iconst = _value_. These parameters must satisfy tdc = iconst.
	6: On entry, rss must not be less than or equal to 0.0: rss = _value_.  On entry, df must not be less than or equal to 0.0: df = _value_.
	438: matrix c does not give a model of full rank.
	73: Memory allocation failed.
		
	successful call of the nag_regsn_mult_linear_tran_model function.
	
*/
  
  int nag_regsn_mult_linear_tran_model(
    int ip,		// the number of terms in the linear model. 
    int iconst,	// the number of constraints to be imposed on the parameters.
    double p[],	// p as returned by (g02dac) and (g02ddc).
    double c[],	// the iconst constraints stored by column.
    int tdc,	// the last dimension of the array c.
    double b[],	// Input: the parameter estimates computed by using the singular value decomposition.  Output: the parameter estimates of the parameters with the constraints imposed.
    double rss, // the residual sum of squares as returned by (g02dac) or (g02ddc).
    double df,	// the degrees of freedom associated with the residual sum of squares as returned by (g02dac) or (g02ddc). 
    double se[], // the standard error of the parameter estimates in b.
    double cov[] // the upper triangular part of the variance-covariance matrix of the ip parameter estimates given in b.
  ); // Parameter estimates, Variance-covariance matrix, Residual sum of squares...
  
/** g02dnc
		gives the estimate of an estimable function along with its standard error.
	
Example:
	Data from an experiment with four treatments and three observations per treatment 
	are read in.  This example uses three different linear regression functions.  
	Depend on what you want, you can assign each coefficient for each function.

void test_nag_regsn_mult_linear_est_func()
{
	
	double rss, sestat, stat, t, tol;
	int i, ip, rank, j, m, n, nestern;
	double df;
	Boolean est, svd;
	Nag_IncludeMean mean;
	char meanc, weight;
	double b[5], cov[15], f[5], h[12], p[35];
	matrix q;
	q.SetSize(12,6);
	double  res[12], se[5], com_ar[45], wt[12];
	double *wtptr;
	int success;

	n = 12;
	m = 4;
	weight = 'u';
	meanc = 'm';
	
	if (meanc=='m')
		mean = Nag_MeanInclude;
	else
		mean = Nag_MeanZero;
	if (weight=='w')
		wtptr = wt;
	else

	wtptr = NULL;
	matrix x = {{1.0, 0.0, 0.0, 0.0, 0}, {0.0, 0.0, 0.0, 1.0, 0},
				{0.0, 1.0, 0.0, 0.0, 0}, {0.0, 0.0, 1.0, 0.0, 0},
				{0.0, 0.0, 0.0, 1.0, 0}, {0.0, 1.0, 0.0, 0.0, 0},
				{0.0, 0.0, 0.0, 1.0, 0}, {1.0, 0.0, 0.0, 0.0, 0},
				{0.0, 0.0, 1.0, 0.0, 0}, {1.0, 0.0, 0.0, 0.0, 0},
				{0.0, 0.0, 1.0, 0.0, 0}, {0.0, 1.0, 0.0, 0.0, 0}};		

	double y[12] = {33.63, 39.62, 38.18, 41.46, 38.02, 35.83,
					35.99, 36.58, 42.92, 37.80, 40.43, 37.89};
	
	int sx[5] =	{1, 1, 1, 1};
	
	ip = 5;
	tol = 0.00001e0;
	nag_regsn_mult_linear(mean, n, x, 5, m, sx, ip, y, wtptr, &rss, &df, b, se, cov, 
						res, h, q, 6, &svd, &rank, p, tol, com_ar);
	printf("\n");
	printf("Estimates from g02dac\n\n");
	printf("Residual sum of squares = %12.4e\n", rss);
	printf("Degrees of freedom = %3.1f\n\n", df);
	
	printf("Variable       Parameter          Estimate Standard error\n\n");
	for (j=0; j<ip; j++)
		printf("%6ld%20.4e%20.4e\n", j+1, b[j], se[j]);
	printf("\n");
	nestern = 3;	
	matrix fl = {{1.0, 1.0, 0.0, 0.0, 0.0}, {0.0, 1.0,-1.0, 0.0, 0.0},
				 {0.0, 1.0, 0.0, 0.0, 0.0}};
				
	for (i = 0; i < nestern; ++i)
	{
		for(j = 0; j < 5; j++)
			f[j] = fl[i][j];
		success = nag_regsn_mult_linear_est_func(ip, rank, b, cov, p, f, &est, &stat, &sestat, &t, tol);
		if (success == 0 || success == NE_RANK_EQ_IP)
		{
			printf("\n");
			printf("Function %ld\n\n", i+1);
			for (j=0; j<ip; ++j)
				printf("%8.2f%c", f[j], (j%5==4 || j==ip-1) ? '\n' : ' ');
			printf("\n");
		if(est)
			printf(" stat = %10.4f se = %10.4f t = %10.4f\n", stat, sestat, t);
		else
			printf("Function not estimable\n");
		}
	}
}

	The output is as follows:
	
	Estimates from g02dac

	Residual sum of squares =  2.2227e+001
	Degrees of freedom = 8.0
	
	Variable       Parameter          Estimate Standard error

    	 1         3.0557e+001         3.8494e-001
    	 2         5.4467e+000         8.3896e-001
    	 3         6.7433e+000         8.3896e-001
    	 4         1.1047e+001         8.3896e-001
    	 5         7.3200e+000         8.3896e-001


	Function 1

    	1.00    1.00    0.00    0.00    0.00
 	 stat =    36.0033 se =     0.9623 t =    37.4119

	Function 2

     	0.00    1.00   -1.00    0.00    0.00
  	 stat =    -1.2967 se =     1.3610 t =    -0.9528

	Function 3

    	0.00    1.00    0.00    0.00    0.00
	Function not estimable


Return:
	The function returns NAG error code or 0 if no error.
	
	11: On entry, ip must not be less than 1: ip = _value_.  On entry, rank must not be less than 1: rank = _value_.
	19: On entry ip = _value_ while rank = _value_. These parameters must satisfy rank = ip.
	440: On entry, rank = ip. In this case, the boolean variable est is returned as TRUE and all statistics are calculated.
	439: se(F) = 0.0 probably due to rounding error or due to incorrectly specified inputs cov and f.
	73: Memory allocation failed.
	
	successful call of the nag_regsn_mult_linear_est_func function.
	
*/

  int nag_regsn_mult_linear_est_func(
    int ip,		// the number of terms in the linear model.
    int rank,	// the rank of the independent variables.
    double b[], // the ip values of the estimates of the parameters of the model.
    double cov[], // the upper triangular part of the variance-covariance matrix of the ip parameter estimates given in b. 
    double p[], // p as returned by (g02dac) or (g02ddc).
    double f[], // the linear function to be estimated.
    Boolean *est,	// indicates if the function was estimable.
    double *stat, // if est = TRUE, stat contains the estimate of the function.
    double *sestat, // if est = TRUE, sestat contains the standard error of the estimate of the function.
    double *t,	// if est = TRUE, t contains the t-statistic for the test of the function being equal to zero.
    double tol // the tolerance value used in the check for estimability.
  ); // Estimate of an estimable function along with its standard error.

/** g02fac
		calculates two types of standardised residuals and two measures of influence for a linear regression.
	
Example:
	A set of 24 residuals and h(i) values from an 11 parameter model fitted to the cloud 
	seeding data considered in Cook and Weisberg (1982) are input and the standardised 
	residuals etc calculated and printed for the first 10 observations.
	
void test_nag_regsn_std_resid_influence()
{

	double rms;
	int i, ip, j, n, nres;
	matrix sres;
	sres.SetSize(24,4);
	n = 24;
	ip = 11;
	nres = 10;
	rms = 0.5798;
	
	double res[24] = { 	0.2660,-0.1387,-0.2971, 0.5926,-0.4013, 0.1396,
					   -1.3173, 1.1226, 0.0321,-0.7111, 0.3439,-0.4379,
					   	0.0633,-0.0936, 0.9968, 0.0209,-0.4056, 0.1396,
					  	0.0327, 0.2970,-0.2277, 0.5180, 0.5301,-1.0650};

	double h[24] = {	0.5519, 0.9746, 0.6256, 0.3144, 0.4106, 0.6268,
						0.5479, 0.2325, 0.4115, 0.3577, 0.3342, 0.1673,
						0.3874, 0.1705, 0.3466, 0.3743, 0.7527, 0.9069,
						0.2610, 0.6256, 0.2485, 0.3072, 0.5848, 0.4794};
						
	printf("The input data are as follows:\n");
	printf("n = 24, ip = 11, nres = 10,	rms = 0.5798,\n");
	printf("res\n");
	for(i = 0; i < n; i ++)
	{
		printf("%5.4f, ",res[i]);
		if(( i + 1) % 6 == 0)
			printf("\n");
	}	
	
	printf("\nh\n");
	for(i = 0; i < n; i ++)
	{
		printf("%5.4f, ",h[i]);
		if(( i + 1) % 6 == 0)
			printf("\n");
	}	
	printf("\nThe results are as follows:\n");
		
	nag_regsn_std_resid_influence(n, ip, nres, res, h, rms, sres);
	printf("\n");
	printf("Internally Externally\n");
	printf("Obs.       studentized   studentized   Cook's D     Atkinson's T\n");
	printf(" \t\t\t\t\t\t   residuals    residuals\n\n");
	for (i=0; i<nres; i++)
	{
		printf("%2ld", i+1);
		for (j=0; j<4; j++)
		printf("%14.3f", sres[i][j]);
		printf("\n");
	}
}

	The output is as follows:
		
	The input data are as follows:
	n = 24, ip = 11, nres = 10,	rms = 0.5798,
	res
	0.2660, -0.1387, -0.2971, 0.5926, -0.4013, 0.1396, 
	-1.3173, 1.1226, 0.0321, -0.7111, 0.3439, -0.4379, 
	0.0633, -0.0936, 0.9968, 0.0209, -0.4056, 0.1396, 
	0.0327, 0.2970, -0.2277, 0.5180, 0.5301, -1.0650, 

	h
	0.5519, 0.9746, 0.6256, 0.3144, 0.4106, 0.6268, 
	0.5479, 0.2325, 0.4115, 0.3577, 0.3342, 0.1673, 
	0.3874, 0.1705, 0.3466, 0.3743, 0.7527, 0.9069, 
	0.2610, 0.6256, 0.2485, 0.3072, 0.5848, 0.4794, 

	The results are as follows:

	Internally Externally
	Obs.       studentized   studentized   Cook's D     Atkinson's T
 						   				   residuals    residuals

  	 1         0.522         0.507         0.030         0.611
 	 2        -1.143        -1.158         4.557        -7.797
	 3        -0.638        -0.622         0.062        -0.875
	 4         0.940         0.935         0.037         0.689
	 5        -0.686        -0.672         0.030        -0.610
	 6         0.300         0.289         0.014         0.408
	 7        -2.573        -3.529         0.729        -4.223
	 8         1.683         1.828         0.078         1.094
	 9         0.055         0.053         0.000         0.048
	10        -1.165        -1.183         0.069        -0.960


Return:
	The function returns NAG error code or 0 if no error.
	
	11: On entry, ip must not be less than 1: ip = _value_.  On entry, nres must not be less than 1: nres = _value_.
	18: On entry, n = _value_ while ip + 1 = _value_. These parameters must satisfy n > ip + 1.
	17: On entry, nres = _value_ while n = _value_. These parameters must satisfy nres = n.
	6: On entry, rms must not be less than or equal to 0.0: rms = _value_.  On entry, h[_value_] must not be less than or equal to 0.0: h[_value_] = _value_.
	8: On entry, h[_value_] must not be greater than or equal to 1.0: h[_value_] = _value_.
	433: On entry, the value of a residual is too large for the given value of rms: res[_value_] = _value_, rms = _value_.
	
	successful call of the nag_regsn_std_resid_influence function.
	
*/

  int nag_regsn_std_resid_influence(
    int n,		// number of observations included in the regression.
    int ip,		// the number of linear parameters estimated in the regression model.
    int nres,	// the number of residuals.
    double res[], // the residuals.
    double h[],	// the diagonal elements of H, h(i), corresponding to the residuals in res.
    double rms,	// the estimate of s2 based on all n observations, s2, i.e., the residual mean square.
    double sres[] // the standardised residuals and influence statistics.
  ); // Internal studentized statistic, External studentized statistic, Cook's D residual, Atkinson's T residual.  
  
/** g02gac
		fits a generalized linear model with normal errors.
	
Example:
	The model: [y = ß1 + ß2x + epsilon] for a sample of 5 observations.

void test_nag_glm_normal()
{
	char linkc, meanc, weightc;
	Nag_Link link;
	Nag_IncludeMean mean;
	int i, j, m, n, ip;
	int max_iter, print_iter, rank;
	double ex_power, scale, tol, eps, rss, df;
	int sx[2];
	double b[2];
	matrix v;
	v.SetSize(5,8);
	double wt[5];
	matrix x;
	x.SetSize(5,2);
	
	double y[5], se[2], cov[3];
	double *wtptr;
	double *offsetptr;
	offsetptr = NULL;
	int success;
	linkc = 'r';
	meanc = 'm';
	weightc = 'n';
	n = 5;
	m = 1;
	print_iter = 0;
	scale = 0.0;
	
	set_enum(linkc, &link, meanc, &mean);

	wtptr = NULL;
	x[0][0] = 1.0;
	y[0] = 25.0;
	x[1][0] = 2.0;
	y[1] = 10.0;
	x[2][0] = 3.0;
	y[2] = 6.0;
	x[3][0] = 4.0;
	y[3] = 4.0;
	x[4][0] = 5.0;
	y[4] = 3.0;	
	
	sx[0] = 1;
	ip = 0;
	for (j=0; j<m; j++)
		if (sx[j]>0) ip += 1;
	if (mean == Nag_MeanInclude)
		ip += 1;
	if (link == Nag_Expo)
		printf("there is some problemtttttttt");
	else
		ex_power = 0.0;
	printf("ipip = %d",ip);
	max_iter = 10;
	tol = 5e-5;
	eps = 1e-6;
	String str = "c:\\c\\test.txt";
	success = nag_glm_normal(link, mean, n, x, 2, m, sx, ip, y, wtptr, 
								offsetptr, &scale, ex_power, &rss, &df, b, &rank,se, cov,
									v, 8, tol, max_iter,	print_iter, str, eps);
	
	if (success == NE_NOERROR || success == NE_LSQ_ITER_NOT_CONV || success == NE_RANK_CHANGED || success == NE_ZERO_DOF_ERROR)
	{
		printf("\nResidual sum of squares = %12.4e\n", rss);
		printf("Degrees of freedom = %3.1f\n\n", df);
		printf("       Estimate    Standard error\n\n");
		for (i=0; i<ip; i++)
			printf("%14.4f%14.4f\n", b[i], se[i]);
		printf("\n");
		printf("   y     fitted value     Residual    Leverage\n\n");
		for (i = 0; i < n; ++i)
		{
	
			printf("%7.1f%10.2f%12.4f%10.3f\n", y[i], v[i][1], v[i][4], v[i][5]);
		}		
		
	}
	else
	{
		printf("there are some problem");
	}
}

int set_enum(char linkc, Nag_Link *link, char meanc, Nag_IncludeMean *mean)
{
	if(linkc == 'E' ||linkc =='e')
		*link = Nag_Expo;
	else if(linkc == 'I' ||linkc =='i')
		*link = Nag_Iden;
	else if(linkc == 'L' ||linkc =='l')
		*link = Nag_Log;
	else if(linkc == 'S' ||linkc =='s')
		*link = Nag_Sqrt;
	else if(linkc == 'R' ||linkc =='r')
		*link = Nag_Reci;
	else
		return 1;
	if(meanc == 'M'|| meanc == 'm')
		*mean = Nag_MeanInclude;
	else if(meanc == 'Z' || meanc == 'z')
		*mean = Nag_MeanZero;
	else
		return 1;
		return 0;
}

	The output is as follows:
		
	ip = 2
	Residual sum of squares =  3.8717e-001
	Degrees of freedom = 3.0

	Estimate       Standard error

    -0.0239        0.0028
     0.0638        0.0026

	y    fitted value   Residual   Leverage

	25.0     25.04     -0.0387     0.995
	10.0      9.64      0.3613     0.458
 	 6.0      5.97      0.0320     0.268
 	 4.0      4.32     -0.3221     0.167
 	 3.0      3.39     -0.3878     0.112

Return:
	The function returns NAG error code or 0 if no error.
	
	70: On entry, parameter link had an illegal value.  On entry, parameter mean had an illegal value.
	11: On entry, n must not be less than 2: n = _value	.  On entry, m must not be less than 1: m = _value	.  On entry, ip must not be less than 1: ip = _value	.  On entry, max_iter must not be less than 0: max_iter = _value	.  On entry, sx[_value	] must not be less than 0: sx[_value	] = _value	.
	5: On entry, scale must not be less than 0.0: scale = _value	.  On entry, tol must not be less than 0.0: tol = _value	.  On entry, eps must not be less than 0.0: eps = _value	.  On entry, wt[_value	] must not be less than 0.0: wt[_value	] = _value	.
	17: On entry, tdx = _value	while m = _value	. These parameters must satisfy tdx = m.  On entry, tdv = _value	while ip = _value	. These parameters must satisfy tdv = ip+6.
	37: On entry, ex power = 0.0 , link = Nag Expo.  These parameters must satisfy link == Nag Expo && ex power _= 0.0.
	457: Parameter ip is incompatible with parameters mean and sx.
	458: Parameter ip is greater than the effective number of observations.
	459: A fitted value is at a boundary. This will only occur with link = Nag Expo, Nag Log or Nag Reci. This may occur if there are small values of y and the model is not suitable for the data. The model should be reformulated with, perhaps, some observations dropped.
	427: The singular value decomposition has failed to converge.
	463: The iterative weighted least-squares has failed to converge in max_iter = _value	 iterations.  The value of max_iter could be increased but it may be advantageous to examine the convergence using the print iter option. This may indicate that the convergence is slow because the solution is at a boundary in which case it may be better to reformulate the model.  
	464: The rank of the model has changed during the weighted least-squares iterations. The estimate for ß returned may be reasonable, but the user should check how the deviance has changed during iterations.
	465: The degrees of freedom for error are 0. A saturated model has been fitted.
	78: Cannot open file _string for appending.
	79: Cannot close file _string	.
	73: Memory allocation failed.
		
	successful call of the nag_glm_normal function.
	
*/

int nag_glm_normal(
    Nag_Link link,
    Nag_IncludeMean mean,
    int n,  		// the number of observations.
    double x[],		// contain the ith observation for the jth independent variable.
    int tdx,		// the second dimension of the array x.
    int m,			// the total number of independent variables.
    int sx[],		// indicates which independent variables are to be included in the model. 
    int ip,			// the number p of independent variables in the model, including the mean or intercept if present.
    double y[],		// observations on the dependent variable.
    double wt[],	// an (optional) weight is specified to be used with the model.
    double offset[], // an optional offset is specified. 
    double *scale,	// indicates the scale parameter for the model.
    double ex_power, // if link = Nag Expo then ex power must contain the power a of the exponential.
    double *rss, 	// the residual sum of squares for the fitted model. 
    double *df,		// the degrees of freedom associated with the residual sum of squares for the fitted model.
    double b[],		// the estimates of the parameters of the generalized linear model.
    int *rank,		// the rank of the independent variables.
    double se[],	// the standard errors of the linear parameters.
    double cov[],	// the estimate variance-covariance of estimate parameters.
    double v[],		// auxiliary information on the fitted model.
    int tdv,		// the second dimension of the array v.
    double tol,		// indicates the accuracy required for the fit of the model.
    int max_iter,	// the maximum number of iterations for the iterative weighted least-squares.
    int print_iter, // indicates if the printing of information on the iterations is required which printing is produced. 
    const char *outfile, // a null terminated character string giving the name of the file to which results should be printed.
    double eps		// the value of eps is used to decide if the independent variables are of full rank 
  ); // Linear Model with Normal Errors, Residual, Leverage.
  
/** g02gbc
		fits a generalized linear model with binomial errors.
	
Example:
	A linear trend (x = -1, 0, 1) is fitted to data relating the incidence of carriers of 
	Streptococcus pyogenes to size of tonsils. The data is described in Cox (1983).

void test_nag_glm_binomial()
{
	out_str("You must create a Text.txt file under C:\C\ folder first.");
	out_str("Otherwise, it will say this function nag_glm_binomial has some problem");
	char linkc, meanc, weightc;
	Nag_Link link;
	Nag_IncludeMean mean;
	int i, j, m, n, nvar;
	double beta[2];
	matrix v;
	v.SetSize(3,8);
	double wt[3];
	double *wtptr;
	double *offsetptr;
	int max_iter, print_iter;
	double tol, eps;
	int rank;
	double df, dev;
	double se[2], cov[3];
	String str;	
	int success;
	linkc = 'g';
	meanc = 'm';
	weightc = 'n';
	offsetptr = NULL;
	n = 3;
	m = 1;
	print_iter = 0;
	wtptr = NULL;

	link = Nag_Logistic;
	mean = Nag_MeanInclude;
	matrix x =  {
				{1.0, 0},
				{0.0, 0},
				{-1.0, 0}
				}
	double y[3] = {19.0, 29.0, 24.0};
	double binom[3] = {516.0, 560.0, 293.0};
	int ivar[2] = {1, 0};	

	nvar = 0;
	for (i=0; i<m; i++)
		if (ivar[i]>0)
			 nvar = nvar + 1;
	if (mean == Nag_MeanInclude)
		nvar =nvar + 1;
	max_iter = 10;
	tol = 5e-5;
	eps = 1e-6;
	str = "c:\\c\\c";
	success = nag_glm_binomial(link, mean, n, x, 2, m, ivar, nvar, y, binom, wtptr, 
								offsetptr, &dev, &df, beta, &rank, se, cov, v, 8, tol, max_iter, 
								print_iter, str, eps);
	if(success == NE_NOERROR || success == NE_SVD_NOT_CONV ||
		success == NE_LSQ_ITER_NOT_CONV ||
		success == NE_RANK_CHANGED || success == NE_ZERO_DOF_ERROR)
	{
		printf("\nDeviance = %12.4e\n", dev);	
		printf("Degrees of freedom = %3.1f\n\n", df);
		
		printf("Estimate       Standard error\n\n");
		for (i=0; i<nvar; i++)
			printf("%7.4f%14.4f\n", beta[i], se[i]);

		printf("\n");
		printf("binom   y   fitted value    Residual  Leverage\n\n");
		for (i = 0; i < n; ++i)
			printf("%5.1f%7.1f%10.2f%12.4f%10.3f\n", binom[i], y[i],v[i][1], v[i][4], v[i][5]);
	}
	else
		printf("This function nag_glm_binomial has some problem\n");	
}
	
	The output is as follows:
		
	You must create a Text.txt file under C:\C\ folder first.
	Otherwise, it will say this function nag_glm_binomial has some problem

	Deviance =  7.3539e-002
	Degrees of freedom = 1.0

	Estimate       Standard error

	-2.8682        0.1217
	-0.4264        0.1598

	binom   y   fitted value    Residual  Leverage

	516.0   19.0     18.45      0.1296     0.769
	560.0   29.0     30.10     -0.2070     0.422
	293.0   24.0     23.45      0.1178     0.809

Return:
	The function returns NAG error code or 0 if no error.
	
	70:	On entry parameter link had an illegal value. On entry parameter mean had an illegal value.
	11:	On entry, n must not be less than 2: n = _value_.  On entry, m must not be less than 1: m = _value_.  On entry, ip must not be less than 1: ip = _value_.  On entry, max_iter must not be less than 0: max_iter = _value_.  On entry, sx[_value_] must not be less than 0: sx[_value_] = _value_.
	5: On entry, tol must not be less than 0.0: tol = _value_.  On entry, eps must not be less than 0.0: eps = _value_.  On entry, wt[_value_] must not be less than 0.0: wt[_value_] = _value_.  On entry, binom t[_value_] must not be less than 0.0: binom t[_value_] = _value_.  On entry, y[_value_] must not be less than 0.0: y[_value_] = _value_.
	17: On entry tdx = _value_ while m = _value_. These parameters must satisfy tdx = m.  On entry tdv = _value_ while ip = _value_. These parameters must satisfy tdv = ip+6.
	25: On entry y[_value_] = _value_ while binom t[_value_] = _value_. These parameters must satisfy y[_value_] = binom t[_value_] .
	457: Parameter ip is incompatible with parameters mean and sx.
	458: Parameter ip is greater than the effective number of observations.
	460: A fitted value is at a boundary i.e., 0.0 or 1.0. This may occur if there are y values of 0.0 or binom t and the model is too complex for the data. The model should be reformulated with, perhaps, some observations dropped.
	73: Memory allocation failed.
	427: The singular value decomposition has failed to converge.
	463: The iterative weighted least-squares has failed to converge in max_iter = _value_ iterations.  The value of max_iter could be increased but it may be advantageous to examine the convergence using the print iter option. This may indicate that the convergence is slow because the solution is at a boundary in which case it may be better to reformulate the model.
	464: The rank of the model has changed during the weighted least-squares iterations. The estimate for ß returned may be reasonable, but the user should check how the deviance has changed during iterations.
	465: The degrees of freedom for error are 0. A saturated model has been fitted.
	78: Cannot open file _string_ for appending.
	79: Cannot close file _string_.
	
	successful call of the nag_glm_binomial function.
	
*/

  int nag_glm_binomial(
    Nag_Link link,
    Nag_IncludeMean mean, // indicates if a mean term is to be included.
    int n, 				// the number of observations
    double x[],			// contain the ith observation for the jth independent variable.
    int tdx,			// the second dimension of the array x.
    int m,				// the total number of independent variables.
    int sx[],			// indicates which independent variables are to be included in the model.
    int ip, 			// the number p of independent variables in the model, including the mean or intercept if present.
    double y[],			// observations on the dependent variable.
    double binom_t[],	// the binomial denominator
    double wt[],		// an (optional) weight is specified to be used with the model.
    double offset[],	// if an offset is required then offset must contain the values of the offset. 
    double *dev,		// the deviance for the fitted model.
    double *df,			// the degrees of freedom associated with the deviance for the fitted model.
    double b[],			// the coefficient of the variable contained in column j of x.
    int *rank,			// the rank of the independent variables.
    double se[],		// the standard errors of the linear parameters.
    double cov[],		// the estimate variance-covariance matrix of parameter estimates.
    double v[],			// auxiliary information on the fitted model.
    int tdv,			// the second dimension of the array v.
    double tol,			// indicates the accuracy required for the fit of the model.
    int max_iter,		// the maximum number of iterations for the iterative weighted least-squares.
    int print_iter,		// indicates if the printing of information on the iterations is required.
    const char *outfile, // a null terminated character string giving the name of the file to which results should be printed.
    double eps 		//the value of eps is used to decide if the independent variables are of full rank.
  ); // Binomial, Residual, Leverage...
  
/** g02gcc
		fits a generalized linear model with Poisson errors.
	
  
Example:
	A 3 by 5 contingency table given by Plackett (1974) is analysed by fitting terms 
	for rows and columns. The table is:

	141		67		114		79		39
	131 	66 		143 	72 		35
	36 		14 		38 		28 		16
	
void test_nag_glm_poisson()
{

	char linkc, meanc, weightc;
	Nag_Link link;
	Nag_IncludeMean mean;
	int i, j, m, n, ip;
	double ex_power, scale;
	int  sx[9];
	double b[9];
	matrix v;
	v.SetSize(15,15);
	double wt[15];
	matrix x;
	x.SetSize(15,9);
	double y[15];
	double *wtptr;
	double  *offsetptr;
	int max_iter;
	int print_iter;
	double eps;
	double tol;
	double df;
	double dev;
	int rank;
	double se[9], cov[45];
	int success;
	offsetptr = NULL;
	linkc = 'l';
	meanc = 'm';
	weightc = 'n';
	n = 15;
	m = 8;
	print_iter = 0;
	
	set_enum(linkc, &link, meanc, &mean);

	wtptr = NULL;
	
	x[0][0] = 1.0;
	x[0][1] = 0.0;
	x[0][2] = 0.0;
	x[0][3] = 1.0;
	x[0][4] = 0.0;
	x[0][5] = 0.0;
	x[0][6] = 0.0;
	x[0][7] = 0.0;
	y[0] = 141.0;
	
	x[1][0] = 1.0;
	x[1][1] = 0;
	x[1][2] = 0;
	x[1][3] = 0;
	x[1][4] = 1;
	x[1][5] = 0;
	x[1][6] = 0;
	x[1][7] = 0;
	y[1] = 67;

	x[2][0] = 1;
	x[2][1] = 0;
	x[2][2] = 0;
	x[2][3] = 0;
	x[2][4] = 0;
	x[2][5] = 1;
	x[2][6] = 0;
	x[2][7] = 0;
	y[2] = 114;
	
	x[3][0] = 1;
	x[3][1] = 0;
	x[3][2] = 0;
	x[3][3] = 0;
	x[3][4] = 0;
	x[3][5] = 0;
	x[3][6] = 1;;
	x[3][7] = 0;
	y[3] = 79;
	
	x[4][0] = 1.0;
	x[4][1] = 0;
	x[4][2] = 0;
	x[4][3] = 0;
	x[4][4] = 0;
	x[4][5] = 0;
	x[4][6] = 0;
	x[4][7] = 1;
	y[4] =  39;
		 
	x[5][0] = 0;
	x[5][1] = 1;
	x[5][2] = 0;
	x[5][3] = 1;
	x[5][4] = 0;
	x[5][5] = 0;
	x[5][6] = 0;
	x[5][7] = 0;
	y[5] =131;


	x[6][0] = 0;
	x[6][1] = 1;
	x[6][2] = 0;
	x[6][3] = 0;
	x[6][4] = 1;
	x[0][5] = 0;
	x[6][6] = 0;
	x[6][7] = 0;
	y[6] = 66;

	x[7][0] = 0;
	x[7][1] = 1;
	x[7][2] = 0;
	x[7][3] = 0;
	x[7][4] = 0;
	x[7][5] = 1;
	x[7][6] = 0;
	x[7][7] = 0;
	y[7] = 143;
	
	x[8][0] = 0;
	x[8][1] = 1;
	x[8][2] = 0;
	x[8][3] = 0;
	x[8][4] = 0;
	x[8][5] = 0;
	x[8][6] = 1;
	x[8][7] = 0;
	y[8] = 72;
	
	x[9][0] = 0;
	x[9][1] = 1;
	x[9][2] = 0;
	x[9][3] = 0;
	x[9][4] = 0;
	x[9][5] = 0;
	x[9][6] = 0;
	x[9][7] = 1;
	y[9] = 35;

	x[10][0] = 0;
	x[10][1] = 0;
	x[10][2] = 1;
	x[10][3] = 1;
	x[10][4] = 0;
	x[10][5] = 0;
	x[10][6] = 0;
	x[10][7] = 0;
	y[10] = 36;	

	x[11][0] = 0;
	x[11][1] = 0;
	x[11][2] = 1;
	x[11][3] = 0;
	x[11][4] = 1;
	x[11][5] = 0;
	x[11][6] = 0;
	x[11][7] = 0;
	y[11] = 14;
	
	x[12][0] = 0;
	x[12][1] = 0;
	x[12][2] = 1;
	x[12][3] = 0;
	x[12][4] = 0;
	x[12][5] = 1;
	x[12][6] = 0;
	x[12][7] = 0;
	y[12] =38;

	x[13][0] = 0;
	x[13][1] = 0;
	x[13][2] = 1;
	x[13][3] = 0;
	x[13][4] = 0;
	x[13][5] = 0;
	x[13][6] = 1;
	x[13][7] = 0;
	y[13] = 28;

	x[14][0] = 0; 
	x[14][1] = 0;
	x[14][2] = 1;
	x[14][3] = 0;
	x[14][4] = 0;
	x[14][5] = 0;
	x[14][6] = 0;
	x[14][7] = 1;
	y[14] =16;
	
	for (j=0; j<m; j++)
		sx[j] = 1;
	ip = 0;
	for (j=0; j<m; j++)
		if (sx[j]>0) ip =ip + 1;
	if (mean == Nag_MeanInclude)
		ip = ip + 1;
	
	ex_power = 0.0;
	max_iter = 10;
	tol = 5e-5;
	eps = 1e-6;
	success = nag_glm_poisson(link, mean, n, x, 9, m,sx, ip, y,
	wtptr, offsetptr, ex_power, &dev, &df, b, &rank, se, cov,
	v, 15, tol, max_iter, print_iter, "c:\\c\\tre.txt", eps);
	if (success == NE_NOERROR || success == NE_LSQ_ITER_NOT_CONV ||
	success == NE_RANK_CHANGED || success == NE_ZERO_DOF_ERROR)
	{
		printf("\nDeviance = %12.4e\n", dev);
		printf("Degrees of freedom = %3.1f\n\n", df);
		printf(" Estimate Standard error\n\n");
		for (i=0; i<ip; i++)
			printf("%14.4f%14.4f\n", b[i], se[i]);
		printf("\n");
		printf("   y    fitted value   Residual   Leverage\n\n");
		for (i = 0; i < n; ++i)
			printf("%7.1f%10.2f%12.4f%10.3f\n", y[i], v[i][1], v[i][4],	v[i][5]);
		
	}
	else
		printf("the function is OK");
}

int set_enum(char linkc, Nag_Link *link, char meanc, Nag_IncludeMean *mean)
{
	if(linkc == 'E' ||linkc =='e')
		*link = Nag_Expo;
	else if(linkc == 'I' ||linkc =='i')
		*link = Nag_Iden;
	else if(linkc == 'L' ||linkc =='l')
		*link = Nag_Log;
	else if(linkc == 'S' ||linkc =='s')
		*link = Nag_Sqrt;
	else if(linkc == 'R' ||linkc =='r')
		*link = Nag_Reci;
	else
		return 1;
	if(meanc == 'M'|| meanc == 'm')
		*mean = Nag_MeanInclude;
	else if(meanc == 'Z' || meanc == 'z')
		*mean = Nag_MeanZero;
	else
		return 1;
		return 0;
}

	The output is as follows:
		
	Deviance =  9.0379e+000
	Degrees of freedom = 8.0

	Estimate Standard error

    	    2.5977        0.0258
        	1.2619        0.0438
        	1.2777        0.0436
        	0.0580        0.0668
        	1.0307        0.0551
        	0.2910        0.0732
        	0.9876        0.0559
        	0.4880        0.0675
       	   -0.1996        0.0904

	y    fitted value    Residual   Leverage

	141.0    132.99      0.6875     0.604
 	 67.0     63.47      0.4386     0.514
	114.0    127.38     -1.2072     0.596
 	 79.0     77.29      0.1936     0.532
 	 39.0     38.86      0.0222     0.482
	131.0    135.11     -0.3553     0.608
 	 66.0     64.48      0.1881     0.520
	143.0    129.41      1.1749     0.601
 	 72.0     78.52     -0.7465     0.537
 	 35.0     39.48     -0.7271     0.488
 	 36.0     39.90     -0.6276     0.393
 	 14.0     19.04     -1.2131     0.255
 	 38.0     38.21     -0.0346     0.382
 	 28.0     23.19      0.9675     0.282
 	 16.0     11.66      1.2028     0.206

Return:
	The function returns NAG error code or 0 if no error.
	
	70: On entry, parameter link had an illegal value.  On entry, parameter mean had an illegal value.
	11: On entry, n must not be less than 2: n = _value_.  On entry, m must not be less than 1: m = _value_.  On entry, ip must not be less than 1: ip = _value_.  On entry, sx[_value_] must not be less than 0: sx[_value_] = _value_.  On entry, max_iter must not be less than 0: max_iter = _value_.
	17: On entry, tdx = _value	 while m = _value_. These parameters must satisfy tdx = m.  On entry, tdv = _value_ while ip = _value_. These parameters must satisfy tdv = ip+6.
	5: On entry, tol must not be less than 0.0: tol = _value_.  On entry, eps must not be less than 0.0: eps = _value_.  On entry, wt[_value_] must not be less than 0.0: wt[_value_] = _value_.  On entry, y[_value_] must not be less than 0.0: y[_value_] = _value_.
	37: On entry ex power = 0.0, link = Nag Expo. These parameters must satisfy link == Nag Expo && ex power _= 0.0.
	73: Memory allocation failed.
	457: Parameter ip is incompatible with mean and sx.		
	458: Parameter ip is greater than the effective number of observations.
	427: The singular value decomposition has failed to converge.
	461: A fitted value is at a boundary, i.e., ˆµ = 0.0. This may occur if there are y values of 0.0 and the model is too complex for the data. The model should be reformulated with, perhaps, some observations dropped.
	463: The iterative weighted least-squares has failed to converge in max_iter = _value_ iterations.  The value of max_iter could be increased but it may be advantageous to examine the convergence using the print iter option. This may indicate that the convergence is slow because the solution is at a boundary in which case it may be better to reformulate the model.
	464: The ranko f the model has changed during the weighted least-squares iterations. The estimate for ß returned may be reasonable, but the user should checkho w the deviance has changed during iterations.
	465: The degrees of freedom for error are 0. A saturated model has been fitted.
	78: Cannot open file _string for appending.
	79: Cannot close file _string	.
	
	successful call of the nag_glm_poisson function.
	
*/

  int nag_glm_poisson(
    Nag_Link link,
    Nag_IncludeMean mean, // indicates if a mean term is to be included.
    int n,			// the number of observations
    double x[],		// contain the ith observation for the jth independent variable.
    int tdx,		// the second dimension of the array x.
    int m,			// the total number of independent variables.
    int sx[],		// indicates which independent variables are to be included in the model.
    int ip,			// the number p of independent variables in the model, including the mean or intercept if present. 
    double y[],		// observations on the dependent variable.
    double wt[],	// an (optional) weight is specified to be used with the model.
    double offset[], // if an offset is required then offset must contain the values of the offset. 
    double ex_power, // if link = Nag Expo then ex power must contain the power a of the exponential.
    double *dev,	// the deviance for the fitted model.
    double *df,		// the degrees of freedom associated with the deviance for the fitted model.
    double b[],		// the estimates of the parameters of the generalized linear model.
    int *rank,		// the rank of the independent variables.
    double se[],	// the standard errors of the linear parameters.
    double cov[],	// the estimate variance-covariance of parameter estimates.
    double v[],		// auxiliary information on the fitted model.
    int tdv,		// the second dimension of the array v.
    double tol,		// indicates the accuracy required for the fit of the model.
    int max_iter,	// the maximum number of iterations for the iterative weighted least-squares.
    int print_iter, // indicates if the printing of information on the iterations is required.
    const char *outfile, // a null terminated character string giving the name of the file to which results should be printed. 
    double eps 	// the value of eps is used to decide if the independent variables are of full rank.
  ); // Poisson, Parameter Estimates, Standard Error, Residual, Leverage...
  
/** g02gdc
		fits a generalized linear model with gamma errors.
	

Example:
	A set of 10 observations from two groups are input and a model for the two groups is fitted.
	
	
void test_nag_glm_gamma()
{

	char linkc, meanc, weightc;
	Nag_Link link;
	Nag_IncludeMean mean;
	int i, j, m, n, ip;
	double ex_power, scale;
	int sx[2];
	double b[2];
	matrix v;
	v.SetSize(10,8);
	double wt[10];
	double *wtptr;
	double *offsetptr;
	int max_iter;
	int print_iter;
	double eps;
	double tol;
	double df;
	double dev;
	int rank;
	double se[2], cov[3];
	int success;
	offsetptr = NULL;
	linkc = 'r';
	meanc = 'm';
	weightc = 'u';
	n = 10;
	m = 1;
	print_iter = 0;
	scale = 0.0;	
	set_enum(linkc, &link, meanc, &mean);

	wtptr = NULL;

	matrix x =  {{1.0, 0.0}, {1.0, 0.0}, {1.0, 0.0},
				 {1.0, 0.0}, {1.0, 0.0}, {0.0, 0.0},
				 {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0},
				 {0.0, 0.0}};
	
	double y[10] = {1.0, 0.3, 10.5, 9.7, 10.9, 
					0.62, 0.12, 0.09, 0.50, 2.14};
	
	for (j=0; j<m; j++)
		sx[j] = 1;
		
	ip = 0;
	for (j=0; j<m; j++)
		if (sx[j]>0) ip = ip + 1;
	if (mean == Nag_MeanInclude)
		ip = ip + 1;
	ex_power = 0.0;
	max_iter = 10;
	tol = 5e-5;
	eps = 1e-6;
	success = nag_glm_gamma(link, mean, n, x, 2, m, sx, ip, y,	wtptr, offsetptr, &scale, ex_power, &dev, &df, b, &rank,
	se, cov, v, 8, tol, max_iter, print_iter, "c:\\c\\test.txt", eps);
	if (success == NE_NOERROR || success == NE_LSQ_ITER_NOT_CONV ||
	success == NE_RANK_CHANGED || success == NE_ZERO_DOF_ERROR)
	{
		printf("\nDeviance = %12.4e\n", dev);		
		printf("Degrees of freedom = %3.1f\n\n", df);
		printf("     Estimate     Standard error\n\n");
		for (i=0; i<ip; i++)
			printf("%14.4f%14.4f\n", b[i], se[i]);
		printf("\n");
		
		printf("y    fitted value   Residual   Leverage\n\n");
		for (i = 0; i < n; ++i)
			printf("%4.1f%10.2f%12.4f%10.3f\n", y[i], v[i][1], v[i][4], v[i][5]);
	}
	else
	printf("This function has some problem\n");
}

int set_enum(char linkc, Nag_Link *link, char meanc, Nag_IncludeMean *mean)
{
	if(linkc == 'E' ||linkc =='e')
		*link = Nag_Expo;
	else if(linkc == 'I' ||linkc =='i')
		*link = Nag_Iden;
	else if(linkc == 'L' ||linkc =='l')
		*link = Nag_Log;
	else if(linkc == 'S' ||linkc =='s')
		*link = Nag_Sqrt;
	else if(linkc == 'R' ||linkc =='r')
		*link = Nag_Reci;
	else
		return 1;
	if(meanc == 'M'|| meanc == 'm')
		*mean = Nag_MeanInclude;
	else if(meanc == 'Z' || meanc == 'z')
		*mean = Nag_MeanZero;
	else
		return 1;
		return 0;
}
	
	The output is as follows:

	Deviance =  3.5034e+001
	Degrees of freedom = 8.0

	Estimate	   Standard error
	
     1.4408        0.6678
    -1.2865        0.6717

	y    fitted value   Residual   Leverage

 	 1.0      6.48     -1.3909     0.200
 	 0.3      6.48     -1.9228     0.200
	10.5      6.48      0.5236     0.200
 	 9.7      6.48      0.4318     0.200
	10.9      6.48      0.5678     0.200
 	 0.6      0.69     -0.1107     0.200
 	 0.1      0.69     -1.3287     0.200
 	 0.1      0.69     -1.4815     0.200
 	 0.5      0.69     -0.3106     0.200
 	 2.1      0.69      1.3665     0.200
	 
Return:
	The function returns NAG error code or 0 if no error.
	
	70: On entry, parameter link had an illegal value.  On entry, parameter mean had an illegal value.
	11: On entry, n must not be less than 2: n = _value_.  On entry, m must not be less than 1: m = _value_.  On entry, ip must not be less than 1: ip = _value_.  On entry, sx[_value_] must not be less than 0: sx[_value_] = _value_.  On entry, max_iter must not be less than 0: max_iter = _value_.
	17: On entry, tdx = _value_ while m = _value_. These parameters must satisfy tdx = m.  On entry, tdv = _value_ while ip = _value_. These parameters must satisfy tdv = ip+6.
	5: On entry, scale must not be less than 0.0: scale = _value_.  On entry, tol must not be less than 0.0: tol = _value_.  On entry, eps must not be less than 0.0: eps = _value_.  On entry, wt[_value_] must not be less than 0.0: wt[_value_] = _value_.  On entry, y[_value_] must not be less than 0.0: y[_value_] = _value_.
	37: On entry ex power = 0.0 , link = Nag Expo.  These parameters must satisfy link == Nag Expo && ex power _= 0.0.
	73: Memory allocation failed.
	457: ip is incompatible with mean and sx.
	458: ip is greater than the effective number of observations.
	427: The singular value decomposition has failed to converge.
	462: A fitted value is at a boundary, i.e., ˆµ = 0.0. This may occur if there are small values of y and the model is not suitable for the data. The model should be reformulated with, perhaps, some observations dropped.
	463: The iterative weighted least-squares has failed to converge in max_iter = _value_ iterations.  The value of max_iter could be increased but it may be advantageous to examine the convergence using the print iter option. This may indicate that the convergence is slow because the solution is at a boundary in which case it may be better to reformulate the model.
	464: The rank of the model has changed during the weighted least-squares iterations. The estimate for ß returned may be reasonable, but the user should check how the deviance has changed during iterations.
	465: The degrees of freedom for error are 0. A saturated model has been fitted.
	78:Cannot open file _string	 for appending.
	79: Cannot close file _string	.
	
	successful call of the nag_glm_gamma function.
	
	 
*/

  int nag_glm_gamma(
    Nag_Link link,
    Nag_IncludeMean mean,
    int n,		// the number of observations.
    double x[], // contain the ith observation for the jth independent variable.
    int tdx,	// the second dimension of the array x.
    int m,		// the total number of independent variables.
    int sx[],	// indicates which independent variables are to be included in the model.
    int ip,		// the number p of independent variables in the model, including the mean or intercept if present.
    double y[], // observations on the dependent variable.
    double wt[], // an (optional) weight is specified to be used with the model.
    double offset[], // if an offset is required then offset must contain the values of the offset 
    double *scale, // the scale parameter for the gamma model.
    double ex_power, // if link = Nag Expo then ex power must contain the power a of the exponential.
    double *dev, // the deviance for the fitted model.
    double *df, // the degrees of freedom associated with the deviance for the fitted model.
    double b[], // the estimates of the parameters of the generalized linear model.
    int *rank, 	// the rank of the independent variables.
    double se[], // the standard errors of the linear parameters.
    double cov[], // the variance-covariance matrix of the parameter estimates.
    double v[], // auxiliary information on the fitted model.
    int tdv,	// the second dimension of the array v.
    double tol, // indicates the accuracy required for the fit of the model.
    int max_iter, // the maximum number of iterations for the iterative weighted least-squares.
    int print_iter, // indicates if the printing of information on the iterations is required.
    const char *outfile, // a null terminated character string giving the name of the file to which results should be printed.
    double eps  // the value of eps is used to decide if the independent variables are of full rank.
  ); // Gamma, Parameter Estimates, Standard Error, Residual, Leverage...
  
/** g02gkc
		calculates the estimates of the parameters of a generalized linear model for given 
		constraints from the singular value decomposition results.
	
Example:
	A loglinear model is fitted to a 3 by 5 contingency table by nag_glm_poisson (g02gcc). 
	The model consists of terms for rows and columns. The table is:

	141   67 	114 	79 	  	39
	131   66 	143 	72 		35
	 36   14 	 38 	28 		16

void test_nag_glm_tran_model()
{

	int i, nclin, ip, j, m, n;
	double ex_power;
	int sx[9];
	double b[9];
	matrix c;
	c.SetSize(15,9);
	double cov[45], se[9];
	matrix v;
	v.SetSize(15,15);
	double *wtptr;
	int max_iter;
	int print_iter;
	double eps;
	double tol;
	double df;
	double dev;
	int rank;
	int success;
	
	n = 15;
	m = 8;
	print_iter = 0;
	wtptr = NULL;

	double y[15] = {141.0, 67, 114, 79, 39, 131, 66, 
					143, 72, 35, 36, 14, 38, 28, 16};

	matrix x =  {{1, 0, 0, 1, 0, 0, 0, 0, 0}, {1, 0, 0, 0, 1, 0, 0, 0, 0},
				 {1, 0, 0, 0, 0, 1, 0, 0, 0}, {1, 0, 0, 0, 0, 0, 1, 0, 0},
				 {1, 0, 0, 0, 0, 0, 0, 1, 0}, {0, 1, 0, 1, 0, 0, 0, 0, 0},
				 {0, 1, 0, 0, 1, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 1, 0, 0, 0},
				 {0, 1, 0, 0, 0, 0, 1, 0, 0}, {0, 1, 0, 0, 0, 0, 0, 1, 0},
			 	 {0, 0, 1, 1, 0, 0, 0, 0, 0}, {0, 0, 1, 0, 1, 0, 0, 0, 0},
				 {0, 0, 1, 0, 0, 1, 0, 0, 0}, {0, 0, 1, 0, 0, 0, 1, 0, 0},
				 {0, 0, 1, 0, 0, 0, 0, 1, 0}};
	
	double c1[18] = { 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0,
					1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0 };
			
	
	for (j=0; j<m; j++)
		sx[j]= 1;
	ip = 9;
	
	max_iter = 10;
	tol = 5e-5;
	eps = 1e-6;
	ex_power = 0.0;
	success = nag_glm_poisson(Nag_Log, Nag_MeanInclude, n, x, 9, m, sx, ip, y, wtptr, NULL, ex_power, &dev, 
					&df, b, &rank, se, cov,	v, 15, tol, max_iter, print_iter, "c:\\c\\test.txt", eps);
	if (success == NE_NOERROR || success == NE_LSQ_ITER_NOT_CONV ||
		success == NE_RANK_CHANGED || success == NE_ZERO_DOF_ERROR)
	{
		printf("\nDeviance = %12.4e\n", dev);
		printf("Degrees of freedom = %3.1f\n\n", df);
		nclin = ip - rank;
			
		for (i=0; i<ip; ++i)
			for (j=0; j<nclin; ++j)
				c[i][j] = c1[i * nclin + j];
		success = nag_glm_tran_model(ip, nclin, v, 15, c, 9, b, 1.0e0, se, cov);
		if (success == NE_NOERROR)
		{
			printf("Estimate        Standard error\n\n");
			for (i=0; i<ip; i++)
				printf(" %7.4f%14.4f\n", b[i], se[i]);
			printf("\n");
		}
		else
			printf("some problem here\n");	
	}
	else
		printf("some problem here\n");
}

	The output is as follows:
	
	Deviance =  9.0379e+000
	Degrees of freedom = 8.0

	Estimate        Standard error

  	  3.9831        0.0396
  	  0.3961        0.0458
  	  0.4118        0.0457
 	 -0.8079        0.0622
  	  0.5112        0.0562
 	 -0.2285        0.0727
  	  0.4680        0.0569
 	 -0.0316        0.0675
 	 -0.7191        0.0887
Return:
	The function returns NAG error code or 0 if no error.
	
	11: On entry, ip must not be less than 1: ip = _value_.
	12: On entry, nclin must not be less than or equal to 0: nclin = _value_.
	20: On entry, nclin = _value_ while ip = _value_. These parameters must satisfy nclin < ip.
	17: On entry, tdv = _value_ while ip = _value_. These parameters must satisfy tdv = ip+6.  On entry, tdc = _value_ while nclin = _value_. These parameters must satisfy tdc = nclin.
	6: On entry, scale must not be less than or equal to 0.0: scale = _value_.
	438: matrix c does not give a model of full rank.
	73: Memory allocation failed.

	successful call of the nag_glm_tran_model function.
	

*/

  int nag_glm_tran_model(
    int ip, 	// the number of terms in the linear model.
    int nclin,	// the number of constraints to be imposed on the parameters.
    double v[], // v as returned by (g02gac), (g02gbc), (g02gcc) or (g02gdc).
    int tdv,	// the second dimension of the array v 
    double c[], // the nclin constraints stored by column.
    int tdc,	// the last dimension of the array c.
    double b[], // Input: the parameter estimates computed by using the singular value decomposition.
				// Output: the parameter estimates of the parameters with the constraints imposed.
    double scale, // the estimate of the scale parameter.
    double se[], // the standard error of the parameter estimates in b.
    double cov[]  // the estimate variance-covariance of the parameter estimates.
  ); // Constraint, Estimate, Standard Error. 
  
/** g02gnc
		gives the estimate of an estimable function along with its standard error 
		from the results from fitting a generalized linear model.
	
Example:
	A loglinear model is fitted to a 3 by 5 contingency table by (g02gcc). The model consists of 
	terms for rows and columns. The table is:

	141   67 	114 	79 	  	39
	131   66 	143 	72 		35
	 36   14 	 38 	28 		16

	The number of functions to be tested is read in, then the linear functions themselves are 
	read in and tested with nag_glm_est_func. The results of nag_glm_est_func are printed.
	
	
void  test_nag_glm_est_func()
{
	Boolean est;
	int i, ip, j, m, n, nestfn, tdx, tdv;
	double ex_power;
	double sestat, stat, z;
	int sx[9];
	double b[9], f[9];
	matrix v;
	v.SetSize(15,15);
	double  wt[15];
	double *wtptr;
	int max_iter;
	int print_iter;
	double eps;
	double tol;
	double df;
	double dev;
	int rank;
	double cov[45], se[9];
	int success;
	n = 15;
	m = 8;
	print_iter = 0;
	
	double y[15] = {141.0, 67, 114, 79, 39, 131, 66, 143, 72, 35, 36, 14, 38, 28, 16};
	
	matrix x =  {{1, 0, 0, 1, 0, 0, 0, 0, 0}, {1, 0, 0, 0, 1, 0, 0, 0, 0},
				 {1, 0, 0, 0, 0, 1, 0, 0, 0}, {1, 0, 0, 0, 0, 0, 1, 0, 0},
				 {1, 0, 0, 0, 0, 0, 0, 1, 0}, {0, 1, 0, 1, 0, 0, 0, 0, 0},
				 {0, 1, 0, 0, 1, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 1, 0, 0, 0},
				 {0, 1, 0, 0, 0, 0, 1, 0, 0}, {0, 1, 0, 0, 0, 0, 0, 1, 0},
				 {0, 0, 1, 1, 0, 0, 0, 0, 0}, {0, 0, 1, 0, 1, 0, 0, 0, 0},
			 	 {0, 0, 1, 0, 0, 1, 0, 0, 0}, {0, 0, 1, 0, 0, 0, 1, 0, 0},
				 {0, 0, 1, 0, 0, 0, 0, 1, 0}};
	
	matrix fl = {{1, 1, 0, 0, 1, 0, 0, 0, 0}, {0, 1,-1, 0, 0, 0, 0, 0, 0},
				 {0, 1, 0, 0, 0, 0, 0, 0, 0}};
				
	wtptr = NULL;
	for (j=0; j<m; j++)
		sx[j]=1;
	ip = 9;
	max_iter = 10;
	tol = 5e-5;
	eps = 1e-6;
	ex_power = 0.0;
	success = nag_glm_poisson(Nag_Log, Nag_MeanInclude, n, x, 9,
					m, sx, ip, y, wtptr, NULL, ex_power, &dev, &df, b, &rank, se, cov,
					v, 15, tol, max_iter, print_iter, "c:\\c\\test.txt", eps);

	if (success == NE_NOERROR || success == NE_LSQ_ITER_NOT_CONV ||
			success == NE_RANK_CHANGED || success == NE_ZERO_DOF_ERROR)
	{
			printf("\nDeviance = %12.4e\n", dev);
			printf("Degrees of freedom = %3.1f\n\n", df);
			printf("Estimate       Standard error\n\n");
			for (i=0; i<ip; i++)
				printf("%7.4f%14.4f\n", b[i], se[i]);

			printf("\n");
			nestfn = 3;			  
			for(i = 0; i < 3; i++)
			{
				for(j = 0; j < 9; j++)
					f[j] = fl[i][j];
			
				success = nag_glm_est_func(ip, rank, b, cov, v, 15, f, &est, &stat, &sestat, &z, tol);
				if (success == NE_NOERROR || success == NE_RANK_EQ_IP)
				{
					printf("\n");
					printf("Function %ld\n\n", i+1);
					for (j=0; j<ip; ++j)
					printf("%8.2f%c", f[j], (j%5==4 || j==ip-1) ? '\n' : ' ');
					printf("\n");
					if (est)
						printf("stat = %10.4f\nsestat = %10.4f\nz = %10.4f\n", stat, sestat, z);
					else
					printf("Function not estimable\n");
				}
			}
	}
	else
	{
		printf("This function has some problem\n");
	}
}

	
	The output is as follows:
	
	Deviance =  9.0379e+000
	Degrees of freedom = 8.0

	Estimate       Standard error

 	 2.5977        0.0258
 	 1.2619        0.0438
 	 1.2777        0.0436
 	 0.0580        0.0668
 	 1.0307        0.0551
 	 0.2910        0.0732
 	 0.9876        0.0559
 	 0.4880        0.0675
	-0.1996        0.0904


	Function 1

    	1.00    1.00    0.00    0.00    1.00    0.00    0.00    0.00    0.00
	stat =     4.8903
	sestat =     0.0674
	z =    72.5934

	Function 2

    	0.00    1.00   -1.00    0.00    0.00    0.00    0.00    0.00    0.00
	stat =    -0.0158
	sestat =     0.0672
	z =    -0.2350

	Function 3

    	0.00    1.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00
	Function not estimable

Return:
	The function return NAG error code or 0 if no error.
	
	11: On entry, ip must not be less than 1: ip = _value_.  On entry, rank must not be less than 1: rank = _value_.
	19: On entry, ip = _value_ while rank = _value_. These parameters must satisfy rank = ip.
	17: On entry, tdv = _value_ while ip = _value_. These parameters must satisfy tdv = ip+6.
	440: On entry, rank = ip. In this case, the boolean variable est is returned as TRUE and all statistics are calculated.
	439: sestat, the standard error of the estimate of the function, se(F) = 0.0; probably due to rounding error or due to incorrectly specified input values of cov and f.
	73: Memory allocation failed.
	
	successful call of the nag_glm_est_func function.
	
*/

  int nag_glm_est_func(
    int ip,		// the number of terms in the linear model.
    int rank, 	// the rank of the independent variables.
    double b[], // the ip values of the estimates of the parameters of the model.
    double cov[], // the estimate variance-covariance of the parameter estimates.
    double v[],	// v as returned by (g02gac), (g02gbc), (g02gcc) or (g02gdc).
    int tdv,	// the second dimension of the array v.
    double f[], // the linear function to be estimated.
    Boolean *est, 	// indicates if the function was estimable.
    double *stat, // if est = TRUE, stat contains the estimate of the function.
    double *sestat, // if est = TRUE, sestat contains the standard error of the estimate of the function.
    double *z,	// if est = TRUE, z contains the z statistic for the test of the function being equal to zero.
    double tol // tol is the tolerance value used in the check for estimability.
  ); // Estimable Function, Parameter Estimates, Standard Errors...
  
/**  g02hac
		performs bounded influence regression (M-estimates). Several standard methods are available.
	
Example:
	The number of observations and the number of x variables are read in followed byt he data. The
	option parameters are then read in (in this case giving: Schweppe type regression with Hampel's
	Psi function and Huber's Chi function and then using the 'replace expected byo bserved' option in
	calculating the covariances). Finally a set of values for the constants are read in.  After a 
	call to nag_robust_m_regsn_estim, theta-hat, its standard error and sigma-hat are printed. In addition the
	weight and residual for each observation is printed.


void test_nag_robust_m_regsn_estim()
{	
	char regtypec, psifunc, sigestc, covmat_estc;
	Nag_RegType regtype;
	Nag_PsiFun psifun;
	Nag_SigmaEst sigma_est;
	Nag_CovMatrixEst covmat_est;
	int i, j, m, n;
	double sigma;
	matrix c;
	c.SetSize(3,3);
	double theta[3];
	double cpsi;
	double hpsi[3];
	double cucv;
	double dchi;
	double rs[8];
	double wt[8];
	double tol;
	int max_iter;
	int print_iter;
	double info[4];
	int success;
    String str;
	m = 3;
	n = 8;
	
	matrix x = 	{{1, -1, -1}, {1, -1, 1}, {1, 1, -1},
				 {1, 1, 1}, {1, -2, 0}, {1, 0, -2},
				 {1, 2, 0}, {1, 0, 2}};

	double y[8] = {2.1, 3.6, 4.5, 6.1, 1.3, 1.9, 6.7, 5.5};

	regtypec = 'S';
	psifunc = 'P';
	sigestc = 'X';	

	set_enum1(regtypec, &regtype, psifunc, &psifun, sigestc, &sigma_est);

	if (regtype != Nag_HuberReg)
	{
		covmat_estc = 'R';
		cucv = 3.0;
	}	
	if (covmat_estc == 'A')
		covmat_est = Nag_CovMatAve;
	else
		covmat_est = Nag_CovMatObs;
	if (psifun != Nag_Lsq)
	{
		if (psifun == Nag_HuberFun)
			cpsi = 1.5;
		else
			cpsi = 1.5;
		if (psifun == Nag_HampelFun)
		{	
			hpsi[0] = 1.5;
			hpsi[1] = 3.0;
			hpsi[2] = 4.5; 
		}
		if (sigma_est == Nag_SigmaChi)
			dchi= 1.5;
	}


	printf("dchi =%f \n", dchi);

	tol = 5e-5;
	max_iter = 50;
	print_iter = 1;
	sigma = 1.0e0;
	for (i = 0; i < m; ++i)
		theta[i] = 0.0e0;
	str = "c:\\c\\teax.txt";
	
	success = nag_robust_m_regsn_estim(regtype, psifun, sigma_est, covmat_est, n, m, x,
										3, y, cpsi, hpsi, cucv, dchi, theta, &sigma, c, 3, rs,
									    wt, tol, max_iter, print_iter, str, info);

	if(success == 0)
	{
	}
	else
	{
		printf("THE PRINTF");
	}


	if((success == NE_NOERROR) || (success == NE_THETA_ITER_EXCEEDED) ||
	(success== NE_LSQ_FAIL_CONV) || (success == NE_MAT_SINGULAR) ||
	(success == NE_WT_LSQ_NOT_FULL_RANK) ||
	(success == NE_REG_MAT_SINGULAR) ||
	(success == NE_COV_MAT_FACTOR_ZERO) ||
	(success == NE_VAR_THETA_LEQ_ZERO) ||
	(success == NE_VAR_THETA_LEQ_ZERO) ||
	(success == NE_ERR_DOF_LEQ_ZERO) ||
	(success == NE_ESTIM_SIGMA_ZERO))
	{
		if (success != NE_NOERROR)
		{
			
			printf(" Some of the following results may be unreliable\n");
		}
		printf("Sigma = %10.4f\n\n", sigma);
		printf("Theta         Standard errors\n\n");
		for (j = 0; j < m; ++j)
			printf("%6.4f %13.4f\n", theta[j], c[j][j]);

		printf("\nWeights       Residuals\n\n");
		for (i = 0; i < n; ++i)
			printf("%6.4f %13.4f\n", wt[i], rs[i]);
	}
}

int set_enum1(char regtypec, Nag_RegType *regtype, char psifunc, Nag_PsiFun *psifun,
			 char sigestc, Nag_SigmaEst *sigma_est)
{
	if(regtypec == 'h'||regtypec == 'H')
		*regtype = Nag_HuberReg;
	else if(regtypec == 'M'||regtypec == 'm')
		*regtype = Nag_MallowsReg;
	else if(regtypec == 'S'||regtypec == 's')
		*regtype = Nag_SchweppeReg;	
	else
		return 1;
	
	
	if(psifunc == 'L' || psifunc =='l')
		*psifun = Nag_Lsq;
	else if(psifunc == 'H' || psifunc =='h')
		*psifun = Nag_HuberFun;
		
	else if(psifunc == 'p' || psifunc =='P')
		*psifun = Nag_HampelFun;	
	else if(psifunc == 'A' || psifunc =='a')
		*psifun = Nag_AndrewFun;
	else if(psifunc == 'T' || psifunc =='t')
		*psifun = Nag_TukeyFun;
	else
		return 1;
		
	if(sigestc == 'R' || sigestc == 'r')
		*sigma_est = Nag_SigmaRes;
	else if(sigestc == 'C' || sigestc == 'c')
		*sigma_est = Nag_SigmaConst;
	else if(sigestc == 'X' || sigestc == 'x')
		*sigma_est = Nag_SigmaChi;
	else 
		return 1;
	return 0;
		
	
}

	
	The output is as follows:
	
	dchi =1.500000 
	Sigma =     0.2026

	Theta         Standard errors

	4.0423        0.0384
	1.3083        0.0272
	0.7519        0.0311

	Weights       Residuals

	0.5783        0.1179
	0.5783        0.1141
	0.5783       -0.0987
	0.5783       -0.0026
	0.4603       -0.1256
	0.4603       -0.6385
	0.4603        0.0410
	0.4603       -0.0462

	
Return:
	The function returns NAG error code or 0 if no error.
	
	11: On entry, n must not be less than 2: n =_value_.  On entry, m must not be less than 1: m =_value_.
	20: On entry, m =_value_ while n =_value_. These parameters must satisfy m < n.
	17: On entry, tdx =_value_ while m =_value_. These parameters must satisfy tdx = m.  On entry, tdc =_value_ while m =_value_. These parameters must satisfy tdc = m.
	70: On entry, parameter regtype had an illegal value.  On entry, parameter psifun had an illegal value.  On entry, parameter sigma est had an illegal value.  On entry, parameter covmat est had an illegal value.
	12: On entry, max_iter must not be less than or equal to 0: max_iter =_value_.
	6: On entry, sigma must not be less than or equal to 0.0: sigma =_value_.  On entry, tol must not be less than or equal to 0.0: tol = value_.
	473: On entry, psifun = Nag HuberFun and cpsi =_value_. For this value of psifun, cpsi must be > 0.0.
	475: On entry, regtype = Nag MallowsReg, cucv =_value_ and m =_value_. For this value of regtype, cucv must be >= m.
	476: On entry, regtype = Nag_SchweppeReg, cucv =_value_ and m =_value_. For this value of regtype, cucv must be >= m^(1/2).
	474: On entry, psifun _= Nag Lsq, sigma est = Nag SigmaChi and dchi =_value_. For these values of psifun and sigma est, dchi must be > 0.0.
	477: On entry, psifun = Nag HampelFun and hpsi[0] =_value_, hpsi[1] =_value_ and hpsi[2] =_value_. For this value of psifun, the elements of hpsi must satisfyt he condition 0.0 = hpsi[0] = hpsi[1] = hpsi[2] and hpsi[2] > 0.0.
	478: The number of iterations required to calculate the weights exceeds max_iter.  This is only applicable if regtype _= Nag_HuberReg.
	479: The number of iterations required to calculate ß1 exceeds max_iter.  This is onlya pplicable if regtype = Nag MallowsReg and sigma est = Nag SigmaRes.
	480: The number of iterations required to calculate theta and sigma exceeds max_iter. In this case, info[2] = max_iter on exit.
	481: The iterations to solve the weighted least-squares equations failed to converge.
	482: The weighted least-squares equations are not of full rank.
	483: Failure to invert matrix while calculating covariance.  If regtype = Nag HuberReg, then (XTX) is almost singular.  If regtype _= Nag HuberReg, then S1 is singular or almost singular. This maybe due to too many diagonal elements of the matrix being zero, see Section 6.
	484: In calculating the correlation factor for the asymptotic variance-covariance matrix, the factor for covariance matrix = 0.  For this error, see Section 6. In this case c is returned as inv(t(X)*X).  (This is only applicable if regtype = Nag HuberReg).
	485: The estimated variance for an element of ? = 0. In this case the diagonal element of c will contain the negative variance and the above diagonal elements in the row and the column corresponding to the element will be returned as zero.  This error maybe caused byrounding errors or too many of the diagonal elements of p being zero. See Section 6.
	486: n = value_, rank of x =_value_. The degrees of freedom for error, n - (rank of x) must be > 0.0.
	487: The estimated value of s was 0.0 during an iteration.
	78: Cannot open file string_ for appending.
	79: Cannot close file string_.
	73: Memory allocation failed.
	
	successful call of the nag_robust_m_regsn_estim function.
	

*/
  
  int nag_robust_m_regsn_estim(
    Nag_RegType regtype,
    Nag_PsiFun psifun,
    Nag_SigmaEst sigma_est,
    Nag_CovMatrixEst covmat_est,
    int n,		// the number of observations.
    int m,		// the number m, of independent variables.
    double x[], // the values of the X matrix, i.e., the independent variables.
    int tdx,	// the second dimension of the array x.
    double y[], // the data values of the dependent variable.
    double cpsi, // if psifun = Nag HuberFun, cpsi must specify the parameter, c, of Huber's psi function.
    double hpsi[], // if psifun = Nag_HampelFun then hpsi[0], hpsi[1] and hpsi[2] must specify the parameters h1, h2, and h3, of Hampel's piecewise linear psi function. 
    double cucv, // if regtype = Nag_MallowsReg then cucv must specify the value of the constant, c, of the function u for Maronna's proposed weights. If regtype = Nag_SchweppeReg then cucv must specify the value of the function u for the Krasker-Welsch weights.
    double dchi, // the constant, d, of the Chi function.
    double theta[], // starting values of the parameter vector theta. These maybe obtained from least-squares regression.
    double *sigma, // a starting value for the estimation of sigma.
    double c[],	// the diagonal elements of c contain the estimated asymptotic standard errors of the estimates of theta,
    int tdc, 	// the second dimension of the array c.
    double rs[], // contains the residuals from the model evaluated at final value of theta.
    double wt[], // contains the vector of weights.
    double tol,	// the relative precision for the calculation of A.
    int max_iter, // the maximum number of iterations that should be used in the calculation of A.
    int print_iter, // the amount of information that is printed on each iteration.
    const char *outfile, // a null terminated character string giving the name of the .le to which results should be printed.
    double info[]	// elements of info contain the values of parameters, rank, and number of iterations.
  ); // Sigma, Theta, Weight, Standard Error....
  
/** g02hkc
		computes a robust estimate of the covariance matrix for an expected fraction of gross errors.
	
Example:
	A sample of 10ob servations on three variables is read in and the robust estimate of the covariance 
	matrix is computed assuming 10% gross errors are to be expected and file c;\c\test.txt exists. 

void test_nag_robust_corr_estim()
{
	int i, j, k, m, n;
	int ifail;
	matrix x;
	x.SetSize(20,10);
	double theta[10];
	int tdx=10;
	int max_iter, l1, l2;
	int print_iter;
	double eps, cov[15];
	int iter;
	double tol;
	x[0][0] = 3.4;
	x[0][1] =  6.9 ;
	x[0][2] = 12.2;
	x[1][0] = 6.4 ;
	x[1][1] = 2.5;
	x[1][2] = 15.1;
	x[2][0] = 4.9;
	x[2][1] = 5.5;
	x[2][2] = 14.2;
	x[3][0] = 7.3;
	x[3][1] = 1.9;
	x[3][2] = 18.2;
	x[4][0] = 8.8;
	x[4][1] = 3.6;
	x[4][2] = 11.7;
	x[5][0] = 8.4;
	x[5][1] =  1.3;
	x[5][2] =  17.9;

	x[6][0] = 5.3;
	x[6][1] = 3.1;
	x[6][2] = 15.0;
	x[7][0] = 2.7; 
	x[7][1] = 8.1;
	x[7][2] = 7.7;
	x[8][0] = 6.1; 
	x[8][1] = 3.0;
	x[8][2] = 21.9;
	x[9][0] = 5.3;
	x[9][1] = 2.2;
	x[9][2] = 13.9;

	n = 10;
	m = 3;
	eps = 0.1;
	max_iter = 100;
	tol = 5e-5;
	print_iter = 1;
	printf("The input data are as follows\n");
	printf("\nn = 10, m = 3, eps = 0.1, max_iter = 100, tol = 5e-5, print_iter = 1\n");
	printf("\nthe ouput file name = c:\\c\\test.txt\n");
	String	strFile = "c:\\c\\test.txt";
	nag_robust_corr_estim(n, m, x, tdx, eps, cov, theta, max_iter, print_iter,
							strFile, tol, &iter);
	printf("\n\ng02hkc required %ld iterations to converge\n\n", iter);
	printf("Covariance matrix\n");
	l2 = 0;
	for (j = 1; j <= m; ++j)
	{
		l1 = l2 + 1;
		l2 += j;
		for (k = l1; k <= l2; ++k)
		{
			printf("%6.3f\t", cov[k - 1]);
		}
		printf("\n");
	}
	printf("\n\ntheta\n");
	for (j = 1; j <= m; ++j)
	{
		printf("%6.3f\n", theta[j - 1]);
	}
	
}


  The output is as follows

	The input data are as follows

	n = 10, m = 3, eps = 0.1, max_iter = 100, tol = 5e-5, print_iter = 1

	the ouput file name = c:\c\test.txt


	g02hkc required 23 iterations to converge

	Covariance matrix
	  3.461	
     -3.681	   5.348	
      4.682	  -6.645	14.439


	theta
	 5.818
	 3.681
	15.037
	
Return:
	The function returns NAG error code or 0 if no error.
	
	11: On entry, n must not be less than 2: n = _value_.  On entry, m must not be less than 1: m = _value_.
	19: On entry, m = _value_ while n = _value_. These parameters must satisfy m = n.
	17: On entry, tdx = _value_ while m = _value_	. These parameters must satisfy tdx = m.
	12: On entry, max_iter must not be less than or equal to 0: max_iter = _value_.
	5: On entry, eps must not be less than 0.0: eps = _value_.
	8: On entry, eps must be not be greater than or equal to 1.0: eps = _value_.
	6: On entry, tol must not be less than or equal to 0.0: tol = _value_.
	497: On entry, column _value_ of array x has constant value.
	3: Too many iterations(_value_).  The iterative procedure to find the co-variance matrix C, has failed to converge in max_iter iterations.
	490: The iterative procedure to find C has become unstable. This may happen if the value of eps is too large.
	78: Cannot open file _string	 for appending.
	79: Cannot close file _string	.
	73: Memory allocation failed.
	
	successful call of the nag_robust_corr_estim function.
	
*/
  
  int nag_robust_corr_estim(
    int n,		// the number of observations.
    int m,		// the number of independent variables
    double x[],	// must contain the ith observation for the jth variable.
    int tdx,	// the second dimension of the array x.
    double eps, // the expected fraction of gross errors expected in the sample.
    double cov[], // the covariance matrix.
    double theta[], // the robust estimate of the location parameters.
    int max_iter, // the maximum number of iterations that will be used during the calculation of the covariance matrix.
    int print_iter, // indicates if the printing of information on the iterations is required.
    const char *outfile, // a null terminated character string giving the name of the file to which results should be printed. 
    double tol, 	// the relative precision for the final estimates of the covariance matrix.
    int *iter		// the number of iterations performed.
  ); // Robust Covariance matrix, Robust Parameter Estimates.
  /* end proto */
#ifdef __cplusplus
}
#endif
#endif /* not NAGG02 */