Newer
Older
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
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
;+
; NAME:
; MPFIT
;
; 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:
; Perform Levenberg-Marquardt least-squares minimization (MINPACK-1)
;
; MAJOR TOPICS:
; Curve and Surface Fitting
;
; CALLING SEQUENCE:
; parms = MPFIT(MYFUNCT, start_parms, FUNCTARGS=fcnargs, NFEV=nfev,
; MAXITER=maxiter, ERRMSG=errmsg, NPRINT=nprint, QUIET=quiet,
; FTOL=ftol, XTOL=xtol, GTOL=gtol, NITER=niter,
; STATUS=status, ITERPROC=iterproc, ITERARGS=iterargs,
; COVAR=covar, PERROR=perror, BESTNORM=bestnorm,
; PARINFO=parinfo)
;
; DESCRIPTION:
;
; MPFIT uses the Levenberg-Marquardt technique to solve the
; least-squares problem. 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 MINPACK-1 (LMDIF.F) by More' 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.
;
; The Levenberg-Marquardt technique is a particular strategy for
; iteratively searching for the best fit. This particular
; implementation is drawn from MINPACK-1 (see NETLIB), and seems to
; be more robust than routines provided with IDL. This version
; allows upper and lower bounding constraints to be placed on each
; parameter, or the parameter can be held fixed.
;
; The IDL user-supplied function should return an array of weighted
; deviations between model and data. 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:
;
; DEVIATES = (Y - F(X)) / ERR
;
; where F is the function representing the model. You are
; recommended to use the convenience functions MPFITFUN and
; MPFITEXPR, which are driver functions that calculate the deviates
; for you. If ERR are the 1-sigma uncertainties in Y, then
;
; TOTAL( DEVIATES^2 )
;
; will be the total chi-squared value. MPFIT will minimize the
; chi-square value. The values of X, Y and ERR are passed through
; MPFIT to the user-supplied function via the FUNCTARGS keyword.
;
; Simple constraints can be placed on parameter values by using the
; PARINFO keyword to MPFIT. See below for a description of this
; keyword.
;
; MPFIT does not perform more general optimization tasks. See TNMIN
; instead. MPFIT is customized, based on MINPACK-1, to the
; least-squares minimization problem.
;
; USER FUNCTION
;
; The user must define a function which returns the appropriate
; values as specified above. The function should return the weighted
; deviations between the model and the data. For applications which
; use finite-difference derivatives -- the default -- the user
; function should be declared in the following way:
;
; FUNCTION MYFUNCT, p, X=x, Y=y, ERR=err
; ; Parameter values are passed in "p"
; model = F(x, p)
; return, (y-model)/err
; END
;
; See below for applications with explicit derivatives.
;
; The keyword parameters X, Y, and ERR in the example above are
; suggestive but not required. Any parameters can be passed to
; MYFUNCT by using the FUNCTARGS keyword to MPFIT. Use MPFITFUN and
; MPFITEXPR if you need ideas on how to do that. The function *must*
; accept a parameter list, P.
;
; In general there are no restrictions on the number of dimensions in
; X, Y or ERR. However the deviates *must* be returned in a
; one-dimensional array, and must have the same type (float or
; double) as the input arrays.
;
; See below for error reporting mechanisms.
;
;
; CHECKING STATUS AND HANNDLING ERRORS
;
; Upon return, MPFIT will report the status of the fitting operation
; in the STATUS and ERRMSG keywords. The STATUS keyword will contain
; a numerical code which indicates the success or failure status.
; Generally speaking, any value 1 or greater indicates success, while
; a value of 0 or less indicates a possible failure. The ERRMSG
; keyword will contain a text string which should describe the error
; condition more fully.
;
; By default, MPFIT will trap fatal errors and report them to the
; caller gracefully. However, during the debugging process, it is
; often useful to halt execution where the error occurred. When you
; set the NOCATCH keyword, MPFIT will not do any special error
; trapping, and execution will stop whereever the error occurred.
;
; MPFIT does not explicitly change the !ERROR_STATE variable
; (although it may be changed implicitly if MPFIT calls MESSAGE). It
; is the caller's responsibility to call MESSAGE, /RESET to ensure
; that the error state is initialized before calling MPFIT.
;
; User functions may also indicate non-fatal error conditions using
; the ERROR_CODE common block variable, as described below under the
; MPFIT_ERROR common block definition (by setting ERROR_CODE to a
; number between -15 and -1). When the user function sets an error
; condition via ERROR_CODE, MPFIT will gracefully exit immediately
; and report this condition to the caller. The ERROR_CODE is
; returned in the STATUS keyword in that case.
;
;
; EXPLICIT DERIVATIVES
;
; In the search for the best-fit solution, MPFIT by default
; calculates derivatives numerically via a finite difference
; approximation. The user-supplied function need not calculate the
; derivatives explicitly. However, the user function *may* calculate
; the derivatives if desired, but only if the model function is
; declared with an additional position parameter, DP, as described
; below. If the user function does not accept this additional
; parameter, MPFIT will report an error. As a practical matter, it
; is often sufficient and even faster to allow MPFIT to calculate the
; derivatives numerically, but this option is available for users who
; wish more control over the fitting process.
;
; There are two ways to enable explicit derivatives. First, the user
; can set the keyword AUTODERIVATIVE=0, which is a global switch for
; all parameters. In this case, MPFIT will request explicit
; derivatives for every free parameter.
;
; Second, the user may request explicit derivatives for specifically
; selected parameters using the PARINFO.MPSIDE=3 (see "CONSTRAINING
; PARAMETER VALUES WITH THE PARINFO KEYWORD" below). In this
; strategy, the user picks and chooses which parameter derivatives
; are computed explicitly versus numerically. When PARINFO[i].MPSIDE
; EQ 3, then the ith parameter derivative is computed explicitly.
;
; The keyword setting AUTODERIVATIVE=0 always globally overrides the
; individual values of PARINFO.MPSIDE. Setting AUTODERIVATIVE=0 is
; equivalent to resetting PARINFO.MPSIDE=3 for all parameters.
;
; Even if the user requests explicit derivatives for some or all
; parameters, MPFIT will not always request explicit derivatives on
; every user function call.
;
; EXPLICIT DERIVATIVES - CALLING INTERFACE
;
; When AUTODERIVATIVE=0, the user function is responsible for
; calculating the derivatives of the *residuals* with respect to each
; parameter. The user function should be declared as follows:
;
; ;
; ; MYFUNCT - example user function
; ; P - input parameter values (N-element array)
; ; DP - upon input, an N-vector indicating which parameters
; ; to compute derivatives for;
; ; upon output, the user function must return
; ; an ARRAY(M,N) of derivatives in this keyword
; ; (keywords) - any other keywords specified by FUNCTARGS
; ; RETURNS - residual values
; ;
; FUNCTION MYFUNCT, p, dp, X=x, Y=y, ERR=err
; model = F(x, p) ;; Model function
; resid = (y - model)/err ;; Residual calculation (for example)
;
; if n_params() GT 1 then begin
; ; Create derivative and compute derivative array
; requested = dp ; Save original value of DP
; dp = make_array(n_elements(x), n_elements(p), value=x[0]*0)
;
; ; Compute derivative if requested by caller
; for i = 0, n_elements(p)-1 do if requested(i) NE 0 then $
; dp(*,i) = FGRAD(x, p, i) / err
; endif
;
Loading full blame...