Fitting an Image (Choosing Fit Statistics and Solvers)

In the example code below, we assume that an instance of the fitting.Imfit class has already been created and supplied with the necessary ModelDescription.

For example,

imfit_fitter = pyimfit.Imfit(someModelDescription)

Finding the best-fit solution by minimizing a fit statistic

Fitting a model to an image in PyImfit involves calculating a model image and then comparing it pixel-by-pixel with a data image to derive a summary “fit statistic” (based on some total likelihood value). The goal is to minimize the fit statistic (which corresponds to maximizing the likelihood). This is done by iteratively adjusting the parameters of the model and recomputing the fit statistic until convergence is achieved; the algorithm which oversees this process is called a “minimizer” or “solver”.

Fit statistics (chi^2 and all that)

Which fit statistic to use depends in part on your data source, and so is determined as part of the loadData method of the Imfit class.

  1. Chi^2 Statistics
    1. Sigmas estimated from the data values
    2. Sigmas estimated from the model values
    3. Sigmas from a pre-existing error/uncertainty/variance image

To specify model-based chi^2 as the fit statistic

imfit_fitter.loadData(imageData, ..., use_model_for_errors=True)

To supply your own error/uncertainty map

imfit_fitter.loadData(imageData, error=errorImageData, ...)
imfit_fitter.loadData(imageData, error=errorImageData, error_type="variance", ...)
  1. Pure Poisson Statistics

To specify Poisson-MLR as the fit statistic

imfit_fitter.loadData(imageData, ..., use_poisson_mlr=True)
  1. Possible specification of image A/D gain, exposure time, etc.

For any case except using a pre-existing error map, you may need to supply information about how the values in the data image can be converted to detected counts (e.g., detected photoelectrons), since the underlying statistical models assume the latter. For example, if the per-pixel values were converted to ADUs via an A/D gain, you should supply the gain value (in electrons/ADU); if the values are counts/second, you should also supply the total intgration time. If there was a significant read noise term, this should also be described. The relevant keywords for the loadData and fit methods are: gain (A/D gain in electrons/ADU), read_noise (Gaussian read noise in electrons), exp_time (seconds, if the data values are ADU/sec), and n_combined (number of combined exposures). An example:

imfit_fitter.loadData(imageData, ..., gain=3.1, exp_time=800, read_noise=7.5)

Note: If you are using Poisson-MLR as the fit statistic, then read_noise should not be used (the Poisson MLR statistical model cannot handle a Gaussian read-noise term).

  1. Possible pre-subtracted background level

In some cases, it may be convenient to work with data images where the sky background has been removed. The fitting process needs to know about this, since otherwise there will be problems with data pixels having values near or below zero. You can specify a constant background level that has already been subtracted from the image, using the original_sky keyword for the loadData and fit methods; the value should be in the same units as the data pixels (e.g., ADU, ADU/sec, etc.). An example:

imfit_fitter.loadData(imageData, ..., original_sky=244.9)

Minimizers/Solvers

To actually find the best fit, you tell the Imfit object to find the minimum fit-statistic value, using a particular minimization algorithm.

PyImfit has three minimizers (a.k.a. solvers):

  • Levenberg-Marquardt: This is the default, and is by far the fastest; it has the drawback of being the most prone to being trapped in local minima in the fit-statistic landscape. It requires initial guesses for each parameter value.

    imfit_fitter.DoFit()
imfit_fitter.DoFit(solver="LM")

  • Nelder-Mead Simplex: This is a slower algorithm, generally held to be less likely to be trapped in local fit-statistic minima. Like Levenberg-Marquardt, it requires initial guesses for each parameter value.

imfit_fitter.DoFit(solver="NM")
  • Differential Evolution: This is a genetic-algorithms approach, and is the slowest of all the algorithms. Unlike the other two solvers, it requires lower and upper parameter limits for all non-fixed parameters. Initial guesses for the non-fixed parameter values are ignored.
imfit_fitter.DoFit(solver="DE")

Roughly speaking, the Nelder-Mead simplex minimizer is about an order of magnitude slower than Levenberg-Marquardt, and Differential Evolution is itself about an order of magnitude slower than Nelder-Mead simplex.

More Information

For more information: