mirror of
https://gitlab.rlp.net/pgp/pgp1-python-einfuehrung
synced 2024-11-16 13:48:11 +00:00
1078 lines
37 KiB
Text
1078 lines
37 KiB
Text
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "5064e2e2",
|
|
"metadata": {},
|
|
"source": [
|
|
"# Fortgeschrittenes Beispiel\n",
|
|
"In diesem Abschnitt wollen wir uns mit einem komplexeren Beispiel beschäftigen, um weitere Methoden von `iminuit` kennzulernen.\n",
|
|
"Hierzu betrachten wir ein Zählexperiment, z.B. ein Teilchendetektor, bei dem ein Energiespektrum aufgenommen wird. Für jedes Energieintervall (Bin) wird die Anzahl der registrierten Ereignisse bestimmt. Hierbei können wir annehmen, dass die Verteilung der gemessenen Anzahl durch eine Poisson-Verteilung beschrieben wird. Dann entspricht der Fehler in jedem Bin gerade $\\sqrt n$. \n",
|
|
"Dieses Spektrum soll aus zwei gauß-förmigen Peaks über einem exponentiellen Untergrund bestehen und wird mit Hilfe eines Zufallszahlengenerator \"erzeugt\"."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "100a4fe4-a5c4-4be3-a7f7-13337b97a194",
|
|
"metadata": {},
|
|
"source": [
|
|
"Nun wollen wir die Messdaten mit Hilfe von `iminuit` fitten. Hierzu müssen wir zunächste zwei Module des packages importieren und eine Funktion für die Entladekurve des Kondensators definieren:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 1,
|
|
"id": "520f4973",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Diese Zelle nur auf JupyterHub des ZDV ausführen um `iminuit` zu installieren!\n",
|
|
"# import sys\n",
|
|
"# import subprocess\n",
|
|
"# subprocess.check_call([\n",
|
|
"# sys.executable, \n",
|
|
"# '-m',\n",
|
|
"# 'pip',\n",
|
|
"# 'install',\n",
|
|
"# '--proxy',\n",
|
|
"# 'http://webproxy.zdv.uni-mainz.de:3128',\n",
|
|
"# 'iminuit'\n",
|
|
"# ])"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 2,
|
|
"id": "2ffe340b-cd0f-45ec-b5b8-42e7a0349d4c",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"from iminuit import Minuit, cost\n",
|
|
"import matplotlib.pyplot as plt\n",
|
|
"import numpy as np"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "143a2a23-0a62-439f-9d28-9f555ae85589",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"rnd_bkd = np.random.exponential(39.7, 5000)\n",
|
|
"rnd_bkd += 40\n",
|
|
"\n",
|
|
"peak1 = np.random.normal(53.3, 2.1, 5000)\n",
|
|
"peak2 = np.random.normal(60.5, 2.78, 12000)\n",
|
|
"data = np.concatenate([rnd_bkd, peak1, peak2])\n",
|
|
"\n",
|
|
"entries, edges = np.histogram(data, bins=120, range=(40, 80))\n",
|
|
"center = edges[:-1] + np.diff(edges)/2\n",
|
|
"\n",
|
|
"plt.errorbar(center, entries, np.sqrt(entries), ls='', marker='.')\n",
|
|
"plt.xlabel('Energy [keV]')\n",
|
|
"plt.ylabel('Number of counts per bin')"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "b582615c-9251-409d-bcfc-d19fd579e161",
|
|
"metadata": {},
|
|
"source": [
|
|
"Zunächst wollen wir das Fitmodel in der Form\n",
|
|
"\n",
|
|
"$$f(x) = A_1 \\cdot \\exp \\bigg\\{\\frac{-(x - \\mu_1)^2}{2 \\cdot \\sigma_1^2}\\bigg\\} + A_2 \\cdot \\exp \\bigg\\{\\frac{-(x - \\mu_2)^2}{2 \\cdot \\sigma_2^2}\\bigg\\} + A_3 \\exp\\{-x/\\tau\\}$$\n",
|
|
"\n",
|
|
"definieren. Hier lohnt es sich, erst Funktionen für die einzelnen Komponenten zu definieren und dann das Gesamtmodel. Hierdurch lassen sich später die einzelnen Komponenten besser darstellen."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 4,
|
|
"id": "f84d7527-c0d2-475d-966d-5363a8e09369",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def peak(x, A, mu, sigma):\n",
|
|
" return A*np.exp(-(x-mu)**2/(2*sigma**2))\n",
|
|
"\n",
|
|
"def bkg(x, A, tau):\n",
|
|
" return A*np.exp(-x/tau)\n",
|
|
"\n",
|
|
"def fit_model(x, A_p1, A_p2, mu_p1, mu_p2, sigma_p1, sigma_p2, A_bkg, tau_bkg):\n",
|
|
" return peak(x, A_p1, mu_p1, sigma_p1) + peak(x, A_p2, mu_p2, sigma_p2) + bkg(x, A_bkg, tau_bkg)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "32014861-316c-4692-9d52-48f2fb71321c",
|
|
"metadata": {},
|
|
"source": [
|
|
"Nun wollen wir wieder die Kostenfunktion und die Minimierungsfunktion definieren. Startwerte können wir anhand unseres Plots ablesen, lediglich $\\tau$ lässt sich auf diese Weise nicht gut bestimmen."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 5,
|
|
"id": "a31901cf-a0ce-4db8-a072-a661fbbb7296",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"ls = cost.LeastSquares(center, entries, np.sqrt(entries), fit_model)\n",
|
|
"\n",
|
|
"mi = Minuit(ls, \n",
|
|
" A_p1 = 400, \n",
|
|
" A_p2 = 700,\n",
|
|
" mu_p1 = 54,\n",
|
|
" mu_p2 = 60,\n",
|
|
" sigma_p1 = 2,\n",
|
|
" sigma_p2 = 2,\n",
|
|
" A_bkg = 100,\n",
|
|
" tau_bkg = 10, \n",
|
|
" )\n",
|
|
"mi.limits['tau_bkg'] = (0, None)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "1e69a046-770f-4c38-9b91-0176bb0686a1",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"plt.errorbar(center, entries, np.sqrt(entries), ls='', marker='.')\n",
|
|
"plt.xlabel('Energy [keV]')\n",
|
|
"plt.ylabel('Number of counts per bin')\n",
|
|
"\n",
|
|
"x = np.arange(40, 80, 0.1)\n",
|
|
"plt.plot(x, fit_model(x, *mi.values), color='k', label='Initial guess')\n",
|
|
"plt.legend()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "89f755f4-b780-43a6-a923-49662c4c701a",
|
|
"metadata": {},
|
|
"source": [
|
|
"Unsere Startparameter sind bereits nicht schlecht, aber weichen noch stark von den Daten ab. Bei komplexeren Daten und Fitmodellen lohnt es sich, den Fit schrittweise durchzuführen. Bevor wir uns den beiden Peaks widmen, welche uns eigentlich interessieren, sollten wir versuchen, den Untergrund etwas besser zu beschreiben. Um den Untergrund besser fitten zu können, sollten wir erst den Fitbereich auf einen Energiebereich limitieren, in welchem der Untergrund dominiert. Dem Plot können wir entnehmen, dass dies für alle Werte unterhalb von 45 keV und oberhalb von 70 keV der Fall ist. Im Allgemeinen können wir Wertebereiche in Python mit Hilfe von „Masken“ selektieren. Eine Maske lässt sich wie folgt erstellen:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 505,
|
|
"id": "d53e8386-ea7f-43fa-b4fe-65229308a2ec",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"mask_outside_of_peaks = (center < 45) | (center >= 70)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "84cef7a6-13a0-4ba8-ac40-eb86a54411dc",
|
|
"metadata": {},
|
|
"source": [
|
|
"Die Maske hat hierbei die Selbe länge wie unseren Daten…"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "d1d06116-d726-4163-b414-6ccde6a19027",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"len(mask_outside_of_peaks), len(mask_outside_of_peaks)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "80db0ae0-5cbd-4db9-b184-610d77bf1c58",
|
|
"metadata": {},
|
|
"source": [
|
|
"… und beinhaltet Wahrheitswerte `True` und `False`, bzw. 1 und 0, mit welchen wir unsere Daten selektieren können:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "f24d19d8-3483-45b5-aee9-1d3f8755da22",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"mask_outside_of_peaks, center[mask_outside_of_peaks]"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "5b5c07e7-1865-48f2-bd9e-0540661fd71e",
|
|
"metadata": {},
|
|
"source": [
|
|
"Unsere Selektion können wir an unsere Kostenfunktion direkt übergeben."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 508,
|
|
"id": "3034bb22-0b96-498d-9736-ed9bb2189460",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"ls.mask = (center < 45) | (center >= 70)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "77a664fd-513e-4c89-ba52-945b6f68512f",
|
|
"metadata": {},
|
|
"source": [
|
|
"Nun können wir nochmal unsere Funktion und Messwerte für den ausgewählten Bereich plotten…"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "81232354-a7b8-4e2a-9ac0-159ce0a03da4",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"plt.errorbar(center[ls.mask], entries[ls.mask], np.sqrt(entries[ls.mask]), ls='', marker='.', label='Not masked')\n",
|
|
"plt.errorbar(center[~ls.mask], entries[~ls.mask], np.sqrt(entries[~ls.mask]), ls='', marker='.', label='Masked')\n",
|
|
"plt.xlabel('Energy [keV]')\n",
|
|
"plt.ylabel('Number of counts per bin')\n",
|
|
"\n",
|
|
"x = np.arange(40, 80, 0.1)\n",
|
|
"plt.plot(x, fit_model(x, *mi.values), color='k', label='Initial guess')\n",
|
|
"plt.legend()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "ec675b22",
|
|
"metadata": {},
|
|
"source": [
|
|
"Außerdem müssen wir noch alle Fitparameter, welche nicht zum Untergrund beitragen, als konstant festhalten"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "4a93a1c2-17df-46c2-b38e-9a509fe16fc7",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"mi.fixed[:] = True\n",
|
|
"mi.fixed[['tau_bkg', 'A_bkg']] = False\n",
|
|
"print (mi.fixed)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "c5a8d247-5b71-42ae-9706-d16192374686",
|
|
"metadata": {},
|
|
"source": [
|
|
"bevor wir die Minmierung starten und das Resultat darstellen."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "3e90c2ed-c282-47c2-b0fe-3063f9545639",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"mi.migrad()\n",
|
|
"mi.hesse()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "0b435af3-73ea-42de-9ab7-6a16ae9dbceb",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"plt.errorbar(center, entries, np.sqrt(entries), ls='', marker='.')\n",
|
|
"plt.xlabel('Energy [keV]')\n",
|
|
"plt.ylabel('Number of counts per bin')\n",
|
|
"\n",
|
|
"x = np.arange(40, 80, 0.1)\n",
|
|
"plt.plot(x, fit_model(x, *mi.values), color='k', label='Initial guess')\n",
|
|
"plt.legend()\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "6def3e2b-5edf-48bb-99b8-2b7fdaae51c5",
|
|
"metadata": {},
|
|
"source": [
|
|
"Das Resultat sieht bereits sehr gut aus. Nun können wir uns den eigentlichen Peaks widmen und starten im Folgenden mit dem kleineren der beiden. Zunächst sollten wir den maskierten Bereich entweder neu definieren oder komplett entfernen."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 513,
|
|
"id": "ebd77c40-6fcd-4881-bc1d-e3ca8ae0bf3b",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"ls.mask = None"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "7850ae53-ae2d-49aa-ac7b-dcef60a2dab7",
|
|
"metadata": {},
|
|
"source": [
|
|
"Außerdem können wir dem Plot entnehmen, dass durch den höheren Untergrund unsere Anfangsstartwerte nicht mehr ganz so gut passen. Diese können wir wie folgt aktualisieren:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 514,
|
|
"id": "823e05a0-516c-4d30-8dc7-5381e0e2e617",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"mi.values['A_p1'] = 700\n",
|
|
"mi.values['sigma_p1'] = 3"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "8648bf00-901e-40dc-ada2-9a6b684e8f31",
|
|
"metadata": {},
|
|
"source": [
|
|
"Nun sollten wir alle Parameter wieder festhalten und nur die Parameter des ersten Peaks freigeben."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "3c83690c-103e-47ff-b18f-13ac763ee87d",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"mi.fixed[:] = True\n",
|
|
"mi.fixed[['A_p1', 'mu_p1', 'sigma_p1']] = False\n",
|
|
"mi.migrad()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "34df75bf-3750-4186-ae12-4f6bb9e49931",
|
|
"metadata": {},
|
|
"source": [
|
|
"Jetzt wiederholen wir das ganze für den zweiten Peak…"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "264a9891-423c-479a-8906-c048aac2fd2e",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"mi.fixed[:] = True\n",
|
|
"mi.fixed[['A_p2', 'mu_p2', 'sigma_p2']] = False\n",
|
|
"mi.migrad()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "32d67543-870f-4bd9-bba4-2d01086c671a",
|
|
"metadata": {},
|
|
"source": [
|
|
"Zum Schluss geben wir wieder alle Parameter frei und führen einen letzten Fit durch. "
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "72d43004-cd80-418a-996a-f1e7a7133ce9",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"mi.fixed[:] = False\n",
|
|
"mi.migrad()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "067fbf6f-14c4-4a46-afb3-71753d06af23",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"plt.errorbar(center, entries, np.sqrt(entries), ls='', marker='.')\n",
|
|
"plt.xlabel('Energy [keV]')\n",
|
|
"plt.ylabel('Number of counts per bin')\n",
|
|
"\n",
|
|
"x = np.arange(40, 80, 0.1)\n",
|
|
"plt.plot(x, fit_model(x, *mi.values), color='k', label='Best fit')\n",
|
|
"plt.plot(x, peak(x, *mi.values['A_p1', 'mu_p1', 'sigma_p1']), color='gray', ls='--', label='Peak 1')\n",
|
|
"plt.plot(x, peak(x, *mi.values['A_p2', 'mu_p2', 'sigma_p2']), color='gray', ls='-.', label='Peak 2')\n",
|
|
"plt.plot(x, bkg(x, *mi.values['A_bkg', 'tau_bkg']), color='gray', label='Background')\n",
|
|
"plt.legend()\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "7ef19633-0947-4568-b537-a1c69e42b7c2",
|
|
"metadata": {},
|
|
"source": [
|
|
"Das Ergebnis sieht sehr gut aus. Alle Kacheln sind grün und die Daten scheinen durch die Funktion gut beschrieben zu werden. Natürlich können wir das gesamte Fitverfahren auch etwas kompakter in einer Zelle darstellen:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "2311f135-8410-4f35-8d58-b9bcef0fed53",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"ls = cost.LeastSquares(center, entries, np.sqrt(entries), fit_model)\n",
|
|
"\n",
|
|
"mi = Minuit(ls, \n",
|
|
" A_p1 = 800, \n",
|
|
" A_p2 = 1400,\n",
|
|
" mu_p1 = 54,\n",
|
|
" mu_p2 = 60,\n",
|
|
" sigma_p1 = 2,\n",
|
|
" sigma_p2 = 2,\n",
|
|
" A_bkg = 100,\n",
|
|
" tau_bkg = 10, \n",
|
|
" )\n",
|
|
"mi.limits['tau_bkg'] = (0, None)\n",
|
|
"mi.fixed[:] = True\n",
|
|
"ls.mask = (center < 45) | (center >= 70)\n",
|
|
"mi.fixed[['tau_bkg', 'A_bkg']] = False\n",
|
|
"mi.migrad()\n",
|
|
"ls.mask = None\n",
|
|
"mi.values['A_p1'] = 700\n",
|
|
"mi.values['sigma_p1'] = 3\n",
|
|
"mi.fixed[:] = True\n",
|
|
"mi.fixed[['A_p1', 'mu_p1', 'sigma_p1']] = False\n",
|
|
"mi.migrad()\n",
|
|
"mi.fixed[:] = True\n",
|
|
"mi.fixed[['A_p2', 'mu_p2', 'sigma_p2']] = False\n",
|
|
"mi.migrad()\n",
|
|
"mi.fixed[:] = False\n",
|
|
"mi.migrad()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "b2d4c8e9-da2c-489e-9b2f-de24f042c341",
|
|
"metadata": {},
|
|
"source": [
|
|
" # Wann fittet ein Fit?\n",
|
|
"Nach dem wir nun unser Model an unsere Daten angepasst haben, stellt sich die Frage: „Spiegelt unser Model unsere Daten gut wider?“. Um diese Frage beantworten zu können, gibt es verschiedene Möglichkeiten, welche wir im Folgenden etwas näher betrachten wollen. \n",
|
|
"## Fit Residual: \n",
|
|
"Schauen wir uns zunächst noch einmal an, wie das Chi-Quadrat definiert ist:\n",
|
|
"$$ \\chi^2 = \\sum_i \\frac{(y_i - \\lambda_i)^2}{\\Delta y_i^2} $$\n",
|
|
"Wir minimieren den Abstand zwischen einem Messwert und unserem Model und gewichten diesen mit den Unsicherheiten unserer Messwerte. Fitresiduen spiegeln genau dies wider. Sie sind definiert als \n",
|
|
"$$ \\frac{(y_i - \\lambda_i)}{\\Delta y_i} $$\n",
|
|
"Für unseren Fit sehen sie wie folgt aus.\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "30cafddc-ea17-4158-82cc-f132dee2c8de",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"residuals = (entries - fit_model(center, *mi.values))/np.sqrt(entries)\n",
|
|
"\n",
|
|
"plt.plot(center, residuals, ls='', marker='.')\n",
|
|
"plt.xlabel('Energy [keV]')\n",
|
|
"plt.ylabel('Residuals [$\\sigma$]')"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "d0ef61ca-afc5-472d-8e8e-b4726ef2a3dd",
|
|
"metadata": {},
|
|
"source": [
|
|
"Als einzelner Plot sind sie noch nicht sehr informativ. Hilfreicher ist es bereits, wenn wir die Residuen zusammen mit unseren Daten und Fitmodel darstellen. "
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "d9fbe83b-3146-4d72-89a4-084c29752e24",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"fig_fit = plt.figure(constrained_layout=True)\n",
|
|
"gs = fig_fit.add_gridspec(5, 5, hspace=0)\n",
|
|
"\n",
|
|
"\n",
|
|
"main_axis = fig_fit.add_subplot(gs[:4, :])\n",
|
|
"res_axis = fig_fit.add_subplot(gs[4:, :], sharex=main_axis)\n",
|
|
"fig_fit.tight_layout()\n",
|
|
"\n",
|
|
"\n",
|
|
"main_axis.errorbar(center, entries, np.sqrt(entries), ls='', marker='.', color='k')\n",
|
|
"\n",
|
|
"main_axis.plot(x, peak(x, *mi.values['A_p1', 'mu_p1', 'sigma_p1']), color='gray', ls='--')\n",
|
|
"main_axis.plot(x, peak(x, *mi.values['A_p2', 'mu_p2', 'sigma_p2']), color='gray', ls='-.')\n",
|
|
"main_axis.plot(x, bkg(x, *mi.values['A_bkg', 'tau_bkg']), color='gray')\n",
|
|
"\n",
|
|
"x = np.arange(40, 80, 0.1)\n",
|
|
"main_axis.plot(x, fit_model(x, *mi.values), color='purple', label='Best fit')\n",
|
|
"main_axis.legend()\n",
|
|
"main_axis.set_ylabel('Number of entries per bin')\n",
|
|
"main_axis.xaxis.set_tick_params(direction='inout')\n",
|
|
"main_axis.tick_params(axis='x', labelcolor=(0, 0, 0, 0))\n",
|
|
"main_axis.set_xlim(40, 80)\n",
|
|
"\n",
|
|
"res_axis.set_xlabel('Energy [keV]')\n",
|
|
"res_axis.set_ylabel('Res [$\\sigma$]')\n",
|
|
"res_axis.set_ylim(-3, 3)\n",
|
|
"res_axis.set_yticks([-2, 0, 2])\n",
|
|
"res_axis.fill_between((40, 80), -1, 1, alpha=0.3, color='purple')\n",
|
|
"res_axis.fill_between((40, 80), -2, 2, alpha=0.3, color='purple')\n",
|
|
"res_axis.axhline(0, color='purple')\n",
|
|
"res_axis.set_xlim(40, 80)\n",
|
|
"res_axis.plot(center, \n",
|
|
" residuals,\n",
|
|
" color='k', marker='.', ls=''\n",
|
|
" )\n",
|
|
"plt.tight_layout()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "dbe65a21-572e-4618-bcd8-78f13e945e8a",
|
|
"metadata": {},
|
|
"source": [
|
|
"Sofern unser Fitmodel unsere Daten gut beschreibt, erwarten wir, dass die Residuen sich Gaußförmig zufällig um den Wert 0 herum verteilen. Dies folgt direkt aus der Annahme, dass sich die Unsicherheiten unserer Messwerte durch eine Gaußverteilung darstellen lassen. Dies können wir direkt überprüfen, sofern wir unsere Residuen in ein Histogramm eintragen. "
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "05e24224-66f7-45ed-99c6-f6d257e2c779",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"plt.hist(residuals, bins=10, range=(-3, 3), histtype='step')\n",
|
|
"plt.xlabel('Residual [$\\sigma$]')\n",
|
|
"plt.ylabel('#Entries per bin')\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "24ce04cc-5234-4326-9a28-7624b9c7d23e",
|
|
"metadata": {},
|
|
"source": [
|
|
"Bzw. den Anteil an Residuen berechnen, welcher innerhalb der 1 $\\sigma$ Umgebung liegt."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "39009321-41f4-49f4-820a-717be277b1b0",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"np.sum(np.abs(residuals) < 1)/len(residuals)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "08579cdf-3b28-4ea2-9c61-6ae62974af51",
|
|
"metadata": {},
|
|
"source": [
|
|
"Zeigen unsere Residuen eine Struktur oder ein systematisches Verhalten, deutet dies auf einen ungenauen Fit oder ein falsches Fitmodel hin. Dies ist im Folgenden gezeigt. "
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "850870af-e546-4d95-b9de-8a4e7b61c241",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"pseudo_data = np.random.normal(0, 2, 5000)\n",
|
|
"\n",
|
|
"fig_fit = plt.figure(constrained_layout=True)\n",
|
|
"gs = fig_fit.add_gridspec(5, 5, hspace=0)\n",
|
|
"\n",
|
|
"main_axis = fig_fit.add_subplot(gs[:4, :])\n",
|
|
"res_axis = fig_fit.add_subplot(gs[4:, :], sharex=main_axis)\n",
|
|
"fig_fit.tight_layout()\n",
|
|
"\n",
|
|
"entries1, edges1, _ = main_axis.hist(pseudo_data, bins=25, range=(-5,5), histtype='step', color='k')\n",
|
|
"center1 = edges1[:-1] + np.diff(edges1)/2\n",
|
|
"\n",
|
|
"residuals1 = (entries1 - peak(center1, 400, 0.2, 2))/np.sqrt(entries1)\n",
|
|
"\n",
|
|
"x = np.arange(-5, 5, 0.1)\n",
|
|
"\n",
|
|
"main_axis.plot(x, peak(x, 400, 0.2, 2), color='purple')\n",
|
|
"main_axis.set_ylabel('Number of entries per bin')\n",
|
|
"main_axis.xaxis.set_tick_params(direction='inout')\n",
|
|
"main_axis.tick_params(axis='x', labelcolor=(0, 0, 0, 0))\n",
|
|
"main_axis.set_xlim(-5, 5)\n",
|
|
"\n",
|
|
"res_axis.set_xlabel('Energy [keV]')\n",
|
|
"res_axis.set_ylabel('Res [$\\sigma$]')\n",
|
|
"res_axis.set_ylim(-3, 3)\n",
|
|
"res_axis.set_yticks([-2, 0, 2])\n",
|
|
"res_axis.fill_between((-5, 5), -1, 1, alpha=0.3, color='purple')\n",
|
|
"res_axis.fill_between((-5, 5), -2, 2, alpha=0.3, color='purple')\n",
|
|
"res_axis.axhline(0, color='purple')\n",
|
|
"res_axis.set_xlim(-5, 5)\n",
|
|
"res_axis.plot(center1, \n",
|
|
" residuals1,\n",
|
|
" color='k', marker='.', ls=''\n",
|
|
" )\n",
|
|
"plt.tight_layout()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "48e95a88-0742-4221-a716-17dacfc02823",
|
|
"metadata": {},
|
|
"source": [
|
|
"Zusätzlich zu den Fit-Residuen bietet das $\\chi^2$ selbst einen Weg, um die „goodness-of-fit“ unseres Model bestimmen zu können ...\n",
|
|
"\n",
|
|
"### $\\chi^2$:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "fe1789cf-7ed3-4db3-a0ae-9e563a9dc85e",
|
|
"metadata": {},
|
|
"source": [
|
|
"Wie gut fittet unsere obige Funktion unsere Messdaten? Sehr gut? Gut? Befriedigend? Oder doch eher schlecht? Wäre es nicht gut, ein Maß für die Güte des Fits zu haben? Wie könnte ein solches Maß aussehen?\n",
|
|
"\n",
|
|
"Sie haben das entscheidende Kriterium bereits kennengelernt: bei der Methode der kleinsten Quadrate geht es darum, das $\\chi^2$ zu minimieren. Gucken wir uns hierzu erst noch einmal an, wie sich das $\\chi^2$ berechnet:\n",
|
|
"\n",
|
|
"$$ \\chi(\\phi_1 ... \\phi_N)^2 = \\sum_{i = 1}^{N} \\frac{ (y_i - \\lambda(x_i; \\phi))^2}{\\Delta y_i^2}$$\n",
|
|
"\n",
|
|
"Bei der Minimierung werden dabei Werte mit geringerer Unsicherheit bevorzugt, d.h. stärker gewichtet (s. Bild unten).\n",
|
|
"\n",
|
|
"<figure class=\"image\">\n",
|
|
"<img src=\"images/MaterialPythonkurs092018/LeastSquare.png\" alt=\"{{ Least Square Beispiel }}\" width=80%>\n",
|
|
"</figure>\n",
|
|
"\n",
|
|
"Damit man für einen gegebenen Datensatz nicht hunderte von verschiedenen Funktionen durchprobieren muss, gibt es für das $\\chi^2$ eine allgemeine Faustregel, welche den berechneten $\\chi^2$-Wert mit der Anzahl unserer Freiheitsgrade vergleicht. Die Anzahl an Freiheitsgrade ist gemeinhin gegeben als *Anzahl der Messwerte - Anzahl der Funktionsparameter* ($m - n$).\n",
|
|
"\n",
|
|
"1. Sofern $\\chi^2/\\text{ndof} >> 1$: sollte die Hypothese bzw. die Fitfunktion angezweifelt werden. Sie beschreibt in diesem Fall die Messdaten nur unzureichend. (Bzw. sollte $\\chi^2/\\text{ndof} > 1$ kann dies auch bedeuten, dass die Unsicherheiten unterschätzt sind)\n",
|
|
"2. Sofern $\\chi^2/\\text{ndof} \\approx 1$: beschreibt die Hypothese bzw. die Fitfunktion die Daten wie erwartet und wird nicht abgelehnt. \n",
|
|
"3. Falls $\\chi^2/\\text{ndof} << 1$ beschreibt die Hypothese bzw. die Fitfunktion die Daten wesentlich besser als erwartet. In diesem Fall heißt das nicht automatisch, dass unsere Hypothese falsch ist, aber man sollte überprüfen, ob die gemessenen Fehler nicht überschätzt worden sind (oder eine Korrelation zwischen den Messfehlern vorliegt). \n",
|
|
"\n",
|
|
"Sofern Sie eine Arbeit schreiben und Ihre **Goodness-of-the-Fit** ($\\chi^2/\\text{ndof}$) angeben wollen, so geben Sie immer beides an, das $\\chi^2$ und die Anzahl an Freiheitsgraden *ndof*. Beide Werte getrennt haben einen größeren Informationsgehalt als der resultierende Quotient (Genaueres lernen Sie z.B. in der Vorlesung *Statistik, Datenanalyse und Simulationen* im Master).\n",
|
|
"\n",
|
|
"Sehen wir uns hierzu nochmal unseren Doppelpeakfit etwas genauer an. `iminuit` berechnet hier für uns bereits das reduzierete $\\chi^2$."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "fa85a19a-f066-4567-abb0-6283ae1bc90b",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"mi"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "9f464246-d333-4143-baf0-aa2a632c5be4",
|
|
"metadata": {},
|
|
"source": [
|
|
"Eine eigene Abschätzung für das $\\chi^2$ ergibt:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "b0ad46ce-f541-40bb-898c-154ad5f94787",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def chi_square_ndof(x_values, y_values, dy_values, fit_model, minuit):\n",
|
|
" ndof = len(x_values) - len(minuit.values)\n",
|
|
" chi2 = np.sum((y_values - fit_model(x_values, *minuit.values))**2/dy_values**2)\n",
|
|
" return chi2, ndof\n",
|
|
"\n",
|
|
"\n",
|
|
"chi_square, ndof = chi_square_ndof(center, entries, np.sqrt(entries), fit_model, mi)\n",
|
|
"print(chi_square, ndof, chi_square/ndof)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "295031f4-6d18-411c-b5dd-a62ed97da7f1",
|
|
"metadata": {},
|
|
"source": [
|
|
"### Hypothesen-Test mittels $\\chi^2$\n",
|
|
"Wie schon im vorherigen Abschnitt erwähnt, kann man das $\\chi^2$ auch dazu verwenden, die Gültigkeit des gewählten Models zu prüfen.\n",
|
|
"Hierzu schauen wir uns die $\\chi^2$-Verteilung an. Der einzige freie Parameter ist die Anzahl der Freiheitsgrade. Die Anzahl der Freiheitsgrade ist auch gleichzeitig der Erwartungswert der $\\chi^2$-Verteilung. In unserem Beispiel oben ist die Anzahl der Freiheitsgrade 112 und die entsprechende Verteilung sieht wie folgt aus..."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 527,
|
|
"id": "8c11bc85-4e25-4d40-8397-257414d48a1f",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"from scipy.stats import chi2\n",
|
|
"# chi_distribution = lambda x, ndof: chi2.pdf(x, ndof)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "76836863-109c-4e7c-989e-04b62ec4ca9d",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"x = np.arange(40., 180.)\n",
|
|
"# plt.plot(x, chi_distribution(x, 112))\n",
|
|
"plt.plot(x,chi2.pdf(x, 112))\n",
|
|
"x = np.arange(chi_square, 180, 0.1)\n",
|
|
"plt.fill_between(x, chi2.pdf(x, 112), alpha=0.3)\n",
|
|
"plt.ylim(0, None)\n",
|
|
"plt.xlabel('x')\n",
|
|
"plt.ylabel('$\\chi^2(x)$')\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "b30829b3-a9e8-4d93-8895-9fd9f67ab9dc",
|
|
"metadata": {},
|
|
"source": [
|
|
"Der erste Schritt für den Hypothesen-Test ist die Berechnung des $P$-Werts\n",
|
|
"$$ P = \\int_{\\chi^2}^{\\infty} f(z,n_d)dz $$\n",
|
|
"wobei $f(z,n_d)$ die $\\chi^2$-Verteilung und $n_d$ die Anzahl der Freiheitsgrade ist.\n",
|
|
"Im Bild oben entspricht dies der ausgefüllten Fläche.\n",
|
|
"\n",
|
|
"Die praktische Berechnung erfolgt mittels der kumulativen Verteilungsfunktion via\n",
|
|
"$$ P = 1 - \\chi^2_{CDF}(x, n_d) $$\n",
|
|
"wobei für $x$ das im Fit bestimmte $\\chi^2$ eingesetzt wird. Die praktische Bedeutung des $P$-Werts ist die Wahrscheinlichkeit bei einer Wiederholung des Experiments in größeres $\\chi^2$ zu erhalten, wenn unser Model die Daten richtig beschreibt und die ermittelten Fitparameter den wahren Werten entsprechen."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "cfa9d88a-eada-49dd-8cb3-73c7dd345c08",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"p_value = lambda x, ndof: 1 - chi2.cdf(x, ndof)\n",
|
|
"p_value(chi_square, ndof), p_value(chi_square*10, ndof*10), p_value(ndof, ndof)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "9cba146a-6309-42d1-92cb-8bdde2da42a2",
|
|
"metadata": {},
|
|
"source": [
|
|
"Kehren wir zu unserem Doppelpeak-Spektrum zurück und änderen das Fitmodell, indem wir statt eines exponentiellen einen konstanten Untergrund annehmen."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "9b91ee55-ac17-4dd6-9827-48677f772096",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def alternative_fit_model(x, A_p1, A_p2, mu_p1, mu_p2, sigma_p1, sigma_p2, c):\n",
|
|
" return peak(x, A_p1, mu_p1, sigma_p1) + peak(x, A_p2, mu_p2, sigma_p2) + c\n",
|
|
"\n",
|
|
"ls = cost.LeastSquares(center, entries, np.sqrt(entries), alternative_fit_model)\n",
|
|
"\n",
|
|
"mi = Minuit(ls, \n",
|
|
" A_p1 = 800, \n",
|
|
" A_p2 = 1400,\n",
|
|
" mu_p1 = 54,\n",
|
|
" mu_p2 = 60,\n",
|
|
" sigma_p1 = 2,\n",
|
|
" sigma_p2 = 2,\n",
|
|
" c = 100, \n",
|
|
" )\n",
|
|
"mi.limits['c'] = (0, None)\n",
|
|
"mi.fixed[:] = True\n",
|
|
"ls.mask = (center < 45) | (center >= 70)\n",
|
|
"mi.fixed[['c']] = False\n",
|
|
"mi.migrad()\n",
|
|
"ls.mask = None\n",
|
|
"mi.values['A_p1'] = 700\n",
|
|
"mi.values['sigma_p1'] = 3\n",
|
|
"mi.fixed[:] = True\n",
|
|
"mi.fixed[['A_p1', 'mu_p1', 'sigma_p1']] = False\n",
|
|
"mi.migrad()\n",
|
|
"mi.fixed[:] = True\n",
|
|
"mi.fixed[['A_p2', 'mu_p2', 'sigma_p2']] = False\n",
|
|
"mi.migrad()\n",
|
|
"mi.fixed[:] = False\n",
|
|
"mi.migrad()\n",
|
|
"mi.hesse()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "c9fbbebc",
|
|
"metadata": {},
|
|
"source": [
|
|
"Diese Änderung ist gering und der Fit scheint die Daten weiterhin zu beschreiben. Allerdings gibt bei kleinen Energien eine deutlich sichtbare Diskrepanz. Dies zeigt sich auch in einem größeren $\\chi^2$-Wert. Wie wirkt sich dies auf den $P$-Wert aus?"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "4aa0f3d9-1d0b-4b4c-b816-2a0cb9ae9793",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"chi_square, ndof = chi_square_ndof(center, entries, np.sqrt(entries), alternative_fit_model, mi)\n",
|
|
"print(chi_square, ndof, chi_square/ndof)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "607ddd33",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"p_value = lambda x, ndof: 1 - chi2.cdf(x, ndof)\n",
|
|
"print(chi_square, ndof)\n",
|
|
"p_value(chi_square, ndof)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "bcb62098-1e8b-4c9f-8aa3-048037f0d21e",
|
|
"metadata": {},
|
|
"source": [
|
|
"Der Fit ist offensichtlich viel schlechter und der $P$-Wert liegt nahe bei null, so dass man dieses Model ausschließen sollte.\n",
|
|
"\n",
|
|
"Was aber, wenn die Änderung nicht so dramatisch ist? Ist ein $P$-Wert von 0,4 besser als 0,2? Nein, das kann man so nicht beantworten. Aber für einen Hypothesen-Test sollten man vorher eine Schwelle festlegen für die Akzeptanz oder Ablehnung des Models.\n",
|
|
"\n",
|
|
"Wie ein solcher Hypothesen-Test aussehen kann, wollen wir im Folgenden betrachten. Hierbei benutzen wir\n",
|
|
"1. ein korrektes Model (Normalverteilung),\n",
|
|
"2. ein korrektes Model mit überschätztem Fehler (10% größer),\n",
|
|
"3. und ein falsches Model (Lorentzverteilung)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 264,
|
|
"id": "c3f1f1d4-4b84-45a1-9d23-4cbb8ba32c8c",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def lorentzian( x, x0, a, gam ):\n",
|
|
" return a * gam**2 / ( gam**2 + ( x - x0 )**2)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "0e3fcfd5",
|
|
"metadata": {},
|
|
"source": [
|
|
"Den Fit der drei Modelle und die Bestimmung des entsprechenden $P$-Werts wiederholen wir 5000-mal um eine ausreichende Statistik zu erhalten."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 7,
|
|
"id": "9667c766",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Diese Zelle nur auf JupyterHub des ZDV ausführen um `tqdm` zu installieren falls es nicht vorhanden sein sollte!\n",
|
|
"# import sys\n",
|
|
"# import subprocess\n",
|
|
"# subprocess.check_call([\n",
|
|
"# sys.executable, \n",
|
|
"# '-m',\n",
|
|
"# 'pip',\n",
|
|
"# 'install',\n",
|
|
"# '--proxy',\n",
|
|
"# 'http://webproxy.zdv.uni-mainz.de:3128',\n",
|
|
"# 'tqdm'\n",
|
|
"# ])"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "c3b58808-f155-4194-b02e-e5f649cb86aa",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"from tqdm.notebook import tqdm\n",
|
|
"\n",
|
|
"res_good_model = []\n",
|
|
"res_overfitting = []\n",
|
|
"res_wrong_model = []\n",
|
|
"\n",
|
|
"def peak(x, A, mu, sigma):\n",
|
|
" return A*np.exp(-(x-mu)**2/(2*sigma**2))\n",
|
|
"\n",
|
|
"def lorentzian( x, x0, a, gam ):\n",
|
|
" return a * gam**2 / ( gam**2 + ( x - x0 )**2)\n",
|
|
"\n",
|
|
"\n",
|
|
"for i in tqdm(range(5000)):\n",
|
|
" \n",
|
|
" test_data = np.random.normal(0, 2, 5000)\n",
|
|
" \n",
|
|
" entries, edges = np.histogram(test_data, bins=25, range=(-4,4))\n",
|
|
" center = edges[:-1] + np.diff(edges)/2\n",
|
|
" \n",
|
|
" ls = cost.LeastSquares(center, entries, np.sqrt(entries), peak)\n",
|
|
" mi = Minuit(ls, \n",
|
|
" mu=0.1,\n",
|
|
" sigma=1.5,\n",
|
|
" A = 300\n",
|
|
" )\n",
|
|
" mi.migrad()\n",
|
|
" \n",
|
|
" chi, ndof = chi_square_ndof(center, entries, np.sqrt(entries), peak, mi)\n",
|
|
" res_good_model.append(p_value(chi, ndof))\n",
|
|
"\n",
|
|
"\n",
|
|
" ls = cost.LeastSquares(center, entries, np.sqrt(entries)*1.1, peak)\n",
|
|
" mi = Minuit(ls, \n",
|
|
" mu=0.1,\n",
|
|
" sigma=1.5,\n",
|
|
" A = 300\n",
|
|
" )\n",
|
|
" mi.migrad()\n",
|
|
" \n",
|
|
" chi, ndof = chi_square_ndof(center, entries, np.sqrt(entries)*1.1, peak, mi)\n",
|
|
" res_overfitting.append(p_value(chi, ndof))\n",
|
|
"\n",
|
|
"\n",
|
|
" ls = cost.LeastSquares(center, entries, np.sqrt(entries), lorentzian)\n",
|
|
" mi = Minuit(ls, \n",
|
|
" x0=0,\n",
|
|
" gam=3,\n",
|
|
" a = 300,\n",
|
|
" )\n",
|
|
" mi.migrad()\n",
|
|
" \n",
|
|
" chi, ndof = chi_square_ndof(center, entries, np.sqrt(entries), lorentzian, mi)\n",
|
|
" res_wrong_model.append(p_value(chi, ndof))\n",
|
|
"\n",
|
|
"res_wrong_model = np.array(res_wrong_model)\n",
|
|
"res_good_model = np.array(res_good_model)\n",
|
|
"res_overfit_model = np.array(res_overfitting)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "7ec4cf79",
|
|
"metadata": {},
|
|
"source": [
|
|
"Die Schwelle des $P$-Werts für den Hypothesen-Test setzen wir auf 0,1, d.h. Ergebnisse mit eine, $P$-Wert $<$ 0,1 werden verworfen, alle anderen akzeptiert."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "f41e0b38-56b6-4f2a-bc75-075f622a2068",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"fig, axes = plt.subplots()\n",
|
|
"axes.hist(res_good_model, bins=25, range=(0, 1), histtype='step', color='purple', label='Good model')\n",
|
|
"axes.hist(res_wrong_model, bins=25, range=(0, 1), histtype='step', color='orange', label='Wrong model')\n",
|
|
"axes.hist(res_overfitting, bins=25, range=(0, 1), histtype='step', color='firebrick', label='Too large uncertainties (10 %)')\n",
|
|
"axes.set_xlabel('p-value')\n",
|
|
"axes.set_ylabel('Number of fits')\n",
|
|
"axes.legend()\n",
|
|
"axes.axvline(0.1, color='k')\n",
|
|
"axes2 = plt.twiny()\n",
|
|
"axes2.set_xlabel('Red. $\\chi^2$')\n",
|
|
"axes2.set_xticks([0.2, 0.5, 0.8], ['> 1', '1', '< 1'])\n",
|
|
"plt.show()\n",
|
|
"\n",
|
|
"axes.set_yscale('log')\n",
|
|
"fig"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "86237b4b",
|
|
"metadata": {},
|
|
"source": [
|
|
"Wie man sieht, wird das falsche Modell nahezu immer verworfen während das richtige Modell meistens nicht verworfen wird. Das Modell mit dem überschätzten Fehler wird sogar häufiger akzeptiert, so dass man hier keine Unterscheidung vornehmen kann."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "fc58ee5c-308c-4479-9236-751d7f158fe5",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"print(f'Fraction of wrong model fits rejected: {np.sum(res_wrong_model<0.1)/len(res_wrong_model):.4f}')\n",
|
|
"print(f'Fraction of good model fits rejected: {np.sum(res_good_model<0.1)/len(res_good_model):.4f}')\n",
|
|
"print(f'Fraction of overfitting model fits rejected: {np.sum(res_overfit_model<0.1)/len(res_overfit_model):.4f}')"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "392f4ef2",
|
|
"metadata": {},
|
|
"source": [
|
|
"Wenn man das Limit für den Hypothesen-Test auf 0,05 festlegt, ändern die Ergebnisse wie folgt:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "d5f5efbe-ef8f-48b0-b27b-166f21cb5a06",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"print(f'Fraction of wrong model fits rejected: {np.sum(res_wrong_model<0.05)/len(res_wrong_model):.4f}')\n",
|
|
"print(f'Fraction of good model fits rejected: {np.sum(res_good_model<0.05)/len(res_good_model):.4f}')\n",
|
|
"print(f'Fraction of overfitting model fits rejected: {np.sum(res_overfit_model<0.05)/len(res_overfit_model):.4f}')"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "de9861f6-7870-4dd8-8366-15e0c7dd5125",
|
|
"metadata": {},
|
|
"source": [
|
|
"Der Hypothesen-Test kann das Modell nicht ablehnen, statt es zu bestätigen!"
|
|
]
|
|
}
|
|
],
|
|
"metadata": {
|
|
"kernelspec": {
|
|
"display_name": "jupyter",
|
|
"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.11"
|
|
}
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 5
|
|
}
|