Week2/LikelihoodFit/LikelihoodFit_original.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Principle of Maximum Likelihood\n",
    "\n",
    "\n",
    "## Description:\n",
    "\n",
    "Python script for illustrating the principle of maximum likelihood and a likelihood fit.\n",
    "\n",
    "__This is both an exercise, but also an attempt to illustrate four things:__\n",
    "  1. How to make a (binned and unbinned) Likelihood function/fit.\n",
    "  2. The difference and a comparison between a Chi-square and a (binned) Likelihood.\n",
    "  3. The difference and a comparison between a binned and unbinned Likelihood.\n",
    "  4. What goes on behind the scenes in Minuit, when it is asked to fit something.\n",
    "\n",
    "In this respect, the exercise is more of an illustration rather than something to be used directly, which is why it is followed later by another exercise, where you can test if you have understood the differences, and how and when to apply which fit method.\n",
    "\n",
    "The example uses 50 exponentially distributed random times, with the goal of finding the best estimate of the lifetime (data is generated with lifetime, tau = 1). Three estimates are considered:\n",
    "  1. Chi-square fit (chi2)\n",
    "  2. Binned Likelihood fit (bllh)\n",
    "  3. Unbinned Likelihood fit (ullh)\n",
    "\n",
    "The three methods are based on a scan of values for tau in the range [0.5, 2.0]. For each value of tau, the chi2, bllh, and ullh are calculated. In the two likelihood cases, it is actually -2*log(likelihood) which is calculated, which you should (by now) understand why.\n",
    " \n",
    "Note that the unbinned likelihood is in principle the \"optimal\" fit, but also the most difficult for several reasons (convergence, numerical problems, implementation, speed, etc.). However, all three methods/constructions essentially yield the same results, when there is enough statistics (i.e. errors are Gaussian), though the $\\chi^2$ also gives a fit quality.\n",
    " \n",
    "The problem is explicitly chosen to have only one fit parameter, such that simple 1D graphs can show what goes on. In this case, the analytical solution (simple mean) is actually prefered (see Barlow). Real world problems will almost surely be more complex.\n",
    "\n",
    "Also, the exercise is mostly for illustration. In reality, one would hardly ever calculate and plot the Chi-square or Likelihood values, but rather do the minimization using an algorithm (Minuit) to do the hard work.\n",
    "\n",
    "### Reference:\n",
    "- Barlow, chapter 5 (5.1-5.7)\n",
    "- Cowan, chapter 6\n",
    "\n",
    "### Authors: \n",
    "- Troels C. Petersen (Niels Bohr Institute, petersen@nbi.dk)\n",
    "\n",
    "### Date:    \n",
    "- 23-11-2023 (latest update)\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np                                     # Matlab like syntax for linear algebra and functions\n",
    "import matplotlib.pyplot as plt                        # Plots and figures like you know them from Matlab\n",
    "import seaborn as sns                                  # Make the plots nicer to look at\n",
    "from iminuit import Minuit                             # The actual fitting tool, better than scipy's\n",
    "import sys                                             # Module to see files and folders in directories\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.path.append('../../External_Functions')\n",
    "from ExternalFunctions import Chi2Regression, BinnedLH, UnbinnedLH\n",
    "from ExternalFunctions import nice_string_output, add_text_to_ax # useful functions to print fit results on figure\n",
    "\n",
    "plt.rcParams['font.size'] = 16 # set some basic plotting parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Program settings:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "save_plots = False     # Determining if plots are saved or not\n",
    "verbose = True         # Should the program print or not?\n",
    "veryverbose = False    # Should the program print a lot or not?\n",
    "\n",
    "ScanChi2 = True        # In addition to fit for minimum, do a scan...\n",
    "\n",
    "# Parameters of the problem:\n",
    "Ntimes = 50           # Number of time measurements.\n",
    "tau_truth = 1.0;       # We choose (like Gods!) the lifetime.\n",
    "\n",
    "# Binning:\n",
    "Nbins = 50             # Number of bins in histogram\n",
    "tmax = 10.0            # Maximum time in histogram\n",
    "binwidth = tmax / Nbins      # Size of bins (s)\n",
    "\n",
    "# General settings:\n",
    "r = np.random          # Random numbers\n",
    "r.seed(42)             # We set the numbers to be random, but the same for each run"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Generate data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Produce array of exponentially distributed times and put them in a histogram:\n",
    "t = r.exponential(tau_truth, Ntimes) # Exponential with lifetime tau.\n",
    "yExp, xExp_edges = np.histogram(t, bins=Nbins, range=(0, tmax))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Is the data plotted like we wouls like to? Let's check..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# In case you want to check that the numbers really come out as you want to (very healthy to do at first):\n",
    "if (veryverbose) :\n",
    "    for index, time in enumerate(t) :\n",
    "        print(f\"  {index:2d}:   t = {time:5.3f}\")\n",
    "        if index > 10: \n",
    "            break # let's restrain ourselves"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks like values are coming int, but are they actually giving an exponential? Remember the importance of __plotting your data before hand__!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_center = xExp_edges[:-1] + (xExp_edges[1]-xExp_edges[0])/2.0   # Get the value of the histogram bin centers\n",
    "plt.plot(X_center,yExp,'o')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check that it looks like you are producing the data that you want. If this is the case, move on (and possibly comment out the plot!)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analyse data:\n",
    "The following is \"a manual fit\", i.e. scanning over possible values of the fitting parameter(s) - here luckely only one, tau - and seeing what value of chi2, bllh, and ullh it yields. When plotting these, one should find a <b>parabola</b>, the minimum value of which is the optimal fitting parameter of tau. The rate of increase around this minimum represents the uncertainty of the fitting parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Define the number of tau values and their range to test in Chi2 and LLH:\n",
    "# As we know the \"truth\", namely tau = 1, the range [0.5, 1.5] seems fitting for the mean.\n",
    "# The number of bins can be increased at will, but for now 50 seems fitting.\n",
    "Ntau_steps = 50\n",
    "min_tau = 0.5\n",
    "max_tau = 1.5\n",
    "delta_tau = (max_tau-min_tau) / Ntau_steps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Loop over hypothesis for the value of tau and calculate Chi2 and (B)LLH:\n",
    "chi2_minval = 999999.9   # Minimal Chi2 value found\n",
    "chi2_minpos = 0.0        # Position (i.e. time) of minimal Chi2 value\n",
    "bllh_minval = 999999.9\n",
    "bllh_minpos = 0.0\n",
    "ullh_minval = 999999.9\n",
    "ullh_minpos = 0.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tau  = np.zeros(Ntau_steps+1)\n",
    "chi2 = np.zeros(Ntau_steps+1)\n",
    "bllh = np.zeros(Ntau_steps+1)\n",
    "ullh = np.zeros(Ntau_steps+1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now loop of POSSIBLE tau estimates:\n",
    "for itau in range(Ntau_steps+1):\n",
    "    tau_hypo = min_tau + itau*delta_tau         # Scan in values of tau\n",
    "    tau[itau] = tau_hypo\n",
    "\n",
    "    # Calculate Chi2 and binned likelihood (from loop over bins in histogram):\n",
    "    chi2[itau] = 0.0\n",
    "    bllh[itau] = 0.0\n",
    "    for ibin in range (Nbins) :\n",
    "        # Note: The number of EXPECTED events is the intergral over the bin!\n",
    "        xlow_bin = xExp_edges[ibin]\n",
    "        xhigh_bin = xExp_edges[ibin+1]\n",
    "        # Given the start and end of the bin, we calculate the INTEGRAL over the bin,\n",
    "        # to get the expected number of events in that bin:\n",
    "        nexp = Ntimes * (np.exp(-xlow_bin/tau_hypo) - np.exp(-xhigh_bin/tau_hypo))\n",
    "        # The observed number of events... that is just the data!\n",
    "        nobs = yExp[ibin]\n",
    "\n",
    "        if (nobs > 0):      # For ChiSquare but not LLH, we need to require Nobs > 0, as we divide by this:\n",
    "            chi2[itau] += (nobs-nexp)**2 / nobs                           # Chi2 summation/function\n",
    "        bllh[itau] += -2.0*np.log(stats.poisson.pmf(int(nobs), nexp))     # Binned LLH function\n",
    "\n",
    "        if (veryverbose and itau == 0) :\n",
    "            print(f\"  Nexp: {nexp:10.7f}   Nobs: {nobs:3.0f}     Chi2: {chi2[itau]:5.1f}    BLLH: {bllh[itau]:5.1f}\")\n",
    "\n",
    "    # Calculate Unbinned likelihood (from loop over events):\n",
    "    ullh[itau] = 0.0\n",
    "    for time in t :     # i.e. for every data point generated...\n",
    "        ullh[itau] += -2.0*np.log(1.0/tau_hypo*np.exp(-time/tau_hypo))   # Unbinned LLH function\n",
    "    \n",
    "    if (verbose) :\n",
    "        print(f\" {itau:3d}:  tau = {tau_hypo:4.2f}   chi2 = {chi2[itau]:6.2f}   log(bllh) = {bllh[itau]:6.2f}   log(ullh) = {ullh[itau]:6.2f}\")\n",
    "\n",
    "    # Search for minimum values of chi2, bllh, and ullh:\n",
    "    if (chi2[itau] < chi2_minval) :\n",
    "        chi2_minval = chi2[itau]\n",
    "        chi2_minpos = tau_hypo\n",
    "    if (bllh[itau] < bllh_minval) :\n",
    "        bllh_minval = bllh[itau]\n",
    "        bllh_minpos = tau_hypo\n",
    "    if (ullh[itau] < ullh_minval) :\n",
    "        ullh_minval = ullh[itau]\n",
    "        ullh_minpos = tau_hypo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"  Decay time of minimum found:   chi2: {chi2_minpos:7.4f}s    bllh: {bllh_minpos:7.4f}s    ullh: {ullh_minpos:7.4f}s\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"  Chi2 value at minimum:   chi2 = {chi2_minval:.1f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot and fit results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Define range around minimum to be fitted:\n",
    "min_fit = 0.15\n",
    "max_fit = 0.20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2, 2, figsize=(16, 12))\n",
    "\n",
    "ax_chi2 = axes[0,0]\n",
    "ax_bllh = axes[1,0]\n",
    "ax_ullh = axes[0,1]\n",
    "# A fourth plot is available for plotting whatever you want :)\n",
    "\n",
    "# ChiSquare:\n",
    "# ----------\n",
    "ax_chi2.plot(tau, chi2, 'k.', label='chi2')\n",
    "ax_chi2.set_xlim(chi2_minpos-2*min_fit, chi2_minpos+2*max_fit)\n",
    "ax_chi2.set_title(\"ChiSquare\")\n",
    "ax_chi2.set_xlabel(r\"Value of $\\tau$\")\n",
    "ax_chi2.set_ylabel(\"Value of ChiSquare\")\n",
    "\n",
    "# Binned Likelihood:\n",
    "# ----------\n",
    "ax_bllh.plot(tau, bllh,'bo')\n",
    "ax_bllh.set_xlim(bllh_minpos-2*min_fit, bllh_minpos+2*max_fit)\n",
    "ax_bllh.set_title(\"Binned Likelihood\")\n",
    "ax_bllh.set_xlabel(r\"Value of $\\tau$\")\n",
    "ax_bllh.set_ylabel(r\"Value of $\\ln{LLH}$\")\n",
    "\n",
    "# Unbinned Likelihood:\n",
    "# ----------\n",
    "ax_ullh.plot(tau, ullh, 'g.')\n",
    "ax_ullh.set_xlim(ullh_minpos-2*min_fit, ullh_minpos+2*max_fit)\n",
    "ax_ullh.set_title(\"Unbinned Likelihood\")\n",
    "ax_ullh.set_xlabel(r\"Value of $\\tau$\")\n",
    "ax_ullh.set_ylabel(r\"Value of $\\ln{LLH}$\")\n",
    "\n",
    "fig;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Parabola function\n",
    "Note that the parabola is defined differently than normally. The parameters are:\n",
    "   * `minval`:    Minimum value (i.e. constant)\n",
    "   * `minpos`:    Minimum position (i.e. x of minimum)\n",
    "   * `quadratic`: Quadratic term."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func_para(x, minval, minpos, quadratic) :\n",
    "    return minval + quadratic*(x-minpos)**2\n",
    "func_para_vec = np.vectorize(func_para)           # Note: This line makes it possible to send vectors through the function! "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Double parabola with different slopes on each side of the minimum:\n",
    "In case the uncertainties are asymmetric, the parabola will also be so, and hence needs to be fitted with two separate parabolas meeting at the top point. Parameters are now as follows:\n",
    "   * `minval`:   Minimum value (i.e. constant)\n",
    "   * `minpos`:   Minimum position (i.e. x of minimum)\n",
    "   * `quadlow`:  Quadratic term on lower side\n",
    "   * `quadhigh`: Quadratic term on higher side"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func_asympara(x, minval, minpos, quadlow, quadhigh) :\n",
    "    if (x < minpos) :\n",
    "        return minval + quadlow*(x-minpos)**2\n",
    "    else :\n",
    "        return minval + quadhigh*(x-minpos)**2\n",
    "func_asympara_vec = np.vectorize(func_asympara)   # Note: This line makes it possible to send vectors through the function! "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Perform both fits:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fit chi2 values with our parabola:\n",
    "indexes = (tau>chi2_minpos-min_fit) & (tau<chi2_minpos+max_fit)\n",
    "\n",
    "# Fit with parabola:\n",
    "chi2_object_chi2 = Chi2Regression(func_para, tau[indexes], chi2[indexes])\n",
    "minuit_chi2 = Minuit(chi2_object_chi2, minval=chi2_minval, minpos=chi2_minpos, quadratic=20.0)\n",
    "minuit_chi2.errordef = 1.0\n",
    "minuit_chi2.migrad()\n",
    "\n",
    "# Fit with double parabola:\n",
    "chi2_object_chi2_doublep = Chi2Regression(func_asympara, tau[indexes], chi2[indexes])\n",
    "minuit_chi2_doublep = Minuit(chi2_object_chi2_doublep, minval=chi2_minval, minpos=chi2_minpos, quadlow=20.0, quadhigh=20.0)\n",
    "minuit_chi2_doublep.errordef = 1.0\n",
    "minuit_chi2_doublep.migrad();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot (simple) fit:\n",
    "minval, minpos, quadratic = minuit_chi2.values       # Note how one can \"extract\" the three values from the object.\n",
    "print(minval)\n",
    "minval_2p, minpos_2p, quadlow_2p, quadhigh_2p = minuit_chi2_doublep.values\n",
    "print(minval_2p)\n",
    "x_fit = np.linspace(chi2_minpos-min_fit, chi2_minpos+max_fit, 1000)\n",
    "y_fit_simple = func_para_vec(x_fit, minval, minpos, quadratic)\n",
    "ax_chi2.plot(x_fit, y_fit_simple, 'b-')\n",
    "d = {'Chi2 value':        minval,\n",
    "     'Fitted tau (s)':    minpos,\n",
    "     'quadratic':         quadratic}\n",
    "\n",
    "text = nice_string_output(d, extra_spacing=3, decimals=3)\n",
    "add_text_to_ax(0.02, 0.95, text, ax_chi2, fontsize=14)\n",
    "fig.tight_layout()\n",
    "if save_plots: \n",
    "    fig.savefig(\"FitMinimum.pdf\", dpi=600)\n",
    "fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Given the parabolic fit, we can now extract the uncertainty on tau (think about why the below formula works!):\n",
    "err = 1.0 / np.sqrt(quadratic)\n",
    "# For comparison, I give one extra decimal, than I would normally do:\n",
    "print(f\"  Chi2 fit gives:    tau = {minpos:.3f} +- {err:.3f}\")\n",
    "\n",
    "# For the asymmetric case, there are naturally two errors to calculate.\n",
    "#err_lower = 1.0 / np.sqrt(quadlow)\n",
    "#err_upper = 1.0 / np.sqrt(quadhigh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Go through tau values to find minimum and +-1 sigma:\n",
    "# This assumes knowing the minimum value, and Chi2s above Chi2_min+1\n",
    "if (ScanChi2) :\n",
    "    if (((chi2[0] - chi2_minval) > 1.0) and ((chi2[Ntau_steps] - chi2_minval) > 1.0)) :\n",
    "        found_lower = False\n",
    "        found_upper = False\n",
    "        for itau in range (Ntau_steps+1) :\n",
    "            if ((not found_lower) and ((chi2[itau] - chi2_minval) < 1.0)) :\n",
    "                tau_lower = tau[itau]\n",
    "                found_lower = True\n",
    "                \n",
    "            if ((found_lower) and (not found_upper) and ((chi2[itau] - chi2_minval) > 1.0)) :\n",
    "                tau_upper = tau[itau]\n",
    "                found_upper = True\n",
    "      \n",
    "    \n",
    "        print(f\"  Chi2 scan gives:   tau = {chi2_minpos:6.4f} + {tau_upper-chi2_minpos:6.4f} - {chi2_minpos-tau_lower:6.4f}\")\n",
    "    else :\n",
    "        print(f\"  Error: Chi2 values do not fulfill requirements for finding minimum and errors!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Discussion:\n",
    "One could here of course have chosen a finer binning, but that is still not very satisfactory, and in any case very slow. That is why we of course want to use e.g. iMinuit to perform the fit, and extract all the relevant fitting parameters in a nice, fast, numerically stable, etc. way."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "# Fit the data using iminuit (both chi2 and binned likelihood fits)\n",
    "\n",
    "Now we want to see, what a \"real\" fit gives, in order to compare our result with the one provided by Minuit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the function to fit with:\n",
    "def func_exp(x, N0, tau) :\n",
    "    return N0 * binwidth / tau * np.exp(-x/tau)\n",
    "\n",
    "# Define the function to fit with:\n",
    "def func_exp2(x, N0, tau) :\n",
    "    return N0 / tau * np.exp(-x/tau)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $\\chi^2$ fit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prepare figure\n",
    "fig_fit, ax_fit = plt.subplots(figsize=(8, 6))\n",
    "ax_fit.set_title(\"tau values directly fitted with iminuit\")\n",
    "ax_fit.set_xlabel(\"Lifetimes [s]\")\n",
    "ax_fit.set_ylabel(\"Frequency [ev/0.1s]\")\n",
    "\n",
    "# Plot our tau values\n",
    "indexes = yExp>0 # only bins with values!\n",
    "xExp = (xExp_edges[1:] + xExp_edges[:-1])/2    # Move from bins edges to bin centers\n",
    "syExp = np.sqrt(yExp)                          # Uncertainties\n",
    "ax_fit.errorbar(xExp[indexes], yExp[indexes], syExp[indexes], fmt='k_', ecolor='k', elinewidth=1, capsize=2, capthick=1)\n",
    "\n",
    "# Chisquare-fit tau values with our function:\n",
    "chi2_object_fit = Chi2Regression(func_exp, xExp[indexes], yExp[indexes], syExp[indexes])\n",
    "# NOTE: The constant for normalization is NOT left free in order to have only ONE parameter!\n",
    "\n",
    "minuit_fit_chi2 = Minuit(chi2_object_fit, N0=Ntimes, tau=tau_truth)\n",
    "minuit_fit_chi2.fixed[\"N0\"] = True\n",
    "minuit_fit_chi2.errordef = 1.0\n",
    "minuit_fit_chi2.migrad()\n",
    "\n",
    "# Plot fit\n",
    "x_fit = np.linspace(0, 10, 1000)\n",
    "y_fit_simple = func_exp(x_fit, *minuit_fit_chi2.values)\n",
    "ax_fit.plot(x_fit, y_fit_simple, 'b-', label=\"ChiSquare fit\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print the obtained fit results:\n",
    "# print(minuit_fit_chi2.values[\"tau\"], minuit_fit_chi2.errors[\"tau\"])\n",
    "tau_fit = minuit_fit_chi2.values[\"tau\"]\n",
    "etau_fit = minuit_fit_chi2.errors[\"tau\"]\n",
    "\n",
    "print(f\"  Decay time of minimum found:   chi2: {tau_fit:.3f} +- {etau_fit:.3f}s\")\n",
    "print(f\"  Chi2 value at minimum:         chi2 = {minuit_fit_chi2.fval:.1f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Alternatively to the above, one can in iMinuit actually ask for the Chi2 curve to be plotted by one command:\n",
    "minuit_fit_chi2.draw_mnprofile('tau')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "### Binned likelihood fit:\n",
    "\n",
    "Below is an example of a binned likelihood fit. Try to write an unbinned likelihood fit yourself!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Binned likelihood-fit tau values with our function\n",
    "# extended=True because we have our own normalization in our fit function\n",
    "bllh_object_fit = BinnedLH(func_exp2, t, bins=Nbins, bound=(0, tmax), extended=True)\n",
    "minuit_fit_bllh = Minuit(bllh_object_fit, tau=tau_truth, N0=Ntimes)\n",
    "\n",
    "minuit_fit_bllh.errordef = 0.5   # Value for likelihood fit\n",
    "minuit_fit_bllh.migrad()\n",
    "\n",
    "# Plot fit\n",
    "N0, tau_fit_bllh = minuit_fit_bllh.values\n",
    "x_fit = np.linspace(0, 10, 1000)\n",
    "y_fit_simple = func_exp(x_fit, N0, tau_fit_bllh)\n",
    "ax_fit.plot(x_fit, y_fit_simple, 'r-', label=\"Binned Likelihood fit\")\n",
    "\n",
    "# Define the ranges:\n",
    "ax_fit.set_xlim(0, 5)\n",
    "ax_fit.set_ylim(bottom=0)     # We don't want to see values below this!\n",
    "fig_fit.legend(loc=[0.45, 0.75])\n",
    "fig_fit.tight_layout()\n",
    "fig_fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if (save_plots) :\n",
    "    fig_fit.savefig(\"ExponentialDist_Fitted.pdf\", dpi=600)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "# Questions:\n",
    "\n",
    "1) Consider the four plots (bottom right one empty) showing chi2, bllh, and ullh as a function of lifetime, tau. Do the four curves resemble each other in shape? Are they identical in shape? Do the three methods give similar results, or are they different? Do you see the relation between the curves and the fit result? This question requires that you also fit a parabola to the other two cases. Remember to consider both central value and uncertainty of tau.\n",
    "\n",
    "2) Now consider the two (chi2 and bllh) fits by iMinuit. How alike results do they obtain? Again, consider both the central values and the uncertainty.\n",
    "\n",
    "3) Try to decrease the number of exponential numbers you consider to say 10, and see how things change. Does the difference between Chi2, bllh, and ullh get bigger or not?\n",
    "\n",
    "4) Try to increase the number of exponential numbers you consider to say 10000, and see what happens to the difference between Chi2 and BLLH? Also, does the errors become more symetric? Perhaps you will need to consider a shorter range of the fit around the mimimal value, and have to also increase the number of points you calculate the chi2/bllh/ullh (or decrease the range you search!), and possibly change the ranges of your plotting.\n",
    "\n",
    "\n",
    "### Advanced Questions:\n",
    "\n",
    "5) Make (perhaps in a new program) a loop over the production of random data,\n",
    "   and try to see, if you can print (or plot) the Chi2 and BLLH results for each\n",
    "   turn. Can you spot any general trends? I.e. is the Chi2 uncertainty always\n",
    "   lower or higher than the (B/U)LLH? And are any of the estimators biased?\n",
    "\n",
    "6) Make a copy of the program and put in a different PDF (i.e. not the exponential).\n",
    "   Run it, and see if the errors are still asymetric. For the function, try either\n",
    "   e.g. a Polynomial or a Gaussian."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Learning points:\n",
    "\n",
    "Make sure that you understand how the likelihood is different from the ChiSquare,\n",
    "and how the binned likelihood is different from the unbinned. If you don't do it,\n",
    "this exercise, and much of the course and statistics in general will be a bit lost\n",
    "on you! :-)\n",
    "\n",
    "You should understand that the binned likelihood fit resembels the ChiSquare fit a bit,\n",
    "only the evaluation in each bin is different, especially if the number of events in the bin is low,\n",
    "as the PDF considered (Poisson for the LLH, Gaussian for the ChiSquare) is then different.\n",
    "At high statistics, they give the same result, but the ChiSquare fit quality can be evaluated.\n",
    "\n",
    "Also know that the unbinned likelihood uses each single event, and is thus different at its core.\n",
    "This can make a difference, if there are only few events and/or if each event has several\n",
    "attributes (e.g. uncertainties in x and y), which can't be summarized in a simple histogram with bins.\n",
    "\n",
    "Finally, uncertainties can be asymmetric (from low statistics and/or asymmetric fit), but usually\n",
    "they are not.\n",
    "\n",
    "\n",
    "## Conclusion:\n",
    "Fitting \"manually\" is damn hard, cumbersome, and not a thing that one wants to do.\n",
    "Always let a well tested program (e.g. iMinuit or SciPy Optimize) do it, and instead\n",
    "take the inspired position of checking that the fitting program actually is doing\n",
    "what it is supposed to do, and that everything comes out reasonable.\n",
    "\n",
    "The art of fitting is multiple. **Very importantly, a fit requires good input parameters**,\n",
    "as it will otherwise not converge. Also, the Chi-square fit is more robust, so it is often\n",
    "a good idea to start with this, and if the fit converges, one can use the fitting parameters\n",
    "as input values for subsequent (likelihood) fits. Finally, one needs to consider the binning\n",
    "and fitting range carefully, and make good use of the p-value from the Chi-square fit."
   ]
  }
 ],
 "metadata": {
  "executable": "/usr/bin/env python",
  "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.11"
  },
  "main_language": "python"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}