{ "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 04 \n", "Abgabe: Montag 4.12.2023 bzw. Dienstag 5.12.2023 \n", "\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Auf dem vierten Übungsblatt beschäftigen Sie sich vertieft mit der Parameterschätzung und wenden einen Hypothesentest an, um die Gemeinsamkeit zweier Datensätze zu überprüfen. \n", "\n", "Auf dem vergangenen Übungsblatt haben Sie die Parameterschätzung selbst implementiert, indem Sie eine Kostenfunktion (*Likelihood*) geschrieben, diese an diskreten Stellen eines vorgegebenen Parameterraumes ausgewertet und das Minimum an einer der Stützstellen berechnet haben. Wie Sie aus der Vorlesung oder Ihrer eigenen Erfahrung kennen, weist das Vorgehen an vielen Stellen Optimierungsbedarf auf. Sie könnten viel mehr Stützstellen verwenden, um einen besseren Scan der Likelihood zu erhalten, dann müssen Sie aber viel mehr rechnen, zudem müssen Sie Unsicherheiten berücksichtigen, das Verfahren für alle zukünftigen Herausforderungen verallgemeinern und alles so niederschreiben, dass der Code lesbar und einfach verständlich ist. An dieser Stelle können wir Sie beruhigen, denn all diese Aufgaben haben bereits eine Vielzahl an Menschen für Sie erledigt.\n", "\n", "Im Folgenden beschäftigen Sie sich mit einem Minimierungsverfahren, welches Ihnen erlaubt deutlich effizienter ein Minimum in einer Parameterlandschaft zu finden, anschließend sind Sie eingeladen mithilfe von [*kafe2*](http://etpwww.etp.kit.edu/~quast/kafe2/htmldoc/) die Parameterschätzung vollständig durchzuführen. Zum Abschluss lernen Sie den *Student'schen-t-Test* kennen." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "# Aufgabe 1: Simulated Annealing \n", "---\n", "\n", "Wie Sie in Ihrem Studium möglicherweise bereits erkannt haben, ist das Finden eines Optimums einer hochdimensionalen Kostenfunktion $F(x)$ eine häufige Fragestellung. Damit Sie eine Vorstellung der Lösungsansätze erhalten, sollen Sie in dieser Aufgabe selbst einen Algorithmus schreiben, der Ihnen dabei hilft, jenes Optimum zu bestimmen. Der Algorithmus ist bei Weitem nicht der beste oder schnellste, verdeutlich Ihnen aber den Einsatz von Monte-Carlo-Methoden in der Modernen Physik.
\n", "Als Beispiel betrachten Sie die zweidimensionale, modifizierte Rosenbrock Funktion\n", "\n", "$$\n", "F(x,y)=\\left(x^2+y-a\\right)^2+\\left(x+y^2-b\\right)^2+c\\cdot\\left(x+y\\right)^2\n", "$$\n", "\n", "mit drei Parametern $a, b$ und $c$. Die Funktion weist einige lokale Minima, aber immer nur ein globales Minimum auf. Die Bearbeitung der Aufgabe erfolgt für die zufällige Wahl an Startparametern: $a=11,\\, b=7,\\, c=0.1$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## a) Implementierung des Algorithmus\n", "In der ersten Teilaufgabe erstellen Sie das Grundgerüst der Minimierung, indem Sie zunächst die modifizierte Rosenbrock-Funktion und anschließend das simulierte Abkühlen als Funktionen definieren. Zur Regulierung der Temperatur wird eine konstante Kühlrate verwendet, die als Parameter des Algorithmus betrachtet wird. Zusätzlich wird eine Schrittweite als Parameter eingeführt, welche die Distanz eines neuen Zustandes zum alten beeinflusst." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import matplotlib.animation\n", "from IPython.display import Image\n", "# initialize random number generator\n", "rng = np.random.default_rng( 42 )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# ToDo: implement modified rosenbrock function\n", "def modified_rosenbrock_function(x, par):\n", " xValue = x[0] \n", " yValue = x[1] \n", " a = par[0] \n", " b = par[1] \n", " c = par[2] \n", " \n", " return xValue**2+yValue**2" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def simulated_annealing(\n", " init_vals=[0,0], \n", " rosenbrock_pars=[0, 0, 0],\n", " init_temp=100, \n", " final_temp=1, \n", " cool_speed=1, \n", " step_size=1\n", "):\n", " \"\"\"\n", " Minimize the modified Rosenbrock function using simulated annealing.\n", " \n", " params:\n", " init_vals: Initial x and y values.\n", " rosenbrock_pars: Parameters of the modified Rosenbrock function.\n", " init_temp: Initial temperature the cooling starts from.\n", " final_temp: Final temperature of the cooling.\n", " cool_speed: Cooling speed in percent of the current temperature.\n", " step_size: Step size used in the cooling procedure.\n", " \n", " returns:\n", " min_pars: List of floats.\n", " List of the x and y values at the found minimum.\n", " listOfPoints: List of floats.\n", " List of the visited points during the minimization process.\n", " \"\"\"\n", " nParameter = 2 # 2 parameters: x and y\n", " if len(init_vals) != nParameter:\n", " raise Exception(\"Number of function parameters does not correspond to given number of initial values.\"\n", " f\"Expected number of initial values: {nParameter}\")\n", " \n", " # =======================================================================================================\n", " # Presetting and book keeping\n", " # =======================================================================================================\n", "\n", " # parameters of the algorithm\n", " initialTemperature = init_temp\n", " finalTemperature = final_temp\n", " coolingSpeed = cool_speed # in percent of current temperature --> defines number of iterations\n", " stepSize = step_size \n", " \n", " # current parameters and cost\n", " currentParameters = init_vals \n", " currentFunctionValue = modified_rosenbrock_function(currentParameters, rosenbrock_pars)\n", "\n", " # keep reference of best parameters\n", " bestParameters = currentParameters\n", " bestFunctionValue = currentFunctionValue\n", "\n", " # log keeping: keep track of path\n", " listOfPoints= []\n", " \n", " # =======================================================================================================\n", " # =======================================================================================================\n", "\n", " # Heat the system\n", " temperature = initialTemperature\n", "\n", " # Start to slowly cool the system\n", " while (temperature > finalTemperature): \n", "\n", " # ToDo: Change parameters\n", " newParameters = rng.standard_normal(nParameter)\n", " \n", " # Get the new value of the function\n", " newFunctionValue = modified_rosenbrock_function(newParameters, rosenbrock_pars)\n", "\n", " # ToDo: Compute Boltzman probability\n", "\n", " # ToDo: Acceptance rules\n", "\n", " \n", " # Delete this as soon as acceptance rules are implemented\n", " bestParameters = currentParameters\n", " listOfPoints.append(currentParameters)\n", " \n", " # Cool the system\n", " temperature *= (1 - coolingSpeed/100.)\n", "\n", "\n", " # end of cooling loop\n", " return bestParameters, listOfPoints" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## b) Finden des globalen Minimums\n", "Testen Sie zunächst Ihren Algorithmus aus, indem Sie Startwerte für $x$ und $y$ in der Nähe des globalen Minimums wählen $(x,y)=(-2, 3)$. Machen Sie sich dabei mit den Eigenschaften der Anfangs- und Endtemperaturen, der Abkühlrate und der Schrittgröße vertraut.\n", "\n", "Sobald Sie mit Ihrem Algorithmus zufrieden sind, dieser also zuverlässig gegen das globale Minimum konvergiert, wählen Sie als Startpunkt einen Punkt in einem der Nebenminima (z.B. $(x,y)=(3, -2)$) aus und passen Sie die Parameter des Algorithmus so an, dass er wieder zuverlässig gegen das globale Minimum konvergiert.
\n", "Zur korrekten Einstellung der Parameter bietet es sich an, entweder einen Bereich an Parametern systematisch zu untersuchen oder den Einfluss eines jeden Parameters graphisch zu untersuchen. Sie finden einen Vorschlag zur graphischen Untersuchung in den folgenden Zellen. In beiden Fällen bietet es sich jedoch an, sich die Eigenschaften der Algorithmusparameter vorher zu überlegen.\n", "\n", "#### Hinweise und Überlegungen:\n", "* Die Anfangs- und Endtemperatur sollen die Kostenfunktionswerte der Anfangsbedingung und des gesuchten Minimums widerspiegeln.\n", "* Die Differenz zwischen Anfangs- und Endtemperatur und die Kühlrate beeinflussen die Anzahl an Iterationen.\n", "* Die Temperaturskala beeinflusst die Wahrscheinlichkeit für Sprünge.\n", "* Die Schrittgröße sollte adäquat im Bezug zum Abstand der Minima gewählt werden. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "a = 11\n", "b = 7\n", "c = 0.1\n", "\n", "\n", "initial_values = [-2, 3]\n", "initial_temp, final_temp, cool_speed, step_size = 100, 1, 1, 1\n", "\n", "bestParameters, listOfPoints = simulated_annealing(init_vals=initial_values, \n", " rosenbrock_pars=[a, b, c],\n", " init_temp=initial_temp, \n", " final_temp=final_temp, \n", " cool_speed=cool_speed, \n", " step_size=step_size\n", " )\n", "\n", "minValue = modified_rosenbrock_function(bestParameters, [a, b, c])\n", "\n", "print(f\"Resultat für das Minimum: f(x,y)={minValue} bei x={bestParameters[0]}, y={bestParameters[1]}\")\n", " \n", "# check if minimum is the global one (then the mimimum value should be around 0.01): \n", "if minValue < 0.1:\n", " print(\"Ja, das ist das globale Minimum!\")\n", "else:\n", " print(\"Oh nein, das ist nicht das globale Minimum, sondern nur ein lokales!\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%matplotlib notebook" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Points = np.array(listOfPoints).T\n", "frames = 30\n", "total = len(Points[0])\n", "\n", "\n", "def init():\n", " #p1.set_zscale('log')\n", " p1.set_xlim(-5, 5)\n", " p1.set_ylim(-5, 5)\n", " p1.set_zlim(0.0001, 1000)\n", " \n", " p1.set_xlabel('X')\n", " p1.set_ylabel('Y')\n", " p1.set_zlabel('Z')\n", " \n", " # plot function\n", " X = np.linspace(-5, 5, 1000)\n", " Y = np.linspace(-5, 5, 1000)\n", " X, Y = np.meshgrid(X, Y)\n", "\n", " Z = modified_rosenbrock_function([X, Y], [a, b, c])\n", " p1.plot_surface(X, Y, Z, cmap='coolwarm', alpha=0.4)\n", " p1.scatter(bestParameters[0], bestParameters[1], modified_rosenbrock_function(bestParameters, [a, b, c]), c=\"k\", label=\"foo\")\n", "\n", "\n", "def update(num, data, plotcom): \n", " low = int(num/frames * len(Points[0]))\n", " up = int((num+1)/frames * len(Points[0]))\n", " \n", " plotcom.set_data(data[0:2, low:up])\n", " plotcom.set_3d_properties(data[2, low:up])\n", "\n", "\n", "\n", "fig=plt.figure(figsize=(10, 10))\n", "p1=fig.add_subplot(1, 1, 1, projection='3d')\n", "\n", "data = np.array([Points[0], Points[1], modified_rosenbrock_function(Points, [a,b,c])])\n", "plotcommand = p1.plot(data[0], data[1], data[2], color='blue', marker=\".\", linestyle=\"\")[0]\n", "\n", "# animate\n", "animation = matplotlib.animation.FuncAnimation(fig, update, frames=frames, init_func=init, fargs=(data, plotcommand), interval=250, blit=False)\n", "\n", "animation.save('annealing.gif', writer=matplotlib.animation.PillowWriter(fps=5))\n", "\n", "from IPython.display import HTML\n", "HTML(animation.to_jshtml())\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "# Aufgabe 2: Anpassung einer Geraden \n", "---\n", "\n", "Die zweite Aufgabe auf diesem Übungsblatt wird Ihnen sehr bekannt vorkommen, da Sie im Anfängerpraktikum garantiert eine Messung dieser Art durchgeführt haben. In der Datei *StromSpannungsKennlinie.csv* finden Sie 40 Strom- und Spannungsmesswerte. Diese Daten könnten aus einer beliebigen Kalibrationsmessung stammen. Ihre Aufgabe ist, aus diesen Strom- und Spannungswerten den elektrischen Widerstand des untersuchten Bauteils zu bestimmen.\n", "\n", "Der Kernaspekt der Aufgabe liegt darin, die Anpassung korrekt durchzuführen und Ihnen zu verdeutlich, wie einfach eine Anpassung, unter der Berücksichtgung aller typischen Unsicherheiten, ist. Es wird Ihnen nahegelegt, hier *kafe2* als Werkzeug zu verwenden, da dieses *Python*-Paket viele wichtige Funktionen beinhaltet, die sonst durch mühseelige Handarbeit in die Anpassung eingebaut werden müssten.\n", "\n", "Im Folgenden betrachten Sie einen Datensatz aus je 20 gemessenen Strom- und Spannungswerten aus zwei Messbereichen: Der erste Messbereich umfasst die Spannung $U\\in[0,2]$V und den Strom $I\\in[0,20]$mA, der zweite die Bereiche $U\\in[2,20]$V und $I\\in[20,200]$mA. Für die Unsicherheiten der Messungen können Sie sich auf die Anzeige des Messinstruments und die Herstellerangaben beziehen:\n", "* Jeder Messbereich umfasst eine Angabe von fünf Stellen mit einer **absoluten Unsicherheit** von **zwei Stellen**. Messen Sie also im ersten Messbereich die Spannung, wird Ihnen der Bereich *x,xxxx*V angezeigt, die Unsicherheit bezieht sich dann auf die letzten zwei Stellen. \n", "* Alle Strom- und Spannungsmessungen unterliegen einem **konstanten Untergrundrauschen**, bei der Spannung gibt der Hersteller des verwendeten Messgerätes **20mV**, beim Strom **1mA** an. \n", "* Zusätzlich wird eine **relative, vollständig korrelierte Unsicherheit** von **5%** auf alle Messwerte angegeben." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## a) Vorbereitung der Auswertung\n", "\n", "Bereiten Sie im ersten Aufgabenteil die Auswertung der Messdaten vor, indem Sie \n", "* eine lineare Funktion definieren, die Ihnen als Modellfunktion bei der Anpassung dienen wird,\n", "* die Messdaten korrekt einlesen\n", "* und alle beschriebenen Unsicherheiten korrekt aufschreiben und gegebenenfalls verrechnen.\n", "\n", "Beachten Sie dabei die Aufteilung der Unsicherheiten auf die Messbereiche. Es bietet sich beispielsweise an, Unsicherheiten mithilfe eines Arrays zu definieren. In der folgenden Teilaufgabe haben Sie die Möglichkeit, entweder Fleißkommazahlen (*float*) oder Arrays der Länge der Messdaten anzugeben." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import PhyPraKit as ppk\n", "from kafe2 import XYContainer, Fit, Plot, ContoursProfiler" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# linear model\n", "\n", "\n", "# uncertainty components\n", "\n", "\n", "# data read in\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## b) Durchführung der Anpassung\n", "\n", "Führen Sie in dieser Teilaufgabe die Anpassung der Geraden durch. Beispiele, wie die Anpassung mithilfe von *kafe2* funktioniert, finden Sie in der [Dokumentation zu *kafe2*](http://etpwww.etp.kit.edu/~quast/kafe2/htmldoc/index.html).\n", "\n", "Der Einfluss der relativen, korrelierten Unsicherheit wird zu einer Verzerrung der Anpassung führen. Führen Sie dafür die Anpassung einmal mit dem Bezug der relativen Unsicherheit auf die Daten und einmal mit Bezug auf das Modell durch und vergleichen Sie beide Ergebnisse.\n", "\n", "#### Hinweise:\n", "* Die Daten übergeben Sie am einfachsten an *kafe2* mithilfe der Klasse *XYContainer(x, y)*.\n", "* Der Fit lässt sich mithilfe der Klasse *Fit(container)* erstellen.\n", "* Die Unsicherheiten können Sie mithilfe der Methode [*add_error(err_val)*](http://etpwww.etp.kit.edu/~quast/kafe2/htmldoc/parts/api_documentation/index.html?highlight=add_error#kafe2.fit.indexed.IndexedContainer.add_error) der Klasse *Fit* definieren.\n", "* Die Argumente *relative*, *correlation* und *reference* (*'data'* oder *'model'*) sind notwendig zur korrekten Fehlerbehandlung in der Methode *add_error()*.\n", "* Verbessern Sie die Lesbarkeit des Plots mithilfe von geschickt ausgewählten Labeln und Farben im Plot.\n", "* Beim Plotten mithilfe der Klasse *Plot()* bietet es sich an, *asymmetric_parameter_errors=True* zu setzen, um asymmetrische Unsicherheiten zu erhalten. Dies ist äquivalent zu dem Vorgehen, mit dem Sie auf dem vorigen Blatt Konfidenzintervalle berechnet haben.\n", "\n", "##### Zusatzaufgabe (ohne Wertung)\n", "* Betrachten Sie die Likelihood-Konturen der Parameter in beiden Fällen der Anpassung. Wo sehen Sie den Einfluss der unterschiedlichen Behandlung der relativen Unsicherheiten?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# prepare data for fit\n", "\n", "\n", "# perform first fit: rel. uncertainty to data\n", "\n", "\n", "# add uncertainties to first fit\n", "\n", "\n", "# perform first fit\n", "\n", "\n", "# perform second fit: rel. uncertainty to model\n", "\n", "\n", "# add uncertainties to second fit\n", "\n", "\n", "# perform second fit\n", "\n", "\n", "# plot fits\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "# Aufgabe 3: Vergleich von zwei Verteilungen durch *Hypothesentests* \n", "---\n", "\n", "In der letzten Aufgaben machen Sie sich an einem einfach Beispiel klar, wie ein sehr simpler Hypothesentest Ihnen dabei helfen kann, eine von Ihnen geforderte Entscheidung zu treffen. Als Beispiel betrachten Sie zwei Datensätze, *sample1.dat* und *sample2.dat*. Sie sollen die Aussage treffen, ob es möglich ist, dass die Datensätze aus einem gemeinsamen Datensatz stammen. Als Mittel der Überprüfung können Sie einen sogenannten verallgemeinerten [*Student'schen-t-Test*](https://de.wikipedia.org/wiki/Zweistichproben-t-Test#Welch-Test) durchführen.\n", "\n", "Sie gehen von der Nullhypothese $H_0$ **Die Stichproben haben gleiche Erwartungswerte $\\mu_1=\\mu_2$** aus. Zur Quantifizierung des Unterschieds wird eine Prüfgröße (engl.: test statistics) $d$ berechnet, die auf einen möglichen Unterschied der beiden Proben empfindlich ist. Für einen solchen Test wird der Betrag der Differenz der Erwartungswerte verwendet, die auf die erwartete Unsicherheit normiert werden:\n", "\n", "\n", "\\begin{align*}\n", "d &= \\frac{\\left|\\mu_1-\\mu_2\\right|}{\\sigma}, \\\\\n", "\\sigma &= \\sqrt{\\frac{\\sigma_1^2}{N_1}+\\frac{\\sigma_2^2}{N_2}}.\n", "\\end{align*}\n", "\n", "\n", "Dabei geben $\\sigma_{1,2}$ bzw. $N_{1,2}$ die Standardabweichung bzw. Größe der Stichproben an." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## a) Vorbereitung des Hypothesentests\n", "\n", "Implementieren Sie eine Funktion zur Berechnung der Größe $d$ und der Anzahl der Freiheitsgrade $n_\\mathrm{dof}=N_1+N_2-2$. Lesen Sie die Stichproben ein und histogrammieren Sie diese." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Implement test statistic\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# data read in\n", "\n", "\n", "# plot both distributions\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## b) *Studentscher-t-Test*\n", "\n", "Berechnen Sie nun die beobachtete Größe $d_\\mathrm{obs}$. Stellen Sie die *t*-Verteilung graphisch dar und zeichnen Sie den beobachteten Wert ein. Stellen Sie ebenfalls fest, ob Sie die Nullhypothese mit einer Irrtumswahrscheinlichkeit von $\\alpha=0.045 (2\\sigma)$ verwerfen können.\n", "\n", "#### Hinweise:\n", "* Die *SciPy*-Methode *[scipy.stats.t](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.t.html).pdf()* ist eine große Hilfe bei der Darstellung der *t*-Verteilung.\n", "* Die Überprüfung der Nullhypothese kann entweder über die Bestimmung des $p$-Wertes der beobachteten Größe $d_\\mathrm{obs}$ oder graphisch durch Einzeichnung von $d_0$ geschehen. In beiden Fällen bietet *[scipy.stats.t](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.t.html)* Methoden (*cdf*, *sf*) an, die bei der Berechnung hilfreich sind." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import scipy.stats as stat\n", "import scipy.integrate as integrate" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# calculate observed quantity\n", "\n", "\n", "# accept or reject 0 hypothesis\n", "alpha = 0.045\n", "\n", "\n", "# plot t distribution\n", "\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.11.5" } }, "nbformat": 4, "nbformat_minor": 4 }