(*:Mathematica Version: 3.0 *)

(*:Package Version: 1.2 *)

(*:Context: Statistics`LinearRegression` *)

(*:Name: Statistics`LinearRegression` *)

(*:Title: Linear Regression Analysis *)

(*:Author: Wolfram Research, Inc. *)

(*:Copyright: Copyright 1990-1999, Wolfram Research, Inc. *)

(*:Summary:
This package provides least squares or weighted least squares linear
regression for data whose errors are assumed to be normally and
independently distributed.  It supports a number of commonly used statistics
such as RSquared, EstimatedVariance, and an ANOVATable.  It also provides
diagnostics for leverage, influence, collinearity, and correlation.
*)

(*:History:
  Version 1.0, Wolfram Research, Inc., 1990. 
  Version 1.1, ECM (Wolfram Research), 1994: mostly rewritten,
  Replaced options OutputList, OutputControl, BestFitCoefficients, and
     ConfidenceIntervalTable (and option value NoPrint) with options
     RegressionReport, BestFitParameters, and ParameterCITable (and option
     value SummaryReport);
  Added options BestFitParametersDelta, CatcherMatrix, CoefficientOfVariation,
     CookD, StudentizedResiduals, JackknifedVariance,
     CovarianceMatrixDetRatio, DurbinWatsonD, EigenstructureTable, HatDiagonal,
     MeanPredictionCITable, ParameterConfidenceRegion,
     PartialSumOfSquares, PredictedResponseDelta,
     SequentialSumOfSquares, SinglePredictionCITable, StandardizedResiduals,
     and VarianceInflation.
  Added ability for DesignedRegress to accept either an unweighted design
     matrix or the singular value decomposition of an unweighted design
     matrix.  Added IncludeConstant option to DesignMatrix.
  Added RegressionReportValues so the possible values of the option
     RegressionReport are easily discovered.
  Version 1.2, ECM (Wolfram Research), 1998: modified Weights option so it could
        be set to a pure function of the entire observation vector... not
        just the response.  A pure function of a single argument is still
	applied to the response as before.  Some speedups implemented based
	on suggestions of D. Withoff and experiments by ECM.  Modified
	StandardizedResiduals and StudentizedResiduals to take account of
	Weights.  (Note that the definitions for StandardizedResiduals and
	StudentizedResiduals for the unweighted case follow the definitions
	on p. 157 of "Modern Applied Statistics with S-Plus", Venables and
	Ripley, Springer-Verlag, 1994.  There is not uniform agreement on
	what to call these diagnostics.)
*)

(* :Discussion:

(1) Weights -> vector

  When Weights -> {w1, ..., wn}, then the
        weighted residual sum of squares is minimized:
                Sum[wi (yi - g[params, xi1, ..., xik])^2, {i, 1, n}]
        (where g[params, ] represents the linear model).

(2) Weights -> pureFcn

  (i) In Regress[{y1, ..., yn}, basisFcns, x, Weights -> (w[#] &)],
	the weight function is applied to the response y, and the weighted sum 
		Sum[w[yi] (yi - g[params, i])^2, {i, 1, n}]
	is minimized.
      In Regress[{y1, ..., yn}, basisFcns, x, Weights -> (w[#1, #2] &)],
	the weight function is applied to the response y and the implied
	independent variable x, and the weighted sum
		Sum[w[yi, i] (yi - g[params, i])^2, {i, 1, n}]
	is minimized.

  (ii) In Regress[{{x1, y1}, ..., {xn, yn}}, basisFcns, x, Weights -> (w[#] &)],
	 the weight function is applied to the response y, and the weighted sum
		Sum[w[yi] (yi - g[params, xi])^2, {i, 1, n}]
	 is minimized.
       In Regress[{{x1, y1}, ..., {xn, yn}}, basisFcns, x,
	 					Weights -> (w[#1, #2] &)],
	 the weight function is applied to the response y and the explicit
	 independent variable x, and the weighted sum
		Sum[w[yi, xi] (yi - g[params, xi])^2, {i, 1, n}]
	 is minimized.

  (iii) In Regress[{{x11, ..., x1k, y1}, ..., {xn1, ..., xnk, yn}}, basisFcns,
	  				{x1, ..., xk}, Weights -> (w[#] &)],
	  the weight function is applied to the response y, and the weighted sum
                Sum[w[yi] (yi - g[params, xi, ..., xk])^2, {i, 1, n}]
          is minimized.
	In Regress[{{x11, ..., x1k, y1}, ..., {xn1, ..., xnk, yn}}, basisFcns,
                         {x1, ..., xk}, Weights -> (w[#1, ..., #(k+1)] &)],
	  the weight function is applied to the response y and the explicit
          independent variables x1, ..., xk, and the weighted sum
		Sum[w[yi, x1i, ..., xki] (yi - g[params, x1i, ..., xki])^2,
								 {i, 1, n}]
	  is minimized.
	

  Note that while Regress will apply a pure function to the ith observation
	{yi, xi1, ..., xik}, DesignedRegress does not have
	access to the independent variables {xi1, ..., xik}, so it will
	apply the pure function to the ith response yi joined with the ith
	entry in the design matrix:
	{y1, f1[xi1, ..., xik], ..., fp[xi1, ..., xik]}.

  Warning:  Note that while DesignedRegress[mat, response] or
	DesignedRegress[svd, response] is supposed to
	save time because the design matrix or the singular value
	decomposition of the design matrix is already computed (and
	may be used with different response vectors), there is
	a sanity check done in vDR to calculate constantPresent that
	can eliminate the potential time savings.  This sanity check is
	not needed if you use Regress instead of DesignedRegress.
*)

(*:Keywords: linear regression, design matrix  *)

(*:Requirements: No special system requirements. *)

(*:Warning:
Options OutputList, OutputControl, BestFitCoefficients, and
ConfidenceIntervalTable and option value NoPrint are obsolete.

V2.2 syntax:
  Regress[d, f, v, OutputList -> Null, OutputControl -> Automatic]
  Regress[d, f, v, OutputList -> {BestFitCoefficients, ConfidenceIntervalTable},
	 OutputControl -> Automatic]
  Regress[d, f, v, OutputList -> {BestFitCoefficients, ConfidenceIntervalTable},
	 OutputControl -> NoPrint]
V3.0 syntax:
  Regress[d, f, v, RegressionReport -> SummaryReport]
  Regress[d, f, v, RegressionReport -> {BestFitParameters, ParameterCITable,
	 SummaryReport}]
  Regress[d, f, v, RegressionReport -> {BestFitParameters, ParameterCITable}]
*)

(*:Sources:
F. Mosteller & J. W. Tukey, Data Analysis and Regression, 1977, Addison-Wesley.
P. F. Velleman & R. E. Welsch, Efficient Computing of Regression Diagnostics,
	The American Statistician, 1981, Vol. 35, No. 4.
D. A. Belsey, E. Kuh, & R. E. Welsch, Regression Diagnostics, 1980, Wiley.
Weisberg, S., Applied Linear Regression, 1985, Wiley.
*)


BeginPackage["Statistics`LinearRegression`",
	      (* needed for FRatioDistribution and StudentTDistribution: *)	
             "Statistics`NormalDistribution`",
	      (* needed for CDF: *)
	     "Statistics`Common`DistributionsCommon`",
             "Statistics`ConfidenceIntervals`",
              (* needed for Mean and CoefficientOfVariation: *)
             "Statistics`DescriptiveStatistics`",
	      (* needed for Weights, RegressionReport,
		 RegressionReportValues, and all symbols allowed by
		 RegressionReportValues[Regress] and
		 RegressionReportValues[DesignedRegress] *)
	     "Statistics`Common`RegressionCommon`",
	      (* needed for CovarianceMatrix and CorrelationMatrix: *)
	     "Statistics`Common`MultivariateCommon`"]

Regress::usage =
"Regress[data, funs, vars] finds a least-squares fit to a list of data as
a linear combination of the functions funs of variables vars. The data
can have the form {{x1, y1, ..., f1}, {x2, y2, ..., f2}, ...}, where the
number of coordinates x, y, ... is equal to the number of variables in
the list vars.  The data can also be of the form {f1, f2, ...}, with a
single coordinate assumed to take values 1, 2, .... The argument funs can
be any list of functions that depend only on the objects vars.  The result
is a list of rules identifying summary statistics about the fit, or other
statistics as specified by the option RegressionReport."

Options[Regress] = {RegressionReport -> SummaryReport,
 IncludeConstant -> True, BasisNames -> Automatic, Weights -> Automatic} ~Join~
   Options[SingularValues] ~Join~ Options[StudentTCI]

DesignedRegress::usage =
"DesignedRegress[designmatrix, response] finds a least-squares fit given a 
design matrix and a vector of the response data.  The result is a list of
rules identifying summary statistics about the fit, or other statistics as
specified by the option RegressionReport."

Options[DesignedRegress] = {RegressionReport -> SummaryReport,
   BasisNames -> Automatic, Weights -> Automatic} ~Join~
   Options[SingularValues] ~Join~ Options[StudentTCI]

RegressionReportValues[DesignedRegress] =
{AdjustedRSquared, ANOVATable, BestFitParameters, 
BestFitParametersDelta, CatcherMatrix, CoefficientOfVariation,
CookD, CorrelationMatrix, CovarianceMatrix,
CovarianceMatrixDetRatio, DurbinWatsonD, EigenstructureTable,
EstimatedVariance, FitResiduals, HatDiagonal, JackknifedVariance,
MeanPredictionCITable, ParameterCITable, ParameterConfidenceRegion,
ParameterTable, PartialSumOfSquares, PredictedResponse, PredictedResponseDelta,
RSquared, SequentialSumOfSquares, SinglePredictionCITable,
StandardizedResiduals, StudentizedResiduals, SummaryReport, VarianceInflation}

RegressionReportValues[Regress] =
	Sort[ RegressionReportValues[DesignedRegress] ~Join~ {BestFit} ]

DesignMatrix::usage = 
"DesignMatrix[data, funs, vars] gives the design matrix for modeling the
data as a linear combination of the functions. The data can have the 
form {{x1, y1, ..., f1}, {x2, y2, ..., f2}, ...}, where the
number of coordinates x, y, ... is equal to the number of variables in
the list vars.  The data can also be of the form {f1, f2, ...}, with a
single coordinate assumed to take values 1, 2, .... The argument funs can
be any list of functions that depend only on the objects vars."

Options[DesignMatrix] = {IncludeConstant -> True}


	(* options unique to Regress and DesignedRegress *)

IncludeConstant::usage =
"IncludeConstant is an option to Regress and DesignMatrix and specifies whether
a constant term is automatically included in the model.  If the set of basis
functions does not include a constant term and IncludeConstant->True, then a
constant term is added.  If the set of basis functions includes a constant term
and IncludeConstant -> False, then the set of basis functions is left
unchanged."


	(* support for obsolete options and option values *)

OutputList::usage =
"OutputList is an obsolete symbol, superceded by RegressionReport.
The obsolete options OutputList -> Null and OutputControl -> Automatic
can be achieved by RegressionReport -> SummaryReport."

OutputControl::usage = 
"OutputControl is an obsolete option to regression functions.  The
obsolete option OutputControl -> Automatic can be achieved by including
SummaryReport in the value of the RegressionReport option.  The obsolete
option OutputControl -> NoPrint can be achieved by omitting SummaryReport
from the value of the RegressionReport option."

NoPrint::usage =
"NoPrint is an obsolete value of the obsolete option OutputControl.
The obsolete option OutputControl -> NoPrint can be achieved by omitting
SummaryReport from the value of the RegressionReport option." 

BestFitCoefficients::usage =
"BestFitCoefficients is an obsolete symbol, superceded by BestFitParameters."

ConfidenceIntervalTable::usage =
"ConfidenceIntervalTable is an obsolete symbol, superceded by ParameterCITable."


	(* updating usage of shared symbols *)

If[StringPosition[CovarianceMatrix::usage, "regression"] === {},

(* CovarianceMatrix exists, but regression usage does not. *)
CovarianceMatrix::usage = CovarianceMatrix::usage <> " " <>
"CovarianceMatrix is used in the output of regression functions
to identify the estimated covariance matrix of the fit parameters.";
CorrelationMatrix::usage = CorrelationMatrix::usage <> " " <>
"CorrelationMatrix is used in the output of regression functions
to identify the estimated correlation matrix of the fit parameters.",

(* CovarianceMatrix does not exist. *)
CovarianceMatrix::usage =
"CovarianceMatrix is used in the output of regression functions
to identify the estimated covariance matrix of the fit parameters.";
CorrelationMatrix::usage =
"CorrelationMatrix is used in the output of regression functions
to identify the estimated correlation matrix of the fit parameters."

]

If[StringPosition[SummaryReport::usage, "For linear regression"] === {},

(* SummaryReport exists, but linear regression usage does not. *)
SummaryReport::usage = SummaryReport::usage <> " " <>
"For linear regression functions, SummaryReport specifies a summary including
ParameterTable, RSquared, AdjustedRSquared, EstimatedVariance, and ANOVATable.",

(* SummaryReport does not exist. *)
SummaryReport::usage =
"SummaryReport is a value of the option RegressionReport specifying that
linear regression functions should return a summary including ParameterTable,
RSquared, AdjustedRSquared, EstimatedVariance, and ANOVATable."

]

	(* updating usage of DescriptiveStatistics.m symbols *)

If[StringPosition[CoefficientOfVariation::usage, 
	 "linear regression"] === {},

(* CoefficientOfVariation exists, but linear regression usage does not. *)
CoefficientOfVariation::usage = CoefficientOfVariation::usage <> " " <>
"CoefficientOfVariation is used in the output of linear regression functions
to identify the ratio of the estimated error standard deviation (residual
root mean square) to the mean of the response variable."

]


Begin["`Private`"]

(* =========================== obsolete symbols =========================== *)

BestFitCoefficients = BestFitParameters
ConfidenceIntervalTable = ParameterCITable


(* ============================== Regress ================================= *)

Regress[data_, basis_, vars_, options___Rule] :=
  Block[{answer, varlist, varcount, coordcount},
    answer /; (varlist = If[VectorQ[vars], vars, {vars}];
	       varcount = Length[varlist];
	       answer = Which[
		  Apply[Or, Map[NumberQ, N[varlist]]],
		     Message[Regress::crdvl, varlist];
		     $Failed,	
 	          MatrixQ[data], (* data is a matrix *)
		     If[varcount <= (coordcount = Length[data[[1]]]-1),
		        vRegress[data, basis, varlist, options],
		        Message[Regress::coord, coordcount, varcount];
		        $Failed],
		  True, (* data is a vector *)
		     If[varcount === 1,
		        vRegress[Transpose[{Range[Length[data]],data}],
				basis, varlist, options],
		        Message[Regress::coord, 1, varcount];
		        $Failed]
	       ];
	       answer =!= $Failed) 	
  ] /; If[!(VectorQ[data] || MatrixQ[data]),
	    Message[Regress::notdata]; False, True] &&
       If[!VectorQ[basis], 
	    Message[Regress::notfcn]; False, True] &&
       Head[vars] =!= Rule

vRegress[data_, basis_, varlist_, options___] :=
    Module[{response = Map[Last, data], cbasis = basis,
	    tempdata, tempweights,  
	    designMatrix, optionlist = {options},
	    dmPrecision, nResponse, nDesignMatrix,
	    report, includeC, names, weights, tol, cLevel, outlist, control},
	(* establish the values of these options according to the defaults
		of Regress *)
	{includeC, names, weights, tol, cLevel} =
		{IncludeConstant, BasisNames, Weights,
		 Tolerance, ConfidenceLevel} /. optionlist /. Options[Regress];
	If[!(Or @@ Apply[And, Outer[FreeQ, basis, varlist], {1}]) && includeC,
            PrependTo[cbasis, 1]
        ];

	(* check and set Weights *)
        (* This must be done before calling vDesignedRegress, because
		if Weights -> pureFcn, then the Weights setting must be
		converted to Weight -> vector while the data is still 
		available.  vDesignedRegress can't do this conversion,
		because it has access to the design matrix,
		not the raw data. *)
	Which[
         weights === Equal || weights === Automatic,
                weights = Table[1, {Length[response]}],
         Head[weights] === Function,
		(* Rearrange data so that response is first... if weights
			has only one arg, then it is applied to the response.
		   Note that as of May 1998, the ability to use a weights
			function of more than one arg is not documented. *)
		tempdata = Map[Join[{Last[#]}, Drop[#, -1]]&, data];
                tempweights = Map[Apply[weights, #]&, tempdata];
                If[!VectorQ[N[tempweights], NumberQ],
                        Message[Regress::bdwghts, weights];
                        weights = Table[1, {Length[response]}],
                        weights = tempweights
                ],
         VectorQ[weights],
                If[!VectorQ[N[weights], NumberQ] ||
                   Length[weights] =!= Length[response],
                        Message[Regress::bdwghts, weights];
                        weights = Table[1, {Length[response]}]
                ],
         True,
                Message[Regress::bdwghts, weights];
                weights = Table[1, {Length[response]}]
        ];

	(* calculate designMatrix *)
	If[!FreeQ[designMatrix =
	   DesignMatrix[data, cbasis, varlist, IncludeConstant -> includeC],
	   DesignMatrix],
	   Message[Regress::nodm];
	   Return[$Failed]];
	If[Length[response] <= Length[cbasis],
            Message[Regress::mindata]
        ];
	If[names === Automatic, 
		names = cbasis,
		If[!VectorQ[names],
		   Message[Regress::invnam]; names = cbasis,
		   If[Length[names] =!= Length[cbasis],
		      Message[Regress::namnum]; names = cbasis
		   ]
		] 
	];


	dmPrecision = Precision[designMatrix];
	(* make nDesignMatrix suitable for SingularValues *)
	{nDesignMatrix, nResponse} = If[dmPrecision === Infinity,
                      N[{designMatrix, response}],
			{designMatrix, response} ];

	(* make sure that the default options of Regress are used instead of
		the default options of DesignedRegress *)
	If[MemberQ[optionlist, OutputControl -> _] ||
	   MemberQ[optionlist, OutputList -> _],
	   (* invoke vDesignedRegress with obsolete options *)
	   {control, outlist} = {OutputControl, OutputList} /. optionlist /.
		{OutputList -> Null, OutputControl -> Automatic};
	   vDesignedRegress[nDesignMatrix, nResponse, cbasis,
                OutputControl -> control, OutputList -> outlist,
		IncludeConstant -> includeC,
                BasisNames -> names, Weights -> weights, Tolerance -> tol,
                ConfidenceLevel -> cLevel],
	   (* invoke vDesignedRegress with new options *)
	   report = RegressionReport /. optionlist /. Options[Regress];
           vDesignedRegress[nDesignMatrix, nResponse, cbasis,
		RegressionReport -> report, IncludeConstant -> includeC,
		BasisNames -> names, Weights -> weights, Tolerance -> tol,
		ConfidenceLevel -> cLevel]
	]

    ] (* end vRegress *)


(* =========================== DesignMatrix ============================== *)

DesignMatrix[data_, basis_, vars_, opts___Rule] :=
    Block[{answer, varlist, varcount, coordcount},
        answer /; (varlist = If[VectorQ[vars], vars, {vars}];
	           varcount = Length[varlist];
               answer = Which[
		  Apply[Or, Map[NumberQ, N[varlist]]],
                     Message[DesignMatrix::crdvl, varlist];
                     $Failed,
                  MatrixQ[data], (* data is a matrix *)
                     If[varcount <= (coordcount = Length[data[[1]]]-1),
                        vDesignMatrix[data, basis, varlist, opts],
                        Message[DesignMatrix::coord, coordcount, varcount];
                        $Failed],
                  True, (* data is a vector *)
                     If[varcount === 1,
                        vDesignMatrix[Transpose[{Range[Length[data]],data}],
                                basis, varlist, opts],
                        Message[DesignMatrix::coord, 1, varcount];
                        $Failed]
               ];
	       answer =!= $Failed)
    ] /; If[!(VectorQ[data] || MatrixQ[data]),
	    Message[DesignMatrix::notdata]; False, True] &&
	 If[!VectorQ[basis],
            Message[DesignMatrix::notfcn]; False, True]	 &&
	 Head[vars] =!= Rule

vDesignMatrix[data_, basis_, varlist_, opts___] :=
    Module[{cpresent = Or @@ Apply[And, Outer[FreeQ, basis, varlist], {1}],
	    cbasis = basis, pureFunction, newbasis, optionlist = {opts},
        x, newvars, vartrans},
	If[optionlist =!= {} && !MemberQ[optionlist, IncludeConstant->_],
	   Message[DesignMatrix::invopt, optionlist]];
	(* If there is not already a constant in the basis functions, then
		add one if IncludeConstant is True. *)
	If[!cpresent &&
	   (IncludeConstant /. optionlist /. Options[DesignMatrix]),
	   PrependTo[cbasis, 1]
        ];
    newvars = Table[Unique[x], {Length[varlist]}];
    vartrans = Thread[varlist -> newvars];
	If[!FreeQ[cbasis,Function],
	   (* pure fcns in the basis list are supported, but not documented *)
           Message[DesignMatrix::pure];
           newbasis = Map[Module[{fcn = #},
			    If[FreeQ[fcn, Function],
			       fcn,
			       fcn @@ varlist]]&, cbasis];
	   pureFunction := Function[Evaluate[newvars], Evaluate[newbasis/.vartrans]],
	   (* normally the basis list does not include pure fcns *)
	   pureFunction := Function[Evaluate[newvars], Evaluate[cbasis/.vartrans]]
        ];
	Apply[pureFunction, data, 1]
    ]


(* ================= DesignedRegress[designMatrix, ...] =================== *)

DesignedRegress[designMatrix_?MatrixQ, response_?VectorQ, options___Rule] :=
  Block[{answer, dmPrecision, nDesignMatrix, nResponse, optionlist,
	 report, pos}, (* ttt does not belong here *)
    (
    answer
    ) /; (dmPrecision = Precision[{designMatrix, response}];
	  (* make nDesignMatrix suitable for SingularValues *)
	  {nDesignMatrix, nResponse} = If[dmPrecision === Infinity,
          	      N[{designMatrix, response}],
          	      {designMatrix, response} ];
	  optionlist = {options};
	  If[MemberQ[optionlist, IncludeConstant->_],
	     (* IncludeConstant should be an option of DesignMatrix,
			not DesignedRegress. *)
	     Message[DesignedRegress::invopt, IncludeConstant];
	     optionlist = Delete[optionlist, Position[optionlist,
		IncludeConstant->_]]	
	  ];
	  (* filter out invalid report items *)
          report = RegressionReport /. optionlist /. Options[DesignedRegress];
	  If[ListQ[report],
	     If[MemberQ[report, BestFit],
		(* we can't return a symbolic expression when the input is
			entirely numeric *)
		Message[DesignedRegress::invrep, BestFit];
		report = Delete[report, Position[report, BestFit]]
  	     ];
	     pos = Position[optionlist, RegressionReport -> _];
	     optionlist = Delete[optionlist, pos];	
	     optionlist = Insert[optionlist, RegressionReport->report, pos]
	  ]; 
	  (* ttt array is just a place holder *)
	  answer = vDesignedRegress[nDesignMatrix, nResponse,
		 Array[ttt, Length[nDesignMatrix[[1]]]],
		 Apply[Sequence, optionlist] ];
	  answer =!= $Failed)	
  ] /; If[MatrixQ[N[designMatrix], NumberQ],
	  True,
	  Message[DesignedRegress::invdm]; False] &&
       If[Length[response] == Length[designMatrix],
	  True,
	  Message[DesignedRegress::invresp, Length[response],
		 Length[designMatrix]];
	  False] (* end of DesignedRegress[designMatrix, etc.] *)


vDesignedRegress[designMatrix_?MatrixQ, response_,
	 basis_, (* this is either a valid list of basis functions
			or a place holder {ttt[1], ttt[2], ...} *)
	 options___] :=
    Module[{temp, constantQ, n, p, dmPrecision,
	    optionlist = {options}, tol, weights,
	    tempdata, tempweights,	
	    sqrtweights, weightedDesignMatrix, svd, u, w, v},

	If[FreeQ[basis, ttt],
           (* basis list is not just a place-holder *)	
	   With[{temp = Apply[Plus, Map[If[NumberQ[#], 1, 0]&, basis]]},
		Which[temp == 1,
		      constantQ = True,
		      temp == 0,
		      constantQ = False,
		      temp > 1,
		      Message[DesignedRegress::twocon];
		      Return[$Failed],
		      True,
		      Message[DesignedRegress::uncon];
		      Return[$Failed]
		]	
	   ],
	   (* It is unknown whether constantQ = True or False, that is
	      whether there is a single constant basis fcn represented in the
	      design matrix or not. *)
	   constantQ = UNKNOWN];

	{n, p} = Dimensions[designMatrix];

        {tol, weights} = {Tolerance, Weights} /.
		 optionlist /. Options[DesignedRegress];
	(* set and check weights *)
	Which[
         weights === Equal || weights === Automatic,
                weights = Table[1, {Length[response]}],
         Head[weights] === Function,
		(* Rearrange data so that response is first... if weights
                        has only one arg, then it is applied to the response.
                   Note that as of May 1998, the ability to use a weights
                        function of more than one arg is not documented. *)
		(* vDesignedRegress has no access to the independent
		   variables, so a pure function for computing
		   weights must be applied to the response vector. *)
		tempdata = Map[Prepend[designMatrix[[#]], response[[#]]]&,
                                Range[n]];
                tempweights = Map[Apply[weights, #]&, tempdata];
                If[!VectorQ[N[tempweights], NumberQ],
                        Message[DesignedRegress::bdwghts, weights];
                        weights = Table[1, {Length[response]}],
                        weights = tempweights
                ],
         VectorQ[weights],
                If[!VectorQ[N[weights], NumberQ] ||
                   Length[weights] =!= Length[response],
                        Message[DesignedRegress::bdwghts, weights];
                        weights = Table[1, {Length[response]}]
                ],
         True,
                Message[DesignedRegress::bdwghts, weights];
                weights = Table[1, {Length[response]}]
        ];
	dmPrecision = Precision[designMatrix];
	If[Precision[weights] === Infinity,
		 weights = N[weights, Max[dmPrecision, $MachinePrecision]]
	];
	sqrtweights = Sqrt[weights];


	weightedDesignMatrix = sqrtweights designMatrix;

        svd = SingularValues[weightedDesignMatrix, Tolerance -> tol];
        If[Head[svd] === List && Length[svd] === 3,
            {u, w, v} = svd,
            Message[DesignedRegress::nosvd];
            Return[$Failed]
        ];
        If[Length[w] < p, Message[DesignedRegress::rank,
		 Length[w], p, If[Length[w]==1, "is", "are"] ] ]; 
	vDR[weightedDesignMatrix, {u, w, v}, sqrtweights response, basis, 
		sqrtweights, constantQ, options]
		
    ] (* end of vDesignedRegress[designMatrix, etc] *)


(* =================== DesignedRegress[{u, w, v}, ...] ===================== *)

DesignedRegress[{u_?MatrixQ, w_?VectorQ, v_?MatrixQ}, response_?VectorQ,
	 options___Rule] :=
  Block[{answer, optionlist, report, pos}, (* ttt does not belong here *)
    (
    answer
    ) /; (optionlist = {options};
	  If[MemberQ[optionlist, IncludeConstant->_],
	     (* IncludeConstant should be an option of DesignMatrix,
			not DesignedRegress. *)
	     Message[DesignedRegress::invopt, IncludeConstant];
	     optionlist = Delete[optionlist, Position[optionlist,
		IncludeConstant->_]]	
	  ];
	  (* filter out invalid report items *)
          report = RegressionReport /. optionlist /. Options[DesignedRegress];
	  If[ListQ[report],
	     If[MemberQ[report, BestFit],
		(* we can't return a symbolic expression when the input is
			entirely numeric *)
		Message[DesignedRegress::invrep, BestFit];
		report = Delete[report, Position[report, BestFit]]
  	     ];
	     pos = Position[optionlist, RegressionReport -> _];
	     optionlist = Delete[optionlist, pos];	
	     optionlist = Insert[optionlist, RegressionReport->report, pos]
	  ]; 
	  (* ttt matrix array is just a place holder *)
	  answer = vDesignedRegress[{u, w, v}, response, Array[ttt, Length[w]],
		 Apply[Sequence, optionlist] ];
	  answer =!= $Failed)	
  ] /; If[MatrixQ[u, NumberQ] && VectorQ[w, NumberQ] && MatrixQ[v, NumberQ],
	  True,
	  Message[DesignedRegress::invsvd]; False] &&
       If[Length[response] == Length[u[[1]]],
	  True,
	  Message[DesignedRegress::invresp, Length[response], Length[u[[1]]]];
	  False] (* end of DesignedRegress[{u, w, v}, etc.] *)
	

vDesignedRegress[{u_?MatrixQ, w_?VectorQ, v_?MatrixQ}, response_,
	 basis_, (* this is a place holder {ttt[1], ttt[2], ...} *)
	 options___] :=
    Module[{n, p, optionlist = {options}, tol, weights, tempdata, tempweights,
		wPrecision,  sqrtweights,
	  	designMatrix = Transpose[u].DiagonalMatrix[w].v,
	  	weightedDesignMatrix, svd, uprime, wprime, vprime},
	{n, p} = Dimensions[designMatrix];

        {tol, weights} = {Tolerance, Weights} /.
		 optionlist /. Options[DesignedRegress];
	(* set and check weights *)
        Which[
         weights === Equal || weights === Automatic,
                weights = Table[1, {Length[response]}],
         Head[weights] === Function,
                (* Rearrange data so that response is first... if weights
                        has only one arg, then it is applied to the response.
                   Note that as of May 1998, the ability to use a weights
                        function of more than one arg is not documented. *)
                (* vDesignedRegress has no access to the independent
                   variables, so a pure function for computing
                   weights must be applied to the response vector. *)
                tempdata = Map[Prepend[designMatrix[[#]], response[[#]]]&,
                                Range[n]];
                tempweights = Map[Apply[weights, #]&, tempdata];
                If[!VectorQ[N[tempweights], NumberQ],
                        Message[DesignedRegress::bdwghts, weights];
                        weights = Table[1, {Length[response]}],
                        weights = tempweights
                ],
         VectorQ[weights],
                If[!VectorQ[N[weights], NumberQ] ||
                   Length[weights] =!= Length[response],
                        Message[DesignedRegress::bdwghts, weights];
                        weights = Table[1, {Length[response]}]
                ],
         True,
                Message[DesignedRegress::bdwghts, weights];
                weights = Table[1, {Length[response]}]
        ];
	wPrecision = Precision[w];
	If[Precision[weights] === Infinity,
		 weights = N[weights, Max[wPrecision, $MachinePrecision]]
	];
	sqrtweights = Sqrt[weights];


	weightedDesignMatrix = sqrtweights designMatrix;

	(* Here it is unknown whether constantQ = True or False, that is
	   whether there is a single constant basis fcn represented in the
	   design matrix or not.  The problem with this is that vDR
	   will check the designMatrix and assign a value to constantQ,
	   and this can be time-consuming.  Perhaps all the time saved
	   in reusing the singular value decomposition of the design matrix
	   is lost in making this sanity check. *)
	constantQ = UNKNOWN;
	If[TrueQ[Apply[Equal, sqrtweights]],
	   (* use existing svd *)
	   vDR[weightedDesignMatrix, {u, (sqrtweights[[1]])^2 w, v},
		sqrtweights response, basis, sqrtweights,
		constantQ, options],
	   (* compute new svd *)
	   Message[DesignedRegress::newsvd];
           svd = SingularValues[weightedDesignMatrix, Tolerance -> tol];
           If[Head[svd] === List && Length[svd] === 3,
               {uprime, wprime, vprime} = svd,
               Message[DesignedRegress::nosvd];
               Return[$Failed]
           ];
           If[Length[wprime] < p,
		 Message[DesignedRegress::rank, Length[wprime], p]]; 
	   vDR[weightedDesignMatrix, {uprime, wprime, vprime},
		sqrtweights response, basis, sqrtweights,
		constantQ, options]
	]

    ] (* end of vDesignedRegress[{u, w, v}, etc] *)


(* ================================ vDR ================================= *)

vDR[weightedDesignMatrix_, {u_, w_, v_}, weightedResponse_, basis_,
	 sqrtweights_, constantQ_, options___] :=
    (* weightedDesignMatrix:   weighted design matrix for regression         *)
    (* {u, w, v}:      singular value decomposition of weightedDesignMatrix *)
    (* weightedResponse:       vector of response values                     *)
    (* basis:								     *)
    (*   symbolic basis functions if vDesignedRegress is called by vRegress,  *)
    (*   place-holding array if vDesignedRegress is called by DesignedRegress *)
    (* optionlist:     list of options                                    *)
    Module[{n, p, optionlist = {options}, constantPresent = constantQ,
	  meanWeightedResponse,
          modelDOF, totalDOF, errorDOF, modelSS, totalSS, errorSS,
          pValue, fRatio,
	  weights = sqrtweights^2,
	  designMatrix = weightedDesignMatrix/sqrtweights,
	  response = weightedResponse/sqrtweights,
          tol,			(* tolerance for SingularValues *)
	  names, 		(* basis fcn names for table headings *)
	  catcherMatrix, 	(* terminology of Mosteller & Tukey *)
          fit, 			(* best fit parameters *)
          responseMean, predictedResponse, fitResiduals,
	  rSquared, adjRSquared, coeffvar,
          totalResiduals, 	(* residuals from fitting a constant *)
          responseVariance, 	(* estimated variance from best fit *)
          report, 		(* output accumulation *)
          designMatrixT, headings, cLevel, singlepredci, meanpredci, varinf,
	  sosDiagnostics	(* flag indicating whether design matrix is
				   suitable for computing PartialSumOfSquares
				   or SequentialSumOfSquares *)
	 },

	{n, p} = Dimensions[designMatrix];


     (* READ IN OPTIONS *)

       	optionlist = {options};
        {tol, names, cLevel} =
		 {Tolerance, BasisNames, ConfidenceLevel} /.
		 optionlist /. Options[DesignedRegress];
 

    (*  COMPUTE PARAMETERS *) 

	(* If X = design matrix and Y = response, then
		X = Transpose[u].DiagonalMatrix[w].v      
		fit = Inverse[Transpose[X].X].Transpose[X].Y
		    = Transpose[v].DiagonalMatrix[1/w].u.Y
		    = catcherMatrix.Y = Transpose[v].((u.Y)/w)		*)

	(* NOTE: The following manner of calculating catcherMatrix unnecessarily
		constructs DiagonalMatrix[sqrtweights] (an n x n matrix) and is
		slow.
	   catcherMatrix =
	     Transpose[v].DiagonalMatrix[1/w].u.DiagonalMatrix[sqrtweights]; *)

	catcherMatrix = Map[sqrtweights # &, u];
	catcherMatrix = Transpose[v].DiagonalMatrix[1/w].catcherMatrix;
        fit = catcherMatrix . response;

	(* NOTE: Calculate hatMatrix later when you know that the HatDiagonal
		diagnostic is needed.    hatMatrix is also an n x n matrix.
		 hatMatrix = designMatrix . catcherMatrix; *)

        predictedResponse = designMatrix . fit; (* hatMatrix . response *)
        fitResiduals = response - predictedResponse;
        errorDOF = n - p;
        errorSS = weights . (fitResiduals^2);

	designMatrixT = Transpose[designMatrix];

     (*  COMPARING MODELS  *)

	If[!(constantPresent === True || constantPresent === False),
	   (* Examine design matrix to determine whether the model included
		a single constant basis function.  This is time-consuming
		for a large sample. *)
	   With[{temp = Apply[Plus,
			 Map[If[Apply[Equal, #], 1, 0]&, designMatrixT]]},
		Which[temp == 1,
		      constantPresent = True,
		      temp == 0,
		      constantPresent = False,
		      temp > 1,
		      Message[DesignedRegress::twocon];
		      Return[$Failed],
		      True,
		      Message[DesignedRegress::uncon];
		      Return[$Failed]
		]	
	   ]
	];
	
        If[constantPresent,
	   meanWeightedResponse =
	      sqrtweights (weightedResponse.sqrtweights)/Apply[Plus,weights];
	   totalResiduals = weightedResponse - meanWeightedResponse;
           totalDOF = n - 1,
           (* else if !constantPresent *)
           totalResiduals = weightedResponse;
           totalDOF = n
        ];
        totalSS = totalResiduals . totalResiduals;
        modelDOF = totalDOF - errorDOF;
        modelSS = totalSS - errorSS;
	If[modelSS < 0,
	   Message[DesignedRegress::badfit];
	   rSquared = $Failed;  adjRSquared = $Failed,
	   rSquared = modelSS/totalSS;
	   adjRSquared = 1-(1-rSquared)(totalDOF/errorDOF)];

        responseVariance = errorSS/errorDOF;


     (*  REPORT GENERATION *)

        report = GetOutputList[optionlist];


	(* ================= compute ANOVA table ================ *)

 	If[MemberQ[report, ANOVATable],
	   If[modelSS >= 0,
	      If[!constantPresent,
	         (* Use message to remind user of the changed interpretation 
		    of total sum of squares when there is no constant in
		    the model. *)
		 Message[DesignedRegress::tsos];
	         (* Use "U Total" sum of squares heading to remind user of 
		    the changed interpretation when there is no constant term
		    in the model. *)
		 headings = {{"Model", "Error", "U Total"},
			{"DF", "SumOfSq", "MeanSq", "FRatio", "PValue"}},
		 headings = {{"Model", "Error", "Total"},
			{"DF", "SumOfSq", "MeanSq", "FRatio", "PValue"}}
	      ];
              fRatio = (modelSS/modelDOF)/responseVariance;
              pValue = 1 - CDF[FRatioDistribution[modelDOF, errorDOF], fRatio];
              report = report /. ANOVATable -> (ANOVATable -> TableForm[
              	{{modelDOF, modelSS, modelSS/modelDOF, fRatio, pValue},
               	{errorDOF, errorSS, errorSS/errorDOF},
               	{totalDOF, totalSS}}, TableHeadings->headings ]),
	      (* modelSS < 0 *)
	      If[!constantPresent,
	         Message[DesignedRegress::tsos];
		 headings = {{"Error", "U Total"}, {"DF", "SumOfSq", "MeanSq"}},
		 headings = {{"Error", "Total"}, {"DF", "SumOfSq", "MeanSq"}}
	      ]; 
	      report = report /. ANOVATable -> (ANOVATable -> TableForm[
		{{errorDOF, errorSS, errorSS/errorDOF},
		 {totalDOF, totalSS}}, TableHeadings->headings ])	
	   ] (* end If modelSS >= 0 *)
        ]; 

        (* ========= compute things involving the covariance matrix ======== *)

        If[Or @@ Map[MemberQ[report, #]&, {CovarianceMatrix,
             CorrelationMatrix, MeanPredictionCITable,
	     ParameterTable, ParameterCITable,
	     ParameterConfidenceRegion, ParameterConfidenceRegion[_?VectorQ],
	     PartialSumOfSquares, SinglePredictionCITable}],
           Block[{unscaledcovMatrix, covMatrix, corMatrix, parmci, predci,
		   infy, indet,
                   parmse, (* list of standard errors of parm estimates *)
                   predse, (* list of standard errors of predicted response *)
		   tRatio, tPValues},
	    (* expression for covariance matrix does not depend on whether the
		regression is weighted or not, except that v and w are from
		the svd of a weighted design matrix *)
	    unscaledcovMatrix = Transpose[v] . (v / w^2);
	    covMatrix = responseVariance unscaledcovMatrix;
	    report = report /. CovarianceMatrix ->
                   (CovarianceMatrix -> MatrixForm[covMatrix]);
            parmse = Sqrt[DiagonalElements[covMatrix]];
	    If[MemberQ[report, CorrelationMatrix],
               infy = (Head[Power::infy] === $Off);  Off[Power::infy];	
	       indet = (Head[Infinity::indet] === $Off);  Off[Infinity::indet]; 
               corMatrix = Transpose[Transpose[covMatrix/parmse] / parmse];
	       If[!infy, On[Power::infy]];  If[!indet, On[Infinity::indet]];
 	       report = report /. CorrelationMatrix ->
                   (CorrelationMatrix -> MatrixForm[corMatrix]);
	    ]; 
	    If[MemberQ[report, ParameterTable],
	        infy = (Head[Power::infy] === $Off);  Off[Power::infy];
	        indet = (Head[Infinity::indet] === $Off);  Off[Infinity::indet];
            	tRatio = fit / parmse; 
	        If[!infy, On[Power::infy]];  If[!indet, On[Infinity::indet]];
		tPValues = Map[If[# === Indeterminate, Indeterminate,
			CDF[StudentTDistribution[errorDOF],#]]&, tRatio];
		tPValues = Map[If[# === Indeterminate, Indeterminate,
				  2 If[#>.5, 1-#, #]]&, tPValues];
	        (* Note: when names is not defined, the default row headings
			are used (i.e., 1, 2, 3, ....) *)
	    	report = report /. ParameterTable ->
             		(ParameterTable -> TableForm[
			   Transpose[{fit, parmse, tRatio, tPValues}],
              		TableHeadings ->
			 {names, {"Estimate", "SE", "TStat", "PValue"}}])
	    ];
	    If[MemberQ[report, PartialSumOfSquares],
 	       sosDiagnostics = p > 1 &&
			TrueQ[Apply[Equal, designMatrixT[[1]]]] &&
			Apply[And, Map[!TrueQ[Apply[Equal, #]]&,
			 	Drop[designMatrixT, 1] ]];
	       If[sosDiagnostics,
		        partialsos =
			 Drop[fit^2/DiagonalElements[unscaledcovMatrix], 1];
			report = report /. {PartialSumOfSquares ->
			  (PartialSumOfSquares -> partialsos)},
			Message[DesignedRegress::parsos];
			report = report /. {PartialSumOfSquares ->
                                (PartialSumOfSquares -> $Failed)}	]
  	    ];
	    If[MemberQ[report, ParameterCITable],
            	parmci = Apply[N[StudentTCI[#1, #2, errorDOF,
                                   ConfidenceLevel -> cLevel]]&,
                           Transpose[{fit, parmse}], 1];        
                report = report /. ParameterCITable ->
             		(ParameterCITable -> TableForm[
				Transpose[{fit, parmse, parmci}],
			 TableDepth -> 2,
	      		TableHeadings -> {names, {"Estimate", "SE", "CI"}}])
	    ];
	    If[MemberQ[report, MeanPredictionCITable] ||
	       MemberQ[report, SinglePredictionCITable],
	       	predse = Sqrt[Map[(#.covMatrix.#)&, designMatrix]] ];	
            If[MemberQ[report, MeanPredictionCITable],
	        predci = Apply[N[StudentTCI[#1, #2, n-p,
				ConfidenceLevel -> cLevel]]&,
				Transpose[{predictedResponse, predse}], 1];
	        report = report /. MeanPredictionCITable ->
	     	   (MeanPredictionCITable -> TableForm[
		    Transpose[{response, predictedResponse, predse, predci}],
			TableDepth -> 2, TableHeadings ->
		 	{None, {"Observed", "Predicted", "SE", "CI"}}])
            ];
            If[MemberQ[report, SinglePredictionCITable],
		predse = Sqrt[responseVariance + predse^2];
	        predci = Apply[N[StudentTCI[#1, #2, n-p,
				ConfidenceLevel -> cLevel]]&,
				Transpose[{predictedResponse, predse}], 1];
	        report = report /. SinglePredictionCITable ->
	     	   (SinglePredictionCITable -> TableForm[
		    Transpose[{response, predictedResponse, predse, predci}],
			TableDepth -> 2, TableHeadings ->
		 	{None, {"Observed", "Predicted", "SE", "CI"}}])
            ];
	    If[MemberQ[report, ParameterConfidenceRegion] ||
	       MemberQ[report, ParameterConfidenceRegion[_?VectorQ]],
               (* Weisberg, p. 97 *)
	       Block[{ellipsoid, arglist, pos, 
		      basisSubset, fitSubset, covMatrixSubset, pSubset},
	          If[MemberQ[report, ParameterConfidenceRegion],
		     ellipsoid = ellipsoidalLocus[fit, covMatrix * p *
			Quantile[FRatioDistribution[p, n-p], cLevel]];
	             report = report /. ParameterConfidenceRegion ->
                  	(ParameterConfidenceRegion -> ellipsoid)
	          ];
	          If[MemberQ[report, ParameterConfidenceRegion[_?VectorQ]],
 		     arglist = Map[First, Select[report,
			 MatchQ[#, ParameterConfidenceRegion[l_List] /;
			 	 Apply[And, Map[MemberQ[basis, #]&, l]]]& ]];
	             Scan[(basisSubset = #;
		           pos = Flatten[Map[Position[basis, #, 1]&,
				 basisSubset]];
		   	   fitSubset = fit[[pos]];
		           covMatrixSubset =
				 Transpose[Transpose[covMatrix[[pos]]][[pos]]];
		   	   pSubset = Length[basisSubset];
		   	   ellipsoid = ellipsoidalLocus[fitSubset,
			     covMatrixSubset * pSubset * Quantile[
			     FRatioDistribution[pSubset, n-p], cLevel]];
		   	   report = report /.
				 (ParameterConfidenceRegion[basisSubset] ->
		     	         (ParameterConfidenceRegion[basisSubset] ->
				 ellipsoid))
		          )&, arglist]
	          ]  (* end If *)
	       ] (* end Block *)
	    ] (* end If *)


	   ]  (* end Block *)
        ];	(* end of covariance matrix calculations *)


        (* =========== compute things involving the hat matrix ========== *)

        If[Or @@ Map[MemberQ[report, #]&, {HatDiagonal,
		StandardizedResiduals, StudentizedResiduals,
		PredictedResponseDelta, CookD, CovarianceMatrixDetRatio,
		BestFitParametersDelta, JackknifedVariance}],
	   Block[{hatMatrix, (* yhat = hatMatrix . y *)
		  hatdiag, standResiduals, cookD, jackvar, studResiduals,
	          dfit, covratio, parmInfo, obsrInfo, dbeta},
	     hatMatrix = designMatrix . catcherMatrix; 
	     hatdiag = DiagonalElements[hatMatrix];
	     (* 8/98: changed standResiduals from
		fitResiduals/Sqrt[responseVariance (1-hatdiag)]
		to
		fitResiduals/Sqrt[responseVariance (1-hatdiag)/weights]
		*)
	     standResiduals = fitResiduals/
		Sqrt[responseVariance (1-hatdiag)/weights];
	     cookD = hatdiag/(1 - hatdiag) standResiduals^2 / p;
	     (* Exhibit 2B.1, p. 75, Belsley, Kuh, Welsch *)
	     jackvar = ((n-p) responseVariance -
		fitResiduals^2/(1 - hatdiag)) / (n-p-1);
	     (* 8/98: changed studResiduals from
	        studResiduals = fitResiduals/Sqrt[jackvar (1-hatdiag)];
		to
	        studResiduals = fitResiduals/Sqrt[jackvar (1-hatdiag)/weights];
		*)
	     studResiduals = fitResiduals/Sqrt[jackvar (1-hatdiag)/weights];
	     dfit = Sqrt[hatdiag/(1 - hatdiag)] studResiduals;
	     covratio = (n-p)^p/( (n-p-1+studResiduals^2)^p (1-hatdiag) ); 	
	     If[MemberQ[report, BestFitParametersDelta],
		parmInfo = 1/Sqrt[Map[(Apply[Plus, #^2])&, catcherMatrix]];
		obsrInfo = dfit / Sqrt[hatdiag];
		dbeta = DiagonalMatrix[obsrInfo].Transpose[catcherMatrix].
				DiagonalMatrix[parmInfo]
	     ];
	     report = report /. {
		HatDiagonal -> (HatDiagonal -> hatdiag),
		StandardizedResiduals -> (StandardizedResiduals ->
			 standResiduals),
		CookD -> (CookD -> cookD),
		JackknifedVariance -> (JackknifedVariance -> jackvar),
		StudentizedResiduals ->
		 (StudentizedResiduals -> studResiduals),
		PredictedResponseDelta -> (PredictedResponseDelta -> dfit),
		CovarianceMatrixDetRatio -> (CovarianceMatrixDetRatio ->
			 covratio),
		BestFitParametersDelta -> (BestFitParametersDelta -> dbeta) }
	   ]
	];
	
 	(* ==================== compute DurbinWatsonD =================== *)	
	
	If[MemberQ[report, DurbinWatsonD],
	   dw = Apply[Plus,
		 (Drop[fitResiduals-RotateRight[fitResiduals], 1])^2] /
		Apply[Plus, fitResiduals^2];
	   report = report /. DurbinWatsonD -> (DurbinWatsonD -> dw)
	];

	(* ================= compute sequential sum of squares ============ *)

	If[MemberQ[report, SequentialSumOfSquares],
	   If[!ValueQ[sosDiagnostics],
	      sosDiagnostics = p > 1 &&
			TrueQ[Apply[Equal, designMatrixT[[1]]]] &&
	                Apply[And, Map[!TrueQ[Apply[Equal, #]]&,
			 	Drop[designMatrixT, 1] ]]	];
	   If[sosDiagnostics,
                 With[{modelSS = ModelSS[designMatrix, response, weights, tol]},
		   (* compute vector of model sum of squares partition *)
	      	   report = report /. {SequentialSumOfSquares ->
			(SequentialSumOfSquares -> modelSS)} 
		 ],
	         Message[DesignedRegress::seqsos];
	         report = report /. {SequentialSumOfSquares ->
		 	(SequentialSumOfSquares -> $Failed)}	]
	]; (* end SequentialSumOfSquares *)


	(* ================= add predictedResponse to report ============ *)

	If[MemberQ[report, PredictedResponse],
		report = report /. PredictedResponse ->
        	 	(PredictedResponse -> predictedResponse)
	];


	(* ================= add fitResiduals to report ================= *)

	If[MemberQ[report, FitResiduals],
		report = report /. FitResiduals ->
        	 	(FitResiduals -> fitResiduals)
	];


	(* ================= add catcherMatrix to report ================ *)

	If[MemberQ[report, CatcherMatrix],
	   report = report /. CatcherMatrix ->
		 (CatcherMatrix -> MatrixForm[catcherMatrix])
	];


	(* ================= compute variance inflation factors ============ *)

 	If[MemberQ[report, VarianceInflation],
           varinf = Map[(Apply[Plus, #^2])&, catcherMatrix] *
	      If[constantPresent, 	
	  	Map[(Apply[Plus, (#-Mean[#])^2])&, designMatrixT],
	  	Map[(Apply[Plus, #^2])&, designMatrixT]
	      ];
	   report = report /. VarianceInflation -> (VarianceInflation -> varinf)
	];

	(* ================= compute eigenstructure table ============ *)

 	If[MemberQ[report, EigenstructureTable],
	   Block[{weightedDesignMatrixT, pos, carrierscale, scaledDesignMatrix,
		svd, uu, ww, vv, eval, condnum, partition, proportion,
		localNames},
             localNames = names;	
	     (* weightDesignMatrix is n x p, where p is the # of basis fcns,
			including the constant term *)
	     weightedDesignMatrixT = Transpose[weightedDesignMatrix];
	     (* Don't assume that the constant term corresponds to the
			first column of designMatrix. *)
	     If[constantPresent,
	           pos = Position[Transpose[designMatrix],
			 x_List /; Apply[Equal, x]];
	           localNames = Delete[names, pos];	
	           weightedDesignMatrixT = Delete[weightedDesignMatrixT, pos]
	     ];
	     (* now weightedDesignMatrixT is (p-1) x n *)
	     localNames = Map[If[Head[#]=!=String, InputForm[#], #]&,
		 	localNames];
	     weightedDesignMatrixT -= Map[Mean, weightedDesignMatrixT];
	     carrierscale = Map[Sqrt[#.#]&, weightedDesignMatrixT];	
	     scaledDesignMatrix = Transpose[weightedDesignMatrixT /
		 	carrierscale];
	     (* scaledDesignMatrix is n x (p-1); each column is scaled *)	
   	     (* Note the singular values found are identical to the eigenvalues
		of CorrelationMatrix[Transpose[weightedDesignMatrixT]],
		where CorrelationMatrix is defined in
		MultiDescriptiveStatistics.m *)	
	     svd = SingularValues[scaledDesignMatrix];	
	     If[!(Head[svd] === List && Length[svd] === 3),
		    report = report /. EigenstructureTable ->
			(EigenstructureTable -> $Failed),
		   {uu, ww, vv} = svd;   eval = ww^2;   condnum = First[ww]/ww;
		   (* rows of vv are the eigenvectors of Transpose[X].X,
			where X is scaledDesignMatrix *)
		   partition = Map[(#^2/eval)&, Transpose[vv]];
                   (* each row of
                      ( responseVariance partition / carrierscale^2 )
                      gives the partition of the variance of the corresponding
                      parameter *)
		   proportion = Map[(#/Apply[Plus,#])&, partition];
		   (* each row of proportion gives the proportion of the 
		      variance attributable to the various eigenvalues for 
		      the corresponding parameter *)
                   report = report /. EigenstructureTable ->
		     (EigenstructureTable -> TableForm[
		     Transpose[
			Join[{eval, condnum}, Chop[proportion, 10^(-4)]] ],
		     TableHeadings ->
			 {None, Join[{"EigenV", "Index"}, localNames]}
			])
	     ]
	   ] (* end Block *)
        ]; 
  
	(* ================= compute coefficient of variation  ============ *)

 	If[MemberQ[report, CoefficientOfVariation],
	  coeffvar = If[And @@ Map[Positive, weightedResponse],
		      responseMean =
		  	Apply[Plus, weights response] / Apply[Plus, weights];
		      Sqrt[responseVariance] / responseMean,
		      $Failed];
	  report = report /. CoefficientOfVariation ->
	         (CoefficientOfVariation -> coeffvar)
	];

	(* ========== warn user of the changed interpretation of =========== *)
	(* === RSquared and AdjustedRSquared when there is no constant === *)

	If[!constantPresent,
	   If[(MemberQ[report, RSquared] || MemberQ[report, AdjustedRSquared]),
	      Message[DesignedRegress::rsqr] ];
	   If[MemberQ[report, VarianceInflation],
              Message[DesignedRegress::vif] ]
	];	
	

	(* ================================================================ *)
	report = report /. {
		 EstimatedVariance -> (EstimatedVariance -> responseVariance),
        	 BestFitParameters -> (BestFitParameters -> fit),
		 RSquared -> (RSquared -> rSquared),
		 AdjustedRSquared -> (AdjustedRSquared -> adjRSquared),
		 (* Note: BestFit is not a valid output of DesignedRegress
			since it involves basis functions *)
		 BestFit -> (BestFit -> fit.basis)
		};

        report
    ] (* end of vDR *)


(* ======================================================================= *)

(* used by ParameterConfidenceRegion option *)
ellipsoidalLocus[mu_, sigma_] :=
  Module[{esystem, esystemT, sorted, r, dir},
    (* (x-mu).Inverse[sigma].(x-mu) == 1 *)
    If[!FreeQ[esystem = Eigensystem[sigma], Eigensystem], Return[$Failed]];
    (* radii are square roots of eigenvalues *)
    esystemT = Transpose[MapAt[Sqrt[Chop[#]]&, esystem, 1]];
    (* sort semi-axes from largest to smallest *)
    sorted = Sort[esystemT, #1[[1]] > #2[[1]]&];
    {r, dir} = Transpose[sorted];
    Ellipsoid[mu, r, dir]
  ]

(* used to process BasisNames option *)
ChangeValue[list_, name_, newval_]:=
   If[FreeQ[list,name],
      Prepend[list,(name->newval)],
      Map[If[!FreeQ[#,name], Replace[name->_,(name->_)->(name->newval)],#]&,
	  list]
   ]

DiagonalElements[matrix_] :=
    Module[{k, size = Length[matrix]},
        Table[matrix[[k, k]], {k, 1, size}]
    ]

GetOutputList[optionlist_] :=
    Module[{localReport, control, outlist, regReport, summary},
      summary = {ParameterTable, RSquared, AdjustedRSquared, EstimatedVariance,
	ANOVATable};
      If[!(FreeQ[optionlist, OutputControl] && FreeQ[optionlist, OutputList]),
	(* ============= support obsolete options ============ *)
	If[!FreeQ[optionlist, OutputControl],
          Message[DesignedRegress::obs, OutputControl];
      	  control = OutputControl /. optionlist /. {OutputControl -> Automatic}
	];
	If[!FreeQ[optionlist, OutputList],
	  Message[DesignedRegress::obs, OutputList];
	  outlist = OutputList /. optionlist /. {OutputList -> Null}
	];
	If[outlist === Null, localReport = {}, localReport = outlist];
        Which[control === Automatic,
          	localReport = Join[localReport,summary],
              control =!= NoPrint,
          	Message[DesignedRegress::optx,
			 control, OutputControl -> control]
	],
	(* ======= ignore obsolete options and support current option ====== *)
	regReport = RegressionReport /. optionlist;
        (* fixing mistakes user made in RegressionReport spec *)
        If[Head[regReport] === List && MemberQ[regReport, RegressionReport] &&
                regReport =!= {RegressionReport},
           Message[DesignedRegress::regrep];
           regReport = DeleteCases[regReport, RegressionReport]
        ];
        regReport = regReport /. Options[DesignedRegress];
	If[Head[regReport] =!= List,
	   regReport = {regReport}];
        If[MemberQ[regReport, SummaryReport],
	   regReport = DeleteCases[regReport, SummaryReport];
	   (* delete from `summary' items that are individually listed
		in `regReport' *)	
	   Scan[(summary = DeleteCases[summary, #])&,
		Intersection[regReport,
			{ParameterTable, RSquared, AdjustedRSquared,
			EstimatedVariance, ANOVATable}]];
           localReport = Join[regReport, summary],
	   localReport = regReport]
      ]; 
      localReport
    ]

ModelSS[designMatrix_, response_, weights_, tol_] :=        
  Module[{includeConstant, n, p, mean, sqrtweights, weightedDesignMatrixT,
	  weightedResponse, totalSS, partialWeightedDesignMatrix,
 	  svd, u, w, v, pred, scan, modelSS = {}},
    includeConstant = Apply[Equal, Map[#[[1]]&, designMatrix]];
    {n, p} = Dimensions[designMatrix];
    mean = Apply[Plus, response]/n;
    If[VectorQ[weights] && Length[weights]==n &&
		Apply[And, Map[NumberQ, N[weights]]],
       (* numericalize weights *)
       If[Precision[weights] == Infinity, weights = N[weights]];
       sqrtweights = Sqrt[weights],
       sqrtweights = Table[1., {n}]
    ];
    weightedDesignMatrixT = Transpose[sqrtweights designMatrix];
    weightedResponse = sqrtweights response;
    If[includeConstant,
	 p--;
	 totalSS = (response-mean).(response-mean),
	 totalSS = response.response
    ];
    scan = Scan[
         (partialWeightedDesignMatrix =
		 Transpose[Drop[weightedDesignMatrixT, -#]];
          svd = SingularValues[partialWeightedDesignMatrix, Tolerance -> tol];
          If[Head[svd] === List && Length[svd] === 3,
            	{u, w, v} = svd,
                Return[$Failed]
          ];
	  pred = partialWeightedDesignMatrix . (Transpose[v] .
			 DiagonalMatrix[1/w] . u) . response;
          modelSS = AppendTo[modelSS, totalSS-((response-pred).(response-pred))]
         )&, Reverse[Range[p]-1]];
    If[scan === $Failed, Return[$Failed]];
    (* leave the model sum of squares as is for the first nonconstant basis
       function included in model; for each additional basis function,
       compute additional model sum of squares that this function contributes *)
    Scan[(modelSS[[#+1]] = modelSS[[#+1]]-modelSS[[#]])&, Reverse[Range[p-1]]];
    modelSS
  ] 


(* ======================================================================== *)
DesignMatrix::notfcn =
"The second argument to DesignMatrix is not a list of functions."

DesignMatrix::notdata =
"The first argument to DesignMatrix is not a list or rectangular array of data."

DesignMatrix::coord = "The number of independent coordinates (`1`) in the data 
is less than the number of variables (`2`) in the model."

DesignMatrix::crdvl = "One or more of the coordinates `1` has a value.
The coordinates must be symbols."

DesignMatrix::pure =
"Warning: The basis functions include pure functions that will be applied to
the variables from the variable list."

DesignMatrix::invopt =
"Warning: The list `` includes options that are invalid for DesignMatrix."

(* ======================================================================== *)
DesignedRegress::invdm =
"The design matrix is not a numerical matrix."
(* if Regress called DesignedRegress, this typically means that the basis
 functions yielded non-numerical results when evaluated at the data points *)

DesignedRegress::invsvd =
"The singular valid decomposition argument is invalid."

DesignedRegress::invresp =
"The response vector has length `1` but should have length `2`."

DesignedRegress::invrep =
"The item `` is an invalid RegressionReport item for DesignedRegress.
Deleting item from RegressionReport."

DesignedRegress::invopt =
"The option `` is an invalid option for DesignedRegress."

DesignedRegress::bdwghts =
"Warning: Value of option Weights -> `1` is not Automatic, a pure function
mapping a response to a numerical weight, or a numerical vector having
the same length as the data.  Setting all weights to 1 (Weights -> Automatic)."

DesignedRegress::newsvd =
"Warning: DesignedRegress was unable to use the input singular value
decomposition of the design matrix due to unequal weights.  Decomposing
the weighted design matrix."

DesignedRegress::nosvd =
"DesignedRegress was unable to obtain the singular value decomposition for the
design matrix of this problem."

DesignedRegress::rank =
"Warning: the rank of the design matrix is `1`, less than full rank `2`.
Only `1` of the `2` basis functions `3` needed to provide this fit.
Try using a different model or greater precision."

DesignedRegress::badfit =
"Warning: unable to find a fit that is better than the mean response."

DesignedRegress::twocon =
"DesignedRegress determined that two or more constant terms are present 
in the model." (* (This should not happen under normal conditions.) *)

DesignedRegress::uncon =
"DesignedRegress was unable to determine whether a single constant term is
present in the model." (* (This should not happen under normal conditions.) *)

DesignedRegress::seqsos =
"Basis functions are unsuitable for sequential sum of squares analysis.  
Constant should be listed first if included in basis function list."

DesignedRegress::parsos =
"Basis functions are unsuitable for partial sum of squares analysis.  
Constant should be listed first if included in basis function list."

DesignedRegress::obs =
"Warning: the option `` is an obsolete option of Regress and DesignedRegress.
It is superceded by RegressionReport."

DesignedRegress::optx = "Unknown option `1` in `2`."

DesignedRegress::regrep =
"RegressionReport may not be included in the list specified by the
RegressionReport option.  Deleting from list."

DesignedRegress::tsos =
"Warning: the total sum of squares in the ANOVATable is uncorrected (not
centered on the response mean) when there is no constant term in the model;
it is designated U Total."

DesignedRegress::rsqr =
"Warning: the RSquared and AdjustedRSquared diagnostics are redefined when
there is no constant term in the model."

DesignedRegress::vif =
"Warning: the VarianceInflation collinearity diagnostics are redefined when
there is no constant term in the model."

(* ======================================================================== *)
Regress::bdwghts =
"Warning: Value of option Weights -> `1` is not Automatic, a pure function
mapping an observation to a numerical weight, or a numerical vector having
the same length as the data.  Setting all weights to 1 (Weights -> Automatic)."

Regress::notfcn =
"The second argument to Regress is not a list of functions."

Regress::notdata =
"The first argument to Regress is not a list or rectangular array of data."

Regress::coord = "The number of independent coordinates (`1`) in the data is
less than the number of variables (`2`) in the model."

Regress::crdvl = "One or more of the coordinates `1` has a value.  The 
coordinates must be symbols."

Regress::nodm = "Unable to compute design matrix."

Regress::namnum =
"Warning: number of basis function names conflicts with number of basis
functions.  Using default names."

Regress::invnam =
"Warning: invalid BasisNames option value.  Using default names."

Regress::mindata =  "The number of parameters to be estimated is greater than
or equal to the number of data points.  Subsequent results may be misleading."


        End[]


SetAttributes[ Regress ,ReadProtected];
SetAttributes[ DesignedRegress ,ReadProtected];
SetAttributes[ DesignMatrix ,ReadProtected];

Protect[ Regress, DesignedRegress, DesignMatrix];


        EndPackage[]

(* :Example:
x: number of days material is stored
y: amount of methane gas produced

data = {{1, 8.2}, {1, 6.6}, {1, 9.8}, {2, 19.7}, {2, 15.7}, {2, 16.0},
        {3, 28.6}, {3, 25.0}, {3, 31.9}, {4, 30.8}, {4, 37.8}, {4, 40.2},
        {5, 40.3}, {5, 42.9}, {5, 32.6}}

Regress[data, {1, x, x^2}, x, Weights ->
        1/{1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5}^2 ]

{ParameterTable ->

>          Estimate    SE         TStat      PValue     ,
      1    -3.62995    2.53839    -1.43002   0.178228

      x    12.421      2.67432    4.64453    0.000565756

       2
      x    -0.717542   0.520765   -1.37786   0.193398

>    RSquared -> 0.948947, AdjustedRSquared -> 0.940438,

>    EstimatedVariance -> 1.52011,

>    ANOVATable ->         DF   SumOfSq   MeanSq    FRatio    PValue      }
                                                                        -8
                   Model   2    339.059   169.53    111.525   1.77066 10

                   Error   12   18.2413   1.52011

                   Total   14   357.3


Example of undocumented feature... pure function pointed to by Weights
may be a function of independent variable x....


Regress[data, {1, x, x^2}, x, Weights -> (1/#2^2 &)]

	same as above

Regress[data, {1, x, x^2}, {x}, Weights -> (1/#2^2 &)]

	same as above

*)


(* :Example:
x1, x2, y: simulated data

data = Table[(x1 = i;
	      x2 = Prime[i];
	      {x1, x2, 3 x1 + 2 x2 + Random[NormalDistribution[0, 
			Sqrt[x1+x2]/2]]}), {i, 20}];
	
data = {{1, 2, 6.71949}, {2, 3, 11.5972}, {3, 5, 18.0687}, {4, 7, 30.2149}, 
    {5, 11, 36.7307}, {6, 13, 42.6693}, {7, 17, 55.28}, {8, 19, 61.2664}, 
    {9, 23, 72.7406}, {10, 29, 81.9494}, {11, 31, 96.9331}, 
    {12, 37, 108.132}, {13, 41, 119.574}, {14, 43, 130.867}, 
    {15, 47, 138.833}, {16, 53, 157.414}, {17, 59, 177.146}, 
    {18, 61, 175.105}, {19, 67, 191.956}, {20, 71, 212.905}}

Example of undocumented feature... pure function pointed to by Weights
may be a function of independent variables x1 and x2....

Regress[data, {1, x1, x2}, {x1, x2}, Weights -> ((4/(#2 + #3))&)]

{ParameterTable -> 
 
>          Estimate   SE         TStat        PValue      , 
      1    -0.03908   1.05484    -0.0370483   0.970878

      x1   2.70723    0.692708   3.90818      0.00113143

                                                        -9
      x2   2.12577    0.203008   10.4714      7.84933 10
 
>    RSquared -> 0.998064, AdjustedRSquared -> 0.997836, 
 
>    EstimatedVariance -> 1.17182, 
 
>    ANOVATable ->         DF   SumOfSq   MeanSq    FRatio   PValue}
                   Model   2    10268.4   5134.2    4381.4   0.

                   Error   17   19.9209   1.17182

                   Total   19   10288.3


Regress[data, {1, x1, x2}, {x1, x2},
	 Weights -> Table[4/(i + Prime[i]), {i, 20}]]

	same as above
	
*)
