The class acts as a base class for L2-type (least-squares) fitting.
Actual classes (se e.g. LinearFit and
NonLinearFit.
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.
Create a fitter: the normal way to generate a fitter object. Necessary
data will be deduced from the Functional provided with
setFunction()
Destructor
Sets the function to be fitted. Upon entry, the argument function object
is cloned. The cloned copy is used in the later fitting process.
A valid function should be an instance of the
Function class,
so that derivatives with respect to the adjustable parameters
can be calculated. The current values of the "available" parameters
of the function are taken as the initial guess for the non-linear fitting.
Set sigma values to be interpreted as weight (i.e. 1/sigma/sigma).
A value of zero or -1 will be skipped. The switch will stay in effect
until set False again explicitly. Default is False.
Return a pointer to the function being fitted. Should
never delete this pointer.
Return the nth constraint equation derived from SVD
Note that the number present will be given by getDeficiency()
Fit the function to the data. If no sigma provided, all ones assumed.
In the case of no x,y,sigma the fitting equations are supposed to be
generated by previous calls to buildNormalMatrix. Note that the ones
with a scalar sigma will assume sigma=1 (overloading problem). The mask
assumes that if present, points with False will be skipped.
Obtain the chi squared. It has already been calculated during the
fitting process.
Get the errors on the solved values
Get covariance matrix
Generate the normal equations by one or more calls to the
buildNormalMatrix(), before calling a fit() without arguments.
The arguments are the same as for the fit(arguments) function.
A False is returned if the Array sizes are unmatched.
Prerequisite
Etymology
A set of data point is fit with some functional equation.
The class acts as a generic base class for L2 type
fits.
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 GenericL2Fit
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 tGenericL2FitSVD 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(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;
// Linear combination of functions describing 1 + x + x*x
combination.setCoefficient(0, 1.0); // 1
combination.setCoefficient(1, 1.0); // x
combination.setCoefficient(2, 1.0); // x^2
// Get the solution
fitter.setFunction(combination);
Vector<Double> solution = fitter.fit(x, primesTable, sigma);
// Try with a function with automatic derivatives (note that default
// polynomial has zero first guess)
LinearFit<AutoDiffA<Double> > fitad;
Polynomial<AutoDiffA<Double> > sqre(2);
fitad.setFunction(sqre);
solution = fitad.fit(x, primesTable, sigma);
In the test program examples are given on how to get the other
information, and other examples.
Member Description
GenericL2Fit()
GenericL2Fit(const GenericL2Fit &other)
Copy constructor (deep copy)
GenericL2Fit &operator=(const GenericL2Fit &other)
Assignment (deep copy)
virtual ~GenericL2Fit()
template <class U> void setFunction(const U<U,U> &function)
template <class U> Bool setConstraint(const uInt n, const Function<U,U> &function, const Vector<typename FunctionTraits<T>::BaseType> &x, const typename FunctionTraits<T>::BaseType y= typename FunctionTraits<T>::BaseType(0))
Set the possible constraint functions. The addConstraint
will add one; the setConstraint will [re-]set the
nth constraint. If unsucessful, False returned.
Bool setConstraint(const uInt n, const Vector<typename FunctionTraits<T>::BaseType> &x, const typename FunctionTraits<T>::BaseType y= typename FunctionTraits<T>::BaseType(0))
Bool setConstraint(const uInt n, const typename FunctionTraits<T>::BaseType y= typename FunctionTraits<T>::BaseType(0))
Bool addConstraint(const Function<typename FunctionTraits<T>::DiffType, typename FunctionTraits<T>::DiffType> &function, const Vector<typename FunctionTraits<T>::BaseType> &x, const typename FunctionTraits<T>::BaseType y= typename FunctionTraits<T>::BaseType(0))
Bool addConstraint(const Vector<typename FunctionTraits<T>::BaseType> &x, const typename FunctionTraits<T>::BaseType y= typename FunctionTraits<T>::BaseType(0))
Bool addConstraint(const typename FunctionTraits<T>::BaseType y= typename FunctionTraits<T>::BaseType(0))
Constraint functional can only be set when the function to be fitted
has been set. It should have the same number of parameters as the function
to be fitted. The x should have the correct dimension.
void setCollinearity(const Double cln)
Set the collinearity factor as the square of the sine of the
minimum angle allowed between input vectors (default zero for non-SVD,
1e-8 for SVD)
void asWeight(const Bool aswgt)
void asSVD(const Bool svd)
Set the use of SVD or not (default). When set the default collinearity
is set as well.
Function<typename FunctionTraits<T>::DiffType, typename FunctionTraits<T>::DiffType> *fittedFunction()
const Function<typename FunctionTraits<T>::DiffType, typename FunctionTraits<T>::DiffType> *const fittedFunction() const
uInt fittedNumber() const
Return the number of fitted parameters
uInt NConstraints()
Return the number of constraints, and pointers to constraint functions.
A 0-pointer will be returned if no such constraint present.
This pointer should never be destroyed.
Function<typename FunctionTraits<T>::DiffType, typename FunctionTraits<T>::DiffType> *getConstraint(const uInt n)
Vector<typename LSQTraits<typename LSQTraits<T>:: BaseType>::base> getSVDConstraint(uInt n)
void setParameterValues (const Vector<typename FunctionTraits<T>::BaseType> &parms)
Set the parameter values. The input is a vector of parameters; all
or only the masked ones' values will be set, using the input values
void setMaskedParameterValues (const Vector<typename FunctionTraits<T>::BaseType> &parms)
Vector<typename FunctionTraits<T>::BaseType> fit(const Vector<typename FunctionTraits<T>::BaseType> &x, const Vector<typename FunctionTraits<T>::BaseType> &y, const Vector<typename FunctionTraits<T>::BaseType> &sigma, const Vector<Bool> *const mask=0)
Vector<typename FunctionTraits<T>::BaseType> fit(const T<typename FunctionTraits<T>::BaseType> &x, const Vector<typename FunctionTraits<T>::BaseType> &y, const Vector<typename FunctionTraits<T>::BaseType> &sigma, const Vector<Bool> *const mask=0)
Vector<typename FunctionTraits<T>::BaseType> fit(const Vector<typename FunctionTraits<T>::BaseType> &x, const Vector<typename FunctionTraits<T>::BaseType> &y, const Vector<Bool> *const mask=0)
Vector<typename FunctionTraits<T>::BaseType> fit(const T<typename FunctionTraits<T>::BaseType> &x, const Vector<typename FunctionTraits<T>::BaseType> &y, const Vector<Bool> *const mask=0)
Vector<typename FunctionTraits<T>::BaseType> fit(const Vector<Bool> *const mask=0)
Bool fit(Vector<typename FunctionTraits<T>::BaseType> &sol, const Vector<typename FunctionTraits<T>::BaseType> &x, const Vector<typename FunctionTraits<T>::BaseType> &y, const Vector<typename FunctionTraits<T>::BaseType> &sigma, const Vector<Bool> *const mask=0)
Bool fit(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> &sigma, const Vector<Bool> *const mask=0)
Bool fit(Vector<typename FunctionTraits<T>::BaseType> &sol, const Vector<typename FunctionTraits<T>::BaseType> &x, const Vector<typename FunctionTraits<T>::BaseType> &y, const typename FunctionTraits<T>::BaseType &sigma, const Vector<Bool> *const mask=0)
Bool fit(Vector<typename FunctionTraits<T>::BaseType> &sol, const T<typename FunctionTraits<T>::BaseType> &x, const Vector<typename FunctionTraits<T>::BaseType> &y, const typename FunctionTraits<T>::BaseType &sigma, const Vector<Bool> *const mask=0)
Bool fit(Vector<typename FunctionTraits<T>::BaseType> &sol, const Vector<Bool> *const mask=0)
Thrown Exceptions
Double chiSquare() const
const Vector<typename FunctionTraits<T>::BaseType> &errors() const
Bool errors(Vector<typename FunctionTraits<T>::BaseType> &err) const
Thrown Exceptions
void compuCovariance(Matrix<Double> &cov)
void buildNormalMatrix (const Vector<typename FunctionTraits<T>::BaseType> &x, const Vector<typename FunctionTraits<T>::BaseType> &y, const Vector<typename FunctionTraits<T>::BaseType> &sigma, const Vector<Bool> *const mask=0)
void buildNormalMatrix (const Matrix<typename FunctionTraits<T>::BaseType> &x, const T<typename FunctionTraits<T>::BaseType> &y, const T<typename FunctionTraits<T>::BaseType> &sigma, const T<Bool> *const mask=0)
void buildNormalMatrix (const Vector<typename FunctionTraits<T>::BaseType> &x, const Vector<typename FunctionTraits<T>::BaseType> &y, const Vector<Bool> *const mask=0)
void buildNormalMatrix (const Matrix<typename FunctionTraits<T>::BaseType> &x, const T<typename FunctionTraits<T>::BaseType> &y, const T<Bool> *const mask=0)
Bool residual(Vector<typename FunctionTraits<T>::BaseType> &y, const T<typename FunctionTraits<T>::BaseType> &x, const Vector<typename FunctionTraits<T>::BaseType> &sol, const Bool model=False)
Return the residual after a fit in y. x can
be a vector (if 1D function) or a matrix (ND functional), as in the
fit() methods. If sol is given, it is the solution derived from
a fit and its value will be used; otherwise only the parameters
in the fitted functional will be used.
If model is given as True, the model, rather
the residual <data>-<model> will be returned in y.
False is returned if residuals cannot be calculated.
Bool residual(Vector<typename FunctionTraits<T>::BaseType> &y, const T<typename FunctionTraits<T>::BaseType> &x, const Bool model=False)
Thrown Exceptions
uInt getRank() const
Get the rank of the solution (or zero of no fit() done yet). A
valid solution will have the same rank as the number of unknowns (or
double that number in the complex case). For SVD solutions the
rank could be less.
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) = 0
Generalised fitter
void buildMatrix(const Array<typename FunctionTraits<T>::BaseType> &x, const T<typename FunctionTraits<T>::BaseType> &y, const T<typename FunctionTraits<T>::BaseType> *const sigma, const T<Bool> *const mask=0)
Build the normal matrix
void buildConstraint()
Build the constraint equations
void fillSVDConstraints()
Get the SVD constraints
Bool buildResidual(Vector<typename FunctionTraits<T>::BaseType> &y, const T<typename FunctionTraits<T>::BaseType> &x, const Vector<typename FunctionTraits<T>::BaseType> *const sol, const Bool model=False)
Calculate residuals
typename FunctionTraits<T>::BaseType getVal_p(const T<typename FunctionTraits<T>::BaseType> &x, uInt j, uInt i) const
Function to get evaluated functional value
void initfit_p(uInt parcnt)
Initialise the fitter with number of solvable parameters
uInt testInput_p (const Array<typename FunctionTraits<T>::BaseType> &x, const T<typename FunctionTraits<T>::BaseType> &y, const T<typename FunctionTraits<T>::BaseType> *const sigma)
Return number of condition equations and check sizes x, y, sigma
Thrown Exceptions
void resetFunction()
Reset all the input
void setFunctionEx()
Set function properties
Bool setConstraintEx(const uInt n, const Vector<typename FunctionTraits<T>::BaseType> &x, const typename FunctionTraits<T>::BaseType y)
Set Constraint properties