@cindex fitting
@cindex least squares fit
@cindex regression, least squares
@cindex weighted linear fits
@cindex unweighted linear fits
This chapter describes routines for performing least squares fits to
experimental data using linear combinations of functions.  The data
may be weighted or unweighted, i.e. with known or unknown errors.  For
weighted data the functions compute the best fit parameters and their
associated covariance matrix.  For unweighted data the covariance
matrix is estimated from the scatter of the points, giving a
variance-covariance matrix.

The functions are divided into separate versions for simple one- or
two-parameter regression and multiple-parameter fits.  The functions
are declared in the header file @file{gsl_fit.h}.

@menu
* Fitting Overview::            
* Linear regression::           
* Linear fitting without a constant term::
* Multi-parameter fitting::
* Robust linear regression::
* Troubleshooting::
* Fitting Examples::            
* Fitting References and Further Reading::  
@end menu

@node Fitting Overview
@section Overview

Least-squares fits are found by minimizing @math{\chi^2}
(chi-squared), the weighted sum of squared residuals over @math{n}
experimental datapoints @math{(x_i, y_i)} for the model @math{Y(c,x)},
@tex
\beforedisplay
$$
\chi^2 = \sum_i w_i (y_i - Y(c, x_i))^2
$$
\afterdisplay
@end tex
@ifinfo

@example
\chi^2 = \sum_i w_i (y_i - Y(c, x_i))^2
@end example

@end ifinfo
@noindent
The @math{p} parameters of the model are @c{$c = \{c_0, c_1, \dots\}$}
@math{c = @{c_0, c_1, @dots{}@}}.  The
weight factors @math{w_i} are given by @math{w_i = 1/\sigma_i^2},
where @math{\sigma_i} is the experimental error on the data-point
@math{y_i}.  The errors are assumed to be
Gaussian and uncorrelated. 
For unweighted data the chi-squared sum is computed without any weight factors. 

The fitting routines return the best-fit parameters @math{c} and their
@math{p \times p} covariance matrix.  The covariance matrix measures the
statistical errors on the best-fit parameters resulting from the 
errors on the data, @math{\sigma_i}, and is defined
@cindex covariance matrix, linear fits
as @c{$C_{ab} = \langle \delta c_a \delta c_b \rangle$}
@math{C_@{ab@} = <\delta c_a \delta c_b>} where @c{$\langle \, \rangle$}
@math{< >} denotes an average over the Gaussian error distributions of the underlying datapoints.

The covariance matrix is calculated by error propagation from the data
errors @math{\sigma_i}.  The change in a fitted parameter @math{\delta
c_a} caused by a small change in the data @math{\delta y_i} is given
by
@tex
\beforedisplay
$$
\delta c_a = \sum_i {\partial c_a \over \partial y_i} \delta y_i
$$
\afterdisplay
@end tex
@ifinfo

@example
\delta c_a = \sum_i (dc_a/dy_i) \delta y_i
@end example

@end ifinfo
@noindent
allowing the covariance matrix to be written in terms of the errors on the data,
@tex
\beforedisplay
$$
C_{ab} =  \sum_{i,j} {\partial c_a \over \partial y_i}
                     {\partial c_b \over \partial y_j} 
                     \langle \delta y_i \delta y_j \rangle
$$
\afterdisplay
@end tex
@ifinfo

@example
C_@{ab@} = \sum_@{i,j@} (dc_a/dy_i) (dc_b/dy_j) <\delta y_i \delta y_j>
@end example

@end ifinfo
@noindent
For uncorrelated data the fluctuations of the underlying datapoints satisfy
@c{$\langle \delta y_i \delta y_j \rangle = \sigma_i^2 \delta_{ij}$}
@math{<\delta y_i \delta y_j> = \sigma_i^2 \delta_@{ij@}}, giving a 
corresponding parameter covariance matrix of
@tex
\beforedisplay
$$
C_{ab} = \sum_{i} {1 \over w_i} {\partial c_a \over \partial y_i} {\partial c_b \over \partial y_i} 
$$
\afterdisplay
@end tex
@ifinfo

@example
C_@{ab@} = \sum_i (1/w_i) (dc_a/dy_i) (dc_b/dy_i) 
@end example

@end ifinfo
@noindent
When computing the covariance matrix for unweighted data, i.e. data with unknown errors, 
the weight factors @math{w_i} in this sum are replaced by the single estimate @math{w =
1/\sigma^2}, where @math{\sigma^2} is the computed variance of the
residuals about the best-fit model, @math{\sigma^2 = \sum (y_i - Y(c,x_i))^2 / (n-p)}.  
This is referred to as the @dfn{variance-covariance matrix}.
@cindex variance-covariance matrix, linear fits

The standard deviations of the best-fit parameters are given by the
square root of the corresponding diagonal elements of
the covariance matrix, @c{$\sigma_{c_a} = \sqrt{C_{aa}}$}
@math{\sigma_@{c_a@} = \sqrt@{C_@{aa@}@}}.
The correlation coefficient of the fit parameters @math{c_a} and @math{c_b}
is given by @c{$\rho_{ab} = C_{ab} / \sqrt{C_{aa} C_{bb}}$}
@math{\rho_@{ab@} = C_@{ab@} / \sqrt@{C_@{aa@} C_@{bb@}@}}.

@node   Linear regression
@section Linear regression
@cindex linear regression

The functions described in this section can be used to perform
least-squares fits to a straight line model, @math{Y(c,x) = c_0 + c_1 x}.

@cindex covariance matrix, from linear regression
@deftypefun int gsl_fit_linear (const double * @var{x}, const size_t @var{xstride}, const double * @var{y}, const size_t @var{ystride}, size_t @var{n}, double * @var{c0}, double * @var{c1}, double * @var{cov00}, double * @var{cov01}, double * @var{cov11}, double * @var{sumsq})
This function computes the best-fit linear regression coefficients
(@var{c0},@var{c1}) of the model @math{Y = c_0 + c_1 X} for the dataset
(@var{x}, @var{y}), two vectors of length @var{n} with strides
@var{xstride} and @var{ystride}.  The errors on @var{y} are assumed unknown so 
the variance-covariance matrix for the
parameters (@var{c0}, @var{c1}) is estimated from the scatter of the
points around the best-fit line and returned via the parameters
(@var{cov00}, @var{cov01}, @var{cov11}).   
The sum of squares of the residuals from the best-fit line is returned
in @var{sumsq}.  Note: the correlation coefficient of the data can be computed using @code{gsl_stats_correlation} (@pxref{Correlation}), it does not depend on the fit.
@end deftypefun

@deftypefun int gsl_fit_wlinear (const double * @var{x}, const size_t @var{xstride}, const double * @var{w}, const size_t @var{wstride}, const double * @var{y}, const size_t @var{ystride}, size_t @var{n}, double * @var{c0}, double * @var{c1}, double * @var{cov00}, double * @var{cov01}, double * @var{cov11}, double * @var{chisq})
This function computes the best-fit linear regression coefficients
(@var{c0},@var{c1}) of the model @math{Y = c_0 + c_1 X} for the weighted
dataset (@var{x}, @var{y}), two vectors of length @var{n} with strides
@var{xstride} and @var{ystride}.  The vector @var{w}, of length @var{n}
and stride @var{wstride}, specifies the weight of each datapoint. The
weight is the reciprocal of the variance for each datapoint in @var{y}.

The covariance matrix for the parameters (@var{c0}, @var{c1}) is
computed using the weights and returned via the parameters
(@var{cov00}, @var{cov01}, @var{cov11}).  The weighted sum of squares
of the residuals from the best-fit line, @math{\chi^2}, is returned in
@var{chisq}.
@end deftypefun

@deftypefun int gsl_fit_linear_est (double @var{x}, double @var{c0}, double @var{c1}, double @var{cov00}, double @var{cov01}, double @var{cov11}, double * @var{y}, double * @var{y_err})
This function uses the best-fit linear regression coefficients
@var{c0}, @var{c1} and their covariance
@var{cov00}, @var{cov01}, @var{cov11} to compute the fitted function
@var{y} and its standard deviation @var{y_err} for the model @math{Y =
c_0 + c_1 X} at the point @var{x}.
@end deftypefun

@node Linear fitting without a constant term
@section Linear fitting without a constant term

The functions described in this section can be used to perform
least-squares fits to a straight line model without a constant term,
@math{Y = c_1 X}.

@deftypefun int gsl_fit_mul (const double * @var{x}, const size_t @var{xstride}, const double * @var{y}, const size_t @var{ystride}, size_t @var{n}, double * @var{c1}, double * @var{cov11}, double * @var{sumsq})
This function computes the best-fit linear regression coefficient
@var{c1} of the model @math{Y = c_1 X} for the datasets (@var{x},
@var{y}), two vectors of length @var{n} with strides @var{xstride} and
@var{ystride}.  The errors on @var{y} are assumed unknown so the 
variance of the parameter @var{c1} is estimated from
the scatter of the points around the best-fit line and returned via the
parameter @var{cov11}.  The sum of squares of the residuals from the
best-fit line is returned in @var{sumsq}.
@end deftypefun

@deftypefun int gsl_fit_wmul (const double * @var{x}, const size_t @var{xstride}, const double * @var{w}, const size_t @var{wstride}, const double * @var{y}, const size_t @var{ystride}, size_t @var{n}, double * @var{c1}, double * @var{cov11}, double * @var{sumsq})
This function computes the best-fit linear regression coefficient
@var{c1} of the model @math{Y = c_1 X} for the weighted datasets
(@var{x}, @var{y}), two vectors of length @var{n} with strides
@var{xstride} and @var{ystride}.  The vector @var{w}, of length @var{n}
and stride @var{wstride}, specifies the weight of each datapoint. The
weight is the reciprocal of the variance for each datapoint in @var{y}.

The variance of the parameter @var{c1} is computed using the weights
and returned via the parameter @var{cov11}.  The weighted sum of
squares of the residuals from the best-fit line, @math{\chi^2}, is
returned in @var{chisq}.
@end deftypefun

@deftypefun int gsl_fit_mul_est (double @var{x}, double @var{c1}, double @var{cov11}, double * @var{y}, double * @var{y_err})
This function uses the best-fit linear regression coefficient @var{c1}
and its covariance @var{cov11} to compute the fitted function
@var{y} and its standard deviation @var{y_err} for the model @math{Y =
c_1 X} at the point @var{x}.
@end deftypefun

@node Multi-parameter fitting
@section Multi-parameter fitting
@cindex multi-parameter regression
@cindex fits, multi-parameter linear
The functions described in this section perform least-squares fits to a
general linear model, @math{y = X c} where @math{y} is a vector of
@math{n} observations, @math{X} is an @math{n} by @math{p} matrix of
predictor variables, and the elements of the vector @math{c} are the @math{p} unknown best-fit parameters which are to be estimated.  The chi-squared value is given by @c{$\chi^2 = \sum_i w_i (y_i - \sum_j X_{ij} c_j)^2$}
@math{\chi^2 = \sum_i w_i (y_i - \sum_j X_@{ij@} c_j)^2}.

This formulation can be used for fits to any number of functions and/or
variables by preparing the @math{n}-by-@math{p} matrix @math{X}
appropriately.  For example, to fit to a @math{p}-th order polynomial in
@var{x}, use the following matrix,
@tex
\beforedisplay
$$
X_{ij} = x_i^j
$$
\afterdisplay
@end tex
@ifinfo

@example
X_@{ij@} = x_i^j
@end example

@end ifinfo
@noindent
where the index @math{i} runs over the observations and the index
@math{j} runs from 0 to @math{p-1}.

To fit to a set of @math{p} sinusoidal functions with fixed frequencies
@math{\omega_1}, @math{\omega_2}, @dots{}, @math{\omega_p}, use,
@tex
\beforedisplay
$$
X_{ij} = \sin(\omega_j x_i)
$$
\afterdisplay
@end tex
@ifinfo

@example
X_@{ij@} = sin(\omega_j x_i)
@end example

@end ifinfo
@noindent
To fit to @math{p} independent variables @math{x_1}, @math{x_2}, @dots{},
@math{x_p}, use,
@tex
\beforedisplay
$$
X_{ij} = x_j(i)
$$
\afterdisplay
@end tex
@ifinfo

@example
X_@{ij@} = x_j(i)
@end example

@end ifinfo
@noindent
where @math{x_j(i)} is the @math{i}-th value of the predictor variable
@math{x_j}.

The functions described in this section are declared in the header file
@file{gsl_multifit.h}.

The solution of the general linear least-squares system requires an
additional working space for intermediate results, such as the singular
value decomposition of the matrix @math{X}.

@deftypefun {gsl_multifit_linear_workspace *} gsl_multifit_linear_alloc (size_t @var{n}, size_t @var{p})
@tpindex gsl_multifit_linear_workspace
This function allocates a workspace for fitting a model to @var{n}
observations using @var{p} parameters.
@end deftypefun

@deftypefun void gsl_multifit_linear_free (gsl_multifit_linear_workspace * @var{work})
This function frees the memory associated with the workspace @var{w}.
@end deftypefun

@deftypefun int gsl_multifit_linear (const gsl_matrix * @var{X}, const gsl_vector * @var{y}, gsl_vector * @var{c}, gsl_matrix * @var{cov}, double * @var{chisq}, gsl_multifit_linear_workspace * @var{work})
This function computes the best-fit parameters @var{c} of the model
@math{y = X c} for the observations @var{y} and the matrix of
predictor variables @var{X}, using the preallocated workspace provided
in @var{work}.  The @math{p}-by-@math{p} variance-covariance matrix of the model parameters
@var{cov} is set to @math{\sigma^2 (X^T X)^{-1}}, where @math{\sigma} is
the standard deviation of the fit residuals.
The sum of squares of the residuals from the best-fit,
@math{\chi^2}, is returned in @var{chisq}. If the coefficient of
determination is desired, it can be computed from the expression
@math{R^2 = 1 - \chi^2 / TSS}, where the total sum of squares (TSS) of
the observations @var{y} may be computed from @code{gsl_stats_tss}.

The best-fit is found by singular value decomposition of the matrix
@var{X} using the modified Golub-Reinsch SVD algorithm, with column
scaling to improve the accuracy of the singular values. Any components
which have zero singular value (to machine precision) are discarded
from the fit.
@end deftypefun

@deftypefun int gsl_multifit_wlinear (const gsl_matrix * @var{X}, const gsl_vector * @var{w}, const gsl_vector * @var{y}, gsl_vector * @var{c}, gsl_matrix * @var{cov}, double * @var{chisq}, gsl_multifit_linear_workspace * @var{work})
This function computes the best-fit parameters @var{c} of the weighted
model @math{y = X c} for the observations @var{y} with weights @var{w}
and the matrix of predictor variables @var{X}, using the preallocated
workspace provided in @var{work}.  The @math{p}-by-@math{p} covariance matrix of the model
parameters @var{cov} is computed as @math{(X^T W X)^{-1}}. The weighted
sum of squares of the residuals from the best-fit, @math{\chi^2}, is
returned in @var{chisq}. If the coefficient of determination is
desired, it can be computed from the expression @math{R^2 = 1 - \chi^2
/ WTSS}, where the weighted total sum of squares (WTSS) of the
observations @var{y} may be computed from @code{gsl_stats_wtss}.
@end deftypefun

@deftypefun int gsl_multifit_linear_svd (const gsl_matrix * @var{X}, const gsl_vector * @var{y}, double @var{tol}, size_t * @var{rank}, gsl_vector * @var{c}, gsl_matrix * @var{cov}, double * @var{chisq}, gsl_multifit_linear_workspace * @var{work})
@deftypefunx int gsl_multifit_wlinear_svd (const gsl_matrix * @var{X}, const gsl_vector * @var{w}, const gsl_vector * @var{y}, double @var{tol}, size_t * @var{rank}, gsl_vector * @var{c}, gsl_matrix * @var{cov}, double * @var{chisq}, gsl_multifit_linear_workspace * @var{work})
In these functions components of the fit are discarded if the ratio of
singular values @math{s_i/s_0} falls below the user-specified
tolerance @var{tol}, and the effective rank is returned in @var{rank}.
@end deftypefun

@deftypefun int gsl_multifit_linear_usvd (const gsl_matrix * @var{X}, const gsl_vector * @var{y}, double @var{tol}, size_t * @var{rank}, gsl_vector * @var{c}, gsl_matrix * @var{cov}, double * @var{chisq}, gsl_multifit_linear_workspace * @var{work})
@deftypefunx int gsl_multifit_wlinear_usvd (const gsl_matrix * @var{X}, const gsl_vector * @var{w}, const gsl_vector * @var{y}, double @var{tol}, size_t * @var{rank}, gsl_vector * @var{c}, gsl_matrix * @var{cov}, double * @var{chisq}, gsl_multifit_linear_workspace * @var{work})
These functions compute the fit using an SVD without column scaling.
@end deftypefun

@deftypefun int gsl_multifit_linear_est (const gsl_vector * @var{x}, const gsl_vector * @var{c}, const gsl_matrix * @var{cov}, double * @var{y}, double * @var{y_err})
This function uses the best-fit multilinear regression coefficients
@var{c} and their covariance matrix
@var{cov} to compute the fitted function value
@var{y} and its standard deviation @var{y_err} for the model @math{y = x.c} 
at the point @var{x}.
@end deftypefun

@deftypefun int gsl_multifit_linear_residuals (const gsl_matrix * @var{X}, const gsl_vector * @var{y}, const gsl_vector * @var{c}, gsl_vector * @var{r})
This function computes the vector of residuals @math{r = y - X c} for
the observations @var{y}, coefficients @var{c} and matrix of predictor
variables @var{X}.
@end deftypefun

@node Robust linear regression
@section Robust linear regression
@cindex robust regression

Ordinary least squares (OLS) models are often heavily influenced by the presence of outliers.
Outliers are data points which do not follow the general trend of the other observations,
although there is strictly no precise definition of an outlier. Robust linear regression
refers to regression algorithms which are robust to outliers. The most common type of
robust regression is M-estimation. The general M-estimator minimizes the objective function
@tex
\beforedisplay
$$
\sum_i \rho(e_i) = \sum_i \rho (y_i - Y(c, x_i))
$$
\afterdisplay
@end tex
@ifinfo

@example
\sum_i \rho(e_i) = \sum_i \rho (y_i - Y(c, x_i))
@end example

@end ifinfo
where @math{e_i = y_i - Y(c, x_i)} is the residual of the ith data point, and
@math{\rho(e_i)} is a function which should have the following properties:
@itemize @w{}
@item @math{\rho(e) \ge 0}
@item @math{\rho(0) = 0}
@item @math{\rho(-e) = \rho(e)}
@item @math{\rho(e_1) > \rho(e_2)} for @math{|e_1| > |e_2|}
@end itemize
@noindent
The special case of ordinary least squares is given by @math{\rho(e_i) = e_i^2}.
Letting @math{\psi = \rho'} be the derivative of @math{\rho}, differentiating
the objective function with respect to the coefficients @math{c}
and setting the partial derivatives to zero produces the system of equations
@tex
\beforedisplay
$$
\sum_i \psi(e_i) X_i = 0
$$
\afterdisplay
@end tex
@ifinfo

@example
\sum_i \psi(e_i) X_i = 0
@end example

@end ifinfo
where @math{X_i} is a vector containing row @math{i} of the design matrix @math{X}.
Next, we define a weight function @math{w(e) = \psi(e)/e}, and let
@math{w_i = w(e_i)}:
@tex
\beforedisplay
$$
\sum_i w_i e_i X_i = 0
$$
\afterdisplay
@end tex
@ifinfo

@example
\sum_i w_i e_i X_i = 0
@end example

@end ifinfo
This system of equations is equivalent to solving a weighted ordinary least squares
problem, minimizing @math{\chi^2 = \sum_i w_i e_i^2}. The weights however, depend
on the residuals @math{e_i}, which depend on the coefficients @math{c}, which depend
on the weights. Therefore, an iterative solution is used, called Iteratively Reweighted
Least Squares (IRLS).
@enumerate
@item Compute initial estimates of the coefficients @math{c^{(0)}} using ordinary least squares

@item For iteration @math{k}, form the residuals @math{e_i^{(k)} = (y_i - X_i c^{(k-1)})/(t \sigma^{(k)} \sqrt{1 - h_i})},
where @math{t} is a tuning constant depending on the choice of @math{\psi}, and @math{h_i} are the
statistical leverages (diagonal elements of the matrix @math{X (X^T X)^{-1} X^T}). Including @math{t}
and @math{h_i} in the residual calculation has been shown to improve the convergence of the method.
The residual standard deviation is approximated as @math{\sigma^{(k)} = MAD / 0.6745}, where MAD is the
Median-Absolute-Deviation of the @math{n-p} largest residuals from the previous iteration.

@item Compute new weights @math{w_i^{(k)} = \psi(e_i^{(k)})/e_i^{(k)}}.

@item Compute new coefficients @math{c^{(k)}} by solving the weighted least squares problem with
weights @math{w_i^{(k)}}.

@item Steps 2 through 4 are iterated until the coefficients converge or until some maximum iteration
limit is reached. Coefficients are tested for convergence using the critera:
@tex
\beforedisplay
$$
|c_i^{(k)} - c_i^{(k-1)}| \le \epsilon \times \hbox{max}(|c_i^{(k)}|, |c_i^{(k-1)}|)
$$
\afterdisplay
@end tex
@ifinfo

@example
|c_i^(k) - c_i^(k-1)| \le \epsilon \times max(|c_i^(k)|, |c_i^(k-1)|)
@end example

@end ifinfo
for all @math{0 \le i < p} where @math{\epsilon} is a small tolerance factor.
@end enumerate
@noindent
The key to this method lies in selecting the function @math{\psi(e_i)} to assign
smaller weights to large residuals, and larger weights to smaller residuals. As
the iteration proceeds, outliers are assigned smaller and smaller weights, eventually
having very little or no effect on the fitted model.

@deftypefun {gsl_multifit_robust_workspace *} gsl_multifit_robust_alloc (const gsl_multifit_robust_type * @var{T}, const size_t @var{n}, const size_t @var{p})
@tpindex gsl_multifit_robust_workspace
This function allocates a workspace for fitting a model to @var{n}
observations using @var{p} parameters. The type @var{T} specifies the
function @math{\psi} and can be selected from the following choices.
@deffn {Robust type} gsl_multifit_robust_default
This specifies the @code{gsl_multifit_robust_bisquare} type (see below) and is a good
general purpose choice for robust regression.
@end deffn

@deffn {Robust type} gsl_multifit_robust_bisquare
This is Tukey's biweight (bisquare) function and is a good general purpose choice for
robust regression. The weight function is given by
@tex
\beforedisplay
$$
w(e) = \left\{ \matrix{ (1 - e^2)^2, & |e| \le 1\cr
                        0, & |e| > 1} \right.
$$
\afterdisplay
@end tex
@ifinfo

@example
w(e) = (1 - e^2)^2
@end example

@end ifinfo
and the default tuning constant is @math{t = 4.685}.
@end deffn

@deffn {Robust type} gsl_multifit_robust_cauchy
This is Cauchy's function, also known as the Lorentzian function.
This function does not guarantee a unique solution,
meaning different choices of the coefficient vector @var{c}
could minimize the objective function. Therefore this option should
be used with care. The weight function is given by
@tex
\beforedisplay
$$
w(e) = {1 \over 1 + e^2}
$$
\afterdisplay
@end tex
@ifinfo

@example
w(e) = 1 / (1 + e^2)
@end example

@end ifinfo
and the default tuning constant is @math{t = 2.385}.
@end deffn

@deffn {Robust type} gsl_multifit_robust_fair
This is the fair @math{\rho} function, which guarantees a unique solution and
has continuous derivatives to three orders. The weight function is given by
@tex
\beforedisplay
$$
w(e) = {1 \over 1 + |e|}
$$
\afterdisplay
@end tex
@ifinfo

@example
w(e) = 1 / (1 + |e|)
@end example

@end ifinfo
and the default tuning constant is @math{t = 1.400}.
@end deffn

@deffn {Robust type} gsl_multifit_robust_huber
This specifies Huber's @math{\rho} function, which is a parabola in the vicinity of zero and
increases linearly for a given threshold @math{|e| > t}. This function is also considered
an excellent general purpose robust estimator, however, occasional difficulties can
be encountered due to the discontinuous first derivative of the @math{\psi} function.
The weight function is given by
@tex
\beforedisplay
$$
w(e) = \left\{ \matrix{ 1, & |e| \le 1\cr
                        {1 \over |e|}, & |e| > 1} \right.
$$
\afterdisplay
@end tex
@ifinfo

@example
w(e) = 1/max(1,|e|)
@end example

@end ifinfo
and the default tuning constant is @math{t = 1.345}.
@end deffn

@deffn {Robust type} gsl_multifit_robust_ols
This specifies the ordinary least squares solution, which can be useful for quickly
checking the difference between the various robust and OLS solutions. The weight function is given by
@tex
\beforedisplay
$$
w(e) = 1
$$
\afterdisplay
@end tex
@ifinfo

@example
w(e) = 1
@end example

@end ifinfo
and the default tuning constant is @math{t = 1}.
@end deffn

@deffn {Robust type} gsl_multifit_robust_welsch
This specifies the Welsch function which can perform well in cases where the residuals have an
exponential distribution. The weight function is given by
@tex
\beforedisplay
$$
w(e) = \exp{(-e^2)}
$$
\afterdisplay
@end tex
@ifinfo

@example
w(e) = \exp(-e^2)
@end example

@end ifinfo
and the default tuning constant is @math{t = 2.985}.
@end deffn
@end deftypefun

@deftypefun void gsl_multifit_robust_free (gsl_multifit_robust_workspace * @var{w})
This function frees the memory associated with the workspace @var{w}.
@end deftypefun

@deftypefun {const char *} gsl_multifit_robust_name (const gsl_multifit_robust_workspace * @var{w})
This function returns the name of the robust type @var{T} specified to @code{gsl_multifit_robust_alloc}.
@end deftypefun

@deftypefun int gsl_multifit_robust_tune (const double @var{tune}, gsl_multifit_robust_workspace * @var{w})
This function sets the tuning constant @math{t} used to adjust the residuals at each iteration to @var{tune}.
Decreasing the tuning constant increases the downweight assigned to large residuals, while increasing
the tuning constant decreases the downweight assigned to large residuals.
@end deftypefun

@deftypefun int gsl_multifit_robust (const gsl_matrix * @var{X}, const gsl_vector * @var{y}, gsl_vector * @var{c}, gsl_matrix * @var{cov}, gsl_multifit_robust_workspace * @var{w})
This function computes the best-fit parameters @var{c} of the model
@math{y = X c} for the observations @var{y} and the matrix of
predictor variables @var{X}, attemping to reduce the influence
of outliers using the algorithm outlined above.
The @math{p}-by-@math{p} variance-covariance matrix of the model parameters
@var{cov} is estimated as @math{\sigma^2 (X^T X)^{-1}}, where @math{\sigma} is
an approximation of the residual standard deviation using the theory of robust
regression. Special care must be taken when estimating @math{\sigma} and
other statistics such as @math{R^2}, and so these
are computed internally and are available by calling the function
@code{gsl_multifit_robust_statistics}.
@end deftypefun

@deftypefun int gsl_multifit_robust_est (const gsl_vector * @var{x}, const gsl_vector * @var{c}, const gsl_matrix * @var{cov}, double * @var{y}, double * @var{y_err})
This function uses the best-fit robust regression coefficients
@var{c} and their covariance matrix
@var{cov} to compute the fitted function value
@var{y} and its standard deviation @var{y_err} for the model @math{y = x.c} 
at the point @var{x}.
@end deftypefun

@deftypefun gsl_multifit_robust_stats gsl_multifit_robust_statistics (const gsl_multifit_robust_workspace * @var{w})
This function returns a structure containing relevant statistics from a robust regression. The function
@code{gsl_multifit_robust} must be called first to perform the regression and calculate these statistics.
The returned @code{gsl_multifit_robust_stats} structure contains the following fields.
@itemize @w{}
@item double @code{sigma_ols} This contains the standard deviation of the residuals as computed from ordinary least squares (OLS).

@item double @code{sigma_mad} This contains an estimate of the standard deviation of the final residuals using the Median-Absolute-Deviation statistic

@item double @code{sigma_rob} This contains an estimate of the standard deviation of the final residuals from the theory of robust regression (see Street et al, 1988).

@item double @code{sigma} This contains an estimate of the standard deviation of the final residuals by attemping to reconcile @code{sigma_rob} and @code{sigma_ols}
in a reasonable way.

@item double @code{Rsq} This contains the @math{R^2} coefficient of determination statistic using the estimate @code{sigma}.

@item double @code{adj_Rsq} This contains the adjusted @math{R^2} coefficient of determination statistic using the estimate @code{sigma}.

@item double @code{rmse} This contains the root mean squared error of the final residuals

@item double @code{sse} This contains the residual sum of squares taking into account the robust covariance matrix.

@item size_t @code{dof} This contains the number of degrees of freedom @math{n - p}

@item size_t @code{numit} Upon successful convergence, this contains the number of iterations performed

@item gsl_vector * @code{weights} This contains the final weight vector of length @var{n}

@item gsl_vector * @code{r} This contains the final residual vector of length @var{n}, @math{r = y - X c}
@end itemize
@end deftypefun

@node Troubleshooting
@section Troubleshooting
@cindex least squares troubleshooting

When using models based on polynomials, care should be taken when constructing the design matrix
@math{X}. If the @math{x} values are large, then the matrix @math{X} could be ill-conditioned
since its columns are powers of @math{x}, leading to unstable least-squares solutions.
In this case it can often help to center and scale the @math{x} values using the mean and standard deviation:
@tex
\beforedisplay
$$
x' = {x - \mu(x) \over \sigma(x)}
$$
\afterdisplay
@end tex
@ifinfo

@example
x' = (x - mu)/sigma
@end example

@end ifinfo
@noindent
and then construct the @math{X} matrix using the transformed values @math{x'}.

@node Fitting Examples
@section Examples

The following program computes a least squares straight-line fit to a
simple dataset, and outputs the best-fit line and its
associated one standard-deviation error bars.

@example
@verbatiminclude examples/fitting.c
@end example

@noindent
The following commands extract the data from the output of the program
and display it using the @sc{gnu} plotutils @code{graph} utility, 

@example
$ ./demo > tmp
$ more tmp
# best fit: Y = -106.6 + 0.06 X
# covariance matrix:
# [ 39602, -19.9
#   -19.9, 0.01]
# chisq = 0.8

$ for n in data fit hi lo ; 
   do 
     grep "^$n" tmp | cut -d: -f2 > $n ; 
   done
$ graph -T X -X x -Y y -y 0 20 -m 0 -S 2 -Ie data 
     -S 0 -I a -m 1 fit -m 2 hi -m 2 lo
@end example

@iftex
@sp 1
@center @image{fit-wlinear,3.0in}
@end iftex

The next program performs a quadratic fit @math{y = c_0 + c_1 x + c_2
x^2} to a weighted dataset using the generalised linear fitting function
@code{gsl_multifit_wlinear}.  The model matrix @math{X} for a quadratic
fit is given by,
@tex
\beforedisplay
$$
X=\pmatrix{1&x_0&x_0^2\cr
1&x_1&x_1^2\cr
1&x_2&x_2^2\cr
\dots&\dots&\dots\cr}
$$
\afterdisplay
@end tex
@ifinfo

@example
X = [ 1   , x_0  , x_0^2 ;
      1   , x_1  , x_1^2 ;
      1   , x_2  , x_2^2 ;
      ... , ...  , ...   ]
@end example

@end ifinfo
@noindent
where the column of ones corresponds to the constant term @math{c_0}.
The two remaining columns corresponds to the terms @math{c_1 x} and
@math{c_2 x^2}.

The program reads @var{n} lines of data in the format (@var{x}, @var{y},
@var{err}) where @var{err} is the error (standard deviation) in the
value @var{y}.

@example
@verbatiminclude examples/fitting2.c
@end example

@noindent
A suitable set of data for fitting can be generated using the following
program.  It outputs a set of points with gaussian errors from the curve
@math{y = e^x} in the region @math{0 < x < 2}.

@example
@verbatiminclude examples/fitting3.c
@end example

@noindent
The data can be prepared by running the resulting executable program,

@example
$ GSL_RNG_TYPE=mt19937_1999 ./generate > exp.dat
$ more exp.dat
0.1 0.97935 0.110517
0.2 1.3359 0.12214
0.3 1.52573 0.134986
0.4 1.60318 0.149182
0.5 1.81731 0.164872
0.6 1.92475 0.182212
....
@end example

@noindent
To fit the data use the previous program, with the number of data points
given as the first argument.  In this case there are 19 data points.

@example
$ ./fit 19 < exp.dat
0.1 0.97935 +/- 0.110517
0.2 1.3359 +/- 0.12214
...
# best fit: Y = 1.02318 + 0.956201 X + 0.876796 X^2
# covariance matrix:
[ +1.25612e-02, -3.64387e-02, +1.94389e-02  
  -3.64387e-02, +1.42339e-01, -8.48761e-02  
  +1.94389e-02, -8.48761e-02, +5.60243e-02 ]
# chisq = 23.0987
@end example

@noindent
The parameters of the quadratic fit match the coefficients of the
expansion of @math{e^x}, taking into account the errors on the
parameters and the @math{O(x^3)} difference between the exponential and
quadratic functions for the larger values of @math{x}.  The errors on
the parameters are given by the square-root of the corresponding
diagonal elements of the covariance matrix.  The chi-squared per degree
of freedom is 1.4, indicating a reasonable fit to the data.

@iftex
@sp 1
@center @image{fit-wlinear2,3.0in}
@end iftex

The next program demonstrates the advantage of robust least squares on
a dataset with outliers. The program generates linear @math{(x,y)}
data pairs on the line @math{y = 1.45 x + 3.88}, adds some random
noise, and inserts 3 outliers into the dataset. Both the robust
and ordinary least squares (OLS) coefficients are computed for
comparison.

@example
@verbatiminclude examples/robfit.c
@end example

The output from the program is shown in the following plot.

@iftex
@sp 1
@center @image{robust,3.0in}
@end iftex

@node Fitting References and Further Reading
@section References and Further Reading

A summary of formulas and techniques for least squares fitting can be
found in the ``Statistics'' chapter of the Annual Review of Particle
Physics prepared by the Particle Data Group,

@itemize @w{}
@item
@cite{Review of Particle Properties},
R.M. Barnett et al., Physical Review D54, 1 (1996)
@uref{http://pdg.lbl.gov/}
@end itemize

@noindent
The Review of Particle Physics is available online at the website given
above.

@cindex NIST Statistical Reference Datasets
@cindex Statistical Reference Datasets (StRD)
The tests used to prepare these routines are based on the NIST
Statistical Reference Datasets. The datasets and their documentation are
available from NIST at the following website,

@center @uref{http://www.nist.gov/itl/div898/strd/index.html}.

@noindent
The GSL implementation of robust linear regression closely follows the publications

@itemize @w{}
@item DuMouchel, W. and F. O'Brien (1989), "Integrating a robust
option into a multiple regression computing environment,"
Computer Science and Statistics:  Proceedings of the 21st
Symposium on the Interface, American Statistical Association

@item Street, J.O., R.J. Carroll, and D. Ruppert (1988), "A note on
computing robust regression estimates via iteratively
reweighted least squares," The American Statistician, v. 42, 
pp. 152-154.
@end itemize