Fitting Models to Data

This module provides wrappers, called Fitters, around some Numpy and Scipy fitting functions. All Fitters can be called as functions. They take an instance of FittableModel as input and modify its parameters attribute. The idea is to make this extensible and allow users to easily add other fitters.

Linear fitting is done using Numpy’s numpy.linalg.lstsq function. There are currently two non-linear fitters which use scipy.optimize.leastsq and scipy.optimize.fmin_slsqp.

The rules for passing input to fitters are:

• Non-linear fitters currently work only with single models (not model sets).

• The linear fitter can fit a single input to multiple model sets creating multiple fitted models. This may require specifying the model_set_axis argument just as used when evaluating models; this may be required for the fitter to know how to broadcast the input data.

• The LinearLSQFitter currently works only with simple (not compound) models.

• The current fitters work only with models that have a single output (including bivariate functions such as Chebyshev2D but not compound models that map x, y -> x', y').

Simple 1-D model fitting

In this section, we look at a simple example of fitting a Gaussian to a simulated dataset. We use the Gaussian1D and Trapezoid1D models and the LevMarLSQFitter fitter to fit the data:

import numpy as np
import matplotlib.pyplot as plt
from astropy.modeling import models, fitting

# Generate fake data
np.random.seed(0)
x = np.linspace(-5., 5., 200)
y = 3 * np.exp(-0.5 * (x - 1.3)**2 / 0.8**2)
y += np.random.normal(0., 0.2, x.shape)

# Fit the data using a box model.
# Bounds are not really needed but included here to demonstrate usage.
t_init = models.Trapezoid1D(amplitude=1., x_0=0., width=1., slope=0.5,
                            bounds={"x_0": (-5., 5.)})
fit_t = fitting.LevMarLSQFitter()
t = fit_t(t_init, x, y)

# Fit the data using a Gaussian
g_init = models.Gaussian1D(amplitude=1., mean=0, stddev=1.)
fit_g = fitting.LevMarLSQFitter()
g = fit_g(g_init, x, y)

# Plot the data with the best-fit model
plt.figure(figsize=(8,5))
plt.plot(x, y, 'ko')
plt.plot(x, t(x), label='Trapezoid')
plt.plot(x, g(x), label='Gaussian')
plt.xlabel('Position')
plt.ylabel('Flux')
plt.legend(loc=2)

(png, svg, pdf)

As shown above, once instantiated, the fitter class can be used as a function that takes the initial model (t_init or g_init) and the data values (x and y), and returns a fitted model (t or g).

Simple 2-D model fitting

Similarly to the 1-D example, we can create a simulated 2-D data dataset, and fit a polynomial model to it. This could be used for example to fit the background in an image.

import warnings
import numpy as np
import matplotlib.pyplot as plt
from astropy.modeling import models, fitting

# Generate fake data
np.random.seed(0)
y, x = np.mgrid[:128, :128]
z = 2. * x ** 2 - 0.5 * x ** 2 + 1.5 * x * y - 1.
z += np.random.normal(0., 0.1, z.shape) * 50000.

# Fit the data using astropy.modeling
p_init = models.Polynomial2D(degree=2)
fit_p = fitting.LevMarLSQFitter()

with warnings.catch_warnings():
    # Ignore model linearity warning from the fitter
    warnings.simplefilter('ignore')
    p = fit_p(p_init, x, y, z)

# Plot the data with the best-fit model
plt.figure(figsize=(8, 2.5))
plt.subplot(1, 3, 1)
plt.imshow(z, origin='lower', interpolation='nearest', vmin=-1e4, vmax=5e4)
plt.title("Data")
plt.subplot(1, 3, 2)
plt.imshow(p(x, y), origin='lower', interpolation='nearest', vmin=-1e4,
           vmax=5e4)
plt.title("Model")
plt.subplot(1, 3, 3)
plt.imshow(z - p(x, y), origin='lower', interpolation='nearest', vmin=-1e4,
           vmax=5e4)
plt.title("Residual")

(png, svg, pdf)