{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "---\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 05
\n", "Abgabe: Montag 11.12.2023 bzw. Dienstag 12.12.2023\n", "\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Auf dem letzten Übungsblatt im Statistik-Block sollen viele der Konzepte und Techniken aus dem Block der Veranstaltung \"Rechnernutzung in der Physik\" in einer Projektaufgabe angewendet und vertieft werden. In der Projektaufgabe sollen mit Hilfe eines Ensemble-Tests aus vielen Pseudoexperimenten (auch:Toy-MC) die Eigenschaften einer (einfachen) Parameteranpassung an einem nahezu realistischen Beispiel untersucht werden. Zum Abschluss des Statistik-Blocks erwarten wir, dass Sie aus den Beispielen in Vorlesung und Tutorium, aus den weiteren hier angegebenen Tutorials und aus der Online-Dokumentation der verwendeten Software-Werkzeuge eigenständig Elemente zur Lösung der Aufgaben zusammenstellen.\n", "\n", "Das Blatt darf innerhalb einer Tutoriumsgruppe zu maximal dritt bearbeitet werden. Die Abfrage erfolgt jedoch gleichmäßig. Jede teilnehmende Person muss in der Lage sein jeden Programmschritt vollständig erklären zu können.\n", "**Anmerkung:** die Bearbeitung dieses Blattes ist **verpflichtend** für den Erhalt des Übungsscheins." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "# Projektaufgabe \"Ensemble-Test\"\n", "---\n", "\n", "Gemessen werden soll die Lebensdauern von Teilchen, die nach Durchgang durch einen Detektor in einem Absorber gestoppt werden. Die Zerfälle werden über die Zerfallsprodukte wieder vom selben Detektor registriert. Die theoretische Lebensdauer der Teilchen beträgt $\\tau=$2µs, dies entspricht in etwa der Lebensdauer des Myons. Wegen der Überlappung der Detektorsignale können Lebensdauern kleiner als $t_{\\rm min}=$1µs nicht zuverlässig gemessen werden. Die Messelektronik ist nur bis zum Zeitpunkt $t_{\\rm max}=$10µs aktiv. Die Zahl der so registrierten Zerfälle ist mit $N=50$ registrierten Ereignissen nur sehr klein, so dass in einem ungebinner Maximum-Likelihood-Fit (d.h. alle Datenpunkte werden in der Likelihood-Funktion berücksichtigt, nicht nur die Einträge in Bins eines Histogramms) eine Exponentialfunktion an die im Intervall $[t_{\\rm min}, t_{\\rm max}]$ gemessenen Lebensdauern angepasst werden soll.\n", "\n", "Es stellen sich zwei Fragen:\n", "1. Ist der Schätzwert für die Lebensdauer $\\hat{\\tau}$ *erwartungstreu* (auch: unverzerrt, engl.: unbiased)? \n", "2. Wie gut ist die *Abdeckung* (engl.: coverage) des Konfidenzintervalls für $\\hat{\\tau}$?\n", "\n", "**Zur Erläuterung:** Die Unsicherheit auf die gemessene Lebensdauer soll aus einem Scan der negativen Log-Likelihood (NLL) gewonnen werden. Das Intervall $[\\hat{\\tau}-\\Delta^-, \\hat{\\tau}+\\Delta^+]$ mit den asymmetrischen Unsicherheiten $\\Delta^+$ und $\\Delta^-$ heißt Konfidenzintervall. Zur Erinnerung: in der frequentistischen Statistik liegt der wahre Wert eines Parameters in einem Bruchteil $\\alpha$ aller aus Daten konstruierten Konfidenzintervalle. Eine häufige Wahl ist $\\alpha\\approx 0,683$, dies entspricht einer Standardabweichung der Gaußverteilung. Ein wichtiger Test bei der Bestimmung von Konfidenzintervallen mithilfe von Ensemble-Tests ist die Überprüfung der Abdeckung, also die Frage, ob der wahre Wert wirklich in in einem Bruchteil $\\alpha$ aller Konfidenzintervalle liegt.\n", "\n", "**Hinweise:** Lesen Sie alle Aufgaben zuerst durch, planen Sie die Programmstruktur sorgfältig und achten Sie auf einen modularen Aufbau der einzelen Teile, damit Sie die notwendigen Schritte in einer Schleife ausführen und die gesamte Studien ggf. mit unterschiedlichen Parametereinstellungen wiederholen können. In den Codebeispielen zur Vorlesung finden Sie eine Reihe von Lösungsansätzen. Auch in den Tutorials zu Likelihood-Anpassungen in Jupyter [*negLogLFits.ipynb*](https://git.scc.kit.edu/yh5078/datenanalyse/-/blob/master/jupyter/negLogLFits.ipynb) und in Beispiel-Python-Skripten zu PhyPraKit wie [*PhyPraKit/examples/toyMC_Fit.py*](https://git.scc.kit.edu/yh5078/PhyPraKit/-/blob/master/examples/toyMC_Fit.py) von Prof. Quast finden Sie eine Vorlage zum Aufsetzen der Programmstruktur. Wenn Sie die Maximum-Likelihood-Anpassung nicht selbst programmieren möchten, können Sie die Funktion *PhyPraKit.phyFit.mFit()* dazu verwenden, die in der [PhyPraKit-Dokumentation](https://etpwww.etp.kit.edu/~quast/PhyPraKit/htmldoc/) beschrieben ist." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# notwendige Importe: numpy, matplotlib.pyplot, Fit Pakete\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "# -> eigenen Code hier einfügen " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Aufgaben" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1) Schreiben Sie eine Funktion, die 50 exponentiell verteilte Zufallszahlen im sensitiven Detektorintervall $[t_{\\rm min}, t_{\\rm max}]$ erzeugt." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Parameter für die Ausführung der Studie als globale Variable\n", "Nexp = 3000\n", "N = 50\n", "tau = 2.\n", "tmin = 1.\n", "tmax = 10.\n", "\n", "npar = 1 # Zahl der angepassten Parameter\n", "pnams = [\"tau\"] # Liste mit Namen der Parameter\n", "true_vals= np.array([tau, tmin, tmax]) # die \"Wahren Werte\"" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# definieren Sie hier die Verteilungsdichte \n", "## def exponentialDecayPDF(t, tau= tau, tmin=tmin, tmax=tmax):\n", "\n", "# -> eigenen Code hier einfügen \n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Funktion zur Erzeugung der Daten\n", "## def generateExpData(N, tau, tmin, tmax):\n", "\n", "# -> eigenen Code hier einfügen \n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "2) Setzen Sie eine Anpassung mit einer ungebinnten negativen log-Likelihood-Funktion auf. Dazu können Sie das Beispiel aus der Vorlesung und aus einem früheren Übungsblatt nutzen. Alternativ können Sie die Funktion `mfit()` aus dem Paket `PhyPrakit.phyFit` nutzen." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Schleife zur wiederholten Ausführung von Datenerzeugung, Anpassung und zum Speichern der Ergebnisse\n", "## def MC_loop():\n", "\n", "# -> eigenen Code hier einfügen \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "3) Führen Sie 1) und 2) in einer Schleife aus; das Ziel ist es, 3000 Pseudoexperimente zu simulieren. Beginnen Sie zum Testen des Codes aber zunächst mit einer kleineren Anzahl! Denken Sie daran, in der Schleife die zur Bestimmung von Erwartungstreue und Abdeckung notwendigen Daten in einem Array zu speichern. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# run MC loop\n", "# -> eigenen Code hier einfügen \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "4) Analysieren Sie die in der Monte Carlo-Schleife gewonnenen Daten und geben Sie Erwartungstreue und Abdeckung für die geschätzte Lebensdauer $\\hat{\\tau}$ und deren Unsicherheitsintervall $[\\hat{\\tau}-\\Delta^-, \\hat{\\tau}+\\Delta^+]$ an. Denken Sie daran, dass Sie auch die statistische Signifikanz der so bestimmten Eigenschaften sicherstellen müssen. Schätzen Sie daher die Unsicherheiten der ermittelten Erwartungstreue und Abdeckung ab. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Berechnung und Ausgabe der Ergebnisse\n", "\n", "# -> eigenen Code hier einfügen \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "5) Bewerten Sie Ihre Ergebnisse: \n", " a) Wird eine statistisch signifikante Verzerrung nachgewiesen? \n", " b) Ist sie ggf. im Bezug auf die Unsicherheiten der Einzelmessungen relevant? \n", " c) Benötigen Sie auf Grund der festgestellten Unter- oder Überabdeckung eine Korrektur des für eine Einzelmessung relevanten Konfidenzintervalls ? " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Antworten\n", "*Hier Antworten eintippen*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "6) **Freiwillig**: In der Vorlesung haben Sie die \"Bootstrap-Methode\" kennen gelernt, um die Eigenschaften von Funtionen von Zufallszahlen zu untersuchen. Diese Methode lässt sich auch auf das hier behandelte Problem anwenden. Schreiben sie dazu eine Funktion, die aus einem einzigen initialen, mit der unter 1) geschriebenen Funktion erzeugten Daten neue Datensätze durch \"Ziehen mit Zurücklegen\" erzeugt. Dazu können Sie die Methode `rng.choice(data, size=N)` verwenden. Führen Sie nun wiederum die Schritte 2) - 4) durch und vergeleichen Sie das Ergebnis. Sie sollten diese Studie mit mehreren verschiedenen initialen Datensätzen wiederholen, da im Einzelfall, abhängig vom initialen Datensazt, größere Fluktuationen auftreten." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Funktion zum Re-Sampling der Daten (für Bootstrap)\n", "#def resampleData(N, dT0):\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Einen Datensatz erzeugen\n", "# run MC loop\n", "\n", "# Daten mit Resampling (\"bootstrapping\")\n" ] } ], "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 }