- Almost all my fitting programs fit a generic function using the
CURFIT routine from Bevington's book "Data reduction
and Error Analysis for the Physical Sciences" (McGraw-Hill 1969)
The few exceptions are simple regression programs using a plain
least square routine, some polynomial fits using Bevington's POLFIT
and some surface fitting using a routine by S.McKechnie of WKS Leiden.
- CURFIT is part of the class of "generalized least square" or
Levenberg-Marquardt methods, as described e.g. in the
Numerical Recipes book, and is widely used by astronomers
(try a Google search !).
- The typical CURFIT based program must provide a Fortran function
FUNC(X,I,A,NT) where Y=FUNC(X,I,A,NT)
evaluates y=f(xi) where f(x) is the function to be fitted and
the NT elements of array A are the parameters to be fitted.
- What matters are the parameters to be fitted. If you want to parametrize
the function on other parameters but not to fit them, one handy way is to
pass them through a COMMON block, or write a wrapper. A wrapper
is also needed to hide a routine which has a different calling interface.
- On the other hand it does not matter what are the X, although
often and in many places of the code X is the independent variable
of the data to be fitted.
One could call FUNC(W,1,A,NT) to compute y=f(w) at a single scalar
One can also use X as an arbitrary array of indices. This is customary
in fitting X-ray spectra, where one is not actually fitting versus energy, but
versus channel number. Here usually X(I)=I or X(i)=I+K.
One is not even bound to fit unidimensional data. I fitted two-dimensional
functions z=f(x,y) to image data using as X a one-dimensional
index built from the pixel addresses I,J
basic CURFIT program
- The basic CURFIT program assumes that FUNC has such name and calls
it in the finite difference derivative routine, and in the chi-square computation
routine (while FUNC is passed as an EXTERNAL argument to
If desired one can replace the derivative routine with a specific routine
computing the derivative analytically (for speed), or with some other
optimized variant (see below).
If desired one could propagate the named EXTERNAL argument
to the other routines as well (I've never done this).
- The CURFIT routine does not do any fit ! It simply performs a
single fit iteration which leads from an initial guess of the
parameters to a better one which has a smaller chi-square.
(Actually there is a loop, potentially an infinite loop, inside CURFIT
in case the adjustment does not lead to a smaller chi-square, changing a
tuning variable FLAMBDA)
- Therefore the typical, plain fitting program must provide something like
- code reading NPTS data points X,Y and eventual errors SY
- code asking (or computing) the initial guess for the NT parameters to be fitted A
- initializing the finite difference increments DA
used for the computation of the partial derivatives of FUNC
- initializing FLAMBDA and MODE
- DO loop on maximum number of allowed iterations
- CALL CURFIT(X,Y,SY,NPTS,NT,MODE,A,DA,SA,FLAMBDA,YFIT,CHI)
- eventually displaying partial results
- leaving the loop if a convergence criterion is achieved
- END DO
About MODE and FLAMBDA see Bevington's code.
The convergence criterion is care of the user. Generally one checks that the
difference between the chi-square CHI of the current iteration and
the one (saved) of the previous one is below a predefined threshold. But
one could also check on the parameter values themselves to converge.
On exit A contains the fitted parameters, SA a rough
estimate of their errors, and YFIT the values of fitted function
- I tend not to trust errors computed by CURFIT this way, but prefer
to use confidence limits computed according to standard prescriptions (
Lampton, Margon & Bowyer 1976 or
selecting 1, 2 or more interesting parameters
stepping them over an equispaced grid
for each step do a fit on all remaining parameters, and save the
look for the points where the chi square differs from the minimum chi
square (the best fit) by less than a given delta chi square (which depends
on the confidence level and number of interesting parameters)
- Several of my ad-hoc programs, all the times I need to
fit a predefined model and am just interested in the parameters, are written
as the basic programs described above, i.e. use an ad-hoc FUNC
routine, and call CURFIT within an iteration loop.
- A typical variant may include the selection of a range of data to be fitted,
narrower than the original data array.
Another variant may use an ad-hoc optimized version of the FDER or
FDERIV derivative routine.
- By spectrum I mean any distribution of spectral flux vs an energy
variable (actually this can be a wavelength, a frequency or an energy).
In practice my library of spectral models supports only three "standard"
- Flambda vs lambda : erg/cm2/s/Å vs Å
- dN/dE vs E : photon/cm2/s/keV vs keV
- Fnu vs nu : erg/cm2/s/Hz vs Hz
The first is used by the UV fitting programs, the second by the X-ray
fitting programs and the third is currently not used.
- The library (see report R7, entry
for details) defines a number of spectral functions allowing to select,
setting appropriate COMMON flags via dedicated calls :
- the standard unit representation (one of the three above)
- the unit for temperature parameters (when applicable) : keV or Kelvin
- whether the user supplied normalization is comprehensive of
any physical constants (c, h, pi) or not. Current programs assume
that physical constants are included and display in their
log file the source-specific part of the normalization when
applicable (e.g. the emission measure or the part proportional
to r²/d² for a Bremsstrahlung or blackbody)
The conversions (when applicable) are either taken care by branches in
the code of each individual spectral form routine, or, when the model
is natively computed in dN/dE, by a conversion routine called MOAMED
(according to the saying
"if the mountain will not come to Mohammed, Mohammed
will go to the mountain")
As a technicality, all spectral routines are defined with 2 (or 3)
entry points as follows :
- a numbered function entry point Y=SPEnn(X,A)
- an optional alias entry point Y=alias(X,A)
- an info entry point CALL INFOnn(M,NAME,PARA)
where the former ones compute the flux Y at energy/lambda/frequency
X, while the latter information call returns the number
of parameters M. Array A(M) contains the user supplied
values of the parameters, character array PARA(M) contains their
names (and NAME the spectral form name).
- One can see that the above function has not the calling sequence
expected for FUNC. In fact this is not the function to be fitted,
and there are several reasons for this :
- the above calling sequence is handier for non-fitting usages of the
- even with single component models, the actual spectrum is usually the
product of an unsabsorbed model by (interstellar) absorption,
be it UV reddening or X-ray photoelectric absorption
- in most cases one desires a multiple component model (the present
version of the software supports only additive components, i.e.
the sum of several models, all absorbed through a common absorption
(a former version supported separate by-model absorptions) )
- finally not all parameters of a model are necessarily to be fitted
- and last but not least one is not necessarily fitting a model (for
X-rays one is fitting the convolution of a model with a response
- The library provides a set of single model wrapper routines (not
used for fitting) namely SPECUV and SPECSX (which handle
the UV reddening or the X-ray absorption, and call a common dispatcher
SPECS which uses a computed GOTO to call the desired spectral model;
and a set of multiple model wrapper routines (namely SPEUV2
and SPEX2 which handle absorption and call a common dispatcher
- The dispatcher SPESUM takes care of calling each model with its
array of parameters, but first builds the array of parameters mixing in
the correct way (using flags) the parameters fixed by the user with those
to be fitted.
The calling wrappers act in a similar way for the absorption parameter
(AV or NH) which may or may not be fitted.
- The UV form of FUNC offers the choice of fitting plain flux vs
wavelength (i.e. uses "raw" SPEUV2 at x(i)) or of fitting the
integral of the model over the data bin (by numeric integration of
SPEUV2 on 20 steps over the bin). In the latter the integral
is computed once for all bins, and the stored value returned at later calls.
Similarly the derivative FDER routines also computes derivatives
once at each parameter change, and returns the stored value at other calls.
- the X-ray form of FUNC also computes the model once at parameter
change only, using the stored values at other calls. However it performs
an important additional task, which is the convolution of the spectrum
with the response matrix for the I-th channel.
The FDER routine in the melib library is also used,
which computes derivatives only at fitted parameter change.
- One technical adjustment done quite recently to the X-ray fitting program
to speed them up a bit, is that, instead of reading the entire response
matrix (which in XAS is structured as an image), only the values pre-analysed
to be above a given threshold are saved internally and used for convolution.
- A further, more recent addition is the native support to cross-normalizations
when fitting data from more than one instrument. In XAS (unofficial add-on),
inspired in this by my Exosat software, spectra (and matrices)
from different instruments are physically joined in new files.
In the case of multi-instrument files, one can fit the cross-normalization of
all instruments (except one which is kept as reference). This is done saving
for each instrument :
For stepped parameters the index array (see below) p-th location contains -2,J to indicate
that the p-th stepped parameter is a cross-norm (-2) for the J-th instrument.
- a normalization value
- a normalization flag (0=not fitted 1=fitted 2=stepped)
- a normalization index (N means that it is fitted parameter A(N))
- a couple of arrays with the start and end channel of each instrument
The FUNC code takes care of this and recovers the normalizations of
individual spectral bins from the relevant fitted value (otherwise they are
initialized to the user-supplied fixed cross-normalization) and divides the
computed count spectrum by them.
- This is similar to but separate from what is done to decide whether a parameter
is fitted, stepped or fixed. There are the following arrays :
- a matrix (indexed by component and parameter) of values, which record
the user supplied value of the j-th parameter of the i-th component
- a matrix of flags which tell whether
the j-th parameter of the i-th component is fixed (0), fitted (1) or
- a matrix of indices (N means that the j-th parameter of the i-th component
is fitted parameter A(N))
There is a further index array for the stepped case. It has two couple of
elements (for the first and second stepped parameter). The p-th location
contains I,J to indicate that the p-th stepped parameter is the j-th parameter
of the i-th component.
The absorption parameter (AV or NH) is handled similarly,
being flagged as -1,-1 in the index array, while it does not use the matrix
of values, flags and indices but a separate scalar value, flag and index.