Data-fitting

Header file datafit.h

Functions

See also methods to fit surfaces to data.

Polynomial fitting

lsqline
Finds best-fitting straight line and standard deviation, with intercept calculated or specified. Equation is of the form y = a.x + b. 
int lsqline(
   // Data in:
   const Vector& x, // ordinates
   const Vector& y, // y(x) data values
   const int fixb, // 1: define bhat / 0: fit it
   const double eps, // numerical bandwidth to prevent overflow on division
   // Data in/out:
   double& bhat, // best-fitting (or supplied) intercept
   // Data out:
   double& ahat, // best-fitting gradient
   double& sd // standard deviation
)
Returns 1 on success; 0 if failed because gradient or intercept calculation would involve division by small number.
polyfitg
Fit polynomial to data. Solves overdetermined system of equations by Givens rotations. 
Array polyfitg(
   const Array& x, // Array of x-ordinates
   const Array& y, // Array of y-values
   const Array& w, // Array of weights
   const int n, // number of polynomial coefficients, i.e. order + 1
   double& res, // residual of l2 fit (returned value)
   const double eps = 1.0e-8 // numerical bandwidth for spotting singular matrix
)
Returns polynomial coefficients (order of increasing powers of x), or an Array of length zero if it fails.
polyfitn
Fit polynomial to data. Solves normal equations. This function is included for interest: polyfitg is the recommended function. 
Array polyfitn(
   const Array& x, // Array of x-ordinates
   const Array& y, // Array of y-values
   const Array& w, // Array of weights
   const int n, // number of polynomial coefficients, i.e. order + 1
   const double eps = 1.0e-8 // numerical bandwidth for spotting singular matrix
)
Returns polynomial coefficients (order of increasing powers of x), or an Array of length zero if it fails.

Multilinear fitting

Least squares fit in n linear parameters
Best n-dimensional least squares fit to weighted data using Givens rotations to solve overdetermined system of equations 
Array linfit(
   const Matrix& xf, // xf[i][j]: function j's value at data point i
   const Vector& y, // data values
   const Vector& w, // weight of each data value
   double& res, // residual of fit
   const double eps = 1.0e-8 // numerical bandwidth: singular matrix
)
natural fit to data defined at the corners of a rectangle
Array frectfit(
   const Array& x, // x ordinates of corners (2)
   const Array& y, // y ordinates of corners (2)
   const matrix& f, // function values (2x2)
   const eps = 1.0e-8 // numerical bandwidth
)
Returns coefficients for the functions {1,x,y,xy} respectively, or an Array of no elements if the fitting fails.

Utilities

Work out which run of points to use for fitting tabular data
int fitstart(
   const Array& r, // Array of ordinates
   const double rr, // value around which fit is to be performed
   const int nr // number of points required for the fit
)
(returns index of first point)

Generic classes

Class datafit acts as a wrapper for the fitting of general functions (especially ones nonlinear in the adjustable parameters) to data. This class provides alternative parameter-fitting methods. To use, a program defines a subclass inheriting from datafit, and providing instantiations of essential methods describing the fitting functions and their derivatives. The exact functions required depend on the method used.

Specifications of the functions are: 

// evaluate function given coefficient array and index of data point
double func(const Array& co, const int ipt) const

// evaluate row of Jacobian matrix: derivative of function with respect
// to coefficients, at specified data point
Vector jacob(const Array& co, const int ipt) const

// evaluate Hessian matrix: second derivatives of function with respect
// to coefficients, at specified data point
Matrix jacob(const Array& co, const int ipt) const

Fitting methods are described below:

Gauss-Newton
Least squares fit in n non-linear parameters. Uses a Gauss-Newton scheme, and tries both directions along steepest gradient. Requires function evaluations and the Jacobian matrix. Returns 1 on success, 0 on failure. 
   int gaussnewton(
      const Array& y, // data values
      const Array& w, // weight of each data value
      Array& co, // coefficients: guess in, optimised out
      int& hitmax, // reached maximum iterations before finding solution
      int& sing, // Jacobian was close to singular: could not proceed
      int& stuck, // could not find valid parameter change
      int& nits, // iterations taken
      double& el2norm, // least squares norm in final answer
      const double acc = 1.0e-6, // accuracy
      const double eps = 1.0e-8, // bandwidth: singular matrix or tiny steps
      const double dampf = 0.1, // damping factor
      const double dampfac = 0.5, // change factor for damping factor
      const int maxits = 0, // maximum iterations (0 => no limit)
      const int verbose = 0 // diagnostic output
   )
Levenberg-Marquardt
Least squares fit in n non-linear parameters. Uses a Levenberg-Marquardt scheme, and tries both directions along steepest gradient. Requires function evaluations and the Jacobian matrix. Returns 1 on success, 0 on failure. 
   int levmar(
      const Array& y, // data values
      const Array& w, // weight of each data value
      Array& co, // coefficients: guess in, optimised out
      int& hitmax, // reached maximum iterations before finding solution
      int& sing, // Jacobian was close to singular: could not proceed
      int& stuck, // could not find valid parameter change
      int& nits, // iterations taken
      double& el2norm, // least squares norm in final answer
      const double tregin, // initial trust region
      const double gamstep, // initial jump in growth parameter
      const double dfac = 2, // change factor for growth and trust region
      const double acc = 1.0e-6, // accuracy
      const double eps = 1.0e-8, // bandwidth: singular matrix or tiny steps
      const int bothways = 1, // to alternate direction of search for robustness      const int maxits = 0, // maximum iterations (0 => no limit)
      const int verbose = 0 // diagnostic output
   )
Gauss-Newton
Least squares fit in n non-linear parameters Uses Gauss-Newton scheme with the Hessian matrix of second-order corrections, and tries both directions along steepest gradient. Requires function evaluations, the Jacobian matrix and the Hessian matrix. Returns 1 on success, 0 on failure. 
   int gaussnewtonhess(
      const Array& y, // data values
      const Array& w, // weight of each data value
      Array& co, // coefficients: guess in, optimised out
      int& hitmax, // reached maximum iterations before finding solution
      int& sing, // Jacobian was close to singular: could not proceed
      int& stuck, // could not find valid parameter change
      int& nits, // iterations taken
      double& el2norm, // least squares norm in final answer
      const double cutoff = 0.1, // error below which Hessian is not calculated
      const double acc = 1.0e-6, // accuracy
      const double eps = 1.0e-8, // bandwidth: singular matrix or tiny steps
      const double dampf = 0.1, // damping factor
      const double dampfac = 0.5, // change factor for damping factor
      const int maxits = 0, // maximum iterations (0 => no limit)
      const int verbose = 0 // diagnostic output
   )