JupyterHub Screenshot aktualisiert

Einfuehrung_iminuit.ipynb

@@ -0,0 +1,594 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "c9a4045f-f389-40f5-9f19-32f8bbebc75d",
+   "metadata": {},
+   "source": [
+    "# Methode der kleinsten Quadrate\n",
+    "\n",
+    "Im folgenden wollen wir die **Methode der kleinsten Quadrate (Least Squares)** näher beleuchten. Diese Methode wird oft benutzt, um eine Funktion $\\lambda(x; \\ $**$\\phi$**$)$ mit den Funktionsparametern $\\mathbf{\\phi}$ an die gemessenen Punkte **$(x,y)$** anzupassen. Um jedoch die **Methode der kleinsten Quadrate** zu verstehen, wollen wir sie erst einmal anschaulich und halb-mathematisch herleiten. Dabei stüzen wir uns im Folgenden auf eine Herleitung aus dem Buch **\"Statistical Data Analysis\"**  von **Glen Cowan**."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5a8f4ef9-222c-440f-8621-e612a8988fd0",
+   "metadata": {},
+   "source": [
+    "Bevor wir dies jedoch tun, schauen wir uns das Problem des Fittens doch erst einmal anschaulich an. \n",
+    "\n",
+    "<figure class=\"image\">\n",
+    "<img src=\"images/SketchLeastSquares.png\"  alt=\"{{ Beispiel PDF }}\" width=70%>\n",
+    "</figure>\n",
+    "\n",
+    "Beim Fitten, zum Beispiel einer Geraden (lila) an eine Reihe von Messpunkten (schwarz), wollen wir den Abstand zwischen den einzelnen Messpunkten und der Geraden (orange) möglichst klein halten. Sprich die Summe über alle $\\Delta Y_i$ \n",
+    "\n",
+    "$$\\sum_i \\Delta Y_i $$\n",
+    "\n",
+    "sollte möglichst klein sein, wobei $\\Delta Y_i$ durch \n",
+    "\n",
+    "$$ \\Delta Y_i = y_i – f(x_i, \\vec{\\theta})$$\n",
+    "\n",
+    "gegeben ist und $f(x, \\vec{\\theta})$ unsere Fitfunktion repräsentiert. Hierbei symbolisiert $\\vec{\\theta}$ die Parameter unserer Funktion. Sprich im Fall einer Geraden die Steigung $m$ und den Offset $y_0$ ($\\vec{\\theta}=(m, y_0)$). \n",
+    "\n",
+    "Darüber hinaus sollte die Richtung des Abstandes, sprich ob ein Messpunkt unterhalb oder oberhalb der Fitfunktion liegt, keine Rolle spielen. Daher quadrieren wir das Ganze und erhalten somit\n",
+    "\n",
+    "$$ LS = \\sum_i = (y_i – f(x_i, \\theta))^2$$\n",
+    "\n",
+    "Dies ist die allgemeinste Form der Methode der kleinsten Quadrate. Sie besagt, dass die Funktion, welche die Messpunkte am besten beschreibt, sprich die optimalen Werte für $\\vec{\\theta}$ aufweist, den Ausdruck LS minimiert. \n",
+    "\n",
+    "Nun weisen unsere Messpunkte nicht nur Werte für X und Y aus, sondern sind noch zusätzlich durch einen Messunsicherheit (Messfehler) charakterisiert. Diese sollten wir natürlich bei der Bestimmung unserer Parameter $\\vec{\\theta}$ berücksichtigen. Sprich Messwerte mit einer großen Unsicherheit sollten weniger stark berücksichtigt werden wie Messwerte mit einer kleinen Unsicherheit. Dies können wir gewährleisten, sofern wir die Distanzen $\\Delta Y_i$ mit den jeweiligen Unsicherheiten $\\Delta y_i$ gewichten, sprich  \n",
+    "\n",
+    "$$ \\chi^2 = \\sum_i =\\frac{(y_i – f(x_i, \\theta))^2}{\\Delta y_i^2}$$\n",
+    "\n",
+    "berechnen. Das Quadrieren der Unsicherheiten sorgt dafür, dass der Ausdruck dimensionslos wird. Diese besondere Form der kleinsten Quadrate nennt man auch oft $\\chi^2$-Fit. Wir werden später noch einmal genauer beleuchten warum. "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1153a474-8afe-44ae-8511-b403a4ad861d",
+   "metadata": {},
+   "source": [
+    "Nun wollen wir uns erst einmal ein Beispiel anschauen, wie dies in der Praxis aussieht. In der nachfolgenden Animation wird ein Ohm’schwer Widerstand an eine Reihe von Spannungs- und Strommessungen gefittet. Dies entspricht unserem obigen Geradenbeispiel. \n",
+    "<figure class=\"image\">\n",
+    "<img src=\"images/MaterialPythonkurs092018/Fitting_gif.gif\"  alt=\"{{ Least Square Beispiel }}\" width=100%>\n",
+    "</figure>\n",
+    "\n",
+    "Wie die Animation zeigt, werden so lange verschiedene Widerstände ausprobiert, bis ein Wert gefunden wurde, bei dem das $\\chi^2$ minimal wird. Dieses Variieren der Widerstandswerte passiert nicht zufällig, sondern basiert auf einem Algorithmus für Optimierungsverfahren. \n",
+    "\n",
+    "Es gibt verschiedene Arten von Algorithmen, um Minimierungsprobleme zu lösen. Wie diese genau aufgebaut sind, lernen Sie in anderen Programmierveranstaltungen, wie zum Beispiel *Programmieren für Physiker* oder *Computer in der Wissenschaft*. Zum Glück haben uns in Python bereits andere Menschen diese Arbeit abgenommen. Im folgenden wollen wir uns das package `iminuit` etwas genauer anschauen, welches bereits ein sehr umfangreiches und mächtiges Fittingtool darstellt. \n",
+    "\n",
+    "[iminuit](https://iminuit.readthedocs.io/en/stable/tutorials.html) verfügt auch über eine exzellente Dokumentation, mit Hilfe derer Sie auch komplexere Probleme lösen können."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "702b9a12-539d-43b5-8314-7f8dcdf93cda",
+   "metadata": {},
+   "source": [
+    "Um mit Hilfe von `imnuit` etwas zu fitten brauchen wir zunächst einmal ein paar Messdaten und ein Fitmodel. Im Folgenden wollen wir die Entladekurve eines Kondensators mit der Kapazität $C$ über einen Widerstand $R$ bestimmen. Die Entladekurve ist durch eine einfache Exponentialfunktion der Form \n",
+    "\n",
+    "$$ I = I_0 \\exp\\{-t/RC\\}$$\n",
+    "\n",
+    "gegeben. Die Messdaten befinden sich in einer CSV-Datei im Ordner `data`. Die CSV-Datei kann mit Hilfe des `pandas` package eingelesen werden. [pandas](https://pandas.pydata.org/) ist ähnlich wie `numpy` ein package welches eine Fülle an Funktionen zum Verarbeiten und Verwalten von Daten bereitstellt. Es gehört ähnlich wie auch `numpy`, `scipy` und `matplotlib` zu den Standardbibliotheken, welche sehr häufig in der Wissenschaft verwendet werden. Aufgrund der zeitlichen Limitierung des Versuchstages können wir leider nicht auf alle Funktionen von `pandas` eingehen und wollen uns im Folgenden lediglich auf die Grundlagen beschränken. Für ihre zukünftigen Praktika lohnt es sich jedoch, noch mehr über `pandas` in Ihrer Eigenstudienzeit zu lernen."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5215840a-1276-49c1-9181-274cd8a2b4bf",
+   "metadata": {},
+   "source": [
+    "CSV-Datein können wie folgt eingelesen werden"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 485,
+   "id": "f8ef1be0-a42d-4a11-b674-c2ed099fefcb",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import pandas as pd\n",
+    "\n",
+    "data_frame = pd.read_csv('data/discharge_data.csv')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "15800aa8-8a7f-4d59-ab06-3edc6bb1e443",
+   "metadata": {},
+   "source": [
+    "Dabei gibt pandas die Daten als so genannten DataFrames zurück. Dies sind Objekte, welche ähnlich wie strukturierte `numpy.arrays` zu behandeln sind. DataFrames werden allgemein als Tabellen dargestellt."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f14ca80f-e0d7-4447-9335-b3744f7a028f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "data_frame"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "fc24d5fa-d3c0-4866-b18a-9dd07768a222",
+   "metadata": {},
+   "source": [
+    "Um die Daten aus einer bestimmente Spalte zu bekommen, können diese einfach mit dem Spaltennamen aufgerufen werden:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "93b7cbb1-1095-4a53-83d9-7b32f068daea",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "data_frame['time']"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4476302c-eb65-409a-b1aa-2342ecbd9c88",
+   "metadata": {},
+   "source": [
+    "oder"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "969d8afa-5d52-4e01-8b64-ddab090891b8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "data_frame['current']"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "32249263-ec9f-44de-81b7-7a6c69e23332",
+   "metadata": {},
+   "source": [
+    "Einzelne Messwerte lassen sich mit Hilfe von `.loc` bestimmen."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e4b44637-8e25-46c1-863d-3cd7604f52dd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "data_frame.loc[5, 'current']"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a2c0c04b-be37-482d-aabc-802bfa2965d2",
+   "metadata": {},
+   "source": [
+    "Sollten Sie eine Spalte von Messdaten in ein `numpy.array` umwandeln wollen, so können Sie dies wie folgt erreichen"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f246f55e-5fc8-427c-990d-3e97799b5aeb",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "data_frame['current'].values"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3d2bd9ed-852d-4448-a051-a6f677ea891d",
+   "metadata": {},
+   "source": [
+    "Die Messdaten können Sie auch wie gewohnt mit Hilfe von `matplotlib` darstellen."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e3898686-3926-48a0-be4c-4d460a1792f3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "plt.errorbar(\n",
+    "    data_frame['time'], \n",
+    "    data_frame['current'], \n",
+    "    data_frame['delta_current'], \n",
+    "    ls='', \n",
+    "    marker='.', \n",
+    "    color='k'\n",
+    ")\n",
+    "plt.xlabel('Time []')\n",
+    "plt.ylabel('Current []')\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "100a4fe4-a5c4-4be3-a7f7-13337b97a194",
+   "metadata": {},
+   "source": [
+    "Nun wollen wir die Messdaten mit Hilfe von `iminuit` fitten. Hierzu müssen wir zunächste zwei Module des packages importieren und eine Funktion für die Entladekurve des Kondensators definieren:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "520f4973",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Diese Zelle nur auf JupyterHub des ZDV ausführen um `iminuit` zu installieren!\n",
+    "# import sys\n",
+    "# import subprocess\n",
+    "# subprocess.check_call([\n",
+    "#     sys.executable, \n",
+    "#     '-m',\n",
+    "#     'pip',\n",
+    "#     'install',\n",
+    "#     '--proxy',\n",
+    "#     'http://webproxy.zdv.uni-mainz.de:3128',\n",
+    "#     'iminuit'\n",
+    "# ])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 492,
+   "id": "2ffe340b-cd0f-45ec-b5b8-42e7a0349d4c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "from iminuit import Minuit, cost\n",
+    "import numpy as np\n",
+    "\n",
+    "def discharge_current(t, I0, R, C):\n",
+    "    return I0 * np.exp(-t/(R*C))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ef87da8f-7af9-4e3f-af63-c28a2e1d9830",
+   "metadata": {},
+   "source": [
+    "Nun können wir den Fit selbst durchführen. Hierzu muss zuerst mittels dem `cost` Modul eine sogenannte Kostenfunktion erstellt werden. Die Kostenfunktion ist im Grunde unsere $\\chi^2$ Funktion\n",
+    "\n",
+    "$$ \\chi^2 = \\sum_i =\\frac{(y_i – f(x_i, \\theta))^2}{\\Delta y_i^2}$$\n",
+    "\n",
+    "welche minimiert werden soll. Dies ist bereits bei `iminuit` für uns vordefiniert. Anschließend können wir die genutzt Kostenfunktion über `Minuit` minimieren lassen. Hierzu müssen wir zunächst geeignete Startwerte für den Minimierungsprozess vorgeben. Diese sollten im Idealfall nicht allzu weit von den wahren Werten entfernt liegen. Wir werden an einem späteren Beispiel noch einmal genauer zeigen, wie man hier vorgehen kann. Um den Minimierungsprozess zu starten muss noch am Ende `migrad()` aufgerufen werden."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "bf36b7b9-fb20-47b7-8538-479026b48fb2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Kostenfunktion:\n",
+    "ls = cost.LeastSquares(\n",
+    "    data_frame['time'],\n",
+    "    data_frame['current'], \n",
+    "    data_frame['delta_current'], \n",
+    "    discharge_current,\n",
+    ")\n",
+    "\n",
+    "# Minimierung\n",
+    "mi = Minuit(ls,        # Kostenfunktion \n",
+    "            I0=0.9,    # Startwerte\n",
+    "            R=10*10**3, \n",
+    "            C=10**-6\n",
+    "           )\n",
+    "mi.migrad()\n",
+    "mi.hesse()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b5df2d60-8284-4757-96c8-7f26afc17942",
+   "metadata": {},
+   "source": [
+    "Wie ihr seht gibt euch minuit euch vier verschiedene Objekte zurück. Für euch am wichtigsten ist die erste Tabelle, welche euch zeigt, ob euer Fit funktioniert hat. Im Allgemeinen gilt sind hier alle Felder grün hat euer Fit funktioniert, gelbe Felder können ein Problem andeuten müssen sie aber nicht und lila Felder bedeuten, dass etwas mit eurem Fit nicht in Ordnung ist. Die Bedeutungen der einzelnen Felder für unseren obigen Fit sind auch nochmal in der nachfolgenden Abbildung einzeln erklärt. Die Bedeutung der meisten Felder werden wir noch im laufe des Kurses kennen lernen. \n",
+    "\n",
+    "<figure class=\"image\">\n",
+    "<img src=\"images/FitPerformance.png\"  alt=\"{{Fit Performance }}\" width=100%>\n",
+    "</figure>\n",
+    "\n",
+    "Wie wir unserer Tabelle entnehmen können, gibt es also ein Problem mit unserem Fit. Um besser verstehen zu können, was das Problem sein könnte, wollen wir uns auch noch die anderen Outputs ansehen.\n",
+    "\n",
+    "Die zweite Tabelle zeigt uns die bestimmten Werte für die Parameter in der Spalte `Value` und deren Unsicherheiten in der Spalte `Hess error`. Hierbei fällt auf, dass für unseren obigen Fit die Unsicherheiten der Parameter $R$ und $C$ größer sind als die bestimmten Werte selbst. \n",
+    "\n",
+    "Die dritte Tabelle ist die sogennnante **Kovarianzmatrix**. Die Kovarianzmatrix hat als Einträge auf ihrer **Hauptdiagonalen** die **Varianzen der entsprechenden Parameter** auf der **Nebendiagonalen** stehen die **Kovarianzen**. Der Wert in Klammern gibt die **Korrelation** zwischen den entspechenden Parametern an. Sind zwei Parameter stark **korreliert**, wird das entsprechende Feld **blau** dargestellt, bei einer **Antikorrelation** ist das Feld **rot**. \n",
+    "\n",
+    "Die letzte Ausgabe ist ein Plot unserer Messwerte zusammen mit der Fitfunktion basierend auf den Parametern des besten Fits. (Nur für neuere Version von `iminuit`)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "72665daa-1d74-41da-8b9a-1e4c427eed07",
+   "metadata": {},
+   "source": [
+    "Obwohl underser Fit unsere Messdaten gut widerspiegelt, scheint es ein Problem mit der Bestimmung einiger Parameter zu geben. Die große Unsicherheit in $R$ und $C$ deutet an, dass hier das Problem liegt. Um dies zu bestätigen, können wir uns einmal das reduzierte $\\chi^2(x, I_0, R, C)$ als Funktion des entsprechenden Parameters von `iminuit` plotten lassen, während wir die anderen Parameter, so wie die x-Werte, konstant lassen. \n",
+    "\n",
+    "Für $I_0$ sieht das entsprechende Profil so aus:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d3230cb6-fbe3-4093-ba09-5271dc168a4d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mi.draw_profile('I0')\n",
+    "plt.ylabel('$\\chi^2(I_0, x, R, C)/ndof$')\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b837e542-d3c9-4f61-a8d1-4f22db7d5137",
+   "metadata": {},
+   "source": [
+    "Bei den anderen beiden Parametern ist dies nicht der Fall:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "af339c6e-f0e7-40cd-a2cf-61aaaa4df1e4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mi.draw_profile('R')\n",
+    "plt.ylabel('$\\chi^2(R, x, I_0, C)/ndof$')\n",
+    "plt.show()\n",
+    "\n",
+    "mi.draw_profile('C')\n",
+    "plt.ylabel('$\\chi^2(C, x, I_0 R)/ndof$')\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "499ffb9a-00a7-4c06-ada4-1d3a41a7f1d4",
+   "metadata": {},
+   "source": [
+    "Das liegt daran, dass $R$ und $C$ vollständig korreliert sind. Reduziert `iminuit` $C$ um ein Faktor zwei, so wird dies dadurch kompensiert, dass das optimale Minimum verlangt, dass $R$ um einen Faktor zwei größer sein muss. Das heißt, es ist ohne weitere Information nicht möglich, $R$ und $C$ näher zu bestimmen, sondern lediglich das Produkt der beiden Größen.\n",
+    "\n",
+    "Deshalb müssen wir in unserer Fitfunktion $R$ und $C$ durch die Zeitkonstante $\\tau$ ersetzen und schreiben\n",
+    "\n",
+    "$$ I = I_0 \\exp\\{-t/\\tau\\}$$\n",
+    "\n",
+    "mit $\\tau = R \\cdot C$.\n",
+    "\n",
+    "Führen wir nun erneut den Fit durch, so erhalten wir ein fehlerfreies Ergebnis..."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "847419a7-d77b-4207-8607-44af9d615ffc",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#from iminuit import Minuit, cost\n",
+    "\n",
+    "def discharge_current2(t, I0, tau):\n",
+    "    return I0 * np.exp(-t/tau)\n",
+    "\n",
+    "ls = cost.LeastSquares(\n",
+    "    data_frame['time'],\n",
+    "    data_frame['current'], \n",
+    "    data_frame['delta_current'], \n",
+    "    discharge_current2\n",
+    ")\n",
+    "mi = Minuit(ls, I0=0.9, tau=0.3)\n",
+    "mi.migrad()\n",
+    "mi.hesse()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a46c76ec-5b00-48f4-9a46-2ea083ca5dba",
+   "metadata": {},
+   "source": [
+    "... und die Werte und Fehler lassen sich über ..."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "69f540a5-e89b-4c24-aa7e-b03eaedb28d1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mi.values['I0']"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "66733c05-692d-46e3-ae82-6f84d66ef28c",
+   "metadata": {},
+   "source": [
+    "... bzw. ..."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "66e6da5b-ff32-4982-a3aa-5b9b93262073",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mi.errors['I0']"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c670cd3f-fcfb-4cfc-a8d4-eded75e9a669",
+   "metadata": {},
+   "source": [
+    "... für jeden Parameter auslesen. Dies lässt sich nun auch nutzen, um unsere Messwerte samt Fit in einem etwas schöneren Plot mit Achsenbeschriftungen darzustellen. Hierbei können wir ausnutzen, dass `iminuit` die Parameter in der Reihenfolge der Argumente unser definierten Fitfunktion speichert."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "45fcf856-c58e-424d-8fd7-15037cb6698e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.errorbar(data_frame['time'], \n",
+    "             data_frame['current'], \n",
+    "             xerr=data_frame['delta_time'], \n",
+    "             yerr=data_frame['delta_current'], \n",
+    "             ls='', \n",
+    "             marker='.',\n",
+    "             label='Data'\n",
+    "            )\n",
+    "x = np.arange(0, 1, 0.01)\n",
+    "plt.plot(x, \n",
+    "         discharge_current2(x, *mi.values), # Sternchen operator zum entpacken der Werte\n",
+    "         color='k',\n",
+    "         label='Best fit'\n",
+    "        )\n",
+    "fit_info = [\n",
+    "    f\"$\\\\chi^2$/$n_\\\\mathrm{{dof}}$ = {mi.fval:.1f} / {mi.ndof:.0f} = {mi.fmin.reduced_chi2:.1f}\",\n",
+    "]\n",
+    "for p, v, e in zip(mi.parameters, mi.values, mi.errors):\n",
+    "    fit_info.append(f\"{p} = ${v:.3f} \\pm {e:.3f}$\")\n",
+    "\n",
+    "plt.legend(title=\"\\n\".join(fit_info))\n",
+    "plt.ylabel('Current [mA]')\n",
+    "plt.xlabel('Time [s]')\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1cd73609-8593-4725-a7c4-317d6a48a72f",
+   "metadata": {},
+   "source": [
+    "# Mathematisch motivierte Herleitung des $\\chi^2$-Fits:\n",
+    "\n",
+    "Nach diesen anfänglichen Beispielen wollen wir uns eine semi-mathematische Herleitung des $\\chi^2$-Fits angucken um etwas besser zu verstehen, warum diese Methode für uns in der Physik so wichtig ist. In unserem Grundpraktikum haben wir bereits gelernt, dass Messwerte durch Zufallszahlen $x_i$ repräsentiert werden und einer gewissen **Wahrscheinlichkeitsdichtefunktion (probability density function)** $f(x)$ unterliegen.\n",
+    "\n",
+    "<figure class=\"image\">\n",
+    "<img src=\"images/MaterialPythonkurs092018/PorbDensFun.png\"  alt=\"{{ Beispiel PDF }}\" width=70%>\n",
+    "</figure>\n",
+    "\n",
+    "\n",
+    "Eine **pdf** gibt an, mit welcher **Wahrscheinlichkeit ein Wert $x_i$** innerhalb eines **infinitesimalen Intervals $\\text{d}x_i$** zu finden ist. Des Weiteren gilt, dass die Gesamtwahrscheinlichkeit gegeben ist durch $\\int_S f(x) dx = 1$. \n",
+    "\n",
+    "Nun betrachten wir folgendes Beispiel: In unserem Labor messen wir genau drei mal die Raumtemperartur T. Auch hier gilt, dass unsere Messungen der einzelnen $T_i$ einer gewissen **Wahrscheinlichkeitsdichtefunktion** folgen. Betrachten Sie nun das folgende Bild; Welche **Wahrscheinlichkeitsdichtefunktion** passt besser zu den gezeigten Daten und **Warum?**\n",
+    "\n",
+    "<figure class=\"image\">\n",
+    "<img src=\"images/MaterialPythonkurs092018/ProbMaxTemp.png\"  alt=\"{{ Beispiel PDF }}\" width=100%>\n",
+    "</figure>\n",
+    "\n",
+    "Die rechte Verteilung spiegelt unsere Messdaten besser wider. Dies können wir auch mathematisch ausdrücken. Für $N$ voreinander unabhängige Zufallszahlen bzw. Messpunkte (in unserem Beispiel $N = 3$) ist die Gesamtwahrscheinlichkeit gegeben durch das Produkt der einzelnen Wahrscheinlichkeitsdichten $f(x_i, \\theta)$ multipliziert mit dem jeweiligen infinitesimalen Element $dx_i$\n",
+    "\n",
+    "$$\\prod_{i = 1}^{N} f(x_i,\\theta) \\  dx_i \\text{   für alle } x_i \\text{ in } [x_i, x_i + dx_i]$$\n",
+    "\n",
+    "wobei $x_i$ in unserem Beispiel den Messpunkten $T_i$ und $f(x_i,\\theta)$ unserer Gaussverteilung mit $\\theta = (\\mu, \\sigma)$ entspricht. Sofern unsere Werte gut von der jeweiligen **Wahrscheinlichkeitsdichtefunktion** repräsentiert werden, d.h. wir die richtigen Parameter $\\theta$ gewählt haben (wie im rechten oberen Plot), gilt \n",
+    "\n",
+    "$$ \\prod_{i = 1}^{N} f(x_i,\\theta)  dx_i \\ \\ \\text{ist} \\ \\textbf{maximal.}$$\n",
+    "\n",
+    "Da die einzelnen $dx_i$ von unseren Parametern $\\theta$ unabhängig sind, gilt die gleiche Argumentation auch für \n",
+    "\n",
+    "$$ \\mathcal{L}(x_1 ... x_N; \\theta_1 ... \\theta_N) = \\prod_{i = 1}^{N} f(x_i,\\theta)$$ \n",
+    "\n",
+    "wobei $\\mathcal{L}(x_1 ... x_N; \\theta_1 ... \\theta_N)$ die sogenannte **\"likelihood\"** function darstellt.\n",
+    "\n",
+    "Wie kommen wir nun von der **likelihood function** auf unsere **Methode der kleinsten Quadrate** und das Fitten einer Funktion $\\lambda(x; \\ $**$\\phi$**$)$ an die gemessenen Punkte **$(x,y)$**? Dazu brauchen wir noch einen Zwischenschritt. Oftmals ist es einfacher, statt die **likelihood function** zu maximieren, die so genannte **log likelihood function**\n",
+    "\n",
+    "$$ \\log( \\mathcal{L}(x_1 ... x_N; \\theta_1 ... \\theta_N)) = \\sum_{i = 1}^{N} \\log(f(x_i,\\theta))$$\n",
+    "\n",
+    "zu maximieren. Dies ist im Grunde das Gleiche, da der Logarithmus eine monoton-steigende Funktion ist. Auch in unserem Fall der **Methode der kleinsten Quadrate** benötigen wir die **log likelihood function**. \n",
+    "\n",
+    "Stellen Sie sich nun vor, wir haben eine Messung mit $N$ voneinander unabhängigen Messpunkten (x,y). Des Weiteren nehmen wir an, dass alle $x_i$ ohne Fehler sind und dass unsere $y_i$ gaußförmig um einen unbekannten wahren Wert $\\lambda_i$ (sprich $\\lambda_i$ entspricht dem Erwartungswert $\\mu_i$ unserer Gaußverteilung) mit einer bekannten Varianz $\\Delta y_i^2$ verteilt sind (Diese Annahme lässt sich mit dem zentralen Grenzwertsatz begründen, so lange der Fehler sich aus der Summe kleiner Fehler zusammensetzt). Die dazugehörige  **likelihood function** ist dann gegeben durch:\n",
+    "\n",
+    "$$ \\mathcal{L}(y_1 ... y_N; \\lambda_1 ... \\lambda_N, \\Delta y_1 ... \\Delta y_N)) = \\prod_{i = 1}^{N}\\frac{1}{\\sqrt{2 \\pi \\Delta y_i^2}} \\cdot \\exp \\bigg( \\frac{ -(y_i - \\lambda_i)^2}{2 \\cdot \\Delta y_i^2}\\bigg)$$\n",
+    "\n",
+    "Beziehungsweise die **log likelihood function** mit $\\lambda_i = \\lambda(x_i; \\phi)$ ergibt sich zu\n",
+    "\n",
+    "$$ \\log(\\mathcal{L}(y, \\theta)) \\approx -\\frac{1}{2} \\sum_{i = 1}^{N}\\bigg( \\frac{ (y_i - \\lambda(x_i; \\phi))^2}{\\Delta y_i^2}\\bigg)$$\n",
+    "\n",
+    "wobei die konstanten Terme, welche nicht von unserer Funktion $\\lambda(x_i; \\phi)$ abhängen, vernachlässigt worden sind. Durch den Faktor $-\\frac{1}{2}$ ist das Maximieren dieser **log likelihood function** gleich dem Minimieren von\n",
+    "\n",
+    "$$ \\chi(\\phi_1 ... \\phi_N)^2 = \\sum_{i = 1}^{N} \\frac{ (y_i - \\lambda(x_i; \\phi))^2}{\\Delta y_i^2}$$\n",
+    "\n",
+    "Diese Funktion ist unsere gesuchte **Methode der kleinsten Quadrate**. Mit ihrer Hilfe kann eine beliebige Funktion $\\lambda(x; \\phi)$, welche linear in ihren Parametern $\\phi$ ist, an unsere Messdaten $(x,y\\pm\\Delta y)$ gefittet werden. Dabei stellt der Fitprozess selbst lediglich ein Minimierungsproblem dar. Im Folgenden sind unsere Annahmen noch einmal grafisch in einem Beispiel dargestellt.\n",
+    "\n",
+    "<figure class=\"image\">\n",
+    "<img src=\"images/MaterialPythonkurs092018/LeastSquare.png\"  alt=\"{{ Least Square Beispiel }}\" width=100%>\n",
+    "</figure>"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "be4a8d21-29db-4866-9117-8746b80d5945",
+   "metadata": {},
+   "source": [
+    "Wie ein Algorithmus bei der Minimierung vorgeht, sprengt den Rahmen dieses Vorversuchs. Hier sei auf entsprechende Vorlesungen verwiesen. Aber um einen kleinen Einblick zu erhalten, kann man sich die Werte der Parameter und von $\\chi^2$ für jeden Schritt ausgeben lassen. Dazu wird der Parameter `verbose` auf 1 gesetzt.\n",
+    "\n",
+    "Man erkannt, dass für jeden Parameter zunächst separat geprüft wird, welche Änderung (größer oder kleiner) die Kostenfunktion minimiert. Danach beginnt die eigentliche Minimierung der Kostenfunktion durch den Algorithmus."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "43bfd15e-7b68-4b70-bc06-0b23f89f7bff",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ls = cost.LeastSquares(data_frame['time'], data_frame['current'], data_frame['delta_current'], discharge_current, verbose=1)\n",
+    "\n",
+    "mi = Minuit(ls, I0=0.9, R=10*10**3, C=10**-6)\n",
+    "mi.migrad()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c664e4f0-0226-4be0-91f4-7c28aee16a4a",
+   "metadata": {},
+   "source": [
+    "__Zur Erinnerung__: die Methode der kleinsten Quadrate berücksichtigt nur Fehler in `y` aber nicht in `x`. Hierfür sind komplexere Methoden notwendig, die wir hier nicht betrachten wollen.\n",
+    "Ebenfalls wichtig ist, die Statusmeldungen von `iminuit` zu prüfen, d.h. eine <span style=\"color:yellow\"> **gelbe** </span> Box zeigt an, das man sich Gedanken über das Ergebnis machen sollte (in unserem Fall, dass die Variablen $R$ und $C$ korreliert sind) und eine ,<span style=\"color:#7878FA\"> **violette** </span> Box, dass der Fit nicht konvergiert ist und das Ergebnis nicht verwendet werden kann."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "jupyter",
+   "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.10.11"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

Herleitung_Methode_der_kleinsten_Quadrate.ipynb

@@ -1,100 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "id": "f3b16f1a-521e-4caf-828e-85251faf2c6c",
-   "metadata": {},
-   "source": [
-    "Dieses Dokument beinhaltet die Herleitung der Methode der kleinsten Quadrate, die im Kapitel 1 beim Fitten benutzt wird."
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "a144903b-e362-40de-8af0-7da7b39fe260",
-   "metadata": {},
-   "source": [
-    "### Methode der kleinsten Quadrate\n",
-    "\n",
-    "Im folgenden wolllen wir die **Methode der kleinsten Quadrate (Least Squares)** näher beleuchten. Diese Methode wird oft benutzt, um eine Funktion $\\lambda(x; \\ $**$\\phi$**$)$ mit den Funktionsparametern $\\mathbf{\\phi}$ an die gemessenen Punkte **$(x,y)$** anzupassen. Um jedoch die **Methode der kleinsten Quadrate** zu verstehen, wollen wir sie erst einmal anschaulich und mathematisch herleiten. Dabei stüzen wir uns im Folgenden auf eine Herleitung aus dem Buch **\"Statistical Data Analysis\"**  von **Glen Cowan**.\n",
-    "\n",
-    "In unserem Grundpraktikum haben wir bereits gelernt, dass Messwerte durch Zufallszahlen $x_i$ representiert werden und einer gewissen **Wahrscheinlichkeitsdichtefunktion (probability density function)** $f(x)$ unterliegen.\n",
-    "\n",
-    "<figure class=\"image\">\n",
-    "<img src=\"images/MaterialPythonkurs092018/PorbDensFun.png\"  alt=\"{{ Beispiel PDF }}\" width=70%>\n",
-    "</figure>\n",
-    "\n",
-    "\n",
-    "Eine **pdf** gibt an, mit welcher **Wahrscheinlichkeit ein Wert $x_i$** innerhalb eines **infinitesimalen Intervals $\\text{d}x_i$** zu finden ist. Des Weitren gilt, dass die Gesamtwahrscheinlichkeit gegeben ist durch $\\int_S f(x) dx = 1$. \n",
-    "\n",
-    "Nun betrachten wir folgendes Beispiel: In unserem Labor messen wir genau drei mal die Raumtemperartur T. Auch hier gilt, dass unsere Messung der einzelnen $T_i$ einer gewissen **Wahrscheinlichkeitsdichtefunktion** folgen. Betrachten Sie nun das folgende Bild; Welche **Wahrscheinlichkeitsdichtefunktion** passt besser zu den gezeigten Daten und **Warum?**\n",
-    "\n",
-    "<figure class=\"image\">\n",
-    "<img src=\"images/MaterialPythonkurs092018/ProbMaxTemp.png\"  alt=\"{{ Beispiel PDF }}\" width=100%>\n",
-    "</figure>\n",
-    "\n",
-    "Die rechte Verteilung spiegelt unsere Messdaten besser wieder. Dies können wir auch mathematisch ausdrücken. Für $N$ voreinander unabhängige Zufallszahlen bzw. Messpunkte (in unserem Beispiel $N = 3$) ist die Gesamtwahrscheinlichkeit gegeben durch das Produkt der einzelnen Wahrscheinlichkeitsdichten $f(x_i, \\theta)$ multipliziert mit dem jeweiligen infinitesimalen element $dx_i$\n",
-    "\n",
-    "$$\\prod_{i = 1}^{N} f(x_i,\\theta) \\  dx_i \\text{   für alle } x_i \\text{ in } [x_i, x_i + dx_i]$$\n",
-    "\n",
-    "wobei $x_i$ in unserem Beispiel den Messpunkten $T_i$ und $f(x_i,\\theta)$ unserer Gausverteilung mit $\\theta = (\\mu, \\sigma)$ entspricht. Sprich sofern unsere Werte gut von der jeweiligen **Wahrscheinlichkeitsdichtefunktion** repräsentiert werden, d.h. wir die richtigen Parameter $\\theta$ gewählt haben (wie im rechten oberen Plot), gilt \n",
-    "\n",
-    "$$ \\prod_{i = 1}^{N} f(x_i,\\theta)  dx_i$$ \n",
-    "\n",
-    "ist **maximal**. Da die einzelnen $dx_i$ von unseren Parametern $\\theta$ unabhängig sind, gilt die gleiche Argumentation auch für \n",
-    "\n",
-    "$$ \\mathcal{L}(x_1 ... x_N; \\theta_1 ... \\theta_N) = \\prod_{i = 1}^{N} f(x_i,\\theta)$$ \n",
-    "\n",
-    "wobei $\\mathcal{L}(x_1 ... x_N; \\theta_1 ... \\theta_N)$ die sogenannte **likely hood function** darstellt.\n",
-    "\n",
-    "Wie kommen wir nun von der **likely hood function** auf unsere **Methode der kleinsten Quadrate** und das Fitten einer Funktion $\\lambda(x; \\ $**$\\phi$**$)$ an die gemessenen Punkte **$(x,y)$**? Dazu brauche wir noch einen Zwischenschritt. Oftmals ist es einfacher, statt die **likely hood function** zu maximieren, die so genannte **log likely hood function**\n",
-    "\n",
-    "$$ \\log( \\mathcal{L}(x_1 ... x_N; \\theta_1 ... \\theta_N)) = \\sum_{i = 1}^{N} \\log(f(x_i,\\theta))$$\n",
-    "\n",
-    "zu maximieren. Dies ist im Grunde das Gleiche, da der Logarithmus eine monoton-steigende Funktion ist. Auch in unserem Fall der **Methode der kleinsten Quadrate** benötigen wir die **log likely hood function**. \n",
-    "\n",
-    "Stellen Sie sich nun vor, wir haben eine Messung mit $N$ voneinander unabhängigen Messpunkten (x,y). Des Weiteren nehmen wir an, dass alle $x_i$ ohne Fehler sind und dass unsere $y_i$ gaußförmig um einen unbekannten Wahrenwert $\\lambda_i$ (sprich $\\lambda_i$ entspricht dem Erwartungswert $\\mu_i$ unserer Gaußverteilung) mit einer bekannten Varianz $\\Delta y_i^2$ verteilt sind (Diese Annahme lässt sich mit dem zentralen Grenzwertsatz begründen, so lange der Fehler sich aus der Summe kleiner Fehler zusammensetzt). Die dazugehörige  **likely hood function** ist dann gegeben durch:\n",
-    "\n",
-    "$$ \\mathcal{L}(y_1 ... y_N; \\lambda_1 ... \\lambda_N, \\Delta y_1 ... \\Delta y_N)) = \\prod_{i = 1}^{N}\\frac{1}{\\sqrt{2 \\pi \\Delta y_i^2}} \\cdot \\exp \\bigg( \\frac{ -(y_i - \\lambda_i)^2}{2 \\cdot \\Delta y_i^2}\\bigg)$$\n",
-    "\n",
-    "Beziehungsweise die **log likely hood function** mit $\\lambda_i = \\lambda(x_i; \\phi)$ ergibt sich zu\n",
-    "\n",
-    "$$ \\log(\\mathcal{L}(y, \\theta)) \\approx -\\frac{1}{2} \\sum_{i = 1}^{N}\\bigg( \\frac{ (y_i - \\lambda(x_i; \\phi))^2}{\\Delta y_i^2}\\bigg)$$\n",
-    "\n",
-    "wobei die konstanten Terme, welche nicht von unserer Funktion $\\lambda(x_i; \\phi)$ abhängen, vernachlässigt worden sind. Durch den Faktor $-\\frac{1}{2}$ ist das Maximieren dieser **log likely hood function** gleich dem Minimieren von\n",
-    "\n",
-    "$$ \\chi(\\phi_1 ... \\phi_N)^2 = \\sum_{i = 1}^{N} \\frac{ (y_i - \\lambda(x_i; \\phi))^2}{\\Delta y_i^2}$$\n",
-    "\n",
-    "Diese Funktion ist unsere gesuchte **Methode der kleinsten Quadrate**. Mit ihrer Hilfe kann eine beliebige Funktion $\\lambda(x; \\phi)$, welche liniear in ihren Parametern $\\phi$ ist, an unsere Messdaten $(x,y\\pm\\Delta y)$ gefittet werden. Dabei stellt der Fitprozess selbst lediglich ein Minimierungsproblem dar. Im Folgenden sind unsere Annahmen noch einmal grafisch in einem Beispiel dargestellt.\n",
-    "\n",
-    "<figure class=\"image\">\n",
-    "<img src=\"images/MaterialPythonkurs092018/LeastSquare.png\"  alt=\"{{ Least Square Beispiel }}\" width=100%>\n",
-    "</figure>\n",
-    "\n",
-    "Es gibt verschiedene Arten von Algorithmen um Minimierungsprobleme zu lösen. Wie diese genau aufgebaut sind, lernen Sie in anderen Progrmmierveranstaltungen wie zum Beispiel *Programmieren für Physiker* oder *Computer in der Wissenschaft*. Zum Glück haben uns bereits in Python andere Menschen diese Arbeit abgenommen und wir können aus dem Package `scipy.optimize` die Funktion `curve_fit` verwenden.\n",
-    "\n",
-    "Hierbei stellt `curve_fit` eine Methode dar, Fit-Funktionen nach der obigen vorgestellten Methode der *kleinsten Quadraten* zu bestimmen. Dies hat zur Folge, dass lediglich die y-Fehler der Messwerte für den Fit verwendet werden können."
-   ]
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
-   "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.9.6"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 5
-}

README.md

@@ -1,17 +1,17 @@
-# PGP1 Python Einfuehrung
+# PGP2 Vorversuch *Erweiterte Statistik*
 
 Im folgenden finden Sie die Installationsanweisungen für Python und das Kursmaterial. Sofern Sie Fragen haben, bzw. Hilfe benötigen, können Sie an diese Personen eine Email schreiben:
 
-Mo (Mohamad Bitar): mobitar@students.uni-mainz.de
+Matthias Hoek: matthias.hoek@uni-mainz.de
 
 ## Kursmaterial installieren:
 
-Die installation des Kursmaterials geschieht vollständig automatisch. Einfach den Anweisungen im PGP-Skript folgen und diesen Link verwenden: [Kursmaterial Installieren](https://jupyterhub.zdv.uni-mainz.de/hub/user-redirect/git-pull?repo=https%3A%2F%2Fgitlab.rlp.net%2Fpgp%2Fpgp1-python-einfuehrung&urlpath=tree%2Fpgp1-python-einfuehrung%2F&branch=master)
+Die installation des Kursmaterials geschieht vollständig automatisch. Einfach den Anweisungen im PGP-Skript folgen und diesen Link verwenden: [Kursmaterial Installieren](https://jupyterhub.zdv.uni-mainz.de/hub/user-redirect/git-pull?repo=https%3A%2F%2Fgitlab.rlp.net%2Fpgp%2Fpgp1-python-einfuehrung&urlpath=tree%2Fpgp1-python-einfuehrung%2F&branch=add_pgp2_lecture)
 
 ## Einloggen auf dem Jupyterhub:
 
 Sie können auch den durch die ZDV angebotenen JupyterHub (https://jupyterhub.zdv.uni-mainz.de) zur Bearbeitung Ihrer Notebooks verwenden. Hier müssen Sie sich zunächst mit Ihrem Uni-Account anmelden. Danach erscheint eine Auswahlseite, auf der Sie die Art der Jupyter Umgebung auswählen. Für das Praktikum ist die Standardumgebung die richtige Wahl, s. Bild unten.
 
-<img src = "./images/Screenshot_ZDV_JupyterHub.png" >
+<img src = "./images/Screenshot_ZDV_JupyterHub_neu.png" >
 
 Klicken Sie auf die Schaltfläche **Spawn**, dann öffnet sich das Notebook Dashboard.

custom_set_up.py

@@ -3,7 +3,7 @@
 # 
 # Since the directory structure for the custom.css file changed 
 # this set-up only works for IPython and Jupyter versions after 
-# the so called the "big splitt" (version 4 and higher).
+# the so called the "big split" (version 4 and higher).
 import os
 import jupyter_core
 import IPython