1********************** 2Fitting Models to Data 3********************** 4 5This module provides wrappers, called Fitters, around some Numpy and Scipy 6fitting functions. All Fitters can be called as functions. They take an 7instance of `~astropy.modeling.FittableModel` as input and modify its 8``parameters`` attribute. The idea is to make this extensible and allow 9users to easily add other fitters. 10 11Linear fitting is done using Numpy's `numpy.linalg.lstsq` function. There are 12currently two non-linear fitters which use `scipy.optimize.leastsq` and 13`scipy.optimize.fmin_slsqp`. 14 15The rules for passing input to fitters are: 16 17* Non-linear fitters currently work only with single models (not model sets). 18 19* The linear fitter can fit a single input to multiple model sets creating 20 multiple fitted models. This may require specifying the ``model_set_axis`` 21 argument just as used when evaluating models; this may be required for the 22 fitter to know how to broadcast the input data. 23 24* The `~astropy.modeling.fitting.LinearLSQFitter` currently works only with 25 simple (not compound) models. 26 27* The current fitters work only with models that have a single output 28 (including bivariate functions such as 29 `~astropy.modeling.polynomial.Chebyshev2D` but not compound models that map 30 ``x, y -> x', y'``). 31 32.. _modeling-getting-started-1d-fitting: 33 34Simple 1-D model fitting 35------------------------ 36 37In this section, we look at a simple example of fitting a Gaussian to a 38simulated dataset. We use the `~astropy.modeling.functional_models.Gaussian1D` 39and `~astropy.modeling.functional_models.Trapezoid1D` models and the 40`~astropy.modeling.fitting.LevMarLSQFitter` fitter to fit the data: 41 42.. plot:: 43 :include-source: 44 45 import numpy as np 46 import matplotlib.pyplot as plt 47 from astropy.modeling import models, fitting 48 49 # Generate fake data 50 np.random.seed(0) 51 x = np.linspace(-5., 5., 200) 52 y = 3 * np.exp(-0.5 * (x - 1.3)**2 / 0.8**2) 53 y += np.random.normal(0., 0.2, x.shape) 54 55 # Fit the data using a box model. 56 # Bounds are not really needed but included here to demonstrate usage. 57 t_init = models.Trapezoid1D(amplitude=1., x_0=0., width=1., slope=0.5, 58 bounds={"x_0": (-5., 5.)}) 59 fit_t = fitting.LevMarLSQFitter() 60 t = fit_t(t_init, x, y) 61 62 # Fit the data using a Gaussian 63 g_init = models.Gaussian1D(amplitude=1., mean=0, stddev=1.) 64 fit_g = fitting.LevMarLSQFitter() 65 g = fit_g(g_init, x, y) 66 67 # Plot the data with the best-fit model 68 plt.figure(figsize=(8,5)) 69 plt.plot(x, y, 'ko') 70 plt.plot(x, t(x), label='Trapezoid') 71 plt.plot(x, g(x), label='Gaussian') 72 plt.xlabel('Position') 73 plt.ylabel('Flux') 74 plt.legend(loc=2) 75 76As shown above, once instantiated, the fitter class can be used as a function 77that takes the initial model (``t_init`` or ``g_init``) and the data values 78(``x`` and ``y``), and returns a fitted model (``t`` or ``g``). 79 80.. _modeling-getting-started-2d-fitting: 81 82Simple 2-D model fitting 83------------------------ 84 85Similarly to the 1-D example, we can create a simulated 2-D data dataset, and 86fit a polynomial model to it. This could be used for example to fit the 87background in an image. 88 89.. plot:: 90 :include-source: 91 92 import warnings 93 import numpy as np 94 import matplotlib.pyplot as plt 95 from astropy.modeling import models, fitting 96 97 # Generate fake data 98 np.random.seed(0) 99 y, x = np.mgrid[:128, :128] 100 z = 2. * x ** 2 - 0.5 * x ** 2 + 1.5 * x * y - 1. 101 z += np.random.normal(0., 0.1, z.shape) * 50000. 102 103 # Fit the data using astropy.modeling 104 p_init = models.Polynomial2D(degree=2) 105 fit_p = fitting.LevMarLSQFitter() 106 107 with warnings.catch_warnings(): 108 # Ignore model linearity warning from the fitter 109 warnings.simplefilter('ignore') 110 p = fit_p(p_init, x, y, z) 111 112 # Plot the data with the best-fit model 113 plt.figure(figsize=(8, 2.5)) 114 plt.subplot(1, 3, 1) 115 plt.imshow(z, origin='lower', interpolation='nearest', vmin=-1e4, vmax=5e4) 116 plt.title("Data") 117 plt.subplot(1, 3, 2) 118 plt.imshow(p(x, y), origin='lower', interpolation='nearest', vmin=-1e4, 119 vmax=5e4) 120 plt.title("Model") 121 plt.subplot(1, 3, 3) 122 plt.imshow(z - p(x, y), origin='lower', interpolation='nearest', vmin=-1e4, 123 vmax=5e4) 124 plt.title("Residual") 125