Least-Squares Fitting ↑

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.

Contents:

  1. Overview

  2. Linear regression

    1. Module functions for linear regression

  3. Linear fitting without a constant term

  4. Multi-parameter fitting

    1. GSL::MultiFit::Workspace class

    2. Module functions

    3. Higer level interface

    4. NDLINEAR: multi-linear, multi-parameter least squares fitting (GSL extension)

  5. Examples

    1. Linear regression

    2. Exponential fitting

    3. Multi-parameter fitting

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), The p parameters of the model are 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} = <delta c_a delta c_b> where < > 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 allowing the covariance matrix to be written in terms of the errors on the data, For uncorrelated data the fluctuations of the underlying datapoints satisfy <delta y_i delta y_j> = sigma_i^2 delta_{ij}, giving a corresponding parameter covariance matrix of 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} = sqrt{C_{aa}}. The correlation coefficient of the fit parameters c_a and c_b is given by 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_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 - (c0 + c1 x_i))^2

for the parameters c0, c1. For unweighted data the sum is computed with w_i = 1.

Module functions for linear regression ↑

Linear fitting without a constant term ↑

Multi-parameter fitting ↑

Workspace class ↑

Module functions ↑

Higer level interface ↑

Examples ↑

Linear regression ↑

#!/usr/bin/env ruby
require("gsl")
include GSL::Fit

n = 4
x = Vector.alloc(1970, 1980, 1990, 2000)
y = Vector.alloc(12, 11, 14, 13)
w = Vector.alloc(0.1, 0.2, 0.3, 0.4)

#for i in 0...n do
#   printf("%e %e %e\n", x[i], y[i], 1.0/Math::sqrt(w[i]))
#end

c0, c1, cov00, cov01, cov11, chisq = wlinear(x, w, y)

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

Exponential fitting ↑

#!/usr/bin/env ruby
require("gsl")

# Create data
r = Rng.alloc("knuthran")
a = 2.0
b = -1.0
sigma = 0.01
N = 10
x = Vector.linspace(0, 5, N)
y = a*Sf::exp(b*x) + sigma*r.gaussian

# Fitting
a2, b2, = Fit.linear(x, Sf::log(y))
x2 = Vector.linspace(0, 5, 20)
A = Sf::exp(a2)
printf("Expect: a = %f, b = %f\n", a, b)
printf("Result: a = %f, b = %f\n", A, b2)
graph([x, y], [x2, A*Sf::exp(b2*x2)], "-C -g 3 -S 4")

Multi-parameter fitting ↑

#!/usr/bin/env ruby
require("gsl")
include GSL::MultiFit

Rng.env_setup()

r = GSL::Rng.alloc(Rng::DEFAULT)
n = 19
dim = 3
X = Matrix.alloc(n, dim)
y = Vector.alloc(n)
w = Vector.alloc(n)

a = 0.1
for i in 0...n
  y0 = Math::exp(a)
  sigma = 0.1*y0
  val = r.gaussian(sigma)
  X.set(i, 0, 1.0)
  X.set(i, 1, a)
  X.set(i, 2, a*a)
  y[i] = y0 + val
  w[i] = 1.0/(sigma*sigma)
  #printf("%g %g %g\n", a, y[i], sigma)
  a += 0.1
end

c, cov, chisq, status = MultiFit.wlinear(X, w, y)

prev next

Reference index top