Concepts in my fitting software

CURFIT basics

• 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(x_i) 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 point x=w.

  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 CURFIT itself.

  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 the following:

  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 y=f(x).

• 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 Avni 1976)

  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 chisquare 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.

  Spectral fitting

• 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" representations :
  • Flambda vs lambda : erg/cm²/s/Å vs Å
  • dN/dE vs E : photon/cm²/s/keV vs keV
  • Fnu vs nu : erg/cm²/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 85 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 spectral model
  • 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 matrix)

• 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 SPESUM).

• 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 (A_V or N_H) 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 :
  • 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
  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.

  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 stepped (2)
  • 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 (A_V or N_H) 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.