GAUFIT: least-square fits

Fitting with the least-squares algorithm proceeds much as before:

% gaufit

ESP GAUFIT running.

[...]

So far, the interaction with GAUFIT is exactly as before, but now we choose to use the new algorithm.

FWHM - Work in FWHM, rather than sigmas? /TRUE/ > 
LSQFIT - Use the non-linear least-squares method? /FALSE/ > y
CALCSD - Calculate and display uncertainties? /TRUE/ > 
BACK - Background count value ( < 0 to have it fitted) /760/ > -1
MAXITER - Maximum number of iterations (-1 for default) /-1/ > 

MODEL - Output NDF filename /@ic3374c-out/ > 
MODTYP - Whole image model (W)/Residuals (R)/reGression diag. (G) /'w'/ >

You choose to use the new algorithm by setting the parameter LSQFIT to true. Around half of the function-evaluations in the fitting algorithm are used to calculate the uncertainties in the fit parameters, so if you have no need for these (for some reason), you will save a significant amount of run time by not requesting them.

This algorithm can fit the background count, so you do not need to use another routine to determine this. If you do give a positive value for the BACK parameter, then the routine will use that instead of fitting it (there is absolutely no speed advantage to this).

Unless you suspect the algorithm is somehow misbehaving, you should not give a positive value for the iteration count: the default is 150, and the fit should converge well before that, so that this parameter merely acts as a check on any pathological cases which cause the routine to somehow run away with itself.

For the LSQ algorithm, there is a third option for MODTYP. MODTYP=G gives a `regression diagnostic', which is an image in which the value at each point is the change in the residual function if the corresponding point in the data were deleted. The residual function is half the sum of the squares of the differences between the model and the data. This is expensive to create, so you should not select it unless you actually wish to examine it.

First estimates of source data

Source      X           Y         Angle   FWHMa/as     FWHMb/as         Peak
  1       93.67       95.15       21.00     59.          57.        0.933E+03

   IT    NF     DRIFT    NL'D RESID
    0     1    0.00        0.00E+00
    1     2   0.110E-04    0.42E+03
    2     4   0.120E-02    0.41E+03
    3     5   0.383E-01    0.41E+03
    4     6   0.287        0.39E+03
    5     8   0.103        0.38E+03
    6    10   0.230        0.32E+03
    7    12   0.168        0.28E+03
    8    14   0.190        0.15E+03
    9    15   0.198        0.89E+02
   10    16   0.199        0.75E+02
   11    17   0.199        0.74E+02
   12    18   0.199        0.74E+02
   13    19   0.199        0.74E+02
   14    20   0.199        0.74E+02
   15    21   0.199        0.74E+02

GAUFIT2: algorithm performance
         Effective data s.d:     12.      Check reasonable
           Condition number:    0.21E+05  poor -- uncertainties plausible
        Optimisation metric:    0.91      Acceptable


Fitted parameter values:
Source      X           Y         Angle   FWHMa/as     FWHMb/as         Peak
  1       93.32       96.11       20.73     10.          6.1        0.189E+03

Parameter uncertainties: (-ve values indicate no estimate made)
Source      X           Y         Angle   FWHMa/as     FWHMb/as         Peak
  1        0.04        0.07        0.02    0.15         0.85E-01   -0.100E+01

Background (fitted) = 769.0236

OUT - Text file for parameter /@ic3374c-txt6/ >

There is a good deal of information in the text which GAUFIT produces while it is working.

For each iteration, GAUFIT prints the number of function evaluations so far, the `drift' from the initial guessed position, and a residual. The `drift' is a scaled estimate of how far the current solution $(x,y,\sigma_a,\sigma_b,\theta)$ is from the initial one, with the different components of the solution appropriately weighted. If the drift reaches 1, the routine will conclude that it has somehow got lost, and give up with an error message to that effect. The scale for this calculation is set when you `Indicate the outer limit of the source' at the time you pick the source positions - you might want to select a small circle if there are several sources which might become confused with each other. The `normalised residual' is related to how much the model deviates from the data, but the numerical value is unimportant - it should decrease as the calculation goes on.

When the routine calculates the parameter uncertainties, it assumes a value for the data standard deviation based on the size of the residual and the number of parameters and data points, and it reports this for your information. If this value seems unreasonable, because it is substantially different from an alternative estimate you have of the standard deviation, please let me know.

The routine also reports an estimated condition number (the ratio between the highest and lowest eigenvalues) for the Hessian used in the uncertainty calculation. If this is huge (more than about 10⁶), you should not place too much reliance on the uncertainties produced. In this case, the condition number is larger than we'd like, but the reported uncertainties should not be too far out. A condition number less than 10³ would suggest good reliable uncertainties.

In adapting the least-squares algorithm for this particular application, I made one or two optimisations, and the `optimisation metric' indicates how well this is performing. It's scaled so that a value towards 1 indicates that the optimisation is working as expected, zero means I needn't have bothered, and negative values suggest it's actually creating more work than the unimproved case. If you have a non-pathological case which gives negative values here, I'd like to hear about it.

Finally, the routine reports that it has converged successfully (if it hasn't, it should give a more-or-less useful explanation here).

The routine displays its results, then their corresponding uncertainties, with negative uncertainties indicating that no estimate was made. The routine does not (in this release) make an estimate of the peak flux uncertainty.