General.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Basic statistical concepts\n",
    "\n",
    "\n",
    "**Course: \"Introduction to data analysis methods for HEP\"**\n",
    "\n",
    "*Author: Jose A. Hernando*, January 2017\n",
    "\n",
    "*Particle Physics Deparment. Universidade de Santiago de Compostela, Spain.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## General information about the course\n",
    "\n",
    "These is a series of lectures on Data Analysis for High Energy Physics. Emphasis is put on 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",
    "The most common data analysis packages in HEP are ROOT, RooFit and TMVA [9]. But we will use the Python scientific toolkits, Matplotlib, Numpy, Scipy, Scikit-learn [10] that are distributed with Anaconda Python [10]. \n",
    "\n",
    "Structure of the course:\n",
    "* General Concepts: probability density functions and likelihoods. Frequentist and bayesian statistics\n",
    "* Hypothesis testings: test statistic. Continuous hypothesis. CLs method.\n",
    "* Confidence Interval: Classical and F.C. intervals.\n",
    "* MVA: Classification.\n",
    "* Regression: Likelihood estimation.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Chapter : Probability, Bayes, probability density funcions and likelihoods.\n",
    "\n",
    "Structure of the chapter:\n",
    "1. Introduction\n",
    "2. Probability and Bayes' theorem\n",
    "3. Frequentist and Bayes: the great divide\n",
    "4. probability density functions: moments, likelihood and posterior probabity\n",
    "5. common pdfs\n",
    "6. non-parametric pdfs\n",
    "7. Bibliography"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# general imports\n",
    "%matplotlib inline\n",
    "\n",
    "# numpy and matplotlib\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as stats\n",
    "\n",
    "# possible styles: ggplot (simplicity), bmh (scientify data), \n",
    "matplotlib.style.use('ggplot')\n",
    "\n",
    "# ploting images\n",
    "from IPython.display import display, Image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Last Execution  Fri Feb  1 17:50:59 2019\n"
     ]
    }
   ],
   "source": [
    "import time\n",
    "\n",
    "print(' Last Execution ', time.asctime() )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 1. Introduction\n",
    "\n",
    "This notes present some **statististical data analysis methods** common to High Energy Physics (HEP). \n",
    "\n",
    "These methods help us to **extract information from complex data**. \n",
    "\n",
    "In addition, Monte Carlo (MC) techniques allow us to **simulate experiments**.\n",
    "\n",
    "But we must understand **how and when to use a method, and what its result means**.\n",
    "\n",
    "In HEP we collect a **dataset** of events (**samples**), each one has a list of variables (**features**). \n",
    "\n",
    "\n",
    "\n",
    "For example a *dataset* is a collection of interactions of the same type and under the same conditions in our detector, and the variables that summarize each interaction, (i.e number of particles, total energy ,etc) are the *features*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Different tasks of data analysis:\n",
    "\n",
    "* **Classification** of the data. Data belongs to different categories (**classes**). For example a particle can be identified as a pion or a kaon. What is the most likely class for a given sample?\n",
    "* **Hypothesis testing**. In some cases we have competing hypothesis. Usually we have a default (or null) hypothesis and the alternative one. When do we reject the null hypothesis and accept the alternative one? When to discard the alternative hypothesis at a given confidence level? For example: is there a new particle, is the Higgs?\n",
    "* Setting **limits or confidence intervals** for a relevant parameter at a given **confidence level**. If a certain particle has not been found we set a upper limit on a cross-section, otherwise we provide a confidence interval that covers with a certain confidence level the range of the cross-section. \n",
    "* **Regression** (fits): estimate the value of a parameter (or parameters) and their uncertainties, i.e, to estimate the Higgs mass.\n",
    "\n",
    "We will briefly review some of the methods used in each of the tasks. But before, we should revisit basic concepts from probability and statistics."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "**Example:**\n",
    "\n",
    "Consider the case of a dice. There are six possible events $\\Omega = \\{1,2,3,4,5,6\\}$, each one has a probability $P(E)=1/6 \\; E=\\{1,2,3,4,5,6\\}$. \n",
    "\n",
    "If you are a frequentist you could compute the probability for each event just rolling the dice $n$, large, number of times, and computing the frequency of that event.\n",
    "If you are a bayesian, most likely, you will asing a prior probability of 1/6 for each event."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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) = \\frac{P(B|A)\\,P(A)}{P(B)}\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.\n",
    "\n"
   ]
  },
  {
   "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": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nsamples 100000\n",
      "p(A) = 0.138\n",
      "p(B) = 0.304\n",
      "p(B|A) = 0.398\n",
      "p(A|B) =  0.18\n",
      "p(B|A) p(A) = 0.055\n",
      "p(A|B) p(B) = 0.055\n"
     ]
    }
   ],
   "source": [
    "nsamples = 100000\n",
    "xs = np.random.choice(np.arange(6), size=(nsamples, 2))\n",
    "\n",
    "print('nsamples {:d}'.format(nsamples))\n",
    "As = np.array([xi for xi in xs if np.sum(xi) == 6])\n",
    "Bs = np.array([xi for xi in xs if len(np.argwhere(xi == 4)) > 0])\n",
    "pA = 1.*len(As)/(1.*nsamples)\n",
    "pB = 1.*len(Bs)/(1.*nsamples)\n",
    "print('p(A) = {:5.3f}'.format(pA))\n",
    "print('p(B) = {:5.3f}'.format(pB))\n",
    "\n",
    "BAs = np.array([xi for xi in As if len(np.argwhere(xi == 4) > 0)])\n",
    "ABs = np.array([xi for xi in Bs if np.sum(xi) == 6])\n",
    "\n",
    "pBA = 1.*len(BAs)/(1.*len(As))\n",
    "pAB = 1.*len(ABs)/(1.*len(Bs))\n",
    "print('p(B|A) = {:5.3f}'.format(pBA))\n",
    "print('p(A|B) = {:5.2f}'.format(pAB))\n",
    "\n",
    "print('p(B|A) p(A) = {:5.3f}'.format(pBA*pA))\n",
    "print('p(A|B) p(B) = {:5.3f}'.format(pAB*pB))"
   ]
  },
  {
   "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 +?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We want to estimate $p(HIV+ | +)$!\n",
    "\n",
    "$$\n",
    "\\begin{array}{l}\n",
    "p(HIV+) = 0.0001, \\; p(HIV-) = 0.9999, \\; \\\\ p(+|HIV+) = 0.99, \\, p(+|HIV-)=0.02 \\\\\n",
    "p(HIV+|+) \\, p(+)= \\, p(+| HIV+) \\, p(HIV+) \\\\\n",
    "p(+) = p(+|HIV+)\\, p(HIV+) \\, + \\, p(+|HIV-) \\, p(HIV-) \\\\ \\\\\n",
    "p(HIV+|+) = \\frac{p(+|HIV+) \\, p(HIV+)}{p(+)} \\\\\n",
    "p(HIV+|+) = \\frac{p(+|HIV+) \\, p(HIV+)}{p(+|HIV+)\\, p(HIV+) \\, + \\, p(+|HIV-) \\, p(HIV-)} \\\\\n",
    "p(HIV+|+) = \\frac{0.99 \\, 0.0001}{0.99 \\, 0.0001 + 0.02 \\, 0.999} \\\\\n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "probability to test + is p(+) = 0.000\n",
      "probability to be HIV+ is tested + is p(HIV+|+) = 0.331\n"
     ]
    }
   ],
   "source": [
    "pH, pHnot   = 0.0001, 0.9999\n",
    "ppH, ppHnot = 0.99  , 0.0002\n",
    "pp          = ppH*pH + ppHnot*pHnot\n",
    "pHp         = ppH*pH / pp\n",
    "\n",
    "print('probability to test + is p(+) = {:5.3f}'.format(pp))\n",
    "print('probability to be HIV+ is tested + is p(HIV+|+) = {:5.3f}'.format(pHp))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Exercise:** What is the probability if repeating the test, it still gives positive?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "probability to test + is p(+) = 0.328\n",
      "probability to be HIV+ is tested + is p(HIV+|+) = 1.000\n"
     ]
    }
   ],
   "source": [
    "pH    = pHp\n",
    "pHnot = 1.-pH\n",
    "\n",
    "pp  = ppH*pH + ppHnot*pHnot\n",
    "pHp = ppH*pH /pp\n",
    "\n",
    "print('probability to test + is p(+) = {:5.3f}'.format(pp))\n",
    "print('probability to be HIV+ is tested + is p(HIV+|+) = {:5.3f}'.format(pHp))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**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": "fragment"
    }
   },
   "source": [
    "$p(!R| !D)$?\n",
    "\n",
    "$$\n",
    "\\begin{array}{l} \n",
    "p(!D |R) = 0.05, \\; p(R) = 0.75 \\\\\n",
    "p(!R|!D) \\, p(!D) = p(!D|!R)\\, P(!R) \\\\\n",
    "p(!R|!D) = p(!R)/p(!D) \\\\\n",
    "P(!D) = p(!R) \\, + \\, P(!D|R) \\,p(R) \\\\\n",
    "p(!R|!D) = \\frac{p(!R)}{p(!R) \\, + \\, p(!D|R) \\, p(R)}\n",
    "\\end{array}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "probability that the reactor is off p(!R|!D) = 0.984\n"
     ]
    }
   ],
   "source": [
    "pnotR  = 0.75\n",
    "pR     = 1.-pnotR\n",
    "pnotDR = 0.05\n",
    "\n",
    "pnotRnotD = pnotR / (pnotR + pnotDR*pR)\n",
    "\n",
    "print('probability that the reactor is off p(!R|!D) = {:5.3f}'.format(pnotRnotD))"
   ]
  },
  {
   "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. They measure a **likelihood function**:\n",
    "\n",
    "$$\n",
    "f(x | \\mu)\n",
    "$$\n",
    "where $\\mu$ is the theory and $x$ the data. \n",
    "\n",
    "**Bayesians are deductive**. *They compute the probability using the Bayes' theorem. They need to assign first a prior probability based in a 'reasonable' initial guess*. They measure the credibility of a theory based on the data. They measure a **posterior probability**\n",
    "\n",
    "$$\n",
    "p(\\mu | x)\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "\n",
    "<img src=\"imgs/bayes_cartoon.png\">\n",
    "\n",
    "Most of physicist we are bayesians, but we think of ourselves as frequentists. 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, $f(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",
    "From the probability axions we have: \n",
    "\n",
    "$$\n",
    "\\int_{-\\infty}^{\\infty} f(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} f(x) \\mathrm{d}x\n",
    "$$\n",
    "\n",
    "In the case the input variables ${\\bf x}$ is a rv vector, the $f({\\bf x})$ is a scalar n-dimensional function. If the variables are uncorrelated them, the pdf is the product of the individual pdfs.\n",
    "\n",
    "$$\n",
    "f({\\bf x}) = \\Pi_{i=1,n} f_i(x_i)\n",
    "$$\n",
    "\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",
    "f(x_1,\\dots,x_{n-1}) = \\int f({\\bf x}) \\mathrm{d}x_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",
    "f(x_j) = \\int f({\\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",
    "f({\\bf n}|s,b) = \\Pi_{i=1,n} \\frac{(s_i+b_i)^{n_i}}{n_i!} e^{-(s_i+b_i)}, \\;\\;\\;\n",
    "f({\\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": [
    "### Likelihood \n",
    "\n",
    "The **probability density function**, *pdf*, is the probability to measure $x$ given some parameters $\\mu$ of the distribution, $f(x|\\mu)$. \n",
    "\n",
    "A **likelihood** is the probability evaluted with the observed data, $f(x|\\mu)$. Notice that now $x$ is a data, a measurement! \n",
    "\n",
    "For n, ${\\bf x}$, independent measurements of a distribution with a pdf, the likelihood is the product of the likelihood of each measurement $x_i$:\n",
    "\n",
    "$$\n",
    "f({\\bf x} | \\mu) = \\Pi_{i=1}^n f(x_i | \\mu)\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "\n",
    "The *pdf* of gaussian distribution with mean 0 and sigma 1, a normal distribution, is:\n",
    "\n",
    "$$\n",
    "f(x|\\mu=0, \\sigma=1) = \\frac{1}{\\sqrt{2\\pi}} e^{-\\frac{x^2}{2}}\n",
    "$$\n",
    "\n",
    "Consisder two measurements, $x_0=1., x_1=2.$, from a gaussian pdf expected distribution, with known sigma, equal to 1, and unknow mean, $\\mu$, the *likelihood* is:\n",
    "\n",
    "$$\n",
    "f(\\{1., 2.\\} | \\mu, \\, \\sigma = 1) = \\frac{1}{2\\pi} e^{-(1-\\mu)^2/2} e^{-(2 - \\mu)^2/2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "\n",
    "A **likelihood function** gives the likelihood a given measurement ${\\bf x}$ as a function of the distribution parameters $\\mu$.\n",
    "\n",
    "$$\n",
    "\\mathcal{L}(\\mu) \\equiv f(x | \\mu),\n",
    "$$\n",
    "\n",
    "\n",
    "One of the common methods to estimate a parameter is to find the best parameter value that maximizes the likelihood function (see the regression chapter). \n",
    "\n",
    "The likelihood can be a very small number, for that reason, it is common to take the logarithm, and called **log likelihood function**:\n",
    "\n",
    "$$\n",
    "-2 \\log \\mathcal{L}(x |\\mu)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Exercise:** 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": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbAAAAEPCAYAAAAj0pGKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmclNWZ8P3fqb27AQEbkG52WRQQkEYWiSyKG4qOGo/RmH3GSTLJJBPf98kkJsZxJnnMk48ZnXmSeYbX5HGMRucYN1QUFBSQfd8FBAQbBGyQxV5qPe8fd4FN03RXVVfXXVV9fT+f+th1171cx4K+OPd9znWUtRYhhBCi0HjcDkAIIYTIhCQwIYQQBUkSmBBCiIIkCUwIIURBkgQmhBCiIEkCE0IIUZAkgQkhhChIksCEEEIUJJ/bAQghhMgvWus/AjcDR4wxI5PbfgPMAiLAbuAbxpjj7kVZ+D0wKy95yUte8sro1ZIngRuabHsLGGmMGQXsBH7SyjnaXcH3wA4ePJjRceXl5dTU1GQ5GndIW/JPsbQDpC35qi1tqaioaPFzY8xirfWAJtvmN3q7AvhiRhfPokLvgQkhhMi9bwJvuB1EwffAysvLMzrO5/NlfGy+kbbkn2JpB0hb8lVb26K1XtPo7WxjzOwUj3sAiAHPZHzxLCn4BJZpF1puJeSnYmlLsbQDpC35qq23EI0x49I9Tmv9NZzBHdcYY1p7jtbuCj6BCSGEaH9a6xuAHwNTjTF1bscDksCEEEI0obV+FpgGlGutq4Ff4Iw6DAJvaa0BVhhjvu1akOQogTU3p6DJ5wp4HJgJ1AFfN8asy0VsQgghzmaMubuZzX/IeSCtyNUoxCc5d05BYzcCQ5Kv+4D/aPeIZCVqIYQoaDlJYMaYxcCxFna5FXjKGGONMSuArlrr3u0SjLUEli/H82//hqqtbZdLCCGEaH/5Mg+sEvio0fvq5LbsUwrPkSOokyfxb9zYLpcQQgjR/vJlEIdqZluz9/i01vfh3GbEGJPZPIhrrsHzzDN03bWLxMyZ4MuX/w2Zkbkt+Seb7fjnf/Y2u/3nP49n5fytKZbvBKQtxSZffnNXA30bve8DNFsjKjnZ7vSEO5vRPIiyMnr26EH93r3Uv/cesZHnjCspKDK3Jf9ksx11dZ2b3V5Tcyor529NsXwnIG05rbVSUoUiXxLYHOB7WuvngAnACWPMx+12NaWwEyfC3r0E1q4lNmIEqOY6gUIIIfJVrobRNzenwA9gjPk/wFycIfQf4Ayj/0Z7x2SHDydRVoa3pgbv/v3E+/dv70sKIYTIopwksPPMKWj8uQX+LhexnOHzER09muCyZQTWraNeEpgQQhSUfBmF6Iro6NFYrxffnj2oYy2N8hdCCJFvOnQCs6WlRIcPByCwTgp/CCFEIenQCQwgWlUFgH/rVlR9vcvRCCGESFWHT2CJCy8kNnAgKhaTic1CCFFAOnwCA4ic7oVt2ACxmMvRCCGESIUkMCDerx/xHj3w1Nbie/99t8MRQgiRAklgAEqd6YUF1qyRSvVCCFEAJIElxS65xJnYfPQo3g8/dDscIYQQrZAEdprXS3TsWCDZCxNCCJHXJIE1Ehk1Cuv349u/H8/hw26HI4QQogWSwBoLhYiOGgVAYPVql4MRQgjREklgTUTGjsV6PPh27UKdOOF2OEIIIc5DElgTtksXYsOGoRIJAmvXuh2OEEKI85AE1ozIuHEA+DdvBikvJYQQeUkSWDMSPXsS698fFYsRkPJSQgiRlySBnUfkiisA8K9fD9Goy9EIIYRoShLYecT79SPeqxeeujr8W7e6HY4QQogmJIGdj1JnemGBNWsgkXA5ICGEEI353A4gn8WGDCHRtSue48fx7dxJ7JJL3A5JCCHandb6j8DNwBFjzMjktu7AfwMDgA8BbYz51K0YQXpgLfN4Pu+FrVolRX6FEB3Fk8ANTbb9I7DAGDMEWJB87ypJYK2IDh/uFPn95BO8+/a5HY4QQrQ7Y8xi4FiTzbcC/5X8+b+Av8ppUM0o+FuI5eXlGR3n8/lSPlZNn45n4UJKtm8nkZwjlk/SaUu+K5a2ZLMdP/yht9nt5eXBrJy/NcXynYC0pTGtdeOq5bONMbNbOaSXMeZjAGPMx1rrnhlfPEsKPoHV1NRkdFx5eXnqx/bvT6dYDLVtG7WbNpGoqMjomu0lrbbkuWJpSzbb8dhjnZvdfv/9p7Jy/tYUy3cC0pbTKioqMMbk37/G0yS3EFMRChEZMwaA4KpVLgcjhBCuOKy17g2Q/O8Rl+ORBJaqaFUV1uvFt3s3niL5F5wQQqRhDvC15M9fA15xMRZAEljKbGkp0csuA5IjEoUQokhprZ8FlgPDtNbVWutvAY8A12qtdwHXJt+7quCfgeVSZNw4/Js24Xv/fdSVV2K7dnU7JCGEyDpjzN3n+eianAbSCumBpcFecAGxSy9FWSsLXgohhMskgaUpMn48Vin8W7agTuVmFJgQQohzSQJLU6J7d2JDhjgLXq5Z0/oBQggh2oUksAxEJkwAwL9pE6quzuVohBCiY5IEloFEz57ELr4YFYvhl16YEEK4ImejELXWNwCPA17gCWPMI00+74dTX6trcp9/NMbMzVV86QqPH49v924CGzY4BX9LStwOSQghOpSc9MC01l7gd8CNwHDgbq318Ca7/QwwxpjLgS8Bv89FbJlKVFQQ698fFY0SWL/e7XCEEKLDydUtxPHAB8aYPcaYCPAcTmXjxizQJfnzBcDBHMWWscjEiQAE1q2DhgaXoxFCiI4lVwmsEvio0fvq5LbGHgLu1VpXA3OB7+cmtMzF+/Qh1qcPKhwmsGGD2+EIIUSHkqtnYKqZbU1Xh7wbeNIY86jWehLwJ631SGNMovFOWuv7gPsAjDE5WU6lRTNn4n36aUref5/EtddCMDdLXDQmS0Tkn2y2o7RUllPJFmlLcclVAqsG+jZ634dzbxF+i+QKoMaY5VrrEFBOk4rHyTVrTq9bY3OynEpLOnWipFs3fAcOEF6wgMj48W0/Z5pkiYj8k8121NU1v5xKTY0sp5IuaYujIs+WhMpUrm4hrgaGaK0Haq0DOIM05jTZZz/JOlta60uBEPBJjuLLnFKfPwtbswYiEZcDEkKIjiEnCcwYEwO+B8wDtjubzFat9cNa61uSu90P/I3WeiPwLPB1Y0zT24x5Kd6/P/HevVH19QQ2bnQ7HCGE6BByNg8sOadrbpNtDzb6eRswOVfxZJVShCdNovTFF/GvWUNk9GgIBNyOSgghippU4siS+IABxHv3xlNXJ70wIYTIAUlg2ZLshQHOUivyLEwIIdqVJLAsig8YQLyiwnkWJvPChBCiXUkCy6amvbBw2OWAhBCieEkCy7J4//7EKitRDQ1SI1EIIdqRJLBsU4rI6V7Y2rVSI1EIIZSaiFIPodSbKLUJpXah1HKUehKlvoFS3TI5rSSwdhDv18+pkdjQ4BT6FUKIjkipr6HUZmAZ8EOgFNgFrAQ+BSYATwAHkslsYDqnlwTWHpQiMtmZ0hZYuxZVX+9yQEIIkWNKbQQewZn/WwV0w9opWHsH1t6LtTOx9lKgO/A3QE9gK0rdleolJIG1k3ifPs56YZEI/tWr3Q5HCCFy7f8CA7H2x1i7Hmubr6xk7QmsfQZrZwKTgOOpXkASWDsKn+6FrV+Pqq11ORohhMghax/D2vQGAVi7EWvnpbq7JLB2lOjdm9jFF6NiMQKrVrkdjhBCuEOpr6DUSJQ6d20gpTJeFyhntRA7qvCVV+LbvRv/xo1EqqqwXbq0fpAQQrhMa/0PwF/jrN24GfiGMSbTYdX/lTxPGKW2AuuTrw3AOJS6C2u/kO5JpQfWzhI9exIdOhQVjxNYscLtcIQQolVa60rg74FxxpiRgBdnGaxMdQdmAD8DduAUbv934D3gMWBEJieVBJYD4cmTsR4P/q1bUceOuR2OEEKkwgeUaK19OMPfmy5CnDprj2PtO1j72+QIxBFAb+BXQD1wdyanVecbGFIgbCTDork+n49YLJblcM5Pvfoqng0bSIwYgb399qyeO9dtaU/F0pZstuPw4XMfGwD06hXPyvlbUyzfCUhbTgsEAmit1zbaNDu52v0ZWusfAL/ESTDzjTFfzjjYlij1O6AT1n4t7UMLPYEdPJjZPwpyvbS4OnmSsj/8AZVIUHvvvSR69crauWWZ9PyTzXY8+mjnZrfff/+prJy/NcXynYC05bSKigoAdb7PtdbdgBeAu3CGtT8P/MUY83RGF2yJUrcCT2HtBekeKrcQc8R26UJ0zBgAgkuXuhyNEEK0aAaw1xjziTEmCrwIXJnx2ZT6GUrNRKnezXzaE8honpGMQsyhyIQJ+Ddvxrd3L97qauJ9+rgdkhBCNGc/MFFrXYpzC/EaYE0bzvc/gE6ARakjwDqcUYhxnJGOP87kpNIDyyFbWkpk3DgAgosXQ2HfvhVCFCljzErgLziJZjNOrpjd4kEtuwC4FLgXeAYoAb4L/BxnMMdvUWoeSv0KpVIeJCA9sByLVFXh37gR78cf49u9m9jgwW6HJIQQ5zDG/AL4RVZO5gy22JF8PXtmu1IXA2OTryqcmog/xhm23yrpgeVaMEhk4kQAAu+9B4mEywEJIYRLrN2Ntc9j7U+w9jqs7QGkXJFeemAuiI4aRWDtWrxHj+LfupXoZZe5HZIQQmSXUl9tw7Fg7VOt7SYJzA1eL+HJkymZO5fAsmVEL7kE/H63oxJCiGx6sg3HWqDVBCa3EF0Su+QS4j174vnsMwLr17sdjhBCZJe1nja85BlYXlOK8JQpAARWrkTV1bkckBBCFBZJYC6K9+9PbMAAVCQihX6FEMVLqbHtcVpJYC4LT5mCVQr/xo2oTz91OxwhhGgP76DU9GyfVBKYyxI9ehAbMQKVSBBcssTtcIQQoj38GZiLUnec84lSX0Cp9zI5qSSwPBC+8kqsz4d/1y48Bw64HY4QQmSXtd8B/ifwHEp9GwClLkOpV4HFQLdMTisJLA/Yzp2JVFUBEFq0SEpMCSGKj7UPA98G/g2lFuHUQhwJfBPIaDJsyvPAtNbTgQ+NMXu11r2BR3AKMf7UGHMok4uLz0XGj8e/ZYtTYmrHDmKXXOJ2SEIIkT1KdQeG4uSNq4BlwDSszXiBtnR6YL9PXhjgUcCPM9msLQUexWmBAJHJkwGcZ2FFsuieEEKg1C+APcDf4eSPbwLjgN+25bTpVOKoNMbsTy4vfT3QH4jQlmWmxVmiI0bgX7cOb00NgXXriIwf73ZIQgiRDQ8ATwD/hLWHAVBqP/ASSvUC7sXaaLonTSeBndRa98K5Z7nNGPOZ1jqA0xNrldb6BuBxnCrDTxhjHmlmHw08hNOz22iMuSeN+Aqfx0N42jRK//IXAitXEh0xAltW5nZUQgjRVpdi7e6ztli7MDm0fi7wJs6aY2lJ5xbivwOrcdZy+V1y22Tg/dYO1Fp7k8fcCAwH7tZaD2+yzxDgJ8BkY8wI4IdpxFY04v37Exs0yJncLCs3CyGa4d2zp7AGezVNXp9vXwd8ARiQyWlTTmDGmF/jLDM92RjzXHLzAZzVNFszHvjAGLPHGBMBngNubbLP3wC/M8Z8mrzekVRjKzbhqVOxHg/+LVvwHOmw/xuEEM3w7tlD6Usv4XnmmcJKYudj7QfAlZkcmu4w+r1Ahdb6ruT7AzgP5lpTCXzU6H11cltjQ4GhWuulWusVyVuOHVKie3eiY8agrCX4zjvF8YdUCNF28Tihd98FwA4e7Cw7kq+UegWlLk9pX2sPo1QIpX50Zp5YCtIZRn8ZMAcIA32A/wamAl8D7mrhUIDm/i83/a3sA4YA05LnX6K1HmmMOd4kjvuA+wCMMZSXl6fahLMv5vNlfGxO3Hwznn37UMeOUVZTA5deet5d874taSiWtmSzHaWlzRfmLi8PZuX8rSmW7wQKvy1q+XI84TC2b188kyZRnt//uN0PrECpDTiPnt4DNp01bF6pCpw7dLOA23E6Rd9M9QLpDOL4D+BBY8yftNani/YtAv6/FI6tBvo2et+Hc0cvVgMrjDFRYK/WegdOQlvdeCdjzGw+H7pva2pq0mjC58rLy8n02FzxjxlDaMECEi+/TG23buBr/usqhLakqljaks121NV1bnZ7Tc2prJy/NcXynUBht0XV1lL25puoSIS666+nm7UZt6WioiLL0TXD2u+j1GM44xkeAi4ALEqdxOkIdcMZBKiAVcn9/oS1KS9Tn04CGwE8fTo0AGNMrda6JIVjVwNDtNYDcTLsl4CmIwxfBu4GntRal+PcUkzl9mTRio4ahX/jRmdY/Zo1RCZOdDskIYRLAkuXoiIRYgMHEh80yO1wUuMM3vg+St0PTAImABVACDiKMwhwMdbuy+T06TwD+xCoarxBaz0e+KC1A40xMeB7wDxgu7PJbNVaP6y1viW52zzgqNZ6G/AO8P8aY46mEV/x8XgIT3cKOAdWrkSdPOlyQEIIN3gOH8a/ZQvW46Fh2jS3w0mftRGsXYS1/wtrf4i138baB7D2T5kmL0ivB/Zz4HWt9f8BAlrrn+DUtfqbVA42xszFGe/feNuDjX62wI+SL5EU79eP6NCh+HfuJLhoEQ2zZrkdkhAil6wltHAhyloiY8diu3d3O6LsUcp5IGltRvdC0xlG/xrOPK4eOM+++gO3G2PmZ3Jhkbrw1KlOtfqdO/Hu3+92OEKIHPJt3Yr34EESZWWEJ01yO5z0KKVQ6lqUmnaePY4BP0KpjArLp9MDwxizDvhuJhcSmbNduhCZMIHg0qUEFy6k7qtfBY8sJCBE0WtoIPies1RW+KqrIJibkadZ9Cjw94BCqVlYe9ZdOKxNoNQcnLt5v0/35C0mMK31w6mcpPGtQNE+IuPGOdXqjx7Fv3490aqq1g8SQhS04IoVeGpriVdUEBs+vPUD8s8g4GLgF8DO8+yzHvhPsp3AOHvo+/nk9USEouHz0TBtGqWvvEJw2TJiw4ZhO3VyOyohRDvx1NTgX78eqxQNV1+d80nLWuuuOAV4R+L8nv+mMWZ5mqdZCfiwtqW5XeOA4Sjlxdp4C/udo8UEZoz5RjonE+0rPngwsUGD8O3ZQ3DxYhpmznQ7JCFEe7CW4IIFqESCyKhRJHr1ciOKx4E3jTFfTBZuL83gHL8H5qLUT7F20Xn2+Q6QSL7SktYzsGTBXY0zjv8gznD4XeleVGSuYfp0yvbtw799O9FRo4j36eN2SEKILPNt346vuhpbUkL4C1/I+fW11l2AKcDXAZI1bCNpn8jaEyj1MPAaSn0GLAUWA29h7fbkXn2A/dj0y4qkPBJAa30Pzr3KUUAtzhLQ65LbRY7Yrl3PrBMWfPttiKfV4xZC5LuGBoKLnM5Kw5QpUJJKrYisGwR8AvxfrfV6rfUTWuvM1naydh7OHOL3gJuBx4AtKLUXpX4KHAF2ZHLqdHpg/wLMNMYsPr1Ba30V8Cfgz5lcPBuKthZiS2bOdOokHj9O2Z49eCsrC7ctTRT099JINtvxwx9KLcRsKYS2qDfewAPYYcMITp163mdfbW2L1npNo7ezk2X6zpweGAt83xizUmv9OPCPOPOB02ftTuBOlOqEs3zKlOTrweS1XsvktOkksM5A0wd4KwBXV1ws5lqILfFOmEDpSy9h33iD0PDh1ETTXsw0LxX693JaNtvx2GPN10K8/36phZiufG+L5/BhSpcsAaWomzCBxNHzFyNqS1sqKiowxoxrYZdqoNoYszL5/i84CaxtrP0MZ/HKNwFQKgRcC/xtJqdLZzLRb4Ffaa1DAMkaiL9Mbhc5Fh80iOiQIahoFM+8eW6HI4Roq0SC0FtvoawlevnlJHr0cC0UY8wh4COt9bDkpmuAbVm/kLUNWPsqGSbHdHpg3wUuAn6QrEbfDaeK8Mda6++c3skY0y+TQET6wtOn49u3D7VjB75+/YgNHux2SEKIDPk3bsR7+DCJTp3ypeLG94FnkiMQ9wDtNyrd2i2ZHJZOArs3kwuI9mM7dyY8eTJlK1cSXLiQWL9+EAi4HZYQIk3q1KnPK25cc01eVNwwxmzAmaPVdkotbOHTBHACWAv8AWsPp3ralBOYMeZ8Y/iFi6JjxmD37cOzZw/BZcsIF2KlaiE6uOA776AiEaKDBxfrnRSFs0RWb2AvcBjoBQwEPk6+nwn8A0pNxdqUblemsyKzD2e9rsuBs0pAGGPuS/U8Iss8HhI33YT93/8b/7p1RC+91K1Jj0KIDHh378a/axfW7z+zfFIR+i3O8PkqrF1/ZqtSVYAB/gmnBzYfZ2zFbamcNJ1BHE/jPGhL4GTLxi/hpooKomPHoqwlNH8+JNKe0C6EcEM4TGjBAufHK6/EdunickDt5l+Ah85KXgDWrsVJXv+CtdXAb3CG16cknWdgNwB9jTG5Gbsr0hK+8kp8O3fiPXIE/9q1RK+4wu2QhBCtCC5diufUKeK9ehEdO9btcNrTUOB8Y/4/AU7fN91NGlOz0umBbQOKaCW1IhMI0DBjBgDBZctQx4+7HJAQoiWegwfxb9jgrLJ83XXFvkTSh8Bfn+ez+5KfA5QD55/81kS6oxCf0FrPp8ltQ2PMU2mcR7ST+KBBRIcNw79jB6G336b+jjtyXsFaCJGCeJzQ/PkoawmPG0eiZ0+3I2pvDwNPo9Qm4AWc8lE9gTtwqt2fLkk4A6eCfUrSSWBfB67Cmf9V32i7BSSB5YnTc8N8+/bh27qV2MiRbockhGgisGoV3qNHSXTtSiQ/5ny1L2ufRakanOddPwX8QBRYA1yHtW8n9/wRkHKB13QS2A+Ay40x21vdU7jGlpXRMG0aJW++Sejdd6kdMEDWDRMij3hqagisWAFAw7XXgt/vckQ5Yu1bwFso5cG5VViDtYkm+zSkc8p0broeBvanc3Lhjtjw4cQGDkSFwwSTI5yEEHkgkSA0b96Zdb7i/Tpg4SJrE1h75JzklYF0emD/ilNW5BGc+5dnGGP2tDUQkUVK0TBjBmVPPon/gw+I7dxJbOhQt6MSosPzr12L99Ahp1zUlJRHixcHpXoD9wNTcQYEHgXeBX6LtYcyOWU6PbDfAbcAy4APGr1kQcs8ZLt0OfMXJLhgAaquzuWIhOjY1LFjBJcuBZK3DvOgXFTOKDUU2AD8PfAZsApnXckfABtQakgmp02nlFRRj/EsRtHRo/Ht3Invo48ILlhAw6xZbockRMeUSDijDuNxoiNGEB80yO2Icu3XwElgAtZ+eGarUv1xqm/8Grg93ZNKUipmStFw3XVYvx//zp34dmS06KkQoo3869bhO3CARFkZDVOnuh2OG6YDPz8reQFYuw94KPl52tKthfhdnPuX5TjFGQEwxnSwm7mFw3btSnjKFEILFhBasIDavn2xpaVuhyVEh+E5duxMpfmGa6+FkhKXI3JFADhfFadTyc/Tlk4P7F9xVs1cDFThTEbrCbRUJl/kgejo0cT69kXV1xN8+22w1u2QhOgYEglCb7zh3DocOZL4xRe7HZFbNgDfTw6h/5xSCqdjtCGTk6aTwG4HbjTGPA7Ekv/9KzLs+okcUoqG66/HBgL4d+3Ct12m8gmRC4HVq8+MOuygtw5PexinysZ2lHoYpb6DUv8EbAWuxZngnLZ0Elgp8FHy53qtdakx5n2c5VVEnrMXXEBDcqmG0MKFqFNSk1mI9uQ5coTAsmUANFx/PYRCLkfkImvfBG7GuV34AM6o9p/hjEi8GWvnZ3LadBLYduB0ifM1wENa658BBzK5sMi92IgRxC6+GBUOE3rzTbmVKER7icUIzZ3rTFgePZr4gAFuR+Q+a9/E2nFAZ6Av0Blrx2PtvExPmU4C+wFO7Spw6lWNxcmosphloVCKhmuvxZaU4Nu/H/+GjG47CyFaEVy61Kl12K1bx5uw3Bpr67D2ANa2eXJqOgmsE3B6jY7PgIM4S6zIROYCYsvKPl92ZdEiPEdTXrlACJEC7/79+NeuxXo81N9wAwQyGmBX+JRKoFQ8xVcsk0ukU0rq98D1yZ8fTf43CszGqdDRIq31DcDjgBd4whjzyHn2+yLwPHCFMWZNGvGJFMWGDiU6YgT+rVsJzZ1L3T33gNfrdlhCFL6GBkJvvukskzJxIomKCrcjctPDOKuVtJt0ElilMWZ/cj7Y9UB/IILTE2uR1tqL89DuWqAaWK21nmOM2dZkv844pUZSXg9GZKZh+nS81dV4jxwhuHSp3OYQoq2sJfT2284KyxddRGTCBLcjcpe1D7X3JdK5hXhSa90LZyLzNmPMZ8ntqawFMB74wBizxxgTAZ4Dbm1mv38G/heQVkl9kYFgkPobb8QqhX/NGrz7ZaEBIdrCt3Ur/h07sH4/9TNnyl2NHEgngf07sBp4Bqc3BTAZeD+FYyv5fAg+OL2wysY7aK0vB/oaY15LIybRBonKSiITJqCsdSZb1te3fpAQ4hzq008JLXRqOjRcfTW2WzeXI+oY0inm+2ut9UtA3BizO7n5APDXKRze3Lr2Z+6Naq09OJU+vt7aibTW95Ec+WiMoby8PIXLn8vn82V8bL5pU1tmzcJz7Biqupqy5ctJ3HknqOa+rtwolu8lm+0oLW3+X/Ll5bmpZl4s3wm0U1vicTxz5qD8fhKjRxOcOjUnf4eK6XvJVDrPwDDG7GzpfQuqccb9n9aHs5+ddQZGAu9qrQEuAuZorW9pOpDDGDMbZ+AIgK2pqUm9AY2Ul5eT6bH5pq1tUVOmUPbUU6j162no3p3omDFZjC49xfK9ZLMddXWdm91eU5ObyejF8p1A+7QlsHgxwV27SHTuTO348ZCjkb1taUtFCoNLkmMX1gAHjDE3Z3ShdpZWAmuD1cAQrfVAnF7bl4B7Tn9ojDmBUyAYAK31u8D/I6MQc8NecAEN111HyWuvEXz3XeKVlSR69HA7LCHynvfDDwmuXo1Vioabbiq2Qr0/wClg0aVNZ1FqGs4jo+1Yu66ZzyuBb2Htw+meOifLqRhjYsD3gHk4/0OMMWar1vphrXWrQ/BF+4sNG0bksstQ8Tih11+HSMTtkITIa6q2ltAbbwAQufJK4pWVrRxROLTWfYCbgCcyPolSnVBqGbAA+BOwGqXeRKmm3b8+wC8ckgQqAAAYcElEQVQyuUSuemAYY+YCc5tse/A8+07LRUzibOFp0/AePIj36FFCCxfScMMNbockRH5KJAjNnYunro5Y375Exo93O6Jsewz4HziPdzL1U+BSnLENq4FpOEV7V6LU9Vi77fyHpkbZwq6HZyMZ9hR8Ph+xWEaTv/NOVtty5AieP/4RFY2SmDULm+PnYcXyvWSzHYcPNz+Io1eveFbO35pi+U4ge21RixfjWbQIW1JC4m//Fjq35fd8ZtrSlkAggNZ6baNNs5PjC9Ba3wzMNMZ8V2s9DedxTvrPwJR6H/g91v5bo22VwCvAAOBGrF2NUhOAZVib9ryDgk9gBw+2Oo+6WfJg+vx8W7ZQMm8e1uej7stfJpHDkU7F8r1ksx2PPtr8L8f775dBHOnKRlu8+/ZR8sILANTffrtrhXqzMIij2aGSWuv/CXwFiAEhnGdgLxpj7k3rIkrVAddj7ZIm28uAV4FxOFWc6skwgeXkGZgoLLERI4gOH46KxQi9+qo8DxMiSX32mVNl3loiEyYUZZV5Y8xPjDF9jDEDcAbcLUw7eTmO4DzfOpu1tcCNwBKcx0o3ZRqrJDBxLqVouOYa4hdeiPfYMULz58vSK0IkEoRef/3z516TJrkdUb5bQ/MVl8DacPKz13DWBcuIJDDRvECAhltuwfr9+HfswL9+vdsRCeGq4JIl+KqrSZSVOUPmPcX/69MY824b5oA9C/RHqQub/dTaGHAX8J9ARrXsiv8bEBlLdO/urCRLcumVA7J2qeiYfDt3ElizBuvx0HDzzdiyMrdDyn/WvoC1k7D2/DO7rbVY+x2sHZjJJSSBiRbFhg0jUlWFSiQoee01VG2t2yEJkVOeo0edFcyB8JQpxPuc+1hHuEMSmGhV+KqriFVW4vnsM0KvvQbx3AzfFsJ14TChOXNQ0SjRYcOIjh3rdkSFT6mFKJWVfwVIAhOt83ppmDWLRFkZvupqgu++63ZEQrQ/ayl54w28x44Rv/BCGq67ztVC10VkGlCajRNJAhMpsWVl1N9yC9bjIbBhA74tW9wOSYh2FVi2DN/u3dhgkPpbb4VAwO2QRBOSwETKEhUVhGfMACD01lt4MpxELkS+8+3aRXDFCqxS1N90k6zvlackgYm0RC+7jMjo0c6gjjlzUKdyUw1CiFzx1NR8XqT3qquID8xogJzIAUlgIm3h6dOJ9e2Lp7aWkpdfhmjU7ZCEyApVV0fJSy85gzYuuYTIuHFuhyRaIAlMpO/0oI6uXfEeOeIMMZZKHaLQxWKE5szBc/Ik8YsukkEbBUASmMiILSmh/q/+ChsI4N+5k8CyZW6HJETmrCX01lv4Dhwg0amTM2jD73c7KtEKSWAiY4kLL3QecCtFcMUKfFu3uh2SEBkJrFqFf9s2rM/n/MOsUye3QxIpkAQm2iQ+aBDh6dMBCM2fj/ejj1yOSIj0+LZvJ/jee1ilaJg5k0SvXm6HJFIkCUy0WfTyy4mMHXtmZKLn2DG3QxIiJd4DB84qExUbMsTliDqEa8mweG9TksBEVoSnTiV28cWohgZKXnxRaiaKvKeOHaPk5ZdRiQSRMWOIVlW5HVLHYO0CrG3IxqkkgYns8HionzmT+EUX4TlxgpKXXpKFMEXeUp99RukLL6AaGoidvg0uIw4LjiQwkT2BAPW33eYMrz98mJJXX5XCvyL/hMOUvPjimeHy9R1kba9iJN+ayCpbWkrd7bdjS0rwffihrOYs8kssRsmcOXg/+YRE167U33ab1DgsYJLARNbZbt2ou+02ZzXnbdsILlokSUy4L5FAvfIKvv37SZSVUXfHHdjSrBRFFy6RBCbaRaJ378+r169dS2DVKrdDEh2ZtQTffhvPtm3Y5K1u27Wr21GJNpIEJtpNfMAAGmbOdCY6v/ce/o0b3Q5JdFCBJUsIbN7sTFS+7TaZ61UkJIGJdhUbNozwNdcAEFywQKp1iJwLrFpFcPVqrMdD4o47iPfJymLAIg9IAhPtLjp6NOEpU1DWEpo3D9/OnW6HJDoI/9q1BJcscapsXH89DB3qdkgiiySBiZyIXHEF4UmTnCT2+ut4d+92OyRR5PwbNxJ6910AwjNmEBs+3N2ARNZJAhM5E5k0ici4cWdKTnn37HE7JFGkfFu2EHr7bQAarr6a6KhRLkck2oMkMJE7ShGeMoXI5Zc7SeyVVySJiazzbdnizD/EqW8YvfxylyMS7cXndgCig1HqTPX6wPr1lMyZQ/0ttxAfNMjlwEQx8G/e/HnyuuoqIldc4XJEhUlr3Rd4CrgISACzjTGPuxvVuaQHJnIvmcQiY8ag4nFKXnkF3wcfuB2VKHD+TZvO6nlFxo93OaKCFgPuN8ZcCkwE/k5rnXcPESWBCXcoRfjqq88swxJ69VV8O3a4HZUoUP41awi99RYADVOmSM+rjYwxHxtj1iV/PgVsByrdjepcObuFqLW+AXgc8AJPGGMeafL5j4C/xsn8nwDfNMbsy1V8wgVKEZ42Dev1Ely9mtDrr6M6dYLKvPt7IvKVtQRWrCC4bBmQHLAhz7yySms9ALgcWOlyKOdQNgc16rTWXmAnzkJm1cBq4G5jzLZG+0wHVhpj6rTW3wGmGWPuauXUNpLhkh0+n49YLJbRsfmm4NtiLWrJEjyLFuHxeIjNmIGdMMHtqNokm9/J4cPeZrf36pWbSv95++fLWtSCBXiWL8cqhb35ZuyYMS0ekrdtyUBb2hIIBNBar220abYxZnbT/bTWnYBFwC+NMS9mFmn7yVUPbDzwgTFmD4DW+jngVuBMAjPGvNNo/xXAvamcuKamJqOAysvLMz423xRFW4YPx19fT/dVq2h4+WXCH39MZPLkgl2jKZvfyWOPdW52+/33n8rK+VuTl3++EglC8+fj37oV6/HQMHMmsT59oJU487ItGWpLWyoqKjDGjGtpH621H3gBeCYfkxfkLoFVAh81el8NtPRP7G8BbzT3gdb6PuA+AGMM5eXlGQXk8/kyPjbfFE1brr8e1acPJS+/TOnmzSR8PuzMmQW5VlM2v5PS0uZ7YOXlwaycvzV59+crGsXzwguovXuxF1xA4otfJDR4cEqH5l1b2qA926K1VsAfgO3GmN+2y0WyIFcJrLl/Rjd771JrfS8wDpja3OfJbu7prq6VHliRtWXECD49dYqS115DLV1K7NAh6mfOLLg1m7L5ndTVNd8Dq6npeD0wVV9Pycsv4z14EBsKUTdzJomuXVvteZ2WT21pq7b2wFoxGfgKsFlrvSG57afGmLkZXbCd5CqBVQN9G73vAxxsupPWegbwADDVGBPOUWwiz8QHD6b+zjspefllfLt3U/r8887yF7J2U4emPv2U0hdfxHP8OInOnam/4w4SF17odlhFyRjzHs13PPJKru7NrAaGaK0Haq0DwJeAOY130FpfDvwncIsx5kiO4hJ5Kl5ZSd3dd5Po0gXvoUOU/vnPeI4edTss4RLPgQOUPfssnuPHiffoQd0990jyErlJYMaYGPA9YB7OfAJjjNmqtX5Ya31LcrffAJ2A57XWG7TWc85zOtFBJLp3p+6ee4j36oXnxAlKn30W74cfuh2WyDHftm2UPv88qr6e2MCB1N11F7ZTJ7fDEnkgZ/PAkvdO5zbZ9mCjn2fkKhZROGxZGXV33UXojTfw79pFyUsvEZ42jeiYMQU7QlGkyFoCS5YQXL0agMjo0YSvvrogB/WI9iF/EkT+8/tpmDWL8PjxTtWOhQsJvvUWFMl8HtGMcJjQnDlnFqJsuOYawjNmSPISZ5FivqIwKEXkqqtIXHghofnzCWzejLemhvpZs7Cdmx+lJwqT5+hRQnPm4D12DBsMUj9rFvH+/d0OS+Qh+eeMKCix4cOdwR2dO+P9+GNKn3kG70cftX6gKAi+nTud7/TYMeLl5dR++cuSvMR5SQITBSfRqxd1995LrG9fPLW1lDz/PIGVKyEHZdFEO4nHCb7zDiWvvoqKRokOG0bd3Xdju3VzOzKRxySBiYJkS0up/+IXnedi1hJ87z1KXnoJVVfndmgiTer4cUqffZbAunXO865p02i46aaCm7wuck8SmChcHg+Rq66i7rbbsKEQvr17KX3qKRlqX0B827dT9vTTeA8fJtGlC3V33UW0qkpGmIqUSAITBS8+aBC1X/kKsT598NTWUvrCCwQXLZJRivmsoYHQ3LmUzJ2LCoeJDh5M7b33kmi9xJEQZ8goRFEUbJcu1N95J4FVqwgsX05gzRq8e/fScMMNJC66yO3wRCPe/fsJzZuH5+RJrM9HePp0opddJr0ukTZJYKJ4eDxEJk4k1q8fJW++iffoUUqffZbI+PFEJk4Eb/NV3UWOhMMEFy8msGkTAPFevaifORPbvbvLgYlCJQlMFJ1ERQW1X/kKwffeI7BuHcEVK/Dt2kXDddfJLSqXePfuJfT2206vy+MhMmkSkSuukH9UiDaRBCaKk99PePp0YoMHE5o/3+mNPfcc0TFjCE+eDMHcrKXV0anaWoLvvIN/xw4A4j17Ord1e/RwOTJRDCSBiaIW79uX2q9+1XkutnYtgfXr8e3cSXjKFGKXXirPXdpLIoF/wwaCy5ahwmGsz0fkyiuJjB0rvS6RNZLARPHz+4lMmULskksILViA9+BBSt54g9imTYSnTyfRq5fbERYV7/79BBcuxJtc/iY2cCAN11yDveAClyMTxUYSmOgwEj17UvelL+HbupXgkiX4DhzA9/TTRIcPJzx5MrZLF7dDLGieo0ed/6+7dwOQuOACwtOmEbv4YunpinYhCUx0LEoRGzmS2ODBBFeuxL9uHf5t2/Dt2EF07FgiV1yBLSlxO8qCok6dIrB8Of4tW1DWYv1+IhMmEKmqAp/8ihHtR/50iY4pFCI8dSqR0aMJLlmCf+dOAqtX49+wgUhVlfOsRhJZi9SpUwRWrcK/aRMqkXBGF44aRWTSJGxZmdvhiQ5AEpjo0GzXrjTMmkXk0CGCS5fi+/BDgitWEFi7luiYMUSqquSXcRPq008JrF3r9LjicaxSRIcOJTJ5MgmZ0yVySBKYEEDioouov+MOvAcOEFi2DN/+/U6PbN06osOHEx07lkR5udthuspz8CCBtWvx7dqFSlb+jw4dSmTSpA7//0a4QxKYEI3EKyupv/NOPB9/7Nwe++ADAps3E9i8mVi/fkTHjCE2aFDHGQoei6E2baL0nXfwHjoEgPV4iA4fTmTcOElcwlWSwIRoRqJ3bxpuvZXIsWOfD/TYvx/f/v0kysqcXtnIkcVZBslaPIcP49+yBf/77+PxevHW1WFDIaKXXUbk8stlFWyRFySBCdGCRPfuhGfMIPyFLzi/0LdswXv0KMHVqwmuXk28Vy+il1xCbNiwgv+l7jl2DN+OHfh27DgzhwvAXnwx9YMHExs2DPx+FyMU4mySwIRIRShEdNw4olVVzrOgLVvw7dyJ9/BhvIcPw6JFxHv1IjZ4sDPv6cIL3Y64dYkEnkOH8O3ejW/PHrw1NWc+siUlRC+9lOiIEQSHDyfW6DMh8oUkMCHSoRSJykoaKivh6qvx7d2Lb/t2fB9+eCaZBZcuxdOjB6HycmL9+hGvqMB26+b+ZF5r8Rw7hvejj/Du34+vuhpVX//5x4EAsSFDiA4dSrx//47znE80S2t9A/A44AWeMMY84nJI55AEJkSm/H5iQ4cSGzoUolG8+/bh370b7969qNpa/J98gn/7dgASpaUkevcm3rMniZ49iZeXO5U/PO20pmwigefTT/HU1OD55BO8hw7hPXQIFQ6fvVuXLsQuvpjYoEHE+/SRiccCAK21F/gdcC1QDazWWs8xxmxzN7KzyZ9WIbLB7yc+eDDxwYPBWkqspWH9erzV1XgPHsRTV4dn9+4zZZbAGc2X6NYN27Uric6dsZ06kSgrcyqBhELYYBDr9RIKJ0goLwqLsgmUTeCLR/BWH4NwGE9dHaq2FvXZZ3hOnMBz4gTqxAlUInFOmIlOnYhXVhLv149Y377Yrl3d7xmKfDQe+MAYswdAa/0ccCsgCUyIoqYU9OhBtKqKaFUVWIv69FPnFuMnn+A5csTpGdXWOoMlGg2YaM6UtYFmt5f+d6TF4xJdupAoLydeXk6iVy/ivXsX/EATkTOVwEeN3lcDE1yK5bwKPoGVZzgPxefzZXxsvpG25J9z2tGjBwwdevZO4TAcOwbHj6NOnoSTJ6G21nku1dDgvGIxxl9rUfE4KIX1eEB5sMEgJeV+CAaxpaXQqRN06uRUfO/eHbp1y9qIwWL5TkDa0pjWek2jt7ONMbMbvW+uW24zvlg7KfgEVpPh6Kjy8vKMj8030pb8k3I7/H4nubWwwOOjjzbTa/oM7v/6qfOf98SJFKJMTbF8JyBtOa2iogJjzLgWdqkG+jZ63wc4mNHF2lHBJzAhhBBZtxoYorUeCBwAvgTc425I52qnIVBCCCEKlTEmBnwPmAdsdzaZre5GdS7pgQkhhDiHMWYuMNftOFoiPTAhhBAFKWc9sNZmdWutg8BTQBVwFLjLGPNhruITQghRWHLSA2s0q/tGYDhwt9Z6eJPdvgV8aowZDPwr8OtcxCaEEKIw5eoW4plZ3caYCHB6VndjtwL/lfz5L8A1WmspESCEEKJZuUpgzc3qrjzfPskRMCeAAijpLYQQwg25egaWyqzulGZ+a63vA+4DMMZQUVGRcVBtOTbfSFvyT7ba8eij5/skd2WhiuU7AWlLMclVDyyVWd1n9tFa+4ALgGNNT2SMmW2MGZecRa4yfWmt17bl+Hx6SVvy71Us7ZC25O8rC20peLnqgaUyq3sO8DVgOfBFYKExJu9qbwkhhMgPOemBnW9Wt9b6Ya31Lcnd/gBcqLX+APgR8I+5iE0IIURhytk8sOZmdRtjHmz0cwNwZ67iAWa3vkvBkLbkn2JpB0hb8lUxtSUjylq5SyeEEKLwSCkpIYQQBalDF/PVWv8zzgTqBHAE+LoxJu/WvEmF1vo3wCwgAuwGvmGMOe5uVOnTWt8JPARcCow3xqxp+Yj801rZtEKhtf4jcDNwxBgz0u14MqW17otTpu4inL/rs40xj7sbVWa01iFgMRDE+f39F2PML9yNyj0dvQf2G2PMKGPMGOA14MHWDshjbwEjjTGjgJ3AT1yOJ1NbgNtx/pIWnBTLphWKJ4Eb3A4iC2LA/caYS4GJwN8V8HcSBq42xowGxgA3aK0nuhyTazp0D8wYc7LR2zLycMnsVBlj5jd6uwJnKkLBMcZsB9Baux1Kps6UTQPQWp8um7bN1agyYIxZrLUe4HYcbWWM+Rj4OPnzKa31dpzKP4X4nVjgs+Rbf/JVsL+32qpDJzAArfUvga/ilK6a7nI42fJN4L/dDqKDaq5s2gSXYhFNJBPy5cBKl0PJWLKXvxYYDPzOGFOwbWmrok9gWuu3ce59N/WAMeYVY8wDwANa65/gzFXL2/vJrbUluc8DOLdMnsllbOlIpR0FrLkKBx32X8j5RGvdCXgB+GGTuy8FxRgTB8ZorbsCL2mtRxpjtrgdlxuKPoEZY2akuOufgdfJ4wTWWlu01l/Deeh+TT5XMUnjOylEqZRNEzmmtfbjJK9njDEvuh1PNhhjjmut38V5TtkhE1iHHsShtR7S6O0twPtuxdJWyZFvPwZuMcbUuR1PB3ambJrWOoBTNm2OyzF1aMllmf4AbDfG/NbteNpCa90j2fNCa10CzKCAf2+1VYeeyKy1fgEYhjO0dh/wbWPMAXejykyyBFcQZzVrgBXGmG+7GFJGtNa3Af8O9ACOAxuMMde7G1V6tNYzgcdwhtH/0RjzS5dDyojW+llgGlAOHAZ+YYz5g6tBZUBr/QVgCbAZ5+86wE+T1YEKitZ6FM66iV6cDogxxjzsblTu6dAJTAghROHq0LcQhRBCFC5JYEIIIQqSJDAhhBAFSRKYEEKIgiQJTAghREGSBCaEEKIgSQITQghRkCSBCSGEKEhFXwtRiFzQWg8GtgO/arzAoNb6P4B7gemFuDinEPlMemBCZIEx5gPgCeAftNblAFrrB3GWtrlNkpcQ2SelpITIEq31RcBu4Pc4BVZnA3cbY4yrgQlRpCSBCZFFyQVS78e5Pf8DY8zvXA5JiKIltxCFyK5dOKsCLJfkJUT7kgQmRJZora8G/hNYDkzWWo92OSQhipokMCGyQGs9FngZZyDHNGA/8Cs3YxKi2EkCE6KNkkPo3wDmA983xkSAfwJmaq2nuBqcEEVMEpgQbZAceTgfZw7Yl40xp1f8fQpnJOIjbsUmRLGTUYhCCCEKkvTAhBBCFCRJYEIIIQqSJDAhhBAFSRKYEEKIgiQJTAghREGSBCaEEKIgSQITQghRkCSBCSGEKEiSwIQQQhSk/x/DbsWnWGvpwQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n  = 1\n",
    "xs = stats.norm.rvs(size=n)\n",
    "\n",
    "def ll(x, xs):\n",
    "    dx = (x-xs)\n",
    "    return np.sum(dx*dx)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.hist(xs, bins=50, alpha=0.5, range=(-3., 3.), color='blue')\n",
    "axb = ax.twinx()\n",
    "\n",
    "mus = np.linspace(-3., 3., 100)\n",
    "lls = np.array([ll(mu, xs) for mu in mus])\n",
    "axb.plot(mus, lls, alpha=0.5, color='red', lw=2)\n",
    "ax.set_xlabel('$x$', fontsize=16); ax.set_ylabel('samples');\n",
    "axb.set_ylabel('$-2 \\, \\log \\mathcal{L}(x | \\mu)$', color='red', fontsize=16)\n",
    "fig.savefig('figs/gen_likelihood_example.png');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "\n",
    "\n",
    "### 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 hypothesis $h$, with a prior probability $\\pi(h)$. Given some observation, $x$, Bayes' theorem allow us to compute the **posterior probability**, $\\pi'(h)$. The probability of $h$ given the data $x$, $p(h|x)$, is: \n",
    "\n",
    "$$\n",
    "p(h|x) = \\frac{p(x|h) \\pi(h)}{\\int p(x|h) \\pi(h) \\mathrm{d}h}\n",
    "$$\n",
    "\n",
    "As the total probability for the ensemble of hypothesis $h$ should be one, we can obtain the posterior probability $\\pi'(h)$, just normalizing the conditional probability of each hypothesis:\n",
    "\n",
    "$$\n",
    "\\pi'(h) = \\frac{p(h|x)}{ \\int p(h|x) \\mathrm{d}h}\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}$\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": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "posterior integral 0.995\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEQCAYAAABMXyhMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl4U2Xax/FvmnSjLUtpEQE3FBUHEUTBbQQVQRBB1HkEUcQFxhF0FB03HHVQ0HEdVFwQFXHDR0FEB0VfFRzHBURFKYsiMgqo0NIFui95/zgBSymQtElPmv4+15WrzenJyf0kbe4+u8fv9yMiIhKMOLcDEBGRxkNJQ0REgqakISIiQVPSEBGRoClpiIhI0JQ0REQkaEoaIiISNCUNEREJmpKGiIgEzed2ABGgKe4iIqHzBHNSLCYNNm7cWKfHZWRkkJ2dHeZo3BErZYmVcoDKEo1ipRxQv7K0a9cu6HPVPCUiIkFT0hARkaApaYiISNCUNEREJGhKGiIiEjTXRk8ZY54BBgGbrLVdavm5B5gCDASKgFHW2i8bNkoRcVtlZSUlJSUAeDxBjQoN2m+//UZpaWlYr+mW3ZXF7/fj9XpJSkoKy/O4OeR2BvAoMHM3Px8AdArcegGPB76KSBNRWVlJcXExKSkpYU8YAD6fD6/XG/brumFPZSkpKaG8vJz4+Ph6P49rzVPW2o+ALXs4ZQgw01rrt9Z+BrQ0xuzbMNGJSDQoKSmJWMJoShITEykrKwvLtaK5T6M98HO1++sDx0SkCVHCqL9wvobRPCO8tlLWukSIMWYMMAbAWktGRkadntDn89X5sdEmVsoSK+WA8JUlMTFxx/dutcc31PtSWVmJzxfZj6lIX78h7aksSUlJYXnPovnVWg/sV+1+B6DW9UGstdOAaYG7/rpOpdeSAtEnVsoBkSmLW69NQ70vpaWlEe1z8Pl8VFRUROz6DWlvZSktLd3texbKMiLRnDTmAeOMMbNwOsDzrbW/uByTiEiT5uaQ25eBPkCGMWY9cDsQD2CtfQKYjzPcdg3OkNtL3IlURES2cy1pWGuH7+XnfmBsA4UjIuK6vLw8evfuzRtvvMGBBx4Y9ONGjx5Nr169uPzyyyMXXEA0j54SEYlpDzzwANOmTdtx/5FHHuHUU08NKWEAjB8/noceeoiCgoIwR7grJQ0REZe8/fbbHH/88QAUFxfz8ssvM3z4HhthatW5c2cOOOAA5syZE+4Qd6GkISISAVu2bKF9+/ZMmzaNgQMH0rFjR0466SQWLVoEOMt+pKSkcOSRRwLw/vvv4/F4OPbYY3e6zhdffEGHDh0oLCzccSwvL4/27duTlZW141j//v2ZO3duxMulpCEiEgHLly8HYMaMGdx888289957dO7cmbFjx1JcXMyiRYu45JLfx/csXryYrl277jIRLysri4MOOoiUlJSdrp2QkMChhx6641j37t35+uuvKS4ujmi5lDRERCIgKysLr9fLzJkz+eMf/8jBBx/MhAkTyM3N5YcffmDdunWceeaZO85fv349bdq0qfU622sj2y1fvpxOnTrttJbUPvvsQ3l5Ob/99lvkCkV0z9MQEdlJ+/bBT0ILpw0bap1XvEdZWVmcfvrpHHLIITuOVf+Qv+GGG3Y6v6SkpNYZ2ytWrGDgwIE7HVu+fDl/+MMfdjqWnJy84zqRpJqGiEgEZGVl0aXLzrs+LFu2jKSkJA4++OBdzk9PTyc/P3+nY1VVVaxatarW69RMGrm5uQC0bt06HOHvlmoaItJo1OU//j2J1DIiJSUl/PDDD1RVVe10fPr06QwePHhHraC6Ll26YK3d6dgPP/xAcXExbdu23XFs5cqVrF27dpeksWrVKtq2bUtmZmYYS7Ir1TRERMJs1apVAMydO5fPP/+cNWvWcNVVV7Fu3TpuvvnmWh/Tu3dvvv/+e7Zs+X3HiO2jo5599lnWrl3LRx99xJVXXgmwy1Lnn3/+OX369IlAaXampCEiEmbbRzxdf/31XHnllfTv35/CwkL+/e9/19rZDc5ci27dujFv3rydrnPyySfzyy+/cNppp/GPf/yDm266iZYtW/LMM8/sOK+kpIT58+dzwQUXRLxsap4SEQmzrKwsDj/8cIYMGcKQIUOCftz48eO57bbbuOiii/B6vWRlZdG1a9ddaifV52cAzJo1i6OPPpoePXqEJf49UU1DRCTMsrKy6Ny5c8iPO+WUUxg1ahS//OIs6L1ixYqgruPz+Zg8eXLIz1cXShoiImHk9/tZtWoVRxxxRJ0ef9lll9GhQwc2b97M5s2bg0oaF1544U5DeyNJzVMiImHk8XhYvXp1va+TmZnJhg0bwhBReKmmISIiQVPSEBGRoClpiIhI0JQ0REQkaEoaIiISNCUNEREJmpKGiIgETUlDRESCpqQhIiJBU9IQEZGgaRkREWlU2rdv36DPF41LebhJNQ0REQmakoaISAS8+eabHHTQQaxfv37Hsdtuu40TTjiBzZs3uxhZ/ShpiIhEwKBBgzj88MOZMmUKAE888QRz587lhRdeiPg+3pGkPg0RkQjweDzcdNNNjBw5kgMOOICHH34Yay0dO3Z0O7R6UdIQEYmQ3r17c9RRR3HvvfcyY8YMunXr5nZI9abmKRGRCPn4449ZsWIFfr+/UTdJVaekISISAVlZWYwePZq77rqLM844g7vvvtvtkMLCteYpY8wZwBTAC0y31t5T4+f7A88BLQPn3GStnd/ggYqIhGj9+vWMHDmSMWPGMGzYMLp160bfvn355JNPOOGEE9wOr15cqWkYY7zAVGAAcAQw3BhTcxf2WwFrre0ODAMea9goRURCl5uby4gRIzjttNO49tprATj88MMZNGhQTNQ23Kpp9ATWWGvXAhhjZgFDgBXVzvEDzQPftwA2NmiEIhKVwjlD2+fzUVFREbbrAbRq1YpFixbtcvyJJ54I6/O4xa2k0R74udr99UCvGufcAbxrjLkKSAH6NkxoIiKyO24lDU8tx/w17g8HZlhrHzDGHA88b4zpYq2tqvlAY8wYYAyAtZaMjIw6BeXz+er82GgTK2WJlXJAZMri1mvTUO/Lb7/9hs8X2Y+pSF+/Ie2pLImJiWF5z9x6tdYD+1W734Fdm58uA84AsNZ+aoxJAjKATTUvZq2dBkwL3PVnZ2fXKaiMjAzq+thoEytliZVyQGTK4tZr01DvS2lpKV6vN2LXj0TzlFv2VpbS0tLdvmft2rUL/nlCjiw8lgCdjDEHARtwOrovqHHOT8BpwAxjTGcgCWi8C7aIiMQAV0ZPWWsrgHHAAmClc8hmGWMmGmMGB067DhhtjFkGvAyMstbWbMISEZEG5PH7Y+5z2L9xY90GWqkpJPrESjkgfGWpvp+EW3s9NNT7UlhYSEpKSsSu35Sap/b0Wgaap2rra96FZoSLSNTyer2UlJQQg//cNqiKigri4sLzcR87wwZEJOYkJSVRXl5OUVER4KwcG06JiYmUlpaG9Zpu2V1Z/H4/cXFxJCUlheV5lDREJKrFx8cTHx8fkWur+TN0ap4SEZGgKWmIiEjQlDRERCRoShoiIhI0JQ0REQmakoaIiARNSUNERIIW9DwNY8w+QD/gKJwtWPOAZcB71tpfIxOeiIhEk70mjcAKs3cCpwBLcRYY/BVIAy4C/mWM+RC4zVq7YrcXEhGRRi+YmsYM4D5ghLV2lznqxpgEnK1anwaOD2t0IiISVfaaNKy1O7ZhNcYkWGvLavy8DHg1cBMRkRgW6tpT24wxq3D6Mr4OfF0HTLDWXhLm2EREJMqEmjTaAN0Ct6NwNlLaH1gb5rhERCQKhZQ0rLV5wMLADQBjzF3AL2GNSkREolI45mncBdwUhuuIiEiUC6mmYYyZitOX8TXwrbW2BNgXKIpAbCIiEmVC7dPYAJwK/BXoaIzZgNPPscAYMxRnDsf31trK8IYpIiLRINQ+jcnbvzfGxANHAEcGbqMDXzOB8OwrKCIiUaXO271aa8txhtwuq37cGNOyvkGJiEh02mtHuDHmamNM4l7OSTTGXA07RliJiEgMCqam0RZYY4yZDywCVgNbcdaeOhToAwwAZkYoRhERiRJ7rWlYa28BugPfA5cBbwPLgfnApcAqoLu19tYIxikiIlEgqD4Na202cH/gJiIiTVQwfRrNjDFpDRGMiIhEt2BmhJ8KLDXGTDfGHBvpgEREJHoF06fxlrX2UJw9NY4xxlwf+bBERCQaBT1Pw1q7GmfkVK37aoiISOwLde2p/sATwP7GmG04I6luDyQUERGJcaGucjsNeBRn7sYJOMNtFxljeu3xUSIiEhNCXUYkAfhXYEHCzUCWMWYJ8AjQM9zBiYhIdAk1aSwCzgZmVzv2DvByqE9sjDkDmAJ4genW2ntqOccAdwB+YJm19oJQn0dERMIn1KSxHzDTGNMNJ3FsBi4A/i+UixhjvMBU4HRgPbDEGDPPWrui2jmdgJuBE621ucaYNiHGKiIiYRZq0piKs/x5d2AkThKpAF40xvwVZz+NVdban/ZynZ7AGmvtWgBjzCxgCLCi2jmjganW2lwAa+2mEGMVEZEwC3U/jZeq3zfGtOD3/TSOBM4D/gCk7+VS7YGfq91fD9TsTD808Bz/xWnCusNa+04o8YqISHjVeT8NAGttPvBx4BYKTy3H/DXu+4BOOKvodgD+Y4zpUtvS68aYMcCYQExkZGSEGE7gCX2+Oj822sRKWWKlHBCZsrj12sTK+xIr5YCGK0u9kkY9rMdp2tquA7CxlnM+C2z29KMxZjVOEllS82LW2mk4w4EB/NnZ2XUKKiMjg7o+NtrESllipRwQmbK49drEyvsSK+WA+pWlXbt2QZ/rVtJYAnQyxhyEs+/4MJwO9ermAsOBGcaYDJzmqrUNGqWIiOwk1Ml9YWGtrQDGAQtwOs+ttTbLGDPRGDM4cNoCIMcYswL4EPibtTbHjXhFRMTh8ftrdiU0ev6NG2u2dAVHVdXoEyvlgPCVpX379ju+37BhQ72vVxex8r7ESjkgLM1TtfU176LezVPGmGdwOsKfC8wUFxGRGBWO5ikPTn/EN2G4loiIRLGQkoYxZpfzrbWXWGv7At3CFpWIiESloJNGYOmPQmNMYm0/DwyNFRGRGBZ00gj0V3wHtI5cOCIiEs1C7Qh/EXjLGDMFZ/LdjqFX1toPwhmYiIhEn1CTxl8CX++ocdwPdKx3NCIiEtVCXbDwoEgFIiIi0S/keRqBfS6G46xUuwF42Vr7fbgDExGR6BPqkNuzgKXA4cAW4DDgi2pLf4iISAwLtaYxGRhirf1w+wFjTB/gUWBeGOMSEZEoFOqM8A7Af2oc+zhwXEREYlyoSeNr4Loax8YHjouISIwLtXnqSmBeYD/wn4H9gW2A+jRERJqAkGoa1tqVQGfAAA8AfwKOCBwXEZEYt9ekYYw5udr3pwInAwlAduDrHwPHRUQkxgXTPPUY0CXw/dO7OUczwkVEmoC9Jg1rbZdqdw/RRksiIk1XqEujb9vd0ugiIhL7tDS6iIgETUuji4hI0LQ0uoiIBE1Lo4uISNBCXUZERESasLrsp3E6zn4amdbas4wxxwDN1achIhL7Qt1P4yrgcZxRVNtnihcDd4U5LhERiUKhNk9dA/S11t4DVAWOrcLZjElERGJcqEkjDWd1W/h9uG08UBa2iEREJGqFmjQ+Am6qcexq4MNazhURkRgTakf4VcCbxpjRQJoxZjVQAJwV9shERCTqhJo0RgDHAj1xNmD6GViM09fxYHhDExGRaBNq0rjNWns/8HngBoAx5laUNEREYl5QSaPaJkteY8wpgKfajzsCW0N9YmPMGcAUwAtMD4zIqu2884BXgWOttV+E+jwiIhI+wdY0tm++lAQ8U+24H/gVp68jaIFl1qcCp+MsfLjEGDPPWruixnlpOB3tn+96FRERaWhBJY3ta04ZY2Zaa0eG4Xl7AmustWsD150FDAFW1DjvTuBe4PowPKeIiNRTqAsWhiNhALTn9/ke4NQ2elU/wRjTHdjPWvuWMUZJQ5qkvDwPq1fH8913Pr7/3sfmzV7gHaAl4OXss1uTlAStWlWx//4V7L9/JYccUkHXrmUkJ7scvMSkkJJGoD9jnbX2R2PMvsA9QCVwi7X21xAu5anl2I69OYwxccBDwKgg4xoDjAGw1pKRkRFCKL/z+Xx1fmy0iZWyxEo5ILiylJfDBx94eP/9OBYu9LBsWW1Tqfrv+G7JktqvEx/vp3t3Pyed5Gfw4Cp69fITF8blSWPlfYmVckDDlSXU0VOP8ftv7AOBrxXANGBwCNdZD+xX7X4HYGO1+2lAF2ChMQagLTDPGDO4ts5wa+20QAwA/uzs7BBC+V1GRgZ1fWy0iZWyxEo5YPdl8fthyZIE5sxJ5q23ksjN9e74WWKin86dy+nUqYJOnSpo376SsWMvAPKACmbPfpuSEg85OXH89JOXdet8rFgRz8qVPhYvjmPxYnjwQS9t21YycGAxF1xQROfOFRErS2MTK+WA+pWlXbt2QZ8batJob639yRjjw0keB+AsIbJxzw/bxRKgkzHmIGADMAy4YPsPrbX5wI6UaYxZCFyv0VMSS0pLYe7cZJ56KpWVK+N3HD/ssHL69SvhpJNK6dFj12amsWPf2fH9ccfVvoLP1q0evvwygQ8/TGT+/CQ2bPDxzDOpPPNMKj17ljJqVCFnnlmCL+R1rqWpC/VXpsAYsw9OLWCFtXabMSYBZ/2poFlrK4wx44AFOENun7HWZhljJgJfWGvnhRiXSKNRUgLPP5/C1KmpgT4KyMysxJgihg4tDktNIC3NT+/epfTuXcrttxfwzTfxWNuM115LZvHiRBYvTuTAAyu45pqtDB1arOQhQQv1V+URnFpCAs4scIATcVa6DYm1dj4wv8ax23Zzbp9Qry8SbSor4eWXm/Hgg6ls3Oj86R1xRDmjR29jyJBiEhMj87weDxx1VDlHHZXPzTcXMHt2MtOmpbJunY9rrmnFlClpTJhQwBlnlOCprbdRpJqQusastf8E+gInWmtnBQ5vAC4Pd2AisWTp0niOP97H9de3ZONGH507lzNjRg7vvrsZYyKXMGpKTfVz8cVFLFq0iYceyuXAAyv48Ucfl1+ezrBhrVm1SlUO2bO6/Ib4gRHGmPY4CWOWtfbb8IYlEhtycz1Mntycl15KAaB9+wpuvnkrQ4YUh3U0U6h8PjCmmHPOKeaFF5px333N+fjjRPr1y+TPf97G+PFbNWRXahXqzn1nAUuBw4EtOJsvLTHGhDJySqRJeO+9RE49tQ0vvZRCfLyfG26oZOHCzQwd6m7CqM7ng1GjivjPf37j4osL8fvhscfS6NevDUuWJLgdnkShUGsak4Eh1tod+2cYY/oAjwLqvBYBCgo83HFHC155pRkAPXuWct99+Rx3XEuys/17ebQ70tP9TJ6cz7nnFnHddS35/vt4hg5tzRVXFHLDDQUkKH9IQKj/73QA/lPj2MeB4yJN3rJl8ZxxRiavvNKMxEQ/t92Wz2uv5XDIIfUfEdUQevQoZ8GCzVx99Vbi4uDxx1M5++wM1q717v3B0iSEmjS+Bq6rcWx84LhIk+X3w/TpKQwZksH//uejS5cyFizYzJ//XIi3kX3eJibCjTduZc6cbDp0qGDZsgT698/k9dfVySGhN0/9BWfnvr/irB21P7CN0GaDi8SUwkIP11zTkvnznQ/VUaMK+fvf80lKcjmwejrmmHLefXczN93Uknnzkhk3rhVffRXP3/9e4HZo4qJQh9yuAjoDBrgfOA84wlq7MgKxiUS9deu8nHVWBvPnJ5OWVsWTT25h0qTGnzC2a9HCz2OP5XL33XnEx/t5+ulUjGnNL7+4HZm4JdTRUwnAbcCzwIzA7TZjTIz8iYgEb9GiRAYOzGT16ngOOaSct97azKBBJW6HFXYeD4wcWcTs2dm0bVvJ4sWJnHhiPMuXa05HUxRqn8bjwKk4GyMdG/jaG2chQ5Em47nnmnHRRenk58fRr18xb72VzSGHVLodVkT16FHOO+9spmfPUjZs8HD22Rm8/bb+X2xqQk0aZwODrLVvW2tXWGvfDhw7O/yhiUSfykq4447m3HJLSyorPVx99VaefjqXtLToHEobbpmZVcyalcNFF1VSXBzH5ZenM3VqKv6mUXwh9KTxK9CsxrFkQC2cEvOKiz2MHt2Kp55KJT7ez0MP5XLjjVujZqJeQ0lMhKeequSWWwrwePxMntycW25pQWVsV7QkINRGyeeBd4wxj/D7nhhjgZnGmFO3n2St/SB8IYq4b8sWD6NGtWbp0gRatqxi+vQtHH987cuSNwUeD4wdu43996/gr39txcyZKfz2WxxTp+aRnKxqRywLNWn8OfD1lhrHrwjcwFmbqmN9ghKJJuvXexkxIp01a+Jp376CF1/cQqdOjWOyXqSddVYJbdrkcOml6SxYkMz553t57rkcWrVS4ohVoe4RflCkAhGJRqtX+7jggtb8+quXzp3Lef75HPbdt8rtsKJKr15lzJ2bzYgR6SxdmsC552bw4ot6nWJVE2uNFQneV1/Fc845Gfz6q5fjjitl9uxsfRDuRqdOFcydm02nTuWsXh3P2Wdn8MMPjWwqvARFSUOkFh9/nMD557cmLy+Ovn1LeOGFHFq0UJPLnrRrV8WcOdl0717G+vU+zjkngxUrNJcj1ihpiNTw7ruJjBzZmsLCOM45p4jp07dob4kgpaf7eeWVHP74x1Kys72cd14GX34Z0m7QEuWUNESqmTcvidGj0ykt9XDxxYVMmZJHvD7zQpKS4mfGjBz69y8mPz+OYcNa88knWls9VihpiARYm8zYsa2oqPDwl79sY9Kk/CY3ByNckpLgySdzOeecIgoL47jootYsXNhAe9pKROlPQgR44YVmXHttK6qqPFx/fQETJhTg8bgdVeMWHw9TpuQxYkQhJSUeLrkknffeU+Jo7JQ0pMl79tlm3HhjSwBuvTWfa6/dpoQRJnFx8M9/5nPppdsoK/Nw+eXp/PvfWq+qMVPSkCbtySdTuPVWJ2FMnJjPX/5S6HJEscfjgYkTC7jiim2Bpr9WvPGGEkdjpaQhTdbUqalMnNgCgLvvzuOyy5QwIsXjgVtvLeDqq7dSWelh3LhWzJmjIWmNkQZRS5M0ZUoq997bHI/Hz3335TN8eJHbIcU8jwduuGErPh88+GAaV1/dkooKMKbY7dAkBEoa0uQ8+GAqDzzgJIwHH8zTh1YD8njguuu24vX6ue++5owf35KqKhg2TO9BY6GkIU2G3w/335/Gv/6VRlycn3/9K49zz9WHlRuuuWYbPh/cfXdzrruuFZWVHkaMUG2vMVCfhjQJfj/8859OwvB6/TzyiBKG28aN28att+YDcMMNLZk5s+ZWPRKNVNOQmOf3w913pzF1qpMwpk7N5ayzYm8v78boL38pxOuFf/yjBTff7DRVjRqlGkc0U01DYprfD3fe2ZypU9Pw+fw8/rgSRrQZM6aQiROdGseECS2ZPj3F5YhkT5Q0JGb5/XD77c158klne9YnnsjlzDOVMKLRZZcVMmlSHgC3396CJ59U4ohWShoSk/x+uPXWFjz9dCoJCX6mTdvCgAFKGNFs1Kgi7rnHSRwTJ7Zg6tRUlyOS2ihpSMypqoIbb2zBjBkpJCb6mT59C/36lbodlgThoouKuP/+PDweP5MnN+ehh5Q4oo1rHeHGmDOAKYAXmG6tvafGz8cDlwMVwGbgUmvt/xo8UGlUKivhuuta8uqrzUhK8vP001vo00cJozEZPrwIn8/P+PEtuf/+5pSXe/jb37ZqPbAo4UpNwxjjBaYCA4AjgOHGmCNqnPYVcIy1tivwGnBvw0YpjU1FBVx9tZMwkpOrmDkzRwmjkfrTn4p5+OE8vF4/U6akcdddzfFr48So4FZNoyewxlq7FsAYMwsYAqzYfoK19sNq538GXNigEUqjUloKY8e24u23k0lNreL557fQs2eZ22FJPQwdWkx8vJ+xY1vxxBOpFBd7uOsu7XHiNreSRnvg52r31wO99nD+ZcDbu/uhMWYMMAbAWktGRkadgvL5fHV+bLSJlbIEU46iIrjkEh/vvhtHy5Z+3nyzkp49mzdQhMGLxHvi1nvcUL9fo0ZBZmYFw4f7eO65FCCJxx+vxOsNz/Vj5e8EGq4sbiWN2lona618GmMuBI4Beu/uYtbaacC07dfJzs6uU1AZGRnU9bHRJlbKsrdybNvmYdSodD79NI7WrSt5+eUcOnasIBqLHon3xK33uCF/v3r1ghkzErj00nSee87Lli1lPPxwLglh2EE2Vv5OoH5ladeuXdDnulXRWw/sV+1+B2BjzZOMMX2BCcBga60ap2UnW7Z4OP/81nz6aSJt21Yye3YOf/hDhdthSQScfHIZL720hbS0Kt58M5nLLkunuFg9425wK2ksAToZYw4yxiQAw4B51U8wxnQHnsRJGJtciFGi2C+/xHHuuRl8/XUC++9fwZw52XTqpIQRy3r2LOPVV3NIT6/kgw+SGDEinYICJY6G5krSsNZWAOOABcBK55DNMsZMNMYMDpx2H5AKvGqM+doYM283l5Mm5scfvZx9dgbffRfPYYeV8/rr2RxwQKXbYUkDOPLIcubMyaFt20o+/zyR887LYPNm9Yw3JI8/9sax+Tdu3KWlKyhq34w+NcvxzTfxXHhhOjk5Xrp3L2PmzBzS0xvH73C43pP27dvv+H7Dhg31vl5duP379fPPXoYPb82PP/o48MAKXnopp07/OLhdjnAKQ59GUNU2pWhpND76KJHzzmtNTo6X3r1LeOWVxpMwJLz226+SuXOzOfLIMtat83H22RksX65FuxuCkoY0Cq+/nszIkekUFsZxzjlFzJixhZQUJYymLCOjildfzeHEE0vZtMnLuedm8NFHiW6HFfOUNCSq+f3w8MOpjBvXivJyD2PGbGPKlLywDLeUxi8tzc/zz+cwZEgR27bFcdFF6bzySrLbYcU0JQ2JWuXlcOWVXv75T2c/7zvuyOf22ws0I1h2kpgIjz6ax5VXbqWiwsP48a144IE0LTsSIfrzk6iUn+9h5Mh0nnnGS1KSn6eeymX06EK3w5IoFRcHEyZsZdLQVQyZAAAO4UlEQVQkZ4XcBx9MY9y4lpRoNfywU9KQqPPjj17OOiuDjz5KIjPTj7XZ2gtDgjJqVBHPPruFlJQq5s5thjEZZGfrYy6c9GpKVPnvfxMYNCiTH36Ip3Pncv7733J69Ch3OyxpRE4/vZS5c7Np166CpUsTGDhQI6vCSUlDooLfD9OnpzB8eGvy8uI4/fQS5s7N5oAD3I5MGqMjjqjgrbey6d69jA0bfAwZksEbbyS5HVZMUNIQ1xUXe7j66pbcfnsLKis9jBu3laef3kJqqnoype722aeK117L5vzziygpiePKK9O5667mVGi1mXpR0hBXrV3rZfDgDObMaUazZlU8+eQWbr55a9iWvpamLSkJHnggj7vucjZ0evzxVIYNa82mTfroqyu9cuKaN99MYsCATFasiOfAA53mhEGD1OEt4eXxwCWXFGFtDm3aVPLpp4n075/JJ59osk9dKGlIgysuhgkTWnDFFels2xbHoEHFvPPOZg47TO0GEjnHHVfGggWbOf54Zwb5+ee3ZuJEr5qrQqSkIQ1q5UofgwZlMmNGCvHxfiZNyuOJJ3JJS1P/hURemzZVzJqVw1VXbcXvh0mTvPzpT63ZsEHtocFS0pAGUVXljI4688xMVq2Kp2PHCubNy2bUqCI82hJBGpDPBzfdtJVZs3LYd18/ixcn0rdvJq+9lqxZ5EFQ0pCI+9//vBjTmttvb0FpqYcRIwpZsGAzXbtq/oW456STyvjii3JOP72EgoI4/vrXVowe3YqcHH0s7oleHYmYqiqYMaMZfftm8umniWRkVPLMM1u49958mjXTv3TivowMePbZLTzwQC6pqVW8/XYyffpk8vrrqnXsjpKGRMSKFc6EqgkTWlJUFMeQIUV8+OFm+vfX6CiJLh4PDBtWzPvvb+bEE0vZssXLuHGtGDkynfXr1ddRk5KGhFVhoYdJk9IYMCCTL79MYJ99KnnqqS089lge6elVbocnslsdOlTyyis53H9/Hi1aVPHBB0n06ZPJI4+kUlrqdnTRQ0lDwsLvh9mzkzn55DY89lgalZUwalQhCxduYuBA1S6kcfB4YPjwIhYu3MSgQcUUF8dxzz3NOe20Nnz4oTZ4AiUNCYMlSxIYMiSDq69uxa+/eunWrYx587KZNCmf5s3VMCyNT5s2VTz5ZC4vv5zDwQeX8+OPPi68sDUjRqSzcmXTXvxQSUPq7LvvfFx6aSvOPjuDpUsTyMys5MEHc3nzzWyOPlojo6TxO/nkUv7v/zbz97/nk5ZWxcKFSfTrl8l117Vosv0dShoSsjVrfFx1VUtOOy2TBQuSSU6u4pprtvKf/2zi/POLtbOexJSEBLjiikI++WQTl166jbg4mDUrhZNOasMtt7Tgl1+a1i980yqt1MuKFT7Gjm1Jnz6ZzJnTjLg4uOgi54/pb3/bqlndEtPS06u4884CPvxwE0OHFlFRAc89l8KJJ+7DjTe24Mcfm0bNQ0lD9sjvh48+SuSCC9I5/fQ2zJ3bDJ8PLrywkI8/3sQ99+TTpo1GRUnT0bFjJY8+mscHH2xm0KBiSks9vPBCCief3IYrrmjF0qXxMT3Ho2n36MhuFRR4mD07meeeS+H77+MBaNasiuHDi/jznwtp377S5QhF3HXooRU8+WQua9Zs5fHHU5g9uxlvvpnMm28mc9RRZVxySSFnnVVMUozt/aSkITv4/bB4cQLWJjNvXjJFRU5FdJ99KrnkkkIuvLCQVq1i+F8okTo45JAKHnggn+uu28qMGSm8+GIKy5YlcM01CdxxRwuGDi1i2LAiunSJjeV0lTSE777zMW9eMq+/nsy6db//SpxwQimjRhXSr18J8fEuBijSCLRrV8Utt2zl2mu38sYbTi39m28SePbZVJ59NpXOncsZOrSYIUOK6dCh8dbUlTSaIL8fsrJ8LFiQzPz5Saxa9XtGaNu2kvPOK8KYIg4+uPH+You4JTnZWZZk2LBili/3MWtWM15/vRkrV8azcmU8kyc3p3v3MgYMKOGMM4ob3d+ZkkYTUVDg4b//TWTRokQ++CCRDRt+f+tbtKhiwIBiBg8u4aSTSrXVqkiYdOlSwV13FfD3vxewaFEir7/ejHffTeSrrxL46qsEJk9uTseOFZxySgmnnFLKcceVkpzsdtR7pqQRo/Lz4YMPEvnsswQ+/TSRZcviqaz8feOKffappG/fEvr1K+Hkk0tJ0M6XIhGTmAj9+pXSr18pRUUeFi5M5O23k3j//STWrvWxdm0qTz+dSkKCn6OPLuOEE8o47rhSuncvj7oVoZU0YkBJCaxaFc+338bz9dfxfPVVAt9958Pvb73jHK/XT8+epfTuXUqfPqV07VquSXgiLmjWzM/AgSUMHFhCRQUsXZrAhx8m8uGHiWRlxfPZZ4l89lkikIbX6+eII8rp3r2crl3L6dq1jEMPrXC1j1FJoxEpLoZ163ysXevju+98rF4dz3ff+VizxrdTLQIgIcFPly5lHH98KccdV8axx5Zp8p1IlPH5oFevMnr1KuOmm7aSm+vhs88S+eSTBL74IoGsrHi+/TaBb7/9vSkgIcHPwQdX0LlzOYcdVkGnThUcfHA5LVo0UMwN8zS7MsacAUwBvMB0a+09NX6eCMwEegA5wPnW2nUNHWdD8fth2zYPv/0Wxy+/eNm40bn9/LOPn37y8tNPzn2/f9e9UePi/Bx2WDl/+IPz30iPHmX07t2CrVuzXSiJiNRVq1Z+BgwoYcAAZ2XooiIPX38dz7Jl8XzzTQLffBPPunW+HZ3q1fl8ft54I55u3SK77psrScMY4wWmAqcD64Elxph51toV1U67DMi11h5ijBkG/BM4v+GjDY3fD6WlsG1bHNu2edi6NY6CAg8FBXHk53vIzY0jLy+O3Nw4srPjyMnxkp0dx6ZNcTvmReyOz+dn//0rOOgg57+Lww4r5/DDne+Tk3euRSQmwtatkSypiERas2Z+TjjB6eOAQsD553L1ah+rVjktDT/84Nx+/tnbIEN53app9ATWWGvXAhhjZgFDgOpJYwhwR+D714BHjTEea23Y21h++y2O119PJjExjry8VCorPZSXQ3m587W01ENZmYeyMigr81Ba6qGkxLkVFzu3khIPRUUeCgs9uzQVBSspqYp99qli330radfOuXXoUMkBB1Sy334VdOhQqfkSIk1caqqfHj3K6dFj5xpFSkoG27ZFfkkft5JGe+DnavfXA712d461tsIYkw+0BsLe5rJxo5c779zeINi83tdLSPCTmlpFaqqf1FQ/LVpU0bx5Fc2b+2nVqoqWLZ1bRoZza926ijZtKklL8+OpW74RkSYuORkKCyP/PG4ljdo+GmvWIII5BwBjzBhgDIC1lnbt2oUUTLt2hHmBMQ9OV427Qn0dolWslAPCUxZ/lKyGFyvvS6yUAxqmLG4NulwP7Fftfgdg4+7OMcb4gBbAltouZq2dZq09xlp7DM4ndp1uxpil9Xl8NN1ipSyxUg6VJTpvsVKOMJUlKG7VNJYAnYwxBwEbgGHABTXOmQdcDHwKnAd8EIn+DBERCZ4rNQ1rbQUwDlgArHQO2SxjzERjzODAaU8DrY0xa4DxwE1uxCoiIr9zbZ6GtXY+ML/GsduqfV8C/KmBw5rWwM8XSbFSllgpB6gs0ShWygENVBZPtHSqiYhI9NPqQyIiEjStPVWDMeZOnImFVcAmYJS1tubIrqhnjLkPOAsoA34ALrHW5rkbVd0YY/6EM9GzM9DTWvuFuxGFZm9L5jQmxphngEHAJmttF7fjqStjzH44yxS1xflbn2atneJuVHVjjEkCPgIScT7TX7PW3h6p51NNY1f3WWu7Wmu7AW8Bt+3tAVHqPaCLtbYr8B1ws8vx1Mdy4BycP4xGpdqSOQOAI4Dhxpgj3I2qXmYAZ7gdRBhUANdZazsDxwFjG/H7Ugqcaq09CugGnGGMOS5ST6aaRg3W2oJqd1PYzYTCaGetfbfa3c9whi03StbalQDGGLdDqYtglsxpNKy1HxljDnQ7jvqy1v4C/BL4fqsxZiXOKhSN7n0JTEXYFrgbH7hF7HNLSaMWxphJwEggHzjF5XDC4VLgFbeDaKKCWTJHXBRIgt2Bz10Opc4CNdqlwCHAVGttxMrSJJOGMeb/cNoya5pgrX3DWjsBmGCMuRlnPknE2gfrY2/lCJwzAacq/mJDxhaqYMrSSNU207ZR1l5jkTEmFZgNXFOjlaFRsdZWAt2MMS2B140xXay1yyPxXE0yaVhr+wZ56kvAv4nSpLG3chhjLsbptDwt2mfTh/CeNDbBLJkjLjDGxOMkjBettXPcjiccrLV5xpiFOP1OEUka6givwRjTqdrdwcAqt2Kpj8CInRuBwdbaIrfjacJ2LJljjEnAWTJnnssxNXnGGA/OqhMrrbUPuh1PfRhjMgM1DIwxyUBfIvi5pcl9NRhjZgOH4QzD+x9whbV2g7tRhS6w/Eoizq6HAJ9Za69wMaQ6M8YMBR4BMoE84GtrbX93owqeMWYg8C+cIbfPWGsnuRxSnRljXgb6ABnAb8Dt1tqnXQ2qDowxJwH/Ab7F+VsHuCWwUkWjYozpCjyH8/sVh7Ms08RIPZ+ShoiIBE3NUyIiEjQlDRERCZqShoiIBE1JQ0REgqakISIiQVPSEBGRoClpiIhI0JrkMiIiDcEYE4ez6OV91SdbGWNaAVtw9mp5zq34ROpCNQ2RyDkUSAW+qnG8e+BrzeMiUU9JQyRyjg58/bLG8e44G+esbNhwROpPSUMkcnrgbItac+2yo4Esa225CzGJ1IuShkjkHM2utQxwahpqmpJGSUlDJAICS293o0ZyMMa0wVlFWUlDGiUlDZHIOBhoCVTWOH4Vzt/d1w0ekUgYaMitSGRs7wS/3BjzM7AJ6AdcHDh+jDHmS2ttsSvRidSRahoikXE0zlyMm4B/AM8DacCfgALgfCUMaYxU0xCJjB7Al4HJezUn8LVwIR6RsFBNQyQyugNL3Q5CJNyUNETCzBhzANAaJQ2JQdojXEREgqaahoiIBE1JQ0REgqakISIiQVPSEBGRoClpiIhI0JQ0REQkaEoaIiISNCUNEREJ2v8Dis1ZlUfFLOUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def posterior(x, xs):\n",
    "    \"\"\" the posterior of a gaussian likelihood with known sigma and mu mean, and flat prior on mu \n",
    "    is a gaussian with mean the average of the n-measurements and sigma sigma/n\n",
    "    \"\"\"\n",
    "    mu = np.mean(xs)\n",
    "    msig = np.sqrt(1./(1.*len(xs)))\n",
    "    return stats.norm.pdf(x, mu, msig)\n",
    "fig, ax = plt.subplots()\n",
    "ax.hist(xs, bins=100, range=(-3., 3.), color='black', label='$x$')\n",
    "mus = np.linspace(-3., 3., 100)\n",
    "pos = np.array([posterior(mu, xs) for mu in mus])\n",
    "dx  = mus[1]-mus[0]\n",
    "print('posterior integral {:5.3f}'.format(np.sum(pos)*dx))\n",
    "ax.plot(mus, pos, color='blue', lw=2, label=\"$p'(\\mu)$\"); ax.legend(fontsize=14)\n",
    "ax.set_xlabel('$\\mu$', fontsize=16); ax.set_ylabel(\"posterior, $p'(\\mu)$\");                        "
   ]
  },
  {
   "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": "fragment"
    }
   },
   "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": 18,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.46153846 0.30769231 0.15384615 0.07692308]\n",
      "[0.63157895 0.28070175 0.07017544 0.01754386]\n",
      "[0.         0.87671233 0.10958904 0.01369863]\n",
      "[0.         0.93772894 0.05860806 0.003663  ]\n"
     ]
    }
   ],
   "source": [
    "def posterior(xi, hpriors, faces = [4, 6, 12, 24]):\n",
    "    \n",
    "    def pxh(x, xface):\n",
    "        if (x <= xface): return 1./(1.*xface)\n",
    "        return 0.\n",
    "\n",
    "    pxhs = np.array([pxh(xi, faces[i])*hpriors[i] for i in range(4)])\n",
    "    px = np.sum(pxhs)\n",
    "    \n",
    "    hpos = pxhs/px\n",
    "    hpos = hpos/np.sum(hpos)\n",
    "    return hpos\n",
    "\n",
    "measurements = [1, 4, 5, 2]\n",
    "hpriors      = [0.25, 0.25, 0.25, 0.25]\n",
    "\n",
    "for meas in measurements:\n",
    "    hpriors = posterior(meas, hpriors)\n",
    "    print (hpriors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Expectation values\n",
    "\n",
    "Given a $x$ rv that follows a pdf, $\\, f(x)$, and a function $g(x)$ on $x$, we define the expected value of $g(x)$ as:\n",
    "\n",
    "$$\n",
    "E[g(x)] \\equiv \\int g(x) \\, f(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 \\, f(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 \\, f(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 $f(x)$.\n",
    "\n",
    "The **median**, is the $x_{med}$ value that divides the distribution in half, \n",
    "\n",
    "$$\n",
    "\\int_0^{x_{med}} f(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 \\, f({\\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": [
    "## 5.  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",
    "f(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",
    "f(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",
    "f(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",
    "f(x|\\tau) = \\frac{1}{\\tau} e^{-x/\\tau}\n",
    "$$\n",
    "\n",
    "* **Gaussian** (or normal). Associated with the distributions of measurements.\n",
    "\n",
    "$$ \n",
    "f(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",
    "f(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",
    "f(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",
    "f(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",
    "f(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": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEPCAYAAACKplkeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X+cVGXd//HXsLv88me6GfJDQwRzES1B0EIlTEJLiTv8gNqtKUp2i19v5S7ptsBQvEFBpOSrAv7AsuijRjcphpD31x+hRkpYqBVwqyzoFqEG8sNddr9/nJltGGZghr1mzgDv5+PBgznnXNd1Pntmdj57nXOd6ySampoQEREJoVXcAYiIyL5DSUVERIJRUhERkWCUVEREJBglFRERCUZJRUREglFSERGRYJRUREQkGCUVEREJpjLuAEpIUweIiOyZRL4F96ekwrp164K3WV1dzfr164O321LlGheUb2zlGheUb2yKq3DlGluuuDp27FhQOzr9JSIiwSipiIhIMEoqIiISjJKKiIgEo6QiIiLB7Fejv0SkuLZv387WrVsBSCTyHoUaXF1dHdu2bYtt/7tSLrGlHtDYtm1bKioqgrWrpCIiQWzfvp0tW7ZwwAEHxJpQACorK4N+UYZUTrE1NTXx4Ycf0q5du2Bt6vSXiASxdevWskgokr9EIsEBBxzQ3LsMQUlFRIJRQtn7hH7PdPpLJMPUqQcVVH7MmI1FimTvooSy9wr53qmnIiIiwainIvucQnsaoN7G/mzYsGEcd9xxTJw4cY+2l1KhsaxZs4ZTTz2VBQsWcNJJJxU5uoiSiogU1Z4k+ZYI/QfCrFmzqKqqCtrmniqnWHKJNamY2WBgOlABzHb3SRnbzwDuBE4ERrj7oxnbDwZeB+a5++jSRC0i+5OPfexjcYfQrJxiySW2aypmVgHMAM4BaoALzawmo9jbwNeBn+Ro5mbgmWLFKCL7h4aGBsaNG0dNTQ01NTXcfPPNNDY2AtEppxtvvLG5bL9+/bjzzjv59re/zXHHHUfv3r25++67d2hv7dq1jBw5kh49etCjRw+uuOKKHR69MXXqVAYOHIi7069fP7p37851113HRx99xIMPPkifPn3o2bMnN910U3Mc2WJ57LHHOPfcc+nRowcnnngio0aN4p133inWYcpLnBfq+wIr3X21u38EzAWGpBdw9zfd/VWgMbOymfUGPgE8VYpgRWTfNW/ePBobG5k/fz6TJk3i4YcfZtasWTnLz5o1i+OPP55f/epXXH311dxyyy387ne/A6IbCi+//HL+9re/4e488sgjvPvuu4wcObL5LnaIrnc89dRTzJkzh1mzZvH4449z+eWXs3z5cn76058yZcoUHnjgAZ588smccdTX1zNmzBgWLVrEnDlz2LBhA1dffXW4A7MH4jz91QlYk7ZcC/TLp6KZtQKmAv8KnBU+NBHZnxxxxBHcfPPNJBIJjj32WFavXs3MmTP5xje+kbX8mWeeyWWXXQZA165due+++3j++efp06cPzz33HK+99hpLliyhS5cuAMyYMYPPfe5zPPfccwwcOBCAxsZG7rjjDg4++GA+9alPMWDAAF588UVefvllWrduTffu3TnllFNYsmQJX/rSl7LGMWLEiObXRx99NJMmTeLMM89k3bp1BT9cK5Q4k0q2gdH5PvL334AF7r7GzHIWMrNRwCgAd6e6urrgIHensrKyKO22VLnGBcWP7UtLbym4TnX1d5vjKrR+dfV3C95focr1/UyPq66ujsrKnb9SWrUq7QmRVAzZYskmkUjQp0+fHS6A9+3bl9tvv50tW7aQSCRo1apVc3uJRIKePXvu0P6RRx7Jhg0bqKysZNWqVXTo0IGuXbs2b+/WrRsdOnRg1apVDBw4kFatWtGpUycOO+yw5jKf+MQn6NatG+3bt29e9/GPf7y53dS+02N59dVXmTJlCitWrOC9995r7gnV1dVx1FFHNU8HU1FRscvj0aZNm2CfsTiTSi3QJW25M5Dv835PA043s38DDgRam9kmdx+bXsjdZwIzk4tNxXiE5972aNByUOzY6uvrC66zfv365rgKrV+K41yu72d6XNu2bcs6p1X6NYFSaGhooLKykoaGhrzKNzU10djYuEP57du3N7eVub2pqYmKioqd2t++fTsNDQ071M3cT+pLv7GxcacYm5qassbd0NCww75TsWzevJnhw4dz+umnM336dKqrq9mwYQNDhw5ly5YtO8SSii2Xbdu20dDQEORxwnEmlaVAdzPrCqwFRgAX5VPR3S9OvTazrwN9MhOKiEi+li1bRlNTU/Od5a+88godOnTgoIMKHw7do0cP3n33XdasWdN8+uutt96irq6O7t27B4t55cqVbNiwgbFjx3LUUUcBsGDBgmDt76nYLtS7ewMwGlhINCzY3X2FmU0ws/MBzOwUM6sFLgDuNbMVccUrIvuuuro6xo8fz8qVK3n88ce55557uPLKK/eordNPP52amhpGjx7Nq6++yvLlyxk9ejS9evWif//+wWLu1KkTbdq04YEHHuCtt95i8eLF3H777cHa31Ox3qfi7guABRnrxqW9Xkp0WmxXbTwIPFiE8ERkPzF06FC2b9/OeeedRyKRYMSIEXucVBKJBPfffz/f+973GDZsGBAlmtRAgFAOP/xwpk2bxuTJk5kzZw7HH38848eP5+KLL9595SJKpA9x28c1pY8TD2VvONddbood2x+++sOC6/R67JrmuAqt3+uxawreX6HK9f1Mj2vz5s07XGSOUyHXVEqtHGPbvHkzRx111K6uqeSdDTWhpIiIBKOkIiIiwSipiIhIMEoqIiISjJKKiIgEo6QiIiLBKKmIiEgwSioiIhKMkoqIiASjZ9SLBFTo89hDP09dSm/JkiVccMEF/OEPf9hhKvv9lZKKiBTVQVOnlnR/G8eMKen++vTpw7Jly/aK58eXgpKKiEgLtG7dmiOOOCLuMMqGrqmIyH5t2LBh3HDDDYwbN46amhpqamq4+eabmx8u9v7773PttddSU1NDt27dGD58OH/605+a6y9ZsoROnTqxYcMGAP7xj39wzTXXcOKJJ3LMMcdw2mmn7fC8+4ceeoj+/ftzzDHH0KtXLy666KLmCSYbGxuZNm0affr0oWvXrpx11lksXLiwue6aNWvo1KkTTzzxBCNGjKBbt24MGDCAZ599thSHKi9KKiKy35s3bx6NjY3Mnz+fSZMm8fDDDzcnguuuu45ly5Zx//338/jjj9OuXTsuvvhitmzZkrWt2267jTfeeIM5c+bwzDPPMHXqVDp06ADA8uXLGTt2LNdffz3PPvssc+fOZcCAAc11Z8+ezT333MONN97I4sWLGTx4MFdccQV//OMfd9jH5MmTufzyy1m0aBEnnXQS3/zmN/nwww+Lc3AKpNNfIrLfO+KII5qfd3LssceyevVqZs6cydlnn81TTz3FY489xqmnngrAD37wA/r27cu8efO46KKdH1ZbW1vLCSecwGc+8xmA5qc/Aqxdu5b27dszaNAgDjzwQDp37kzPnj2bt997771cddVVDB06FIBvfetbvPTSS9x777388If/fCTDlVdeyaBBgwAYO3Ysjz76KCtWrKBv377hD06B1FMRkf3eySefvMMDtHr37s27777LypUradWqFb17927edvDBB/OpT32KP//5z1nbuuSSS/jlL3/JF77wBSZMmMALL7zQvO2MM86gc+fOnHrqqYwePRp3Z9OmTQBs3LiRd999lz59+uzQ3imnnLLTvo4//vjm16leULk8b0c9FZGABi25tbAKY4r/gC/Zc7t6iGGupzgOHDiQ3/72tzz99NM8//zzXHLJJXz5y19m2rRpHHjggSxevJjnn3+e5557jrvuuovJkyfzxBNPcMABB+RsN3NdVVXVTttS14DipqQiZUn3e0gpLVu2jKampuYv6FdeeYUOHTrQvXt3Ghsbefnll5tPf23cuJE33niD4cOH52zvsMMOY9iwYQwbNozPf/7zXH311UyaNIk2bdpQWVlJ//796d+/P2PGjOGkk05i8eLFfO1rX6NDhw4sXbp0h2fZL126lO7duxf3AAQUa1Ixs8HAdKACmO3ukzK2nwHcCZwIjHD3R5PrPw3cDRwMbAcmuvvPShm7iOw76urqGD9+PJdccglvvPEG99xzD9deey3HHHMMX/ziF7nhhhu47bbbOPjgg5k8eTIHHnggX/nKV7K2dfvtt9OrVy969OjB9u3befLJJzn66KNp06YNixYtYs2aNfTt25dDDz2U3/zmN2zatKk5aVx11VVMmTKFrl270qtXL37+85/z0ksv8eSTT5bycLRIbEnFzCqAGcDZQC2w1Mzmu/tracXeBr4O/EdG9c3AJe7+FzPrCLxsZgvd/f0ShC4loNNIUkpDhw5l+/btnHfeeSQSCUaMGMGVV14JwB133MH48eO57LLL2LZtG3369OHhhx+mXbt2Wdtq3bo1kydP5u2336Zt27acfPLJPPjggwAccsghzJw5k6lTp7J161aOPvpopkyZQr9+/QAYOXIkmzZt4pZbbmH9+vV069aNWbNmccIJJ5TkOIQQZ0+lL7DS3VcDmNlcYAjQnFTc/c3kth1OFrr7n9NerzOzvwIfB5RURMpMqe9w3xOVlZVMnDiRiRMn7rTt0EMPZfr06Tnrfvazn2Xt2rXNy9deey3XXntt1rKpUWOp+1IytWrViuuuu47rrrsu6/YuXbrssK+UbOviEufor07AmrTl2uS6gphZX6A1sCpQXCIisofi7KlkGzqRe6hFFmZ2JPAj4FJ332nog5mNAkYBuDvV1dV7EucuVVZWFqXdlirXuCC/2NJHt+Qjvb1C66bqp+Iq5b7zfY/K9f1Mj6uuro7KyvIZ+5NvLIlEglatWpU09nI6TkDzAIIQn7E4f7JaoEvacmdgXb6Vzexg4Angu+7+YrYy7j4TmJlcbCrGOO7q6uqyGR+erlzjgvxiq6+vL6jN9PYKrZuqn4qrlPvO9z0q1/czPa5t27ZRUVERc0SRysrKnKeYMj3yyCMAeZdvqUJiK5Vt27bR0NCQ9TPWsWPHgtqKM6ksBbqbWVdgLTAC2Pn21CzMrDUwD3jI3R8pXogiIlKI2K6puHsDMBpYCLwerfIVZjbBzM4HMLNTzKwWuAC418xWJKsbcAbwdTP7ffLfp2P4MUREJE2sJ/bcfQGwIGPduLTXS4lOi2XW+zHw46IHKCJ529Xd51LeQr53mvtLRIJRYtn7hH7PlFREJIi2bdvy4YcfKrHsRZqamvjwww9p27ZtsDbLa1ybiOy1KioqaNeuHZs3bwZyT7hYCm3atGHbtm2x7X9XyiW2VPJv165d0FF7SioiEkxFRUXzbLtxKtch2FDesYWg018iIhKMkoqIiASjpCIiIsEoqYiISDBKKiIiEoySioiIBKOkIiIiwSipiIhIMEoqIiISjJKKiIgEo6QiIiLBKKmIiEgwSioiIhKMkoqIiASjpCIiIsEoqYiISDCxPqTLzAYD04EKYLa7T8rYfgZwJ3AiMMLdH03bdinw3eTiLe4+pzRRixTH1KkH5dzWvn0FmzfvvH3MmI3FDEmkYLH1VMysApgBnAPUABeaWU1GsbeBrwM/yah7GDAe6Af0Bcab2ceKHbOIiOxanKe/+gIr3X21u38EzAWGpBdw9zfd/VWgMaPuF4FF7r7B3d8DFgGDSxG0iIjkFufpr07AmrTlWqKex57W7ZRZyMxGAaMA3J3q6uo9i3QXKisri9JuS5VrXJBfbFVVVQW1md5eoXVT9VNxlXLf6XXbt6/IWa5Vq1a0b98+S/02Be0vtHL9nJVrXFC+sYWKK86kksiyrilkXXefCcxMbV+/fn2ezeevurqaYrTbUuUaF+QXW319fUFtprdXaN1U/VRcpdx3et1s10xS2rdvz+bNm7PUj/eaSrl+zso1Lijf2HLF1bFjx4LaiTOp1AJd0pY7A+sKqDsgo+7/CxKVSEwGLbk157aqqqrsCWvMNUWMSKRwcSaVpUB3M+sKrAVGABflWXchcGvaxflBwHfChygiIoWI7UK9uzcAo4kSxOvRKl9hZhPM7HwAMzvFzGqBC4B7zWxFsu4G4GaixLQUmJBcJyIiMYr1PhV3XwAsyFg3Lu31UqJTW9nq3g/cX9QARUSkILqjXkREglFSERGRYJRUREQkGCUVEREJRklFRESCUVIREZFgYh1SLPu2g6ZOzbq+on17Dsoy5cjGMWOKHZKIFJl6KiIiEoySioiIBKPTX1I0S5a0zrq+qipBff3O23rp7JfIXk89FRERCUZJRUREglFSERGRYJRUREQkGCUVEREJRklFRESCUVIREZFgCrpPxcxOBQYDpwIdgXbAeuBPwDPAL9z9vdBBiojI3iGvpGJmlwL/AfQE/gG8CvwF2AIcBvQD/hWYYWYOfN/d/7coEYuISNnabVIxs+XAEcBDwCXA7929KUu5Q4AvAxcDK8zsMnf/2W7aHgxMByqA2e4+KWN7m+R+ewN/B4a7+5tmVgXMBk5O/gwPuft/7e5nERGR4srnmsoDQFd3v8Hdl2VLKADu/oG7P+zu5wKnAe/vqlEzqwBmAOcANcCFZlaTUWwk8J67HwtMAyYn118AtHH3XkQJ5xtm9sk8fhYRESmi3fZU3P3OQht19+XA8t0U6wusdPfVAGY2FxgCvJZWZghwU/L1o8BdZpYAmoADzKyS6LrOR0Sn5UREJEYFjf4ys5MD7rsTsCZtuTa5LmsZd28APgAOJ0owHwLvAG8DU9x9Q8DYRERkDxQ6S/H/mNlX3P1/Auw7kWVd5qm1XGX6AtuJRqB9DHjOzBanej0pZjYKGAXg7lRXV7c46EyVlZVFabelyiGuqqqqrOsTiUTWbenx5qqbS0vqpuqnjlkp951v3XyOWRzK4XOWTbnGBeUbW6i4Ck0qPwEWmNnX3P2x9A1m1h+Y5O7982yrFuiSttwZWJejTG3yVNchwAbgIuBX7l4P/NXMfgP0AXZIKu4+E5iZXGxav359nqHlr7q6mmK021LlEFd9fX3W9VVVVVm3pcebq24uLambqp86ZqXcd7518zlmcSiHz1k25RoXlG9sueLq2LFjQe0UdPrL3b8J/Bcw18yuAjCzXmb2S+BZol5DvpYC3c2sq5m1BkYA8zPKzAcuTb4eBjydHCjwNjDQzBJmdgDRfTNvFPKziIhIeAXfUe/uE4CrgB+Y2TPAMuAE4HKgVwHtNACjgYXA69EqX2FmE8zs/GSx+4DDzWwlcD0wNrl+BnAg8Eei5PSAu79a6M8iIiJhFfzkRzM7DOhBdE3jdGAJMCCZJAri7guABRnrxqW93ko0fDiz3qZs60VEJF6Fjv4aT3Td4mpgKlHvpA9wR/jQRERkb1NoT+VGojvZv+/udQBm9jYwz8w+AXwtefFcRET2Q4UmlePdfVX6Cnd/2sw+T3Qa61fAWaGCExGRvUuho79W5Vj/CtAf+GSAmEREZC8V7Hkq7r4S+Gyo9kREZO+TzyzF/w3c5O7LdlfW3evMrC3wb8Bmd78nQIwikoepUw8qqPyYMRuLFInsz/K5pvIW8KKZ/Z7ojvrngFfThxCbWUeiqVPOA/4FWEs0MkxESmTQklsLqzDmmuIEIvu1fJLKR8DngQuB8URTpTSZ2T+AbUR30VcRzdP1W+DfgR+5e2NRIhYRkbKVT1L5d6K73a8xs41EI7xOA44E2hI9POsN4Fl3f6tokYqISNnLJ6ls4J9zet1A9Bz6ybsoLyIi+6l8ksrzwBQz+zjRKa6sT34UERHJJ6mMBuYk/zUBi83sVaKJJFP/VuhOehER2e19Ku6+zt3PJnoKYwL4GdETFwcTTdnyMrDRzF4xs/uKGayIiJS3vG9+dPd3gXnANHc3d+9BNBLsTOBbwO+BkI8bFhGRvUxBc3+5+1czljcRXXN5PmRQIiKydwo2TYuIiIiSioiIBKOkIiIiwSipiIhIMEoqIiISTKFPfgzKzAYD04EKYLa7T8rY3gZ4COhNNMfYcHd/M7ntROBe4GCgETjF3beWLnoREckUW0/FzCqAGcA5QA1woZnVZBQbCbzn7scC04DJybqVwI+Bq9y9JzAA0B39IiIxi7On0hdY6e6rAcxsLjAEeC2tzBDgpuTrR4G7zCwBDCJ6pstyAHf/e6mCFhGR3OJMKp2ANWnLtUC/XGXcvcHMPgAOB3oQPdNlIfBxYK6735a5AzMbBYxK1qe6ujr4D1FZWVmUdluqHOKqqqrKuj6RSGTdlh5vrrq5tKRuqn7qmJVy3/nWLfYx21Pl8DnLplzjgvKNLVRccSaVRJZ1mTMg5ypTCfQHTgE2A782s5fd/dfpBd19JjAzVW/9+vUtiziL6upqitFuS5VDXPX12c9IVlVVZd2WHm+uurm0pG6qfuqYlXLf+dYt9jHbU+XwOcumXOOC8o0tV1wdO3YsqJ04R3/VAl3SljsD63KVSV5HOYTo+S61wDPuvt7dNwML0LxjIiKxi7OnshTobmZdiZ5pPwK4KKPMfOBS4AVgGPC0u6dOe33bzNoTPe74TKIL+SIiEqPYkkryGsloYCHRkOL73X2FmU0Afufu84H7gB+Z2UqiHsqIZN33zOwOosTUBCxw9ydi+UH2YQdNnVpQ+Y1jxhQpEhHZW8R6n4q7LyA6dZW+blza663ABTnq/phoWLGIiJQJ3VEvIiLBKKmIiEgwSioiIhJMrNdUpLwtWdK6oPK9dJ1eZL+nnoqIiASjpCIiIsEoqYiISDBKKiIiEoySioiIBKOkIiIiwSipiIhIMLpPRUQ0eagEo56KiIgEo56KiGj2BAlGPRUREQlGSUVERIJRUhERkWCUVEREJBglFRERCSbW0V9mNhiYDlQAs919Usb2NsBDQG/g78Bwd38zbftRwGvATe4+pVRxi4hIdrH1VMysApgBnAPUABeaWU1GsZHAe+5+LDANmJyxfRrwZLFjFRGR/MR5+qsvsNLdV7v7R8BcYEhGmSHAnOTrR4GzzCwBYGZfAVYDK0oUr4iI7EacSaUTsCZtuTa5LmsZd28APgAON7MDgBuA75cgThERyVOc11QSWdY15Vnm+8A0d99kZjl3YGajgFEA7k51dfUehppbZWVlUdptqRBxVVVVFVQ+c3+56icSiazb0uu3ZN+F1k3VTx2zUu4737rldsxSdfflz3+xlGtsoeKKM6nUAl3SljsD63KUqTWzSuAQYAPQDxhmZrcBhwKNZrbV3e9Kr+zuM4GZycWm9evXB/8hqqurKUa7LRUirvr6+oLKZ+4vV/2qqqqs29Lrt2TfhdZN1U8ds1LuO9+65XbMUnX35c9/sZRrbLni6tixY0HtxJlUlgLdzawrsBYYAVyUUWY+cCnwAjAMeNrdm4DTUwXM7CZgU2ZCERGR0ovtmkryGsloYCHwerTKV5jZBDM7P1nsPqJrKCuB64Gx8UQrIiL5iPU+FXdfACzIWDcu7fVW4ILdtHFTUYITEZGC6Y56EREJRklFRESCUVIREZFglFRERCQYJRUREQlGSUVERIJRUhERkWCUVEREJBglFRERCUZJRUREgol1mhYR2TccNHVqQeU3jhlTpEgkbkoq+zj9sotIKSmpiEiLLVnSuqDyvfS3yz5L11RERCQYJRUREQlGSUVERILRNZV9nM51i0gpqaciIiLBKKmIiEgwSioiIhJMrNdUzGwwMB2oAGa7+6SM7W2Ah4DewN+B4e7+ppmdDUwCWgMfAd9y96dLGryIiOwktp6KmVUAM4BzgBrgQjOrySg2EnjP3Y8FpgGTk+vXA+e5ey/gUuBHpYlaRER2Jc7TX32Ble6+2t0/AuYCQzLKDAHmJF8/CpxlZgl3X+bu65LrVwBtk70aERGJUZynvzoBa9KWa4F+ucq4e4OZfQAcTtRTSfkqsMzdt2XuwMxGAaOS9amurg4XfVJlZWVR2m2pVFxVVVUF1Uv/WVpSd1f1E4lE1m2h9l1o3VT9uI/ZruqW2zFL1Q1xzIqhXH8voXxjCxVXnEklkWVdUyFlzKwn0SmxQdl24O4zgZmpeuvXr89WrEWqq6spRrstlYqrvr6+oHrpP0tL6u6qflVVVdZtofZdaN1U/biP2a7qltsxS9UNccyKoVx/L6F8Y8sVV8eOHQtqJ87TX7VAl7TlzsC6XGXMrBI4BNiQXO4MzAMucfdVRY9WRER2K86eylKgu5l1BdYCI4CLMsrMJ7oQ/wIwDHja3ZvM7FDgCeA77v6bEsYsIiK7EFtPxd0bgNHAQuD1aJWvMLMJZnZ+sth9wOFmthK4HhibXD8aOBb4npn9PvnviBL/CCIikiHW+1TcfQGwIGPduLTXW4ELstS7Bbil6AGKSNHpQXL7Ft1RLyIiwSipiIhIMJr6XkRipccz7FvUUxERkWCUVEREJBglFRERCUZJRUREglFSERGRYDT6q8wVemMY6OYwEYmPeioiIhKMeiplrtAx/KBx/LL/yNWTr2jfnoM2b866TT354lJPRUREglFPRUT2Wrl68lVVCerrs29TT7641FMREZFglFRERCQYJRUREQlG11REZL+lB4SFp56KiIgEo56KiOy39CyX8GJNKmY2GJgOVACz3X1SxvY2wENAb+DvwHB3fzO57TvASGA78H/cfWEJQy+Iutgisr+ILamYWQUwAzgbqAWWmtl8d38trdhI4D13P9bMRgCTgeFmVgOMAHoCHYHFZtbD3beX9qfIj/4aEtn36I/F7OLsqfQFVrr7agAzmwsMAdKTyhDgpuTrR4G7zCyRXD/X3bcB/2tmK5PtvVCi2EVkP9eSPxb35Yli40wqnYA1acu1QL9cZdy9wcw+AA5Prn8xo26n4oUKf/jqD7Our6qqor6+fqf1vR67ppjhiMherKVz+uX6PspZt4TfR4mmpqaS7SydmV0AfNHdr0gu/yvQ192vSSuzIlmmNrm8iqhHMgF4wd1/nFx/H7DA3R/L2McoYBSAu/cu/k8lIrJPSuRbMM4hxbVAl7TlzsC6XGXMrBI4BNiQZ13cfaa793H3PkQHJfg/M3u5WG3vi3GVc2zlGlc5x6a49p3YdhNX3uI8/bUU6G5mXYGb7gRrAAAGtUlEQVS1RBfeL8ooMx+4lOhayTDgaXdvMrP5wE/M7A6iC/Xdgd+WLHIREckqtp6KuzcAo4GFwOvRKl9hZhPM7PxksfuAw5MX4q8HxibrrgCc6KL+r4Cry3Xkl4jI/iTW+1TcfQGwIGPduLTXW4ELctSdCEwsaoD5mRl3ADmUa1xQvrGVa1xQvrEprsKVa2xB4ortQr2IiOx7NPeXiIgEo7m/8tSSKWWKGFOX5D47AI3ATHefnlFmAPDfwP8mV/3c3ScUM660fb8JbCSaSqchOQovfXuC6JieC2wGvu7urxQ5puOAn6WtOgYY5+53ppUZQImOmZndD3wZ+Ku7n5Bcd1gyxk8CbwLm7u9lqXsp8N3k4i3uPqfIcd0OnAd8BKwCLnP397PUfZNdvO9FiOsm4Ergb8li/5k8tZ5Zd5e/w0WK7WfAcckihwLvu/uns9R9k+Ids6zfE8X6nKmnkoe0KWXOAWqAC5NTxaRrnlIGmEY0pUyxNQBj3P144FTg6ixxATzn7p9O/itJQknz+eR+s/2SnEM0cq870f1Edxc7GHf/U+pYEP0BsBmYl6VoqY7Zg8DgjHVjgV+7e3fg18nlHSS/EMYT3TDcFxhvZh8rclyLgBPc/UTgz8B3dlF/V+976LgApqW9X9kSSj6/w8Fjc/fhaZ+3x4Cf76J+sY5Zru+JonzOlFTy0zyljLt/BKSmlEk3BEhl8EeBs5J/iReNu7+T+sve3TcSjaIr6swCgQ0BHnL3Jnd/ETjUzI4s4f7PAla5+1sl3OcO3P1Zonuv0qV/luYAX8lS9YvAInffkPzrchHZv2yDxeXuTyVHbUI0o0XnUPtrSVx5yud3uGixJb8LDPhpyH3mYxffE0X5nCmp5CfblDKZX947TCkDpKaUKQkz+yTwGeClLJtPM7PlZvakmfUsVUxAE/CUmb2cnN0gUz7HtZhGkPuXPK5jBvAJd38Hoi8E4IgsZeI+dpcDT+bYtrv3vRhGm9mrZnZ/jr+k4z5epwN17v6XHNtLcswyvieK8jlTUslPth5H5rC5fMoUhZkdSNS1/nd3/0fG5leAo939JOCHwC9KEVPS59z9ZKJTDleb2RkZ2+M8Zq2B84FHsmyO85jlK85jdyPRKZWHcxTZ3fse2t1AN+DTwDtAttkaYzteSRey615K0Y/Zbr4ncin4uCmp5KclU8oUlZlVEX1QHnb3nc7Xuvs/3H1T8vUCoMrMqosdV3J/65L//5XoukXfjCJ5TbdTJOcAr7h7XeaGOI9ZUl3qNGDy/79mKRPLsUtetP0ycLG7Z/1yyeN9D8rd69x9u7s3ArNy7C+2z1ry++Bf2HGAyA6KfcxyfE8U5XOm0V/52eMpZYoZVPI87X3A6+5+R44yHYi63U1m1pfoD4m/FzOu5H4PAFq5+8bk60FEE4Gmm0902mIu0YXAD1Ld8RLI+ZdjXMcsTeqzNCn5/39nKbMQuDXtVM8gdn3hvMWSo6duAM509805yuTzvoeO68i0z81Q4I9ZiuXzO1wsXwDe8OTEuJmKfcx28T1RlM+Zbn7Mk5mdC9xJNBzxfnefaGYTgN+5+3wzawv8iOh85QZghCefFVPEmPoDzwF/IBoqCPCfwFEA7n6PmY0Gvkl0umILcL27LylmXMnYjuGfo6oqgZ8kj9lVabElgLuILvxtJhqi+rsSxNae6DzxMe7+QXJdelwlO2Zm9lNgAFAN1BGNtPkF0TRERwFvAxe4+wYz6wNc5f+c2ftyovcbYKK7P1DkuL4DtOGfCfZFd7/KzDoSDdE9N9f7XuS4BhCd+moiGhr7DXd/Jz2uZN2dfodDxZUrNne/z8weJDpW96SVLeUxy/U98RJF+JwpqYiISDC6piIiIsEoqYiISDBKKiIiEoySioiIBKOkIiIiwSipiIhIMEoqIiISjJKKiIgEo2laRGJgZscSTUF+q7uPT1t/N/A1omdrFH12AZHQ1FMRiYG7rwRmA9elJqs0s3FEU8oPVUKRvZWmaRGJSXLiylXA/wXeAGYCF7q7xxqYSAsoqYjEyMwmAmOITkVf6+4zYg5JpEV0+kskXn8hmvn3BSUU2RcoqYjExMwGAvcSPYPnc2Z2UswhibSYkopIDMzsZKLnpswmegbH28CtccYkEoKSikiJJYcTPwk8BVzj7h8B3wfOLcHz3EWKSklFpISSI76eIrpH5eLkc9UBHiIaATYprthEQtDoLxERCUY9FRERCUZJRUREglFSERGRYJRUREQkGCUVEREJRklFRESCUVIREZFglFRERCQYJRUREQnm/wMjvKtsUzj/YQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "N, p = 100, 0.08\n",
    "\n",
    "ns = np.arange(20)\n",
    "fig, ax = plt.subplots()\n",
    "ax.bar(ns, stats.binom.pmf(ns, N, p), color='blue', alpha=0.5, label='binomial')\n",
    "ax.bar(ns, stats.poisson.pmf(ns, N*p), color='red', alpha=0.5, label='poisson')\n",
    "ax.set_xlabel('$x$', fontsize=16)\n",
    "ax.set_ylabel('$f(x$)', fontsize=16)\n",
    "ax.legend(fontsize=14)\n",
    "fig.savefig('figs/gen_binomial_poisson.png');"
   ]
  },
  {
   "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$ samples, each one with $n$ measurements, $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$ samples, each with $n$ measurements, $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$.\n",
    "\n",
    "**Exercise:** Get the possion 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": [
    "## 6. Non parametric pdfs\n",
    "\n",
    "In some cases we do not know the parametric pdf that describes our sample, but we can try to infer a **non-parametrical pdf** using the data itself.\n",
    "\n",
    "One possibility is to use a kernel function around each point of the sample. Given a kernel function, $k(u)$, we can construct the pdf as:\n",
    "$$\n",
    "f(x) = \\frac{1}{n} \\sum_{i=1}^{n} k(u_i)\n",
    "$$\n",
    "where $u = |x-x_i|/w$, and $w$ the width of the kernel. The kernel is also a symmetric function on $u$. There are many kernel functions that we can use: \n",
    "$$\n",
    "k(u) = (1-|u|), \\;\\;\\; k(u) = \\frac{3}{4} (1-u^2), \\;\\;\\; \\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{u^2}{2}}\n",
    "$$\n",
    "for the two first kernels $|u|\\le 1$. They correspond to a triangular, a Epanechnikov and gaussian shapes, respectively. There is also the possibility to use \"adaptative\" kernels, which width depend on the distance to the nearest neighbors. There is a large industry about the kernel functions to use!\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Example:** Generate n (small) samples with two gaussian with sigma 1 and means 0, 5.\n",
    "Plot the normalized histogram. Use different kernels to obtain a non-parametric pdf."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAENCAYAAAAYIIIKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFqZJREFUeJzt3X+0ZWV52PHvzL0UcwIocmrTO0DEShLxR6EC2lLBBDSTRmAl1YeRRYoRO8s2SFPGZWRpsCVqETJW1gomzBppyQpKHzG1kxYltGDMaiQOFiMLEAtI4HJRuaJI10Vw4PaPvWc4nLl3Zr937tnnzDnfz1qzuPvd73v3865zmGfe/b773WsWFxeRJKmptcMOQJK0fzFxSJKKmDgkSUVMHJKkIiYOSVIRE4ckqYiJQ5JUxMQhSSpi4pAkFZkedgAD4uPwkrQya/ZWYVwTB3NzcwP9/d1ul/n5+YFeY9TY58kxif22zzAzM9OonbeqJElFTBySpCImDklSEROHJKmIiUOSVMTEIUkqYuKQJBUxcUiSipg4JElFxvbJcUmjY/Pmg4dy3U2bnhjKdcedIw5JUhEThySpiIlDklTExCFJKmLikCQVMXFIkoqYOCRJRUwckqQiJg5JUhEThySpiIlDklTExCFJKmLikCQVMXFIkoqYOCRJRUwckqQiJg5JUhEThySpSGuvjo2I9cAVwBSwNTMv7Tt/IfAuYAfwKPDOzPzb+twzwB111Qcz84y24pYkPV8riSMipoArgTcBs8D2iNiWmXf1VLsdOD4zFyLiXwGXAWfV557MzGPbiFWStGdtjThOBO7NzPsBIuI64ExgV+LIzFt66t8KnNNSbJKkAm0ljnXAQz3Hs8Dr9lD/POALPccviIjbqG5jXZqZn1/9ECVJTbSVONYsUba4VMWIOAc4Hjilp/jIzJyLiJcBN0fEHZl5X1+7jcBGgMyk2+2uTuTLmJ6eHvg1Ro19nhyr3e9OZ2rVfleJbvfAxnUn8bNeaZ/bShyzwBE9x4cDc/2VIuI04APAKZn51M7yzJyr/3t/RHwJOA54XuLIzC3AlvpwcX5+fjXj302322XQ1xg19nlyrHa/FxYOXrXfVWJ+/onGdSfxs+7v88zMTKN2bSWO7cDREXEU8DCwATi7t0JEHAdcBazPzO/1lB8KLGTmUxHRBU6imjiXJA1BK89xZOYO4HzgRuDuqijvjIhLImLn0trLgYOAz0bE1yNiW13+CuC2iPgb4BaqOY67kCQNxZrFxSWnGvZ3i3Nzu90JW1UOayfDJPYZVr/fmzcP51bVpk3eqtqTZW5VLTUn/Tw+OS5JKmLikCQVMXFIkoqYOCRJRUwckqQiJg5JUhEThySpiIlDklTExCFJKmLikCQVMXFIkoqYOCRJRUwckqQiJg5JUhEThySpiIlDklTExCFJKmLikCQVMXFIkopMDzsASfuHYb03XKPHEYckqYiJQ5JUxMQhSSpi4pAkFTFxSJKKuKpK0tgqWQnW6UyxsPBc/U2bnhhESGPBEYckqYiJQ5JUpLVbVRGxHrgCmAK2ZualfecvBN4F7AAeBd6ZmX9bnzsX+GBd9cOZeU1bcUuSnq+VEUdETAFXAr8CHAO8PSKO6at2O3B8Zr4GuB64rG77YuBDwOuAE4EPRcShbcQtSdpdWyOOE4F7M/N+gIi4DjgTuGtnhcy8paf+rcA59c+/DNyUmY/VbW8C1gOfaSFuSVKftuY41gEP9RzP1mXLOQ/4wgrbSpIGqK0Rx5olyhaXqhgR5wDHA6eUtI2IjcBGgMyk2+2uLNKGpqenB36NUWOfJ8dS/e50poYUTTvWrl1Lp9PZddztHjjEaNqx0u93W4ljFjii5/hwYK6/UkScBnwAOCUzn+pp+8a+tl/qb5uZW4At9eHi/Pz8Pge9J91ul0FfY9TY58mxVL97n3EYR51Oh4WFhV3H8/Pj/xxH/+c8MzPTqF1biWM7cHREHAU8DGwAzu6tEBHHAVcB6zPzez2nbgQ+2jMh/mbgosGHLElaSitzHJm5AzifKgncXRXlnRFxSUScUVe7HDgI+GxEfD0ittVtHwN+jyr5bAcu2TlRLklq35rFxSWnGvZ3i3Nzu90JW1WTeAvDPk+Opfo97i9y6r9VNQlbjixzq2qpeeXn8clxSVIRE4ckqYiJQ5JUpChxRMTaiPj7gwpGkjT6Gi3HjYgXAZ8E3gr8BPjpejXUiZn5wT02liSNlaYjjj8CHgd+Fni6LvsKcNYggpIkja6mieNU4ILMfIR6u4/MfBR4yaACkySNpqaJ43HgeRuaRMSRwCOrHpEkaaQ1TRxbgc9FxC8CayPiHwPXUN3CkiRNkKZ7VX0M+DHVy5gOAK6m2lfqigHFJUkaUU0Tx9/LzE8An+gtjIifAb6z6lFJkkZW01tV31qm/K5lyiVJY6pp4tht06uIOAR4dnXDkSSNuj3eqoqIh6iW3/5URDzYd/owfO+3JE2cvc1xnEM12rgB+I2e8kXgu5l5z6ACkySNpj0mjsz8C4CI6Gbmwp7qSpImQ6NVVZm5EBHHAm+gehBwTc+5iwcUmyRpBDWaHI+IjcD/Bn4J+B3g1cAm4OWDC02SNIqarqp6H7A+M38NeLL+786dciVJE6Rp4nhJZv5l/fOzEbE2M78AnD6guCRJI6pp4piNiJfWP38LODMi3sBzW6xLkiZE0y1HLgNeATwAXAJcD/wd4ILBhCVJGlV7HXFExBrgy8BNAPUtqkOBQzPzDwcbniRp1Ow1cWTmInAHPduLZObTmfn/BhmYJGk0NZ3juB34uUEGIknaPzSd4/gS8MWI+M/Azv2rAMjMq1c/LEnSqGqaOE4Cvg2c0le+SPVSJ0nShGi65cgvDjoQSdL+oemIY59FxHqqV81OAVsz89K+8ydTvWHwNcCGzLy+59wzVBP0AA9m5hntRC1J6tdK4oiIKar3lb8JmAW2R8S2zOx9g+CDwDuA9y7xK57MzGMHHqgkaa/aGnGcCNybmfcDRMR1wJn0vHo2Mx+oz/lWQUkaYW0ljnVUq7F2mgVeV9D+BRFxG7ADuDQzP99fod7BdyNAZtLtdvch3L2bnp4e+DVGjX2eHEv1u9OZGlI07Vi7di2dTmfXcbd74BCjacdKv9+NEkdEPJaZL16i/HuZ+ZIGv2K3d5bTs6S3gSMzcy4iXgbcHBF3ZOZ9vRUycwuwZefvnp+fL/j15brdLoO+xqixz5NjqX4vLBw8pGja0el0WFh47n118/NPDDGadvR/zjMzM43aNX0A8ID+gog4gGqiu4lZ4Iie48OBuYZtycy5+r/3Uz1TclzTtpKk1bXHEUdE/CXVyOAFEfHlvtOHA3/V8DrbgaMj4ijgYWADcHaThhFxKLCQmU9FRJfqmZLLGl5XkrTK9naraivVbaYTgE/1lC8C3wVubnKRzNwREecDN1KNUq7OzDsj4hLgtszcFhEnAP+VagPF0yPi32fmK6l25b2qnjRfSzXHcdcyl5IkDdiaxcW9TzVExC9k5jdbiGe1LM7NNb4TtiKTeO/bPk+Opfq9efNkzXFs2jSxcxxLzUk/T9Mnx78ZEW8GjgUO6jt3cVGkkqT9WtNVVX8ABHALsNBzqmRllCRpDDR9juPtwLGZ+dBea0qSxlrT5bjfB344yEAkSfuHpiOOzcC1EfEfqFZT7bJzGxFJ0mRomjh2vlv8LX3lizR/CFCSNAaarqpqektLkjTmijY5jIgjgHWZeeuA4tGQNFmj3+lMLblf0SSsd5f0nKbLcY8EPkP1HMcicFBEvBVYn5nvGmB8kqQR0/QW1FXA/wAOBn5Sl91E9WImSdIEaXqr6kTgVzPz2YhYBMjMxyPihYMLTdJqKtkyZLnbkhI0H3F8F3h5b0FEHEP1uldJ0gRpmjh+H/jvEfGbwHREvB34L8DHBhaZJGkkNUocmXk18D7gbVSvgP0XwO9m5rUDjE2SNIKarqqaqt/zvdu7viVJk6XprarvRMQnI+KkgUYjSRp5TVdVvZlqh9zP1G/i+wzw6cy8Y2CRSdrNuL9MSfuHpnMct2fm+zLzSOBcqte7/q+I+MZAo5MkjZyV7EF1D3A31ST5S1c1GknSyGs6Of4i4J8DZwOvB/6cainutsGFJkkaRU3nOOaAvwI+Dfx6Zj4+uJAkSaOsaeL4B5n5yEAjkSTtF5q+j+ORiHgTsAF4SWaeHhHHA4dk5s0DjVCSNFIaTY5HxHuo3gL4f4GT6+IngQ8PKC5J0ohquqrqt4HTMvNS4Nm67JvAzw8kKknSyGqaOA6mWn4L1YucAA4Anl71iCRJI61p4vgy8P6+sguAW1Y3HEnSqGu6quo9wJ9FxL8EDo6Ie4AfAacPLDJJ0kgqWVV1AnAC8LNUt62+mpnP7rnlcyJiPXAFMAVsredLes+fDHwCeA2wITOv7zl3LvDB+vDDmXlN0+tKklZX0xEHmbkIfLX+UyQipoArqd5RPgtsj4htmXlXT7UHgXcA7+1r+2LgQ8DxVPMrX6vb/qA0DknSvlvJXlUrcSJwb2ben5lPA9cBZ/ZWyMwHMvMbPLdqa6dfBm7KzMfqZHETsL6NoCVJu2s84thH63huVRZUo47X7UPbdf2VImIjsBEgM+l2uyuLtKHp6emBX6NNnc7UXuusXbuWTqezW3m3e+AgQhoJo/Y5N/mcVsNyn/U46+/zOH+vd1rp97utxLFmibLFJcpW3DYztwBbdp6fn59v+OtXptvtMuhrtGlhYe/veeh0OiwsLOxWPj//xCBCGgmj9jk3+ZxWw3Kf9Tjr7/M4f6936v9+z8zMNGrX1q2qWeCInuPDqTZOHHRbSdIqa2vEsR04OiKOAh6m2vPq7IZtbwQ+GhGH1sdvBi5a/RAlSU20MuLIzB3A+VRJ4O6qKO+MiEsi4gyAiDghImaBtwFXRcSdddvHgN+jSj7bgUvqMknSELQ14iAzbwBu6Cu7uOfn7VS3oZZqezVw9UADlCQ10tYchyRpTJg4JElFTBySpCKtzXFofG3e3M6zBf02bRr/dfbSKHLEIUkqYuKQJBUxcUiSipg4JElFTBySpCKuqhoTw1rZJI2rff1/apxX/TnikCQVMXFIkoqYOCRJRUwckqQiJg5JUhEThySpiIlDklTExCFJKmLikCQVMXFIkoqYOCRJRUwckqQiJg5JUhEThySpiIlDklTExCFJKuKLnFaZL1SSNO5aSxwRsR64ApgCtmbmpX3nDwT+GHgt8H3grMx8ICJeCtwN3FNXvTUz391W3JKk52slcUTEFHAl8CZgFtgeEdsy866eaucBP8jMl0fEBuBjwFn1ufsy89g2YpUk7VlbcxwnAvdm5v2Z+TRwHXBmX50zgWvqn68HTo2INS3FJ0lqqK3EsQ54qOd4ti5bsk5m7gAeBw6rzx0VEbdHxF9ExBsGHawkaXltzXEsNXJYbFjnEeDIzPx+RLwW+HxEvDIzf9RbMSI2AhsBMpNut7sKYS9venp6yWt0OlMDve4wrV27lk6nM+wwdul2Dxz4NZb7nIelre/XqH3WbVjtPrfx/dxXK/1+t5U4ZoEjeo4PB+aWqTMbEdPAC4HHMnMReAogM78WEfcBPwfc1ts4M7cAW+rDxfn5+VXvRK9ut8tS11hYGN9VVZ1Oh4WFhWGHsctFF6287aZNTzSqt9znPCxtfb9G7bNuw2r3eX6+2XdsmPq/3zMzM43atZU4tgNHR8RRwMPABuDsvjrbgHOBrwBvBW7OzMWI+LtUCeSZiHgZcDRwf0txS5L6tDLHUc9ZnA/cSLW0NjPzzoi4JCLOqKt9CjgsIu4FLgTeX5efDHwjIv6GatL83Zn5WBtxS5J219pzHJl5A3BDX9nFPT//GHjbEu0+B3xu4AFKkhpxyxFJUhEThySpiHtVSdIA7Mu+dU1X/Q2LIw5JUhEThySpiIlDklTExCFJKmLikCQVMXFIkoqYOCRJRUwckqQiJg5JUhEThySpiIlDklTEvao0kZruI9TpTC351r192UtoX/YwkkaBIw5JUhEThySpiLeqltDkVsJytzAkaV+N+pbsjjgkSUVMHJKkIiYOSVIRE4ckqYiJQ5JUxFVV0gr4EJ8mmSMOSVIRE4ckqYiJQ5JUxMQhSSrS2uR4RKwHrgCmgK2ZeWnf+QOBPwZeC3wfOCszH6jPXQScBzwDXJCZN7YVtyTp+VoZcUTEFHAl8CvAMcDbI+KYvmrnAT/IzJcD/xH4WN32GGAD8EpgPfDJ+vdJkoagrVtVJwL3Zub9mfk0cB1wZl+dM4Fr6p+vB06NiDV1+XWZ+VRmfhu4t/59kqQhaCtxrAMe6jmercuWrJOZO4DHgcMatpUktaStOY41S5QtNqzTpC0RsRHYCJCZzMzMlMa4y+bNTWtO4kNg9nlyTGK/x6HPZX1Yyd+VbY04ZoEjeo4PB+aWqxMR08ALgccatiUzt2Tm8Zl5PFWyGeifiPhaG9cZpT/2eXL+TGK/7fOuP3vV1ohjO3B0RBwFPEw12X12X51twLnAV4C3Ajdn5mJEbAM+HREfB2aAo4GvthS3JKlPKyOOes7ifOBG4O6qKO+MiEsi4oy62qeAwyLiXuBC4P112zuBBO4Cvgj8VmY+00bckqTdtfYcR2beANzQV3Zxz88/Bt62TNuPAB8ZaIDltgw7gCGwz5NjEvttnxtas7i42zyzJEnLcssRSVIR38exDyLicuB04GngPuA3M/OHw41qMPa2Zcy4iYgjqLbA+RngWWBLZl4x3KjaUe/McBvwcGa+ZdjxDFpEvAjYCryKaqn/OzPzK8ONavAi4t8C76Lq8x1Uf3/9uElbRxz75ibgVZn5GuBbwEVDjmcgGm4ZM252AJsy8xXA64HfmoA+7/RvqBaxTIorgC9m5i8A/5AJ6HtErAMuAI7PzFdR/YNwQ9P2jjj2QWb+ec/hrVTLiMfRri1jACJi55Yxdw01qgHKzEeAR+qfn4iIu6l2LBjbPgNExOHAr1ItRrlwyOEMXEQcApwMvAOg3hLp6WHG1KJp4Kci4idAhyWej1uOI47V807gC8MOYkAmetuXiHgpcBzw10MOpQ2fAN5HdXtuErwMeBT4TxFxe0RsjYifHnZQg5aZDwO/DzxI9Q+kx/v+IbxHjjj2IiL+J9V97n4fyMz/Vtf5ANWtjWvbjK1FSz1NOhHL8SLiIOBzwG9n5o+GHc8gRcRbgO9l5tci4o3Djqcl08A/At6TmX8dEVdQPUP2u8MNa7Ai4lCquwZHAT8EPhsR52TmnzRpb+LYi8w8bU/nI+Jc4C3AqZk5rn+ZNtr2ZdxExAFUSePazPzTYcfTgpOAMyLinwEvAA6JiD/JzHOGHNcgzQKzmblzNHk99cPHY+404NuZ+ShARPwp8E8AE8eg1SuNfgc4JTMXhh3PADXZMmas1Fv6fwq4OzM/Pux42pCZF1Ev8KhHHO8d86RBZn4nIh6KiJ/PzHuAUxnzeazag8DrI6IDPEnV79uaNnaOY9/8AdVWlDdFxNcj4o+GHdAgLLdlzHCjGriTgN8Afqn+bL9e/0tc4+c9wLUR8Q3gWOCjQ45n4OoR1vXA/6FairuWgqfIfXJcklTEEYckqYiJQ5JUxMQhSSpi4pAkFTFxSJKKmDgkSUVMHJKkIiYOSVIRtxyRWhARlwFHZ+av1ceXU22utz4zfzLU4KRCPjkutSAiDqN6S+QbqV8MBfzTzHx8mHFJK2HikFoSEf8O+HXghVRJ46E9t5BGk3McUntuB14NXGTS0P7MxCG1ICJeDfwhcA3V2yKl/ZaJQxqwiFgH/BnwbuBfA6+eoDfsaQyZOKQBiohDgBuAj2fmtvqFX5cDHxluZNLKOTkuSSriiEOSVMTEIUkqYuKQJBUxcUiSipg4JElFTBySpCImDklSEROHJKmIiUOSVOT/A0qRe61GS7WjAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# generate X matrix (n-samples, n-features = 1)\n",
    "N = 500\n",
    "xs = np.concatenate((np.random.normal(0, 1, int(0.3 * N)),\n",
    "                    np.random.normal(5, 1, int(0.7 * N))))\n",
    "X = xs.reshape(len(xs,), 1)\n",
    "\n",
    "X_plot = np.linspace(-5, 10, 1000).reshape(1000, 1)\n",
    "fig, ax = plt.subplots();\n",
    "\n",
    "ax.hist(X[:, 0], bins=20, density=True, color='blue', alpha=0.5);\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel('event rate');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'sklearn'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-22-79ab8c0c2691>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0;32mfrom\u001b[0m \u001b[0msklearn\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mneighbors\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mKernelDensity\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mfig\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0max\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msubplots\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfigsize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m6\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m6\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mhist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbins\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m20\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdensity\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcolor\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'blue'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0malpha\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m;\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'sklearn'"
     ]
    }
   ],
   "source": [
    "from sklearn.neighbors import KernelDensity\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(6, 6))\n",
    "\n",
    "ax.hist(X[:, 0], bins=20, density=True, color='blue', alpha=0.2);\n",
    "\n",
    "colors = ['red', 'blue', 'yellow', 'black']\n",
    "kernels= ['gaussian', 'tophat', 'epanechnikov', 'linear']\n",
    "kernelwidth = 0.75\n",
    "for i, kernel in enumerate(kernels):\n",
    "    kde = KernelDensity(kernel=kernel, bandwidth=kernelwidth)\n",
    "    # fit the kde dataset\n",
    "    kde.fit(X)\n",
    "    # score_samples returns the log for x-points, get back the pdf\n",
    "    pdfs = np.exp(kde.score_samples(X_plot))\n",
    "    integral = np.sum(pdfs)*(X_plot[1]-X_plot[0])[0]\n",
    "    print('pdf {:s} integral {:5.3f}'.format(kernel, integral))\n",
    "    ax.fill(X_plot[:, 0], pdfs, color=colors[i], alpha=0.3, label=kernel)\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel('$f(x)$')\n",
    "ax.legend(loc=2)\n",
    "fig.savefig('figs/gen_nonparam_pdf.png');"
   ]
  },
  {
   "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\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"
   ]
  },
  {
   "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}