Author: Hans Terlouw <gipsy@astro.rug.nl>
This module provides the class Fitter, which uses the implementation in C of MPFIT, Craig Markwardt’s non-linear least squares curve fitting routines for IDL. MPFIT uses the Levenberg-Marquardt technique to solve the least-squares problem, which is a particular strategy for iteratively searching for the best fit. In its typical use, MPFIT will be used to fit a user-supplied function (the “model”) to user-supplied data points (the “data”) by adjusting a set of parameters. MPFIT is based upon the robust routine MINPACK-1 (LMDIF.F) by Moré and collaborators.
For example, a researcher may think that a set of observed data points is best modelled with a Gaussian curve. A Gaussian curve is parameterized by its mean, standard deviation and normalization. MPFIT will, within certain constraints, find the set of parameters which best fits the data. The fit is “best” in the least-squares sense; that is, the sum of the weighted squared differences between the model and data is minimized.
See also
Tutorial with background information for this module and practical examples.
Parameters: |
|
---|
Objects of this class are callable and return the fitted parameters when called.
Residuals function
The residuals function must return a NumPy (dtype=’d’) array with weighted deviations between the model and the data. It takes two arguments: a NumPy array containing the parameter values and a reference to the attribute data which can be any object containing information about the data to be fitted. E.g., a tuple like (xvalues, yvalues, errors).
In a typical scientific problem the residuals should be weighted so that each deviate has a Gaussian sigma of 1.0. If x represents values of the independent variable, y represents a measurement for each value of x, and err represents the error in the measurements, then the deviates could be calculated as follows:
where f is the analytical function representing the model. If err are the 1-sigma uncertainties in y, then
will be the total chi-squared value. Fitter will minimize this value. As described above, the values of x, y and err are passed through Fitter to the residuals function via the attribute data.
Derivatives function
The optional derivates function can be used to compute weighted function derivatives, which are used in the minimization process. This can be useful to save time, or when the derivative is tricky to evaluate numerically.
The function takes three arguments: a NumPy array containing the parameter values, a reference to the attribute data and a list with boolean values corresponding with the parameters. If a boolean in the list is True, the derivative with respect to the corresponding parameter should be computed, otherwise it may be ignored. Fitter determines these flags depending on how derivatives are specified in item side of the attribute parinfo, or whether the parameter is fixed.
The function must return a NumPy array with partial derivatives with respect to each parameter. It must have shape (n,m), where n is the number of parameters and m the number of data points.
Configuration attributes
The following attributes can be set by the user to specify a Fitter object’s behaviour.
Required attribute. A NumPy array, a tuple or a list with the initial parameters values.
Required attribute. Python object with information for the residuals function and the derivatives function. See above.
A list of dictionaries with parameter contraints, one dictionary per parameter, or None if not given. Each dictionary can have zero or more items with the following keys and values:
'fixed': a boolean value, whether the parameter is to be held fixed or not. Default: not fixed.
'limits': a two-element tuple or list with upper end lower parameter limits or None, which indicates that the parameter is not bounded on this side. Default: no limits.
'step': the step size to be used in calculating the numerical derivatives. Default: step size is computed automatically.
'side': the sidedness of the finite difference when computing numerical derivatives. This item can take four values:
0 - one-sided derivative computed automatically (default)
1 - one-sided derivative \((f(x+h) - f(x) )/h\)
-1 - one-sided derivative \((f(x) - f(x-h))/h\)
2 - two-sided derivative \((f(x+h) - f(x-h))/2h\)
3 - user-computed explicit derivatives
where \(h\) is the value of the parameter 'step' described above. The “automatic” one-sided derivative method will chose a direction for the finite difference which does not violate any constraints. The other methods do not perform this check. The two-sided method is in principle more precise, but requires twice as many function evaluations. Default: 0.
'deriv_debug': boolean to specify console debug logging of user-computed derivatives. True: enable debugging. If debugging is enabled, then 'side' should be set to 0, 1, -1 or 2, depending on which numerical derivative you wish to compare to. Default: False.
As an example, consider a function with four parameters of which the first parameter should be fixed and for the third parameter explicit derivatives should be used. In this case, parinfo should have the value [{'fixed': True}, None, {'side': 3}, None] or [{'fixed': True}, {}, {'side': 3}, {}].
Relative \(\chi^2\) convergence criterium. Default: 1e-10
Relative parameter convergence criterium. Default: 1e-10
Orthogonality convergence criterium. Default: 1e-10
Finite derivative step size. Default: 2.2204460e-16 (MACHEP0)
Initial step bound. Default: 100.0
Range tolerance for covariance calculation. Default: 1e-14
Maximum number of iterations. Default: 200
Maximum number of function evaluations. Default: 0 (no limit)
Result attributes
After calling the method fit(), the following attributes are available to the user:
A NumPy array, list or tuple with the fitted parameters. This attribute has the same type as params0.
Final parameter uncertainties (\(1 \sigma\))
Final parameter covariance (NumPy-) matrix.
Final \(\chi^2\).
Starting value of \(\chi^2\).
Minimum reduced \(\chi^2\).
Standard errors.
Number of parameters.
Number of free parameters.
Number of pegged parameters.
Number of degrees of freedom.
Final residuals.
Number of iterations.
Number of function evaluations.
mpfit.c’s version string.
Fitting status code.
Message string.
Methods:
Perform a fit with the current values of parameters and other attributes.
Optional argument params0: initial fitting parameters. (Default: previous initial values are used.)
This method requires SciPy. After the method fit() has been called, this method calculates the upper and lower value of the confidence interval for all elements of the NumPy array x. The model values and the arrays with confidence limits are returned and can be used to plot confidence bands.
Parameters: |
|
---|---|
Returns: | A tuple with the following elements, each one is a Numpy array:
|
Note
If parameters were fixed in the fit, the corresponding error is 0 and there is no contribution to the confidence interval.
In principle, the process of computing explicit derivatives should be straightforward. In practice, the computation can be error prone, often being wrong by a sign or a scale factor.
In order to be sure that the explicit derivatives are correct, for debugging purposes the user can set the attribute parinfo[‘deriv_debug’] = True for any parameter. This will cause Fitter.fit() to print both explicit derivatives and numerical derivatives to the console so that the user can compare the results.
When debugging derivatives, it is important to set parinfo[‘side’] to the kind of numerical derivative to compare with: it should be set to 0, 1, -1, or 2, and not set to 3. When parinfo[‘deriv_debug’] is set for a parameter, then Fitter.fit() automatically understands to request user-computed derivatives.
The console output will be sent to the standard output, and will appear as a block of ASCII text like this:
1 2 3 4 5 6 7 8 | FJAC DEBUG BEGIN
# IPNT FUNC DERIV_U DERIV_N DIFF_ABS DIFF_REL
FJAC PARM 1
.... derivative data for parameter 1 ....
FJAC PARM 2
.... derivative data for parameter 2 ....
.... and so on ....
FJAC DEBUG END
|
which is to say, debugging data will be bracketed by pairs of “FJAC DEBUG” BEGIN/END phrases. Derivative data for individual parameter i will be labeled by “FJAC PARM i”. The columns are, in order,
IPNT - data point number \(j\)
FUNC - residuals function evaluated at \(x_j\)
DERIV_U - user-calculated derivative \({\partial f(x_j)}/{\partial p_i}\)
DERIV_N - numerically calculated derivative according to the value of parinfo[‘side’]
DIFF_ABS - difference between DERIV_U and DERIV_N: fabs(DERIV_U-DERIV_N)
DIFF_REL - relative difference: fabs(DERIV_U-DERIV_N)/DERIV_U
Since individual numerical derivative values may contain significant round-off errors, it is up to the user to critically compare DERIV_U and DERIV_N, using DIFF_ABS and DIFF_REL as a guide.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 | #!/usr/bin/env python
import numpy
from kapteyn import kmpfit
def residuals(p, d):
a, b, c = p
x, y, w = d
return (y - (a*x*x+b*x+c))/w
x = numpy.arange(-50,50,0.2)
y = 2*x*x + 3*x - 3 + 2*numpy.random.standard_normal(x.shape)
w = numpy.ones(x.shape)
a = [x, y, w]
f = kmpfit.Fitter(residuals, params0=[1, 2, 0], data=a)
f.fit() # call fit method
print f.params
print f.message
# result:
# [2.0001022845514451, 3.0014019147386, -3.0096629062273133]
# mpfit (potential) success: Convergence in chi-square value (1)
a[1] = 3*x*x - 2*x - 5 + 0.5*numpy.random.standard_normal(x.shape)
print f(params0=[2, 0, -1]) # call Fitter object
# result:
# [3.0000324686457871, -1.999896340813663, -5.0060187435412962]
|
Simple interface to Fitter.
Parameters: |
|
---|---|
Returns: | a Fitter object from which the fit results can be extracted. |