Class for linear least-squares fit. More...

#include <LinearFit.h>

Inheritance diagram for casa::LinearFit< T >:

Public Member Functions
	LinearFit ()
	Create a fitter: the normal way to generate a fitter object.
	LinearFit (const LinearFit &other)
	Copy constructor (deep copy)
LinearFit &	operator= (const LinearFit &other)
	Assignment (deep copy)
virtual	~LinearFit ()
	Destructor.
Protected Member Functions
virtual Bool	fitIt (Vector< typename FunctionTraits< T >::BaseType > &sol, const Array< typename FunctionTraits< T >::BaseType > &x, const Vector< typename FunctionTraits< T >::BaseType > &y, const Vector< typename FunctionTraits< T >::BaseType > const sigma, const Vector< Bool > const mask=0)
	Generalised fitter.

Detailed Description

template<class T>
class casa::LinearFit< T >

Class for linear least-squares fit.

Review Status

Reviewed By:: wbrouw
Date Reviewed:: 2004/06/15
Test programs:: tLinearFitSVD

Prerequisite

Etymology

A set of data point is fit with some functional equation. The equations solved are linear equations. The functions themselves however can be wildly nonlinear.

Synopsis

NOTE: Constraints added. Documentation out of date at moment, check the tLinearFitSVD and tNonLinearFirLM programs for examples.

The following is a brief summary of the linear least-squares fit problem. See module header, Fitting , for a more complete description.

Given a set of N data points (measurements), (x(i), y(i)) i = 0,...,N-1, along with a set of standard deviations, sigma(i), for the data points, and M specified functions, f(j)(x) j = 0,...,M-1, we form a linear combination of the functions:

    z(i) = a(0)f(0)(x(i)) + a(1)f(1)(x(i)) + ... + a(M-1)f(M-1)(x(i)),

where a(j) j = 0,...,M-1 are a set of parameters to be determined. The linear least-squares fit tries to minimize

    chi-square = [(y(0)-z(0))/sigma(0)]^2 + [(y(1)-z(1))/sigma(1)]^2 + ..\.
                 + [(y(N-1)-z(N-1))/sigma(N-1)]^2.

by adjusting {a(j)} in the equation.

For complex numbers, [(y(i)-z(i))/sigma(i)]^2 in chi-square is replaced by [(y(i)-z(i))/sigma(i)]*conjugate([(y(i)-z(i))/sigma(i)])

For multidimensional functions, x(i) is a vector, and

    f(j)(x(i)) = f(j)(x(i,0), x(i,1), x(i,2), ...)

Normally, it is necessary that N > M for the solutions to be valid, since there must be more data points than model parameters to be solved.

If the measurement errors (standard deviation sigma) are not known at all, they can all be set to one initially. In this case, we assume all measurements have the same standard deviation, after minimizing chi-square, we recompute

    sigma^2 = {(y(0)-z(0))^2 + (y(1)-z(1))^2 + ..\.
              + (y(N-1)-z(N-1))^2}/(N-M) = chi-square/(N-M).

A statistic weight can also be assigned to each measurement if the standard deviation is not available. sigma can be calculated from

    sigma = 1/ sqrt(weight)

Alternatively a 'weight' switch can be set with asWeight(). For best arithmetic performance, weight should be normalized to a maximum value of one. Having a large weight value can sometimes lead to overflow problems.

The function to be fitted to the data can be given as an instance of the Function class. One can also form a sum of functions using the CompoundFunction .

For small datasets the usage of the calls is:

Create a functional description of the parameters
Create a fitter: LinearFit<T> fitter();
Set the functional representation: fitter.setFunction()
Do the fit to the data: fitter.fit(x, data, sigma) (or do a number of calls to buildNormalMatrix(x, data, sigma) and finish of with fitter.fit() or fitter.sol())
if needed the covariance; residuals; chiSquared, parameter errors can all be obtained

Note that the fitter is reusable. An example is given in the following.

The solution of a fit always produces the total number of parameters given to the fitter. I.e. including any parameters that were fixed. In the latter case the solution returned will be the fixed value.

Template Type Argument Requirements (T)

Float
Double
Complex
DComplex

If there are a large number of unknowns or a large number of data points machine memory limits (or timing reasons) may not allow a complete in-core fitting to be performed. In this case one can incrementally build the normal equation (see buildNormalMatrix()).

The normal operation of the class tests for real inversion problems only. If tests are needed for almost collinear columns in the solution matrix, the collinearity can be set as the square of the sine of the minimum angle allowed.

Singular Value Decomposition is supported by the LinearFitSVD class, which has a behaviour completely identical to this class (apart from a default collinearity of 1e-8).

Other information (see a.o. LSQFit ) can be set and obtained as well.

Motivation

The creation of this class was driven by the need to write code to perform baseline fitting or continuum subtraction.

Example

In the following a polynomial is fitted through the first 20 prime numbers. The data is given in the x vector (1 to 20) and in the primesTable (2, 3, ..., 71) (see tLinearFitSVD test program). In the following all four methods to calculate a polynomial through the data is used

        // The list of coordinate x-values
        Vector<Double> x(nPrimes);
        indgen((Array<Double>&)x, 1.0);  // 1, 2, ..\.
        Vector<Double> primesTable(nPrimes);
        for (uInt i=1; i < nPrimes; i++) {
           primesTable(i) =
           Primes::nextLargerPrimeThan(Int(primesTable(i-1)+0.01));
         };   
        Vector<Double> sigma(nPrimes);
        sigma = 1.0;
        // The fitter
        LinearFit<Double> fitter;
        Polynomial<AutoDiff<Double> > combination(2);
        // Get the solution
        fitter.setFunction(combination);
        Vector<Double> solution = fitter.fit(x, primesTable, sigma);
        // create a special function (should probably at beginning)
        static void myfnc(Vector<Double> &y, const Double x) {
        y(0) = 1; for (uInt i=1; i<y.nelements(); i++) y(i) = y(i-1)*x; };
        fitter.setFunction(3, &myfnc);
        solution = fitter.fit(x, primesTable, sigma);
        // Create the direct coefficients table
        fitter.setFunction(3);
        Matrix<Double> xx(nPrimes, 3);
        for (uInt i=0; i<nPrimes; i++) {
           xx(i,0) = 1;
           for (uInt j=1; j<3; j++) xx(i,j) = xx(i,j-1)*Double(i+1);
         };
         solution = fitter.fit(xx, primesTable, sigma);

In the test program examples are given on how to get the other information, and other examples.

Definition at line 207 of file LinearFit.h.

Constructor & Destructor Documentation

template<class T>

casa::LinearFit< T >::LinearFit ( )

Create a fitter: the normal way to generate a fitter object.

Necessary data will be deduced from the Functional provided with setFunction()

template<class T>

casa::LinearFit< T >::LinearFit ( const LinearFit< T > & other )

Copy constructor (deep copy)

template<class T>

virtual casa::LinearFit< T >::~LinearFit ( ) [virtual]

Destructor.

Member Function Documentation

template<class T>

virtual Bool casa::LinearFit< T >::fitIt	(	Vector< typename FunctionTraits< T >::BaseType > &	sol,
		const Array< typename FunctionTraits< T >::BaseType > &	x,
		const Vector< typename FunctionTraits< T >::BaseType > &	y,
		const Vector< typename FunctionTraits< T >::BaseType > *const	sigma,
		const Vector< Bool > *const	mask = `0`
	)		`[protected, virtual]`

Generalised fitter.

Implements casa::GenericL2Fit< T >.

template<class T>

LinearFit& casa::LinearFit< T >::operator= ( const LinearFit< T > & other )

Assignment (deep copy)

The documentation for this class was generated from the following file:

casacore/scimath/Fitting/LinearFit.h

Public Member Functions

Protected Member Functions

Detailed Description

template<class T> class casa::LinearFit< T >

Review Status

Prerequisite

Etymology

Synopsis

Template Type Argument Requirements (T)

Motivation

Example

Constructor & Destructor Documentation

Member Function Documentation

template<class T>
class casa::LinearFit< T >