{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "# Rechnernutzung in der Physik\n", "**Institut für Experimentelle Teilchenphysik**
\n", "**Institut für Theoretische Teilchenphysik**
\n", "Prof. G. Quast, Prof. M. Steinhauser
\n", "Dr. A. Mildenberger, Dr. Th. Chwalek
\n", "[Ilias Seite zum Kurs](https://ilias.studium.kit.edu/ilias.php?ref_id=2212147&cmdClass=ilrepositorygui&cmdNode=x1&baseClass=ilrepositorygui)
\n", "WS 2023/24 – Blatt 03
\n", "Abgabe: Montag 27.11.2023 bzw. Dienstag 28.11.2023\n", "\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Auf dem dritten Übungsblatt beschäftigen Sie sich mit einem der wichtigsten Grundbegriffe der modernen Parameteranpassung - der Maximum-Likelihood-Methode. Sie haben mit großer Wahrscheinlichkeit im Laufe Ihres Studiums oder Privatlebens bereits die Maximum-Likelihood-Methode (oder ihre Abwandlungen) verwendet, um die Steigung einer Geraden oder andere Parameter einer Messung zu bestimmen. Mit diesem Übungsblatt beschäftigen Sie sich noch mal ganz genau mit der Methodik und nehmen hoffentlich mit, auf welchem Prinzip die unterschiedlichen eingebauten Fitfunktionen in *SciPy*, *kafe2*, *probfit* und vielen mehr (ja, auch Excel) aufgebaut sind. Sie werden merken, dass sich hinter all den angesprochenen Methoden keine schwarze Magie befindet, sondern die Umsetzung und Überführung auf eigene Programme ganz einfach ist. Mit einem guten Verständis der Maximum-Likelihood-Methode können Sie dann nicht nur die Methode der kleinsten Quadrate verstehen, sondern auch schon bald Anpassungen wie diese selbst durchführen." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "# Aufgabe 1: Poisson-Likelihood \n", "---\n", "\n", "In der ersten Aufgabe beschäftigen Sie sich mit einem vermutlich bereits bekanntem Beispiel aus dem Praktikum: Sie haben eine radioaktive Quelle erhalten und messen elektronisch die Klicks eines Geiger-Müller-Zählrohrs in aufeinanderfolgenden, gleich langen Zeitintervallen. Die Messdaten haben Sie in der Datei *Geiger_Zaehler.txt* abgespeichert und wollen nun für die weitere Auswertung die Anzahl der radioaktiven Zerfälle pro Zeitintervall bestimmen. \n", "\n", "Bei der Auswertung der Daten erinnern Sie sich an die Vorlesung *Rechnernutzung* und die Parameterschätzung mithilfe der Maximum Likelihood Methode und beschließen, statt *kafe2* zu nutzen, die gesuchte Zahl selbst abzuschätzen. Sie schauen zurück in die Folien der Vorlesung und entdecken die Funktion\n", "\n", "$$ \n", "L\\left(\\left\\{x_i\\right\\}|a\\right)=\\prod_{i\\leq n}f\\left(x_i,a\\right).\n", "$$\n", "\n", "Sie erinnern sich daran, dass $\\left\\{x_i\\right\\}$ ihre Messungen sind, $f(x,a)$ die zugrunde liegende Wahrscheinlichkeitsdichte und $a$ der gesuchte Parameter ist." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## a) Likelihoodfunktion\n", "\n", "Im ersten Aufgabenteil überprüfen Sie graphisch die Form der Likelihoodfunktion.\n", "> Frage: Warum ist es gerechtfertig die Poissonverteilung als zugrunde liegende Wahrscheinlichkeitsdichte für den Ausgang des Zählexperiments zu wählen? Welche Abweichungen von der Annahme gibt es beim Aufbau und der Messung?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "> Antwort:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Implementieren Sie die Likelihoodfunktion $L\\left(\\left\\{x_i\\right\\}|a\\right)$. Es bietet sich an, die Wahrscheinlichkeitsdichte als Argument der Funktion zu übergeben. Denken Sie daran, dass Sie einen gewissen Parameterereich abdecken müssen. Es bietet sich also an nicht nur $x$, sondern auch $a$ als Array zu behandeln." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import scipy.special as sp" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Poisson distribution\n", "def fPoisson(x, nu_P):\n", " k=np.around(x)\n", " return (nu_P**k) / np.exp(nu_P) / sp.gamma(k+1.)\n", "\n", "# Likelihood function\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Betrachten Sie nun die Likelihoodfunktion für eine steigende Zahl an Messergebnissen $L\\left(\\left\\{x_i\\right\\}|a\\right)$. Stellen Sie dafür die Likelihoodfunktion für $a$ im Intervall $[2,9]$ dar für\n", "* nur den ersten Messwert\n", "* den ersten und den zweiten Messwert\n", "* die ersten 10 Messwerte\n", "* die ersten 100 Messwerte.\n", "\n", "Sie können mithilfe von *np.loadtxt(\"filename\")* die Messdaten aus dem Verzeichnis laden, indem sie das Notebook ausführen. Sollten Sie also das Notebook von ihrem eigenen PC in [jupytermachine.etp.kit.edu](jupytermachine.etp.kit.edu) laden oder lokal auf ihrem Rechner arbeiten, denken Sie daran, zuerst die Datei *Geiger_Zaehler.txt* herunterzuladen." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# load observations\n", "\n", "\n", "# creater parameter intervall\n", "\n", "\n", "# list for desired observation numbers\n", "\n", "\n", "#build plots\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "> Frage: Erläutern Sie das Verhalten des Maximums der Likelihoodfunktion. Schätzen Sie den Wert des gesuchten Paramters ab." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "> Antwort: " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## b) Negative Log Likelihood (NLL)\n", "\n", "Beachten Sie die $y$-Skala bzw. den Wertebereich der Likelihoodfunktion für die vier verschiedenen Anzahlen an Messpunkten. Für eine große Anzahl an Messpunkten nimmt die Funktion immer kleinere Werte an. Eine Logarithmierung des Ausdruck vereinfacht nicht nur die mathematische Analyse, sondern hilft auch numerisch, da das Produkt einer großen Anzahl kleiner Wahrscheinlichkeiten die numerische Genauigkeit des Computers für Fließkommaoperationen leicht unterschreiten kann, und dies wird gelöst, indem stattdessen die effizientere Summe der natürlichen Logarithmen der Wahrscheinlichkeiten (\"Log Likelihood\") berechnet wird. Statt nach dem Maximum der Likelihoodfunktion zu suchen, wird häufig nach dem Minimum der negativen Log Likelihood (NLL) gesucht. \n", "\n", "$$ \n", "-\\ln\\left(L\\left(\\left\\{x_i\\right\\}|a\\right)\\right)=-\\sum_{i\\leq n}\\ln\\left(f\\left(x_i,a\\right)\\right).\n", "$$\n", "\n", "Wiederholen Sie den ersten Aufgabenteil mit der NLL Methode. Implementieren Sie dazu eine neue Funktion, die die NLL ausgibt. Betrachten Sie den gleichen Parameterbereich aber diesmal für \n", "* nur den ersten Messwert\n", "* den ersten und den zweiten Messwert\n", "* die ersten 10 Messwerte\n", "* **alle Messwerte** statt nur 100 Messwerte\n", "\n", "Stellen Sie Ihr Ergebnis geeignet dar." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Negative Log Likelihood function\n", "\n", "\n", "#build plots\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## c) Parameterschätzung\n", "\n", "Bestimmen Sie nun den gesuchten Parameter und schätzen Sie dessen $1\\sigma$ Konfidenzintervall ab, indem Sie die *asymmetrischen* Grenzen für die Unsicherheit angeben.\n", "\n", "Für den ersten Schritt kann die *NumPy*-Methode *np.argmin()* hilfreich sein. Sie gibt den Index des Arrays am Minimum aus.
\n", "Den zweiten Schritt erhalten Sie, indem Sie diejenigen Punkte finden, deren NLL Wert um $1/2$ größer sind als der des Maximums. Das $1\\sigma$ Konfidenzintervall entspricht dann dem Wertebereich $[\\hat{a}-\\Delta^-,\\hat{a}+\\Delta^+]$. Die üblich Angabe für den besten Schätzwert erfolgt dann über folgende Darstellung:\n", "$$\\hat{a}^{+\\Delta^+}_{-\\Delta^-}.$$\n", "Hier eine schnelle Illustration:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = np.linspace(0, 10, 1000)\n", "y = (x - 3)**2 + 0.2\n", "\n", "dots_x = [3 - 1/np.sqrt(2), 3, 3 + 1/np.sqrt(2)]\n", "dots_y = [0.7, 0.2, 0.7]\n", "\n", "plt.plot(x, y, 'b-', label='NLL')\n", "plt.hlines(0.7, 1, 5, color='green', label=r'NLL$[\\hat{a}]$+0.5')\n", "plt.plot(dots_x, dots_y, 'ro')\n", "\n", "colors = ['gold', 'orange', 'tomato']\n", "labels = [r'Unterschranke $\\hat{a}-\\Delta^-$', r'Minimum $\\hat{a}$', r'Oberschranke $\\hat{a}+\\Delta^+$']\n", "for i in range(3):\n", " plt.vlines(dots_x[i], 0, dots_y[i], linestyle='dashed', color=colors[i], label=labels[i])\n", "\n", "plt.ylim(-0.1, 1)\n", "plt.xlim(1, 5)\n", "plt.xlabel(r'$a$')\n", "plt.ylabel(r'NLL $-\\ln L\\left(\\left\\{x_i\\right\\}|a\\right)$')\n", "plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Es bietet sich an die Zwischenschritte als Funktion zu implementieren, damit sie erneut verwendet werden können." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# define funtcion for estimator and CL calculation\n", "\n", "\n", "# print results\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "# Aufgabe 2: Profile Likelihood \n", "---\n", "In der zweiten Aufgabe schauen Sie sich ein vielleicht bereits bekanntes Beispiel aus dem Praktikum an: Der Myon-Zerfall. Das Myon hat dieselben Eigenschaften wie das Elektron, nur ist es ca. 200 Mal schwerer. Myonen entstehen beispielsweise als Zerfallsprodukte kosmischer Strahlung und können dann am Boden z. B. durch ein praktisch zeitgleiches Signal (\"Koinzidenz\") in mehreren Lagen von Szintillationszählern (Detektoren, die bei Teilchendurchgang Licht emittieren, das mit Photondetektoren nachgewiesen werden kann) nachgewiesen werden. Dabei zerfallen sie größtenteils gemäß dem folgenden Prozess, illustriert als \"Feynman-Diagramm\", bei dem das Myon in ein Elektron, ein Elektron-Antineutrino und ein Myon-Neutrino zerfällt und ein virtuelles W-Boson als \"Kraftteilchen\" fungiert:\n", "\n", "![MuonDecay](https://atlas.physicsmasterclasses.org/wpath_files/img/Feynman/Paarerzeugung.png)\n", "\n", "Bei physikalischen Experimenten treten häufig Störparameter auf, die für das physikalische Ergebnis irrelevant sind, aber dessen Unsicherheit beeinflussen. Die Berechnung und Analyse der Profile-Likelihood für den jeweils interessierenden physikalischen Parameter ist eine einfache Möglichkeit, die durch Störparameter verursachten Unsicherheiten zu berücksichtigen.\n", "\n", "In der Datei *tau_mu.dat* finden Sie 150 Werte gemessener Lebensdauern von Myonen, die im oder in der Nähe eines Myondetektors gestoppt wurden. Die Zeitdifferenz zwischen dem Nachweis des einfallenden Myons und dem Eintreffen des Elektrons aus dem Zerfall ist die Messgröße. Bedingt durch das Messverfahren sind nur Zeiten zwischen $\\Delta t_{min}=1.0\\,\\mu$s und $\\Delta t_{max}=11.5\\,\\mu$s erfasst. Als Störparameter bei den Messungen tritt ein Untergrund aus Zufallskoinzidenzen von Signalen auf, die im Bereich von $[\\Delta t_{min},\\Delta t_{max}]$ als flach verteilt angenommen werden können.\n", "\n", "In den folgenden Teilaufgaben soll die Myonenlebensdauer $\\tau_{\\mu}$ bestimmt werden. Sie werden dabei bemerken, wie einfach es ist, mit dem Wissen aus der Vorlesung und der [ersten Aufgabe](#Aufgabe1) eine reelle Größe abzuschätzen." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## a) Verteilungsdichte der Zeitdifferenzen\n", "\n", "Zunächst soll die zugrunde liegende Wahrscheinlichkeitsdichte der Messdaten erstellt werden. Aus dem Aufbau der Messung ist zu entnehmen, dass die Verteilung der gemessenen Zeitdifferenzen aus der Summe einer Gleichverteilung und einer Exponentialverteilung besteht. \n", "* Die **Gleichverteilung** für den Untergrund aus Zufallskoinzidenzen ist durch den relativen Anteil $f_b$ an der Gesamtanzahl gemessener Zeitdifferenz als Störparameter skaliert.\n", "* Die **Exponentialverteilung** für den Myonenzerfall ist charakterisiert durch die Myonenlebensdauer $\\tau_{\\mu}$ und beinhaltet die noch zu bestimmende Normierungskonstante $C_n$. Außerdem wird der Störparameter $f_b\\in[0,1]$ berücksichtigt durch den Faktor $(1-f_b)$, so dass die relativen Anteile in der Summe Eins ergeben.\n", "\n", "\n", "$$\n", "\\mathrm{PDF}(t)=C_n(1-f_b)\\exp\\left\\{-\\frac{t}{\\tau_{\\mu}}\\right\\}+\\frac{f_b}{\\Delta t_{max}-\\Delta t_{min}}\n", "$$\n", "\n", "Zur Vorbereitung der Auswertung führen Sie folgende Schritte aus:\n", "1. Implementieren Sie eine Funktion, die die Wahrscheinlichkeitsdichte ausgibt. Beachten Sie, dass die PDF auf Eins normiert sein muss. Bestimmen Sie dazu die Normierungskonstante $C_n$.\n", "2. Implementieren Sie eine Funktion, die die NLL der PDF ausgibt, falls diejenige, die Sie in der [ersten Aufgabe](#Aufgabe1) geschrieben haben, nicht ausreichen sollte.\n", "3. Lesen Sie die Daten aus der gegebenen Datei ein und histogrammieren Sie diese, damit Sie eine grobe Vorstellung des Problems erhalten." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# pdf for myuon decay\n", "\n", "\n", "# NLL if not implemented\n", "\n", "\n", "# load observations\n", "\n", "\n", "# plot obsevration\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## b) Bestimmung der Profile Likelihood\n", "\n", "Zur Bestimmung der Profile Likelihood wählen Sie 1000 Werte $\\tau_i$ für die Lebensdauer $\\tau_{\\mu}$ im Intervall zwischen $1.5\\,\\mu$s und $2.9\\,\\mu$s. Minimieren Sie für jeden Wert die NLL $-\\ln L(\\Delta t;\\tau_i,f_b)$ bezüglich des Parameters $f_b$. Die Untergrundrate $f_b$ können Sie durch 200 Werte im Intervall $[0,0.3]$ wählen – dies entspricht einem kleinen Untergrundanteil von ca. 0-30%.
\n", "Stellen Sie die Profile Likelihood geeignet dar. Zeichnen Sie ebenfalls die Linie ein, die die Abschätzung des $1\\sigma$ Konfidenzintervalls zeigt." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# create array for tau and f_b\n", "\n", "\n", "# creat profile likelihood\n", "\n", " \n", "# plot profile likelihood\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## c) Parameterschätzung\n", "\n", "Bestimmen Sie als letzten Schritt \n", "* den bestmöglichen Schätzer für die Myonenlebenszeit $\\hat{\\tau}_{\\mu}$, \n", "* dessen asymmetrischen Grenzen ür die Unsicherheiten $+\\Delta^+,-\\Delta^-$\n", "* und den bestmöglichen Schätzer für die Untergrundrate $\\hat{f}_b$ unter der Annahme $\\hat{\\tau}_{\\mu}$.\n", "\n", "Stellen Sie die histogrammierten Zeitdifferenzen, die Verteilung des Signals unter der Annahme der bestmöglichen Schätzer und alle Ergebnisse in einem Plot dar." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# calculate best estimate for tau and f_b\n", "\n", "\n", "# plot final result\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [] } ], "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": 4 }