{ "cells": [ { "cell_type": "markdown", "id": "3795466b", "metadata": {}, "source": [ "\n", "\n", "\n", "\n", "\n", "\n", "
Institut für Theoretische Physik
Universität zu Köln 		 Prof. Dr. Simon Trebst
Christoph Berke
\n", "
\n", "

Statistische Physik

\n", "

Übungsblatt 4

\n", "
\n", "

Wintersemester 23/24

\n", "\n", "\n", "**Website:** https://www.thp.uni-koeln.de/trebst/Lectures/2023-StatPhys.shtml\n", "\n", "**Abgabe**: 06.11.2023, 10:00 Uhr **Besprechung**: 07.11.2023 \n", "\n", "**Name**: Bitte geben Sie Ihren Namen an. " ] }, { "cell_type": "markdown", "id": "fa28e11f", "metadata": {}, "source": [ "\n", "\n", "
Aufgabe 14: Wiederholung und Zentraler Grenzwertsatz [7 Punkte]
" ] }, { "cell_type": "markdown", "id": "5aa8bfcf", "metadata": {}, "source": [ "In der Vorlesung haben Sie den **Zentralen Grenzwertsatz** als einen der mathematischen Kerne der statistischen Physik kennengelernt, der erklärt, warum wir für viele makroskopische Größen statistische, aber gleichzeitig *präzise* Aussagen machen können. In dieser Teilaufgabe sollen Sie die wesentlichen Aussagen des Zentralen Grenzwertsatzes numerisch demonstrieren." ] }, { "cell_type": "markdown", "id": "5efad1f0", "metadata": {}, "source": [ "### Vorbereitungen ###\n", "In den folgenden Zellen werden zunächst die benötigten Pakete geladen, sowie einige Voreinstellungen für die spätere Visualisierung vorgenommen. Stellen Sie sicher, dass Sie die Zellen ohne Fehlermeldung ausführen können und installieren Sie ggf. Pakete nach." ] }, { "cell_type": "code", "execution_count": null, "id": "7eeb23a7", "metadata": {}, "outputs": [], "source": [ "using CairoMakie, QuadGK, Optim, Statistics, Distributions, Colors\n", "import Base.rand" ] }, { "cell_type": "code", "execution_count": null, "id": "edce87ae", "metadata": {}, "outputs": [], "source": [ "# Define some Makie themes.\n", "ax_theme = Theme(Axis = (xticksvisible = false, yticksvisible = false, xticklabelsvisible = false, yticklabelsvisible = false, xautolimitmargin = (0,0)))\n", "ln_theme = Theme(Lines = (linewidth = 4, cycle = nothing, color = :red3))\n", "hi_theme = Theme(Hist = (cycle = nothing, color = :gray70))\n", "lb_theme = Theme(Label = (fontsize = 16, valign = :top, halign = :left, tellwidth = false, tellheight = false, padding = (10,10,10,10)))\n", "lg_theme = Theme(Legend = (valign = :top, halign = :right, tellwidth = false, tellheight = false, margin = (10,10,10,10)))\n", "fs_theme = Theme(fontsize = 24)\n", "customtheme1 = merge(ax_theme, ln_theme, fs_theme, hi_theme, lb_theme);" ] }, { "cell_type": "markdown", "id": "395a483b", "metadata": {}, "source": [ "### Zufallszahlen und Verteilungen ###\n", "In dieser Aufgabe sollen Sie die Gültigkeit des Zentralen Grenzwertsatzes für Zufallszahlen, die verschiedenen Verteilungen gehorchen, überprüfen. Hierfür generieren wir Zufallszahlen gemäß einiger Standardverteilungen, die vom Paket ```Distributions``` zur Verfügung gestellt werden. Weiterhin sollen Sie auch Zufallszahlen aus dem Intervall $\\left[0,1 \\right]$ die einer *beliebigen* Funktion $f(x)$ folgen generieren. In den beiden nächsten Zeilen werden zunächst die Funktionen definiert, mit denen eine, bzw. $n$ $f$-verteilte Zufallszahlen generiert werden können. Anschließend werden eine Testfunktion $f$ und mehrere Verteilungen aus dem Paket ```Distributions``` definiert." ] }, { "cell_type": "code", "execution_count": null, "id": "080dcea6", "metadata": {}, "outputs": [], "source": [ "# Funktion erzeugt eine f-verteilte Zufallszahl mittels Verwerfungsmethode.\n", "# f muss nicht normiert sein.\n", "function rand(f::Function)\n", "\n", " integral, error = quadgk(f, 0, 1)\n", " λ = 1.05 * (-minimum(optimize(x->-f(x)/integral, 0, 1)))\n", " \n", " x = rand()\n", " while rand() > f(x)/integral/λ\n", " x = rand()\n", " end \n", " return x\n", "end\n", "\n", "# Funktion erzeugt n f-verteilte Zufallszahlen mittels Verwerfungsmethode.\n", "function rand(f::Function, n::Integer)\n", " \n", " integral, error = quadgk(f, 0, 1)\n", " λ = 1.05 * (-minimum(optimize(x->-f(x)/integral, 0, 1))) \n", " x = zeros(n)\n", " for i in 1:n\n", " y = rand()\n", " while rand() > f(y)/integral/λ\n", " y = rand()\n", " end \n", " x[i] = y\n", " end \n", " \n", " return x\n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "2cf06736", "metadata": {}, "outputs": [], "source": [ "# beliebige Verteilung (es muss f(x)>0 in [0,1] gelten,\n", "# aber f muss nicht normiert auf [0,1] sein.\n", "f(x) = x^4 - x^2 + 1/2 # Testen Sie hier weitere Verteilungen.\n", "\n", "# Verteilungen aus dem Paket Distributions.\n", "unif = Uniform() # Uniforme Verteilung\n", "gauss = Normal() # Normalverteilung\n", "cauchy = Cauchy() # Cauchy-Verteilung\n", "arcsin = Arcsine(0,1) # Arcsine-Verteilung\n", "bindist = Binomial(100,0.7) # Binomialverteilung" ] }, { "cell_type": "markdown", "id": "a5bbb040", "metadata": {}, "source": [ "### Testen der Verteilungen ###\n", "Mit den oben definierten Funktionen können Sie jetzt Zufallszahlen nach beliebigen Funktionen oder nach den ausgewählten Standardverteilungen generieren. Genauer erzeugt ```rand(d, n)``` $n$ $d$-verteilte Zufallszahlen. Dabei kann $d$ entweder eine von Ihnen definierte Funktion sein, oder eine der Verteilungen aus ```Distributions```. Die nachfolgende Zelle erzeugt für alle oben definierten Möglichkeiten $n = 100{.}000$ Zufallszahlen und stellt diese gemeinsam mit der erwarteten Verteilung in einem Histogramm dar." ] }, { "cell_type": "code", "execution_count": null, "id": "325f306c", "metadata": {}, "outputs": [], "source": [ "n = 100_000 # Anzahl der Zufallszahlen\n", "fig = Figure()\n", "set_theme!(customtheme1)\n", "ga = fig[1,1] = GridLayout()\n", "\n", "# Uniform verteilte Zufallszahlen.\n", "ax11 = Axis(ga[1,1], ylabel = L\"p(x)\", limits = (nothing,nothing,0,nothing))\n", "hist!(ax11, rand(unif, n), normalization = :pdf, bins = 40)\n", "lines!(ax11, [Point2f(0,1), Point2f(1,1)])\n", "ylims!(ax11,0,1.2)\n", "Label(ga[1,1], \"uniform\")\n", "\n", "# Normalverteilte Zufallszahlen.\n", "ax12 = Axis(ga[1,2], limits = (nothing,nothing,0,nothing))\n", "hist!(ax12, rand(gauss, n), normalization = :pdf, bins = 40)\n", "lines!(ax12, -4:0.01:4, pdf.(gauss, -4:0.01:4))\n", "Label(ga[1,2], \"Gauß\")\n", "\n", "# Arcsin verteilte Zufallszahlen.\n", "ax13 = Axis(ga[1,3], limits = (nothing,nothing,0,nothing))\n", "hist!(ax13, rand(arcsin, n), normalization = :pdf, bins = 40)\n", "lines!(ax13, 0:0.01:1, pdf.(arcsin, 0:0.01:1))\n", "Label(ga[1,3], \"Arcsin\")\n", "\n", "# Binomialverteilte Zufallszahlen.\n", "ax21 = Axis(ga[2,1], ylabel = L\"p(x)\", xlabel = L\"x\", limits = (nothing,nothing,0,nothing))\n", "hist!(ax21, rand(bindist, n), normalization = :pdf, bins = 0:100)\n", "lines!(ax21, 0:100, pdf.(bindist, 0:100))\n", "Label(ga[2,1], \"Binomial\")\n", "\n", "# Cauchyverteilte Zufallszahlen.\n", "ax22 = Axis(ga[2,2], xlabel = L\"x\", limits = (nothing,nothing,0,nothing))\n", "hist!(ax22, rand(cauchy, n), normalization = :pdf, bins = -10:0.5:10)\n", "lines!(ax22, -10:0.02:10, pdf.(cauchy, -10:0.02:10))\n", "Label(ga[2,2], \"Cauchy\")\n", "\n", "# Zufallszahlen die beli\n", "ax23 = Axis(ga[2,3], xlabel = L\"x\", limits = (nothing,nothing,0,nothing))\n", "hist!(ax23, rand(f, n), normalization = :pdf, bins = 40)\n", "normierung, err = quadgk(f, 0, 1)\n", "lines!(ax23, 0..1, x->f(x)/normierung)\n", "Label(ga[2,3], \"benutzerdefiniert\")\n", "\n", "colgap!(ga, 5)\n", "rowgap!(ga, 5)\n", "\n", "fig" ] }, { "cell_type": "markdown", "id": "a83a9305", "metadata": {}, "source": [ "## Jetzt sind Sie dran! ##\n", "Der Zentrale Grenzwertsatz macht eine Aussage über den Mittelwert von $N$ $p$-verteilten, unabhängigen Zufallszahlen. Das besondere ist dabei die Universalität: Die konkrete Verteilung $p$ spielt für die Aussage (fast) keine Rolle. Sie sollen jetzt die Kernaussage des Satzes numerisch demonstrieren.\n", "\n", "## 1. Schärfe der Verteilung der Mittelwerte ##\n", "Wählen Sie eine Verteilung $p$ aus den obigen Beispielen aus. Implementieren Sie anschließend die folgenden Schritte:\n", "1. Erzeugen Sie $N$ $p$-verteilte Zufallszahlen $x_i$.\n", "2. Berechnen Sie den Mittelwert dieser $N$ Zufallszahlen, also $X_N = \\frac{1}{N} \\sum \\limits_{i=1}^N x_i$.\n", "3. Wiederholen Sie die Schritte 1. und 2. um mindestens $k = 10.000$ verschiedene Mittelwerte $X_N$ zu berechnen.\n", "4. Stellen Sie die Verteilung der Mittelwerte $X_N$ als Histogramm dar.\n", "5. Wiederholen Sie Schritt 1 bis 4 für verschiedene Werte von $N$. Stellen Sie die Histogramme für die verschiedenen $N$ gemeinsam dar.\n", "6. Berechnen Sie für jedes $N$ die Standardabweichung der $10.000$ Mittelwerte. Stellen Sie in einem zweiten Plot diese Standardabweichung als Funktion von $N$ dar.\n", "\n", "Führen Sie diese Schritte für (i) die Cauchy-Verteilung, (ii) eine von Ihnen beliebig gewählte Verteilung $f$ und (iii) eine weitere der obendefinierten Verteilungen aus dem Paket ```Distributions``` durch. Diskutieren Sie Ihre Ergebnisse für die Standardabweichung in einem kurzen Kommentar am Ende Ihres Codes.\n", "\n", "*Hinweis:* Varianz und Standardabweichung der Zahlen eines Arrays ```a``` können Sie mit den Funktionen ```var(a)``` und ```std(a)``` aus dem Paket ```Statistics``` berechnen." ] }, { "cell_type": "code", "execution_count": null, "id": "305a2d76", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "55a6205d", "metadata": {}, "source": [ "## 2. Konvergenz gegen eine Normalverteilung ##\n", "Im 1. Aufgabenteil sollten Sie beobachtet haben, dass die Verteilungen der Mittelwerte mit zunehmendem $N$ (fast) immer schärfer wird. Jetzt beschäftigen wir uns mit der konkreten Form der Verteilung. Dazu gehen wir analog zum 1. Aufgabenteil vor, berechnen also erneut viele Mittelwerte $X_N$ von jeweils $N$ $p$-verteilten Zufallszahlen $x_i$. Dieses Mal sollen sie nicht direkt die Mittelwerte als Histogramm darstellen, sondern die Verteilung der Größe\n", "\n", "\\begin{equation}\n", " \\sqrt{N} \\cdot \\frac{X_N - \\langle x \\rangle_p}{ \\sigma_p} \\,,\n", "\\end{equation}\n", "\n", "wobei $\\langle x \\rangle_p$ und $\\sigma_p$ Mittelwert und Standardabweichung der Verteilung $p$ sind, die den Zufallszahlen $x_i$ zu Grunde liegt. Stellen Sie die Histogramme für verschiedene Werte von $N$ dar, etwa $N \\in \\left[1,2,5,10 \\right]$. Gegen welche Verteilung konvergiert das Histogramm für große $N$? Zeichnen Sie diese Funktion ebenfalls in Ihre Abbildung ein.\n", "\n", "\n", "*Hinweis:* Die Werte für $\\langle x \\rangle_p$ und $\\sigma_p$ geben wir Ihnen in der folgenden Zelle für die Verteilungen aus dem Paket ```Distributions``` vor. Für Ihre eigene Verteilung $f$ können Sie $\\langle x \\rangle_f$ und $\\sigma_f$ nähern, indem Sie Mittelwert und Standardabweichung einer großen Anzahl an Samples ($\\approx 1{.}000.{000}$) berechnen. " ] }, { "cell_type": "code", "execution_count": null, "id": "098115ab", "metadata": {}, "outputs": [], "source": [ "mean_asin, std_asin = 0.5, sqrt(0.125)\n", "mean_cauchy, std_cauchy = cauchy.μ, cauchy.σ\n", "mean_gauss, std_gauss = gauss.μ, gauss.σ\n", "mean_binom, std_binom = bindist.p*bindist.n, sqrt(bindist.p*bindist.n*(1-bindist.p))\n", "mean_unif, std_unif = (maximum(unif) + minimum(unif)) / 2, (maximum(unif) - minimum(unif)) / sqrt(12)" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.9.3", "language": "julia", "name": "julia-1.9" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.9.3" } }, "nbformat": 4, "nbformat_minor": 5 }