AEPSHEP2022/notebook2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# **Statistics Lecture 2 : exercises notebook**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the companion notebook to lecture 3 in the _Practical Statistics_ course, covering the following topics:\n",
    "\n",
    "1. Maximum likelihood for a single Poisson\n",
    "2. Maximum likelihood for binned data\n",
    "3. Hypothesis testing basics\n",
    "4. Testing for discovery\n",
    "\n",
    "First perform the usual imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Maximum likelihood for a single Poisson"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To illustrate Maximum Likelihood in practice, we consider an example with binned data. We take a simple model with signal and background:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2\n"
     ]
    }
   ],
   "source": [
    "s = 2\n",
    "b = 1\n",
    "\n",
    "# Now generate some random data for the s+b hypothesis\n",
    "np.random.seed(1) # Set the random seed to make sure we generate the same data every time\n",
    "n = np.random.poisson(s + b)\n",
    "\n",
    "print(n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The situation above shows a rather standard case where we count events with coming from 2 sources: signal and background. Since we are in a counting situation, the observed event count $n$ is described by a Poisson distribution.\n",
    "\n",
    "We now want to estimate the best-fit signal in this (pseudo) experiment. We know the answer already:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "naïve best-fit =  1\n"
     ]
    }
   ],
   "source": [
    "best_s = n - b\n",
    "print(\"naïve best-fit = \", best_s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the best-fit value is lower than the true value $s = 2$ that was used to generate, since we had a slight negative fluctuation when generating the data.\n",
    "\n",
    "However the goal here is to exercise the maximum-likelihood (ML) estimation techniques that will be useful in more complex examples later. So we are going to compute the ML estimate $\\hat{s}$ of $s$, which (hopefully!) will match what we have just above.\n",
    "\n",
    "We will define the likelihood $L(s)$, and then move to the quantity $\\lambda = -2 \\log L(s)$ that will be more useful later.\n",
    "We can use either, since maximizing $L(s)$ is equivalent to minimizing $\\lambda(s)$, and the ML estimate $\\hat{s}$ at the maximum/minimum is the same either way,\n",
    "\n",
    "Let's now compute the ML estimate (best fit) $\\hat{s}$ using $\\lambda(s)$ :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a scan over s values:\n",
    "s_values = np.arange(-0.5, 3.5, 0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathbf{\\color{red}{\\rightarrow\\text{Compute }\\lambda\\text{ for each value of }s \\text{, plot the result and find the best-fit value among the values tested.}}}$\n",
    "<span style=\"color:red\">Note that the likelihood $L(s)$ is just a Poisson PDF, and recall that in scipy you evaluate a Poisson PDF using ```scipy.stats.poisson.pmf(nobs, mean)``` where ```mean``` is the Poisson parameters (i.e. the mean and variance of the distribution) and ```nobs``` is the observed number of events.</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathbf{\\color{green}{Answer}}$ : "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "hat{s}, ML = (-0.5, -999)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEGCAYAAABCa2PoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAARHUlEQVR4nO3de6xlZX3G8e8jg3jhYguDwAx2MKAteIVdgtG2pE5gqlEQIY5NCq00UwlNNTY2UIyX1KTSWi/YipmIFrQVCS0FBaTgpTQNBc8QlBlujopyMqSMQrmURgR//WMv7PGwz5zLvPusM5zvJ9k5a73vetf67Zcz8+x1mU2qCkmSdtYz+i5AkvT0YKBIkpowUCRJTRgokqQmDBRJUhMr+i6gL/vtt1+tWbOm7zIkaZeyadOmH1XVylF9yzZQ1qxZw8TERN9lSNIuJckPZurzkpckqQkDRZLUhIEiSWrCQJEkNWGgSJKaMFAkSU0YKJKkJgwUSVITBookqQkDRZLUhIEiSWrCQJEkNWGgSJKaMFAkSU0YKJKkJgwUSVITBookqQkDRZLUhIEiSWrCQJEkNWGgSJKaMFAkSU0YKJKkJgwUSVITBookqQkDRZLUhIEiSWpiyQVKkpcnuSHJrUm+lGTvKX0v6/q2dP3P6tqP6ta3JjkvSfp7B5K0PC25QAE+DZxVVS8FLgPeDZBkBfB54O1VdQRwLPDTbsz5wAbgsO61bpFrlqRlbykGyouB67vla4E3d8vHAd+uqm8BVNWPq+qJJAcCe1fVDVVVwEXAiYtcsyQte0sxUDYDb+yWTwEO7pZfBFSSa5LcnOTPuvZVwOSU8ZNd21Mk2ZBkIsnE9u3bx1C6JC1fK/o4aJLrgANGdJ0DvA04L8l7gSuAx7q+FcBrgF8HHgW+mmQT8NCI/dSo41bVRmAjwGAwGLmNJGlhegmUqlo7yybHASR5EfD6rm0S+Leq+lHXdxVwJMP7KqunjF0NbGtasCRpVkvukleS/bufzwDeA3yq67oGeFmS53Q36H8LuK2q7gUeTnJM93TXqcDlPZQuScvakgsU4K1J7gLuYHim8VmAqnoA+AjwTeAW4OaqurIbcwbDp8O2At8Frl7kmiVp2cvwwajlZzAY1MTERN9lSNIuJcmmqhqM6luKZyiSpF2QgSJJasJAkSQ1YaBIkpowUCRJTRgokqQmDBRJUhMGiiSpCQNFktSEgSJJasJAkSQ1YaBIkpowUCRJTRgokqQmDBRJUhMGiiSpCQNFktSEgSJJasJAkSQ1YaBIkpowUCRJTRgokqQmDBRJUhMGiiSpCQNFktSEgSJJasJAkSQ1YaBIkpowUCRJTRgokqQmllygJHl5khuS3JrkS0n27tp3T3Jh1357krOnjDmqa9+a5Lwk6e8dSNLytOQCBfg0cFZVvRS4DHh3134KsEfXfhTwR0nWdH3nAxuAw7rXukWtWJK0JAPlxcD13fK1wJu75QKem2QF8GzgMeChJAcCe1fVDVVVwEXAiYtbsiRpKQbKZuCN3fIpwMHd8qXA/wD3Aj8EPlxV9wOrgMkp4ye7tqdIsiHJRJKJ7du3j6N2SVq2egmUJNcl2TzidQLwNuDMJJuAvRieiQAcDTwBHAQcAvxpkhcCo+6X1KjjVtXGqhpU1WDlypXN35ckLWcr+jhoVa2dZZPjAJK8CHh91/a7wFeq6qfAfUn+AxgA/w6snjJ2NbCtbcWSpNksuUteSfbvfj4DeA/wqa7rh8BvZ+i5wDHAHVV1L/BwkmO6p7tOBS7voXRJWtaWXKAAb01yF3AHwzONz3btfwfsyfAeyzeBz1bVt7u+Mxg+HbYV+C5w9aJWLEnq55LXjlTVx4GPj2h/hOFN+lFjJoCXjLk0SdIOLMUzFEnSLshAkSQ1YaBIkpowUCRJTRgokqQmDBRJUhMGiiSpCQNFktSEgSJJasJAkSQ1YaBIkpowUCRJTRgokqQmDBRJUhMGiiSpCQNFktSEgSJJasJAkSQ1YaBIkpowUCRJTRgokqQmDBRJUhMGiiSpCQNFktSEgSJJasJAkSQ1YaBIkpowUCRJTRgokqQmDBRJUhMGiiSpiRV9HDTJKcD7gV8Djq6qiSl9ZwOnA08Af1JV13TtRwF/DzwbuAp4R1VVkj2Ai4CjgB8Db6mqu8dV+we+tIXbtj00rt1L0tgdftDevO8NRzTf74LPUJK8cyeOuxk4Cbh+2j4PB9YDRwDrgE8m2a3rPh/YABzWvdZ17acDD1TVocBHgXN3oi5J0gLtzBnKu4CPLWRgVd0OkGR61wnAxVX1E+D7SbYCRye5G9i7qm7oxl0EnAhc3Y15fzf+UuBvk6SqaiG1zWYcqS5JTwc7cw/lKWnQwCrgninrk13bqm55evsvjKmqx4EHgX1H7TzJhiQTSSa2b9/euHRJWt525gxlh2cASa4DDhjRdU5VXT7TsBmOM1P7jsY8tbFqI7ARYDAYjOUMRpKWqx0GSpKHGf2XcxjeHJ9RVa1dQD2TwMFT1lcD27r21SPap46ZTLIC2Ae4fwHHliTthB1e8qqqvapq7xGvvapqHE+IXQGsT7JHkkMY3ny/qaruBR5OckyGN15OBS6fMua0bvlk4Gvjun8iSZpZX48Nvwn4BLASuDLJLVV1fFVtSXIJcBvwOHBmVT3RDTuD/39s+OruBXAB8LnuBv79DJ8SkyQtsizXD/ODwaAmJiZm31CS9HNJNlXVYFSf/1JektSEgSJJasJAkSQ1YaBIkpowUCRJTRgokqQmDBRJUhMGiiSpCQNFktSEgSJJasJAkSQ1YaBIkpowUCRJTRgokqQmDBRJUhMGiiSpCQNFktSEgSJJasJAkSQ1YaBIkpowUCRJTRgokqQmDBRJUhMGiiSpCQNFktSEgSJJasJAkSQ1YaBIkpowUCRJTRgokqQmegmUJKck2ZLkZ0kG0/rOTrI1yZ1Jju/anpPkyiR3dOM+NGX7PZJ8sRtzY5I1i/x2JEn0d4ayGTgJuH5qY5LDgfXAEcA64JNJduu6P1xVvwq8Enh1kt/p2k8HHqiqQ4GPAucuQv2SpGl6CZSqur2q7hzRdQJwcVX9pKq+D2wFjq6qR6vq693Yx4CbgdVTxlzYLV8KvDZJxvsOJEnTLbV7KKuAe6asT3ZtP5fkecAbgK9OH1NVjwMPAvuOu1BJ0i9aMa4dJ7kOOGBE1zlVdflMw0a01ZR9rgC+AJxXVd+by5hpNW0ANgC84AUvmKEESdJCjC1QqmrtAoZNAgdPWV8NbJuyvhH4TlV9bMSYyS5w9gHun6Gmjd0+GAwGI0NHkrQwS+2S1xXA+u7JrUOAw4CbAJJ8kGFYvHPEmNO65ZOBr1WVYSFJi6yvx4bflGQSeBVwZZJrAKpqC3AJcBvwFeDMqnoiyWrgHOBw4OYktyT5w253FwD7JtkKvAs4a5HfjiQJyHL9MD8YDGpiYqLvMiRpl5JkU1UNRvUttUtekqRdlIEiSWrCQJEkNWGgSJKaMFAkSU0YKJKkJgwUSVITBookqQkDRZLUhIEiSWrCQJEkNWGgSJKaMFAkSU0YKJKkJgwUSVITBookqQkDRZLUhIEiSWrCQJEkNWGgSJKaMFAkSU0YKJKkJgwUSVITBookqQkDRZLUhIEiSWrCQJEkNWGgSJKaMFAkSU0YKJKkJgwUSVITvQRKklOSbEnysySDaX1nJ9ma5M4kx48Ye0WSzVPW90jyxW7MjUnWLMJbkCRN09cZymbgJOD6qY1JDgfWA0cA64BPJtltSv9JwCPT9nU68EBVHQp8FDh3jHVLkmbQS6BU1e1VdeeIrhOAi6vqJ1X1fWArcDRAkj2BdwEfHDHmwm75UuC1STKeyiVJM1lq91BWAfdMWZ/s2gD+Avgb4NGZxlTV48CDwL6jdp5kQ5KJJBPbt29vWbckLXtjC5Qk1yXZPOJ1wo6GjWirJK8ADq2qy+Y6ZtTOq2pjVQ2qarBy5crZ34Qkac5WjGvHVbV2AcMmgYOnrK8GtgGvAo5KcjfDmvdP8o2qOnbKmMkkK4B9gPt3onRJ0gIstUteVwDruye3DgEOA26qqvOr6qCqWgO8BrirC5Mnx5zWLZ8MfK2qRp6hSJLGZ2xnKDuS5E3AJ4CVwJVJbqmq46tqS5JLgNuAx4Ezq+qJWXZ3AfC5JFsZnpmsH2ftkqTRslw/zA8Gg5qYmOi7DEnapSTZVFWDUX1L7ZKXJGkXZaBIkpowUCRJTRgokqQmDBRJUhMGiiSpCQNFktSEgSJJasJAkSQ1YaBIkpowUCRJTRgokqQmDBRJUhMGiiSpCQNFktSEgSJJasJAkSQ1YaBIkpowUCRJTRgokqQmDBRJUhMGiiSpiVRV3zX0Isl24AcLHL4f8KOG5bRkbQtjbQtjbQuzK9f2K1W1clTHsg2UnZFkoqoGfdcxirUtjLUtjLUtzNO1Ni95SZKaMFAkSU0YKAuzse8CdsDaFsbaFsbaFuZpWZv3UCRJTXiGIklqwkCRJDVhoMxBkl9Ocm2S73Q/f2mG7e5OcmuSW5JMjLmmdUnuTLI1yVkj+pPkvK7/20mOHGc986zt2CQPdvN0S5L3LlJdn0lyX5LNM/T3NmdzrK+veTs4ydeT3J5kS5J3jNiml7mbY219zduzktyU5FtdbR8YsU1f8zaX2uY/b1Xla5YX8FfAWd3yWcC5M2x3N7DfItSzG/Bd4IXAM4FvAYdP2+Z1wNVAgGOAGxdpruZS27HAl3v47/ibwJHA5hn6e5mzedTX17wdCBzZLe8F3LWEft/mUltf8xZgz255d+BG4JglMm9zqW3e8+YZytycAFzYLV8InNhfKQAcDWytqu9V1WPAxQxrnOoE4KIa+k/geUkOXCK19aKqrgfu38Emfc0ZMKf6elFV91bVzd3yw8DtwKppm/Uyd3OsrRfdXDzSre7evaY/BdXXvM2ltnkzUObm+VV1Lwx/gYH9Z9iugH9NsinJhjHWswq4Z8r6JE/9QzSXbcZhrsd9VXe6fXWSIxahrrnoa87mo9d5S7IGeCXDT7RT9T53O6gNepq3JLsluQW4D7i2qpbMvM2hNpjnvK1oXeSuKsl1wAEjus6Zx25eXVXbkuwPXJvkju5TZ2sZ0Tb908VcthmHuRz3ZobfB/RIktcB/wIcNu7C5qCvOZurXuctyZ7APwHvrKqHpnePGLJoczdLbb3NW1U9AbwiyfOAy5K8pKqm3iPrbd7mUNu8580zlE5Vra2ql4x4XQ7815Onod3P+2bYx7bu533AZQwv/4zDJHDwlPXVwLYFbDMOsx63qh568nS7qq4Cdk+y3yLUNpu+5mxO+py3JLsz/Av7H6rqn0ds0tvczVbbUvh9q6r/Br4BrJvW1fvv3Ey1LWTeDJS5uQI4rVs+Dbh8+gZJnptkryeXgeOAkU/rNPBN4LAkhyR5JrC+q3F6zad2T5EcAzz45GW7MZu1tiQHJEm3fDTD38MfL0Jts+lrzuakr3nrjnkBcHtVfWSGzXqZu7nU1uO8rew+/ZPk2cBa4I5pm/U1b7PWtpB585LX3HwIuCTJ6cAPgVMAkhwEfLqqXgc8n+FpIwzn9R+r6ivjKKaqHk/yx8A1DJ+q+kxVbUny9q7/U8BVDJ8g2Qo8CvzBOGpZYG0nA2ckeRz4X2B9dY+VjFOSLzB8cmW/JJPA+xjejOx1zuZRXy/zBrwa+D3g1u6aO8CfAy+YUltfczeX2vqatwOBC5PsxvAv40uq6stL4c/pHGub97z51SuSpCa85CVJasJAkSQ1YaBIkpowUCRJTRgokqQmDBRJUhMGiiSpCQNFWiK6b1u4svsyvs1J3tJ3TdJ8+C/lpaVjHbCtql4PkGSfnuuR5sUzFGnpuBVYm+TcJL9RVQ/2XZA0HwaKtERU1V3AUQyD5S+zSP+rWqkVL3lJS0T3ZaP3V9XnkzwC/H7PJUnzYqBIS8dLgb9O8jPgp8AZPdcjzYvfNixJasJ7KJKkJgwUSVITBookqQkDRZLUhIEiSWrCQJEkNWGgSJKa+D+wolgOu1iQogAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def L(s, nobs) :\n",
    "    return -999 # change this to python code that will compute the Poisson likelihood (or PDF) for nobs and mean=s+b\n",
    "\n",
    "def lambda_s(s, nobs) :\n",
    "    return -999 # change this to python code that will compute -2 log L\n",
    "\n",
    "lambda_values = np.array([ lambda_s(s_value, n) for s_value in s_values ]) # compute the likelihood (the Poisson PDF for each s)\n",
    "\n",
    "# Plot the result\n",
    "plt.plot(s_values, lambda_values)\n",
    "plt.xlabel('s')\n",
    "plt.ylabel('lambda(s)')\n",
    "\n",
    "# Find the maximum likelihood value in the list : merge the 2 lists and find the max based on the L value\n",
    "from operator import itemgetter\n",
    "print('hat{s}, ML =', min( [ (s,l) for s,l in zip(s_values, lambda_values) ], key=itemgetter(1)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If all goes well, you should find a best-fit $\\hat{s}$ of about 1.\n",
    "\n",
    "You also get the value of $\\lambda(\\hat{s})$, but this is not very interesting. Since $L(s)$ was defined only up to a multiplicative factor, $\\lambda(s)$ is defined up to an additive constant, so the value $\\lambda(\\hat{s})$ at the minimum isn't meaningful. On the other hand, the difference between values of $\\lambda(s)$ at different $s$ is meaningful, and will be used later in hypothesis testing.\n",
    "\n",
    "Of course the $\\hat{s}$ value we got was a bit imprecise due to the step size of 0.1 in the scan, but we can do better:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "exact hat{s}: 2.6180339603380443\n"
     ]
    }
   ],
   "source": [
    "from scipy.optimize import minimize_scalar\n",
    "# Compute the value at minimum (.x selects the position of the minimum, .fun would return the function value)\n",
    "s_hat = minimize_scalar(lambda s : lambda_s(s, n), bounds=(0, 100)).x\n",
    "print('exact hat{s}:', s_hat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should find a best-fit value of exactly $\\hat{s} = 1$, which is the correct value as we know since the beginning."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Maximum likelihood for binned data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now a consider a more realistic example with binned data. Since there are multiple bins (here 10), there's no way to guess the answer and we'll really need the ML to get the best-fit signal. We again take a simple model with both a signal and a background component:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbfklEQVR4nO3de5RcZZnv8e8vlyG0hItJhiOEpCJCCCfkguEyoBC5yLDCCQNhKdBRxDDtIEdQBwQmKlGmR1mgI2uWMiaAgLaCB5GLMKgIHUANkIQkXGJGwSRGcAhBQi4mpM1z/qjqpJN0d3Z31967u/bvs1at1N61L09d8tTut973eRURmJlZcfTLOwAzM8uWE7+ZWcE48ZuZFYwTv5lZwTjxm5kVzIC8A0hi6NChUSqV8g7DzKxPWbBgwesRMWzn9X0i8ZdKJebPn593GGZmfYqkFe2td1OPmVnBOPGbmRWME7+ZWcH0iTZ+M7P2bNmyhVWrVrFp06a8Q8nVoEGDGD58OAMHDky0fWqJX9KtwBnAaxExtrLuncBdQAlYDnwoIv6cVgxmVttWrVrF4MGDKZVKSMo7nFxEBGvWrGHVqlWMGjUq0T5pNvXcBvz9TuuuAn4REYcAv6gsm5l1y6ZNmxgyZEhhkz6AJIYMGdKlv3pSS/wR8Tjwxk6rzwRur9y/HfiHtM5vZsVQ5KTfqquvQdY/7u4fEa8CVP792442lNQgab6k+atXr84sQDOzWtdre/VExOyImBQRk4YN22XgmZlZrzRr1ixuuOGGDh+/9957efHFFzOMaFdZJ/7/kfQugMq/r2V8frNeY/LkyUyePDnvMCxjRUz89wMXVO5fANyX8fnNrMCampoolUr069ePUqlEU1NTVY7b2NjI6NGjOeWUU1i2bBkAc+bM4aijjmL8+PFMmzaNjRs38qtf/Yr777+fK664ggkTJvDSSy+1u13aUkv8kn4A/BoYLWmVpBnAV4FTJf0WOLWybGaWuqamJhoaGlixYgURwYoVK2hoaOhx8l+wYAF33nknzz77LPfccw/PPPMMAGeffTbPPPMMixcvZsyYMdxyyy0cd9xxTJ06leuvv55FixZx8MEHt7td2lLrxx8R53Xw0MlpndPMrCMzZ87c5Wp648aNzJw5k/r6+m4f94knnuCss86irq4OgKlTpwLw/PPP8/nPf54333yT9evXc9ppp7W7f9Ltqskjd82sEFauXNml9V3RXnfKj33sY9x7772MHz+e2267jebm5nb3TbpdNfXaXj1mZtU0YsSILq1P6oQTTuDHP/4xf/nLX1i3bh0PPPAAAOvWreNd73oXW7Zs2aE5afDgwaxbt27bckfbpcmJ38wKobGxcVtzTKu6ujoaGxt7dNwjjzySD3/4w0yYMIFp06bx/ve/H4Brr72WY445hlNPPZXDDjts2/bnnnsu119/PRMnTuSll17qcLs0KSIyOVFPTJo0KTwRi9Wa1q6cWfxpX6uWLl3KmDFjEm/f1NTEzJkzWblyJSNGjKCxsbFH7fu9SXuvhaQFETFp523dxm9mhVFfX18zib4n3NRjZlYwTvxmZgXjxG9mVjBO/GZmBePEb2ZWMO7VY2a14/tVnpTl/HS6u8+aNYu99tqLyy+/PJXj746v+M3MUlQqlfIOYRdO/GY5aGpqYt68ecydO7eq5YEtexs2bGDKlCmMHz+esWPHctdddyXab/HixZx00kkccsghzJkzJ+Uod+SmHrOMtZYH3rx5M8C28sCABxf1QQ8//DAHHHAADz74IABr165NtN+SJUuYN28eGzZsYOLEiUyZMoUDDjggzVC38RW/WcY6Kw9sfc8RRxzBI488wpVXXskTTzzBPvvswyWXXMKECROYMGECr7zyyrb7besCnXnmmey5554MHTqUD3zgAzz99NOZxewr/p24foqlLc3ywJa9Qw89lAULFvDQQw9x9dVX88EPfpBvfvOb2x4vlUosWrRol/12LuXcXmnntPiK3yxjaZUHtny88sor1NXVMX36dC6//HIWLlyYaL/77ruPTZs2sWbNGpqbmznqqKNSjnQ7X/GbZayxsZGGhoYdmnuqUR7YSK37ZWeee+45rrjiCvr168fAgQO56aabEu139NFHM2XKFFauXMkXvvCFzNr3wYnfLHOtP+DOmDGDzZs3M3LkyJoqD1w0p512WqfTJS5fvnyXdbNmzUovoASc+M1yUF9fv60Ln39Psqy5jd/MrGCc+M3MCsaJ38ysYJz4e5HJkydvG0dgZpYWJ34zs4Jx4jezmiFV97Y7y5cvZ+zYsVV/Hs3NzZxxxhlVP24rJ34zs4Jx4jcz64GWlhYuuOACxo0bxznnnLNLAT6ARYsWceyxxzJu3DjOOuss/vznPwPl3/WuvPJKjj76aA499FCeeOKJHfbbunUrhxxyCKtXr962/J73vIfXX3+9RzE78ZuZ9cCyZctoaGhgyZIl7L333nzrW9/aZZuPfvSjXHfddSxZsoQjjjiCL33pS9sea2lp4emnn+Yb3/jGDusB+vXrx/Tp07fN1/DII48wfvx4hg4d2qOYnfjNzHrgoIMO4vjjjwdg+vTpPPnkkzs8vnbtWt58801OPPFEAC644AIef/zxbY+fffbZALz3ve9tt7zDxz/+ce644w4Abr31Vi688MIex+zEb2bWAz0tr7zHHnsA0L9/f1paWnZ5/KCDDmL//ffn0Ucf5amnnuL000/vfrAVTvwGeAyBWXetXLmSX//61wD84Ac/4H3ve98Oj++zzz7st99+29rvv/vd7267+k/qoosuYvr06XzoQx+if//+PY7Zid/MakZEdW9JjBkzhttvv51x48bxxhtvcPHFF++yze23384VV1zBuHHjWLRoEV/84he79LymTp3K+vXrq9LMAzlV55T0GeAiIIDngAsjYlMesfQWrZNvb968mVKp5DK9Zn1AqVTixRdf3O12EyZMYN68ebusb1uZdejQodva+Hf+C3zx4sWMHz+eww47rKchAzlc8Us6ELgUmBQRY4H+wLlZx9GbdDT5dusv+WZWXF/96leZNm0aX/nKV6p2zLyaegYAe0oaANQBr+QUR6/gybfNrCNXXXUVK1as2OW3g57IPPFHxB+BG4CVwKvA2oj42c7bSWqQNF/S/NbBC7XKk2+bdV8kbYyvYV19DfJo6tkPOBMYBRwAvEPS9J23i4jZETEpIiYNGzYs6zAz5cm3zbpn0KBBrFmzptDJPyJYs2YNgwYNSrxPHj/ungL8PiJWA0i6BzgO+F4OsfQKnnzbrHuGDx/OqlWrqPVWgd0ZNGgQw4cPT7x9Hol/JXCspDrgL8DJwPwc4ug1PPm2WfcMHDiQUaNG5R1Gn5N54o+IpyTdDSwEWoBngdlZx9HbFHny7dZua3k87zzPbZaXXPrxR8Q1wDV5nNvMrOg8ctfMrGCc+M3MCsaJ38ysYJz4zcwKxonfzKxgOuzVI+kBytUz2xURU1OJyMzMUtVZd84bKv+eDfwvto+sPQ9YnmJMZmaWog4Tf0TMBZB0bUSc0OahByQ93sFuZmbWyyVp4x8m6d2tC5JGAbVdNc3MrIYlGbn7GaBZ0suV5RLwidQiMjOzVO028UfEw5IOAVrn/PpNRGxONyyzPuj76tr2r3Vjv/OLW37YqqezXj1nd/DQwZKIiHtSisnMzFLU2RX//+nksQCc+GuEJ3o3K5bOevVcmGUglo+OJnoHnPzNatRue/VI2l/SLZL+q7J8uKQZ6YdmWfBE72bFk6Q7523ATynPjwvw38CnU4rHMuaJ3otp8uTJ2yahseJJkviHRsQPga0AEdEC/DXVqCwznujdrHiSJP4NkoZQqdsj6VhgbapRWWYaGxupq6vbYZ0nejerbUkGcH0WuJ9yN85fUh61e06qUVlmPNG7WfEkGcC1UNKJwGhAwLKI2JJ6ZJaZIk/0blZEnQ3gOikiHm1nINehHsBlZtZ3dXbFfyLwKO0P5PIALjOzPqqzAVzXVO5eFBHuxWNmViOS9Or5naTrJR2eejRmZpa6JIl/HOVBWzdLmiepQdLeKcdlZmYp2W3ij4h1ETEnIo4DPgdcA7wq6XZJ70k9QjMzq6rddueU1B+YAlxIeRKWrwFNwPuBh4BDU4yvx9TFEund2S9cIt3M+pAkA7h+CzwGXB8Rv2qz/m5JJ3Swj5mZ9VJJEv+4iFjf3gMRcWmV4zEzs5QlaeNvN+mbmVnflKRXj5mZ1RAnfjPLlOcCyF9ntXo+29mOEfH16odjZmZp6+yKf3DlNgm4GDiwcvsnoEejeCXtK+luSb+RtFTS3/XkeL2V1LXb3LnlW1f26etaJ3qfO3cupVKJpqamvEOqeX7NrbNaPV8CkPQz4MiIWFdZngX8vx6e90bg4Yg4R9LfAHW728Fqjyd6z55fc4NkbfwjgLfbLL9NeSBXt1TKPZwA3AIQEW9HxJvdPZ71XZ7oPXt+zQ2S9eP/LvC0pB9TLsd8FnBHD875bmA18B1J44EFwGURsaHtRpIagAbw/K+1yhO9Z8+vuUGyfvyNlMs1/Bl4E7gwIv6tB+ccABwJ3BQRE4ENwFXtnHd2REyKiEnDhg3rwemst/JE79nza26QvDtnHfBWRNwIrJI0qgfnXAWsioinKst3U/4isILxRO/Z82tukCDxS7oGuBK4urJqIPC97p4wIv4E/EHS6Mqqk4EXu3s867vq6+uZPXs2e+yxBwAjR45k9uzZ/pExRUV/zT2GoCxJG/9ZwERgIUBEvCJpcA/P+ymgqdKj52XKTUlWRX2lKqknes+eX3NLkvjfjoiQFACS3tHTk0bEIsrjA8zMLGNJ2vh/KOnbwL6S/hF4BLg53bDMzCwtu73ij4gbJJ0KvAWMBr4YET9PPTIzM0tFkhm4rouIK4Gft7POzMz6mCRNPae2s+70agdiZmbZ6Kw658XAJ4GDJS1p89Bg4Fft72XWd3oUmRVVZ0093wf+C/gKO46sXRcRb6QalVl3fb+L3zqvdWO/8/2tY31bh009EbE2IpZTrqT5RkSsiIgVwBZJx2QVoJlZLehNg8eStPHfBLSdd3dDZZ2ZmfVBSRK/Ira3qEbEVpIN/DIzs14oSeJ/WdKlkgZWbpdRLrNgZmZ9UJLE/0/AccAfKVfWPIZKnXwzM+t7kozcfQ04N4NYzHpM9V3tcTO5sl9z4j3i/C6ewqyXSVKW+VBJv5D0fGV5nKTPpx+amZmlIUlTzxzKtfi3AETEEvwXgJlZn5Uk8ddFxNM7rWtJIxgzM0tfksT/uqSDKU+0jqRzgFdTjcrMzFKTpD/+JcBs4DBJfwR+DxRjnjYzsxqUpFfPy8AplZm3+kXEuvTDMjOztCSpxz8EuAZ4HxCSngS+HBFr0g7OrC/JvStpVwrNuThdoSVp478TWA1MA86p3L8rzaDMzCw9Sdr43xkR17ZZ/ldJ/5BSPGZmlrIkV/yPSTpXUr/K7UPAg2kHZmZm6UiS+D9BeVKWzZXbncBnJa2T9FaawZmZWfUl6dUzOItAzMwsG0lq9czYabm/pGvSC8msKJorN7NsJflx92RJ04AZwFDgVmBuqlGZWZd1rTvp5Mo+zYn3cFXS2pGkqed8SR8GngM2AudFxC9Tj8zMzFKRpKnnEOAy4EfAcuAjkupSjsvMzFKSpFfPA8AXI+ITwInAb4FnUo3KzMxSk6SN/+iIeAugMun61yTdn25YZmaWliRX/HtKukXSwwCSDgdOSDcsMzNLS5Ir/tuA7wAzK8v/TblWzy0pxWRmfYy6UOutO/uE68NVVZIr/qER8UNgK0BEtAB/7emJK+MBnpX0k54eq3qagHmUe6uWKstWu/x+WzaampqYN28ec+fOpVQq0dSU72ctSeLfUCnN3DoD17HA2iqc+zJgaRWOUyVNQAPlqhQAKyrLTga1ye+3ZaOpqYmGhgY2by5/1lasWEFDQ0OuyT9J4v8scD9wsKRfAncAn+rJSSUNB6YAN/fkONU1k/IwhbY2sr2Fy2qL32/LxsyZM9m4ccfP2saNG5k5M7/PWpIBXAslnQiMBgQsi4gtPTzvN4DPAR3WAZLUQPkSjBEjRvTwdEms7OJ669v8fls2Vq5s/zPV0fosJLniJyJaIuKFiHi+p0lf0hnAaxGxYDfnnB0RkyJi0rBhw3pyyoQ6+nLJ4kvHsuf327LR0YVrNhe07UuU+KvseGCqpOWUSzyfJOl7OcSxk0Zg5wHJdZX1Vnv8fls2Ghsbqavb8bNWV1dHY2N+n7Uk3TmrKiKuBq4GkDQZuDwipmcdx67qK//OpPzn/gjKSaC+wz2sL2t9X2dQ/oF3JH6/+6be3pW0vr78mZoxYwabN29m5MiRNDY2blufhw4Tv6QjO9sxIhZWP5y81eP/+EVSD8yp3G/OMQ6rdfX19cyZU/6sNTc35xsMnV/xf62TxwI4qacnj4hm/D/OzCxTHSb+iPhAloGYmVk2ErXxSxoLHA4Mal0XEXekFZSZmaVnt4m/Ms3iZMqJ/yHgdOBJygO5rKqa8w7AzAogSXfOc4CTgT9FxIXAeGCPVKMyM7PUJGnq+UtEbJXUImlv4DXg3SnHZWbW66XdlRTSqUyaJPHPl7Qv5X5vC4D1wNPVD8XMzLKQpFbPJyt3/7MyGcveEbEk3bDMzCwtSSZb/0Xr/YhYHhFL2q6zWtGMf1w2K4bORu4Ooly8ZKik/ShX5gTYGzggg9jMzCwFnTX1fAL4NOUk37Y8w1vAN1OMyczMUtTZyN0bgRslfSoi/iPDmMzMLEVJevV8W9KlwAmV5Wbg21WYjMXMzHKQJPF/CxhY+RfgI8BNwEVpBWVF05x3AFYITcA8ymW4SxS5DHdnP+4OiIgW4KiIGN/moUclLU4/NLMsNOcdgGWiifJMrpsryysqy1DE5N9Zd87WQVp/lXRw60pJ7wb+mmpUZmZVNRPYuNO6jZX1xdNZU09r983LgcckvVxZLgEXphmUmVl1dTSxeX4Tnueps8Q/TNJnK/e/DfQHNlAuzTwReCzl2MzMqmQE5ead9tYXT2dNPf2BvYDBlL8gVFkeUFlnZtZHNFIej9pWXWV98XR2xf9qRHw5s0jMzFLT+gPuDMo/8I7EvXra142Co2ZmvVU95SLDUPTeXJ019ZycWRRmZpaZDhN/RLyRZSBmZpaNJFMvmplZDUlSssHMak5z3gFYjnzFb2ZWME78ZmYF48RvZlYwTvxmZgXjxG9mVjBO/GZmBePEb2ZWME78ZmYF48RvZlYwmSd+SQdJekzSUkkvSLos6xjMLC9NwHLgccqT+TXlGUxh5VGyoQX454hYKGkwsEDSzyPixRxiMbPMtE543jr3bbEnPM9T5lf8EfFqRCys3F8HLAUOzDoOM8uaJzzvLXJt45dUojx/71PtPNYgab6k+atXr848NjOrNk943lvklvgl7QX8CPh0RLy18+MRMTsiJkXEpGHDhmUfoJlVWUcTmxdzwvM85ZL4JQ2knPSbIuKePGIws6z1hgnPm3FJ6nx69Qi4BVgaEV/P+vxmlpd6YDblic5V+Xc2/mE3e3n06jke+AjwnKRFlXX/EhEP5RCLmWWqHif6/GWe+CPiScpf92ZmlgOP3DUzKxjPuWtmlonmvAPYxlf8ZmYF48RvZlYwTvxmZgXjxG9mVjBO/GZmBePEb2ZWME78ZmYF48RvZlYwTvxmZgXjxG9mVjBO/GZmBePEb2ZWME78ZmYF48RvZlYwTvxmZgXjxG9mVjBO/GZmBePEb2ZWME78ZmYF48RvZlYwTvxmZgXjxG9mVjBO/GZmBePEb2ZWME78ZmYF48RvZlYwTvxmZgXjxG9mVjBO/GZmBePEb2ZWME78ZmYFk0vil/T3kpZJ+p2kq/KIwcysqDJP/JL6A98ETgcOB86TdHjWcZiZFVUeV/xHA7+LiJcj4m3gTuDMHOIwMyukATmc80DgD22WVwHH7LyRpAagobK4XtKyDGJrNRR4PenGUoqRZHtuP+/sz90lft5V0WeedxXOP7K9lXkk/vaeRuyyImI2MDv9cHYlaX5ETMrj3Hny8y4WP+/iyqOpZxVwUJvl4cArOcRhZlZIeST+Z4BDJI2S9DfAucD9OcRhZlZImTf1RESLpP8L/BToD9waES9kHcdu5NLE1Av4eReLn3dBKWKX5nUzM6thHrlrZlYwTvxmZgXjxN9GEUtJSDpI0mOSlkp6QdJleceUJUn9JT0r6Sd5x5IVSftKulvSbyrv+9/lHVMWJH2m8hl/XtIPJA3KO6a8OPFXFLiURAvwzxExBjgWuKQgz7vVZcDSvIPI2I3AwxFxGDCeAjx/SQcClwKTImIs5Y4l5+YbVX6c+LcrZCmJiHg1IhZW7q+jnAQOzDeqbEgaDkwBbs47lqxI2hs4AbgFICLejog3cw0qOwOAPSUNAOoo8PghJ/7t2islUYgE2EpSCZgIPJVzKFn5BvA5YGvOcWTp3cBq4DuVJq6bJb0j76DSFhF/BG4AVgKvAmsj4mf5RpUfJ/7tEpWSqFWS9gJ+BHw6It7KO560SToDeC0iFuQdS8YGAEcCN0XERGADUPO/Z0naj/Jf8KOAA4B3SJqeb1T5ceLfrrClJCQNpJz0myLinrzjycjxwFRJyyk3650k6Xv5hpSJVcCqiGj9q+5uyl8Ete4U4PcRsToitgD3AMflHFNunPi3K2QpCUmi3N67NCK+nnc8WYmIqyNieESUKL/Xj0ZEzV8BRsSfgD9IGl1ZdTLwYo4hZWUlcKykuspn/mQK8KN2R/Koztkr9ZFSEmk4HvgI8JykRZV1/xIRD+UXkqXsU0BT5QLnZeDCnONJXUQ8JeluYCHlnmzPUuDSDS7ZYGZWMG7qMTMrGCd+M7OCceI3MysYJ34zs4Jx4jczKxgnfqspkkqSnu/gsZt7QwG6zmI0y4L78VthRMRFecdQDZIGRERL3nFY3+UrfqtFAyTdLmlJpe58HYCkZkmTKvfXS2qUtFjSPEn773wQSbMk3VrZ72VJl1bW73DFLulySbPanOPfJT1eqXV/lKR7JP1W0r8miPG9kuZKWiDpp5Le1ea4/yZpLuVS0mbd5sRvtWg0MDsixgFvAZ9sZ5t3APMiYjzwOPCPHRzrMOA0ymW7r6nUNdqdtyPiBOA/gfuAS4CxwMckDekoxsqx/wM4JyLeC9wKNLY57r4RcWJEfC1BDGYdcuK3WvSHiPhl5f73gPe1s83bQOusWwuAUgfHejAiNkfE68BrwC5/GbSjtcbTc8ALlTkPNlMuj9BaCLC9GEdT/oL4eaV8xucpFwtsdVeCc5vtltv4rRbtXIekvbokW2J7vZK/0vH/hc1t7rdu18KOF007T+HXus/Wnfbf2uY87cUoyl8UHU2FuKGD9WZd4it+q0Uj2swjex7wZJWP/z/A30oaImkP4IxuHKO9GJcBw1rXSxoo6X9XJWKzNpz4rRYtBS6QtAR4J3BTNQ9eqef+Zcozlf0E+E03DrNLjJUpP88BrpO0GFhEgWvGW3pcndPMrGB8xW9mVjBO/GZmBePEb2ZWME78ZmYF48RvZlYwTvxmZgXjxG9mVjD/H+PFsT6Z6lClAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Define the binning\n",
    "nbins = 10\n",
    "x = np.linspace(0, nbins - 1, nbins)\n",
    "\n",
    "# The background follows a linear shape\n",
    "b_yields = np.array([ (1 - i/2/nbins) for i in range(0, nbins) ])\n",
    "b_fracs = b_yields/np.sum(b_yields)\n",
    "\n",
    "# The signal shape is a peak\n",
    "s_fracs = np.zeros(nbins)\n",
    "s_fracs[4:7] = [ 0.1, 0.8, 0.1 ]\n",
    "\n",
    "# Define the signal and background\n",
    "s = 3\n",
    "b = 50\n",
    "s_and_b = s*s_fracs + b*b_fracs\n",
    "b_only = b*b_fracs\n",
    "\n",
    "# Now generate some data for the s+b model\n",
    "np.random.seed(14) # Set the random seed to make sure we generate the same data every time\n",
    "data = [ np.random.poisson(s*s_frac + b*b_frac) for s_frac, b_frac in zip(s_fracs, b_fracs) ]\n",
    "\n",
    "# Make a plot of the result\n",
    "plt.bar(x, s_and_b, color='orange', yerr=np.sqrt(s_and_b), label='s+b')\n",
    "plt.bar(x, b_only, color='b', label='b only')\n",
    "plt.scatter(x, data, zorder=10, color='k', label='data')\n",
    "plt.xlabel('bin number')\n",
    "plt.ylabel('Total expected yield')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot above shows the typical case of a binned analysis: we have 2 *templates*, i.e. reference shapes for signal and background. These are given by bin fractions $f_{S,i}$ and $f_{B,i}$ which are normalized to 1. They are then scaled by the yields $s$ and $b$.\n",
    "\n",
    "The ML estimate $\\hat{s}$ corresponds to the \"best fit\" of the templates to the data, so one can more or less read it off the plot: naively, one expects something a bit larger than the $s=3$ value that is shown in the plot, perhaps $\\hat{s} \\approx 5$ or so.\n",
    "\n",
    "We now check this by performing the ML estimation. A small change compared to the previous case is that instead of $L(s)$, we'll be using the quantity $\\lambda = -2 \\log L(s)$. The reason is as follows: first, the total PDF for all the bins together is a product of Poisson terms, one for each bin:\n",
    "$$\n",
    "L(s) = \\prod_{i=1}^{n_{\\text{bins}}} P(n_i, s f_{s,i} + b f_{b,i})\n",
    "$$\n",
    "Technically it's easier to deal with sums than products, and the log does this for us. As we'll see in a moment, this is also useful to deal with Gaussian PDFs, since in this case the log just extracts the squared term in the exponent. This is also the reason for the $-2$ term in the formula.\n",
    "The only consequence in practical terms is that instead of $L(s)$ we compute $\\lambda(s)$, and that we now want to minimize it instead of maximizing. But the ML estimate $\\hat{s}$ remains the same, since the value of $s$ that maximizes $L(s)$ also minimizes $\\lambda(s)$.\n",
    "\n",
    "Let's now compute the ML estimate (best fit) $\\hat{s}$ using $\\lambda(s)$ :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the s scan\n",
    "s_values = np.arange(-2, 15, 0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathbf{\\color{red}{\\rightarrow\\text{Compute }\\lambda\\text{ for each value of }s \\text{, plot the result and find the best-fit value among the values tested}}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "hat{s}, ML = (-2.0, -9990)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZoAAAEGCAYAAABcolNbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAATAElEQVR4nO3de6zk5X3f8ffHXofGNmAuS7ks7mKB3UKdNGaKSeO2KKCF2I3XdkBZt5K3MtLWlqPGcpsWhBVcY6mluRDhxFjUxAUn8kW0lHUxJYuJS1URzFnEZZfrEpNwutReexHgWrK9+Ns/5jl4OMzZPZzdZ2YW3i9pdH7zPL/nN9/fc2bnM78Lh1QVkiT18qppFyBJenkzaCRJXRk0kqSuDBpJUlcGjSSpq1XTLmDWHH300bV27dpplyFJB5WtW7d+t6pWj+szaBZZu3Ytc3Nz0y5Dkg4qSf5qqT5PnUmSujJoJEldGTSSpK4MGklSVwaNJKkrg0aS1JVBI0nqyqCRJHVl0EiSujJoJEldGTSSpK4MGklSVwaNJKkrg0aS1JVBI0nqyqCRJHVl0EiSujJoJEldGTSSpK4MGklSVwaNJKkrg0aS1JVBI0nqyqCRJHVl0EiSujJoJEldGTSSpK5mLmiS/HySO5Lcn+SrSQ4b6fu51re99f+N1n56e74jyZVJ0toPSfLl1n5nkrVT2i1JesWauaABPgdcVFVvBW4AfgsgySrgT4APVdVpwFnAj9uYq4BNwCntcV5rvxB4qqpOBq4ALp/QPkiSmlkMmrcAt7flLcCvteV1wH1VdS9AVX2vqp5LchxwWFXdUVUFXAe8p41ZD1zblq8Hzl442pEkTcYsBs024N1t+QLgxLb8ZqCS3JLk7iT/prWfAMyPjJ9vbQt9TwBU1R7gaeCoxS+YZFOSuSRzu3btOqA7I0mvdKum8aJJbgWOHdN1CfBB4Mokvw1sBn7U+lYB7wD+PvAD4OtJtgLPjNlOLbzUXvp+2lB1NXA1wGAweFG/JGnlphI0VXXOPlZZB5DkzcC7Wts88D+r6rut72vA2xhet1kzMnYNsHNkzInAfLvGcziw+0DsgyRpeWbu1FmSY9rPVwEfBz7bum4Bfi7Ja1to/GPggap6Eng2yZnt+ssHgBvbmM3AxrZ8PnBbu44jSZqQmQsa4P1JHgEeYnhk8nmAqnoK+H3gLuAe4O6quqmN+TDDu9V2AI8BN7f2a4CjkuwAPgZcNKF9kCQ18Qv+Cw0Gg5qbm5t2GZJ0UEmytaoG4/pm8YhGkvQyYtBIkroyaCRJXRk0kqSuDBpJUlcGjSSpK4NGktSVQSNJ6sqgkSR1ZdBIkroyaCRJXRk0kqSuDBpJUlcGjSSpK4NGktSVQSNJ6sqgkSR1ZdBIkroyaCRJXRk0kqSuDBpJUlcGjSSpK4NGktSVQSNJ6sqgkSR1ZdBIkroyaCRJXRk0kqSuDBpJUlcGjSSpq5kLmiQ/n+SOJPcn+WqSw1r7a5Jc29ofTHLxyJjTW/uOJFcmSWs/JMmXW/udSdZOabck6RVr5oIG+BxwUVW9FbgB+K3WfgFwSGs/HfgXI8FxFbAJOKU9zmvtFwJPVdXJwBXA5RPZA0nS82YxaN4C3N6WtwC/1pYLeF2SVcDPAj8CnklyHHBYVd1RVQVcB7ynjVkPXNuWrwfOXjjakSRNxiwGzTbg3W35AuDEtnw98P+AJ4G/Bn63qnYDJwDzI+PnWxvt5xMAVbUHeBo4avELJtmUZC7J3K5duw7s3kjSK9xUgibJrUm2jXmsBz4IfCTJVuBQhkcuAGcAzwHHAycB/yrJm4BxRyi18FJ76ftpQ9XVVTWoqsHq1av3c+8kSaNWTeNFq+qcfayyDiDJm4F3tbZ/CvyPqvox8J0k/xsYAP8LWDMydg2wsy3PMzwimm+n3A4Hdh+QnZAkLcvMnTpLckz7+Srg48BnW9dfA7+codcBZwIPVdWTwLNJzmzXXz4A3NjGbAY2tuXzgdvadRxJ0oTMXNAA70/yCPAQwyOTz7f2PwJez/Aazl3A56vqvtb3YYZ3q+0AHgNubu3XAEcl2QF8DLhoInsgSXpe/IL/QoPBoObm5qZdhiQdVJJsrarBuL5ZPKKRJL2MGDSSpK4MGklSVwaNJKkrg0aS1JVBI0nqyqCRJHVl0EiSujJoJEldGTSSpK4MGklSVwaNJKkrg0aS1JVBI0nqyqCRJHVl0EiSujJoJEldGTSSpK4MGklSVwaNJKkrg0aS1JVBI0nqyqCRJHVl0EiSujJoJEldGTSSpK4MGklSVwaNJKkrg0aS1JVBI0nqyqCRJHW1ahovmuQC4BPA3wHOqKq5kb6LgQuB54B/WVW3tPbTgf8M/CzwNeA3q6qSHAJcB5wOfA/49ap6vI3ZCHy8bfpTVXVtz/36d1/dzgM7n+n5EpLUzanHH8alv3raAd/uio9oknx0P153G/A+4PZF2zwV2ACcBpwHfCbJq1v3VcAm4JT2OK+1Xwg8VVUnA1cAl7dtHQlcCrwdOAO4NMkR+1GzJGkF9ueI5mPAH6xkYFU9CJBkcdd64EtV9UPgW0l2AGckeRw4rKruaOOuA94D3NzGfKKNvx74www3fC6wpap2tzFbGIbTF1dS83L0+CYgSQe7/blG86KUOABOAJ4YeT7f2k5oy4vbXzCmqvYATwNH7WVbL5JkU5K5JHO7du06ALshSVqwP0c0tbfOJLcCx47puqSqblxq2BKvs1T7Sse8sLHqauBqgMFgsNf9kiS9NHsNmiTPMv7DOQwvyi+pqs5ZQT3zwIkjz9cAO1v7mjHto2Pmk6wCDgd2t/azFo35xgpqkiTth72eOquqQ6vqsDGPQ6uqxx1rm4ENSQ5JchLDi/7frKongWeTnNmuv3wAuHFkzMa2fD5wW1UVcAuwLskR7SaAda1NkjRB07q9+b3Ap4HVwE1J7qmqc6tqe5KvAA8Ae4CPVNVzbdiH+entzTe3B8A1wBfajQO7Gd61RlXtTnIZcFdb75MLNwZIkiYnwy//WjAYDGpubm7fK0qSnpdka1UNxvX5lwEkSV0ZNJKkrgwaSVJXBo0kqSuDRpLUlUEjSerKoJEkdWXQSJK6MmgkSV0ZNJKkrgwaSVJXBo0kqSuDRpLUlUEjSerKoJEkdWXQSJK6MmgkSV0ZNJKkrgwaSVJXBo0kqSuDRpLUlUEjSerKoJEkdWXQSJK6MmgkSV0ZNJKkrgwaSVJXBo0kqSuDRpLUlUEjSepqKkGT5IIk25P8JMlgUd/FSXYkeTjJua3ttUluSvJQG/cfRtY/JMmX25g7k6wd6duY5NH22DixHZQkPW9aRzTbgPcBt482JjkV2ACcBpwHfCbJq1v371bV3wZ+AfilJL/S2i8Enqqqk4ErgMvbto4ELgXeDpwBXJrkiK57JUl6kakETVU9WFUPj+laD3ypqn5YVd8CdgBnVNUPqurP29gfAXcDa0bGXNuWrwfOThLgXGBLVe2uqqeALQzDS5I0QbN2jeYE4ImR5/Ot7XlJ3gD8KvD1xWOqag/wNHDUcrYlSepvVa8NJ7kVOHZM1yVVdeNSw8a01cg2VwFfBK6sqr/cx5i9bmtRrZuATQBvfOMblyhNkrQS3YKmqs5ZwbB54MSR52uAnSPPrwYerao/GDNmvgXR4cDu1n7Wom19Y4lar27bZjAYjA0jSdLKzNqps83AhnYn2UnAKcA3AZJ8imGIfHTMmIU7ys4HbquqAm4B1iU5ot0EsK61SZImqNsRzd4keS/waWA1cFOSe6rq3KranuQrwAPAHuAjVfVckjXAJcBDwN3Da/38YVV9DrgG+EKSHQyPZDYAVNXuJJcBd7WX/WRV7Z7gbkqSgAy//GvBYDCoubm5aZchSQeVJFurajCub9ZOnUmSXmYMGklSVwaNJKkrg0aS1JVBI0nqyqCRJHVl0EiSujJoJEldGTSSpK4MGklSVwaNJKkrg0aS1JVBI0nqyqCRJHVl0EiSujJoJEldGTSSpK4MGklSVwaNJKkrg0aS1JVBI0nqyqCRJHVl0EiSujJoJEldGTSSpK4MGklSVwaNJKkrg0aS1JVBI0nqyqCRJHVl0EiSuppK0CS5IMn2JD9JMljUd3GSHUkeTnLumLGbk2wbeX5Iki+3MXcmWTvStzHJo+2xsetOSZLGmtYRzTbgfcDto41JTgU2AKcB5wGfSfLqkf73Ad9ftK0Lgaeq6mTgCuDytu6RwKXA24EzgEuTHNFlbyRJS5pK0FTVg1X18Jiu9cCXquqHVfUtYAfDkCDJ64GPAZ8aM+batnw9cHaSAOcCW6pqd1U9BWxhGF6SpAmatWs0JwBPjDyfb20AlwG/B/xgqTFVtQd4GjhqH9t6gSSbkswlmdu1a9f+7oMkaUS3oElya5JtYx7r9zZsTFsl+XvAyVV1w3LH7KX9xY1VV1fVoKoGq1ev3kt5kqSXalWvDVfVOSsYNg+cOPJ8DbAT+EXg9CSPM6z5mCTfqKqzRsbMJ1kFHA7sbu1nLdrWN1ZQkyRpP8zaqbPNwIZ2J9lJwCnAN6vqqqo6vqrWAu8AHmkhszBm4Y6y84HbqqqAW4B1SY5oNwGsa22SpAnqdkSzN0neC3waWA3clOSeqjq3qrYn+QrwALAH+EhVPbePzV0DfCHJDoZHMhsAqmp3ksuAu9p6n6yq3T32R5K0tAy//GvBYDCoubm5aZchSQeVJFurajCub9ZOnUmSXmYMGklSVwaNJKkrg0aS1JVBI0nqyqCRJHVl0EiSujJoJEldGTSSpK4MGklSVwaNJKkrg0aS1JVBI0nqyqCRJHVl0EiSujJoJEldGTSSpK4MGklSVwaNJKkrg0aS1JVBI0nqyqCRJHWVqpp2DTMlyS7gr5boPhr47gTLORCseTKseTKseTJWUvPfqqrV4zoMmpcgyVxVDaZdx0thzZNhzZNhzZNxoGv21JkkqSuDRpLUlUHz0lw97QJWwJonw5onw5on44DW7DUaSVJXHtFIkroyaCRJXRk0e5Hkd5I8lOS+JDckecMS6z2e5P4k9ySZm3CZJDkvycNJdiS5aEx/klzZ+u9L8rZJ17ionhOT/HmSB5NsT/KbY9Y5K8nTbU7vSfLb06h1sX39rmdwrt8yMof3JHkmyUcXrTP1uU7yx0m+k2TbSNuRSbYkebT9PGKJsXt9/0+45pn+zFii5k8k+T8jv/93LjF25fNcVT6WeADrgFVt+XLg8iXWexw4eko1vhp4DHgT8DPAvcCpi9Z5J3AzEOBM4M4pz+txwNva8qHAI2NqPgv479N+D7zU3/WszfWY98r/Zfgf1s3UXAP/CHgbsG2k7T8CF7Xli8b9+1vO+3/CNc/0Z8YSNX8C+NfLeO+seJ49otmLqvqzqtrTnv4FsGaa9SzhDGBHVf1lVf0I+BKwftE664HraugvgDckOW7ShS6oqier6u62/CzwIHDCtOo5wGZqrhc5G3isqpb6yxdTU1W3A7sXNa8Hrm3L1wLvGTN0Oe//LsbVPOufGUvM83Ls1zwbNMv3QYbfVMcp4M+SbE2yaYI1wfAD+omR5/O8+EN7OetMRZK1wC8Ad47p/sUk9ya5Oclpk61sSfv6Xc/sXAMbgC8u0TeLc/03q+pJGH45AY4Zs84sz/esfmaM8xvtdN8fL3GKcr/medX+VnewS3IrcOyYrkuq6sa2ziXAHuBPl9jML1XVziTHAFuSPNS+OUxCxrQtvmd9OetMXJLXA/8F+GhVPbOo+26Gp3i+384Z/zfglAmXOM6+ftezOtc/A7wbuHhM96zO9XLM6nzP8mfGYlcBlzGct8uA32MYkqP2a55f8Uc0VXVOVf3dMY+FkNkI/BPgn1U7WTlmGzvbz+8ANzA8zJyUeeDEkedrgJ0rWGeikryGYcj8aVX918X9VfVMVX2/LX8NeE2Soydc5oss43c9c3Pd/Apwd1V9e3HHrM418O2F047t53fGrDNz830QfGYsruXbVfVcVf0E+E9L1LJf8/yKD5q9SXIe8G+Bd1fVD5ZY53VJDl1YZngxcNu4dTu5CzglyUntW+sGYPOidTYDH2h3RJ0JPL1wSmIakgS4Bniwqn5/iXWObeuR5AyG79XvTa7KsTUt53c9U3M94v0scdpsFue62QxsbMsbgRvHrLOc9//EHCSfGYvrGb2G+N4latm/eZ70XQ8H0wPYwfC85D3t8dnWfjzwtbb8JoZ3YNwLbGd4ym3Sdb6T4Z1bjy28PvAh4ENtOcAftf77gcGU5/UdDA+77xuZ23cuqvk32nzey/Ci6j+YgffD2N/1LM91q+m1DIPj8JG2mZprhiH4JPBjht+eLwSOAr4OPNp+HtnWff7fX3v+ovf/FGue6c+MJWr+Qnuv3scwPI470PPsn6CRJHXlqTNJUlcGjSSpK4NGktSVQSNJ6sqgkSR1ZdBIkroyaCRJXRk00kGg/dfkN7U/fLktya9PuyZpuV7xf1RTOkicB+ysqncBJDl8yvVIy+YRjXRwuB84J8nlSf5hVT097YKk5TJopINAVT0CnM4wcP79NP53y9JKeepMOggkOR7YXVV/kuT7wD+fcknSshk00sHhrcDvJPkJw7+8++Ep1yMtm3+9WZLUlddoJEldGTSSpK4MGklSVwaNJKkrg0aS1JVBI0nqyqCRJHX1/wFpUj4ZwbsZWwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def lambda_s(s, data) :\n",
    "    return sum( [ -999 for n, s_frac, b_frac in zip(data, s_fracs, b_fracs) ] ) # replace -999 with the correct expression\n",
    "\n",
    "lambda_values = np.array([ lambda_s(s_value, data) for s_value in s_values ]) # compute the likelihood (the Poisson PDF for each s)\n",
    "\n",
    "# Plot the result\n",
    "plt.plot(s_values, lambda_values)\n",
    "plt.xlabel('s')\n",
    "plt.ylabel('lambda(s)')\n",
    "\n",
    "# Find the maximum likelihood value in the list : merge the 2 lists and find the max based on the L value\n",
    "from operator import itemgetter\n",
    "print('hat{s}, ML =', min( [ (s,l) for s,l in zip(s_values, lambda_values) ], key=itemgetter(1)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If all goes well, you should find a best-fit $\\hat{s}$ of about 5. We can again compute the exact best-fit value using `minimize_scalar`:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathbf{\\color{red}{\\rightarrow\\text{Use the minimize_scalar minimization routine to find the precise value of } \\hat{s}}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import minimize_scalar\n",
    "# ==> Use the minimize_scalar function to find s_hat\n",
    "#     The syntax is minimize_scalar(func, (min, max)).x  (the .x returns the position of the minimum)\n",
    "#     and 'func' should be the lambda_s function you defined above\n",
    "s_hat = -999 # replace with real code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should find a best-fit value of about $\\hat{s} = 5.06$. \n",
    "\n",
    "Note than since $L(s)$ was defined only up to a multiplicative factor, $\\lambda(s)$ is defined up to an additive constant, so the value $\\lambda(\\hat{s})$ at the minimum isn't meaningful. On the other hand, the difference between values of $\\lambda(s)$ at different $s$ is meaningful, and will be used later in hypothesis testing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Hypothesis testing basics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far we've focused on estimating parameters, but to get realistic results (including setting uncertainties on the best-fit values), we need an additional ingredient, *hypothesis testing*. The procedure for this is as follows:\n",
    "* **Define the hypotheses to test**. We need two of them: a baseline *null* hypothesis, and an *alternate* which is tested against it\n",
    "* **Define the *test* to perform** -- meaning the observables to use, and procedure to decide which hypothesis to accept.\n",
    "* **Compute the result of the test for the observed data**, in particular the *p-value* of the test, and decide whether to accept the null or the alternate hypothesis.\n",
    "\n",
    "To start with, we consider a simple Poisson counting process\n",
    "* Define two different signal hypotheses ($n=0$ and $n=5$)\n",
    "* Define the test in terms of the observed number of events: if it is less than some threshold, accept the null $n=0$, otherwise $n=5$.\n",
    "\n",
    "For the study below, we consider a single-bin counting experiment. In this case the disciminant is just the observed count $n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the  hypotheses\n",
    "b = 10 # background of 10 events\n",
    "s_null = 0 # the null hypothesis: no signal, only background\n",
    "s_alt  = 5 # the alternate hypothesis: \n",
    "\n",
    "# consider various values of the observed n (our discriminant)\n",
    "ns = np.arange(1, 30, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathbf{\\color{red}{\\rightarrow\\text{Compute and plot the PDFs for the two hypotheses. These are Poisson PDFs with parameters s_null + b and s_alt + b}}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note that the Poisson PDF is implemented in scipy as `scipy.stats.poisson.pmf`\n",
    "pass # real code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As for any test, there are two ways to go wrong:\n",
    "* Rejecting the null although it is true (false positive): the fraction of such cases is the *Type-I error rate* or *p-value*.\n",
    "* Accepting the null although it is false (false negative): the fraction of such cases is the *Type-II error rate*, also (1 - Power)\n",
    "\n",
    "A more stringent test reduces the Type-I error rate, but at the expense of an increased Type-II rate and vice versa. This is illustrated by the *ROC curve* which shows each error rate as a function of the other. The threshold for the test should be chosen depending on the desired balance between these error rates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAACQCAYAAAAWRCxZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAn6klEQVR4nO2deZhU1bHAfzULDMg+IBhQh6CJLA8GYQBFZVQeakQRF5ZIVIwgKhj3HSUxJjwxGjdwZXGDKETBLUGNqEQWQRAREJB1RHbZt1nO+6NuM80wzPQM3XP79tTv+87Xt+/St07f7lv3VNWpEucchmEYhnEkkvwWwDAMw4hvTFEYhmEYJWKKwjAMwygRUxSGYRhGiZiiMAzDMErEFIVhGIZRIil+CxBN6tev7zIyMvwWwzAMIzDMnTt3s3OuQUn7JJSiyMjIYM6cOX6LYRiGERhEZHVp+5jpyTAMwyiRhBpRGHHIN99ATk5k+6amQpcuULVqbGUyDKNMmKIwYoNz8MgjMHRo2Y479VR4801o1iw2chmGUWZMURjRZ88e+P3vYcIE6NcPhgwBkdKPW7oUBg9WZTF6NFx2WexlNRKG3NxccnJy2Ldvn9+ixCVpaWk0adKE1NTUMh9risKILj/+CD16wNdfw/DhcNddkSkJgKws6NwZeveGyy9XBTNihJmijIjIycmhZs2aZGRkIJH+5ioJzjm2bNlCTk4OTZs2LfPx5sw2osfs2Xqz//57mDwZ7r47ciURIiMDvvgCbr0Vnn5aFceKFTER10gs9u3bR3p6uimJYhAR0tPTyz3aMkVhRIc33oCzzoK0NJgxAy66qPyfVaUKPP44vP02/PCDmqL++c/oyWokLKYkjszRfDemKIyjo6AA7r0XrrwSOnXSUUWrVtH57EsuURPWr36l/oo//AH274/OZxtGHDBt2jS6d+8OwNixYxk8eLDPEhWPKQqj/OzcCT17qi9i4ECYOhXq14/uOZo2henTVUk89RSceSasXBndcxiGUSKmKIzysXIlnH46vP+++hKee05NRrGgShX4+9/V/LR0KbRtq2Ypw4gzVq1aRfPmzRkwYAAtW7akW7du7N27l+zs7INZIzZv3kzQUg1Z1JNRdr74Ai69FPLy4F//gq5dK+a8PXtCZib06qXnf/VVDb81jKLccgvMnx/dz8zM1AeWUli2bBnjx4/nxRdfpFevXkyaNCm6cviAjSiMsrF6NVx4IaSnw6xZFackQoRMUV26wIABMG9exZ7fMEqhadOmZGZmAtCuXTtWrVrlqzzRwEYURuQ4pxPpnNORhF/D56pVdfZ2u3Y6ypgzJ/q+ESPYRPDkHyuqhs37SU5OZu/evaSkpFBQUAAQyAmBNqIwIuf55+GTT+Cxx/xTEiGOPVZ9FuvX6wS9vDx/5TGMEsjIyGDu3LkATJw40Wdpyo4pCiMyVq6EO+5QU9PAgX5Lo2RlqRP9P/+Be+7xWxrDOCJ33HEHo0aN4vTTT2fz5s1+i1NmxDnntwxRo3379s7qUcSAggI491yYOxcWLoQTTvBbokMZMgSeeQZefx1++1u/pTF8YvHixTRv3txvMeKa4r4jEZnrnGtf0nE2ojBKZ+RImDYNnngi/pQE6CzuM8+E666LfqSLYRixVRQicr6IfC8iy0XkMNuAiJwiIjNEZL+I3FGWY40KYvlyzdl0wQVw7bV+S1M8qanw1ltQr546t7dsqZDTFhRE3hJo4G5UQmKmKEQkGXgWuABoAfQVkRZFdtsK3Aw8Vo5jjVhTUAD9++uN+MUXy57gryJp2FCd2+vWQZ8+MXNu5+bCa69p+qnk5Mhb8+Y6sXz79piIZRgxJZbhsR2A5c65FQAiMgHoASwK7eCc2whsFJELy3qsUQE89ZTOWRg3Dho39lua0unQAUaN0hDe++6DRx+N2kdv2wYvvKBfyY8/QosWWpMpJYJ/UEEB/PvfmoXk/vvhqqvgppv0MwwjCMRSUTQG1oa9zwE6RvtYERkIDAQ4IR7t50Fl6VJN9nfRRfC73/ktTeRce6063UeM0Mf+Pn2O6uNWr4Ynn9QB1a5dcM45unzeeZBUhvH4sGE63ePZZ+Hll9Xtc8456ofv3j0yhWMYfhFLH0VxdopILbURH+uce8E51945175BgwYRC2eUQH4+XHMNVKumcyfi2eRUHE88oXUsrr1Wa3aXgzlzVMc0a6aprEK1mD75RN01ZVESIdq3hzFjYO1a+OtfYdkydak0a6Z5FQMYNWlUEmKpKHKA48PeNwHWVcCxxtHyxBNaU+KZZ+C44/yWpuxUqQITJ0Ldunon3ro1osMKCuDddzU7SFYWfPih1k9asUL9Em3bRke8Bg102seKFZrb8OSTdfDWpIm6hJYujc55jPghIyODzZs3s23bNkaOHOm3OGUmloriK+BkEWkqIlWAPsCUCjjWOBoWL4YHHtAbbN++fktTfho1gkmT1KHQt6+Okkpgxw41AV18MaxapRG3a9eqBev440s8tNykpGjJjY8/hu++U9fKW2+pxezVV2NzTsNfgqoocM7FrAG/AZYCPwD3e+sGAYO85Ubo6GEHsM1brnWkY0tr7dq1c8ZRkJvrXIcOzqWnO7d+vd/SRIcXXnAOnLvuuiPusmKFcy1bOpeS4txTT+nX4Bc5Oc6ddZaK3L+/c7t3+ydL0Fi0aJHfIjjnnOvRo4c79dRTXYsWLdzzzz/vnHPuxBNPdJs2bXK9e/d2aWlprk2bNu6OO+6ocNmK+46AOa6Ue6vNzDYKGT5cbSD/+Iem8k4EnNOIrZ9+0lQfZ599yOYvv9Sn+txctVade64/YoaTlwd//CM88oiG1b71lkVIRUL4rGMfs4yzdetW6tWrx969e8nKyuKzzz6jXbt2zJkzh127dtG9e3cWLlwYXeEixGZmG0fHwoXw0ENwxRWJoyRAHfFt2ujyxRdDWObO119XvVG7NsycGR9KAtQk9fDDWjBw82Z1go8d67dURqQ89dRTtGnThk6dOrF27VqWLVvmt0hHjQXlGfrUPWCA3jGffdZvaaJPKO3zrl3QsycF73/IsGF6M+7SRV0Z6em+SlgsXbvqU/GVV6qT+9NP9fLUqOG3ZPGPX1nGp02bxscff8yMGTOoXr062dnZgUwrXhQbURhqapo5E/7v/zQkJ4HZ86/P6NN5DQ8/rNGzU6fGp5IIcdxx8NFHOg/j1Vc1Guvbb/2WyjgS27dvp27dulSvXp0lS5Ywc+bMQ7bXrFmTnTt3+iRd+TFFUdnZt09jNTMz4eqr/ZYmpvxEI7KZxsSZTRjxp7289FLsynxHk+RktQp+8onOEO/QAV56yfJHxSPnn38+eXl5tG7dmqFDh9KpU6dDtqenp9O5c2datWrFnXfe6ZOUZcdMT5WdJ5/U6cejR5dvFllAmE8bLuJdtlKPt+lJj8k/wtBgBT6cfbaaon73O7UUfvqpluOoWdNvyYwQVatW5cMPPzxsfXg51DfeeKMCJYoOiXtnMEpn40b4y180Tcc55/gtTcyY/FMHzmA6ANM5gx5M0TQfTz7ps2Rlp2FDrUL75z/DhAl62WxGtxFrTFFUZoYNg927o5o8L94YMwZ6zr6H5ixmNh1oy/zCjbffrrPqAkZSkiYXnDxZg9XOOkvnFRpGrDBFUVlZvFjTod5wA5xyit/SxIR33tFaRimSz808yXGsP3SH/HzIzg6ssb97dx1d5OTAGWdo6RDDiAWmKCord96pcZYPPeS3JDHh8881qV9KCuS5JIrPM4kmXLr33gqVLZp06aK+ip07VVksWOC3REYiYoqiMvLxx/D++2q/qF/fb2mizjff6Nw60ER/xyQfKPmARx8NdMxpu3bwxReqFLt00XyOhhFNTFFUNvLz1TbftKkWQ0gwVqyA88+HvXv1fUSF7pzTadmlJA6MZ5o31xpT9evrRL2PPvJbIiORMEVR2Rg7Vu0Tw4dDWprf0kSVDRugWzctmS0C+/eX4eBNmzR9a4DJyNCRxUknwYUX6oxzI/iMHTuWwYMH+yqDKYrKxK5dmkL8tNM0p1MCsX27jiRWr9aooDIpiRDjxqnBP8A0agSffaYzuHv10ukxlZZGjfSJIVqtUSO/e+QbpigqEyNGwPr1WmwhaFXrSmDfPs0Au2CBzmIul5IIcfHFhXargFKnjqYm6dpVB0mPP+63RD6xYUOFf96qVas45ZRTuPrqq2ndujWXX345e/bsObi9oKCAjIwMtm3bdnDdSSedxIYNG3j33Xfp2LEjbdu2pWvXrmwo5nzXXHMNEydOPPi+RljirxEjRpCVlUXr1q15KMpBKiUqChGZGrYc3NAQQ2MoR4yA3r2hSFqBIJOfD7/9LUybps7co1ISoKOuCy+Mhmi+cswxMGWKDhxvvx2GDg1sFHDg+P777xk4cCALFiygVq1ahxQqSkpKokePHrz99tsAzJo1i4yMDBo2bMgZZ5zBzJkzmTdvHn369OHRMsxvmjp1KsuWLWP27NnMnz+fuXPn8vnnn0etT6WNKMIzxCWWraKy8cADGgI0fLjfkkQN53QayNtvQ2oqHCgluCliPv0URo2K0of5R9WqMH68jir+/Ge47TZTFhXB8ccfT+fOnQHo168f06dPP2R77969+cc//gHAhAkT6N27NwA5OTmcd955/M///A8jRozgu+++i/icU6dOZerUqbRt25ZTTz2VJUuWRDW9eWmKwn5WicDXX8Mrr8Af/qAezwThwQfhxRd1JJGbG+UPHzxYQ6gCTnKyfkc336ypt2++2ZRFrJEiZt3t27eTmZlJZmYmU6ZM4bTTTmP58uVs2rSJd955h0svvRSAIUOGMHjwYL799luef/75YtOTp6SkUFBQAGh10gPe05FzjnvvvZf58+czf/58li9fzu+jGJxRmqL4pYhMEZF3w5YPtqhJYcQO59T2kJ4O993ntzRR46mn9Ck5JSXCENiyUlAAZ54Z6JDZECKqJG6/HZ55Bm68UbtnxIY1a9Yww5vMMn78eLp3737wBn7xxRcjIvTs2ZPbbruN5s2bk+7lud++fTuNGzcGYNy4ccV+dkZGBnPnzgVg8uTJ5HpPSOeddx6jR49m165dAPz4449s3Lgxan0qLXtsj7Dlx6J2VqPiePddNeA/+6wWJkoAJk7UwVFqagxGEuGsW6cVg155JYYnqRhE1EWVmqrWx7w8eP75hE4Y7BvNmzdn3LhxXH/99Zx88snccMMNh+3Tu3dvsrKyGBtWunDYsGFcccUVNG7cmE6dOrFy5crDjhswYAA9evSgQ4cOnHvuuRxzzDEAdOvWjcWLF3PaaacB6uR+7bXXOPbYY6PTqdKKaoca6q9oEOn+frR27dqVr+J4onLggHO/+pVzp5ziXG6u39JEhRkznEtLc65KFed0uFR6q5G8x71Cv8gPKNomT/a721GjoMC5Bx7Qbl19tXN5eX5LFD0WLVp06IqGDct/zYtrDRuWKsPKlStdy5YtY9TDo+ew78g5B8xxpdxbSxxRiBrbHgSGoMlykkQkD3jaOfen6KgqI2aMGgVLl+qoIiX4pUdWroQePdRsUqGmk169dHRRr14FnjQ2iGgJ2NRUTfOVl6dzMBPg53E469eXvo8REaUNPG8BzgCynHPpzrm6QEegs4jcGmvhjKNg40b19nbtmhDhntu2aTd+/llvdjHxSxyJ/fs1iVICeYEffFB9PK+/roWQKvT7TGAyMjJYuHCh32JEndIUxVVAX+fcQWOZc24F0M/bZsQrd98Ne/bA008HfnJdbi5cfjl8/30UJtSVl4ULdTJCAnH//VomfcIEzbQbU3+PEWhKUxSpzrnD6mc55zYBqbERyThqvvxS7Qm33Rb4WhPOaZTOJ5+oeaSYiMGK4y9/ga++8lGA6HPXXTpze9IktbBFbS6KkVCUpihK+tnYTyoeycuDm26CJk10kl3AGTECXnpJlYTvNzHn1JTnq7aKPrfequHG77yjIzdfRmxGXFOaC6uNiOygsOpLyEgrQGKlHk0UnnsO5s+HN9/UwkQBZuJEtaDFbK5EedixQ2uMJ1ge7yFD9Hu+8Ubo2VNHGNWq+S2VES+UOKJwziU752o552p6rVbYezM9xRsbN+ooomtXfTQMMLNmqZO1SpU4UhIhPv5YJyEkGDfcoNVx//UvzcS7fbvfEgWLLVu2HJyB3ahRIxo3bnzw/YEoD4fDkwNmZ2czZ86cqH5+UUpLCpgmIreIyDMiMlBEyhREJyLni8j3IrJcRO4pZruIyFPe9gUicmrYtlUi8q2IzBeR2H4LiUKCOLBXrdIkrs7F8Qzim27SeN0EY8AAeOMNdXNlZ0c/AWtFUtFZxtPT0w/OwB40aBC33nrrwfdVqlSpmE7HiNJ8FOOA9sC3wG+Av0X6wSKSDDwLXAC0APqKSIsiu10AnOy1gUDRTGxnO+cynXPtIz1vpSVBHNjbtxeGwUIcjiZC5OdrkeoEDBXq00en3ixdql0Mqj70Icv4Iezdu5emTZseTLOxY8cOMjIyyM3NJTs7m1tuuYXTTz+dVq1aMXv2bAB2797NtddeS1ZWFm3btmXy5MnR7UQ5KU1RtHDO9XPOPQ9cDpxZhs/uACx3zq1wzh0AJnBoShC89694EwRnAnVE5LgynMOAhHFg5+ZqWuwlS46i+FBFsm6dllBNoPkVIc4/Xy1sW7ZA584aHWyUjWrVqpGdnc37778PaKbYyy67jNRUtdrv3r2bL7/8kpEjR3LttdcC8Mgjj3DOOefw1Vdf8emnn3LnnXeye/du3/oQojRFcfBxyTlX1me7xsDasPc53rpI93HAVBGZKyIDj3QSzyQ2R0TmbNq0qYwiJgghB/bjjwfWge2cJmz96CMf50qUhy++CHwJ1SNx2mnaPRHNj/jll35LFDyuu+46xowZA8CYMWPo37//wW19+/YF4KyzzmLHjh1s27aNqVOnMnz4cDIzM8nOzmbfvn2sWbPGF9nDiTTqCTTSqVpYFJRzztUq4djijORFH71K2qezc26diBwLfCQiS5xzh1XicM69ALwA0L59+8R7tCuNBHFgDxumjtSYpAyPNWPGwK9+Bfcc5oYLPC1bwn//C//7v/oTmzQJLrjAb6mCQ+fOnVm1ahWfffYZ+fn5tGrV6uC2ounIRQTnHJMmTeLXv/71Idv69+/PvHnz+MUvfsEHH3xQIbKHE2nUUyjSKSVsuSQlATo6OD7sfRNgXaT7OOdCrxuBt1FTllGUBHBg//3v8Kc/xVkYbFm5916N501AMjJg+nR1fV18sRZDMiLnqquuom/fvoeMJoCDxYumT59O7dq1qV27Nueddx5PP/10KBEr8+bNA3Q0Mn/+fF+UBMS2ZvZXwMki0lREqgB9gKI1LKYAV3nRT52A7c65n0TkGBGpCSAixwDdALOSFiUBHNhjxuiEr0AriRC9eyfczO0QDRtq4b/OneHKK7WuhREZV155JT///PNBU1OIunXrcvrppzNo0CBefvllAIYOHUpubi6tW7emVatWDI2XtDGlpZc9moZGSi0FfgDu99YNAgZ5y4JGRv2ARla199b/EvjGa9+Fji2tVao047m5zmVmOtekiXM7d/otTbmYNMm5pCTnUlOjmw066mnGy9KqVnVuzRq/v9qYsXevcz16aFcfekjTlscLRVNo+5BlvFjeeust169fv0PWdenSxX311Vfl7Wq5iUma8SgooQ+AD4qsey5s2QE3FXPcCqBNLGULPAGfgf3xx9C3b0B9EiWxfz+0basxpTVr+i1N1ElLUwvbwIHwxz/C5s1qOozHNOXxkGV8yJAhfPjhh76ZjKJFHF5eo1Q2bAi0A3vmTLjkEl3Oz0/A6NItWyArS2NK4/EOepSkpMDLL0P9+pqLa8kS9Vs0aOC3ZPHH008/Xez6adOmVawgR4kVQgwi99wTWAf2t99q1ExeniqIBChJXTzff68dTTgtqIjAo4+qj2n6dGjfHmKcRcLwEVMUQeOLLwLrwF6+HLp1Ux3nXIKZnIrj4481y14Cc801Gj4LOot79GhfxTkYLWQcztF8N6YogsT69WrYz8gI3AzsH39US9mWLfo06nvK8IriuefgiSf8liKmtGsHc+eqovj972HQIH8mTKalpbFlyxZTFsXgnGPLli2kpZUv6XfiGVATlQMHNL/F1q0wY0agHNibN+uErXXrApKaI9rcdhs0a6aTEBKU+vU16+wDD2jVvG++Uad346K5GGJIkyZNyMnJodJmaCiFtLQ0mjRpUq5jTVEEhVtvVWPw+PHQJjgBYTt2qKl+6VJ1glY6JRHikktg3DjNnZ6gpKTA8OHqx7/mGjj1VHjrLTjrrIo5f2pqKk2bNq2Yk1UyzPQUBEaPhpEj4Y47NLVnQNi+XR+iv/66kiuJEFddpYWqE9w0ctllWk+kTh045xx48smE73LCY4oi3pk9WyvKdO0Kf/2r39JEzOrVOov3iy8CluQvVoTulH/5i4Y0B34aesm0aKE/3e7d4ZZboF8/iIMkqEY5MUURz6xfD5deqobeCRMCE5M/ezZ07KjmptTUShDdVFb++U+1z+zc6bckMaV2be3qI4+oxbRjR51DYwQPUxTxSrjz+u23IT3db4kiYtIkrYz288+V1HEdKfPna8bZtWtL3TXIJCXBffepo3vbNjj9dC2dYmVWg4Upingl5LwePToQzmvndJbu5ZcXzrY2JVEK69fDr3+tQ7AEp1s3WLwYhgyBUaOgeXONijLfRTAwRRGPBMx5nZsL118Pd92l1rEDB8zcFDH79mmFIC/ldCJTs6Y6tmfN0my0V1yhwQ5xUJfHKAVTFPFGwJzXoRrXL76YIKnCK5pQctI+fTTLXiUgK0uzsT/2GPznP+r4fuIJ++3EM6Yo4omAOa9XrVKb8yefqNPa/ujlJGR/GTZMZ94nbAKsQlJS4Pbb4bvvoEsXnZPYsaPO8DbiD1MU8ULAnNezZukfe9kyi2yKKhMmQKtWOrW5EpCRAe+9p9ny162DDh3UPZfgAWGBwxRFvBAQ57Vzei/LztYoFotsigFLlkBmps7i3rrVb2lijog+Iy1erL6uJ5+Ek07SKSfbtvktnQGmKPwnNxfuvDMQzusfftAJVH37QkGBNlMSMUIEXnsNmjTRdPKVwBxVp47+Db78UtN/3H8/nHCCBkmsW+e3dJUbUxR+smaNGmgfe0wd2HHqvN6zBx58EFq2hKlTCyObzCcRQ0J+i/374eabdc7Fl1/6K1MF0akTfPghzJungRJ/+xs0bQoDBugkTqPiMUXhF++9pyUzFy7U0MiRI+POee0cTJ6sCuLhh/WhVsQURIVSUKCvK1dqTpRLLomPGp8VQGamzuheulTTl7/6qpZgueIKK5JU0ZiiqGhCpqaLLtJx9ddfQ69efkt1GMuXq5npkkt02J+WpgrCnNY+ERphTJ6sv5tHHqk0RT2aNdPnqNWrtbjj1KkaYtu1q9aGskl7sccURUUSbmq68UatK3HSSX5LdQhFzUypqXo/2rfPb8mMg+Tna+GHjAytdlhJsu01bKgO7rVrtQzrd99pnZMWLfQ3u2CBKY1YYYqionjvPR1Lh0xNzz6rj+lxQsjM1KKFmpkKCtTMZCOIOCRkjlq/Hvr311DqK6+Ezz8v3JbA1Kqlg/KVK+Gll6BRIx1gtWmjGVHuu0/9G6Y0oocpilgTbmo68cS4MzVt365/tjPPVDPTTz+ZmSkwhO6EeXnwxhs6Wv3FL+Chh/QumuCkpanv4tNP1Tz63HNqlXv0UY2aOukkuPtu9WeY0jhKnHMJ09q1a+fiitWrnevUSZM03Hijc3v3+i2Rc865Awece+8953r1ci4tTcVLTdVWmFMicVqN5D3uFfr5L0hFtOTkwuWsLOfGjnVu506/f3IVysaNzr34onPnnedcSop+FSee6Nzttzv373879/PPfksYXwBzXCn3VtH9EoP27du7OX6HQzin+ZrGj9fSl/n5+sju8yjCOR3MvPqqPnxu2qRBVgUF6odI5PkQNZL3MjJ/IL/jNb9FqTiSkgrVR5Uqmhqme3dNQNi0qdoVKwFbt8KUKVqS9aOPdJQsotlrO3UqbC1aaIGtyoiIzHXOtS9xH1MUUWLRIr0Djx8PK1bon7N7d60076PDeu1aeP11VRCLFun9A1S8yuKgrpSKIpzk5EMn7NWpo6G2Z5yhiiMrC6pX9028imLHDn2GmzmzsG3Zottq1tT0ISHF0bEjNGjgr7wVhSmKWLN6teazGD9ec/MkJcG55+rU5Z499Q9ZgezYoWLMm1fYQpEgKSmVtyRppVcU4SQl6Q8h3AGVlKQTFLKzNctjhw466oizeT3RxjnNNjBjRqHi+OabQp3aoIE+4zVrpq/hy+npiTMo811RiMj5wJNAMvCSc254ke3ibf8NsAe4xjn3dSTHFkeFKIoNG7SM2xtvwH//q+s6dVLl0KuXhmBUABs3FiqDr7/W9sMPhduTk9WslJKif4jKPEnOFEUppKbq3TE8YiopCY49VhXGSSdpIEZGRuHr8cdD1ap+SRwz9uzRDLazZulEvx9+0DlFa9ce6hCvXbtQaTRrpn/7Y49V5XLssdrS04OhayNRFDHrhogkA88C/wvkAF+JyBTn3KKw3S4ATvZaR2AU0DHCY6NPQYGG/axZo6OF1asPX96xQ/dt0UJj8vr0gV/+Miqnz8vTofDGjYe2TZsOfb9qlYoZIiVF/+dVquhTzv79hU9FFrlklErRH0lSkv6otmzRB6MZMw4/RgTq19cwowYNoG5dbXXqHL4ceq1VC6pVUwUTpw6B6tU1AvDMMw9dv2+f/u+WL9cWUiDz5mld8OIexESgXr1C5dGggbaaNbXVqFHy6zHH6FcV+l/7SSz1XQdguXNuBYCITAB6AOE3+x7AK57nfaaI1BGR44CMCI6NDgUFPHXKSA5s2QnbtuG8pyqHXhmX1gDqnoyrUw9a1sHVqUtBs5MpaNCQglyhYFzhw1hBweHL+/frjyy87d17aNu3T59ktm8/9KklnORk3VZQoP/jqlV1OTe38EdaSSbqGrGmoKDkH1NysiqSbdtg82b9YYbuZJFaKFJS9A5YtarGuaalqRKpVk3v1lWr6j7hLXTeoi0p6dAmcuT3ITlFCltJ7z3SRDgFOCW07QTgRIFz9ev6eV81Nu6qzqZd+rpxV3U27a7uLVdjU051Fi2pzubd1di5vwr78sp2662SnEfVlPzDW7K+Nqyzn8nfZJTpM8tCLBVFYyC8cnwOOmoobZ/GER4LgIgMBAYCnHDCCWWXMimJe3+4jj0FR5j8tg/4yWs+EnrIS0k59OmiWjX/ZAoKB3Ylk5SaAtVq+S1K4pKfry0vLzI7Z2i/PXtiL1uMSQLSvdY8wmPySGYXNdhFDXZSs9jXXdRgP1W15Xttf9XCdWEtdc1O2FFPR20xIJaKorjBUtHHjSPtE8mxutK5F4AXQH0UZREwxNpNaTz7rA4l4dCHi/AHjPBt4Q85If9g6KElfDklRW/moYemtLTChyWjYpANWzm71uVQ7RK/RTEiJVzpFG3FbQsNt0OtuPeh4X7R2Sch30xx78Nf4VA/Tvj2SEdS3n4pQB2vFc8er0VIy5ZQ49LI9y8jsbxd5QDHh71vAhTNKn+kfapEcGzUqFcPhg6N1acb/tMIuNBvIQwjsMQyhcdXwMki0lREqgB9gClF9pkCXCVKJ2C7c+6nCI81DMMwKoCYjSicc3kiMhj4NxriOto5952IDPK2Pwd8gIbGLkfHWf1LOjZWshqGYRhHJqEm3InIJmB12Kr6wGafxIklidovSNy+Wb+CR6L2rWi/TnTOlTgPPaEURVFEZE5pE0mCSKL2CxK3b9av4JGofStPvyzNuGEYhlEipigMwzCMEkl0RfGC3wLEiETtFyRu36xfwSNR+1bmfiW0j8IwDMM4ehJ9RGEYhmEcJQmpKETkfBH5XkSWi8g9fssTTURklYh8KyLzRcTncn7lR0RGi8hGEVkYtq6eiHwkIsu817p+ylhejtC3YSLyo3fd5ovIb/yUsTyIyPEi8qmILBaR70TkD976QF+3EvoV6GsmImkiMltEvvH69UdvfZmvV8KZnrwU5UsJS1EO9I15ivIKQkRWAe2dc4GO7xaRs4BdaPbgVt66R4GtzrnhnoKv65y72085y8MR+jYM2OWce8xP2Y4GL7Pzcc65r0WkJjAXuAS4hgBftxL61YsAXzOv3s8xzrldIpIKTAf+AFxKGa9XIo4oDqY3d84dAEIpyo04wjn3ObC1yOoewDhveRz6Zw0cR+hb4HHO/RQqLOac2wksRjM9B/q6ldCvQOOUXd7bVK85ynG9ElFRHCl1eaLggKkiMtdLsZ5INPRyfeG9HuuzPNFmsIgs8ExTgTLPFEVEMoC2wCwS6LoV6RcE/JqJSLKIzAc2Ah8558p1vRJRUUScojygdHbOnYpWB7zJM3MY8c8ooBmQiVY3+Zuv0hwFIlIDmATc4pzb4bc80aKYfgX+mjnn8p1zmWgG7g4i0qo8n5OIiiKS9OaBxTm3znvdCLyNmtoShQ2evThkN97oszxRwzm3wfvTFgAvEtDr5tm6JwGvO+f+6a0O/HUrrl+Jcs0AnHPbgGnA+ZTjeiWiokjYFOUicoznbENEjgG6AQtLPipQTAGu9pavBib7KEtUCf0xPXoSwOvmOUdfBhY75x4P2xTo63akfgX9molIAxGp4y1XA7oCSyjH9Uq4qCcAL4zt7xSmKH/EX4mig4j8Eh1FgKaIfyOofROR8UA2mslyA/AQ8A7wJlqReA1whXMucE7hI/QtGzVhOGAVcH3IThwUROQM4AvgWyBU6u0+1J4f2OtWQr/6EuBrJiKtUWd1MjooeNM59ycRSaeM1yshFYVhGIYRPRLR9GQYhmFEEVMUhmEYRomYojAMwzBKxBSFYRiGUSKmKAzDMIwSMUVhGGGISEZ41td4QUSmiUjC1W82goEpCsOIMSKS4rcMhnE0mKIwKi0icpuILPTaLWGbUkRknJcMbqKIVPf2Hy4ii7z1j3nrGojIJBH5ymudvfXDROQFEZkKvCIis0SkZdi5p4lIO2+2/Wjv2Hki0sPbXk1EJnjn+gdQ7Qh9WCUifxSRr0XrlJwSm2/LqMyYojAqJSLSDugPdAQ6AQNEpK23+dfAC8651sAO4EYRqYemcWjprf+zt++TwBPOuSzgMuClsNO0A3o4536Lprvv5Z37OOAXzrm5wP3Af7zjzwZGeOlZbgD2eOd6xPusI7HZSxQ5Crij3F+KYRwBUxRGZeUM4G3n3G4vZ/8/gTO9bWudc//1ll/z9t0B7ANeEpFLgT3e9q7AM14q5ylArVA+LmCKc26vt/wmcIW33At4y1vuBtzjHT8NSENTK5zlnRvn3AJgQQl9CSXnmwtkRNZ9w4gcs50alZXi0tGHKJrXxjnn8kSkA3AummhyMHAO+rB1WphC0A8XAdgd9gE/isgWL/9Ob+D6MDkuc859X8zxkebX2e+95mP/aSMG2IjCqKx8DlwiItU9U09PNDEcwAkicpq33BeY7tUqqO2c+wC4BU0WBzAVVRoAiEhofXFMAO7yPudbb92/gSFeBlPCzF+fA1d661oBrcvXTcM4ekxRGJUSr/TlWGA2mv30JefcPG/zYuBqEVkA1ENt/zWB97x1nwG3evveDLT3nM6LgEElnHYiOhp5M2zdw2iJygVeWO7D3vpRQA3vfHd5chqGL1j2WMMwDKNEbERhGIZhlIgpCsMwDKNETFEYhmEYJWKKwjAMwygRUxSGYRhGiZiiMAzDMErEFIVhGIZRIqYoDMMwjBL5f3D7euUTLC7uAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.subplot(211)\n",
    "plt.plot(ns, scipy.stats.poisson.pmf(ns, s_null + b), color='r', label='null')\n",
    "plt.plot(ns, scipy.stats.poisson.pmf(ns, s_alt  + b), color='b', label='alt')\n",
    "\n",
    "# First a loose test with a threshold at n=12 (i.e s=2)\n",
    "threshold1 = 12\n",
    "lo1 = np.arange( 1, threshold1 + 1, 1)\n",
    "hi1 = np.arange(threshold1, 30, 1)\n",
    "plt.fill_between(hi1, scipy.stats.poisson.pmf( hi1, s_null + b ), color='r', label='p-value')\n",
    "plt.fill_between(lo1, scipy.stats.poisson.pmf( lo1, s_alt  + b ), color='b', label='Type-II')\n",
    "plt.xlabel('observed n')\n",
    "plt.ylabel('PDF')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathbf{\\color{red}{\\rightarrow\\text{Compute the Type-I error probability (p-value) and the Type-II error probability for this test.}}}$\n",
    "\n",
    "Recall that\n",
    "- the p-value is the integral of the null hypothesis PDF above the threshold, so given by `scipy.stats.poisson.sf(threshold, s_null+b)`\n",
    "- the Type-II error is the integral of the alt. hypothesis PDF below the threshold,  so given by `scipy.stats.poisson.cdf(threshold, s_alt+b)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Type-I error (p-value) = -999\n",
      "Type-II error = -999\n"
     ]
    }
   ],
   "source": [
    "print('Type-I error (p-value) = %g' % -999 ) # replace with the correct expression\n",
    "print('Type-II error = %g' % -999 )          # replace with the correct expression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now for a more stringent test with a threshold at n=16 (i.e s=6)\n",
    "threshold2 = 16"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathbf{\\color{red}{\\rightarrow\\text{Compute the Type-I error probability (p-value) and the Type-II error probability for this new test.}}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Type-I error (p-value) = -999\n",
      "Type-II error = -999\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAACCCAYAAABFCHMXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAe8UlEQVR4nO2deXhUVdLG30pCEgjRhBCWASHIoIDOjEBkkZFFI4LCICqbOiguiAIuAyK44igji4DgKAwioH4CMiKLo6MggqKymLApIiEmRINAQEwgIJCk6/ujuidtyHITuvv0vV2/5zlPb/d21+nb/d5z69SpImaGoiiK4nzCTBugKIqiBAYVfEVRlBBBBV9RFCVEUMFXFEUJEVTwFUVRQgQVfEVRlBAhwrQBZVG3bl1OSkoybYaiKIptSEtLO8LMiRVtE5SCn5SUhNTUVNNmKIqi2AYiyq5sG3XpKIqihAhBOcJXgpAdO4CcHGvb1qgBdO0KREX51yZFUaqECr5SMczAxInAk09Wbb+2bYGlS4Hmzf1jl6IoVUYFXymfkyeBu+4CliwBbrsNGDUKIKp8v/R0YORIEf3584GbbvK/rYqiVIoKvlI2+/cDffsCW7cCkyYBY8daE3sAuPxyoHNnYOBA4Oab5UQxdaq6eBTFMDppq5zNli0i2nv2ACtXAo8+al3sPSQlARs2AA8/DLz0kpwAMjP9Yq6iKNZQwVd+y6JFQJcuQHQ0sHEj0KdP9d8rMhKYPh1Yvhz4/ntx8bz7ru9sVRSlSqjgK4LLBYwfD9x6K9Cxo4zyL73UN+99ww3iGrroIvHnP/ggcPq0b95bURTLqOArwPHjQL9+4qsfNgxYvRqoW9e3n9GsGfD55yL2s2YBV14JZGX59jMURakQFfxQJysLuOIK4P33xdc+Z464YvxBZCTw4ovi1klPB9q0EXePoigBQQU/lNmwAWjfXhZUffihhFJWdXK2OvTrB2zbBrRoAdx4I/B//+f/z1QURQU/ZMnOBq6/HkhIADZvBlJSAvv5HhdP167APffICUBRFL+igh+KMMuCKmYZ2V90kRk7oqJkNW7dujLqP3LEjB2KEiKo4Ici//oXsHYt8MILEi9vknr1xKd/8KAs1CoqMmuPojgYFfxQIysLGDNGXDjDhpm2Rrj8cpks/uQTYNw409YoimPR1AqhhMsF3HknEBYGvPZaYCZorXLHHUBaGjBtmizQuuUW0xYpiuPQEX4o8corwPr1wIwZQJMmpq05m+nTJT7/7ruB7dtNW6MojsOS4BNRTyLaQ0QZRHTWNTcRtSSijUR0mojGVGVfJUBkZEhOnF69ZJQfjNSoAfz730CdOjKJ+/PPAflYl8t6Yw6ISYriFyoVfCIKB/AygF4AWgMYTEStS212FMADAF6oxr6Kv3G5gKFDRVBffTW4XDmlqV9fJnF/+gkYNMhvk7iFhRL+37YtEB5uvbVqJQuF8/P9Ypai+BUrPvz2ADKYORMAiGgJgL4AvvVswMy5AHKJ6Pqq7qsEgFmzJOb99deBRo1MW1M57dsDs2dL6OhjjwFTpvjsrfPygLlz5SvZvx9o3Vpqu0RY+Ce4XMBHH0l2iMcfB4YMAUaMkPdQFDtgRfAbAfjR63EOgA4W39/yvkQ0DMAwAGgSjP5lu5KeLknR+vQB/vpX09ZY5847ZRJ36lQZhg8adE5vl50NzJwpFzgFBcBVV8n9a6+VOWyrTJgApKYCL78s896vvCLvNWoU0Lu3tROHopjCyk+9rOt/q55My/sy81xmTmbm5MTERItvr1RIcbFEv9SsKbH3wezKKYsZMySP/p13Sk3dapCaKueK5s0lVZCnpsvatTKdURWx95CcDCxYAPz4I/D888DevTLl0Ly55J/T9WNKsGLl554D4AKvx40B/GTx/c9lX+VcmTFDctr/859Aw4amrak6kZHAO+8A8fGiqEePWtrN5QLee0+yNlx+OfDf/0odlsxM8du3aeMb8xITZdlAZqbkgGvRQi6mGjeWKZP0dN98jqL4DGausEHcPpkAmgGIBLADwCXlbDsBwJjq7Ovd2rVrx8o58u23zFFRzP36Mbtcpq05NzZuZI6MZO7Rg7moqMJN8/OZe/ViBpibNGGePl2eCxS7djHffz9zTIy0N94I3GcroQ2AVK5MzyvbQN4H1wFIB/A9gMfdzw0HMNx9vwFkNH8MQJ77/nnl7VtZU8E/RwoLmdu3Z05IYD540LQ1vmHuXPm5TplS7iaZmcyXXMIcEcE8a5Z8DabIyWHu0kVMHjqU+cQJc7YooYHPBD/QTQX/HHn+eTm0b79t2hLf4XLJ1UpkJPOOHWe9/MUXzImJzHFxzB9/bMC+MigsZH7iCWYi5tatZfSvKP7CiuDrSlun8c03wNNPA/37AwMGmLbGdxDJxHN8vEQbeZVIfOstoHt34PzzgU2bgKuvNminFxERwLPPSgGxI0dksnfhQtNWKaGMCr6TYJbc8uefL3GDTiMxEZg3D9i5E3j6abhcwFNPAbfdBnTqJGJ/8cWmjTyblBTJFNGxo0zm3n67hIYqSqBRwXcSb78tqjd5soijE+ndG7jnHpyc/BIGpRzGs89K1Obq1VLLJVhp2BBYs0bi+N98U6KHvv7atFVKqEEchMlBkpOTOTU11bQZ9uLUKaBlS3F5pKVVL8DcJhzIOIG+l2Yg9fQfMOW5Mxj9WLStlhisWyfJQPPyZG3AXXfZb4mEEnwQURozJ1e0jXNVIdSYOVOWk06b5mix374daN89BrvoUiynmzAma6TtxLJ7d+nHlVeKB+6224Djx01bpYQCzlWGUCI3F/jHPyR9wlVXmbbGb6xcCfz5z3L/8y/D0ffRlpLfYNUqs4ZVg/r1pbrkc88BS5bIYdMVuoq/UcF3AhMmACdO+DTJWLCxYIEstm3dGtiyxb1a9plngD/9SYbJubmmTawyYWGShG3lSgmu6tJFEropir9Qwbc7u3dL+sf77hMfvgNZsUJqoqSkSP2W/2WJiIyUXAl5eVKuMQjno6zQu7eM9nNy5AomI8O0RYpTUcG3O488AtSuLbH3DuSzzyT5WXKypMmvVavUBpdeKu6slSttHeTetatM5h4/LqK/c6dpixQnooJvZz7+GHj/ffEL1K1r2hqfs2MH8Je/AM2aSTdr1y5nw4cfFsV88EFg375AmuhT2rUDNmyQBVtdu0reO0XxJSr4dqW4GBg9WtRw1CjT1viczEygZ08gNlaKjlR4PgsLk+IugKxqKi4OiI3+oFUrqVVTt664sNasMW2R4iRU8O3KwoVy3T9pEhAdbdoan3LoENCjB3DmjIi9pXo4TZtKGavPPpO00DYmKUlG+r//PXD99cCyZaYtUpyCCr4dKSgAnnhC8gn072/aGp+Sny8j+wMHxI1TpfKBt98O3HCDuLhs7gRv0AD49FNZkTtgADB/vmmLFCeggm9Hpk4FDh4Epk931BLNU6dEr7/5RuqedOxYxTcgkoiluDhJsHbqlB+sDBxxcZIyIiVFVuNOn27aIsXuqODbjZwcEfyBA6uhiMFLcbGkG1i/XrxVvXpV840SE2Ux1s6dUo7K5sTEyLqy/v1lyubJJ20bfaoEAVpy2W488YTU8Js0ybQlPoNZlhEsXw68+CJw663n+Ia9ewMPPCDpJlJS5LGNiYoCFi8GzjtPVuYWFDju4k4JECr4dmLrVuCNNyT2PinJtDU+46mngFdfBR57TCIrfcLkyeIEHzpU4jt/9zsfvbEZwsPlO4qJkZNiUZHMUavoK1VBXTp2gVmu6RMSRBkdwqxZMmq9+2659RnR0ZKk5uRJ8efbOFTTA5GI/ejRUpf+/vvlYk9RrKKCbxfee08c3M88IwVOHMA778iIvl8/YPZsP4xWW7aUM8onn8i8hwMgkq6MGwfMmQPce6+KvmIdzYdvBwoLJYVAWJhUzYiwvydu0yZJE9y2LbB2rR+XEjBLboZly2RFk0MmupnFFfbccxKN+tpr4vZRQhcr+fDtrxyhwOzZQHq6jPIdIPZZWUDfvuJWX7HCz+vGPLVwN2+WMKBt2xxxhUQk9XJr1JA0SkVFEt3kgJ+H4kfUpRPs5ObKUC4lRZZd2py8POnGmTOysCoglRjj4oBFi4AffgCGD3dUXKNnlP/WWzJVUVRk2iIlmFHBD3YefVQmHl96yfYhGYWFwM03A3v3SubLgGZzvuIKmf9YsqQk745DePxxCUpaskS8V4WFpi1SghUV/GDmyy/lOv1vf7N9rntmiSpZu1bCC7t3N2DEuHFAt27AyJHAnj0GDPAfY8dKbP6yZZKK4cwZ0xYpwYgKfrBSVASMGAE0biyLrWzO1KnAvHkyGr3jDkNGhIdLwZToaGDwYOD0aUOG+IeHH5agpBUr5ErKYd1TfIAKfrAyZ45Uup4+vYJE8PbgnXfEMzVwIPD3vxs2plEjqZe4bRswfrxhY3zPqFHAK6/I/H6/fsCvv5q2SAkmVPCDkdxcGdWnpMhQzcZs3iyTiZ06iXcqLBh+cX36iDLOmAF88IFpa3zOffdJDrkPP5TMo/n5pi1SggVLfz8i6klEe4gog4jOykhFwiz36zuJqK3Xa/uI6Gsi2k5EGlxvBYdM1O7bJxWrGjaUCoRBlbZ/yhTgj38U/9KBA6at8Tn33COBSV9+KdMWhw6ZtkgJBioVfCIKB/AygF4AWgMYTESls5T3AtDC3YYBmF3q9e7MfFlliwIUOGaiNj9fwi9Pnw5g+GVV8KReKCiQbG0ODG0ZNEhcO+npUic3K8u0RYpprIzw2wPIYOZMZj4DYAmAvqW26QvgDRY2AYgjooY+ttX5OGSitrBQ0vmmp0v4ZatWpi0qh1atZK5k3TrJUeCg+HwPPXtK6eOffwY6d5ZaA0roYkXwGwH40etxjvs5q9swgNVElEZEw8r7ECIaRkSpRJR6+PBhC2Y5EAdM1DJL1OOaNbLA9aqrTFtUCUOGyOqlBQt8nL0teOjUSUomEgFXXikXkUpoYkXwy3Iilx4KVbRNZ2ZuC3H7jCCiLmV9CDPPZeZkZk5ODLrr/wDgkInaCRNkwnD8eODOO01bY5EJE0qE/403TFvjFy65BPjii5Li6P/9r2mLFBNYEfwcABd4PW4M4Cer2zCz5zYXwHKIi0gpjQMmal98UcIuhw4FJk40bU0VIJLVYFddJbUEP/nEtEV+ISlJ8se1bCmT6YsXm7ZICTRWBP8rAC2IqBkRRQIYBGBVqW1WARjijtbpCCCfmQ8QUQwRxQIAEcUA6AFAvYilccBE7YIFsvDnxhtlhG+7c1ZkpCxTvfhi6cSuXaYt8gv168uURefOMlf9z3+atkgJJJUKPjMXARgJ4CMAuwEsZeZdRDSciIa7N/sAQCaADACvArjf/Xx9AJ8T0Q4AWwC8z8wf+rgP9sYBE7XvvisFTK65RkIBbZuxMS5O4vJr1QKuu86R4ZqAJAv98EMZ5Y8aJR4tB85XK2XBzEHX2rVrxyHDSy8xA8xLl5q2pFqsWcMcGcncsSPz8eOmrfERW7cyx8Qwt23roE6dTWEh89Ch8vMbMUIeK/YFQCpXoq3BsO4xdDl0yNYTtZs2ATfcIF6QDz6wbWDR2bRpAyxdKrVwBw50bM7hiAgpnPLII8DLL0sIZ6gGyIUKKvgmGTfOthO1X38N9Oolq2hXrwbi401b5GOuu05U8IMPJM7UoT4PIll0vGCBTOgmJwNabM65qOCbYsMG207UZmQAPXoAMTESb9+ggWmL/MS998pJ+V//ElV0MHfcIWGbgKzKnT/fqDmKn1DBN8HBg5KeNynJdhO1+/eLB6qwUMQ+Kcm0RX5m4kTJUTBunKRicDDt2gFpaSL4d90lxcE0xbKzUMEPNGfOSN6Bo0clcbmNHN9HjkgkztGjEuURtCkTfElYmFyJXXmlVAv/7DPTFvmVunXl2D76qFzYdOsmJ3nFGajgB5qHHxZn6fz5wJ/+ZNoayxw7Jj77rCxJyJUcSmnwoqLk5HzhhTKzuXy5aYv8SkQEMGmS1DH45hugbVvHn+dCBhX8QDJ/vlSnGDNG3AQ2IT9fYra3bwf+/W+ga1fTFhmgTh1g/XpJqXzTTZLvyKETuR5uuknqGcTFySLkmTMd32XHo4IfKLZskcoUKSnA88+btsYy2dmyKvOLL4A33wR69zZtkUE8y1RvvBEYPVqidxwasumhdWv56fbuDTz0EHDbbcCJE6atUqqLCn4gOHhQRKJRI5n4s8lS1C1bgA4dgJwc4KOPbHVR4j9q1pQY/bFj5WrtL38Bjh83bZVfOf98WU09caLk3+nQQdZgKPZDBd/feE/SLl8OJCSYtsgSy5bJhF2tWsDGjTZIcxxIwsKAyZNlVnP1aglryckxbZVfCQsDHntMJnTz8oArrpCMIFo+0V6o4Psbm03SMgNTp8rC38suk5FcSETjVIdhw6ScV1aWDHu3bTNtkd/p0QPYvVty8MyeLb+Nd95R375dUMH3JzabpC0slNjrsWMlo8DatUC9eqatCnKuvVYmOMLDJXTzP/8xbZHfiY2VCdzNm2Vao39/8Wz98INpy5TKUMH3FzabpPXUoJ07Vy7dFy0Sd7VigT/8QdSvZUugb19JlRECXH458NVXwAsvSAmB1q2BGTMcP49ta1Tw/YHNJmn37ZNInHXr5KJk4kTx2SpVoGFD4NNPgT59gAcekJCW4mLTVvmdiAgJWNq1S8J1//Y38W6lpZm2TCkL/Vv7GptN0noicfbvl0icoUNNW2RjYmJktvuhh8Tn0a2bZNwMAZKSxJu1dCnw009A+/YyfeXwACbboYLva2w0Sfv22zIqi4nRSByfER4ufo2FC4HvvpNlqiNGyADA4RDJWGf3bsk7N3Mm8PvfA//4h0T2KEFAZQnzTTRbFkA5c4Z5zBipJjFmjGlrKiQjg/n668XUTp2Yc3NNW+RQjh5lHjmSOSyMOSGBec4c5qIi01YFjI0bmXv2lN9ZbCzzI48w799v2irnAgsFUIyLe1nNdoKfnS3KCTAPHx60pYNOnmR+6inmqCjm2rWZp02T85TiZ3bsYO7SRX4fbdsyf/mlaYsCyrZtzIMGyXkvMpL57ruZ9+wxbZXzUMEPBO+9x1ynjgxhliwxbU2ZuFzMK1cyJyXJEb/lFh1pBRyXi3nxYuZGjeQgDBnCfOCAaasCSkYG8333yYCDiPnmm5m/+sq0Vc7BiuCrD7+6FBZKbbg+fYAmTSQsYeBA01adxfffSx6Uvn3FV79uHfDWW8DvfmfashCDSNZifPcdMH68RG9ddBEwbZpM9IcAzZvLspTsbCkvsHq1hHampAAff6yLtwJCZWcEEy3oR/jeLpz77mP+9VfTFp3FiRPMTz4po6nYWHXfBB1795ZMpLRsyfz668wFBaatCij5+cxTpjA3aFDyNTz5pHjAXC7T1tkPqEvHD7z3HnN8fNC6cFwu5hUrmJs2VfeNLfjPf5gvvlgOVu3azEOHMn/6aUgp3q+/Ms+bx9ytm/j5AeYWLZjHj2feujWkvopzQgXfl3hH4Vx2GXN6ummLfkNeHvOrrzJ37iwmXnIJ8/r1pq1SLOFyMW/YwHzXXSL6APOFFzI/8wxzVpZp6wLKwYMSzHT11czh4SVfxdix4u9X8S8fFXxfkZ3N3LGjfF333x80LpwzZ2SAOGAAc3S0mHfRRcwzZ6r7xrYUFDC/+aYoHpEc1O7dQ9Llk5srg5hrr2WOiJCvomlT5tGjmT/6iPmXX0xbGFxYEXyS7YKL5ORkTk1NNWsEsyxDXbwYeP11WSY/bx4wYIBxs7ZulWIkixYBhw/LYt5Bg4AhQ2QSjMioiYqvyM6WA71wocy+164taUyvuQbo1EmWt4bIwT56FFi1SiqurVkjMRNEkq2zY8eS1rq1rH0LRYgojZkrLD6qgl+ab78VJV28GMjMBCIjJcxl8mRZNmiIH3+U6Jo33xQTIyMlQOivf5Vas5GRxkxT/A2zZORcuFByF3jyFdSvLyrXqZO05GQpYOBwjh2TsdimTSXt55/ltdhYSevgOQF06AAkJpq1N1Co4FslO1vC5BYvltwnYWHA1VcDgwcD/fpJUc8AcuyYmLFtW0nbuVP+91dcISP5AQOA+PiAmqUEA0VFUll840ZpmzYBe/fKa+Hhks7DcwJo3x5o1izok/edK8xyAeT5OjZtkv+PJ3ddYqKM1Zo3l1vv+wkJzrlI8pngE1FPADMBhAOYx8yTSr1O7tevA3ASwB3MvNXKvmUREME/dEgSXS1aJKMnQIYEgweLmjZo4N/Pd5ObWyLqW7fKbUZGyev16gFt2ojQ33KL0YsMJVg5ckRUzqN4W7YABQXyWlgY0Lgx0LSpuIBK315wARAVZdB4/3DypCyN2bwZSE+XE0JGhlwpe0ve+eeXiH/z5vK3r1dPThL16klLSLDHOdMngk9E4QDSAVwDIAfAVwAGM/O3XttcB2AURPA7AJjJzB2s7FsW5yz4LpekKM7OlqoM2dm/bT/8IMNoQJx+t94qTvALL6z+Z3pRVCSXmLm5v22HD//2cXa2ZBb00KyZiLt3a9jQOSMQJUAUF8tVQFqaVOPKzpYc2NnZUorR5SrZlkh+ZBdcANSpI1ez8fElt973PbfnnSfFEqKibOcwP3VKvoqMDGmeE8H338tXVVYufyL5ajwngcREabGx0mrXrvg2Jka+qshI//6XrQi+lfNWewAZzJzpftMlAPoC8BbtvgDecM8UbyKiOCJqCCDJwr6+weXCrFav4MzPx4G8fLD7eo4h3zBHJwLxLcBxdYBL4sBx8XA1bwFXYn24Cgku97ysyyWt9P3Tp+XH4t1+/fXs506cAH75pexVg+HhQN26JSOHlBS5Am/TRsoJqotG8Qke105Z2VoLC0X0vU8C+/bJc4cPy3A4L0+alXz+NWoA0dFyAoiO/m2rWVNULiLity08/OznIiLkasTTiCq+71FOopJW0WM30URoCaCl57UmAJoScLX81385VRO5BbVwuEBucwtq4fCJWu77NXE4pxa+/a4WjpyoieOnI3GqqGpD/8jwIkRFFJfdwotRP+40Vu5IqtJ7VgUr1jYC8KPX4xzIKL6ybRpZ3BcAQETDAAwDgCZNmlgwqxRhYRj//d04WRxd9uunABxwt/LfAuHhJb8t7/ulf8uelpBw9m/c+3LQu8XHa2ERxTA1asilZLNmFW/HLG6hX34R8fe+PXbs7JFOeaOgggI5cRQVVdwKC0VxmUtGWuXd9xNhABLczWoZ5yKEowC1UYDaOI7YMm8LUBunESWt2N1OR5U859VqHCqEjJP9gxXBL+sipPT4tbxtrOwrTzLPBTAXEJeOBbvO4uAvIvbeJ/eKbr1FXd0miuIFUYnPojoDMH/iuXyW0PzKH5fer/R7lPW6RRsiAMS5W6X2WsHPQmRF8HMAXOD1uDGAnyxuE2lhX58RG+uvd1YUJWgoPXJTLGPFwfAVgBZE1IyIIgEMArCq1DarAAwhoSOAfGY+YHFfRVEUJQBUOsJn5iIiGgngI0ho5Xxm3kVEw92vzwHwASRCJwMSljm0on390hNFURSlQoJy4RURHQaQ7fVUXQBHDJnjT5zaL8C5fdN+2Q+n9q10v5oyc4XrioNS8EtDRKmVxZfaEaf2C3Bu37Rf9sOpfatOvzRIUFEUJURQwVcURQkR7CL4c00b4Cec2i/AuX3TftkPp/atyv2yhQ9fURRFOXfsMsJXFEVRzpGgFnwi6klEe4gog4jGmbbHlxDRPiL6moi2E5Hh8l7Vh4jmE1EuEX3j9VwdIlpDRHvdt7ZMC1dO3yYQ0X73cdvuzhRrK4joAiJaR0S7iWgXET3oft7Wx62Cftn6mBFRNBFtIaId7n49436+yscraF061U2tbBeIaB+AZGa2dXwwEXUBUADJlnqp+7kpAI4y8yT3iTqemR81aWd1KKdvEwAUMPMLJm07F9yZbBsy81YiigWQBuAGAHfAxsetgn4NgI2PmbveSAwzFxBRDQCfA3gQwI2o4vEK5hH+/9IyM/MZAJ7UykoQwcyfATha6um+AF53338d8qezHeX0zfYw8wFPgSJmPg5gNySzra2PWwX9sjXuGuXuijao4W6MahyvYBb88lIuOwUGsJqI0typoZ1EfXcuJbhv6xm2x9eMJKKdbpePrdwepSGiJABtAGyGg45bqX4BNj9mRBRORNsB5AJYw8zVOl7BLPiWUyvblM7M3BZALwAj3O4DJfiZDaA5gMsg1RWmGbXmHCCi2gCWAXiImY+ZtsdXlNEv2x8zZi5m5ssgGYfbE9Gl1XmfYBZ8K2mZbQsz/+S+zQWwHOLCcgqH3P5Uj18117A9PoOZD7n/fC4Ar8Kmx83tC14G4C1mftf9tO2PW1n9csoxAwBmzgOwHkBPVON4BbPgOza1MhHFuCeVQEQxAHoA+KbivWzFKgC3u+/fDmClQVt8iucP5qYfbHjc3JOArwHYzczTvV6y9XErr192P2ZElEhEce77NQGkAPgO1TheQRulA/yvOPqLKEmtPNGsRb6BiC6EjOoBSVG9yK59I6LFALpBMvcdAvA0gBUAlkIqhv4AoD8z227ys5y+dYO4BhjAPgD3evyodoGI/gxgA4CvAXhqBj4G8Xfb9rhV0K/BsPExI6I/QiZlwyGD9KXM/HciSkAVj1dQC76iKIriO4LZpaMoiqL4EBV8RVGUEEEFX1EUJURQwVcURQkRVPAVRVFCBBV8RVGUEEEFX1EUJURQwVcURQkR/h+R4vN128KrbgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.subplot(212)\n",
    "plt.plot(ns, scipy.stats.poisson.pmf(ns, s_null + b), color='r', label='null')\n",
    "plt.plot(ns, scipy.stats.poisson.pmf(ns, s_alt  + b), color='b', label='alt')\n",
    "pass # draw the shaded regions for the test at threshold2\n",
    "print('Type-I error (p-value) = %g' % -999 ) # replace with the correct expression\n",
    "print('Type-II error = %g' % -999 )          # replace with the correct expression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathbf{\\color{red}{\\rightarrow\\text{Now compute the Type-I error probability (p-value) and the Type-II error probability for the range of thresholds below, and draw the ROC curve.}}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "thresholds = np.arange(1, 30, 1)\n",
    "pass # replace with code to plot the ROC curve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The threshold for discovery in physics is often set at $5\\sigma$, which corresponds to a p-value of about $3 \\cdot 10^{-7}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathbf{\\color{red}{\\rightarrow\\text{What threshold would one use to get so small a Type-I error rate ?}}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pvalue= 2.866515718791933e-07 threshold= -999 power= -999\n"
     ]
    }
   ],
   "source": [
    "pfive = scipy.stats.norm.sf(5)\n",
    "threshold_five = -999 # replace with the correct expression\n",
    "print('pvalue=', pfive, 'threshold=', threshold_five, 'power=', -999) # replace -999 with the expression for the power"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should find $n_{\\text{obs}} \\ge 29$. If we do osberve such a large n, then we can indeed reject the null at with a p-value corresponding to a $5\\sigma$ threshold. This corresponds to testing for discovery: we reject the no-signal hypothesis in favor of an hypothesis with non-zero signal.\n",
    "\n",
    "However the power of the test is also very low (you should have found about $4\\, 10^{-4}$ above), which means that even if the $s=5$ signal is there, we are very unlikely to observe an event count of $n_{\\text{obs}} \\ge 29$. The underlying problem here is the fact that the difference between the alternate hypothesis $s=5$ and the null $s=0$ is quite small on the scale of the measurement uncertainties, so the two hypotheses are difficult to separate.\n",
    "\n",
    "If the alternate hypothesis was larger, say $20$ signal events, then things would be easier. The threshold for discovery would still be 29 events (this only depends on the null, not the alt), but now the power would be larger. Let's compute it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_alt2 = 20"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathbf{\\color{red}{\\rightarrow\\text{Compute the power for this hypothesis}}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pvalue= 2.866515718791933e-07 threshold= -999 power= -999\n"
     ]
    }
   ],
   "source": [
    "pfive = scipy.stats.norm.sf(5)\n",
    "threshold_five = -999 # replace with the correct expression (same as above)\n",
    "print('pvalue=', pfive, 'threshold=', threshold_five, 'power=', -999) # replace -999 with the expression for the power (not the same as above)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should find that the power is about 50% -- i.e. this means that if $s_{\\text{alt}} = 20$ is indeed true, then we have an approximately 50% chance of obtaining $n \\ge 29$ and therefore of making the discovery."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Discovery testing as a hypothesis test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Discovery testing can be seen as a hypothesis test like the one we have studied above. The setup is as follows:\n",
    "* **Null hypothesis**: no signal\n",
    "* **Alternate hypothesis**: any positive signal\n",
    "* **Discriminant**: the likelihood ratio -- or in simple cases the measured yield.\n",
    "* **p-value**: the tail integral seen above (false positive: probability to get an outcome that is larger than the threshold, even though the null is true)\n",
    "\n",
    "We need to add one more element though: since the alternate can be *any* positive signal, how do we decide what value to use ? The natural solution is to use the best-fit value $\\hat{s}$, so that the discriminant is finally\n",
    "$$\n",
    "q_0 = -2 \\log \\frac{L(s=0)}{L(\\hat{s})}.\n",
    "$$\n",
    "the numerator is for the null, and uses $s=0$. The denominator uses $\\hat{s}$ for the alternate.\n",
    "\n",
    "\n",
    "First, let's check that this works for single bin counting. Recall that this is all for a Gaussian likelihood:\n",
    "$$\n",
    "L(s) = \\exp\\left[-\\frac{1}{2} \\left(\\frac{(n - (s+b)}{\\sqrt{s+b}}\\right)^2 \\right]\n",
    "$$\n",
    "First, what is the best-fit value $\\hat{s}$ for a given $n$, the value that maximizes $L$ ? \n",
    "\n",
    "We need the exponent of the Gaussian to be zero, so $n = \\hat{s} + b$  and $\\hat{s} = n - b$.\n",
    "\n",
    "So what is the value of $q_0$ ? From a simple computation,\n",
    "\n",
    "$$\n",
    "q_0 = \\left(\\frac{n - b}{\\sqrt{b}}\\right)^2 = \\left(\\frac{\\hat{s}}{\\sqrt{b}}\\right)^2\n",
    "$$\n",
    "\n",
    "So it is equal to the square of the significance! \n",
    "\n",
    "In other words, once we have $q_0$,\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "q0 = 24"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can get the significance as\n",
    "$$\n",
    "z = \\sqrt{t_0}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.898979485566356"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sqrt(q0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and the p-value as"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.816785043215482e-07"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scipy.stats.norm.sf(np.sqrt(q0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course this all is a bit overkill here, since we already had the answer using simpler methods. But now that we validated the technique on a simple example, we can apply it to a less trivial one, in the next section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Testing for discovery in a histogram\n",
    "\n",
    "Now that we have the general technique, let's apply all this to the full binned distribution above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAdnUlEQVR4nO3de5RcZZnv8e8vl0nbSgSShhNI0h0wghASLgGjCOQIWcoZJoAwGk5HYwzT6HAE74aJa/Ay7YEFesRZA2MHgWh6UESEMOOgGGguKjCBgZCLkVsSWiIJUSAmEhLynD9qd1npVHdXX3bt6q7fZ61eVfvde9d+CpI8vd/33c+riMDMzAxgWNYBmJlZ5XBSMDOzPCcFMzPLc1IwM7M8JwUzM8sbkXUA/TF27NhoaGjIOgwzs0Hl0UcffSki6ortG9RJoaGhgRUrVmQdhpnZoCJpQ1f73H1kZmZ5TgpmZpbnpGBmZnmDekzBzKwru3btor29nddeey3rUDJTU1PD+PHjGTlyZMnnOCmY2ZDU3t7OfvvtR0NDA5KyDqfsIoKtW7fS3t7OpEmTSj7P3UdmNiS99tprjBkzpioTAoAkxowZ0+s7JSeFErS2ttLQ0MCwYcNoaGigtbU165DMrATVmhA69OX7u/uoB62trTQ1NbFjxw4ANmzYQFNTEwCNjY1ZhmZmNuB8p9CDRYsW5RNChx07drBo0aKMIjKzwejLX/4yV199dZf7b7/9dtasWVPGiIpzUujBxo0be9VuZtYXTgqDxMSJE3vVbmaDUxpjh83NzRxxxBGcccYZrFu3DoDFixdz4oknMm3aNM477zx27NjBr371K5YtW8bnP/95jj32WJ555pmix5VFRAzanxNOOCHStnTp0qitrQ0g/1NbWxtLly5N/dpm1ndr1qwp+dg0/p6vWLEipkyZEtu3b49XXnklDj/88LjqqqvipZdeyh+zaNGi+Pa3vx0REfPmzYsf/ehH+X1dHddbxf47ACuii39XfafQg8bGRlpaWqivr0cS9fX1tLS0eJDZbAhJY+zwgQce4Nxzz6W2tpbRo0cze/ZsAFatWsUpp5zCMcccQ2trK6tXry56fqnHDTTPPipBY2Ojk4DZEJbW2GGxKaEf/ehHuf3225k2bRo33XQTbW1tRc8t9biBltqdgqQbJG2WtKrIvs9JCkljC9ouk/S0pHWS3pdWXGZmnaUxdnjqqafyk5/8hD//+c9s27aNO++8E4Bt27Yxbtw4du3atde4xX777ce2bdvy210dl7Y0u49uAt7fuVHSBGAWsLGg7ShgDnB0cs61koanGJuZWV5zczO1tbV7tdXW1tLc3Nznzzz++OP50Ic+xLHHHst5553HKaecAsDXvvY13vnOdzJr1iyOPPLI/PFz5szhqquu4rjjjuOZZ57p8rjUdTXYMBA/QAOwqlPbrcA0YD0wNmm7DLis4JifAe/q6fPLMdBsZoNTbwaaI3KDzfX19SEp6uvrh8xkkt4ONJd1TEHSbOB3EfFEp762Q4GHCrbbk7Zin9EENEF5p4XOnDkToGz9emZWXh47zCnb7CNJtcAi4B+L7S7SFsU+JyJaImJ6REyvqyu6xKiZmfVROe8UDgcmAR13CeOBxySdRO7OYELBseOBF8oYm5mZUcY7hYh4MiIOioiGiGgglwiOj4jfA8uAOZJGSZoETAYeKVdsZmaWk+aU1JuBXwNHSGqXtKCrYyNiNXALsAa4C7g4It5IKzYzMysute6jiLigh/0Nnbabgb7P/zIzs35zmYsK5wV+zAaGNLA/PVm/fj1TpkwZ8O/R1tbGWWedNeCf28FlLiqYF/gxs3LznUIF8wI/ZoPb7t27mTdvHlOnTuX8888vWv768ccfZ8aMGUydOpVzzz2XP/7xj0Du2agvfvGLnHTSSbz97W/ngQce2Ou8PXv2MHnyZLZs2ZLfftvb3sZLL73Ur5idFCqYF/gxG9zWrVtHU1MTK1euZPTo0Vx77bX7HPORj3yEK6+8kpUrV3LMMcfwla98Jb9v9+7dPPLII3zrW9/aqx1g2LBhzJ07N9+l/Itf/IJp06YxduxY+sNJoYJ5gR+zwW3ChAmcfPLJAMydO5cHH3xwr/2vvPIKL7/8MqeddhoA8+bN4/7778/v/8AHPgDACSecwPr16/f5/I997GN873vfA+CGG25g/vz5/Y7ZSaGCpVGky8zKp3Pp7GKltLszatQoAIYPH87u3bv32T9hwgQOPvhg7rnnHh5++GHOPPPMvgebcFKoYB0L/HT8wfACP2aDy8aNG/n1r38NwM0338x73vOevfa/9a1v5YADDsiPF3z/+9/P3zWU6sILL2Tu3Ll88IMfZPjw/heX9uyjCtfY2MjixYsBF+Mz648oWk0tXe94xztYsmQJF110EZMnT+YTn/jEPscsWbKEj3/84+zYsYPDDjuMG2+8sVfXmD17NvPnzx+QriNwUjAzS0VDQwNr1qzp8bhjjz2Whx56aJ/2wl8Cx44dmx9TmDlzZr5qM8ATTzzBtGnTBmzNBScFM7NB6oorruC6664b0Idaqzop9HLMp9fnZHG7ambVY+HChSxcuHBAP9MDzWaduLTI0BFV/ptZX75/Vd8pmHXm0iJDR01NDVu3bmXMmDG9ngo6FEQEW7dupaamplfnOSlYt6ptGdLuSos4KQwu48ePp729PV8GohrV1NQwfvz4Xp3jpGBWwKVFho6RI0cyadKkrMMYdDymYFYg69IiHs+wrDkpmBXIsrRIx3jGhg0biIj8eIYTg5WTk4JZgY7SIvX19Ugqa2kRl0q3SuAxBbNOGhsbMxlU9niGVYLU7hQk3SBps6RVBW1XSfqNpJWSfiJp/4J9l0l6WtI6Se9LK65K0Ztl/+67L/czkEsFWuXJejzDDNLtProJeH+ntruBKRExFfgtcBmApKOAOcDRyTnXSup/uT+zQcSl0q0SpJYUIuJ+4A+d2n4eER1FwR8COibQng38ICJ2RsRzwNPASWnFZlaJshzPMOuQ5UDzx4D/TN4fCjxfsK89abMq1draSk1NDZKqampmY2Mj69evZ8+ePaxfv94Jwcouk6QgaRGwG+j4m16sF7xo0Q5JTZJWSFpRzU8qDmUdUzN37twJ4KmZZmVU9qQgaR5wFtAYf6nW1A5MKDhsPPBCsfMjoiUipkfE9Lq6unSDtUx4aqZZdsqaFCS9H/giMDsiCv/WLwPmSBolaRIwGXiknLFZ5fDUTLPspPacgqSbgZnAWEntwOXkZhuNAu5OqhY+FBEfj4jVkm4B1pDrVro4It5IK7a+acs6gKoxceJENmzYULTdzNKVWlKIiAuKNH+3m+ObAc+9M5qbm/cqXw2emmlWLn6iuQqlveIc9G/VuY4ZNwsWLGDnzp3U19fT3NzsmThmZeCkYBWpsbGRxYsXA9WzloNZJXBBPDMzy3NSMDOzPCcFsyJmzpyZX4rUrJo4KZiZWZ6TgpmZ5Xn20aDQlnUAZlYlfKdgZmZ5Xd4pSLqTLiqVAkTE7FQisiGt0h+cM6t23XUfXZ28fgD4H8DSZPsCYH2KMZmZWUa6TAoRcR+ApK9FxKkFu+6UdH/qkZmZWdmVMqZQJ+mwjo2ktLUXMjAzG4JKmX30aaBN0rPJdgNwUWoRmZlZZnpMChFxl6TJwJFJ028iYme6YZlVr44nqV0I0LLQ3eyjD3Sx63BJRMRtKcVkFaUt6wDMrIy6u1P4m272BeCkYGY2xHQ3+2h+OQMxM7Ps9TimIOlg4OvAIRFxpqSjgHdFRJdLa5pVIj84Z9azUqak3gT8DDgk2f4t8KmU4jEzswyVkhTGRsQtwB6AiNgNvNHTSZJukLRZ0qqCtgMl3S3pqeT1gIJ9l0l6WtI6Se/rw3exAddKbgbysOS1NctgzKwMSkkK2yWNIamDJGkG8EoJ590EvL9T20JgeURMBpYn2yRdUnOAo5NzrpU0vJQvYGlpBZqADeT+129Itp0YzIayUpLCZ4Bl5Kai/hL4HvDJnk6KiPuBP3RqPhtYkrxfApxT0P6DiNgZEc8BTwMnlRCbpWYRsKNT246k3cyGqh6TQkQ8BpwGvJvck8xHR8TKPl7v4IjYlHzuJuCgpP1Q4PmC49qTtn1IapK0QtKKLVu29DEM69nGXrbbYNfa2kpNTQ2SaGhooLXVd4XVqLuH194bEfcUeYjt7Sk8vFZsjkfReRwR0QK0AEyfPt1zPVIzkVyXUbH2cmkr47WqW2trK01NTezcmStWsGHDBpqamgBobGzMMjQrs+7uFE5LXv+myM9Zfbzei5LGASSvm5P2dmBCwXHjgRf6eA0bEM1Abae22qTdhppFixaxY8fe3YU7duxg0SJ3F1ab7h5euzx5e2FE9DjbqETLgHnAFcnrHQXt/ybpm+Smvk4GHhmga1qfdPx2uIhcl9FEcgnBvzUORRs3Fu8W7Krdhq5SBpqflnRVMkOoZJJuBn4NHCGpXdICcslglqSngFnJNhGxGrgFWAPcBVw8gInI+qyR3HpKe5JXJ4ShauLE4t2CXbXb0FVK6eyp5KaLXi9pGHADuZlCr3Z3UkRc0MWu07s4vhn3TdgQVelPUzc3N9PU1LRXF1JtbS3Nzf4rWW1KmX20LSIWR8S7gS8AlwObJC2R9LbUIzSz1DU2NtLS0sKoUaMAqK+vp6WlxYPMVajHpCBpuKTZkn4CXAN8AzgMuBP4acrxmVmZNDY2MmPGDE477TTWr19f1oTQ2tpKQ0MDw4YN83TYjJXSffQUcC9wVUT8qqD9VkmndnGOmVlJOqbDdnRdeTpstkoZaJ4aEQs6JQQAIuKSFGIysyri6bCVpZTlOP9UjkDMKktb1gFUDU+HrSyl3CmYmaXG02Eri5OCmWWqubmZ2tq9n573dNjsdJkUJH2mu59yBmlmQ1fW02E982lv3Y0p7Je8HgGcSK4UBeRqH92fZlBmNnAq/cE5yCWGxYsXA9DW1ta/D+sFz3zaV5d3ChHxlYj4CjAWOD4iPhsRnwVOIFewzsxsUPPMp32VMqYwEXi9YPt1cmszmpkNap75tK9SHl77PvBI8kRzAOeSW33NzGxQmzhxIhs27LtuSDXPfCql9lEzMB/4I/AyMD8ivp5yXGZVrA0/J1Eenvm0r1KnpNYCr0bENUC7pEkpxmRmVhYdM5/q6+uR5EKAlNB9JOlyYDq5WUg3AiOBpcDJ6YZmZpa+jplPDQ0NZZ35VKlKuVM4F5gNbAeIiBf4y3RVMzMbQkpJCq9HRJAbZEbSm9MNyczMslLK7KNbJH0H2F/S3wEfA65PNywzGwoGw4NztrdSqqReLWkW8Cq5cYV/jIi7U4/MzMzKrpSV166MiLsj4vMR8bmIuFvSlf25qKRPS1otaZWkmyXVSDpQ0t2SnkpeD+jPNczMrPdKGVOYVaTtzL5eUNKhwCXA9IiYAgwH5gALgeURMRlYnmybmVkZdVcl9ROSngSOlLSy4Oc54Ml+XncE8CZJI8g9A/ECcDawJNm/BDinn9cws15rww/OlU8lVmjtbkzh34D/BP4ve//Wvi0i/tDXC0bE7yRdDWwE/gz8PCJ+LungiNiUHLNJ0kHFzpfUBDRBdT+KbmaDW6VWaO2uSuorEbEeuAb4Q0RsiIgNwC5J7+zrBZOxgrOBScAhwJslzS31/IhoiYjpETG9rq6ur2GYmWWqUiu0ljKmcB1QuE7z9qStr84AnouILRGxC7gNeDfwoqRxAMnr5n5cw8ysolVqhdZSkoKSh9cAiIg9lPZ8Q1c2AjMk1UoScDqwltwiPvOSY+YBd/TjGmZmFa1S16YuJSk8K+kSSSOTn0uBZ/t6wYh4GLgVeIzcgPUwoAW4Apgl6SlyM56u6Os1zMyk0n/uuy/305tz+vJgXqFKrdBaym/8Hwe+DXyJXKmL5SQDvX0VEZcDl3dq3knursHMbMjrGExesGABO3fupL6+nubm5swrtJbyRPNmcs8RmJnZAMpqberulPJE89slLZe0KtmeKulL6YdmZmblVsqYwmLgMmAXQESsxHcOZjbg2vCDc9krJSnURsQjndp2pxGMmZllq5Sk8JKkw/nLegrnA5tSjcrMzDJRyuyji8lNGT1S0u+A54DqXcDUzGwIK2X20bPAGcmKa8MiYlv6YZmZWRZ6TAqSxpB7puA9QEh6EPhqRGxNOzgzs8FoMK84V8qYwg+ALcB5wPnJ+x+mE46ZmWWplDGFAyPiawXb/yTpnJTiMTOzDJVyp3CvpDmShiU/HwT+I+3AzMzKpw0/I5FTSlK4iNyCOzuTnx8An5G0TdKraQZnZmblVcrso/3KEYiZmWWvlNpHCzptD5fUucKpmZkNAaV0H50u6aeSxkk6BngI8N2DmdkQVEr30f+W9CFyC+LsAC6IiF+mHpmZmZVdKd1Hk4FLgR8D64EPS6rt9iQzMxuUSuk+uhP4x4i4CDgNeAr4r1SjMjOzTJTy8NpJEfEqQEQE8A1Jy9INy8zMslDKncKbJH1X0l0Ako4CTk03LDOzatFGJT04V0pSuAn4GTAu2f4t8Kn+XFTS/pJulfQbSWslvUvSgZLulvRU8npAf65hZma9V0pSGBsRtwB7ACJiN/BGP697DXBXRBwJTAPWAguB5RExGViebJuZWRmVkhS2J+WzO1ZemwG80tcLShpNrvvpuwAR8XpEvAycDSxJDlsCnNPXa5iZWd+UMtD8GWAZcLikXwJ15Epo99Vh5Mpv3yhpGvAouSmvB0fEJoCI2CTpoGInS2oCmgAmTpzYjzDMzKyzHu8UIuIxclNR302uON7REbGyH9ccARwPXBcRxwHb6UVXUUS0RMT0iJheV1fXjzDMzKyzUrqPiIjdEbE6IlZFxK5+XrMdaI+Ih5PtW8kliRcljQNIXjf38zpmZtZLJSWFgRQRvweel3RE0nQ6sIZcF9W8pG0ecEe5YzMzq3aljCmk4ZNAq6S/Ap4F5pNLULckVVk3An+bUWxmZlWry6Qg6fjuTkzGGvokIh4HphfZdXpfP9PMzPqvuzuFb3SzL4D3DnAsZmaWsS6TQkT8z3IGYmZm2StpTEHSFOAooKajLSK+l1ZQZmaWjR6TQrL05kxySeGnwJnAg4CTgpnZEFPKlNTzyQ0A/z4i5pOrVTQq1ajMzCwTpSSFP0fEHmB3UrdoM7lSFWZmNsSUMqawQtL+wGJydYr+BDySZlBmZpaNHpNCRPx98vZfk4V2Rvez9pGZmVWoHruPJC3veB8R6yNiZWGbmZkNHd090VwD1AJjk1XQlOwaDRxShtjMzKzMuus+uojcspuHAIUlLV4F/iXFmMzMLCPdPdF8DXCNpE9GxD+XMSYzM8tIKbOPviPpEnJLaAK0Ad8ZgHUVzMyswpSSFK4FRiavAB8GrgMuTCsoMzPLRncDzSMiYjdwYkRMK9h1j6Qn0g/NzMzKrbspqR0PqL0h6fCORkmHAW+kGpWZmWWiu+6jjimonwPulfRsst1AbqU0MzMbYrpLCnWSPpO8/w4wHNhOrnz2ccC9KcdmZmZl1l1SGA68hb/cMZBsA+yXWkRmZpaZ7pLCpoj4aloXljQcWAH8LiLOknQg8ENy3VPrgQ9GxB/Tur6Zme2ru4FmdbNvIFwKrC3YXggsj4jJwPJk28zMyqi7pHB6WheVNB74a+D6guazgSXJ+yXAOWld38zMiusyKUTEH1K87reALwB7CtoOjohNybU3AQcVO1FSk6QVklZs2bIlxRDNzKpPKSuvDShJZwGbI+LRvpwfES0RMT0iptfV1Q1wdGZm1a2UMhcD7WRgtqT/RW5662hJS4EXJY2LiE2SxpFb9tPMzMqo7HcKEXFZRIyPiAZgDnBPRMwFlgHzksPmAXeUOzYzs2pX9qTQjSuAWZKeAmYl22ZmVkZZdB/lRUQbuVLcRMRWUpzxZGZmPaukOwUzM8uYk4KZmeU5KZiZWZ6TgpmZ5TkpmJlZnpOCmZnlOSmYmVmek4KZmeU5KZiZWZ6TgpmZ5TkpmJlZnpOCmZnlOSmYmVmek4KZmeU5KZiZWZ6TgpmZ5TkpmJlZnpOCmZnlOSmYmVle2ZOCpAmS7pW0VtJqSZcm7QdKulvSU8nrAeWOzcys2mVxp7Ab+GxEvAOYAVws6ShgIbA8IiYDy5NtMzMro7InhYjYFBGPJe+3AWuBQ4GzgSXJYUuAc8odm5lZtct0TEFSA3Ac8DBwcERsglziAA7q4pwmSSskrdiyZUvZYjUzqwaZJQVJbwF+DHwqIl4t9byIaImI6RExva6uLr0AzcyqUCZJQdJIcgmhNSJuS5pflDQu2T8O2JxFbGZm1SyL2UcCvgusjYhvFuxaBsxL3s8D7ih3bGZm1W5EBtc8Gfgw8KSkx5O2fwCuAG6RtADYCPxtBrGZmVW1sieFiHgQUBe7Ty9nLGZmtjc/0WxmZnlOCmZmluekYGZmeU4KZmaW56RgZmZ5TgpmZpbnpGBmZnlOCmZmluekYGZmeU4KZmaW56RgZmZ5TgpmZpbnpGBmZnlOCmZmluekYGZmeU4KZmaW56RgZmZ5TgpmZpbnpGBmZnlOCmZmlldxSUHS+yWtk/S0pIVZx2NmVk0qKilIGg78C3AmcBRwgaSjso3KzKx6VFRSAE4Cno6IZyPideAHwNkZx2RmVjVGZB1AJ4cCzxdstwPvLDxAUhPQlGz+SdK6MsUGMBZ4qdSDpRQjKe+1e/W9U7h+rwzgtf29y3/tXvP37pP6rnZUWlIo9jVjr42IFqClPOHsTdKKiJiexbWz5O9dXfy9q1uldR+1AxMKtscDL2QUi5lZ1am0pPBfwGRJkyT9FTAHWJZxTGZmVaOiuo8iYrek/wP8DBgO3BARqzMOq1Am3VYVwN+7uvh7VzFFRM9HmZlZVai07iMzM8uQk4KZmeU5KZSgWktvSJog6V5JayWtlnRp1jGVi6Thkv5b0r9nHUs5Sdpf0q2SfpP8f39X1jGVg6RPJ3/GV0m6WVJN1jFlxUmhB1VeemM38NmIeAcwA7i4ir77pcDarIPIwDXAXRFxJDCNKvhvIOlQ4BJgekRMITfJZU62UWXHSaFnVVt6IyI2RcRjyftt5P6BODTbqNInaTzw18D1WcdSTpJGA6cC3wWIiNcj4uVMgyqfEcCbJI0Aaqni56OcFHpWrPTGkP+HsTNJDcBxwMMZh1IO3wK+AOzJOI5yOwzYAtyYdJ1dL+nNWQeVtoj4HXA1sBHYBLwSET/PNqrsOCn0rMfSG0OdpLcAPwY+FRGvZh1PmiSdBWyOiEezjiUDI4Djgesi4jhgOzDkx9AkHUDu7n8ScAjwZklzs40qO04KPavq0huSRpJLCK0RcVvW8ZTBycBsSevJdRW+V9LSbEMqm3agPSI67gZvJZckhrozgOciYktE7AJuA96dcUyZcVLoWdWW3pAkcv3LayPim1nHUw4RcVlEjI+IBnL/r++JiKr4rTEifg88L+mIpOl0YE2GIZXLRmCGpNrkz/zpVMEAe1cqqsxFJRoEpTfSdDLwYeBJSY8nbf8QET/NLiRL2SeB1uQXoGeB+RnHk7qIeFjSrcBj5Gbc/TdVXPLCZS7MzCzP3UdmZpbnpGBmZnlOCmZmluekYGZmeU4KZmaW56RgVUFSg6RVXey7vhIK/XUXo1m5+DkFq3oRcWHWMQwESSMiYnfWcdjg5jsFqyYjJC2RtDJZM6AWQFKbpOnJ+z9Japb0hKSHJB3c+UMkfVnSDcl5z0q6JGnf6zd9SZ+T9OWCa/w/Sfcn6xScKOk2SU9J+qcSYjxB0n2SHpX0M0njCj7365LuI1fu26xfnBSsmhwBtETEVOBV4O+LHPNm4KGImAbcD/xdF591JPA+cqXVL09qRPXk9Yg4FfhX4A7gYmAK8FFJY7qKMfnsfwbOj4gTgBuA5oLP3T8iTouIb5QQg1m3nBSsmjwfEb9M3i8F3lPkmNeBjtXWHgUauvis/4iInRHxErAZ2OeOooiOmllPAquT9Sp2kisn0VF0sViMR5BLHncn5Ua+RK4wY4cflnBts5J4TMGqSeeaLsVqvOyKv9R+eYOu/47sLHjfcdxu9v5Fq/OSjh3n7Ol0/p6C6xSLUeSSSFdLY27vot2s13ynYNVkYsGawxcADw7w578IHCRpjKRRwFl9+IxiMa4D6jraJY2UdPSARGzWiZOCVZO1wDxJK4EDgesG8sOTWvxfJbc63b8Dv+nDx+wTY7IM7PnAlZKeAB6niuv9W7pcJdXMzPJ8p2BmZnlOCmZmluekYGZmeU4KZmaW56RgZmZ5TgpmZpbnpGBmZnn/H8AcleX8P1XkAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Define the binning\n",
    "nbins = 10\n",
    "x = np.linspace(0, nbins - 1, nbins)\n",
    "\n",
    "# The background follows a linear shape\n",
    "b_yields = np.array([ (1 - i/2/nbins) for i in range(0, nbins) ])\n",
    "b_fracs = b_yields/np.sum(b_yields)\n",
    "\n",
    "# The signal shape is a peak\n",
    "s_fracs = np.zeros(nbins)\n",
    "s_fracs[3:7] = [ 0.1, 0.4, 0.4, 0.1 ]\n",
    "\n",
    "# Specify signal and background yields\n",
    "s = 70\n",
    "b = 1000\n",
    "s_and_b = s*s_fracs + b*b_fracs\n",
    "b_only = b*b_fracs\n",
    "\n",
    "# Generate a dataset\n",
    "np.random.seed(6) # make sure we always generate the same\n",
    "data = [ np.random.poisson(s*s_frac + b*b_frac) for s_frac, b_frac in zip(s_fracs, b_fracs) ]\n",
    "plt.bar(x, b_only, yerr=np.sqrt(b_only), color='b', label='b only')\n",
    "plt.scatter(x, data, zorder=10, color='k', label='data')\n",
    "plt.xlabel('bin number')\n",
    "plt.ylabel('Total expected yield')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now compute the best-fit $\\hat{s}$ and from it the value of $q_0$, and use this to compute the significance of the signal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathbf{\\color{red}{\\rightarrow\\text{Compute the values of } \\hat{s} \\text {, } q_0 \\text { and the significance}}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "exact hat{s}: -999\n",
      "q0 = -999 z = -999 p-value = -999\n"
     ]
    }
   ],
   "source": [
    "from scipy.optimize import minimize_scalar\n",
    "s_hat = -999 # replace -999 with the correct expression\n",
    "print('exact hat{s}:', s_hat)\n",
    "\n",
    "# Compute and print q0, z, and the p-value\n",
    "q0 =  -999 # replace -999 with the correct expression\n",
    "print('q0 =', q0, 'z =', -999, 'p-value =', -999) # replace the -999 with the correct expressions\n",
    "\n",
    "pass # here plot the best-fit s+b distribution together with the data"
   ]
  },
  {
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}