06_Scikit-HEP_fitting.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center><h1><b><span style=\"color:blue\">Fitting and model building</span></b></h1></center>\n",
    "\n",
    "#### **Quick intro to the following packages**\n",
    "- The core package `iminuit`.\n",
    "- Model building and a word on the Scikit-HEP affiliated project/package `zfit`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&nbsp;\n",
    "<center>\n",
    "    <img src=\"images/logo_iminuit.png\" style=\"width:150px;\"/>\n",
    "    <h2><b><span style=\"color:green\">Python wrapper to Minuit2 minimization and error computation package</span></b></h2>\n",
    "</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `iminuit` package provides Python bindings for the [C++ Minuit2 library](https://root.cern.ch/guides/minuit2-manual) maintained at CERN, which is effectively the only fitting engine used in HEP. The package has no external dependency apart from NumPy.\n",
    "\n",
    "`iminuit` gives the user full power of the engine internals. That can be really useful, but it often demands a certain level of expertise. It is hence no wonder that many fitting libraries around build atop it. That's in particular the case for Astronomy.\n",
    "\n",
    "Note: feel free to complement the introduction below with the several *great tutorials* available from the [GitHub repository](https://github.com/scikit-hep/iminuit). Indeed what is shown here is just an appetiser."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **1. A very simple example**\n",
    "\n",
    "Minimisation of a function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from iminuit import Minuit\n",
    "\n",
    "def fcn(x, y, z):\n",
    "    return (x - 2) ** 2 + (y - 3) ** 2 + (z - 4) ** 2\n",
    "\n",
    "fcn.errordef = Minuit.LEAST_SQUARES\n",
    "\n",
    "m = Minuit(fcn, x=0, y=0, z=0)\n",
    "\n",
    "m.migrad()  # run optimiser"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.hesse()  # run covariance estimator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **2. A more evolved example**\n",
    "\n",
    "Let's look at a little sample of track information generated with a toy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "\n",
    "<b>Acknowledgements</b>\n",
    "\n",
    "This mini-tutorial is kindly provided by Hans Dembinski (TU Dortmund), with minor modifications.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import uproot\n",
    "\n",
    "f = uproot.open(\"data/sample_tracks.root\")\n",
    "\n",
    "event = f[\"event\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "event.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Legend :-):\n",
    "* mc_trk_len: number of true tracks in event\n",
    "* mc_trk_px: x-component of true momentum of particle (variable-length array)\n",
    "* mc_trk_py: y-component of true momentum of particle (variable-length array)\n",
    "* mc_trk_pz: z-component of true momentum of particle (variable-length array)\n",
    "* trk_len: number of reconstructed tracks in event\n",
    "* trk_px: x-component of momentum of reconstructed track (variable-length array)\n",
    "* trk_py: y-component of momentum of reconstructed track (variable-length array)\n",
    "* trk_pz: z-component of momentum of reconstructed track (variable-length array)\n",
    "* trk_imc: index of matched true particle or -1 (variable-length array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the content of branches as arrays (don't use this in large trees - you will exhaust the computer memory):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trk_len = event[\"trk_len\"].array()\n",
    "mc_trk_len = event[\"mc_trk_len\"].array()\n",
    "#trk_mom = event.arrays([\"trk_px\", \"trk_py\", \"trk_pz\"], library=\"pd\")\n",
    "trk_px = event[\"trk_px\"].array(library=\"pd\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first ten entries, this is a normal numpy array\n",
    "trk_len[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "print(f\"{np.sum(trk_len == 0)} events with zero tracks\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.hist(trk_len);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first ten entries, trk_px is a special jagged array\n",
    "trk_px[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(trk_px);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Fits**\n",
    "\n",
    "* Typical analysis work flow:\n",
    "    1. Pre-select data and make compact data trees\n",
    "    2. Make histograms and profiles from tree data\n",
    "    3. Fit histograms and profiles to extract physical parameters\n",
    "* Many specialized fitting tools for individual purposes, e.g.:\n",
    "    - [scipy.optimize.curve_fit](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html#scipy.optimize.curve_fit)\n",
    "    - [RooFit](https://root.cern.ch/roofit)\n",
    "* Generic method\n",
    "    - Select mathematical model (PDF) which describes data\n",
    "    - Use maximum-likelihood method to adapt model to data\n",
    "* Specialised methods give fast results for some types of problems\n",
    "* Generic method allows one to do advanced things not implemented in specialised methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's fit the $p_x$ distribution with a normal distribution to extract the parameters $\\mu$ and $\\sigma$.\n",
    "\n",
    "- To apply a maximum-likelihood method, we need a statistical model that describes the data\n",
    "- Assumption 1: original data before histogramming is normal distributed; pdf is $\\mathcal{N}(\\mu, \\sigma)$ with parameters $\\mu$ and $\\sigma$\n",
    "- Assumption 2: count in histogram cell is Poisson distributed $P(n_i, \\lambda_i)$\n",
    "- Expected content in a histogram cell is $\\lambda_i = N \\int_{x_i}^{x_{i+1}} \\mathcal{N}(\\mu, \\sigma) \\, \\text{d}x$, where $N$ is total number of events\n",
    "- Likelihood is joint probability of data under model\n",
    "  $L = \\prod_i P(n_i, \\lambda_i)$, need to maximize this by varying model parameters $\\mu$ and $\\sigma$\n",
    "- Technical step to achieve this: Minimize score $S(\\mu,\\sigma) = -2\\ln L(\\mu, \\sigma) = -2\\sum_i \\ln P(n_i; \\lambda_i)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import boost_histogram as bh\n",
    "\n",
    "xaxis = bh.axis.Regular(50, -2, 2)\n",
    "h_px = bh.Histogram(xaxis)\n",
    "\n",
    "h_px.fill(trk_px)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get data from before\n",
    "px_axis = h_px.axes[0]\n",
    "cx = px_axis.centers\n",
    "dx = px_axis.widths\n",
    "xe = px_axis.edges\n",
    "n = h_px.view()\n",
    "\n",
    "plt.errorbar(cx, n, n**0.5, dx, fmt=\"o\", label=\"data\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# scipy has efficient and correct implementations for most statistical distributions\n",
    "from scipy.stats import norm, poisson\n",
    "\n",
    "n_total = np.sum(n)\n",
    "\n",
    "def score(mu, sigma):\n",
    "    cdf = norm(mu, sigma).cdf\n",
    "    lambdas = n_total * (cdf(xe[1:]) - cdf(xe[:-1]))\n",
    "    probs = poisson.pmf(n, lambdas)\n",
    "    return -2 * np.sum(np.log(probs + 1e-100)) # avoid taking log of zero"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.ticker import LogLocator\n",
    "\n",
    "mus = np.linspace(-1, 1, 20)\n",
    "sigmas = np.linspace(1e-10, 2, 20)\n",
    "\n",
    "g_mu, g_sigma = np.meshgrid(mus, sigmas)\n",
    "g_score = np.vectorize(score)(g_mu, g_sigma)\n",
    "\n",
    "plt.contourf(g_mu, g_sigma, g_score)\n",
    "plt.xlabel(\"$\\mu$\")\n",
    "plt.ylabel(\"$\\sigma$\")\n",
    "plt.colorbar();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = Minuit(score, mu=0, sigma=1)\n",
    "m.errordef = Minuit.LEAST_SQUARES\n",
    "m.migrad()\n",
    "\n",
    "mu, sigma = m.values\n",
    "s_mu, s_sigma = m.errors\n",
    "\n",
    "plt.errorbar(cx, n, n ** 0.5, dx, fmt=\"o\", label=\"data\");\n",
    "plt.plot(cx, norm(mu, sigma).pdf(cx) * n_total * dx, label=\"fit\")\n",
    "plt.title(f\"$\\mu = {mu:.3f} \\pm {s_mu:.3f}, \\quad \\sigma = {sigma:.3f} \\pm {s_sigma:.3f}$\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check whether the fit is good:\n",
    "- by looking at *pull distribution*\n",
    "    - $(n_i - \\lambda_i) / \\lambda_i$ for Poisson distribute data\n",
    "- by checking the $\\chi^2$ value against the degrees of freedom\n",
    "    - Simple check: $\\chi^2/n_\\text{dof}$ should be about 1\n",
    "    - Better check: chance probability $\\int_{\\chi^2_\\text{observed}}^{-\\infty} P(\\chi^2; n_\\text{dof}) \\, \\text{d}\\chi^2$ to obtain a higher value than the observed should not be too small"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cdf = norm(mu, sigma).cdf\n",
    "\n",
    "n_pred = (cdf(xe[1:]) - cdf(xe[:-1])) * n_total\n",
    "n_sigma = n_pred ** 0.5 # for Poisson-distributed data\n",
    "\n",
    "pull = (n - n_pred) / n_sigma\n",
    "\n",
    "plt.errorbar(cx, pull, np.ones_like(pull), fmt=\"o\")\n",
    "plt.axhline(0, ls=\"--\")\n",
    "\n",
    "# degrees of freedom: number of fitted bins minus number of fitted parameters\n",
    "n_dof = len(n) - 2 # need to subtract two fitted parameters\n",
    "\n",
    "chi2_obs = np.sum(pull ** 2)\n",
    "\n",
    "print(f\"chi2/ndof = {chi2_obs} / {n_dof} = {chi2_obs / n_dof}\")\n",
    "\n",
    "from scipy.stats import chi2\n",
    "\n",
    "chance_prob = 1 - chi2(n_dof).cdf(chi2_obs)\n",
    "\n",
    "print(f\"{chance_prob}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **3. Towards realistic HEP data use-cases**\n",
    "\n",
    "HEP analyses typically involve far more sophisticated fitting work. It is often about **data model building** and performing **(un)binned maximum likelihood fits** to describe experimental distributions. Also, fitting is often performed to extract decay mode yields atop the model parameters.\n",
    "\n",
    "The `iminuit` package has functionality to performed so-called *simultaneous fits*,\n",
    "where the models being fitted to two or more datasets share at least a parameter.\n",
    "But it is not meant as a package providing commonly-used goodies and a \"framework\" for (complex) model building. Imagine how cumbersome it would be to verbosely describe the model below for an amplitude analysis ...!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **A simple binned extended maximum-likelihood fit example**\n",
    "\n",
    "(Taken straight from one of the excellent `iminuit` tutorials.)\n",
    "\n",
    "The binned extended maximum-likelihood fit is strictly the binned equivalent of the corresponding unbinned fit. One sums the logarithm of Poisson probabilities for the observed counts as a function of the predicted counts in this case (times -1 to turn maximization into minimization).\n",
    "\n",
    "Instead of a density, you need to provide a cdf of the density in this case (which must be vectorized). There is no need to separately return the total integral like the unbinned case. The parameters are the same as in the unbinned extended fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from iminuit import cost\n",
    "from scipy.stats import norm, uniform\n",
    "\n",
    "xrange = -1, 1\n",
    "\n",
    "rng = np.random.default_rng(1)\n",
    "\n",
    "xdata = rng.normal(0, 0.1, size=400)\n",
    "xdata = np.append(xdata, rng.uniform(*xrange, size=1000))\n",
    "\n",
    "def model_pdf(x, z, mu, sigma):\n",
    "    return (z * norm.pdf(x, mu, sigma) + \n",
    "            (1 - z) * uniform.pdf(x, xrange[0], xrange[1] - xrange[0]))\n",
    "\n",
    "c = cost.UnbinnedNLL(xdata, model_pdf)\n",
    "\n",
    "m = Minuit(c, z=0.4, mu=0, sigma=0.2)\n",
    "m.limits[\"sigma\"] = (0, None)\n",
    "m.limits[\"z\"] = (0, 1)\n",
    "m.limits[\"mu\"] = (-1, 1)\n",
    "\n",
    "m.migrad()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n, xe = np.histogram(xdata, bins=50, range=xrange)\n",
    "cx = 0.5 * (xe[1:] + xe[:-1])\n",
    "dx = np.diff(xe)\n",
    "\n",
    "plt.errorbar(cx, n, n ** 0.5, fmt=\"ok\")\n",
    "xm = np.linspace(xe[0], xe[-1])\n",
    "plt.plot(xm, model_pdf(xm, *[p.value for p in m.init_params]) * len(xdata) * dx[0],\n",
    "         ls=\":\", label=\"init\")\n",
    "plt.plot(xm, model_pdf(xm, *m.values) * len(xdata) * dx[0], label=\"fit\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Scaling up in complexity**\n",
    "\n",
    "In spite of the example above being a *simple* example of an extended maximum-likelihood fit to binned data, the code is, honestly, very verbose.\n",
    "Clearly much of the above could be simplified for the general user. And now imagine that the (realistic example!) function to fit is not a uniform+normal distribution but rather this:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center><img src=\"images/fitting_AmAn-example.png\" width=\"55%\"/></center>\n",
    "[Example taken from a presentation by Anton Poluektov.]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ROOT framework provides `RooFit` as a model fitting library. The PyHEP ecosystem in particular has the `zfit` package, a model fitting library based on TensorFlow and optimised for simple and direct manipulation of probability density functions ..."
   ]
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}