lecture_intro.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Statistical Methods in Python For Rare Event Searches\n",
    "\n",
    "*Summer Course at Canfranc Underground Laboratory (LSC) July 2019*\n",
    "\n",
    "### Lecture I\n",
    "\n",
    "## Introduction\n",
    "\n",
    "*Author: Jose A. Hernando*, June 2019\n",
    "\n",
    "*Instituto Galego de Altas Enerxías. Universidade de Santiago de Compostela, Spain.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Last Execution  Mon Jul  8 16:20:20 2019\n"
     ]
    }
   ],
   "source": [
    "import time\n",
    "print( ' Last Execution ', time.asctime())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**About**\n",
    "\n",
    "These lectures are about statistical method for rare event searches in Particle Physics using Python. They cover Hypothesis Testing and Confidence Intervals. They are based on the excellent lectures on statistics by Prosper [1], Cowan [2] and Cranmer[3] given at CERN Academic Training.\n",
    "\n",
    "We will use the Python scientic toolkits, Matplotlib, Numpy, Scipy, that are distributed with Anaconda Python [9]. \n",
    "\n",
    "A disclaimer: I am not a statistician not a programmer! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# general imports\n",
    "%matplotlib inline\n",
    "%reload_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "# numpy and matplotlib\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats       as stats\n",
    "import scipy.optimize     as optimize \n",
    "\n",
    "import htintro_examples  as htexam\n",
    "\n",
    "plt.style.context('seaborn-colorblind');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 1. Introduction\n",
    "\n",
    "We are making an experiment to discover a new particle. \n",
    "\n",
    "We are expecting few events if the particle exits in Nature (it is usually the case at any discovery!). This is a **rare search**.\n",
    "\n",
    "Rare events usually follow poissonian disbritubions. Unfortunately, as the statistics are nicely and friendly in the \"Gaussian domain\".\n",
    "\n",
    "And second,\n",
    "\n",
    "When we could clain an **observation** or a **discovery** of the new particle? And if not, what is the **limit** in a given observable (i.e. the half lifetime) we could impose?\n",
    "\n",
    "In fact, *what does it mean discovery, observation, a limit, a confidence interval?*\n",
    "And, *how do we compute them from data?*\n",
    "\n",
    "These are the question we try to answer in these lectures."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We will discovery that we need to cover three issues:\n",
    "- Hypothesis Testing\n",
    "- Confidence Interval\n",
    "- Regression \n",
    "\n",
    "They seem different but they are not!\n",
    "\n",
    "We start for the a **simple hypothesis**, $H_1$ (i.e BSM) that confronts the current hypothesis, $H_0$, (called null, i.e SM). \n",
    "\n",
    "Soon we realize that in general the alternative hypothesis depends on a **strength parameter** $\\mu$, that is, $H_1(\\mu)$. For example the half-lifetime of $\\beta\\beta0\\nu$ decay. This is called **composite hypothesis**.\n",
    "\n",
    "And then we are back to estimate the region of $\\mu$ \"compatible\" with data $x$, that is to define some **confidence intervals**, which in most cases implies to estimate the parameter, $\\hat{\\mu}$, that is **regression**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The starting point is a bifurcation: either we follow a **Bayes** or a **Frequentist** path.\n",
    "\n",
    "Bayesian compute *credible intervals* while Frequentist talk about *coverage intervals* (CI) and they do mean very different things! \n",
    "\n",
    "Be a Bayesian usually implies do **integration** (sometime complicated integrals!).\n",
    "\n",
    "Be a frequentist implies usually either do regression (fits!) or do **simulations**. \n",
    "But thanks to the current computer power, we can play the **frequentist game**!\n",
    "\n",
    "In this lectures we will use Python (in fact, scipy, numpy and mathplotlib) for our (mostly) frequentist journey!\n",
    "\n",
    "But let's start with refreshing some basic ideas about probability density functions, likelihoods and posterior probabilities."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 2. Probability and Bayes' theorem\n",
    "\n",
    "Nature is probabilistic and our measurements are too. \n",
    "\n",
    "The starting point of probability are the **Kolmogorov axioms**:\n",
    "\n",
    "* the probability for an event, $E$, is non negative $P(E)\\ge 0$,\n",
    "* the probability for the entire space of possibilities,$\\Omega$, is one, $P(\\Omega) = 1$,\n",
    "* the probability for disjoint events, $E_1,\\dots,E_n$, is additive, $P(E_1  \\dots  \\cup E_n) = \\sum_{i=1,n} P(E_i)$.\n",
    "\n",
    "From there we obtain the following corollaries:\n",
    "* $P(A \\cup B) = P(A) + P(B) - P(A \\cap B) $\n",
    "* $P(!E) = 1- P(E)$ \n",
    "where $!E$ is probability that $E$ does not happen."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Bayes' theorem\n",
    "\n",
    "The probability $P(A)$ that an event happens, for example that rolling the dice we get 3, is called **marginal probability**. The probability that an event $A$ happens if another one $B$ has happened is called **conditional probability**, $P(A|B)$. For example, the probability that rolling your dice is 3 if you know that the number obtained was odd, in that case $P(3|\\mathrm{odd})=1/3$. The Bayes' theorem relates both probabilities. \n",
    "\n",
    "If $A,B$ are two events with marginal probabilities $P(A),\\, P(B)$, the conditional probabilities $P(A|B), \\, P(B|A)$ are:\n",
    "\n",
    "$$\n",
    "P(A|B) \\, P(B) = P(B|A)\\,P(A)\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Example:** consider the experiment of rolling two dices, $A$ is the case when the sum of both dices is 6, and $B$ when one of them is 4. Verify the Baye's theorem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The probability to get six adding both dices is $P(A)= 5/36$, possible outcomes are $\\{(1,5),(2,4),(3,3),(4,2),(5,1)\\}$, the probability to get four when the addition is six is therefore $P(B|A)=2/5$. Now, the probability to get one dice with four is $P(B)=11/36$, and the probability that both add six if one is four is $2/11$, therefore:\n",
    "\n",
    "$$\n",
    "\\frac{11}{36} \\, \\frac{2}{11} = \\frac{2}{5} \\, \\frac{5}{36} = \\frac{1}{18}\n",
    "$$\n",
    "\n",
    "what verifies the Bayes' theorem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Exercise:** The probability of an individual to have HIV+ is 1 in 10000. Consider a medical test of HIV. If the patient is HIV+, the test is positive in 99 % of the times, but if the person is HIV- person, the test can wrongly be positive in 2 % of the times. After a test, a patient is diagnosed with HIV+, what is the probability that he has HIV +? Should the doctor repeat the test?\n",
    "\n",
    "**Exercise:** A neutrino experiment has a problem with the DAQ, and 5% of the runing time it does not observe interactions. The neutrinos are produced by a reactor nearby that operates 75% of the time. In this moment, the experiment does not observe neutrinos, what is the probability that the reactor is off?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 3. Frequentist vs Bayesians\n",
    "\n",
    "The is a great divide in statistics: frequentist vs Bayesians. \n",
    "\n",
    "**Frequentist are inductive**. *They compute the probability of an observation by repeating the same experiment many times*. They test is the data is compatible with the theory. \n",
    "\n",
    "They measure a **likelihood function**:\n",
    "$$\n",
    "\\mathcal{L} (x | \\mu)\n",
    "$$\n",
    "\n",
    "where $\\mu$ is the theory and $x$ the data. \n",
    "\n",
    "**Bayesians are deductive**. *They compute the probability using the Bayes' theorem. \n",
    "They measure the credibility of a theory based on the data, but they need to assign first a prior probability based in a 'reasonable' initial guess*. \n",
    "\n",
    "They measure a **posterior probability**:\n",
    "$$\n",
    "p(\\mu | x)\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "A conversation of a bayesian and a frequentist about rolling a dice with 6 faces!\n",
    "\n",
    "<img src=\"figs/bayes_cartoon.png\" width=\"300\">\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Most of physicist we are bayesians, but we think of ourselves as frequentists. \n",
    "\n",
    "Consider the discovery of the Higgs. If you are frequentist you claim: \"LHC data strongly disagrees with the hypothesis of the SM without the Higgs and agrees with a Higgs in the SM.\", while a Bayesian claims: \"We have discovered the Higgs!\"\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 4. Probability density functions\n",
    "\n",
    "\n",
    "In physics the variables that we measure are of probabilistic nature, therefore are **random variables** (rv). The probability distribution, $g(x)$ that follow a rv is called **probability density function**, pdf. If $x$ takes discrete values, the name is **probability mass function** (pmf). \n",
    "\n",
    "In most of the cases the pdfs depends on a parameter(s) $\\mu$, indicated as $g(x | \\mu)$.\n",
    "\n",
    "From the probability axions we have: \n",
    "\n",
    "$$\n",
    "\\int_{-\\infty}^{\\infty} g(x) \\mathrm{d}x = 1\n",
    "$$\n",
    "\n",
    "We give the cumutative density function (cdf), $F(x)$, as:\n",
    "\n",
    "$$\n",
    "F(x) = \\int_{-\\infty}^{x} g(x) \\mathrm{d}x\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "We call **marginal** pdf* when one or more rvs are integrated out, for example if we integrate $x_n$:\n",
    "$$\n",
    "g(x_1,\\dots,x_{n-1}) = \\int g({\\bf x}) \\mathrm{d}x_n\n",
    "$$\n",
    "\n",
    "We can \"project\" the pdf in one axis, that is get the marginal pdf for each variable $x_i$ individually:\n",
    "$$\n",
    "g(x_j) = \\int g({\\bf x}) \\, \\Pi_{i=1,n; i \\neq j} \\, \\mathrm{d}x_i\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Example:** Consider the case in which we try to observe a signal event in $n$ independent channels. In each one we expect $b_i$ background events and, if there is signal, $s_i$ events. What is the total PDF?\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "As the channels are $n$ independent, the pdf is the product of the pdf of each channel. In each channel we can observe $n_i$ events, the total observation is ${\\bf n} = (n_1,\\dots,n_n)$. The pdfs for the $b$, background only, and $b+s$, background and signal hypothesis are distributed with a Poisson:\n",
    "$$\n",
    "g({\\bf n}|s,b) = \\Pi_{i=1,n} \\frac{(s_i+b_i)^{n_i}}{n_i!} e^{-(s_i+b_i)}, \\;\\;\\;\n",
    "g({\\bf n}|b) = \\Pi_{i=1,n} \\frac{b_i^{n_i}}{n_i!} e^{-b_i}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 4. Likelihood and posterior  probability - An example\n",
    "\n",
    "Let's consider the basic measurement case seen by a frequentist and by a bayesian.\n",
    "\n",
    "We just measure a quantity $\\mu$ using a set of data $x$. We know that $x$ are gaussian distributed with known sigma, $\\sigma$. We want to estimate $\\hat{\\mu}$.\n",
    "\n",
    "For example, we have a set of measurement from a gaussian distribution with mean zero and sigma one. How we obtaine $\\hat{mu}$. \n",
    "\n",
    "Obvious, it is the average! And the uncertanty in the average is $\\sigma/\\sqrt{n}$, where $n$ is the number of events in our sample!\n",
    "\n",
    "But let's see how a frequentist and a bayesian proceed in detail."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Likelihood \n",
    "A frequentist evaluate the likelihood of its data.\n",
    "\n",
    "The **probability density function**, *pdf*, $g(x | \\mu)$is the probability to measure $x$, which depends via the parameter $\\mu$.\n",
    "\n",
    "A **likelihood** is the probability evaluted on a observation data, $x$, (on data!)\n",
    "\n",
    "For n, ${\\bf x}$, independent events that follow a pdf, $g(x | \\mu)$, the likelihood is the product of the likelihood of each measurement $x_i$:\n",
    "\n",
    "$$\n",
    "\\mathcal{L}({\\bf x} | \\mu) = \\Pi_{i=1}^n g(x_i | \\mu)\n",
    "$$\n",
    "\n",
    "\n",
    "The likelihood can be a very small number, for that reason, it is common to take the logarithm, and called **log likelihood**:\n",
    "\n",
    "$$\n",
    "\\ln \\mathcal{L}(x |\\mu) = \\sum_{i=1}^n \\ln \\left(g(x_i |\\mu) \\right)\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "\n",
    "Given a dat $x$, and with $\\mu$ inknown, the likelihood is a function that depends on the parameter $\\mu$, $\\mathcal{L}(\\mu)$. In most cases a frequentist estimate the best-parameters, $\\hat{\\mu}$ as the parameters that maximizes the likelihood (this is the maximum likelihood estimation).\n",
    "\n",
    "But for convenience, use instead:\n",
    "$$\n",
    "-2 \\log \\mathcal{L}(x |\\mu)\n",
    "$$\n",
    "\n",
    "In the following example, you can generate a n-size sample with $x$ values generated random from a normal gaussian. The likelihood function is drawn as a function $\\mu$. The best-estimate $\\hat{\\mu}$ is the value where the likelihood is maximum.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Explore:** Generate $n$ samples of a normal. Compute the likelihood as funciton of the mean (fix sigma of the gaussian to 1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mu mean : -0.06629666962817765 , mu std : 0.1\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbAAAAEKCAYAAABzHwA5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzt3XmcVNW16PHf6omGZp6UyW5UZJCZVkAZZYgSFcQ5OGtI7stokheT3CQmJt4k3htjyE3yNNFEIxKioogaFRFknmdBBJGWBkTmoRvoab0/9il6bqq6q+rUsL6fT326+tSpc1Z1V/eqvc/ea4uqYowxxsSbFL8DMMYYY+rDEpgxxpi4ZAnMGGNMXLIEZowxJi5ZAjPGGBOXLIEZY4yJS5bAjDHGxCVLYMYYY+KSJTBjjDFxKc3vAIKRkpKijRs39jsMY4yJK4WFhaqqCdtQiYsE1rhxYwoKCvwOwxhj4oqInPI7hkhK2MxsjDEmsVkCM8YYE5csgRljjIlLcXENzBgTe4qLi8nPz+f06dN+h5L0MjMz6dy5M+np6X6HElWWwIwx9ZKfn0+zZs3IyclBRPwOJ2mpKocOHSI/P5+uXbv6HU5URawLUUSeEZHPRWRzhW2tRWSuiGz3vraK1PmNMZF1+vRp2rRpY8nLZyJCmzZtkrIlHMlrYH8Hrq6y7QfAPFXtBszzvjfGxClLXhGkGvSuyfp7iFgCU9WFwOEqmycCz3r3nwUmRer8xhgT106dggMHIAlbVsGK9ijE81R1H4D3tX2Uz2+MMfGhoACKi6GszO9IYlbMDqMXkakislpEVpeUlPgdjjExJzB4ouotJyfH79BMQxUXu5sIZGb6HU3MinYC2y8iHQC8r5/XtqOqPqWquaqam5ZmgyWNqSovLw9VrXbLy8vzO7SYsXv3bkaPHk3Pnj259NJL+f3vf19tn6ZNm0bk3KdOnWLkyJGUlpbWud+oUaPYtWtX5Y2Fhe5r48aQkkJRUREjRozAPsxXFu0E9hpwt3f/bmB2lM9vjEkiaWlp/Pa3v2Xr1q0sX76cP/7xj2zZsiXs55k5cyZTp07lscceO7vtmWeeYfLkyaSmpoZ2sLKy8gSWlQVARkYGY8aMYebMmeEKOSFErGkjIjOAUUBbEckHHgZ+DfxLRO4HPgVujtT5jTFR9LOfxeRxO3ToQIcOHQBo1qwZPXv2ZM+ePfTq1avG/R9//HGeeeYZAB544AG+/e1vA/CLX/yC6dOn06VLF9q2bcugQYP43ve+B8C0adM4evQo9957Lw8++CDf//73AZg+fTovvPDC2WOPHj2aH/3oR4wbN44f//jHHD9+nGnTplU6/6hRo3jyiSfo3r49h06cYOT48Wze7GYiTZo0iR/+8IdMmTKlQT+TRBKxBKaqt9fy0JhIndMYY2qza9cu1q1bx+DBg2t8fM2aNfztb39jxYoVqCqDBw8+2wX48ssvs27dOkpKShg4cCCDBg0CoLCwkEceeYTt27dTUlLCc889B0BRURE7d+6sdD3y5z//OT/96U/5/PPPWbduHa+99lq1GHbs2EG3jh2hpISNH39Mnz59zj7Wu3dvVq1aFcafSPyzi0vGmIaLVAssCGPHjuWzzz6rtv3RRx9l4sSJAJw8eZIbb7yRJ554gubNm9d4nMWLF3PDDTeQ5XXbTZ48mUWLFlFWVsbEiRMJrEl43XXXnX3O/Pnzyc7OplUrV5OhXbt2ABw8eJCWLVtWOv6IESNQVR5//HEWLFhQrWsxLy+PTh07klJSAiJs/Ogj+vbte/bx1NRUMjIyOHHiBM2aNQvpZ5SoLIEZY+Lau+++W+fjxcXF3HjjjUyZMoXJkyfXup/WMnG4tu0ACxcuZPz48dW2N27cuFpljE2bNrFv3z7atm1bYwJav349fXv2DByANWvXcuutt1ba58yZM2TaqMSzYnYYvTHGNJSqcv/999OzZ0++853v1LnviBEjePXVVyksLKSgoIBXXnmF4cOHM2zYMObMmcPp06c5efIkb7zxxtnnLF26lMGDB3PgwAEAjh8/zpYtW2jVqhWlpaVnk9i+ffuYMmUKs2fPJisri7fffrva+TesX89pb+He7Xv2MHv27EpdiIcOHaJdu3ZJV7C3LpbAjDEJa8mSJfzjH//gvffeo3///vTv358333yzxn0HDhzIPffcw+WXX87gwYN54IEHGDBgAJdddhnXX389/fr1Y/LkyeTm5tKiRQsWLFhAeno6jz76KBdddBHXXHMNzz//PD169ABg/PjxLF68mMLCQiZPnsxvf/tbevbsyU9+8hN+VkOX6/q1aykrK6PfuHE88utf07NnT5599tmzj8+fP58JEyZE5OcUt2qaRxJrtyZNmqgxpjL35xv89nDbsmVLVM4TC06cOKGqqgUFBTpo0CBds2bNOZ+zdu1aveOOO86538iRI/WTTz7Ri3Jy9Pi2baonT9a43w033KAffvhhrcep6fcBFGgM/A+P1M1aYMYYcw5Tp06lf//+DBw4kBtvvJGBAwee8zkDBgxg9OjR55zIDHDi8GFSRNy1MW+wSEVFRUVMmjSJ7t271yv+2oSyaog400Rkh4hsFJGBFZ5zt7f/dhG5u6ZzRYJoCBWP/ZKVlaUFXt+wMcYRkRoHGNS2Pdy2bt1Kz8CgA1Nvf//735k0ejQt09PdxOUWLep1nJp+HyJSqKpZtT1HREYAJ4HnVLW3t+0x4LCq/lpEfgC0UtWHRGQC8A1gAjAY+L2qDhaR1sBqIBdQYA0wSFWP1OuFhMBaYMYY46N77rqLlhkZ7psmTaJ6bg1t1ZCJuESnqrocaOmVBPwCMFdVD3tJay7Vl9KKCEtgxhjjp1On3NpfGRkQGyMMa1s1pBOwu8J++d622rZHXGInsLIyOHjQ7yiMMaZmquV1DyPT+koLrOrh3aY24Fg1rZqpdWyPuMSdyHziBPz1r1BaCg8+CKEW1DTGmEgLLJuSklLj4I0wKFHV3BCfs19EOqjqviqrhuQDXSrs1xnY620fVWX7gvqFG5rEbYE1bQqNGsHJk7B1q9/RGGNMdYHBaU2auLW/YkNtq4a8BtzljUYcAhzzuhjfBsaLSCtvxOJ4b1vEJW4CE4HLL3f3rQCmMSbWlJa6618Q9cEbAd6qIcuA7iKS760U8mtgnIhsB8Z53wO8CewEdgB/Af4PgKoeBn4BrPJuj3jbIi5xuxAB+vSBuXMhLw/274fzzvM7ImOMcQLJq1Ej8GnRXg1h1RBvYvTXajnOM8AzYQwtKInbAgP3xujXz923VpgxJlaolncfZtU6TcucQ2InMIDLLnNfN26EM2f8jcUYExGvvPIKIsKHH34Y1P6nTp06u9ZXXYqKihgxYgQlJSX1ju2+++6jffv29O7du3zjmTOuCzE11X3QNvWS+AmsXTvIyYGiItiwwe9ojDERMGPGDHJzc/nnP/9Z7bGZM2cydepUHnvssbPbnnnmGSZPnlxtTa6qMjIyGDNmDDNnzqx3bPfccw9vvfVW5Y2xOXgj7iR+AoPyVtiqVa7pboxJGCdPnuT999/n6aefZsaMGZUemzZtGtu2bePee+9l1qxZZ7dPnz797GKXAKNGjWLbtm2AW7akYmtp0qRJTJ8+vd7xjRgxgtatW5dvKCkp7w3yafBGokjsQRwBPXpAs2Zw4ADs2gVdu/odkTEJ59Ynl1Xbdm3fDtw5NIdTRaXc87eV1R6/aVBnbs7twuGCIv7j+TWVHpv5laFBnffVV19l7Nix9O3bl6ysLNauXcvAgQMpLCzkkUceYfv27ZSUlPDcc88Brltw586d5OTknD3Gjh076NatGwAbN26stA5X7969WVXlGvrw4cM5ceJEtVj+53/+h7Fjx9YdcGDicuPGNj+1gZIjgaWmQm4uzJ8PK1daAjMmgcyYMYOpU12BiVtuuYUZM2YwcOBA5s+fT3Z2Nq1atQKgXbt2ABw8eJCWLVuefX5eXh6dOnUiJcV1SG3cuJG+ffuefTw1NZWMjAxOnDhxdiXlRYsW1S/YsjIbvBFGyZHAAAYNgoUL4cMP4ehRqPAGNsY0XF0tpsYZqXU+3jorI+gWV0WHDh1i5cqVZ7sHb731VkaOHMljjz3GwoULGT9+fPVYGjc+u1IywPr16yslrDVr1nDrrbdWes6ZM2fIzMw8+329W2CBuofp6a72oWmQ5ElgTZtCr16waROsXg3nauYbY2LeSy+9xIQJE2jkjeTr2rUr559/PosXL2bp0qV897vf5cCBA7Rr147jx4+Tn59Pr169KC0t5fTp02RmZrJhw4azCW379u3Mnj2bX/7yl2fPcejQIdq1a0d6hUK79W6BWesrrJJjEEfA4MHu69q1rv6YMSauzZgxgzlz5pCTk3P2tnXrVqZOnUp6ejqPPvooF110Eddccw3PP/88PXr0AGD8+PEsXrwYcC2wsrIy+vXrxyOPPELPnj159tlnz55j/vz5TJgwod4x3n777QwdOpRt27bRuV8/nv7nPyNV9zDpJE8LDKBTJ+jYEfbuhc2bYcAAvyMyxjTAggUL6vW8r3/96zz++OOMHTuWjRs3sm7durPXt6p64YUX+NWvflXvGM+OjDx8GE6fdr1BNnQ+LJKrBVaxPuLKlTak3pgkNWDAAEaPHs2xY8dISUmpNXkVFRUxadIkunfv3rATlpa65AU2dD6MkiuBAfTu7d5A+/bB7t3n3t8Yk5Duu+8+WrRowUcffVTrPhkZGdx1110NP1ng2ldmpm91DxNR8iWwtDQ3pB5gxQp/YzHGJL6Ki1ba4I2wSr4EBq4yR0qKWyfs2DG/ozHGJLJTp9z8Lxs6H3bJmcCaNYNLL3VvKqtSb0y9qV1HrpuqW1QXXOsrQoM3kvX3kJwJDMqH1K9ZY0PqjamHzMxMDh06lLT/PINSVORqH6akRGzovKpy6NChShOtk0XyXk3s3Nnd8vPdUiuDBvkdkTFxpXPnzuTn53PgwAG/Q4ldBQXuA3JmJhw5ErHTZGZm0rlz54gdP1YlbwID1wrLz4fly2HgQJubYUwI0tPT6Wp1RWt35AhMm+b+rzz4oLt0YcLKly5EEXlQRD4Qkc0iMkNE/Gn79upVXqV+505fQjDGJKjAXNPevS15RUjUE5iIdAK+CeSqam8gFbgt2nEArkp9YGLz8uW+hGCMSUBnzriSdQBDQy9SbILj1yCONKCxiKQBTYC9PsXhrn2lp8P27a4lZowxDbVunUti2dnQoYPf0SSsqCcwVd0D/A/wKbAPOKaq70Q7jrOaNIF+/dx9a4UZYxqqrKy8SIK1viLKjy7EVsBEoCvQEcgSkTtq2G+qiKwWkdUlJSWRDWrIEPd1w4byGfPGGFMf27a5ARytWsEll/gdTULzowtxLPCJqh5Q1WJgFnBF1Z1U9SlVzVXV3LRI1w5r29a90UpK3FphxhhTX8uWua9Dhrj5XyZi/PjpfgoMEZEmIiLAGGCrD3FUFmiFrVzpEpkxxoRqzx749FNo1Aj69/c7moTnxzWwFcBLwFpgkxfDU9GOo5quXeG881zZl82b/Y7GGBOPAq2v3FyXxExE+dK+VdWHVbWHqvZW1TtV9YwfcVQiAld4PZlLl9paYcaY0Bw9Ch984LoNA6XqTERZB21FgQmHn38OH3/sdzTGmHiyfHn5xOXmzf2OJilYAqsoNbX8WtjSpf7GYoyJH6dOlU9cvqLamDQTIZbAqho0yK3Zs3OnW7XZGGPOZc0aV3n+wgvh/PP9jiZpWAKrKjOzvDJ94IKsMcbUprS0fOKytb6iyhJYTQYPdhdiN2+2FZuNMXXbtAlOnID27eGii/yOJqlYAqtJy5blKzZbeSljTG1UYckSd/+KK2xJpiizBFabQFfAmjXuAq0xxlQVKALevDn06eN3NCGraWkrEekqIitEZLuIzBSRDG/fRt73O7zHc/yN3hJY7Tp0cN0BRUWwapXf0RhjYlGg9TV0qBvFHEfqWNrqN8DvVLUbcAS433vK/cARVb0Y+J23n68sgdXlyivd1xUr3LLgxhgTsHs35OW5gV8DB/odTX1VXdpqH3AVrloSwLPAJO/+RO97vMfHeOUAfWMJrC5du0LHjlBQ4CrVG2NMQKD1ddllcVk2qqalrYA1wFFVDRSEzQc6efc7Abu955Z4+7eJZsxVWQKri0h5K2zpUjeowxhjDh50y6akpcV62ai0wLJU3m1q4IGalrYCrqnhGIG6ejW1tnytuRfhdUoSQM+e0Lo1HD4MW7e60YnGmOS2eLEbgdivHzRt6nc0dSlR1dxaHju7tBWAiASWtmopImleK6szsNfbPx/oAuR7XY4tgMMRjf4crAV2Likp5SMSFy2yIr/GJLtjx2Djxso9NPGppqWttgDzgZu8fe4GZnv3X/O+x3v8PVV//yFaAgtG//7uU9Znn8GOHX5HY4zxU+BywqWXut6ZOFXH0lYPAd8RkR24a1xPe095Gmjjbf8O8IOoB12F+JxAg5KVlaUFBQX+BrFkCcydC9nZcO+9/sZiDCAi1PT3W9t2EwYFBfDEE25U8n/8h1tDMIaJSKGqZvkdR6RYCyxYubluuGxenltx1RiTfAJTai65JOaTVzKwBBasRo3KRxstWuRvLMZEQU5ODiJS7ZaTk+N3aP44cwZWrnT3hw3zNxYDWAILzeDBkJ7uysd89pnf0RgTUXl5eahqtVteXp7foflj1So4fdpdRrjgAr+jMVgCC02TJq4rEawVZkwyKS4uX15p+HB/YzFnWQILVaDm2ZYtroinMSbxrV3rBnB07GhLpsQQS2Chat4cBgxw88GsFWZM4ispKS8bNWKELZkSQyyB1cewYW6C86ZNrkKHMSZxbdgAx4+7UYfdu/sdjanAElh9tGzpSsioupIyxpjEVFpa3tMyfLi1viJARLJEpF5r0VgCq69hw9ybef16OHrU72iMMZGwaZP7+27TBnr18juahCAiKSLyJRF5Q0Q+Bz4E9nkLa/63iHQL9liWwOqrTRvo3duVlAn0jxtjEkdZGSxc6O6PGOEuG5hwmA9cBPwQOF9Vu6hqe2A4sBz4tYjcEcyBrBp9Q4wYAZs3uxFKw4ZBixZ+R2SMCZfANe7WraFPH7+jSSRjVbXaCsGqehh4GXhZRNKDOZB9pGiIdu1cQc/SUrsWZkwisdZXxNSUvOqzD1gCa7jAsNq1a91IJWNM/Nu8GQ4dcq2vvn39jiYhicgmEZkuIg+JyDUi0llE/jOUY1gCa6j27a0VZkwiqdj6Gj7cWl+RMxL4C3AKuA3YDHwxlAPYbyYcAq2wNWusFWZMvNu8GQ4ehFatrPUVQap6WFUXqOo0Vb0buAzYHsoxLIGFQ/v2bohtxTkjxpj4U1YG77/v7o8Y4crGmYioOlxeVbcDIX1isAQWLiNHll8LO3bM72iMMfWxcWP5ta9+/fyOJtE9JSKfisgyEXlSRJ4FNotIk2APYAksXNq3d/PCSkvL+8+NMfGjtLS89TVypF37ijBVHa2qFwC3Am8AO4AmwAYR+TCYY/gyD0xEWgJ/BXoDCtynqsv8iCWsRo50/efr1rl5Ya1a+R2RMSZYGzbAkSPQtq3N+4ogEalpMbX13u1ZQIBmItJcVescVODXRObfA2+p6k0ikoHLuvGvbVt30XfDBvdJbtIkvyMyxgSjYutr1ChrfUXWs7iGS9XCklrl+78Dz9V1oKgnMBFpDowA7gFQ1SKgKNpxRMzIkW4G/4YNrhXWtq3fERljziVw7TowIMtEjKqODtex/PiYcSFwAPibiKwTkb+KSJYPcURG69bQv7+rVB/4RGeMiV3FxeXXra31FTUi8kRDjxHyb8qrJNy8AedMAwYCf1bVAUAB8IMazjNVRFaLyOqSkpIGnM4HI0e64bebN8P+/X5HY4ypy6pVcOIEdOgAPXv6HU0yOSkicwINGBEZLyIhVUYPKoGJyAsi0tw70RZgm4j839DjBSAfyFfVFd73L+ESWiWq+pSq5qpqblpanNUcbtECcnNdK+y99/yOxhhTmzNnyudujhlj631Fkar+GJgBLBCRxcB3qaExU5dgW2C9vNEgk4A3gQuAO0M5UYCqfgbsFpHA0qZjcEkxsQwfDhkZsG0b7N7tdzTGmJosWwanTkF2Nlx0kd/RJBURGQN8GdcL1w74pqqGVAki2ASW7pW3nwTMDrZScB2+AUwXkY1Af+C/Gni82NO0KQwZ4u7Pm+daY8aY2FFYCEuXuvtXXWWtr+j7T+AnqjoKuAmYKSJXhXKAYPvmngR2ARuAhSKSDdS73ISqrgdy6/v8uHHFFa5/fdcu2LnTPuEZE0sWLYKiIrj4YtcCM1GlqldVuL9JRK7BrQd2RbDHCLYFNkdVO6nqBFVV4FPgvpCiTUaZmXDlle6+tcKMiR3HjsHKle7+mDH+xpJkRGpu6qrqPtwlpVr3qSrYBPZylRMp8M8gn5vcBg+GZs1g717YkniX+oyJS/Pnu8nLvXu70YcmmuaLyDeqVuTwiloM9Woi3h3MgersQhSRHsClQAsRmVzhoeZAZmgxJ6n0dDe3ZM4c1wrr0cMqXBvjp88/d4UGUlLctS8TbVfjevBmiEhX4Cgun6QC7wC/8y4zndO5roF1B64FWgLXVdh+Ajd6xARjwAA32ungQTfj/7LL/I7ImOQV6M6/7DJXeMBElaqeBv4E/MkbHNgWOKWqR0M9Vp0JTFVnA7NFZGhCFNv1S0qK62efORMWLHDLNGRk+B2VMcnn00/d1JaMDLfel/GVN6J9H7geP1UNqgp9QLDXwHaIyI9E5CkReSZwCzXYpNajB3TuDAUFrjVmjIkuVZg7190fOtRNdTG+EJH/KyJLReTiCpv3iMhXQzlOsAlsNtACeBe3bkvgZoIlAuPGuftLlsDJk/7GY0yy2brVFRXIynJTXAwi0lJEXhKRD0Vkq4gMFZHWIjJXRLZ7X1t5+4qITBORHSKyUUSqVVAKwcXAg7jrXwCo6gkqX6o6p2ATWBNVfUhV/6WqLwduoZzI4OaadO/u5p4sWOB3NMYkj9JSePddd3/UKGjUyNdwYkhgaaseQD9gK66c0zxV7QbMo7y80zVAN+82FfhzA847DxgHnC2KISJtgStDOUiwCex1EZkQyoFNLcaOddfE1qyBAwf8jsaY5LB6NRw+7JY3GtiQhkPiqLC01dPglrbyBlJMxK3Zhfc1sLDhROA5dZYDLUWkXnMQVPVfuGIYO0RklYg8ipvAvC2U4wSbwL6FS2KnReS4iJwQkTpXyjS1aNcOBg2q3B9vjImc06fLlzYaO9amsZSrbWmr87xJxYHJxe29/TsBFQu75nvb6kVV/4Crq/swbgj993Aj3IMWVAJT1WaqmqKqmara3Pu+IUuqJLdAF8ZHH7kSU8aYyFm82NU9DHThJ5e0wLJU3m1qxccIYmmrCmqqjtGg8kKqekpV31TVH6jqCOAXoTw/qFqIXlmPKUBXVf2FiHQBOqjqytBDNmRludWa582Dd96BqVNtET1jIuHIkfJRv+PHJ2PB3hJVra3ubE1LW/0A2C8iHVR1n9dF+HmF/btUeH5nYG99AxOR79Sw+ZiIHAt2InOw/zX/BAwFvuR9fxL4Y5DPNTUZMgSaN4fPPnNVAYwx4ffuu24AR9++0KnevV0JqY6lrV6jvJTT3bhR6Hjb7/JGIw4BjgW6GuspF/gqrhuyE25gyCjgLyLy/WAOEGw1+sGqOlBE1gGo6hGvbpWpr/R01x8/a5ZrifXqZSOjjAmnTz+FDz4o/1szNQksbZUB7ATuxTVs/iUi9+MKt9/s7fsmMAHYARR6+zZEG2Cgqp4EEJGHca3AEcAa4LFzHSDYBFYsIql4/Z0i0g4oq0/EpoI+fVxF7Px8109vVbGNCQ9VeOstd/+KK1xvh6mmjqWtqv0z8oq4fy2Mp78AKKrwfTGQraqnRORMMAcItgtxGvAK0N4b7riYRFyEMtpE4Oqr3f1ly+BoyKXA4l5OTg4iEvQtJycnrs9romTjRrcCRLNm5UsamVjzArBcRB4WkZ8BS3AFfrNwXZnnJBrkGlVeZfoxuJEo81R1a71CroesrCwtKCiI1umi7+WXYdMmuPRSuPnmc++fQESEYN+D9dk/1s4bTrXFFOmfUSz+LCopKoI//AFOnIAbbnC1R5OUiBSqapbfcdRGRAYBw3B5ZbGqrg7l+UG1wETk90BrVf2jqv5vNJNXUhg71vXTf/AB5OX5HY0x8W3RIpe8OnVygzdMLCvBXY4qoUJVjmAF24W4FvixVwPrv0WktmGZpj5atCjv5njzTSizy4vG1Mvhw7B0qbt/zTXJOGw+bojIt4DpuOVU2gPPi8g3QjlGsBOZn1XVCcDlwEfAb0Rke4jxmrpceSW0bAn797syU8aY0L39ths237+/W/3BxLL7cSPcH1bVnwJDCHGdyVBnz14M9ABygJDWbTHnkJ4OX/iCu//ee65ygDEmeDt2lK/1ZSN644EApRW+L6Xmah+1CvYa2G9E5CPgEWAzMEhVQyp7b4LQowdceCGcOgXz5/sdjTHxo7S0fNj8yJFu9KGJdX8DVojIz7xRiMvxCgsHK9gWWB7wS2Cpqv4NaC4il4dyIhOEwLD6lBRXPXtfQya5G5NEli2DgwehTRsYPNjvaEwQVPVx4D7gMHAEuFdVnwjlGMEmsD7AYOB27/sTWCmpyGjf3v0BqsIbb7ivxpjaHTtWXm1+wgRIC7Y+g/Gbqq5R1Wmq+ntVXRfq84NNYINV9WvAae+kRwArJRUpo0a5LpD8fFgX8u/UmOTy9ttQXOzKsV10kd/RmHMILMdVwy3kZbqCTWBWSiqaGjVylbPBFSO1AR3G1Ozjj2HLlsqDoExMCyzHVcMt5GW6rJRUrOrdG7p2dclr3jy/ozEm9pSUuHmT4AZutGjhbzwm6oKdBzYd+D7wK2AfMElVX4xkYElPxPXnp6bC2rWuO9EYU27JEjh0CNq2haFD/Y7G+CDoeWCq+qGVkoqydu3cH6YqzJnjhgobY1ziWrTI3b/2WvdBz8QNEekhImNEpGmV7VeHchxbBjjWjRwJrVq5Ch0rVpx7f2MSXWCEbkmJq7hhKwXEFRH5Jm6RzG8Am0VkYoWHQ7o0ZQks1qWnwxdL3uFTAAAYD0lEQVS/6O7Pn5+US64YU8mmTbBzJzRuXD7YycSTL+OKYUzCrcD8E68uIkSiEofx2cUXu0EdxcU2N8wkt1On3LB5cMmrSRN/4zH1kRpYhVlVd+GS2DUi8jjxksBEJFVE1onI637FEFeuvhoyM2H7djds2Jhk9M47UFAA2dmu+9DEo89E5Owvz0tm1+Kq0vcJ5UB+tsC+BdhgkGA1berWDQM3dPjUKX/jMSbadu50E/vT0uC662yplPh1F/BZxQ2qWqKqdwEjQjmQLwlMRDoDXwT+6sf549agQe6TZ0FBeTeKMcmguNiNxAU3sKltW3/jMfWmqvmqejaBich1FR5bEsqx/GqBPYGbV2bVPEIhAtdf7z6Brl/vqhAYkwwWLIAjR+C88+CKK/yOxoTXo/V9YtQTmIhcC3yuqnWu2igiU0VktYisLikpiVJ0caBNG/cJFNwn0qKiqIeQk5ODiFS75dQynLm2/UWE7OxsX2IycWTvXldtPvABzuZ8JZp69wWLRnlEm4j8CrgTKAEygebALFW9o7bnZGVlaUFBQZQijAOlpfCXv8Bnn8GQIW6ARxSJCDW9b0Ld7ue5Q40pnK8hXCL9847G7/OcSkrgqafg88/dpH6rdxgSESlU1Sy/46iLiGxU1b71eW7UW2Cq+kNV7ayqOcBtwHt1JS9Tg9RUmDjRrRu2YgXk5fkdkTGRsXChS16tW8NVV/kdjYkxNg8sXnXoAMOHuzlhs2e7i9zGJJJ9+2DxYtd1OHGim9RvTAW+JjBVXaCq1/oZQ1wbMcJd1D582CrWm8RSWuo+mJWVweWXu9G3JlHtr+8TrQUWz1JTYdIk60o0ief999013latYMwYv6MxEaSq4+r7XEtg8a5DBxg2zHUlvvoqnDnjd0TGNMyePeVdh5MmQYYt/m5qZgksEYwcCeef7+bJvPOO39EYU3/FxfDKK67rcOhQ6zo0dbIElghSU2HyZPd1zRr46CO/IzKmft59Fw4edGvh2ahDcw6WwBJF+/bl1wpeew0KC/2Nx5hQ7dzpruWmpMANN7iKM8bUwRJYIhkyxHW5nDzpqnTE2MRbY2pVWOi6DsF1iXfs6G88Ji5YAkskKSnuonejRrB1q6vcbUysU4XXX4cTJ6BLFze/0ZggWAJLNK1ala/g/O9/u+sJxsSy9evdGneNGrlruSn2bymaqq7NKCJdRWSFiGwXkZkikuFtb+R9v8N7PMfPuMESWGLq29fdioth1iw3KdSYWHT4sPugBTBhgvsAZqKt6tqMvwF+p6rdgCPA/d72+4Ejqnox8DtvP19ZAktUEyZAy5aukvd77/kdjTHVlZTAiy+6FRV693YfukxUVV2bUUQEuAp4ydvlWWCSd3+i9z3e42O8/X1jCSxRZWbCjTe67pglS2DHDr8jMqayd9919Q5btoRrr7UVlv1RdW3GNsBRVQ2sYZUPdPLudwJ2g1tBGTjm7e8bS2CJrEsXGD3a3Z81y10kNyYWbNsGy5e7D1g33eQ+cJlISAusq+jdpgYeqGVtxpo+RWgQj/nCJlokumHDYNcut3rzrFlw5512kdz46/hxV6gX3NzFzp39jSexlahqbi2PXQlcLyITKF+b8QmgpYikea2szsBeb/98oAuQLyJpQAvgcESjPwf7T5boRNyk0KZN4ZNP3PpKxviltNRd9yoshIsvhiuu8DuipFXL2oxTgPnATd5udwPepw1e877He/w99XmVV0tgyaBpUzc8WcRV+f74Y78jMslq3jzYvRuaN3cfrOy6Vyx6CPiOiOzAXeN62tv+NNDG2/4d4Ac+xXeWxNoy6TXJysrSgoICv8OIf++/D/PnQ5Mm8JWvQIsW9TpMqEvNh3MJ+nCdO9SYwvkawiXSP++wH3/rVpg503Vh33MPXHBBg2M0dRORQlXN8juOSLEWWDIZMcJ12xQWum6cc8wPy8nJQUSq3bKTsEJ4dnZ2jT+LUG85OTk1Hr+2n3VdN79+D6H+LHJyctx8r8B1r3Hj4IILan3Ntf2MjKnKWmDJprAQnnwSjh1zK91OmFDrruFqpSRCCyxc/PwZ+XWcDBGK/vQn2L8fevaEW24BL1nF0u8mEVkLzCSWJk3g5pvd0isrV7oyPsZEiirXg0tebdrAxIl23cuEjSWwZNS5c3nL6/XXXbUOYyJh+XJ6g1tV+bbbbL6XCStLYMlq0CDIzXXlfGbOBOuiNeH2yScwd667f8MNbpFKY8LIElgyu/pqV63j2DH417+s6K8Jn8OH3XuqrIxF4K59GRNmlsCSWVqau6DerBnk5cGbb9oimKbhTp+GGTPg1Cm45BLm+x2PSViWwJJds2bu2kRaGqxZ45Z0N6a+ysrg5ZfhwAHXZXjjjf4WyzMJzRKYgU6d3ErOAG+/Ddu3+xuPiV9z57r3T+PG8KUvuUUqjYkQS2DG6d0bRo50XYgvveSGPRsTilWrYNkyV2njlltscUoTcZbATLlRo1wiO3MGpk+nmd/xmPixfbu7hgpw/fXQtau/8ZikYAnMlBNxXYkXXADHj/MlcMnMmLp89pkrTabqypX17+93RCZJWAIzlaWluUEdbdpwPgRVM9EksaNHYfp0KCqCPn3KF1A1JgosgZnqmjSBKVMoANixwxVhteH1pqrCQvjHP9xK3zk5VibKRJ0lMFOz1q15AVwJoI0b4Z13LImZckVFruV16BCcf375VAxjosgSmKnVXoBbb3WFf5ctg6VL/Q7JxILSUldlY88eaNkSpkyxGofGF5bATN0uuqh8jtjcuW6ys0legYnKO3a4ruY77nCT4Y3xQdQTmIh0EZH5IrJVRD4QkW9FOwYToj59Klev37jR33iMP1RhzhzYssW1uO68E9q29Tsqk8T8aIGVAN9V1Z7AEOBrItLLhzhMKC6/HMaOdf/EXn3VLQ9vkocqvPUWrFsH6emuykaHDn5HZZJc1BOYqu5T1bXe/RPAVqBTtOMw9TBsmJvnU1bmqnV89JHfEZloUHUlxlascNdDb7vNzRU0xme+XgMTkRxgAFCtgqyITBWR1SKyuqSkJGox5eTkIN5y5xVvOTk5UYshpo0eDUOHugv5M2dWSmLZ2dk1/uyys7PDdvpQzxGNmOJFbe/t2m7Z2dkueb3zDixf7pLXrbe666LGxABRn4ZGi0hT4H3gUVWdVde+WVlZWhClBRdFhJp+JrVtT2S1vubAJ/KK/9QuuST6AcahaLy/wnYOVTdwZ+lS93u+5Rbo3j3i8STj31qkiEihqmb5HUek+NICE5F04GVg+rmSl4lBIvCFL8CQIeUtMbsmllhU4d//dskrJQVuvrleycuYSPJjFKIATwNbVfXxaJ/fhEkgiQW6E1980UYnJoqyMnjtNVi5sryF3aOH31EZU40fLbArgTuBq0RkvXeb4EMcpqFEYPz48oEdr7wCq1f7HZVpiNJSmDWr8mhDa3mZGBX12i+quhiwgmmJQgSuusotXDh3rpsnVlgIw4dbXbx4U1TkKmzs2OFKiE2ZAkk42MXEDyteZsLjyitdEnvjDXjvPTh5Eq6+2l0/MbGvoABeeMGVh8rKcsmrY0e/ozKmTpbATPjk5rryQi+/7K6fnDwJkydbkddYd/hweWHeli1dhY02bfyOyphzsv8sJrx69XJJbMYMV3LoxAk38TUrYUfyxrfdu93vqrAQzjvPahuauOLbPLBQ2DwwfzToNe/f77qkjh1zn+q/9CVo3z68AcahmJoHtnmzKwtWUgIXX+yGyjdqFJYY6hXPObab0CX6PDBLYFXYH1W5Br/mkyfdp/s9e9w/xptugm7dwhdgHIqJBKYK8+fDwoXugdxcV6w5QtcrLYH5xxJYDLAE5o+wvObiYvcp/4MPykcsDhuWtCMUfU9gp0+7YfLbtpVPgxgyJKK/D0tg/rEEFgMsgfkjbK9ZFRYtcp/6VaFnT7fGWAS6q2KdnwmsnQgH/vAHOHgQGjd2LeIo1DW0BOafRE9gNsbZRJ6Im+x8++0uaW3dCk89BZ995ndkyWPDBqaCS17t28OXv2xFeZNcbWszikhrEZkrItu9r6287SIi00Rkh4hsFJGB/r4CS2Ammi65BKZOdaPdDh2Cv/7VVe6wT9uRU1zsFqF85RXSAfr2hQcegNat/Y7M+K+2tRl/AMxT1W7APO97gGuAbt5tKvDn6IdcmSUwE11t2rh/oIMGudFvr7/uqj8UFvodWeLZtw+efBLWrIG0NF4DuOEGV2XDJL061macCDzr7fYsMMm7PxF4Tp3lQEsR8XVVU0tgJvrS0+G66+DGG8u7FP/0J1fCyDSYACxZ4lq4Bw9Cu3bwwAOsg6QdPGPqVmVtxvNUdR+4JAcE5r90AnZXeFo+Pi9GbBOZjX/69IHOnV0R4E8/heefd0O6x41LygEeYXHwIPeCq0sJcPnl7ueZnu5nVMY/aSJSscL2U6r6VMUdvLUZXwa+rarHpfYPOTU94Gv/vyUw469WreCee9y6U/Pnu2tiH30E115ri2SGoqzM/QwXLKALuGoa11+f9PPuDCWqmlvbg7WszbhfRDqo6j6vi/Bzb3s+uLeXpzOwNxJBB8u6EI3/UlLc3LCvfAU6dYLjx10VjxdfdPdN3T791F3revddKClxXYVf+5olL1OnOtZmfA2427t/NzC7wva7vNGIQ4Bjga5Gv9g8sCpsbko5X15zWRmsWOEq2hcXuwEHI0e6ybapqdGNJQLC+v46edIlrfXr3fetWsEXv4h06xZT72GbB+afuuaBicgwYBGwCSjzNv8Idx3sX8AFwKfAzap62Et4/wtcDRQC96qqrwsAWhdikLKzs6mjbzjkY+3atavBx8nJySEvL6/Bx6/rOFGXkuJWee7ZE95+2w3wmDvXjaQbO9Ztj+OBCLW9j0L6WRcXw7JlsHixW8MrNdW1YIcNg/T08JwjjMIZT23v1brOXdPfQrj+duLZOdZmHFPD/gp8LaJBhchaYFVE49NfuM4Rrk+wMf2J9+OP4d//dqPpwA36GDcuORdaLC2FDRtgwYLyrtXu3V05qARa/qSu92O43tvJ0vpL9EoclsCqsAQWg0pLYe1a94878D7o2tV1Lebk+BlZdJSVucS1cCEcOeK2dejgElfXrv7GFgGWwMLHElgMsAQW2nESLoEFnDnjus6WL4fTp9227Gy44go3YjGOuxZrVFTkEvfy5XD0qNvWtq1L3L17J97r9VgCCx9LYDHAElhox0nYBBZw+rT7p14xkbVpA4MHu1JJmZn+xtdQhw+7a35r18KpU25bmzaunmSfPhFb9iRWWAILH0tgMcASWGjHSfgEFnDmjPsnv2JFeQslPR0uvRQGDoQuXeKnlVJc7JY4WbfOXfcL6NIFrrzSXeuKl9fSQJbAwscSWAywBBbacZImgQWUlbnRiqtXwyeflG9v0cJ1tfXuDeefH3sJoKTExfvBBy7+M2fc9rQ0F3Nurhu0kmQsgYWPJbAYYAkstOMkXQKr6NAh14rZuLHyJOjmzd11sm7d3HUzv7oZjx51SWv7dlf7saio/LFOnVwXaN++br2uJGUJLHwsgcUAS2ChHSepE1iAqqtQsXkzfPghnDhR/piIa5FlZ0PHjm5EX5s24b+2VFwM+/fD3r3ulpdXPoow4PzzoUcPd20rgYbCN4QlsPCxBBYDLIGFdhxLYFWousUzt22DnTthzx43NL+i9HS3RlabNu5rs2bQtKm7NWrkKoJUXIZE1bWezpxxA0lOnnRJ8vhxNwjj4EHX2qr6c83MdEP/L7zQtQhbtoz4y483lsDCJ9ETWMJX4rj1yWXVtl3btwN3Ds3hVFEp9/xtZaXHzrv9V7y4ejc353bhcEER//H8mmrPv2NINtf168jeo6d4cOb6ao9/efiFjO11Hh8fOMmPZm2q9nhmdj8APth7jEfmbKn2+Pev7s6g7NasyTvMY29tq/b4T6/rxaUdW5CZ3a/G15fW2q1w8O6W/fxl0c5qj//u1v50bNmYORv28vzyPM67/VeVjvPnOwbROiuDF1fv5qU1+dWe//d7L6dxRir/WLaL1zdWL4U28ytDAXhq4cfM2/p5pccy01N59r7LAZg2bztLdhys9HirJhn8vzsHAfCbtz5kbV7lFkuHFpk8cdsAAH4+5wO27K1cK/HCdln8anJfAH44ayM7D1T84NOIXh0H8/Cdl0B+Pt+e8xH7jp9xSehMCZyEgfnHeSjd/U6+WnQRR7Tyn8iVqcf5Zpp7zXcXdeO0Vm61jUk9xtQ0t9L0rUU9XGLMyHBJsFEjrh10AXde0bXG9x7ATYM6R/S9942rujGsW9sGv/cWbz/IH97bXu3x/5rch4vaNQ36vVfVn+9wv/va3nuS5lYpCPa9V/G9XfG91+KK26r97bRqUv4BJfzvPejVsTkPX3cpAN/+5zr2HTtd6fGB2a146OoeAHz1H2s4Ulh09vWYmiV8AjOmmvR0NwG4/VFo5P0TKSuDkmJokQrZ3VyLamMJFJW5x1Td1/QMyGrqnnMiAyTVdUmmpkJaKnQ8Dy4b6eZrzdpOtUo9sTaQxJg4Zl2IVVgXojH+si7E8En0LsTEnhFpjDEmYVkCM8YYE5csgRljjIlLlsCMMcbEJV8SmIhcLSLbRGSHiPzAjxiMMcbEt6gnMBFJBf4IXAP0Am4XkV7RjsMYY0x886MFdjmwQ1V3qmoR8E9gog9xGGOMiWN+JLBOwO4K3+d724wxxpig+VGJo6ZSBNVmDorIVGBq4HEROVXP86UBJaE8QaJQLaGe56j2Wmo7TqjHj8ZrriDk30kMs9cSAXW9H4N8r559LaH+jUT5b+FcGvo7SehlDfxIYPlAlwrfdwb2Vt1JVZ8CnmroyURktarmNvQ4sSBRXkuivA6w1xKrEuW1JMrriBQ/uhBXAd1EpKuIZAC3Aa/5EIcxxpg4FvUWmKqWiMjXgbeBVOAZVf0g2nEYY4yJb75Uo1fVN4E3o3S6BndDxpBEeS2J8jrAXkusSpTXkiivIyLiohq9McYYU5WVkjLGGBOXkiKBicgvRGSjiKwXkXdEpKPfMdWXiPy3iHzovZ5XRCQu16QXkZtF5AMRKRORuBxllSgl0UTkGRH5XEQ2+x1LQ4hIFxGZLyJbvffWt/yOqb5EJFNEVorIBu+1/NzvmGJRUnQhikhzVT3u3f8m0EtVv+pzWPUiIuOB97zBML8BUNWHfA4rZCLSEygDngS+p6qrfQ4pJF5JtI+AcbipIauA21V1i6+B1YOIjABOAs+pam+/46kvEekAdFDVtSLSDFgDTIrT34kAWap6UkTSgcXAt1R1uc+hxZSkaIEFkpcnixomTscLVX1HVQMTG5fj5tHFHVXdqqrb/I6jARKmJJqqLgQO+x1HQ6nqPlVd690/AWwlTqv8qHPS+zbdu8Xt/61ISYoEBiAij4rIbmAK8FO/4wmT+4B/+x1EkrKSaDFMRHKAAcAKfyOpPxFJFZH1wOfAXFWN29cSKQmTwETkXRHZXMNtIoCq/qeqdgGmA1/3N9q6neu1ePv8J67EzHT/Iq1bMK8jjgVVEs1En4g0BV4Gvl2l9yWuqGqpqvbH9bJcLiJx270bKb7MA4sEVR0b5K4vAG8AD0cwnAY512sRkbuBa4ExGsMXMUP4ncSjoEqimejyrhe9DExX1Vl+xxMOqnpURBYAVwNxPdAm3BKmBVYXEelW4dvrgQ/9iqWhRORq4CHgelUt9DueJGYl0WKMN/DhaWCrqj7udzwNISLtAiOMRaQxMJY4/r8VKckyCvFloDtu1Fse8FVV3eNvVPUjIjuARsAhb9PyeBxRKSI3AH8A2gFHgfWq+gV/owqNiEwAnqC8JNqjPodULyIyAxgFtAX2Aw+r6tO+BlUPIjIMWARswv2tA/zIq/wTV0SkL/As7r2VAvxLVR/xN6rYkxQJzBhjTOJJii5EY4wxiccSmDHGmLhkCcwYY0xcsgRmjDEmLlkCM8YYE5csgRljjIlLlsCMMcbEJUtgxtSDt+7UOO/+L0Vkmt8xGZNsEqYWojFR9jDwiIi0x1U9v97neIxJOlaJw5h6EpH3gabAKG/9KWNMFFkXojH1ICJ9gA7AGUtexvjDEpgxIfKWrp+OW4G5QETiqgixMYnCEpgxIRCRJsAs4LuquhX4BfAzX4MyJknZNTBjjDFxyVpgxhhj4pIlMGOMMXHJEpgxxpi4ZAnMGGNMXLIEZowxJi5ZAjPGGBOXLIEZY4yJS5bAjDHGxKX/D357mBCKmZaYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "nsize = 100\n",
    "xs = stats.norm(0., 1.).rvs(nsize)\n",
    "htexam.normal_likelihood(xs);\n",
    "#plt.yscale('log');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Posterior probability in bayesian statistics\n",
    "\n",
    "With bayesian statistic one can improve a prior \"knowledge\" using the data.\n",
    "\n",
    "Suppose that we have an ensemble of possible hypotheses $\\mu$, with a prior probability $\\pi(\\mu)$. Given some observation, $x$, Bayes' theorem allow us to compute the **posterior probability**, $\\pi'(\\mu)$. The probability of $\\mu$ given the data $x$, $p(\\mu|x)$, is: \n",
    "\n",
    "$$\n",
    "p(\\mu|x) = \\frac{p(x|\\mu) \\, \\pi(\\mu)}{\\int p(x|\\mu) \\pi(\\mu) \\, \\mathrm{d}\\mu}\n",
    "$$\n",
    "\n",
    "\n",
    "As the total probability for the ensemble of hypothesis $\\mu$ should be one, we can compute the posterior, $\\pi'(\\mu)$\n",
    "\n",
    "$$\n",
    "\\pi'(\\mu) = \\frac{p(\\mu | x)}{\\int p(\\mu | x) \\mathrm{d}\\mu}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Exercise:** Show that the posterior of n $x$ normal distributed measurements, $p(x | \\mu, \\sigma)$, with a flat prior on $\\mu$, and known $\\sigma$, is a gaussian with mean $\\mu = \\bar{x}$, the average of the $n$ $x$ measurements, and the sigma $\\sigma/\\sqrt{n}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Example:** Consider the example above. An sample of $n$ measurements from a normal distribution. Obtain the posterior probability of the mean, if the sigma is know. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "posterior integral 1.000\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAEMCAYAAAAWDss+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzt3XuUHWWZ7/Hvr7uTdGgSLkk8YhLTOQQRQSQaGVFRBGQQWYTj5YhXGHVlOEsUb8cRx4Ual8tRZyGjOHiicIQxMri4aEbxICqorCWXEJKQkCABbWhBgQAhkAs0ec4fVTvsNLu796VqX2r/PmvV6tq7alc9lezOk+d933pLEYGZmVkeelodgJmZFZeTjJmZ5cZJxszMcuMkY2ZmuXGSMTOz3DjJmJlZbpxkzMwsN04yZmaWGycZMzPLTV8zT9bT0xNTp05t5inNzDretm3bIiLGLQok9QIrgb9ExMmjtk0BLgVeBWwG3hURf84p3D00NclMnTqVp556qpmnNDPreJK2V7Hb2cAGYHqFbR8CHouIBZJOA74GvCvDEMfk5jIzsw4naQ7wVuD7Y+yyGLgkXb8COE6SmhGbk4yZWec7H/gMsGuM7bOB+wEiYgTYAsxoRmBOMmZm7a9P0sqyZUlpg6STgYci4rZxPl+pamnKFPxN7ZMxM6vVM888w/DwMDt27Gh1KLnr7+9nzpw5TJo0afSmkYhYNMbHXgecIukkoB+YLumHEfG+sn2GgbnAsKQ+YB/g0YzDr0jNfJ7MwMBAuOPfzGrxpz/9iWnTpjFjxgya1I3QEhHB5s2b2bp1K/Pnz99jm6RtETEw0TEkHQN8usLoso8AL4+IM9OO/7dFxP/MMPwxTdhcJuliSQ9JWlf23v6SrpN0d/pzv3zDNLNutWPHjsInGABJzJgxI7OKTdJSSaekLy8CZkjaBHwS+GwmJ6lCNX0yPwBOHPXeZ4FfR8RBwK9pYsBm1n2KnmBKGr3OiLihVMVExLkRsSJd3xER74yIBRFxZETcm0G4VZkwyUTE73h+2135cLhLgFMzjsssH488Au99L9x8c6sjMesK9Xb8/7eIeBAgIh6U9IIMYzLLz1VXwY9+BGvWwB13QJf8D9msVXIfwixpSWnY3cjISN6nMxvf448nP9evh1/8orWxpAYHB5H0vGVwcLDVoVkVrr766laH0NbqTTJ/k3QAQPrzobF2jIhlEbEoIhb19XnEtLXYE088t/71r7cujjJDQ0NExPOWoaGhVodmE7jnnnv41a9+BcD27dt54xvfyLPPPjvm/k8//TRveMMb6Kb/cNebZFYAp6frpwM/zSYcs5yVJ5nf/hZuuaV1sVjHu/baa/noRz8KwMUXX8zb3vY2ent7x9x/8uTJHHfccVx++eXNCrHlqhnCfBnwB+BgScOSPgT8C/BmSXcDb05fm7W/UpKZMyf5+Y1vtC4W6xibNm1i1qxZDA4OcsQRR7D//vtz4IEHMjIywktf+lIAli9fzuLFi3d/5phjjuGuu+4CYPPmzRx22GEAnHrqqSxfvrz5F9EiE7ZfRcS7x9h0XMaxmOWvlGQ+8xn41KeSgQCbNsGCBa2Ny6qT10CNCW5KX7BgAa9//ev55Cc/ydFHH80xxxzDt7/9bV7+8pcDSTPYvffeu0c/2qZNmzjooIMAWLt27e59DzvsMG699dZ8rqMNee4y6y6lJPPSlyZDmXftggsuaG1M1hHWr1+/uxrZuHEjBx988O5tjzzyCPvuu+/u10NDQ8yePZuenuSf2LVr13L44YcD0Nvby+TJk9m6dWsTo28dJxnrLqUkM306vP3tyXrapGEdICKfZQLbt29nx44d7Lffftx///3MmDGDyZMn794+derUPe7UX7169e6kAnDbbbft8Xrnzp309/dn9IfS3pxkrLuUJ5l99knWu+R/lFa/O++8k0MOOQSADRs27F4v2W+//Xj22Wd3J5o1a9bsXr/77rv56U9/uru5bPPmzcyaNavSJJiF5CRj3aU8yey9d7L+5JOti8c6QnlT2dSpU1m1ahUbN27cY58TTjiBG2+8EUgqmV27dvGKV7yCpUuXcsghh3DJJckkKddffz0nnXRScy+ghXzjinWXLVuSn9Onw/b0ibauZGwCH/jAB3avH3300dx77/On/jrrrLM477zzOP7441m7di23334706ZNe95+P/rRj/jqV7+aa7ztxJWMdY+REdi2LRmhNDDgSsYytXDhQt70pjexZcsWenp6KiaYp59+mlNPPXWPQQNF50rGukepYpk2DXp6nksyrmQsIx/84AcB+OMf/1hx++TJk/eoirqBKxnrHuX9MQB77ZVUNdu3wzhTgZhZ/ZxkrHuMTjI9PUmzGYCf2GqWCycZ6x6jkwwkTWfgJrM218zHxLdSEa/TSca6RynJlO6PAXf+d4D+/n42b95cyH+Ay0UEmzdvrvkmTUn9km6RtEbSeklfqrDPGZIelrQ6XT6cWeATcMe/dQ9XMh1pzpw5DA8P8/DDD7c6lNz19/czpzR5a/V2AsdGxJOSJgE3SvpFRNw0ar/LI+KsTAKtgZOMdY9KScaVTNubNGkS8+fPb3UYbSuSEq/0BZ6ULm1T9rm5zLqHKxnrXH2lJwyny5LyjZJ6Ja0meYDkdRFxc4VjvF3SWklXSJrblKhxJWPdxJWMda6RiFg01saIeBY4QtK+wNWSDouIdWW7/BdwWUTslHQmcAlwbL4hJ1zJWPcYL8m4krECiIjHgRuAE0e9vzkidqYvvwe8qlkxOclY9xivucyVjHUoSbPSCgZJU4HjgY2j9jmg7OUpwIZmxefmMusermSsmA4ALpHUS1I4/DgifiZpKbAyIlYAH5N0CjACPAqc0azgnGSse7iSsQKKiLXAwgrvn1u2fg5wTjPjKnFzmXUPd/ybNZ2TjHWP8mfJlHgIs1munGSse7iSMWs6JxnrHr4Z06zpnGSsO+zatedDy0pcyZjlyknGukMpiQwMQG/vc+97CLNZrpxkrDtUaioDD2E2y5mTjHWHsZKMKxmzXDnJWHeo9MAy2LOSKfhDscxawUnGusNYlczkyTBpEoyMwM6dz/+cmTXESca6w1hJBtwvY5YjJxnrDuMlGQ9jNsuNk4x1h2oqGXf+m2WuoSQj6ROS1ktaJ+kySf1ZBWaWKVcyZi1Rd5KRNBv4GLAoIg4DeoHTsgrMLFPVJBlXMmaZa7S5rA+YKqkP2At4oPGQzHLgjn+zlqg7yUTEX4B/Be4DHgS2RMQvswrMLFOuZKygJPVLukXSmrT74ksV9pki6XJJmyTdLGmwWfE10ly2H7AYmA+8CBiQ9L4K+y2RtFLSypGRkfojNWtEpWfJlLiSsc62Ezg2Il4BHAGcKOk1o/b5EPBYRCwAvgl8rVnBNdJcdjzwp4h4OCKeAa4CXjt6p4hYFhGLImJRX5+f9mwt4krGCioSpf8hTUqX0dNXLAYuSdevAI6TpGbE10iSuQ94jaS90mCPAzZkE5ZZxtwnYwUmqVfSauAh4LqIuHnULrOB+wEiYgTYAsxoRmyN9MncTJIRVwF3pMdallFcZtlyJWOdra/U7ZAuS8o3RsSzEXEEMAc4UtJhoz5fqWppymR9DbVfRcQXgC9kFItZflzJWGcbiYhFE+0UEY9LugE4EVhXtmkYmAsMp6OB9wEezSPQ0XzHvxVfxHNJpvypmCW+GdM6mKRZkvZN16eS9JdvHLXbCuD0dP0dwG8imjPtuHvirfi2bUsev9zfn8y6PJqby6yzHQBcIqmXpHD4cUT8TNJSYGVErAAuAv5D0iaSCqZpN847yVjxjfUsmRI3l1kHi4i1wMIK759btr4DeGcz4ypxc5kV33j9MeBKxixHTjJWfBMlGVcyZrlxkrHicyVj1jJOMlZ8rmTMWsZJxopvoiQzMJD8fOqpZBSamWXGScaKb6Ik09OzZ6Ixs8w4yVjxlZJMqe+lEvfLmOXCScaKr1SdVJNk3C9jliknGSu+UpIpNYlV4s5/s1w4yVjxVZNk3FxmlgsnGSs+VzJmLeMkY8XnSsasZZxkrPhcyZi1jJOMFV+XVzKDg4NIet4yODjY6tCsC3iqfyu+Lq9khoaGqPR8KqnSE3nNsuVKxoqvyysZs1ZykrHiqyXJFLCSMWslJxkrvlqay1zJWIeRNFfS9ZI2SFov6ewK+xwjaYuk1elybqVj5cF9MlZsEbBtW7LuSsaKaQT4VESskjQNuE3SdRFx56j9fh8RJzc7OFcyVmzbtyeJZsoU6O0de78Cd/xbsUXEgxGxKl3fCmwAZrc2quc4yVixVdNUBu74t3bXJ2ll2bKk0k6SBoGFwM0VNh8laY2kX0g6NMdY9+DmMiu2apOMnydj7W0kIhaNt4OkvYErgY9HxBOjNq8C5kXEk5JOAn4CHJRPqHtyJWPFVmsl4yRjHUjSJJIEszwirhq9PSKeiIgn0/VrgEmSZjYjNicZK7ZaKxn3yViHUXJX7UXAhog4b4x9Xpjuh6QjSf7t39yM+NxcZsXmSsaK73XA+4E7JK1O3/sc8GKAiPgu8A7gf0kaAbYDp0WlaSBy4CRjxVZtktlrr+Tntm2waxf0uMi3zhARNwLjzhEUERcAFzQnoj35N8mKrdok09OzZ6Ixs0w4yVixVZtkyvdxv4xZZpxkrNhqSTLulzHLnJOMFZsrGbOWaijJSNpX0hWSNqaTsx2VVWBmmSglDFcyZi3R6OiyfwP+X0S8Q9JkYK8MYjLLjisZs5aqO8lImg68ATgDICKeBp7OJiyzjLhPxqylGmku++/Aw8D/lXS7pO9LquI32ayJSgmjlEDG40rGLHONJJk+4JXAhRGxEHgK+OzonSQtKc0cOjIy0sDpzOrgSsaspRpJMsPAcESUppS+giTp7CEilkXEoohY1NfnCQasydwnY9ZSdSeZiPgrcL+kg9O3jgNGP4nNrLVcyZi1VKOlxUeB5enIsnuBf2g8JLMMuZIxa6mGkkxErAbGfZCOWUu5kjFrKd/xb8XmSsasYZIGJPXW81knGSs2VzJmNZPUI+k9kn4u6SFgI/CgpPWSviGp6kc3O8lYcUW4kjGrz/XAgcA5wAsjYm5EvAA4GrgJ+BdJ76vmQB5TbMW1c2fyALJJk5JlIq5krANJmgtcCrwQ2AUsi4h/G7WPSKYBOwnYBpwREavGOezxEfHM6Dcj4lHgSuBKSVX8UjnJWJHVUsWU7+dKxjrLCPCpiFglaRpwm6TrIqL8lpK3AAely98BF6Y/KyolGEl3AGvLljuA0yPiK5WSUCVuLrPiqjXJuJKxDhQRD5aqkojYCmwAZo/abTFwaSRuAvaVdEAVh38j8D1gO3AasA54ay3xuZKx4nIlY11G0iCwELh51KbZwP1lr4fT9x4c73hp89gN6ULa4f/5WmJyJWPF5UrGiqOvNAdkuiwZvYOkvUn6Sz4eEU+M3lzhmDHRSUePIouIu4HDa4jblYwVWL2VzFNPJSPTVOn30qwlRiJizBvf0074K4HlEXFVhV2Ggbllr+cAD1Rx3mWSDgT+QtIn0w+sk7RXRGyrJnBXMlZctSaZvj6YMiUZkbZjR35xmWUoHTl2EbAhIs4bY7cVwAeUeA2wJSLGbSoDiIg3RcSLgXcBPwc2kTycco2kjdXE50rGiqvWJFPad+fO5LNTp+YTl1m2Xge8H7hD0ur0vc8BLwaIiO8C15AMX95EMoR53HkmJSkidjenRcR9wH0kyaq0z/RqgnOSseKqJ8nsvTc8+mjS+T9zZj5xmWUoIm6kcp9L+T4BfKSGw14v6Urgp2mCASCdDPn1wOkkN2z+YKIDOclYcdVbyZR/1qw7nQh8ELhM0nzgcWAqSRfLL4FvphMkT8hJxoqr3koGPIzZulpE7AD+Hfj3dFDBTGB7RDxe67Hc8W/F5UrGrGaSTpX0bUkvgeTu//SGz5oTDDjJWEYGBweRVPUyODiY+3n/delSAD6zdGn153UlY10uIn4CfBM4TtJZjR7PScYyMTQ0RERUvQwNDeV+3k+feSYAX7/ggurP60rGjIi4NyIuBM6W9DJJu7tWShVOtdwnY8XlPhmzRu0PfAd4qaRHSOZFeyWwoNoDOMlYcblPxqxR90XEmwAkzQEOIZkss2pOMlZcrmTMGrWPpNcC6yNimGR6mpo4yVhxuZIxa9TewKeBQyVNIWkuWxcR/7vaAzjJWHG5kjFr1FERcQ+ApH7gZcChtRzAScaKy5WMWUNKCSZd3wGsSpeqeQizFVepGnElY9YyTjJWXK5kzFrOScaKy30yZi3nJGPF9PTTMDICPT3Jg8iq5UrGLFNOMlZM5VVMLY9RLiUZVzJmmXCSsWIqJZlS81e1Svu7kjHLhJOMFVM9/THl+7uSsQ4i6WJJD0laN8b2YyRtkbQ6Xc5tVmy+T8aKqd4k40rGOtMPgAuAS8fZ5/cRcXJzwnmOKxkrpnqTzOTJ0NsLzzyTDB4w6wAR8Tvg0VbHUUnDSUZSr6TbJf0si4DMMlFvkpFczVg76pO0smxZUscxjpK0RtIvJNU0NUwjsmguO5tk0rTpGRzLLBv1JpnSZ7ZsSfpl9tsv27jM6jMSEYsa+PwqYF5EPCnpJOAnwEHZhDa+hiqZ9PkCbwW+n004ZhlpJMm4krGCiYgnIuLJdP0aYJKkmc04d6PNZecDnwF2ZRCLWXYarWTKj2HW4SS9UEpuGJN0JMm//Zubce66k4ykk4GHIuK2CfZbUmpHHBkZqfd0loPBwUEkPW8ZHBysaX9JzJs3ryUxjSmLSsbDmK1DSLoM+ANwsKRhSR+SdKakM9Nd3gGsk7QG+BZwWkREM2JrpE/mdcApaftePzBd0g8j4n3lO0XEMmAZwMDAQFMuyqozNDREpe+ZxrhDfqz9WxnTmFzJWBeJiHdPsP0CkiHOTVd3JRMR50TEnIgYBE4DfjM6wZi1jCsZs7bg+2SsmB57LPm57761f9aVjFlmMrnjPyJuAG7I4lhmmSglmXqGILuSMcuMKxkrpkaSjCsZs8w4yVgxuZIxawtOMlZMrmTM2oKTjBVTKcnsv3/tn3UlY5YZJxkrnmefTeYeA9hnn9o/70rGLDNOMlY85Qmmt7f2z7uSMcuMk4wVTyP9MeBKxixDTjJWPI0mGVcyZplxkrHieTR9QKCTjFnLOclY8TRayZQ+92hbPs3WrKM4yVjxNJpkZsxIfm7eDM2ZDd2ssJxkrHgaTTL9/Unn/zPPwNat2cVl1oWcZKx4Gk0ysGc1Y2Z1c5Kx4skiycxMH3/+yCONx2OWM0kXS3pI0roxtkvStyRtkrRW0iubFZuTjBVPlknGlYx1hh8AJ46z/S3AQemyBLiwCTEBTjJWRFk2l7mSsQ4QEb8DxhsOuRi4NBI3AftKOqAZsTnJdIHBwUEkPW+ZN29eq0PLxzhJZt68eRX/LEYv37rsMgA+/v73V9w+ODhY8dRj/VmPt7Tq76HaP4t6r3ms/a0lZgP3l70eTt/LXSZPxrT2NjQ0RHTTUNxxksyf//zn6o7xpS/BF7/I+Z//POd/+cvP2yyp4sc66c+66j+LVK3XPNb+Vpc+SSvLXi+LiGU1fL7SX0ZTvqhOMlY8jUzzX+KOf2svIxGxqIHPDwNzy17PAR5oLKTquLnMiqXRaf5LPITZimUF8IF0lNlrgC0R8WAzTuxKxoql0Wn+S1zJWAeRdBlwDDBT0jDwBWASQER8F7gGOAnYBGwD/qFZsTnJWLFkMbIMPITZOkpEvHuC7QF8pEnh7MHNZVYsWSUZD2E2y4STjBVL1knGk2SaNcRJxoolqySz114wdSrs3OknZJo1wEnGiiWrJAPu/DfLgJOMFUuWScbDmM0a5iRjxeJKxqytOMlYseSRZFzJmNXNScaK5dF0Itosm8tcyZjVzUnGisWVjFlbcZKxYsmj49+VjFndnGSsWFzJmLWVupOMpLmSrpe0QdJ6SWdnGZhZXbKY5r/ElYxZwxqZIHME+FRErJI0DbhN0nURcWdGsZnVJqtp/ks8hNmsYXVXMhHxYESsSte3Ahto0uM8zSrKapr/Et+MadawTPpkJA0CC4GbK2xbImmlpJUjIyNZnK4qfu54/cZ69nuWz6Kv9RxV7Z9lfwy0bSUz1nd7rCXLvzezWqnR55FL2hv4LfCViLhqvH0HBgbiqSZNNihpzOeOd8oz2LPSNde8ciW8+tWwcCGsWtX48SKSiTJ37Egmydxrr92bmvH9arfvcK3xdM33rgkkbYuIgVbHUY+GKhlJk4ArgeUTJRiz3GVdyUju/DdrUCOjywRcBGyIiPOyC8msTlknGfAwZusIkk6UdJekTZI+W2H7GZIelrQ6XT7crNgaGV32OuD9wB2SVqfvfS4irmk8LLM65JFkXMlYm5PUC3wHeDMwDNwqaUWFkb6XR8RZzY6v7iQTETcCyjAWs8bkWck4yVj7OhLYFBH3Akj6T2Ax0Ba3k/iOfyuOPCsZN5dZa/WVRummy5KybbOB+8teD1P5dpK3S1or6QpJc3ONtkwjzWVm7aWUCFzJWPGMRMSiMbZValEaPazvv4DLImKnpDOBS4BjswxwLK5krDj++Mfk54EHZndMd/xb+xsGyiuTOcAD5TtExOaI2Jm+/B7wqibF5iRjBREB69cn6y97WXbHdce/tb9bgYMkzZc0GTgNWFG+g6QDyl6eQjJDS1O4ucyK4aGHkgeW7bMPvOhF2R3XlYy1uYgYkXQWcC3QC1wcEeslLQVWRsQK4GOSTiGZc/JR4IxmxeckY8VQqmIOPTS5iTIrrmSsA6S3jlwz6r1zy9bPAc5pdlzg5jIrijyayuC5Subhh7M9rlmXcJKxYrgzvSXg0EOzPe7s2TBlCgwPw+OPZ3tssy7gJGPFUN5clqVJk+CII5L1227L9thmXcBJxjpfXiPLSl796uTnrbdmf2yzgnOSsc6X18iyEicZs7o5yVjny2tkWcmRRyY/b7kl+2ObFZyTjHW+PJvKAF7yEpg+Pen8/+tf8zmHWUE5yVjny2tkWUlPD7wqnYXDTWZmNXGSsc6X18iycu6XMatL1yWZefPmISmTZXBwMJOYBgcHMzn+WMeZN29eJnG2pfKRZXkmmVH9MmN9j7L8s27GOVoVz1jf1Vp/17L63bH8KGL0jND5GRgYiKeeeqop55JE3teW1TnGOk6tx2/GNbedv/0NXvjCZGTZY4/l0/EPcN99MG8e7L9/MsVMXucpkPG+j1l9t7P63Wl3krZFxECr46hH11UyVjB5jywrmTsXXvCCZKj0n/6U33nMCsZJxjpbM5rKIElg7pcxq5mTjHW2O+5IfuY1fLmc75cxq5mTjHWue+6BH/4wWX/ta/M/nysZs5o5yVhnioB//EfYvh3e857nqow8lZLMLbfA7bfnfz6zAnCSsc506aXw618nDxU7//zmnHPmTDjjDNi5E04+OZkBwKwNSDpR0l2SNkn6bIXtUyRdnm6/WdJgs2JzkrHO87e/wSc+kax/85swa1bzzv3d78LRR8MDDySJZuvW5p3brAJJvcB3gLcALwPeLWl0J+WHgMciYgHwTeBrzYqvMx6/XMdTCWfW+blWnGOs49R6/GZccy5G388Q8dx7IyOwY0ey3HUX/PKXcM01yT0xb34zvO99zY11yhS4+mo46ihYswaOPRYWL4aFC2HBApg8OXkGTV9fMiKtNKy6C++rGe/7mNV3O6vfnabYe2+YOjWPIx8JbIqIewEk/SewGLizbJ/FwBfT9SuACyQpmnAzUWfcjDkwANu2ZR+Qda7585PmsvnzW3P+TZuSRPPII605v3WeCy+EM8+s66Pj3Ywp6R3AiRHx4fT1+4G/i4izyvZZl+4znL6+J90n9y9wZ1QyM2fWnGQefuQRZpWez56TrM4x1nFqPX4zrjk3o/+nX6oCenuT//319yc3Qx5/PJxwQlI59Pa2JlZIqpZ16+Daa2H16mQgwPAwPPNMsoyMPFeNFejO81qM933M6rud1e9OU/T3N/LpPkkry14vi4hl6XqlMnn0l66afXLRGZVMHTytjFlreVqZ7ExQyRwFfDEi/j59fQ5ARHy1bJ9r033+IKkP+CswqxnNZe74NzPrbLcCB0maL2kycBqwYtQ+K4DT0/V3AL9pRoKBTmkuMzOziiJiRNJZwLVAL3BxRKyXtBRYGRErgIuA/5C0CXiUJBE1hZvL2uAcbi6zInJzWXY8C7OZmVkFDSWZie4yNTOz7lZ3kqnyLlMzM+tijVQyu+8yjYingdJdpmZmZkBjSWY2cH/Z6+H0PTMzM6CxIcxV3UEqaQmwpLRd0vY6z9cHjNTyATVhvqg6z/G8axnrOLUevxnXXKbmv5M25mvJwXjfxyq/q7uvpdbfkSb/Lkyk0b+TXCY9a4ZGkswwMLfs9RzggdE7pVMfLBv9fq0krYyIRY0epx0U5VqKch3ga2lXRbmWolxHPRppLqvmLlMzM+tidVcyY91lmllkZmbW8RqaViYirgGuySiWiTTc5NZGinItRbkO8LW0q6JcS1Guo2ZNnVbGzMy6i6eVMTOz3HRUkpH0ZUlrJa2W9EtJL2p1TPWS9A1JG9PruVrSvq2OqR6S3ilpvaRdkjpy9ExRpkeSdLGkh9KnIHYsSXMlXS9pQ/rdOrvVMdVLUr+kWyStSa/lS62Oqdk6qrlM0vSIeCJd/xjwsoio73mmLSbpBJJnOoxI+hpARPxTi8OqmaRDgF3A/wE+HRErJ/hIW0mnR/oj8GaSYfm3Au+OiDvH/WAbkvQG4Eng0og4rNXx1EvSAcABEbFK0jTgNuDUDv07ETAQEU9KmgTcCJwdETe1OLSm6ahKppRgUgM06fGheYiIX0ZE6easm0juM+o4EbEhIu5qdRwNKMz0SBHxO5JnhXS0iHgwIlal61uBDXTobCKReDJ9OSldOvbfrXp0VJIBkPQVSfcD7wXObXU8Gfkg8ItWB9GlPD1SG5M0CCwEbm5tJPWT1CtpNfAQcF1EdOy11KPtkoykX0laV2FZDBAR/xwRc4HlwFmtjXZ8E11Lus8/k0w3sbx1kY6vmuvoYFVNj2TNJ2lv4Erg46NaMTpKRDwbEUeQtFYcKaljmzLr0XaPX46I46vc9UfAz4Ev5Bh+fTKJAAAB2klEQVROQya6FkmnAycDxzXredv1qOHvpBNVNT2SNVfaf3ElsDwirmp1PFmIiMcl3QCcCHT04IxatF0lMx5JB5W9PAXY2KpYGiXpROCfgFMiYlur4+linh6pzaSd5RcBGyLivFbH0whJs0ojRyVNBY6ng//dqkenjS67EjiYZDTTEHBmRPyltVHVR9ImYAqwOX3rpk4cKSfpfwDfBmYBjwOrI+LvWxtVbSSdBJzPc9MjfaXFIdVF0mXAMcBM4G/AFyLiopYGVQdJrwd+D9xB8rsO8Ll0hpGOIulw4BKS71YP8OOIWNraqJqro5KMmZl1lo5qLjMzs87iJGNmZrlxkjEzs9w4yZiZWW6cZMzMLDdOMmZmlhsnGTMzy42TjHU1STdIOjhdn9Hpz2IxazdOMtbtFgB3p+uHk9xlbmYZcZKxriVpHvCXiChNXXI4sLaFIZkVjpOMdbMj2DOpvAonGbNMOclYN3sF0A+7Z/hejJvLzDLlJGPd7AigR9IakqesbgBOb21IZsXiWZita6WPW1iYPkfezHLgSsa6kqRpwC4nGLN8uZIxM7PcuJIxM7PcOMmYmVlunGTMzCw3TjJmZpYbJxkzM8uNk4yZmeXGScbMzHLjJGNmZrn5/+lHyQ4LliQiAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "htexam.normal_posterior(xs);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mu mean : -0.20607850253635096 , mu std : 0.31622776601683794\n",
      "posterior integral 1.000\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAEMCAYAAACLA8K2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzt3XmYFNW5+PHvOzMMMwz7oqzOIBBWkUVBXFgEjeLCogaNJuISrrnXaJKbmz0x6vXnTXJjjInxaiKKioCKiFsUUVDQoAgCIktA1mERGLZhmH3e3x+nmxlmYbp7uruqe97P89TT3VXVVad6eurtc+qtc0RVMcYYY7yW4nUBjDHGGLCAZIwxxicsIBljjPEFC0jGGGN8wQKSMcYYX7CAZIwxxhcsIBljjPEFC0jGGGN8wQKSMcYYX0jzugBVpaSkaGZmptfFMMaYhHH8+HFV1aSoXMQsIIlIN+AZoCNQATyhqn861XsyMzMpKCiIVZGMMSbpiEih12WIlljWkMqA/1TVlSLSAlghIu+o6roY7tMYY0yCilk1T1X3qOrKwPN8YD3QJVb7M8aYhJSfD7m5XpfCF+LS7igiOcBg4ON47M8YYxJCURE89xzMmAHbtnldGs9JrIefEJHmwPvAA6r6ci3LpwHTANLT04cWFxfHtDym4XJycti+fXuN+dnZ2Wyzf6qI1PWZ1iWen3VpaSm5ubkUFRXFZX+NhioUFEBZGaSmQlYWpNRdR8jIyKBr1640adLkpPkiclxVs2Jd3HiIaUASkSbA68DbqvpQfetnZWWpJTX4n4hQ2/emrvmmfuF+dvH8rLdu3UqLFi1o164dIhKXfSY9VTh4EIqLXRBq3x7S6r6kr6rk5eWRn59P9+7dT1qWTAEpZk124r65TwLrQwlGxhh/KioqsmAUTapw5IgLRiLQrt0pgxG4HyDt2rVL+lpqLK8hXQB8C7hYRFYFpvEx3J8xJkYsGEWJKhw9CsePVwajak1wdWkMf4OYpX2r6lIg+T9BY4wJ1bFj7roRQJs2kJ4el92KyHTgSmCfqg4IzGsLzAFygG3AN1T1UKB160/AeOA4MDWYMR1rSXF3rzHG+F5BgUvxBheMMjLiufengcuqzfsp8K6q9gLeDbwGuBzoFZimAY/FqYwWkIwxJuaOH3fXjQBatYI4d5Gmqh8AB6vNngDMCDyfAUysMv8ZdZYBrUWkUzzKaQHJGJPwdu7cyZgxY+jbty/9+/fnT3+q2UtZ8+bNY7LvwsJCRo0aRXl5eV0rwOHDjL72WrYdOuTSu+tQUlLCyJEjKSsri0lZqzldVfeA68gAOC0wvwuws8p6ucSpUwMLSMaYhJeWlsYf/vAH1q9fz7Jly3j00UdZty76vZTNmTOHadOm8bvf/e7EvOnTpzN58mRSU1NrvqGoCA4dcs+D9xqdQnp6OmPHjmXOnDnhFCtNRD6tMk0L5821qO3af1zuMfBVb9/GGJ/7zW98ud1OnTrRqZNrVWrRogV9+/Zl165d9OvXr9b1H3roIaZPnw7A7bffzve//30A7r//fmbOnEm3bt1o3749Q4cO5Uc/+hEAjzzyCIcPH+aWW27hBz/4AT/+8Y8BmDlzJs8///yJbY8ZM4af//znXDJyJL/8yU84euwYjzz0kAtIAaNHj+bxxx+nd+/e5OXlMWrUKNauXQvAxIkT+dnPfsaNN94Y6uGXqeo5IX9Ylb4SkU6quifQJLcvMD8X6FZlva7A7gi2HzYLSMaYpLJt2zY+++wzhg8fXuvyFStW8NRTT/Hxxx+jqgwfPvxEk9vcuXP57LPPKCsrY8iQIQwdOhSA48ePc99997Fp0ybKysp45plnANfEtmXLFnJyck5s/9577+XXv/wl+6ZM4bO1a3n1hRegZcuTyrB582Z69eoFwJo1azjrrLNOLBswYADLly+P5kdSl1eBm4H/CTzOrzL/ThGZDQwHjgSb9mLNApIxJnSxqiGFYNy4cezdu7fG/AceeIAJEyYAcOzYMa655hoefvhhWlYLAkFLly5l0qRJZAWazyZPnsySJUuoqKhgwoQJBMdku+qqq068Z9GiRWRnZ9OmTRsAOnToAMCBAwdo3br1Sdsfed55aGkpDz3+OItff53Utm3dPUcB27dvp0uXLqQEuglas2YNAwcOPLE8NTWV9PR08vPzadGiRXgfUh1EZBYwGmgvIrnAPbhA9IKI3AbsAK4LrP4mLuV7My7t+5aoFCIEFpCMMQlh4cKFp1xeWlrKNddcw4033sjkyZPrXK+uLpdO1RXTBx98wKWXXlpjfmZm5sm9J5SU8PnSpez56ivat29Pi27dTgpGAKtWrTopAK1YsYIpU6actE5xcTEZUUwLV9Ub6lg0tpZ1FfiPqO08DJbUYIxJeKrKbbfdRt++ffnhD394ynVHjhzJK6+8wvHjxykoKGDevHlcdNFFXHjhhbz22msUFRVx7Ngx3njjjRPv+eijjxg+fDj79+8H4OjRo6xbt442bdpQXl7uglJJCXvWrePGO+9k/syZZLVqxdsLFtTY/+rVq08EsU2bNjF//vyTmuzy8vLo0KFDjU5UGwMLSMaYhPfhhx/y7LPP8t577zFo0CAGDRrEm2++Weu6Q4YMYerUqQwbNozhw4dz++23M3jwYM4991yuvvpqzj77bCZPnsw555xDq1atWLx4MU2aNOGBBx6gR48eXH755Tz33HP06dMHgEsvvZSlixdzfOdOJt92G3+4/376nncev/rVr/hNLU2cq1atoqKigrPPPpv77ruPvn37MmPGjBPLFy1axPjxjbSXNVX1zdSsWTM1/ue+NqHPN/UL97OL52e9bt26uO3La/n5+aqqWlBQoEOHDtUVK1bU+56VH3+sN11zjequXaoHDqhWVNS63qhRo3Tr1q3ao0cPPXr0aJ3bmzRpkm7YsKHWZbX9LYAC9cH52xWFLCA10vfbNSRjjAmYNm0a69ato6ioiJtvvpkhQ4ac+g2lpQzu1o0x559PeVpajQSG6vLz80lJSakzWaGkpISJEyfSu3fvhhxG3IhICnA9cCNwLlAMNBWR/bjkiCdUdVPI21MfjV9j4yElBhsPKfr8PB7S+vXr6du3b1z2lVBKSiAvz/Xg3bQp1BOMnn76aSZOnFgjKy8ctf0tvBwPSUTeBxbiUsbXqmpFYH5bYAzwTWCeqj4XyvashmSMMeEKMxgBTJ06NT5li69xqlpafaaqHgTmAnMDA7WGxAKSMcaEo2owyshwPXc3grGKalNbMIpknaDEz7KrqID58+HLL70uiTEm2RUXWzCqhYh8LiIzReQnInK5iHQVkV+Eu53ED0hr1sBnn8Hzz8OGDV6XxhiTrIqK4OBBF4wyMy0YnWwU8DegEJfksBa4ItyNJH5AOvtsGD4cysvhhRfg88+9LpExJtkUFlYGo2bNoHVrC0ZVqOpBVV2sqo+o6s24jLuQs+uCEj8gicBll8FFF7nmu5dfhpVxGW3XGNMYHD9eOYREVpYbYM+C0UlEpFfV14FU74F1rF6n5EhqEIGxY9349O++C6++6qrX55/vdcmMMYmsoKBypNfmzaFFCwtGtXtCRHoAu4A1QAawVkSaqerxUDeSHAEp6KKLXArmm2/CggWumn3xxfYFMsaERxWOHYP8fPe6ZUsXkEytVHUMgIicAQwCzg48rhaRclXtE8p2Er/Jrrphw2DyZEhJgSVLXHCqqPC6VMaYKJg3bx4iwoYQE5jqHV484KShw1Xh6NHKYNSqVUjB6NZbb+W0005jwIABIZUtGYjIGVWnwOxVwAzgB8A44FoRqX0skGqSLyABDBwIU6ZAWhosXw5z50J8xqg3xsTQrFmzOOecc5g9e3aNZWEPL17FiaHDZ8+Gw4ddUx24TLp6hh0Pmjp1Km+99VboB5McZgBPBx6rTk8HpqeAR4CJoWwsOQMSQO/ecNNNrgnviy9g5kx3D4ExJiEdO3aM999/nyeffJJZs2adtOyRRx5h48aN3HLLLbz88ssn5s+cOfPE4H3ghg7fuHEj4IZ5qFqbmXj11cycMcM19Yu43hcCg/WFYuTIkbRt2zbSw0tIqjpGVS8OPFadLq42PRPK9pLrGlJ1OTlwyy3w3HOwdSs8/TTceKO1BRvTAFMe/2eNeVcO7MS3RuRQWFLO1Kc+qbH82qFdue6cbhwsKOG7z604admcfxsR0n5feeUVxo0bx8CBA8nKymLlypUMGTIkrOHF6xw6vLycAZ06sXzlStfc37YtpKdz0UUXkR9suqvif//3fxk3blxI5W4MRORhVf1+Q7cTdkAK9O7aXFWPNnTncdGxI9x2Gzz7LOzZA08+6WpO7dp5XTJjTBhmzZrFtGnTAPjGN77BrFmzGDJkSMjDi9c5dHhZGeTlkVpR4YYOT0+nRXo6AEuWLInnISayYyLyGnC9qhaIyKXAPap6QTgbCSkgicjzwB1AObACaCUiD6nq78MttSfatHFBaeZM2L3bBaVvfhO6dvW6ZMYknFPVaDLTU0+5vG1Wesg1oqry8vL45JNPTjTHTZkyhVGjRvG73/0u5OHFax06fPJkOHDAJT6lpVFcWkpGlaEhrIYUGlX9pYh8E1gsIsVAAfDTcLcT6jWkfoEa0UTcGBdnAN8Kd2eeysqCqVOhVy93o9uMGdbVkDEJ4qWXXmL8+PE0bdoUgO7du9OxY0eWLl0a2vDi1DF0eJcuLhilp5MnUmPo8CVLlrBq1aoakwWjk4nIWOA7uEDUAbhLVcOuXoYakJoEuhCfCMwPp/dWX0lPhxtugCFDoLQU5syBZcu8LpUxph6zZs3itddeIycn58S0fv16pk2bVv/w4kuXAtWGDv/1r+nbowczXnjBJS60a8ei999v0NDhN9xwAyNGjGDjxo107dqVJ598MirHniB+AfxKVUcD1wJzROTicDcS6jWkx4FtwGrgAxHJBo6EuzNfSEmBq65yzXjvvgtvveW6Bfn6190yY4zvLF68OKL33XnnnTz00EOMGzeONWvW8NnKlbQAd9MrnNT7wvPPP8+DDz4YcRmrZ/41Jqp6cZXnn4vI5bjxkMLqLifUM/BrqtpFVccHxnDfAdwazo58RcT16nDNNZCaCh9/DLNnW1q4MUlm8ODBjBkzhiNHjrihw8vKKoNR69auBwaRhBs63C9Eau8GR1X3AGNPtU6t2wtlGGQRWamqQ6rNW6GqQ0PdUSg8GcJ8+3YXjAoL4fTTXbJDq1bxLUOCsSHMo8+GMI+x8nLXW3dpaeU9RoHrUYnEh0OYL8bVhOar6o4q89OBC4GbgUWq+nQo2ztlk52I9AH647LqJldZ1BLXeV7iy86G73zHZeB99RX87W9w/fWWgWdMsigtdYPqVVS4FpG2baFJyKNqm1O7DNdaNktEugOHcbEhFVgA/FFVV4W6sfquIfUGrgRaA1dVmZ+Py6hIDm3bwu23u/GUgjfQXn2164LIGJO4CgtdV0CqLqmpTRsXlExUqGoR8Ffgr4HEt/ZAoaoejmR7pwxIqjofmC8iI1S15u3ZySQz090w+49/wKefunGV9u1zw1pYb+GmkVNVwrgU4L3qvXVnZib8oHp+bw4PZF/vAde6pqph31cTalLDZhH5uYg8ISLTg1O4O/O91FS44goYP95l3C1dCrNmubGVjGmkMjIyyMvL8/0J8YSKCpc5W3XoiCQIRnl5eWRk+O9KiYj8l4h8JCI9q8zeJSJ3hL2tEJMaPgKW4HppONGPu6rODXeHp+JJUkNdtmyBF190Vf527dz9S+3be10qX7Ckhujzc1JDaWkpubm5J/V64FsVFa6n7vJyF4CaNUua60UZGRl07dr1pBt3wdukhsD+HwemA1+q6oEq899Q1SvC2laIAWmVqg4Ku6Rh8lVAAvcra/Zsl+zQtKkbZ8nSQi0gxYCfA1LC2LwZXnrJtWi0b++SkxrBj0gfBKRvAF8D/qyqRwLz2gObVbX1Kd9cTahNdq+LSFi3MAea9faJyNpw3ucrwT7w+vd39yjNmgXvvWcD/hnjJ6rwwQcuU7aoCPr0cZmzjSAY+YGqvoDrKGGziCwXkQdwN8RuDHdbodaQ8oEsoCQwiSuH1jkKoIiMBI4Bz6hqSEMo+q6GFKQKH30ECxe65z17utpSs2Zel8wTVkOKPqshRaiwEObNg3/9yzXRjR4NI0cm9PWicHldQ6pSjkxgDDASF5BKVDWsTv9CCkiREpEc4PWED0hBW7e660rHj7ubZ6+7rlHer2QBKfosIEVgzx53q8ahQy6L7ppr3I/FRsYvAak6ERmlqu+H9Z4Qa0gC3Ah0V9X7RaQb0ElVa47EdfL7ckimgARw5IgLSrm5Livv0kth2LDG9ovMAlKUWUAKgyqsWOH6oSwrg06d4BvfcE3sjZAfApKI/LCW2UeAFeHcGBvqNaS/AiOAbwZeHwMeDXUnpyIi00TkUxH5tKysLBqbjK1WrdwotOed5zJ5/vEPF6ASIQMpTDk5OYhIjSk7O7vW9bOzs2tdv+qInfFW1zF4VaZwP1NTTUmJa6J7/XUXjM45x13nbaTByEfOwY2Z1yUwTQNGA38TkR+HupGw+rITkc9UdXBg3mpVPbue9+WQbDWkqr74Al591SU8tGkD114LXbp4XaqoidavcC9/zfutNpcMn6ln9u51P/7y8lyvC1deab2pEFoNSUS24XrYKQfKVPUcEWkLzAFycKM5fENVD0VYhreBa1T1WOB1c+AlYBKultQvlO2EWkMqFZFUQAM76wBYqln//vBv/+aaDA4dgunTXfJDYztRGBNLqrB8Ofz97y4YnX66y6KzYBSuMao6SFXPCbz+KfCuqvYC3iWCEV6rOAOX8BZUCmSraiEQ8jAKoQakR4B5wGmBlL6lwP871RtEZBbwT6C3iOSKyG2hFiqhtG3rmgyGD3dNeAsWuPTTYBf3xpjIHT/u7gV8443KJrrbb4cOHbwuWTKYAMwIPJ+BG4A1Us8Dy0TkHhH5DfAhrsPVLGBdqBsJOcsu0PP3WFzK97uquj7sItcj4Zrsqtu4EebPd/9EWVkwcaIbMj1BJUPzkjXZJbAtW9z1ovx8yMhwA2v27+91qXxHREqAz6vMekJVn6i2zlbgEK6V63FVfUJEDle9cVVEDqlqxBfjRGQobsgJAZaq6qdhbyPEa0h/Auao6kdhlzIMCR+QAI4edf9EW7e61+ee6zLxErD7kmQ4eVpASkBlZW40538G+nM+4wx331/rsG76bzRCvIbUWVV3i8hpwDvA94BXoxyQzsbdg6TAElVdHe42Qm2yWwn8UkQ2i8jvReScet/RWLVsCd/6FlxyiUsLX74cHn8cdu/2umTG+N/evfDEEy4YpaTAmDEwdaoFowZS1d2Bx324yy/DgK9EpBNA4HFfpNsXkbuBmbjhJ04DnhOR74W9nTDvfWgLXANcD5wRuBgWNUlRQ6pqzx43jMX+/e6fa+RIN3R6gozHkgy/5q2GlCAqKuDDD2HxYncttl07VytKoqzVWKmvhhS4jpOiqvmB5+8A9+EuweSp6v+IyE+Btqoacop2tX2sAUaoakGVff5TVcPKPKlvgL7qegJ9cGmCIV+oarQ6dYJp01zzw7Jl7p9t40aYNAlOO83r0hnjDwcOwCuvuJvNwTVzX3KJS+020XA6MM/1b0Aa8LyqviUiy4EXAglnO4DrGrAPocpIEIHnYfcWEOo1pN/i8sm34PLW50U6IuCpJF0NqaqtW13Cw+HDroY0ahRccIGva0vJ8Gveakg+VlHhfqi99567btSyJUyYAD16eF2yhOKjnhpuxjUHgsvYe1pVHw5rOyEGpH/H9c6Qo6r3icgZQMf6ug4KV1IHJHA30C5Y4Lo9AVeDmjABOnb0tlx1SIaTpwUkn9q3z/1A27XLvR40CC67zGXTmbD4ISAFyjEUuABXM/pAVT8LexshBqTHcDfCXqyqfUWkDbBAVc8Nd4enkvQBKWjLFtfDw+HD7trSBRe4GlNauC2osZUMJ08LSD5TVuauFX3wgbtW1LKlS+dO4NsjvOaXgBQNMe06KFyNJiCB65Nr4UKXhafqLuJedRV42O9bdclw8rSA5CM7dsBrr7kkH4ChQ921IqsVNYiXAUnc0ES1fRHrHaKoNqH+JLeug6ItPR3Gj4ezznK1pf374emnXdPFJZe4G2uNSQaFhS6x59PAfZI+/PFlIqOqLaK5vVBrSDcCU4AhuC4mrgV+qaovRrMwjaqGVFWwGWPJEvc8MxPGjYMhQzwd1iIZfs1bDclDqrBmjbtuWlDgmqcvvNDd/uCz5ulE1uia7MC6DoqLvDzXZ9eWLe51ly5wxRXQubMnxUmGk6cFJI/s3euGZtm+3b3OyXHfZeuDLuoaZUCKh0YfkMD9qvziC3j7bdeHl4irKV18cdyb8ZLh5GkBKc4KC2HRospro1lZruusgQM9re0nM68DUqCy0gX4ODj8RGD+Zar6Vljb8tOX2gJSFcXF8P777j6Nigp34XfUKDc6bZzuXUqGk6cFpDipqHDXiBYtckEpJcV9V0ePtqSFGPM4qeEu4D+A9cAg4G5VnR9YtlJVh4S1PT99qS0g1WL/fjdU85dfutft27ukh699Lea/OJPh5GkBKQ42b3bXifYFukLLyYHLL3fjFpmY8zggfY7rMuiYuAFZXwKeVdU/Vc3KDnl7vvlSYwGpTqqwaZNrxsvLc/NyclxTSAyvLyXDydMCUgzt3QvvvFP5Y6lNG/ed7NPHmufiyOOAtE6rjAYrlSPFrsPdtzoorO15/qWuwgJSPcrLXdv8+++7ZhGAAQPc9aW2baO+u2Q4eVpAioHDh13T3Jo17sdSRobrNHj4cMue84DHAek94IequqrKvDRgOnCjqoZ1fcECUiIqLHQp4h9/7IJUSoq7yXDkSGgRvdsCkuHkaQEpigoK3Pdu+XL3vUtNdR2hjhwJzZrFtyzmBI8DUlegTFX31rLsAlX9MKztWUBKYEeOuF+qq1e7X6ppae5C8gUXRCUjL6FPnvXs2wJSGI4fh48+gk8+cT2MiLgbukePjknN3ITH6yy7qkTkKlV9LeL3W0BKAvv2ucC0PnBrWHq6C0wjRjQoMCXkyTPEfVtACkFhocvyXLbMZX2CS6YZO9YSFnzEZwFpjYY5BtJJ77eAlET27HGB6V//cq+bNHFNKiNGRNSUl1AnzzD3bQHpFAoK3IitwRoRQM+erkbUtWts9mki5rOA9LmqnhXx+y0gJaFdu1ziQzAwpaW5PvLOPz+sJpaEOHlGuG8LSLU4fNg1zX32GZSWunk9e7prRGecEd19majxWUCyGpKpw549rpv/DRvcNSYR6NfPBaYQhob29cmzgfu2gFTFnj2uRrR2rbvBFVzT3KhRNoR4ArCAFCMWkGLkwAHXeevq1ZUnnOxsOO886N3bZenVwpcnzyjtu9EHpOC9bR99BNu2uXkpKS5Z4YIL4LTTGlxGEx8WkGLEAlKMHT3qUsVXrICiIjevdWuXADF4sOtlvArfnDxjsO9GG5CKilyT3PLlcPCgm9e0qesvcfhw930wCcVnAekdVb0k4vdbQGqEiovdSemTTypPSmlp7ibbc85xzTQi3p88o8ACUsDu3a6vubVrKxMVgj9Ghgyx/uYSmJ8CUkNZQGrMgs02H39c2f0LQMeOMHgwmeedR6EFpLiUJybbKSpyAWjlSheQgrp3d7Whr32tzuZakziSKSBZPx+NmYg7KX3ta66m9OmnsGrVibFs/hPgpZfg7LOhRw87eSWCigrYutVdL1y3zg34CK45dtAg16NH+/beltGYOlhAMk7btq5jzIsvho0bYeVKUsH9wl67Fpo3dxe8BwxwHbpa55n+oep+RHz+uZvy8yuXnXmmC0R9+7r70ozxMWuyM3VqLcLh9993v7aDvYyDC179+7sU8o4d6w1O1mQX/f2KCLp3r6sFrV1b8+8zcKCr2bZp0+B9GX9LpiY7C0imTidOnqruZttgbenYscqV2rRxv7779HF38dfSrGcBKUr7rahwf4eNG7nrkkt45J57KpdlZbkfCAMHur+D1WAbjWQKSNZkZ+on4k5yXbu6Zr3t290v8/Xr4dAhdy/LRx+5Hp979XJTjx410shNBAoLYcsWl3yyaZPr1gdoC+7z7tPHNaPm5Ng1PpPwrIZk6lTvr/mKCti50/UEsXFjZQq5e7NLHz/zTHLGjmVbaaknY+UkXA2pvBxyc11iwpdfuudV12/dGnr3Jmf8eLYFhx4xjZrVkIwBdzLMznbTpZe6HiGCv+R37HAn09xcpgL89reuhpWd7fpF69rV9Ure2JWWus9pxw5X89y5s7IfOaj8jHv2dNmQHTqACNuDy4xJIhaQTHSIuJNlhw6ur7ziYneC3bKFr8CdZLdudVNw/dNPd4Gpc2c3dejgBn1LVhUVnAYutX7XLheIvvqqsjunoNNOc/cKBaemTb0orTFxZ012pk5RzQg7dqyyFrBjh0tTrn4iTktzQen0090UDHAtW0Z8kd6TJjtV101TXp4bq+qrr9y0bx+/+dWv+E3VZAQRl6l4xhmVtcfmzevdhadDmBtfsSY7Y8KVleWy8fr2da9LS10v07m5rheB3bvdNag9e9xUVXq6y+Zr29Y9tmkDrVq5qUULlzwRz6wyVdcLwtGjbtTeo0ddcsehQ+4YDh6s7J6nmoPgUuY7d3a1w06drOnSmICkCUhTHv9njXlXDuzEt0bkUFhSztSnPqmx/NqhXbnunG4cLCjhu8+tqLH8pvOyuerszuw+XMgP5qyqsfw7F53JuH6n8+X+Y/z85c9rLP/exb24sFd7vth9hPteW1dj+Y8v683Q7Las2H6Q3721scbyX1/Vj/6dW7F00wH+/N6mGsv/3+Sz6NGhOQvXfcXflmypsfyPUwbRuXUmr63ezXPLttdY/thNQ2mblc6Ln+7kpRW5NZZLmmsqevaf23h9zZ4ay+f82wgAnvjgS95dv++kZRlNUplx6zAAWp1/fY2/T5tm6fzft4bCGWfw27c2sLJpFzi9wgWqkhI6SQkPtzsABw5w76E2rMtvBjsqgDwgjzNTiniwiTumn5V1Z4s0c819KSmQmkq/Finc068pFwLff3Qhe4oqQFJAAITH6bE2AAAfc0lEQVQzR97sEjFUueO9PRwqLgfFBRtVLmibwl3ZAsXF3Ly6jKLSClejKy+H8nLGphxmWtpeAKYU965yZK2AVlzZrIBvdUmhsN1pTN2UDulNoEk679+QwVcHu3Jt965clx347j1V87tb33cvs4f7bJP1u/f0LcPITE9t8HfvkXc38eHmAyctP/HdA/fd237opOWdWmXw8PWDAbj3tS9Yt/voScvP7JDFg5Ndh9Y/e3kNW/af3KrTr3NL7rmqPwDfn/0Ze44UnbR8SHYbfnJZHwDueHYFh46XnDiexi5pApJJAikp7npJ06bQKgOuv9LNn7cacg+7YFVW5qbUptC8g+uVIL8Cyssqu8kBOH4cDu5kLHB4127Qk3sp6Acwa5Z7UdIDtNq/wuGjsCNwIizpBVotgSAttbI5cWcL19wYnJo0gUFdYEQOlJTD3po/howxNdk1JFOnhOntu7TU3ax7/HjlVFQEJSWMGzmShW+84YJVaamr4ahyw/XXM2vmTNfUJ+KCYVqaq2WlpbmgmJ7uHjMzK6esLHeNJ8JueBLmMzUJI5muIcU0IInIZcCfgFTg76r6P6da3wKSvyTDyTPh7kOK83ZM4kumgBSzGxlEJBV4FLgc10Jyg4j0i9X+jDHGJLZY3lk3DNisqltUtQSYDUyI4f6MMcYksFgmNXQBdlZ5nQsMr76SiEwDpgVeqogURri/NKCs3rUSg2+ORRqWTn3iOBq4nQapa99hlilqf5NofRYN2I5vvl8NlCzHAQ07lqTpNDKWAam2/5Yajd6q+gTwRIN3JvKpqp7T0O34QbIcS7IcB9ix+FGyHAck17E0RCyb7HKBblVedwV217GuMcaYRi6WAWk50EtEuotIOnA98GoM92eMMSaBxazJTlXLRORO4G1c2vd0Vf0iVvsjCs1+PpIsx5IsxwF2LH6ULMcByXUsEfPVjbHGGGMaLxtQxRhjjC9YQDLGGOMLSRWQROR+EVkjIqtEZIGIdPa6TJEQkd+LyIbAscwTkdZelylSInKdiHwhIhUiknBprSJymYhsFJHNIvJTr8vTECIyXUT2ichar8vSECLSTUQWicj6wHfrbq/LFCkRyRCRT0RkdeBY7vW6TF5KqmtIItJSVY8Gnt8F9FPVOzwuVthE5FLgvUBiyG8BVPUnHhcrIiLSF6gAHgd+pKqfelykkAW6v/oXcAnuNoblwA2qWnM8hwQgIiOBY8AzqjrA6/JESkQ6AZ1UdaWItABWABMT8e8i7u7mLFU9JiJNgKXA3aq6zOOieSKpakjBYBSQRS034iYCVV2gqsG7tpfh7uFKSKq6XlVrDriTGJKq+ytV/YDAGIGJTFX3qOrKwPN8YD2uZ5iEo86xwMsmgSkhz1vRkFQBCUBEHhCRncCNwK+9Lk8U3Ar8w+tCNFK1dX+VkCe+ZCUiOcBg4GNvSxI5EUkVkVXAPuAdVU3YY2mohAtIIrJQRNbWMk0AUNVfqGo3YCZwp7elrVt9xxFY5xe4/q1melfS+oVyLAkqpO6vjDdEpDkwF/h+tdaRhKKq5ao6CNcSMkxEErY5taESbsRYVR0X4qrPA28A98SwOBGr7zhE5GbgSmCs+vxCXxh/k0Rj3V/5VOB6y1xgpqq+7HV5okFVD4vIYuAyIKETTyKVcDWkUxGRXlVeXg1s8KosDREY2PAnwNWqetzr8jRi1v2VDwUSAZ4E1qvqQ16XpyFEpEMwi1ZEMoFxJOh5KxqSLctuLtAbl9W1HbhDVXd5W6rwichmoCmQF5i1LBGzBQFEZBLwZ6ADcBhYpapf97ZUoROR8cDDVHZ/9YDHRYqYiMwCRgPtga+Ae1T1SU8LFQERuRBYAnyO+18H+LmqvuldqSIjIgOBGbjvVwrwgqre522pvJNUAckYY0ziSqomO2OMMYnLApIxxhhfsIBkjDHGF3yV9p2SkqKZmUkzPLwxxsTc8ePHVVWTonIRs4AkIt2AZ4COuEyYJ1T1T6d6T2ZmJgUFBbEqkjHGJB0RKfS6DNESyxpSGfCfVTtAFJF3ErEDRGOMMbEXs2peMnWAaJLQu+/C738PuxLuNjVjklZc2h2ToQNEkyT274cbb4Rx4+DHP4acHJgyBZY1yt7+jfGVmN8YG+gA8X3ggdr6nBKRacA0gPT09KHFxcUxLY9puJycHLZv315jfnZ2Ntu2bYt/gUL1yitw221w8CBkZsKoUfDOO1Be7pa//DJMmuRJ0er6TOvi+886TkpLS8nNzaWoqMjrosRcRkYGXbt2pUmTJifNF5HjqprlUbGiKqYBKdAB4uvA26H0OZWVlaWW1OB/IkJt35u65vvC1q3Qrx8UFbna0f/9H/ToAbm58OCD8Ne/Qrt2sHYtdOwY9+KF+9n5+rOOo61bt9KiRQvatWuH6+IuOakqeXl55Ofn071795OWJVNAilmTXTJ1gGgSnCp873suGF1/PSxY4IIRQNeu8Oc/wyWXQF4e3H67W98khKKioqQPRuB+gLRr1y7pa4KxvIZ0AfAt4GIRWRWYxsdwf8bUbv58eOMNaNkS/vhHqH7ySkmBp56C1q3den//uzflNBFJ9mAU1BiOM5ZZdktVVVR1oKoOCkwJ1xuvSXDHjsFdd7nnDzxQd3Ncly6u2Q7gBz9wTXzGmLhKirt7janTfffBzp0wdCh897unXveGG+Daa6GgAP73f+NTPpPU5s2b53UREooFJJO8DhyARx5xTXSPPQapqfW/59e/do8zZsDhw7Etn0lqX375JQsXLgSgsLCQUaNGUR7M6KxFSUkJI0eOpKysLF5F9B0LSCZ5Pf00FBfD5ZfDueeG9p6zzoKLL3a1pOnTY1o8k9zefvttvve97wEwffp0Jk+eTOopfhSlp6czduxY5syZE68i+o4FJJOcKirg8cfd8zvCHGz37rvd45//XHmPkjF12Lx5Mx06dCAnJ4dBgwbRtm1bevToQVlZGX369AFg5syZTJgw4cR7Ro8ezcaNGwHIy8tjwIABAEycOJGZM2fG/yB8wgKSSU7vvgubN0O3bjA+zOTOK66AM8+EbdvgtddiUjwTAyKxmerRs2dPLrzwQp599llWrVrFwIEDeeWVV7grkExTUlLCli1byMnJOfGezZs306tXLwDWrFnDWWedBcCAAQNYvnx59D+bBGEBySSn//s/9zhtWmjXjqpKTXX3LQH86ZQd1BsDwBdffHGilrNhwwZ69+59YtmBAwdo3br1idfbt2+nS5cupKS40++aNWsYOHAgAKmpqaSnp5Ofnx/H0vuHBSSTfHbvdvcepaa6roIiccst0Lw5LF4Ma9ZEtXgmRlRjM9WjsLCQoqIi2rRpw86dO2nXrh3p6eknlmdmZp50Q2uwFhW0YsWKk14XFxeTkZERpQ8lsVhAMsnnySfdtZ+JE6FTp8i20aoVTJ3qnj/1VNSKZpLPunXr6Nu3LwDr168/8TyoTZs2lJeXnwhKq1evPvF806ZNzJ8//0STXV5eHh06dKjRX11jYQHJJJeKCvjb39zzcJMZqrvpJvf44otuu8bUompzXWZmJitXrmTDhg0nrXPppZeydOlSwNWQKioqOPvss7nvvvvo27cvM2bMAGDRokWMD/eaZzJRVd9MzZo1U+N/7msT+vy4+uAD19CSk6NaXt6wbVVUuO2A224MhfvZ+eKz9oF169Z5XYSQrFy5Um+66SZVVe3Ro4cePXq01vUmTZqkGzZsqHM7tR0vUKD1nFuB6cA+YG0dywV4BNgMrAGG1LfNWExWQzLJ5aWX3ON117k+6hpCBL7xDfd89uyGbcs0aoMHD2bMmDEcOXKElJQUWrRoUWOdkpISJk6ceFJCRBQ9DVx2iuWXA70C0zTgsVgUoj4WkEzyqKiAuXPd82uvjc42p0xxjy+9BI34DnrTcLfeeiutWrXiX//6V63L09PT+fa3vx2TfavqB8DBU6wyAXgmUOlaBrQWkQgvwEbOApJJHsuWuSHJzzgj9J4Z6jN4MPTsCfv2wfvvR2ebxvhPF2Bnlde5gXlxZQHJJI9gc92114Z0Q2NIRCprSY24Sxfja2ki8mmVaVoE26jtHybuA4NZQDLJQfXkgBRNwYA0dy6UlkZ328Y0XJmqnlNleiKCbeQC3aq87grsjk7xQmcBySSH5cvdMBNdusDw4dHd9oAB0LcvHDzouiQyvqKNZITfGB/nq8C3xTkPOKKqe2K5w9pYQDLJ4cUX3eO11zY8u666qtl2waQJ4wsZGRnk5eUlfVBSVfLy8iLuwUFEZgH/BHqLSK6I3CYid4hI8Ga9N4EtuLTvvwH/Ho1yh11OP/0hs7KytKCgwOtimHqISK0ngLrmx5xqZWeoS5bAhRdGfx+rVrkEh44dXeJElINeuJ+dZ5+1z5SWlpKbm3tS1zzJKiMjg65du9boxUFEjqtqlkfFiioLSCZsvgtIMQ4WgAt62dmuWfDjj2HYsKhu3gKSiVQyBSRrsjOJ7/XX3eNVV8UmGIFrtrvqKvfchqQwJiYsIJnEFwwQV14Z2/1cfbV7fPXV2O7HmEbKmuxM2HzVZLd3r+vROyMD8vKgWbPY7au4GNq3h2PH3PWq7Oyobdqa7EykrMnOGL944w33OHZsbIMRQNOm8PWvu+fWbGdM1FlAMoktGBiC13diLdhsZwHJmKizJjsTNt802RUVQbt2cPy4y37r2jX2+zxwAE4/3Y1Ge+AAtGwZlc1ak52JlDXZGeMHixa5YDR4cHyCEbhrSOef77oQWrAgPvs0ppGwgGQSV7yb64Is286YmLCAZBKTauX9R7FO964uuL9//APKy+O7b2OSmAUkk5jWrHHXjTp2hKFD47vvPn1cV0UHDsAnn8R338YkMQtIJjFVrR3FqneGuohU1pKCaefGmAazgGQSk1fNdUFXXHFyOYxppEQkS0RSo7ItP6WOWtp3YvA87XvfPtdU16SJ652hefPY77O64mKXcl5QADt2QLdu9b/nFCzt20Qq3mnfIpICXA/cCJwLFANNgf24YSyeUNVNkWzbakgm8fzjHy6pYcwYb4IRuF4bLrnEPX/zTW/KYIw3FgE9gJ8BHVW1m6qeBlwELAP+R0RuimTDFpBM4vG6uS4ouH9rtjONyzhVvV9V16hqRXCmqh5U1bmqeg0wJ5INW5OdCZunTXYlJe7m1Px82LIFuneP7f5OZc8e6NwZMjNd02FmZsSbsiY7EymvemoQkc+BNVWmz4GbVfWBSLdpNSSTWJYudcGof39vgxG4XsaHDoXCQtdrhDGNyyjccOeFuGtKa4ErGrJBC0gmsfiluS4oWA7rbNU0MoEmusWq+oiq3oxLcIgomSHIApJJLH4LSFVHkbUmNNOIiEivqq8DmXUDG7LNmAUkEZkuIvtEZG2s9mEamY0bYdMmaNsWzjvP69I4Q4a4jl137YIVK7wujTHx9ISI7BCRf4rI4yIyA1grIhEPTBbLGtLTwGUx3L5pbObPd4+XXw5pad6WJUiksrPVV17xtizGxJGqjlHVM4ApwBvAZqAZsFpENkSyzZgFJFX9ADgYq+2bRmjePPc4aZK35ahuwgT3GAyYxiQxEZGqr1V1h6q+GkgFv0ZVewHDItm2XUMyiWH3bli2DDIy4DKfVbxHj3YD9a1d61LRjUlui0TkeyJyRtWZIpIuIhcHmu4mR7JhzwOSiEwTkU9F5NOysjKvi2OqyMnJQURqTNnZ2bWun52dXev6OTk5DS9McOyhSy+FrNBvuajrGKJSpqD0dBg/3j2vp5YU7mdqjA9dBpQDs0Rkt4isE5GtuAy7G4A/qurTkWw4pjfGikgO8LqqDghlfbsx1l+idfNlVLbz9a+7EVqnT4dbbmnwvqN+Y+ns2XDDDTByJLz/ftjlCZfdGGuCvBzCXESaAO2BQlU93ODtWUAydfHNyfPwYejQASoq4KuvXE8NDdx31E/oR464MpaXn7KMvvlMTdKIZ0ASkYnAWODPqvqvaG8/lmnfs4B/Ar1FJFdEbovVvkySe/NNKCtztY8wglFctWrlriVVVNgYScZ3ROQyEdkoIptF5Ke1LJ8qIvtFZFVgur227ajqK8AfgbEicme0yxnLLLsbVLWTqjZR1a6q+mSs9mWSXDCdeuJEb8tRn2D5gtmAxvhAYKyiR4HLgX7ADSLSr5ZV56jqoMD097q2p6pbVPUx4G4R6SciJ+7BEJGvNaSsnic1GHNKRUVuuAlIjIAkAm+95ZrwjPGHYcDmQCApAWYDE6Kw3ba4QLdTRD4XkRdw4yFFzAKS8bcFC+DYMdcjgt8z0Tp3hlGj3OB9dpOs8Y8uwM4qr3MD86q7RkTWiMhLIhLKiJM7AjfHdsLVvv4GTG1IQS0gGX97/nn3eN113pYjVNdf7x5nz/a2HKYxSQveOhOYplVbLrW8p3pGzGtAjqoOBBYCM0LYbysROV9EWqlqrqq+o6pLIyh/ZUH9lKljWXb+4nlGWH4+nH66G95h27aIakhxH7vpwAE3vDrA3r01kjA8/0xN0qkvy05ERgC/UdWvB17/DEBVH6xj/VTgoKq2qme/+4ClQH/cEObrgbWq+l8RHQhWQzJ+9sorLhhddJH/m+uC2rd3Q5uXl8PcuV6XxhiA5UAvEekuIum4sYterbqCiHSq8vJqXHCpzwhVnayqvYE+wC9wA/VFzAKS8a+ZM93jN7/pbTnCZc12xkdUtQy4E3gbF2heUNUvROQ+EQn0DMxdIvKFiKwG7iKEa0Gq+mWV50WqulJVn21IWa3JztTJ0+alr76CLl1c1trevdCuXVT3HdMmryNHXFNjSQns3OmOI8r7tSY7E+RlTw3RZjUk408vvOCavS6/POJg5JlWrVzfdqrw4otel8aYhGEByfhTojbXBQWb7YLHYYyplzXZmTp51rz05ZfQsyc0b+6a7ppFPAClN0124JIxOnd2/fCtWOHuo4rifq3JzgRZk50xsTR9unucNKlBwchTmZlw883u+eOPe1sWYxKE1ZBMnTz5NV9UBN26uft5PvwQzj8/JvuOSw1j/Xro18+N37R7N7RsaTUkE3VWQzImVl580QWjwYNhxAivS9Mwffu6roQKCuxakjEhsIBk/OUvf3GPd97pUr4T3R13uMfHHnNZd8aYOlmTnalT3JuXPvkEhg+Htm0hN9ddh4nRvuPW5FVc7Jog9++Hjz5Czj/fmuxMVFmTnTGx8Oij7vG226ISjHyhadPKIdcfe8zbshjjc1ZDMnWKaw1p/37o2hVKS13ad/fuDd7vqfYd1xrGl19Cr16QlsYZpaXssBqSiSKrIRkTbX/5i+tq54orohaMfKNHD3ejbGkpP/O6LMb4mNWQTJ3iVkPavx/OPNMNxLdkCVx4YYP3Wd++417DWL8e+venRJX0CIfSqMpqSCbIakjGRNODD7pgNH58VIORr/TtCzfcQDq44zXG1GA1JFOnuNSQdu5011eKi2HlSnf/URT5poYEsGED5X37ktqkCWza1KBaktWQTJDVkIyJlnvvdcFoypSoByPf6dOHWeASNx54wOvSGOM7VkMydYp5DWnjRujf3z1fv97VlKLMVzUkoLcIG1NS3E2/H38MQ4dGtB2rIZkgqyEZ01AVFfDd77oxj269NSbByI/+BXD33ZXHXVLidZGM8Q0LSMYbjzwCixZBhw7w3//tdWni6/77XVbhmjXw2996XRpjfMOa7EydYtZkt26dGx+ouBjmz4err27wPkLedz3zY+3Eft97D8aOhSZN4LPPKpsuw92OafSsyc6YSJWWwre/7YLRrbfGNBj52sUXw3e+4z6PW25xn4cxjZwFJBM/qvCjH7kRVHNy4I9/9LpE3vr9713Hq8uXu8H8Kiq8LpExnrKAZOLn3nvdtaMmTeCZZ6BlS69L5K1WreC116BFC5gzB378Y69LZIynLCCZ+PjDH1xASkmBWbPgoou8LpE/nH02vPwypKW5z+jhh70ukTGesYBkYkuV/wLXVAfw1FNwzTVelsh/xo1znwvAD34AP/kJlJV5WyZjPGABycTOgQNw1VX8Lvj6L39xCQ2mpptugr/+FVJT4Xe/g8svd5+fMY2IBSQTfarw+uswaBC88QYHAebNg//4D69L5m/f/S4sXOjuzVq40KXGP/+8JTuYRsMCkokeVXj1VRg2DK66CnbtggsuYBDAxIlely4xjB7tshCHDXMdz954o+te6I03XO8OxiQxC0imYcrLYelS+M//dAPRTZgAn34Kp5/uLtIvXsxOr8uYaLp1gw8/hCefdKPorloFV17p5n//+/DRRzTxuowmoYjIZSKyUUQ2i8hPa1neVETmBJZ/LCI58S+l9dRgTuFEbwAlJXD4MOzZ436179jhelv47DNYvRqq/s06d3bpy9/5DjRrdvJ2vDyGEOd7VZ46FRbCo4/CY4/Bli0nZhcBGeedB+eeCz17ulF2s7Ndc1+7dpCeHv3CG1+qr6cGEUnFdaN4CZALLAduUNV1Vdb5d2Cgqt4hItcDk1R1SoyLXrOsCR+QNm2qzOBqbOr621WdX/15cKqoOHkqLa2ciouhsJA9W7fSqVkzOH781OXo0cM1yU2cCCNGuAvzVVhAisJ+Vd0NtLNnw5tvup7ST6VZM2jeHLKyIDMTmjZ1QapJE/f3SU11KfjBnseDj66Qlc+Dr+s+oPCPxdTu0UddjThMIQSkEcBvVPXrgdc/A1DVB6us83ZgnX+KSBqwF+igcf4nSYvnzmLi8GF33cJEXSdwwSg1Fdq0cb++zzjDNR317Okuug8eDO3be13U5CfirisNGwYPPUQbEQ4tWOBqqVu3wrZtruZ64AAcPOj+bvX9kDD+EnlHu2ki8mmV10+o6hNVXneBk1rOc4Hh1bZxYh1VLRORI0A7IK6pnjENSCJyGfAnIBX4u6r+T9R30rMnvPJK1DebMOr6hVrXL9yqv4ZTUlywEXG/nJs0cTdoZmRAZiade/Zkd36++5Vtv4R95TDAJZe4qTpVNyR8QYELSgUFrtm1tNQ9lpdXTlVrzMH31lXDrm0/Jnq6dIn0nWWqes4pltf2z1v9jxfKOjEXs4AUaLd8lCrtliLyatV2y6ho08ZdSDdRtwdcs49JLCKuO6IWLbwuifGHXKBbldddgd11rJMbaLJrBe6OjXiKZZbdMGCzqm5R1RJgNmCRwxhj4ms50EtEuotIOnA9UP06x6vAzYHn1wLvxfv6EcQ2INXWbhlxndQYY0z4VLUMuBN4G1gPvKCqX4jIfSISHP/lSaCdiGwGfgjUSA2Ph1heQwqpTVJEpgHTgstFpDDC/aUBydIBmG+ORRp27ejEcTRwOw1S177DLFPU/ibR+iwasB3ffL8aKFmOAxp2LJn1raCqbwJvVpv36yrPi4DrItx/1MQyIIXSbkkgG+SJ6vPDJSKf1nNhL2Eky7Eky3GAHYsfJctxQHIdS0PEsskulHZLY4wxBohhDSmQyx5st0wFpqvqF7HanzHGmMQW0/uQamu3jKEGN/v5SLIcS7IcB9ix+FGyHAck17FEzFddBxljjGm8rLdvY4wxvpBUAUlE7heRNSKySkQWiEhnr8sUCRH5vYhsCBzLPBFp7XWZIiUi14nIFyJSISIJl0VUX7f9iUREpovIPhFZ63VZGkJEuonIIhFZH/hu3e11mSIlIhki8omIrA4cy71el8lLSdVkJyItVfVo4PldQD9VvcPjYoVNRC7F3SldJiK/BVDVn3hcrIiISF+gAngc+JGqflrPW3wjlG77E4mIjASOAc+o6gCvyxMpEekEdFLVlSLSAlgBTEzEv4u4m8myVPWYiDQBlgJ3q+oyj4vmiaSqIQWDUUAWHnQOGA2quiBwdzXAMtw9XAlJVderaj1jJfhWUnV/paof4EH/ZNGmqntUdWXgeT6u94GE7AVGnWOBl00CU0Ket6IhqQISgIg8ICI7gRuBX9e3fgK4FfiH14VopKz7K58LjGw6GPjY25JETkRSRWQVsA94R1UT9lgaKuECkogsFJG1tUwTAFT1F6raDZiJ67/Jl+o7jsA6v8B1JzLTu5LWL5RjSVC+6JLf1E5EmgNzge9Xax1JKKparqqDcC0hw0QkYZtTGyrhBuhT1XEhrvo88AZwTwyLE7H6jkNEbgauBMZ60etuOML4mySakLq/MvEXuN4yF5ipqi97XZ5oUNXDIrIYuAxI6MSTSCVcDelURKRXlZdXAxu8KktDBAY2/AlwtarasJ/ese6vfCiQCPAksF5VH/K6PA0hIh2CWbQikgmMI0HPW9GQbFl2c4HeuKyu7cAdqrrL21KFL9AFfFMgLzBrWSJmCwKIyCTgz0AH3ECnq1T1696WKnQiMh54mMrurx7wuEgRE5FZwGigPfAVcI+qPulpoSIgIhcCS4DPcf/rAD8P9AyTUERkIDAD9/1KwQ0NcZ+3pfJOUgUkY4wxiSupmuyMMcYkLgtIxhhjfMECkjHGGF+wgGSMMcYXLCAZY4zxBQtIxhhjfMECkjHGGF+wgGTMKYjIYhHpHXjeLtHHEjLGzywgGXNqPYFNgecDcb0DGGNiwAKSMXUQkWxgl6oGu6cZCKzxsEjGJDULSMbUbRAnB6ChWEAyJmYsIBlTt7OBDDjRk/wErMnOmJixgGRM3QYBKSKyGjf68HrgZm+LZEzyst6+jalDYBiQwaqa73VZjGkMrIZkTC1EpAVQYcHImPixGpIxxhhfsBqSMcYYX7CAZIwxxhcsIBljjPEFC0jGGGN8wQKSMcYYX7CAZIwxxhcsIBljjPEFC0jGGGN84f8DhdGjqFldka0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "nsize = 10\n",
    "xs    = stats.norm(0., 1.).rvs(nsize)\n",
    "plt.subplot(2, 1, 1); htexam.normal_likelihood(xs)\n",
    "plt.subplot(2, 1, 2); htexam.normal_posterior(xs);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Exercise:** Show that the posterior of a normal prior, $\\pi(\\mu | \\mu_0, \\sigma_{\\mu_0})$,\n",
    "and a normal likelihood, $p(x |\\mu, \\sigma)$, for n-measurements, $x$, is a normal distribution with mean, $\\mu'$, and sigma, $\\sigma_{\\mu'}$:\n",
    "\n",
    "$$\n",
    "\\frac{1}{\\sigma^2_{\\mu'}} = \\frac{n}{\\sigma^2} + \\frac{1}{\\sigma^2_{\\mu_0}}, \n",
    "\\;\\;\\;\n",
    "\\mu' = \\sigma^2_{\\mu'} \\left( \\frac{\\mu_0}{\\sigma^2_{\\mu_0}} + \\frac{\\sum_i x_i}{\\sigma^2} \\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Exercise:** There are several dices in a box, with 4, 6, 12 and 24 sides. We pick one randomly and we roll it four times, with outcomes, $1,4,5,2$, what is the posterior probability that the selected dice has 4, 6, 12 or 24 sides?, what is the posterior probability if we roll it twice again and wet get 6 and 1? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data : 1  posteriors:  [0.375      0.25       0.16666667 0.20833333]\n",
      "data : 4  posteriors:  [0.59340659 0.26373626 0.08791209 0.05494505]\n",
      "data : 5  posteriors:  [0.         0.82051282 0.13675214 0.04273504]\n",
      "data : 2  posteriors:  [0.         0.91211401 0.0760095  0.01187648]\n"
     ]
    }
   ],
   "source": [
    "measurements = [1, 4, 5, 2]\n",
    "hpriors      = [0.15, 0.15, 0.20, 0.50]\n",
    "for meas in measurements:\n",
    "    hpriors = htexam.dice_posterior(meas, hpriors) \n",
    "    print('data :', meas, ' posteriors: ', hpriors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 5. Expectation values\n",
    "\n",
    "Given a $x$ rv that follows a pdf, $\\, g(x)$, and a function $f(x)$ on $x$, we define the expected value of $g(x)$ as:\n",
    "\n",
    "$$\n",
    "E[g(x)] \\equiv \\int f(x) \\, g(x) \\, \\mathrm{d}x\n",
    "$$\n",
    "\n",
    "The **mean**, or average value, is the expected value of $x$:\n",
    "\n",
    "$$\n",
    "E[x] = \\mu \\equiv  \\int x \\, g(x) \\, \\mathrm{d}x\n",
    "$$\n",
    "\n",
    "The **variance** is the expected value of $(x-\\mu)^2$:\n",
    "\n",
    "$$\n",
    "V[x] = \\sigma^2 = E[(x-\\mu)^2]= E[x^2]-\\mu^2 \\equiv \\int (x-\\mu)^2 \\, g(x) \\, \\mathrm{d}x\n",
    "$$\n",
    "\n",
    "We call **standard deviation** to:\n",
    "\n",
    "$$\n",
    "\\sigma = \\sqrt{V[x]}\n",
    "$$\n",
    "\n",
    "The standard deviation is related with the spread of the pdf, while the mean is related with the center of density of the pdf. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "\n",
    "The **expected value** is the $x$ value with the highest probability $g(x)$.\n",
    "\n",
    "The **median**, is the $x_{med}$ value that divides the distribution in half, \n",
    "\n",
    "$$\n",
    "\\int_0^{x_{med}} g(x) \\mathrm{d}x = 0.5\n",
    "$$\n",
    "\n",
    "For a symmetric pdf, the mean and median are the same.\n",
    "\n",
    "\n",
    "If ${\\bf x}$ is a vector, the pdf is a n-dimensional function. In that case, we define the covariance element between $i$ and $j$ elements as:\n",
    "\n",
    "$$\n",
    "\\mathrm{cov}[x_i,x_j] = \\int x_i \\, x_j \\, g({\\bf x}) \\; \\Pi_{i=1,n} \\mathrm{d}x_i\n",
    "$$\n",
    "\n",
    "If the variables $x_i, \\, x_j$ variables are independent then the covariance is zero. But the contrary is not necessarily true.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 6.  Common probability density functions \n",
    "\n",
    "The most commont distributions in HEP are:\n",
    "\n",
    "* **Binomial**. A event can happen with probability $p$. The probability to get $n$ events if we try $N$ times is given by the binomial pdf.\n",
    "\n",
    "$$\n",
    "g(n|N,p) = \\frac{N!}{n! (N-n)!} p^n (1-p)^{N-n}\n",
    "$$\n",
    "\n",
    "* **Poisson**. Related with the number of events if we expect $\\nu$. The number of interactions in a crossing of the LHC beams is modeled with a poisson.\n",
    "\n",
    "$$ \n",
    "g(n|\\nu) = \\frac{\\nu^n}{n!} e^{-\\nu}\n",
    "$$ \n",
    "\n",
    "* **Uniform**. Equal probability to get $x$ in an interval $[a, b]$.\n",
    "\n",
    "$$\n",
    "g(x|a,b) = \\frac{1}{(b-a)}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "\n",
    "* **Exponential**. An event can happen in $x$ with probability $1/\\tau$. It describes the time of particles decays.\n",
    "\n",
    "$$\n",
    "g(x|\\tau) = \\frac{1}{\\tau} e^{-x/\\tau}\n",
    "$$\n",
    "\n",
    "* **Gaussian** (or normal). Associated with the distributions of measurements.\n",
    "\n",
    "$$ \n",
    "g(x|\\mu,\\sigma) = \\frac{1}{\\sqrt{2\\pi} \\sigma} e^{-\\frac{(x-\\mu)^2}{2\\sigma^2}}\n",
    "$$ \n",
    "\n",
    "* **$\\chi^2_n$** with $n$ degrees of freedom. Associated with goodness of fit.\n",
    "\n",
    "$$\n",
    "g(x|n) = \\frac{1}{2^{n/2} \\Gamma(n/2)} x^{n/2-1}e^{-x/2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "\n",
    "* **Breit-Wigner**. Describes the distribution of masses in resonances. \n",
    "\n",
    "$$\n",
    "g(x|\\Gamma,x_0) = \\frac{1}{\\pi} \\frac{\\Gamma/2}{\\Gamma^2/4 + (x-x_0)^2} \n",
    "$$\n",
    "\n",
    "* **Beta**. Associated with the measurement of probabilities. Used in Bayesian statistics.\n",
    "\n",
    "$$\n",
    "g(x | \\alpha,\\beta) = \\frac{\\Gamma(\\alpha+\\beta)}{\\Gamma(\\alpha)\\Gamma(\\beta)} x^{\\alpha-1} (1-x)^{\\beta-1}\n",
    "$$\n",
    "\n",
    "where $x, \\alpha, \\beta$ are positive.\n",
    "* **Gamma**. Associated with bayesian statistics with poisson distributions. \n",
    "\n",
    "$$\n",
    "g(x | \\alpha,\\beta) = \\frac{\\beta^\\alpha}{\\Gamma(\\alpha)} x^{\\alpha-1} e^{-x \\beta}\n",
    "$$\n",
    "\n",
    "where $x, \\alpha, \\beta$ are positive.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "\n",
    "Some of the pdfs are nicely related:\n",
    "\n",
    "In the case that an event is rare but we do a large number of trials, their distribution follow a binomial with $p \\to  0$, $N \\to \\infty$, if $ N \\, p = \\nu$ is constant, (this happens in nuclei decays), the binomial transform into Poisson with mean $\\nu$. \n",
    "\n",
    "When $\\nu$ is \"large\", 12 is enough, the Poisson transforms into Gaussian with $\\mu=\\nu$.\n",
    "\n",
    "When we take n values of $x$, gaussian distributed, with mean $\\mu$ and sigma $\\sigma$, and compute $\\chi^2 = \\sum_{i=1,n} \\frac{(x_i -\\mu)^2}{\\sigma^2}$, it follows a chi-squared distribution with $n$ degrees of freedom.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Exercise:** Compare a poisson distribution, with $\\nu=p \\, N$ with a binomial when $N$ is large and $p$ is small."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEKCAYAAADuEgmxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAG2hJREFUeJzt3XuUlPWd5/H3B+SyRI0oxEkEBBI0gHREWsJGgyhq8Ioz0RXjBUQP0ZUko5Nzxjnuagc3idEzesIuRklC1hgVE7MZOxsMYzQkx41EGsQGVIaWaaXF9QIG13ht+O4f9TSnKKqbX0E/Xd3yeZ3Tp5/L71f15anq+vDcfqWIwMzMLEWvahdgZmY9h0PDzMySOTTMzCyZQ8PMzJI5NMzMLJlDw8zMkjk0zMwsmUPDzMySOTTMzCzZAdUuoLMMGjQohg8fXu0yzMx6lJUrV74REYNT239kQmP48OE0NDRUuwwzsx5F0ouVtPfhKTMzS+bQMDOzZA4NMzNL9pE5p2Fm+4cPP/yQlpYW3nvvvWqX0qP079+fIUOG0KdPn316HIeGmfUoLS0tHHTQQQwfPhxJ1S6nR4gItmzZQktLCyNGjNinx/LhKTPrUd577z0OO+wwB0YFJHHYYYd1yt5ZrqEhaZqk9ZKaJF1fZv1VktZIWi3pCUljitb9U9ZvvaQv5VmnmfUsDozKddY2yy00JPUGFgBnAGOAi4pDIXN/RIyLiGOBW4Hbs75jgBnAWGAacGf2eGZmVkV5ntOYCDRFxEYASYuB6cCzbQ0i4q2i9h8D2r6wfDqwOCLeB/5dUlP2eE/mWK+Z9UB1y+o69/Gm7PnxmpubOfvss1m7du0uy6+88kquu+46xowp/f9x5znzzDO5//77OeSQQ9pt03az86BBgzr9+fMMjSOATUXzLcDnSxtJuga4DugLnFLUd3lJ3yPyKdP2N5V+yKR8iJgB/OhHP8r9OZYsWZL7c3Qkz3Ma5Q6gxW4LIhZExKeBfwT+SyV9Jc2R1CCp4fXXX9+nYs3MKtHa2srMmTOpqanh/PPP55133mHKlCk7hzM68MADueGGG/jc5z7HpEmTePXVVwF48cUXmTp1KjU1NUydOpWXXnoJgFmzZnH11Vdz8sknM3LkSP7whz8we/ZsRo8ezaxZs3Y+7/Dhw3njjTcAOO+885gwYQJjx45l4cKFXfLvzjM0WoChRfNDgM0dtF8MnFdJ34hYGBG1EVE7eHDyeFtmZvts/fr1zJkzh8bGRg4++GDuvPPOXdb/9a9/ZdKkSTzzzDNMnjyZH/7whwDMnTuXyy67jMbGRi6++GK+/vWv7+zz5ptv8vjjj3PHHXdwzjnncO2117Ju3TrWrFnD6tWrd6th0aJFrFy5koaGBubPn8+WLVvy/UeTb2isAEZJGiGpL4UT2/XFDSSNKpo9C9iQTdcDMyT1kzQCGAU8lWOtZmYVGTp0KCeccAIAl1xyCU888cQu6/v27cvZZ58NwIQJE2hubgbgySef5Ctf+QoAl1566S79zjnnHCQxbtw4Dj/8cMaNG0evXr0YO3bszv7F5s+fv3NPZtOmTWzYsGG3Np0tt3MaEdEqaS6wFOgNLIqIdZLmAQ0RUQ/MlXQq8CHwJjAz67tO0s8pnDRvBa6JiO151WpmVqnSS1hL5/v06bNzWe/evWltbd3j4/Tr1w+AXr167Zxumy/tv2zZMn73u9/x5JNPMmDAAKZMmdIld8nnep9GRCyJiKMi4tMR8e1s2Y1ZYBAR34iIsRFxbEScHBHrivp+O+t3dEQ8kmedZmaVeumll3jyycIFnQ888AAnnnhiUr8vfOELLF68GID77rsvuV+pbdu2MXDgQAYMGMDzzz/P8uXL99ypE3gYETPr0ap1ddvo0aO55557+OpXv8qoUaO4+uqr+fWvf73HfvPnz2f27NncdtttDB48mJ/85Cd79fzTpk3jrrvuoqamhqOPPppJkybt1eNUShG7XZTUI9XW1oa/hMlS+JLbnu25555j9OjR1S6jRyq37SStjIja1Mfw2FNmZpbMoWFmZskcGmZmlsyhYWZmyRwaZmaWzKFhZmbJfJ+GmfVsjXWd+3g1nfx4dM2Q6V3FoWFmlrOuGDK9q/jwlJlZhZqbm/nsZz+729Dojz32GOPHj2fcuHHMnj2b999/H2DnkOnbt29n1qxZHHPMMYwbN4477rgDKNwlPmbMGGpqapgxYwYAW7du5bzzzqOmpoZJkybR2NgIQF1dHbNnz2bKlCmMHDmS+fPnd+m/3aFhZrYXSodGv/3225k1axYPPvgga9asobW1lR/84Ae79Fm9ejUvv/wya9euZc2aNVx++eUA3HLLLTz99NM0NjZy1113AXDTTTcxfvx4Ghsb+c53vsNll12283Gef/55li5dylNPPcW3vvUtPvzwwy77dzs0zMz2QunQ6I899hgjRozgqKOOAmDmzJn88Y9/3KXPyJEj2bhxI1/72tf47W9/y8EHHwxATU0NF198MT/72c844IDCWYMnnniCSy+9FIBTTjmFLVu2sG3bNgDOOuss+vXrx6BBg/jEJz6x8wueuoJDw8xsL5QOhZ5i4MCBPPPMM0yZMoUFCxZw5ZVXAvCb3/yGa665hpUrVzJhwgRaW1spNy5g23MWD5ve0bDreXBomJnthdKh0U899VSam5tpamoC4N577+Wkk07apc8bb7zBjh07+PKXv8zNN9/MqlWr2LFjB5s2beLkk0/m1ltv5S9/+Qtvv/02kydP5r777gMK350xaNCgnXsm1eSrp8ysZ8vhEtkUpUOjf//732fSpElccMEFtLa2cvzxx3PVVVft0ufll1/m8ssvZ8eOHQB897vfZfv27VxyySVs27aNiODaa6/lkEMOoa6ujssvv5yamhoGDBjAPffcU41/5m4cGmZme6FXr147T1q3mTp1Kk8//fRubZctW7ZzetWqVbutL/2qWIBDDz2Uhx9+eLfldXV1u8yvXbs2seLO4dAwq4C/i8P2dz6nYWZWoeHDh3f5//C7C4eGmfU4H5VvHO1KnbXNHBpm1qP079+fLVu2ODgqEBFs2bKF/v377/Nj+ZyGmfUoQ4YMoaWlhddff73apfQo/fv3Z8iQIfv8OA4NM+tR+vTpw4gRI6pdxn7Lh6fMzCyZQ8PMzJI5NMzMLFmuoSFpmqT1kpokXV9m/XWSnpXUKOkxSUcWrdsuaXX2U59nnWZmlia3E+GSegMLgNOAFmCFpPqIeLao2dNAbUS8I+lq4FbgwmzduxFxbF71mZlZ5fLc05gINEXExoj4AFgMTC9uEBG/j4h3stnlwL5fD2ZmZrnJMzSOADYVzbdky9pzBfBI0Xx/SQ2Slks6r1wHSXOyNg2+ZtvMLH953qdR7htKyt7CKekSoBYoHnx+WERsljQSeFzSmoh4YZcHi1gILASora317aFmZjnLc0+jBRhaND8E2FzaSNKpwA3AuRHxftvyiNic/d4ILAPG51irmZklyDM0VgCjJI2Q1BeYAexyFZSk8cDdFALjtaLlAyX1y6YHAScAxSfQzcysCnI7PBURrZLmAkuB3sCiiFgnaR7QEBH1wG3AgcAvsu++fSkizgVGA3dL2kEh2G4puerKzMyqINexpyJiCbCkZNmNRdOnttPvT8C4PGszM7PK+Y5wMzNL5tAwM7NkDg0zM0vm0DAzs2QODTMzS+bQMDOzZA4NMzNL5tAwM7NkDg0zM0vm0DAzs2QODTMzS+bQMDOzZA4NMzNL5tAwM7NkDg0zM0vm0DAzs2QODTMzS+bQMDOzZA4NMzNL5tAwM7NkDg0zM0vm0DAzs2QODTMzS+bQMDOzZA4NMzNLlmtoSJomab2kJknXl1l/naRnJTVKekzSkUXrZkrakP3MzLNOMzNLk1toSOoNLADOAMYAF0kaU9LsaaA2ImqAh4Bbs76HAjcBnwcmAjdJGphXrWZmlibPPY2JQFNEbIyID4DFwPTiBhHx+4h4J5tdDgzJpr8EPBoRWyPiTeBRYFqOtZqZWYIDcnzsI4BNRfMtFPYc2nMF8EgHfY/o1Oqsx6pbVld5nymV9zGz3eUZGiqzLMo2lC4BaoGTKukraQ4wB2DYsGF7V6WZmSXL8/BUCzC0aH4IsLm0kaRTgRuAcyPi/Ur6RsTCiKiNiNrBgwd3WuFmZlZenqGxAhglaYSkvsAMoL64gaTxwN0UAuO1olVLgdMlDcxOgJ+eLTMzsyrK7fBURLRKmkvhw743sCgi1kmaBzRERD1wG3Ag8AtJAC9FxLkRsVXSzRSCB2BeRGzNq1YzM0uT5zkNImIJsKRk2Y1F06d20HcRsCi/6szMrFK+I9zMzJI5NMzMLJlDw8zMkjk0zMwsmUPDzMySOTTMzCyZQ8PMzJI5NMzMLJlDw8zMkjk0zMwsmUPDzMySOTTMzCyZQ8PMzJI5NMzMLJlDw8zMkiV/n0b2DXqfAt4FmiNiR25VmZlZt9RhaEj6OHANcBHQF3gd6A8cLmk5cGdE/D73Ks3MrFvY057GQ8BPgS9GxF+KV0iaAFwqaWRE/DivAs3MrPvoMDQi4rQO1q0EVnZ6RWZm1m0lnQiXdEXJfG9JN+VTkpmZdVepV09NlbRE0iclHQMsBw7KsS4zM+uGkq6eioivSLoQWAO8A1wUEf8n18rMzKzbST08NQr4BvBLoJnCCfABOdZlZmbdUOrhqV8D/zUivgqcBGwAVuRWlZmZdUupN/dNjIi3ACIigH+WVJ9fWWZm1h11uKch6USAtsAoFhEbJB2cnRhvr/80SeslNUm6vsz6yZJWSWqVdH7Juu2SVmc/Digzs25gT3saX5Z0K/BbCvdktN0R/hngZOBI4B/KdZTUG1gAnAa0ACsk1UfEs0XNXgJmAd8s8xDvRsSx6f8UMzPL255u7rs2G3PqfOAC4G8ojD31HHDXHq6gmgg0RcRGAEmLgenAztCIiOZsncexMjPrAfZ4Ijwi3gQOBhqBR4EngDeAz0rqaE/gCGBT0XxLtixVf0kNkpZLOq+CfmZmlpPUE+ETgFqgHhBwFoWrp66S9IuIuLVMH5VZFhXUNiwiNksaCTwuaU1EvLDLE0hzgDkAw4YNq+Chzcxsb6RecnsYcFxEfDMi/oFCgAwGJlM4J1FOCzC0aH4IsDm1sIjYnP3eCCwDxpdpszAiaiOidvDgwakPbWZmeyk1NIYBHxTNfwgcGRHvAu+302cFMErSCEl9gRkU9lT2SNJASf2y6UHACRSdCzEzs+pIPTx1P7Bc0sPZ/DnAA5I+Rjsf5hHRKmkusBToDSyKiHWS5gENEVEv6XjgV8BA4BxJ34qIscBo4O7sBHkv4JaSq67MzKwKUseeulnSEuBECucqroqIhmz1xR30WwIsKVl2Y9H0CgqHrUr7/QkYl1KbmZl1neSve/X3Z5iZWeo5DTMzs/Q9DbPuYspby6ra32x/5j0NMzNL5tAwM7NkDg0zM0vm0DAzs2QODTMzS+bQMDOzZA4NMzNL5tAwM7NkDg0zM0vm0DAzs2QeRsSsAh6CxPZ3Dg2rjsa6ytrXVNjezHLhw1NmZpbMoWFmZskcGmZmlsyhYWZmyRwaZmaWzKFhZmbJHBpmZpbMoWFmZskcGmZmlsx3hJt1kbpldZX3mVJ5H7M8eU/DzMyS5RoakqZJWi+pSdL1ZdZPlrRKUquk80vWzZS0IfuZmWedZmaWJrfQkNQbWACcAYwBLpI0pqTZS8As4P6SvocCNwGfByYCN0kamFetZmaWJs89jYlAU0RsjIgPgMXA9OIGEdEcEY3AjpK+XwIejYitEfEm8CgwLcdazcwsQZ6hcQSwqWi+JVvWaX0lzZHUIKnh9ddf3+tCzcwsTZ5XT6nMsujMvhGxEFgIUFtbm/rYZlXhL3Cyj4I89zRagKFF80OAzV3Q18zMcpJnaKwARkkaIakvMAOoT+y7FDhd0sDsBPjp2TIzM6ui3EIjIlqBuRQ+7J8Dfh4R6yTNk3QugKTjJbUAFwB3S1qX9d0K3EwheFYA87JlZmZWRbneER4RS4AlJctuLJpeQeHQU7m+i4BFedZnZmaV8TAitvca6yprX1NhezPrdjyMiJmZJfOehlXFsuZlFbWfUpNPHWZWGe9pmJlZMoeGmZklc2iYmVkyh4aZmSVzaJiZWTKHhpmZJXNomJlZMoeGmZklc2iYmVkyh4aZmSXzMCJmPYUHiLRuwHsaZmaWzKFhZmbJHBpmZpbMoWFmZskcGmZmlsyhYWZmyRwaZmaWzKFhZmbJHBpmZpbMoWFmZsk8jMj+zMNSmFmFct3TkDRN0npJTZKuL7O+n6QHs/V/ljQ8Wz5c0ruSVmc/d+VZp5mZpcltT0NSb2ABcBrQAqyQVB8RzxY1uwJ4MyI+I2kG8D3gwmzdCxFxbF71mZlZ5fLc05gINEXExoj4AFgMTC9pMx24J5t+CJgqSTnWZGZm+yDP0DgC2FQ035ItK9smIlqBbcBh2boRkp6W9AdJX8yxTjMzS5TnifByewyR2OYVYFhEbJE0AfgXSWMj4q1dOktzgDkAw4YN64SSzcysI3mGRgswtGh+CLC5nTYtkg4APg5sjYgA3geIiJWSXgCOAhqKO0fEQmAhQG1tbWkgmVkbXylnnSTPw1MrgFGSRkjqC8wA6kva1AMzs+nzgccjIiQNzk6kI2kkMArYmGOtZmaWILc9jYholTQXWAr0BhZFxDpJ84CGiKgHfgzcK6kJ2EohWAAmA/MktQLbgasiYmtetZqZWZpcb+6LiCXAkpJlNxZNvwdcUKbfL4Ff5lmbmZlVzsOImJlZMoeGmZklc2iYmVkyh4aZmSVzaJiZWTKHhpmZJXNomJlZMoeGmZkl8zf3mdmeeewqyzg0ejr/MZtZF/LhKTMzS+Y9Ddtry5qXVdR+Sk0+dZhZ1/GehpmZJXNomJlZMoeGmZklc2iYmVkyh4aZmSXz1VNmli/fS/SR4j0NMzNL5tAwM7NkDg0zM0vmcxrVVunxXvAxXzOrGu9pmJlZMu9pmFn35T3xbsd7GmZmlizXPQ1J04DvA72BH0XELSXr+wE/BSYAW4ALI6I5W/dPwBXAduDrEbE0z1r3Rx6l1swqlVtoSOoNLABOA1qAFZLqI+LZomZXAG9GxGckzQC+B1woaQwwAxgLfAr4naSjImJ7XvWa2UeQbyzsdHnuaUwEmiJiI4CkxcB0oDg0pgN12fRDwP+QpGz54oh4H/h3SU3Z4z2ZY717z29Ms48e/12XlWdoHAFsKppvAT7fXpuIaJW0DTgsW768pO8R+ZXac/kQk1k39BE+gZ9naKjMskhsk9IXSXOAOdns25LWV1RhmkHAG537kN/qjP57WVe5TdvpfduprUueu6P++9k225f32c6+e7nNOuW5O9JBXbk/9576VmGbJWmvriMreZA8Q6MFGFo0PwTY3E6bFkkHAB8Htib2JSIWAgs7sebdSGqIiNo8n2NvdNe6oPvW1l3rgu5bm+uqXHetrbPqyvOS2xXAKEkjJPWlcGK7vqRNPTAzmz4feDwiIls+Q1I/SSOAUcBTOdZqZmYJctvTyM5RzAWWUrjkdlFErJM0D2iIiHrgx8C92YnurRSChazdzymcNG8FrvGVU2Zm1ZfrfRoRsQRYUrLsxqLp94AL2un7beDbedaXKNfDX/ugu9YF3be27loXdN/aXFflumttnVKXCkeDzMzM9szDiJiZWTKHRkbSNEnrJTVJur7M+n6SHszW/1nS8C6oaaik30t6TtI6Sd8o02aKpG2SVmc/N5Z7rJzqa5a0JnvehjLrJWl+ts0aJR3XBTUdXbQtVkt6S9Lfl7Tpsm0maZGk1yStLVp2qKRHJW3Ifg9sp+/MrM0GSTPLtenkum6T9Hz2Wv1K0iHt9O3wdc+hrjpJLxe9Xme207fDv+GcanuwqK5mSavb6ZvnNiv7OZHb+ywi9vsfCifqXwBGAn2BZ4AxJW3+M3BXNj0DeLAL6vokcFw2fRDwb2XqmgL87yptt2ZgUAfrzwQeoXCzwSTgz1V4Xf8vcGS1thkwGTgOWFu07Fbg+mz6euB7ZfodCmzMfg/MpgfmXNfpwAHZ9PfK1ZXyuudQVx3wzYTXusO/4TxqK1n/z8CNVdhmZT8n8nqfeU+jYOeQJxHxAdA25Emx6cA92fRDwFRJ+3rHWYci4pWIWJVN/z/gOXrWnfHTgZ9GwXLgEEmf7MLnnwq8EBEvduFz7iIi/kjhysBixe+le4DzynT9EvBoRGyNiDeBR4FpedYVEf8aEa3Z7HIK90d1qXa2V4qUv+Hcass+C/4T8EBnPmeKDj4ncnmfOTQKyg15UvrhvMuQJ0DbkCddIjscNh74c5nV/1HSM5IekTS2q2qicJf+v0paqcLd+aVStmueZtD+H3G1thnA4RHxChT+4IFPlGlT7W03m8JeYjl7et3zMDc7bLaoncMs1d5eXwRejYgN7azvkm1W8jmRy/vMoVGwL0Oe5E7SgcAvgb+PiLdKVq+icPjlc8B/B/6lK2rKnBARxwFnANdImlyyvprbrC9wLvCLMquruc1SVXPb3UDh/qj72mmyp9e9s/0A+DRwLPAKhcNApaq2vTIX0fFeRu7bbA+fE+12K7Osw+3m0CioZMgTtOuQJ7mS1IfCG+G+iPhfpesj4q2IeDubXgL0kTQo77qy59uc/X4N+BWFQwTFkoaDyckZwKqIeLV0RTW3WebVtsN02e/XyrSpyrbLToSeDVwc2UHvUgmve6eKiFcjYntE7AB+2M7zVe29ln0e/B3wYHtt8t5m7XxO5PI+c2gU7MuQJ7nJjpP+GHguIm5vp83ftJ1bkTSRwmu6Jc+6suf6mKSD2qYpnERdW9KsHrhMBZOAbW27y12g3f/5VWubFSl+L80EHi7TZilwuqSB2eGY07NluVHhS9P+ETg3It5pp03K697ZdRWfB/vbdp4v5W84L6cCz0dES7mVeW+zDj4n8nmf5XE2vyf+ULjS598oXIFxQ7ZsHoU/IID+FA51NFEYB2tkF9R0IoVdxUZgdfZzJnAVcFXWZi6wjsLVIsuBL3TR9hqZPecz2fO3bbPi2kThi7heANYAtV1U2wAKIfDxomVV2WYUgusV4EMK/6u7gsK5sMeADdnvQ7O2tRS+4bKt7+zs/dYEXN4FdTVROL7d9l5ru1rwU8CSjl73nOu6N3v/NFL4IPxkaV3Z/G5/w3nXli3/n23vraK2XbnN2vucyOV95jvCzcwsmQ9PmZlZMoeGmZklc2iYmVkyh4aZmSVzaJiZWTKHhpmZJXNomJlZMoeGWU4kHZ8Nstc/uyt4naRjql2X2b7wzX1mOZL03yiMJvAfgJaI+G6VSzLbJw4Nsxxl4yCtAN6jMFzJ9iqXZLZPfHjKLF+HAgdS+Ea1/lWuxWyfeU/DLEeS6il8i9wICgPtza1ySWb75IBqF2D2USXpMqA1Iu6X1Bv4k6RTIuLxatdmtre8p2FmZsl8TsPMzJI5NMzMLJlDw8zMkjk0zMwsmUPDzMySOTTMzCyZQ8PMzJI5NMzMLNn/B127y8c9V/FDAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "N, p = 10, 0.8; ns = np.arange(20)\n",
    "plt.bar(ns, stats.binom.pmf(ns, N, p) , color='green', alpha=0.5, label='binomial')\n",
    "plt.bar(ns, stats.poisson.pmf(ns, N*p), color='orange', alpha=0.5, label='poisson')\n",
    "plt.xlabel('x'); plt.ylabel('g(x)'); plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Exercise:** Check that from a \"large\" $\\nu$ the poisson distribution is equivalent to a gaussian distribution.\n",
    "\n",
    "**Exercise:** Generate $m$ datasets, each one with $n$ data, $x_i$ distributed random in $[0, 1]$ interval, sum the $n$ numbers, $\\sum_i^n x_i$ and obtain the distribution of the $m$ samples.\n",
    "\n",
    "**Exercise:** Generate $m$ datasets, each with $n$ data, $x_i$, normal distributed, compute its distance squared $\\chi^2 = \\sum_{i=1}^n x^2_i$, what is the distribution of $\\chi^2$?\n",
    "\n",
    "**Exercise:** Generate $n$ $x_i$-values, each one gaussian distributed with mean $\\mu_i$ and sigma $\\sigma_i$, show that the sum $\\sum_i x_i$ is gaussian distributed with mean, $\\mu = \\sum_i \\mu_i$, and sigma $\\sigma^2 = \\sum_i \\sigma^2_i$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "\n",
    "**Exercise:** Get the poisson distribution from a Gamma distribution.\n",
    "\n",
    "**Exercise:** Get an exponential distribution from a Gamma distribution.\n",
    "\n",
    "**Exercise:** Get the $\\chi^2$ distribution from a Gamma distribution.\n",
    "\n",
    "**Exercise:** Study the Beta distribution for different parameters of $\\alpha, \\beta$ (start $\\alpha=\\beta=1$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "***\n",
    "\n",
    "## Bibliography\n",
    "\n",
    "[1] \"Practical Statistic for LHC physicist,\" H. B. Prosper, CERN Academic Training Lectures (2015). https://indico.cern.ch/event/358542/ https://arxiv.org/pdf/1504.00945.pdf\n",
    "\n",
    "[2] \"Statistic for HEP,\" G. Cowan. CERN Academic Training Lectures (2012). http://indico.cern.ch/event/173726/\n",
    "\n",
    "[3] \"Statistics for Particle Physics,\" K. Cranmer, CERN Academic Training Lectures (2009). \n",
    "https://indico.cern.ch/event/48425/\n",
    "\n",
    "[4] \"Unified approach to the classical statistical analysis of small signals, \"G. J. Feldman and R. D. Cousins, Phys. Rev. D57 (1998) 3873. http://journals.aps.org/prd/abstract/10.1103/PhysRevD.57.3873\n",
    "\n",
    "[5] “Asymptotic formulae for likelihood-based tests of new physics,” Glen Cowan, Kyle Cranmer, Eilam Gross, Ofer Vitells. Eur. Phys. J. C71 1554 (2011). https://arxiv.org/abs/1007.1727\n",
    "\n",
    "[6] \"Incorporating systematic uncertainties into an upper limit,\" R.D. Cousins and V.L. Highland. Nucl. Instrum. Meth. A320, 331 (1992). http://www.sciencedirect.com/science/article/pii/0168900292907945\n",
    "\n",
    "[7] \"Confidence Level Computation for Combining Searches with Small Statistics,\" T. Junk, Nucl. Instrum. Meth. A434, 435 (1999). https://arxiv.org/abs/hep-ex/9902006  \n",
    "\n",
    "[8] \"How good are your fits? Unbinned multivariate goodness-of-fit tests in high energy physics,\" M. Willians, https://arxiv.org/abs/1006.3019\n",
    "\n",
    "[9] ROOT https://root.cern.ch, TMVA http://tmva.sourceforge.net, RooFit https://root.cern.ch/roofit\n",
    "\n",
    "[10] Anaconda https://anaconda.org, SciPy https://www.scipy.org, NumPy http://www.numpy.org, Scikit-learn http://scikit-learn.org/stable/, Matplotlib http://matplotlib.org\n",
    "\n",
    "[11] \"Lectures on Statistics in Theory: Prelude to Statistics in Practice\" B. Cousins, https://arxiv.org/abs/1807.05996\n",
    "\n",
    "[12] \"Statistics for physics\". D. Tonelli. Invisible School 2019, Canfranc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}