{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "# Jupyter-Notebook-Tutorial: \n", "\n", "# Anpassung von Modellen an Messdaten mit der Maximum-Likelihood-Methode\n", "\n", " Günter Quast, Juni 2021\n", "***\n", "## Zusammenfassung\n", "\n", "Dieses Tutorial führt in die Grundlage der Anwendung der Maximum-Likelihood-Methode\n", "zur Anpassung an Messdaten ein. Für viele in der Praxis häufig auftretende Problemstellungen\n", "ist die in zahlreichen Anwendungen und Programmpaketen implementierte Methode der \"kleinsten\n", "Fehlerquadrate\" (\"least squares\", LSQ) nicht adäquat. Es sollte statt dessen die auf soliden\n", "mathematischen Grundlagen beruhende Maximum-Likelihood-Methode (\"Maximum Likelihood Estimation\",\n", "MSE) zur Schätzung der Modellparameter und ihrer Vertrauensbereiche angewandt werden. \n", "\n", "Dieses Tutorial in Form eines Jupyter-Notebooks bietet neben kurzen Erläuterungen der \n", "mathematischen Grundlagen insbesondere viele Code-Beispiele in der Sprache *Python*, \n", "die die Maximum-Likelihood-Methode an typischen Anwendungsbeispielen erläutern und zur \n", "Verwendung in eigenen Anpassungsproblemen modifiziert und erweitert werden können.\n", "\n", "Der Beispielcode nutzt das *Python*-Paket [iminuit](https://iminuit.readthedocs.io/en/stable/)\n", "zur numerischen Optimierung und Analyse der Unsicherheiten sowie die im wissenschaftlichen\n", "Rechnen mit *Python* weit verbreiteten Bibliotheken *numpy*, *scipy* und *matplotlib*.\n", "Als anwendungsorientierte Zwischenebene wird das Paket \n", "[PhyPraKit.phyFit](https://phyprakit.readthedocs.io/en/latest/) verwendet.\n", "\n", "\n", "\n", "## Inhaltsübersicht:\n", "\n", "1. Anpassung von Modellen an Messdaten\n", "\n", "2. Mathematische Grundlagen der Parameterschätzung\n", "\n", " 2.1 Parameterschätzung mit dem Maximum-Likelihood-Verfahren\n", "\n", " 2.2 Beispiel 1: Anpassung einer Normalverteilung\n", "\n", " 2.3 Bestimmung der Parameterunsicherheiten \n", "\n", " 2.4 Eigenschaften der Maximum-Likilihood-Schätzung\n", " \n", " 2.5 Intervallschätzung\n", "\n", " - Bestimmung der Parameterunsicherheiten mit der Profil-Likelihood\n", "\n", " - Konfidenzkonturen \n", "\n", " - Beispiel 2: log-Likelihood Anpassung einer Exponentialfunktion\n", "\n", "3. Anpassung von Verteilungsdichten an Histogramme\n", "\n", " 3.1 Beispiel \n", " \n", " 3.2 Signifikanz eines Signals und Ausschlussgrenze \n", "\n", "4. Maximum-Likelihood und Methode der kleinsten Fehlerquadrate\n", "\n", " 4.1 Konstruktion der Kovarianzmatrix\n", "\n", " 4.2 Beispiel zur Anpassung an x/y-Daten mit komplexen Unsicherheiten \n", "\n", "5. Beurteilung der Qualität von Anpassungsverfahren\n", "\n", " 5.1 Überprüfung mittels Monte-Carlo-Methode (\"toy-MC\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## Grundsätzliches zu 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.\n", "Die Eingabe von Formeln im *LaTeX*-Format wird ebenfalls unterstützt.\n", "\n", "Eine Zusammenstellung der wichtigsten Befehle zur Verwendung von *Jupyter* als Arbeitsumgebung\n", "findet sich im Notebook\n", "[*JupyterCheatsheet.ipynb*](https://git.scc.kit.edu/yh5078/datenanalyse/-/blob/master/jupyter/JupyterCheatsheet.ipynb).\n", "Grundlagen zur statistischen Datenauswertung finden sich in den Notebooks \n", "[*IntroStatistik.ipynb*](https://git.scc.kit.edu/yh5078/datenanalyse/-/blob/master/jupyter/IntroStatistik.ipynb)\n", "und\n", "[*Fehlerrechnung.ipynb*](https://git.scc.kit.edu/yh5078/datenanalyse/-/blob/master/jupyter/Fehlerrechnung.ipynb).\n", "\n", "In *Jupyter* werden Code und Text in jeweils einzelne Zellen eingegeben. \n", "Aktive Zellen werden durch einen blauen Balken am Rand angezeigt.\n", "Zellen können sich in zwei Zuständen befinden: im Edit-Mode ist das Eingabefeld weiß, \n", "im Command-Mode ist es ausgegraut.\n", "Durch Klicken in den Randbereich wird der Command-Mode gewählt, ein Klick in das Textfeld einer\n", "Code-Zelle schaltet in den Edit-Mode.\n", "Die Taste `esc` kann ebenfalls verwendet werden, um den Edit-Mode zu verlassen.\n", "\n", "Die 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.\n", "Die Eingabe von `m` im Command-Mode setzt den Typ Markdown, Eingabe von `y` wählt den Typ Code.\n", "\n", "Prozessiert - also Text gesetzt oder Code ausgeführt - wird der Zelleninhalt durch Eingabe von\n", "`shift+return`, oder auch `alt+return` wenn zusätzlich eine neue, leere Zelle erzeugt werden soll.\n", "\n", "Die hier genannten Einstellungen sowie das Einfügen, Löschen oder Ausführen von Zellen sind\n", "auch über das PullDown-Menü am oberen Rand verfügbar.\n", "\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Abhängigkeiten\n", "\n", "Dieses Notebook nutzt Module aus dem Paket \n", "[PhyPraKit](https://github.com/GuenterQuast/PhyPraKit) vers. >=1.3.0, \n", "das von den Pakten *numpy*, *scipy*, *matplotlib* und *iminuit* abhängt.\n", "\n", "---\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Likelihood-Anpassung von Modellen an Messdaten\n", "***\n", "Die Physik als \"Theorie-geleitete Erfahrungswissenschaft\" lebt vom\n", "Zusammenspiel theoretischer Hypothesenbildung und deren stringenter\n", "Überprüfung im Experiment. Bei der Interpretation von Messdaten kommt\n", "es daher zunächst darauf an, ob sie überhaupt einem erwarteten Modell \n", "entsprechen. Erst, wenn diese Frage mit \"ja\" beantwortet werden kann, \n", "ist die Extraktion vom Modellparametern sinnvoll.\n", "\n", "Leider tragen die meisten der häufig verwendeten Werkzeuge in der \n", "Voreinstellung dieser Anforderung nicht Rechnung, sondern zielen \n", "auf eine möglichst genaue Parametrisierung der Daten ab. \n", "Bei Vorgabe von unpassenden Modellen werden dazu die Parameterfehler so \n", "weit vergrößert, dass die Parametrisierung trotzdem möglichst genau passt.\n", "\n", "Um den besonderen Bedingungen der Datenauswertung in der Physik\n", "gerecht zu werden, ist die Verwendung von eigenem Programmcode, \n", "idealerweise in der Programmiersprache *Python* in einem \n", "Jupyter-Notebook, häufig notwendig, da die meisten mittels einer\n", "grafischen Oberfläche zu bedienenden Werkzeuge zur statistischen\n", "Datenauswertung den Anforderungen nicht oder nur zum Teil genügen. \n", "\n", "Wichtige Anforderungen an eine Methode zur Modellanpassung sind:\n", "\n", "- Alle Messunsicherheiten sollten direkt bei der Modellanpassung berücksichtigt\n", " und entsprechend modelliert werden, um sie vollständig und korrekt auf die\n", " Unsicherheiten der Parameter zu propagieren. Eine später von Hand\n", " durchgeführte Fehlerfortpflanzung ist aufwändig und funktioniert nur\n", " in einfachsten Fällen. \n", "\n", "- Häufig treten auch relative Unsicherheiten als Bruchteil des wahren \n", " Messwerts auf, z.$\\,$B. bei Kalibrationsunsicherheiten o.$\\,$Ä..\n", " Solche relativen Unsicherheiten korrekt zu behandeln, gelingt mit der\n", " in vielen Anpassungsprogrammen implementierten Methode der kleinsten \n", " Fehlerquadrate nur noch eingeschränkt. \n", " Bessere Parameterschätzungen nutzen die Maximum-Likelihood-Methode \n", " (\"MLE\", Maximum Likelihood Estimation) mit an das Problem angepassten\n", " Likelihood-Funktionen. \n", "\n", "- Sehr oft, insb. in der Quanten-, Kern- und Teilchenphysik, treten als \n", " Messgrößen Zählraten in Abhängigkeit von Energie, Streuwinkel oder anderen\n", " Größen auf, die üblicherweise als Häufigkeitsverteilung (Histogramm)\n", " dargestellt werden. Die Einträge in einem Intervall eines Histogramms\n", " (einem sog. `Bin`) sind aber nicht Gauß-, sondern Poisson-verteilt. \n", " Auch näherungsweise können Bins mit Null Einträgen mit einer - oft immer\n", " noch verwendeten - angepassten Methode der kleinsten Fehlerquadrate\n", " nicht berücksichtigt werden. Zur Anpassung von Modellen an Histogramme\n", " müssen also ebenfalls maximum-likelihood Methoden angewendet werden.\n", " \n", "- Bei der Analyse von sehr kleinen Zählraten bei nur selten eintretenden\n", " Ereignissen kann schon die Notwendigkeit zur Aufteilung in Bins zur einem\n", " Informationsverlust und zu einer Verzerrung des Ergebnisses führen.\n", " In solchen Fällen sollte die Möglichkeit zur Durchführung eines ungebinnten\n", " Maximum-Likelihod-Fits zur Verfügung stehen.\n", " \n", "- Bei der Anpassung an zweidimensionale Datenpunkte $(x_i, y_i)$ sollten \n", " die verwendeten Werkzeuge Unsicherheiten in Abszissenrichtung (entlang\n", " der \"x-Achse\") zusätzlich zu denen in Ordinatenrichtung behandeln können. \n", " Die Unterstützung von korrelierten Unsicherheiten ist ebenfalls notwendig, \n", " um typische Charakteristika von Messgeräten, wie korrelierte \n", " Kalibrationsunsicherheiten zusätzlich zu den unabhängigen \n", " Digitalisierungsunsicherheiten und Rauschbeiträgen, darstellen zu können.\n", "\n", "Leider gibt es kaum Anpassungswerkzeuge, die alle genannten Anforderungen\n", "gleichzeitig erfüllen. Deshalb wurde am ETP in zahlreichen Bachelor-Arbeiten\n", "ein quelloffenes Anpassungswerkzeug, [kafe2](https://github.com/dsavoiu/kafe2)\n", "entwickelt, das die genannten Aspekte abdeckt. Allerdings ist der Programmcode\n", "wegen der großen Zahl an Optionen und Methoden sehr komplex und für Einsteiger\n", "unübersichtlich.\n", "\n", "Zur Illustration der Vorgehensweise wird daher in diesem Tutorial ein\n", "einziges, sehr flexibles Werkzeug zur numerischen Optimierung und \n", "Analyse der Parameterunsicherheiten verwendet, nämlich das am CERN \n", "entwickelte Paket `Minuit` bzw. dessen Python-Interface \n", "[iminuit](https://iminuit.readthedocs.io/en/stable/). \n", "Eine schlanke Implementierung von Anwendungsbeispielen bietet das Paket \n", "[PhyPraKit.phyFit](https://phyprakit.readthedocs.io/en/latest/), das \n", "hier verwendet wird, um die Grundlagen der Anpassung von Modellen an \n", "Messdaten mit der Maximum-Likelihood Methode darzustellen. \n", "\n", "***" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Mathematische Grundlagen der Parameterschätzung\n", "***\n", "\n", "In diesem Tutorial werden ausgehend vom Maximum-Likelihood-Prinzip \n", "zunächst typische Anwendung zur Anpassung von Verteilungsdichten an\n", "Messdaten diskutiert und dann die Verbindung zur altbekannten Methode\n", "der kleinsten Quadrate zur Anpassung von parameterbehafteten \n", "Modellfunktionen $y_i = f(x_i; *par )$ an zweidimensionale Datenpunkte\n", "$(x_i, y_i)$ dargestellt. \n", "\n", "Ausgangspunkt der Überlegungen ist der auf der Maximum-Likeliood Methode\n", "basierende Formalismus zur Parameterschätzung von Verteilungsdichten, \n", "die den als statistische Stichprobe betrachteten Messergebnissen zu \n", "Grunde liegen. \n", "Manchmal sind die Parameter der Verteilungsdichte selbst interessant, \n", "z.$\\,$B. bei der Lebensdauerverteilung von quantenmechanischen Zuständen \n", "oder von (Elementar-)Teilchen.\n", "Häufig treten Parameter von Verteilungsdichten aber auch nur als \n", "\"Störparameter\" auf, die statistische Fluktuationen um einen angenommenen\n", "wahren Wert beschreiben. Die Behandlung von Messunsicherheiten von \n", "Datenpunkten $(x_i, y_i)$ ist dafür ein Beispiel. \n", "\n", "\n", "## 2.1 Parameterschätzung mit dem Maximum-Likelihood Verfahren \n", "\n", "Beginnen wir an dieser Stelle mit der einfachsten Problemstellung, \n", "der Schätzung der Parameter der Verteilungsdichte einer Grundgesamtheit,\n", "aus der eine endliche Stichprobe (=\"Daten\") gezogen werden. \n", "\n", "Wir gehen von einer Verteilungsdichte ${\\rm pdf}(x, \\vec p)$ aus,\n", "die von einer Anzahl von Parametern $p_j$ abhängt. Ein Datensatz\n", "$\\vec x$ besteht aus ${n_d}$ voneinander unabhängigen\n", "Einzelmessungen $x_i$.\n", "\n", "Nach dem Maximum Likelihood-Prinzip ist der beste Schätzwert der\n", "Parameter ${\\hat{ \\vec {p}}}$ gegeben durch die Werte, die die\n", "sog. `Likelihood`, das Produkt der Werte der $pdf$ für die\n", "Einzelmessungen, ${\\rm pdf}(x_i, \\vec p)$, maximieren:\n", "\n", "\\begin{equation}\\label{equ:Likelihood1}\n", " {\\cal L}(\\vec x, \\vec p) =\n", " \\displaystyle\\prod_{i=1}^{n_d} \\, {\\rm pdf}(x_i, \\vec p) \\,.\n", "\\end{equation} \n", "\n", "Um unhandlich kleine Werte des Produkts zu vermeiden, aus Gründen\n", "der numerischen Stabilität und auch aus mathematischen Gründen\n", "(s. \"Fisher-Information\") verwendet man den negativen Logarithmus\n", "von $\\cal L$, die \"negative log-Likelihood\" $nl\\cal L$, gegeben durch\n", "\\begin{equation}\\label{equ:Likelihood2}\n", " {nl\\cal L}(\\vec x, \\vec p) =\n", " -\\displaystyle\\sum_{i=1}^{n_d} \\, \\ln\\left({\\rm pdf}(x_i, \\vec p)\\right) \\,.\n", "\\end{equation} \n", "\n", "Die negative log-Likelihood-Funktion wird als Funktion der Parameter\n", "aufgefasst und bezüglich der Parameter minimiert. Üblicherweise verwendet\n", "man sie als sogenannte \"Kostenfunktion\" in numerischen Optimierungsverfahren.\n", "\n", "In der Physik ist es auch üblich, einen mit Zwei multiplizierten Wert zu \n", "verwenden, ${n2l\\cal L}(\\vec x; \\vec p) = 2 \\cdot {nl\\cal L}(\\vec x, \\vec p)$.\n", "Wie wir unten sehen werden, entspricht dies in Spezialfällen dem Wert der\n", "Residuensumme, $S$, der nach dem von Gauß vorgeschlagenen Verfahren der \n", "kleinsten Fehlerquadrate häufig als Kostenfunktion verwendet wird. \n", "\n", "Es lässt sich zeigen, dass das log-Likelihood-Verfahren unter allen \n", "Verfahren zur Parameterschätzung optimal ist, also die kleinste Varianz \n", "der Parameterwerte liefert. Auch dies ist ein Grund, dem Likelihood-Verfahren\n", "den Vorzug zu geben, wenn man maximalen Informationsgewinn aus bisweilen\n", "extrem aufwändigen und damit teuren Experimenten ziehen möchte.\n", "\n", "\n", "#### Likelihood und Wahrscheinlichkeitsdichte. \n", "Die Likelihoodfunktion ist im klassischen Sinn der frequentistischen Statistik \n", "keine Wahrscheinlichkeitsdichte der Parameter $\\vec p$, denn erstens ist\n", "sie nicht normiert, und zweitens setzt sie sich aus Produkten von \n", "Wahrscheinlichkeiten zusammen, einen Wert $x_i$ zu beobachten wenn\n", "die Parameterwerte vorgegeben sind. \n", "In der Interpretation der Bayes'schen Statistik\n", "kann man aber unter Verwendung einer A-priori-Wahrscheinlichkeit, eines\n", "sogenannten \"Bayes-Priors\" eine solche Wahrscheinlichkeitsdichte angeben.\n", "Dazu wird das Bayes'sche Theorem verwendet:\n", "\n", "$\\displaystyle {\\cal P}( p | \\vec x) = \\frac \n", " { {\\cal L}(\\vec x | p) \\, \\cdot \\, {\\cal P}_p(p)} \n", " {\\cal P_x(\\vec x)} \\, $.\n", "\n", "In diesem Fall ist ${\\cal P_x(\\vec x)}$ eine Normierungskonstante, die dafür\n", "sorgt, dass ${\\cal P}( p | \\vec x)$, die Wahrscheinlichkeit einen Parameterwert\n", "$p$ zu beobachten, wenn die Daten festliegen. Wenn man z.B. annehmen kann, dass \n", "der Parameterbereich beschränkt ist und die Prior-Verteilung innerhalb dieses \n", "Bereichs flach ist, dann ist ${\\cal P}_p(p)$ eine Gleichverteilung, also eine Konstante.\n", "Es ist aber auch möglich und ein häufig genutztes Verfahren, Vorwissen über \n", "die Verteilung des Parameters $p$ aus anderen Untersuchungen als Verteilungsdichte\n", "zu berücksichtigen. \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2.1 Beispiel 1: Anpassung einer Normalverteilung\n", "\n", "Wir betrachten als erstes, einfache Beispiel Daten $(x_i)$, an die eine Gauß'sche \n", "Normalverteilung angepasst werden soll. \n", "\n", "Die Gaußverteilung mit Erwartungswert $\\mu$ und Standardabweichung $\\sigma$ ist \n", "\n", "${\\cal N}(x; \\mu, \\sigma) \\,=\\, \\frac{1}{\\sqrt{2\\pi} \\, \\sigma} \\exp -{\\frac{(x-\\mu)^2}{2 \\sigma^2}}\\,$.\n", "\n", "Das Zweifache des negativen natürlichen Logarithmus der Likelihood davon ist\n", "\n", "$-2\\, \\ln({\\cal N}(x; \\mu, \\sigma))\\,=\\, \\left(\\frac{x-\\mu}{\\sigma}\\right)^2 + 2\\ln(\\sigma) + \\ln (2\\pi)\\,$.\n", "\n", "Den von allen Variablen unabhängigen Term $\\ln(2\\pi)$ kann man weglassen, \n", "weil er auf die Lage des Minimums keinen Einfluss hat. \n", "\n", "Die Kostenfunktion $n2l{\\cal L}$ ist die Summe solcher Ausdrücke über alle $n_d$ Werte\n", "einer Messreihe:\n", "\n", "$n2l{\\cal L} = 2 n_d \\ln(\\sigma) + \n", " \\displaystyle \\sum_i^{n_d} \\left(\\frac{x_i-\\mu}{\\sigma}\\right)^2\\,$.\n", "\n", "Durch partielles Differenzieren dieser Funktion nach den Parametern $\\mu$ und $\\sigma$\n", "und Nullsetzen der Ableitungen kann man analytisch die optimalen Parameterwerte bestimmen.\n", "\n", "Wir wollen zunächst den Schätzwert für den Parameter $\\mu$ bestimmen und $\\sigma$ als bekannt\n", "und fest voraussetzen; $\\sigma$ könnte z.B. die Auflösung einer Messapparatur sein, mit der\n", "wir mehrfach Messungen der Größen $x_i$ durchführen. \n", "Man erhält nach kurzer Rechnung (die man in jedem Buch zur Statistik findet) für den\n", "Maximum-Likelihood Schätzwert $\\hat{\\mu}$:\n", "\n", "$\\hat{\\mu} = \\displaystyle \\frac{1}{n_d} \\sum_i^{n_d} x_i \\,$ . \n", " \n", "Das ist die erwartete Formel für den Mittelwert der Stichprobe, die also dem MLE-Schätzwert entspricht.\n", "\n", "Wir fassen nun die Likelihood als Funktion von $\\mu$ auf. Durch Einsetzen des Schätzwerts $\\hat{\\mu}$\n", "und nach kurzer Umformung erhält man:\n", "\n", "$\\displaystyle n2l{\\cal L}(\\mu) = \n", " 2 n_d \\ln(\\sigma) + \\frac{1}{n_d} \\sum_i^{n_d} ({x_i}^2 - \\hat{\\mu}^2) + \n", " \\frac{n_d}{\\sigma^2} (\\mu -\\hat{\\mu})^2 \\,=\\, C \\,+\\, \\frac{n_d}{\\sigma^2} (\\mu -\\hat{\\mu})^2\\,, $\n", " \n", "wobei mit $C$ alle Terme abgekürzt sind, die nicht vom Parameter $\\mu$ abhängen.\n", "Die erhaltene Darstellung ist eine Parabel mit Scheitelpunkt bei $\\mu=\\hat{\\mu}$ und\n", "Krümmung $\\frac{2 n_d}{\\sigma^2}$, deren Minimum den Wert $n2l{\\cal L}(\\hat{\\mu}) = C$ hat.\n", "Den Vorfaktor des quadratischen Terms identifizieren wir als das Inverse der Varianz des\n", "Mittelwerts, $V_\\hat{\\mu} = {\\sigma_\\hat{\\mu}}^2 = \\sigma^2/n_d$. Damit kann man auch \n", "schreiben \n", "$\\displaystyle \n", "n2l{\\cal L}(\\mu) =\n", "C \\,+\\, \\frac{1}{{\\sigma_\\hat{\\mu}}^2} (\\mu -\\hat{\\mu})^2$. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Wenn man auch noch die Auflösung $\\sigma$ als unbekannt annimmt und den Schätzwert aus den\n", "Daten bestimmen möchte, so ist das Mimimum von $n2l{\\cal L}(x_i; \\mu, \\sigma)$ durch Nullsetzen\n", "beider partieller Ableitungen nach $\\mu$ und $\\sigma$ gegeben. Für den Schätzwert $\\hat{\\sigma}$ \n", "erhält man: \n", "\n", "$\\hat{\\sigma^2} = \\displaystyle \\frac{1}{n_d} \\sum_i^{n_d} (x_i -\\hat{\\mu})^2 = \\sum_i^{n_d} {x_i}^2 - {\\hat{\\mu}}^2 \\,.$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Numerisches Verfahren zur Bestimmung des Optimums\n", "\n", "Weil es in diesem Tutorial um numerische Methoden zu Behandlung realistischer \n", "Problemstellungen gehen soll, die rein analytisch nicht mehr zugänglich sind, \n", "nutzen wir dieses bekannte und einfache Beispiel, um die numerische Herangehensweise \n", "zu illustrieren. Die dazu verwendeten Methoden können auch recht leicht auf \n", "allgemeinere und komplexere Probleme mit komplizierteren Likelihood-Funktionen\n", "angewandt werden. \n", "\n", "Mit Hilfe des folgen Codefragments werden zunächst (Pseudo-)Daten erzeugt:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "# generate Gaussian-distributed data\n", "\n", "mu0 = 2.0\n", "sig0 = 0.5\n", "np.random.seed(314159) # initialize random generator\n", "data = mu0 + sig0 * np.random.randn(100)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Wie bei der Analyse von Daten aller Art üblich, sollte man sich zunächst einen Überblick\n", "über die Daten verschaffen. Dazu gehören die Anzahl der Datenpunkte, Mittelwert, Median\n", "und Standardabweichung.\n", "Zur Visualisierung eignent sich die Darstellung als sog. \"BoxPlot\", der den Median, das\n", "zentrale 50%-Quantil als Rechteck und den Bereich der Daten als \"Antennen\" zeigt. \n", "Extremwerte, die weiter vom Rand des Rechtecks entfernt liegen als der Länge des Rechtecks \n", "entspricht, werden als als individuelle Datenpunkte (\"Ausreißer\") angezeigt.\n", "Im Code-Beispiel unten werden zusätzlich die einzelnen Datenwerte als senkrechte Linien sowie \n", "Mittelwert und Standardabweichung als \"Fehlerbalken\" dargestellt. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "\n", "# determine characteristics of the data\n", "nd = len(data)\n", "mn = min(data)\n", "mx = max(data)\n", "mean = data.mean()\n", "med = np.median(data)\n", "std = data.std()\n", "\n", "print(\"\\n*==* Data characteristics\")\n", "print(\" number of data points:\", nd)\n", "print(\" minimum: :\", mn)\n", "print(\" maximum: :\", mx)\n", "print(\" mean of data : {:.3g}\".format(mean))\n", "print(\" median : {:.3g}\".format(med))\n", "print(\" standard deviation : {:.3g}\".format(std))\n", "\n", "print(\"\\nBoxPlot of data\")\n", "\n", "plt.boxplot(data, vert=False, whis=1.0)\n", "plt.vlines(data, 0.0, 0.8, lw=1, color=\"grey\", alpha=0.5, label=\"x\")\n", "plt.errorbar(mean, 0.4, marker=\"|\", xerr=std)\n", "tx = plt.xlabel(\"x\", size=\"x-large\")\n", "plt.gca().yaxis.set_visible(False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Wenn die Stichprobe groß genug ist, wie in unserem Fall, dann ist eine \n", "Häufigkeitsverteilung ebenfalls aussagekräftig, wie im folgenden Beispiel gezeigt. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot histogram of the data\n", "if nd < 30:\n", " print(\"\\nHistogram of data not meaningful\")\n", "else:\n", " print(\"\\nHistogram of data\")\n", " nb = min(nd // 10, 50)\n", " bc, be, _ = plt.hist(data, nb, color=\"lightgrey\", rwidth=0.95)\n", " plt.errorbar(mean, 0.6065 * max(bc), marker=\"|\", xerr=std)\n", " tx = plt.xlabel(\"x\", size=\"x-large\")\n", " ty = plt.ylabel(\n", " \"frequency\",\n", " size=\"x-large\",\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Zur Anpassung einer Normalverteilung an die Stichprobendaten verwenden wir \n", "nun ein numerisches Verfahren, das auf der Maximierung der Likelihood bzw.\n", "Minimierung des negativen natürlichen Logarithmus der Likelihood beruht. \n", "In der numerischen Optimierung ist es üblich, eine solche zu optimierende\n", "skalare Funktion als \"Kostenfunktion\" zu bezeichnen. \n", "\n", "Die Implementierung von $n2l{\\cal L}$ als zu minimierende Kostenfunktion \n", "in *Python* ist einfach:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# define cost function: 2 * negative log likelihood of Gauß distribution\n", "def myCost(mu=1.0, sigma=1.0):\n", " # simple -2*log-likelihood of a 1-d Gauss distribution\n", " r = (data - mu) / sigma\n", " return np.sum(r * r) + 2.0 * len(data) * np.log(sigma)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In der folgenden Codezelle wird das Paket MINUIT mit Hilfe des *Python*-basierten \n", "Interfaces *iminuit* zum Auffinden des Minimums der Kostenfunktion verwendet. \n", "Neben der Minimierung einer skalaren, von Parametern abhängigen Kostenfunktion \n", "bietet *iminuit* auch umfangreiche Methoden zur Steuerung der Anpassung und zur \n", "Analyse der Parameterunsicherheiten. \n", "\n", "Für unser einfaches Problem sieht der Code zur Verwendung von *iminuit* z.B. so aus:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# code for iminuit vers. >2.0\n", "from iminuit import Minuit\n", "\n", "# initialize Minuit object\n", "m = Minuit(myCost, mu=1.0, sigma=1.0)\n", "m.errordef = 1.0 # internal parameter, needed to control uncertainty analysis\n", "\n", "# perform optimization\n", "m.migrad()\n", "\n", "# print results\n", "print(\"parameter names: \", m.parameters)\n", "print(\"best-fit values: \", tuple(m.values))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Technisch einfacher wird die Anpassung mit Hilfe einer zusätzlichen Software-Ebene, \n", "die *Minuit* und das mächtige, aber damit auch notwendigerweise komplizierte Interface \n", "kapselt und so typische Anwendungsfälle vereinfacht, wie z.$\\,$B. das \n", "Werkzeug [*kafe2*](https://github.com/dsavoiu/kafe2), \n", "oder das schlanke, im Paket *PhyPraKit* enthaltene Paket *phyFit*. \n", "Für Anpassungen mit *iminuit* bei vom Anwender vorgegebener Kostenfunktion ist die \n", "Interface-Funktion *mFit()* geeignet. Übergeben wird lediglich die Kostenfunktion; \n", "Rückgabewerte sind die aus den Keyword-Argumenten der Kostenfunktion bestimmten Namen \n", "und die Schätzwerte der Parameter. Zurückgegeben wird weiter eine Schätzung der \n", "Parameterunsicherheiten und eine Größe zur Einschätzung der Qualität der Anpassung. \n", "Über diese letztgenannten Größen wird weiter unten noch zu sprechen sein. \n", "\n", "Das Anpassungsbeispiel mit *mFit* benötigt klarerweise weniger eigenen Code\n", "und sieht wie folgt aus: " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from PhyPraKit.phyFit import mFit\n", "\n", "fit_results = mFit(myCost)\n", "\n", "# Print results\n", "print(\"\\n*==* user-defined cost: Fit Result:\")\n", "print(\" parameter names: \", fit_results[\"parameter names\"])\n", "print(\" parameter values: \", fit_results[\"parameter values\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Die erhaltenen Werte sind identisch zu denen im Beispiel vorhin. \n", "Sie stimmen übrigens perfekt mit den erwarteten Ergebnissen aus den \n", "analytischen Betrachtungen von oben überein: \n", "$\\mu$ ist ist der Erwartungswert und $\\sigma$ ist die Standardabweichung\n", "der Werte der Stichprobe. \n", "\n", "Die Parameterwerte sind allerdings nicht identisch zu denen, die für die Erzeugung\n", "der Pseudo-Daten verwendet wurden. Die Abweichungen sind natürlich durch die endliche \n", "Größe der Stichprobe erklärbar." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Als letzen Schritt sollte man das erhaltene Ergebnis visuell überprüfen.\n", "Dazu verwenden wir die Histogrammdarstellung und überlagern die angepasste\n", "Funktion." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot histogram of the data and best-fit curve\n", "from scipy.stats import norm\n", "\n", "if nd < 30:\n", " print(\"\\nHistogram of data not meaningful\")\n", " bw = 1\n", "else:\n", " print(\"\\nHistogram of data\")\n", " nb = min(nd // 10, 50)\n", " bc, be, _ = plt.hist(data, nb, color=\"lightgrey\", rwidth=0.95, label=\"data\")\n", " bw = be[1] - be[0]\n", " plt.errorbar(mean, 0.6065 * max(bc), marker=\"|\", xerr=std)\n", "\n", "# best-fit parameters\n", "mu = fit_results[\"parameter values\"][0]\n", "sigma = fit_results[\"parameter values\"][1]\n", "# plot noraml distribution with these parameters\n", "xp = np.linspace(mn, mx, 150)\n", "plt.plot(xp, nd * bw * norm.pdf(xp, loc=mu, scale=sigma), color=\"orange\", label=\"fit\")\n", "plt.legend(loc=\"best\")\n", "tx = plt.xlabel(\"x\", size=\"x-large\")\n", "ty = plt.ylabel(\n", " \"frequency\",\n", " size=\"x-large\",\n", ")" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "## 2.3 Bestimmung der Parameterunsicherheiten\n", "\n", "In diesem Abschnitt werden Methoden beschrieben und untersucht, mit \n", "deren Hilfe die statistischen Unsicherheiten der Schätzungen der\n", "optimalen Parameterwerte untersucht werden können. Die Parameterschätzung\n", "beruht ja auf einem Algorithmus, der als Eingabe Zufallsgrößen nutzt.\n", "Daher sind die Ergebnisse der Schätzung selbst Zufallsgrößen, die eine \n", "nicht verschwindende Varianz haben. Das Ergebnis ist also erst vollständig\n", "und als Ergebnis einer physikalischen Studie überhaupt erst vewendbar, wenn\n", "auch die Unsicherheiten der Parameterwerte bekannt sind. \n", "\n", "Die Varianz der geschätzten Parameterwerte, oder besser die Kovarianzmatrix\n", "der Unsicherheiten aller Parameter, lässt sich aus dem Verlauf der\n", "log-Likelihood am Minimum bestimmen. Anschaulich ist ein Parameter\n", "um so genauer bestimmt, je \"schärfer\" das Minimum ausgeprägt ist. \n", "Diesen Zusammenhang hatten wir oben bei der Darstellung der negativen\n", "log-Likelihood Funktion für den Parameter $\\mu$ einer Normalverteilung\n", "schon vermutet. Einen mathematischen Beweis und eine allgemeingültige\n", "Grenze für die Varianz haben unabhängig voneinander H. Cramér und C.R. Rao \n", "sowie später auch andere hergeleitet.\n", "In der für uns passenden Form geschrieben, gilt für die Varianz einer\n", "unverzerrten Schätzung $\\hat{p}$ eines Parameters $p$:\n", "\n", "${\\rm Var}_{\\hat p} = {\\sigma_\\hat{p}}^2 \\ge \n", "\\left[ { - \\left. \\frac { {\\partial}^2 \\ln{\\cal{L} }}\n", "{{\\partial p}^2 } \\right|}_{\\hat p}\\right]^{-1}\\,. $\n", "\n", "Die kleinstmögliche Varianz eines Schätzwerts ist also durch das Inverse\n", "der Krümmung der negativen log-Likelihood Funktion am Minimum gegeben.\n", "\n", "Eine Modifikation ist notwendig, wenn es sich um eine verzerrte\n", "Schätzung handelt, bei der eine Verzerrung (ein \"Bias\") \n", "$b := E[\\hat{p}] - p_0$ auftritt, also eine Abweichung des Erwartungswerts \n", "des Schätzwerts vom wahren Wert $p_0$. \n", "Für Likelihood-Schätzungen geht die Verzerrung für große\n", "Stichproben gegen Null. Ein Beispiel hatten wir schon oben bei\n", "der Diskussion der Schätzung der Standardabweichung einer Normalverteilung\n", "gesehen. Da man die Verzerrung im Allgemeinen nicht kennt und sie nur\n", "für kleine Stichproben relevant ist, wird sie in der Praxis weggelassen\n", "und in der Cramér-Rao-Grenze für maximum-Likelihood-Schätzungen das \n", "Gleichheitszeichen angenommen, weil solche Schätzverfahren optimal\n", "effizient sind.\n", "\n", "Unter diesen Annahmen lässt sich das Verfahren zur Schätzung der \n", "Parameterunsicherheiten erweitern, um eine Schätzung der \n", "vollständigen Kovarianzmatrix aller Parameter zu erhalten:\n", "\n", "\\begin{equation}\n", "\\left( {\\bf V}\\right)_{ij} =\n", "\\left( { - \\left. \\frac {{\\partial}^2 \\ln{\\cal{L}}}\n", "{\\partial p_i \\, \\partial p_j} \n", "\\right|}_{\\hat p_i \\hat p_j}\\right)^{-1}\\,.\n", "\\end{equation}\n", "\n", "Diese Methode zur Bestimmung der Parameterunsicherheiten wird in praktisch allen\n", "Applikationen, Programmen und Bibliotheksfunktionen zur Parameteranpassung verwendet. \n", "\n", "Wenn die Zuverlässigkeit der so bestimmten Grenzen für die Parameterwerte\n", "wichtig ist, sollten Ensembletests mit einer großen Anzahl an Pseudodatensätzen \n", "durchgeführt werden, wie sie weiter unten beschrieben werden. \n", "\n", "Die Bestimmung der vollständigen Kovarianzmatrix der Parameterunsicherheiten \n", "ist wichtig, weil Korrelationen zwischen den Parametern auch bei unkorrelierten \n", "Daten sehr häufig auftreten. Wenn die Korrelationen besonders groß sind, kann \n", "dies auch ein Hinweis auf eine ungeschickte Parametrisierung des Problems sein. \n", "Bei Korrelationen größer als 10% sollten sie als Teil des Ergebnisses unbedingt \n", "kommuniziert werden, weil die Form der innerhalb eines vorgegebenen \n", "Konfidenzniveaus noch zulässigen Parameterbereiche davon abhängt. \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Analyse der Parameterunsicherheiten in Beispiel 1\n", "\n", "In unserem Beispiel mit *iminuit* von oben wurde die Analyse der Unsicherheiten\n", "nach der beschriebenen Methode bereits durchgeführt, und wir können uns die \n", "Ergebnisse einfach ausgeben lassen. \n", "Das um entsprechenden Zeilen zur Ausgabe sehen so aus:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# print results of fit with iminuit\n", "print(\"parameter names :\", m.parameters)\n", "print(\"best-fit values :\", tuple(m.values))\n", "print(\"uncertainties :\", tuple(m.errors))\n", "print(\"covariace matrix:\\n\", np.asarray(m.covariance))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Beim Aufruf der Funktion *mFit()* wurde (natürlich) auch gleich eine vollständige \n", "Analyse der Unsicherheiten durchgeführt.\n", "Das um die Ausgabe von Unsicherheiten und Korrelationen ergänzte Beispiel \n", "sieht folgendermaßen aus:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Print results\n", "print(\"\\n*==* user-defined cost: Fit Result:\")\n", "print(\" parameter names: \", fit_results[\"parameter names\"])\n", "print(\" parameter values: \", fit_results[\"parameter values\"])\n", "print(\" parameter uncertainties:\\n\", fit_results[\"confidence intervals\"])\n", "print(\" parameter correlations:\\n\", fit_results[\"correlation matrix\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Bitte beachten, dass die Analyse der Unsicherheiten in diesem Fall nach \n", "einem anderen Verfahren erfolgt ist, das unten beschrieben wird und das \n", "auch mit asymmetrischen Unsicherheiten umgehen kann, die relevant werden,\n", "wenn die log-Likelihood am Minimum stark von einer Parabel abweicht. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Grafische Analyse der negativen Log-Likelihood Funktion aus Beispiel 1\n", "\n", "Für unser erstes Beispiel wollen wir zur Illustration den Verlauf von \n", "$n2l{\\cal L}$ grafisch darstellen. In diesem einfachen Fall mit nur zwei\n", "freien Parametern stellt dies bereits die Lösung des Minimierungsproblems dar. \n", "\n", "In der Codezelle unten ist dafür eine nützliche Hilfsfunktionen definiert.\n", "Die Codierung der Kostenfunktion muss zu deren Verwendung noch angepasst \n", "werden, um eine optimale vektorbasierte Berechnung für Arrays von \n", "Eingabeparametern zu ermöglichen.\n", "Für die Minimierung mit *iminuit* wurde die Kostenfunktion jeweils für\n", "einzelnen Parameterwerte berechnet - zur grafischen Darstellung wollen wir\n", "sie aber effizient in einem einzigen Aufruf auf einem zweidimensionalen \n", "Gitternetz auswerten. \n", "\n", "Die für diesen Fall angepasste Kodierung der Kostenfunktion ist ebenfalls \n", "in der Codezelle unten gezeigt." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from matplotlib import cm\n", "\n", "\n", "def get_3Dfunction_data(\n", " xmin=1.0,\n", " xmax=1.0,\n", " ymin=-1.0,\n", " ymax=1.0,\n", " func=None,\n", " fkwargs={},\n", " delta_x=0.1,\n", " delta_y=0.1,\n", "):\n", " \"\"\"\n", " get function values on 2D-grid\n", "\n", " Args:\n", " float xmin: minimum x-range\n", " float xmax: maximum x-range\n", " float ymin: minimum y-range\n", " float ymys: maximum y-range\n", " funtion func: funtion f(x,y **kwargs)\n", " dict fkwargs: prameters to pass to function\n", " delta: resolution of grid\n", "\n", " Returns:\n", " tuple X, Y, Z\n", " \"\"\"\n", " x = np.arange(xmin, xmax, delta_x)\n", " y = np.arange(ymin, ymax, delta_y)\n", " X, Y = np.meshgrid(x, y)\n", " Z = func(X, Y, **fkwargs)\n", " return X, Y, Z\n", "\n", "\n", "# define efficient version of cost function for arrays mu and sigma:\n", "def plot_myCost(mu, sigma):\n", " # -2*log-likelihood of a 1-d Gauss distribution\n", " # mu, sigma: numpy-arrays of parameter values\n", " n2lL = np.zeros_like(mu)\n", " for i in range(len(data)):\n", " n2lL += (data[i] - mu) * (data[i] - mu)\n", " n2lL /= sigma * sigma\n", " return n2lL + 2.0 * len(data) * np.log(sigma)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Einen recht guten quantitativen Überblick über eine skalare Funktion von zwei\n", "Variablen erhält man, wenn eine Konturdarstellung, analog zu Höhenlinien\n", "in der Geografie, gewählt wird. \n", "In dem für die Darstellung verwendeten Array kann man auch nach dem minimalen \n", "Wert und den zugehörigen Parameterwerten suchen und sie ausgeben.\n", "\n", "Hier der Code dazu:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot cost function near minimum\n", "# - ranges close to minimum\n", "minx = 1.9\n", "maxx = 2.2\n", "miny = 0.4\n", "maxy = 0.64\n", "dx = (maxx - minx) / 125\n", "dy = (maxy - miny) / 125\n", "# get Z data for grid\n", "X, Y, Z = get_3Dfunction_data(\n", " xmin=minx, xmax=maxx, ymin=miny, ymax=maxy, func=plot_myCost, delta_x=dx, delta_y=dy\n", ")\n", "\n", "# find and print x/y of minimum\n", "minz = np.amin(Z)\n", "i_xmn, i_ymn = np.unravel_index(np.argmin(Z), np.shape(Z))\n", "x_opt = X[i_xmn, i_ymn]\n", "y_opt = Y[i_xmn, i_ymn]\n", "print(\"Minimum z =\", minz, \" at (x,y) =(\", x_opt, \",\", y_opt, \")\")\n", "\n", "# create figure ...\n", "fig = plt.figure(figsize=(10.0, 7.5))\n", "# ... plot (filled) contour lines ...\n", "cont = plt.contourf(X, Y, Z, cmap=\"GnBu\", levels=30)\n", "# ... and plot (x,y) of minimum ...\n", "plt.plot(x_opt, y_opt, marker=\"x\", color=\"red\", markersize=12)\n", "plt.vlines(x_opt, miny, maxy, linestyles=\"dashed\", color=\"orange\")\n", "plt.hlines(y_opt, minx, maxx, linestyles=\"dashed\", color=\"orange\")\n", "# ... and add a legend\n", "cbar = plt.colorbar(cont, shrink=0.5, aspect=10, label=\"-2ln(L)\")\n", "xl = plt.xlabel(\"par_x\", size=\"x-large\")\n", "yl = plt.ylabel(\"par_y\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Einen guten Eindruck des Verhaltens am Minimum gewinnt man, wenn man den Verlauf der\n", "Funktion am Minimum entlang der Achsenrichtungen darstellt, die in der Grafik oben als\n", "orange, gestrichelte Linien eingezeichnet sind. Hier der Code für die x-Achse:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(\"x-projection at minimum y\")\n", "\n", "# extract data for plot from X and Z arrays\n", "pxvals = X[i_ymn, :]\n", "proj_n2lLx = Z[i_ymn, :]\n", "\n", "fig_x = plt.plot(pxvals, proj_n2lLx)\n", "plt.xlabel(\"par_x\", size=\"x-large\")\n", "plt.ylabel(\"-2ln(L)\", size=\"x-large\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Und noch einmal fast das gleiche für die y-Achse:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(\"y-projection at minimum x\")\n", "\n", "# extract data for plot from Y and Z arrays\n", "pyvals = Y[:, i_xmn]\n", "proj_n2lLy = Z[:, i_xmn]\n", "\n", "fig_y = plt.plot(pyvals, proj_n2lLy)\n", "plt.xlabel(\"par_y\", size=\"x-large\")\n", "plt.ylabel(\"-2ln(L)\", size=\"x-large\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Man sieht, dass schon in diesem sehr einfachen Beispiel die negative log-Likelihood \n", "nicht für alle Parameter parabelfömig ist, sondern dass dies für den Parameter \n", "*mu=sigma* nur näherungsweise zutrifft. Wenn die Stichprobe vergrößert wird (bitte ausprobieren!),\n", "stimmt die parabolische Näherung in der Nähe des Minimums immer besser, \n", "und das Minimum wird schärfer. Der nichtparabolische Verlauf ist der Grund dafür, \n", "dass die mit *mFit()* oben bestimmte Unsicherheit für diesen Parameter leicht asymmetrisch\n", "war, nämlich $(-0.033, +0.038)$, während die symmetrischen Unsicherheiten als \n", "$\\Delta\\sigma=\\pm 0.035$ ausgegeben wurden. \n", "\n", "Die mit dieser einfachen, auf der feingranularen grafischen Darstellung der \n", "negativen Log-Likelihood-Funktion basierende Methode erhaltenen Ergebnisse \n", "stimmen sehr gut mit den oben mit Hilfe der Fuktion *mFit()* und dem Paket \n", "*iminuit* erhaltenen Resultaten überein.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2.4 Eigenschaften der Maximum-Likilihood-Schätzung\n", "\n", "Zur Charakterisierung von Verfahren zur Parameterschätzung werden einige grundlegende\n", "Eigenschaften definiert. \n", "\n", "Da die Parameter $\\hat{p}$ und ihre Unsicherheiten aus zufälligen Größen, \n", "den Daten, bestimmt werden, sind sie auch Zufallsgrößen mit Erwartungswerten \n", "$E[\\hat{p}]$ und und Varianz $V[\\hat{p}]$. \n", "\n", "Unabdingbar für ein brauchbares Verfahren ist die Forderung von **Konsistenz**,\n", "d.h. dass mindestens im Grenzfall großer Stichproben für den Schätzwert $\\hat{p}$ \n", "eines Parameters $p$ mit \"wahrem Wert\" $p_0$ gelten muss:\n", "\n", " - $\\displaystyle \\lim_{n \\to \\infty} \\hat{p} = p_0$ .\n", " \n", "\n", "Wünschenswert ist es, wenn für beliebige Stichprobengrößen gilt:\n", " \n", " - $E[\\hat{p}] = p_0\\, ;$\n", " diese Eigenchaft nennt man **Erwartungstreue**.\n", " \n", "\n", "Als Maß für die Abweichung des Erwartungswerts einer Schätzung vom wahren Wert dient \n", "die **Verzerrung** (engl. \"Bias\"), \n", "\n", " - $b = E[\\hat{p}] - p_0 = 0\\,; $ \n", " für große Stichproben muss wegen der geforderten Konsistenz gelten \n", " $\\displaystyle \\lim_{n \\to \\infty} b = 0$.\n", "\n", "Eine weitere wünschenswerte Eigenschaft ist die Forderung nach einer möglichst kleinen\n", "Varianz der Schätzung $\\hat{p}$:\n", "\n", " - $V[\\hat{p}_i]$ minimal; Verfahren, die Cramér-Rao-Grenze erreichen nennt man\n", " **effizient** bzw. **optimal**.\n", "\n", "\n", "Für viel Anwendungen, insbesondere bei der Verwendung von Schätzparametern in Hypothesentests,\n", "ist die Angabe eines korrekten Konfidenzbereichs für den Parameter wichtiger\n", "als die genaue Lage des Optimums, es geht in diesen Fällen also um die Schätzung eines\n", "(Konfidenz-)Intervalls. Als Maß nutzt man die **Überdeckungswahrscheinlichkeit** \n", "(engl. \"coverage\"). Dies ist die Wahrscheinlichkeit, mit der das geschätzte \n", "$\\pm 1 \\sigma$-Intervall den wahren Parameterwert enthält:\n", "\n", " - $W(p_0\\in [\\hat{p} - \\sigma_{\\hat{p}}, \\hat{p} +\\sigma_{\\hat{p}} ])\\,$;\n", " für das $\\pm 1 \\sigma$-Intervall sollte dieser Wert nahe bei 0.683 liegen, \n", " dem $\\pm 1 \\sigma$-Quantil der Gaußverteilung.\n", "\n", "In der Praxis brauchbare Verfahren sollten auch möglichst unbeeinflusst von \"Ausreißern\" \n", "oder (leichten) Fehlmodellierungen der Daten sein. Diese Eigenschaft nennt man \n", "**Robustheit**. Sie ist allerdings quantitativ schwer zu fassen. \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Beispiel: Schätzung der Parameter der Gaußverteilung**, $\\hat{\\mu}$ und $\\hat{\\sigma}$\n", "\n", "Beispielhaft sollen diese Eigenschaften anhand des hier diskutierten Beispiels illustriert werden.\n", "\n", "Die Schätzung des Erwartungswerts $\\mu$ ist für beliebig kleine Stichproben \n", "erwartungstreu und effizient:\n", "\n", " - $b_\\hat{\\mu} = E\\left[ \\frac{1}{n_d} \\sum_i^{n_d} x_i \\right] - \\mu =\n", " \\frac{1}{n_d} \\sum_i^{n_d}E[x_i] - \\mu = \\frac{1}{n_d} \\cdot n_d \\mu - \\mu = 0\\,$, \n", " d.h. $\\hat{\\mu}$ ist unverzerrt und damit erwartungstreu. \n", "\n", "\n", " - $V_{\\hat{\\mu}} = V\\left[\\frac{1}{n_d} \\sum_i^{n_d} x_i\\right]= \n", " \\frac{1}{{n_d}^2} \\sum_i^{n_d} V\\left[ x_i \\right]=\n", " \\frac{1}{{n_d}^2} \\sum_i^{n_d} \\sigma^2 = \\frac{\\sigma^2}{n_d}\\,$; \n", " dieses Ergebnis ist identisch zum oben aus der Anwendung der Cramér-Rao-Grenze \n", " bestimmten Wert, $\\hat{\\mu}$ ist also eine effiziente Schätzung. \n", "\n", "Etwas anders sieht es für die Schätzung der Varianz $v=\\sigma^2$ aus. Aus der\n", "grafischen Darstellung und auch aus der analytischen Form der Likelihoodfunktion \n", "ist ersichtlich, dass die negative log-Likelihood eine Asymmetrie um das Minimum\n", "aufweist und von der Parabelform abweicht. \n", "\n", "\n", " - $b_\\hat{v} = E[\\hat{v^2}] - \\sigma^2 =\n", " E\\left[ \\frac{1}{n_d} \\sum_i^{n_d} {x_i}^2 - {\\hat{\\mu}}^2 \\right ] - \\sigma^2 = \n", " \\frac{1}{n_d} \\sum_i^{n_d} E[{x_i}^2] - {\\hat{\\mu}}^2 - \\sigma^2 $. \n", " Durch intelligentes Einsetzen einer Null, $\\mu^2 - \\mu^2$, erhält man \n", " $b_{\\hat{v^2}} = \n", " \\frac{1}{n_d} (\\sum_i^{n_d} E[{x_i}^2] -\\mu^2) - ({\\hat{\\mu}}^2 - \\mu ^2) - \\sigma^2 =\n", " \\sigma^2 - {\\sigma_{\\hat{\\mu}}}^2 - \\sigma^2 = \n", " {\\sigma_{\\hat{\\mu}}}^2 = \\frac{\\sigma^2}{n_d} \\ne 0 \\,. $\n", " \n", " Diese Schätzung ist also verzerrt! Allerdings geht die Verzerrung für große Stichproben \n", " $n_d \\to \\infty$ gegen Null, die Schätzung ist also \"asymptotisch unverzerrt\". \n", " Oft ist die Verzerrung unerwünscht, und deshalb wendet man die sog. \"Bessel-Korrektur\" an, \n", " um eine unverzerrte Schätzung der Varianz einer Stichprobe zu erhalten:\n", " ${\\hat{\\sigma}^*}^2 = \\displaystyle \\frac{n_d}{n_d-1} \\, {\\hat{\\sigma}}^2 = \n", " \\frac{1}{n_d-1} \\sum_i^{n_d} {x_i}^2 -\\hat{\\mu}^2 \\,$.\n", " \n", " - Als Cramér-Rao-Grenze für die Varianz von $v=\\sigma^2$ erhält man unter Verwendung der\n", " oben angegebenen Formel, die die Verzerrung nicht berücksichtigt: \n", " $V_{\\hat{v}} \\ge \\frac{\\sigma^4}{n_d/2}\\,$. \n", " Eine direkte Berechnung der Varianz aus der oben erhaltenen Formel für den\n", " Schätzwert von $\\sigma^2$, $\\hat{v}=\\frac{1}{n_d} \\sum_i^{n_d} {x_i}^2 - {\\hat{\\mu}}^2$ ergibt:\n", " \n", " $V_{\\hat{v}} = \\frac{n_d-1}{n_d}\\frac{\\sigma^4}{n_d/2}\\,$; \n", " dieser Wert ist sogar kleiner als der Grenzwert. Der Grund dafür ist die \n", " unberücksichtigt gebliebene Verzerrung.\n", " \n", " Aus der **vollständigen Formel für die Cramér-Rao-Bedingung** unter Berücksichtigung der\n", " Verzerrung $b$,\n", " ${\\rm Var}_{\\hat p} = {\\sigma_\\hat{p}}^2 \\ge \n", " \\left[ { - \\left. \\frac { {\\partial}^2 \\ln{\\cal{L} }} \n", " {{\\partial p}^2 } \\right|}_{\\hat p}\\right]^{-1} \\cdot\n", " \\left[ 1 + \\left. \\frac{\\partial b}{\\partial p}\\right|_\\hat{p} \\right]^2 $ \n", " ergibt sich allerdings \n", " $V_{\\hat{v}} \\ge \\left( \\frac{n_d-1}{{n_d}}\\right)^2 \\frac{\\sigma^4}{n_d/2}\\,$; \n", " d.h. die Varianz der Schätzung ist also tatsächlich größer als durch die (korrekt angewendete)\n", " Cramér-Rao-Grenze gegeben. Weil die Grenze erst für unendlich großer Stichproben erreicht\n", " wird, nennt man die Schätzung \"asymptotisch effizient\".\n", "\n", " !!! Übrigens: Programme zur Anpassung basieren ausnahmslos auf der einfachen Variante der \n", " Cramér-Rao-Grenze und geben in diesem Fall also eine zu große Varianz aus! \n", " \n", "**Resumé:**\n", "Die schon in diesem Beispiel sichtbaren Eigenschaften von Maximum-Likelihood-Schätzungen \n", "sind verallgemeinerbar. In einigen Fällen sind ML-Schätzungen effizient und unverzerrt, \n", "im Grenzfall großer Stichproben sind sie generell \"asymptotisch effizient\" und \"asymptotisch\n", "unverzerrt\". Für Verteilungen aus der \"Exponentialfamilie\" werden diese Eigenschaften\n", "schon unabhängig von der Stichprobengröße erreicht. \n", "\n", "**Literaturempfehlung:**\n", "Eine sehr lesenswerte Übersicht über die mathematischen Grundlagen findet sich z.B.\n", "unter dem Link \n", "[Statistical Methods: Likelihood, Bayes and Regression\n", "](http://strimmerlab.org/publications/lecture-notes/MATH20802/).\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2.5 Intervallschätzung\n", "\n", "Wenn die negative log-Likelihood einen parabelförmigen Verlauf hat, ist es ausreichend,\n", "nur die Schätzung des optimalen Werts eines Parameters $\\hat{p}$ und seiner Varianz $V_p$\n", "anzugeben. Dann wird als Konfidenzintervall für den Parameterwert ein symmetrischer Bereich\n", "der Breite $2\\cdot\\sqrt{V_p}$ um den Schätzwert angenommen. \n", "ver\n", "Bei kleinen Stichproben ist der Verlauf am Minimum aber häufig nicht mehr gut durch \n", "eine Parabel anzunähern, und die Angabe eines aus den zweiten Ableitungen am Punkt \n", "des Minimums berechneten symmetrischen Intervalls für die Unsicherheiten ist dann nicht\n", "mehr ausreichend. In solchen Fällen verwendet man zur Abschätzung des Vertrauensintervalls\n", "der Parameter die sogenannte **Profil-Likelihood** und leitet daraus ein Vertrauensintervall\n", "mit vorgegebenem Konfidenzniveau ab - üblich ist es dabei, den $\\pm 1 \\sigma$-Bereich einer\n", "Gaußverteilung zu Grunde zu legen, also ein Konfidenzniveau von 68.3%. Dies entspricht\n", "der Wahrscheinlichkeit, einen Wert einer Gauß-verteilten Größe $x$ innerhalb des \n", "Intervalls $[\\mu-\\sigma, \\mu+\\sigma]$ zu finden. \n", "\n", "### Bestimmung der Parameterunsicherheiten mit der Profil-Likelihood\n", "\n", "Zur Bestimmung der Profil-Likelihood wird für einen vorgegebenen Parameter $p_i$\n", "eine Anzahl von dicht in der Nähe des Minimums liegenden Punkten betrachtet und \n", "die negative log-Likelihood bezüglich aller übrigen Parameter $p_j, j\\ne i$, minimiert. \n", "Auf diese Weise werden die Korrelationen von $p_i$ mit allen anderen Parametern $p_j$\n", "berücksichtigt, und die Porfil-Likelihood ist eine Funktion nur eines einzigen Parameters,\n", "${\\cal L}_p(\\vec x, p_i)$.\n", "\n", "Das Vertrauensintervall, das dem Bereich plus/minus einer Standardabweichung einer\n", "Gaußverteilung um den Parameterwert, also einem Konfidenzniveau von 68.3% entspricht, \n", "erhält man an den Stellen, an denen die negative log-Likelihood um den Wert\n", "$\\frac{1}{2}$ größer ist als am Minimum (bzw. um Eins größer, wenn \n", "man $-2\\ln(\\cal L)$ verwendet). \n", "\n", "Für parabolische Verläufe, wenn der negative Logarithmus der Profil-Likeliood\n", "also die Form $0.5\\cdot ((p_i - \\hat{p}_i)/\\sigma_{p_i})^2$ hat, ist das \n", "Verfahren zur auf den zweiten Ableitungen basierenden Methode äquivalent. \n", "Ein Vorteil der Profil-Likelihood wird offensichtlich, wenn man Transformationen \n", "von Parametern betrachtet, also z.B. \n", "$nl{\\cal L}(\\vec{x};p) \\to nl{\\cal L}(\\vec{x};\\tilde{p})$ \n", "mit $\\tilde{p} = \\tilde{p}(p)$. \n", "Hat man für den Parameter $p$ einen Wert $p_{\\Delta}$ bestimmt bei dem $nl{\\cal L}$ \n", "um den Wert $\\Delta$ ansteigt, so ergibt sich der Wert für $\\tilde{p}$ durch \n", "einfaches Einsetzen: $\\tilde{p}_{\\Delta} = \\tilde{p}(p_{\\Delta})$. \n", "Da man durch eine geeignete Transformation $p \\to {p_\\cal N} = p_{\\cal N}(p)$ \n", "beliebige Verteilungsdichten auf die Gaußform transformieren kann, \n", "lässt sich das mit der Methode bestimmte Konfidenzintervall als das \n", "äquivalente $\\pm 1 \\sigma$-Konfidenzintervall einer Gaußverteilung verstehen. \n", "\n", "Die durch die zweiten Ableitungen am Optimum festgelegte Varianz bzw. die \n", "Kovarianzmatrix der Parameter hat diese vorteilhafte Transformationseigenschaft\n", "übrigens nicht, denn eine Parametertransformation führt zum Auftreten von Ableitungen\n", "$\\partial \\tilde{p}/\\partial p$; man bestimmt so also die Varianz der Parameter,\n", "kennt aber die genaue Verteilungsdichte nicht! \n", "Der Scan der Profil-Likelihood liefert dagegen die Schätzung eines Konfidenzbereichs \n", "für die Parameterwerte - eine **Intervallschätzung**.\n", "\n", "Durch Variation des Wertes von $\\Delta$ lassen sich mit der gleichen Methode\n", "auch andere Konfidenzintervalle bestimmen: \n", "\n", "| | ΔnlL | Δn2lL |\n", "| --: | :--: | :--: |\n", "| 1σ | 0.5 | 1 |\n", "| 2σ | 2 | 4 |\n", "| 3σ | 4.5 | 9 |\n", "| zσ | z²/2 | z² |\n", "\n", "\n", "Den Wert $z$ nennt man - in Anlehnung an den $p$-Wert eines Hypothesentests -\n", "auch \"$z$-Value\". Er drückt ein Konfidenzintervall als Anzahl der äquivalenten\n", "$\\sigma$s einer Gaußverteilung aus. \n", "\n", "Bei nicht-parabolischen Verläufen ergibt sich aus der Analyse der Profil-Likelihood \n", "ein asymmetrisches Konfidenzintervall um den Parameterwert am Minimum, \n", "das in solchen Fällen an Stelle der aus den zweiten Ableitungen am Minimum bestimmen \n", "Schätzungen der Standardabweichungen der Parameter angegeben werden sollte.\n", "\n", "Unter Verwendung des Bayes'schen Theorems und Vorgabe einer A-priori-Wahrscheinlichkeit \n", "${\\cal P}(p_i)$ für den jeweils interessierenden Parameter $p_i$ kann aus der \n", "Profil-Likelihood eine Wahrscheinlichkeitsdichte abgeleitet werden:\n", "\n", "$\\displaystyle {\\cal P}( p_i | \\vec x) \\propto \n", " { {\\cal L}_p(\\vec x | p_i) \\, \\cdot \\, {\\cal P}_p(p_i)} \\, $.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Grafische Darstellung der Profil-Likelihood zu Beispiel 1\n", "\n", "Profil-Likelihoods lassen sich für unser Beispiel mit zwei Parametern direkt aus dem schon\n", "für die grafische Darstellung genutzten Array erzeugen. Dazu muss der Code zur Darstellung\n", "der Projektionen der negativen Log-Likelihood nur leicht ergänzt werden. Anstelle der Position\n", "des globalen Minimums in $x$ bzw. $y$ wird nun das Minimum bzgl. $y$ bzw. $x$ für vorgegebene\n", "Werte $x$ bzw. $y$ bestimmt. \n", "\n", "Zur Analyse der Profil-Likelihood ist es zweckmäßig, die an den Stützstellen der Arrays \n", "X und Y bestimmten Werte mit Hilfe von kubischen Splines (*scipy.interpolate.CubicSpline*) \n", "zu interpolieren, die zweimal stetig differenzierbar sind. \n", "Im Code unten ist dazu eine allgemein verwendbare Funktion *analyze_n2lLp()* definiert, \n", "die eine eindimensionale Profil-Likelihood auswertet und grafisch darstellt. \n", "Zunächst wird die dazu an den Stützstellen bestimme Profil-Likelihood durch einen kubischen \n", "Spline repräsentiert. Bestimmt werden daraus der optimale Parameterwert durch Nullsetzen \n", "der ersten Ableitung und die symmetrische Unsicherheit aus der zweiten Ableitung am\n", "Optimalwert.\n", "Zeichnet man die entsprechenden horizontale Linie bei $-\\Delta 2\\ln{\\cal L}_{prof} = 1$ \n", "ein und projiziert die Schnittpunkte mit der Profil-Likelihood auf die $x$-Achse, erhält man \n", "das 68%-Konfidenzintervall für die Parameterunsicherheit. Die Schnittpunkte werden (natürlich) \n", "numerisch mit Hilfe des Newton-Verfahrens bestimmt. Üblicherweise zieht man bei der grafischen \n", "Darstellung den Minimalwert ab, weil nur die Änderung im Bezug auf das Minimum relevant ist. \n", "\n", "Hier der Code für den Logarithmus der Profil-Likelihood des in $x$-Richtung \n", "aufgetragenen Parameters, d.h. den Erwartungswert $\\mu$ der Normalverteilung:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# extract data Z array\n", "prof_n2lLx = np.zeros_like(pxvals)\n", "for i, x in enumerate(pxvals):\n", " prof_n2lLx[i] = np.amin(Z[:, i])\n", "\n", "\n", "def analyze_n2lLp(pvals, n2plL, name=\"p\"):\n", " from scipy.interpolate import CubicSpline\n", " from scipy.optimize import newton\n", "\n", " \"\"\"Analyze profile likelihood\n", "\n", " Input:\n", " - parameter values \n", " - corresponding values of -2 *log(L_p) \n", " - parameter name \n", " \n", " Returns:\n", " - popt : optimum (= min of profile likelihood) \n", " - d_popt: uncertainty from 2nd derivative\n", " - d_pm: negative uncertainty\n", " - d_pp: positve uncertainty from profile likelihood scan\n", " \"\"\"\n", "\n", " minp = min(pvals)\n", " maxp = max(pvals)\n", " # define interpolating cubic Spline\n", " n2plL = CubicSpline(pvals, n2plL)\n", "\n", " # calculate derivatives\n", " deriv = n2plL.derivative(nu=1)\n", " deriv2 = n2plL.derivative(nu=2)\n", "\n", " # find minimum using Newton method on 1st derivative\n", " popt = newton(deriv, x0=minp, x1=maxp)\n", " min_n2plL = n2plL(popt)\n", " # determin uncertainty from Cramer-Rao\n", " errdef = 1.0 # for 2negLogL, 0.5 else\n", " d_popt = np.sqrt(2.0 * errdef / deriv2(popt))\n", "\n", " # determine confidence interval of parameter from profile likelihood\n", " def f(x):\n", " return n2plL(x) - min_n2plL - errdef\n", "\n", " p_neg = newton(f, x0=minp, x1=popt)\n", " p_pos = newton(f, x0=popt, x1=maxp)\n", " d_pm = popt - p_neg\n", " d_pp = p_pos - popt\n", "\n", " # output\n", " print(\"*==* Analysis of Profile Likelihood for parameter\", name)\n", " print(\n", " \"optimum: {:.4g}+/-{:.2g} [-{:.2g}, +{:.2g}]\".format(popt, d_popt, d_pm, d_pp)\n", " )\n", "\n", " # plot profile likelihood\n", " fig_p = plt.plot(pvals, n2plL(pvals) - min_n2plL)\n", " # plot lines for uncertainty analysis\n", " plt.hlines(1.0, minp, maxp, color=\"orange\", linestyle=\"dashed\")\n", " plt.arrow(p_neg, 1.0, 0.0, -1.0, head_length=0.33, color=\"orange\")\n", " plt.arrow(p_pos, 1.0, 0.0, -1.0, head_length=0.33, color=\"orange\")\n", "\n", " plt.title(\"Profile Likelihood\")\n", " plt.xlabel(name, size=\"x-large\")\n", " plt.ylabel(\"-2$\\Delta$pln(L)\", size=\"x-large\")\n", "\n", " return popt, d_popt, d_pm, d_pp\n", "\n", "\n", "val, dval, dm, dp = analyze_n2lLp(pxvals, prof_n2lLx, name=r\"$\\mu$\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Die Kurve sieht in diesem Fall fast genau so aus wie die einfache Projektion.\n", "Das liegt daran, dass die Korrelation zwischen den Parametern klein ist und \n", "deshalb eine Änderung des $x$-Parameters kaum Einfluss auf den $y$-Parameter hat.\n", "\n", "Hier noch der Code für die Standardabweichung $\\sigma$ der Normalverteilung, die\n", "wir in der grafischen Darstellung in $y$-Richtung aufgetragen hatten. \n", "Zur Analyse und grafischen Darstellung der Profil-Likelihood verwenden wir wieder \n", "die eben schon definierte Funktion *analyze_n2lLp()*:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# extract nlL values from Z array\n", "prof_n2lLy = np.zeros_like(pyvals)\n", "for i, y in enumerate(pyvals):\n", " prof_n2lLy[i] = np.amin(Z[i, :])\n", "\n", "val, dval, dm, dp = analyze_n2lLp(pyvals, prof_n2lLy, name=r\"$\\sigma$\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Diese der Illustration des Verfahrens dienende Vorhehensweise funktioniert\n", "nur für einfache Problemstellungen mit maximal zwei Parametern und ist numerisch\n", "nicht effizient, wenn gar keine grafische Darstellung benötigt wird.\n", "\n", "Wenn mehr Parameter im Spiel sind oder die Effizienz der Berechnungen von Bedeutung\n", "ist, sollten *iminuit* oder, komfortabler, die Interface-Funktionen aus den Paketen \n", "*phyFit* oder *kafe2* verwendet werden.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Profil-Likelihood in zwei Dimensionen: Konfidenzkonturen\n", "\n", "Das Prinzip der Profil-Likelihood lässt sich auch in zwei Dimensionen,\n", "also für Paare von Parametern, anwenden. Man erhält dann Konturlinien,\n", "die Konfidenzbereichen analog zu ein- oder auch zwei-$\\sigma$-Konturen von \n", "Gauß-förmigen Verteilungen entsprechen. Bei starker Abweichung der Konturen \n", "von der Ellipsenform sollten sie zusammen mit dem Ergebnis berichtet werden.\n", "\n", "Da das numerische Verfahren aufwändig ist, sollten zur Darstellung von \n", "Konfidenzkonturen die Implementierungen in den Paketen *iminiut* oder\n", "*kafe2* verwendet werden. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Profil-Likelihood und Konfidenzkonturen mit Hilfe der Funktion *phyFit.mFit()*\n", "\n", "Profile-Likelihood-Kurven und Konfidenzkonturen können leicht auch mit der \n", "Funktion *mFit()* erzeugt werden, wenn die zusätzliche Option `plot_cor=True` \n", "angegeben wird. Ansonsten ist der Code identisch zum bereits oben besprochenen\n", "ersten Beispiel mit *mFit()*:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fit_results = mFit(myCost, plot_cor=True)\n", "\n", "# Print results\n", "pvals, perrs, cor, gof, pnams = fit_results.values()\n", "print(\"\\n*==* user-defined cost: Fit Result:\")\n", "print(\" parameter names: \", pnams)\n", "print(\" parameter values: \", pvals)\n", "print(\" neg. parameter errors: \", perrs[:, 0])\n", "print(\" pos. parameter errors: \", perrs[:, 1])\n", "print(\" correlations : \\n\", cor)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In der mit dem Code oben erzeugten Matrix aus Grafiken werden auf der Diagonalen die\n", "Verläufe der Profil-Likelihood angezeigt; in den unteren Nebendiagonalelementen sind\n", "die ein- und zwei-$\\sigma$ Konturen dargestellt. Die Fehlerbalken sind die aus der Analyse\n", "der zweiten Ableitungen gewonnenen Unsicherheiten, die für $\\Delta\\mu$ perfekt und \n", "für $\\Delta\\sigma$ noch recht gut mit den mit Hilfe der Profil-Likelihood bestimmten Konfidenzintervallen übereinstimmen. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Direkte Anpassung einer Verteilungsdichte an Daten mit der Funktion *phyFit.mFit()*\n", "\n", "Die Funktion *mFit()* bietet ein weiters Interface zur Durchführung von Anpassungen mit der ML-Methode an, bei denen als Eingaben die Daten die Verteilungsdichte angegeben werden. Damit kann *mFit()* auch die Darstellung der Daten und der angepassten Verteilungsdichte inklusive eines Konfidenzbands um die Modellfunktion übernehmen (Option `plot=True`).\n", "Die zusätzliche Option `plot_cor=True` ist ebenfalls verfügbar und erzeugt die gleiche\n", "Ausgabe wie im vorigen Beispiel. Zusätzlich sind im Beispiel in der Codezelle unten noch einige weitere Optionen zur Beschriftung der Achsen und der Legende angegeben. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Definition der Gaußvereteilung\n", "def fGauss(x, mu=10.0, sigma=0.5):\n", " return np.exp(-((x - mu) ** 2) / 2.0 / sigma**2) / np.sqrt(2.0 * np.pi) / sigma\n", "\n", "\n", "# Anpassung Verteilungsdichte an Daten\n", "rdict = mFit(\n", " fGauss,\n", " data=data, # data - if not None, normalised PDF is assumed as model\n", " plot=True, # plot data and model\n", " # plot_cor=True, # plot profiles likelihood and contours\n", " axis_labels=[\"Daten x\", \"Dichte pdf(x; *p)\"],\n", " data_legend=\"Daten\",\n", " model_legend=\"Gaußverteilung\",\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Beispiel 2: log-Likelihood Anpassung einer Exponentialfunktion\n", "\n", "Wir wollen nun ein etwas komplexeres, realistisches Beispiel zur Anpassung von\n", "Parametern an Daten mit Hilfe des Maximum-Likelihood Verfahrens anschauen. \n", "\n", "Dazu verwenden wir echte Messungen der Lebensdauern von gestoppten \n", "kosmischen Myonen, die in einem Wasser-Cherenkov-Detektor (in diesem Fall \n", "einer Kaffeekanne mit Photomultiplier) nachgewiesen wurden. Beim Durchgang \n", "durch das Wasser erzeugen die Myonen ein Lichtsignal; wenn sie im Wasser\n", "oder im Boden darunter gestoppt werden, können Elektronen aus den zerfallenden\n", "Myonen wieder in den Detektor gelangen und ebenfalls nachgewiesen werden. \n", "Die Daten entsprechen den Zeitdifferenzen (in µs) solcher Doppelpulse. Zufällig\n", "eintreffende Myonen, Umgebungsstrahlung und evtl. Detektorrauschen erzeugen \n", "ebenfalls Doppelpulse, die im Gegensatz zu den Lebensdauern aber flach verteilt\n", "sind und einen Untergrund bei der Lebensdauermessung darstellen.\n", "\n", "Die relevante, auf Eins normierte Verteilungsdichte und die Daten finden \n", "sich in der folgenden Code-Zelle:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def exponentialDecayPDF(t, tau=2.0, fbg=0.2, a=1.0, b=11.5):\n", " \"\"\"Probability density function\n", "\n", " for an exponential decay with flat background. The pdf is normed for the interval\n", " [a=1µs, b=11.5µs); these parameters a and b must be fixed in the fit!\n", " \"\"\"\n", "\n", " # 1.) Exponential distribution of life time measuremtens in the intervall [a, b]\n", " pdf1 = np.exp(-t / tau) / tau / (np.exp(-a / tau) - np.exp(-b / tau))\n", "\n", " # 2. flat distribution in the interval [a, b]\n", " pdf2 = 1.0 / (b - a)\n", "\n", " # the full distribution is a combination of both, with a background fraction fbg\n", " return (1 - fbg) * pdf1 + fbg * pdf2\n", "\n", "\n", "# real data from measurement with a Water Cherenkov detector (\"Kamiokanne\")\n", "dT = [7.42, 3.773, 5.968, 4.924, 1.468, 4.664, 1.745, 2.144, 3.836, 3.132, 1.568, 2.352,\n", " 2.132, 9.381, 1.484, 1.181, 5.004, 3.06, 4.582, 2.076, 1.88, 1.337, 3.092, 2.265,\n", " 1.208, 2.753, 4.457, 3.499, 8.192, 5.101, 1.572, 5.152, 4.181, 3.52, 1.344, 10.29,\n", " 1.152, 2.348, 2.228, 2.172, 7.448, 1.108, 4.344, 2.042, 5.088, 1.02, 1.051, 1.987,\n", " 1.935, 3.773, 4.092, 1.628, 1.688, 4.502, 4.687, 6.755, 2.56, 1.208, 2.649, 1.012,\n", " 1.73, 2.164, 1.728, 4.646, 2.916, 1.101, 2.54, 1.02, 1.176, 4.716, 9.671, 1.692,\n", " 9.292, 10.72, 2.164, 2.084, 2.616, 1.584, 5.236, 3.663, 3.624, 1.051, 1.544, 1.496,\n", " 1.883, 1.92, 5.968, 5.89, 2.896, 2.76, 1.475, 2.644, 3.6, 5.324, 8.361, 3.052,\n", " 7.703, 3.83, 1.444, 1.343, 4.736, 8.7, 6.192, 5.796, 1.4, 3.392, 7.808, 6.344,\n", " 1.884, 2.332, 1.76, 4.344, 2.988, 7.44, 5.804, 9.5, 9.904, 3.196, 3.012, 6.056,\n", " 6.328, 9.064, 3.068, 9.352, 1.936, 1.08, 1.984, 1.792, 9.384, 10.15, 4.756, 1.52,\n", " 3.912, 1.712, 10.57, 5.304, 2.968, 9.632, 7.116, 1.212, 8.532, 3.000, 4.792, 2.512,\n", " 1.352, 2.168, 4.344, 1.316, 1.468, 1.152, 6.024, 3.272, 4.96, 10.16, 2.14, 2.856,\n", " 10.01, 1.232, 2.668, 9.176]\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Die Anpassung mit Hilfe des log-Likelihood-Verfahrens funktioniert genau so \n", "wie oben. Dieses Mal werden allerdings an Stelle der Kostenfunktion die Daten\n", "und die Verteilungsdichte sowie weitere Optionen zur grafischen Darstellung \n", "als Parameter angegeben. Damit die Anpassung funktioniert, muss das\n", "Anpassungswerkzeug einige besondere Optionen unterstützen: \n", "- die Begrenzung von Parametern auf einen sinnvollen Bereich \n", " (Option `limits=('fbg', 0., 1.)`) und \n", "- die Fixierung von Parametern, die in der Anpassung nicht variiert werden. \n", " In diesem Fall ist dies das Intervall, in dem die Messungen verlässlich sind\n", " und auf das die Verteilungsdichte normiert ist \n", " (Option `fixPars = ['a', 'b']`).\n", " \n", "Hier nun der vollständige Programmcode zur Anpassung einer Verteilungsdichte\n", "an (ungebinnte) Daten:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "results = mFit(\n", " exponentialDecayPDF,\n", " data=dT, # data - if not None, a normalised PDF is assumed as model\n", " limits=(\"fbg\", 0.0, 1.0), # parameter limits\n", " fixPars=[\"a\", \"b\"], # fix parameter(s)\n", " neg2logL=True, # use -2 * ln(L)\n", " plot=True, # plot data and model\n", " plot_band=True, # plot model confidence-band\n", " plot_cor=True, # plot profiles likelihood and contours\n", " axis_labels=[\n", " \"life time \" + \"$\\Delta$t ($\\mu$s)\",\n", " \"Probability Density pdf($\\Delta$t; *p)\",\n", " ],\n", " data_legend=\"$\\mu$ lifetime data\",\n", " model_legend=\"exponential decay + flat background\",\n", ")\n", "\n", "# Print results\n", "pvals, perrs, cor, gof, pnams = results.values()\n", "print(\"\\n*==* unbinned ML Fit Result:\")\n", "print(\" parameter names: \", pnams)\n", "print(\" parameter values: \", pvals)\n", "print(\" neg. parameter errors: \", perrs[:, 0])\n", "print(\" pos. parameter errors: \", perrs[:, 1])\n", "print(\" correlations : \\n\", cor)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Anders als bei der Anpassung der Normalverteilung an Gauß-verteilte Daten zeigen die\n", "Profil-Likelihoods starke Abweichungen von der Parabelform, und die Kontur weicht stark\n", "von der Form einer Ellipse ab. \n", "In einem solchen Fall ist die Bereitstellung der Grafiken zusätzlich zur Angabe der numerischen Ergebnisse notwendig. Auch die stark asymmetrischen Unsicherheiten müssen angegeben werden. \n", "\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. Anpassung von Verteilungsdichten an Histogramme\n", "***\n", "\n", "Ein weiteres typisches Beispiel für die Anwendung der log-Likelihood-Methode\n", "ist die Behandlung von Problemen, bei denen Poisson-verteilte Größen\n", "auftreten. Dies sind z.$\\,$B. die Anzahlen von Einträgen in einzelnen\n", "Intervallen von Häufigkeitsverteilungen (Histogramme).\n", "\n", "Die Poissonverteilung von Anzahlen $n_i$ mit Erwartungswerten $\\mu_i$\n", "ist gegeben durch\n", "\n", "\\begin{equation}\\label{equ-Poisson}\n", " P(n_i;\\mu_i)=\\frac{{\\mu_i}^{n_i}} {n_i\\,!} \\, {\\rm e}^{-\\mu_i} \\,.\n", "\\end{equation}\n", "\n", "Durch Bilden des Produkts über alle (als statistisch unabhängig\n", "angenommenen) $n_b$ Bins eines Histogramms erhält man die Likelihood \n", "\n", "\\begin{equation}\\label{equ-LPoisson}\n", " {\\cal{L}}_{Poisson}=\\displaystyle\\prod_{i=1}^{n_b}\n", " {\\rm P}\\left(n_i;\\mu_i(\\vec p)\\right)\\, ,\n", "\\end{equation}\n", "\n", "und schließlich durch Logarithmieren und Weglassen des nicht von den \n", "Parametern abhängigen Terms $\\ln(n!)$ die negative log-Likelihood der\n", "Poisson-verteilten Anzahlen der Einträge in den Bins eines Histogramms\n", "\\begin{equation}\\label{equ-nlLPoisson}\n", "nl{\\cal{L}}_{Poisson}= -\\ln {\\cal{L}}_{Poisson} =\n", "\\displaystyle\\sum_{i=1}^{n_b} - n_i \\cdot \\ln(\\mu_i(\\vec p))\\,+\\,\\mu_i(\\vec p) \\,.\n", "\\end{equation}\n", "\n", "Bei Verwendung dieser Kostenfunktion im Minimierungsprozess können\n", "Anpassungen von Verteilungsdichten an Histogramme ganz analog wie bei\n", "der Verwendung der oben bereits diskutierten Likelihood zur Anpassung\n", "von Verteilungsdichten an Messdaten durchgeführt werden, die in diesem\n", "Fall aber als Histogramm vorliegen. \n", "Der Rechenaufwand hängt von der Zahl der Bins und nicht - wie im ersten Beispiel - \n", "von der Größe des Datensatzes ab. \n", "\n", "Die Anzahl der auf Grund der Verteilungsdichte erwarteten Einträge $n_i$ \n", "bestimmt man sinnvollerweise durch Bildung des Integrals der Verteilung\n", "über die Binbreite - eine numerische Näherung (z.B. Simpson 2. Ordnung)\n", "ist meist ausreichend. \n", "\n", "In manchen Fällen wird auch eine Gauß-Näherung der Poisson-Verteilung\n", "verwendet, $P(n_i;\\mu_i) = {\\cal N}(n_i; \\mu_i, \\sigma=\\sqrt{\\mu_i})$.\n", "Bei großen Anzahlen an Bin-Einträgen ist diese Näherung akzeptabel und\n", "erlaubt häufig die Anwendung von Programmen zur Anpassung, die auf\n", "der Methode der kleinsten Fehlerquadrate beruhen und daher nur eine\n", "auf der Gaußverteilung beruhende Kostenfunktion erlauben. \n", "\n", "Wichtig ist es, die Schätzung der Unsicherheiten $\\sigma_i = \\sqrt{\\mu_i}$\n", "aus der Modellerwartung und nicht aus den beobachteten Werten in den\n", "Daten zu bestimmen, da letzteres zu einer verzerrten Parameterschätzung\n", "führen würde und Bins mit Null Einträgen sogar weggelassen werden müssten. \n", "Im Notfall kann man sich behelfen, indem zunächst eine Anpassung mit den aus\n", "der Beobachtung gewonnenen Unsicherheiten durchgeführt wird, die in \n", "einer wiederholten Anpassung durch die aus der dann näherungsweise bekannten \n", "Modellerwartung gewonnenen Werte ersetzt werden.\n", "Bei den empfohlenen Werkzeugen *kafe2* und *phyFit* sind solche Umstände\n", "allerdings nicht notwendig, da sie als Voreinstellung das korrekte, oben \n", "skizzierte Likelihood-Verfahren für die Anpassung an Histogrammdaten nutzen. \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.1 Beispiel: Anpassung an Histogramm-Daten\n", "\n", "Die Anpassung eines Modells an histogrammierte Daten wird durch die Funktion\n", "*hFit()* aus dem Paket *PhyPraKit.phyFit* unterstützt. Die Vorgehensweise ist ganz\n", "analog zur Anpassung von ungebinnten Daten. Von Anwenderseite müssen lediglich\n", "die Daten und die Verteilungsdichte bereit gestellt sowie Angaben zu den \n", "Ausgabeoptionen gemacht werden. Die passenden Kostenfunktionen, also das zweifache\n", "des negativen natürlichen Logarithmus der Poissonverteilung - oder, oft in \n", "ausreichender Näherung der Gaußverteilung - sind in *hFit()* implementiert. \n", "\n", "Hier das vollständige Programmbeispiel:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from PhyPraKit.phyFit import hFit\n", "\n", "\n", "# the model function to fit\n", "def SplusB_model(x, mu=6.987, sigma=0.5, s=0.2):\n", " \"\"\"(normalized) pdf of a Gaussian signal on top of flat background\"\"\"\n", " normal = np.exp(-0.5 * ((x - mu) / sigma) ** 2) / np.sqrt(2.0 * np.pi * sigma**2)\n", " linear = 1.0 / (xmx - xmn)\n", " return s * normal + (1 - s) * linear\n", "\n", "\n", "# the histogram data\n", "nbins = 40\n", "xmn = 1\n", "xmx = 10\n", "bedges = np.linspace(xmn, xmx, nbins + 1)\n", "bcontents = np.array([1, 1, 1, 2, 2, 2, 4, 1, 0, 3, 1, 1, 0, 2, 3, 3, 1, 1, 0, 2,\n", " 3, 2, 3, 1, 1, 8, 6, 7, 9, 1, 0, 1, 2, 1, 3, 2, 1, 3, 2, 4])\n", "#\n", "# --- perform fit\n", "#\n", "hFit_results = hFit(\n", " SplusB_model,\n", " bcontents,\n", " bedges, # bin entries and bin edges\n", " limits=[(\"s\", 0.0, None), (\"sigma\", 0.0, None)], # parameter limits\n", " plot=True, # plot data and model\n", " plot_band=True, # plot model confidence-band\n", " plot_cor=True, # plot profiles likelihood and contours\n", " axis_labels=[\"x\", \"y \\ f(x, *par)\"],\n", " data_legend=\"peudo-data\",\n", " model_legend=\"signal + background model\",\n", ")\n", "\n", "# Print results\n", "# pvals, perrs, cor, chi2, pnams = hFit_results\n", "print(\"\\n*==* histogram fit Result:\")\n", "print(\" parameter names: \", hFit_results[\"parameter names\"])\n", "print(\" goodness-of-fit: {:.3g}\".format(hFit_results[\"goodness-of-fit\"]))\n", "print(\" parameter values: \", hFit_results[\"parameter values\"])\n", "print(\" neg. parameter errors: \", hFit_results[\"confidence intervals\"][:, 0])\n", "print(\" pos. parameter errors: \", hFit_results[\"confidence intervals\"][:, 1])\n", "print(\" correlations : \\n\", hFit_results[\"correlation matrix\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Die grafische Ausgabe zeigt die angepasste Verteilungsdichte sowie deren Unsicherheit, die durch \n", "Propagation der Parameterunsicherheiten auf die Modellvorhersage bestimmt wird. \n", "Die rechteckigen Flächen repräsentieren die Zahl der in den Daten beobachteten Einträge und \n", "entsprechen der klassischen Histogramm-Darstellung. Die um die Modellfunktion eingezeichneten\n", "Fehlerbalken zeigen den 68.3%-Konfidenzbereich für die unter Zugrundelegung der Poission-Verteilung\n", "erwartete Zahl an Einträgen in jedem Bin. Die Größe \"g.o.f\" (= \"goodness of fit\") ist die \n", "Differenz des beobachteten Werts von $n2l{\\cal L}$ und dem bestmöglichen Wert dieser Größe, \n", "den man erhält, wenn alle Daten auf der Modellvorhersage liegen, dem sog. \"Saturated Model\"; \n", "g.o.f. konvergiert für große Datensätze gegen die aus der Anpassung mit der Methode der kleinsten\n", "Fehlerquadrate bekannte Größe $\\chi^2$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ein **Vergleich mit der Gauß-Näherung** ist durch Angabe einer weiteren Option\n", "`use_GaussApprox=True` beim Aufruf von *hFit* leicht möglich. Hier das gleiche\n", "Beispiel, bei dem im Vergleich zu vorhin nur diese Option geändert wurde. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "#\n", "# --- perform 2nd Fit with Gaussian cost function\n", "#\n", "results_ga = hFit(\n", " SplusB_model,\n", " bcontents,\n", " bedges, # bin entries and bin edges\n", " limits=[(\"s\", 0.0, None), (\"sigma\", 0.0, None)], # parameter limits\n", " use_GaussApprox=True, # use Gaussian approxmiation of Poisson distr.\n", " plot=True, # plot data and model\n", " plot_band=True, # plot model confidence-band\n", " plot_cor=False, # plot profile likelihoods and contours\n", " axis_labels=[\"x\", \"y \\ f(x, *par)\"],\n", " data_legend=\"peudo-data\",\n", " model_legend=\"s+b, Gauss Approx.\",\n", ")\n", "\n", "# Print results\n", "# pvals, perrs, cor, chi2, pnams = results_ga\n", "print(\"\\n*==* histogram fit with Gaussian cost:\")\n", "print(\" parameter names: \", results_ga[\"parameter names\"])\n", "print(\" goodness-of-fit: {:.3g}\".format(results_ga[\"goodness-of-fit\"]))\n", "print(\" parameter values: \", results_ga[\"parameter values\"])\n", "print(\" neg. parameter errors: \", results_ga[\"confidence intervals\"][:, 0])\n", "print(\" pos. parameter errors: \", results_ga[\"confidence intervals\"][:, 1])\n", "print(\" correlations : \\n\", results_ga[\"correlation matrix\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "Die Ergebnisse sind recht ähnlich, allerdings zeigt sich ein signifikanter Unterschied \n", "beim angepassten Signalanteil, der in der Gauß-Näherung um fast die Hälfte seiner Unsicherheit größer als in der korrekten, auf der Poissonverteilung beruhenden \n", "Anpassung geschätzt wird. \n", "Der Grund dafür sind die asymmetrischen Unsicherheitsintervalle der Poissonverteilung, die \n", "im Bereich des Signals kleinere Werte der angepassten Verteilungsdichte bevorzugen. In \n", "der Praxis würde man in diesem Fall auf die Anwendung der Gauß-Näherung verzichten. \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.2 Signifikanz eines Signals und Ausschlussgrenze \n", "\n", "Mit Hilfe der Profil-Likelihood lassen sich noch viel weitergehende Aussagen über \n", "einen Parameter treffen.\n", "Ein oft anzutreffendes Beispiel ist die Suche nach einem Signal-Peak über einem\n", "Untergrund, wie wir sie oben schon behandelt haben. Die Größe des Signals, also der\n", "Wert des Parameters $s$, war dabei sehr deutlich von Null verschieden. Bei \n", "kleineren Signalen stellt sich die Frage, ob sie überhaupt \"signifikant\" sind \n", "oder nicht auch das Ergebnis einer Fluktuation des Untergrunds sein könnten.\n", "Ein häufig verwendetes Mass für die Signifikanz $S$ einer Signalstärke \"s\" ist \n", "der auf seine Unsicherheit $\\sigma_s$ normierte Abstand des Signalparameters von \n", "Null - das entpricht der Aussage:\n", "\"Wir haben ein Signal mit einer Signifikanz von z Sigma beobachtet\". Je größer\n", "dieser Abstand, desto unwahrscheinlicher ist eine statistische Fluktuation des\n", "Untergrunds als Ursache für den lokalen, als Signal fehlinterpretierten \n", "Überschusses an Histogrammeinträgen. \n", "\n", "Wenn der Verlauf der Likelihood am Minimum von der Parabelform abweicht, \n", "ist das kein gutes Maß mehr. \n", "Man spezifiziert dann den Abstand von Null als den Logarithmus des Verhältnisses der \n", "Profil-Likelihood für Signalgröße $s=0$ und am Optimum, $s=\\hat{s}$. \n", "Die Wurzel aus diesem Verhältnis stellt den sogenannten $z$-Wert dar, den wir schon\n", "ober eingeführt hatten und der die äquivalente Anzahl an Standardabweichungen einer\n", "Gaußverteilung ausdrückt, um die der Erwartungswert von Null abweicht. \n", "Es gilt\n", "\n", "\\begin{equation}\n", " z^2 = -2 \\ln \\left( \\frac{{\\cal L}_{prof}(0)} {{\\cal L}_{prof}(\\hat{s})} \\right)\n", " = \\left( -2\\ln{\\cal L}_p(s=0) \\right) - \\left( -2\\ln{\\cal L}_p(s=\\hat{s}) \\right)\n", " = -2\\Delta\\ln{\\cal L}_p(s) \\,.\n", "\\end{equation}\n", "\n", "$S=z$ ist die Signifikanz in \"Anzahl $\\sigma$s\", wie wir es oben etwas\n", "unpräzise ausgedrückt hatten. \n", "Der Term $-2\\Delta\\ln{\\cal L}$ entpricht der auf den Wert am Minimum bezogenen \n", "Kostenfunktion, die zur Besimmung der Signifikanz an den Stellen $s=0$ und \n", "$s=\\hat{s}$ bestimmt werden muss.\n", "\n", "Anders als häufig in der Physik gibt man die Signifikanz in der professionellen \n", "Statistik als $p$-Wert an, den man durch Berechnung des entsprechenden Wertes der \n", "kumulativen Verteilungsdichte der Standardnormalverteilung erhält:\n", "$p_s$ = 1-*scipy.stats.norm.cdf*($z$). \n", "Üblicherweise wählt man den $p$-Wert bzgl. der Nullhypothese, $p_0 = 1 - p_s$, \n", "der die Übereinstimmung der Beobachtung mit der Null-Hypothese quantifiziert. \n", "Je kleiner der Wert von $p_0$, desto schlechter passt die Null-Hypothese und \n", "desto gesicherter ist damit die Existenz der beobachteten, als Signal \n", "bezeichneten Abweichung.\n", "\n", "Wenn das Signal nicht hinreichende signifikant ist, sollte man eine Obergrenze\n", "für die Signalgröße angeben; üblich ist hier eine Ausschlussgrenze mit einem \n", "Konfidenznievau von cl=95%. Übersetzt in eine Änderung der Profil-Likelihood\n", "relativ zu ihrem Wert am Optimum wird diese Grenze wieder durch die Verwendung\n", "des entsprechenden Quantils der Standardnormalverteilung bestimmt, das über die \n", "Funktion $z_{cl}$ = scipy.statsnorm.ppf(cl)* (ppf = percent point function, \n", "die Inverse der kumulativen Verteilungsfunkton cdf) berechnet werden kann.\n", "Die zugehörige Änderung des Logarithmus der Profil-Likelihood ist $z_{cl}^2$\n", "und entspricht für cl=95% dem Wert 2.71. \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Um solche Untersuchungen durchführen zu können, braucht man Zugriff auf die \n", "Profil-Likelihood. Mit Hilfe der Funktion *mFit()* gelingt das, wenn man nicht \n", "die Ergebnisse der Anpassung, sondern die entsprechende Instanz der Klasse *mnFit* \n", "als Ausgabe anfordert, indem man den Parameter `return_fitObject = True` angibt. \n", "Damit hat man Zugriff auf alle Daten und Methoden der Klasse *mnFit*.\n", "\n", "Um die Signalsuche etwas realistischer zu gestalten, wurde in den Daten im Beispiel \n", "von vorhin das Signal um einen Faktor Zwei verkleinert. \n", "Damit die Anpassung auch für kleine angenommene Signalstärken noch konvergiert, \n", "sollten noch Position und Breite des gesuchten Signals fixiert werden und \n", "vorher geeignete Werte gesetzt werden. \n", "\n", "Der leicht modifizierte Code für unser bereits bekanntes Beispiel sieht nun so aus: " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# the modified histogram data\n", "nbins = 40\n", "xmn = 1\n", "xmx = 10\n", "bedges = np.linspace(xmn, xmx, nbins + 1)\n", "bcontents2 = np.array(\n", " [\n", " 1,\n", " 1,\n", " 1,\n", " 2,\n", " 2,\n", " 2,\n", " 4,\n", " 1,\n", " 0,\n", " 3,\n", " 1,\n", " 1,\n", " 0,\n", " 2,\n", " 3,\n", " 3,\n", " 1,\n", " 1,\n", " 0,\n", " 2,\n", " 3,\n", " 2,\n", " 3,\n", " 1,\n", " 1,\n", " 4,\n", " 3,\n", " 4,\n", " 4,\n", " 1,\n", " 0,\n", " 1,\n", " 2,\n", " 1,\n", " 3,\n", " 2,\n", " 1,\n", " 3,\n", " 2,\n", " 4,\n", " ]\n", ")\n", "# --- perform fit\n", "hFitObj = hFit(\n", " SplusB_model,\n", " bcontents2,\n", " bedges, # bin entries and bin edges\n", " p0=[7.0, 0.25, 0.2], # set (initial) parameter values\n", " fixPars=[\"mu\", \"sigma\"], # fix signal position and shape\n", " limits=[(\"s\", 0.0, None)], # parameter limits\n", " return_fitObject=True,\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Das *Python-Dictionary* mit den Ergebnissen erhält man mit der Methode *getResult()*\n", "der Klasse *mnFit*: \n", "```\n", "resultDict = hFit.getResult()\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# get and print results\n", "pvals, perrs, cor, gof, pnams = hFitObj.getResult().values()\n", "print(\"\\n*==* histogram fit with small signal:\")\n", "print(\" parameter names: \", pnams)\n", "print(\" parameter values: \", pvals)\n", "print(\" neg. parameter errors: \", perrs[:, 0])\n", "print(\" pos. parameter errors: \", perrs[:, 1])\n", "print(\" correlations : \\n\", cor)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Die Profil-Likelihood für einen Parameter mit dem Namen *name* wird mit dem Aufruf \n", "der Methode *getProfile()* der Klasse *mnFit* berechnet und zurückgegeben:\n", "```\n", "pvals, nlLp, status = hFitObj.getProfile(name, range, npoints)\n", "```\n", "\n", "Ähnlich wie im grafische Beispiel in Kap. 2 verwedenden wir ein Hilfsfunktion\n", "*analyse_signal()*, um die Profil-Likelihood auzuwerten. Dieses Mal interessiert\n", "uns aber nicht der Zentralbereich der Kurve in der Nähe des Minimums, sondern\n", "die Ränder. Hier der Code zur Erzeugung und Auswertung der Profil-Likelihood,\n", "der auch eine anschauliche grafische Darstellung der in der Beschreibung oben\n", "skizzierten Zusammenhänge erzeugt: " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# get profile likelihood (in range [0., mx]\n", "mx = 0.225\n", "vals, n2lLp, status = hFitObj.getProfile(\"s\", range=(0.0, mx), npoints=50)\n", "\n", "\n", "# analyze profile likelihood\n", "def analyze_signal(vals, n2lLp):\n", " from scipy.stats import norm\n", " from scipy.interpolate import CubicSpline\n", " from scipy.optimize import newton\n", "\n", " n2lLp_spline = CubicSpline(vals, n2lLp)\n", "\n", " # Determine significance of parameter\n", " z = np.sqrt(n2lLp_spline(0.0))\n", " p_0 = 1.0 - norm.cdf(z)\n", "\n", " # Determin upper limit on parameter\n", " z_95 = norm.ppf(0.95) # z-value for 95%\n", "\n", " def f(x):\n", " return n2lLp_spline(x) - z_95**2\n", "\n", " s95 = newton(f, x0=pvals[0], x1=mx)\n", "\n", " # - plot profile likelihood\n", " plt.plot(vals, n2lLp)\n", " plt.xlabel(\"Signal fraction s\", size=\"x-large\")\n", " plt.ylabel(\"$-2 \\Delta \\ln L_p$\", size=\"x-large\")\n", "\n", " # - plot significance\n", " z2 = z * z\n", " plt.arrow(0, 0, 0, z2, head_width=0.01, color=\"orange\")\n", " plt.text(-mx / 15, z2, \" $z^2$\", size=\"large\", color=\"sienna\")\n", "\n", " # - plot limit\n", " plt.hlines(z_95**2, 0.66 * mx, mx, color=\"maroon\", linestyle=\"dotted\")\n", " plt.vlines(s95, 0, n2lLp_spline(mx), color=\"maroon\")\n", " plt.text(\n", " mx,\n", " z_95**2,\n", " \" ${{z_{{95}}}}^2=$\" + str(\"{:.3g}\".format(z_95**2)),\n", " size=\"large\",\n", " color=\"darkred\",\n", " )\n", " plt.text(s95, 0.0, \" 95% limit\", size=\"large\", color=\"darkblue\")\n", "\n", " return z, p_0, s95\n", "\n", "\n", "# perform analysis of signal properties using function\n", "z, p_0, s95 = analyze_signal(vals, n2lLp)\n", "print(\"\\n\")\n", "print(\" == Analysis of profile likelihood\")\n", "print(\" Significance, z-value: {:.3g}\".format(z))\n", "print(\" p-value of Null-hypothesis: {:.3g}\".format(p_0))\n", "print(\" 95% limit: \" + pnams[0] + \" < {:.3g}\".format(s95))\n", "print(\"\\n\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Das Signal hat also eine relativ schwache Signifikanz von nur\n", "$z=2.14$ (oder \"$2.14 \\sigma$\"); dies entspricht einem $p$-Wert \n", "für die Null-Hypothese von $1.6\\%$. Üblicherweise verlangt man\n", "$z$-Werte von mindestens $3$, häufig sogar $z \\ge 5$, um sich gegen\n", "falsche Entdeckungen abzusichern, die auftreten, wenn man an sehr\n", "viele Stellen sucht, also mit verschiedenen Methoden in vielen\n", "Verteilungen an verschiedenen, vorher nicht bekannten Stellen\n", "im Parameterraum (der sog. \"look-elsewhere effect\"). \n", "\n", "Im vorliegenden Fall würde man also nicht von der Entdeckung eines\n", "Signals sprechen, sondern eine Obergrenze auf einen möglichen Signalanteil\n", "von $s < 0.188$ bei einem Konfidenzniveau von 95% setzen." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "## 4. Maximum-Likelihood und Methode der kleinsten Fehlerquadrate\n", "***\n", "\n", "Die Methode der kleinsten Fehlerquadrate zur Anpassung von Modellen ist\n", "ein Spezialfall des Maximum-Likelihood Verfahrens.\n", "Um das zu sehen, leiten wir allgemein ein auf der log-Likelihood Methode \n", "basierendes Anpassungsverfahren für Modellfunktionen an Messdaten her.\n", "\n", "Wir bezeichnen die zufällige Abweichung eines Messwertes vom wahren\n", "Wert mit $z$, die durch eine Verteilungsdichte $f_z$\n", "beschrieben wird. Im Falle von Messunsicherheiten ist $f_z$ häufig\n", "die Normalverteilung mit Erwartungswert Null und Standardabweichung $\\sigma$,\n", "${\\cal N}(z; \\sigma)$. Im mehrdimensionalen Fall für $n_d$ nicht\n", "notwendigerweise unabhängige Datenpunkte ist die multivariate\n", "Gaußverteilung ${\\cal N}(\\vec z; V)$ mit der Kovarianzmatrix $V$ relevant,\n", "\n", "\\begin{equation}\\label{equ-mStandardGauss}\n", " {\\cal{N}}\\left(\\vec{z}, V) \\right) = \n", "\\frac{1} {\\sqrt{(2\\pi)^{n_d} \\det(V)} }\n", "\\cdot\n", "\\exp \\left( -\\frac{1}{2} \\vec{z}^T V^{-1} \\vec{z}\n", "\\right) \\,.\n", "\\end{equation}\n", "\n", "Die Zufallsgröße $z$ mit Erwartungswert Null entspricht den Fluktuationen \n", "um den wahren Wert. Wenn es für den wahren Wert eine exakte theoretische \n", "Erwartung in Form eines parameterbehafteten Modells $f_i(\\vec p)$ mit einem \n", "Satz an Parametern $\\vec p$ gibt, lässt sich ein Messwert $y_i$ schreiben als\n", "\\begin{equation}\\label{equ-zplusw}\n", "y_i = f_i(\\vec p) + z_i \\,.\n", "\\end{equation}\n", "\n", "Üblicherweise führt man zum Test von parameterbehafteten, durch\n", "Funktionen beschriebenen Modellen mehrere Messungen an verschiedenen\n", "Stützstellen $x_i=(\\vec x)_i$ durch, man betrachtet also eine Modellfunktion\n", "$f(\\vec x; \\vec p)$.\n", "Die Verteilungsdichte, die alle Messungen beschreibt, sieht dann so aus:\n", "\n", "\\begin{equation}\\label{equ-mGauss}\n", " {\\cal{N}}\\left(\\vec{y}, V, f(\\vec{x}, \\vec{p}) \\right) = \n", " \\frac{1} {\\sqrt{(2\\pi)^{n_d} \\det(V)} } \\cdot\n", " \\exp\\left( -\\frac{1}{2} (\\vec{y}-\\vec{f})^T V^{-1} (\\vec{y}-\\vec{f}) \\right) \\,.\n", "\\end{equation}\n", "\n", "Nach dem Maximum-Likelihood-Prinzip ist der beste Parametersatz durch\n", "den Punkt ${\\hat{\\vec{ p}}}$ im Parameterraum gegeben, für den die Likelihood\n", "maximal wird. Wie schon oben verwendet man das Zweifache des negativen natürlichen\n", "Logarithmus der Likelihood, $n2l{\\cal L}_{Gauß}$, der dann minimiert wird.\n", "\n", "\\begin{equation}\\label{equ-nlLGauss}\n", "-2\\, \\ln{\\cal{L_{\\rm Gauß}}}\\left( \\vec y, V, \\vec{f}(\\vec x, \\vec {p}) \\right)\n", "\\,=\\, \\left(\\vec y - \\vec f(\\vec x; \\vec p ) \\right)^T V^{-1}\n", " \\left(\\vec y - \\vec f(\\vec x; \\vec p ) \\right)\n", " + \\ln(\\det(V)) + n_d\\,\\ln(2\\pi) \\,.\n", "\\end{equation}\n", "\n", "Für die Bestimmung des Minimums von $n2l{\\cal L}_{Gauß}$ sind bzgl. der\n", "Parameter konstante Terme nicht relevant, man kann sie daher weglassen.\n", "Wenn die Normierung der Gaußverteilung, also der Term $\\det(V)$, nicht\n", "von den Parametern abhängt, vereinfacht sich der obige Ausdruck\n", "zum altbekannten Ausdruck für die quadratische Residuensumme mit \n", "Kovarianzmatrix, $S$:\n", "\n", "\\begin{equation}\\label{equ-chi2}\n", "S\\left(\\vec y, V, \\vec{f}(\\vec x, \\vec {p})\\right) \\,=\\,\n", " \\left(\\vec y - \\vec f(\\vec x; \\vec p ) \\right)^T V^{-1}\n", " \\left(\\vec y - \\vec f(\\vec x; \\vec p ) \\right) \\,.\n", "\\end{equation}\n", "\n", "Unter speziellen Bedingungen sind also die Minimierung von \n", "$n2l{\\cal{L}}$ und der quadratischen Summe der Residuen $S$ \n", "äquivalent. \n", "Dieser einfache Fall ist aber nicht mehr gegeben, wenn relative, \n", "auf den Modellwert bezogene Unsicherheiten auftreten, oder wenn \n", "Unsicherheiten in Abszissen-Richtung behandelt werden sollen, die \n", "mit Hilfe einer Taylor-Entwicklung erster Ordnung von der Abszisse \n", "auf die Ordinate übertragen werden und damit von der Ableitung des \n", "Modells nach $x$ und so auch von den Parameterwerten abhängen.\n", "Die praktische Regel lautet daher, möglichst immer die Likelihood\n", "zu verwenden, und nur in gut begründeten, berechtigten Fällen\n", "auf die Methode der kleinsten Fehlerquadrate zurück zu greifen.\n", "Allerdings ist man bei der Wahl der verfügbaren numerischen\n", "Werkzeuge stark eingeschränkt, wenn man diese Empfehlung\n", "umsetzen möchte.\n", "\n", "Sollen Problemstellungen mit nicht Gauß-förmigen Verteilungen\n", "behandelt werden, ist die Verwendung von log-Likelihoodverfahren\n", "unumgänglich. Das Aufstellen entsprechender Likelihood-Funktionen\n", "zur Behandlung spezieller Problemstellungen gehört in der\n", "wissenschaftlichen Praxis heute zum Standard. Dank sehr leistungsfähiger\n", "Algorithmen zur numerischen Minimierung in hoch-dimensionalen\n", "Parameterräumen und auch Dank moderner Computertechnik stellt\n", "die Verwendung korrekter Likelihood-Verfahren kein\n", "unüberwindliches Problem mehr dar." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Praktische Hinweise\n", "\n", "In der Praxis sind die oben beschriebenen Verfahren nur in Kombination mit \n", "numerischen Verfahren zur Minimierung von skalaren Funktionen in mehr- oder\n", "sogar hoch-dimensionalen Räumen und Programmcode zur Verwaltung der Daten\n", "und ihrer Unsicherheiten durchführbar. Obwohl in seltenen Fällen auch analytische\n", "Lösungen existieren, nutzt man in der Praxis fast ausschließlich Programmpakete\n", "zur Durchführung von Anpassungen; analytische (Teil-)Lösungen müssen nur\n", "eingesetzt werden, wenn es um zeitkritische Problemstellungen geht. \n", "\n", "#### Konstruktion der Kovarianzmatrix\n", "\n", "Die Kovarianzmatrix der Daten bildet die Unsicherheiten Gauß-verteilter \n", "Eingabedaten vollständig ab. Wird sie bei der Anpassung berücksichtigt, \n", "werden alle zwischen den Daten unabhängigen und korrelierten Unsicherheiten\n", "beim Anpassungsprozess in das Endergebnis propagiert - eine klassische \n", "Fehlerrechnung von Hand ist dann nicht mehr notwendig. \n", "\n", "Fassen wir zunächst kurz die wesentlichen Eigenschaften der\n", "Kovarianzmatrix ${\\bf V}$ zusammen:\n", "\n", "- ${\\bf V} = \\left(\\mathrm{V}_{ij}\\right)$ ist eine symmetrische Matrix;\n", "- sie hat die Dimension $n_d$, 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", "- die Kovarianzmatrix-Elemente für voneinander unabhängige Unsicherheiten\n", " werden addiert.\n", "\n", "Gerade der letzte Punkt ist von entscheidender Bedeutung, denn er erlaubt es,\n", "die vollständige Kovarianzmatrix sukzessive aus einzelnen Beiträgen zur\n", "Unsicherheit aufzubauen, mit der Konstruktionsvorschrift:\n", "\n", "- Unsicherheiten der Messwerte werden nach verschiedenen, unabhängigen\n", " Quellen aufgeschlüsselt.\n", "- Unabhängige Unsicherheiten jeder einzelnen Messung werden quadratisch in\n", " den Diagonalelementen aufaddiert.\n", "- Allen Messwerten oder Gruppen von Messungen gemeinsame absolute oder\n", " relative Unsicherheiten werden quadratisch in den betreffenden\n", " Diagonal- und Nebendiagonalelementen $V_{ii}$, $V_{jj}$ und $V_{ij}$ aufaddiert.\n", "\n", "Unsicherheiten in Abszissenrichtung, gegeben durch eine Kovarianzmatrix\n", "${\\bf V}^x$, können berücksichtigt werden, in dem man sie mit Hilfe der\n", "Modellvorhersage $f(\\vec x, \\vec p)$ per Taylor-Entwicklung in erster Ordnung\n", "in y-Richtung transformiert und dann zur Kovarianz-Matrix der Datenpunkte in\n", "y-Richtung, ${\\bf V}^y$, addiert:\n", "\n", "\\begin{equation}\\label{equ-xyCovariance}\n", " V_{ij} = (V^y)_{ij} \n", " + \\frac{\\partial f}{\\partial x_i} \\frac{\\partial f}{\\partial x_j} (V^x)_{ij} \n", "\\end{equation}\n", "\n", "Insgesamt ergeben sich so acht Arten von Unsicherheiten, nämlich \n", "unabhängige und / oder korrelierte\n", "absolute und / oder relative Unsicherheiten in\n", "$x$- und / oder $y$-Richtung.\n", "\n", "Zum Bau der Kovarianzmatrix setzt man idealerweise Programmcode ein. Die\n", "wenigsten gängigen Anpassungsprogramme bringen dafür direkt Optionen mit,\n", "sondern es muss die vollständige Kovarianzmatrix als Parameter übergeben\n", "werden.\n", "Das Paket *kafe2* enthält die Methode \n", "`add_error(err_val, axis=?, correlation=?, relative=?, reference=?)` \n", "mit deren Hilfe einzelne Komponenten der Unsicherheit hinzugefügt werden\n", "können. Die Interfaces zu verschiedenen Anpassungsprogrammen\n", "im Paket *PhyPraKit* sehen ebenfalls die Angabe einzelner Komponenten\n", "der Unsicherheit als Parameter vor, wenn die entsprechenden Anpassungspakete\n", "dies unterstützen.\n", "Alle acht Arten von Unsicherheiten können direkt nur mit *kafe2* oder \n", "mit *PhyPraKit.phyFit* behandelt werden. \n", "\n", "#### Berücksichtigung externer und eingeschränkter Parameter\n", "\n", "Häufig hängen Modellfunktionen von externen, mit Unsicherheiten behafteten \n", "Parametern ab, die z.$\\,$B. in einer Hilfsmessung bestimmt wurden oder aus der\n", "Literatur stammen. Dies können die Ergebnisse von Kalibrationsmessungen sein,\n", "oder auch Natur- oder Apparatekonstanten. Statt die Effekte der Parameterunsicherheiten \n", "mit Hilfe einer händischen Fehlerrechnung auf das Endergebnis zu propagieren, \n", "können sie auch als eingeschränkte Parameter (\"constrained parameter\") direkt \n", "in der Anpassung berücksichtigt werden. \n", "Dazu werden solche Parameter gleichzeitig als freie Parameter in der \n", "Anpassung und als Messgrößen eingeführt - d.$\\,$h. ein entsprechender Term \n", "zur log-Likelihood hinzugefügt. \n", "\n", "Notwendig ist bisweilen auch eine Methode, mit deren Hilfe Parameter auf \n", "feste Werte fixiert werden können.\n", "Der Einfluss externer Parameter kann damit auch untersucht werden, indem\n", "man sie nacheinander auf ihren Erwartungswert und die jeweiligen oberen bzw. \n", "untere Grenzen ihres Konfidenzbereichs fixiert und die Veränderungen des \n", "Anpassungsergebnisses beobachtet.\n", "\n", "Das temporäre Fixieren und wieder frei geben von Parametern kann auch hilfreich \n", "sein, wenn eine komplexe Anpassung nicht zum globalen Minimum konvergiert. Man \n", "kann dann einen oder einige Parameter in der Nähe des erwarteten Wertes fixieren \n", "und eine Minimierung bezüglich der übrigen Parameter vornehmen. Wenn man dann an \n", "diesem temporären Minimum die fixierten Parameter wieder frei gibt, sollte \n", "die Anpassung zum globalen Minimum konvergieren. \n", "\n", "Ein Problem während des Anpassungsprozesses stellen oft auch Parameter dar,\n", "die temporär Werte in mathematisch nicht definierten Wertebereichen oder in \n", "unphysikalischen Bereichen annehmen (negative Werte unter Wurzeln, negative\n", "Massen etc). Das Setzen von Limits zum Vermeidung solcher Parameterbereiche\n", "ist eine wichtige Option für jedes Anpassungspaket. \n", "\n", "Bei nichtlinearen Problemstellungen gibt es häufig neben dem globalen Minimum\n", "weitere Nebenminima - oder sogar mehrere oder viele gleichwertige Lösungen. \n", "In solchen Fällen werden Startwerte für die Parameter in der Nähe\n", "einer \"vernünftigen\" Lösung benötigt. Auch diese Möglichkeit muss ein \n", "Anpassungspaket bieten. Da im Laufe des besseren Verständnisse aus manchem \n", "einfachen Problem später ein nichtlineares Problem werden kann, ist es eine \n", "gute Angewohnheit, grundsätzlich immer Startwerte zu setzen.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ein Beispiel mit allen acht oben genannten Arten von Unsicherheiten ist in der Code-Zelle\n", "unten gezeigt. Alle in der Funkion *xyFit()* verfügbaren Optionen sind ebenfalls angegeben,\n", "auch wenn für die meisten \"vernünftige\" Voreinstellungen gelten und man sich die explizite\n", "Angabe häufig ersparen kann. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Beispiel zur Anpassung an x/y-Daten\n", "\n", "Ein Beispiel, bei dem alle acht Fehlertypen vorkommen, ist im folgenden Code-Beispiel zur \n", "Anwendung der Funktion *xyFit*" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ " from PhyPraKit.phyFit import xyFit\n", "\n", "#\n", "# *** Example of an application of phyFit.mFit()\n", "#\n", "# define the model function to fit\n", "def exp_model(x, A=1., x0=1.):\n", " return A*np.exp(-x/x0)\n", "\n", "# another model function\n", "def poly2_model(x, a=0.1, b=1., c=1.):\n", " return a*x**2 + b*x + c\n", "\n", "# set model to use in fit\n", "#fitmodel=exp_model # also try poly2_model !\n", "fitmodel=poly2_model # also try exp_model!\n", " \n", "# the data ...\n", "data_x = [0.0, 0.2, 0.4, 0.6, 0.8, 1., 1.2,\n", " 1.4, 1.6, 1.8, 2.0, 2.2, 2.4, 2.6]\n", "data_y = [1.149, 0.712, 0.803, 0.464, 0.398, 0.354, 0.148,\n", " 0.328, 0.181, 0.140, 0.065, 0.005,-0.005, 0.116]\n", "# ... and components of the uncertaity \n", "sabsy = 0.07 # independent y\n", "srely = 0.05 # 5% of model value\n", "cabsy = 0.04 # correlated\n", "crely = 0.03 # 3% of model value correlated\n", "sabsx = 0.05 # independent x\n", "srelx = 0.04 # 4% of x\n", "cabsx = 0.03 # correlated x\n", "crelx = 0.02 # 2% of x correlated\n", "\n", "# perform fit to data with function xyFit using class mnFit\n", "results = xyFit(fitmodel, # the model function\n", " data_x, data_y, # x and y data\n", " sx=sabsx, # the error components defined above\n", " sy=sabsy,\n", " srelx=srelx,\n", " srely=srely,\n", " xabscor=cabsx,\n", " xrelcor=crelx,\n", " yabscor=cabsy,\n", " yrelcor=crely,\n", "# p0=(1., 0.5), # start values of parameters (taken from model if not given)\n", "# constraints=['A', 1., 0.03], # constraint(s)\n", "# constraints=[0, 1., 0.03] # alternative specificaiton of constraints\n", "# limits=('A', 0., None), # parameter limits\n", "# fixPars = ['A'], # parameter(s) to be fixed to start value \n", " use_negLogL=True, # use full n2lL (default)\n", " plot=True, # plot data and model \n", " plot_band=True, # plot model uncerctainty\n", " plot_cor=True, # plot profiels and contours\n", " quiet=True, # silent output it True\n", " axis_labels=['x', 'y \\ f(x, *par)'], \n", " data_legend = 'pseudo data', \n", " model_legend = 'model')\n", "\n", "# Print results \n", "print('\\n*==* xyFit Result:')\n", "for key in results:\n", " print(\"{}\\n\".format(key), results[key])\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "## 5. Beurteilung der Qualität von Anpassungsverfahren\n", "***\n", "\n", "Die oben schon genannten Qualitätsanforderungen an ein Schätzverfahren können \n", "je nach gewähltem Verfahren, aber auch abhängig von der Größe einer Stichprobe \n", "bzw. der gewählten Parametrisieung der Modellfunktion unterschiedlich gut erfüllt\n", "sein. \n", "\n", "Hier nich einmal die Zusammenstellung der relevanten Eigenschaften,die eine\n", "Parameterschätzung erfüllen sollte:\n", "\n", " - konsistent: $\\displaystyle \\lim_{n \\to \\infty} \\hat{p}_i = p_i$\n", " \n", " - möglichst erwartungstreu schon bei kleinen Stichproben: $E[\\hat{p}_i] = p_i$ \n", " \n", " - möglichst unverzerrt: Verzerrung $b$ (=bias) $b_i=E[\\hat{p}_i] - p_i$ = 0, \n", " für große Stichproben muss wegen der geforderten Konsistenz gelten\n", " $\\displaystyle \\lim_{n \\to \\infty} b = 0$\n", " \n", " - optimal: $V[\\hat{p}_i]$ minimal \n", " \n", " - robust: unbeeinflusst von \"Ausreißern\" oder leichten Fehlmodellierungen der Daten\n", " \n", " - korrekte Überdeckungswahrscheinlichkeit (engl. \"coverage\"): \n", " $W(p_i \\in [\\hat{p}_i - \\Delta p_i, \\hat{p}_i +\\Delta p_i])= 68.3\\%$ \n", "\n", "Die hier diskutierte Maximum-Likelihood Methode (MLE) ist optimal, allerdings eher nicht\n", "robust, weil falsche Annahmen über das Modell einen empfindlichen Einfluss auf die\n", "Parameterwerte haben können. MLE-Schätzungen sind in der Regel auch nicht unverzerrt; \n", "wenn nicht unbedingt der Punkt der besten Anpassung (\"Punktschätzung\"), sondern, wie \n", "häufig in der Physik ein möglichst korrektes Konfidenzintervall (\"Intervallschätzung\") \n", "für den Schätzwert von Belang ist, spielt das allerdings eine untergeordnete Rolle. \n", "Wie wir schon oben gesehen hatten, existiert für die Bestimmung der Unsicherheiten nur \n", "eine Ungleichung (Cramér-Rao-Fréchet Grenze). Wie gut bestimme Verfahren typischerweise\n", "nach den oben genannten Kriterien funktionieren, ist oft empirisch bekannt. Mann kann\n", "zeigen, dass z.$\\,$B. die Methode der kleinsten Fehlerquadrate für lineare Problemstellungen\n", "und Gauß-förmige Messunsicherheiten optimal und unverzerrt ist. Die hier skizzierten\n", "Likelihood-Verfahren erfüllen die oben genannten Bedingungen schon für recht kleine \n", "Stichproben ebenfalls gut. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Überprüfung mittels Monte-Carlo - Methode (\"toy-MC\")\n", "\n", "Zur Überprüfung der Eigenschaften eines Schätzverfahrens und des Einflusses der Stichprobengröße\n", "werden in der Praxis häufig Ensemble-Tests durchgeführt - d.h. man erzeugt zufällig eine große \n", "Anzahl gleichwertiger Stichproben und führt jeweils das Anpassungsverfahren durch.\n", "Aus einer statistischen Analyse der Ergebnisse lassen sich dann die Eigenschaften des\n", "Schätzverfahrens bestimmen.\n", "\n", "Zur Simulation der Stichproben verwendet man die Monte-Carlo Methode; weil die Simulation\n", "für diesen Zweck sehr einfach sein kann, nennt man solche Studien auch \"toy-Monte Carlo\".\n", "\n", "Beispielcode zur Erzeugung künstlicher \"Pseudo-Daten\" ist in den folgenden Codezellen gezeigt.\n", "Die durchgeführte Anpassung entspricht dem schon oben diskutierten Beispiel zur Anpassung\n", "einer Modellfunktion an x/y-Daten mit verschiedenen Typen von Unsicherheiten. \n", "\n", "\n", "Hier der Code-Teil zur Definition aller Parameter zur Datenerzeugung und zur Anpassung:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from PhyPraKit import generateXYdata\n", "from PhyPraKit.phyFit import xyFit, get_functionSignature\n", "import matplotlib.pyplot as plt\n", "from scipy import stats\n", "\n", "\n", "# define some models to test\n", "def exp_model(x, A=1.0, x0=1.0):\n", " return A * np.exp(-x / x0)\n", "\n", "\n", "# another model function\n", "def poly2_model(x, a=0.75, b=1.0, c=1.0):\n", " return a * x**2 + b * x + c\n", "\n", "\n", "# set model to use ...\n", "model = poly2_model\n", "# model = exp_model\n", "\n", "# ... and number of pseudo-experiments to perform\n", "Nexp = 200\n", "\n", "# parameters of pseudo-data\n", "nd = 120 # number of data points\n", "xmin = 0.0 # x-range\n", "xmax = 2.5\n", "data_x = np.linspace(xmin, xmax, nd) # x of data points\n", "\n", "# extract model paramerters from model function\n", "mpardict = get_functionSignature(model)[1] # keyword arguments of model\n", "true_vals = np.asarray(list(mpardict.values()))\n", "npar = len(mpardict)\n", "parnams = list(mpardict.keys())\n", "\n", "# set uncertainties for x/y pseudo-data\n", "# - independent absolute uncertainties\n", "sabsx = 0.05 # 0.05\n", "sabsy = 0.07 # 0.07\n", "# - indepentent relative uncertainties\n", "srelx = 0.04 # 4%\n", "srely = 0.03 # 5% of model value\n", "# correlated uncertainties in x and y direction\n", "cabsx = 0.03 # 0.03\n", "crelx = 0.02 # 2%\n", "cabsy = 0.04 # 0.04\n", "crely = 0.03 # 3% of model value\n", "\n", "# calulate total uncertainties\n", "sigx = np.sqrt(sabsx * sabsx + (srelx * data_x) ** 2)\n", "sigy = np.sqrt(sabsy * sabsy + (srely * model(data_x, **mpardict)) ** 2)\n", "\n", "# set keyword arguments for fit function xyFit\n", "kw_uncertainties = { # - components of uncertainty\n", " \"sx\": sabsx, # indep x\n", " \"sy\": sabsy, # indel y\n", " \"srelx\": srelx, # indep. rel. x\n", " \"srely\": srely, # indep. rel. y\n", " \"xabscor\": cabsx, # correlated x\n", " \"xrelcor\": crelx, # correlated rel. x\n", " \"yabscor\": cabsy, # correlated y\n", " \"yrelcor\": crely, # correlated rel. y\n", "}\n", "\n", "kw_options = { # options for Fit\n", " \"ref_to_model\": True, # reference of rel. uncert. to model\n", " \"use_negLogL\": True, # standard chi^2 if false\n", "}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Die Schleife zur Erzeugung der Daten, Modellanpassung und die Analyse der Ergebnisse \n", "ist in der folgenden Codezelle gezeigt. Für die Analyse der individuellen Ergebnisse \n", "der Anpassung sind die folgenden Zeilen zuständig. Sie speichern die individuellen \n", "Ergebnisse, berechnen die $\\chi^2$-Wahrscheinlichkeit, die Verzerrung und eine Bool'sche \n", "Variable für die Auswertung der Überdeckung zur späteren statistischen Auswertung:\n", "\n", "```\n", " # calculate chi2 probability\n", " chi2prb = 1.- stats.chi2.cdf(chi2, nd-len(pvals))\n", " c2prb.append(chi2prb)\n", " \n", " # analyze bias and coverage\n", " for i in range(npar):\n", " biases[i].append( pvals[i] - true_vals[i] ) # bias\n", " # get parameter confidence interval (CI)\n", " pmn = pvals[i] + perrs[i,0]\n", " pmx = pvals[i] + perrs[i,1]\n", " # coverage: count how often true value is in CI\n", " if (true_vals[i] >= pmn and true_vals[i]<=pmx): N_coverage[i] +=1\n", "```\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# this is the loop part of the MC simulation\n", "\n", "\n", "def MC_loop():\n", " global nfail, biases, c2prb, N_coverage\n", " # initialize arrays for statistical analysis in loop\n", " nfail = 0 # counter for failed fits\n", " biases = [[] for i in range(npar)] # bias\n", " c2prb = [] # chi^2 probatility\n", " N_coverage = npar * [0] # counters for coverage\n", "\n", " #\n", " # start toy-MC loop --------------------------------------------\n", " #\n", " print(\"generating data - patience please!\")\n", " for iexp in range(Nexp):\n", " # generate a sample of peudo-data\n", " np.random.seed(\n", " 314159 + (iexp * 2718281) % 100000\n", " ) # initialize random generator\n", " xt, yt, data_y = generateXYdata(\n", " data_x,\n", " model,\n", " sigx,\n", " sigy,\n", " xabscor=cabsx,\n", " xrelcor=crelx,\n", " yabscor=cabsy,\n", " yrelcor=crely,\n", " mpar=true_vals,\n", " )\n", " # perform fit to pseudo data\n", " try:\n", " plot = True if iexp < 1 else 0\n", " result = xyFit(\n", " model,\n", " data_x,\n", " data_y, # model & data\n", " **kw_uncertainties, # all uncertainties\n", " **kw_options, # fit options\n", " plot=plot, # plot data and model\n", " showplots=False, # call plt.show() in user code if False\n", " )\n", "\n", " # save fit results in local variables\n", " pvals, perrs, cor, chi2, pnams = result.values()\n", "\n", " if plot:\n", " plt.suptitle(\"Result of fit {:d}\".format(iexp))\n", "\n", " if iexp < 1:\n", " # Print results to illustrate how to use output\n", " np.set_printoptions(precision=6)\n", " print(\"\\n*==* Fit {:d} Result:\".format(iexp))\n", " print(f\" chi2: {chi2:.4g}\")\n", " print(f\" parameter names: \", pnams)\n", " print(f\" parameter values: \", pvals)\n", " np.set_printoptions(precision=3)\n", " print(f\" parameter errors: \", perrs)\n", " np.set_printoptions(precision=3)\n", " print(f\" correlations : \\n\", cor)\n", " np.set_printoptions(precision=8) # default output precision\n", " print()\n", " else:\n", " # show simple progress bar\n", " if (100 * iexp) % Nexp == 0:\n", " print(\"*\", end=\"\")\n", "\n", " # calculate chi2 probability\n", " chi2prb = 1.0 - stats.chi2.cdf(chi2, nd - len(pvals))\n", " c2prb.append(chi2prb)\n", "\n", " # analyze bias and coverage\n", " for i in range(npar):\n", " biases[i].append(pvals[i] - true_vals[i]) # bias\n", " # get parameter confidence interval (CI)\n", " pmn = pvals[i] + perrs[i, 0]\n", " pmx = pvals[i] + perrs[i, 1]\n", " # coverage: count how often true value is in CI\n", " if true_vals[i] >= pmn and true_vals[i] <= pmx:\n", " N_coverage[i] += 1\n", "\n", " except Exception as e:\n", " nfail += 1\n", " print(\"!!! fit failed \", nfail)\n", " print(e)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Auswertung der Ergebnisse\n", "\n", "Die Code zur Analyse der in der Schleife aufgezeichneten Daten gibt zunächst die für\n", "jeden der Parameter beobachtete Verzerrung aus, zusammen mit deren durch die endliche\n", "Statistik bedingten Parameterunsicherheiten. Statistisch signifikant sind nur Verzerrungen, \n", "die deutlich größer als die Unsicherheit auf den Mittelwert der Abweichungen sind. Zur Bewertung\n", "der Relevanz muss man sie in Bezug auf die Unsicherheit der Parameterwerte, also der Breite\n", "der Verteilung der Abweichungen, setzen. Erst wenn sie eine bedeutsame Größe davon erreichen, \n", "typischerweise größer als 5% - 10%, sollten die Verzerrungen angegeben werden. \n", "\n", "Als zweite Größe wird die Coverage ausgegeben, die dem Anteil der Fälle entspricht, in\n", "denen der wahre Parameterwert innerhalb des Unsicherheitsintervalls liegt." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# run MC loop\n", "MC_loop()\n", "print(\"\\n\\n*==* sucessfully analyzed {:d} data sets\".format(Nexp - nfail))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def print_result():\n", " global nfail, biases, c2prb, N_coverage\n", "\n", " # - convert output to numpy arrays\n", " for i in range(npar):\n", " biases[i] = np.array(biases[i])\n", " c2prb = np.array(c2prb)\n", "\n", " ## plt.show() # blocks at this stage until figure deleted\n", " # results as printout\n", " N_succ = Nexp - nfail\n", " print(\"\\n\\n*==* \", N_succ, \" successful fits done:\")\n", " print(\" * parameter names:\")\n", " for i in range(npar):\n", " print(\" {:d}: {:s}\".format(i, pnams[i]))\n", "\n", " print(\" * biases:\")\n", " for i in range(npar):\n", " # bias = deviation of parameters from their true values\n", " b = biases[i].mean()\n", " s = biases[i].std()\n", " e = s / np.sqrt(Nexp)\n", " print(\" {:d}: {:.3g} +\\- {:.2g}, std {:.3g}\".format(i, b, e, s))\n", "\n", " print(\" * coverage:\")\n", " for i in range(npar):\n", " # coverage: fraction of true val in confidence interval\n", " p_coverage = N_coverage[i] / (N_succ) * 100 // 0.682689492\n", " print(\" {:d}: {:.3g}%\".format(i, p_coverage))\n", "\n", "\n", "print_result()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In diesem Beispiel gibt es keine signifikanten Verzerrungen, und die Überdeckung liegt\n", "in der Nähe des für ein $\\pm 1\\sigma$-Intervall erwarteten Werts, d.h. 100% von 68.3%." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Grafische Darstellung der Ergebnisse\n", "\n", "Eine einfache grafische Darstellung der Ergebnisse zeigt das Code-Fragment unten. \n", "Verwendung findet eine allgemeine Funktion zur Darstellung von Verteilungen und Korrelationen \n", "eines Datensatzes mit mehreren Variablen aus dem Paket *PhyPraKit*. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from PhyPraKit import plotCorrelations\n", "\n", "\n", "def plot_correlations():\n", " global nfail, biases, c2prb, N_coverage\n", " # plot results as array of sub-figures\n", " names = [r\"$\\Delta$\" + pnams[i] for i in range(len(pnams))]\n", " fig, axarr = plotCorrelations(biases, names)\n", " fig.suptitle(\"Biases and Correlations\", size=\"xx-large\")\n", " ax = axarr[0][1]\n", " ax.text(\n", " 0.1,\n", " 0.45,\n", " \"$\\\\Delta$: fit - true \\n \\n\"\n", " + \"$\\\\mu$: mean \\n\"\n", " + \"$\\\\sigma$: standard deviation \\n\"\n", " + \"$\\\\sigma_\\\\mu$: error on mean \\n \\n\"\n", " + \"$\\\\rho$: correlation coefficient\",\n", " transform=ax.transAxes,\n", " )\n", " plt.show()\n", "\n", "\n", "plot_correlations()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Die beobachteten Verteilungen bestätigen die Schätzungen der Unsicherheiten und der Korrelationen durch die in *xyFit()* verwendeten Methoden. Wenn sich deutliche Unterschiede ergeben hätten, würde man die aus der Ensemble-Studie abgeleiteten Unsicherheiten und Korrelationen als Ergebnis angeben." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Verlässlichkeit des Hypothesentests\n", "\n", "Als letztes Beispiel ist unten die Überprüfung der Qualität der als Hypothesentest eingesetzten\n", "und aus dem in der Anpassung erhaltenen Wert von $\\chi^2$ am Optimum und der Zahl der Freiheitsgrade\n", "berechneten $\\chi^2$-Wahrscheinlichkeit gezeigt. Diese sollte flach sein, wenn die Größe zur Bestimmung\n", "der Goodness-of-Fit tatsächlich einer $\\chi^2$-Verteilung folgt. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def analyse_chi2():\n", " global nfail, biases, c2prb, N_coverage\n", " # analysis of chi2 probability\n", " figc2prb = plt.figure(figsize=(7.5, 5.0))\n", " ax = figc2prb.add_subplot(1, 1, 1)\n", " nbins = int(min(50, Nexp / 20))\n", " binc, bine, _ = ax.hist(c2prb, bins=nbins, ec=\"grey\")\n", " plt.xlabel(\"chi2 probability\")\n", " plt.ylabel(\"entries/bin\")\n", " # - check compatibility with flat distribution\n", " mean = (Nexp - nfail) / nbins\n", " chi2flat = np.sum((binc - mean) ** 2 / mean)\n", " prb = 1.0 - stats.chi2.cdf(chi2flat, nbins)\n", " print(\n", " \"compatibility of chi2 probability with flat distribution: {:f}%\".format(\n", " prb * 100\n", " )\n", " )\n", "\n", " plt.show()\n", "\n", "\n", "analyse_chi2()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In diesem Fall ist die beobachtete Verteilung tatsächlich im Rahmen eines ebenfalls mit der\n", "$\\chi^2$-Methode durchgeführten Tests auf Übereinstimmung mit einer Rechteckverteilung flach. Allerdings zeigen sich bei Erhöhung der Zahl der untersuchten Stichproben statistisch\n", "signifikante Abweichungen, die aber noch hinreichend klein sind, so dass man den Wert von\n", "$\\chi^2$ als Qualitätskriterium für die Anpassung verwenden kann. Mit Hilfe der nun\n", "aus dem Ensemble-Test bekannten Verteilung könnte man auch die kritische Grenze für das \n", "Verwerfen der Modellhypothese entsprechend anpassen. \n", "\n", "Übrigens: der Hauptgrund für die Abweichung der $\\chi^2$-Wahrscheinlichkeit von der \n", "Rechteckverteilung sind in diesem Beispiel die relativen, auf den Modellwert bezogenen\n", "Unsicherheiten der Messdaten. Schaltet man diese in der Simulation und im Modell aus,\n", "ergibt sich eine deutlich bessere Übereinstimmung mit der Rechteckverteilung. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Überprüfung der Robustheit\n", "\n", "Die Überprüfung der Robustheit einer Anpassung ist eines der schwierigsten Probleme. \n", "Auch hier helfen Ensemble-Tests, bei denen dann aber Stichproben der Daten erzeugt werden,\n", "die anders aussehen als im Modell angenommen. So können Effekte von falsch abgeschätzten\n", "Unsicherheiten, leicht variierten Modellfunktionen oder von seltenen Ausreißern in \n", "Verteilungen studiert werden. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Studie 1: Behandlung relativer Unsicherheiten\n", "\n", "Mit dem oben definierten Set-up sind nun leicht weitere Studien der Auswirkungen bestimmter \n", "Optionen der Anpassung leicht durchzuführen. Dazu müssen lediglich die entsprechenden Optionen\n", "verändert und die Sequenz der oben durchgeführten Schritte wiederholt werden. \n", "Gezeigt ist das im Beispiel in der nächsten Zelle für den Fall, dass relative Unsicherheiten auf die \n", "Daten statt auf den Modellwert bezogen werden. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# modify fit optoins\n", "kw_options = { # options for Fit\n", " \"ref_to_model\": False, # !!! rel. uncert. w.r.t to data\n", " \"use_negLogL\": True, # standard chi^2 if false\n", "}\n", "# execute MC loop and analysis ...\n", "MC_loop()\n", "# ... and print results\n", "print(\"\\n\\n*==* sucessfully analyzed {:d} data sets\".format(Nexp - nfail))\n", "print_result()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Wie die Ausgabe zeigt, sind die Parameter nun signifikant verzerrt, und im Vergleich zur \n", "Standardabweichung sind diese Verzerrungen auch relevant. die Überdeckung ist ebenfalls \n", "sehr viel schlechter geworden, d.h. die Angaben der Unsicherheiten sind in diesem Fall\n", "nicht zuverlässig. Die grafische Darstellung der Ergebnisse des ersten Pseudodatensatzes\n", "zeigt, was passiert: Datenpunkten, die zu kleineren Werten fluktuieren, werden kleinere \n", "relative Unsicherheiten zugewiesen, und dementsprechend größere Unsicherheiten für Datenpunkte,\n", "die zu größeren Werten fluktuieren. Im Ergebnis wird die angepasste Modellfunktion dadurch hin zu\n", "kleineren $y$-Werten gezogen. Nicht in allen Fällen ist der Effekt so deutlich in der Ergebnisgrafik\n", "sichtbar wie in diesem Fall! Den größten Beitrag zur offensichtlichen Abweichung des angepassten\n", "Modells an die Daten liefert übrigens die relative, korrelierte Unsicherheit von 3% auf die $y$-Werte." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Studie 2: Vergleich mit der Methode der kleinsten Quadrate\n", "\n", "Verwendet man anstelle der Likelihood-Methode die klassische $\\chi^2$-Methode, so ergeben \n", "sich ebenfalls starke Verzerrungen und schlechte Werte für die Abdeckung des Konfidenzintervalls\n", "der Parameter. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# modify fit optoins\n", "kw_options = { # options for Fit\n", " \"ref_to_model\": True, # rel. uncert. w.r.t model\n", " \"use_negLogL\": False, # !!! use standard chi^2\n", "}\n", "# execute MC loop and analysis ...\n", "MC_loop()\n", "print(\"\\n\\n*==* sucessfully analyzed {:d} data sets\".format(Nexp - nfail))\n", "# ... and print results\n", "print_result()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In diesem Fall wird die angepasste Modellfunktion hin zu größeren $y$-Werten gezogen. Der Grund\n", "ist die Vernachlässigung der parameterabhängigen Unsicherheiten in der Kostenfunktion. \n", "Wenn die Modellfunktion größere Werte vorhersagt, werden auch die relativen Unsicherheiten \n", "entsprechend größer und damit der $\\chi^2$-Wert entsprechend kleiner. Der Term $\\ln(\\det(V))$ \n", "in der vollen log-Likelihood wirkt diesem Effekt entgegen. Wieder ist in vielen Fällen das Problem\n", "nicht so klar sichtbar wie in diesem Beispeil mit recht großen Anteilen relativer Unsicherheiten.\n", "Der Test \"mit dem Auge\", d.h. scheinbar gute Übereinstimmung von Datenpunkten und Modell, \n", "täuscht oft." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Resumé:** \n", "Als Ergebnis dieser Studie ist festzuhalten, dass relative Unsicherheiten auf die Modellwerte \n", "bezogen werden müssen und dazu eine Likelielihood als Kostenfunktion notwendig ist.\n", "In den Paketen *kafe2* und *PhyPraKit.phyFit* ist das korrekte Verhalten voreingestellt. \n", "Bei Verwendung anderer Anpassungsprogramme ist Vorsicht geboten, wenn relative \n", "Unsicherheiten im Spiel sind. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "# Abschließende Anmerkungen\n", "\n", "Mit den in diesem Tutorial beschriebenen Verfahren zur Anpassung von Modellen an Messdaten\n", "sind Sie für Anwendungen in den Fortgeschrittenenpraktika oder in Ihren Abschlussarbeiten\n", "bestens gerüstet. Die beschriebenen Methoden entsprechen den in der Wissenschaft üblichen,\n", "und Sie werden in wissenschaftlichen Veröffentlichungen einigen davon wieder begegnen. \n", "\n", "Die hier gezeigten Beispiele basieren mit Absicht auf dem schlanken Anpassungswerkzeug \n", "*phyFit*, um die Transparenz der Vorgehensweise sicher zu stellen und Ihnen durch \n", "Modifikation der Beispiele eigene Studien und Anpassungen an spezielle Problemstellungen\n", "zu ermöglichen. Da *phyFit* nur eine schlanke Interface-Ebene darstellt und im Übrigen auf dem \n", "im wissenschaftlichen Umfeld seit langem etablierten und anerkannten Paket MINUIT zur Optimierung\n", "in hochdimensionalen Parameterräumen und zur Schätzung der Parameterunsicherheiten basiert, sind \n", "Ihre damit erzielten Ergebnisse fachlich sehr gut abgesichert und in eigenen Arbeiten \"zitierfähig\".\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "***" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.12" } }, "nbformat": 4, "nbformat_minor": 4 }