{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Introduction to Scipy: Fitting data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We have talked about the Numpy and Matplotlib libraries, but there is a third library that is invaluable for Scientific Analysis: [Scipy](http://www.scipy.org). Scipy is basically a very large library of functions that you can use for scientific analysis. A good place to start to find out about the top-level scientific functionality in Scipy is the [Documentation](http://docs.scipy.org/doc/scipy/reference/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Examples of the functionality include:\n",
    "\n",
    "* Integration (scipy.integrate)\n",
    "* Optimization/Fitting (scipy.optimize)\n",
    "* Interpolation (scipy.interpolate)\n",
    "* Fourier Transforms (scipy.fftpack)\n",
    "* Signal Processing (scipy.signal)\n",
    "* Linear Algebra (scipy.linalg)   **[contains e.g. scipy.linalg.inv() as used in chapter 9 with np.linalg.inv()\n",
    "for systems of linear equation]**\n",
    "* Spatial data structures and algorithms (scipy.spatial)\n",
    "* Statistics (scipy.stats)\n",
    "* Multi-dimensional image processing (scipy.ndimage)\n",
    "\n",
    "and so on."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "In this section, we will take a look at how to fit models to data. When analyzing scientific data, fitting models to data allows us to determine the parameters of a physical system (assuming the model is correct)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "There are a number of routines in Scipy to help with fitting, but we will use the simplest one, ``curve_fit()``, which is imported as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import curve_fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The full documentation for the ``curve_fit()`` is available [here](http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html#scipy.optimize.curve_fit), and we will look at a simple example here, which involves fitting a straight line to a dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We first create a fake/mock dataset with some random noise:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.random.uniform(0., 100., 100)\n",
    "y = 3. * x + 2. + np.random.normal(0., 10., 100)\n",
    "plt.plot(x, y, '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's now imagine that this is real data, and we want to determine the slope and intercept of the best-fit line to the data. We start off by definining a function representing the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def line(x, a, b):\n",
    "    return a * x + b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The arguments to the function should be ``x``, followed by the parameters. We can now call ``curve_fit()`` to find the best-fit parameters using a least-squares fit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "popt, pcov = curve_fit(line, x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The ``curve_fit()`` function returns two items, which we call ``popt`` and ``pcov``. The ``popt`` argument are the best-fit paramters for ``a`` and ``b``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "popt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which is close to the initial values of ``3`` and ``2`` used in the definition of ``y``.\n",
    "\n",
    "The reason the values are not exact is because there are only a limited number of random samples, so the best-fit slope is not going to be exactly those used in the definition of ``y``. The ``pcov`` variable contains the *covariance* matrix, which indicates the uncertainties and correlations between parameters. This is mostly useful when the data has uncertainties."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Important Note:**  ``curve_fit()`` by default uses sigma just as weights, i.e., as hints how important a point should be in the least-squares computation.  To compute the covariance matrix, it then estimates the actual error essentially from the scatter in the data. If, however, your data points come with (reliable) error estimates, that actually throws away a lot of information, and the covariance matrix is less expressive than it could be.  To make ``curve_fit()`` treat its sigma argument as actual error estimates and use their absolute (rather than relative, as with weights) values in pcov calculation, one needs to pass ``absolute_sigma=True``. This is done here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's now try and fit the data assuming each point has a vertical error (standard deviation) of +/-10:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "e = np.repeat(10., 100)\n",
    "plt.errorbar(x, y, yerr=e, fmt=\"none\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "popt, pcov = curve_fit(line, x, y, sigma=e, absolute_sigma=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "popt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Now ``pcov`` will contain the true variance and covariance of the parameters, so that the best-fit parameters are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"a =\", popt[0], \"+/-\", pcov[0,0]**0.5)\n",
    "print(\"b =\", popt[1], \"+/-\", pcov[1,1]**0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We can now plot the best-fit line:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.errorbar(x, y, yerr=e, fmt=\"none\")\n",
    "xfine = np.linspace(0., 100., 100)  # define values to plot the function for\n",
    "plt.plot(xfine, line(xfine, popt[0], popt[1]), 'r-')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should now be able to fit simple models to datasets! Note that for more complex models, more sophisticated techniques may be required for fitting, but ``curve_fit()`` will be good enough for most simple cases."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Note that there is a way to simplify the call to the function with the best-fit parameters, which is:\n",
    "\n",
    "    line(x, *popt)\n",
    "\n",
    "The * notation will expand a list of values into the arguments of the function. This is useful if your function has more than one or two parameters. Hence, you can do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.errorbar(x, y, yerr=e, fmt=\"none\")\n",
    "plt.plot(xfine, line(xfine, *popt), 'r-')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the following code, we generate some random data points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.random.uniform(0., 10., 100)\n",
    "y = np.polyval([1, 2, -3], x) + np.random.normal(0., 10., 100)\n",
    "e = np.random.uniform(5, 10, 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit a line and a parabola to it and overplot the two models on top of the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As before, we use the [data/munich_temperatures_average_with_bad_data.txt](data/munich_temperatures_average_with_bad_data.txt) file, which gives the temperature in Munich every day for several years:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The following code reads in the file and removes bad values\n",
    "import numpy as np\n",
    "date, temperature = np.loadtxt('data/munich_temperatures_average_with_bad_data.txt', unpack=True)\n",
    "keep = np.abs(temperature) < 90\n",
    "date = date[keep]\n",
    "temperature = temperature[keep]\n",
    "plt.plot(date,temperature)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit the following function to the data:\n",
    "\n",
    "$$f(t) = a~\\cos{(2\\pi t + b)} + c$$\n",
    "\n",
    "where $t$ is the time in years.\n",
    "\n",
    "Make a plot of the data and the best-fit model in the range 2008 to 2012.\n",
    "\n",
    "What are the best-fit values of the parameters? What is the overall average temperature in Munich, and what are the typical daily average values predicted by the model for the coldest and hottest time of year?\n",
    "\n",
    "What is the meaning of the ``b`` parameter, and does its value make sense?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}