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:
For complex numbers,
For multidimensional functions, x(i) is a vector, and
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
A statistic weight can also be assigned to each measurement if the
standard deviation is not available. sigma can be calculated from
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:
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.
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
Other information (see a.o. LSQFit) can
be set and obtained as well.
Destructor
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.
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.
[(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)])
f(j)(x(i)) = f(j)(x(i,0), x(i,1), x(i,2), ...)
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).
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.
Note that the fitter is reusable. An example is given in the following.
Template Type Argument Requirements (T)
The following data types can be used to instantiate the LinearFit
templated class:
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.
Member Description
LinearFit()
Create a fitter: the normal way to generate a fitter object. Necessary
data will be deduced from the Functional provided with
setFunction()
LinearFit(const LinearFit &other)
Copy constructor (deep copy)
LinearFit &operator=(const LinearFit &other)
Assignment (deep copy)
virtual ~LinearFit()
virtual Bool fitIt (Vector<typename FunctionTraits<T>::BaseType> &sol, const T<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