fitting.texi

This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY without ev

@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::     
* 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 @math{\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 @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

@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

@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

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

@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}.
@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{c00}, double @var{c01}, double @var{c11}, 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{c11}, 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.

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