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