likelihood.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parameter Estimation - Likelihood\n",
    "=================================\n",
    "\n",
    "Take a random variable $x$ described by a pdf $f(x)$: the **sample\n",
    "space** is defined to be the set of all possible values of $x$. The set\n",
    "of $n$ independent measurements of the random variable $x$, $\\{x_i\\}$ is\n",
    "called a **sample of size** $n$.\n",
    "\n",
    "In theory of probability, from the pdf $f(x)$ we can compute all sorts\n",
    "of quantities (mean, moments, etc\\...). In statistical inference we are\n",
    "concerned about the opposite problem: take the distribution of a\n",
    "quantity measured in data and infer its parent pdf $f(x)$. In the\n",
    "simplest case of data distributed following a known pdf which depends on\n",
    "a parameter $\\theta$ (i.e. $f(x,\\theta)$), the statistical inference is\n",
    "reduced to the extraction of the best estimate of the parameter from\n",
    "data.\n",
    "\n",
    "Often when talking about an estimate we use the adjectives \"accurate\"\n",
    "and \"precise\". In what follows we mean (see {numref}`fig-PrecAcc`):\n",
    "\n",
    "-   *accuracy:* how close is the estimated value to the true reference\n",
    "    value\n",
    "\n",
    "-   *precision:* how reproducible the measurements are.\n",
    "\n",
    "This means for instance that a poorly calibrated device can give you\n",
    "high precision but poor accuracy.\n",
    "\n",
    "```{figure} ./images/ch6/PrecAcc.png\n",
    "---\n",
    "width: 400px\n",
    "align: center\n",
    "name: fig-PrecAcc\n",
    "---\n",
    "FIXME: Meaning of \"accuracy\" and \"precision\".\n",
    "```\n",
    "\n",
    "Any function of the observed measurements is called a **statistic**. A\n",
    "statistic used to estimate some parameter of a distribution is called\n",
    "**estimator**. We will generally denote an estimator of a parameter by\n",
    "adding a circumflex ($\\;\\hat{ }\\;$) to the symbol of the parameter:\n",
    "$\\hat{\\theta}$.\n",
    "\n",
    "You can build several estimators for any parameter. As an example take\n",
    "the estimation of the average height of all students enrolled at ETH. Be\n",
    "$h_i$ the outcome of the measurement of each of the $N$ students, then\n",
    "any of the following procedures would produce an estimate:\n",
    "\n",
    "-   add all $h_i$ and divide by N\n",
    "\n",
    "-   add only the first 15, divide by 15; ignore the rest\n",
    "\n",
    "-   add all $h_i$ and divide by N-1\n",
    "\n",
    "-   just quote it to be 1.82 m\n",
    "\n",
    "-   multiply all $h_i$ and take $N^{th}$-root of result\n",
    "\n",
    "-   choose the most popular height (mode)\n",
    "\n",
    "-   take shortest and tallest and divide by 2\n",
    "\n",
    "-   add 2$^{nd}$, 4$^{th}$, 6$^{th}$,\\... and divide by N/2 \\[ or\n",
    "    (N-1)/2 if N odd\\]\n",
    "\n",
    "-   take only $h_i$ of students with brown hair, divide by M\n",
    "\n",
    "All these are by definition estimators. Some appear to be clearly better\n",
    "than others, but how do we define what is a better/worse estimator? To\n",
    "answer this question we define some general properties of the\n",
    "estimators: bias, consistency, efficiency and robustness."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Properties of the estimators\n",
    "----------------------------\n",
    "\n",
    "The estimator $\\hat{\\theta}$ being a function of random variables (data)\n",
    "is itself a random variable and it will be distributed according to a\n",
    "pdf $g(\\hat{\\theta}|\\theta)$, which will clearly depend on the parameter\n",
    "$\\theta$. We define the following properties for an estimator (see {numref}`fig-estProp` )\n",
    "\n",
    "```{figure} ./images/ch6/estProp.png\n",
    "---\n",
    "width: 400px\n",
    "align: center\n",
    "name: fig-estProp\n",
    "---\n",
    "FIXME: Some estimator properties.\n",
    "```\n",
    "\n",
    "-   An estimator is called *unbiased* if its expectation value is equal\n",
    "    to the true value: $<\\hat{\\theta}>=\\theta$. Thus an estimator is\n",
    "    biased if $b_n = <\\hat{\\theta}> - \\theta \\ne 0$. The number $b_n$ is\n",
    "    called the bias of the estimator. We include the subscript $n$ in\n",
    "    this definition since we will see that some estimators are unbiased\n",
    "    only asymptotically, i.e. only for $n\\to \\infty$. An example of an\n",
    "    unbiased estimator is the mean ($\\langle \\bar{\\mu} \\rangle = \\mu$);\n",
    "    the third in the list of the previous section is asymptotically\n",
    "    unbiased ($\\langle \\hat{\\mu} \\rangle = n/(n-1)\\mu$) and so\n",
    "    $b_n(\\hat{\\mu}) \\to 0$ for $n\\to \\infty$. If we know the bias, we\n",
    "    can construct an unbiased estimator by correcting it.\n",
    "\n",
    "-   An estimator is called *consistent* if $\\forall \\epsilon > 0$,\n",
    "    $\\lim_{n \\to \\infty} P(|\\hat{\\theta}-\\theta| \\ge \\epsilon) = 0$. For\n",
    "    instance if $\\hat{\\theta}$ is the average of data distributed\n",
    "    according to a p.d.f. where we can apply the CLT, then\n",
    "    $\\hat{\\theta}$ is a consistent estimator because\n",
    "    $N(\\bar{x}; \\mu, \\sigma^2/n)$ tends to a delta function for\n",
    "    $n \\to \\infty$. In the list of the previous section for example, the\n",
    "    first and the third are consistent the second is not.\n",
    "\n",
    "-   An estimator is called *efficient* if it has the smallest possible\n",
    "    variance of $\\hat{\\theta}$ (see later in this section the minimum\n",
    "    variance bound). The efficiency $\\epsilon$ is defined as\n",
    "    $\\epsilon=\\frac{{\\rm minimal\\, Variance\\, of\\,} \\hat{\\theta}}{{\\rm Variance\\, of\\,} \\hat{\\theta}}$.\n",
    "\n",
    "-   An estimator is called *robust* if it is insensitive to wrong data\n",
    "    or wrong assumptions, especially in the tails of a distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Estimation of the Mean\n",
    "----------------------\n",
    "\n",
    "The estimator for the mean $\\mu$ obtained from $n$ independent\n",
    "measurements $x_{i}$ is: \n",
    "\n",
    "$$\\label{estimate_mean}\n",
    "\\hat{\\mu}=\\frac{1}{n}\\sum_i x_i.$$\n",
    "\n",
    "This estimator is unbiased, i.e.\n",
    "$\\langle\\hat{\\mu}\\rangle=\\langle \\frac{1}{n}\\sum_i x_i\\rangle=\\frac{1}{n}\\sum_i\\langle x_i\\rangle = \\mu$.\n",
    "\n",
    "Furthermore it is consistent because of the CLT. Its variance is given\n",
    "by \n",
    "\n",
    "$$V(\\hat{\\mu})=\\frac{1}{n}\\sigma^2.$$\n",
    "\n",
    "Whether this estimator is\n",
    "efficient or not depends on the p.d.f. of the parent distribution. For\n",
    "instance, given a uniform distribution the mean is not the most\n",
    "efficient estimator; the estimator $\\hat{\\mu}=0.5(x_{max}+x_{min})$ has\n",
    "a smaller variance. The robustness for the sample mean is increased if\n",
    "the truncated mean is used. This means that the largest and smallest\n",
    "values are trimmed (truncated). This more robust mean is less sensitive\n",
    "to outliers, but unless the parent distribution is symmetric it will be\n",
    "biased. An example for a truncated mean can be found in sports rating\n",
    "when only 4 out of 6 grades are used to form the final grade."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Estimation of the Variance\n",
    "--------------------------\n",
    "\n",
    "An estimator for the variance of a parent distribution $\\sigma^2$, when\n",
    "we know the true mean $\\langle x \\rangle = \\mu$ is:\n",
    "\n",
    "$$\n",
    "s_1^2 = \\frac{1}{n}\\sum_{i=1}^n (x_i - \\mu)^2.\n",
    "$$ \n",
    "\n",
    "This estimator is unbiased \n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "<s_1^2> &=& \\frac{1}{n}\\langle\\sum(x_i - \\mu)^2\\rangle  \\\\\n",
    "        &=& \\frac{1}{n}\\left(\\langle \\sum x_i^2 \\rangle - 2\\mu\\langle \\sum x_i \\rangle + n \\mu^2\\right) \\\\\n",
    "        &=& \\frac{1}{n}\\left( n \\langle  x^2 \\rangle - 2n\\mu\\langle x \\rangle + n \\mu^2\\right)  \\qquad (independent~x_i~so~: \\langle\\sum x_i^2\\rangle= n\\langle x^2 \\rangle)\\\\\n",
    "        &=& \\langle x^2 \\rangle - 2\\mu^2 + \\mu^2 \\\\\n",
    "        &=& \\sigma^2 -\\mu^2 + \\mu^2  \\qquad (\\sigma^2 = \\langle x^2 \\rangle - \\mu^2) \\\\\n",
    "        &=& \\sigma^2        \\\\\n",
    "\\end{aligned}\n",
    "$$ \n",
    "        \n",
    "So $s_1^2$ is an unbiased\n",
    "estimator of the variance of the parent p.d.f $\\sigma^2$, when $\\mu$ is\n",
    "known.\n",
    "\n",
    "When it is not known, we use the estimate $\\bar{x} = \\hat{\\mu}$ and\n",
    "define:\n",
    "\n",
    "$$\n",
    "s_x^2 = \\frac{1}{n}\\sum(x_i-\\bar{x})^2 = \\bar{x^2} - \\bar{x}^2\n",
    "$$\n",
    "\n",
    "The expectation value of $s_x^2$ is: \n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\langle s_x^2 \\rangle= \\frac{1}{n}\\left( \\langle\\sum x_i^2\\rangle- \\frac{1}{n}\\langle\\left( \\sum x_i  \\right)^2  \\rangle\\right).\\end{aligned}\n",
    "$$\n",
    "\n",
    "Substituting: \n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    " \\langle\\sum x_i^2\\rangle&=& n\\langle x^2 \\rangle\\\\\n",
    " \\sigma^2 &=& \\langle x^2\\rangle- \\mu^2\\\\\n",
    " V\\left( \\sum x_i \\right) &=& \\langle\\left( \\sum x_i\\right)^2\\rangle- \\left( \\langle\\sum x_i\\rangle\\right)^2\n",
    " \\end{aligned}\n",
    "$$\n",
    "\n",
    "we get\n",
    "\n",
    "$$\n",
    "\\langle s_x^2\\rangle= \\frac{1}{n}\\left( n(\\sigma^2+\\mu^2) - \\frac{1}{n}\\left( \\langle\\left( \\sum x_i\\right)^2 \\rangle+ \\langle\\sum x_i\\rangle^2 \\right)  \\right)\n",
    "$$\n",
    "\n",
    "and using \n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    " V(\\langle\\sum x_i\\rangle) &=& \\sum V(x_i) = nV(x) = n\\sigma^2 \\\\\n",
    "\\langle\\sum x_i\\rangle&=& n \\langle x \\rangle= n \\mu\n",
    "\\end{aligned}\n",
    "$$ \n",
    "\n",
    "we finally have\n",
    "\n",
    "$$\n",
    "\\langle s_x^2 \\rangle= \\frac{1}{n}\\left( n\\sigma^2 + n \\mu^2 - \\frac{1}{n}(n\\sigma^2 + (n\\mu)^2 ) \\right) = \\frac{1}{n}(n-1) \\sigma^2 .\n",
    "$$\n",
    "\n",
    "This means that $s_x^2$ is a biased estimator of $\\sigma^2$. The reason\n",
    "is that we used the sample mean $\\bar{x}$ as an estimator of the true\n",
    "mean $\\mu$. The spread of the data around the sample mean is less than\n",
    "the spread around the true mean and since the variance is the spread\n",
    "around the true mean $s_x^2$ underestimates the true variance.\\\n",
    "The formula for the variance we are used to is the one with the\n",
    "corrected bias:\n",
    "\n",
    "$$\n",
    "s^2 = \\frac{n}{n-1} s_x^2 = \\frac{n}{n-1}(<x^2> - <x>^2) = \\frac{1}{n-1}\\sum (x_i -<x>)^2.\n",
    "$$\n",
    "\n",
    "In the same way as we got to the unbiased estimator of the variance, we\n",
    "obtain the expression for the unbiased estimator of the covariance\n",
    "$V_{xy}$ of two random variables $x$ and $y$ with unknown (but\n",
    "estimated) means\n",
    "\n",
    "$$\n",
    "\\hat{V}_{xy}=\\frac{1}{n-1}\\sum(x_i-<x>)(y_i-<y>)=\\frac{n}{n-1}(<xy>-<x><y>).\n",
    "$$\n",
    "\n",
    "The correlation coefficient is then given by\n",
    "\n",
    "$$\n",
    "\\rho_{xy}=\\frac{\\hat{V}_{xy}}{s_xs_y}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maximum Likelihood Method\n",
    "-------------------------------------------------------------------------------------\n",
    "\n",
    "Assume we have $n$ measurements of a random variable $x$, distributed\n",
    "according to a known probability density function $f(x|\\theta)$. Where\n",
    "$\\theta$ stands for one parameter of the p.d.f. (e.g. the mean). Note\n",
    "that here you assume you know what is the correct pdf to fit and you are\n",
    "\\\"only\\\" interested in finding the value of the parameter $\\theta$ that\n",
    "allows the model to best fit the data.\n",
    "\n",
    "Can we find a *general* method to build an estimator for $\\theta$ ( i.e.\n",
    "$\\hat{\\theta}$ ) ? The procedure we present here goes under the name of\n",
    "**maximum likelihood method** and it is the most intuitive way to set up\n",
    "such an estimator.\n",
    "\n",
    "To understand the maximum likelihood method for parameter estimation\n",
    "(sometimes abbreviated as ML method) we start from the probability\n",
    "$f(x|\\theta)dx$ (i.e. the probability to observe $x\\in (x,x+dx)$ given\n",
    "$\\theta$). With this we can compute the probability to observe a certain\n",
    "set of data $\\{x_i\\}$ *given* the parameter $\\theta$, as the joint\n",
    "probability:\n",
    "\n",
    "$$\n",
    "P= f(x_1|\\theta)dx_1\\cdot f(x_2|\\theta)dx_2\\cdot \\ldots \\cdot f(x_n|\\theta)dx_n\n",
    "$$\n",
    "\n",
    "```{maergin}\n",
    "Note that likelihood and probability are both translated in German\n",
    "as \"Wahrscheinlichkeit\". Nevertheless there is a fundamental\n",
    "difference between a true analytical probability density function\n",
    "and a likelihood function. The latter is a function of a sample and\n",
    "therefore also a random variable.\n",
    "```\n",
    "The **likelihood function** is then defined as:\n",
    "\n",
    "$$\n",
    "L(\\theta)=f(x_1|\\theta)\\cdot f(x_2|\\theta)\\cdots f(x_n|\\theta)=\\prod_{i = 1}^{Nevts} f(x_i|\\theta).\n",
    "$$\n",
    "\n",
    "where the product runs over all the events in the data sample $\\{x_i\\}$\n",
    "and $L(\\theta)$ is normalized to 1 for all values of $\\theta$.\n",
    "\n",
    "The function $L(\\theta)$ is, for a given data sample, a function of only\n",
    "the parameter $\\theta$ and it gives us the probability to get with this\n",
    "choice of the parameter $\\theta$ the measured values $\\{x_{i}\\}$. The\n",
    "likelihood function is *not* a probability density function in the\n",
    "parameters $\\theta$ (if that was the case we would calculate explicitly\n",
    "the expectation value of $\\theta$ and all its higher moments).\n",
    "\n",
    "The **ML principle** states that the best estimate of $\\theta$ is given\n",
    "by the estimator $\\hat{\\theta}$ which maximizes $L(\\theta)$, i.e. the\n",
    "value which maximizes the probability to obtain the observed set of\n",
    "observed data $\\{x_{i}\\}$ given $\\theta$:\n",
    "\n",
    "$$\n",
    "\\hat{\\theta} = max_\\theta~L(\\theta)\n",
    "$$\n",
    "\n",
    "The maximum is then computed by differentiating $dL(\\theta) / d\\theta = 0$.\n",
    "\n",
    "The concept is trivially generalized to several parameters $a_{k}$\n",
    "requiring: $\\partial L/\\partial a_k=0 \\ \\ \\forall k$.\n",
    "\n",
    "In practice we often work with the (natural) logarithm of the likelihood\n",
    "function $l(\\theta) = \\ln L(\\theta)$ (slang: the *log-likelihood*). The\n",
    "reason for this is that likelihoods are usually calculated with\n",
    "computers and the product of several probabilities (i.e. numbers smaller\n",
    "than 1) will hit the numerical precision of the machine. The logarithm\n",
    "transforms the product into a sum which does not create precision\n",
    "problems. Since the logarithm is a monotonic rising function, the value\n",
    "that maximizes $L(\\theta)$ also maximizes $\\ln L(\\theta)$, and our\n",
    "condition becomes:\n",
    "\n",
    "$$\n",
    "l(\\theta)=\\ln L(\\theta)=\\sum_{i=1}^{Nevts} \\ln f(x_i|\\theta) = \\mbox{Maximum}.\n",
    "$$\n",
    "\n",
    "Any monotonic transformation of the likelihood will leave its maximum\n",
    "unchanged.\n",
    "\n",
    "It has to be stressed again that the ML estimation yields a value\n",
    "$\\hat{\\theta}$, for which the observed data are the most \"likely\\\"\n",
    "(compared to other parameter values), and not vice-versa! It is not to\n",
    "be confused with the statement that the parameter $\\hat{\\theta}$ is the\n",
    "most probable value.\n",
    "\n",
    "```{admonition} Example:\n",
    ":class: tip\n",
    "Take the probability density given by\n",
    "$f(x|\\theta) = 1 + \\theta(x-0.5)$ with $x$ between 0 and 1. The provided\n",
    "sample data $\\{x_{i}\\}= \\{0.89, 0.03, 0.5, 0.36, 0.49\\}$. The\n",
    "log-likelihood function is then given by \n",
    "\n",
    "$$\n",
    "l(\\theta)=\\sum_{i=1}^5 \\ln(1+\\theta(x_i-0.5))\n",
    "$$\n",
    "\n",
    "The maximum of the log-likelihood function can be determined graphically\n",
    "to be $\\hat{\\theta} = -0.6$.\n",
    "```\n",
    "\n",
    "```{figure} ./images/ch6/LogLikelihoodJPG.png\n",
    "---\n",
    "width: 400px\n",
    "align: center\n",
    "name: fig-LogLikelihoodJPG\n",
    "---\n",
    "FIXME: The log-likelihood function from the example $l(\\theta)=\\sum_{i=1}^5 \\ln(1+\\theta(x_i-0.5))$.\n",
    "```\n",
    "\n",
    "\n",
    "The numerical libraries in use to solve optimization problems search for\n",
    "minima, not maxima. To use a minimization program to find a maximum you\n",
    "just need to flip the sign of the likelihood: the problem is trivially\n",
    "moved from a maximization to a minimization. In the following we will\n",
    "typically work with the $-\\ln L(\\theta)$ (slang: negative log-likelihood\n",
    "or NLL). \n",
    "\n",
    "```{admonition} Example:\n",
    ":class: tip\n",
    "*Exponential decay.* Let's derive the ML estimator for a\n",
    "particle's lifetime. The proper decay time of an unstable particle with\n",
    "lifetime $\\tau$ follows the exponential distribution:\n",
    "\n",
    "$$\n",
    "f(t;\\tau)=\\frac{1}{\\tau}e^{-t/\\tau}\n",
    "$$ \n",
    "\n",
    "Given a set of $n$ measurements\n",
    "$\\{t_i\\}$ of the proper-decay time, we can write the log-likelihood as:\n",
    "\n",
    "$$\n",
    "l(\\tau)=\\ln L(\\tau)= \\ln \\prod_i f(t_i;\\tau) = \\sum_i \\ln f(t_i;\\tau)=\\sum_i\\left(\\ln\\frac{1}{\\tau}-\\frac{t_i}{\\tau}\\right).\n",
    "$$\n",
    "\n",
    "Maximizing the log-likelihood with respect to $\\tau$:\n",
    "\n",
    "$$\n",
    "\\frac{\\partial l}{\\partial \\tau}= \\sum_i \\left(-\\frac{1}{\\tau} + \\frac{t_i}{\\tau^2} \\right) = -\\frac{n}{\\tau}+\\frac{1}{\\tau^2}\\sum t_i = 0\n",
    "$$\n",
    "\n",
    "we obtain the ML estimator $\\hat{\\tau}$:\n",
    "\n",
    "$$\n",
    "\\hat{\\tau}=\\frac{1}{n}\\sum_i t_i\n",
    "$$\n",
    "Hence we get the mean as the ML\n",
    "estimator of the lifetime! Furthermore it shows that the ML estimator is\n",
    "asymptotically unbiased: increasing the number of measurements the bias\n",
    "decreases.\n",
    "```\n",
    "\n",
    "```{admonition} Example:\n",
    ":class: tip\n",
    "*Gaussian distribution.* The Gaussian probability\n",
    "distribution function is given by\n",
    "\n",
    "$$\n",
    "f(x_i;\\mu)=\\frac{1}{\\sqrt{2\\pi}}\\frac{1}{\\sigma_i}\\cdot e^{-\\frac{1}{2}\\left(\\frac{x_i-\\mu}{\\sigma_i}\\right)^2}.\n",
    "$$\n",
    "\n",
    "To get a ML estimator for the mean $\\hat{\\mu}$ we construct again the\n",
    "log-likelihood function: \n",
    "\n",
    "$$\n",
    "\\label{eq_gauss_llh}\n",
    "l(\\mu)=\\ln L(\\mu)= \\ln \\prod_i f(t_i;\\mu) = \\sum_i \\ln f(t_i;\\mu) = \n",
    "\\sum_i\\left( \\ln\\frac{1}{\\sqrt{2\\pi}} - \\ln \\sigma_i - \\frac{1}{2} \\left( \\frac{x_i-\\mu}{\\sigma_i} \\right)^2 \\right)\n",
    "$$\n",
    "\n",
    "Differentiating with respect to $\\mu$, the determination of the maximum\n",
    "yields: \n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\frac{dl(\\mu)}{d\\mu}&=& \\frac{d}{d\\mu} \\sum_i -\\frac{1}{2}\\left( \\frac{x_i-\\mu}{\\sigma_i}\\right)^2   = \\sum_i\\frac{x_i-\\mu}{\\sigma_i^2}=0 \\\\\n",
    "\\hat{\\mu}&=&\\frac{\\sum_i x_i/\\sigma_i^2 }{\\sum_i 1/\\sigma_i^2}\\end{aligned}\n",
    "$$\n",
    "\n",
    "Which is the weighted mean of the sample $\\{x_i\\}$ and it simplifies to\n",
    "$\\hat{\\mu} = \\frac{1}{n} \\sum_{i}x_{i}$ if all the $x_{i}$ have the same\n",
    "$\\sigma_{i}$. In this case ($\\sigma_{i} = \\sigma \\, \\forall i$) we can\n",
    "use the likelihood method to get an estimate for the variance\n",
    "$\\hat{\\sigma}^{2}$. The ML method yields\n",
    "\n",
    "$$\n",
    "\\hat{\\sigma}^{2}=\\frac{1}{n}\\sum_i(x_i-\\hat{\\mu})^2.\n",
    "$$\n",
    "which is, as\n",
    "already discussed, asymptotically unbiased.\n",
    "```\n",
    "\n",
    "```{admonition} Example:\n",
    ":class: tip\n",
    "*Poisson distribution:* Consider a set of data $\\{r_i\\}$\n",
    "which we assume to be distributed according to a Poisson distribution\n",
    "with parameter $\\lambda$. We use the ML method to find an estimator for\n",
    "$\\lambda$. The log-likelihood function for the Poisson distribution is\n",
    "given by\n",
    "\n",
    "$$\n",
    "l(\\lambda) = \\sum_i\\ln \\frac{\\lambda^{r_i}}{r_i!}e^{-\\lambda} = \\sum_i \\ln \\lambda^{r_i} - n\\lambda - \\sum_i \\ln r_i! =\n",
    "\\ln\\lambda\\cdot\\sum_i r_i-n\\lambda - \\sum_i \\ln r_i!\n",
    "$$ \n",
    "\n",
    "Differentiating\n",
    "$l(\\lambda)$ w.r.t. $\\lambda$ and equating it to zero gives as estimator\n",
    "for the mean of a Poisson distribution\n",
    "$\\hat{\\lambda} = \\frac{1}{n} \\sum_{i} r_{i}$, which is again the mean of\n",
    "the sample.\n",
    "```\n",
    "\n",
    "```{admonition} Example:\n",
    ":class: tip\n",
    "*Binomial distribution:* As for the Poisson case above,\n",
    "consider a set of data $\\{s_i, r_i\\}$ (for the measurement $i$ we obtain\n",
    "$s_i$ successes and $r_i$ failures, $n_i = s_i +r_i$) which we assume to\n",
    "be distributed according to a binomial distribution\n",
    "$B(p) = {n \\choose s} p^{s} (1-p)^{r}$ with $s + r = n$. We use the ML\n",
    "method to find an estimator for $p$. The log-likelihood function for the\n",
    "binomial distribution is given by:\n",
    "\n",
    "$$\n",
    "l(p)=\\ln B(p) = \\ln{n \\choose s} + s\\ln p + r\\ln(1-p)\n",
    "$$ \n",
    "\n",
    "The\n",
    "requirement $\\frac{\\partial l(p)}{\\partial p} = 0$ yields\n",
    "$\\frac{s}{p} - \\frac{r}{1-p}=0$ and hence $\\hat{p}=s/n$, which is the\n",
    "fraction of successes given n trials.\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Minimum Variance Bound\n",
    "----------------------\n",
    "\n",
    "### Information\n",
    "\n",
    "We introduce here the concept of information following Fisher's\n",
    "definition. Any information definition should fulfil the following\n",
    "criteria:\n",
    "\n",
    "-   the information should increase if we make more observations - add\n",
    "    more data\n",
    "\n",
    "-   data, which are irrelevant to the estimation of the parameters we\n",
    "    wish to estimate or to the hypothesis we wish to test, should\n",
    "    contain no information\n",
    "\n",
    "-   the precision of the estimation should be greater if we have more\n",
    "    information\n",
    "\n",
    "The **Fisher information** (information for short in the following) on a\n",
    "parameter $\\theta$ given by a data set $\\{\\vec{x}\\}$ of the random\n",
    "variable $x$ is defined as the expectation value:\n",
    "\n",
    "$$\n",
    "I_{\\vec{x}}(\\theta) = \\langle\\left( \\frac{\\partial \\ln L(\\vec{x};\\theta)}{\\partial \\theta}  \\right)^2  \\rangle=\n",
    "$$\n",
    "\n",
    "$$= \\langle\\left( \\frac{\\partial l}{ \\partial \\theta }  \\right)^2  \\rangle= \\int \\left( \\frac{\\partial \\ln L(\\vec{x};\\theta)}{\\partial \\theta} \\right)^2 L(\\vec{x}; \\theta) d\\vec{x}\n",
    "$$\n",
    "\n",
    "To have a more compact notation we define the **score** of one\n",
    "measurement as the random variable:\n",
    "\n",
    "$$\n",
    "S_1 = \\frac{\\partial}{\\partial \\theta}\\ln f(x;\\theta)\n",
    "$$\n",
    "\n",
    "The score of a sample is the sum of the score of each measurement:\n",
    "\n",
    "$$\n",
    "S(\\vec{x},\\theta) = \\sum_{i=1}^{n} S_1(x_i;\\theta)\n",
    "$$\n",
    "\n",
    "and it is equal\n",
    "to the derivative of the log-likelihood w.r.t to the parameter of\n",
    "interest:\n",
    "\n",
    "$$\n",
    "S(\\vec{x},\\theta) = \\frac{\\partial \\ln L(\\vec{x},\\theta)}{\\partial \\theta}\n",
    "$$\n",
    "\n",
    "So the definition of the information of the sample $\\vec{x}$ on the\n",
    "parameter $\\theta$, can be rewritten as the expectation value of the\n",
    "square of the score:\n",
    "\n",
    "$$\n",
    "I_{\\vec{x}}(\\theta) = \\langle S^2(\\vec{x};\\theta) \\rangle\n",
    "$$ \n",
    "\n",
    "If $\\ln L(\\vec{x},\\theta)$ is twice differentiable w.r.t. $\\theta$, then\n",
    "the Fisher information can be rewritten as:\n",
    "\n",
    "$$\n",
    "I_{\\vec{x}}(\\theta)  = \\langle\\; \\left( \\frac{\\partial}{\\partial \\theta} \\ln L \\right)^2 \\;\\rangle\n",
    "= -\\langle\\; \\frac{\\partial^2}{\\partial \\theta^2} \\ln L \\;\\rangle\n",
    "$$\n",
    "\n",
    "*Proof*: \n",
    "\n",
    "$$\\frac{\\partial^2}{\\partial\\theta^2}\\ln L =\n",
    "\\frac{\\partial}{\\partial\\theta}\\left( \\frac{1}{L} \\frac{\\partial L}{\\partial \\theta} \\right) =\n",
    "\\frac{-(\\partial L / \\partial \\theta )^2}{L^2} + \\frac{1}{L}\\frac{\\partial^2 L}{\\partial\\theta^2}\n",
    "$$\n",
    "\n",
    "$$\n",
    "-\\langle\\frac{\\partial^2 \\ln L}{\\partial\\theta^2}\\rangle=\n",
    "\\langle\\left(\\frac{\\partial L / \\partial \\theta }{L} \\right)^2\\rangle- \\langle\\frac{ \\partial^2  L / \\partial\\theta^2}{L} \\rangle=\n",
    "\\langle\\left(\\frac{\\partial \\ln L}{\\partial \\theta}\\right)^2 \\rangle-\n",
    "\\int \\left( \\frac{\\partial^2 L }{\\partial \\theta^2} \\right) \\frac{1}{L} L dx =\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\langle\\left(\\frac{\\partial \\ln L}{\\partial \\theta}\\right)^2 \\rangle= I_{\\vec{x}} (\\theta)\n",
    "$$\n",
    "\n",
    "the last is because\n",
    "$\\int \\left( \\frac{\\partial^2 L }{\\partial \\theta^2} \\right) \\frac{1}{L} L dx = \\frac{\\partial^2}{\\partial \\theta^2} \\int L dx = \\frac{\\partial^2}{\\partial \\theta^2} 1 = 0$.\n",
    "\n",
    "With these definitions we can check that the Fisher information fulfils\n",
    "the requirements shown above.\n",
    "\n",
    "$\\rightarrow 1)$ The information should increase if we make more\n",
    "observations. For n measurements:\n",
    "\n",
    "$$\n",
    "I(\\theta) = \\langle\\left( \\sum_{i=1}^n S_1 (x_i ; \\theta \\right)^2 \\rangle\n",
    "$$\n",
    "\n",
    "$$\n",
    "I(\\theta) = \\langle\\left( \\sum_{i=1}^n S_1(x_i; \\theta)  \\right)^2\\rangle\n",
    "= V\\left( \\sum_{i=1}^n S_1(x_i; \\theta) \\right) + \\langle\\sum_{i=1}^n S_1(x_i; \\theta)  \\rangle^2\n",
    "$$\n",
    "\n",
    "where we used $V(a) = \\langle a^2\\rangle- \\langle a \\rangle^2$. Assuming\n",
    "that the single measurements $x_i$ are independent, the variance of the\n",
    "sum is the sum of the variances. And since all the measurements are\n",
    "taken from the same p.d.f., the variance is the same for all $i$. A\n",
    "similar argument applies to the second term. So:\n",
    "\n",
    "$$\n",
    "I(\\theta) = n V\\left( S_1(x; \\theta) \\right) + n^2 \\langle S_1(x; \\theta)  \\rangle^2\n",
    "$$\n",
    "\n",
    "which shows that the information increases with the number of\n",
    "observations.\n",
    "\n",
    "$\\rightarrow 2)$ Irrelevant data carry no information. For irrelevant\n",
    "data the p.d.f. will not depend on $\\theta$; the score will be 0 adding\n",
    "no information.\n",
    "\n",
    "$\\rightarrow 3)$ The precision should be greater if we have more\n",
    "information. This comes from the definition of Fisher's information:\n",
    "\n",
    "$$\n",
    "\\sigma^2(\\hat{\\theta}) = \\frac{1}{d^2 L/d\\theta^2|_{\\theta = \\hat{\\theta}}} = \\frac{1}{I(\\hat{\\theta})}\n",
    "$$\n",
    "\n",
    "The variance is the inverse of the second derivative of the likelihood,\n",
    "i.e. the inverse of the information. The larger the information the\n",
    "smaller the variance. Another way to look at it: think about the second\n",
    "derivative computed at the best estimate of the parameter\n",
    "($\\hat{\\theta}$) as the curvature of the likelihood at that point. The\n",
    "larger the curvature, the more pronounced the minimum, the larger the\n",
    "information in the data set (the smaller the uncertainty on the\n",
    "parameter see [ML uncertainty](likelihood.html#uncertainty-for-ml-estimators-sec-mlunc)). \n",
    "Now go to the other extreme: a likelihood that\n",
    "does not depend on a parameter will be flat with respect to it, so the\n",
    "curvature and the information will be zero, and the variance infinite."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rao-Cramér-Frechet inequality\n",
    "\n",
    "The \\\"Rao-Cramér-Frechet inequality\\\" (1945) tells that any estimator\n",
    "will never have a variance smaller than a given number that depends on\n",
    "the information contained in the dataset and the bias of the estimator.\n",
    "\n",
    "To see how this is possible, let's compute the covariance among the\n",
    "score and the MLE (both are random variables):\n",
    "\n",
    "$$\n",
    "\\mbox{cov}[S(\\vec{x}; \\hat{\\theta}(\\vec{x})), \\hat{\\theta}(\\vec{x})] = \\langle S(\\vec{x}; \\hat{\\theta}(\\vec{x}))\\hat{\\theta}(\\vec{x})\\rangle-\\langle S(\\vec{x}; \\hat{\\theta}(\\vec{x}))\\rangle\\langle\\hat{\\theta}(\\vec{x})\\rangle\n",
    "$$\n",
    "\n",
    "Let's compute it taking $\\hat{\\theta}$ as an estimator of $\\theta$ with\n",
    "bias $b_n(\\hat{\\theta}) = \\langle\\hat{\\theta} \\rangle- \\theta$ and\n",
    "assume that its variance is finite and that the range of $x$ does not\n",
    "depend on $\\theta$. Then we can write the first term as:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\langle\\hat{\\theta} S(\\vec{x},\\theta) \\rangle&=& \\int\\ldots\\int \\hat{\\theta}\\left( \\frac{\\partial}{\\partial \\theta} \\ln L(\\vec{x},\\theta) \\right) L(\\vec{x},\\theta) dx_1\\ldots dx_n\\\\\n",
    "&=& \\int\\ldots\\int \\hat{\\theta}\\left( \\frac{1}{L(\\vec{x},\\theta)}\\frac{\\partial}{\\partial \\theta} L(\\vec{x},\\theta) \\right) L(\\vec{x},\\theta) dx_1\\ldots dx_n\\\\\n",
    "&=& \\int\\ldots\\int \\hat{\\theta}\\left( \\frac{\\partial}{\\partial \\theta} L(\\vec{x},\\theta) \\right)dx_1\\ldots dx_n\\\\\n",
    "&=& \\int\\ldots\\int \\hat{\\theta}\\frac{\\partial}{\\partial \\theta}\\left( \\prod_{i=1}^nf(x_i;\\theta)dx_i \\right)\\\\\n",
    "&=& \\int\\ldots\\int \\frac{\\partial}{\\partial \\theta}\\left(\\hat{\\theta} \\prod_{i=1}^nf(x_i;\\theta)dx_i \\right)\\end{aligned}\n",
    "$$\n",
    "\n",
    "the last step follows because $\\hat{\\theta}$ is a statistic (a function\n",
    "of the data only) and therefore does not depend on $\\theta$. Then we\n",
    "change the order of integration and differentiation: \n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\langle\\hat{\\theta} S(\\vec{x},\\theta) \\rangle&=&  \\frac{\\partial}{\\partial \\theta}\\int\\ldots\\int \\hat{\\theta} \\prod_{i=1}^n f(x_i;\\theta)dx_i \\\\\n",
    "&=& \\frac{\\partial}{\\partial \\theta} \\langle\\hat{\\theta}\\rangle= \\frac{\\partial}{\\partial \\theta} (\\theta + b_n(\\hat{\\theta}))\\\\\n",
    "&=& 1 + \\frac{\\partial}{\\partial \\theta} b_n(\\hat{\\theta})\\end{aligned}\n",
    "$$\n",
    "\n",
    "The second term is zero because\n",
    "$\\langle S(\\vec{x};\\theta)\\rangle= \\sum \\langle S_1 (x_i;\\theta)\\rangle$\n",
    "\n",
    "and \n",
    "\n",
    "$$\\begin{aligned}\n",
    " \\langle S_1(x;\\theta)  \\rangle&=& \\langle\\frac{\\partial}{\\partial \\theta} \\ln f(x;\\theta) \\rangle\\\\\n",
    " &=& \\int \\left(  \\frac{\\partial}{\\partial \\theta} \\ln f(x;\\theta)  \\right) f(x; \\theta) dx \\\\\n",
    " &=& \\int \\frac{1}{f(x;\\theta)} \\left( \\frac{\\partial}{\\partial \\theta} f(x;\\theta) \\right) f(x;\\theta) dx \\\\\n",
    " &=&  \\int \\frac{\\partial}{\\partial \\theta} f(x;\\theta) dx\\\\\n",
    " \\end{aligned}\n",
    "$$\n",
    "\n",
    "interchanging the order of integration and\n",
    "differentiation (this usually holds for smooth distributions encountered\n",
    "in physics):\n",
    "\n",
    "$$\n",
    "\\langle S_1(x;\\theta) \\rangle=  \\frac{\\partial}{\\partial \\theta} \\int f(x;\\theta) dx = \\frac{\\partial}{\\partial \\theta} 1 = 0\n",
    "$$\n",
    "\n",
    "since $f(x;\\theta)$ is normalized for all values of $\\theta$.\n",
    "\n",
    "Putting everything together:\n",
    "\n",
    "$$\n",
    "\\mbox{cov}[S(\\vec{x},\\hat{\\theta}),\\hat{\\theta}(\\vec{x})] = 1 + \\frac{\\partial}{\\partial \\theta} b_n(\\hat{\\theta}).\n",
    "$$\n",
    "\n",
    "Their correlation coefficient is:\n",
    "\n",
    "$$\n",
    "\\rho^2 = \\frac{(\\mbox{cov}[S,\\hat{\\theta}])^2}{V(S)V(\\hat{\\theta})} = \\frac{\\left( 1 + \\frac{\\partial}{\\partial \\theta} b_n(\\hat{\\theta})\\right)^2}{I(\\theta)V(\\hat{\\theta})}\n",
    "$$\n",
    "\n",
    "and since $\\rho^2 \\le 1$, we have\n",
    "\n",
    "$$\n",
    "V(\\hat{\\theta}) \\ge \\frac{\\left( 1 + \\frac{\\partial}{\\partial \\theta} b_n(\\hat{\\theta})\\right)^2}{I(\\theta)}\n",
    "$$\n",
    "\n",
    "This is the so called \"**Rao-Cramér-Frechet inequality**\" (RCF) or\n",
    "\"Information inequality\".\n",
    "\n",
    "This inequality means that there is a lower bound on the variance of the\n",
    "estimator; i.e. given a certain amount of information (a data set) we\n",
    "can never find an estimator with lower variance than this bound. To\n",
    "reduce the bound we need to get more information or get rid of the bias.\n",
    "For an unbiased estimator the bound becomes\n",
    "$V(\\hat{\\theta}) = 1/I(\\theta)$.\n",
    "\n",
    "Now that we know what is the minimum variance of an estimator we can\n",
    "also define the **efficiency** of the estimator as:\n",
    "\n",
    "$$\n",
    "\\epsilon(\\hat{\\theta}) = \\frac{V_{min}(\\hat{\\theta})}{V(\\hat{\\theta})} \\le 1\n",
    "$$\n",
    "\n",
    "which for an unbiased estimator is\n",
    "\n",
    "$$\n",
    "\\epsilon(\\hat{\\theta}) = \\frac{1}{V(\\hat{\\theta})I(\\theta)} \\le 1\n",
    "$$\n",
    "\n",
    "An estimator with $\\epsilon = 1$ is called *efficient*. It is not always\n",
    "possible to find an *efficient* estimator, but it can be shown that:\n",
    "\n",
    "-   if an efficient estimator for a given problem exist, it will be\n",
    "    found using the ML method\n",
    "\n",
    "-   ML estimators are efficient in the large sample limit.\n",
    "\n",
    "In simple words, the maximum likelihood estimator is the best you can\n",
    "get\\...\n",
    "\n",
    "*Theorem*: An efficient estimator can be found if and only if it belongs\n",
    "to the exponential family:\n",
    "\n",
    "$$\n",
    "f(x;\\theta) = \\exp[A(\\theta)\\hat{\\theta}(x) + B(\\theta) + C(x) ]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Uncertainty for ML estimators\n",
    "-----------------------------\n",
    "\n",
    "Let's take the simplest case of a likelihood with only one parameter in\n",
    "the large sample limit (i.e. the estimator is asymptotically unbiased,\n",
    "efficient and the RCF is valid as an equality). Expand its NLL function\n",
    "around $\\theta = \\hat{\\theta}$:\n",
    "\n",
    "$$\n",
    "F(\\theta) = -\\ln L(\\theta)=F(\\hat{\\theta})+\\frac{1}{2}\\frac{d^2F}{d\\theta^2}\\lvert_{\\theta=\\hat{\\theta}} (\\theta-\\hat{\\theta})^2+\\cdots\n",
    "$$\n",
    "\n",
    "(the first derivative vanishes by construction because of the ML\n",
    "principle). Then let's approximate the likelihood with a Gaussian in the\n",
    "neighborhood of its maximum: \n",
    "\n",
    "$$\n",
    "L(\\theta)\n",
    "        \\sim const\\cdot \\exp\\left(-\\frac{1}{2}\\cdot\\frac{d^2F}{d\\theta^2}\\lvert_{\\theta=\\hat{\\theta}}(\\theta-\\hat{\\theta})^2\\right) \n",
    "        := const\\cdot \\exp\\left(-\\frac{(\\theta-\\hat{\\theta})^2}{2\\sigma^2}\\right).\n",
    "$$\n",
    "\n",
    "By comparing the exponents we find:\n",
    "\n",
    "$$\n",
    "\\sigma^2(\\hat{\\theta})=\\frac{1}{\\left.d^2F/d\\theta^2\\right|_{\\theta=\\hat{\\theta}}} = \\frac{1}{I(\\hat{\\theta})}.\n",
    "$$\n",
    "\n",
    "The variance is the inverse of the second derivative of the\n",
    "log-likelihood at $\\theta = \\hat{\\theta}$, i.e. the inverse of the\n",
    "information.\n",
    "\n",
    "The difference $F(\\theta) - F(\\hat{\\theta})$ calculated at\n",
    "$\\theta = \\hat{\\theta} \\pm n \\cdot \\sigma(\\hat{\\theta})$, using the\n",
    "equations above is:\n",
    "\n",
    "$$\n",
    "F(\\hat{\\theta} \\pm n\\sigma) - F(\\hat{\\theta}) = \\frac{1}{2}\\frac{d^2F}{d\\theta^2}\\lvert_{\\theta=\\hat{\\theta}} (\\hat{\\theta} \\pm n\\sigma - \\hat{\\theta})^2 = \\frac{1}{2}\\frac{1}{\\sigma^2} ( n\\sigma )^2  = \\frac{1}{2}n^2\n",
    "%F(\\hat{\\theta}\\pm n\\cdot\\sigma) - F(\\hat{\\theta}) = \\frac{1}{2}n^2.\n",
    "$$\n",
    "\n",
    "This enables us to find the uncertainty of an estimator $\\hat{\\theta}$\n",
    "easily by looking at the graph for the log-likelihood function. When the\n",
    "log-likelihood has decreased from the maximum by $0.5$ you are at\n",
    "$\\pm 1 \\sigma$, by $2$ you are at $\\pm2 \\sigma$, by $4.5$ you are at\n",
    "$\\pm 3\\sigma$ and so on.\n",
    "\n",
    "If the log-likelihood function is not parabolic at the maximum then you\n",
    "can try with a non-linear transformation ($\\theta$ goes into\n",
    "$z = z(\\theta)$) such that $F(z)$ shows the desired parabolic behavior:\n",
    "the best estimator is then $\\hat{z} = z(\\hat{\\theta})$ and the standard\n",
    "deviation $\\sigma_{z}$ of $z$ can then be determined as above.\n",
    "\n",
    "If a transformation cannot be found (which is the typical case in any\n",
    "realistic application), you can always proceed numerically and find the\n",
    "values for which the likelihood crosses $1/2 n^2$.\n",
    "\n",
    "Monte Carlo techniques can also be used to estimate the standard\n",
    "deviation or the variance of a parameter. One can simulate a large\n",
    "amount pseudo-experiments (slang: toy data or toys) and for each of them\n",
    "compute the ML estimator: the distribution of the ML estimators is then\n",
    "used to compute the variance. To generate the toy data, one can choose\n",
    "as \"true\" value of the parameter the one from the real experiment and as\n",
    "the size of the sample the number of events of the real experiment.\n",
    "Finally the value of the variance can be computed from\n",
    "$s^2 = 1/(n-1) \\sum (x_i - \\bar{x})^2$ (where $x_i$ are the ML estimates\n",
    "and $i$ runs over the toy datasets) and give this as the statistical\n",
    "error of the parameter estimated from the real measurement.\n",
    "\n",
    "In the case of several parameters\n",
    "$\\theta_{1},\\theta_{2}, \\ldots , \\theta_{m}$ the likelihood function is\n",
    "generalized to\n",
    "\n",
    "$$\n",
    "L(\\theta_1,\\theta_2,\\ldots,\\theta_m)=\\prod_{i=1}^n f(x_i;\\theta_1,\\theta_2,\\ldots,\\theta_m).\n",
    "$$\n",
    "\n",
    "Expanding the NLL function around its minimum at $\\hat{\\theta}$, we\n",
    "obtain (the first derivative vanishes - ML principle): \n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "F(\\theta_1,\\ldots,\\theta_m)&=&F(\\hat{\\theta}_1,\\ldots,\\hat{\\theta}_m)+\\frac{1}{2}\\sum_{i,k}\\frac{\\partial^2F}{\\partial \\theta_i\\partial \\theta_k}(\\theta_i-\\hat{\\theta}_i)(\\theta_k-\\hat{\\theta}_k)+\\cdots \\\\ &=&F(\\hat{\\theta}_1,\\ldots,\\hat{\\theta}_m)+\\frac{1}{2}\\sum_{i,k}G_{ik}(\\theta_i-\\hat{\\theta}_i)(\\theta_k-\\hat{\\theta}_k)+\\cdots\\end{aligned}\n",
    "$$\n",
    "\n",
    "where $G$ is given by: \n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "V&=&G^{-1} \\\\\n",
    "G_{ik}&=&\\frac{\\partial^2 F}{\\partial \\theta_i\\partial \\theta_k },\\end{aligned}\n",
    "$$\n",
    "\n",
    "evaluated at the minimum $\\hat{\\theta}$. In the case of only two\n",
    "parameters the contour lines are drawn as lines with the same likelihood\n",
    "values $F({\\bf \\theta}) = F({\\bf \\hat{\\theta}}) + 1/2 r^{2}$, which\n",
    "correspond to ellipses (see [Error matrix](errors.html#a-closer-look-at-the-error-matrix)\n",
    ")."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Binned Maximum Likelihood\n",
    "-------------------------\n",
    "\n",
    "The likelihood function as we described it in the previous chapter is\n",
    "\"unbinned\". This means that it is constructed out of all available data\n",
    "points $x_{i}$ and therefore no information is lost due to binning. For\n",
    "large data sets, using each single point might very time consuming (in\n",
    "the minimization of the NLL at each variation of the parameters you need\n",
    "to loop over all data points) and it might be practical to bin the data\n",
    "and represent it in histograms. We assume that the random variables\n",
    "$x_{i}$ are distributed according to a p.d.f. $f(x;\\theta)$ and that the\n",
    "expectation value $\\nu = (\\nu_{1}, \\ldots, \\nu_{N})$ for the number of\n",
    "entries per bin $i$ is given by \n",
    "\n",
    "```{math}\n",
    ":label:eq-nuBinned\n",
    "\\nu_i=\\int_{x_i^{min}}^{x_i^{max}}f(x;\\theta) dx.\n",
    "$$ \n",
    "\n",
    "The boundaries of\n",
    "bin $i$ are denoted by $x_{i}^{min}$ and $x_{i}^{max}$, respectively. We\n",
    "can now think of the histogram as some sort of single measurement of a\n",
    "$N$-dimensional random vector for which the combined probability density\n",
    "is given by a multinomial distribution. This means we are asking for the\n",
    "joint probability to observe $n_i$ entries in bin $i$ when the expected\n",
    "is $\\nu_i$. Normalizing by $n_{tot} = \\sum n_{i}$ we get:\n",
    "\n",
    "$$\n",
    "f_{comb}(\\vec{n};\\vec{\\nu})=\\frac{n_{tot}!}{n_1!\\cdots n_N!}\\left(\\frac{\\nu_1}{n_{tot}}\\right)^{n_1}\\cdots \\left(\\frac{\\nu_N}{n_{tot}}\\right)^{n_N}\n",
    "$$\n",
    "\n",
    "Remember that the dependence on the parameter $\\theta$ is embedded in\n",
    "the $\\nu_i$ as in Eq. {eq}`eq-nuBinned`. The negative logarithm of the joint\n",
    "probability yields now the binned NLL function (all uninteresting terms\n",
    "are dropped): \n",
    "\n",
    "$$\n",
    "l(\\theta) = \\sum_{i=1}^N \\ln\\nu_i(\\theta)^{n_i}\n",
    "$$\n",
    "\n",
    "The estimations for $\\hat{\\theta}$ are found as in the unbinned case by\n",
    "minimizing the NLL.\n",
    "\n",
    "Taking the number of bins $N\\to \\infty$ brings back the unbinned\n",
    "likelihood case. Provided that the expected number of entries in a bin\n",
    "is not zero ($\\nu_i(\\theta) > 0$) the binned ML is usable even when some\n",
    "bins have zero entries observed (in contrast with the least square\n",
    "method that we will discussed in the next chapter)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extended Maximum Likelihood Method \n",
    "----------------------------------\n",
    "\n",
    "We applied up to now the ML to a fixed number of events $N$. We can\n",
    "easily extend the ML to the case where the total number of events is\n",
    "itself not known and it is treated as a parameter to be estimated. To do\n",
    "this we can multiply the previous expression of the likelihood by a\n",
    "Poisson p.d.f. which represents the probability to observe $n$ events\n",
    "when the expected number of events is $\\nu$:\n",
    "\n",
    "$$\n",
    "L(x;\\vec{\\theta}) = \\prod_{i=1}^n f(x_i;{\\bf \\theta}) \\to L_E(x;{\\bf \\theta}, \\nu) = \\frac{e^{-\\nu}\\nu^n}{n!} \\prod_{i=1}^n f(x_i;{\\bf \\theta})\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "```{math}\n",
    ":label: eq-eml\n",
    "l(\\vec{x};\\vec{\\theta}) = \\sum_{i=1}^n \\ln f(x_i;{\\bf \\theta}) \\to l_E(\\vec{x};{\\bf \\theta},\\nu)  = \\sum_{i=1}^n \\ln \\nu f(x_i;{\\bf \\theta}) - \\nu + const\\label{eq_eml}\n",
    "```\n",
    "\n",
    "where $\\ln\\nu f(x;{\\bf \\theta})$ is now normalized to $\\nu$ instead of 1\n",
    "and where we dropped the constant term $\\ln n!$ which is irrelevant in\n",
    "the minimization. This new likelihood is called the\n",
    "*extended-maximum-likelihood* or EML.\n",
    "\n",
    "We can now distinguish two cases:\n",
    "\n",
    "Case 1: the parameter $\\nu$ depends on ${\\bf \\theta}$. The EML\n",
    "log-likelihood function can be written as \n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\ln L({\\bf \\theta})&=&n\\ln \\nu({\\bf \\theta})-\\nu({\\bf \\theta})+\\sum_{i=1}^{n}\\ln f(x_i;{\\bf \\theta})\\\\\n",
    "&=&-\\nu({\\bf \\theta})+\\sum_{i=1}^{n}\\ln(\\nu({\\bf \\theta}) f(x_i;{\\bf \\theta}))\\end{aligned}\n",
    "$$\n",
    "\n",
    "where the additive terms not depending on ${\\bf \\theta}$ are dropped. By\n",
    "taking the Poisson term into consideration in the EML function, the\n",
    "resulting variance is usually smaller, because when estimating\n",
    "${\\bf \\hat{\\theta}}$, we use the extra information brought in by $n$.\n",
    "\n",
    "Case 2: $\\nu$ does not dependent on ${\\bf \\theta}$. Differentiating\n",
    "Eq. {eq}`eq-eml` w.r.t $\\nu$ and equating it to zero yields as estimator simply\n",
    "$\\hat{\\nu} = n$, as expected. We also obtain as estimators the same\n",
    "$\\hat{\\theta}$ of the standard ML. Nevertheless the variance of the\n",
    "$\\hat{\\theta}$ would be bigger because now not only $\\hat{\\theta}$ but\n",
    "also $n$ is a source of statistical uncertainty."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Combination of Measurements with the ML Method\n",
    "----------------------------------------------\n",
    "\n",
    "Suppose we have different measurements of the same parameter $\\theta$ by\n",
    "different experiments and you want to combine them using the ML method.\n",
    "More precisely, suppose we have a set of $n$ measured data points with\n",
    "probability density $f(x;\\theta)$ from one experiment and a second set\n",
    "with $m$ measured data points $y_{i}$, which are distributed according\n",
    "to a probability density $g(y;\\theta)$ from a second experiment. The two\n",
    "probability densities $f(x;\\theta)$ and $g(y;\\theta)$ can have different\n",
    "functional forms, because of the different experimental techniques used\n",
    "to determine $\\theta$. As an example you can think of $\\theta$ being the\n",
    "mass of a particle and $f$ and $g$ the results of two experiments or the\n",
    "results of the mass measurement in two decay modes.\\\n",
    "The two experiments together can be interpreted as one single experiment\n",
    "and the resulting likelihood is just the product of:\n",
    "\n",
    "$$\n",
    "L(\\theta)=\\prod_{i=1}^{n}f(x_i;\\theta)\\cdot \\prod_{i=1}^{m}g(y_i;\\theta)=L_x(\\theta)\\cdot L_y(\\theta)\n",
    "$$\n",
    "\n",
    "This expression becomes clear if you think back at the definition of\n",
    "likelihood. The likelihood is based on the conditional probability that,\n",
    "given a parameter $\\theta$ we observe the data set we have. The product\n",
    "above, just extends the conditional probability further to a larger data\n",
    "set comprising two experiments.\\\n",
    "This way of combining different measurements is only valid in the case\n",
    "where the two likelihood are totally uncorrelated, i.e. the two\n",
    "experiments do not share any common source of uncertainty. If that is\n",
    "not the case then the parameters correlation has to be included in the\n",
    "likelihoods expressions. A real life example can be found in the\n",
    "combination of the Higgs mass measurement performed by the ATLAS and CMS\n",
    "collaborations [FIXME @massHiggs]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Constraining parameters\n",
    "-----------------------\n",
    "\n",
    "It often happens that the parameters to be estimated are constrained,\n",
    "for instance by a physical reason (e.g. mass $> 0$) or by other\n",
    "measurements. Imposing constraints always implies adding some\n",
    "information, and therefore the errors of the parameters are in general\n",
    "reduced.\n",
    "\n",
    "The most efficient method to deal with a constraint is to rewrite the\n",
    "parameters such that the constraints are embedded in their definition.\n",
    "\n",
    "```{admonition} Example:\n",
    ":class: tip\n",
    "Take $\\theta_i$ as fractions subjected to the constraint\n",
    "that they should add to 1:\n",
    "\n",
    "$$\n",
    "0 \\le \\theta \\le 1 \\qquad \\; \\qquad \\sum_{i=1}^n \\theta_i =1.\n",
    "$$\n",
    "Then we can redefine the parameters as: \n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\theta_1 &=& \\psi_1\\\\\n",
    "\\theta_2 &=& (1-\\psi_1)\\psi_2\\\\\n",
    "\\theta_3 &=& (1-\\psi_1)(1-\\psi_2)\\psi_3\\\\\n",
    "&\\ldots&\\\\\n",
    "\\theta_{k-1} &=& (1-\\psi_1)(1-\\psi_2)\\ldots(1-\\psi_{k-2})\\psi_{k-1}\\\\\n",
    "\\theta_k &=& (1-\\psi_1)(1-\\psi_2)\\ldots(1-\\psi_{k-1})\\end{aligned}\n",
    "$$\n",
    "\n",
    "where the $\\psi_i \\; \\forall i$ are bounded to be between 0 and 1.\n",
    "```\n",
    "\n",
    "The most general way to express a constraint is through an implicit\n",
    "equation (or in general a set of equations) of the form:\n",
    "$\\vec{g}(\\vec{\\theta}) =0$ and the general method to implement them is\n",
    "to use the *Lagrange multipliers*. Given a likelihood\n",
    "$L(\\vec{x};\\vec{\\theta})$ and the constraint $\\vec{g}(\\vec{\\theta}) = 0$\n",
    "we will find the maximum of:\n",
    "\n",
    "$$\n",
    "F(\\vec{x};\\vec{\\theta}, \\vec{\\alpha}) = \\ln L(\\vec{x};\\vec{\\theta}) + \\vec{\\alpha} \\vec{g}(\\vec{\\theta})\n",
    "$$\n",
    "\n",
    "with respect to $\\vec{\\theta}$ and $\\vec{\\alpha}$. The estimators of\n",
    "$\\vec{\\theta}$ found in this way satisfy the constraints and also have\n",
    "all the usual properties of maximum likelihood estimators.\n",
    "\n",
    "```{admonition} Example:\n",
    ":class: tip\n",
    "Take the likelihood $L(x; \\theta_1; \\theta_2)$ and say we\n",
    "want to estimate $\\theta_1$ but we know from a different measurement\n",
    "that $\\theta_2$ has a a value $\\bar{\\theta}_2 \\pm \\sigma_{\\theta_2}$. We\n",
    "can introduce the constraint on $\\theta_2$ by simply multiplying the\n",
    "likelihood by a gaussian function centred at $\\bar{\\theta}_2$ with width\n",
    "$\\sigma_{\\theta_2}$ (or adding the equivalent parabolic term to the\n",
    "log-likelihood).\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some general remarks concerning ML estimators\n",
    "---------------------------------------------\n",
    "\n",
    "-   For large data sets (large $n$) the ML estimator $\\hat{\\theta}$ is\n",
    "    unbiased and normally distributed around the true value $\\theta$.\n",
    "    The variance approaches the RCF-boundary, i.e. ML estimators are\n",
    "    efficient. They are furthermore consistent for large $n$. These\n",
    "    properties explain the popularity of the ML method.\n",
    "\n",
    "-   A way to study the bias of an estimator $\\hat{\\theta}$ is through\n",
    "    toy MC. Multiple data sets of the same size $n$ of the original one,\n",
    "    all depending on the parameter $\\theta$, have to be generated and\n",
    "    analysed with the ML method. With the obtained results we can\n",
    "    produce the so called \"pull-distribution\", i.e. the histogram filled\n",
    "    with the value $(\\theta - \\hat{\\theta}) / \\sigma_{\\hat{\\theta}}$.\n",
    "    For an unbiased estimator, this distribution should be centred at\n",
    "    zero and have width 1.\n",
    "\n",
    "-   It is important for the likelihood function and the ML estimator\n",
    "    that the probability density $f(x,\\theta)$ is normalized, i.e.\n",
    "    $\\int f(x,\\theta) dx$ is independent from $\\theta$. If the\n",
    "    maximization/minimization is done numerically, it is important to\n",
    "    ensure the normalization at each step (numerical programs like\n",
    "    ROOFIT do that automatically).\n",
    "\n",
    "-   Errors on the ML estimator: For sufficiently large $n$ the\n",
    "    likelihood function $L(\\theta)$, the error $\\sigma$ of the estimator\n",
    "    can be determined from:\n",
    "\n",
    "    -   the values for $\\hat{\\theta}$ for which $\\ln L$ decreases by 0.5\n",
    "\n",
    "    -   $1/\\sqrt{-\\frac{\\partial^2\\ln L}{\\partial a^2}}$\n",
    "\n",
    "    -   the boundaries of the integral $\\int_{a}^{b} L(\\theta) d\\theta$\n",
    "        that encloses 68% of the whole area.\n",
    "\n",
    "    In the case of small $n$, $\\ln L$ is not necessarily parabolic. The\n",
    "    function does not even have to be symmetric around its maximum.\n",
    "    Nevertheless $\\Delta \\ln L = \\frac{1}{2}$ can be used as a rule of\n",
    "    thumb to get an approximation for the 68% confidence interval. For a\n",
    "    better estimation of the uncertainty you can use Monte Carlo toys."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "References\n",
    "----------\n",
    "\n",
    "-   W. Metzger [@Metzger], \"Statistical Methods in Data Analysis\": Ch.8\n",
    "\n",
    "-   G. Cowan [@CowanBook], \"Statistical Data Analysis\",Ch. 6\n",
    "\n",
    "-   R. Barlow [@Barlow], \" A guide to the use of statistical methods in\n",
    "    the physical sciences\". Ch. 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}