{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## Jupyter Notebook Tutorial zur \n", "# Fehlerrechnung im Anfängerpraktikum\n", "\n", " Günter Quast, März 2020, \n", " aktualisiert Jan. 2022\n", "---\n", "\n", "## Jupyter Notebooks\n", "\n", "Diese Datei vom Typ `.ipynb` enthält ein Tutorial als `jupyter notebook`.\n", "`jupyter` bietet eine Browser-Schnittstelle mit einer (einfachen) Entwicklungsumgebung\n", "für `python`-Code und erklärende Texte im intuitiven *markdown*-Format. Die Eingabe von Formeln im `LaTeX`-Format wird ebenfalls unterstützt. \n", "\n", "Eine Zusammenstellung der wichtigsten Befehle zur Verwendung von *jupyter* als\n", "Arbeitsumgebung findet sich im Notebook\n", "[*JupyterCheatsheet.ipynb*](JupyterCheatsheet.ipynb). \n", "Die kurzen Abschnitte geben einen Überblick über das Notwendigste, um mit dem \n", "Tutorial arbeiten zu können. \n", "\n", "In *jupyter* werden Code und Text in jeweils einzelne Zellen eingegeben. \n", "Aktive Zellen werden durch einen blauen Balken am Rand angezeigt. Sie können sich in zwei Zuständen befinden: im Edit-Mode ist das Eingabefeld weiß, im Command-Mode ist es\n", "ausgegraut. Durch Klicken in den Randbereich wird der Command-Modus gewählt, ein Klick \n", "in das Textfeld einer Code-Zelle schaltet in den Edit-Mode. Die Taste `esc` kann\n", "ebenfalls verwendet werden, um den Edit-Modus zu verlassen. \n", "\n", "Eingabe von `a` im Command-Mode erzeugt eine neue leere Zelle oberhalb der aktiven Zelle,\n", "`b` eine unterhalb. Eingabe von `dd` löscht die betreffende Zelle.\n", "\n", "Zellen können entweder den Typ `Markdown` oder `Code` haben. Die Eingabe von `m` im Command-Modus setzt den Typ Markdown, Eingabe von `y` wählt den Typ Code.\n", "\n", "Prozessiert - also Text gesetzt oder Cocde ausgeführt - wird der Zelleninhalt durch \n", "Eingabe von `shift+return`, oder auch `alt+return` wenn zusätzliche eine neue, leere Zelle erzeugt werden soll. \n", "\n", "Die hier genannten Einstellungen sowie das Einfügen, Löschen oder Ausführen von Zellen können auch über das PullDown-Menü am oberen Rand ausgeführt werden. \n", "\n", "---\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "> Bei diesem Tutorial handelt es sich um einen Ansatz nach dem Prinzip\n", "\"_learning by doing_\".\n", "Es werden einfache Beispiele in der Sprache _Python_ verwendet, um konkrete Anwendungen\n", "zu verdeutlichen, die auch später in eigenen Arbeiten verwendet werden können. \n", "Dieses Tutorial baut auf anderen auf, die _Python_-Grundkenntnisse vermitteln \n", "(*PythonIntro.ipynb* oder *PythonCheatsheet.ipynb*), den Umgang mit den Paketen \n", "_numpy_ und _matplotlib_ zur Visualisierung von Daten (*matplotlibTuotrial.ipynb) \n", "und Grundlagen der Statistik (*IntroStatistik.ipynb*) vertiefen." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Allgemeine Einstellungen und nützliche Pakete" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import sys, os\n", "\n", "#\n", "import numpy as np, matplotlib.pyplot as plt\n", "\n", "# Zeilen mit % oder %% am Anfang sind sogenannte \"magische Kommandos\",\n", "# die den Zellentyp oder die Art der Anzeige von Grafiken angeben\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "---\n", "\n", "# (Moderne) Fehlerrechnung im Praktikum \n", "---\n", "\n", "# Inhalt \n", "\n", "## 1. [Einleitung](#Sec-Einleitung)\n", "\n", "### 1.1 [Bedeutung der Messunsicherheit](#SubSec-Messunsicherheit)\n", "\n", "### 1.2 [Aufnahme einer Messreihe zur Bestimmung der Messunsicherheit](#SubSec-Messreihe)\n", "\n", "### 1.3 [Häufigkeiten und Histogrammdarstellung](#SubSec-Häufigkeit)\n", "\n", "### 1.4 [Nomenklatur und Rechenregeln](#SubSec-Rechenregeln)\n", "\n", "\n", "## 2. [Fehlerfortpflanzung](#Sec-Fehlerfortplanzung)\n", "\n", "### 2.1 [Vorüberlegungen](#Sec-Vorüberlegungen)\n", "\n", "### 2.2 [Das Fehlerfortpflanzungsgesetz bei unabhängigen Unsicherheiten](#SubSec-Unabhängig)\n", "\n", "### 2.3 [Fehlerfortpflanzung mit Computer-Unterstützung](#SubSec-Computer)\n", "\n", "### 2.4 [Das Fehlerfortpflanzungsgesetz bei korrelierten Unsicherheiten](#SubSec-korreliert)\n", "\n", "\n", "## 3. [Anpassung von Modellen an Messdaten und Parameterschätzung](#Sec-Parameterschätzung)\n", "\n", "### 3.2 [Bestimmung der Parameterunsicherheiten](#Sec-Parameterunsicherheiten)\n", "\n", "### 3.3 [Beurteilung der Qualität einer Anpassung](#Sec-Fitqualität)\n", "\n", "### 3.4 [Verfahren und Programmpakete zur Parameteranpassung](#Sec-Programmpakete)\n", "\n", "### 3.5 [Details zum Anpassungspaket *kafe2*](#Sec-kafe)\n", "\n", "### 3.6 [Details zur Modellanpassung mit PhyPraKit/phyFit](#Sec-phyFit)\n", "\n", "\n", "## 4.0 [Spezielle Anwendungsfälle](#Sec-Sepical)\n", "\n", "### 4.1 [Transformation von Parametern und Konfidenzintervallen mittels Anpassung](#Sec-ParTrafo)\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "---\n", "# 1. Einleitung [^](#Inhalt)\n", "---\n", "\n", "Grundlage der Physik sind quantitative Beobachtungen von Naturphänomenen, ausgedrückt als\n", "Messgrößen mit einem Zahlenwert und einer Maßeinheit. Beobachtungen sind immer zufälligen\n", "Einflüssen unterworfen, die sich nicht vermeiden lassen und die das Messergebnis verändern.\n", "Dadurch bedingte Abweichungen vom angenommenen \"wahren Wert\" werden als Messunsicherheiten\n", "bezeichnet. Früher hat man dafür den unzutreffenden Begriff \"Messfehler\" verwendet, weil\n", "man der Ansicht war, durch genügend Sorgfalt und Aufwand beim Messprozess und den Einsatz\n", "bestmöglicher Technik die Messunsicherheit auf Null drücken zu können. Seit der Entwicklung\n", "der Quantenphysik als Grundlage aller Gebiete der Physik wissen wir aber, dass zufällige \n", "Einflüsse auf Beobachtungsgrößen prinzipiell unvermeidlich sind: in der Quantenpysik werden Erwartungswerte für Beobachtungsgrößen vorhergesagt, denen Verteilungsdichten mit nicht verschwindender Varianz zu Grunde liegen. \n", "Beispiele sind die mittlere Lebensdauer eines atomaren Zustands, eines radioaktiv\n", "zerfallenden Kerns oder Elementarteilchens, oder auch die prinzipielle Unbestimmtheit\n", "von Ort und Impuls eines Quantenteilchens. Schon in klassischen Systemen, die durch die\n", "gegenseitigen Wechselwirkungen vieler Teilchen bestimmt sind, ist es praktisch unmöglich,\n", "statistische Fluktuationen von Messergebnissen durch unkontrollierbare Einflüsse zu\n", "vermeiden. Bei immer weiter verbesserten Beobachtungsmethoden stößt man dann letztlich\n", "an die quantenmechanischen Grenzen. \n", "\n", "Die typische Problemstellung beim Messen besteht also darin, den Werte einer interressierenden Größe $y$ aus einer oder mehreren Messgrößen $x_i$ zu bestimmen, \n", "also $y(x_1, \\ldots, x_n)$.\n", "Beispiele sind etwa die Bestimmung eines ohmschen Widerstands aus Messungen von Strom und\n", "Spannung, einer Federkonstanten aus der Auslenkung und einer wirkenden Kraft, die Messung der Erdbeschleunigung durch die Auslenkung einer Feder mit bekannter Federkonstante\n", "durch ein angehängtes Gewicht mit bekannter Masse, oder die Messung der Elementarladung\n", "aus dem Radius einer Kreisbahn, auf der sich ein Elektronenstrahl in einem Magnetfeld bekannter Stärke bewegt, und viele mehr. \n", "\n", "Die Eingangsgrößen $x_i$ sind als Messgrößen mit unvermeidbaren, zufälligen Messunsicherheiten behaftet. Damit ist auch die abgeleitete Größe $y$ eine Größe\n", "mit einer zufälligen Unsicherheit, deren Wert es zu bestimmen gilt. \n", "\n", "Messgrößen sind ein Beispiel für *Zufallsgrößen*, also Größen, die abhängig vom Zufall verschiedene Werte annehmen können. Diese Werte sind nach einer bestimmten *Verteilungsdichte* (für kontinuierliche Zufallsgrößen) bzw. einer bestimmten *Wahrscheinlichkeit* $p(x)$ (für diskrete Zufallsgrößen) um den wahren Wert der Größe angeordnet. Als Maß für die Unsicherheit wird die Standardabweichung dieser Verteilungsdichte um den wahren Wert der Messgröße verwendet. Da dieser wahre Wert (nach einer endlichen Anzahl von Messungen) nur näherungsweise bekannt ist, verwendet man den Messwert als Schätzung für den wahren Wert und gibt die Standardabweichung als\n", "Unsicherheit an:\n", "\n", "$M = (m \\pm \\Delta m) \\cdot {\\rm }$\n", "\n", "Eine übliche grafiche Darstellung einer Messung \n", "\n", "$G=1\\pm0.3\\,{\\rm a.u.}$\n", "\n", "(dabei steht a.u. für \"arbitrary units\", d. h. beliebige Einheiten) erhalten Sie,\n", "wenn Sie die beiden folgenden Zellen mit *python*-Code ausführen \n", "(durch Anklicken und Eingabe von `+`)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# diesen Code durch Eingabe von ausführen\n", "\n", "# Datenpunkt 1.0 +/- 0.3\n", "g_m = [1.0]\n", "g_e = [0.3]\n", "\n", "# Erzeugen einer Grafik\n", "plt.figure(1, figsize=(5.0, 5.0))\n", "plt.errorbar([1], g_m, yerr=g_e, fmt=\"bo\")\n", "\n", "# Achsenbeschriftungen hinzufügen\n", "plt.ylabel(\"G (a.u.)\", size=\"xx-large\")\n", "plt.xlabel(\"Nr. der Messung\", size=\"x-large\")\n", "# Grenzen festlegen\n", "plt.xlim(0.5, 1.5)\n", "plt.ylim(0.0, 1.1 * (max(g_m) + max(g_e)))\n", "\n", "# Text hinzufügen\n", "plt.text(0.75, 0.4, \"G=%.2f +/- %.2f a.u.\" % (g_m[0], g_e[0]), size=\"x-large\")\n", "\n", "# keine Anzeige der Achsenbeschriftung für die x-Achse\n", "plt.gca().axes.get_xaxis().set_visible(False) # keine Anzeige der x-Achsenbeschriftung\n", "\n", "# Graifk anzeigen\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## 1.1 Bedeutung der Messunsicherheit \n", "---\n", "\n", "Die Angabe der Unsicherheit ist von entscheidender Wichtigkeit zur Bewertung einer Messung. \n", "Als Beispiel betrachten wir zwei Messungen der Lichtgeschwindigkeit, die zwar \n", "den gleichen Messwert ergeben, aber unterschiedliche Unsicherheiten haben.\n", "\n", "$$\\displaystyle\n", "\\begin{array}{ll}\n", "{\\rm Messung\\,1:\\,}\\, c = & (3.09 \\pm 0.15)\\cdot 10^8 \\frac{\\rm m}{\\rm s} \\\\\n", "{\\rm Messung\\,2:\\,}\\, c = & (3.09 \\pm 0.01)\\cdot 10^8 \\frac{\\rm m}{\\rm s} \\\\\n", "\\end{array}$$\n", "Führen Sie den Code in der Zeile unten aus [1].\n", "\n", "[1] **Anmerkung**: Stören Sie sich nicht an den Details im Code unten, die zur\n", "\"Verschönerung\" der grafischen Darstellung dienen. Das verwendete Paket \n", "*matplotlib.pyplot* bietet sehr viele Optionen zur Beeinflussung der Grafik,\n", "die man ggf. nachschlägt. Im Laufe der Zeit werden Sie sich die wichtigsten\n", "davon merken. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# diesen Code durch Eingabe von ausführen\n", "\n", "c_m = [3.09e8, 3.09e8]\n", "c_e = [0.15e8, 0.01e8]\n", "c_w = 2.99792458e8\n", "\n", "plt.figure(1, figsize=(6.0, 5.0))\n", "plt.errorbar([1, 2], c_m, yerr=c_e, fmt=\"bo\")\n", "\n", "plt.axhline(c_w, color=\"darkred\", linewidth=3)\n", "plt.text(0.55, 3.005e8, \"wahrer Wert\", color=\"darkred\")\n", "\n", "plt.ylabel(\"c (m/s)\", size=\"x-large\")\n", "plt.xlabel(\"Nr. der Messung\")\n", "plt.title(\"Messungen der Lichtgeschwindigkeit\")\n", "# noch einige Verschönerungen:\n", "plt.xlim(0.5, 2.5)\n", "ax = plt.gca()\n", "ax.locator_params(nbins=3)\n", "plt.grid()\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Man sieht, dass die Messunsicherheit der ersten Messung mit dem wahren Wert der Lichtgeschwindigkeit von 299 792 458 m s−1 überlappt, während der Fehlerbalken der zweiten Messung weit davon entfernt endet. Die zweite Messung\n", "ist damit klar im Widerspruch zum - übrigens im SI-Einheitensystem als exakt\n", "definierten - Literaturwert. \n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "### 1.2 Aufnahme einer Messreihe zur Bestimmung der Messunsicherheit \n", "---\n", "\n", "Zufällige Einflüsse auf den Wert einer Messung lassen sich durch wiederholte Ausführung des Messprozesses bestimmen. Dazu wird ein Messreihe aus mehreren Messungen \n", "$x_i, {\\small i=1, \\ldots, N}$ \n", "durchgeführt. Die einzelnen Messwerte $m_i$ unterscheiden sich, weil der allen \n", "zu Grunde liegende \"wahre Werte\" $w$ von jeweils zufälligen Abweichungen $z_i$ betroffen ist:\n", "\n", "$m_i \\,=\\, w + z_i \\,.$\n", "\n", "Unter recht schwachen Voraussetzung ist die zu Grunde liegende Verteilung in der Praxis häufig eine Normal- (oder auch Gauß-)Verteilung. Dies gilt insbesondere dann,\n", "wenn sehr viele Störungen (die sogar einer beliebigen Verteilungsdichte folgen dürfen)\n", "zufällig zusammen wirken. Das ist die Aussage des \n", "**zentralen Grenzwertsatzes der Wahrscheinlichkeitstheorie**. \n", "\n", "Die **statistische Messunsicherheit** einer Einzelmessung lässt sich als \n", "Standardabweichung der Werte der Messreihe darstellen: \n", "\n", "$\\sigma_{m} = \\sqrt{ \\frac{1}{N-1} \\displaystyle \\sum_{i=1}^N (m_i-\\bar{m})^2} \\,\\,{\\rm mit}\\, \\bar{m} = \\frac{1}{N} \\displaystyle \\sum_{i=1}^N m_i \\,.$\n", "\n", "Dabei ist $\\bar{m}$ der (arithmetische) Mittelwert der aufgenommenen Messwerte.\n", "Das Quadrat der Standardabweichung bezeichnet man als **\"Varianz\"** $V_m = {\\sigma_m}^2$.\n", "In Worten ausgedrückt: \n", "> Die Varianz ist die Summe der quadrierten Abweichungen der Einzelmessungen vom Erwartungswert. \n", " \n", "$\\bar{m}$ und $\\sigma_m$ nach obiger Formel sind die besten unverzerrten Schätzungen der Parameters $w$ und $\\sigma$ der angenommenen Gauß-Verteilung \n", "\n", "$g(m) = \\frac{1}{\\sqrt{2\\pi}{\\sigma}} \\exp {\\frac{(m - w)^2}{2{\\sigma}^2} } \\,.$\n", "\n", "Vielleicht ist Ihnen aufgefallen, dass die mittlere quadratische Abweichung durch Division durch $(N-1)$ statt $N$ berechnet wird. Dies ist die sogenannte \"Bessel-Korrektur\", die dafür sorgt, dass die Schätzung aus den gemessenen Daten unverzerrt ist. Eine Schätzung wird im allgemeinen dann als unverzerrt bezeichnet,\n", "wenn der Erwartungswert des Schätzwerts (siehe [Abschnitt 1.4](#SubSec-Rechenregeln)) dem wahren Wert entspricht. Eine Konsequenz ist, dass man - sinnvollerweise - zur Berechnung der Standardabweichung mindestens zwei Messwerte braucht.\n", "\n", "Abweichungen vom wahren Wert, die durch \"echte\" Fehler im Messprozess oder die Eichung der verwendeten Messgeräte verursacht werden, können durch wiederholte Messungen nicht bestimmt werden. Erst durch Variation des Messverfahrens, die Verwendung jeweils verschiedener Geräte oder durch Vergleich mit genau bekannten Referenzgrößen in einer \"Kalibrationsmessung\" werden solche Unsicherheiten erkennbar. Man nennt sie **systematische Unsicherheiten**. Typisch für systematische Unsicherheiten ist, dass sie zu Abweichungen führen, die allen Messungen gemeinsam, also korreliert sind. Beispielsweise ist bei Widerstandsthermometern die Temperatur nicht perfekt proportional zum elektrischen Widerstand. Alle mit demselben Thermometer durchgeführten Messungen zeigen dann ähnliche Abweichungen vom wahren Wert. Oft werden die zu erwartenden systematischen Unsicherheiten vom Hersteller im Datenblatt eines Messgeräts angegeben. \n", "\n", "Wenn man eine Messreihe mit $N$ Einzelmessungen ausgeführt hat, sollte man als Messwert nicht eine beliebige Einzelmessung, sondern den Mittelwert $\\bar{m}$ aller Messungen der Messreiche zitieren. Der Grund dafür ist, dass die Genauigkeit des Mittelwertes besser ist als die einer Einzelmessung, da sich statistisch bedingte Unsicherheiten teilweise herausmitteln. Somit ergibt sich ein besserer Schätzwert für den wahren Wert. Weiter unten wird gezeigt, dass die **Unsicherheit des Mittelwerts** gegeben ist durch\n", "\n", "$\\sigma_{\\bar{m}} = \\frac{\\sigma_m} {\\sqrt{N}}\\,,$\n", "\n", "also um einen Faktor $1/\\sqrt{N}$ geringer ist als die Unsicherheit eines einzelnen Messwerts.\n", "\n", "Das *python*-Paket *numpy* bietet einfache Möglichkeiten, Mittelwert und Standardabweichung einer Messreihe zu berechnen. Dazu müssen die Daten als *numpy array* vorliegen. Betrachten Sie folgendes Code-Fragment:\n", "\n", "```\n", "# Messdaten\n", "d=[0.846, 1.513, 0.388, 0.782, 0.403, 0.743, 0.976, 0.988, 1.736, 1.128]\n", "n= len(d)\n", "\n", "data=np.array(d)\n", "mean=data.mean()\n", "std=data.std(ddof=1)\n", "\n", "print(\"Mittelwert der Messreihe: \", mean)\n", "print(\"Standardabweichung der Messungen: \", std)\n", "print(\"Unsicherheit des Mittelwerts: \", std/np.sqrt(n) )\n", "\n", "```\n", "\n", "[2] **Anmerkung**: der Parameter \"*ddof=1*\" (*ddof*= \"delta degrees of freedom\", d. h. Differenz in der Zahl der Freiheitsgrade) sorgt dafür, dass die Standardabweichung mit der \"Bessel-Korrektur\" berechnet wird (die Zahl der Freiheitsgrade wird durch die Verwendung des Mittelwerts um eins reduziert)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Übung:\n", "Geben Sie den Code oben in die leere Zelle unten ein und führen Sie ihn aus." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# eigenen Code hier eingeben:\n", "# -->" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "### 1.3 Häufigkeiten und Histogrammdarstellung \n", "---\n", "\n", "\n", "Wenn - häufig automatisiert - große Messreihen aufgenommen wurden, empfiehlt es sich, \n", "sie grafisch darzustellen. Dazu werden die Messdaten $m_1, \\ldots, m_N$ in $k$ Klassen\n", "eingeteilt, z. B. Intervalle $1, \\ldots, k$ mit Intervallgrenzen $M_0, \\ldots, M_k$,\n", "und die Häufigkeit des Vorkommens von Werten in jedem der Intervalle aufgetragen.\n", "Die übliche Darstellung ist ein Balkendiagramm, bei dem die Höhe eines Balkens im\n", "Intervall $j$ die Häufigkeit $H_j$ einer Beobachtung in diesem Intervall angibt. \n", "\n", "Für $N$ Messungen gilt:\n", "\n", "Summe der Häufigkeiten:\n", " $\\displaystyle \\sum_{j=1}^k H_j = N\\,,$\n", "\n", "relative Häufigkeiten: \n", " $h_j=\\frac{H_j}{N} \\displaystyle \\,\\Leftrightarrow\\,\\sum_{j=1}^k h_j = 1\\,,$\n", " \n", "Mittelwert der Messgrößen:\n", " $\\overline{m}= \\displaystyle \\sum_{j=1}^k h_j \\, \\overline{M}_j \\,,$\n", "\n", "und Messunsicherheit: \n", " $\\sigma_x= \\displaystyle \\sum_{j=1}^k h_j (\\overline{M}_j - \\overline{m})^2\\,.$ \n", "$\\overline{M}_j$ ist dabei die Mitte des jeweiligen Intervalls, \n", "$\\overline{M}_j=(M_j - M_{j-1})/2 \\,.$ \n", "\n", "Folgendes Code-Beispiel illustriert eine Histogrammdarstellung:\n", "```\n", "# Daten für Histogramm\n", "Ri = np.array([139.9,139.9,140.1,140.0,139.9,140.0,140.0,139.8,140.0,140.3,140.0,140.0,140.0,140.2,140.], dtype=float)\n", "\n", "bconts, bedges, _p = plt.hist(Ri, bins=np.linspace(139.8, 140.3, 6))\n", "plt.ylabel('Häufigkeit')\n", "plt.xlabel('Widerstand R (k$\\Omega$)')\n", "plt.show()\n", "```\n", "Zur statistischen Auswertung kann das Paket *PhyPraKit* verwendet werden, in dem die Berechnung der obigen Formeln implementiert ist:\n", "```\n", "# Berechnung der statistischen Daten aus dem Histogramm\n", "from PhyPraKit import histstat\n", "mean, sigma, sigma_mean = histstat(bconts, bedges, pr= True)\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Übung: \n", "Führen Sie diesen Beispiel-Code in der Zelle unten aus." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# eigenen Code hier eingeben\n", "# -->" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 1.3.1 Berücksichtigung von systematischen Unsicherheiten \n", "Beim obigen Code-Beispiel handelte es sich offensichtlich um die wiederholte Bestimmung eines ohmschen Widerstands, z. B. mit einem Multimeter. Nehmen Sie nun an, dass die Zuleitungen einen Widerstand von $0.3\\,\\Omega$ haben, der auf $\\pm 0.1\\,\\Omega$ genau bekannt ist.\n", "\n", "Es ist klar, dass der zusätzliche Widerstand einfach korrigierbar ist - der tatsächliche\n", "Widerstand des Bauteils ist um $0.3\\,\\Omega$ kleiner als der aus dem Histogramm bestimmte Mittelwert.\n", "\n", "Im Folgenden wird es darum gehen, wie sich die Unsicherheit auf den ohmschen Widerstand der Zuleitungen auf die Unsicherheit des Bauteilwiderstands auswirkt. Wir werden sehen, dass sich die Varianz, also das Quadrat der Unsicherheit, als Summe der statistischen Varianz und dem Quadrat der systematischen Unsicherheit ergibt, \n", "$\\sigma_{\\rm tot}^2 = \\sigma_{\\rm stat}^2 + \\sigma_{\\rm syst}^2 \\,.$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "### 1.4 Nomenklatur und Rechenregeln \n", "---\n", "\n", "Der (arithmethische) Mittelwert $\\overline{m}$ von Messgrößen wurde bereits in [Abschnitt 1.2](#SubSec-Messreihe) eingeführt. Man kann allgemein für Zufallsgrößen $x$ mit unterschiedlichen Verteilungsdichten zeigen (u. a. für gaußverteilte Zufallsgrößen), dass der Erwartungswert der Verteilungsdichte dem Mittelwert entspricht. Daher gibt es verschiedene Nomenklaturen für den Erwartungswert einer Zufallsgröße $x$, die einer diskreten Wahrscheinlichkeit $p(x)\\equiv h_i$ oder einer kontinuierlichen Wahrscheinlichkeitsdichte $f_x(x)$ folgt. Bisher haben wir für den Mittelwert die Bezeichung $\\bar{x}$ verwendet. Für den Erwartungswert üblich sind die äquivalenten Darstellungsweisen \n", "\n", "$\\bar{x} \\equiv \\rm{E}[x] \\equiv \\left$\n", " \n", "Je nach Art der Daten bzw. deren Verteilung gelten im Einzelfall zur Berechnung \n", "des Erwartungswerts $\\left< x \\right>$ die folgenden Formeln: \n", "\n", "$$\\begin{array}{c|ccc}\n", "{\\rm Erwartungswert} & {\\rm Stichprobe} \\, x_i, {\\small i=1,...,N} &\n", " {\\rm diskrete \\,Verteilung} \\, h_i & {\\rm Verteilungsdichte} \\,f_x(x) \\\\\n", "\\hline % ---------------------------------------------------------------------------- \n", "\\left< x \\right> \\,= & \\frac{1}{N}\\sum_{i=1}^N x_i &\n", " \\sum_{i=1}^N x_i h_i & \\int_{-\\infty}^\\infty xf_x(x)\\,dx \\\\\n", "\\end{array}$$\n", "\n", "Für Erwartungswerte gelten allgemein einfache Rechenregeln, die aus der Linearität der Definition des Erwartungswerts folgen. Im folgenden sind $a$, $b$ Konstanten und $x$, $x_1$ und $x_2$ Zufallsgrößen. Es gilt: \n", "\n", "$\\left< a \\right> = a$, d. h. auch $\\left< \\left< x \\right>\\right> = \\left< x \\right>$\n", "\n", "$\\left< ax \\right> = a \\left< x \\right>$\n", " \n", "$\\left< x_1 + x_2 \\right> = \\left< x_1 \\right> + \\left< x_2 \\right>$\n", "\n", "und damit auch $\\left< a x_1 + b x_2 \\right> = a \\left< x_1 \\right> + b \\left< x_2 \\right>$ \n", " \n", "Der Erwartungswert der quadratischen Abweichung vom Mittelwert, die Varianz $\\mathrm{V}_x$, \n", "lässt allgemein wie folgt schreiben:\n", "\n", "$\\rm{V}_x = \\left< { \\left(x - \\left< x\\right> \\right)}^2 \\right>$ \n", "\n", "Durch Ausmultiplizieren erhält man die äquivalente Darstellung\n", "\n", "$$\\begin{array}{lcl} \n", "\\rm{V}_x & = & \\left< x^2 - 2 x \\left< x \\right> + \\left< x \\right>^2 \\right> \\\\\n", " ~ & = & \\left< x^2 \\right> - \\left<2x \\left< x\\right>\\right> + \\left<\\left< x \\right>^2\\right> \\\\ \n", " ~ & = & \\left< x^2 \\right> - 2 \\left< x \\right>^2 + \\left< x \\right>^2 \\\\\n", " ~ & = & \\left< x^2 \\right> - \\left< x \\right>^2 \\\\\n", "\\end{array}$$\n", "\n", "Beachten Sie: der Erwartungswert eines Produktes ist im allgemeinen nicht gleich dem Produkt der Erwartungswerte,\n", "\n", "> $\\left< x_1 x_2 \\right> \\ne \\left< x_1 \\right> \\left< x_2 \\right>$.\n", "\n", "Gleichheit ($\\left< x_1 x_2 \\right> = \\left< x_1 \\right> \\left< x_2 \\right>$) gilt allerdings, wenn $x_1$ und $x_2$ unabhängige Zufallszahlen sind. \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## 2. Fehlerfortpflanzung [^](#Inhalt)\n", "---\n", "\n", "Wir widmen uns nun der Frage, wie sich die Unsicherheit einer Funktion $y$ von $N$ Zufallsgrößen $y = y(x_1, \\ldots, x_N)$ ergibt. \n", "\n", "### 2.1 Vorüberlegungen \n", "\n", "Beginnen wir mit einem ganz einfachen Beispiel und berechnen die Varianz einer Größe $s$, die sich als Linearkombination zweier Zufallsgrößen $x_1$ und $x_1$ ergibt:\n", "$s = a\\,x_1 + b\\,x_2$. \n", "\n", ">$$\\displaystyle\\begin{eqnarray}\n", " \\rm{V}_s &=& \\left< { \\left(s - \\left< s \\right> \\right)}^2 \\right>\n", " = \\left< s^2 \\right> - \\left< s \\right>^2 \\\\\n", " &=& \\left< a^2 x_1^2 + b^2 x_2^2 + 2 a b x_1 x_2 \\right> \n", " - a^2 \\left< x_1 \\right>^2 - b^2 \\left< x_2 \\right>^2 \n", " - 2 a b \\left< x_1 \\right> \\left< x_2 \\right> \\\\\n", " &=& a^2 \\left( \\left - \\left^2 \\right)\n", " + b^2 \\left( \\left - \\left^2 \\right) \n", " + 2ab \\,\\bigl( \\left - \\left \\left \\bigr)\\\\\n", " &=& a^2 \\rm{V}_1 + b^2 \\rm{V}_2 \n", " + 2ab \\,\\bigl(\\left< x_1 x_2 \\right> - \\left< x_1 \\right>\\left< x_2 \\right>\\bigr)\n", "\\end{eqnarray}$$\n", "\n", "Der letzte Term verschwindet für unabhängige Zufallsvariable $x_1$ und $x_2$. \n", "Er enhält die sogenannte *Kovarianz*, \n", "\n", "$\\mathrm{cov}(x_1, x_2) = \\left< (x_1 - \\left)(x_2 - \\left)\\right>\n", " = \\left< x_1 x_2 \\right> - \\left< x_1 \\right>\\left< x_2 \\right>$\n", "\n", "Die Kovarianz ist eine Verallgemeinerung der Varianz für zwei Variable: Es gilt $\\rm{cov}(x,x) = \\rm{V}_x$, d. h. die Kovarianz einer Variable mit sich selbst entspricht der Varianz.\n", "\n", "Wir halten also insgesamt als Ergebnis für die Varianz einer Summe von Zufallsvariablen fest:\n", "\n", "$ \\rm{V}_s = a^2 \\rm{V}_1 + b^2 \\rm{V}_2 + 2ab \\, \\rm{cov}(x_1, x_2) \\,\\,\\, $\n", " für $s = a\\,x_1 + b\\,x_2\\,$\n", " Gleichung (2.1) .\n", "\n", "\n", "Für eine allgemeine Funktion $f(x)$ einer Zufallsgröße $x$ lässt sich die Frage\n", "nach der Unsicherheit bzw. Varianz $\\mathrm{V}_f$ näherungsweise beantworten, wenn man\n", "die Taylor-Entwicklung von $f$ um den Erwartungswert $\\left$ betrachtet:\n", "\n", "$f(x) = f(\\left) + \\left. {\\frac{\\mathrm{d}f}{\\mathrm{d}x}}\\right| _{x=\\left} \n", " \\left(x-\\left\\right) + \\ldots$\n", "\n", "Da die Unsicherheit typischerweise eine kleine Größe ist und wir uns nur um eine Näherung der Funktion $f$ über den Bereich der Unsicherheit $\\sigma_x=\\sqrt{\\mathrm{V}_x}$ um den Erwartungswert $\\left$ interessieren, ist die Näherung in den meisten Fällen hinreichend genau.\n", "\n", "Unter Verwendung des Ergebnisses von oben erhalten wir:\n", "\n", "$\\mathrm{V}_f = \\left< f(x)^2 \\right> - \\left< f(x) \\right>^2 \n", " \\simeq \\left( \\left. \\frac{\\mathrm{d}f}{\\mathrm{d}x}\\right|_{x=\\left} \\right)^2 \\left< (x-\\left)^2 \\right> $ \n", "\n", "Der letzte Term ist die Varianz von $x$; damit erhalten wir also:\n", "> $ \\mathrm{V}_f \\simeq\\, \\left( \\left. \\frac{\\mathrm{d}f}{\\mathrm{d}x}\\right|_{x=\\left} \\right)^2 \\mathrm{V}_x \\,$;\n", "Gleichheit gilt bei linearen Funktionen $f(x)$. \n", "\n", "Die Varianz der Funktion $f(x)$ entspricht also der Varianz von $x$, multipliziert mit dem Wert der ersten Ableitung der Funktion am Erwartungswert $\\left$ zum Quadrat. \n", "\n", "Bei Funktionen eines Vektors von $n$ unabhängigen Zufallsvariablen $f(x_1, \\ldots, x_n)\\equiv f(\\vec{x})$ ergibt sich die Varianz $V_f$ als Summe von Termen analog zur obigen Gleichung, diesmal mit partiellen Ableitungen:\n", "> $\\mathrm{V}_f \\simeq \\sum\\limits_{i=1}^n \\left( \\left.\\frac{\\partial V}{\\partial x_i}\\right|_{\\vec{x}=\\left<\\vec{x}\\right>}\\right)^2 \\mathrm{V}_{x_i}$ Gleichung (2.2)\n", "\n", "Der in [Gleichung (2.1)](#Equ-2.1) hergeleitete Term, der die Kovarianz beinhaltet, musste für diese Herleitung nicht berücksichtigt werden, da die Zufallszahlen als unabhängig angenommen wurden. Dieser Term wird in [Abschnitt 2.4](#SubSec-korreliert) im Zusammenhang mit korrelierten Unsicherheiten aber wieder aufgegriffen." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "### 2.2 Das Fehlerfortpflanzungsgesetz bei unabhängigen Unsicherheiten \n", "--- \n", "\n", "Zusammenfassend kann man nach diesen Vorüberlegungen den allgemeinen Fall \n", "einer Funktion $y$ von $n$ unabhängigen Zufallsvariablen $y(x_1, \\ldots, x_n)$ \n", "wie folgt schreiben:\n", "\n", "\n", "> $\\sigma^2_{y}= \\displaystyle \\sum_{i=1}^n\n", " \\left( \\frac{\\partial y}{\\partial x_i} \\right)^2 \\sigma^2_{x_i}\\,$ \n", " Gleichung (2.3) ,\n", "> wobei die partiellen Ableitungen wieder am Erwartungswert ausgewertet werden.\n", "\n", "Wichtige Spezialfälle sind Summen, Differenzen, Produkte oder Quotienten von zwei Größen.\n", "\n", "Unter Anwendung des allgemeinen Ergebnisses erhält man:\n", "\n", "$y=x_1 \\pm x_2$\n", "> $\\Rightarrow \\color{blue}{\n", " {\\sigma_y}^2 = {\\sigma_1}^2 + {\\sigma_2}^2 } $\n", "\n", "Die quadrierte Unsicherheit der Summe oder Differenz\n", "zweier Messungen ergibt sich als die quadratische Summe ihrer Unsicherheiten. \n", "\n", "$y=x_1 \\cdot x_2 {\\rm ~oder~} y = \\frac{x_1}{x_2}$\n", "> $\\Rightarrow \\color{blue}{\n", " \\left( \\frac{\\sigma_y}{y} \\right)^2 \\simeq\n", " \\left(\\frac{\\sigma_1}{x_1} \\right)^2 \n", " + \\left(\\frac{\\sigma_2}{x_2} \\right)^2 }$\n", " \n", " Die quadrierte *relative* Unsicherheit (d. h. die Unsicherheit dividiert durch den Wert selbst) \n", " des Produkts oder Verhältnisses zweier Messungen ist die quadratische Summe ihrer relativen Unsicherheiten. \n", " \n", "Interessant und von besonderer Relevanz ist auch der Fall, wenn eine Funktion von einer\n", "Potenz einer Zufallsvariablen abhängt, also \n", " \n", "$y = x^m \\quad\\rightarrow \\frac{\\partial y}{\\partial x} = m x^{m-1}$\n", "\n", "$\\sigma_y^2 = \\left(\\frac{\\partial y}{\\partial x}\\right)^2 \\sigma_x^2 \n", " = m^2 x^{2m-2} \\sigma_x^2 = m^2 \\frac{y^2}{x^2} \\sigma_x^2 \n", " \\quad\\rightarrow \\frac{\\sigma_y^2}{y^2} = m^2 \\frac{\\sigma_x^2}{x^2}$\n", "> $\\Rightarrow \\color{blue}{\n", " \\frac{\\sigma_y}{x} \\simeq |m| \\frac{\\sigma_x}{x}} $\n", " \n", " \n", "Mit Hilfe des Fehlerfortpflanzungsgesetzes kann man auch die weiter oben schon besprochene Reduktion der Unsicherheit durch Mittelwertbildung der Messwerte $x_i$ von $N$ gleichartigen Messungen begründen: \n", "\n", "$y = \\bar{x} = \\frac{1}{N} \\displaystyle \\sum_{i=1}^N x_i $ \n", "$\\Rightarrow {\\sigma_y}^2 = \\frac{1}{N^2} \\displaystyle \\sum_{i=1}^N {\\sigma_i}^2 \n", " = \\frac{1}{N^2} N {\\sigma_x}^2 = \\frac{1}{N} {\\sigma_x}^2 $, also\n", " \n", "> $ \\color{blue}{ {\\sigma_\\bar{y}}^2 = \\frac{1}{N} {\\sigma_x}^2 }$\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Kleine Übung zu 2.2\n", "\n", "Bei der Untersuchung eines Sensors oder der Bewertung eines medizinischen Tests\n", "wird bestimmt, wie häufig er er auf ein Signal anspricht. \n", "Es wird $p$ mal ein positiver und $n$ mal eine negative Ausgang beobachtet. \n", "Die \"Effizienz\" zur Bewertung der Leistungsfähigkeit ist definiert als \n", "$\\epsilon = \\frac{p}{p+n}\\,$.\n", "\n", "\n", "Die statistischen Unsicherheiten der beobachteten Häufigkeiten sind\n", "$\\sigma_p =\\sqrt{p}$ und $\\sigma_n=\\sqrt{n}$. \n", "Wie groß ist die Unsicherheit $\\sigma_\\epsilon$ der Effizienz $\\epsilon$?\n", "\n", "> Lösung: $\\sigma_\\epsilon = \\sqrt{ \\frac{\\epsilon \\, (1-\\epsilon)}{p+n} }$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Hier eigene Berechnung eingeben:** \n", "--> \n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "### 2.3 Fehlerfortpflanzung mit Computer-Unterstützung \n", "---\n", "\n", "\n", "Für *python* gibt es das Paket *uncertainties*, das mit Unsicherheiten von Variablen\n", "umgehen kann.\n", "\n", "Falls nicht schon vorhanden, installieren Sie es mit `pip3 install uncertainties`. \n", "Die vollständige Dokumentation finden Sie unter dem Link\n", "[https://pythonhosted.org/uncertainties](https://pythonhosted.org/uncertainties) .\n", "\n", "Ein Beispiel illustriert die Anwendung:\n", "```\n", "from uncertainties import ufloat, umath\n", "x = ufloat(1.0, 0.1)\n", "y = ufloat(2.0, 0.1)\n", "r = umath.sqrt(x**2 + y**2)\n", "print( r ) \n", "```\n", "\n", "Lassen Sie den Code in der Zelle unten laufen und zeigen Sie, dass das\n", "Ergebnis korrekt ist ! \n", "Bemerkenswert ist noch, dass die *print()*-Funktion die korrekt gerundete\n", "Ausgabe erzeugt, also die Genauigkeit des ausgegebenen Wertes der der auf\n", "zwei Stellen gerundeten Unsicherheit entspricht." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Code hier eingeben:\n", "# -->" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Kurze Begründung, warum Sie dem Ergebnis des Programms trauen** \n", "->\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "---\n", "### 2.4 Das Fehlerfortpflanzungsgesetz bei korrelierten Unsicherheiten \n", "--- \n", "\n", "\n", "Bisher haben wir unabhängige Messungen betrachtet und konnten daher den Kovarianz-Term \n", "in der obigen Formel für die Summe von zwei Zufallsvariablen zu Null setzen. Dies\n", "ändert sich, wenn die Messgrößen gemeinsame Unsicherheiten haben, z. B. eine oben schon\n", "besprochene zusätzliche, allen Messungen gemeinsame systematische Unsicherheit. \n", "\n", "Eine gemeinsame Fehlerkomponente führt dazu, dass Abhängigkeiten zwischen den\n", "beobachteten Werten von Größen entstehen: ist die eine zu groß, so steigt auch\n", "die Wahrscheinlichkeit, dass eine zweite Variable zu groß gemessen wird - man\n", "nennt dies eine \"positive Korrelation\". Eine negative Korrelation bedeutet, dass\n", "für die zweite Variable mit höherer Wahrscheinlichkeit ein zu kleiner Wert gemessen\n", "wird, wenn der Wert der ersten zu groß war. Ein Beispiel dafür ist die Anpassung einer\n", "Ausgleichsgeraden $y=ax+b$ an Messwerte. Wird die Steigung $a$ der Gerade größer, so wird ihr\n", "$y$-Achsenabschnitt $b$ kleiner (und umgekehrt), die beiden Größen sind negativ korreliert. \n", "Bei zwei unabhängigen Variablen kann man aus dem Wert der ersten keine Information über die zweite Variable gewinnen (und umgekehrt), die Korrelation ist exakt Null. \n", "\n", "Zur Behandlung systematischer Unsicherheiten ist ein Verständnis der Kovarianzmatrix\n", "unerlässlich; daher beschäftigen wir uns etwas eingehender damit.\n", "\n", "#### 2.4.1 Kovarianz, Kovarianz- und Korrelationsmatrix \n", "---\n", "\n", "\n", "Als Maß für den Zusammenhang zwischen zwei Zufallsgrößen $x_1$ und $x_2$ dient die\n", "Kovarianz. Sie ist definiert als das Produkt der Abweichungen vom Erwartungswert\n", "der ersten und zweiten Zufallsgröße:\n", "\n", "$\\mathrm{cov}(x_1, x_2) = \\left< (x_1 - \\left< x_1\\right>) \\, (x_2 - \\left< x_2\\right>) \\right>.$\n", "\n", "Die Kovarianz ist entsprechend zur Varianz $\\mathrm{V}$ einer einzelnen Variablen konstruiert. \n", "Es gilt $\\mathrm{cov}(x_1, x_1) = \\mathrm{V}_{x_1}$. Für unabhängige Variablen verschwindet die\n", "Kovarianz; sie kann maximal den Wert $\\sqrt{\\mathrm{V}_1 \\cdot \\mathrm{V}_2} = \\sigma_1 \\cdot \\sigma_2$ annehmen.\n", "Wenn es mehrere Variable gibt, werden paarweise Werte für die Kovarianz berechnet;\n", "dann erhält man die sogenannte Kovarianzmatrix \n", "\n", "${\\bf V} = (\\mathrm{V}_{ij}) = \\left( \\mathrm{cov}(x_i, x_j) \\right).$\n", "\n", "Übersichtlicher ist eine Vektor-Schreibweise mit dem Variablen-Vektor $\\vec{x}$:\n", "\n", "${\\bf V} = \\bigl< \\bigl( \\vec{x} - \\left< \\vec{x} \\right> \\bigr) \\,\n", " \\bigl(\\vec{x} - \\left< \\vec{x} \\right> \\bigr)^T \\bigr>\n", " = \\left< \\vec{x}\\vec{x}^T\\right> \\, - \\, \\left< \\vec{x}\\right> \\left<\\vec{x}\\right>^T $\n", "\n", "Als quadratische Größe ist die Kovarianzmatrix wegen des typischerweise großen Wertebereichs\n", "der Matrixelemente nur schwer zu überblicken. \n", "Man verwendet daher eine geeignete Normierung und erhält die Matrix der Korrelationen, \n", "deren Elemente $\\rho_{ij}$ Werte zwischen $-1$ und $1$ annehmen:\n", "\n", "$\\bigl( \\rho_{ij} \\bigr) \\,=\\, \\Bigl (\\frac{ \\mathrm{V}_{ij} } {\\sigma_i\\,\\sigma_j} \\Bigr)$\n", "\n", "Zum besseren Verständnis mag ein Beispiel helfen. Wir erzeugen dazu Zufallsgrößen mit\n", "unabhängigen und gemeinsamen Fehlern:\n", "\n", "#### 2.4.2 Beispiel und grafische Darstellung von Kovarianz und Korrelation\n", "---\n", "\n", "\n", "Zwei Messungen $m_i$ werden modelliert als Summe aus dem wahren Werte $w_i$, einem *unabhängigen* zufälligen Anteil $z_i$ und einem *gemeinsamen* zufälligen Anteil $z_g$:\n", "\n", "> $m_1 = w_1 + z_1 + z_g$ \n", " $m_2 = w_2 + z_2 + z_g$ \n", "\n", "Zur Veranschaulichung wählt man ein sogenanntes \"Streudiagramm\" (engl.: \"scatter-plot\"), \n", "wie es in *matplotlib.pyplot.scatter()* implementiert ist. Hier einige Zeilen Code dazu:\n", "```\n", "import numpy as np, matplotlib.pyplot as plt\n", "\n", "def scatterplot(x, y):\n", " plt.figure(figsize=(7.5,7.5)) \n", " plt.scatter(x, y)\n", " plt.show()\n", "```\n", "\n", "Nach dem Import der notwendigen Pakete erzeugen wir drei Felder mit je 250\n", "standard-normalverteilten Zufallsgrößen. Durch Multiplikation mit vorgegebenen\n", "Werten für die jeweiligen Unsicherheiten und Addition zu den angenommenen wahren\n", "Werten der Messgrößen mit *m1* und *m2* erhalten wir simulierte Messwerte.\n", "Zusammengefasst in einer *python*-Funktion sieht das so aus:\n", "```\n", "def plot_correlation(w1, w2, sig1, sig2, sigg):\n", "# Zufallszahlen erzeugen\n", " zg = np.random.randn(250)\n", " z1 = np.random.randn(250)\n", " z2 = np.random.randn(250)\n", "# simulierte Messwerte erzeugen\n", " m1 = w1 + sig1*z1 + sigg*zg\n", " m2 = w2 + sig2*z2 + sigg*zg\n", "# graphische Darstellung\n", " scatterplot(m1, m2)\n", "```\n", "\n", "Bitte den obigen Code in einer eigenen Zelle eingeben. \n", "In einer neuen Zelle werden dann die Werte definiert und die grafische Darstellung\n", "aufgerufen:\n", "```\n", "w1 = 5.\n", "sig1 = 1.\n", "w2 = 5. \n", "sig2 = 1.\n", "sigg = 0.\n", "plot_correlation(w1, w2, sig1, sig2, sigg)\n", "```\n", "\n", "Es ist nützlich, den Korrelationskoeffizienten berechnen zu lassen.\n", "Die Elemente der Kovarianzmatrix sind\n", "> $\\mathrm{V}_{11} = \\sigma_1^2 + \\sigma_g^2$ \n", " $\\mathrm{V}_{22} = \\sigma_2^2 + \\sigma_g^2$ \n", " $\\mathrm{V}_{12} = V_{21} = \\sigma_g^2\\,$;\n", "\n", "damit ist der Korrelationskoeffizient $\\rho$ der beiden Messungen\n", "> $\\displaystyle\\rho = \\frac{\\mathrm{V}_{12}}{\\sqrt{\\mathrm{V}_{11}*\\mathrm{V}_{22}}}\\,$.\n", "\n", "Der folgende Code erledigt die Berechnung von $\\rho$:\n", "```\n", "rho = sigg*sigg / np.sqrt( (sig1**2 + sigg**2) * (sig2**2+sigg**2))\n", "print ('Korrelationskoeffizient rho = ', rho)\n", "```\n", "\n", "Es ist hilfreich, diese Zeilen vor der Darstellung der Grafik einzufügen,\n", "um den Wert des Korrelationskoeffizienten mit der Form der Punktewolke \n", "in Verbindung zu setzen. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# eigenen Code hier zur Erzeugung der Zufallszahlen\n", "# und zur Darstellung des Streudiagramms hier eingeben:\n", "# -->" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Geben Sie in der folgenden Code-Zelle verschiedene Werte für w1, w1, sig1, sig2 und sigg vor\n", "(in der Größenordnung von z. B. 0., 0.5, 1.0, 1.5), berechnen Sie jeweils den Korrelationskoeffizienten\n", "$\\rho$ und schauen sich das Streudiagramm für verschiedene Werte von $\\rho$ an." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# hier Werte für w1, w2, sig1, sig2 und den gemeinsamen Fehler sigg eingeben\n", "# -->\n", "\n", "\n", "# hier Formel zur Berechnung des Korrelationskoeffizienten eingeben\n", "# -->" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Dokumentation des Ergebnisses** \n", "Beschreiben Sie hier kurz die Form der Punktewolke im Streudiagramm für\n", "> rho = 0. \n", " rho = 1. \n", " rho = 0.5 \n", " rho = 0.75 \n", " rho = 0.9 \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 2.4.3 Konstruktion der Kovarianzmatrix\n", "---\n", "\n", "\n", "Wir fassen kurz die Eigenschaften der Kovarianzmatrix ${\\bf V}$ zusammen:\n", "\n", " - ${\\bf V} = \\left(\\mathrm{V}_{ij}\\right)$ ist eine symmetrische Matrix\n", " - sie hat die Dimension $N$, die der Anzahl der Messwerte entspricht\n", " - für unabhängige Messwerte ist die Matrix diagonal\n", " - die Nebendiagonalelemente $V_{ij},\\,{\\small i\\ne j}$ lassen sich verstehen als\n", " das Produkt der gemeinsamen Unsicherheiten $\\sigma^g_i$ und $\\sigma^g_j$ der \n", " Messungen $i$ und $j$. \n", "\n", "Der letzte Punkt verdient eine kurze *Diskussion*: \n", ">Wir hatten das Kovarianzelement folgendermaßen definiert: \n", ">\n", ">$\\mathrm{V}_{ij}=\\left< (x_i - \\left< x_i\\right>) \\, (x_j - \\left< x_j\\right>) \\right>\\,.$\n", ">\n", ">Wir denken uns nun den Zufallswert der $x_{i,j}$ zusammengesezt aus einem wahren Wert $x^w_{i,j}$, \n", "einer Zufallsgröße $z_{i,j}$ und einer Zufallsgröße $z_g$, die $x_i$ und $x_j$ gemeinsam ist, \n", "also $x_i = x^w_i + z_i + z_g$ und $x_j = x^w_j + z_j + z_g$. \n", ">Weiter gilt $\\left< x_{i,j}\\right> = x^w_{i,j}$, der Erwartungswert der $x_{i,j}$ entspricht also dem wahren Wert. \n", ">\n", ">Durch Einsetzen in die Formel für $\\mathrm{V}_{ij}$ ergibt sich:\n", ">\n", ">$\\mathrm{V}_{ij}=\\left< (x^w_i + z_i + z_g - x^w_i) \\, (x^w_j + z_j + z_g - x^w_j) \\right>\n", " = \\left< (z_i + z_g) \\, (z_j + z_g) \\right> $\n", "> \n", ">Durch Ausmultiplizieren erhalten wir:\n", ">\n", ">$\\mathrm{V}_{ij}=\\left< z_i z_j + z_i z_g + z_j z_g + z_g^2 \\right>\n", " = \\left< z_i z_j \\right>+\\left+\\left+\\left$\n", "> \n", ">Da $z_i$, $z_j$ und $z_g$ unabhängige Zufallszahlen sind, verschwinden die Erwartungswerte ihrer\n", "Produkte. Lediglich der letzte Term ist ungleich Null. Damit ergibt sich schließlich\n", ">\n", ">$\\mathrm{V}_{ij}\\,=\\,\\left \\,=\\, \\sigma_g^2\\,.$\n", ">\n", "> Bei relativen Unsicherheiten gibt es eine weitere Besonderheit: die vollständig korrelierten\n", " Unsicherheiten können dennoch unterschiedliche Größen haben. Ein Beispiel dafür sind Längenmessungen\n", " mit einer Skalenunsicherheit - vollständige Korrelation bedeutet dann, dass zwei mit demselben \n", " Gerät gemessene Werte um den gleichen relativen Anteil vom wahren Wert abweichen. In der Betrachtung\n", " von eben nehmen wir nun an, das $z_g$ standard-normalverteilte Zufallszahlen sind und\n", " schreiben $z^g_1 = z_g \\cdot \\sigma^g_1$ bzw. $z^g_2 = z_g \\cdot \\sigma^g_2$. Der letzte Term nimmt \n", " dann die folgende Form an:\n", "> \n", ">$\\sigma^g_1 \\, \\sigma^g_2 \\left< z_g^2 \\right> = \\sigma^g_1 \\sigma^g_2 $\n", ">\n", ">Damit ist die Gültigkeit des letzten Punktes in der Liste oben gezeigt.\n", " \n", " Wir sind nun in der Lage, eine Konstruktionsvorschrift für die Kovarianzmatrix anzugeben:\n", " - Unsicherheiten der Messwerte werden nach verschiedenen Quellen aufgeschlüsselt\n", " - unabhängige Unsicherheiten jeder einzelnen Messung\n", " - allen oder Gruppen von Messungen gemeinsame absolute oder relative Unsicherheiten\n", " - die Quadrate der einzelnen Unsicherheiten werden in allen Elementen $\\mathrm{V}_{ij}$ der Kovarianz-Matrix\n", " aufaddiert, die von dieser Unsicherheit betroffen sind\n", " - unabhängige Unsicherheiten stehen nur in den Diagonalelementen $\\mathrm{V}_{ii}$;\n", " - korrelierte Unsicherheiten stehen zusätzlich in den betroffenen Nebendiagonalelementen.\n", " \n", "Ein konkretes Beispiel dazu wird weiter unten untersucht. \n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.4.4 Das Fehlerfortpflanzungsgesetz mit korrelierten Messunsicherheiten\n", "---\n", "\n", "\n", "Bei der Herleitung des Fehlerforpflanzungsgesetzes für unabhängige, d.h. damit auch\n", "unkorrelierte Messungen hatten wir gesehen, dass ein weitere Term auftritt, die Kovarianz der\n", "Messwerte, den wir seinerzeit zu Null gesetzt hatten (siehe [Gleichung (2.1)](#Equ-2.1)). \n", "In den meisten Fällen führen aber allen Messungen einer Messreihe gemeinsame\n", "systematische Unsicherheiten dazu, dass die Messwerte korreliert sind. Wir haben\n", "oben gesehen, wie man diese Korrelationen mit Hilfe der Kovarianz oder des \n", "Korrelationskoeffizienten quantifiziert. Wenn man den Kovarianz-Term in [Gleichung (2.1)](#Equ-2.1)\n", "nicht vernachlässigt, erhält man mit Hilfe der Kovarianzmatrix das Fehlerfortpflanzungsgesetz \n", "in allgemeiner Form:\n", "\n", "$y = y(x_1, \\ldots, x_n)$ \n", "\n", "> $\\Rightarrow \\sigma^2_{y} \\simeq\n", " \\displaystyle \\sum_{i=1}^n \\displaystyle\\sum_{j=1}^n\n", "\\left( \\frac{\\partial y}{\\partial x_i} \\right)\n", " \\mathrm{V}_{ij}\n", "\\left( \\frac{\\partial y}{\\partial x_j} \\right)\\,$\n", " Gleichung (2.4)\n", "\n", "Leichter zu merken ist vielleicht die kompaktere Darstellung in Vektor-Schreibweise\n", "> $ {\\large \\sigma^2_{y} \\simeq\n", " (\\vec{\\nabla}y({\\small \\vec{x}}))^T \\, {\\bf V} \\,\\vec{\\nabla}y({\\small \\vec{x}}) }\\,$\n", " Gleichung (2.5)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Übung** zu 2.4.4 \n", "\n", "Zeigen Sie, dass diese Darstellung für unkorrelierte Messungen, d. h. für eine Kovarianzmatrix mit verschwindenden Nebendiagonalelementen, zur früheren Formel für unabhängige Messungen ([Gleichung (2.3)](#Equ-2.3)) äquivalent ist." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Kurze eigene Begründung:** \n", "-> \n", " \n", " \n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Die schon früher betrachteten Spezialfälle für zwei Messwerte sehen für den \n", "allgemein Fall wie folgt aus:\n", "\n", "$y=x_1 \\pm x_2$ \n", "> $\\Rightarrow \\color{blue}{\n", " {\\sigma_y}^2 = {\\sigma_1}^2 + {\\sigma_2}^2 \\pm 2\\, \\mathrm{cov}(x_1,x_2) }$\n", "\n", "\n", "$y=x_1 \\cdot x_2 {\\rm ~oder~} y = \\frac{x_1}{x_2}$ \n", "> $\\Rightarrow \\color{blue} {\n", " \\left( \\frac{\\sigma_y} {y} \\right)^2 \\simeq \n", " \\left( \\frac{\\sigma_1} {x_1} \\right)^2 \n", " + \\left( \\frac{\\sigma_2} {x_2} \\right)^2\n", " \\pm 2\\frac{\\mathrm{cov}(x_1,x_2)}{x_1 x_2} }$\n", "\n", "Ein **Beispiel** soll die Bedeutung des Terms mit der Kovarianz illustrieren:\n", "\n", "> $ y = x_1 - x_2$ \n", " $\\left< x_1 \\right> = \\left< x_2 \\right> = 10.$ \n", " $\\sigma_1 = \\sigma_2 = 1. $ \n", " \n", "Für $\\rho = 0.$ ergibt sich $\\sigma_y = 1.4$ \n", "Aber mit $\\rho = 1. $ erhält man $\\sigma_y = 0.$ !\n", "\n", "__Anmerkung 1__: Dieses Beispiel ist nicht fiktiv, sondern hat tatsächlich \n", "eine wichtige technische Anwendung bei der Signalübertragung. Zusätzlich\n", "zum eigentlichen Signal $s_+$ wird ein zweites, invertiertes Signal $s_- =-s_+$\n", "übertragen (\"differenzielle Signalübertragung\"). \n", "Alle Störungen am Ende der Übertragungsstrecke sind für beide \n", "Signale zu einem hohen Grad korreliert. \n", "Das Nutzsignal wird durch Bilden der Differenz zurückgewonnen: \n", "$s_o = 0.5\\,(s_+ - s_-)$; wegen der fast vollständigen Korrelation der\n", "beiden Signalkomponenten ist die Unsicherheit des extrahierten Signals \n", "sehr klein, $\\sigma_{s_o} \\simeq 0\\,$.\n", "\n", "__Anmerkung 2__: In der älteren Literatur wird oft zur \"konservativen\" Abschätzung\n", "von Messunsicherheiten die \"Größtfehlermethode\" empfohlen. Dabei werden statt der\n", "Wurzel aus der Summe der Quadrate der einzelnen Unsicherheiten die Absolutbeträge\n", "aufaddiert. Dies entspricht einer angenommenen Korrelation von 100% der einzelnen\n", "Komponenten. In der Praxis tritt eine solch starke Korrelation aller Messwerte \n", "allerdings fast nie auf, daher wird die Unsicherheit mit der \"Größtfehlermethode\" \n", "fast immer zu groß abgeschätzt. Wenn Sie hohe Korrelationen der Unsicherheiten vermuten, sollten Sie diese realistisch abschätzen und die Formeln aus\n", "diesem Kapitel verwenden!\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "#### 2.4.5 Korrelierte Variable im *uncertainties*-Paket\n", "\n", "\n", "Das oben bereits eingeführte *python*-Paket *uncertainties* unterstützt auch die Fehlerfortpflanzung\n", "bei korrelierten Unsicherheiten. Die notwendigen imports sind: \n", "```\n", "import uncertainties as unc, numpy as np\n", "```\n", "Die Eingabe der Messgrößen und der korrelierten Unsicherheiten erfolgt mit der Zeile \n", "`(a,b) = unc.correlated_values([1. ,2.], cov_matrix)`\n", "\n", "Vorher muss eine Kovarianzmatrix erzeugt werden:\n", "```\n", "cov_matrix = np.array([ [0.1**2, 0.2*0.1] ,\n", " [0.1*0.2, 0.2**2] ])\n", "```\n", "In diesem Beispiel sind die Unsicherheiten der beiden Messgrößen zu 100% korreliert. \n", "Die korrekte Eingabe überprüft man am besten durch eine Kontrollausgabe:\n", "```\n", "# Ausgabe zur Kontrolle\n", "print('a=', a)\n", "print('b=', b)\n", "print('correlations \\n', unc.correlation_matrix([a,b]))\n", "```\n", "\n", "Nun werden zwei Funktionen, die Summe und Differenz der Größen $a$ und $b$, berechnet und\n", "ausgegeben:\n", "```\n", "# Summe und Differenz von zwei vollständig korrelierten Variablen\n", "print('sum = ', a + b)\n", "print('dif = ', b - a)\n", "```\n", "Geben Sie den Code in der Zelle unten ein. Begründen Sie kurz das Ergebnis." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Code hier eingeben:\n", "# -->" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Das Ergebnis des Programms ist \n", ">sum = 3.00+/-0.30 \n", " dif = 1.00+/-0.10 \n", "\n", "Die Unsicherheiten ergeben in diesem Fall als Summe bzw. Differenz der Unsicherheiten,\n", "> $\\sigma_s = \\sigma_a + \\sigma_b $ \n", " $\\sigma_d = |\\sigma_a - \\sigma_b|$\n", " \n", "Warum ist das korrekt ? \n", "\n", "**Eigene Begründung hier** \n", "-> \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "#### 2.4.6 Konstruktion einer Kovarianzmatrix an einem aufwändigeren Beispiel\n", "---\n", "\n", "\n", "Wir betrachten folgendes **Übungsbeispiel**: \n", "Sechs Personen in drei Zweiergruppen messen mit jeweils einem Messgerät pro Gruppe eine Größe $m_i, {\\small i=1,\\ldots,6}$. Es wird eine theoretische \n", "Korrektur benötigt, die eine Unsicherheit $\\sigma_t$ hat. Die sechs Messwerte sind $m_i$=(0.82, 0.81, 1.32, 1.44, 0.93, 0.99). \n", "Damit gibt es drei Komponenten der Gesamtunsicherheit einer jeden Einzelmessung: \n", "> 1.) individuelle, unabhängige Messunsicherheiten $\\sigma_u = 0.1$, \n", " 2.) Abweichungen der Messgeräte, die zu einer systematischen Unsicherheit $\\sigma_{s}=0.15$ führen, vollständig korreliert in jeder der drei Zweiergruppen, \n", " 3.) eine Theorieunsicherheit $\\sigma_t=0.05$, die für alle Messungen vollständig korreliert ist. \n", " \n", "Da alle Messungen von gleicher Qualität und Genauigkeit sind, \n", "wird als Endergebnis der Mittelwert $\\bar{m}$ aller Messungen gebildet.\n", "Wie groß ist die Unsicherheit des Mittelwerts $\\sigma_\\bar{m}$?\n", " \n", "Damit wir die oben hergeleitete Formel verwenden können, konstruieren wir zunächst die\n", "Kovarianzmatrix der Messungen $m_i$. \n", "Hier noch einmal zusammengefasst die Regeln:\n", " - die Diagonalelemente $\\mathrm{V}_{ii}$ sind jeweils der quadrierte Gesamtfehler der Messung $i$\n", " - die Nebendiagonalen $\\mathrm{V}_{ij},\\,{\\small i\\ne j}$ enthalten die quadrierten, gemeinsamen\n", " Unsicherheiten der Messungen $i$ und $j$ \n", "\n", "Die Gesamtunsicherheit ergibt sich zu ${\\sigma_g}^2 = {\\sigma_u}^2 + {\\sigma_s}^2 + {\\sigma_t}^2$. \n", "Die Matrix sieht damit so aus:\n", ">$$\\displaystyle{{\\bf V}\\,=}\n", "\\left(\n", "\\begin{array}{cccccc}\n", "\\sigma_g^2 & \\sigma_s^2+\\sigma_t^2 & \\sigma_t^2 & \\sigma_t^2 & \\sigma_t^2 & \\sigma_t^2 \\\\\n", "\\sigma_s^2+\\sigma_t^2 & \\sigma_g^2 & \\sigma_t^2 & \\sigma_t^2 & \\sigma_t^2 & \\sigma_t^2 \\\\\n", "\\sigma_t^2 & \\sigma_t^2 & \\sigma_g^2 & \\sigma_s^2+\\sigma_t^2 & \\sigma_t^2 & \\sigma_t^2 \\\\\n", "\\sigma_t^2 & \\sigma_t^2 & \\sigma_s^2+\\sigma_t^2 & \\sigma_g^2 & \\sigma_t^2 & \\sigma_t^2 \\\\\n", "\\sigma_t^2 & \\sigma_t^2 & \\sigma_t^2 & \\sigma_t^2 & \\sigma_g^2 & \\sigma_s^2+\\sigma_t^2 \\\\\n", "\\sigma_t^2 & \\sigma_t^2 & \\sigma_t^2 & \\sigma_t^2 & \\sigma_s^2+\\sigma_t^2 & \\sigma_g^2 \\\\\n", "\\end{array} \\right) $$\n", "\n", "Man sieht deutlich die drei Blöcke, die den drei verwendeten Messgeräten entsprechen. \n", "Die Komponente $\\sigma_t$ steckt in allen Elementen der Kovarianzmatrix. \n", "\n", "Wir berechnen nun die Unsicherheit des Mittelwerts $\\bar{m}=\\frac{1}{6}\\sum_{j=1}^6 m_j$\n", "mit Hilfe von [Gleichung (2.5)](#Equ-2.5). Die Berechnung des Gradienten von $\\bar{m}$ ergibt\n", "> $(\\vec{\\nabla}\\bar{m})^T = \\frac{1}{6}\\,(1,1,1,1,1,1)$. \n", "\n", "Ausführen der Matrix-Multiplikation mit dem Gradienten von rechts und links ergibt\n", "> $\\sigma_\\bar{m}^2 = \\frac{1}{36}\\sum_{i,j} \\mathrm{V}_{ij}$, also die Summe über alle Elemente der Kovarianzmatrix. \n", "\n", "Wegen der recht einfachen Struktur der Kovarianzmatrix ist diese Summe leicht auszuführen:\n", "$\\sigma_u^2$ kommt 6 mal vor, $\\sigma_s^2$ 12 mal und $\\sigma_t^2$ 36 mal. Damit ergibt sich \n", "> $$\\sigma_\\bar{m}^2 \\,=\\, \\sigma_t^2 \\,+\\, \\sigma_s^2/3 \\,+\\, \\sigma_u^2/6 = (0.108)^2\\,, \\, {\\rm also}\\, \\sigma_\\bar{m} = 0.108$$ \n", " \n", "Insgesamt ergibt sich also für den Mittelwert $\\bar{m}=1.052\\pm 0.108 \\,$.\n", "\n", "__Intuitive Überlegungen__\n", "Dieses Ergebnis kann man auch intuitiv begründen. Wir erinnern uns, dass sich bei mehrmaligen, \n", "unabhängigen Messungen die Unsicherheit um einen Faktor reduziert, der der Wurzel aus der Anzahl\n", "der Messugen entspricht. Damit wird die statistische Unsicherheit um den Faktor $\\sqrt{6}$ kleiner\n", "und die systematische, gerätebedingte Unsicherheit um den Faktor $\\sqrt{3}$; die theoretische\n", "Unsicherheit bleibt unverändert. Addiert man die drei Komponenten quadratisch, \n", "so ergibt sich $\\sigma_\\bar{m} = 0.108\\,$.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Übung** zum Beispiel aus [Abschnitt 2.4.6](#SubSec-BeispielKovarianzmatrix)\n", "\n", "Da die Berechnung \"von Hand\" im allgemeinen mühsam ist, nutzt man zur Konstruktion\n", "der Kovarianzmatrix den Computer und überlässt die Berechnung dann dem \n", "oben eingeführten *python*-Paket *uncertainties*.\n", "\n", "Zur Konstruktion der Kovarianzmatrix helfen folgende Überlegungen:\n", " \n", "Die Elemente der Kovarianzmatrix ${\\bf V}_c$, die vollständig korrelierten Unsicherheiten $\\sigma^c_i$ \n", "der Messungen $i, j,$ entsprechen, erhält man durch Bilden des dyadischen Produkts (manchmal auch: \"äußeres Produkt\", \"Matrixprodukt\") von Vektoren, bei denen nur die Elemente $i$ und $j$ ungleich Null sind:\n", " \n", " $\\vec{\\sigma}_c = (0., 0., \\ldots, \\sigma_i, \\sigma_j, 0., \\ldots, 0.)^T\n", " {\\quad \\rightarrow {\\bf V}_c = \\vec{\\sigma}_c {\\vec{\\sigma}_c}^T} $\n", " (\"Matrixprodukt von Spaltenvektor mit Zeilenvektor\") \n", "\n", "Die Kovarianzmatrix ergibt sich als Summe von solchen dyadischen Produkten für alle Arten von\n", "Unsicherheiten. Im Paket *numpy* ist das dyadische Produkt über die Funktion *numpy.outer(a, b)* verfügbar. \n", "\n", "Die Diagonalelemente der Kovarianzmatrix für die unabhängigen Unsicherheiten \n", "kann man mit der Funktion *numpy.diagflat() setzen. \n", "\n", "Eine Funktion zur komfortablen Konstruktion der Kovarianz-Matrix sieht z.B. so aus:\n", "```\n", "def BuildCovarianceMatrix(sig, sigc=[]):\n", " '''\n", " Construct a covariance matrix from independent and correlated uncertainty components\n", "\n", " Args: \n", " * sig: iterable of independent uncertainties \n", " * sigc: list of iterables of correlated uncertainties\n", " \n", " Returns: \n", " covariance matrix as numpy-array\n", " '''\n", " \n", " # construct a numpy array with diagonal elements \n", " su = np.array(sig)\n", " V = np.diagflat(su*su)\n", "\n", " # if given, add off-diagonal components\n", " if sigc != []:\n", " for s in sigc:\n", " sc = np.array(s)\n", " V += np.outer(sc, sc)\n", " return V\n", "```\n", "Diese Funktion ist dem Paket *PhyPraKit* entnommen; dort gibt es auch Funktionen zur\n", "Konversion einer Kovarianz-Matrix in eine Korrelationsmatrix mit Gesamtunsicherheiten \n", "und umgekehrt.\n", "\n", "Die Anwendung auf die gegebenen Daten sieht etwa wie folgt aus:\n", "```\n", "import uncertainties as unc, numpy as np, PhyPraKit as ppk\n", "\n", "# the data\n", "m = np.array([0.82, 0.81, 1.32, 1.44, 0.93, 0.99])\n", "sig_u = 0.10 # independent uncertainties\n", "sig_c = 0.05 # fully correlated uncertainties \n", "sig_s = 0.15 # systematic uncertainties for pairs of measurements\n", "\n", "# construct vectors with error components\n", "N = len(m)\n", "su = sig_u * np.ones(N)\n", "sc = sig_c * np.ones(N)\n", "ss1 = np.array([sig_s, sig_s, 0., 0., 0., 0.])\n", "ss2 = np.array([0., 0., sig_s, sig_s, 0., 0.])\n", "ss3 = np.array([0., 0., 0., 0., sig_s, sig_s])\n", "V = ppk.BuildCovarianceMatrix(su, [sc, ss1, ss2, ss3])\n", "\n", "print('\\nconstructed covariance matrix: \\n', V)\n", "err, C = ppk.Cov2Cor(V)\n", "print('\\ntotal uncertainties:\\n', err)\n", "print('correlation matrix: \\n', C)\n", "```\n", "\n", "**Aufgabe:** Bauen Sie die Kovarianzmatrix für das Beispiel in \n", "[Abschnitt 2.4.6](#SubSec-BeispielKovarianzmatrix) auf und \n", "nutzen Sie das Paket *uncertainties* zur Berechnung des Mittelwerts und seiner Unsicherheit.\n", "*Tipp:* *uncertainties* unterstützt *numpy*-*arrays*! Man kann also Variable vom Typ\n", "*ufloat* mit *numpy*-Fuktionen nutzen. Hier ein Code-Fragment dazu:\n", "```\n", "import uncertainties as unc, numpy as np\n", "unc_m = np.array( unc.correlated_values(m, V) )\n", "mbar = unc_m.sum()/6.\n", "print('Mittelwert und Unsicherheit: ', mbar)\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# eigenen Code hier eingeben\n", "# -->" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "---\n", "## 3. Anpassung von Modellen an Messdaten und Parameterschätzung [^](#Inhalt)\n", "----\n", "\n", "Der Vergleich von Modellen mit Messungen gehört zu den Standardaufgaben in der\n", "Experimentalphysik. Im einfachsten Fall stellt ein Modell eine Funktion dar, die\n", "Vorhersagen für die in einem Experiment zu erwartende Messdaten liefert. Meist sind\n", "diese Modelle keine absoluten Vorhersagen, sondern hängen von Modellparametern ab. \n", "Neben der Überprüfung, ob das Modell die Daten überhaupt beschreibt, gehört die\n", "Bestimmung der Modellparameter zu den typischen Aufgaben. \n", "\n", ">Ein Beispiel für ein solches Modell ist der lineare Zusammenhang zwischen Stromstärke $I$ und \n", "Spannung $U$ an einem (idealen) ohmschen Widerstand. Der Modellparameter \n", "ist der elektrischer Leitwert $G$, Kehrwert des Widerstands $R$. Das Modell \n", "$I(U) = G\\cdot U$ entspricht geometrisch einer Ursprungsgeraden im $I$-$U$-Diagramm. \n", "Eine solche Gerade wird den gemessenen Stromstärken bei gegebener Spannung angepasst. \n", "Die Steigung entspricht dem Modellparameter $G$.\n", "\n", "Zur Funktionsanpassung wird häufig die Methode der kleinsten Quadrate\n", "verwendet, mit der sich sehr elegant auch kompliziertere Situationen wie korrelierte \n", "Unsicherheiten der Messdaten oder Unsicherheiten in Abszissenrichtung zusätzlich zu denen\n", "in Ordinatenrichtung behandeln lassen. \n", "\n", "Dieses Kapitel gibt einen kurzen Abriss der Methode und enthält praktische Hinweise \n", "und Übungen zur Funktionsanpassung mit Hilfe numerischer Methoden auf dem Computer.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.1 Funktionsanpassung mit der Methode der kleinsten Fehlerquadrate\n", "\n", "\n", "Zunächst muss die Übereinstimmung einer Funktion mit Messdaten quantifiziert werden. \n", "Ausgangspunkt sind $N$ Messpunkte $(x_i, y_i)$ und eine Funktion $f(x_i)$, die die Werte\n", "der $y_i$ möglichst genau vorhersagen soll. Zur Definition eines Abstands der Funktion von\n", "den Messpunkten gibt es einen Vorschlag von Carl Friedrich Gauß: \n", "$d := \\sum (f(x_i) - y_i)^2$ (\"Methode der Gauß'schen Fehlerquadrate\").\n", "\n", "Für Messdaten $y_i$ mit bekannten Unsicherheiten $\\sigma_i$ wird der quadrierte Abstand\n", "auf das Quadrat der Unsicherheit normiert. Beachten Sie, dass nur die Werte der Ordinate $y_i$ eine\n", "Unsicherheit besitzen. Die Ordinatenwerte $x_i$ werden hier als exakt angenommen.\n", "Außerdem soll die Funktion $f$ noch von einem Vektor $\\vec{p} = (p_1, \\ldots, p_k)$ \n", "der $k$ Modellparameter abhängig sein. Dann ergibt sich als \n", "parameterabhängiges Abstandsmaß die sogenannte Residuensumme \n", "\n", "> $ S = \\displaystyle \\sum_{i=1}^N \\left( \\frac{ f(x_i; \\vec{p}) - y_i} {\\sigma_i} \\right)^2\\,$ \n", " Gleichung (3.1)\n", "\n", "Nachdem eine Messung erfolgt ist, liegen die Messpunkte $(x_i, y_i)$ fest. Daher kann man die Summe der Residuen $S$ aus [Gleichung (3.1)](#Equ-3.1) als Funktion im $k$-dimensionalen Raum der Parameter $p_j$\n", "auffassen. Durch Minimierung von $S(\\vec{x}; \\vec{p})$ bezüglich der Parameter $p_j$ erhält man\n", "den Parametersatz $\\hat{\\vec{p}}$, für den die gewählte Funktion die Daten am besten beschreibt. \n", "Das Hutsymbol ^ deutet an, dass es sich um die besten Schätzwerte (nach der Methode der kleinsten\n", "Fehlerquadrate) für die Parameter handelt. \n", "\n", "Diese Vorschrift lässt sich aus einem allgemeinen Prinzip der mathematischen Stochastik,\n", "dem Likelihood-Prinzip ableiten, wenn die Messunsicherheiten einer Gauß-Verteilung unterliegen. \n", "\n", "Die Minimierung von $S$ kann in Spezialfällen analytisch erfolgen.\n", "Dazu werden die partiellen Ableitungen von $S$ nach den Parametern $p_j$ zu Null gesetzt,\n", " $\\frac{\\partial S} {\\partial p_j}=0 \\, , \\,{\\small j=1, \\ldots , k}\\,$.\n", "Für Probleme, die linear in den Parametern sind, also \n", "$f(\\vec{x}; \\vec{p})=\\sum_{j=0}^k g_j(\\vec{x})\\, p_j$ mit beliebigen Funkionen $g_j$, \n", "gibt es eine analytische Lösung. Im allgemeinen verwendet man Methoden der numerischen\n", "Optimierung zur Suche des Minimums im $k$-dimensionalen Parameterraum. In der Praxis\n", "werden heute oft auch lineare Probleme durch numerische Verfahren gelöst.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "#### 3.1.1 Beispiel einer Funktionsanpassung \n", "\n", "\n", "Als einfaches Beispiel betrachten wir zehn Messwerte $y_i$, jeweils mit der gleichen Unsicherheit $\\sigma_i = \\sigma$. \n", "\n", "```\n", "y = [11.55, 9.8, 9.82, 9.15, 10.57, 9.58, 10.44, 10.55, 8.23, 10.93]\n", "sig_y = 1. \n", "```\n", "An diese Daten soll eine konstante Funktion $f(x; c) = c$ angepasst werden - dies entspricht\n", "einer Mittelung der Daten. \n", "\n", "Das Abstandsquadrat scheibt sich in *Python* wie folgt:\n", "```\n", "import numpy as np\n", "y = np.array(y) # konvertieren in Numpy-Array\n", "def S(y, c):\n", " return ( ( ( y - c ) / sig_y )**2).sum()\n", "```\n", "\n", "Die Minimierung kann sehr einfach ausgeführt werden, indem man $S$ für verschiedene\n", "Werte von $c_j$ bestimmt:\n", "```\n", "S_c = []\n", "c_values = np.linspace(8, 12, 100)\n", "for c in c_values:\n", " S_c.append(S(y, c))\n", "```\n", "\n", "Jetzt muss nur noch das Mimimum gefunden werden. Dabei hilft die Funktion `np.argmin()`:\n", "```\n", "id_min = np.argmin(S_c)\n", "print(\"Mimimum bei \", c_values[id_min])\n", "```\n", "Damit ist das Minimierungsproblem numerisch gelöst! \n", "\n", "Illustrativ ist noch eine grafische Darstellung der Abhängigkeit $S(c)$:\n", "```\n", "# graphische Darstellung S(c)\n", "import matplotlib.pyplot as plt\n", "plt.figure()\n", "plt.plot(c_values, S_c)\n", "plt.xlabel('c')\n", "plt.ylabel('S(c)')\n", "plt.show()\n", "```\n", "\n", "**Übung zu 3.1.1** \n", "Fügen Sie die Code-Fragmente zusammen und führen Sie das Programm in der Zelle unten aus.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# eigenen Code hier eingeben:\n", "# -->" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "#### 3.1.2 Analytische Lösung \n", "\n", "\n", "Das eben betrachtete einfache Problem lässt sich natürlich auch analytisch lösen. \n", "\n", "- $S(c)=\\displaystyle \\sum_{i=1}^{N=10} \\frac{(y_i - c)^2} {\\sigma^2}$\n", "- $0 = \\frac{dS}{dc} = \\displaystyle \\sum_{i=1}^{N} \\frac{-2 (y_i -c)}{\\sigma^2} $\n", " > $\\Rightarrow \\displaystyle {\\hat c}=\\frac{1}{N}\\sum_{i=1}^{N} y_i $\n", " \n", " \n", "Es ist wenig erstaunlich, dass wir als Ergebnis genau die Formel für den Mittelwert erhalten. \n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "#### 3.1.3 Mittelwert von Messungen mit unterschiedlichen Unsicherheiten\n", "\n", "Wir nehmen nun wieder ${\\small N=10}$ Messwerte $y_i$, diesmal aber mit unterschiedlichen\n", "Unsicherheiten $\\sigma_i$:\n", "```\n", "y = [11.55, 9.8, 9.82, 9.15, 10.57, 9.58, 10.44, 10.55, 8.23, 10.93]\n", "sig_y = [0.8, 0.9, 0.9, 1.1, 1.0, 1.2, 0.7, 1.1, 1.0, 0.9] \n", "```\n", "Wir erwarten, dass in diesem Fall Messwerte mit kleineren Unsicherheiten \n", "\"wichtiger\" für den Mittelwert sind als Messwerte mit größeren Unsicherheiten.\n", "\n", "Betrachten wir die analytische Vorgehensweise von eben:\n", "\n", " - $S(c) = \\displaystyle\\sum_{i=1}^{N} \\frac{(y_i -c)^2}{{\\sigma_i}^2}$\n", " - $0= \\frac{dS}{dc} = \\displaystyle\\sum_{i=1}^{N} \\frac{2 (y_i -c)}{{\\sigma_i}^2}$\n", " > $\\Rightarrow {\\hat c} = \\frac{1}{\\sum{1/\\sigma_i}^2}\n", " \\displaystyle\\sum_{i=1}^{N}\\frac{1}{{\\sigma_i}^2}y_i $ \n", "\n", "Mit der Definition $ w_i \\equiv \\frac{1} {{\\sigma_i}^2}$ erhalten wir als wichtiges\n", "Ergebnis die Formel für den *gewichteten Mittelwert*:\n", "> ${\\hat c}= \\frac{1} {\\sum {w_i}} \\displaystyle\\sum_{i=1}^{N} {w_i} y_i$ \n", "\n", "Der Mittelwert ist die mit $1/\\sigma_i^2$ gewichtete Summe der Einzelmessungen. \n", "Die Gewichte entsprechen der obigen Erwartung: Messwerte mit den kleinsten $\\sigma_i$\n", "bekommen bei der Mittelwertbildung das größte Gewicht. Am Vorfaktor $1/\\sum w_i$ kann\n", "man ablesen, dass die Zahl der Messwerte $N$ in diesem Fall durch die Summe der Gewichte\n", "$\\sum w_i$ ersetzt werden muss.\n", " \n", "Das Paket *PhyPraKit* enthält eine Funktion zur Berechnung des gewichteten Mittelwerts: *wmean(x, sx)*.\n", "\n", "**Übung zu 3.1.2** Probieren Sie *PhyPraKit.wmean()* in der Code-Zeile unten aus: " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from PhyPraKit import wmean\n", "\n", "# eigenen Code hier eingeben\n", "# -->" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "### 3.2 Bestimmung der Parameterunsicherheiten \n", "---\n", "\n", "Nachdem wir einige Beispiele für analytische und numerische Anpassungen gesehen haben, widmen\n", "wir uns nun der Frage, wie man die Unsicherheiten der Parameter bestimmen kann. \n", "\n", "Wieder soll ein konkretes numerisches Beispiel den Einstieg bilden. Wir kommen dazu \n", "auf [Beispiel 3.1.1](#SubSec-BeispielFunktionsanpassung) zurück,\n", "die Anpassung einer Konstanten an Messdaten. Wir setzen alle Messdaten auf einen festen Wert,\n", "z.B. 10, und studieren den Verlauf von $S(c)$ als Funktion des Werts von $c$ für verschiedene \n", "Unsicherheiten $\\sigma_y$. Damit das wiederholte Aufrufen von Berechnung und grafischer Darstellung\n", "übersichtlicher wird, sind im folgenden Beispiel entsprechende Funktionen definiert:\n", "```\n", "import numpy as np, matplotlib.pyplot as plt\n", "\n", "def calc_S(cs):\n", " ''' Berechnung der Residuen-Summen'''\n", " Sc=[]\n", " for c in cs:\n", " Sc.append((((y - c)/sig_y)**2).sum())\n", " return np.array(Sc)\n", "\n", "def plot_S(x, y):\n", " ''' Grafische Darstellung '''\n", " plt.figure()\n", " plt.plot(x, y)\n", " plt.xlabel('c')\n", " plt.ylabel('S(c)')\n", " plt.ylim(0, 10.)\n", " plt.show()\n", "```\n", "\n", "Im Hauptteil des Programms werden nun noch die Daten gesetzt und\n", "dann die Grafik für verschiedene Werte der Unsicherheit ausgeführt:\n", "```\n", "y = 10*np.ones(10)\n", "c_values = np.linspace(9, 11, 100)\n", "\n", "sig_y = 1.\n", "plot_S(c_values, calc_S(c_values))\n", "\n", "sig_y = 0.5\n", "plot_S(c_values, calc_S(c_values))\n", "\n", "sig_y = 0.2\n", "plot_S(c_values, calc_S(c_values))\n", "```\n", "\n", "Geben Sie den Code in die Zelle unten ein, führen ihn aus und schauen die Grafiken an." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Code von oben hier eingeben:\n", "# -->" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Wir beobachten, dass $S(c)$ einen parabolischen Verlauf hat, wobei das Minimum der Parabel\n", "mit kleiner werdenden Unsicherheiten $\\sigma_y$ immer schärfer wird. Offensichtlich bedeutet\n", "ein schärferes Minimum, d. h. eine größere Krümmung, dass der Parameter genauer bestimmt ist. \n", "Die Krümmung einer Parabel am Minimum ist gegeben durch die 2. Ableitung. Wir vermuten daher,\n", "dass die Krümmung direkt mit dem Kehrwert der Varianz verknüpft ist:\n", "\n", "$\\displaystyle{\\frac{1}{\\mathrm{V}_c} = \\frac{1}{\\sigma_c^2} \\propto \\left. \\frac{\\partial^2 S(c)}{\\partial c^2} \\right|_\\hat{c}}$\n", "\n", "Über den Wert der Proportionalitätskonstanten können wir an dieser Stelle keine allgemeine Aussage treffen.\n", "Wenn wir allerdings die analytische Lösung im Beispiel aus [Abschnitt 3.1.2](#SubSec-AnalytischeLoesung) \n", "zu Rate ziehen und die zweite Ableitung bilden, erhalten wir \n", "\n", "$\\displaystyle\\left. \\frac{\\partial^2 S(c)}{\\partial c^2} \\right|_\\hat{c} = \\frac{2N}{\\sigma_y^2}$. \n", "\n", "Da der Schätzwert $\\hat{c}$ der Mittelwert der Größen $y_i$ ist, können wir auf frühere Betrachtungen zurückgreifen und \n", "feststellen, dass $\\sigma_{\\hat c}^2 = \\frac{\\sigma_y^2}{N}$ ist. \n", "In diesem Beispiel ergibt sich der Proportionalitätsfaktor durch Koeffizientenvergleich also zu $\\frac{1}{2}$. Man kann zeigen, dass dies auch allgemein gilt. \n", "\n", "Wenn es mehrere Parameter $p_j$ gibt, lässt sich das Ergebnis weiter verallgemeinern. Die Matrix der partiellen Ableitungen bestimmt die entsprechenden Elemente der Inversen der Kovarianzmatrix der Parameter:\n", "> $\\displaystyle{\\left( {{\\bf V}_\\hat{p}}^{-1} \\right)_{ij} \n", " = \\frac{1}{2}\\,\\left.\\frac{\\partial^2 S(\\vec{p})}{\\partial p_i\\partial p_j}\n", " \\right|_{\\hat{p_i}\\hat{p_j}}}$\n", " Gleichung (3.2)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "### 3.3 Beurteilung der Qualität einer Anpassung \n", "---\n", "\n", "Bei jeder Anpassung stellt sich die Frage, ob das gewählte Modell überhaupt passt. In der Physik ist die\n", "Form des Modells sehr häufig durch theoretische Überlegungen vorgegeben, oder eine Messung wird explizit\n", "dazu durchgeführt, ein angenommenes Modell zu überprüfen. Daher ist die Frage, \"wie gut\" die Daten durch\n", "das Modell beschrieben werden, von zentraler Bedeutung. \n", "\n", "Wir nähern uns dieser Frage, indem wir die zur Anpassung verwendete Residuensumme, \n", "[Gleichung (3.1)](#Equ-3.1), noch einmal genau anschauen: \n", "$ S = \\displaystyle \\sum_{i=1}^n \\left( \\frac{ y_i - f(x_i; \\vec{p}) } {\\sigma_i} \\right)^2\\,$ \n", "\n", "Wenn die Unsicherheiten gaußförmig sind und das Modell richtig ist, entspricht die Modellvorhersage \n", "$f(x_i; \\vec{p})$ dem wahren Wert $y_i^w$. Der Term in der Summe wird dann zu $((y_i - y_i^w)/\\sigma_i)^2$.\n", "Wenn $y_i$ gaußverteilt ist, dann ist $(y_i -y_i^w)/\\sigma_i$ symmetrisch um Null verteilt und\n", "die Standardaweichung ist Eins; es handelt sich also um eine standard-normalverteilte Größe. \n", "Die Verteilung der Summe der Quadrate von standard-normalverteilten Zufallszahlen nennt man\n", "$\\chi^2$-Verteilung mit $n_f$ Freiheitsgraden. Sie ist analytisch bekannt; ihr Erwartungswert \n", "ist $n_f$ und ihre Varianz 2$n_f$.\n", "\n", "Die Zahl der Freiheitsgrade in einer Anpassung mit $N$ Datenpunkten und $k$ Parametern ist $n_f=N-k$. Der Wert von $S$ am Mimimum, also für die besten Schätzwerte der Parameter, folgt einer $\\chi^2$-Verteilung.\n", "\n", "Mit dieser Kenntnis lässt sich mit Hilfe der $\\chi^2$-Verteilung die Wahrscheinlichkeit angeben, einen \n", "noch größeren Wert von $S$ als den tatsächlich beobachteten zu finden. Ist diese Wahrscheinlichkeit klein,\n", "so ist die Anpassung schlecht; erwartet wird ein typischer Wert von 0.5.\n", "\n", "\n", "**Übung** zu [Abschnitt 3.3](#Sec-Fitqualität)\n", "\n", "Die $\\chi^2$-Verteilung ist im Paket *scipy* enthalten, *scipy.stat.chi2.pdf(S, nf)*, die kumulative Verteilungsfunktion dazu ist *scipy.stats.chi2.cdf(S, nf)*. Mit deren Hilfe berechnet man die $\\chi^2$-Wahrscheinlichkeit mittels `chi2prob = 1. - scipy.stat.chi2.cdf(S, nf)`. Das Paket *PhyPraKit* enthält die Funktion *chi2prob()* für diesen Zweck.\n", "\n", "a) Schauen Sie sich mit Hilfe der genannten Funktionen den Verlauf der $\\chi^2$-Verteilung für verschiedene Anzahlen von Freiheitsgraden $(n_f=2$ bis $n_f=10)$ an. \n", "\n", "Hier eine Hilfestellung zum Code:\n", "```\n", "# grafische Darstellug der Chi2-Funktion\n", "\n", "import numpy as np, matplotlib.pyplot as plt\n", "from scipy.stats import chi2\n", "\n", "s = np.linspace(0, 25 ,100)\n", "nf = 0\n", "plt.figure(figsize=(10.,7.5))\n", "plt.xlabel('S')\n", "plt.ylabel('chi2 pdf')\n", "for nf in range(0,10):\n", " plt.plot(s, chi2.pdf(s,nf),'-')\n", "plt.show()\n", "```\n", "b) Stellen Sie auch die $\\chi^2$-Wahrscheinlichkeit grafisch dar. Eine besonders\n", "aussagekräfige Darstellung erhalten Sie, wenn Sie die $\\chi^2$-Wahrscheinlichkeit gegen $S/n_f$\n", "auftragen.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Code hier eingeben\n", "# -->" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "### 3.4 Verfahren und Programmpakete zur Parameteranpassung \n", "---\n", "\n", "Wir hatten in [Abschnitt 3.1](#Sec-Funktionsanpassung) gesehen, dass durch Minimierung der \n", "Residuensumme und Bestimmung der zweiten Ableitungen Schätzungen der Parameter von angepassten \n", "Modellen und die Bestimmung der Parameterunsicherheiten möglich sind.\n", "\n", "Für einfache Modelle, insbesondere solche, bei denen das Modell linear in den Parametern ist, kann die Minimierung mit Methoden der Analysis ausgeführt werden, d. h. man erhält als Ergebnis fertige Formeln, mit denen die besten Schätzwerte der Parameter und ihrer Unsicherheiten berechnet werden können. \n", "\n", "Für nichtlineare Modelle oder unter Annahme komplexerer Unsicherheiten - z. B. auch Unsicherheiten in Richtung\n", "der Abszisse (\"x-Fehler\") statt lediglich in der Ordinatenrichtung, muss die Minimierung mit numerischen Methoden ausgeführt werden. \n", "\n", "Für beide Vorgehensweisen hatten wir schon einfache Beispiele diskutiert. Die allgemeine Lösung für\n", "lineare Modelle wird in [Abschnitt 3.4.1](#SubSec-LineareModelle) vorgestellt.\n", "In [Abschnitt 3.4.2](#SubSec-Programmpakete) werden dann Programmpakete für numerischer\n", "Verfahren beschrieben.\n", "\n", "Zunächst soll aber noch eine Erweiterung der Residuensumme besprochen werden, die es ermöglicht, auch Anpassungen an korrelierte Messwerte vorzunehmen. Dazu schreiben wir die Residuensumme ([Gleichung (3.1)](#Equ-3.1)) in vektorieller Form, mit dem Vektor $\\vec{y}$ der Messwerte sowie dem Vektor $\\vec{x}$ der Messpunkte und dem Parametervektor $\\vec{p}$. Die Unsicherheiten werden repräsentiert durch eine Kovarianzmatrix ${\\bf V}$ mit den Elementen $\\sigma_i^2$ in der Diagonalen. Damit ergibt sich die allgemeine Darstellung der Residuensumme\n", "\n", "> $S(\\vec{y}, \\vec{p}) = \\left( \\vec{y} - \\vec{f}(\\vec{x},\\vec{p} )\\right)^T \n", " \\, {\\bf V^{-1}} \n", " \\left( \\vec{y} - \\vec{f}(\\vec{x},\\vec{p} )\\right)\\, $ \n", " Gleichung (3.3) \n", "\n", "Dieses zunächst rein formale Umschreiben ist tatsächlich auch die korrekte Darstellung, wenn die Kovarianzmatrix nichtverschwindende Nebendiagonalelemente enthält. Begründbar ist dies durch die Betrachtung der Gaußverteilung für mehrere, auch korrelierte Zufallsvariablen, bei der genau dieser Term im im Argument der Exponentialfunktion auftritt. \n" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "#### 3.4.1 Analytische Lösung zur Anpassung linearer Modelle \n", "---\n", "\n", "\n", "Im Rahmen dieses Hands-On-Tutoriums ist nur eine verkürzte Darstellung des theoretischen Hintergrunds \n", "möglich. Eine detailliertere Diskussion mit ausführlicher Herleitung findet sich im Skript\n", "[Funktionsanpassung mit der $\\chi^2$-Methode](http://www.etp.kit.edu/~quast/Skripte/Chi2Method.pdf)\n", "\n", "Wir betrachten Modelle für die Messwerte der Form\n", "${f_i}=\\sum_{j =1} ^k g_j(x_i)\\cdot p_j \\,, \\,\\, {\\small i=1,\\ldots,N}$.\n", "Dies lässt sich kompakt in Vektor- und Matrixschreibweise\n", "darstellen: \n", " > $A_{ij} := g_j(x_i) \\, \\Rightarrow \\, \\vec{f}( \\vec{x}, \\vec{p} ) = {\\bf A}(\\vec{x})\\vec{p}$\n", " \n", "Unter Verwendung von [Gleichung (3.3)](#Equ-3.3) ergibt sich folgende Darstellung für die Residuensumme: \n", " \n", "> $S(\\vec{p})=(\\vec{y} - {\\bf A} \\vec{p})^T \\, {\\bf V}^{-1} (\\vec{y} - {\\bf A} \\vec{p})$\n", "\n", "Die Minimierungsaufgabe sieht dann wie folgt aus: \n", "\n", "> $\\nabla_p S(\\vec{p})=0$, bzw. etwas ausführlicher \n", "\n", "> $\\nabla_p S(\\vec{p})=\\nabla_p \\left( (\\vec{y} - {\\bf A} \\vec{p})^T \\, {\\bf V}^{-1} (\\vec{y} - {\\bf A}\\vec{p})\\right) = 0$\n", "\n", "Das Ausführen der Ableitung führt auf eine Matrixgleichung für den Parametervektor $\\vec{p}$, die sich \n", "mit etwas linearer Algebra (und recht viel Schreibarbeit) analytisch lösen lässt. Für den Vektor der besten Parameterschätzungen $\\hat{\\vec{p}}$ ergibt sich der kompliziert anmutende Matrixausdruck \n", "\n", "> $\\hat{\\vec{p}} = ({\\bf A}^T {\\bf V}^{-1}{\\bf A})^{-1} \\, {\\bf A}^T {\\bf V}^{-1}\\, \\vec{y}\\,$ \n", " Gleichung (3.4a) \n", "\n", "Letzten Endes handelt es sich im Ergebnis allerdings lediglich um eine Linearkombination der Messwerte $y_i$.\n", "\n", "Die Kovarianzmatrix der Parameter ergibt sich zu\n", "\n", "> ${\\bf V}_{\\hat{p}}=\\left ({\\bf A}^T {\\bf V}^{-1} {\\bf A}\\right)^{-1}\\,$\n", " Gleichung (3.4b) \n", "\n", "\n", "#### Spezialfall 1: Lineare Regression für unabhängige Messungen\n", "\n", "Eine sehr häufig vorkommende Anwendung der obigen Formeln ist die einfache Geradenanpassung (\"lineare Regression\") an unkorrelierte Messdaten. In diesem Fall sind die Matrizen ${\\bf A}$ und ${\\bf V}$ gegeben durch:\n", "\n", "${\\bf A} \\,=\\, \\left( \\begin{array}{cc} \n", " 1 & x_1 \\\\\n", " \\ldots & \\ldots \\\\\n", " 1 & x_N \\\\ \n", " \\end{array} \\right)\\,, \n", "\\quad$\n", "$ {\\bf V}^{-1} \\,=\\, \\left( \\begin{array}{ccccc} \n", " \\frac{1}{\\sigma_1^2} & 0 & \\ldots & 0 & 0 \\\\\n", " 0 & \\ldots & \\frac{1}{\\sigma_i^2} & \\ldots & 0 \\\\ \n", " 0 & 0 & \\ldots & 0 & \\frac{1}{\\sigma_N^2} \\\\\n", " \\end{array} \\right)\n", "$\n", "\n", "Durch Einsetzen in die allgemeine Lösung oben erhält man mit der Einführung einiger Abkürzungen \n", "\n", "$$\\begin{array}{lclcl}\n", "S_{1~} =~\\sum_{i=1}^N \\frac{1}{\\sigma_i^2}\\,, & &\n", "S_{x~} =~ \\sum_{i=1}^N \\frac{x_i}{\\sigma_i^2} \\,,& &\n", "S_{y~} =~ \\sum_{i=1}^N \\frac{y_i}{\\sigma_i^2} \\,,\\\\[0.3pc]\n", "S_{xx} =~ \\sum_{i=1}^N \\frac{x_i^2}{\\sigma_i^2} \\,,& &\n", "S_{xy} =~ \\sum_{i=1}^N \\frac{x_i\\,y_i}{\\sigma_i^2} \\,, & &\n", "D_{~~} = {S_1 S_{xx}-S_x^2}\n", "\\end{array}$$\n", "\n", "die spezielle Lösung \n", "\n", "$$\\begin{array}{ccc c c c c }\n", "\\hat p_1&=&\\frac{S_{xx}S_y-S_xS_{xy}}{D}\\,, & & \n", "{\\sigma_{p_1}}^2=\\frac{S_{xx}}{D}\\,, \\\\[0.3pc]\n", "\\hat p_2&=&\\frac{S_1S_{xy}-S_xS_y}{D}\\,, & &\n", "{\\sigma_{p_2}}^2=\\frac{S_{1}}{D}\\,, & & \n", "V_{12}=\\frac{-S_{x}}{D}\\,. \\\\\n", "\\end{array}$$\n", " \n", "Sicher verstehen Sie nun, warum Anpassungen in der Vergangenheit als schwierig galten.\n", "Über lange Zeit waren die hier gezeigten Formeln die Basis von Modellanpassungen, und\n", "viele historische Methoden beruhten auf einer \"Linearisierung\" des Problems, also\n", "der Anwendung von Funktionen auf die Daten, um sie auf eine Geradenform zu transformieren.\n", "Ein Beispiel dafür ist die logarithmische Darstellung von exponentiellen Zusammenhängen.\n", "Dabei wurde allerdings außer acht gelassen, dass nicht nur die Daten selbst, sondern auch\n", "ihre Unsicherheiten transformiert werden müssten. Nach der Transformation sind die\n", "Unsicherheiten im allgemeinen nicht mehr gaußförmig, und die Voraussetzungen zur \n", "Anwendung der Methode der kleinsten Quadrate sind eigentlich nicht mehr erfüllt. \n", "Wurde aber trotzdem so gemacht ...\n", "\n", "Sie könnten nun versuchen, die obigen Formeln selbst zu implementieren. Sollten Sie eine\n", "einfache lineare Regression benötigen, können Sie auf die Implementierung in \n", "*PhyPraKit.linRegression(x, y, sy)* zurückgreifen. Oder Sie lesen weiter und\n", "nutzen eines der im nächsten Unterkapitel beschriebenen, auf numerischen Methoden\n", "basierenden Programmpakete. \n", "\n", "Zum **Üben** hier ein Beispiel einer Geradenanpassung mit *PhyPraKit.linRegression()*.\n", "Schauen Sie sich zunächst den Code der Funktion *linRegression* an; mit *numpy* ist\n", "die Implementierung doch gar nicht so schwer !?\n", "\n", "```\n", "def linRegression(x, y, sy):\n", " \"\"\"\n", " linear regression y(x) = ax + b \n", "\n", " method: \n", " analytical formula\n", "\n", " Args:\n", " * x: np-array, independent data\n", " * y: np-array, dependent data\n", " * sy: scalar or np-array, uncertainty on y\n", "\n", " Returns:\n", " * float: a slope\n", " * float: b constant\n", " * float: sa sigma on slope\n", " * float: sb sigma on constant\n", " * float: cor correlation\n", " * float: chi2 \\chi-square\n", " \"\"\" \n", " \n", " # calculate auxilary quantities\n", " S1 = sum(1./sy**2)\n", " Sx = sum(x/sy**2)\n", " Sy = sum(y/sy**2)\n", " Sxx = sum(x**2/sy**2)\n", " Sxy = sum(x*y/sy**2)\n", " D = S1*Sxx-Sx**2\n", "\n", " # calculate results:\n", " a = (S1*Sxy-Sx*Sy)/D # slope\n", " b = (Sxx*Sy-Sx*Sxy)/D # constant\n", " sa = np.sqrt(S1/D)\n", " sb = np.sqrt(Sxx/D)\n", " cov = -Sx/D\n", " cor = cov/(sa*sb)\n", " chi2 = sum(((y-(a*x+b))/sy)**2)\n", "\n", " return a, b, sa, sb, cor, chi2\n", "```\n", "\n", "Und hier der Code zum Ausprobieren; das ist schon ein relativ komplettes \n", "Beispiel mit Darstellung der Daten und passenden Beschriftungen.\n", "\n", "```\n", "# Beispiel zur Anwendung von linRegression(x, y, sy)\n", "import numpy as np, matplotlib.pyplot as plt\n", "\n", "# die Daten\n", "x = np.array([ 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.])\n", "y = np.array([1.46, 1.4, 1.87, 2.10, 2.60, 2.87, 3.03 , 3.93, 4.04, 4.08])\n", "sy = np.array([0.07, 0.07 , 0.09, 0.10, 0.12,0.13, 0.15, 0.17, 0.18, 0.20])\n", "\n", "# Fit ausführen und Ergebnis speichern\n", "a, b, sa, sb, cor, chi2 = linRegression(x, y, sy)\n", "\n", "# Daten und Modell anzeigen \n", "plt.errorbar(x, y, yerr=sy, fmt='x',label='Übungsdaten')\n", "plt.plot(x, a*x+b, '-', label='angepasste Gerade')\n", "plt.xlabel('x')\n", "plt.ylabel('measurements y')\n", "plt.legend()\n", "plt.show()\n", "```\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Code von oben hier eingeben:\n", "# -->" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Spezialfall 2: Anpassung einer Konstanten: Mittelwertbildung**\n", "\n", "Die Mittelwertbildung entspricht der Anpassung einer konstanten Funktion $y_i = \\bar{m}$. \n", "In diesem Fall wird die Matrix ${\\bf A}$ in [Gleichung (3.4a)](#Equ-3.4a) zu einem Vektor: \n", "${\\bf A} = {\\small (1, \\ldots , 1)^T }$, und man erhält folgendes Ergebnis:\n", "\n", "> $ \\bar{m} = \\frac{1}{\\sum_{i,j} \\left( {\\bf V}^{-1} \\right)_{ij}}\n", " {\\sum_{i,j} \\left( {\\bf V}^{-1} \\right)_{ij}x_j}\\,$\n", " Gleichung (3.5a) \n", " \n", "Für die Varianz ergibt sich:\n", "> ${\\bf V}_\\bar{m} = \\frac{1}{\\sum_{i,j} \\left( {\\bf V}^{-1} \\right)_{ij}}\\,$\n", " Gleichung (3.5b) \n", " \n", " Glleichungen (3.5a,b) sind die Verallgemeinerung für korrelierte Messungen des schon in \n", " [Abschnitt 3.1.2](#SubSec-AnalytischeLoesung) hergeleiteten Ergebnisses für unabhängige \n", " Messungen. Man kann es sich leicht merken:\n", " in beiden Fällen ist der Mittelwert eine gewichtete Summe, \n", " $\\bar{m}= \\left( \\sum w_i \\right)^{-1} \\sum w_i x_i $, wobei die Gewichte die \n", " Elemente der Inversen der Kovarianzmatrix sind. \n", " \n", " Eine Implementierung findet sich in *PhyPraKit.wmean()*. \n", " \n", "Probieren Sie folgendes **Übungsbeispiel** aus:\n", "```\n", "import numpy as np, PhyPraKit as ppk\n", "\n", "# Daten und Unsicherheiten\n", "m = np.array([0.82, 0.81, 1.32, 1.44, 0.93, 0.99])\n", "sig_u = 0.1 # unabhängig\n", "sig_s = 0.15 # korreliert für Messungen (1,2), (3,4) und (5,6)\n", "sig_t = 0.05 # korreliert für alle Messungen\n", "\n", "# Konstruktion der Komponenten für die Kovarianzmatrix\n", "sys_x= [[sig_t, sig_t, sig_t, sig_t, sig_t, sig_t],\n", " [sig_s, sig_s, 0. ,0. ,0. ,0.], \n", " [0., 0., sig_s, sig_s, 0., 0.], \n", " [0., 0., 0., 0., sig_s, sig_s] ]\n", "stat_x = sig_u * np.ones(6)\n", "\n", "# Zusammenbauen der Kovarianzmatrix incl. der unabhängigen Unsicherheiten\n", "V = ppk.BuildCovarianceMatrix(stat_x, sys_x)\n", "#print(V)\n", "\n", "# Berechnung des Mittelwerts (ohne unabhängige Unsicherheiten, weil schon in V enthalten)\n", "mean, sigm = ppk.wmean(m, None, V, pr=True)\n", "```\n", "Geben Sie das Beispiel in die Code-Zeile unten ein. Verändern Sie dann einige\n", "der Komponenten der Unsicherheit. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Code hier eingeben\n", "# -->" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**_Anmerkung_**: Sie habe vielleicht die Daten im Beispiel wiedererkannt? \n", "Es sind die gleichen, die wir ausführlich in [Abschnitt 2.4.6](#SubSec-BeispielKovarianzmatrix) \n", "besprochen hatten. Die hier gezeigte Vorgehensweise ist im allgemeinen Fall das korrekte Verfahren zur Mittelwertbildung. Wären nämlich die Unsicherheiten der einzelnen Messwerte nicht gleich\n", "groß, würde sich ein anderer Mittelwert ergeben - siehe zum Vergleich [Abschnitt 3.1.2](#SubSec-AnalytischeLoesung) \n", "Mittelwertbildung bei verschiedenen Unsicherheiten. Die einfache Herangehensweise \n", "wie in [Abschnitt 2.4.6](#SubSec-BeispielKovarianzmatrix) wäre dann falsch gewesen.\n", "\n", "__Also__:\n", "> Zur Mittelwertbildung korrelierter Werte eine Anpassung einer Konstanten\n", " unter Berücksichtigung der Kovarianz-Matrix vornehmen!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "#### 3.4.2 Numerische Anpassung und Programmpakete\n", "---\n", "\n", "\n", "Mit der allgemeinen Verfügbarkeit von Computern werden heute Anpassungen mit Methoden der numerischen Optimierung ausgeführt. Die Algorithmen sind schnell genug und funktionieren auch bei nichtlinearen Problemstellungen. Am KIT wurden in den letzten Jahren *Python*-Programmpakete entwickelt, die eine korrekte Anpassung auch für kompliziert korrelierte Datenpunkte erlauben:\n", "- *PhyPraKit* (http://ekpwww.etp.kit.edu/~quast/PhyPraKit/htmldoc/): Sammlung von Funktionen \n", " zum Einlesen von Daten aus diversen Quellen, zur Datenvisualisierung, Signalbearbeitung und\n", " zur statistischen Datenauswertung und Modellanpassung sowie Werkzeuge zur Erzeugung simulierter Daten.\n", "- Das Paket *PhyPraKit* enthält ein schlank gehaltenes Interface zum sehr vielseitig verwendbaren \n", " Optimierungs- und Fehleranalyse-Werkzeug *iminuit* (https://iminuit.readthedocs.io/en/stable/). \n", "- *kafe2* (http://ekpwww.etp.kit.edu/~quast/kafe2/htmldoc/ und [Abschnitt 3.5](#Sec-Kafe)): \n", " Softwareumgebung zur Anpassung von Modellen an Daten, speziell entwickelt für das physikalische Praktikum.\n", "\n", "Ein einfaches Beispiel hatten wir als Übung zu [Abschnitt 3.1.1](#SubSec-BeispielFunktionsanpassung) \n", "schon kennen gelernt: das Abstandsmaß für eine große, recht dicht liegende Punktemenge auftragen\n", "und das Minimum suchen. Dieses Verfahren lässt sich sehr leicht viel effizienter gestalten. \n", "Man startet mit weniger Punkten, z.B. minimal drei, und wählt dann zusätzliche Punkte, die \n", "das Minimum sukzessive immer besser eingrenzen. Moderne Optimierer zur Suche des Minimums einer\n", "skalaren Funktion setzen auf eine Vielzahl von Algorithmen, deren detaillierte Beschreibung hier\n", "zu weit führen würde.\n", "\n", "\n", "**Behandlung von Unsicherheiten bezüglich der Abszisse (\"x-Fehler\")**\n", "\n", "Unsicherheiten der Datenpunkte $(x_i,y_i)$ bezüglich der Abszissenachse, d. h. \n", "Unsicherheiten auf die $x_i$, können mit der $\\chi^2$-Methode recht elegant durch\n", "Iteration berücksichtigt werden: \n", " - zunächst erfolgt eine Anpassung ohne Unsicherheiten auf die Abszissenwerte,\n", " - im zweiten Schritt werden diese dann mit Hilfe der ersten Ableitungen \n", " $f'(x_i)$ der im ersten Schritt angepassten Funktion $f$ in entsprechende Unsicherheiten \n", " der Ordinate umgerechnet und quadratisch zu den betreffenden Unsicherheiten der Ordinate addiert: \n", " ${\\sigma_i}^2={\\sigma_y}_i^2 \\,+\\,(f'(x_i) \\cdot {\\sigma_x}_i)^2\\,$.\n", " Die Größe $\\chi^2$ wird jetzt mit den neuen Unsicherheiten $\\sigma_i$ berechnet\n", " und minimiert. \n", " - Ein dritter Schritt, der der Vorgehensweise beim zweiten Schritt entspricht, dient zur \n", " Verbesserung des Ergebnisses und zur Fehlerkontrolle - der Wert von $\\chi^2$ am Minimum\n", " darf sich vom zweiten zum dritten Schritt nicht signifikant ändern, ansonsten muss\n", " nochmals iteriert werden.\n", "Bei der Anwendung dieses Verfahrens ist es essentiell, dass die Modellfunktion über den\n", "Bereich der Unsicherheiten gut durch eine Gerade approximiert wird. \n", "\n", "Übrigens: Dank hoher Rechenleistung und effizienter numerischer Algorithmen ist es heute auch \n", "kein Problem mehr, die Neuberechung der Unsicherheiten in jedem Schritt der numerischen Optimierung\n", "durchzuführen. Dazu muss das Programmierinterface des verwendeten Optimierers aber den vollen Zugriff\n", "auf die zu minimierende \"Kostenfuktion\" erlauben. Wenn dies nicht der Fall ist, funktioniert das \n", "oben skizzierte Verfahren durch mehrfaches Anwenden mit beliebigen Anpassungsprogrammen. \n", "\n", "Im Paket *scipy.optimize* gibt es die Methode *odr* (\"orthogonal distance regression\"), die\n", "ein anderes Verfahren nutzt. *odr* optimiert den mit den Unsicherheiten gewichteten Abstand \n", "der Datenpunkte von der anzupassenden Kurve im zweidimensionalen Raum.\n", "\n", "\n", "**Funktionsanpassung mittels numerischer Verfahren** \n", "\n", "Aus Anwendersicht ändert sich im Vergleich zum Beispiel aus \n", "[Abschnitt 3.4.1](#SubSec-LineareModelle) nichts wesentliches, wenn man numerische\n", "Verfahren auf dem Computer nutzt.\n", "Statt der Funktion *linRegression()* können Sie einfach die Funktion *odFit()* aus dem \n", "Paket *PhyPraKit* aufrufen. Diese Funktion verwendet die Methoden *scipy.curve_fit* und\n", "*scipy.optimize.odr*. aus dem Paket *scipy*. Mit deren Hilfe ist es möglich, beliebige\n", "in *python* definierte Funktionen an Daten mit Unsicherheiten in Ordinaten- und Abszissenrichtung\n", "anzupassen. Wenn keine Unsicherheiten in Abszissenrichtung benötigt werden, werden sie auf Null gesetzt.\n", "\n", "Ein zur Übung aus [Abschnitt 3.4.1](#SubSec-LineareModelle) analoges Beispiel sieht \n", "also ganz ähnlich aus wie dort; es gibt nur kleinere Anpassungen an das etwas andere \n", "Interface von *odFit* im Vergleich zu *linRegression*:\n", "```\n", "# Beispiel zur Anwendung von PhyPraKit.odFit(f, x, y, sx, sy)\n", "import numpy as np, matplotlib.pyplot as plt\n", "from PhyPraKit import odFit\n", "\n", "# das Modell\n", "def fitf(x, a, b): \n", " ''' linear model ''' \n", " return a*x + b\n", "\n", "# die Daten\n", "x = np.array([ 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.])\n", "y = np.array([1.46, 1.4, 1.87, 2.10, 2.60, 2.87, 3.03 , 3.93, 4.04, 4.08])\n", "sx = np.zeros(10)\n", "sy = np.array([0.07, 0.07 , 0.09, 0.10, 0.12,0.13, 0.15, 0.17, 0.18, 0.20])\n", "\n", "# Fit ausführen und Ergebnis speichern\n", "par, pare, cor, chi2 = odFit(fitf, x, y, sx, sy)\n", "a = par[0]\n", "b = par[1]\n", "# Daten und Modell anzeigen \n", "plt.errorbar(x, y, yerr=sy, fmt='x',label='Übungsdaten')\n", "plt.plot(x, a*x+b, '-', label='angepassete Gerade')\n", "plt.xlabel('x')\n", "plt.ylabel('measurements y')\n", "plt.legend()\n", "plt.show()\n", "```\n", "\n", "Geben Sie als **Übung** den Code in die Zelle unten ein und vergleichen Sie das Ergebnis mit \n", "dem Beispiel von *linRegression*. Geben Sie dazu jeweils auch die Parameterwerte und Unsicherheiten\n", "mit Hilfe eines *print()*-Befehls aus und prüfen Sie, ob beide Methoden im Rahmen der Unsicherheiten\n", "gut übereinstimmen !" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Code von oben hier eingeben\n", "# -->" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Das oben beschriebene Programmpaket reicht für die meisten Belange völlig aus. Allerdings \n", "kann *odFit* nicht mit korrelierten Unsicherheiten umgehen. Die Programmpakete *kafe* \n", "(in der aktuellen Version *kafe2*) oder PhyPraKit/phyFit* unterstützen dagegen korrelierte\n", "Unsicherheiten entlang der Ordinaten- und Abszissen-Richtung und bieten Hilfsfunktionen zur\n", "Konstruktion der Kovarianzmatrix. Die Pakete enthalten ebenfalls eigene grafische Ausgaben,\n", "die neben der Modellfunktion auch die mittels Fehlerfortpflanzung ermittelte Unsicherheit \n", "der angepassten Modellfunktion als Band anzeigt. \n", "In *PhyPraKit* gibt es ein einfach gehaltene Interfaces zur komfortablen Anwendung der \n", "Funktionen dieser Pakete, siehe *k2Fit(f, x, y, sx, sy, ...)* oder *xyFit(f, x, y, sx, sy, ...)*. \n", "\n", "**Übung** Für die gleiche Anpassung wie eben mit *odFit* müssen Sie nur im entsprechenden \n", "Beispiel *odFit* durch *k2Fit* oder *xyFit* ersetzen; den Code zur Anzeige der Grafik können\n", "Sie entfernen, weil diese Fuktionen in der Grundeinstellung eine eigene grafische Ausgabe erzeugen. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Code für k2Fit odf xyFit hier eingeben\n", "# -->" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In *k2Fit* bzw. *xyFit* gibt es zusätzlich zu den aus *odFit* bekannten Parametern \n", "(*func*, *x*, *y*, *sx* und *sy*) eine ganze Reihe an optionalen Parametern, die \n", "viele zusätzliche Möglichkeiten bieten. Die Rückgabeliste als *numpy-array*s oder \n", "*python-dictionary* enthält Informationen, die denen von *odFit* entsrpechen, also\n", "die angepassen Parameterwerte und deren Unsicherheiten, die Matrix der \n", "Parameterkorrelationen und die Qualität der Anpassung.\n", "\n", "Gehen wir die optionalen Argumente und deren Bedeutung der Reihe nach durch:\n", "\n", " * _xabscor_: absolute, correlated error(s) on x, absolute korrelierte Unsicherheiten x-Richtung\n", " * _yabscor_: absolute, correlated error(s) on y, absolute korrelierte Unsicherheiten y-Richtung\n", " * _xrelcor_: relative, correlated error(s) on x, relative korrelierte Unsicherheiten x-Richtung\n", " * _yrelcor_: relative, correlated error(s) on y, relative korrelierte Unsicherheiten y-Richtung\n", " \n", "Optionen zur Kontrolle der Anpassung \n", "\n", " * _ref_to_model_: wenn *True*, werden relative Fehler werden auf den Modellwert, nicht auf den\n", " Messwert bezogen. Verzerrungen der Parameterschätzung durch Abwärtsfluktuationen der Daten \n", " werden so vermieden. \n", " * _p0: array-like_, Startwerte für die Parameteranpassung, die bei nicht-linearen Problemstellungen\n", " mit nicht eindeutigen Minima benötigt werden\n", " * _constrains_: _(name, value, uncertainty)_, externe Einschränkungen an Parameter innerhalb bekannter\n", " Unsicherheiten\n", " * _limits_: _(name, min, max)_, Grenzen für die Parameter, die benötigt werden, um unphysikalische oder\n", " mathemematisch nicht definierte Bereiche, wie den Logarithmus oder die Wurzel aus negativen Zahlen, \n", " bei der Parameteroptimierung auszuschließen\n", " \n", "Optionen zur Steuerung der grafischen Ausgabe in *k2Fit*\n", "\n", " * _plot: boolean flag_, Grafische Ausgabe erzeugen\n", " * _axis_labels: list of strings_, Achsenbeschriftungen x und y\n", " * _data_legend: string_, Name für angezeigte Daten in Legende\n", " * _model_name: string_, Name für angezeigte Modell in der Info-Box \n", " * _model_expression: string_, LaTeX-Ausdruck für die angepasste Funktion\n", " * _model_legend: string_, Name für Modell in Legende \n", " * _model_band: string_, Name für das 1-sigma Band in Legende\n", " * _fit info: boolean_, Info-Box mit Fitergebnisse ein/ausschalten \n", " * _asym_parerrs_: Ausgabe von asymmetrischen Unsicherheiten, d.h. eines 68%-Konfidenz Intervalls\n", " * _plot_cor_: Anzeige von Konfidenzkonturen der Parameter\n", " * _showplots_: Anzeige aller Grafiken auf dem Bildschirm \n", " * _quiet_: Unterdrückung von Textausgaben.\n", "\n", "\n", "\n", "Eine sehr **wichtige Anwendung von _kafe2_** oder **_phyFit_** ist die Behandlung von \n", "korrelierten systematischen Unsicherheiten in der Anpassung. Mit anderen numerischen \n", "Werkzeugen muss zunächst ein Anpassung ohne Berücksichtigung der systematischen \n", "Unsicherheiten durchgeführt und dann deren Effekt mittels Fehlerfortpflanzung auf \n", "die Parameter propagiert werden.\n", "Für absolute, also für alle Messdaten gleiche Unsicherheiten funktioniert das sehr gut. \n", "Bei relativen Unsicherheiten, also z. B. \"$\\pm1\\%$ vom Messwert\", oder in Fällen mit \n", "komplizierter Struktur der Korrelationen gibt es aber keine andere Wahl als sie in der\n", "Anpassung mit zu berücksichtigen, weil solche systematischen Effekte auch die Ergebnisse \n", "für die Parameterwerte ändern. \n", "\n", "Ein Beispiel soll das illustrieren. Hier ein minimales Code-Beispiel zur Anpassung mit *kafe2*\n", "```\n", "# Beispiel zur Anwendung von PhyPraKit.k2Fit(f, x, y, sx, sy)\n", "import numpy as np, matplotlib.pyplot as plt\n", "from PhyPraKit import k2Fit\n", "\n", "# das Modell\n", "def fitf(x, a, b): \n", " ''' linear model ''' \n", " return a*x + b\n", "\n", "# die Daten\n", "x = np.array([ 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.])\n", "y = np.array([1.46, 1.4, 1.87, 2.10, 2.60, 2.87, 3.03 , 3.93, 4.04, 4.08])\n", "sx = np.zeros(10)\n", "sy = np.array([0.07, 0.07 , 0.09, 0.10, 0.12,0.13, 0.15, 0.17, 0.18, 0.20])\n", "\n", "# optionale absolute und relative Unsicherheit \n", "abs_syst = 0.1 # absolute systematic error \n", "rel_syst = 0.1 # 10% relative systematische Unsicherheit\n", "\n", "# Fit ausführen und Ergebnis speichern\n", "par, pare, cor, chi2 = k2Fit(fitf, x, y, sx, sy,\n", " yabscor = 0., \n", " yrelcor= 0.\n", " )\n", "```\n", "Als **Übung** soll zunächst genau dieser Code in der Code-Zelle unten ausgeführt werden.\n", "Bitte das Ergebnis in der Text-Zelle festhalten.\n", "\n", "Jetzt wird eine abosulte systematische Unsicherheit hinzugefügt: *yabscor = abs_syst*.\n", "Bitte den Code modifzieren, ausführen und das Ergebnis sichern. \n", "\n", "Nun die absolute Unsicherheit wieder auf 0. setzen und eine relative systematische Unsicherheit hinzufügen: *yrelcor=rel_syst*, den modifizierten Code ausführen und das Ergebis eintragen. \n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Code k2fit mit Systematik hier eingeben\n", "# -->" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Ergebnis**\n", " 1. ohne systematische Unsicherheit: a = +/- , b = +/- \n", " 2. mit absoluter systematischer Unsicherheit: a = +/- , b = +/- \n", " 3. mit relativer systematischer Unsicherheit: a = +/- , b = +/- \n", "\n", "Für den Fall 2. bitte mit Hilfe des Fehlerfortpflanzungsgesetzes prüfen, \n", " ob der Unterschied zu 1. den Erwartungen entspricht.\n", " \n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Als **weiteres Beispiel** greifen wir noch einmal auf die Daten aus \n", "[Abschnitt 2.4.6](#SubSec-BeispielKovarianzmatrix) zurück und\n", "illustrieren die Mittelwertbildung korrelierter Messwerte mit *kafe2*.\n", "Dank der Fähigkeit von *kafe2*, mit komplizierten Kovarianzmatrizen umgehen zu können, \n", "lässt sich die Aufgabenstellung auch komfortabel mit *k2Fit* lösen. \n", "Hier die notwendigen Code-Teile:\n", " \n", "```\n", "# Mittelung korrelierter Messwerte mit kafe2 / k2Fit\n", "\n", "import numpy as np, matplotlib.pyplot as plt\n", "from PhyPraKit import k2Fit\n", "\n", "# das Modell\n", "def fitf(x, c):\n", " # das Ergebnis muss ein Vektor der Länge von x sein \n", " return c*np.ones(len(x)) \n", "\n", "# Daten und Unsicherheiten\n", "x = np.arange(6) + 1.\n", "m = np.array([0.82, 0.81, 1.32, 1.44, 0.93, 0.99])\n", "sig_u = 0.1 # unabhängig\n", "sig_s = 0.15 # korreliert für Messungen (1,2), (3,4) und (5,6)\n", "sig_t = 0.05 # korreliert für alle Messungen\n", "\n", "# Konstruktion der Komponenten für die Kovarianz-Matrix\n", "sys_y= [[sig_t, sig_t, sig_t, sig_t, sig_t, sig_t],\n", " [sig_s, sig_s, 0. ,0. ,0. ,0.], \n", " [0., 0., sig_s, sig_s, 0., 0.], \n", " [0., 0., 0., 0., sig_s, sig_s] ]\n", "\n", "# Fit ausführen\n", "par, pare, cor, chi2 = k2Fit(fitf, x , m, sx=0., sy=sig_u, yabscor = sys_y)\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Code für Mittelung mit k2Fit() hier eingeben\n", "# -->" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ein häufiges **Beispiel aus der Praxis** stellt die Messung von Strom-Spannungskennlinien dar. \n", "Sowohl die Strom- als auch die auf der Abzisse aufgetragene Spannungsmessung sind von verschiedenen Unsicherheiten betroffen. So gibt es auch bei digitalen Messinstrumenten eine Anzeigeunsicherheit, \n", "die im Datenblatt angegeben ist. Dazu kommt ggf. noch Signalrauschen, dass man aus der Stabilität \n", "der Anzeige bei wiederholten Messungen abschätzt. Diese Unsicherheiten sind für jeden Messwert\n", "unabhängig voneinander. \n", "Zusätzlich gibt es noch eine Kalibrationsunsicherheit, die der Hersteller üblicherweise als\n", "relativen Wert bezogen auf den wahren Messwert in Datenblatt angibt. Die Kalibrationsunsicherheit\n", "ist eine gemeinsame, relative Unsicherheit aller Messwerte. Je nach Kalibrationsverfahen ist\n", "es möglich, dass die Unsicherheiten für in verschiedenen Messbereichen aufgenommene Messungen\n", "jeweils unabhängig sind. \n", "\n", "Hier einige typische Daten für Messgenauigkeiten:\n", "\n", " - **Spannungsmessung**: Anzeige 4000 Counts, Genauigkeit 0.5% $\\pm$ 3 Digits, \n", " Messbereich 2 V, Signalrauschen 0.005 V\n", " - **Strommessung**: Anzeige 2000 Counts, Genauigkeit 1.0% $\\pm$ 3 Digits, \n", " Messbereiche 200 µA, 20 mA und 200 mA, Signalrauschen 0.025 mA\n", "\n", "Als Anwendungsfall betrachten wir Messungen mit diesen Messinstrumenten an einer \n", "Halbleiterdiode: \n", "```\n", "# Messdaten:\n", "# - Spannung im Messbereich 2V\n", " data_x = [0.450, 0.470, 0.490, 0.510, 0.530, \n", " 0.550, 0.560, 0.570, 0.580, 0.590, 0.600, 0.610, 0.620, 0.630,\n", " 0.640, 0.645, 0.650, 0.655, 0.660, 0.665 ]\n", "# - Strommessungen: 2 im Bereich 200µA, 12 mit 20mA und 6 mit 200mA\n", " data_y = [0.056, 0.198, 0.284, 0.404, 0.739, 1.739, 1.962,\n", " 2.849, 3.265, 5.706, 6.474, 7.866, 11.44, 18.98,\n", " 23.35, 27.96, 38.61, 46.73, 49.78, 57.75]\n", "``` \n", "\n", "Angepasst werden soll die Shokley-Gleichung:\n", "\n", "```\n", "def Shockley(U, I_s = 0.5, U0 = 0.03):\n", " \"\"\"Parametrisierung einer Diodenkennlinie\n", "\n", " U/U0 sollte während der Anpassung auf einen Wert <150 \n", " beschränkt werden, um Überscheitungen des mit 64 Bit \n", " Genauigkeit darstellbaren Zahlenraums zu verhindern\n", "\n", " Args:\n", " - U: Spannung (V)\n", " - I_s: Sperrstrom (nA)\n", " - U0: Temperaturspannung (V) * Emissionskoeffizient\n", " \n", " Returns:\n", " - float I: Strom (mA)\n", " \"\"\"\n", " return 1e-6 * I_s * np.exp( (U/U0) - 1.)\n", "```\n", "\n", "Mit Ihrem bisherigen Wissen können Sie dieses Problem recht komfortabel mit der oben beschriebenen \n", "Funktion *k2Fit* oder mit der äquivalenten Funktion *xyFit* lösen, die das Modul *pyhFit* aus\n", "dem Paket *PhyPraKit* nutzt.\n", "\n", "\n", "**Hilfestellung**\n", "\n", "Code zur Berechnung der Unsicherheiten:\n", "```\n", "# Komponenten der Messunsicherheit\n", "# - Genauigkeit Spannungsmessung: 4000 Counts, +/-(0.5% +/- 3 digits)\n", "# - verwendeter Messbereich 2V\n", " crel_U = 0.005\n", " Udigits = 3\n", " Urange = 2\n", " Ucounts = 4000\n", "# - Genauigkeit Strommessung: 2000 Counts, +/-(1.0% +/- 3 digits) \n", "# - verwendete Messbereiche 200µA, 20mA und 200mA \n", " crel_I = 0.010\n", " Idigits = 3\n", " Icounts = 2000\n", " Irange1 = 0.2\n", " Irange2 = 20\n", " Irange3 = 200\n", "# - Rauschanteil (aus Fluktuationen der letzen Stelle)\n", "# - delta U = 0.005 V\n", " deltaU = 0.005\n", "# - delta I = 0.025 mA\n", " deltaI = 0.025\n", "\n", "# - Anzeigegenauigkeit der Spannung (V) \n", " sx = Udigits * Urange / Ucounts\n", " sabsx = np.sqrt(deltaU**2 + sx**2) # Rauschanteil addieren\n", "# - korrelierte Kalibrationsunsicherheit \n", " crelx = crel_U\n", " \n", "# - Anzeigegenauigkeit des Stroms (mA), 3 Messbereiche\n", " sy = np.asarray( 2 * [Idigits * Irange1 / Icounts] + \\\n", " 12 * [Idigits * Irange2 / Icounts] + \\\n", " 6 * [Idigits * Irange3 / Icounts]) \n", " sabsy = np.sqrt(deltaI**2 + sy**2) # Rauschanteil addieren\n", "# - korrelierte Kalibrationsunsicherheit \n", " crely = crel_I\n", "```\n", "\n", "Code für die Anpassung mit *k2Fit* oder *xyFit*:\n", "\n", "```\n", "# thisFit = xyFit ## alternativ \n", " thisFit = k2Fit\n", " model = Shockley\n", "# Anpassung ausführen \n", " fitResult = thisFit(model,\n", " # - data and uncertainties\n", " data_x, data_y, # data x and y coordinates\n", " sx=sabsx, # indep x\n", " sy=sabsy, # indel y\n", " xrelcor=crelx, # correlated rel. x\n", " yrelcor=crely, # correlated rel. y\n", " ref_to_model=True, # reference of rel. uncert. to model\n", " # - fit control\n", " p0=(0.2, 0.05), # initial guess for parameter values \n", " limits=('U0', 0.005, None), # parameter limits\n", "# - output options\n", " plot=True, # plot data and model\n", " plot_cor=False, # plot profiles likelihood and contours\n", " showplots = False, # plt.show() in user code \n", " quiet=False, # suppress informative printout\n", " axis_labels=['U (V)', '$I_D$ (mA) \\ Shockley-Gl.'], \n", " model_legend = 'Shockley-Gleichung'\n", " )\n", " # adjust to different output formats of k2Fit and xyFit\n", " if type(fitResult) is type({}):\n", " parvals, paruncs, cor, chi2, pnames = fitResult.values()\n", " else:\n", " parvals, paruncs, cor, chi2 = fitResult\n", "\n", "```\n", "\n", "Nun sollten wir noch eine Ausgabe in Textformat und als Grafiken vorsehen:\n", "```\n", "# Ausgabe der Ergebnisse in Textform:\n", " print('\\n*==* Fit Result:')\n", " print(\" chi2: {:.3g}\".format(chi2))\n", " print(\" parameter values: \", parvals)\n", " print(\" parameter uncertainties: \\n\", paruncs)\n", " print(\" correlations : \\n\", cor)\n", "\n", "# Anpassung der Optionen für grafische Darstellung\n", " plt.figure(num=1) # activate first figure ...\n", " plt.ylim(-1., 60.) # ... set better y-limit ...\n", " plt.show() # .. and show all figures\n", "```\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Code zur Anpassung der Shockley-Gleichung\n", "# --> hier eingeben" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 3.4.3 Weitere Programmpakete zur Funktionsanpassung\n", "----\n", "\n", "Es gibt auf dem Markt eine ganze Reihe kommerzieller und auch freier Programmpakete zur Funktionsanpassung.\n", "Die bekanntesten sind \n", "- [*origin*](https://www.originlab.com)\n", "- [*gunplot*](http://www.gnuplot.info/)\n", "- [*Grace*](http://plasma-gate.weizmann.ac.il/Grace/)\n", "- [*qtiplot*](https://www.qtiplot.com/)\n", "\n", "Diese Werkzeuge bringen ein eigenes grafisches Interface mit - bisweilen im \"Excel-Look\" - und\n", "sind nach einiger Einarbeitungszeit recht intuitiv zu bedienen. Sie haben allerdings keine\n", "*python*-API (Application Programming Interface); daher ist man bei Datenimport und \n", "Weiterverarbeitung der Ergebnisse auf grafische Methoden wie \"Drag-and-Drop\" angewiesen. \n", "\n", "Eine weitere Schwierigkeit besteht darin, dass diese Werkzeuge nicht auf den Einsatz in \n", "der Physik ausgelegt sind und man die Optionen zur korrekten Behandlung von Unsicherheiten\n", "auf die eingegebenen Daten erst identifizieren und explizit setzen muss! Der Grund ist, dass diese\n", "Programme meist zur Parametrisierung von empirischen Daten eingesetzt werden, wobei folgende \n", "Funktionsweisen auftreten könnten:\n", "- Die anzupassende Modellfunktion muss die Daten innerhalb der ausgegebenen Parameterunsicherheiten richtig beschreiben. Bei der Wahl einer schlecht passenden Modellfunktion werden die Parameterunsicherheiten soweit erhöht, bis eine gute Übereinstimmung erreicht ist - ein Wert vom $\\chi^2/n_f=1$ wird dabei erzwungen!\n", "- Die Unsicherheiten $\\sigma_i$ der Datenpunkte sind unbekannt und werden wie bei der Methode der gaußschen Fehlerquadrate alle auf Eins gesetzt.\n", "\n", "Ein solches möglicherweise als Voreinstellung gesetztes Verhalten muss abgeschaltet werden, um eine Modellanpassung mit Überprüfung der Qualität der Übereinstimmung innerhalb der Messunsicherheiten der Eingabedaten und korrekte Unsicherheiten der Parameter zu erhalten. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## 3.5 Details zum Anpassungspaket *kafe2* \n", "---\n", "\n", "Das in [Abschnitt 3.4.2](#SubSec-Programmpakete) bereits erwähnte Paket *kafe*2 \n", "ist ist eine erweiterte Version des seit dem Jahr\n", "2012 am Instut für Experimentelle Teilchenphysik entwickelten Paktes *kafe* zur\n", "Anpassung von Modellfunktionen an Daten. *kafe2* besitzt mittlerweile einen breitgefächerten\n", "Funktionsumfang, der auch auf komplexe Probleme der Datenauswertung anwendbar ist. \n", "\n", "Unterstützt werden verschiedene Datentypen wie einfache indizierte Daten (also eine Messreihe),\n", "zweidimensionale Datenpunkte (eine Größe *\"x\"* und eine abhängige Größe *\"y\"*) sowie Häufigkeitsverteilungen (Histogramme). \n", "\n", "Bei zweidimensionalen Datenpunkte werden Unsicherheiten sowohl der abhängigen wie\n", "auch der unabhängigen Größe und gegebenenfalls deren Korrelationen unterstützt. \n", "Im Vergleich zu vielen anderen Anpassungswerkzeugen ist diese Möglichkeit \n", "ein Alleinstellungsmerkmal von *kafe(2)*. \n", "\n", "Unterstützt wird die Berücksichtigung von Modellparametern, für die es \n", "Einschränkungen durch externe Messungen gibt. Auch die gleichzeitige Anpassung\n", "mehrerer Modelle mit teilweise gemeinsamen Parametern an verschiedene Datensätze\n", "(\"Multi-Fit\") ist möglich.\n", "\n", "Zur Miminierung des Abstandsmaßes zwischen Daten und Modellfunktion(en)\n", "werden numerische Verfahren angewandt, die aus der quelloffenen, \n", "*Python*-basierten Softwareumgebung *scipy* oder dem am CERN entwickelten\n", "Paket *MINUIT* stammen. Das jeweils minimierte Abstandsmaß (oder auch\n", "die \"Kostenfunktion\") entpricht dem mit einem Faktor Zwei multiplizierten\n", "negativen natürliche Logarithmus der Likelihood-Funktion $(-2\\,\\ln{\\cal L})$ \n", "der Daten für das gegebene Modell. \n", "Für gaußförmige Unsicherheiten der Datenpunkte entspricht dies der Methode\n", "der kleinsten Quadrate (auch \"$\\chi^2$-Methode\"). Andere, auf dem Likelihood-Prinzip \n", "beruhende Kostenfunktionen für die Anpassung von Wahrscheinlichkeitsdichten an\n", "Histogramme oder an indizierte Daten sind ebenfalls verfügbar.\n", "\n", "Zusätlich zur oben diskutierten Methode der \"Punktschätzung\" (also des Auffindens\n", "eines Minimums der Abstandmaßes zwischen Daten und Modellfunktion und Schätzen der\n", "Unsicherheit aus den Werten der zweiten Ableitungen nach den Parametern) bietet\n", "*kafe2* die Möglichkeit, ein- und zweidimensionale Konfidenzintervalle für die\n", "Parameter zu bestimmen, die aus einem Scan des Abstandsmaßes gewonnen werden. \n", "\n", "Trotz dieser vielfältigen Möglichkeiten ist das Interface von *kafe2* recht einfach\n", "gehalten; die generelle Vorgehensweise ist folgende:\n", "\n", " - Definition und Initialisierung eines Daten-Containers für die jeweilige Anpassung \n", " - Erzeugung eines Objekts zur Durchführung der Anpassung, die den\n", " Datencontainer mit einem Modell verbindet \n", " - Durchführung der Anpassung mittels *.do_fit()* und\n", " Ausgabe der Ergebnisse auf der Konsole mittels *.report()* oder\n", " durch direkten Zugriff auf die Ergebnis-Variablen des Fit-Objekts\n", " - Gegebenenfalls Erzeugung und Anzeige von Ergebnisgrafiken mit der generischen\n", " Klasse *Plot*.\n", "\n", "\n", "*kafe2* flexibel über ein *Python*-Interface verwendet werden. Zu *kafe2* gibt es ein\n", "eigenes Tutorial als *jupyter*-Notebook, *kafe2Tutorial.ipynb*, das sich an fortgeschrittene\n", "Nutzer mit komplexeren Problemstellungen wendet.\n", "\n", "Als Anwendungsbeispiel können Sie zunächst die in [Abschnitt 3.4.2](#SubSec-Programmpakete) \n", "schon verwendete Funktion *PhyPraKit.k2Fit()* anschauen, die direkt auf Funktionen aus *kafe2*\n", "zugreift.\n", "\n", "Der folgende Code illustriert die Anpassung einer Exponentialfuktion mit dem Anpassungswerkzeug *kafe2*.\n", "```\n", "from kafe2 import config, XYContainer, XYFit, Plot\n", "import numpy as np, matplotlib.pyplot as plt\n", "\n", "# set better default figure size for kafe2\n", "plt.rcParams['figure.figsize']=[12., 5.] \n", "# !!! must be done after importing kafe2 (will else be overriden)\n", "\n", "\n", "# To define a model function for kafe2 simply write it as a python function\n", "def exponential_model(x, A0=1., x0=5.):\n", " # Our second model is a simple exponential function\n", " # The kwargs in the function header specify parameter defaults.\n", " return A0 * np.exp(x/x0)\n", "\n", "# the data for this exercise\n", "x = [19.8, 3.0, 5.1, 16.1, 8.2, 11.7, 6.2, 10.1]\n", "y = [23.2, 3.2, 4.5, 19.9, 7.1, 12.5, 4.5, 7.2]\n", "data = XYContainer(x_data = x, y_data = y)\n", "data.add_error(axis='x', err_val = 0.3)\n", "data.add_error(axis='y', err_val = 0.15, relative = True)\n", "data.label = 'Datenpunkte'\n", "data.axis_labels=['x-Wert','y-Wert']\n", "\n", "# Create 2 XYFit objects with the same data but with different model functions\n", "exponential_fit = XYFit(data, exponential_model)\n", "exponential_fit.model_label = 'exponentielles Modell'\n", "\n", "# Optional: Assign LaTeX strings to parameters and model functions.\n", "exponential_fit.assign_parameter_latex_names(A0='A_0', x0='x_0')\n", "exponential_fit.assign_model_function_latex_expression(\"{A0} e^{{ {x}/{x0} }}\")\n", "\n", "# Perform the fits.\n", "exponential_fit.do_fit()\n", "\n", "# Optional: Print out a report on the result of each fit.\n", "exponential_fit.report()\n", "\n", "# Optional: Create a plot of the fit results using Plot.\n", "p = Plot(exponential_fit)\n", "\n", "p.plot(fit_info=True)\n", "\n", "# Show the fit results.\n", "plt.show()\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Den Code zur Anpassung einer Exponentialfuktion an Messdaten hier einfügen\n", "# -->" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Um dieses nichtlineare Modell etwas genauer anzuschauen, verwenden wir die Möglichkeit in\n", "*kafe2*, die Likelihood-Kontur im Parameterraum zu scannen:\n", "```\n", "from kafe2 import ContoursProfiler\n", "\n", "# Create contour plot and profiles for the exponential fit\n", "cp = ContoursProfiler(exponential_fit)\n", "cp.plot_profiles_contours_matrix(show_grid_for='contours')\n", "plt.show()\n", "```\n", "Geben Sie diesen Code in die Zelle unten ein - bitte etwas Geduld, die Berechnung der\n", "Likelihood-Konturen braucht Rechenzeit!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# code zum Erstellen von Profil-Likelihood und Kovarianz-Ellipse hier eingeben\n", "# -->" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Die grafische Darstellung zeigt, dass es eine starke Korrelation der Parameter gibt und dass die \n", "Kovarianz-Kontur nur näherungsweise eine Ellipse ist. Auch die Likelihood-Kontur ist im Bereich\n", "des Minimums nicht exakt eine Parabel. Man würde in diesem Fall unsymmetrische Unsicherheiten\n", "der Parameter angeben, die man aus den Schnittpunkten einer horizontalen Line bei Eins mit der\n", "Likelihood-Kurve ablesen kann. Das sich so ergebende Intervall entspricht einem Konfidenzniveau\n", "von 68% (genau wie der Bereich innerhalb von $\\pm 1\\sigma$ der Gaußverteilung).\n", "\n", "Mit dieser Erkenntnis sollte das Ergebnis der Anpassung erneut ausgegeben werden:\n", "Der Aufruf der *report()*-Fuktion mit der Option\n", "*asymmetric_parameter_errors=True* erzeugt eine entsprechende Ausgabe:\n", "```\n", "exponential_fit.report(asymmetric_parameter_errors=True)\n", "p.plot(fit_info=True, asymmetric_parameter_errors= True)\n", "# Show the fit results.\n", "plt.show()\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Code zur Ausgabe asymmetrischer Parameterunsicherheiten hier eingeben:\n", "# -->" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 3.5.1 Simultane Anpassung an mehrere Messreihen mit *kafe2* \n", "\n", "Eine weitere, bisweilen sehr nützliche Fähigkeit von *kafe2* besteht darin, mehrere Messreihen\n", "mit gemeinsamen Parametern gleichzeitig anpassen zu können: das Feature *MultiFit*. Beispiele\n", "sind mehrere Messreihen mit der gleichen Apparatur, um die Eigenschaften von Materialien in \n", "Proben mit unterschiedlicher Geometrie zu bestimmen, wie z. B. die Elastizität oder den \n", "spezifischen Widerstand an Proben mit unterschiedlichen Querschnitten und Längen. \n", "In solchen Fällen sind auf die Apparatur zurückzuführende Unsicherheiten in allen Messreihen \n", "gleich, auch die interessierende Materialeigenschaft ist immer die gleiche, lediglich die \n", "unterschiedlichen Gemoetrie-Parameter und die jeweils bestimmten Werte der Messreihen haben\n", "eigene, unabhängige Unsicherheiten.\n", "\n", "Im folgenden, einfachen Beispiel wird angenommen, dass es einen linearen Zusammenhang zwischen \n", "den Messwerten gibt; die Steigung ist gegeben durch das Produkt aus einer Materialeigenschaft,\n", "$p0$, und einem Geometriefaktor $g1$ bzw. $g2$; außderdem gibt es für jede Messreihe eigene \n", "Offsets $o1$ und $o$, die aus der Geradenanpassung bestimmt werden. \n", "Die beiden Modellfunktionen sehen so aus:\n", "\n", "```\n", "# -- define two model functions with common parameter p0\n", "# remark:\n", "# - p0 might be a material constant,\n", "# e.g. elastic modulus or specific resistance,\n", "# - g might be a geometry factor, like length and/or\n", "# diameter of a sample or a combination of both\n", "# - o might be a nuisance parapeter, e.g. an off-set from noise\n", "# Note that constraints on g1, g2 are needed, i.e. external\n", "# measurents, to give meaningful results\n", "def model1(x, p0=1., g1=1., o1=0):\n", " # p0: Materialeigenschaft\n", " # g1: Geometriefaktor für Probe 1\n", " # o1: offset für Messreihe 1\n", " return g1 * p0 * x + o1 \n", "\n", "def model2(x, p0=1., g2=1., o2=0):\n", " # p0: Materialeigenschaft\n", " # g2: Geometriefaktor für Probe 2\n", " # o2: offset für Messreihe 2\n", " return g2 * p0 * x + o2 \n", "\n", "```\n", "\n", "Natürlich müssen die Geometriefaktoren, ggf. mit Unsicherheiten, bekannt sein oder\n", "separat gemessen werden. In den Anpassungen von *model1* und *model2* werden sie\n", "als eingeschränkte Parameter behandelt. Die Abfolge der Schritte zur Anpassung \n", "ist ansonsten wie üblich, also zunächst das Erzeugen von Daten- und Fitobjekten.\n", "Anders ist lediglich das Zusammenfürhen der Fit-Objekte mit der Zeile \n", "\n", "```# combine the two fit objects to form a MultiFit\n", "multiFit = MultiFit( fit_list=[xyFit1, xyFit2] )\n", "```\n", "Ausführen der Anpassung, Ausgabe der Ergebnisse und die grafische Darstellung\n", "erfolgen dann mit den üblichen Methoden der Instanz des *MultiFit*-Objekts.\n", "\n", "Das vollständige Beispiel ist hier gezeigt:\n", "\n", "```\n", "''' general example for fitting multiple distributions with kafe2\n", " - define models\n", " - set up data objects\n", " - set up fit objects\n", " - set up MultiFit objcet\n", " - perform fit\n", " - show and save output\n", "'''\n", "# Imports #\n", "from kafe2 import XYContainer, Fit, MultiFit, Plot\n", "import numpy as np, matplotlib.pyplot as plt\n", "\n", "\n", "# Workflow #\n", "\n", "# 1. set data \n", "\n", "# data set 1\n", "x1 = [ 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.]\n", "y1 = [ 1.29, 1.78, 3.32, 3.85, 5.27, 6.00, 7.07, 8.57, 8.95, 10.52 ]\n", "e1 = 0.2\n", "e1x = 0.15\n", "c1 = 1.0\n", "ec1 = 0.05\n", "\n", "# data set 2\n", "x2 = [ 1., 2., 3., 4., 5., 6., 7., 8., 9., 10. ]\n", "y2 = [ 0.76, 1.19, 1.71, 2.21, 2.58, 3.01, 3.67, 4.24, 4.69, 4.97 ]\n", "e2 = 0.15\n", "e2x = 0.15\n", "c2 = 0.5\n", "ec2 = 0.05\n", "\n", "# 2. convert to kafe2 data structure and add uncertainties\n", "\n", "xy_d1 = XYContainer(x1, y1)\n", "xy_d1.add_error('y', e1) # independent errors y\n", "xy_d1.add_error('x', e1x) # independent errors x \n", "xy_d2 = XYContainer(x2, y2)\n", "xy_d2.add_error('y', e2) # independent errors y\n", "xy_d2.add_error('x', e2x) # independent errors x\n", "\n", "# set meaningful names \n", "xy_d1.label = 'Beispieldaten (1)'\n", "xy_d1.axis_labels = ['x', 'y (1)']\n", "xy_d2.label = 'Beispieldaten (2)'\n", "xy_d2.axis_labels = ['x', 'y(2) & f(x)']\n", "\n", "# 3. create the individual Fit objects \n", "xyFit1 = Fit(xy_d1, model1)\n", "xyFit2 = Fit(xy_d2, model2)\n", "# set meaningful names for model\n", "xyFit1.model_label = 'Lineares Modell'\n", "xyFit2.model_label = 'Lineares Modell'\n", "# add the parameter constraints\n", "xyFit1.add_parameter_constraint(name='g1', value = c1 , uncertainty = ec1)\n", "xyFit2.add_parameter_constraint(name='g2', value = c2 , uncertainty = ec2)\n", "\n", "# combine the two fit objects to form a MultiFit object\n", "multiFit = MultiFit( fit_list=[xyFit1, xyFit2] )\n", "\n", "# 4. perform the fit\n", "multiFit.do_fit()\n", "\n", "# 5. report fit results\n", "multiFit.report()\n", "\n", "# 6. create and draw plots\n", "multiPlot = Plot(multiFit)\n", "##multiPlot = Plot(multiFit, separate_figures=True)\n", "multiPlot.plot(figsize=(13., 7.))\n", "\n", "# 7. show or save plots #\n", "##for i, fig in enumerate(multiPlot.figures):\n", "## fig.savefig(\"MultiFit-\"+str(i)+\".pdf\")\n", "plt.show()\n", "\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Code für MultiFit-Beispiel hier eingeben\n", "# -->" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "### 3.6 Details zur Modellanpassung mit *PhyPraKit/phyFit* \n", "---\n", "\n", "*phyFit* ist Teil des Pakets *PhyPraKit* (https://readthedocs.org/projects/phyprakit/\n", "und bietet ein so schlank wie möglich gehaltenes Interface zum Anpassungspaket \n", "*iminuit* (https://iminuit.readthedocs.io/en/stable/), das insbesondere auch als Vorlage \n", "zu eigenen Implementierungen dienen kann, die in speziellen Anwendugnsfällen nötig werden. \n", "*iminuit* enthält zahlreiche Methoden zur Minimierung einer in einem hochdimensionalen\n", "Raum definierten skalaren Kostenfunktion, zur Analyse der Unsicherheiten und zur Bestimmung\n", "von Konfidenzintervallen bzw. Konfidenzkonturen. \n", "\n", "Die Basis von *phyFit* bildet die Klasse *mnFit* zum Anpassen eines parameterabhängigen Modells\n", "$f(x, *par)$ an eine Menge von Daten $\\{ x_i\\}$, an Datenpunkte $\\{\\, (x_i, y_i=f(x_i, *par)\\,\\}$ \n", "oder einer Wahrscheinlichkeitsdichtefunktion an Histogram-Daten oder an ungebinnete Daten.\n", "Die Parameterschätzung basiert grundsätzlich auf vorimplementierten\n", "Maximum-Likelihood-Funktionen oder im letzteren Fall auf einer benutzerdefinierten Kostenfunktion. \n", "Klassische Kleinste-Quadrate-Methoden stehen optional zum Vergleich mit anderen Paketen zur Verfügung. \n", "Ein besonderes Merkmal des Pakets ist die Unterstützung verschiedener Arten von Unsicherheiten \n", "für $x-y$-Daten, nämlich unabhängige und/oder korrelierte absolute und/oder relative Unsicherheiten\n", "in $x$- und/oder $y$-Richtung, wie sie auch von *kafe2* unterstützt werden. Die Parameterschätzung\n", "für Dichteverteilungen von histogrammierten Daten basiert auf der verschobenen Poisson-Verteilung, \n", "Poisson$(x - loc, lambda)$ der Anzahl der Einträge in jedem Bin des Histogramms. \n", "\n", "In *phyFit* unterstützt werden eingeschränkte Parameter (zur Einbeziehung von externem Wissen \n", "innerhalb von Gauß'schen Unsicherheiten), Grenzen von Parameterwerten (um problematische Bereiche\n", "im Parameterraum während des Minimierungsprozesses zu vermeiden), und das Fixieren von Parametern.\n", "\n", "Unsicherheiten, die von Modellparametern abhängen, werden korrekt durch dynamische Aktualisierung \n", "des negativen natürlichen Logarithmus der Likelihoodfukntion in jedem Schritt der numerischen\n", "Minimierung behandelt. Datenpunkte mit relativen Fehlern werden als Voreinstellung auf den\n", "Modellwert bezogen, um ansonsten auftretende Verzerrungen zu vermeiden. In Anpassungen an\n", "$x-y$-Daten wird die Ableitung der Modellfunktion nach den Stützstellen $x_i$ verwendet, um \n", "die Kovarianzmatrix der $x_i$ auf die $y$-Achse zu transformieren.\n", "\n", "Die im Paket enthaltenen Wrapperfunktionen *xyFit()*, *hFit()* und *mFit()* veranschaulichen die \n", "Steuerung der Schnittstelle von *mnFit* und erlauben die Anpassungen an (x,y)-Daten, Histogramme \n", "oder unbeginnte Daten mit einem einzigen Funktionsaufruf. \n", " \n", "Ein Beispielskript zur Illustration der einzelnen Schritte bei der Modellanpassung mit der Klasse\n", "*mnFit* aus *PhyPraKit/phyFit* ist hier gezeigt: \n", "```\n", "''' general example for fitting with class mnFit from PhyPraKit.phyFit\n", " - initialize data object\n", " - initialize fit opject \n", " - perform fit (2nd order polynomial)\n", " - show output\n", "'''\n", "# Imports #\n", "from PhyPraKit.phyFit import mnFit\n", "import numpy as np, matplotlib.pyplot as plt\n", "\n", "# fit function definition \n", "def poly2(x, a=1.0, b=0.0, c=0.0):\n", " return a * x**2 + b * x + c\n", "\n", "# set fit function \n", "fitf=poly2 # own definition \n", "\n", "# Workflow #\n", "# 1. load data \n", "x = [0.05, 0.36, 0.68, 0.80, 1.09, 1.46, 1.71, 1.83, 2.44, 2.09, 3.72, 4.36, 4.60]\n", "y = [0.35, 0.26, 0.52, 0.44, 0.48, 0.55, 0.66, 0.48, 0.75, 0.70, 0.75, 0.80, 0.90]\n", "sy = [0.06, 0.07, 0.05, 0.05, 0.07, 0.07, 0.09, 0.10, 0.11, 0.10, 0.11, 0.12, 0.10]\n", "srx = 0.05\n", "# 2. set up Fit object\n", "myFit=mnFit()\n", "# 3. initialize data object\n", "myFit.init_data(x, y, ey=sy, erelx=srx)\n", "# 4. initalize fit object\n", "myFit.init_fit(fitf)\n", "# 5. perform fit\n", "myFit.do_fit()\n", "# 6. retrieve results\n", "FitResult = myFit.getResult()\n", "# 7.make result plot\n", "fig = myFit.plotModel()\n", "# 8.print results\n", "np.set_printoptions(precision=3)\n", "print(\"*==* Result of fit:\\n\", FitResult)\n", "# 9. finally, show plot\n", "plt.show()\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Beispeil für phyFit hier ausprobieren\n", "# -->" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*PhyPraKit* enhält eine Reihe von einfachen Wrapper-Fuktionen, die Anpassungen mit einem\n", "einzigen Funktionsaufruf durchführen und die verschiedenen Anwendsmöglichkeiten von \n", "*phyFit* illustrieren: \n", "\n", " - xyFit() Funktionsanpassung mit (korrelierten) x- und y-Unsicherheiten, \n", " - hFit() Maximum-Likelihood-Anpassung einer Verteilungsdichte an Histogramm-Daten,\n", " - mFit() Anpassung einer Nutzerdefinierten Kostenfunktion oder einer \n", " Verteilungsdichte an ungebinnete Daten mit der maximum-likelood Methode,\n", " - xFit() Anpassung eines Modells an indizierte Daten x_i=x_i(x_j, \\*par).\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4.1 Spezielle Anwendungsfälle von numerischen Anpassungen \n", "\n", "Numerische Optimierungsverfahren sind aus der modernen Datenanalyse nicht mehr wegzudenken.\n", "Einige weitergehende Anwendugsfälle, die auch schon für die Grunglagenpraktika in der Physik\n", "relevant sind, werden im folgenden besprochen. " ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "## 4.1 Transformation von Parametern und Konfidenzkonturen mittels Anpassung \n", "\n", "Eine übliche Problemstellung in der Datenanalyse ist die Weiterverarbeitung oder Interpretation\n", "von an Messdaten angepassten Modellparametern. So müssen häufig Ergebnisse verschiedener Messungen\n", "kombiniert werden, oder die für eine Messung optimierte Wahl der Parameter entspricht nicht\n", "den Parametern, die in fundamentalen theoretischen Modellen verwendet werden. Oft sind auch\n", "Parameter in den Messergebnissen enthalten, die für bestimmte theoretische Studien nicht\n", "relevant sind und daher entfernt werden müssen. So sind die Messgrößen in quantenphysikalischen\n", "Prozessen häufig Raten, denen in der Quantentheorie Summen von Produkten von Kopplungskonstanten\n", "zu Grunde liegen; in solchen Fällen muss \"umparametrisiert\" werden, d.h. es müssen Funktionen\n", "der ursprünglichen Modellparameter berechnet und die Konfidenzbereiche der transformierten\n", "Parameter bestimmt werden.\n", "\n", "In einfachen Fällen reichen dazu die bekannten Formeln zur Fehlerfortpflanzung. \n", "Wenn solche Transformationen von einem Parametersatz $\\{p_i\\}$ auf einen anderen Satz \n", "$\\{p'_j\\}$ nicht-linear sind, dann sind die Konfidenz-Konturen der transformierten Parameter \n", "mit analytischen Methoden allerdings meist nicht leicht zu bestimmen. Die oben besprochenen \n", "Formeln zur linearisierten Fehlerfortpflanzung helfen hier nur bedingt, denn sie liefern eine \n", "Näherung für die Standardabweichung einer unter Umständen von der Gaußverteilung stark abweichenden \n", "Verteilung. Die allgemeine Problemstellung sieht also wie folgt aus:\n", "\n", "Ein oder mehrere Parametersätze $\\{ {p_i}^{(k)} \\}$ sollen durch fundamentalere Parameter \n", "$\\{p'_j\\}$ dargestellt werden. Von Interesse sind dabei nicht nur die transformierten \n", "Parameterwerte selbst, sonderen insbesondere auch deren genaue Konfidenzbereiche. Stellt \n", "man die ursprünglichen Parameter als Funktionen der neuen und ggf. auch der ursprünglichen \n", "Parameter dar, also\n", "$${p_i}^{(k)} = {p_i}^{(k)}\\left( p_i^{(k)}, p'_j \\right)\\, ,$$\n", "so kann man mittels numerischer Anpassung der $p'_j$ die neuen Parameter, deren\n", "Unsicherheitsintervalle und Konfidenzkonturen bestimmen. Wenn es genau so viele\n", "ursprüngliche wie transformierte Parameter gibt, stellt dieses numerische Verfahren\n", "die Lösung eines nicht-linearen Gleichungssystems dar. \n", "\n", "*kafe2* enhält zur Umsetzung dieser Idee ein allgemeines Container-Format, *IndexedContainer*, \n", "das beliebige Eingandsdaten sowie die Berücksichtigung von deren Unsicherheiten und Korrelationen\n", "unterstützt. Damit können als Eingebadaten z.B. auch die Ergebnisse von Anpassungen, also \n", "Ergebnisparameter und deren Kovarianzmatrix, in weiteren Anpassungen verwendet werden, um\n", "mehrere Sätze von Ergebnissen zu mitteln oder um die Parameter zu transformieren. \n", "Dazu wird eine $\\chi^2$-Kostenfunktion minimiert:\n", "$$\\chi^2 = \\left( \\vec{p} - \\vec{p}(\\vec{p}')\\right)^T \n", "{\\bf V}_\\vec{p}^{-1} \n", "\\left(\\vec{p} - \\vec{p}(\\vec{p}')\\right) \\, .$$\n", "Die ursprünglichen Parameter $\\vec{p}$ mit Kovarianzmatrix ${\\bf V}$ werden so durch die \n", "transformierten und deren Kovarianzmatrix V' dargestellt: $\\vec{p'}\\left( \\vec{p} \\right)$. \n", "\n", "Als einfaches Beispiel betrachten wie die Messung eines Ortes in der Ebene in Polarkoordinaten\n", "$(r, \\varphi)$, bestimmt durch Abstandsmessung und Winkelpeilung mit Unsicherheiten.\n", "Wie sieht der $1\\,\\sigma$-Konfidenzbereich, also die entsprechende Kontur, in kartesischen \n", "Koodinaten aus?\n", "\n", "Durch Anpassung von $(x, y)$ an die Messungen $(r\\pm\\sigma_r, \\varphi\\pm\\sigma_\\varphi)$\n", "kann diese Frage beantwortet werden. \n", "Die folgende Funktion sorgt für die Koordinatentransformation: \n", "```\n", "def cartesian_to_polar(x, y):\n", " # determine polar coordinats from cartesian (x,y)\n", " r = np.sqrt(x*x + y*y)\n", " phi = np.arctan2(y, x)\n", " return np.concatenate( (r, phi) )\n", "\n", "```\n", "\n", "Die Daten werden über einen Daten-Container bereit gestellt:\n", "```\n", "from kafe2 import IndexedContainer, Fit, Plot, ContoursProfiler\n", "\n", "# example of parameters: (r, phi) of a space point in polar coordinates\n", "pars = np.array([0.9, 0.755])\n", "puncs = np.array([0.01, 0.15])\n", "indexed_data = IndexedContainer(pars)\n", "indexed_data.add_error(puncs)\n", "```\n", "\n", "Wir wollen die Funktion etwas allgemeiner halten und für den Fall vorsorgen, dass\n", "Messungen mehrerer Punkte übergeben werden, und zwar als aneinandergehängte Arrays\n", "mit $r$ und $\\varphi$-Koordinaten.\n", "```\n", "def cartesian_to_polar(x, y):\n", "\n", " # access data container to find out how many measurements we got\n", " nm = len(indexed_data.data)//2 # expect 2 concatenated arrays (r and phi)\n", "\n", " # determine polar coordinats from cartesian (x,y)\n", " r = np.ones(nm)*np.sqrt(x*x + y*y)\n", " phi = np.ones(nm)*np.arctan2(y, x)\n", "\n", " return np.concatenate( (r, phi) )\n", "```\n", "\n", "Die Durchführung der Anpassung und die Darstellung der Ergebnisse laufen wie üblich: \n", "```\n", "# Set up the fit:\n", "fit = Fit(indexed_data, cartesian_to_polar)\n", "fit.model_label = '$r$-$\\phi$ from x-y'\n", "\n", "# Perform the fit:\n", "fit.do_fit()\n", "\n", "# Report and plot results:\n", "fit.report()\n", "p = Plot(fit)\n", "p.plot()\n", "\n", "contours = ContoursProfiler(fit)\n", "contours.plot_profiles_contours_matrix()\n", "\n", "plt.show()\n", "\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Code für Parameter-Transformation hier eingeben\n", "# -->" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Anzumerken ist noch, dass in diesem Beispiel die Anzahl der Ausgangs- und\n", "Ergebnisparameter gleich ist. Daher ist die Zahl der Freiheitsgrade in diesem\n", "Fall Null, und auch der Wert von $\\chi^2$ am Minimum ist exakt Null. Wenn man\n", "mehrere Messungen der Position in $r-\\varphi$ als Eingabedaten vorgibt, wird\n", "zusätzlich zur Parametertransformation auch eine Mittelung durchgeführt, d.h. es \n", "handelt sich dann um eine echte Parameteranpassung mit einer von Null verschiedenen\n", "Anzahl von Freiheitsgraden. \n", "\n", "\n", "#### Parametertransformation mit *phyFit*\n", "\n", "Auch *phyFit* unterstützt mit der Klasse *xDataContainer* und der Wrapper-Funktion\n", "*xFit()* diese Form von Anpassung an eine allgemeine Menge von Daten $\\{ x_i \\}$\n", "mit durch eine Kovarianz-Matrix $\\bf V$ bestimmten Unsicherheiten. \n", "Von Anwenderseite muss die Darstellung der Daten $\\vec{x}$ als Funktion der Parameter\n", "$\\vec{p}$ bereitgestellt werden. Die Kostenfunktion für die Parameterbestimmung ist\n", "das zweifache des negativen Logarithmus der Likelihood der multivariaten Gaußverteilung,\n", "die für konstante, d.$\\,$h. nicht vom Wert der Daten abhängende Unsicherheiten der\n", "$\\chi^2$-Kostenfunktion entspricht. \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4 Abschließende Anmerkungen \n", "\n", "Mit diesem Ausblick auf fortgeschrittenene Anpassungsprobleme endet das einführende Tutorial \n", "zur Behandlung von statistischen und systematischen Unsicherheiten von physikalischen Messsungen\n", "und zur Bestimmung von Modellparametern aus Messdaten.\n", "\n", "Hoffentlich ist es gelungen, die wesentlichen Grundlagen zu erläutern und in die praktische \n", "Anwendung der vorgestellten Methoden einzuführen und Ihnen damit Hilfestellungen für die\n", "Datenauswertung in den Praktika zur Physik zu geben. \n", "Die Code-Beispiele sind zur freien Anwendung in Ihren eigenen Arbeiten gedacht und \n", "Sie können sie sicher nutzbringend einsetzen. \n", "\n", "---\n", "\n", "Für die Grundlagen-Praktika reichen die in diesem Tutorial beschriebenen\n", "Kenntnisse und Methoden aus. Viele weitere Beispiele, die für Anwendungen \n", "in den Fortgeschrittenenpraktika und in Bachelor-Arbeiten gedacht sind, \n", "werden in den auf diesem Tutorial aufbauenden Tutorials *kafe2Tutorial.ipynb*,\n", "*negLogLFits.ipynb* und *advancedFitting.ipynb* beschrieben. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "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.10.12" } }, "nbformat": 4, "nbformat_minor": 4 }