Go to the first, previous, next, last section, table of contents.


Fitting

This chapter describes functions for performing univariate least squares fits to experimental data. The data may be weighted or unweighted. 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.

Linear regression

The functions described in this section can be used to perform least-squares fits to a straight line model, Y = c_0 + c_1 X. For weighted data the best-fit is found by minimizing the weighted sum of squared residuals, \chi^2,

\chi^2 = \sum_i w_i (y_i - (c_0 + c_1 x_i))^2

for the parameters c_0, c_1. For unweighted data the sum is computed with w_i = 1.

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 datasets (x, y), two vectors of length n with strides xstride and ystride. 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.

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 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 covariance matrix for the parameters (c0, c1) is estimated from weighted data 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 c00, double c01, double c11, double *y, double *y_err)
This function uses the best-fit linear regression coefficients c0,c1 and their estimated 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. For weighted data the best-fit is found by minimizing the weighted sum of squared residuals, \chi^2,

\chi^2 = \sum_i w_i (y_i - c_1 x_i)^2

for the parameter c_1. For unweighted data the sum is computed with w_i = 1.

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 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 estimated from the weighted data and returned via the parameters 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 c11, double *y, double *y_err)
This function uses the best-fit linear regression coefficient c1 and its estimated 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 c are the p unknown best-fit parameters, which are to be estimated.

The best-fit is found by minimizing the weighted sums of squared residuals, \chi^2,

\chi^2 = (y - X c)^T W (y - X c)

with respect to the parameters c. The weights are specified by the diagonal elements of the n by n matrix W. For unweighted data W is replaced by the identity matrix.

This formulation can be used for fits to any number of functions and/or variables by preparing the matrix X appropriately. For example, to fit to a p-th order polynomial in x, use the following matrix

X_{ij} = x_i^j

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. The variance-covariance matrix of the model parameters cov is estimated from the scatter of the observations about the best-fit. The sum of squares of the residuals from the best-fit, \chi^2, is returned in chisq.

The best-fit is found by singular value decomposition of the matrix X using the preallocated workspace provided in work. 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 model y = X c for the observations y and the matrix of predictor variables X. The covariance matrix of the model parameters cov is estimated from the weighted data. The weighted sum of squares of the residuals from the best-fit, \chi^2, is returned in chisq.

The best-fit is found by singular value decomposition of the matrix X using the preallocated workspace provided in work. Any components which have zero singular value (to machine precision) are discarded from the fit.

Examples

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

#include <stdio.h>
#include <gsl/gsl_fit.h>

int
main ()
{
  int i, n = 4;
  double x[4] = { 1970, 1980, 1990, 2000 };
  double y[4] = {   12,   11,   14,   13 };
  double w[4] = {  0.1,  0.2,  0.3,  0.4 };

  double c0, c1, cov00, cov01, cov11, chisq;

  gsl_fit_wlinear (x, 1, w, 1, y, 1, n, 
                   &c0, &c1, &cov00, &cov01, &cov11, 
                   &chisq);

  printf("# best fit: Y = %g + %g X\n", c0, c1);
  printf("# covariance matrix:\n");
  printf("# [ %g, %g\n#   %g, %g]\n", cov00, cov01, cov01, cov11);
  printf("# chisq = %g\n", chisq);

  for (i = 0; i < n; i++)
    printf("data: %g %g %g\n", x[i], y[i], 1/sqrt(w[i]));

  printf("\n");

  for (i = -30 ; i < 130 ; i++)
    {
      double xf = x[0] + (i/100.0) * (x[n-1] - x[0]);
      double yf, yf_err;

      gsl_fit_linear_est (xf, c0, c1, cov00, cov01, cov11, &yf, &yf_err);

      printf("fit: %g %g\n", xf, yf);
      printf("hi : %g %g\n", xf, yf + yf_err);
      printf("lo : %g %g\n", xf, yf - yf_err);
    }
}

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

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.

The Review of Particle Physics is available online at http://pdg.lbl.gov/.

The tests used to prepare these routines are based on the NIST Statistical Reference Datasets. The datasets and their documentation are available from http://www.nist.gov/itl/div898/strd/index.html


Go to the first, previous, next, last section, table of contents.