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

Overview
========

Least-squares fits are found by minimizing \chi^2 (chi-squared), the weighted sum of squared residuals over n experimental datapoints (x_i, y_i) for the model Y(c,x),

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

The p parameters of the model are $c = \c_0, c_1, \dots\$ c = {c_0, c_1, ...}.  The weight factors w_i are given by w_i = 1/\sigma_i^2, where \sigma_i is the experimental error on the data-point 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 c and their p \times p covariance matrix.  The covariance matrix measures the statistical errors on the best-fit parameters resulting from the errors on the data, \sigma_i, and is defined as $C_ab = \langle \delta c_a \delta c_b \rangle$ C_{ab} = <\delta c_a \delta c_b> where $\langle \, \rangle$ < > denotes an average over the Gaussian error distributions of the underlying datapoints.

The covariance matrix is calculated by error propagation from the data errors \sigma_i.  The change in a fitted parameter \delta c_a caused by a small change in the data \delta y_i is given by

     \delta c_a = \sum_i (dc_a/dy_i) \delta y_i

allowing the covariance matrix to be written in terms of the errors on the data,

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

For uncorrelated data the fluctuations of the underlying datapoints satisfy $\langle \delta y_i \delta y_j \rangle = \sigma_i^2 \delta_ij$ <\delta y_i \delta y_j> = \sigma_i^2 \delta_{ij}, giving a corresponding parameter covariance matrix of

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

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

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

Linear regression
=================

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

 -- Function: int gsl_fit_linear (const double * X, const size_t XSTRIDE, const double * Y, const size_t YSTRIDE, size_t N, double * C0, double * C1, double * COV00, double * COV01, double * COV11, double * SUMSQ)
     This function computes the best-fit linear regression coefficients (C0,C1) of the model Y = c_0 + c_1 X for the dataset (X, Y), two vectors of length N with strides XSTRIDE and YSTRIDE.  The errors on Y are assumed unknown so the variance-covariance matrix for the parameters (C0, C1) is estimated from the scatter of the points around the best-fit line and returned via the parameters (COV00, COV01, COV11). The sum of squares of the residuals from the best-fit line is returned in SUMSQ.  Note: the correlation coefficient of the data can be computed using gsl_stats_correlation (see Correlation), it does not depend on the fit.

 -- Function: int gsl_fit_wlinear (const double * X, const size_t XSTRIDE, const double * W, const size_t WSTRIDE, const double * Y, const size_t YSTRIDE, size_t N, double * C0, double * C1, double * COV00, double * COV01, double * COV11, double * CHISQ)
     This function computes the best-fit linear regression coefficients (C0,C1) of the model Y = c_0 + c_1 X for the weighted dataset (X, Y), two vectors of length N with strides XSTRIDE and YSTRIDE.  The vector W, of length N and stride WSTRIDE, specifies the weight of each datapoint. The weight is the reciprocal of the variance for each datapoint in Y.

     The covariance matrix for the parameters (C0, C1) is computed using the weights and returned via the parameters (COV00, COV01, COV11).  The weighted sum of squares of the residuals from the best-fit line, \chi^2, is returned in CHISQ.

 -- Function: int gsl_fit_linear_est (double X, double C0, double C1, double COV00, double COV01, double COV11, double * Y, double * Y_ERR)
     This function uses the best-fit linear regression coefficients C0, C1 and their covariance COV00, COV01, COV11 to compute the fitted function Y and its standard deviation Y_ERR for the model Y = c_0 + c_1 X at the point X.

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, Y = c_1 X.

 -- Function: int gsl_fit_mul (const double * X, const size_t XSTRIDE, const double * Y, const size_t YSTRIDE, size_t N, double * C1, double * COV11, double * SUMSQ)
     This function computes the best-fit linear regression coefficient C1 of the model Y = c_1 X for the datasets (X, Y), two vectors of length N with strides XSTRIDE and YSTRIDE.  The errors on Y are assumed unknown so the variance of the parameter C1 is estimated from the scatter of the points around the best-fit line and returned via the parameter COV11.  The sum of squares of the residuals from the best-fit line is returned in SUMSQ.

 -- Function: int gsl_fit_wmul (const double * X, const size_t XSTRIDE, const double * W, const size_t WSTRIDE, const double * Y, const size_t YSTRIDE, size_t N, double * C1, double * COV11, double * SUMSQ)
     This function computes the best-fit linear regression coefficient C1 of the model Y = c_1 X for the weighted datasets (X, Y), two vectors of length N with strides XSTRIDE and YSTRIDE.  The vector W, of length N and stride WSTRIDE, specifies the weight of each datapoint. The weight is the reciprocal of the variance for each datapoint in Y.

     The variance of the parameter C1 is computed using the weights and returned via the parameter COV11.  The weighted sum of squares of the residuals from the best-fit line, \chi^2, is returned in CHISQ.

 -- Function: int gsl_fit_mul_est (double X, double C1, double COV11, double * Y, double * Y_ERR)
     This function uses the best-fit linear regression coefficient C1 and its covariance COV11 to compute the fitted function Y and its standard deviation Y_ERR for the model Y = c_1 X at the point X.

Multi-parameter fitting
=======================

The functions described in this section perform least-squares fits to a general linear model, y = X c where y is a vector of n observations, X is an n by p matrix of predictor variables, and the elements of the vector c are the p unknown best-fit parameters which are to be estimated.  The chi-squared value is given by $\chi^2 = \sum_i w_i (y_i - \sum_j X_ij c_j)^2$ \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 n-by-p matrix X appropriately.  For example, to fit to a p-th order polynomial in X, use the following matrix,

     X_{ij} = x_i^j

where the index i runs over the observations and the index j runs from 0 to p-1.

To fit to a set of p sinusoidal functions with fixed frequencies \omega_1, \omega_2, ..., \omega_p, use,

     X_{ij} = sin(\omega_j x_i)

To fit to p independent variables x_1, x_2, ..., x_p, use,

     X_{ij} = x_j(i)

where x_j(i) is the i-th value of the predictor variable x_j.

The functions described in this section are declared in the header 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 X.

 -- Function: gsl_multifit_linear_workspace * gsl_multifit_linear_alloc (size_t N, size_t P)
     This function allocates a workspace for fitting a model to N observations using P parameters.

 -- Function: void gsl_multifit_linear_free (gsl_multifit_linear_workspace * WORK)
     This function frees the memory associated with the workspace W.

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

     The best-fit is found by singular value decomposition of the matrix 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.

 -- Function: int gsl_multifit_wlinear (const gsl_matrix * X, const gsl_vector * W, const gsl_vector * Y, gsl_vector * C, gsl_matrix * COV, double * CHISQ, gsl_multifit_linear_workspace * WORK)
     This function computes the best-fit parameters C of the weighted model y = X c for the observations Y with weights W and the matrix of predictor variables X, using the preallocated workspace provided in WORK.  The p-by-p covariance matrix of the model parameters COV is computed as (X^T W X)^-1. The weighted sum of squares of the residuals from the best-fit, \chi^2, is returned in CHISQ. If the coefficient of determination is desired, it can be computed from the expression R^2 = 1 - \chi^2 / WTSS, where the weighted total sum of squares (WTSS) of the observations Y may be computed from gsl_stats_wtss.

 -- Function: int gsl_multifit_linear_svd (const gsl_matrix * X, const gsl_vector * Y, double TOL, size_t * RANK, gsl_vector * C, gsl_matrix * COV, double * CHISQ, gsl_multifit_linear_workspace * WORK)
 -- Function: int gsl_multifit_wlinear_svd (const gsl_matrix * X, const gsl_vector * W, const gsl_vector * Y, double TOL, size_t * RANK, gsl_vector * C, gsl_matrix * COV, double * CHISQ, gsl_multifit_linear_workspace * WORK)
     In these functions components of the fit are discarded if the ratio of singular values s_i/s_0 falls below the user-specified tolerance TOL, and the effective rank is returned in RANK.

 -- Function: int gsl_multifit_linear_usvd (const gsl_matrix * X, const gsl_vector * Y, double TOL, size_t * RANK, gsl_vector * C, gsl_matrix * COV, double * CHISQ, gsl_multifit_linear_workspace * WORK)
 -- Function: int gsl_multifit_wlinear_usvd (const gsl_matrix * X, const gsl_vector * W, const gsl_vector * Y, double TOL, size_t * RANK, gsl_vector * C, gsl_matrix * COV, double * CHISQ, gsl_multifit_linear_workspace * WORK)
     These functions compute the fit using an SVD without column scaling.

 -- Function: int gsl_multifit_linear_est (const gsl_vector * X, const gsl_vector * C, const gsl_matrix * COV, double * Y, double * Y_ERR)
     This function uses the best-fit multilinear regression coefficients C and their covariance matrix COV to compute the fitted function value Y and its standard deviation Y_ERR for the model y = x.c at the point X.

 -- Function: int gsl_multifit_linear_residuals (const gsl_matrix * X, const gsl_vector * Y, const gsl_vector * C, gsl_vector * R)
     This function computes the vector of residuals r = y - X c for the observations Y, coefficients C and matrix of predictor variables X.

Robust linear 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

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

where e_i = y_i - Y(c, x_i) is the residual of the ith data point, and \rho(e_i) is a function which should have the following properties:

  \rho(e) \ge 0

  \rho(0) = 0

  \rho(-e) = \rho(e)

  \rho(e_1) > \rho(e_2) for |e_1| > |e_2|

The special case of ordinary least squares is given by \rho(e_i) = e_i^2. Letting \psi = \rho' be the derivative of \rho, differentiating the objective function with respect to the coefficients c and setting the partial derivatives to zero produces the system of equations

     \sum_i \psi(e_i) X_i = 0

where X_i is a vector containing row i of the design matrix X. Next, we define a weight function w(e) = \psi(e)/e, and let w_i = w(e_i):

     \sum_i w_i e_i X_i = 0

This system of equations is equivalent to solving a weighted ordinary least squares problem, minimizing \chi^2 = \sum_i w_i e_i^2. The weights however, depend on the residuals e_i, which depend on the coefficients c, which depend on the weights. Therefore, an iterative solution is used, called Iteratively Reweighted Least Squares (IRLS).

 1. Compute initial estimates of the coefficients c^(0) using ordinary least squares

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

 3. Compute new weights w_i^(k) = \psi(e_i^(k))/e_i^(k).

 4. Compute new coefficients c^(k) by solving the weighted least squares problem with weights w_i^(k).

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

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

     for all 0 \le i < p where \epsilon is a small tolerance factor.

The key to this method lies in selecting the function \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.

 -- Function: gsl_multifit_robust_workspace * gsl_multifit_robust_alloc (const gsl_multifit_robust_type * T, const size_t N, const size_t P)
     This function allocates a workspace for fitting a model to N observations using P parameters. The type T specifies the function \psi and can be selected from the following choices.
      -- Robust type: gsl_multifit_robust_default
          This specifies the gsl_multifit_robust_bisquare type (see below) and is a good general purpose choice for robust regression.

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

               w(e) = (1 - e^2)^2

          and the default tuning constant is t = 4.685.

      -- 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 C could minimize the objective function. Therefore this option should be used with care. The weight function is given by

               w(e) = 1 / (1 + e^2)

          and the default tuning constant is t = 2.385.

      -- Robust type: gsl_multifit_robust_fair
          This is the fair \rho function, which guarantees a unique solution and has continuous derivatives to three orders. The weight function is given by

               w(e) = 1 / (1 + |e|)

          and the default tuning constant is t = 1.400.

      -- Robust type: gsl_multifit_robust_huber
          This specifies Huber's \rho function, which is a parabola in the vicinity of zero and increases linearly for a given threshold |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 \psi function. The weight function is given by

               w(e) = 1/max(1,|e|)

          and the default tuning constant is t = 1.345.

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

               w(e) = 1

          and the default tuning constant is t = 1.

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

               w(e) = \exp(-e^2)

          and the default tuning constant is t = 2.985.

 -- Function: void gsl_multifit_robust_free (gsl_multifit_robust_workspace * W)
     This function frees the memory associated with the workspace W.

 -- Function: const char * gsl_multifit_robust_name (const gsl_multifit_robust_workspace * W)
     This function returns the name of the robust type T specified to gsl_multifit_robust_alloc.

 -- Function: int gsl_multifit_robust_tune (const double TUNE, gsl_multifit_robust_workspace * W)
     This function sets the tuning constant t used to adjust the residuals at each iteration to TUNE. Decreasing the tuning constant increases the downweight assigned to large residuals, while increasing the tuning constant decreases the downweight assigned to large residuals.

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

 -- Function: int gsl_multifit_robust_est (const gsl_vector * X, const gsl_vector * C, const gsl_matrix * COV, double * Y, double * Y_ERR)
     This function uses the best-fit robust regression coefficients C and their covariance matrix COV to compute the fitted function value Y and its standard deviation Y_ERR for the model y = x.c at the point X.

 -- Function: gsl_multifit_robust_stats gsl_multifit_robust_statistics (const gsl_multifit_robust_workspace * W)
     This function returns a structure containing relevant statistics from a robust regression. The function gsl_multifit_robust must be called first to perform the regression and calculate these statistics. The returned gsl_multifit_robust_stats structure contains the following fields.

       double sigma_ols This contains the standard deviation of the residuals as computed from ordinary least squares (OLS).

       double sigma_mad This contains an estimate of the standard deviation of the final residuals using the Median-Absolute-Deviation statistic

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

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

       double Rsq This contains the R^2 coefficient of determination statistic using the estimate sigma.

       double adj_Rsq This contains the adjusted R^2 coefficient of determination statistic using the estimate sigma.

       double rmse This contains the root mean squared error of the final residuals

       double sse This contains the residual sum of squares taking into account the robust covariance matrix.

       size_t dof This contains the number of degrees of freedom n - p

       size_t numit Upon successful convergence, this contains the number of iterations performed

       gsl_vector * weights This contains the final weight vector of length N

       gsl_vector * r This contains the final residual vector of length N, r = y - X c

Troubleshooting
===============

When using models based on polynomials, care should be taken when constructing the design matrix X. If the x values are large, then the matrix X could be ill-conditioned since its columns are powers of x, leading to unstable least-squares solutions. In this case it can often help to center and scale the x values using the mean and standard deviation:

     x' = (x - mu)/sigma

and then construct the X matrix using the transformed values x'.

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.

          [file fitting.c not included]

The following commands extract the data from the output of the program and display it using the GNU plotutils graph utility,

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

The next program performs a quadratic fit y = c_0 + c_1 x + c_2 x^2 to a weighted dataset using the generalised linear fitting function gsl_multifit_wlinear.  The model matrix X for a quadratic fit is given by,

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

where the column of ones corresponds to the constant term c_0. The two remaining columns corresponds to the terms c_1 x and c_2 x^2.

The program reads N lines of data in the format (X, Y, ERR) where ERR is the error (standard deviation) in the value Y.

          [file fitting2.c not included]

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 y = e^x in the region 0 < x < 2.

          [file fitting3.c not included]

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

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

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.

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

The parameters of the quadratic fit match the coefficients of the expansion of e^x, taking into account the errors on the parameters and the O(x^3) difference between the exponential and quadratic functions for the larger values of 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.

The next program demonstrates the advantage of robust least squares on a dataset with outliers. The program generates linear (x,y) data pairs on the line 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.

          [file robfit.c not included]

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

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,

  Review of Particle Properties, R.M. Barnett et al., Physical Review D54, 1 (1996) http://pdg.lbl.gov/

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

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,

http://www.nist.gov/itl/div898/strd/index.html.

The GSL implementation of robust linear regression closely follows the publications

  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

  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.