(* :Title: Nonlinear Curve Fitting and Regression Analysis *)

(* :Context: Statistics`NonlinearFit` *)

(* :Name: Statistics`NonlinearFit` *)

(* :Author: John M. Novak & E. C. Martin *)

(* :Summary:
This package performs nonlinear least squares curve fitting and
statistical regression analysis.  
*)

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

(* :Package Version: 1.8 *)

(* :Mathematica Version: 4.0 *)

(* :History:
    Version 1.0, October 1991 by John M. Novak.
    Version 1.1, February 1992 by John M. Novak--major revisions.
    Version 1.2, February 1994 by ECM--major revisions and additions:
	added NonlinearRegress to compute a variety of diagnostics for a
	nonlinear fit; revised output of NonlinearFit to give a function
	(as does the output of Fit); replaced use of LinearSolve and
	termination upon encountering an ill-conditioned "alpha" matrix with 
	a more robust procedure for finding the next step in the parameter
	search; eliminated specifying weights in data matrix, because they can 
	be specified using Weights and the data matrix format conflicts with 
	the obvious representation of multiple-response data; tried to make 
	option interpretation more uniform across methods; added two more
	parameter formats: {parameter, {start0, start1}} and
	{parameter, {start0, start1}, minr, max}... currently these are only
	effective at eliminating the need for symbolic derivatives if
	Method -> FindMinimum;  added support for Gradient option, also
	relevent when Method -> FindMinimum.
    Version 1.3, 1.4, October 1996 by ECM.
    Version 1.5, January 1997 by Sirpa Saarinen -- implemented
	LevenbergMarquardt method using System`Private`SLM. 
    Version 1.6, May 1997 by ECM -- reimplemented LevenbergMarquardt method
	(and methods Gradient, Newton, and QuasiNewton) using FindMinimum.
	The setting Method -> FindMinimum is now obsolete since
	the V4.0 function FindMinimum supports several methods (Automatic,
	Gradient, LevenbergMarquardt, Newton, and QuasiNewton) .  To get the
	V3.0 default behavior of Method -> FindMinimum (steepest descent), use
	Method -> Gradient.  
    Version 1.7, May 1998 by ECM -- 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.
    Version 1.8, July 1998 by ECM -- added support for
	ParameterConfidenceRegion (based on AsymptoticCovarianceMatrix).
	Modified StandardizedResiduals to take account of Weights.
*)

(* :Keywords:
    curve fitting, nonlinear regression
*)

(* :Sources:
    Janhunen, Pekka, NonlinearFit`, (a Mathematica package), 1990.
    Press, William, et. al., "Numerical Recipes in Pascal",
        pp. 572-580, Cambridge University Press (Cambridge, 1989).
    Withoff, Dave, DataAnalysis`NonlinearRegression`,
        (a Mathematica package), 1989.
    Bates and Watts, "Relative Curvature Measures of Nonlinearity",
	J.R.Stat.Soc.B, 42 No.1, pp.1-25.
    Ratkowsky, David A., "Nonlinear Regression Modeling", Marcel Dekker, Inc.
	(New York, 1983).
    Kennedy, W. J. and Gentle, J. E., "Statistical Computing", Marcel Dekker,
	Inc. (New York, 1980).
*)

(* :Discussion:

   The methods LevenbergMarquardt (default), Gradient (steepest descent),
	Newton, and QuasiNewton are passed directly to FindMinimum.
	The method Automatic uses FindMinimum[, Method -> LevenbergMarquardt]
	if the model is nonlinear, but if the model is linear, direct linear 
	fitting is employed.  Thus the meaning of the Automatic method is
        different in FindMinimum and NonlinearFit.

   A note about the LevenbergMarquardt method;
	this method starts the search for the best parameter values
	with direction biased
	towards steepest-descent and gradually shifts the bias in favor of
	quadratic minimization.  In practice, this shift in direction can only
	happen if the regularized Hessian is well-conditioned.  The shift
	towards quadratic minimization is allowed only as long as the
	matrix inverse is numerically stable.

  The options AccuracyGoal and PrecisionGoal
	specify the number of digits of absolute and relative error
	allowed in the residual sum of squares, not the error in the
	parameter estimates.

   What the curvature diagnostics in FitCurvatureTable mean:
    A maximum intrinsic curvature that is small compared to the curvature of the
    joint parameter confidence region indicates that the curved solution locus 
    is well approximated by a plane.  A maximum parameter-effects curvature that
    is small compared to the curvature of the joint parameter confidence region
    indicates that the curved parameter lines on the approximating plane are
    well approximated by a grid of straight, parallel, equispaced lines.  The
    validity of inferences made regarding parameters of the nonlinear model 
    depend on how well the `planar' and `uniform coordinate' assumptions are
    satisfied.  Model reparameterization can reduce the parameter-effects
    curvature, but not the intrinsic curvature.  (The confidence region
    curvature, which provides a basis for comparison, is for a region whose
    confidence level is specified by the option ConfidenceLevel.)

   FindMinimum[f, {x, x0}] searches for a local minimum in f, starting from
	the point x=x0.
   FindMinimum[f, {x, {x0, x1}}] searches using x0 and x1 as the first two
	values of x.
   FindMinimum[f, {x, xstart, xmin, xmax}] stops the search if x ever leaves
	the range xmin to xmax.


The Weights option....

(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 nonlinear model).

(2) Weights -> pureFcn

  (i) In NonlinearRegress[{y1, ..., yn}, model, x, parms, 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 NonlinearRegress[{y1, ..., yn}, model, x, parms,
		 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 NonlinearRegress[{{x1, y1}, ..., {xn, yn}}, model, x, parms,
	 					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 NonlinearRegress[{{x1, y1}, ..., {xn, yn}}, model, x, parms,
                                                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 NonlinearRegress[{{x11, ..., x1k, y1}, ..., {xn1, ..., xnk, yn}},
                     model, {x1, ..., xk}, parms, 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 NonlinearRegress[{{x11, ..., x1k, y1}, ..., {xn1, ..., xnk, yn}},
              model, {x1, ..., xk}, parms, 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.




   Unlike some nonlinear fitting procedures, it is not possible to use the
   Weights option to specify a function of a residual, so that the weights can
   be updated each iteration, and iteratively reweighted least squares
   (preferably using more robust criteria than the l2 norm) can be performed.
   For this to be implemented,  it would be necesary for NonlinearFit
   to call FindMinimum in a loop, undating the weights each iteration.
   FindMinimum handles more general functions than sum-of-squares, so
   adding a Weights option to FindMinimum would not make sense as it would
   apply to only some of the functions that FindMinimum can handle.
*)

(* :Warning:
    It is necessary to use NonlinearRegress[data, model, vars, params,
	RegressionReport->BestFitParameters] to get the list of parameter
	replacement rules previously given by NonlinearFit[data, model, vars,
	params] in V2.2.
    Currently there is no way to utilize an estimate for the standard error
	based on replications.  For example, it should be possible to
	calculate both the maximum relative curvatures in FitCurvatureTable
	and the bias given by ParameterBias using a standard error estimate
	input by the user.  Currently the estimate for the standard error
	is always taken to be Sqrt[EstimatedVariance] which is based on the
	fit residuals.  The capability to input and use a separate estimate
	of the standard error should be added in a future version of the
	package.
*)

BeginPackage["Statistics`NonlinearFit`",
	(* needed for FRatioDistribution and StudentTDistribution: *)
       "Statistics`NormalDistribution`",	
        (* needed for CDF: *)
       "Statistics`Common`DistributionsCommon`",	
        (* needed for StudentTCI: *)
       "Statistics`ConfidenceIntervals`",
        (* needed for Mean: *)
       "Statistics`DescriptiveStatistics`",
	(* needed for Weights, RegressionReport,
           RegressionReportValues, and many symbols allowed by
	   RegressionReportValues[NonlinearRegress] *)
       "Statistics`Common`RegressionCommon`",
	(* needed for Ellipsoid *)
	"Statistics`Common`MultivariateCommon`"]


NonlinearFit::usage =
"NonlinearFit[data, model, vars, params] searches for a least-squares 
fit to a list of data according to the model containing the variables vars
and the parameters params.  Parameters may be expressed as a list of symbols 
or a list of lists.  When expressed as a list, a parameter may be specified
with starting value(s) and bounds in one of several different ways:
{symbol, start} or {symbol, min, max} or {symbol, start, min, max}, where
start can be single value or a list of two values.
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 Method option specifies the
LevenbergMarquardt (default), Gradient (steepest descent), Newton,
QuasiNewton or Automatic search methods.  The Automatic method does linear
fitting for linear models and LevenbergMarquardt nonlinear fitting for 
nonlinear models."

Options[NonlinearFit] =
  Sort[Join[
    {AccuracyGoal -> Automatic,
    Compiled -> True,
    Gradient -> Automatic,
    MaxIterations -> 30,
    Method -> Automatic,
    PrecisionGoal -> Automatic,
    ShowProgress->False,
    Weights -> Automatic,
    WorkingPrecision -> $MachinePrecision},
   Options[SingularValues] ]]

NonlinearRegress::usage =
"NonlinearRegress[data, model, vars, params] searches for a least-squares
fit to a list of data according to the model containing the variables vars
and the parameters params.  Parameters may be expressed as a list of symbols
or a list of lists.  When expressed as a list, a parameter may be specified
with starting value(s) and bounds in one of several different ways:
{symbol, start} or {symbol, min, max} or {symbol, start, min, max}, where
start can be single value or a list of two values.
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 Method option specifies the
LevenbergMarquardt (default), Gradient (steepest descent), Newton,
QuasiNewton or Automatic search methods.  The Automatic method does linear
fitting for linear models and LevenbergMarquardt nonlinear fitting for
nonlinear models."

Options[NonlinearRegress] = Sort[Join[
      {RegressionReport -> SummaryReport},
      Options[NonlinearFit],
      Options[StudentTCI] ]]
(* NOTE: Options[SingularValues] = {Tolerance -> Automatic}, but
	Tolerance is used in more than SingularValues.  Tolerance is also
	used in maxcurve.  Also note that FitCurvatureTable uses
	ConfidenceLevel from Options[StudentTCI]. *)

RegressionReportValues[NonlinearRegress] =
{ANOVATable, AsymptoticCorrelationMatrix, AsymptoticCovarianceMatrix,
BestFit, BestFitParameters,
EstimatedVariance, FitCurvatureTable, FitResiduals, 
HatDiagonal, MeanPredictionCITable, ParameterBias, ParameterConfidenceRegion,
ParameterCITable, ParameterTable, PredictedResponse,
SinglePredictionCITable, StartingParameters, StandardizedResiduals,
SummaryReport};

ShowProgress::usage =
"ShowProgress is an option for NonlinearFit and NonlinearRegress.
If True, a summary of intermediate results for each iteration of the search
for the best parameters will be printed.";
  
If[StringPosition[Gradient::usage, "Nonlinear"] === {},

Gradient::usage = Gradient::usage <> "  " <>
"The functions NonlinearFit and NonlinearRegress also accept the
Method->Gradient setting.";

If[Head[LevenbergMarquardt::usage] === String,
LevenbergMarquardt::usage = LevenbergMarquardt::usage <> "  " <>
"LevenbergMarquardt is also an allowed value for the Method option of
NonlinearFit and NonlinearRegress."
];

Newton::usage = Newton::usage <> "  " <>
"Newton is also an allowed value for the Method option of
NonlinearFit and NonlinearRegress.";

QuasiNewton::usage = QuasiNewton::usage <> "  " <>
"QuasiNewton is also an allowed value for the Method option of
NonlinearFit and NonlinearRegress."

]

	(* symbols unique to nonlinear regression *)

AsymptoticCovarianceMatrix::usage = 
"AsymptoticCovarianceMatrix is used in the output of NonlinearRegress to
identify the estimated covariance matrix of the fit parameters, for the
linear model approximating the original nonlinear model.  It is given by
(EstimatedVariance * Inverse[Transpose[J].J]), where J is the n x p Jacobian 
for the nonlinear model (i.e., the design matrix for the linear approximation
to the nonlinear model)."

AsymptoticCorrelationMatrix::usage = 
"AsymptoticCorrelationMatrix is used in the output of NonlinearRegress to
identify the estimated correlation matrix of the fit parameters, for the
linear model approximating the original nonlinear model."   

(* Instead of having an AsymptoticParameterConfidenceRegion symbol based on
    AsymptoticCovarianceMatrix,  NonlinearRegress instead makes use of the
    symbol ParameterConfidenceRegion, just as Regress does.  I thought that
    the symbol AsymptoticParameterConfidenceRegion was simply too long.  ECM *)

FitCurvatureTable::usage =
"FitCurvatureTable is used in the output of NonlinearRegress to identify a
table of information about the maximum relative intrinsic and parameter-effects
curvatures associated with the least squares solution of a nonlinear fit."

ParameterBias::usage =
"ParameterBias is used in the output of NonlinearRegress to identify a vector
giving the bias in the parameter estimates.  It represents the discrepancy
between the parameter estimates and the true parameter values."

StartingParameters::usage =
"StartingParameters is used in the output of NonlinearRegress to identify
the starting point used in the search for parameter values minimizing the
sum of squares."

        (* updating usage of shared symbols *)

If[StringPosition[SummaryReport::usage, "NonlinearRegress"] === {},

(* SummaryReport exists, but nonlinear regression usage does not. *)
SummaryReport::usage = SummaryReport::usage <> " " <>
"For NonlinearRegress, SummaryReport specifies a summary including
BestFitParameters, ParameterCITable, EstimatedVariance, ANOVATable,
AsymptoticCorrelationMatrix, and FitCurvatureTable.",

(* SummaryReport does not exist. *)
SummaryReport::usage =
"SummaryReport is a value of the option RegressionReport specifying
that NonlinearRegress should return a summary including BestFitParameters,
ParameterCITable, EstimatedVariance, ANOVATable, AsymptoticCorrelationMatrix,
and FitCurvatureTable."

];

If[StringPosition[HatDiagonal::usage, "NonlinearRegress"] === {},

(* HatDiagonal exists, but nonlinear regression usage does not. *)
HatDiagonal::usage = HatDiagonal::usage <> "  " <>
"For NonlinearRegress, the `hat' matrix is based on the asymptotic design matrix
(i.e., the design matrix corresponding to the nonlinear model linearized about
the parameter estimate).";

(* StandardizedResiduals exists, but nonlinear regression usage does not. *)
StandardizedResiduals::usage = StandardizedResiduals::usage <> "  " <>
"For NonlinearRegress, the `hat' matrix used to compute the standardized
residuals is based on the asymptotic design matrix (i.e., the design matrix
corresponding to the nonlinear model linearized about the parameter estimate)."

];


Begin["`Private`"]

(* ======================= NonlinearFit warning messages ================= *)

NonlinearFit::fdin = (* i.e., similar to FindMinimum::fdin *)
"Form of parameter specification `` is inconsistent with other parameter
specifications."

NonlinearFit::bddata =
"The data argument of NonlinearFit must be a vector or a matrix.  Weighted
regression is specified using the Weights option.  The data format
{{x11, x12, ..., {y11, ..., y1m}}, {x21, x22, ..., {y21, ..., y2m}}, ...,
{{xn1, xn2, ..., {yn1, ..., ynm}}} is reserved for multiple response data,
which will be supported in the future."

NonlinearFit::bdvar =
"The variable or variable list `1` of NonlinearFit is invalid."

NonlinearFit::bdits =
"Warning: MaxIterations -> `1` is not set to a positive integer;
setting to 30."

NonlinearFit::bdmtd =
"Warning: Value of option Method -> `1` is not Automatic, Gradient,
LevenbergMarquardt, Newton or QuasiNewton.  Using LevenbergMarquardt method."

NonlinearFit::obsmtd =
"Warning: FindMinimum is an obsolete setting for the option Method;  using
Method -> Gradient (steepest descent) instead."

NonlinearFit::sprog = NonlinearRegress::sprog = 
"Warning: ShowProgress will not display with the chosen method."

NonlinearFit::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)."

NonlinearFit::bdwrk =
"Warning: Value of option WorkingPrecision -> `1` is not an integer or is
an integer less than $MachinePrecision.  Setting to $MachinePrecision."

NonlinearFit::bdprec =
"Warning: Value of option PrecisionGoal -> `1` is not a positive integer
or Automatic.  Setting to $MachinePrecision - 10."

NonlinearFit::bdacc =
"Warning: Value of option AccuracyGoal -> `1` is not a positive integer
or Automatic.  Setting to $MachinePrecision - 10."

NonlinearFit::bdtol =
"Warning: Value of option Tolerance -> `1` is not a number between 0 and 1
or Automatic.  Setting to `2`."

NonlinearFit::dof =
"Warning: The data set is smaller than the number of parameters; the returned
parameter estimates are suspect."

(* The following error message is issued only if Method -> Automatic,
the model is linear, and the estimates found by SingularValues do not fall
within the allowable parameter range. *)
NonlinearFit::lnoutrng =
"Warning: model is linear; the least-squares solution for the linear model,
``, fell outside the allowable parameter range.  Using
LevenbergMarquardt nonlinear search."

NonlinearFit::lnnosvd =
"NonlinearFit was unable to obtain the singular value decomposition for the
design matrix of the linear model."

NonlinearFit::lnzerosvd =
"Either u, w, or v of the singular value decomposition for the design matrix of
the linear model is equal to {}.  Try decreasing the setting for the
Tolerance option.  Tolerance should be set to Automatic or a number between
0 and 1."

NonlinearFit::lnnotdep =
"Warning: the fit to the linear model is numerically independent of a linear
combination of the basis functions."

(* ===================== NonlinearRegress warning messages ================ *)

NonlinearRegress::fdin = (* i.e., similar to FindMinimum::fdin *)
"Form of parameter specification `` is inconsistent with other parameter
specifications."

NonlinearRegress::bddata =
"The data argument of NonlinearRegress must be a matrix.  Weighted regression 
is specified using the Weights option.  The data format
{{x11, x12, ..., {y11, ..., y1m}}, {x21, x22, ..., {y21, ..., y2m}}, ...,
{{xn1, xn2, ..., {yn1, ..., ynm}}} is reserved for multiple response data,
which will be supported in the future."

NonlinearRegress::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)."

NonlinearRegress::nosvd =
"NonlinearRegress was unable to obtain the singular value decomposition for the
design matrix corresponding to the linearized problem."  

NonlinearRegress::noqrd =
"NonlinearRegress was unable to obtain the QR decomposition of the design
matrix corresponding to the linearized problem.  Unable to compute curvature
diagnostics or parameter bias."

NonlinearRegress::noinv =
"NonlinearRegress was unable to obtain an inverse or pseudoinverse for the
R matrix in the QR decomposition of the design matrix corresponding to the
linearized problem.  Unable to compute curvature diagnostics or parameter bias."

NonlinearRegress::pefail =
"Maximum parameter-effects curvature search failed in `` iterations to find the
direction vector that when dotted with the gradient vector gives unity to
within a tolerance of ``.  Parameter-effects curvature returned may not be
maximum value."

NonlinearRegress::infail =
"Maximum intrinsic curvature search failed in `` iterations to find the
direction vector that when dotted with the gradient vector gives unity to
within a tolerance of ``.  Intrinsic curvature returned may not be maximum
value."

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

NonlinearRegress::zerovar =
"`` cannot be computed in the case of zero EstimatedVariance."

NonlinearRegress::zerocov =
"`` cannot be computed in the case of zero elements on the diagonal of
AsymptoticCovarianceMatrix."

NonlinearRegress::rank =
"Warning: Parameter solution locus has dimensionality of only `` at the 
parameter estimate, less than full dimensionality ``.  Confidence region
curvature given by 1/Sqrt[Quantile[FRatioDistribution[``, ``], ``]]."

NonlinearRegress::crvtab =
"Warning:  Computing the orthogonal complement of a `` x `` matrix for
quantities needed by FitCurvatureTable and ParameterBias can be lengthy.
To avoid the wait, abort this calculation and repeat the regression,
omitting FitCurvatureTable and ParameterBias from the items specified by
RegressionReport."

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

NonlinearRegress::bdcl =
"Warning: Value of option ConfidenceLevel -> `1` is not a number between 0 and 1.
Setting to `2`."


(* ========= NonlinearFit & NonlinearRegress fatal error messages ========= *)

NonlinearFit::tmny = NonlinearRegress::tmny =
"Too many variables have been provided for the size of the data;
given `1`-tuples for data, given `2` variables."

NonlinearFit::bdpar = NonlinearRegress::bdpar =
"The form of the parameters `1` given to NonlinearFit is not correct.
Parameters should be of form {parameter}, {parameter, startpoint}, 
{parameter, minrange, maxrange}, {parameter, startpoint, minrange, maxrange},
{parameter, {startpoint0, startpoint1}}, or {parameter,
{startpoint0, startpoint1}, minrange, maxrange}."

NonlinearFit::bdgrid = NonlinearRegress::bdgrid =
"Specifying parameters in the form `1` means that a search of
`2` grid points in the `3`-dimensional parameter space `4` is conducted to find 
the best starting point (i.e., the one minimizing the residual sum of squares).
Unfortunately, each of the `2` points makes the model complex valued at one or
more observations.  Try a different starting point for the parameters."
(* NOTE:
 or try specifying a new model (it is possible that one of the parameters
is not meant to take on real values) *)

NonlinearFit::srange = NonlinearRegress::srange =
"The starting point `1` for the least squares search for parameter estimates
is outside the prescribed parameter range `2`.  Try another starting point
for the parameters."

NonlinearFit::unevldr = NonlinearRegress::unevldr =
"One or more derivatives of the model with respect to parameters did not
evaluate."

NonlinearFit::nonnum = NonlinearRegress::nonnum =
"The model is not numerical at ``.  Check that all model parameters are
included in the parameter list."

NonlinearFit::fitfail = NonlinearRegress::fitfail =
"The fitting algorithm failed."

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

(* NonlinearFit arguments:
    data - vector or matrix of numbers;
    model - the model, a symbolic expression;
    vars - a list of symbols, expressing the variables in model;
    params - a list of symbol and {symbol, startingpoint}
        pairs, expressing the parameters and their starting
        points in the model;
    opts - the options
*)

NonlinearFit[args___] := 
	With[{result = Catch[nonlinearFit[NonlinearFit, args]]},
		result /; result =!= $Failed];

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

(* NonlinearRegress arguments are the same as those of NonlinearFit *)

NonlinearRegress[args___] := 
	With[{result = Catch[nonlinearFit[NonlinearRegress, args]]},
		result /; result =!= $Failed];


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

(* nonlinearFit calls
linfit[npts, nresponse, nweights, basis, vars,
       parameters, nminr, nmaxr, tol]
or
nlfit[method, npts, nresponse, nweights, model, vars,
      parameters, nstart, nstart0, nstart1, nminr, nmaxr, 
      compiled, maxits, wp, pg, ag, grad]
*)
nonlinearFit[msghead_, xdata_, model_, xvars_, xparams_, opts___] :=
    Block[{accgoal, compiled, grad, maxits, method, precgoal, tol, weights,
		workprec, progress,
	    params, np, 
            datalen, varlen,
	    min, max, x, s,
	    pts, tempweights, response, l1, scan,
	    start, start0, start1, minr, maxr, 
	    p, parmSpace, prec,
	    npts, nresponse, nweights, nstart, nstart0, nstart1, nminr, nmaxr,
	    tempRule, tempRule0, constant,
	    nptsT, newptsT, newpts, newvars, newresponse, 
		temp, data = xdata, fitflag = SameQ[msghead, NonlinearFit],
		inputparams = xparams, vars = xvars},
    (* Get options *)
	If[Count[{opts}, _?OptionQ] != Length[{opts}],Throw[$Failed]];
	(* Evaluate method first, using only Method -> method and
	   Method :> method, to avoid problems with conflicting rules
	   such as Method->Gradient and Gradient->Automatic. *)
	method = Method /. Select[{opts},
		MatchQ[#, Rule[Method, _] | RuleDelayed[Method, _]]&] /.
			   Select[Options[msghead],
		MatchQ[#, Rule[Method, _] | RuleDelayed[Method, _]]&];
        {accgoal, compiled, grad, maxits, precgoal, tol, weights,
	 workprec, progress} =
            {AccuracyGoal, Compiled, Gradient, MaxIterations, 
		PrecisionGoal, Tolerance, Weights, WorkingPrecision,
		ShowProgress} /. {opts} /. Options[msghead];
    (* Option and Argument checking *)
	(* Check the variables *)
	If[!ListQ[vars], vars = {vars}];
	varlen = Length[vars];
	Map[
		If[!symbolQ[#],
			Message[MessageName[msghead,"bdvar"],xvars];
			Throw[$Failed]],
		vars];
	(* Check parameter specifications which include starting conditions *)
	If[parameterQ[inputparams], inputparams = {inputparams}];
	Map[
		If[!parameterQ[#],
			Message[MessageName[msghead,"bdpar"], xparams];
			Throw[$Failed]],
		inputparams];
    (* - check input data, separate out response and coords *)
	If[!ListQ[data],
		Message[MessageName[msghead,"bddata"]];
		Throw[$Failed]];
	(* try to pack the data if using machine precision *)
	If[workprec == $MachinePrecision,
		data = Developer`ToPackedArray[data, Real],
	(* else *) 
		data = N[data, workprec]
	];
	datalen = Length[data];
	ddims = Dimensions[data];
	Switch[Length[ddims],
		1,
			nresponse = data;
			temp = Range[First[ddims]];
			pts = Developer`ToPackedArray[Transpose[Table[temp,{varlen}]]];
			If[!fitflag, (* Do we really need this ? *)
				response = xdata];
			npts = N[pts, workprec],
		2,
			If[ddims[[2]] - varlen < 1,
				Message[MessageName[msghead, "tmny"], ddims[[2]], varlen];
				Throw[$Failed]];
			(* - assume univariate response *)
			nresponse = data[[All, -1]];
			npts = Take[data, All, {1,varlen}];
			If[!fitflag, (* Do we really need this ? *)
				response = xdata[[All,-1]];
				pts = Take[data, All, {1,varlen}]
			],
		_,
			Message[MessageName[msghead,"bddata"]];
			Throw[$Failed]];
	If[!Developer`PackedArrayQ[nresponse],
		If[!VectorQ[nresponse, NumberQ],
			Message[MessageName[msghead,"bddata"]];
			Throw[$Failed]]];
	If[!Developer`PackedArrayQ[npts],
		If[!MatrixQ[npts, NumberQ],
			Message[MessageName[msghead,"bddata"]];
			Throw[$Failed]]];

    (* - check and set weights *)
        Which[
         weights === Equal || weights === Automatic,
                weights = 1,
         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],
			If[fitflag, 
                           Message[NonlinearFit::bdwghts, weights],
                           Message[NonlinearRegress::bdwghts, weights]
			];
                        weights = 1,
                        weights = tempweights
                ],
         VectorQ[weights],
                If[Length[weights] =!= Length[nresponse] ||
                   !VectorQ[N[weights], NumberQ],
                        If[fitflag, 
                           Message[NonlinearFit::bdwghts, weights],
                           Message[NonlinearRegress::bdwghts, weights]
                        ];
                        weights = 1
                ],
         True,
		If[fitflag,
                   Message[NonlinearFit::bdwghts, weights],
		   Message[NonlinearRegress::bdwghts, weights]
	        ];	
                weights = 1
        ];


    (* - fill in default values for start, start0, start1, minr, and maxr *)
        l1 = mylength[First[inputparams]];
	scan = Scan[If[mylength[#] != l1,
		       If[fitflag,
		          Message[NonlinearFit::fdin, #],
		          Message[NonlinearRegress::fdin, #]
		       ];
		       Return[$Failed]]&, 
	            Rest[inputparams]]; 
	If[scan === $Failed, Return[$Failed]];
        {params, start, start0, start1, minr, maxr} = Transpose[
            Map[
		(
	        Switch[#,
		  (* no start value *)
                  _?symbolQ,
			 {#, {1}, Null, Null, -Infinity, Infinity},
		  (* no start value *)
		  {_?symbolQ},
			 {#[[1]], {1}, Null, Null, -Infinity, Infinity},
		  (* one start value *)
                  {_?symbolQ, _?numberQ},
			 {#[[1]], {#[[2]]}, Null, Null, -Infinity, Infinity},
		  (* min and max values *)
                  {_?symbolQ, _?numberQ, _?numberQ},
                    	 {#[[1]], startposns[#[[2]],#[[3]]],
				 Null, Null, #[[2]], #[[3]]},
		  (* min and max values (max = Infinity) *)
                  {_?symbolQ, _?numberQ, Infinity},
                    	 {#[[1]], {#[[2]] + 1}, Null, Null, #[[2]], Infinity},
		  (* min and max values (min = -Infinity) *)
                  {_?symbolQ, -Infinity, _?numberQ},
                    	 {#[[1]],  {#[[3]] - 1}, Null, Null, -Infinity, #[[3]]},
		  (* min and max values (-Infinity, Infinity) *)
                  {_?symbolQ, -Infinity, Infinity},
                    	 {#[[1]], {1}, Null, Null, -Infinity, Infinity},
		  (* one start value and min and max values *)
                  {_?symbolQ, _?numberQ,
			 (-Infinity | _?numberQ), (_?numberQ | Infinity)},
                    	 {#[[1]], {#[[2]]}, Null, Null, #[[3]], #[[4]]},
		  (* two start values *)
		  {_?symbolQ, {_?numberQ, _?numberQ}},
		    {#[[1]], {Null}, #[[2, 1]], #[[2, 2]], -Infinity, Infinity},
		  (* two start values and min and max values *)
		  {_?symbolQ, {_?numberQ, s1_?numberQ},
			 (-Infinity | _?numberQ), (_?numberQ | Infinity)},
		    {#[[1]], {Null}, #[[2, 1]], #[[2, 2]], #[[3]], #[[4]]},
		  (* anything else *)
                  _,
		    {$Failed,$Failed,$Failed,$Failed,$Failed} 		]
		)&,
            inputparams]
        ];
        If[ !FreeQ[start, $Failed], (* bad parameter list *)    
            If[fitflag,
               Message[NonlinearFit::bdpar, inputparams],
               Message[NonlinearRegress::bdpar, inputparams]
            ]; 
            Return[$Failed]
	];


	(* Use parmSpace to describe parameter space in error messages. *)
	parmSpace = Map[inequality[#[[1]], #[[2]], #[[3]]]&, 
			 Transpose[{params, minr, maxr}] ];

    (* - grid search for starting point for parameter specifications such as
	{{alpha, alphaMin, alphaMax}, {beta, betaMin, betaMax},
		{gamma, gammaMin, gammaMax}} *)
	If[Dimensions[start][[2]] > 1,
	   (* do a grid search of 2^p locations for the best starting point
		(p is the number of parameters) *)
           start = Flatten[Outer[List,##]& @@ start,Length[start] - 1];
	   If[FreeQ[start, Null], 
              If[Length[start] === 1,
              	start = First[start],
              	start = findbeststart[start, npts, response,
                        weights, model, vars, params];
	        If[start === $Failed,
			p = Length[params];
			If[fitflag,
		 	   Message[NonlinearFit::bdgrid, inputparams,
				2^p, p, parmSpace],
		 	   Message[NonlinearRegress::bdgrid, inputparams,
				2^p, p, parmSpace]
		        ];
			Return[$Failed]
	        ]	
              ]
	   ],
	   (* take single start point out of list *)
	   start = Transpose[start][[1]]
	];

    (* - check that either start or start0 and start1 are within bounds; *)
	If[FreeQ[start, Null],
	   (* check start *)
	   If[!TrueQ[Apply[And, Thread[minr <= start <= maxr] ]],

	      If[fitflag,	
	         Message[NonlinearFit::srange, start, parmSpace],
		 Message[NonlinearFit::srange, start, parmSpace] ];
	      Return[$Failed]	
	   ],
	   (* check start0 and start1 *)
 	   select = Select[{start0, start1},
			((!TrueQ[Apply[And, Thread[minr <= # <= maxr]]])&)];
	   If[select =!= {}, 	
	      If[fitflag,       
                 Message[NonlinearFit::srange, select[[1]], parmSpace],
                 Message[NonlinearFit::srange, select[[1]], parmSpace] ];
	      Return[$Failed]	
	   ]
	];

    (* - check degrees of freedom *)
        If[Length[params] > Length[npts],
            Message[NonlinearFit::dof]			
        ];
    (* - check assorted options *)
        If[!(IntegerQ[maxits] && Positive[maxits]),
           Message[NonlinearFit::bdits, maxits];
           maxits = 30];
        If[!(IntegerQ[workprec] && workprec >= $MachinePrecision),
           Message[NonlinearFit::bdwrk, workprec];	
           workprec = $MachinePrecision];
        If[precgoal === Automatic, precgoal = workprec - 10];
        If[accgoal === Automatic, accgoal = workprec - 10];
        If[!(IntegerQ[precgoal] && Positive[precgoal]),
           Message[NonlinearFit::bdprec, precgoal];	
           precgoal = $MachinePrecision - 10];
        If[!(IntegerQ[accgoal] && Positive[accgoal]),
           Message[NonlinearFit::bdacc, accgoal];	
           accgoal = $MachinePrecision - 10];
		temp = method;
		If[ListQ[method] && Length[method] == 1, 
			temp = method[[1]]];
        If[!MemberQ[{Automatic, FindMinimum, Gradient, 
		LevenbergMarquardt, Newton, QuasiNewton}, temp],
           Message[NonlinearFit::bdmtd, method];	
           method = LevenbergMarquardt];
	If[temp === FindMinimum,
	   Message[NonlinearFit::obsmtd];
	   method = If[temp === method,Gradient, {Gradient}]];
        progress = TrueQ[progress];
	If[progress === True,
			If[!MatchQ[method, Automatic | LevenbergMarquardt | QuasiNewton | Newton],
				Message[MessageName[msghead, "sprog"]];
				progress = False],
			progress = False];
        (* Note that tol is used by linfit and MakeRegressionReport.
	   In linfit, tol is used by SingularValues.
	   In MakeRegressionReport, tol is used by SingularValues and
	   maxcurve. *)
	If[!(tol === Automatic ||
		 (NumberQ[tol] && FreeQ[tol, Complex] && 0 < tol < 1)),
	   Message[NonlinearFit::bdtol, tol, Automatic];
	   tol = Automatic];


    (* - set working precision of inputs to nlfit *)
    (* NOTE: Numericalize the arguments to WorkingPrecision. *)
	nstart = N[start, workprec];
	nstart0 = N[start0, workprec];
	nstart1 = N[start1, workprec];
	nminr = N[minr, workprec];
	nmaxr = N[maxr, workprec];
		
    (* - check whether model is numeric *)
	tempRule0 = Thread[Rule[vars, npts[[1]]]];
	If[FreeQ[start, Null],
		(* there is a single starting point: start *)
		tempRule = Join[ Thread[Rule[params, start]], tempRule0];
		If[!NumberQ[N[ model /. tempRule, workprec]],
		   nonnumMessage[fitflag, tempRule];  Return[$Failed]],
		(* there are two starting points:  start0 and start1 *)
		tempRule = Join[ Thread[Rule[params, start0]], tempRule0];
		If[!NumberQ[N[ model /. tempRule, workprec]],
                   nonnumMessage[fitflag, tempRule];  Return[$Failed]];
		tempRule = Join[ Thread[Rule[params, start1]], tempRule0];
		If[!NumberQ[N[ model /. tempRule, workprec]],
                   nonnumMessage[fitflag, tempRule];  Return[$Failed]]
	];
		

    (* - call proper fitting routine: linfit or nlfit *)

	If[FreeQ[start, Null],
	   (* If start is FreeQ of Null, then the first derivatives are
		to be calculated symbolically, rather than estimated. *)
		If[grad === Automatic,
			grad = Map[D[model, #]&, params]];
	   If[!FreeQ[grad, Derivative],
              (* If we need to calculate derivatives symbolically, but we
			can't, then we fail with a message. *)
	      If[fitflag, 
                 Message[NonlinearFit::unevldr],
	         Message[NonlinearRegress::unevldr]
	      ];	
	      Return[$Failed]
	   ];
 	   (* The first derivatives are calculated successfully. *)
           If[method === Automatic && Apply[And, Map[FreeQ[grad, #]&, params]],
	      (* The Automatic method uses direct linear fitting when the
		 model is determined to be linear.  Note that
		linfit returns $Failed if the solution to the linear problem
		does not satisfy the parameter range limits. *)
	      constant = model;	
	      Scan[(constant = Coefficient[constant, #, 0])&, params];
	      (* Now constant is equal to that portion of the model that
			does not depend on parameters, but may depend on
			variables. *)
	      newresponse = nresponse-constant;
              (* Note that newresponse may now be an expression in terms of
			some of the variables. *)
		If[!NumberQ[constant], 
	      newresponse = Table[newresponse[[i]] /.
		Thread[Rule[vars, npts[[i]]]], {i, Length[newresponse]}]];
	      (* Now newresponse is a numerical vector. *)
	      (* Note that some variables may not be relevent and the
			appropriate columns of npts and the appropriate
			entries of vars should be dropped to give newpts
			and newvars. *)
	      newptsT = {}; nptsT = Transpose[npts]; newvars = {};
	      Scan[If[!FreeQ[grad, vars[[#]]],
				AppendTo[newptsT, nptsT[[#]]];
			        AppendTo[newvars, vars[[#]]]
			     ]&,
					Range[Length[vars]] ]; 
	      newpts = Transpose[newptsT];
              (* grad will be the list of basis functions *)
tt = Timing[
	      out = linfit[newpts, newresponse, weights, grad, newvars,
			 params, nminr, nmaxr, tol, workprec];
"linfit"];print[tt];
	      If[out === $Failed,
	         If[fitflag, 
                    Message[NonlinearFit::fitfail],
	            Message[NonlinearRegress::fitfail]
	         ];	
                 Return[$Failed],
		 (* result of linfit OK *)
	         If[fitflag,
		    (* Return model evaluated at parameter estimates. *)
	            Return[model /. Last[out]],
		    (* Return solution list.  Note that the starting point is
			set to Null, because no starting point is needed in
			linear fitting. *) 
				If[NumberQ[weights], weights = Table[weights,{Length[response]}]];
	            Return[MakeRegressionReport[
			out, pts, response, weights, model,
                	vars, params, Table[Null, {Length[params]}],
			tol, maxits, workprec, 
		        opts]]
	         ] (* end If fitflag *)
	      ] (* end If out =!= $Failed *)
	   ] (* end If method === Automatic and linear model *)
        ]; (* end If FreeQ[start, Null] *)

	     
        (* NOTE: eliminated threading over multiple responses, because
	   Fit doesn't do this either.  Also, this behavior conflicts with the
	   approach of letting multivariate responses yield true multivariate
	   regression (i.e., estimating a covariance matrix for the regression 
	   error vector). *)

	(* The Automatic method uses LevenbergMarquardt fitting when the
		model cannot be determined to be linear. *)
	If[method === Automatic, method = LevenbergMarquardt];
tt = Timing[
        out = nlfit[method, npts, nresponse, weights,
		    model, vars, params,
		    nstart, nstart0, nstart1, nminr, nmaxr,
		    compiled, maxits,
		    workprec, precgoal, accgoal, grad, progress];
"nlfit"];print[tt];

        (* NOTE: "out" has the form {errorSS, parameterRules} *)
	
    (* - check output of fitting routine *)
        If[out === $Failed,
	    If[fitflag, 
               Message[NonlinearFit::fitfail],
	       Message[NonlinearRegress::fitfail]
	    ];	
            Return[$Failed]
        ]; 

	If[fitflag,
	   model /. Last[out],
	(* else *)
		If[NumberQ[weights], weights = Table[weights,{Length[response]}]];
	   MakeRegressionReport[
		out, pts, response, weights, model, vars, params, start,
		tol, maxits, workprec, 
                opts]
	]
	
    ] (* end nonlinearFit *)


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

nonnumMessage[fitflag_, tempRule_] :=
	   If[fitflag,
	      Message[NonlinearFit::nonnum, tempRule],
	      Message[NonlinearRegress::nonnum, tempRule]
	   ];
	
(* ================================= linfit =============================== *)

(* This uses SingularValues for finding the least-squares solution for a
   linear model.  Someday this can be incorporated into the nonlinear search
   by separating model parameters into parameters the model is linear in
   (set theta_l) and parameters the model is nonlinear in (set theta_n).
   Evaluating the model at theta_l = theta_l[0], doing a nonlinear search step
   to obtain theta_n = theta_n[1], solving exactly for theta_l = theta_l[1],
   evaluating the model at theta_l = theta_l[1], etc., is probably superior
   to doing nonlinear search for all parameters without regard to how they
   enter into the model. *)

linfit[npts_, nresponse_, nweights_, basis_, vars_, params_,
	nminr_, nmaxr_, tol_, wp_] :=
   Module[{uv, uvars, rules, pureFunction,
	   ndesignMatrix, pP, nsqrtweights = N[Sqrt[nweights],wp],
	   svd, u, w, v, p = Length[params],
	   catcherMatrixT, fit, fitResiduals},
	uvars = Table[Unique[uv], {Length[vars]}];
	rules = Thread[Rule[vars, uvars]];
        pureFunction := Function[Evaluate[uvars], Evaluate[basis /. rules]];
	ndesignMatrix = Apply[pureFunction, npts, 1];
	If[wp == $MachinePrecision, 
		ndesignMatrix = Developer`ToPackedArray[ndesignMatrix]];
	pP = Internal`EffectivePrecision[ndesignMatrix];
	If[pP === Infinity, ndesignMatrix = N[ndesignMatrix, wp]];
	svd = SingularValues[nsqrtweights ndesignMatrix, Tolerance -> tol];
        If[Head[svd] === List && Length[svd] === 3,
            {u, w, v} = svd,
            Message[NonlinearFit::lnnosvd];
            Return[$Failed]
        ];
	If[u === {} || w === {} || v === {},
	   Message[NonlinearFit::lnzerosvd];
	   Return[$Failed]
	];
        If[Length[w] < p, Message[NonlinearFit::lnnotdep]];
	catcherMatrixT = (nsqrtweights Transpose[u]).(v/w);
	fit = nresponse.catcherMatrixT;
	If[!(And @@ Thread[nminr < fit < nmaxr]),
	   Message[NonlinearFit::lnoutrng, parmRules];
	   Return[$Failed]
        ];
	fitResiduals = nsqrtweights (nresponse - (ndesignMatrix . fit));
	{fitResiduals.fitResiduals, Thread[Rule[params, fit]]}
   ] (* end linfit Module *)



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

(* NEW nonlinear fitting algorithm, based on FindMinimum[, Method->method], 
       where method may be LevenbergMarquardt, Gradient, Newton, or
	QuasiNewton.
*)

nlfit[method_, npts_, nresponse_, nweights_, model_, vars_,
      parameters_, nstart_, nstart0_, nstart1_, nminr_, nmaxr_, 
      compiled_, maxits_, wp_, pg_, ag_, grad_, progress_] :=
Block[{ranges, fm},
	If[FreeQ[nstart, Null],
		ranges = Map[If[FreeQ[#,DirectedInfinity[_]],#,Take[#,2]]&,
			Transpose[{parameters, nstart, nminr, nmaxr}] ],
		ranges = Map[If[FreeQ[#,DirectedInfinity[_]],#,Take[#,2]]&,
			Transpose[{parameters, Transpose[{nstart0, nstart1}],
                                nminr, nmaxr}] ]
	];

	If[!ListQ[method], 
		fm = Internal`FindNonlinearFit[model, vars, ranges, nresponse, npts, Sqrt[nweights], compiled, maxits, wp, pg, ag, method, grad, progress],
	(* else constuct function and use FindMinimum *) 
    Module[{u, uvars, rules, modelF, i, n = Length[nresponse],
            chi, chisq, optionlist},
		(* initialize function *)
			uvars = Table[Unique[u], {Length[vars]}];
			rules = Thread[Rule[vars, uvars]];

		(* =========================================================== *)
			modelF = Function[Evaluate[uvars], Evaluate[model /. rules]];
			chi = nresponse - Apply[modelF, npts, {1}];
			chisq = (nweights chi).chi;
			optionlist = {AccuracyGoal -> ag, Compiled -> compiled,
			MaxIterations -> maxits,
			Method -> method[[1]], PrecisionGoal -> pg, WorkingPrecision -> wp};

		precw = (Head[FindMinimum::precw] == $Off);
			fm = Check[FindMinimum[
                     Evaluate[chisq],
                     Evaluate[Apply[Sequence, Join[ranges, optionlist]]]
                   ],
                   $Failed, FindMinimum::fmcv, FindMinimum::fmlim,
                        FindMinimum::fmgl];
		If[!precw, On[FindMinimum::precw]];
	]];

	If[Head[fm] === List,
		 fm,
		 $Failed]
    ] (* end nlfit Module *)




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

GetRegressionReport[options_] :=
    Module[{localReport, control, outlist, regReport, summary},
      (* NOTE: the summary list for Regress:
         summary = {ParameterTable,RSquared,AdjustedRSquared,EstimatedVariance,
         ANOVATable} *)
      summary = {BestFitParameters, ParameterCITable, EstimatedVariance,
		 ANOVATable, AsymptoticCorrelationMatrix, FitCurvatureTable};
      regReport = RegressionReport /. options;
      (* fixing mistakes user made in RegressionReport spec *)
      If[Head[regReport] === List && MemberQ[regReport, RegressionReport] &&
		regReport =!= {RegressionReport},
	 Message[NonlinearRegress::regrep];
         regReport = DeleteCases[regReport, RegressionReport]
      ];
      regReport = regReport /. Options[NonlinearRegress];
      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,
			 {BestFitParameters, ParameterCITable,
			  EstimatedVariance,
			  ANOVATable, AsymptoticCorrelationMatrix,
		          FitCurvatureTable}
			]];
           localReport = Join[regReport, summary],
           localReport = regReport];
      localReport
    ]

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

(* MakeRegressionReport is called by nonlinearFit. *)

MakeRegressionReport[out_, pts_, response_, weights_, model_, vars_,
	params_, start_, tol_, maxits_, wp_, opts___] :=
	Block[{cLevel,
	       report, errorSS, parameterRules, parameterEst, uv, uvars,
	       variableRules, modelF,
	       modelprime, modelFprime, modelprimeprime, modelFprimeprime,
               INcurve, PEcurve, bias, rinv, locusDOF,	
	       predictedResponse, fitResiduals, totalDOF, modelDOF, errorDOF,
	       responseVariance, npts, ndesignMatrix, pP,
	       sqrtweights = Sqrt[weights],
	       weightedResponse, nweightedResponse,
	       meanWeightedResponse, nmeanWeightedResponse,
	       uncorrectedTotalSS, correctedTotalSS, modelSS,
	       arglist, basisSubset},

	  cLevel = ConfidenceLevel /. {opts} /. Options[NonlinearRegress];
	  If[ !(NumberQ[cLevel] && FreeQ[cLevel, Complex] && 0 < cLevel < 1),
	     Message[NonlinearRegress::bdcl, cLevel, 0.95];
	     cLevel = .95];

	  report = GetRegressionReport[{opts}];

	  {errorSS, parameterRules} = out;
	  parameterEst = params /. parameterRules;
	  uvars = Table[Unique[uv], {Length[vars]}];
	  variableRules = Thread[Rule[vars, uvars]];
	  modelF = Function[Evaluate[uvars],
		Evaluate[(model /. parameterRules) /. variableRules]];
	  If[FreeQ[(modelprime = Map[D[model, #]&, params]), Derivative],
	     (* Dimensions[modelFprime[[2]]] = {p} *)
	     modelFprime = Function[Evaluate[uvars],
		  Evaluate[(modelprime /. parameterRules) /. variableRules]];
	     If[FreeQ[(modelprimeprime = Map[D[modelprime, #]&, params]),
			Derivative],
		(* Dimensions[modelFprime[[2]]] = {p, p} *)
		modelFprimeprime = Function[Evaluate[uvars],
			Evaluate[(modelprimeprime /. parameterRules) /.
				variableRules]],
		(* modelprimeprime is NOT FreeQ of Derivative *)
		{modelprimeprime, modelFprimeprime} = {$Failed, $Failed}
	     ],
	     (* modelprime is NOT FreeQ of Derivative *)
	     {modelprime, modelFprime, modelprimeprime, modelFprimeprime} =
		 {$Failed, $Failed, $Failed, $Failed}
	  ];


	  (* Dimensions[predictedResponse] = {n} *)
	  predictedResponse = Apply[modelF, pts, {1}];
	  fitResiduals = response - predictedResponse;
	  {totalDOF, modelDOF} = {Length[pts], Length[params]};
	  errorDOF = totalDOF - modelDOF;
	  responseVariance = errorSS/errorDOF;

	  (* Dimensions[ndesignMatrix] = {n, p} *)
          pP = Precision[pts];
	  npts = If[pP === Infinity, N[pts, wp], pts];
	  If[modelFprime === $Failed,
	     ndesignMatrix = $Failed,
	     ndesignMatrix = Apply[modelFprime, npts, {1}];
	     pP = Precision[ndesignMatrix];
	     If[pP === Infinity, ndesignMatrix = N[ndesignMatrix, wp]]		
	  ];

	  pP = Internal`EffectivePrecision[response];
	  weightedResponse = sqrtweights response;
	  nweightedResponse = If[pP === Infinity, 
		N[weightedResponse, wp], weightedResponse];
	  meanWeightedResponse = 
		sqrtweights (response.weights) / Apply[Plus,weights];
          nmeanWeightedResponse = If[pP === Infinity,
		N[nmeanWeightedResponse, wp], meanWeightedResponse];
	  
          uncorrectedTotalSS = nweightedResponse.nweightedResponse;
	  modelSS = uncorrectedTotalSS - errorSS;
	  If[modelSS < 0,
             Message[NonlinearRegress::bdfit]];

	  correctedTotalSS = (nweightedResponse - nmeanWeightedResponse).
		(nweightedResponse - nmeanWeightedResponse);
	  report = report /. {
		EstimatedVariance ->
		 (EstimatedVariance -> responseVariance),
		BestFit -> (BestFit -> (model /. parameterRules)),
		BestFitParameters ->
	         (BestFitParameters -> parameterRules),
		StartingParameters -> (StartingParameters ->
			If[FreeQ[start, Null],
		 	   (* single start point *)
			   Thread[Rule[params, start]],
			   (* two start points *)
			   {Thread[Rule[params, start0]],
			    Thread[Rule[params, start1]]} ])
	       };

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

	  (* Only add to report if necessary.... depending on the length of
		the data, predictedResponse could be a very long vector. *)
	  If[MemberQ[report, PredictedResponse],
                report = report /. PredictedResponse ->
                        (PredictedResponse -> predictedResponse)
          ];


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

	  (* Only add to report if necessary.... depending on the length of
		the data, fitResiduals could be a very long vector. *)
          If[MemberQ[report, FitResiduals],
                report = report /. FitResiduals ->
                        (FitResiduals -> fitResiduals)
          ];


	  (* ===================== compute ANOVA table ===================== *)
	
          If[MemberQ[report, ANOVATable],
             report = report /. ANOVATable -> (ANOVATable ->
		 TableForm[{{modelDOF, modelSS, modelSS/modelDOF},
               	 	    {errorDOF, errorSS, errorSS/errorDOF},
               	 	    {totalDOF, uncorrectedTotalSS},
	       		    {totalDOF-1, correctedTotalSS}},
              		   TableHeadings->
		{{"Model", "Error", "Uncorrected Total", "Corrected Total"},
                 {"DF", "SumOfSq", "MeanSq"}} ])
          ];


	  (* ================== compute things involving ================== *)
	  (* ============= the SVD of the asymp design matrix ============= *)

	  If[Or @@ Map[MemberQ[report, #]&,
		{AsymptoticCovarianceMatrix, AsymptoticCorrelationMatrix,
	         ParameterConfidenceRegion,
		 ParameterConfidenceRegion[_?VectorQ],
		 HatDiagonal, ParameterTable, ParameterCITable,
		 SinglePredictionCITable, MeanPredictionCITable,
		 StandardizedResiduals}],

	   (* Regardless of whether ndesignMatrix === $Failed or the
	      singular value decomposition of nweightedDesignMatrix works
              or not, we need to know for which subsets of the
              parameters joint confidence regions are required.  If
              ndesignMatrix and the singular value decomposition work,
	      we will return  meaningful confidence regions.  If not, we will
	      return $Failed. *)
           If[MemberQ[report, ParameterConfidenceRegion[_?VectorQ]],
                  arglist = Map[First, Select[report,
                              MatchQ[#,
                                 ParameterConfidenceRegion[l_List] /;
                                 Apply[And, Map[MemberQ[params, #]&, l]]]& ]];
           ];


	   If[ndesignMatrix =!= $Failed,
	     Block[{names, nweightedDesignMatrix, 
		    svd, u, w, v, covMatrix, corMatrix, hatMatrix,
		    hatdiag, standResiduals, parmci, predci,
		    parmse, (* list of asympt. stan. errors of parm estimates *)
		    predse, (* list of asympt. stan. errors of pred response *)
		    tRatio, tPValues},

	       (* names = Map[ToString, params]; *) 	
	       names = params; 
	       nweightedDesignMatrix = sqrtweights ndesignMatrix;
	       pP = Precision[nweightedDesignMatrix];
	       If[pP === Infinity,
		  nweightedDesignMatrix = N[nweightedDesignMatrix, wp]];	
	       svd = SingularValues[nweightedDesignMatrix, Tolerance -> tol];

	       If[!(Head[svd] === List && Length[svd] === 3),
		  
		  (* BAD SVD *)
                  Message[NonlinearRegress::nosvd];
		  (* Need to assign everything that depended on a 
			good singular value decomposition to $Failed. *)
		  report = report /. {
		  	AsymptoticCovarianceMatrix -> 
		    		(AsymptoticCovarianceMatrix -> $Failed),
		  	AsymptoticCorrelationMatrix ->
		    		(AsymptoticCorrelationMatrix -> $Failed),
			HatDiagonal -> (HatDiagonal -> $Failed),
			MeanPredictionCITable ->
				(MeanPredictionCITable -> $Failed),
			ParameterCITable ->
				(ParameterCITable -> $Failed),
		  	ParameterConfidenceRegion ->
				(ParameterConfidenceRegion -> $Failed),
		  	ParameterTable ->
				(ParameterTable -> $Failed),
		  	SinglePredictionCITable ->
				(SinglePredictionCITable -> $Failed),
	          	StandardizedResiduals ->
				(StandardizedResiduals -> $Failed)
			};
                  Scan[(basisSubset = #;
                        report = report /.
                                   (ParameterConfidenceRegion[
					basisSubset] ->
                                   (ParameterConfidenceRegion[
					basisSubset] -> $Failed))
                                )&, arglist],

		  (* GOOD SVD *)
                  {u, w, v} = svd;
	          covMatrix = responseVariance (Transpose[v] . (v / w^2));
	          report = report /. AsymptoticCovarianceMatrix ->
		       (AsymptoticCovarianceMatrix -> MatrixForm[covMatrix]);
	          parmse = Sqrt[DiagonalElements[covMatrix]];
		  If[Apply[Or, Map[TrueQ[#==0]&, parmse]],
		     Message[NonlinearRegress::zerocov,
                         AsymptoticCorrelationMatrix];
		     report = report /. AsymptoticCorrelationMatrix ->
			(AsymptoticCorrelationMatrix -> $Failed),
		     corMatrix = Transpose[Transpose[covMatrix/parmse] /
			 parmse];
		     report = report /. AsymptoticCorrelationMatrix ->
			(AsymptoticCorrelationMatrix -> MatrixForm[corMatrix])
		  ];	
	   

	          (* NOTE: Calculate hatMatrix only when you know that the
			   HatDiagonal diagnostic is needed.
			   hatMatrix is an n x n matrix, which can be
			   quite large. *)
		  If[MemberQ[report, HatDiagonal] ||
		     MemberQ[report, StandardizedResiduals],
		     (* NOTE:  The following manner of calculating hatMatrix
			        unnecessarily constructs
			        DiagonalMatrix[sqrtweights]
			        (an n x n matrix) and is slow.
	               hatMatrix = ndesignMatrix .
				(Transpose[v].DiagonalMatrix[1/w].u.
			 	DiagonalMatrix[sqrtweights])
		     *)
		     hatMatrix = Map[sqrtweights # &, u];
		     hatMatrix =
			 Transpose[v].DiagonalMatrix[1/w].hatMatrix;
		     hatMatrix = ndesignMatrix . hatMatrix;
                     hatdiag = DiagonalElements[hatMatrix];
		     (* 8/98: changed standResiduals from
                        fitResiduals/Sqrt[responseVariance (1-hatdiag)]
                        to
                        fitResiduals/Sqrt[responseVariance (1-hatdiag)/weights]
                     *)
		     standResiduals = If[TrueQ[responseVariance==0],
                          Message[NonlinearRegress::zerovar,
                                 StandardizedResiduals];
                          $Failed,
                          fitResiduals /
				 Sqrt[responseVariance (1-hatdiag)/weights]
                     ]; 
		     report = report /. {HatDiagonal ->
			(HatDiagonal -> hatdiag),
                  			 StandardizedResiduals ->
			 (StandardizedResiduals -> standResiduals)}
		  ]; 	(* end HatDiagonal & StandardizedResiduals *)
		  If[MemberQ[report, ParameterTable],
		        tRatio = Map[If[ TrueQ[#[[2]]==0],
					 Infinity, #[[1]]/#[[2]] ]&,
				     Transpose[{parameterEst, parmse}] ];
                	tPValues = Map[CDF[StudentTDistribution[errorDOF],#]&,
				 tRatio];
                	tPValues = 2 Map[If[#>.5, 1-#, #]&, tPValues];
                	report = report /. ParameterTable ->
                        	(ParameterTable -> TableForm[
                           Transpose[{parameterEst, parmse, tRatio, tPValues}],
                        	TableHeadings -> {names,
                         	{"Estimate", "Asymp. SE", "TStat", "PValue"}}])
            	  ];	(* end ParameterTable *)
	          If[MemberQ[report, ParameterCITable],
                     parmci = Apply[N[StudentTCI[#1, #2, errorDOF,
                                   ConfidenceLevel -> cLevel]]&,
                           Transpose[{parameterEst, parmse}], 1];
                     report = report /. ParameterCITable ->
                        (ParameterCITable -> TableForm[
                          Transpose[{parameterEst, parmse, parmci}], 
			  TableDepth -> 2,
                          TableHeadings -> {names,
				 {"Estimate", "Asymptotic SE", "CI"}}])
                  ];	(* end ParameterCITable *)
		  If[MemberQ[report, MeanPredictionCITable] ||
                     MemberQ[report, SinglePredictionCITable],
		     pP = Precision[response];
		     nresponse = If[pP === Infinity, N[response, wp], response];
                     predse = Sqrt[Map[(#.covMatrix.#)&, ndesignMatrix]] ];
                  If[MemberQ[report, MeanPredictionCITable],
                     predci = Apply[N[StudentTCI[#1, #2, errorDOF,
                                ConfidenceLevel -> cLevel]]&,
                                Transpose[{predictedResponse, predse}], 1];
                     report = report /. MeanPredictionCITable ->
                      (MeanPredictionCITable -> TableForm[
                      Transpose[{nresponse, predictedResponse, predse, predci}],
                      		TableDepth -> 2, TableHeadings -> {None,
                          {"Observed", "Predicted", "Asymptotic SE", "CI"}}])
                  ];	(* end MeanPredictionCITable *)
                  If[MemberQ[report, SinglePredictionCITable],
                     predse = Sqrt[responseVariance + predse^2];
                     predci = Apply[N[StudentTCI[#1, #2, errorDOF,
                                ConfidenceLevel -> cLevel]]&,
                                Transpose[{predictedResponse, predse}], 1];
                     report = report /. SinglePredictionCITable ->
                      (SinglePredictionCITable -> TableForm[
                      Transpose[{nresponse, predictedResponse, predse, predci}],
                      		TableDepth -> 2, TableHeadings -> {None,
                          {"Observed", "Predicted", "Asymptotic SE", "CI"}}])
                  ];	(* end SinglePredictionCITable *)
	          If[MemberQ[report, ParameterConfidenceRegion] ||
		     MemberQ[report,
			 ParameterConfidenceRegion[_?VectorQ]],
                     (* Weisberg, p. 97 *)
                     Block[{ellipsoid, pos,
			 fit = parameterEst, p = modelDOF, n = totalDOF,
                         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]],
                           Scan[(basisSubset = #;
                                 pos = Flatten[Map[Position[params, #, 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 If SVD *)


	     ], (* end Block *)



	     (* ndesignMatrix === $Failed *)
	     report = report /. {
               AsymptoticCovarianceMatrix ->
                (AsymptoticCovarianceMatrix -> $Failed),
               AsymptoticCorrelationMatrix ->
                (AsymptoticCorrelationMatrix -> $Failed),
               HatDiagonal -> (HatDiagonal -> $Failed),
	       MeanPredictionCITable -> (MeanPredictionCITable -> $Failed),	
	       ParameterCITable -> (ParameterCITable -> $Failed),
	       ParameterConfidenceRegion ->
		(ParameterConfidenceRegion -> $Failed),
	       ParameterTable -> (ParameterTable -> $Failed),
	       SinglePredictionCITable -> (SinglePredictionCITable -> $Failed),
               StandardizedResiduals -> (StandardizedResiduals -> $Failed)
				 };
             Scan[(basisSubset = #;
                   report = report /.
                                   (ParameterConfidenceRegion[
					basisSubset] ->
                                   (ParameterConfidenceRegion[
					basisSubset] -> $Failed))
                                )&, arglist]
	   ] (* end If ndesignMatrix =!= $Failed *)
	
	  ];	(* end of SVD of asymptotic design matrix calculations *)
	

	  (* ============= compute curvature diagnostics table =========== *)
	
          If[MemberQ[report, FitCurvatureTable] ||
	     MemberQ[report, ParameterBias],
	   If[!(responseVariance==0 || modelprimeprime===$Failed),
	    If[modelprimeprime === Table[0, {modelDOF}, {modelDOF}],
	      {INcurve, PEcurve} = {0, 0};
              locusDOF = modelDOF;	
	      bias = Table[0, {modelDOF}],
	      (* modelprimeprime is NOT a zero matrix *)
	      Block[{radius, V1, V2, qrd, q, r, rsvd, u, w, v,
		     accel, digits, maxcurvetol},
		radius = Sqrt[modelDOF responseVariance];
	        If[!FreeQ[radius, Sqrt], radius = N[radius, wp]];
	        V1 = radius^(-1) ndesignMatrix;
		V2 = radius^(-1) Apply[modelFprimeprime, npts, {1}];
		pP = Precision[V2];
		If[pP === Infinity, V2 = N[V2, wp]];
		If[Head[qrd = QRDecomposition[V1]] =!= List,
        	   Message[NonlinearRegress::noqrd];
		   rinv === $Failed,
	           (* {p x n, p x p} OR {k x n, k x p} where k < p *)	
		   {q, r} = qrd;
		   If[(locusDOF = Length[r]) == modelDOF,
	              rinv = Inverse[r];
	              If[Head[rinv] === Inverse,
			 Message[NonlinearRegress::noinv];
			 rinv = $Failed
		      ],
		      Message[NonlinearRegress::rank, locusDOF, modelDOF,
			locusDOF, errorDOF, cLevel];
		      If[!(Head[rsvd = SingularValues[r, Tolerance->tol]] ===
			   List && Length[rsvd] == 3),
			 Message[NonlinearRegress::noinv];
			 rinv = $Failed,
		         {u, w, v} = rsvd;
			 rinv = Transpose[v].DiagonalMatrix[1/w].u
		      ]
		   ]
		]; (* end calc of rinv, q, locusDOF *)
		If[rinv =!= $Failed,
		   If[(locusDOF totalDOF - locusDOF(locusDOF+1)/2) > 500,
		      (* NOTE: for large problems,
			 inform user what the wait is about. *)	
		      Message[NonlinearRegress::crvtab, Length[q],
				Length[q[[1]]] ] ];
		   q = Join[q, OrthogonalComplement[q]]; (* n x n *) 	
		   (* {n x p x p} OR {n x k x k} where k < p *)
		   accel =   q.Map[(Transpose[rinv].#.rinv)&, V2  ];
		   digits = Max[Precision[accel], $MachinePrecision];
		   maxcurvetol = If[tol === Automatic, 
		      N[10^(-digits+10), digits-10],
		      tol ];
		   INcurve = maxcurve[Drop[accel, locusDOF],
			 maxcurvetol, maxits];
		   If[Head[INcurve] === List,
	              Message[NonlinearRegress::infail, maxits, maxcurvetol];
		      INcurve = INcurve[[2]]
		   ];
		   PEcurve = maxcurve[Take[accel, locusDOF],
			 maxcurvetol, maxits];
		   If[Head[PEcurve] === List,
        	      Message[NonlinearRegress::pefail, maxits, maxcurvetol];
		      PEcurve = PEcurve[[2]]
		   ];
		   bias = (-1/(2 modelDOF)) rinv . Map[
				Apply[Plus, DiagonalElements[#]]&,
					Take[accel, locusDOF] ]
		]	(* end If rinv =!= $Failed *)
	      ]	(* end Block calc of rinv, INcurve, PEcurve, bias,
			   and locusDOF *)
	    ]; (* end If modelprimeprime === Table[0, {modelDOF}, {modelDOF}] *)
	    If[rinv === $Failed,
	      report = report /. {FitCurvature ->
			(FitCurvatureTable->$Failed),
		   ParameterBias -> (ParameterBias->$Failed)},
	      report = report /. {FitCurvatureTable ->
			(FitCurvatureTable -> TableForm[
			   {{INcurve}, {PEcurve},
		 	    {1/Sqrt[Quantile[FRatioDistribution[locusDOF,
				 errorDOF], cLevel]]}}, 
			 TableDepth -> 2, TableHeadings -> {
		  	{"Max Intrinsic", "Max Parameter-Effects",
		   	ToString[100 cLevel]<>" % Confidence Region"},
		  	{"Curvature"}}]),
		   ParameterBias -> (ParameterBias->bias)}
	    ], (* end If rinv === $Failed *)
	    (* responseVariance==0 || modelprimeprime===$Failed *)
	    report = report /. {FitCurvatureTable ->
			(FitCurvatureTable->$Failed),
			ParameterBias -> (ParameterBias->$Failed)}
	   ] (* end If !(responseVariance==0 || modelprimeprime===$Failed) *)
	  ];	(* end If MemberQ[report, FitCurvatureTable] ||
			  MemberQ[report, ParameterBias] *)

          If[TrueQ[responseVariance==0],
	     If[MemberQ[report, FitCurvatureTable],
                Message[NonlinearRegress::zerovar, FitCurvatureTable];	
	        report = report /. FitCurvatureTable ->
			(FitCurvatureTable -> $Failed)];
	     If[MemberQ[report, ParameterBias],
                Message[NonlinearRegress::zerovar, ParameterBias];	
	        report = report /. ParameterBias ->
			(ParameterBias -> $Failed)]
          ];    (* end If responseVariance==0 *)
	
	  report

	] (* end MakeRegressionReport *)

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

(* maxcurve returns the maximum curvature of an acceleration matrix when
   convergence occurs or {$Failed, approximate maximum curvature} if the
   maximum # of iterations is exceeded. *)

maxcurve[A_, tol_, maxiter_] :=
  Module[{dir, grad, Ad, dAd, thresh = 1 - tol, mc = $Failed},
    (* A is a m x p x p tensor, where p = # of parameters, m = d.o.f. *)

    (* initialize parameter space direction vector *)
    dir = Append[Table[0, {Length[ A[[1,1]] ] - 1}], 1];

    Do[  (* set grad to [dirT A dir]T [A dir] and normalize *)
         Ad = Map[(#.dir)&, A];      (* m x p matrix *)
         dAd = Map[(dir.#)&, Ad];    (* m vector *)
         grad = dAd . Ad;            (* p vector *)
         grad /= Sqrt[grad.grad];    (* normalize *)

         (* at convergence, dir points along grad *)
         If[ Abs[grad.dir] <= thresh,
           dir = 3 grad + dir;       (* update dir *)
           dir /= Sqrt[dir.dir],     (* normalize *)
	   mc = Sqrt[dAd.dAd];
	   Break[]
	 ],

         {maxiter}];

    If[mc === $Failed,
      {$Failed, Sqrt[dAd.dAd]},
      mc] 
  ]

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

(* OrthogonalComplement[vectors] computes the (n-p) length n vectors
        that form the orthogonal complement to the p length n
        orthonormal input vectors.
   The ith input vector is rotated in the plane defined by the ith and jth
   axes, until it lines up with the ith axis.  This same rotation is
   performed on all input vectors, as well as the rotation
   matrix (originally the n x n identity matrix).  A rotation is done for
   entries j=i+1 to j=n in the ith vector until the ith vector is equal to
   the ith axis.  Since the vectors are orthonormal to begin with, any
   rotations will maintain this property.  Once the ith vector is lined up
   with the ith axis, none of the other vectors will have a component in the
   ith direction.
   This same step is performed for each of the p input vectors.
   The resulting rotation matrix gives the inverse of the matrix composed
   of the original p vectors plus the (n-p) orthogonal complement.  In the
   case of orthonormal matrices, Transpose gives the Inverse.
*)
OrthogonalComplement[vectors_] :=
 Module[{p = Length[vectors], n = Length[vectors[[1]]],
        tempvT = Transpose[vectors], rotationT, i, j, theta, cos, sin, rot},

        (* rotation matrix starts out as identity *)
        rotationT = IdentityMatrix[n];

        (* line up each vector in turn: *)
        Do[
          (* zero out each vector element by rotating the vectors: *)
          Do[
	    (
            If[tempvT[[j,i]] != 0,
              (* compute rotation angle to zero jth element in ith vector *)
              theta = ArcTan[tempvT[[i,i]],tempvT[[j,i]]];
		cos = Cos[theta];  sin = Sin[theta];
		rot = {{cos, sin}, {-sin, cos}};
              (* rotate all the vectors by this amount *)
	      tempvT[[{i, j}]] = rot.tempvT[[{i, j}]];
              (* rotate the matrix also *)
	      rotationT[[{i, j}]] = rot.rotationT[[{i, j}]]
            ] ),
          {j,i+1,n}],
        {i,p}];

        Drop[rotationT, p]

 ] (* end OrthogonalComplement *)

(* ====================== miscellaneous utility functions =================== *)

mylength[x_] := If[Head[x]===List,
		   Length[x],
		   (* assume symbolQ *)
		   1]

numberQ[n_] := NumberQ[N[n]]

parameterQ[param_] :=
 	     MatchQ[param, _?symbolQ] ||
	     MatchQ[param, {_?symbolQ, start_?numberQ}] ||
             MatchQ[param, {_?symbolQ,
		 min_?((numberQ[#] || # === -Infinity)&),
	 	 max_?((numberQ[#] || # === Infinity)&)}] ||
	     MatchQ[param, {_?symbolQ, start_?numberQ,
		 min_?((numberQ[#] || # === -Infinity)&),
	 	 max_?((numberQ[#] || # === Infinity)&)}] ||
	     MatchQ[param, {_?symbolQ, {start0_?numberQ, start1_?numberQ}}] ||
	     MatchQ[param, {_?symbolQ, {start0_?numberQ, start1_?numberQ},
		 min_?((numberQ[#] || # === -Infinity)&),
                 max_?((numberQ[#] || # === Infinity)&)}]

symbolQ[x_] := (Head[x] === Symbol) ||
		(!numberQ[x] && !ListQ[x]) (* can't check for head Symbol 
				because what about Subscript[a, 1] ? *)

(* utility to choose starting positions based on a minimum and a maximum
    for the parameter *)

startposns[min_,max_] := 
	(
	{min + (1/3)(max - min), min + (2/3)(max - min)}
	)

(* utility to choose a starting position from those derived with the
    above utility by picking the point with the smallest chi-square. *)

findbeststart[starts_, data_, response_, weights_, model_, vars_, params_] :=
    Module[{up, uparams, paramRules, uv, uvars, varRules, umodel,
	    chis, chisq, pairs},
	(* need to use Unique symbols because Function doesn't like
		variables like a[1] or Subscript[a, 1] *)
	uparams = Table[Unique[up], {Length[params]}];
	paramRules = Thread[Rule[params, uparams]];
	uvars = Table[Unique[uv], {Length[vars]}];
	varRules = Thread[Rule[vars, uvars]];
	umodel = model /. paramRules /. varRules;
        chisq = Function[Evaluate[uparams], 
           Evaluate[Plus @@ MapThread[
             (#3 (#2 - Function[Evaluate[uvars], Evaluate[umodel]] @@ #1)^2) &,
             {data,response,weights}]
           ]
        ];
        pairs = Transpose[{starts,chis = N[Apply[chisq, starts,{1}]]}];
	pairs = Select[pairs, FreeQ[#[[2]], Complex]&];
	If[pairs === {},
		$Failed,
        	First[First[Select[pairs, #[[2]] == Min[chis] &]]]
	]
    ]

(* 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]
  ]


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



inequality[parm_, min_, max_] :=
		Switch[{min, max},
			{-Infinity, Infinity}, -Infinity < parm < Infinity,
			{-Infinity, _}, -Infinity < parm <= max,
			{_, Infinity}, min <= parm < Infinity,
			{_, _}, min <= parm <= max]


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

End[]

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}}

NonlinearRegress[data, b0 + b1 x + b2 x^2, x, {b0, b1, b2}, Weights ->
	1/{1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5}^2 ]

{BestFitParameters -> {b0 -> -3.62995, b1 -> 12.421, b2 -> -0.717542},
 
>    ParameterCITable -> 
 
>          Estimate    Asymptotic SE   CI                  , 
      b0   -3.62995    2.53839         {-9.16063, 1.90072}

      b1   12.421      2.67432         {6.59412, 18.2478}

      b2   -0.717542   0.520765        {-1.85219, 0.417108}
 
>    EstimatedVariance -> 1.52011, 
 
>    ANOVATable ->                     DF   SumOfSq   MeanSq , 
                   Model               3    1115.33   371.775

                   Error               12   18.2413   1.52011

                   Uncorrected Total   15   1133.57

                   Corrected Total     14   357.3
 
>    AsymptoticCorrelationMatrix -> 1.          -0.960244   0.901754 , 

                                    -0.960244   1.          -0.975052

                                    0.901754    -0.975052   1.
 
>    FitCurvatureTable ->                           Curvature}
                          Max Intrinsic             0

                          Max Parameter-Effects     0

                          95. % Confidence Region   0.535265

Compare with result from Regress...

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....

NonlinearRegress[data, b0 + b1 x + b2 x^2, x, {b0, b1, b2}, Weights ->
	(1/#2^2 &)]

*)


(* :Example:
x: time
y: radiation count

data = {{0, 383.}, {14, 373.}, {43, 348.}, {61, 328.}, {69, 324.}, {74, 317.},
	{86, 307.}, {90, 302.}, {92, 298.}, {117, 280.}, {133, 268.},
	{138, 261.}, {165, 244.}, {224, 200.}, {236, 197.}, {253, 185.},
	{265, 180.}, {404, 120.5}, {434, 112.5}};

Regress[data, {1, time, time^2}, time]

{ParameterTable -> 
 
>             Estimate      SE             TStat      PValue       , 
      1       386.918       1.11748        346.241    0.

      time    -1.01988      0.0134187      -76.0041   0.

          2                                                     -15
      time    0.000891277   0.0000303768   29.3407    2.44249 10
 
>    RSquared -> 0.999428, AdjustedRSquared -> 0.999357, 
 
>    EstimatedVariance -> 4.07858, 
 
>    ANOVATable ->         DF   SumOfSq   MeanSq    FRatio    PValue}
                   Model   2    114087.   57043.3   13986.1   0.

                   Error   16   65.2572   4.07858

                   Total   18   114152.


This model is linear in the parameters, so NonlinearRegress is not needed...

NonlinearRegress[data, b0 + b1 time + b2 time^2, time, {b0, b1, b2}]

{BestFitParameters -> 
 
>     {b0 -> 386.918, b1 -> -1.01988, b2 -> 0.000891277}, 
 
>    ParameterCITable -> 
 
>          Estimate      Asymptotic SE   CI                        , 
      b0   386.918       1.11748         {384.549, 389.287}

      b1   -1.01988      0.0134187       {-1.04832, -0.991432}

      b2   0.000891277   0.0000303768    {0.000826881, 0.000955673}
 
>    EstimatedVariance -> 4.07858, 
 
>    ANOVATable ->                     DF   SumOfSq       MeanSq , 
                                                      6
                   Model               3    1.44465 10    481551.

                   Error               16   65.2572       4.07858

                                                      6
                   Uncorrected Total   19   1.44472 10

                   Corrected Total     18   114152.
 
>    AsymptoticCorrelationMatrix -> 1.          -0.860713   0.73743  , 

                                    -0.860713   1.          -0.956687

                                    0.73743     -0.956687   1.
 
>    FitCurvatureTable ->                           Curvature}
                          Max Intrinsic             0

                          Max Parameter-Effects     0

                          95. % Confidence Region   0.555652

*)

(* :Example:
x: fish age
y: fish length 

data = {{14, 590}, {21, 910}, {28, 1305}, {35, 1730}, {42, 2140}, {49, 2725},
	{56, 2890}, {63, 3685}, {70, 3920}, {77, 4325}, {84, 4410}, {91, 4485},
	{98, 4515}, {105, 4480}, {112, 4520}, {119, 4545}, {126, 4525},
	{133, 4560}, {140, 4565}, {147, 4626}, {154, 4566}};

data1 = Map[({age, length} = #;
	{age, age^2, Max[age-80, 0], Max[age-80, 0]^2, length})&, data]

spline regression: known knot (= 80)

If the knot is known, then the model is linear in the parameters...

Regress[data1, {1, x1, x2, x3, x4}, {x1, x2, x3, x4}]

{ParameterTable -> 
 
>          Estimate     SE          TStat      PValue      , 
      1    -320.994     130.734     -2.45532   0.0258957

                                                       -8
      x1   59.6324      6.13763     9.71586    4.096 10

      x2   0.00832024   0.0626651   0.132773   0.896028

                                                         -7
      x3   -61.4221     6.9967      -8.77873   1.62856 10

      x4   0.0167096    0.0623096   0.268171   0.791995
 
>    RSquared -> 0.997143, AdjustedRSquared -> 0.996428, 
 
>    EstimatedVariance -> 6868.4, 
 
>    ANOVATable ->         DF   SumOfSq       MeanSq        FRatio    PValue}
                                          7             6
                   Model   4    3.83521 10    9.58803 10    1395.96   0.

                   Error   16   109894.       6868.4

                                         7
                   Total   20   3.8462 10

spline regression: unknown knot (= k)

If the knot is unknown, then NonlinearFit/NonlinearRegress are required...

Derivative[0, 1][Max][0, x_] := x UnitStep[x]

NonlinearRegress[data, b0 + b1 x + b2 x^2 + b3 Max[x-k, 0] + b4 Max[x-k, 0]^2,
	x, {{b0, -321}, {b1, 60}, {b2, .0083}, {b3, -61}, {b4, .017}, {k, 80}}]
	
{BestFitParameters -> 
 
>     {b0 -> -317.91, b1 -> 59.4066, b2 -> 0.0112144, b3 -> -62.0161, 
 
>      b4 -> 0.0179745, k -> 80.0283}, 
 
>    ParameterCITable -> 
 
>          Estimate    Asymptotic SE   CI                   , 
      b0   -317.91     135.204         {-606.092, -29.7289}

      b1   59.4066     6.3458          {45.8808, 72.9323}

      b2   0.0112144   0.0647726       {-0.126845, 0.149274}

      b3   -62.0161    0.00189227      {-62.0201, -62.0121}

      b4   0.0179745   0.0654434       {-0.121515, 0.157464}

      k    80.0283     0.116694        {79.7796, 80.277}
 
>    EstimatedVariance -> 7351.58, 
 
>    ANOVATable ->                     DF   SumOfSq       MeanSq     , 
                                                      8             7
                   Model               6    2.99233 10    4.98722 10

                   Error               15   110274.       7351.58

                                                      8
                   Uncorrected Total   21   2.99344 10

                                                     7
                   Corrected Total     20   3.8462 10
 
>    AsymptoticCorrelationMatrix -> 
 
>     1.          -0.947631   0.883162    -0.591827   -0.736338   -0.580316, 

      -0.947631   1.          -0.982463   0.742232    0.779545    0.730698

      0.883162    -0.982463   1.          -0.835067   -0.738681   -0.824879

      -0.591827   0.742232    -0.835067   1.          0.288683    0.999814

      -0.736338   0.779545    -0.738681   0.288683    1.          0.270148

      -0.580316   0.730698    -0.824879   0.999814    0.270148    1.
 
>    FitCurvatureTable ->                           Curvature}
                          Max Intrinsic             0.0181284

                          Max Parameter-Effects     4.89117

                          95. % Confidence Region   0.587089



*)




