Models-Quick-Fit.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling 1: Make a quick fit using astropy.modeling\n",
    "\n",
    "## Authors\n",
    "Rocio Kiman, Lia Corrales, Zé Vinícius, Kelle Cruz, Stephanie T. Douglas\n",
    "\n",
    "## Learning Goals\n",
    "* Use `astroquery` to download data from Vizier\n",
    "* Use basic models in `astropy.modeling`\n",
    "* Learn common functions to fit\n",
    "* Generate a quick fit to data\n",
    "* Plot the model with the data\n",
    "* Compare different models and fitters\n",
    "\n",
    "## Keywords\n",
    "modeling, model fitting, astrostatistics, astroquery, Vizier, scipy, matplotlib, error bars, scatter plots\n",
    "\n",
    "## Summary\n",
    "In this tutorial, we will become familiar with the models available in [astropy.modeling](http://docs.astropy.org/en/stable/modeling/ ) and learn how to make a quick fit to our data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy.modeling import models, fitting\n",
    "from astroquery.vizier import Vizier\n",
    "import scipy.optimize\n",
    "# Make plots display in notebooks\n",
    "%matplotlib inline "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1) Fit a Linear model: Three steps to fit data using astropy.modeling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are going to start with a **linear fit to real data**. The data comes from the paper [Bhardwaj et al. 2017](https://ui.adsabs.harvard.edu/?#abs/2017A%26A...605A.100B). This is a catalog of **Type II Cepheids**, which is a type of **variable stars** that pulsate with a period between 1 and 50 days. In this part of the tutorial, we are going to measure the **Cepheids Period-Luminosity** relation using `astropy.modeling`. This relation states that if a star has a longer period, the luminosity we measure is higher.\n",
    "\n",
    "To get it, we are going to import it from [Vizier](http://vizier.u-strasbg.fr/viz-bin/VizieR) using [astroquery](http://astroquery.readthedocs.io/en/latest/vizier/vizier.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "catalog = Vizier.get_catalogs('J/A+A/605/A100')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This catalog has a lot of information, but for this tutorial we are going to work only with periods and magnitudes. Let's grab them using the keywords `'Period'` and `__Ksmag__`.  Note that `'e__Ksmag_'` refers to the error bars in the magnitude measurements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "period = np.array(catalog[0]['Period']) \n",
    "log_period = np.log10(period)\n",
    "k_mag = np.array(catalog[0]['__Ksmag_'])\n",
    "k_mag_err = np.array(catalog[0]['e__Ksmag_'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at the magnitude measurements as a function of period:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.errorbar(log_period, k_mag, k_mag_err, fmt='k.')\n",
    "plt.xlabel(r'$\\log_{10}$(Period [days])')\n",
    "plt.ylabel('Ks')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One could say that there is a linear relationship between log period and magnitudes. To probe it, we want to make a fit to the data. This is where `astropy.modeling` is useful. We are going to understand how in three simple lines we can make any fit we want. We are going to start with the linear fit, but first, let's understand what a model and a fitter are."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Models in Astropy\n",
    "[Models](http://docs.astropy.org/en/stable/modeling/#using-models) in Astropy are known parametrized functions. With this format they are easy to define and to use, given that we do not need to write the function expression every time we want to use a model, just the name. They can be linear or non-linear in the variables. Some examples of models are:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* [Gaussian1D](http://docs.astropy.org/en/stable/api/astropy.modeling.functional_models.Gaussian1D.html#astropy.modeling.functional_models.Gaussian1D)\n",
    "* [Trapezoid1D](http://docs.astropy.org/en/stable/api/astropy.modeling.functional_models.Trapezoid1D.html#astropy.modeling.functional_models.Trapezoid1D)\n",
    "* [Polynomial1D](http://docs.astropy.org/en/stable/api/astropy.modeling.polynomial.Polynomial1D.html#astropy.modeling.polynomial.Polynomial1D)\n",
    "* [Sine1D](http://docs.astropy.org/en/stable/api/astropy.modeling.functional_models.Sine1D.html#astropy.modeling.functional_models.Sine1D)\n",
    "* [Linear1D](http://docs.astropy.org/en/stable/api/astropy.modeling.functional_models.Linear1D.html#astropy.modeling.functional_models.Linear1D)\n",
    "* The [list](http://docs.astropy.org/en/stable/modeling/#module-astropy.modeling.functional_models) continues."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fitters in Astropy\n",
    "Fitters in Astropy are the classes resposable for making the fit. They can be linear or non-linear in the parameters (no the variable, like models). Some examples are:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* [LevMarLSQFitter()](http://docs.astropy.org/en/stable/api/astropy.modeling.fitting.LevMarLSQFitter.html#astropy.modeling.fitting.LevMarLSQFitter)       Levenberg-Marquardt algorithm and least squares statistic.\n",
    "* [LinearLSQFitter()](http://docs.astropy.org/en/stable/api/astropy.modeling.fitting.LinearLSQFitter.html#astropy.modeling.fitting.LinearLSQFitter)       A class performing a linear least square fitting.\n",
    "* [SLSQPLSQFitter()](http://docs.astropy.org/en/stable/api/astropy.modeling.fitting.SLSQPLSQFitter.html#astropy.modeling.fitting.SLSQPLSQFitter)        SLSQP optimization algorithm and least squares statistic.\n",
    "* [SimplexLSQFitter()](http://docs.astropy.org/en/stable/api/astropy.modeling.fitting.SimplexLSQFitter.html#astropy.modeling.fitting.SimplexLSQFitter)      Simplex algorithm and least squares statistic.\n",
    "* More detailles [here](http://docs.astropy.org/en/stable/modeling/#id21)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we continue with our fitting."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 1: Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we need to choose which model we are going to use to fit to our data. As we said before, our data looks like a linear relation, so we are going to use a linear model. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model = models.Linear1D()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 2: Fitter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Second we are going to choose the fitter we want to use. This choice is basically which method we want to use to fit the model to the data. In this case we are going to use the [Linear Least Square Fitting](https://www.mathworks.com/help/curvefit/least-squares-fitting.html). In the next exercise ([Modeling 2: Create a User Defined Model](http://learn.astropy.org/rst-tutorials/User-Defined-Model.html)) we are going to analyze how to choose the fitter. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fitter = fitting.LinearLSQFitter() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 3: Fit Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we give to our **fitter** (method to fit the data) the **model** and the **data** to perform the fit. Note that we are including weights: This means that values with higher error will have smaller weight (less importance) in the fit, and the contrary for data with smaller errors. This way of fitting is called *Weighted Linear Least Squares* and you can find more information about it [here](https://www.mathworks.com/help/curvefit/least-squares-fitting.html) or [here](https://en.wikipedia.org/wiki/Least_squares#Weighted_least_squares)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "best_fit = fitter(model, log_period, k_mag, weights=1.0/k_mag_err**2)\n",
    "print(best_fit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And that's it!\n",
    "\n",
    "We can evaluate the fit at our particular x axis by doing `best_fit(x)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.errorbar(log_period,k_mag,k_mag_err,fmt='k.')\n",
    "plt.plot(log_period, best_fit(log_period), color='g', linewidth=3)  \n",
    "plt.xlabel(r'$\\log_{10}$(Period [days])')\n",
    "plt.ylabel('Ks')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Conclusion:** Remember, you can fit data with three lines of code:\n",
    "\n",
    "1) Choose a [model](http://docs.astropy.org/en/stable/modeling/#module-astropy.modeling.functional_models).\n",
    "\n",
    "2) Choose a [fitter](http://docs.astropy.org/en/stable/modeling/#id21).\n",
    "\n",
    "3) Pass to the fitter the model and the data to perform fit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the model `Polynomial1D(degree=1)` to fit the same data and compare the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2) Fit a Polynomial model: Choose fitter wisely"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For our second example, let's fit a polynomial of degree more than 1. In this case, we are going to create fake data to make the fit. Note that we're adding gaussian noise to the data with the function `np.random.normal(0,2)` which gives a random number from a gaussian distribution with mean 0 and standard deviation 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "N = 100\n",
    "x1 = np.linspace(0, 4, N)  # Makes an array from 0 to 4 of N elements\n",
    "y1 = x1**3 - 6*x1**2 + 12*x1 - 9 \n",
    "# Now we add some noise to the data\n",
    "y1 += np.random.normal(0, 2, size=len(y1)) #One way to add random gaussian noise\n",
    "sigma = 1.5\n",
    "y1_err = np.ones(N)*sigma "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot it to see how it looks:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.errorbar(x1, y1, yerr=y1_err,fmt='k.')\n",
    "plt.xlabel('$x_1$')  \n",
    "plt.ylabel('$y_1$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To fit this data let's remember the three steps: model, fitter and perform fit. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model_poly = models.Polynomial1D(degree=3)\n",
    "fitter_poly = fitting.LinearLSQFitter() \n",
    "best_fit_poly = fitter_poly(model_poly, x1, y1, weights = 1.0/y1_err**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(best_fit_poly)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What would happend if we use a different fitter (method)? Let's use the same model but with `SimplexLSQFitter` as fitter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitter_poly_2 = fitting.SimplexLSQFitter()\n",
    "best_fit_poly_2 = fitter_poly_2(model_poly, x1, y1, weights = 1.0/y1_err**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(best_fit_poly_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we got a warning after using `SimplexLSQFitter` to fit the data. The first line says:\n",
    "\n",
    "`WARNING: Model is linear in parameters; consider using linear fitting methods. [astropy.modeling.fitting]`\n",
    "\n",
    "If we look at the model we chose: $y = c_0 + c_1\\times x + c_2\\times x^2 + c_3\\times x^3$, it is linear in the parameters $c_i$. The warning means that `SimplexLSQFitter` works better with models that are not linear in the parameters, and that we should use a linear fitter like `LinearLSQFitter`. The second line says:\n",
    "\n",
    "`WARNING: The fit may be unsuccessful; Maximum number of iterations reached. [astropy.modeling.optimizers]`\n",
    "\n",
    "So it's not surprising that the results are different, because this means that the fitter is not working properly. Let's discuss a method of choosing between fits and remember to **pay attention** when you choose the **fitter**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Compare results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One way to check which model parameters are a better fit is calculating the [Reduced Chi Square Value](https://en.wikipedia.org/wiki/Reduced_chi-squared_statistic). Let's define a function to do that because we're going to use it several times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def calc_reduced_chi_square(fit, x, y, yerr, N, n_free):\n",
    "    '''\n",
    "    fit (array) values for the fit\n",
    "    x,y,yerr (arrays) data\n",
    "    N total number of points\n",
    "    n_free number of parameters we are fitting\n",
    "    '''\n",
    "    return 1.0/(N-n_free)*sum(((fit - y)/yerr)**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reduced_chi_squared = calc_reduced_chi_square(best_fit_poly(x1), x1, y1, y1_err, N, 4)\n",
    "print('Reduced Chi Squared with LinearLSQFitter: {}'.format(reduced_chi_squared))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reduced_chi_squared = calc_reduced_chi_square(best_fit_poly_2(x1), x1, y1, y1_err, N, 4)\n",
    "print('Reduced Chi Squared with SimplexLSQFitter: {}'.format(reduced_chi_squared))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, the *Reduced Chi Square* for the first fit is closer to one, which means this fit is better. Note that this is what we expected after the discussion of the warnings.\n",
    "\n",
    "We can also compare the two fits visually:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.errorbar(x1, y1, yerr=y1_err,fmt='k.')\n",
    "plt.plot(x1, best_fit_poly(x1), color='r', linewidth=3, label='LinearLSQFitter()')  \n",
    "plt.plot(x1, best_fit_poly_2(x1), color='g', linewidth=3, label='SimplexLSQFitter()')\n",
    "plt.xlabel(r'$\\log_{10}$(Period [days])')\n",
    "plt.ylabel('Ks')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Results are as espected, the fit performed with the linear fitter is better than the second, non linear one. \n",
    "\n",
    "**Conclusion:** Pay attention when you choose the fitter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3) Fit a Gaussian: Let's compare to scipy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scipy has the function [scipy.optimize.curve_fit](https://docs.scipy.org/doc/scipy-1.0.0/reference/generated/scipy.optimize.curve_fit.html) to fit in a similar way that we are doing. Let's compare the two methods with fake data in the shape of a Gaussian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mu, sigma, amplitude = 0.0, 10.0, 10.0\n",
    "N2 = 100\n",
    "x2 = np.linspace(-30, 30, N)\n",
    "y2 = amplitude * np.exp(-(x2-mu)**2 / (2*sigma**2))\n",
    "y2 = np.array([y_point + np.random.normal(0, 1) for y_point in y2])   #Another way to add random gaussian noise\n",
    "sigma = 1\n",
    "y2_err = np.ones(N)*sigma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.errorbar(x2, y2, yerr=y2_err, fmt='k.')\n",
    "plt.xlabel('$x_2$')\n",
    "plt.ylabel('$y_2$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's do our three steps to make the fit we want. For this fit we're going to use a non-linear fitter, `LevMarLSQFitter`, because the model we need (`Gaussian1D`) is non-linear in the parameters. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model_gauss = models.Gaussian1D()\n",
    "fitter_gauss = fitting.LevMarLSQFitter()\n",
    "best_fit_gauss = fitter_gauss(model_gauss, x2, y2, weights=1/y2_err**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(best_fit_gauss)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can get the [covariance matrix](http://mathworld.wolfram.com/CovarianceMatrix.html) from `LevMarLSQFitter`, which provides an error for our fit parameters by doing `fitter.fit_info['param_cov']`. The elements in the diagonal of this matrix are the square of the errors. We can check the order of the parameters using:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_gauss.param_names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cov_diag = np.diag(fitter_gauss.fit_info['param_cov'])\n",
    "print(cov_diag)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print('Amplitude: {} +\\- {}'.format(best_fit_gauss.amplitude.value, np.sqrt(cov_diag[0])))\n",
    "print('Mean: {} +\\- {}'.format(best_fit_gauss.mean.value, np.sqrt(cov_diag[1])))\n",
    "print('Standard Deviation: {} +\\- {}'.format(best_fit_gauss.stddev.value, np.sqrt(cov_diag[2])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can apply the same method with `scipy.optimize.curve_fit`, and compare the results using again the *Reduced Chi Square Value*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def f(x,a,b,c):\n",
    "    return a * np.exp(-(x-b)**2/(2.0*c**2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "p_opt, p_cov = scipy.optimize.curve_fit(f,x2, y2, sigma=y1_err)\n",
    "a,b,c = p_opt\n",
    "best_fit_gauss_2 = f(x2,a,b,c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(p_opt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Amplitude: {} +\\- {}'.format(p_opt[0], np.sqrt(p_cov[0,0])))\n",
    "print('Mean: {} +\\- {}'.format(p_opt[1], np.sqrt(p_cov[1,1])))\n",
    "print('Standard Deviation: {} +\\- {}'.format(p_opt[2], np.sqrt(p_cov[2,2])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Compare results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reduced_chi_squared = calc_reduced_chi_square(best_fit_gauss(x2), x2, y2, y2_err, N2, 3)\n",
    "print('Reduced Chi Squared using astropy.modeling: {}'.format(reduced_chi_squared))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reduced_chi_squared = calc_reduced_chi_square(best_fit_gauss_2, x2, y2, y2_err, N2, 3)\n",
    "print('Reduced Chi Squared using scipy: {}'.format(reduced_chi_squared))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see there is a very small difference in the *Reduced Chi Squared*. This actually needed to happen, because the fitter in `astropy.modeling` uses scipy to fit. The advantage of using `astropy.modeling` is you only need to change the name of the fitter and the model to perform a completely different fit, while scipy require us to remember the expression of the function we wanted to use. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.errorbar(x2, y2, yerr=y2_err, fmt='k.')\n",
    "plt.plot(x2, best_fit_gauss(x2), 'g-', linewidth=6, label='astropy.modeling')\n",
    "plt.plot(x2, best_fit_gauss_2, 'r-', linewidth=2, label='scipy')\n",
    "plt.xlabel('$x_2$')\n",
    "plt.ylabel('$y_2$')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Conclusion:** Choose the method most convenient for every case you need to fit. We recomend `astropy.modeling` because is easier to write the name of the function you want to fit than to remember the expression every time we want to use it. Also, `astropy.modeling` becomes useful with more complicated models like [two gaussians](http://docs.astropy.org/en/stable/modeling/#compound-models) plus a [black body](http://docs.astropy.org/en/stable/modeling/#blackbody-radiation), but that is another tutorial."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary:\n",
    "\n",
    "Let's review the conclusion we got in this tutorial:\n",
    "\n",
    "1. You can fit data with **three lines of code**:\n",
    "    * model\n",
    "    * fitter\n",
    "    * perform fit to data\n",
    "    \n",
    "    \n",
    "2. **Pay attention** when you choose the **fitter**.\n",
    "\n",
    "3. Choose the method most convenient for every case you need to fit. We recomend `astropy.modeling` to make **quick fits of known functions**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4) Exercise: Your turn to choose"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the next data:\n",
    " * Choose model and fitter to fit this data\n",
    " * Compare different options"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "N3 = 100\n",
    "x3 = np.linspace(0, 3, N3)\n",
    "y3 = 5.0 * np.sin(2 * np.pi * x3)\n",
    "y3 = np.array([y_point + np.random.normal(0, 1) for y_point in y3])\n",
    "sigma = 1.5\n",
    "y3_err = np.ones(N)*sigma "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plt.errorbar(x3, y3, yerr=y3_err, fmt='k.')\n",
    "plt.xlabel('$x_3$')\n",
    "plt.ylabel('$y_3$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}