********************** 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 `~astropy.modeling.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 `~astropy.modeling.fitting.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 `~astropy.modeling.polynomial.Chebyshev2D` but not compound models that map ``x, y -> x', y'``). .. _modeling-getting-started-1d-fitting: 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 `~astropy.modeling.functional_models.Gaussian1D` and `~astropy.modeling.functional_models.Trapezoid1D` models and the `~astropy.modeling.fitting.LevMarLSQFitter` fitter to fit the data: .. plot:: :include-source: 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) 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``). .. _modeling-getting-started-2d-fitting: 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. .. plot:: :include-source: 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")