{ "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 02
\n", "Abgabe: Montag 20.11.2023 bzw. Dienstag 21.11.2023\n", "\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Gemäß der Abfolge des Stoffs in der Vorlesung üben Sie auf dem zweiten Blatt den Umgang mit Monte-Carlo-(MC)Methoden und die Erstellung von Monte-Carlo-Simulationen. Die MC-Simulationen sind ein wesentlicher Bestandteil der modernen Physik und kommen in jedem Fachgebiet zum Einsatz – immer dann wenn es nicht möglich oder ineffiezient ist, ein Problem mit vielen gekoppelten Freiheitsgraden analytisch oder mit anderen numerischen Methoden zu lösen. Beispiele für solche Probleme sind über die in der Vorlesung hinaus erwähnten Anwendungsfällen z.B. die Faltung von Proteinen, die Entwicklung von Galaxien oder Teilchenkollisionen\n", "in einem Teilchenbeschleuniger.
\n", "\n", "### Erinnerung: Grundsätzliches Verfahren der Monte Carlo Methoden\n", "1. Erzeuge Folge von Zufallszahlen $r_1, r_2, \\dots, r_m$, gleichverteilt in $[0,1[$.\n", "2. Transformiere diese Foge zur Erzeugung einer anderen Sequenz $x_1, x_2, \\dots, x_n$, die\n", " einer beliebigen anderen Verteilugsdichte $f(x)$ folgen. (Anmerkung: $x, x_i$ können auch Vektoren sein.)\n", "3. Aus den $x_i$ werden dann Eigenschaften von $f(x)$ bestimmt. Zum Beispiel entspricht der Anteil der $x_i$ mit $a\\leq x_i\\leq b$ näherungsweise dem Integral $\\int_{a}^{b}f(x)\\mathrm{d}x$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "# Aufgabe 1: Verwerfungsmethode \n", "---\n", "\n", "Im ersten Aufgabenteil beschäftigen Sie sich mit einem einfachen Verfahren zur Erzeugung einer gewünschten Sequenz $x_1, x_2, \\dots, x_n$ aus gleichverteilten Zufallszahlen, der **Verwerfungsmethode**. In der Vorlesung haben Sie bereits das wohl bekannteste Beispiel, die Berechnung der Kreiszahl $\\pi$, kennengelernt. In dem Beispiel werden Zufallszahlen $(r_1,r_2)$ gleichmäßig in der Einheitsfläche erzeugt (durch *np.random.random(size=2)*) und dann mithilfe der Bedingung, dass das Zufallszahlenpaar innerhalb des Einheitskreises liegt ($r_1^2+r_1^2\\leq 1^2$), die Kreiszahl $\\pi$ approximiert.\n", "\n", "Mit Hilfe der Verwerfungsmethode können recht einfach auch hochdimensionale, bestimmte Integrale nährungsweise numerisch berechnet werden. In dieser Aufgabe sollen Sie die Methodik des Beispiels nutzen, um das Volumen $V_d$ von $d$-dimensionalen Kugeln in Abhängikeit von der Anzahl der Dimensionen $d$ zu berechnen." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## a) Berechnung der Einheitsvolumen\n", "Berechnen Sie für $d\\in[1,10]$ jeweils das Volumen der $d$-dimensionalen Einheitskugel. \n", "\n", "Erzeugen Sie hierfür für jeden Schritt 10'000 Zufallszahlen gleichverteilt im entsprechenden $d$-dimensionalen Einheits-Hyperwürfel. Berechnen Sie mithilfe einer geeigneten Bedingung das Volumen des eingeschlossenen Kugelsegments im Einheits-Hyperwürfel und anschließend das volle Volumen der $d$-dimensionalen Einheitskugel.\n", "\n", "> Hinweis: Vergewissern Sie sich, welche *NumPy* Methode in [*np.random*](https://numpy.org/doc/stable/reference/random/legacy.html#simple-random-data) hier nützlich sein kann." ] }, { "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": [ "# set dimensions\n", "dmax = 10 # maximum\n", "dim = np.arange(1., dmax, dtype=int)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# implement MC method\n", "npoints = 10000 # number of d-dimensional points to test\n", "Vmc = [] # array holding results\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## b) Abschätzung der statistischen Unsicherheit und Darstellung des Ergebnisses\n", "\n", "Bestimmen Sie die statistische Unsicherheit $\\sigma_{V_d}$ der Ergebnisse. Nutzen Sie dazu aus, dass bei $N$ Versuchen die Anzahl $N_{in}$ der innerhalb der Kugel liegenden Punkte einer Binomialverteilung $P(N_{in};p,N)$ folgt. Die Unsicherheit ergibt sich aus der Varianz der gemessenen Größe $\\sigma^2=Var[\\hat{p}]$, sowie der Varianz unserer binomialverteilten Zufallsvariable $N_{in}$ und der Näherung $p\\simeq\\hat{p}=\\frac{N_{in}}{N}$.\n", "\n", "> Hinweis: Konstanten können aus der Varianz $Var[x]$ herausgezogen werden, siehe dazu [Rechenregeln der Varianz](https://de.wikipedia.org/wiki/Varianz_(Stochastik)#Lineare_Transformation).\n", ">\n", "> Die Varianz einer binomialverteilten Zufallsvariable haben Sie auf dem ersten Blatt bereits verwendet.\n", "> \n", "> Vergessen Sie nicht, die Unsicherheit auf $\\hat{p}$ so zu skalieren, dass Sie die Unsicherheit $\\sigma_{V_d}$ erhalten." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# uncertainty on MC results using binomial statistics\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Tragen Sie Ihre Ergebnisse für $V_d$ und die Unsicherheit $\\sigma_{V_d}$ in ein Diagramm ein ([plt.errorbar()](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.errorbar.html) wird hier nützlich sein). Die Funktion *VofSphere(d, r=1)* gibt das Volumen von $d$-dimensionalen Kugeln zurück, wie es sich mit mathematischen Integrationsmethoden ergibt. Tragen Sie auch das \"theoretische\" Ergebnis in die Grafik ein." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def VofSphere(d, r=1):\n", " \"\"\"\n", " calculate volume of a d-dimensional sphere with radius r\n", "\n", " Args:\n", " d: dimension\n", " r: radius\n", "\n", " Returns:\n", " (array of) float: d-dimensional Volume\n", " \"\"\"\n", " V = r**d * np.pi**(d/2) / sp.gamma(d/2 + 1)\n", " return V\n", "\n", "# Volume from mathematical formula\n", "\n", "\n", "# plot results (with matplotlib)\n", "\n", "\n", "\n", "plt.xlim(0, dmax+1.)\n", "plt.ylim(1.5, 6.)\n", "plt.xlabel('Dimension', fontsize='x-large')\n", "plt.ylabel('Volumen der Kugel', fontsize='x-large')\n", "plt.legend(loc='best')\n", "plt.title(\"Volumen der d-dimensionalen Kugel\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "# Aufgabe 2: Transformationsmethode \n", "---\n", "\n", "In der zweiten Aufgaben beschäftigen Sie sich mit der Transformationsmethode. In der ersten Teilaufgabe werden Sie den Zusammenhang zwischen der Verteilung von Zufallszahlen und den dazu gehörigen Wahrscheinlichkeitsdichten vertiefen. Im zweiten Aufgabenteil werden Sie dann selbst korrelierte Zufallszahlen erzeugen, da dies in den üblichen *Python* Bibliotheken nicht zu finden ist, die Notwendigkeit jedoch in der Praxis häufig gegeben ist." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## a) Erzeugung von einfachen Verteilungen\n", "Sie arbeiten mit den gleichverteilten Zufallszahlen $r\\in\\text{U}_{[0,1[}$. Im Folgenden betrachten Sie zwei Fälle:\n", "* Sie kennen die Verteilung $t_i=t(r_i)$ und sollen die zugrunde liegende Wahrscheinlichsdichte $g(t)$ bestimmen.\n", "* Sie kennen die Wahrscheinlichkeitsdichte $g(t)$ und sollen eine Verteilung $t_i=t(r_i)$ berechnen.\n", "\n", "Erzeugen Sie in den folgenden Teilaufgaben immer 10'000 gleichverteilte Zufallszahlen im Intervall $[0,1[$ und bestimmen Sie die fehlende Größe mithilfe der Transformationsmethode. Zeichen Sie anschließend in ein Diagramm die Verteilung $t(r)$ (Histogramm) und die Wahrscheinlichkeitsdichte $g(t(r))$ ein." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### i) Gegeben: $g(t)=\\frac{1}{t}$ – Gesucht: $t(r)$ für $t=1\\dots 2.75$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# generate histogramm entries\n", "\n", "\n", "# compute PDF\n", "\n", "\n", "# plot histogramm and pdf\n", "\n", "\n", "# make plot nice\n", "plt.xlabel('$t(r)$', fontsize='x-large')\n", "plt.ylabel('$g(t(r))$')\n", "plt.legend(loc='upper right')\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### ii) Gegeben: $t(r)=r^2$ – Gesucht: $g(t)$ für $t=0.01\\dots 1$ " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# generate histogramm entries\n", "\n", "\n", "# compute PDF\n", "\n", "\n", "# plot histogramm and pdf\n", "\n", "\n", "# make plot nice\n", "plt.xlabel('$t(r)$')\n", "plt.ylabel('$g(t(r))$')\n", "plt.legend(loc='upper right')\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## b) Erzeugung korrelierter, normalverteilter Zufallszahlen\n", "Dieser Aufgabenteil soll Ihnen bewusst machen, wie Sie aus unkorrelierten normalverteilten Zufallszahlen korrelierte Zufallszahlen erzeugen. Die Herleitung des Prinzips wird in der 3. Vorlesung diskutiert und soll hier auf den einfache Fall von zwei\n", "gaußverteilten Zufallszahlen angewandt werden.\n", "\n", "Erzeugen Sie zwei Sätze aus 10'000 standard-normalverteilten Zufallszahlen $\\vec{u}_{1,2}$ ($\\mu_{1,2}=0$ und $\\sigma_{1,2}=1$) und transfomieren Sie diese zu korrelierten Zufallszahlen $\\vec{x}_{1,2}$ mit einem Korrelationskoeffizienten von $\\rho=0.5$.\n", "> Hinweis: Die Methode *np.dot()* kann bei einem Matrix-Vektor Produkt hilfreich sein." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# implement transformation\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Stellen Sie nun die beiden Sätze an Zufallszahlen geeignet dar. Es bietet sich an, für jedes Paar jeweils ein zweidimensionales Histogramm zu erstellen. Berechnen Sie anschließend den Korrelationskoeffizienten aus dem neu erzeugten Satz an Zufallszahlen, um Ihr Vorgehen zu überprüfen.\n", "> Hinweis: Das Paket *PhyPraKit* enthält einige hilfreiche Methoden zur Bearbeitung der Aufgabe. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot randomnumbers and calculate correlationcoefficient\n", "import PhyPraKit as ppk\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "# Aufgabe 3: Majoranten-Methode \n", "---\n", "\n", "In der letzten Aufgabe kombinieren Sie Ihr erlangtes Wissen aus der Transformations- und der Verwerfungsmethode. Die Kombination aus beiden Methoden wird häufig als Majorantenmethode (*engl.: Importance Sampling*) bezeichnet, welche eine deutlich effizientere Methode zum Erzeugen einer benötigten Verteilung ist. \n", "\n", "Im Folgenden betrachten Sie eine Verteilung, die Ihnen in ähnlicher Form möglicherweise im Fortgeschrittenenpraktikum begegen wird:\n", "$$\n", "f(x)=\\sin^2(\\pi x)e^{-x}\n", "$$\n", "im Intervall $[0,2\\pi]$. \n", "\n", "Dazu erzeugen Sie zunächst mithilfe der Transformationsmethode Zufallszahlen gemäß der Verteilung \n", "der Majorante und wählen dann mit der Verwerfungsmethode davon eine passende Anzahl aus, um die \n", "gewünschten Verteilung zu erhalten. \n", "\n", "## a) Transformationsmethode zum Erhalt der Majorante\n", "\n", "Erzeugen Sie mithilfe der Transformationsmethode exponentiell verteilte Zufallszahlen als Majorante\n", "\n", "$$\n", "m(x)=e^{-x}.\n", "$$\n", "\n", "Überzeugen Sie sich zunächst davon, dass Ihre Majorante die gesuchte Verteilung vollständig umschließt, indem Sie die Majorante $m(x)$, die gewünschte Verteilung $f(x)$ und die erzeugten Zufallszahlen histogrammiert in einer Abbildung darstellen.\n", "> Hinweis: Es bietet sich an, die Verteilungen und Verteilungsdichten als Funktion zu implementieren, da Sie sie häufiger aufrufen und gegebenenfalls im Praktikum wiederverwenden werden können." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# build functions f(x), m(x) and PDF of m(x)\n", "\n", "\n", "# sample majorant\n", "\n", "\n", "# plot functions and random number\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## b) Verwerfungsmethode zum Erhalt der gewünschten Verteilung\n", "\n", "Wenden Sie nun die Verwerfungsmethode an, um die gewünschte Verteilung zu erhalten. Erzeugen Sie für jede der Zufallszahlen $x_i$ aus der Majorante eine gleichverteilte Zufalls Zahl $r_i$ aus dem Interval $[0,1]$. Behalten Sie $x_i$ genau dann, wenn die Bedingung\n", "$$\n", "r_i \\cdot m(x_i) \\leq f(x_i)\n", "$$\n", "erfüllt ist. So erhalten Sie nur die Zufallszahlen, die gemäß f(x) verteilt sind, da Sie für jeden Wert von x nur den Bruchteil \n", "$$\n", "\\frac{m(x)}{f(x)}\n", "$$\n", "der x-Werte behalten. Stellen Sie das Ergebnis, gemäß der vorherigen Teilaufgabe, geeignet dar.\n", "> Tipp: Mithilfe von *plt.hist(\\[hist1, hist2\\], stacked=True,...)* können Sie sowohl die angenommenen als auch die abgelehnten Zufallszahlen in einem Histogramm darstellen." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# sample 2nd set of random numbers for accept-reject-method\n", "\n", "\n", "# accept-reject-method\n", "\n", "\n", "# plot functions and random numbers" ] } ], "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 }