203 - Fits mit Scipy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Python Einführungskurs für das Physikalische Anfängerpraktikum der Universität Heidelberg | [Startseite](index.ipynb)\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 203 - Fits mit Scipy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- [Ein Modell an Daten fitten](#Ein-Modell-an-Daten-fitten)\n",
    "- [Aufgabe 1 - Temperaturverlauf in Heidelberg](#Aufgabe-1---Temperaturverlauf-in-Heidelberg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scipy ist ein sehr großes Modul mit vielen Funktionen, die für wissenschaftliche Analysen unverzichtbar sind. Dazu gehören:\n",
    "\n",
    "- Fits (`scipy.optimize`)\n",
    "- Statistik (`scipy.stats`)\n",
    "- Integration (`scipy.integrate`)\n",
    "- Interpolation (`scipy.interpolate`)\n",
    "\n",
    "Wir beschäftigen uns hier zunächst mit Fits. Die einfachste Funktion für Fits ist `scipy.optimize.curve_fit`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import curve_fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Numpy und Matplotlib brauchen wir immer\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ein Modell an Daten fitten\n",
    "\n",
    "Wir generieren zunächst einige Beispieldaten um dann einen Fit durchzuführen:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Zelle ausführen, um Beispieldaten zu generieren\n",
    "samples = np.random.normal(loc=0, scale=1, size=1000)\n",
    "N, bins, _ = plt.hist(samples, bins=50)\n",
    "x = bins[:-1] + (bins[1:] - bins[:-1]) / 2  # Mitte der Bins\n",
    "dN = np.sqrt(N)  # Fehler von N\n",
    "plt.errorbar(x, N, dN, ls='none')\n",
    "plt.xlabel(r'$x$')\n",
    "plt.ylabel(r'$N$')\n",
    "plt.title('Normalverteilte Beispieldaten')\n",
    "x, N, dN = x[N!=0], N[N!=0], dN[N!=0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Zu den Werten `x` haben wir nun normalverteilte Datenpunkte `N` mit Mittelwert $\\mu=0$ und Standardabweichung $\\sigma=1$. Als Fehler von N wurde $\\sqrt{N}$ angenommen.\n",
    "\n",
    "Diese Daten fitten wir nun an eine Gaußkurve. Wenn alles gut geht, sollten wir den Mittelwert und die Standardabweichung, mit denen die Daten generiert wurden, wieder herausfinden können.\n",
    "\n",
    "### 1. Fit-Funktion definieren\n",
    "\n",
    "Zuerst müssen wir die Fit-Funktion definieren:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import exp, pi, sqrt\n",
    "def gaussian(x, mu, sigma, A):\n",
    "    return A / (sigma * sqrt(2 * pi)) * exp(-(x - mu)**2 / sigma**2 / 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Hinweis:** Das erste Argument der Fit-Funktion muss die unabhängige Variable sein. Alle weiteren Argumente werden als Fit-Parameter interpretiert."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Fit durchführen\n",
    "\n",
    "Anschließend können wir die Fit-Funktion, und die Daten an `scipy.optimize.curve_fit` übergeben:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#curve_fit?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Die Funktion führt einen _least-squares_ Fit durch und gibt die Werte der Parameter zurück, wenn der Fit konvergiert. Als weitere Argumente können wir Startwerte `p0` der Parameter für den Fit sowie Fehler angeben:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "popt, pcov = curve_fit(gaussian, x, N, p0=[ 0, 1, 100 ], sigma=dN)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Der erste Rückgabewert **`popt`** ist eine Liste der gefundenen Werte der Parameter, welche die Daten am besten beschreiben. Sie sind in der Reihenfolge, wie sie in der Fit-Funktion definiert wurden:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(popt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Der zweite Rückgabewert **pcov** ist die Kovarianz-Matrix, welche die Varianzen $\\sigma^2$ der Parameter auf der Diagonalen enthält:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(pcov)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wir erhalten also für den Mittelwert und die Standardabweichung wie erwartet $\\mu=0$ und $\\sigma=1$ im Fehlerbereich:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"mu = {:.3f}+/-{:.3f}\".format(popt[0], np.sqrt(pcov[0][0])))\n",
    "print(\"sigma = {:.3f}+/-{:.3f}\".format(popt[1], np.sqrt(pcov[1][1])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Fit-Ergebnis plotten\n",
    "\n",
    "Die optimalen Werte der Fit-Parameter in `popt` können wir nun in die Fit-Funktion einsetzen und diese plotten:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Datenpunkte plotten\n",
    "plt.errorbar(x, N, dN, ls='none', marker='.', color='grey', alpha=0.5, label='Messwerte')\n",
    "# Fit plotten\n",
    "plt.plot(x, gaussian(x, *popt), label=\"\\n\".join([\"Gauß-Fit mit:\", r'$\\mu={:.3f}\\pm{:.3}$'.format(popt[0], np.sqrt(pcov[0][0])), r'$\\sigma={:.3f}\\pm{:.3f}$'.format(popt[1], np.sqrt(pcov[1][1]))]))\n",
    "plt.legend()\n",
    "plt.xlabel(r'$x$')\n",
    "plt.ylabel(r'$N$')\n",
    "plt.title('Gauß-Fit der Beispieldaten')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Hinweis:** Die Parameter in `popt` sind bereits in der richtigen Reihenfolge für die Fit-Funktion. Mit der folgenden Syntax können wir die Liste daher direkt _entpackt_ der Fit-Funktion übergeben:\n",
    ">\n",
    "> ```python\n",
    "> fit_function(x, *popt) # Entpackt die Liste `popt`\n",
    "> ```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Aufgabe 1 - Temperaturverlauf in Heidelberg\n",
    "\n",
    "Wir fitten nun die Temperaturdaten aus Heidelberg, die wir bereits zuvor verwendet haben. Wir können annehmen, dass sich der Jahres-Temperaturverlauf mit der Funktion $$T(t)=T_0+T_A\\frac{1-\\cos{\\!(2\\pi (t+\\Delta t))}}{2}$$ mit $t\\in [0,1]$ beschreiben lässt.\n",
    "\n",
    "Zuerst lesen wir die Daten ein und bestimmen wie zuvor den Zeitpunkt innerhalb des Jahres $t$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Daten einlesen\n",
    "data = np.loadtxt('data/temperatures.txt')\n",
    "date, T = data[np.abs(data[:,1]) != 99,:].transpose()\n",
    "t = date % 1 # Zeitpunkt innerhalb des Jahres"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "a) Definiere die Fit-Funktion `T_model(t, T_0, T_A, dt)` wie oben angegeben."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "203-1a-sol",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "from numpy import cos, pi\n",
    "### BEGIN SOLUTION\n",
    "T_model = lambda t, T_0, T_A, dt: T_0 + T_A * (1 - cos(2 * pi * (t + dt))) / 2\n",
    "### END SOLUTION"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "203-1a-test",
     "locked": true,
     "points": 1,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "from numpy.testing import assert_array_almost_equal\n",
    "try:\n",
    "    T_model\n",
    "except NameError:\n",
    "    raise NameError(\"Es gibt keine Funktion 'T_model'. Definiere eine Funktion mit diesem Namen.\")\n",
    "assert_array_almost_equal(T_model(np.array([0, 0.5, 1]), -1, 20, -0.04), [ -0.68583161,  18.68583161,  -0.68583161], 3, \"Deine Funktion liefert nicht die richtigen Werte. Überprüfe, ob sie mit der Aufgabenstellung übereinstimmt.\")\n",
    "print(\"Richtig.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "b) Verwende `curve_fit`, um den Fit durchzuführen. Weise drei Variablen `T_0`, `T_A` und `dt` jeweils den durch den Fit bestimmten Wert zu."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "203-1b-sol",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "popt, pcov = curve_fit(T_model, t, T, p0=[20, 0, 0])\n",
    "T_0, T_A, dt = popt[0], popt[1], popt[2]\n",
    "print(popt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy.testing import assert_almost_equal\n",
    "try:\n",
    "    T_0, T_A, dt\n",
    "except NameError:\n",
    "    raise NameError(\"Es gibt nicht alle der Variablen 'T_0', 'T_A' und 'dt'. Weise die Werte Variablen mit diesen Namen zu.\")\n",
    "assert_almost_equal(T_0, -0.91433813, 3, \"Hast du 'T_0' den richtigen Wert zugewiesen? Beachte die Reihenfolge der Parameter in 'popt'.\")\n",
    "assert_almost_equal(T_A, 19.91036717, 3, \"Hast du 'T_A' den richtigen Wert zugewiesen? Beachte die Reihenfolge der Parameter in 'popt'.\")\n",
    "assert_almost_equal(dt, -0.04026058, 3, \"Hast du 'dt' den richtigen Wert zugewiesen? Beachte die Reihenfolge der Parameter in 'popt'.\")\n",
    "print(\"Offenbar ist T_0={T_0:.2f}°C, T_A={T_A:.2f}°C und dt={dt:.1f} Tage. Sieht sinnvoll aus, oder?\".format(T_0=T_0, T_A=T_A, dt=dt * 365))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "c) Plotte das Ergebnis des Fits!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "203-1c-sol",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "plt.scatter(t, T, marker='.', color='gray', alpha=0.2, label='Messwerte')\n",
    "### BEGIN SOLUTION\n",
    "t_model = np.linspace(0, 1, 100)\n",
    "plt.plot(t_model, T_model(t_model, *popt), color='black', lw=2, label=\"\\n\".join([r'Fit $T(t)=T_0+T_A\\frac{-\\cos{\\!(2\\pi (t+\\Delta t))}+1}{2}$ mit:', r'$T_0={:.3f}\\pm{:.3f}'.format(popt[0], pcov[0][0]) + '^\\circ\\mathrm{C}$', r'$T_A={:.3f}\\pm{:.3f}'.format(popt[1], pcov[1][1]) + '^\\circ\\mathrm{C}$', r'$\\Delta t={:.3f}\\pm{:.3f}$'.format(popt[2], pcov[2][2])]))\n",
    "plt.title('Jahres-Temperaturverlauf in Heidelberg')\n",
    "plt.xlim(0, 1)\n",
    "plt.xlabel(r'Zeitpunkt innerhalb des Jahres $t$')\n",
    "plt.ylabel(r'Temperatur $T \\, [^\\circ{}\\mathrm{C}]$')\n",
    "plt.legend(loc='lower center', fontsize='medium')\n",
    "### END SOLUTION"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sieht dein Plot etwa so aus?\n",
    "\n",
    "![Jahres-Temperaturverlauf in Heidelberg mit Fit](plots/203-1c.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setze die Variable auf `True` wenn du mit deinem Plot zufrieden bist:\n",
    "plot_1c_fertig = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nose.tools import assert_true\n",
    "assert_true(plot_1c_fertig, \"Versuche dich an Aufgabe 1c, bis du mit deinem Plot zufrieden bist.\")\n",
    "print(\"😱 Ein Meisterwerk.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "Du kannst jetzt mit Python Daten analysieren, plotten und fitten 👍 🎉. Doch nur, indem du nun selbst damit arbeitest, wirst du die wissenschaftliche Programmierung in Python wirklich lernen. Im Physikalischen Anfängerpraktikum hast du nun die Gelegenheit dazu, sodass du anschließend Experte der wissenschaftlichen Programmierung sein wirst!\n",
    "\n",
    "Für die Auswertung der Praktikumsversuche wirst du Python und das Jupyter Notebook auf deinem eigenen Rechner ausführen müssen:\n",
    "\n",
    "- [>> Python installieren und ausführen](setup.ipynb)\n",
    "\n",
    "Wir hoffen, dass wir dir eine hilfreiche Einführung in die wissenschaftliche Programmierung mit Python bieten konnten. In den folgenden Lektionen findest du noch weitere Hinweise zur praktischen Anwendung im PAP:\n",
    "\n",
    "- [>> 301 - Workflows mit Daten im PAP](301%20-%20Workflows%20mit%20Daten%20im%20PAP.ipynb)\n",
    "\n",
    "Wenn du dein Können in der Programierung in Python unter Beweis stellen möchtest, versuche dich mal an folgenden Aufgaben:\n",
    "\n",
    "- [Challenges](challenges.ipynb)\n",
    "\n",
    "Mit deinem Feedback können wir diesen Kurs und die Verwendung von Python im Physikalischen Anfängerpraktikum für dich weiter verbessern. Wir bitten dich daher, kurz die folgenden Fragen zu beantworten:\n",
    "\n",
    "- [Feedback zum Online-Einführungskurs und Python im PAP](https://goo.gl/forms/nvuPvEOCP1CMrp5X2)\n",
    "\n",
    "Melde dich außerdem jederzeit bei [Nils Fischer](https://nilsleiffischer.de) (Konzeption und Implementierung des Kurses) oder [Jens Wagner](http://www.physi.uni-heidelberg.de/Einrichtungen/AP/) (Praktikumsleitung) mit Anmerkungen oder Vorschlägen. Der Kurs steht zudem als Open Source [auf GitHub](https://github.com/nilsleiffischer/python-course) zur Verfügung.\n",
    "\n",
    "_Viel Erfolg im Physikalischen Anfängerpraktikum!_\n",
    "\n",
    "[Startseite](index.ipynb)"
   ]
  }
 ],
 "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.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}