MPFIT2DPEAK Name
MPFIT2DPEAK Author
Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
craigm@lheamail.gsfc.nasa.gov
UPDATED VERSIONs can be found on my WEB PAGE:
http://cow.physics.wisc.edu/~craigm/idl/idl.html
Purpose
Fit a gaussian, lorentzian or Moffat model to data
Major Topics
Curve and Surface Fitting
Calling Sequence
yfit = MPFIT2DPEAK(Z, A [, X, Y, /TILT ...] )
Description
MPFIT2DPEAK fits a gaussian, lorentzian or Moffat model using the
non-linear least squares fitter MPFIT. MPFIT2DPEAK is meant to be
a drop-in replacement for IDL's GAUSS2DFIT function (and requires
MPFIT and MPFIT2DFUN).
The choice of the fitting function is determined by the keywords
GAUSSIAN, LORENTZIAN and MOFFAT. By default the gaussian model
function is used. [ The Moffat function is a modified Lorentzian
with variable power law index. ] The two-dimensional peak has
independent semimajor and semiminor axes, with an optional
rotation term activated by setting the TILT keyword. The baseline
is assumed to be a constant.
GAUSSIAN A[0] + A[1]*exp(-0.5*u)
LORENTZIAN A[0] + A[1]/(u + 1)
MOFFAT A[0] + A[1]/(u + 1)^A[7]
u = ( (x-A[4])/A[2] )^2 + ( (y-A[5])/A[3] )^2
where x and y are cartesian coordinates in rotated
coordinate system if TILT keyword is set.
The returned parameter array elements have the following meanings:
A[0] Constant baseline level
A[1] Peak value
A[2] Peak half-width (x) -- gaussian sigma or half-width at half-max
A[3] Peak half-width (y) -- gaussian sigma or half-width at half-max
A[4] Peak centroid (x)
A[5] Peak centroid (y)
A[6] Rotation angle (radians) if TILT keyword set
A[7] Moffat power law index if MOFFAT keyword set
By default the initial starting values for the parameters A are
estimated from the data. However, explicit starting values can be
supplied using the ESTIMATES keyword. Also, error or weighting
values can optionally be provided; otherwise the fit is
unweighted. Restrictions
If no starting parameter ESTIMATES are provided, then MPFIT2DPEAK
attempts to estimate them from the data. This is not a perfect
science; however, the author believes that the technique
implemented here is more robust than the one used in IDL's
GAUSS2DFIT. The author has tested cases of strong peaks, noisy
peaks and broad peaks, all with success.
Compatibility
This function is designed to work with IDL 5.0 or greater.
Because TIED parameters rely on the EXECUTE() function, they cannot
be used with the free version of the IDL Virtual Machine.
Inputs
Z - Two dimensional array of "measured" dependent variable values.
Z should be of the same type and dimension as (X # Y).
X - Optional vector of x positions for a single row of Z.
X[i] should provide the x position of Z[i,*]
Default: X values are integer increments from 0 to NX-1
Y - Optional vector of y positions for a single column of Z.
Y[j] should provide the y position of Z[*,j]
Default: Y values are integer increments from 0 to NY-1
Outputs
A - Upon return, an array of best fit parameter values. See the
table above for the meanings of each parameter element.
Returns
Returns the best fitting model function as a 2D array.
Keywords
** NOTE ** Additional keywords such as PARINFO, BESTNORM, and
STATUS are accepted by MPFIT2DPEAK but not documented
here. Please see the documentation for MPFIT for the
description of these advanced options.
CHISQ - the value of the summed squared residuals for the
returned parameter values.
CIRCULAR - if set, then the peak profile is assumed to be
azimuthally symmetric. When set, the parameters A[2)
and A[3) will be identical and the TILT keyword will
have no effect.
DOF - number of degrees of freedom, computed as
DOF = N_ELEMENTS(DEVIATES) - NFREE
Note that this doesn't account for pegged parameters (see
NPEGGED).
ERROR - upon input, the measured 1-sigma uncertainties in the "Z"
values. If no ERROR or WEIGHTS are given, then the fit is
unweighted.
ESTIMATES - Array of starting values for each parameter of the
model. If ESTIMATES is not set, then the starting
values are estimated from the data directly, before
fitting. (This also means that PARINFO.VALUES is
ignored.)
Default: not set - parameter values are estimated from data.
GAUSSIAN - if set, fit a gaussian model function. The Default.
LORENTZIAN - if set, fit a lorentzian model function.
MOFFAT - if set, fit a Moffat model function.
MEASURE_ERRORS - synonym for ERRORS, for consistency with built-in
IDL fitting routines.
NEGATIVE - if set, and ESTIMATES is not provided, then MPFIT2DPEAK
will assume that a negative peak is present -- a
valley. Specifying this keyword is not normally
required, since MPFIT2DPEAK can determine this
automatically.
NFREE - the number of free parameters in the fit. This includes
parameters which are not FIXED and not TIED, but it does
include parameters which are pegged at LIMITS.
PERROR - upon return, the 1-sigma uncertainties of the parameter
values A. These values are only meaningful if the ERRORS
or WEIGHTS keywords are specified properly.
If the fit is unweighted (i.e. no errors were given, or
the weights were uniformly set to unity), then PERROR
will probably not represent the true parameter
uncertainties.
*If* you can assume that the true reduced chi-squared
value is unity -- meaning that the fit is implicitly
assumed to be of good quality -- then the estimated
parameter uncertainties can be computed by scaling PERROR
by the measured chi-squared value.
DOF = N_ELEMENTS(Z) - N_ELEMENTS(A) ; deg of freedom
PCERROR = PERROR * SQRT(BESTNORM / DOF) ; scaled uncertainties
QUIET - if set then diagnostic fitting messages are suppressed.
Default: QUIET=1 (i.e., no diagnostics)
SIGMA - synonym for PERROR (1-sigma parameter uncertainties), for
compatibility with GAUSSFIT. Do not confuse this with the
Gaussian "sigma" width parameter.
TILT - if set, then the major and minor axes of the peak profile
are allowed to rotate with respect to the image axes.
Parameter A[6] will be set to the clockwise rotation angle
of the A[2] axis in radians.
WEIGHTS - Array of weights to be used in calculating the
chi-squared value. If WEIGHTS is specified then the ERR
parameter is ignored. The chi-squared value is computed
as follows:
CHISQ = TOTAL( (Z-MYFUNCT(X,Y,P))^2 * ABS(WEIGHTS) )
Here are common values of WEIGHTS:
1D/ERR^2 - Normal weighting (ERR is the measurement error)
1D/Y - Poisson weighting (counting statistics)
1D - Unweighted
The ERROR keyword takes precedence over any WEIGHTS
keyword values. If no ERROR or WEIGHTS are given, then
the fit is unweighted. Example
; Construct a sample gaussian surface in range [-5,5] centered at [2,-3]
x = findgen(100)*0.1 - 5. & y = x
xx = x # (y*0 + 1)
yy = (x*0 + 1) # y
rr = sqrt((xx-2.)^2 + (yy+3.)^2)
; Gaussian surface with sigma=0.5, peak value of 3, noise with sigma=0.2
z = 3.*exp(-(rr/0.5)^2) + randomn(seed,100,100)*.2
; Fit gaussian parameters A
zfit = mpfit2dpeak(z, a, x, y) References
MINPACK-1, Jorge More', available from netlib (www.netlib.org).
"Optimization Software Guide," Jorge More' and Stephen Wright,
SIAM, *Frontiers in Applied Mathematics*, Number 14.
Modification History
New algorithm for estimating starting values, CM, 31 Oct 1999
Documented, 02 Nov 1999
Small documentation fixes, 02 Nov 1999
Documented PERROR for unweighted fits, 03 Nov 1999, CM
Copying permission terms have been liberalized, 26 Mar 2000, CM
Small cosmetic changes, 21 Sep 2000, CM
Corrected bug introduced by cosmetic changes, 11 Oct 2000, CM :-)
Added POSITIVE keyword, 17 Nov 2000, CM
Removed TILT in common, in favor of FUNCTARGS approach, 23 Nov
2000, CM
Added SYMMETRIC keyword, documentation for TILT, and an example,
24 Nov 2000, CM
Changed SYMMETRIC to CIRCULAR, 17 Dec 2000, CM
Really change SYMMETRIC to CIRCULAR!, 13 Sep 2002, CM
Add error messages for SYMMETRIC and CIRCLE, 08 Nov 2002, CM
Make more consistent with comparable IDL routines, 30 Jun 2003, CM
Defend against users supplying strangely dimensioned X and Y, 29
Jun 2005, CM
Convert to IDL 5 array syntax (!), 16 Jul 2006, CM
Move STRICTARR compile option inside each function/procedure, 9 Oct 2006
Add COMPATIBILITY section, CM, 13 Dec 2007
Clarify documentation regarding what happens when ESTIMATES is not
set, CM, 14 Dec 2008