smd_handson_fitting.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"font-size: 3rem; font-weight: bold;\">SMD Hands-On Schätzen</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/javascript": [
       "$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')\n"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%javascript\n",
    "$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1 id=\"tocheading\">Table of Contents</h1>\n",
    "<div id=\"toc\"></div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>.container { width:100% !important; }</style>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib widget\n",
    "\n",
    "from IPython.core.display import display, HTML\n",
    "display(HTML(\"<style>.container { width:100% !important; }</style>\"))\n",
    "plt.rcParams['figure.figsize'] = (7.5, 5)\n",
    "plt.rcParams['figure.constrained_layout.use'] = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cov_to_corr(cov):\n",
    "    '''Convert covariance to correlation matrix\n",
    "    \n",
    "    Taken from: https://math.stackexchange.com/a/300775/892886\n",
    "    '''\n",
    "    D = np.diag(1 / np.sqrt(np.diag(cov)))\n",
    "    return D @ cov @ D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Methode der kleinsten Quadrate\n",
    "\n",
    "## Analytische Lösung für Linearkombinationen von Funktionen\n",
    "\n",
    "(Analog zu Aufgabe 27 des letzten Übungszettels)\n",
    "\n",
    "Hier werden wir eine Funktion der Form\n",
    "\n",
    "$$\n",
    "f(x) = p_0 + p_1 \\cdot \\sin(x) + p_2 \\cdot \\cos(x)\n",
    "$$\n",
    "\n",
    "an Daten anpassen.\n",
    "\n",
    "Diese Funktion ist eine Linearkombination von Basisfunktionen:\n",
    "$$\n",
    "f(x) = \\sum_{i=0}^2 p_i f_i(x)\n",
    "$$\n",
    "\n",
    "mit \n",
    "\n",
    "$$\n",
    "f_0(x) = 1, \\quad f_1(x) = \\sin(x), \\quad f_2(x) = \\cos(x)\n",
    "$$\n",
    "\n",
    "Für diesen Fall gilt die in der Vorlesung und den Übungen vorgestellte analytische Lösung für das Problem der Kleinsten Quadrate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a8d65c14988b4d3a8ccc60a641504620",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'y')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def linear_combination(x, funcs, parameters):\n",
    "    '''Evaluate a liner combination of base functions\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    x: number or np.ndarray\n",
    "        The point or points at which to evaluate\n",
    "    funcs: iterable of callables\n",
    "        The base functions\n",
    "    parameters: iterable of numbers\n",
    "        The coefficients\n",
    "    '''\n",
    "    if len(funcs) != len(parameters):\n",
    "        raise ValueError('len(funcs) must be equal to len(parameters)')\n",
    "        \n",
    "    return np.sum([p * f(x) for p, f in zip(parameters, funcs)], axis=0)\n",
    "\n",
    "\n",
    "funcs = [np.ones_like, np.sin, np.cos]\n",
    "\n",
    "# create some randomized example data points\n",
    "rng = np.random.default_rng(1337)\n",
    "N = 100\n",
    "true_parameters = np.array([2, 1, 0.5])\n",
    "y_unc = rng.uniform(0.1, 0.4, N)\n",
    "\n",
    "x = np.linspace(0, 4 * np.pi, N)\n",
    "y = linear_combination(x, funcs, true_parameters)\n",
    "\n",
    "# Messunsicherheit simulieren\n",
    "y += rng.normal(0, y_unc)\n",
    "\n",
    "plt.figure()\n",
    "plt.errorbar(x, y, yerr=y_unc, ls='', marker='.')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Aufstellen der Design Matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.        , 0.        , 1.        ],\n",
       "       [1.        , 0.12659245, 0.99195481],\n",
       "       [1.        , 0.25114799, 0.9679487 ],\n",
       "       [1.        , 0.37166246, 0.92836793],\n",
       "       [1.        , 0.48619674, 0.87384938]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.column_stack([f(x) for f in funcs])\n",
    "A[:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Aufstellen der Gewichtsmatrix $\\boldsymbol{W} = \\mathrm{Cov}^{-1}(\\mathbf{y})$ der Messwerte.\n",
    "\n",
    "Wir gehen hier davon aus, dass wir die Unsicherheit jedes Messpunktes kennen und das keine Korrelation zwischen den Punkten existiert. Dies ist eine starke Vereinfachung.\n",
    "\n",
    "Die Methode der kleinsten Quadrate erzeugt verzerrte Ergebnisse, falls die Unsicherheiten aus den Messpunkten selbst geschätzt werden. Mehr dazu weiter unten in dem Myon-Beispiel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "92bb6d6830ff4eff86f02543e85386dd",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7f1237c2e6d0>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# All measurements have a known uncertainty and no correlations\n",
    "W = np.diag(np.full(N, 1 / y_unc**2)) \n",
    "plt.matshow(W)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Berechnung der Lösung und der Kovarianz der Parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([2.00565538, 0.97702516, 0.48411763]),\n",
       " array([[ 3.66470863e-04, -5.58492439e-05, -1.00389068e-05],\n",
       "        [-5.58492439e-05,  6.79396825e-04,  4.87671239e-05],\n",
       "        [-1.00389068e-05,  4.87671239e-05,  7.92958991e-04]]))"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cov = np.linalg.inv(A.T @ W @ A)\n",
    "\n",
    "parameters = cov @ A.T @ W @ y\n",
    "parameters, cov"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Chi-Quadrad-Über-Anzahl-Freiheitsgrade"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9440998449492631"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "residuals = (y - A @ parameters)\n",
    "sum_residuals = (residuals.T @ W @ residuals)\n",
    "\n",
    "ndf = len(y) - len(parameters)\n",
    "chisquare_ndf = sum_residuals / ndf\n",
    "chisquare_ndf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ergebnis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e4b5ed630af24079b264d2c5c2cdd8fb",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f12373bed60>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_fit = np.linspace(0, 4 * np.pi, 1000)\n",
    "y_fit = linear_combination(x_fit, funcs, parameters)\n",
    "y_truth = linear_combination(x_fit, funcs, true_parameters)\n",
    " \n",
    "\n",
    "plt.figure()\n",
    "plt.title(\n",
    "    rf'$f(x) = {parameters[0]:.2f} + {parameters[1]:.2f} \\cdot \\sin(x) + {parameters[2]:.2f} \\cdot \\cos(x)'\n",
    "    rf', \\quad \\chi^2_\\mathrm{{ndf}} = {chisquare_ndf:.2f}$'\n",
    ")\n",
    "plt.errorbar(x, y, yerr=y_unc, ls='', marker='.', label='Daten')\n",
    "plt.plot(x_fit, y_fit, label='Fit-Ergebnis')\n",
    "plt.plot(x_fit, y_truth, label='Wahrheit')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "62b96504e4b44bb99056fa7eaaa05df8",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7f1237287be0>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "corr = cov_to_corr(cov)\n",
    "\n",
    "fig = plt.figure(constrained_layout=True)\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "m = ax.matshow(corr, cmap='RdBu_r', vmin=-1, vmax=1)\n",
    "fig.colorbar(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerische Lösung für allgemeine Funktionen\n",
    "\n",
    "\n",
    "Ist die anzupassende Funktion *keine* Linearkombination einzelner Funktionen, kann die Lösung\n",
    "nur numerisch gefunden werden.\n",
    "\n",
    "Aus dem Praktikum bekannt ist die Funktion `scipy.optimize.curve_fit`, die genau dies tut.\n",
    "\n",
    "Der verwendete Nelder-Mead-Algorithmus hat die Eigenschaft, garantiert die korrekte analytische Lösung zu finden, falls sie existiert.\n",
    "\n",
    "- https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html\n",
    "\n",
    "### Anwendung auf lineares Beispiel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2.00565538 0.97702516 0.48411763]\n",
      "[2.00565538 0.97702516 0.48411763]\n",
      "[[ 3.66470863e-04 -5.58492439e-05 -1.00389068e-05]\n",
      " [-5.58492439e-05  6.79396825e-04  4.87671239e-05]\n",
      " [-1.00389068e-05  4.87671239e-05  7.92958991e-04]]\n",
      "[[ 3.66470865e-04 -5.58492432e-05 -1.00389090e-05]\n",
      " [-5.58492432e-05  6.79396818e-04  4.87671201e-05]\n",
      " [-1.00389090e-05  4.87671201e-05  7.92958998e-04]]\n"
     ]
    }
   ],
   "source": [
    "from scipy.optimize import curve_fit\n",
    "\n",
    "\n",
    "def func(x, p1, p2, p3):\n",
    "    return p1 + p2 * np.sin(x) + p3 * np.cos(x)\n",
    "\n",
    "\n",
    "# absolute_sigma prevents scaling of errors to match χ²/ndf=1\n",
    "parameters_numeric, cov_numeric = curve_fit(\n",
    "    func, x, y,\n",
    "    sigma=np.full(N, y_unc),\n",
    "    absolute_sigma=True,\n",
    ")\n",
    "\n",
    "print(parameters, parameters_numeric, cov, cov_numeric, sep='\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0c6a58c3d2424719b7da1e32af1d6c94",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "(0.0, 15.0)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_fit = np.linspace(0, 4 * np.pi, 1000)\n",
    "y_fit = linear_combination(x_fit, funcs, parameters)\n",
    "y_num = linear_combination(x_fit, funcs, parameters_numeric)\n",
    "y_truth = linear_combination(x_fit, funcs, true_parameters)\n",
    " \n",
    "\n",
    "plt.figure()\n",
    "plt.errorbar(x, y, yerr=y_unc, ls='', marker='.', label='Daten')\n",
    "plt.plot(x_fit, y_fit, label='Fit-Ergebnis Analytisch')\n",
    "plt.plot(x_fit, y_num, label='Fit-Ergebnis Numerisch', linestyle=':')\n",
    "plt.plot(x_fit, y_truth, label='Wahrheit')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.legend(loc='upper right')\n",
    "plt.xlim(0, 15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Warum dann nicht einfach immer die numerische Variante nehmen?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "45.2 µs ± 1.78 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "A = np.column_stack([f(x) for f in funcs])\n",
    "\n",
    "ATW = A.T @ W\n",
    "cov = np.linalg.inv(ATW @ A)\n",
    "parameters = cov @ ATW @ y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "427 µs ± 22 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "parameters_numeric, cov_numeric = curve_fit(func, x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Beispiel mit allgemeiner Funktion \n",
    "\n",
    "Aus dem Praktikum, Spaltfunktion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "70eaf7d6b57a469284beae73e10a5225",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df = pd.read_csv('resources/spalt.csv')\n",
    "\n",
    "\n",
    "LASER_WAVELENGTH_NM = 632.8e-9\n",
    "\n",
    "def theory(phi, A0, b):\n",
    "    return (A0 * b * np.sinc(b * np.sin(phi) / LASER_WAVELENGTH_NM))**2\n",
    "\n",
    "\n",
    "# first try with default initial guess (1 for every parameter)\n",
    "p0 = None\n",
    "\n",
    "# now with an \"educated guess\" based on the data and knowledge of the\n",
    "# order of magnitude of the slit size\n",
    "p0 = [np.sqrt(df['I'].max()) / 1e-4, 1e-4]\n",
    "\n",
    "params, cov = curve_fit(theory, df['phi'], df['I'], p0=p0)\n",
    "\n",
    "\n",
    "plt.figure()\n",
    "x = np.linspace(-0.03, 0.03, 501)\n",
    "plt.plot(x, theory(x, *params), label='Fit')\n",
    "\n",
    "plt.plot(df['phi'], df['I'], '.', label='Daten')\n",
    "\n",
    "plt.xlabel(r'$\\varphi \\,\\, / \\,\\, \\mathrm{rad}$')\n",
    "plt.ylabel(r'$I \\,\\, / \\,\\, \\mathrm{A}$')\n",
    "plt.legend(loc='best')\n",
    "\n",
    "None # to not have mpl objects in the output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Maximimum-Likelihood-Methode\n",
    "\n",
    "## Ungebinnter Fit von Wahrscheinlichkeitsdichten\n",
    "\n",
    "\n",
    "Stark vereinfachtes Beispiel einer CERN-Analyse.\n",
    "\n",
    "Wir suchen den Massenpeak eines Teilchens (normal-verteilt) und haben einen exponential-verteilten Untergrund.\n",
    "\n",
    "Wir erzeugen wieder einen \"Spielzeug\"-Datensatz mit Monte-Carlo Methoden:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "152b02667c654aa9874a75ba48b8f24e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rng = np.random.default_rng(42)\n",
    "\n",
    "E_MIN = 75\n",
    "E_MAX = 175\n",
    "\n",
    "# normally distributed signal\n",
    "higgs_signal = np.random.normal(126, 5, 500)\n",
    "\n",
    "# exponentially distributed background\n",
    "background = np.random.exponential(50, size=20000)\n",
    "\n",
    "# combine signal and background\n",
    "measured = np.append(higgs_signal, background)\n",
    "\n",
    "# remove events outside of \"detector range\"\n",
    "in_range = (E_MIN <= measured) & (measured <= E_MAX)\n",
    "measured = measured[in_range]\n",
    "\n",
    "\n",
    "plt.figure()\n",
    "plt.hist(measured, bins=100)\n",
    "plt.xlabel('$m \\,/\\, \\mathrm{GeV}$')\n",
    "plt.margins(x=0)\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Definition der negativen Log-Likelihood\n",
    "\n",
    "Wir bilden eine Überlagerung von zwei Wahrscheinlichkeitsdichten, die einen Anteil von $p$ bzw. $1-p$ haben.\n",
    "\n",
    "Wir müssen die Wahrscheinlichkeitsdichten auf das gemessene Intervall normieren.\n",
    "In diesem Fall ignorieren wir dies für die Normalverteilung unter der Annahme, dass der Massenpeak vollständig im Intervall enthalten ist.\n",
    "\n",
    "Also:\n",
    "\n",
    "\\begin{align}\n",
    "P_1 &= N(\\mu, \\sigma) \\\\[1ex]\n",
    "P_2 &= \\frac{1}{\\exp(-E_\\mathrm{min} / \\tau) - \\exp(-E_\\mathrm{max} / \\tau)} \\exp(- E / \\tau) \\\\[1ex]\n",
    "P(E | p, \\mu, \\sigma, \\tau) &= p \\cdot P_1(E, \\mu, \\sigma) + (1 - p) P_2(E | \\tau)) \\\\[1ex]\n",
    "\\mathcal{L}(p, \\mu, \\sigma, \\tau) &= \\prod_i P(E_i | p, \\mu, \\sigma, \\tau) \\\\[1ex]\n",
    "-\\log\\mathcal{L}(p, \\mu, \\sigma, \\tau) &= -\\sum_i \\log(P(E_i | p, \\mu, \\sigma, \\tau))\n",
    "\\end{align}\n",
    "\n",
    "Im Code, unter Verwendung der Verteilungs-Klassen aus `scipy.stats`, sieht die PDF wie folgt aus:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import norm, expon\n",
    "\n",
    "\n",
    "def pdf(x, mean, std, tau, p):\n",
    "    expon_norm = np.exp(-E_MIN / tau) - np.exp(-E_MAX / tau)\n",
    "    return (\n",
    "        p * norm.pdf(x, mean, std) \n",
    "        + (1 - p) / expon_norm * expon.pdf(x, scale=tau)\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Lösung mit `scipy.optimize.minimize` und `numdifftools.Hessian`\n",
    "\n",
    "\n",
    "Mit `scipy.optimize.minimize` können allgemeine Funktionen minimiert werden.\n",
    "\n",
    "Diese müssen als erstes Argument ein array der zu variierenden Parameter haben.\n",
    "Weitere Parameter können über das `args` Argument von `minimize` übergeben werden.\n",
    "\n",
    "Naiv sieht dann unsere negative Log-Likelihood für die Minimierung mit scipy so aus:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def neg_log_likelihood(parameters, data):\n",
    "    return -np.sum(np.log(pdf(data, *parameters)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Naiv, weil häufig analytische Umformungen der der Log-Likelihood möglich sind,  \n",
    "die zu weniger Operationen führen und somit sowohl Rechenzeit sparen als auch Rundungsfehler vermeiden.\n",
    "\n",
    "Es macht also in vielen Fällen Sinn, die Log-Likelihood per Hand aufzuschreiben und so weit es geht zu vereinfachen. Zum Beispiel auch durch das Weglassen von Termen die unter Variation der Parameter konstant sind."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wir finden ein Minimum (und hoffentlich das globale) mit `scipy.optimize.minimize`. \n",
    "\n",
    "Wie aus der Vorlesung bekannt, bekommen wir einen Schätzer für die Covarianz unserer Parameter aus der Inversen der Hesse-Matrix der Likelihood, ausgewertet im Minimum.\n",
    "\n",
    "Die Hesse-Matrix muss numerisch bestimmt werden. Leider ist die Hesse-Matrix die `scipy.optimize.minimize` mit dem Ergebnis zusammen zurück gibt nur eine grobe Schätzung.\n",
    "\n",
    "Wir nutzen `numdifftools.Hessian` um numerisch die Hesse-Matrix zu bestimmen.\n",
    "\n",
    "`scipy.optimize.minimize` unterstützt mehrere Minimierungs-Algorithmen mit unterschiedlichen Fähigkeiten und Stärken. \n",
    "Für den Fall, dass man Grenzen für die Parameter (`bounds=`) angibt, ist der default der `L-BFGS-B` Algorithmus.\n",
    "\n",
    "Dokumentation: \n",
    "* https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html\n",
    "* https://docs.scipy.org/doc/scipy/reference/optimize.minimize-lbfgsb.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "      fun: 18997.57951127667\n",
       " hess_inv: <4x4 LbfgsInvHessProduct with dtype=float64>\n",
       "      jac: array([-0.00254659, -0.00400178, -0.0003638 ,  0.08258212])\n",
       "  message: 'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'\n",
       "     nfev: 205\n",
       "      nit: 33\n",
       "     njev: 41\n",
       "   status: 0\n",
       "  success: True\n",
       "        x: array([1.25858411e+02, 4.61753387e+00, 4.80818199e+01, 1.04995478e-01])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.optimize import minimize\n",
    "\n",
    "# a very small number to achive > 0 instead of >= 0\n",
    "eps = np.finfo(np.float64).eps\n",
    "\n",
    "result = minimize(\n",
    "    neg_log_likelihood,\n",
    "    args=(measured, ),\n",
    "    x0=[130, 2, 30, 0.2], # here, the initial guess is required\n",
    "    bounds=[\n",
    "        (0, None),    # mean >= 0\n",
    "        (eps, None),  # std > 0\n",
    "        (eps, None),  # tau > 0\n",
    "        (0, 1),       # 0 <= p <= 1\n",
    "    ],\n",
    ")\n",
    "\n",
    "result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-20-c83882007162>:2: RuntimeWarning: invalid value encountered in log\n",
      "  return -np.sum(np.log(pdf(data, *parameters)))\n"
     ]
    }
   ],
   "source": [
    "from numdifftools import Hessian\n",
    "\n",
    "hesse = Hessian(neg_log_likelihood)\n",
    "cov = np.linalg.inv(hesse(result.x, measured))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Higgs mass is 125.86 ± 0.43 GeV\n"
     ]
    }
   ],
   "source": [
    "higgs_mass = result.x[0]\n",
    "higgs_mass_unc = np.sqrt(cov[0, 0])\n",
    "print(f'Higgs mass is {higgs_mass:.2f} ± {higgs_mass_unc:.2} GeV')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Lösung mit iminuit\n",
    "\n",
    "\n",
    "Iminuit ist ein Python Paket, dass Zugriff auf den Minimierungsalgorithmus \"Minuit\" aus dem ROOT Framework ermöglicht.\n",
    "\n",
    "Hierzu ist keine Installation von ROOT notwendig. \n",
    "\n",
    "\"Minuit\" gilt – zumindest unter Teilchenphysiker:innen – als das \"Non-Plus-Ultra\" der Minimierer.\n",
    "\n",
    "Siehe dazu auch: https://www.babushk.in/posts/adventures-physics-software-minuit.html \n",
    "\n",
    "`iminuit` stellt den Minimierer und Loss-Funktionen zur Verfügung, die es wesentlich einfacher machen,  \n",
    "verschiedene Arten von Fits durchzuführen.\n",
    "\n",
    "Außerdem berechnet es die Unsicherheiten über die Hesse-Matrix und stellt auch über das \"minos\" Interface Likelihood-Scans zur genaueren Bestimmung von Unsicherheiten über Likelihood-Ratio Tests zur Verfügung (dazu mehr in den nächsten Vorlesungen zum Thema Testen). \n",
    "\n",
    "\n",
    "\n",
    "Außerdem stellt es – gerade im Notebook – schöne Übersichten über den Fit-Vorgang zur Verfügung.\n",
    "\n",
    "Für den Beginn einmal das gleiche Problem wie oben, nun gelöst mit `imuit`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "    <tr>\n",
       "        <th colspan=\"5\" style=\"text-align:center\" title=\"Minimizer\"> Migrad </th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td colspan=\"2\" style=\"text-align:left\" title=\"Minimum value of function\"> FCN = 3.8e+04 </td>\n",
       "        <td colspan=\"3\" style=\"text-align:center\" title=\"No. of function evaluations in last call and total number\"> Nfcn = 361 </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td colspan=\"2\" style=\"text-align:left\" title=\"Estimated distance to minimum and goal\"> EDM = 3.71e-06 (Goal: 0.0002) </td>\n",
       "        <td colspan=\"3\" style=\"text-align:center\" title=\"No. of gradient evaluations in last call and total number\">  </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td colspan=\"2\" style=\"text-align:center;background-color:#92CCA6;color:black\"> Valid Minimum </td>\n",
       "        <td colspan=\"3\" style=\"text-align:center;background-color:#92CCA6;color:black\"> No Parameters at limit </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td colspan=\"2\" style=\"text-align:center;background-color:#92CCA6;color:black\"> Below EDM threshold (goal x 10) </td>\n",
       "        <td colspan=\"3\" style=\"text-align:center;background-color:#92CCA6;color:black\"> Below call limit </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td style=\"text-align:center;background-color:#92CCA6;color:black\"> Covariance </td>\n",
       "        <td style=\"text-align:center;background-color:#92CCA6;color:black\"> Hesse ok </td>\n",
       "        <td style=\"text-align:center;background-color:#92CCA6;color:black\" title=\"Is covariance matrix accurate?\"> Accurate </td>\n",
       "        <td style=\"text-align:center;background-color:#92CCA6;color:black\" title=\"Is covariance matrix positive definite?\"> Pos. def. </td>\n",
       "        <td style=\"text-align:center;background-color:#92CCA6;color:black\" title=\"Was positive definiteness enforced by Minuit?\"> Not forced </td>\n",
       "    </tr>\n",
       "</table><table>\n",
       "    <tr>\n",
       "        <td></td>\n",
       "        <th title=\"Variable name\"> Name </th>\n",
       "        <th title=\"Value of parameter\"> Value </th>\n",
       "        <th title=\"Hesse error\"> Hesse Error </th>\n",
       "        <th title=\"Minos lower error\"> Minos Error- </th>\n",
       "        <th title=\"Minos upper error\"> Minos Error+ </th>\n",
       "        <th title=\"Lower limit of the parameter\"> Limit- </th>\n",
       "        <th title=\"Upper limit of the parameter\"> Limit+ </th>\n",
       "        <th title=\"Is the parameter fixed in the fit\"> Fixed </th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> 0 </th>\n",
       "        <td> mean </td>\n",
       "        <td> 125.9 </td>\n",
       "        <td> 0.4 </td>\n",
       "        <td> -0.4 </td>\n",
       "        <td> 0.4 </td>\n",
       "        <td> 0 </td>\n",
       "        <td>  </td>\n",
       "        <td>  </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> 1 </th>\n",
       "        <td> std </td>\n",
       "        <td> 4.6 </td>\n",
       "        <td> 0.5 </td>\n",
       "        <td> -0.4 </td>\n",
       "        <td> 0.5 </td>\n",
       "        <td> 2.22E-16 </td>\n",
       "        <td>  </td>\n",
       "        <td>  </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> 2 </th>\n",
       "        <td> tau </td>\n",
       "        <td> 48.1 </td>\n",
       "        <td> 1.5 </td>\n",
       "        <td> -1.5 </td>\n",
       "        <td> 1.6 </td>\n",
       "        <td> 2.22E-16 </td>\n",
       "        <td>  </td>\n",
       "        <td>  </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> 3 </th>\n",
       "        <td> p </td>\n",
       "        <td> 0.105 </td>\n",
       "        <td> 0.009 </td>\n",
       "        <td> -0.009 </td>\n",
       "        <td> 0.010 </td>\n",
       "        <td> 0 </td>\n",
       "        <td> 1 </td>\n",
       "        <td>  </td>\n",
       "    </tr>\n",
       "</table><table>\n",
       "    <tr>\n",
       "        <td></td>\n",
       "        <th colspan=\"2\" style=\"text-align:center\" title=\"Parameter name\"> mean </th>\n",
       "        <th colspan=\"2\" style=\"text-align:center\" title=\"Parameter name\"> std </th>\n",
       "        <th colspan=\"2\" style=\"text-align:center\" title=\"Parameter name\"> tau </th>\n",
       "        <th colspan=\"2\" style=\"text-align:center\" title=\"Parameter name\"> p </th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th title=\"Lower and upper minos error of the parameter\"> Error </th>\n",
       "        <td> -0.4 </td>\n",
       "        <td> 0.4 </td>\n",
       "        <td> -0.4 </td>\n",
       "        <td> 0.5 </td>\n",
       "        <td> -1.5 </td>\n",
       "        <td> 1.6 </td>\n",
       "        <td> -0.009 </td>\n",
       "        <td> 0.010 </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th title=\"Validity of lower/upper minos error\"> Valid </th>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> True </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> True </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> True </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> True </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> True </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> True </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> True </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> True </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th title=\"Did scan hit limit of any parameter?\"> At Limit </th>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th title=\"Did scan hit function call limit?\"> Max FCN </th>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th title=\"New minimum found when doing scan?\"> New Min </th>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "    </tr>\n",
       "</table><table>\n",
       "    <tr>\n",
       "        <td></td>\n",
       "        <th> mean </th>\n",
       "        <th> std </th>\n",
       "        <th> tau </th>\n",
       "        <th> p </th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> mean </th>\n",
       "        <td> 0.181 </td>\n",
       "        <td style=\"background-color:rgb(239,239,250);color:black\"> -0.016 <strong>(-0.082)</strong> </td>\n",
       "        <td style=\"background-color:rgb(244,244,250);color:black\"> -0.0317 <strong>(-0.049)</strong> </td>\n",
       "        <td style=\"background-color:rgb(239,239,250);color:black\"> -0.000334 <strong>(-0.083)</strong> </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> std </th>\n",
       "        <td style=\"background-color:rgb(239,239,250);color:black\"> -0.016 <strong>(-0.082)</strong> </td>\n",
       "        <td> 0.212 </td>\n",
       "        <td style=\"background-color:rgb(222,222,250);color:black\"> -0.15 <strong>(-0.213)</strong> </td>\n",
       "        <td style=\"background-color:rgb(250,166,166);color:black\"> 0.00243 <strong>(0.559)</strong> </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> tau </th>\n",
       "        <td style=\"background-color:rgb(244,244,250);color:black\"> -0.0317 <strong>(-0.049)</strong> </td>\n",
       "        <td style=\"background-color:rgb(222,222,250);color:black\"> -0.15 <strong>(-0.213)</strong> </td>\n",
       "        <td> 2.34 </td>\n",
       "        <td style=\"background-color:rgb(212,212,250);color:black\"> -0.0042 <strong>(-0.290)</strong> </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> p </th>\n",
       "        <td style=\"background-color:rgb(239,239,250);color:black\"> -0.000334 <strong>(-0.083)</strong> </td>\n",
       "        <td style=\"background-color:rgb(250,166,166);color:black\"> 0.00243 <strong>(0.559)</strong> </td>\n",
       "        <td style=\"background-color:rgb(212,212,250);color:black\"> -0.0042 <strong>(-0.290)</strong> </td>\n",
       "        <td> 8.92e-05 </td>\n",
       "    </tr>\n",
       "</table>"
      ],
      "text/plain": [
       "┌─────────────────────────────────────────────────────────────────────────┐\n",
       "│                                Migrad                                   │\n",
       "├──────────────────────────────────┬──────────────────────────────────────┤\n",
       "│ FCN = 3.8e+04                    │              Nfcn = 361              │\n",
       "│ EDM = 3.71e-06 (Goal: 0.0002)    │                                      │\n",
       "├──────────────────────────────────┼──────────────────────────────────────┤\n",
       "│          Valid Minimum           │        No Parameters at limit        │\n",
       "├──────────────────────────────────┼──────────────────────────────────────┤\n",
       "│ Below EDM threshold (goal x 10)  │           Below call limit           │\n",
       "├───────────────┬──────────────────┼───────────┬─────────────┬────────────┤\n",
       "│  Covariance   │     Hesse ok     │ Accurate  │  Pos. def.  │ Not forced │\n",
       "└───────────────┴──────────────────┴───────────┴─────────────┴────────────┘\n",
       "┌───┬──────┬───────────┬───────────┬────────────┬────────────┬─────────┬─────────┬───────┐\n",
       "│   │ Name │   Value   │ Hesse Err │ Minos Err- │ Minos Err+ │ Limit-  │ Limit+  │ Fixed │\n",
       "├───┼──────┼───────────┼───────────┼────────────┼────────────┼─────────┼─────────┼───────┤\n",
       "│ 0 │ mean │   125.9   │    0.4    │    -0.4    │    0.4     │    0    │         │       │\n",
       "│ 1 │ std  │    4.6    │    0.5    │    -0.4    │    0.5     │2.22045e-16│         │       │\n",
       "│ 2 │ tau  │   48.1    │    1.5    │    -1.5    │    1.6     │2.22045e-16│         │       │\n",
       "│ 3 │ p    │   0.105   │   0.009   │   -0.009   │   0.010    │    0    │    1    │       │\n",
       "└───┴──────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘\n",
       "┌──────────┬───────────────────────┬───────────────────────┬───────────────────────┬───────────────────────┐\n",
       "│          │         mean          │          std          │          tau          │           p           │\n",
       "├──────────┼───────────┬───────────┼───────────┬───────────┼───────────┬───────────┼───────────┬───────────┤\n",
       "│  Error   │   -0.4    │    0.4    │   -0.4    │    0.5    │   -1.5    │    1.6    │  -0.009   │   0.010   │\n",
       "│  Valid   │   True    │   True    │   True    │   True    │   True    │   True    │   True    │   True    │\n",
       "│ At Limit │   False   │   False   │   False   │   False   │   False   │   False   │   False   │   False   │\n",
       "│ Max FCN  │   False   │   False   │   False   │   False   │   False   │   False   │   False   │   False   │\n",
       "│ New Min  │   False   │   False   │   False   │   False   │   False   │   False   │   False   │   False   │\n",
       "└──────────┴───────────┴───────────┴───────────┴───────────┴───────────┴───────────┴───────────┴───────────┘\n",
       "┌──────┬─────────────────────────────────────────┐\n",
       "│      │      mean       std       tau         p │\n",
       "├──────┼─────────────────────────────────────────┤\n",
       "│ mean │     0.181    -0.016   -0.0317 -0.000334 │\n",
       "│  std │    -0.016     0.212     -0.15   0.00243 │\n",
       "│  tau │   -0.0317     -0.15      2.34   -0.0042 │\n",
       "│    p │ -0.000334   0.00243   -0.0042  8.92e-05 │\n",
       "└──────┴─────────────────────────────────────────┘"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from iminuit import Minuit\n",
    "from iminuit.cost import UnbinnedNLL\n",
    "\n",
    "# minuit's UnbinnedNLL takes directly the pdf and the observed data\n",
    "loss = UnbinnedNLL(measured, pdf)\n",
    "\n",
    "m = Minuit(loss, mean=130, std=2, tau=30, p=0.2)\n",
    "\n",
    "# set bounds\n",
    "m.limits['mean'] = (0, None)   # >= 0\n",
    "m.limits['std'] = (eps, None)  # > 0\n",
    "m.limits['tau'] = (eps, None)  # > 0\n",
    "m.limits['p'] = (0, 1)\n",
    "\n",
    "# perform minimization\n",
    "m.migrad()\n",
    "\n",
    "# perform likelihood scan for confidence intervals\n",
    "m.minos()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Higgs mass is 125.86 ± 0.43 GeV\n"
     ]
    }
   ],
   "source": [
    "higgs_mass = m.values['mean']\n",
    "higgs_mass_unc = m.errors['mean']\n",
    "\n",
    "print(f'Higgs mass is {higgs_mass:.2f} ± {higgs_mass_unc:.2} GeV')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Poisson-Likelihood-Fit einer histogrammierten Ereignisverteilung\n",
    "\n",
    "\n",
    "Stehen die einzelnen Beobachtungen nicht zur Verfügung oder ist bei der Analyse von großen Datenmengen die Laufzeit relevant, kann auch ein sogenannter \"gebinnter\" Fit durchgeführt werden.\n",
    "\n",
    "Hier wird die Wahrscheinlichts-Verteilung an das Histogramm angepasst.\n",
    "\n",
    "Da ein Histogramm in jedem Bin ein Zählexperiment darstellt, erwarten wir eine Poisson-Statistik für jedes Bin.\n",
    "\n",
    "Die kumulierte Wahrscheinlichkeitsdichte liefert zusammen mit der Gesamtereignis-Anzahl den Erwartungswert für jedes Histogramm-Bin, abhängig von den zu fittenden Parametern $\\boldsymbol{\\theta}$.\n",
    "\n",
    "$$\n",
    "\\mathcal{L} = \\prod_{i=1}^N \\mathcal{P}(k=H_i, \\lambda=\\lambda_i(\\boldsymbol{\\theta}))\n",
    "$$\n",
    "\n",
    "mit \n",
    "\n",
    "$$\n",
    "\\lambda_i = N_\\mathrm{total} \\cdot (\\mathrm{CDF}(b_i, \\boldsymbol{\\theta}) - \\mathrm{CDF}(a_i, \\boldsymbol{\\theta}))\n",
    "$$\n",
    "\n",
    "Wobei $a_i$ und $b_i$ die Bin-Grenzen des $i$-ten Bins sind. Wir integrieren also die PDF in den einzelnen Bins und multiplizieren mit der Gesamtzahl der Ereignisse.\n",
    "\n",
    "\n",
    "**Beispiel aus Lebensdauer kosmische Myonen (V01)**\n",
    "\n",
    "Im Experiment wird direkt in Hardware ein Histogramm gemessen. \n",
    "\n",
    "Die einzelnen Lebensdauern stehen also nicht zur Verfügung.\n",
    "\n",
    "Die Anleitung fordert einen Least-Squares-Fit an die Bin-Höhen. \n",
    "Dies liefert einen Erwartungstreuen Schätzer, so lange ein ungewichteter Fit durchgeführt wird.\n",
    "\n",
    "In der Vergangenheit wurde empfohlen, einen gewichteten Fit unter der Annahme $\\sigma_i = \\sqrt{H_i}$ durchzuführen, richtig ist $\\sqrt{\\lambda}$.\n",
    "\n",
    "Dies ist falsch! Diese Methode liefert konsistent verzerrte Ergebnisse, da Unterfluktuationen stärker gewichtet werden als Überfluktuationen.\n",
    "\n",
    "Die korrekte Methode ist der gebinnte Poisson-Likelihood Fit wie hier oder ein iteratives Least-Squares-Verfahren wie es in SMD-2 vorgestellt wird.\n",
    "\n",
    "Vergleich der Methoden hier: https://gist.github.com/maxnoe/41730e6ca1fac01fc06f0feab5c3566d\n",
    "\n",
    "In dem Versuch kommt es zusätzlich zu der erwarteten Exponential-Verteilung der Myon-Zerfälle noch zu einem gleichverteilten Untergrund durch koinzidente, nicht-zerfallende Myonen. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "81b04aff59d1442f98a46e6d7ccaeaf7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = np.genfromtxt(\"resources/muon_data.txt\")\n",
    "t = np.arange(len(N) + 1) / 21.48\n",
    "\n",
    "t = t[5:]\n",
    "N = N[5:]\n",
    "\n",
    "plt.figure()\n",
    "plt.stairs(values=N, edges=t)\n",
    "plt.xlabel('t / µs')\n",
    "plt.ylabel('Ereignisanzahl')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Iminuit stellt für diesen Anwendungsfall die `BinnedNLL` Loss-Funktion zur Verfügung,\n",
    "die genau die Poisson-Likelihood für ein gemessenes Histogram darstellt."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "    <tr>\n",
       "        <th colspan=\"5\" style=\"text-align:center\" title=\"Minimizer\"> Migrad </th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td colspan=\"2\" style=\"text-align:left\" title=\"Minimum value of function\"> FCN = 247.5 (chi2/ndof = 1.1) </td>\n",
       "        <td colspan=\"3\" style=\"text-align:center\" title=\"No. of function evaluations in last call and total number\"> Nfcn = 103 </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td colspan=\"2\" style=\"text-align:left\" title=\"Estimated distance to minimum and goal\"> EDM = 1.62e-08 (Goal: 0.0002) </td>\n",
       "        <td colspan=\"3\" style=\"text-align:center\" title=\"No. of gradient evaluations in last call and total number\">  </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td colspan=\"2\" style=\"text-align:center;background-color:#92CCA6;color:black\"> Valid Minimum </td>\n",
       "        <td colspan=\"3\" style=\"text-align:center;background-color:#92CCA6;color:black\"> No Parameters at limit </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td colspan=\"2\" style=\"text-align:center;background-color:#92CCA6;color:black\"> Below EDM threshold (goal x 10) </td>\n",
       "        <td colspan=\"3\" style=\"text-align:center;background-color:#92CCA6;color:black\"> Below call limit </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td style=\"text-align:center;background-color:#92CCA6;color:black\"> Covariance </td>\n",
       "        <td style=\"text-align:center;background-color:#92CCA6;color:black\"> Hesse ok </td>\n",
       "        <td style=\"text-align:center;background-color:#92CCA6;color:black\" title=\"Is covariance matrix accurate?\"> Accurate </td>\n",
       "        <td style=\"text-align:center;background-color:#92CCA6;color:black\" title=\"Is covariance matrix positive definite?\"> Pos. def. </td>\n",
       "        <td style=\"text-align:center;background-color:#92CCA6;color:black\" title=\"Was positive definiteness enforced by Minuit?\"> Not forced </td>\n",
       "    </tr>\n",
       "</table><table>\n",
       "    <tr>\n",
       "        <td></td>\n",
       "        <th title=\"Variable name\"> Name </th>\n",
       "        <th title=\"Value of parameter\"> Value </th>\n",
       "        <th title=\"Hesse error\"> Hesse Error </th>\n",
       "        <th title=\"Minos lower error\"> Minos Error- </th>\n",
       "        <th title=\"Minos upper error\"> Minos Error+ </th>\n",
       "        <th title=\"Lower limit of the parameter\"> Limit- </th>\n",
       "        <th title=\"Upper limit of the parameter\"> Limit+ </th>\n",
       "        <th title=\"Is the parameter fixed in the fit\"> Fixed </th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> 0 </th>\n",
       "        <td> tau </td>\n",
       "        <td> 2.169 </td>\n",
       "        <td> 0.034 </td>\n",
       "        <td> -0.034 </td>\n",
       "        <td> 0.035 </td>\n",
       "        <td> 2.22E-16 </td>\n",
       "        <td>  </td>\n",
       "        <td>  </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> 1 </th>\n",
       "        <td> p </td>\n",
       "        <td> 0.958 </td>\n",
       "        <td> 0.007 </td>\n",
       "        <td> -0.007 </td>\n",
       "        <td> 0.007 </td>\n",
       "        <td> 0 </td>\n",
       "        <td> 1 </td>\n",
       "        <td>  </td>\n",
       "    </tr>\n",
       "</table><table>\n",
       "    <tr>\n",
       "        <td></td>\n",
       "        <th colspan=\"2\" style=\"text-align:center\" title=\"Parameter name\"> tau </th>\n",
       "        <th colspan=\"2\" style=\"text-align:center\" title=\"Parameter name\"> p </th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th title=\"Lower and upper minos error of the parameter\"> Error </th>\n",
       "        <td> -0.034 </td>\n",
       "        <td> 0.035 </td>\n",
       "        <td> -0.007 </td>\n",
       "        <td> 0.007 </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th title=\"Validity of lower/upper minos error\"> Valid </th>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> True </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> True </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> True </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> True </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th title=\"Did scan hit limit of any parameter?\"> At Limit </th>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th title=\"Did scan hit function call limit?\"> Max FCN </th>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th title=\"New minimum found when doing scan?\"> New Min </th>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "        <td style=\"background-color:#92CCA6;color:black\"> False </td>\n",
       "    </tr>\n",
       "</table><table>\n",
       "    <tr>\n",
       "        <td></td>\n",
       "        <th> tau </th>\n",
       "        <th> p </th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> tau </th>\n",
       "        <td> 0.00119 </td>\n",
       "        <td style=\"background-color:rgb(250,131,131);color:black\"> 0.000201 <strong>(0.790)</strong> </td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <th> p </th>\n",
       "        <td style=\"background-color:rgb(250,131,131);color:black\"> 0.000201 <strong>(0.790)</strong> </td>\n",
       "        <td> 5.47e-05 </td>\n",
       "    </tr>\n",
       "</table>"
      ],
      "text/plain": [
       "┌─────────────────────────────────────────────────────────────────────────┐\n",
       "│                                Migrad                                   │\n",
       "├──────────────────────────────────┬──────────────────────────────────────┤\n",
       "│ FCN = 247.5 (chi2/ndof = 1.1)    │              Nfcn = 103              │\n",
       "│ EDM = 1.62e-08 (Goal: 0.0002)    │                                      │\n",
       "├──────────────────────────────────┼──────────────────────────────────────┤\n",
       "│          Valid Minimum           │        No Parameters at limit        │\n",
       "├──────────────────────────────────┼──────────────────────────────────────┤\n",
       "│ Below EDM threshold (goal x 10)  │           Below call limit           │\n",
       "├───────────────┬──────────────────┼───────────┬─────────────┬────────────┤\n",
       "│  Covariance   │     Hesse ok     │ Accurate  │  Pos. def.  │ Not forced │\n",
       "└───────────────┴──────────────────┴───────────┴─────────────┴────────────┘\n",
       "┌───┬──────┬───────────┬───────────┬────────────┬────────────┬─────────┬─────────┬───────┐\n",
       "│   │ Name │   Value   │ Hesse Err │ Minos Err- │ Minos Err+ │ Limit-  │ Limit+  │ Fixed │\n",
       "├───┼──────┼───────────┼───────────┼────────────┼────────────┼─────────┼─────────┼───────┤\n",
       "│ 0 │ tau  │   2.169   │   0.034   │   -0.034   │   0.035    │2.22045e-16│         │       │\n",
       "│ 1 │ p    │   0.958   │   0.007   │   -0.007   │   0.007    │    0    │    1    │       │\n",
       "└───┴──────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘\n",
       "┌──────────┬───────────────────────┬───────────────────────┐\n",
       "│          │          tau          │           p           │\n",
       "├──────────┼───────────┬───────────┼───────────┬───────────┤\n",
       "│  Error   │  -0.034   │   0.035   │  -0.007   │   0.007   │\n",
       "│  Valid   │   True    │   True    │   True    │   True    │\n",
       "│ At Limit │   False   │   False   │   False   │   False   │\n",
       "│ Max FCN  │   False   │   False   │   False   │   False   │\n",
       "│ New Min  │   False   │   False   │   False   │   False   │\n",
       "└──────────┴───────────┴───────────┴───────────┴───────────┘\n",
       "┌─────┬───────────────────┐\n",
       "│     │      tau        p │\n",
       "├─────┼───────────────────┤\n",
       "│ tau │  0.00119 0.000201 │\n",
       "│   p │ 0.000201 5.47e-05 │\n",
       "└─────┴───────────────────┘"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from iminuit.cost import BinnedNLL\n",
    "from scipy.stats import uniform\n",
    "\n",
    "T_MIN, T_MAX = t[0], t[-1]\n",
    "\n",
    "def cdf(x, tau, p):\n",
    "    # normalize to 1 in histogram range\n",
    "    cdf_min, cdf_max = expon.cdf([T_MIN, T_MAX], scale=tau) \n",
    "    norm = 1 / (cdf_max - cdf_min)\n",
    "    \n",
    "    # combine exponential signal with uniform background\n",
    "    return (\n",
    "        p * expon.cdf(x, scale=tau) * norm\n",
    "        + (1 - p) * uniform.cdf(x, T_MIN, T_MAX)\n",
    "    )\n",
    "\n",
    "\n",
    "# histogram counds, histogram edges and cumulative distribution function\n",
    "loss = BinnedNLL(N, t, cdf)\n",
    "\n",
    "m = Minuit(loss, tau=2, p=0.99)\n",
    "m.limits['tau'] = (eps, None)\n",
    "m.limits['p'] = (0, 1)\n",
    "m.migrad()\n",
    "m.minos()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fit: τ = 2.169 ± 0.034 µs\n",
      "Lit: τ = 2.1969811 ± 0.0000022 µs\n"
     ]
    }
   ],
   "source": [
    "muon_lifetime = m.values['tau']\n",
    "muon_lifetime_unc = m.errors['tau']\n",
    "\n",
    "print(f'Fit: τ = {muon_lifetime:.3f} ± {muon_lifetime_unc:.3f} µs')\n",
    "print(f'Lit: τ = 2.1969811 ± 0.0000022 µs')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Likelihood-Scan für die Unsicherheiten:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "tau, ts, valid = m.mnprofile('tau', size=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c7fc27e089ff4f8c806f463d3f50378d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "\n",
    "plt.axvline(m.values['tau'], color='C1')\n",
    "plt.plot(m.values['tau'], m.fval, color='C1', marker='o', zorder=3)\n",
    "\n",
    "plt.axvline(m.values['tau'] + m.merrors['tau'].lower, color='C2')\n",
    "plt.axvline(m.values['tau'] + m.merrors['tau'].upper, color='C2')\n",
    "plt.axhline(m.fval + 1, color='C3')\n",
    "\n",
    "plt.plot(tau, ts)\n",
    "\n",
    "plt.xlabel('τ / µs')\n",
    "plt.ylabel('NLL')\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "252c580714b74d58ae5adbad68d1a6ee",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.contour.ContourSet at 0x7f16f4ac08b0>"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plt.figure()\n",
    "m.draw_mncontour('tau', 'p', size=250)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}