Which MPFIT Files Should I Download and Use?
The easiest solution is to download the entire set of routines as
either a ZIP file or a gzipped TAR file, and then extract them
in your IDL path. However it is possible to selectively download what
you need. Bear in mind that you will always need MPFIT.PRO which is the main fitting
engine.
For standard fitting of 1D curves, where the model is a
compiled function, download MPFITFUN.PRO and MPFIT.PRO. For specialized peak fitting,
use MPFITPEAK.PRO in combination with
these two.
For fitting of 2D images or surfaces, download MPFIT2DFUN.PRO and MPFIT.PRO. For peak fitting, use MPFIT2DPEAK.PRO in combination with
these two.
For a drop-in replacement for the IDL-supplied CURVEFIT in legacy
code, use MPCURVEFIT.PRO along with
MPFIT.PRO.
For cases where you don't have a precompiled function,
either at the command line on within your program, you can use an IDL
expression. Download MPFITEXPR.PRO
and MPFIT.PRO.
Is There a Citation for MPFIT?
I regularly get asked if there is a way for users to cite MPFIT in
their scientific publications. The answer is yes! I presented a
paper about MPFIT at the ADASS XVIII conference in Quebec, Canada, in
Nov, 2008, with proceedings published
by the Astronomical Society of the Pacific (2009). I welcome you to cite
this publication in your own works.
-
Markwardt, C. B. 2009, "Non-Linear Least Squares Fitting in IDL
with MPFIT," in proc. Astronomical Data Analysis Software and
Systems XVIII, Quebec, Canada, ASP Conference Series, Vol. 411, eds.
D. Bohlender, P. Dowler & D. Durand (Astronomical Society of the
Pacific: San Francisco), p. 251-254
(ISBN: 978-1-58381-702-5; Link to ASP title listing)
ADS Bibcode: 2009ASPC..411..251M (click for Bibtex and other citation formats)
Arxiv preprint: arXiv:0902.2850v1
- Bibtex entry:
\bibitem[Markwardt(2009)]{2009ASPC..411..251M} Markwardt, C.~B.\ 2009,
Astronomical Data Analysis Software and Systems XVIII, 411, 251
- Referring to the MPFIT website. Please use:
http://purl.com/net/mpfit
in your citations. Right now this link redirects to my Wisconsin
website, but the Wisconsin website will not exist forever. When I move
the website, I will change the purl.com
redirect as well.
At the same time, I also urge you to cite the work of the original
designer of the MINPACK algorithm, Jorge Moré, especially the
first listed citation below.
-
Moré, J. 1978, "The Levenberg-Marquardt Algorithm:
Implementation and Theory," in Numerical Analysis,
vol. 630, ed. G. A. Watson (Springer-Verlag: Berlin), p. 105
(DOI: 10.1007/BFb0067690; Link to Springer title listing)
-
Moré, J. & Wright, S. 1993, "Optimization Software Guide,"
SIAM, Frontiers in Applied Mathematics, Number 14.
(ISBN: 978-0-898713-22-0; Link to SIAM title listing)
-
MINPACK-1, Jorge More', available from netlib.
http://www.netlib.org/minpack
Are there other versions of MPFIT?
Yes! MPFIT is based on the original MINPACK-1 library, written in
FORTRAN, available from Netlib at http://netlib.org/minpack.
I'm also pleased to make available a C library version of MPFIT.
This library has many of the same capabilities as the IDL version,
including parameter bounds and choice of numerical or explicit
derivatives, in a small, fast library, with simple calling interface.
Of course, some IDL-specific elements such as _EXTRA and "tied"
parameters don't make as much sense in C, and they are not present.
See the C Version of MPFIT page for more
information. The original MINPACK-1 code was translated to C by
Stephen Moshier (http://moshier.net/), from which MPFIT
borrows extensively.
MPFIT was translated to Python's
Numeric library in 2002 by Mark Rivers, and can be found at his
website. In 2009, Sergei Koposov ported that code to use the more
modern Python Numpy numerical
library, and the result can be found in his astrolibpy library.
How can I calculate the best-fit model?
All of the basic functions in the library return the best-fit model
function using the keyword YFIT. Simply pass a named variable with
this keyword and upon return the best-fit model will be in that
variable.
If you change the parameter values manually, it is still simple to
recompute the model function. It differs, depending on what type of
model you are fitting:
yfit = model(x, p) | 1D Model named model |
zfit = model(x, y, p) | 2D Model named model |
yfit = MPEVALEXPR(expr, x, p) | Expression named expr |
Here, p is assumed to be the parameter set.
How do I fit a 2D image?
The simplest answer is to use the MPFIT2DFUN.PRO function. You supply
the image, the X and Y labels for each pixel, and a 2D model
function.
Can I fit a function of several variables?
This is a question when you are fitting a function of several
independent variables, such as this:
y = f(x0,x1,x2, ...; p0, p1, p2, ...)
where the xi are the independent variables and the pi
are fitting parameters.
Most of the MPFIT functions do not care how many independent
variables there are. In fact, the main fitting routine MPFIT does not
accept any independent variables at all! Instead, they are considered
to be implicit to the problem. Functions like MPFITFUN do accept an
independent variable called X. [ Since MPFIT itself doesn't deal with
the independent variables, MPFITFUN creates a common block with that
information so that the model function can gain access to it. ]
Therefore, if you are using plain MPFIT to do your fitting, then
you have the freedom to construct your model in any way you
please.
If you are using MPFITFUN or MPFIT2DFUN, then your model should
still be a function of a single variable X, and you should pass a
single set of independent variables called X to the fitting routine.
However, X can be arbitrarily complicated, so you can have it be a 2D
array, the first row containing "x0" values, the second containing
"x1", and so on. Thus the burden falls to your IDL model function to
decode the contents of X. I recommend the 2D array approach, but you
can also simply concatenate the variables in a 1D vector, or pass them
by COMMON block (not recommended).
The array approach I am advocating would be something like this:
p = mpfitfun([x0, x1, x2], y, yerr, pstart, ...)
This creates a new "X" variable which is the array concatenation of
all of your independent variables. Then the burden would be upon your
user function to extract the relevant quantities from this array.
Help! I can't get my fitting to work!
Generally speaking, this is not something I can help with. I am
probably not an expert in your field, and I am definitely not an
expert in optimization. However, I think there are a few very important
things that you must do.
First, you must make sure that the problem is well defined. Make
sure that you know exactly what function you plan on fitting. The
MPFIT functions are not psychic; they can't figure out how to solve
your problem for you.
I see quite a few people who can't get their model functions right.
The best kind of model function is self-contained, not depending on
any outside data or common blocks. Be sure that your model function
works by itself before throwing MPFIT into the mix. It is very
difficult to debug the problem when the added layer of MPFIT is
disguising everything. Be aware of domain problems. For example, if
you will be taking the square root of a parameter, you had better
constrain it to be positive.
Finally, and this can't be stressed enough, it is crucial to
estimate the starting parameters as best you can. Initializing the
parameters to all zeroes is actually the worst thing you can
do, since then the problem becomes scale-less, and MPFIT has a much
harder time deciding what to do.
Do I need to compute derivatives?
The short answer is, probably not. I have found that in most cases
the automatic finite difference technique used by MPFIT is sufficient.
Computing the derivatives explicitly seems to slow things down, and
can even make a worse final solution. My suggestion is to try it
first without computing the derivatives, and them implement
derivatives as needed (and the AUTODERIVATIVE=0 keyword).
However, see the caveats below.
How can I embed MPFIT-style fitting in a widget application?
My belief is that the best solution is to use the ITERPROC
feature, which in turn would call WIDGET_EVENT manually to dispatch
events. Rob Dimeo has successfully integrated this approach into his
dedicated peak fitting program named PAN.
PAN is the winner of an honorable mention the recent IDL programming
contest. Source code is available (uses objects). Nice job, Rob!
Here is a small screenshot:
Small screen shot of PAN
Why does MPFIT just return my parameters without fitting? / How do I check for processing errors?
Users occcasional report that MPFIT, or one of its drivers, return
without fitting any of the parameters. This is usually an indication
that MPFIT discovered an error in the user-specified parameters and
returned immediately.
Users should check the values of the STATUS and ERRMSG keywords
when MPFIT returns. If STATUS is less than or equal to 0, then a
processing error has occurred. A more detailed description of the
problem can be found by consulting the ERRMSG keyword.
A simple example of error checking using STATUS and ERRMSG might be:
p = MPFITFUN(x, y, dy, p0, STATUS=status, ERRMSG=errmsg)
if status LE 0 then message, errmsg
Of course, more sophisticated error handling approaches can be
implemented as needed.
Why does MPFIT not converge properly?
For certain problem sets, MPFIT will not converge to the "optimal"
solution. It may seem obvious to you that there is a better solution,
but MPFIT does not find it. Why could this be?
Be sure you are checking for error
conditions. MPFIT may be indicating an error condition which
prematurely terminates the fit.
MPFIT normally uses a finite difference approximation to estimate
derivatives. This means that it varies each parameter by a small
amount, and measures the corresponding variation in the chi-square
value. However, if your data or model are discretized for some reason
then this approximation can fail, which then causes the fit not to
converge. You should check for:
- Automatically-selected step size too small. In some cases,
the automatically-selected step sized used to compute derivatives is
too small. This usually happens when the dynamic range of the
parameter is very small compared to the absolute value. Solution: set
the step size manually with PARINFO.STEP or .RELSTEP, to be small
enough to resolve the variations in your fitting function. Another
solution is to re-write your model function to remove large constant
offsets from parameters.
- Discretized parameters. Example: a tabulated model, which is
evaluated on a discrete grid of parameters. Solution: Interpolate the
table smoothly between parameter values; set PARINFO.STEP or .RELSTEP
values to be comparable to the grid spacing.
- Discontinuous functions. Example: use of the ABS() function, or
any function which has a discontinous derivative. Solution: use
analytical derivatives; or use a smoother approximation to your
discontinuous function.
- Undersampled data. Example: fitting of undersampled peaks.
For example, consider data sampled with a grid spacing of 1 unit,
which contains a peak whose ideal full-width at half-max is 0.5 unit.
Since this peak is undersampled, the small parameter variations
attempted by MPFIT may have no measureable effect. Solution: get data
with finer sampling; adjust PARINFO.STEP or .RELSTEP to be comparable
to the grid spacing.
Also, beware that PARINFO.RELSTEP will not work if the parameter
value is zero. Generally speaking, all starting parameter values
should be non-zero.
Following along those lines, it's also worth checking what the
dynamic range of your parameters is. For example, if you are
fitting a model of 'p(0)+p(1)*x' where p(0) is 1 and p(1) is 1d-12,
then the fitting won't do very well. There are matrix operations
inside of MPFIT that may lose precision if there is a large dynamic
range in parameter values. It's better to absorb the 1d-12 factor
into the fitting function, so that all of the parameters are of
equivalent magnitude.
Finally, beware that if your chi-square function has local
minima then MPFIT may become irretrievably stuck. MPFIT is not a
global optimizer.
|