{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 7:  Confidence Intervals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### In the exercises of this sheet we consider a measured sample of events after applying quality cuts to separate signal from background events. The total number of observed events follows a Poisson distribution:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation}\n",
    "  P(n_0|\\nu_t) = \\frac{\\nu_t^{n_0} e^{-\\nu_t}}{n_0!},\n",
    "\\end{equation}\n",
    "$n_0$ is the number of observed events after the quality cuts, and $\\nu_t=\\nu_{t,S}+\\nu_{t,B}$ is the true number of events expected to pass the cuts, where $\\nu_{t,S}$ is the contribution from true signal events and $\\nu_{t,B}$ the background remaining after the cuts.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import numpy as np\n",
    "from scipy.stats import poisson, norm\n",
    "from scipy.interpolate import interp1d as scipy_interp1d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 7.1: Significance (voluntary)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### For a specific experiment, the background is expected to be $\\nu_{t,B}=1$ while the measurement is $n_0=3$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### a) How often do you expect to measure $n_0$ or more events if you expect only background?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### b)  What is the corresponding significance?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### c)  Is the measurement compatible with a statistical fluctuation of the background, or is there an evidence (defined as significance $\\ge 3\\,\\sigma$ )or a discovery of the signal ($\\ge 5\\,\\sigma$)?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### d) How many events would have to be observed to reach the evidence and discovery thresholds?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a function to calculate the probability to observe n_0 or more events, if the expected number is mu\n",
    "# Hint: Have a look at the methods of scipy.stats.poisson and in particular the method poisson.cdf.\n",
    "\n",
    "def pPoisson(n_0, mu):\n",
    "    # TODO: Calculate the probability to observe n_0 or more events\n",
    "    #       and also determine the necessary number of observed events for a discovery / evidence.\n",
    "    # Hint: Take a look at the method scipy.stats.norm.ppf to calculate the significances in terms of standard deviations.\n",
    "    pass\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def exercise7_1():\n",
    "    # TODO: Implement your solution to this exercise\n",
    "    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 7.2: Confidence intervals (obligatory)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### a) Frequentist approach (Neyman construction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This classical approach can be used to determine a confidence interval or lower/upper limit, respectively, at a given confidence level (CL).\n",
    "\n",
    "First, let us assume that the background is negligible: $\\nu_{t,B}=0$. Using the frequentist approach, we will compute a $90\\%$ CL upper \n",
    "limit on $\\nu_t=\\nu_{t,S}$ if the measured number of event is $n_0 = 3$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1) Determine the $90\\,\\%$ CL upper limit $\\nu_{UL}$ such that:\n",
    "  \\begin{eqnarray}\n",
    "    \\label{eqs}\n",
    "    \\sum_{n=0}^{n_0} P(n|\\nu_{UL}) & = & 0.10.\n",
    "  \\end{eqnarray}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Now, consider the realistic case with background, i.e., $\\nu_{t,B} > 0$:\n",
    "   \n",
    "   * Compute the $90\\%$ CL upper limits on $\\nu_{t,S}$ as function of $\\nu_{t,B}$.\n",
    "     By convention, the upper limit on $\\nu_{t,S}$ is the limit which would be obtained for $\\nu_{t,S}$ without background minus the number of expected background events,\n",
    "     i.e., $\\nu_{t,B}$, ($CL_{SB}$-limit).\n",
    "     This is done automatically when calling the function `get_upper_poisson_limit` (defined below) with $\\nu_{t,B} > 0$.\n",
    "     $CL_{SB}$ is the confidence level with respect to the signal-plus-background hypothesis and a measure for the compatibility of the experiment with this hypothesis.\n",
    "     **What makes this procedure inconvenient?**\n",
    "  \n",
    "   * Make a plot of the $CL_{SB}$-limit on $\\nu_{t,S}$ as a function of $\\nu_{t,B}$ varying $n_0$ from 0 to 5 (draw all six curves in the same plot).\n",
    "     \n",
    "   * Utilizing the function `get_upper_poisson_limit_normalized` (also defined below) the limit is calculated using the  $CL_{S}$-method:\n",
    "     $CL_{S} = CL_{SB}/CL_{B}$. $CL_{B}$ is a measure for the compatibility with the background-only hypothesis.\n",
    "     The $CL_{S}$-method provides a useful limit, even if the number of observed events is much smaller than the expected background.\n",
    "    \n",
    "   * Create also a plot of the $CL_{S}$-limit on $\\nu_{t,S}$ as a function of $\\nu_{t,B}$ varying $n_0$ from 0 to 5 (draw all six curves in the same plot).\n",
    "   \n",
    "   \n",
    "**HINT:** Read and understand the functions defined below. Use them to solve the following tasks! They use the following definitions:\n",
    "  \n",
    "\\begin{eqnarray}\n",
    "  CL_{SB} & = & \\sum_{n=0}^{n_0} P(n|\\nu_{t,S}+\\nu_{t,B}), \\\\\n",
    "  CL_{B} & = & \\sum_{n=0}^{n_0} P(n|\\nu_{t,B}), \\\\\n",
    "  CL_{S} & = & CL_{SB}\\, /\\, CL_{B}.\n",
    "\\end{eqnarray}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_x(f, y0, x_min, x_max, kw_params=None, n_int=1000):\n",
    "    \"\"\"\n",
    "    Returns an approximation of the x value of the function f for a given y value (y0).\n",
    "    The paremeter params is a dictionary of keyword-value pairs, which are passed on to the function f.\n",
    "    Requires f beeing invertible (monotonic, ...) on the given invervall.\n",
    "    \"\"\"\n",
    "    \n",
    "    x = np.linspace(x_min, x_max, n_int)\n",
    "\n",
    "    if kw_params is None:\n",
    "        y = np.array([f(x_i) for x_i in x])\n",
    "    else:\n",
    "        y = np.array([f(x_i, **kw_params) for x_i in x])\n",
    "    \n",
    "    f_inv = scipy_interp1d(y, x)\n",
    "    \n",
    "    return f_inv(y0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $90\\,\\%$ CL Upper Limit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def exercise7_2a_90cl_upper_limit(method=0):\n",
    "    # TODO: Use the predefined funciton above to solve the exercise.\n",
    "    pass\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $CL_{SB}$ Limit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_poisson_probability(nu_t_s, nu_t_b, n0):\n",
    "    \"\"\"\n",
    "    Calculates the probability to observe n0 or less signal events for given expected signal and expected background\n",
    "    @param nu_t_s: expected signal\n",
    "    @param nu_t_b: number of expected background events\n",
    "    @param n0: number of observed events\n",
    "    \"\"\"\n",
    "    return poisson.cdf(k=n0, mu=nu_t_s + nu_t_b, loc=0)\n",
    "\n",
    "def get_upper_poisson_limit(n0, nu_t_b, cl):\n",
    "    \"\"\"\n",
    "    returns the upper poisson limit (classical)\n",
    "    @param n0: number of observed events\n",
    "    @param nu_t_b: expected background\n",
    "    @param cl: confidence level\n",
    "    \"\"\"\n",
    "    \n",
    "    # targeted pvalue\n",
    "    p_val = 1 - cl\n",
    "    \n",
    "    # probability to observe n0 or less events\n",
    "    \n",
    "    # interpolation limits\n",
    "    x_min = -1\n",
    "    x_max = 10 * (n0 + 1)\n",
    "    \n",
    "    limit = get_x(\n",
    "        f=get_poisson_probability,\n",
    "        y0=p_val,\n",
    "        x_min=x_min,\n",
    "        x_max=x_max,\n",
    "        kw_params={\"nu_t_b\": nu_t_b, \"n0\": n0},\n",
    "    )\n",
    "    return limit\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def exercise7_2a_Csb():\n",
    "    # TODO: Use the predefined funcitons above to solve the exercise.\n",
    "    # TODO: Create plots for 100 points varying the background from 0 to 8\n",
    "    pass\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $CL_{S}$ Limit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_normalized_probability(nu_t_s, nu_t_b, n0):\n",
    "    \"\"\"\n",
    "    Calculates the probability, to observe n0 or less signal events for given expected signal and expected background\n",
    "    @param nu_t_s: expected signal\n",
    "    @param nu_t_b: number of expected background events\n",
    "    @param n0: number of observed events\n",
    "    \"\"\"\n",
    "    \n",
    "    # p_value for signal+background hypothesis\n",
    "    p1 = poisson.cdf(k=n0, mu=nu_t_b + nu_t_s, loc=0)\n",
    "    \n",
    "    # p_value for background only hypothesis\n",
    "    p2 = poisson.cdf(k=n0, mu=nu_t_b)\n",
    "    \n",
    "    return p1 / p2\n",
    "\n",
    "def get_upper_poisson_limit_normalized(n0, nu_t_b, cl):\n",
    "    \"\"\"\n",
    "    Returns the normalized upper poisson limit (classical)\n",
    "    @param n0: number of observed events\n",
    "    @param nu_t_b: expected background\n",
    "    @param cl: confidence level\n",
    "    \"\"\"\n",
    "    \n",
    "    # targeted pvalue\n",
    "    p_val = 1 - cl\n",
    "    \n",
    "    # interpolation limits\n",
    "    x_min = -1\n",
    "    x_max = 10 * (n0 + 1)\n",
    "    \n",
    "    limit = get_x(\n",
    "        f=get_normalized_probability,\n",
    "        y0=p_val,\n",
    "        x_min=x_min,\n",
    "        x_max=x_max,\n",
    "        kw_params={\"nu_t_b\": nu_t_b, \"n0\": n0},\n",
    "    )\n",
    "    \n",
    "    return limit\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def exercise7_2a_Cs():\n",
    "    # TODO: Use the predefined funcitons above to solve the exercise.\n",
    "    # TODO: Create plots for 100 points varying the background from 0 to 8\n",
    "    pass\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### b) Likelihood approach"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The likelihood function for a Poisson process for one single measurement is:\n",
    "\n",
    "  $ L(n_0|\\nu_t) = \\frac{\\nu_t^{n_0} e^{-\\nu_t}}{n_0!}. $\n",
    "\n",
    "where $n_0=3$ is the number of measured events.\n",
    "\n",
    "- Draw the curve $-2\\ln L$ as function of $\\nu_t$, performing a scan \n",
    "over a meaningful range of values. Where is the minimum, $-2\\ln L_{min}$, of this curve?\n",
    "- Calculate the boundaries of the confidence interval at $68\\%$ CL. These are the points for which $2\\cdot \\Delta \\ln L = 1$, where $\\Delta \\ln L$ is defined as $\\ln L_{min} - \\ln L$.\n",
    "- Determine the $90\\%$ CL upper limit, i.e. the point with $\\nu_t>n_0$ for which $2\\cdot \\Delta \\ln L = (1.28)^{2}$.\n",
    "- Plot the likelihood function and the convidence intervals of interest.\n",
    "\n",
    "\n",
    "**Note:** To translate a CL into the proper $\\Delta \\ln L$ cut for a **one-sided** interval, you can use the equation\n",
    "\n",
    "$ 2 \\cdot \\Delta \\ln L = (\\sqrt{2} \\cdot \\texttt{erfinv}(2 \\cdot CL - 1))^2 $,\n",
    "\n",
    "where `erfinv` is the inverse of the error function as provided for instance by `scipy` as `scipy.special.erfinv`.\n",
    "\n",
    "The term which is passed on as argument to the `erfinv` function is equal to $1 - 2 \\cdot{} (1 - CL) = 2 \\cdot{} CL - 1$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.special import erfinv\n",
    "from scipy.optimize import minimize, root_scalar"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Hint:** You can use the `root_scalar` function provided by `scipy.optimize.root_scalar` to find the intersects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"2 * DeltaL = (\\sqrt(2) * erfinv(0.68...))^2 = ({np.sqrt(2) * erfinv(0.68268949214)})^2\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"2 * DeltaL = (\\sqrt(2) * erfinv(2 * 0.9 - 1))^2 = ({np.sqrt(2) * erfinv(2 * 0.90 -1)})^2\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def negative_log_likelihood(nu_t, n0):\n",
    "    \"\"\"\n",
    "    Returns 2 * negative log likelihood for poisson pdf\n",
    "    @param nu_t: number of expected events\n",
    "    @param n0: number of observed events\n",
    "    @return: NNL\n",
    "    \"\"\"\n",
    "    \n",
    "    return -2. * poisson.logpmf(k=n0, mu=nu_t, loc=0)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def exercise7_2b():\n",
    "    pass\n",
    "    # TODO: Implement your solution:\n",
    "    #       - Get the minimum of the likelihood\n",
    "    #       - Get the 68% confidence level interval (points where the the likelihood is by 1 higher than its minimum value\n",
    "    #       - Get the 90% confidence level upper limit\n",
    "    #       - Plot the likelihood function\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### c) Bayesian approach"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Bayesian posterior probability $P(\\nu_t|n_0)$ is given by the Bayes theorem:\n",
    "\n",
    "$$P(\\nu_t|n_0)=\\frac{L(n_0|\\nu_t)\\ \\Pi(\\nu_t)}{\\int_{\\mathrm{all}\\,\\, \\nu_t} L(n_0|\\nu_t)\\ \\Pi(\\nu_t) d\\nu_t}.$$\n",
    "\n",
    "$\\Pi(\\nu_t)$ is called the prior probability on $\\nu_t$ and describes our prior belief about the distribution of this parameter.\n",
    "\n",
    "Try two different priors and compare the results:\n",
    "\n",
    "   i) $\\Pi(\\nu_t) = const.$ for $\\nu_t>0$ and $0$ otherwise,\n",
    "   \n",
    "   ii) $\\Pi(\\nu_t)$ proportional to $1/\\nu_t$ for $\\nu_t>0$ and $0$ otherwise.\n",
    "\n",
    "\n",
    "\n",
    "- Compute and draw the posterior probability for $n_0=3$ and $\\nu_{t,B}=0$.\n",
    "- What are the $90\\%$ credibility upper limits?\n",
    "  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def likelihood_pdf(nu_t,  n0):\n",
    "    \"\"\"\n",
    "    Returns likelihood for poisson pdf\n",
    "    @param nu_t: number of expected events\n",
    "    @param n0: number of observed events\n",
    "    @return: NNL\n",
    "    \"\"\"\n",
    "    \n",
    "    return poisson.pmf(k=n0, mu=nu_t, loc=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def prior_flat_function(nu_t):\n",
    "    pass\n",
    "    # TODO: Implement\n",
    "\n",
    "    \n",
    "def prior_inverse_function(nu_t):\n",
    "    pass\n",
    "    # TODO: Implement\n",
    "\n",
    "\n",
    "def posterior_function_with_flat_prior(x, par):\n",
    "    pass\n",
    "    # TODO: Implement\n",
    " \n",
    "  \n",
    "def posterior_function_with_inverse_prior(x, par):\n",
    "    pass\n",
    "    # TODO: Implement\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def exercise7_2c():\n",
    "    pass\n",
    "    # TODO: Implement your solution:\n",
    "    #       - Calculate the posterior probabilities\n",
    "    #       - Draw the posterior probabilities\n",
    "    #       - Determine the 90 % credibility upper limits\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Finally, compare the upper limits of the methods a), b), c) for $n_0=3$ and $\\nu_{t,B}=0$**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}