lecture_regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Statistical Methods in Python For Rare Event Searches\n",
    "\n",
    "*Summer Course at Canfranc Underground Laboratory (LSC) July 2019*\n",
    "\n",
    "### Lecture V\n",
    "\n",
    "## Hypothesis test and regression\n",
    "\n",
    "*Author: Jose A. Hernando*, June 2019\n",
    "\n",
    "*Instituto Galego de Altas Enerxías. Universidade de Santiago de Compostela, Spain.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These is a series of lectures on statistical method for rare searches in Particle Physics in Python. 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",
    "We will use the Python scientic toolkits, Matplotlib, Numpy, Scipy, that are distributed with Anaconda Python [9]. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "%reload_ext autoreload\n",
    "%autoreload 2\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",
    "import scipy.optimize as optimize\n",
    "\n",
    "# possible styles: ggplot (simplicity), bmh (scientify data), \n",
    "plt.style.context('seaborn-colorblind')\n",
    "\n",
    "import htcore    as htcore\n",
    "import htcatalog as htcata\n",
    "import httoys    as httoys\n",
    "import htplot    as htplot\n",
    "\n",
    "# ploting images\n",
    "from IPython.display import display, Image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Last Execution  Fri Jul  5 11:03:56 2019\n"
     ]
    }
   ],
   "source": [
    "import time\n",
    "\n",
    "print(' Last Execution ', time.asctime())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "As we have seen in the previous lecture, there is a relation between hypothesis testing, CI and regression. To obtain a CI using FC ordering we need to estimate the parameter $\\mu$.\n",
    "\n",
    "We are going to revisit here some of **regression** techniques in particular the extended maximum likelihood, EML.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "**regression** is the task of estimating the parameters, $\\mu$, of a model using data, $x$.\n",
    "The most common regression methods are the Least Squares (LS) and the Maximum Likelihood (ML).\n",
    "\n",
    "**goodness of the fit** is the task to quantify the agreement between a dataset and a model or between datasets. \n",
    "\n",
    "A typical example is to check if the model obtained by regression is consistent with the original dataset. Another example is to check if two datasets are consistent. \n",
    "\n",
    "The most common goodness of the fit used are the $\\chi^2$ and the Kolmogorov-Smirnov. The later is only valid in 1D problems. \n",
    "\n",
    "Extending goodness of the fit to multivariate datasets is not simple. It is still an open issue"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extended Maximun Likelihood \n",
    "\n",
    "The MLE estimator maximizes:\n",
    "\n",
    "$$\n",
    " \\ln \\mathcal{L}(\\mu) = \\sum_{i=1}^n \\log g(x_i | \\mu).\n",
    "$$\n",
    "\n",
    "where $g(x | \\mu)$ is the pdf of x, and $\\{x\\}$ a n-size sample.\n",
    "\n",
    "\n",
    "**Remember** MLE may not be optimal! and that there is no goodness-of-fit for MLE!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extended Maximum Likelihood\n",
    "\n",
    "The extended maximun likelihood is convenient when we have events of different clases, each one with its one pdf.\n",
    "\n",
    "Consider the case of $n$-size sample $x$, each $x_i$ distributed with a pdf $g(x | \\mu)$, and $n$ is poisson distributed $g(n | \\nu)$. The complete likelihood is:\n",
    "\n",
    "$$\n",
    "\\mathcal{L}(\\nu,\\mu) = \\frac{\\nu^n}{n!} e^{-\\nu} \\Pi_{i=1}^n g(x_i|\\mu)\n",
    "$$\n",
    "\n",
    "And it is called **Extended Maximum Likelihood function**.\n",
    "\n",
    "In the cases where $\\nu$ depends on $\\mu$, we must use EML. If $\\nu$ is independent of $\\mu$, then $\\hat{\\nu} = n$, and we recover the ML result.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Signal and background samples combined\n",
    "\n",
    "Consider two types of samples, background and signal, each one has a given pdf in one variable $x$:\n",
    "\n",
    "$$\n",
    "g_b(x | \\mu_b), \\; g_s(x | \\mu_s)\n",
    "$$ \n",
    "\n",
    "We measure a $n$-size sample of $x$ Each $x_i$ can be from signal or background. And signal and background are poisson distributed with means $\\nu_s, \\, \\nu_b$.\n",
    "\n",
    "The number of events is poisson distributed:\n",
    "\n",
    "$$\n",
    "g_e(n|\\nu_s,\\nu_b) = \\frac{(\\nu_s+\\nu_b)^n}{n!} e^{-(\\nu_s+\\nu_b)}\n",
    "$$\n",
    "\n",
    "And the $x$ variables is distributed with the weighted pdfs:\n",
    "\n",
    "$$\n",
    "g(x|\\nu_s,\\nu_s, \\mu_s, \\mu_b) = \\frac{\\nu_s}{\\nu_s+\\nu_b} g_s(x | \\mu_s) + \\frac{\\nu_b}{\\nu_s+\\nu_b} g_b(x | \\mu_b)\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The complete, extended, pdf is:\n",
    "$$\n",
    "g(n, x | \\nu_s,\\nu_b, \\mu_s, \\mu_b) = g_e(n|\\nu_s,\\nu_b) \\, g(x|\\nu_s,\\nu_b, \\mu_s, \\mu_b)\n",
    "$$\n",
    "\n",
    "That is:\n",
    "\n",
    "$$\n",
    " \\frac{(\\nu_s+\\nu_b)^n}{n!} e^{-(\\nu_s+\\nu_b)} \\; \\left( \\frac{\\nu_s}{\\nu_s+\\nu_b} g_s(x | \\mu_s) + \\frac{\\nu_b}{\\nu_s+\\nu_b} g_b(x | \\mu_b) \\right)\n",
    "$$\n",
    "\n",
    "where the fraction of signal and background events, $f_s, \\, f_b$, are:\n",
    "\n",
    "$$\n",
    "f_s = \\frac{\\nu_s}{\\nu_s+\\nu_b}, \\; f_b = \\frac{\\nu_b}{\\nu_s+\\nu_b}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Consider an dataset with $n$ samples of a measurement $x$. The samples can be either background or signal. For the background $x$ follows and exponential and for the signal a gaussian. Use EML to compute the number of signal and background events in a dataset.\n",
    "\n",
    "Consider the case of the exponential $\\tau = -0.1$, and gaussian with mean $\\mu = 2.5$ and sigma $\\sigma = 0.1$; $\\nu_b = 100$ events for the background and $\\nu_s = 20$ for the signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAADMBJREFUeJzt3W9oXfUdx/HPZ23FOB2BmQ1JvcsGo0wE7biIIzC26LZqxD3ZQEFhw5Enm1QQJD70WR6JezAGQZ0DnSL+gWGdWyFKsbg6q1VaU8FJho1unUjwz4JS992De1uympv7u8k99/TbvF8QvPfmnHO/Ob19ezg5hzoiBADI4wt1DwAA6A3hBoBkCDcAJEO4ASAZwg0AyRBuAEiGcANAMoQbAJIh3ACQzNYqNnrhhRfG2NhYFZsGgLPSwYMH34uIkZJlKwn32NiYXnrppSo2DQBnJdv/KF2WUyUAkAzhBoBkCDcAJEO4ASAZwg0AyXQNt+0dtg+t+PrA9m2DGA4A8HldLweMiDckXS5JtrdIWpT0ZMVzAQA66PVUyVWS/h4RxdcbAgD6q9dw3yDp4SoGAQCUKb5z0vY5kq6XdGeH709JmpKkRqPRl+GqMD4zp8WlZUnS6PCQ9k9P1DwRAPSmlyPuayS9HBH/Wu2bETEbEc2IaI6MFN1uX4vFpWUtzExqYWbyVMABIJNewn2jOE0CALUrCrft8yT9QNIT1Y4DAOim6Bx3RPxH0pcrngUAUIA7JwEgGcINAMkQbgBIhnADQDKEGwCSIdwAkAzhBoBkCDcAJEO4ASAZwg0AyRBuAEiGcANAMoQbAJIh3ACQDOEGgGQINwAkQ7gBIBnCDQDJEG4ASIZwA0AyhBsAkikKt+1h24/ZPmp73vZ3qh4MALC6rYXL/VrSMxHxE9vnSDqvwpkAAGvoGm7bX5L0XUk/k6SI+FTSp9WOBQDopOSI+xuS/i3pd7Yvk3RQ0u6I+HjlQranJE1JUqPR6PecQO3GZ+a0uLQsSRodHtL+6Yk1XweqUnKOe6ukb0v6bUTslPSxpOnTF4qI2YhoRkRzZGSkz2MC9VtcWtbCzKQWZiZPhXqt14GqlIT7mKRjEXGg/fwxtUIOAKhB13BHxD8lvW17R/ulqyS9XulUAICOSq8quVXSQ+0rSt6S9PPqRgIArKUo3BFxSFKz4lkAAAW4cxIAkiHcAJAM4QaAZAg3ACRDuAEgGcINAMkQbgBIhnADQDKEGwCSIdwAkAzhBoBkCDcAJEO4ASAZwg0AyRBuAEiGcANAMoQbAJIh3ACQDOEGgGQINwAkU/SPBdtekPShpM8knYgI/uFgAKhJUbjbvh8R71U2CQCgCKdKACCZ0nCHpL/YPmh7qsqBAABrKz1VMh4R79j+iqS9to9GxL6VC7SDPiVJjUajz2NuzPjMnBaXliVJo8NDNU+Ds8Ho8JDGpvecegwMUlG4I+Kd9n+P235S0hWS9p22zKykWUlqNpvR5zk3ZHFpWQszk3WPgbPI/umJukfAJtb1VIntL9q+4ORjST+UdLjqwQAAqys54v6qpCdtn1z+DxHxTKVTAQA66hruiHhL0mUDmAUAUIDLAQEgGcINAMkQbgBIhnADQDKEGwCSIdwAkAzhBoBkCDcAJEO4ASAZwg0AyRBuAEiGcANAMoQbAJIh3ACQDOEGgGQINwAkQ7gBIBnCDQDJEG4ASIZwA0AyhBsAkikOt+0ttl+x/VSVAwEA1tbLEfduSfNVDQIAKFMUbtvbJU1KurfacQAA3WwtXO4eSXdIuqDTAranJE1JUqPR2PhkksZn5rS4tLzmMqPDQ9o/PdGX9+s2w3rea6PrI5fR4SGNTe859Zg/b1Sha7htXyfpeEQctP29TstFxKykWUlqNpvRj+EWl5a1MDO55jIn/5JUZeUM63mvja6PXFaGmj9vVKXkVMm4pOttL0h6RNKE7QcrnQoA0FHXcEfEnRGxPSLGJN0gaS4ibqp8MgDAqriOGwCSKf3lpCQpIp6T9FwlkwAAinDEDQDJEG4ASIZwA0AyhBsAkiHcAJAM4QaAZAg3ACRDuAEgGcINAMkQbgBIhnADQDKEGwCSIdwAkAzhBoBkCDcAJEO4ASAZwg0AyRBuAEiGcANAMoQbAJLpGm7b59p+0farto/YvmsQgwEAVlfyr7x/ImkiIj6yvU3S87b/FBF/rXg2AMAquoY7IkLSR+2n29pfUeVQAIDOis5x295i+5Ck45L2RsSBascCAHRScqpEEfGZpMttD0t60valEXF45TK2pyRNSVKj0Vj3QOMzc1pcWpYkjQ4PbWjd/dMTlc3Qr/fqZY5e3qPXdQfx82S0kc8jUJWicJ8UEUu2n5O0S9Lh0743K2lWkprN5rpPpSwuLWthZnLD645N71nvCEUz9Ou9qnqPXtcdxM+T0UY+j0BVSq4qGWkfacv2kKSrJR2tejAAwOpKjrgvkvR721vUCv2jEfFUtWMBADopuarkNUk7BzALAKAAd04CQDKEGwCSIdwAkAzhBoBkCDcAJEO4ASAZwg0AyRBuAEiGcANAMoQbAJIh3ACQDOEGgGQINwAkQ7gBIBnCDQDJEG4ASIZwA0AyhBsAkiHcAJAM4QaAZAg3ACTTNdy2L7b9rO1520ds7x7EYACA1W0tWOaEpNsj4mXbF0g6aHtvRLxe8WwAgFV0PeKOiHcj4uX24w8lzUsarXowAMDqejrHbXtM0k5JB6oYBgDQXcmpEkmS7fMlPS7ptoj4YJXvT0makqRGo9G3AbsZHR7S2PSeU4/Xu+7pr/djO6XbGp+Z0+LS8qnl909P9PT+vcxXsv219ku3dav+WTLpdb+vB/t7cyoKt+1takX7oYh4YrVlImJW0qwkNZvN6NuEXWzkg9qvD/lGt7O4tKyFmUlJ6vg/gI1YOV/J9jv9PCXrVv2zZNLrfl8P9vfmVHJViSXdJ2k+Iu6ufiQAwFpKznGPS7pZ0oTtQ+2vayueCwDQQddTJRHxvCQPYBYAQAHunASAZAg3ACRDuAEgGcINAMkQbgBIhnADQDKEGwCSIdwAkAzhBoBkCDcAJEO4ASAZwg0AyRBuAEiGcANAMoQbAJIh3ACQDOEGgGQINwAkQ7gBIBnCDQDJEG4ASKZruG3fb/u47cODGAgAsLaSI+4HJO2qeA4AQKGu4Y6IfZLeH8AsAIACW/u1IdtTkqYkqdFo9Guz6zY6PKSx6T2nHp8JVs50+uvrXff0ZfZPTxTPMz4zp8Wl5eIZ1lp3tfctmXmt9Tfy3iXLl75+Nul1X6C7OvZd38IdEbOSZiWp2WxGv7a7XmfiB28jM5WsWxLJlRaXlrUwM7mueVau2+l9S3/ejcxdsm6n5UteP9v0ui/QXR37jqtKACAZwg0AyZRcDviwpBck7bB9zPYt1Y8FAOik6znuiLhxEIMAAMpwqgQAkiHcAJAM4QaAZAg3ACRDuAEgGcINAMkQbgBIhnADQDKEGwCSIdwAkAzhBoBkCDcAJEO4ASAZwg0AyRBuAEiGcANAMoQbAJIh3ACQDOEGgGQINwAkUxRu27tsv2H7TdvTVQ8FAOisa7htb5H0G0nXSLpE0o22L6l6MADA6kqOuK+Q9GZEvBURn0p6RNKPqx0LANBJSbhHJb294vmx9msAgBo4ItZewP6ppB9FxC/az2+WdEVE3HraclOSptpPd0h6o8+zXijpvT5vMyP2Qwv7oYX90HI27IevRcRIyYJbC5Y5JuniFc+3S3rn9IUiYlbSbNF462D7pYhoVrX9LNgPLeyHFvZDy2bbDyWnSv4m6Zu2v277HEk3SPpjtWMBADrpesQdESds/0rSnyVtkXR/RBypfDIAwKpKTpUoIp6W9HTFs3RT2WmYZNgPLeyHFvZDy6baD11/OQkAOLNwyzsAJJMi3NxyL9m+3/Zx24frnqUuti+2/aztedtHbO+ue6Y62D7X9ou2X23vh7vqnqlOtrfYfsX2U3XPMihnfLi55f6UByTtqnuImp2QdHtEfEvSlZJ+uUk/C59ImoiIyyRdLmmX7StrnqlOuyXN1z3EIJ3x4Ra33EuSImKfpPfrnqNOEfFuRLzcfvyhWn9ZN91dvNHyUfvptvbXpvxlle3tkiYl3Vv3LIOUIdzcco/PsT0maaekA/VOUo/26YFDko5L2hsRm3I/SLpH0h2S/lv3IIOUIdxe5bVNeXSBFtvnS3pc0m0R8UHd89QhIj6LiMvVupP5CtuX1j3ToNm+TtLxiDhY9yyDliHcRbfcY3OwvU2taD8UEU/UPU/dImJJ0nPanL//GJd0ve0FtU6hTth+sN6RBiNDuLnlHpIk25Z0n6T5iLi77nnqYnvE9nD78ZCkqyUdrXeqwYuIOyNie0SMqdWFuYi4qeaxBuKMD3dEnJB08pb7eUmPbsZb7m0/LOkFSTtsH7N9S90z1WBc0s1qHVkdan9dW/dQNbhI0rO2X1PrwGZvRGyaS+HAnZMAkM4Zf8QNAPh/hBsAkiHcAJAM4QaAZAg3ACRDuAEgGcINAMkQbgBI5n+fRLE/higXpwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "nug, nub, tau, gmu, gsig, size = 20., 100., -0.1, 2.5, 0.1, 1\n",
    "par = [nug, nub, tau, gmu, gsig]\n",
    "rv  = htcata.extended_norm_expon\n",
    "xs = rv.rvs(*par, size = size)\n",
    "xs = np.concatenate(xs); x0 = xs.flatten() # add several datasets\n",
    "plt.hist(xs, 100, histtype = 'step');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "parbest  [25.39506292 93.60492481 -0.1         2.5         0.1       ]\n"
     ]
    }
   ],
   "source": [
    "mask, masknu = (True, True, False, False, False), (False, True, False, False, False)\n",
    "ht = htcore.htcomposite(rv, par, mask = mask, masknu = masknu)\n",
    "parbest = ht.parbest(x0)\n",
    "print('parbest ', parbest)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: Check is the case fulfills the Wilk's condition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chi2  1.0728216837730762\n",
      "p-value mean : 0.5175868714039004 ; std * sqrt(12.) : 0.9952688722019466\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAIeCAYAAACyQ9gRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzt3XeYVOXZx/HvvTQRURQwIsVduqAryCogrw0VV6JAohGwoaJrAzQRX0gs2KJoYhcFVEQNxYYGAgrYK4SiUgVXNLpAgooYC1Lv948dfMd1F2Z2ypnZ8/tc117MnDnlfmaH+e15zjnPMXdHRETCJyfoAkREJBgKABGRkFIAiIiElAJARCSkFAAiIiGlABARCSkFgIhISCkARERCSgEgIhJSCgARkZCqHnQBO9OgQQPPzc0NugwRkayyYMGCL9294a7my+gAyM3NZf78+UGXISKSVczsX7HMl5FdQGZ2ipmN/eabb4IuRUSkysrIAHD3ae5etNdeewVdiohIlZWRASAiIqmX0ccARERSacuWLZSUlPDjjz8GXUql7LbbbjRp0oQaNWpUankFgIiEVklJCXXr1iU3NxczC7qcuLg7X331FSUlJeTl5VVqHeoCEpHQ+vHHH6lfv37WffkDmBn169dPaO9FASAioZaNX/47JFq7AkBEJKSq7DGAbiNfYfWGjQA0rlebt4d3D7giEZHYPf/880yfPp1169Zx2WWX0aNHj6RvIyMDwMxOAU5p2bJlpdexesNGPh35awByh09PUmUiIunRp08f+vTpw9dff83QoUNTEgAZ2QWkC8FERErdfPPNXHbZZSlZd0YGgIhIWHz//fc0btyYP/3pTwDMmzePDh06sHHjRoYNG8ZJJ53EoYcempJtKwBERAJUp04dFi1axKRJk/jhhx84//zzGT9+PA899BAvvfQSzzzzDKNHj07JtjPyGEAybNnwbwYPHsyQIUOCLkVEZKfq16/P7rvvzsCBAzn77LPp0KEDHTp0SPn3V5UNALZv4/777yc/Px/YP+hqRCTDXXHFFbz//vtJXWeHDh24++67Y5o3Pz+ftWvXMmHChKTWsDNVtguo+t77U79+fd59992gSxER2akvvviCV199ldNOO42cnPR9LVfZPQAzo0uXLqUB0PvUoMsRkQwX61/qqTBw4EC6d+/OBx98kNbtVtk9AICuXbvy4Ycfsm3jt0GXIiJSrjFjxlC7dm1uu+22tPdYVOkAOOKIIwDYvGZFwJWIiPzSRx99xJ133sno0aNp2rQpjRo1SskFXxWpsl1AAIcddhg5OTlsWvNh0KWIiPxCq1atWLHi//9AnT17dlq3n7Y9ADOrY2aPmdlDZnZmOra5xx57kJ+fz6bVCgARkbISCgAzG2dm68xsSZnphWa2wsyKzWx4ZPJvgWfc/UKgVyLbjUfXrl3ZtHYF27ZtS9cmRUSyQqJ7AOOBwugJZlYNGAWcBLQD+ptZO6AJ8HlktrR9G3ft2hXfvJFly5ala5MiIlkhoQBw9zeA9WUmHw4Uu/sqd98MTAZ6AyWUhkDC241H165dAXQ9gIhIGan4Im7M//+lD6Vf/I2BKcCpZvYgMK2ihc2syMzmm9n8L774IuFiWrRoQU7tPRUAIlIudw+6hEpLtPZUnAVU3j3K3N2/B87b1cLuPhYYC1BQUJDwb8bMqNW4rQJARH5ht91246uvvsrK+wLvuCn8brvtVul1pCIASoCmUc+bAGviWUEybggTrdb+bVnxxuOsX7+effbZJynrFJHs16RJE0pKSkhGb0MQdtttN5o0abLrGSuQigCYB7QyszxgNdAPOCOeFbj7NGBaQUHBhckoqFbjtgC89dZb9OqVthOQRCTD1ahRg7y8vKDLCEyip4FOAt4F2phZiZkNdPetwCBgJrAceMrdlyZeauXV2v9AGjZsyMMPPxxkGSIiGSXRs4D6u3sjd6/h7k3c/ZHI9Bnu3trdW7j7n+Ndr5mdYmZjv/nmm0TK+0mTBnuyuWV3pk37B53+929JWaeISLbLyLGAkn1P4LeHd2fp03dQvXo1Pnrt2aSsU0Qk22VkAKRC48aNOe200/hu0Wy+++67oMsREQlcRgZAsruAdhgyZAi+6XueeOKJpK5XRCQbZWQAJLsLaIcuXbpQc79W3HfffVl98YeISDJkZACkiplRt9MpLF++nJkzZwZdjohIoEIVAACtupxA9b1+Ra+zL6Lrn2cFXY6ISGAyMgBSdQwA4N1rCnnmsTFs+fIzVrzyTNLXLyKSLTIyAFJ1DGCHXr16UVhYyIa3JvCf//wnJdsQEcl0GRkAqWZm3HPPPfjWzQwfPnzXC4iIVEGhDACA1q1bs+dhfRg/fjxvv/120OWIiKRdRgZAKo8BRNvriL4ccMABnH322aR6WyIimSYjAyDVxwB2yKlZm4kTJ/LZZ59xySWX6NoAEQmVjAyAdDriiCO4/vrrmTRpEo8//njQ5YiIpE0q7geQNRrXq03u8On49nxqNT2I8y68mD+99SO5LVrx9vDuQZcnIpJSoQ6A6C/5kiGH0qFDB+q+eRef1bkhwKpERNIjI7uA0nUQOFqTJk2YMmUKq1at4svnR7Jly5a0bVtEJAgZGQDpOghc1lFHHcWYMWP48V/vs0/nPhww7B90G/lKWmsQEUmXjAyAIJ133nlcddVVfPfeDC7/1UpWb9gYdEkiIimhACjHrbfeyqmnnsqVV17Jd4tmB12OiEhKhPogcEWqVavGhAkT+Pbbb5n14n00/M3u1GnTjcb1auvsIBGpMrQHUIFatWoxZcoUjujahW+m/5XRx1RTd5CIVCkZGQBBnAVUnjp16jB9+nTat29P7969+eHjeYHWIyKSTBkZAEGdBVSeevXq8fLLL3PwwQfz5XN/Zt9TryV3+HSdHSQiWU/HAGKwzz778NJLL1FYWMiCqSMZNXEiVy34/9e7jXzlp+4hHScQkWyRkXsAmahevXrMmjWLzp0707dvX759b8ZPr63esJFPR/6aT0f+WscJRCRrKADisOeeezJr1ix69uzJ+lkPcNNNN2kEURHJWuoCitPuu+/Oc889x96HFnLdddexdu1afI/CoMsSEYmbAqASatSoQf2eV3BZzwJuv/12ajefx7fXHkvdunWDLk1EJGbqAqoksxxuu+02Ro8ezcZPFnLkkUdSUlISdFkiIjHTHkCCLrroIq57+T+sevGvHH744fjxQ4FfV3p9OqNIRNIlI/cAMuVCsFjVbt6Jt99+m1q1avHvicP529/+Vul16YwiEUmXjAyATLoQrCI77iaWO3w6jevV5uCDD2bevHnUatyWs88+m6FDh7J169agyxQRqZC6gCqpvK6ZBg0a8KvTb+LkH2Zxxx13sHDhQp588kkaNmwYQIUiIjuXkXsA2cyqVWfUqFGMHz+ed999l06dOvHPf/4z6LJERH5BewApMmDAAPLz8/ntb3/LkUceyZ133smE/7ZlzTc/AjrAKyLBUwCkUMeOHVmwYAHnnHMOgwYNYve2R7J27j/Yc889yR0+PejyRCTk1AWUYqeMfZ/F7S+i3tHn8sPKt+nUqRMLFy4MuiwREQVAqq3esJF/3XYKX7/2KG+89hobN26ka9eu/Hf+VI0jJCKBUgCk0ZFHHsn7779Pjx49+PrlsfTp04cvv/wy6LJEJKQUAGnWoEEDpk6dyt7dL+TFF18kPz+f2bN143kRST8FQJJFXyC24yKxssyMPQ/rzdy5c6lXrx49evRg6NChbNq0KYCKRSSsdBZQksVzameHDh2YP38+V111FXfccQezZs1ic0HRT6/vCJMdj3XaqIgkU9r2AMysuZk9YmbPpGub2WD33Xdn1KhRTJ8+nXXr1rH28d/z17/+lW3btvH28O4aF0hEUiamADCzcWa2zsyWlJleaGYrzKzYzIbvbB3uvsrdByZSbFVStqvo0jecGqffyT5tu3DVVVdx7LHHsmrVqqDLFJEqLNYuoPHA/cDjOyaYWTVgFHACUALMM7OpQDXg1jLLn+/u6xKutgqpqDvHvT9PPPEEgwcPJj8/nzvuuIOioqJy5xURSURMewDu/gawvszkw4HiyF/2m4HJQG93X+zuJ5f50Zd/jMyMc845hyVLltC1a1cuvvhiTjzxRLb+V2+hiCRXIscAGgOfRz0viUwrl5nVN7PRQEcz++NO5isys/lmNv+LL75IoLzs1rRpU2bNmsWDDz7IO++8w5pHLuOhhx76xcVj3Ua+8lM3UreRrwRUrYhko0QCwMqZVuGlre7+lbtf7O4t3L1sF1H0fGPdvcDdC8I+jLKZcfHFF7N48WJq7teKoqIievTowaeffvrTPLqBjIhUViIBUAI0jXreBFiTWDmlsu2OYKmWl5fHoRfdwT49LuXlN96mRZsDGTVqFNu3bw+6NBHJYokEwDyglZnlmVlNoB8wNRlFZcMdwdLtnT8dz1czR/HJyuXUbNyOQYMGcfTRR7PlK92IXkQqJ9bTQCcB7wJtzKzEzAa6+1ZgEDATWA485e5LU1eqABxwwAF0vPB26vf8Pe/Mf5+14wdzyy23sGXLlqBLE5EsE+tZQP3dvZG713D3Ju7+SGT6DHdvHenX/3OyilIX0M6988fj+HL6naz+5CNO7dObq6++mk6dOrFpzYqgSxORLJKRYwGpCyg2++23H08//TTPP/8869ev599PDOXyyy/n22+/Dbo0EckCGguoCujduzfHHnssTY49g/vuu48pU6Zw//3307t3b6D0VNFdnSGksYZEwicjA8DMTgFOadmyZdClZI0999yTfU64hJn3X01RURF9+vThN7/5Dffee+9Pp4rujG5RKRI+6gKqYrp27crChQtpduIFPD9tOs1atMaWzGDr1q2VWp8uNBOpujIyACQxNWrUwDr0oXjFcgqPO4ZPpz9A586dmTdvXtzr0oVmIlVXRgaAzgJKjubNmzN9+nSeeuop1q5dS+fOnbn00kvZsGFD0KWJSAbIyABQF1DymBm/+93v+PDDDxkyZAhjxoyhTZs2PP744xXelD662yf6jmbRQ1irO0gk+2VkAEjy7bnnntx9993Mnz+fvLw8BgwYwNFHH82SJUt+MW90t0/0mUG6QY1I1aIACJmOHTvyzjvv8NBDD7F06VI6dOjAlVdeyfZNPwRdmoikWUYGgI4BpFZOTg4XXHABK1eu5Pzzz+euu+5izcMXM3HixAq7hUSk6snIANAxgPSoX78+Y8eOZc6cOVTboz5nnnkmxxxzDJu/+DTo0kQkDTIyACS9Dj/8cPY7+6+MGTOGJUuWsPbRIVx++eU6W0ikilMACACWU42ioiJWrlzJHh0Kue+++2jdujWPPPKI7jsgUkUpAORn6tevT/0el7JgwQJat27NBRdcQJcuXZg7d27QpYlIkmksIClXx44defPNN5kwYQL/+7//S5cuXRgwYAC33norjRo1int9ZQekS9bgc9Hr1YB2IvHJyD0AHQTODGbGWWedxYoVKxg2bBgTJ06kdevW3H777fjW+G5AE31tQTKvI9BQFSKVl5EBIJmlbt26jBw5kqVLl3LMMccwbNgw1oy7lGnTpum0UZEspgCQmLVq1Ypp06bxwgsvgFWjV69eFBYWsmzZsqBLE5FKUABI3AoLC9n//Pu56667mDt3Lvn5+QwePJj169cHXZqIxCEjDwJL5ewYrG3H41Qua9Wqc8UVV3DmmWdy3XXX8cADDzBhwgRuuOEGLr74YmrUqPGLA7TppIPDIrumAKhCEvmSq+yyDRs25MEHH+SSSy7h97//PUOGDOHBBx/kjjvuYPWG7bu8E1mqRN8FTXc7EylfRnYBaSyg7JOfn89LL73E3//+d7Zs2ULPnj35z1PXsXTp0qBLE5EKZGQA6DTQ7GRm9OrVi6VLl3LHHXewac0KDjnkEC677DK+/PLLoMsTkTIyMgAku9WsWZM//OEPNC4ay0UXXcSYMWNo2bIlf/3rX9m0aVPQ5YlIhAJAUqba7nsxatQoFi1axBFHHMFVV11Fu3btePbZZ3X9gEgG0EFgSbl27doxY8YMZs6cydChQznttNOom3sQjdaupNb+bZJ2lk4iZ0GJhJECQNLmxBNP5LjjjuPRRx/lmmuu4d9PXMkZZ5zB698cn5T161RPkfioC0jSqnr16lx44YUUFxdz9dVXM2XKFFY/dBHDhg1DZ32JpJcCQAJRt25dbr75ZlauXEmdtv/D7bffTsuWLbn//vvZsiW+geZEpHIUABKopk2b0uDkK1mwYAEHH3wwgwcPpn379jz33HM6UCySYhkZALoQLHwOPfRQXn75ZaZNm0b16tX57W9/y1FHHaUb0YikUEYGgC4ECycz4+STT2bRokWMHj2ajz76iC5dutC3b18+/vjjoMsTqXIyMgAk3KpXr85FF13ERx99xHXXXcc//vEPDjzwQK644gpdUSySRAoAyVh169blhhtuoLi4mHPPPZf77ruPFi1aMHLkSDZu1N2/RBKl6wCkUqIvuio7Pd7z8ctewFV2+UaNGrG0eT9+dV5HNrw+nj/+8Y+MGjWKm266Cd++T+UbkYCKhpsue+/jHarakNRVYbjtqtCGRCkApFIq+s9SmaGXo9dV0fKrN2xkzcOXAJew3xkjaVQ8hfPOO48aDXN54dgaFBYWYmZxb7uyKhpuOnp6tKo2JHVVGG67KrQhUeoCkqyzW7ODmTt3Lk899RS+ZRM9e/bkuOOOY/78+UGXJpJVFACSlcyM3/3ud+x/wQPcd999LFmyhMMOO4y+fftSXFwcdHkiWUEBIFnNqtVg0KBBFBcXc+211/50xtCgQYNYt25d0OWJZDQFgFQJe+65JzfeeCPFxcVccMEFjB49mhYtWnD99dezfdMPQZcnkpEUAFKlNGrUiAcffJBly5Zx4okncsMNN7B67IXcf//9bN68OejyRDJK2gLAzPqY2UNm9ncz65Gu7Uo4tW7dmmeeeYY5c+ZQo35TBg8ezIEHHsikSZPYvn170OWJZISYAsDMxpnZOjNbUmZ6oZmtMLNiMxu+s3W4+/PufiFwLtC30hWLxKFz5878qv+tzJgxgz322IMzzjiDTp068eKLL2qwOQm9WPcAxgOF0RPMrBowCjgJaAf0N7N2Znawmf2jzM++UYteE1lOJC3MjJNOOon33nuPJ554gg0bNnDSSSfRvXt3DTYnoRZTALj7G8D6MpMPB4rdfZW7bwYmA73dfbG7n1zmZ52Vug14wd0XJrcZIruWk5PDWWedxYcffsi9997L0qVL6dKlC6eeeirLly8PujyRtEvkGEBj4POo5yWRaRUZDBwPnGZmF1c0k5kVmdl8M5v/xRdfJFCeSPlq1arF4MGD+fjjj7n++uuZNWsWBx10EAMHDuTzzz/f9QpEqohEAqC86+4r7FR193vdvZO7X+zuo3cy31h3L3D3goYNGyZQnsjO1a1blxEjRrBq1SqGDBnC3/72N1q1asXQoUP56quvgi5PJOUSCYASoGnU8ybAmsTKKaUbwkg6NWzYkLvuuouVK1fSv39/7rrrLpo3b87NN9/Md999F3R5IimTSADMA1qZWZ6Z1QT6AVOTUZRuCCNBOOCAA3j00UdZtGgR3bt359prr6VFixbcd999bNq0KejyRJIu1tNAJwHvAm3MrMTMBrr7VmAQMBNYDjzl7ktTV6pIeuy4J/G7775Lu3btGDJkCG3atOHxxx9n27ZtQZcnkjQxDQft7v0rmD4DmJHUiijtAgJOadmyZbJXLSlWdmz/RJcvb9jpeLcRy70LyhsbvkuXLmzqcQ377tudtW88zoABA7jgymupd+TZtDzs2KS3q6I6KiMTxrqP5Z4JsUwPiyDan5H3A3D3acC0goKCC4OuReKT6Ic2lnsDxLuNWO5dUNHY8Gu++ZH/PHkt7tfw7LPPcs0117DiuT/TfM1LvNLZ6N49tlpiadfO6ohXJox1H8s9E2KZHhZBtF9jAYnEwMw47bTTWLJkCY888ghr1qzhuOOO4/jjj2fTmhVBlydSKRkZADoLSDJV9erVOf/88/noo4+4++67WbRoEf9+4kr69OnDkiVLdr0CkQySkQGgs4Ak09WqVYvLL7+cVatWsdeRZ/Hqq6+Sn5/Pl/+4g48//jjo8kRikpEBIJIt9thjD+od0Y9Vq1Zx1VVX8cOKd2jbti0XX3wxq1evDro8kZ1SAIgkQf369bntttvYv2gsF154IY888ggtW7Zk6NChfPnll0GXJ1KujAwAHQOQbFW9bn0eeOABVq5cSd++fbnrrrvIy8tjxIgR6PMsmSYjA0DHACTb5eXlMX78eJYsWUJhYSE33ngjzZs357bbbmP7lh+DLk8EyNAAEKkqDjzwQJ5++mkWLFhA586dGT58OGvGlN6iUsNLSNAUACJpcOihhzJjxgzeeustqu/TmMGDB9O6dWvGjRvH1q1bgy5PQiojA0DHAKSq6tatG7/qfyuzZs1iv/32Y+DAgbRr1073KpZAZGQA6BiAVGVmxgknnMCcOXOYOnUqtWvX5owzzuCQQw7hh5Xv6F7FkjYZGQAiYWBmnHLKKbz33ntMnjyZzZs388Vzt3DYYYfppvWSFgoAkYDl5OTQt29fli5dSv2TLuerr77ipJNO4j8Th/H6668HXZ5UYQoAkQxRvXp19sg/gRUrVvDAAw+wdcO/OeaYY37qLhJJtowMAB0EljCrWbMml1xyCfsXPcRdd93FBx98QNeuXTn55JNZuHBh0OVJFZKRAaCDwCKQU6MWV1xxBatWreLWW2/lnXfeoVOnTpx66qkaeVSSIiMDQET+3x577MHw4cP55JNPuP7665k9ezb5+fn079+fFSt0LwKpPAWASJbYa6+9GDFiBJ988gnDhg1j6tSptGvXjnPPPZdVq1YFXZ5kIQWASJapX78+t956K5988gmXX345Tz75JG3atKGoqIjPPvss6PIkiygARLLUvvvuy5133snHH3/MRRddxPjx42nVqhWDBg1i67dfBV2eZIGMDACdBSQSu/3335/777+f4uJiBgwYwJgxY1gz9kL+8Ic/sG7duqDLkwyWkQGgs4BE4tesWTPGjh3LihUr2L3tUdxzzz3k5eXx9Wvj+eor7RHIL2VkAIhI5TVv3pwGv76C5cuX06dPH/4791lyc3O59tpr+frrr4MuTzKIAkCkimrdujUTJkyg0cBRnHTSSdx8883k5eVx44036u5kAigARKq8mg2a8dRTT/HBBx/QvXt3RowYQV5eHrfeeivffvtt0OVJgBQAIiGRn5/PlClTmD9/Pt26deNPf/oTeXl53H777Xz//fdBlycBUACIhEynTp2YNm0ac+fOpaCggGHDhtG8eXPuvPNOfvjhh6DLkzRSAIiE1OGHH86LL77I22+/TX5+PldeeSUtWrTg3nvv5ccfdeP6MFAAiITcEUccwezZs3n99ddp06YNl19+OS1atGDUqFG6cX0Vl5EBoAvBRNLvqKOO4rXXXuOVV14hLy+PQYMG0bJlS0aPHs3mzZuDLk9SICMDQBeCiQTn2GOP5c0332TWrFk0bdqUSy65hFatWvHQQw+xZcuWoMuTJMrIABCRYO24cf3bb7/Niy++SKNGjSgqKqJ169aMGzcO37Y16BIlCRQAIlIhM+PEE0/k3XffZfr06TRo0ICBAwey5uGLGT9+PFu3KgiymQJARHbJzOjZsyf//Oc/mTp1Kjm16nDeeedx4IEH8t2SlxUEWUoBICIxMzNOOeUU9htwN88//zx77LEHX02/i/bt2/O3v/2Nbdu2BV2ixEEBICJxMzN69+7NggULaPibP7Hbbrtx9tln0759eyZOnKggyBIKABGptJycHHZvfQTvvfcezzzzDDVq1ODMM8/koIMO4vtlrysIMpwCQEQSlpOTw6mnnsoHH3zA008/TbVq1fhy2l84+OCDmTx5Mr5dQZCJFAAikjQ5OTmcdtppLFq0iAa9hmFm9O/fn7XjBvPkk0+yffv2oEuUKAoAEUm6nJwc6hx4JIsWLWLy5MmA069fP/Lz8/l++ZsKggyhABCRlKlWrRp9+/al0fn3M2nSJLZv386XU28jPz+fp59+GncFQZDSFgBmdqCZjTazZ8zsknRtV0SCZznV6NevH4sXL6bBKVexfft2Tj/9dNaOG8TTTz+tPYKAxBQAZjbOzNaZ2ZIy0wvNbIWZFZvZ8J2tw92Xu/vFwOlAQeVLFpFsVa1aNeq0O5rFixczadIkcOf0008v7Rr68C0FQZrFugcwHiiMnmBm1YBRwElAO6C/mbUzs4PN7B9lfvaNLNMLeAt4OWktEJGsU61a6R5Bo/Pv/+m6gS//PlJdQ2kWUwC4+xvA+jKTDweK3X2Vu28GJgO93X2xu59c5mddZD1T3f0I4MxkNkJEspPlVKN///4sWbKEBqdcxbZt29Q1lEaJHANoDHwe9bwkMq1cZnaMmd1rZmOAGTuZr8jM5pvZ/C+++CKB8kQkW+zoGlqyZAkTJ07EI8cI8vPzeeqppxQEKZJIAFg507yimd39NXcf4u4Xufuoncw31t0L3L2gYcOGCZQnItmmWrXSPYL9B45i4sSJbN++nb59+yoIUiSRACgBmkY9bwKsSaycUrojmEi47ega2nGweEcQHHzwwTz55JMaYiJJEgmAeUArM8szs5pAP2BqMorSHcFEBP7/YPHixYsjF5RBv379fhpiQkGQmFhPA50EvAu0MbMSMxvo7luBQcBMYDnwlLsvTV2pIhJWOy4oW7x4MU8++SQ5OTn079+fgw46SKOPJiDWs4D6u3sjd6/h7k3c/ZHI9Bnu3trdW7j7n5NVlLqARKQ8OTk5nH766SxatIinnnqK6tWrc+aZZ9K+fXsmTJigIIhTRg4FoS4gEdmZnJwcfve73/00+mjNmjU566yzaNeuHU888YTuUBajjAwAEZFY7Bh99P333+fZZ5+ldu3anHPOObRr147HHntMQbALGRkA6gISkXjk5OTw29/+loULF/Lcc89Rp04dzj33XNq2bcujjz7Kli1bgi4xI2VkAKgLSEQqIycnhz59+rBw4UL+/ve/s9dee3H++efTpk0bHnnkEQVBGRkZACIiiTDUzqKlAAAfhklEQVQzevXqxfz585k2bRr169fnggsuoHXr1jz00ENs3rw56BIzQkYGgLqARCQZzIyTTz6Zf/7zn0yfPp2GDRtSVFREq1at+Pa9GWzatCnoEgOVkQGgLiARSSYzo2fPnsydO5cXXniB/fffn/WzHqBly5aMGjUK3xrOPYKMDAARkVQwMwoLC3nnnXfY9/SbOOCAAxg0aBCrx1zIfffdx48//hh0iWmlABCR0DEzaud15M033+Sll16ier39GDJkCM2bN+eee+5h48aNQZeYFhkZADoGICLpYGYcd9xx/OqMkbzyyiu0adOGK664gry8PO644w6+//77oEtMqYwMAB0DEJF0MjOOPfZYXn31VV5//XUOOugghg4dSl5eHrfffjvfffdd0CWmREYGgIhIUI466iheeukl3nrrLTp27MiwYcPIzc3llltu4b///W/Q5SWVAkBEpBzdunVj5syZzJkzh86dO3P11VeTm5vLjTfeyIYNG4IuLykUACIiO9G5c2emT5/OvHnz+J//+R9GjBhBbm4uI0aMYP36srdKzy4ZGQA6CCwimaagoICpU6eycOFCunfvzo033khubi5XX301X375ZdDlVUpGBoAOAotIpurYsSNTpkzhgw8+oLCwkFtvvZXc3FyGDRvGunXrgi4vLhkZACIimW7HjeoXL17MKaecwl/+8hdyc3O58sor+fe//x10eTFRAIiIJKB9+/ZMmjSJZcuWcdppp3H33XeTl5fH5ZdfzurVq4Mub6cUACIiSdC2bVsef/xxVqxYwRlnnMEDDzxA8+bNufTSS/nss8+CLq9cCgARkSRq2bIljzzyCCtXrmTAgAE8/PDDtGzZkqKiIj755JOgy/sZBYCISArk5eUxduxYiouLufDCC3nsscdo1aoV5513Hh999FHQ5QEZGgA6DVREqopmzZoxatQoVq1axaBBg5g8eTJt27blrLPOYvny5YHWlpEBoNNARaSqady4MXfffTeffPIJf/jDH3juuedo3749ffv2ZfHixYHUlJEBICJSVe2333785S9/4dNPP2X48OHMmDGD/Px81j33Z95777201qIAEBEJQMOGDbnlllv417/+xXXXXceP/1rEoYceyuDBg9NWQ/W0bUlERH5hn3324YYbbuDR/7anqP5KWrZsybT307NtBYCISAbIqVWHa665BoDh709PzzbTshUREck4CgARkZBSAIiIhFRGBoAuBBMRSb2MDABdCCYiknoZGQAiIpJ6CgARkZBSAIiIhJS5e9A1VMjMvgD+VcnFGwDZeafm/1cV2gBVox1VoQ1QNdqhNuzaAe7ecFczZXQAJMLM5rt7QdB1JKIqtAGqRjuqQhugarRDbUgedQGJiISUAkBEJKSqcgCMDbqAJKgKbYCq0Y6q0AaoGu1QG5Kkyh4DEBGRnavKewAiIrITVTIAzKzQzFaYWbGZDQ+6nniZWVMze9XMlpvZUjO7POiaKsvMqpnZe2b2j6BrqSwzq2dmz5jZh5HfSdega4qXmf0+8llaYmaTzGy3oGuKhZmNM7N1ZrYkato+ZjbbzD6K/Lt3kDXuSgVt+Evk87TIzJ4zs3pB1FblAsDMqgGjgJOAdkB/M2sXbFVx2wpc6e4HAl2Ay7KwDTtcDiwPuogE3QO86O5tgUPIsvaYWWNgCFDg7gcB1YB+wVYVs/FAYZlpw4GX3b0V8HLkeSYbzy/bMBs4yN3zgZXAH9NdFFTBAAAOB4rdfZW7bwYmA70Driku7r7W3RdGHn9L6RdO42Crip+ZNQF+DTwcdC2VZWZ7AkcBjwC4+2Z33xBsVZVSHahtZtWB3YE1AdcTE3d/A1hfZnJv4LHI48eAPmktKk7ltcHdZ7n71sjTOUCTtBdG1QyAxsDnUc9LyMIvzx3MLBfoCMwNtpJKuRv4X2B70IUkoDnwBfBopCvrYTOrE3RR8XD31cBfgc+AtcA37j4r2KoS8it3XwulfywB+wZcT6LOB14IYsNVMQCsnGlZeaqTme0BPAtc4e7/DbqeeJjZycA6d18QdC0Jqg4cCjzo7h2B78n8LoefifSR9wbygP2BOmZ2VrBVCYCZXU1pl++EILZfFQOgBGga9bwJWbK7G83MalD65T/B3acEXU8ldAN6mdmnlHbDdTezvwVbUqWUACXuvmMP7BlKAyGbHA984u5fuPsWYApwRMA1JeI/ZtYIIPLvuoDrqRQzGwCcDJzpAZ2PXxUDYB7QyszyzKwmpQe7pgZcU1zMzCjtc17u7ncGXU9luPsf3b2Ju+dS+jt4xd2z7q9Od/838LmZtYlMOg5YFmBJlfEZ0MXMdo98to4jyw5klzEVGBB5PAD4e4C1VIqZFQLDgF7u/kNQdVS5AIgcWBkEzKT0Q/6Uuy8Ntqq4dQPOpvSv5vcjPz2DLirEBgMTzGwR0AG4JeB64hLZe3kGWAgspvT/fUZciborZjYJeBdoY2YlZjYQGAmcYGYfASdEnmesCtpwP1AXmB35/z06kNp0JbCISDhVuT0AERGJjQJARCSkFAAiIiGlABARCSkFgIhISCkARERCSgEgIhJSCgARkZBSAIiIhJQCQEQkpBQAIiIhpQAQEQkpBYCISEgpAEREQkoBICISUgoAEZGQUgCIiISUAkBEJKQUACIiIaUAEBEJKQWAiEhIKQBEREJKASAiElIKABGRkIo5AMysqZm9ambLzWypmV1ezjxmZveaWbGZLTKzQ6NeG2BmH0V+BiSrASIiUjnm7rHNaNYIaOTuC82sLrAA6OPuy6Lm6QkMBnoCnYF73L2zme0DzAcKAI8s28ndv05qa0REJGYx7wG4+1p3Xxh5/C2wHGhcZrbewONeag5QLxIcJwKz3X195Et/NlCYlBaIiEilVOoYgJnlAh2BuWVeagx8HvW8JDKtoukiIhKQ6vEuYGZ7AM8CV7j7f8u+XM4ivpPp5a2/CCgCqFOnTqe2bdvGW6KISKgtWLDgS3dvuKv54goAM6tB6Zf/BHefUs4sJUDTqOdNgDWR6ceUmf5aedtw97HAWICCggKfP39+PCWKiISemf0rlvniOQvIgEeA5e5+ZwWzTQXOiZwN1AX4xt3XAjOBHma2t5ntDfSITBMRkYDEswfQDTgbWGxm70em/QloBuDuo4EZlJ4BVAz8AJwXeW29md0EzIssd6O7r0+8fBERqayYA8Dd36L8vvzoeRy4rILXxgHj4qpORERSRlcCi4iElAJARCSkFAAiIiGlABARCSkFgIhISCkARERCKu6hIEQku3Qb+QqrN2wEoHG92rw9vHvAFUmmUACIVHGrN2zk05G/BiB3+PSAq5FMoi4gEZGQUgCIiISUAkBEJKQUACIiIaUAEBEJKQWAiEhIKQBEREJKASAiElIKABGRkFIAiIiElAJARCSk4hoLyMzGAScD69z9oHJevwo4M2rdBwINIzeF/xT4FtgGbHX3gkQKFxGRxMS7BzAeKKzoRXf/i7t3cPcOwB+B1919fdQsx0Ze15e/iEjA4goAd38DWL/LGUv1BybFXZGIiKRFSo4BmNnulO4pPBs12YFZZrbAzIpSsV0REYldqu4HcArwdpnun27uvsbM9gVmm9mHkT2Kn4mEQxFAs2bNUlSeiIik6iygfpTp/nH3NZF/1wHPAYeXt6C7j3X3AncvaNiwYYrKExGRpAeAme0FHA38PWpaHTOru+Mx0ANYkuxti4hI7OI9DXQScAzQwMxKgBFADQB3Hx2Z7TfALHf/PmrRXwHPmdmObU509xcTK11ERBIRVwC4e/8Y5hlP6emi0dNWAYfEsy0REUktXQksIhJSCgARkZBSAIiIhJQCQEQkpBQAIiIhpQAQEQkpBYCISEgpAEREQkoBICISUgoAEZGQUgCIiIRUqu4HICKSVt1GvsLqDRsBaFyvNm8P7x5wReXLpDoVACJSJazesJFPR/4agNzh0wOupmKZVKe6gEREQkoBICISUgoAEZGQUgCIiISUAkBEJKQUACIiIRVXAJjZODNbZ2ZLKnj9GDP7xszej/xcF/VaoZmtMLNiMxueaOEiIpKYePcAxgOFu5jnTXfvEPm5EcDMqgGjgJOAdkB/M2sXb7EiIpI8cQWAu78BrK/Edg4Hit19lbtvBiYDvSuxHhERSZJUHAPoamYfmNkLZtY+Mq0x8HnUPCWRaSIiEpBkDwWxEDjA3b8zs57A80ArwMqZ18tbgZkVAUUAzZo1S3J5IrIrmTRWTSapiu9LUvcA3P2/7v5d5PEMoIaZNaD0L/6mUbM2AdZUsI6x7l7g7gUNGzZMZnkiEoMdY9V8OvLXP33hSdV8X5IaAGa2n5lZ5PHhkfV/BcwDWplZnpnVBPoBU5O5bRERiU9cXUBmNgk4BmhgZiXACKAGgLuPBk4DLjGzrcBGoJ+7O7DVzAYBM4FqwDh3X5q0VoiISNziCgB377+L1+8H7q/gtRnAjHi2JyIiqaMrgUVEQkoBICISUgoAEZGQUgCIiISUAkBEJKQUACIiIZXsoSAkDSq6JL0qXKpeFdqwKztrY6rb37hebXKHT0/q+qNrLrutyq6/7DorWlfZ9ysI2fyZVQBkoR2XpAM//Wfe2fRsUhXasCs7a2Oq2x/95ZSs9UfXHC2R9ZddZ0Xrqmjb6ZTNn1l1AYmIhJQCQEQkpBQAIiIhpQAQEQkpBYCISEgpAEREQkoBICISUgoAEZGQUgCIiISUrgSOQzZf8l2eeNuTTUNQJKvWTGyb7FqyhrxI51ATQXzWFABxyOZLvssTb3uyaQiKZNWaiW2TXUvWkBfpHGoiiM9aXF1AZjbOzNaZ2ZIKXj/TzBZFft4xs0OiXvvUzBab2ftmNj/RwkVEJDHxHgMYDxTu5PVPgKPdPR+4CRhb5vVj3b2DuxfEuV0REUmyuLqA3P0NM8vdyevvRD2dAzSpXFkiIpJqqTwLaCDwQtRzB2aZ2QIzK0rhdkVEJAYpOQhsZsdSGgD/EzW5m7uvMbN9gdlm9qG7v1HOskVAEUCzZs1SUZ6IiJCCPQAzywceBnq7+1c7prv7msi/64DngMPLW97dx7p7gbsXNGzYMNnliYhIRFIDwMyaAVOAs919ZdT0OmZWd8djoAdQ7plEIiKSHnF1AZnZJOAYoIGZlQAjgBoA7j4auA6oDzxgZgBbI2f8/Ap4LjKtOjDR3V9MUhtERKQS4j0LqP8uXr8AuKCc6auAQ365hIiIBCXUVwJHX3odLYyX/Cfr0vkwStUl/PH+TjRshcQr1AFQ0WXeYbzkP1mXzodRqi7hj/d3omErJF4aDVREJKQUACIiIaUAEBEJKQWAiEhIKQBEREJKASAiElIKABGRkFIAiIiElAJARCSkQn0lcJDivWy/7PyplCnDQlRURyYPeRBdc9npmVRnotL5eUymdA7/ksmf0x0UAAGJ97L9ioatSIVMGRaiojoyeciDiv6TZ1qdiUrn5zGZ0jn8SyZ/TndQF5CISEgpAEREQkoBICISUgoAEZGQUgCIiISUAkBEJKTiCgAzG2dm68xsSQWvm5nda2bFZrbIzA6Nem2AmX0U+RmQaOEiIpKYePcAxgOFO3n9JKBV5KcIeBDAzPYBRgCdgcOBEWa2d7zFiohI8sQVAO7+BrB+J7P0Bh73UnOAembWCDgRmO3u6939a2A2Ow8SERFJsWRfCdwY+DzqeUlkWkXTf8HMiijde6BZs2ZJKWpnl3/LzpUdjiEbJasNqRgio+zQERXVlynDc6RSrMNLpGIYijC8v+VJdgBYOdN8J9N/OdF9LDAWoKCgoNx54pWtl61ngqrwHyFZbUjFEBmx1pYpw3OkUqz/T1Px/zkM7295kn0WUAnQNOp5E2DNTqaLiEhAkh0AU4FzImcDdQG+cfe1wEygh5ntHTn42yMyTUREAhJXF5CZTQKOARqYWQmlZ/bUAHD30cAMoCdQDPwAnBd5bb2Z3QTMi6zqRnff2cFkERFJsbgCwN377+J1By6r4LVxwLh4ticiIqmjK4FFREJKASAiElIKABGRkFIAiIiElAJARCSkdFP4XUj0svNMG4YiFZfRZ4p425bI5f+ZMkRGsn6f6fxcZNr/iWTKlM9FrBQAu5DoZeeZNgxFptWTTPG2LZHL/zNliIxk/T7T+bmoyp/BTPlcxEpdQCIiIaUAEBEJKQWAiEhIKQBEREJKASAiElIKABGRkFIAiIiElAJARCSkFAAiIiGlABARCSkNBZEEZcdRifdy8OjxQ3Y2T7apaMwX+Pn7lMg4NNk29koyxNLmTHhfEhlraVfrqiqCble89wQuBO4BqgEPu/vIMq/fBRwbebo7sK+714u8tg1YHHntM3fvlUjhmSR6bJN4x5SB7Bs/JFY7G/Ml+n1KZGyYqvre7Uwsbc6E9yWRsZZ2tq6qJOh2xRwAZlYNGAWcAJQA88xsqrsv2zGPu/8+av7BQMeoVWx09w6JlywiIskQzzGAw4Fid1/l7puByUDvnczfH5iUSHEiIpI68QRAY+DzqOclkWm/YGYHAHnAK1GTdzOz+WY2x8z6xF2piIgkVTzHAKycaV7BvP2AZ9x9W9S0Zu6+xsyaA6+Y2WJ3//gXGzErAooAmjVrFkd5IiISj3j2AEqAplHPmwBrKpi3H2W6f9x9TeTfVcBr/Pz4QPR8Y929wN0LGjZsGEd5IiISj3gCYB7QyszyzKwmpV/yU8vOZGZtgL2Bd6Om7W1mtSKPGwDdgGVllxURkfSJuQvI3bea2SBgJqWngY5z96VmdiMw3913hEF/YLK7R3cPHQiMMbPtlIbOyOizh0REJP3iug7A3WcAM8pMu67M8+vLWe4d4OBK1CciIimioSBEREJKQ0FUUtCXcMcjkaEWKpJo+zP5/UvF+xV2Qf6+4912tgy1kQwKgEoK+hLueCQy1EJFEm1/Jr9/qXi/wi7I33e8286WoTaSQV1AIiIhpQAQEQkpBYCISEgpAEREQkoBICISUgoAEZGQUgCIiISUAkBEJKQUACIiIaUrgctRFS7zrgptiFbV2pMJ9J6KAqAcVeEy76rQhmhVrT2ZQO+pqAtIRCSkFAAiIiGlABARCSkFgIhISCkARERCKq4AMLNCM1thZsVmNryc1881sy/M7P3IzwVRrw0ws48iPwOSUbyIiFRezKeBmlk1YBRwAlACzDOzqe6+rMysT7r7oDLL7gOMAAoABxZElv06oepFRKTS4tkDOBwodvdV7r4ZmAz0jnHZE4HZ7r4+8qU/GyiMr1QREUmmeAKgMfB51POSyLSyTjWzRWb2jJk1jXNZERFJk3iuBLZypnmZ59OASe6+ycwuBh4Duse4bOlGzIqAIoBmzZrFUd7PdRv5Cqs3bATSe5l79OX16dh2pl3On2n1xCpb6061su+Lrh6unEz9fMUTACVA06jnTYA10TO4+1dRTx8Cbota9pgyy75W3kbcfSwwFqCgoKDckIjF6g0b+XTkryu7eKWl+z9Ipv2HzLR6YpWtdada9PsS/YeNxCdTP1/xdAHNA1qZWZ6Z1QT6AVOjZzCzRlFPewHLI49nAj3MbG8z2xvoEZkmIiIBiXkPwN23mtkgSr+4qwHj3H2pmd0IzHf3qcAQM+sFbAXWA+dGll1vZjdRGiIAN7r7+iS2Q0RE4hTXaKDuPgOYUWbadVGP/wj8sYJlxwHjKlGjiIikgK4EFhEJKQWAiEhIKQBEREJKASAiElIKABGRkFIAiIiElG4KL5IGmToUQDyqQhvk5xQAImmQqUMBxKMqtEF+Tl1AIiIhpQAQEQkpBYCISEgpAEREQkoBICISUgoAEZGQUgCIiISUAkBEJKQUACIiIaUAEBEJKQWAiEhIxRUAZlZoZivMrNjMhpfz+h/MbJmZLTKzl83sgKjXtpnZ+5GfqckoXkREKi/mweDMrBowCjgBKAHmmdlUd18WNdt7QIG7/2BmlwC3A30jr2109w5JqltERBIUzx7A4UCxu69y983AZKB39Azu/qq7/xB5OgdokpwyRUQk2eIJgMbA51HPSyLTKjIQeCHq+W5mNt/M5phZnzi2KyIiKRDP/QCsnGle7oxmZwEFwNFRk5u5+xozaw68YmaL3f3jcpYtAooAmjVrFkd5IiISj3j2AEqAplHPmwBrys5kZscDVwO93H3Tjunuviby7yrgNaBjeRtx97HuXuDuBQ0bNoyjPBERiUc8ATAPaGVmeWZWE+gH/OxsHjPrCIyh9Mt/XdT0vc2sVuRxA6AbEH3wWERE0izmLiB332pmg4CZQDVgnLsvNbMbgfnuPhX4C7AH8LSZAXzm7r2AA4ExZrad0tAZWebsIRERSbO47gns7jOAGWWmXRf1+PgKlnsHOLgyBYqISGroSmARkZBSAIiIhJQCQEQkpBQAIiIhpQAQEQkpBYCISEgpAEREQkoBICISUgoAEZGQUgCIiISUAkBEJKQUACIiIaUAEBEJKQWAiEhIKQBEREJKASAiElIKABGRkFIAiIiElAJARCSk4goAMys0sxVmVmxmw8t5vZaZPRl5fa6Z5Ua99sfI9BVmdmLipYuISCJiDgAzqwaMAk4C2gH9zaxdmdkGAl+7e0vgLuC2yLLtgH5Ae6AQeCCyPhERCUg8ewCHA8XuvsrdNwOTgd5l5ukNPBZ5/AxwnJlZZPpkd9/k7p8AxZH1iYhIQOIJgMbA51HPSyLTyp3H3bcC3wD1Y1xWRETSqHoc81o50zzGeWJZtnQFZkVAUeTpd2a2IuYKf66B3caXlVw2WzUAtTkE1OYQsNsSavMBscwUTwCUAE2jnjcB1lQwT4mZVQf2AtbHuCwA7j4WGBtHXeUys/nuXpDoerKJ2hwOanM4pKPN8XQBzQNamVmemdWk9KDu1DLzTAUGRB6fBrzi7h6Z3i9yllAe0Ar4Z2Kli4hIImLeA3D3rWY2CJgJVAPGuftSM7sRmO/uU4FHgCfMrJjSv/z7RZZdamZPAcuArcBl7r4tyW0REZE4xNMFhLvPAGaUmXZd1OMfgd9VsOyfgT9XosbKSrgbKQupzeGgNodDyttspT00IiISNhoKQkQkpLI+ABIZniJbxdDmP5jZMjNbZGYvm1lMp4Rlsl21OWq+08zMzSyrzxiJpb1mdnrk97zUzCamu8Zki+Fz3czMXjWz9yKf7Z5B1JlMZjbOzNaZ2ZIKXjczuzfyniwys0OTWoC7Z+0PpQejPwaaAzWBD4B2Zea5FBgdedwPeDLoutPQ5mOB3SOPLwlDmyPz1QXeAOYABUHXneLfcSvgPWDvyPN9g647DW0eC1wSedwO+DToupPQ7qOAQ4ElFbzeE3iB0mupugBzk7n9bN8DSGR4imy1yza7+6vu/kPk6RxKr7vIZrH8ngFuAm4HfkxncSkQS3svBEa5+9cA7r4uzTUmWyxtdmDPyOO9qOBaomzi7m9QesZkRXoDj3upOUA9M2uUrO1newAkMjxFtop3WI2BlP4Fkc122WYz6wg0dfd/pLOwFInld9waaG1mb5vZHDMrTFt1qRFLm68HzjKzEkrPRhycntICldJhdOI6DTQDJTI8RbaKZ1iNs4AC4OiUVpR6O22zmeVQOvrsuekqKMVi+R1Xp7Qb6BhK9/DeNLOD3H1DimtLlVja3B8Y7+53mFlXSq85Osjdt6e+vMCk9Psr2/cA4hmegjLDU2SrmIbVMLPjgauBXu6+KU21pcqu2lwXOAh4zcw+pbSvdGoWHwiO9XP9d3ff4qUj7K6gNBCyVSxtHgg8BeDu7wK7UTpGUFUW8zA6lZHtAZDI8BTZapdtjnSHjKH0yz/b+4ZhF21292/cvYG757p7LqXHPXq5+/xgyk1YLJ/r5yk92I+ZNaC0S2hVWqtMrlja/BlwHICZHUhpAHyR1irTbypwTuRsoC7AN+6+Nlkrz+ouIE9geIpsFWOb/wLsATwdOd79mbv3CqzoBMXY5iojxvbOBHqY2TJgG3CVu38VXNWJibHNVwIPmdnvKe0GOTfL/5jDzCZR2o3XIHJsYwRQA8DdR1N6rKMnpfdQ+QE4L6nbz/L3T0REKinbu4BERKSSFAAiIiGlABARCSkFgIhISCkARERCSgEgIhJSCgARkZBSAIgkiZl9F3QNIvFQAIiIhJQCQELHzCZH7hI318z+ZWa/LvP6bWZ2adTz683sysjj581sQeQuXEXlrDs3+u5OZjbUzK6Pen6Wmf3TzN43szFmVi0ljRSJgQJAwugQYJW7dwbOpHT8lWiTgb5Rz08Hno48Pt/dO1E6zPYQM4v53hKRAcz6At3cvQOlY/icWbkmiCQuqweDE4mXmdWmdAjhGyKTlgF7R8/j7u+Z2b5mtj/QEPja3T+LvDzEzH4TedyU0iGYYx2E7TigEzAvMkhfbaAqjNYqWUoBIGFzEPCRu++4beShwAdmdhmlt1mE0tEXn6F0+PD9KN0jwMyOAY4Hurr7D2b2GqVDEkfbys/3rKNfN+Axd/9j0lojkgB1AUnYHAI0M7PdzKwOpXsCd7n7KHfvEPlZQ+mXfj9KQ+CZyLJ7Ubo38IOZtaX0xjNl/QfY18zqm1kt4OSo114GTjOzfQHMbB8zOyAlrRSJgfYAJGwOASYAr1F6g/Fb3P3tsjNFxqKvC6yOugHHi8DFZraI0jtwzSlnuS2RMeznAp8AH0a9tszMrgFmRW5juQW4DPhXEtsnEjPdD0BCxczeAC509xVB1yISNAWAhIqZrQaaVvEbiYvERAEgIhJSOggsIhJSCgARkZBSAIiIhJQCQEQkpBQAIiIhpQAQEQkpBYCISEgpAEREQur/AJ4Pk9MLoV5iAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x648 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "tmus = ht.tmu_rvs(size = 1000);\n",
    "htplot.chi2_test(tmus);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: Compute the CI at 68 % CL using $t_\\mu(x)$ and the p0-value for the case $\\nu_g = 0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "central CI  : [10.30176774 25.32851881] at  90 % CL\n",
      "p0 value    : 2.102124030400887e-10\n"
     ]
    }
   ],
   "source": [
    "beta = 0.90\n",
    "ci = ht.tmu_cint(x0, parbest = parbest, beta = beta)\n",
    "#ul = ht.qmu_ulim(x0, parbest = parbest, beta = beta)\n",
    "q0 = ht.q0(x0, 0., parbest = parbest)\n",
    "p0 = htcore.q0_pvalue(q0)\n",
    "htplot.prt_wilks_pvalues(beta, ci_central = ci, p0 = p0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEMCAYAAAAxoErWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzt3Xl8VNX5+PHPM5ONLCxZSQj7nrCHRVyDLOKC2AoVbVVat1Zta9Xvz6X9ulWrtdpqv2qrdpG6FC1uiCiyxZUd2ZMQQJZAyAYhCVknc35/ZJKGkJCFZO4sz/v1mhdzZ86d+xxuMk/uOeeeI8YYlFJKqabYrA5AKaWU59IkoZRSqlmaJJRSSjVLk4RSSqlmaZJQSinVLE0SSimlmqVJQimlVLM0SSillGqWJgmllFLNCrA6gLMVHR1t+vXr1659T548SVhYWMcGZBGti+fxlXqA1sVTnU1dNm3aVGCMiWmpnNcniX79+rFx48Z27ZuWlkZqamrHBmQRrYvn8ZV6gNbFU51NXUTkQGvKaXOTUkqpZmmSUEop1Sy3JQkR+YeI5InIjmbeFxH5s4jsEZFtIjLOXbEppZRqmjv7JF4DXgD+1cz7lwKDXY9JwF9c/yql/EB1dTXZ2dlUVFR06nG6detGenp6px7DXVpTl5CQEBITEwkMDGzXMdyWJIwxX4hIvzMUmQ38y9QucLFWRLqLSLwxJsctASqlLJWdnU1ERAT9+vVDRDrtOCUlJURERHTa57tTS3UxxlBYWEh2djb9+/dv1zE8qU+iF3CowXa26zWllB+oqKggKiqqUxOEvxERoqKizurqzJOGwDb1k9HksnkicitwK0BcXBxpaWntOmBpaWm79/U0WhfP4yv1APfUpVu3bpSWlnbqMQBqamooKSnp9OO4Q2vrUlFR0e7z50lJIhvo3WA7ETjSVEFjzCvAKwDjx4837R0nrOOlPVNH1GV/wUkyjhYzc0R8xwTVDnpO2iY9Pd0tzUD+1NxUJyQkhLFjx7brGJ7U3LQYuME1yukc4IT2R6j2+sOyTH76xma2ZRdZHYpSXs2dQ2D/DawBhopItojcJCI/FZGfuoosBfYBe4BXgdvdFZvyLY4aJ19k5QPw+Mfp1I6FUMpa4eHhHfI52dnZvP322x3yWa3hztFN17bwvgHucFM4yodtPlhESYWDC4fE8MXufD7blcslyT2tDkupDrFy5Up27drFNddc45bjeVJzk1IdIi0zjwCb8Pw1YxgUG85Tn2RQ5XBaHZbyEvPmzeOaa65h0qRJ9O3bl48//viU9++77z5eeuml+u1HHnmEZ599FoCrrrqKlJQUkpOTeeWVV0777P379zNixIj67WeeeYZHHnmkfvuNN95g4sSJjBkzhttuu42amppT9v/qq6+4++67WbRoEWPGjGH//v0dUOMz86SOa6U6xOrMfFL69qBHWBC/vmw4P35tA2+uO8CPz2vfOHHlfo9+tJNdR4o79DOTErry8KzkFstt3bqVq666irfffrv+S/nyyy+vf3/evHncdddd3H57bYv4O++8w6effgrAP/7xDyIjIykvL2fChAlcffXVREVFtSq+9PR03n77bb7++msCAwO5/fbbefPNN7nhhhvqy5x//vlMmDCBZ555hhEjRrhllJYmCeVTjp6oID2nmPtmDgMgdWgM5w+K5vmVWVw3qQ/BAXaLI1SerLy8nIKCAh5++GEAkpKSOH78+Cllxo4dS15eHkeOHCE/P58ePXrQp08fAP785z/z/vvvA3Do0CGysrJanSRWrlzJpk2bmDBhQn0ssbGxp5XLzMxk6NCh7a5jW2mSUD7l8915AEwZVjtNvojwk/P78ZPXNrL+u2NcMLjF6fOVB2jNX/ydYceOHQwePJiQkBAANm/ezOjRo3nxxRd59dVXAVi6dClz5sxh0aJFHD16lHnz5gG1w4RXrFjBmjVrCA0NJTU19bSb2AICAnA6/9v02fB9Yww33ngjTz75ZLPxFRYW0q1bt/opNtLT01m9ejX33Xcfd955J08++WSHD+/VPgnlU9Iy8+nZNYShcf/9RZk8IJrgABsr0/MsjEx5g61bt3Lw4EEqKio4efIkDz/8ML/61a+444472LJlC1u2bCEhIYF58+axcOFCFi1axJw5cwA4ceIEPXr0IDQ0lIyMDNauXXva58fFxZGXl0dhYSGVlZUsWbKk/r2pU6eyaNEi8vJqf06PHTvGgQOnLvnw3XffkZCQUL9dl8QAiouLO+X+D00SymdU1zj5KquAKcNiTpnaoUuQnXMHRrE6M0+Hw6oz2rp1Kz/84Q9JTU1lwoQJ/OxnP+O88847rVxycjIlJSX06tWL+PjaGzZnzpyJw+Fg1KhR/O///i/nnHPOafsFBgby0EMPMWnSJK644gqGDRtW/15SUhKPP/44M2bMYNSoUUyfPp2cnFNvFRs2bBgFBQWMGDGCb775hs2bN5OUlMTJkyc7+H/iv7S5SfmMTQeOU1Lp4KIhp7fjXjwsltUf7mRfwUkGxnTMeHXle7Zu3cqrr77K73//+xbLbt++/ZTt4OBgPvnkk2bL10058otf/IJf/OIXTZa55pprzji0NTw8nPXr19dvP/roozzxxBN07dq1/gqko2mSUD5jdWYegXbhvEGndxROGRYLH+5kdUaeJgnVrL179zJ48GCrw2iV6upqIiMjefnllzv1ONrcpHzG55n5jO8bSUTI6fPmJ/YIZUhcOKsytF9CNe/w4cPYbN7xtRgYGNjkvRgdzTv+N5RqwZGicjKOlpA6tPnRSxcPi2P9d8coqah2Y2RKeTdNEsonfL67dq6mKcNO74+oc/GwWBxOw1dZBe4KSymvp0lC+YS0zDwSuoUwOLb5/oZxfbrTrUsgK7XJSalW0yShvF6Vw8nXewq5aGjsGVc1C7DbuHBIDGmZeTidOhRWqdbQJKG83qYDxymtdJyxP6LO1GGxFJRWse3wCTdEppT30yShvF7a7rqhr9Etlr1oSAw2gVXpuW6ITCnvp0lCeb26oa/hwS3f9tMjLIiUvj20X0KpVtIkobxaa4a+NjZ1eBw7jxSTc6K8EyNT3uj5559nxIgRJCcn89xzz9W/fuzYMaZPn87gwYOZPn16/cyw7777LsnJyVxwwQUUFhYCtTfk1U3615TS0lJuu+02Bg4cSHJyMhdeeCHr1q0DOm71uo6kSUJ5tdYMfW1sqqus3linGtqxYwevvvoq69evZ+vWrSxZsoSsrCwAnnrqKaZOnUpWVhZTp07lqaeeAuDZZ59l7dq13HDDDbz11lsA/OY3v+G3v/1ts8e5+eabiYyMJCsri507d/Laa69RUOC5w7I1SSiv1pqhr40Nig2nT2SozgqrTpGens4555xDaGgoAQEBXHTRRfVrQ3z44YfceOONANx444188MEHANhsNiorKykrKyMwMJAvv/yS+Pj4Zqf22Lt3L+vWrePxxx+vv7N7wIABpyxq5Gl07iblteqGvs4anXDGoa+NiQgXD4vl3+sPUl5VQ5cgXYjIE6WmprZY5oorruDee++tLz9//nzmz59PQUFB/RTeddLS0s74WSNGjODXv/41hYWFdOnShaVLlzJ+/HgAcnNz62d7jY+Pr59M7+GHH+aSSy4hISGBN954gx/84AcsXLiw2WPs3LmTMWPGYLd7z8+cXkkor9WWoa+NTRseR6XDydd7PPcyX7nX8OHDue+++5g+fTozZ85k9OjRBASc+e/o6dOns2nTJj766CM++OADLrvsMjIzM5kzZw633HILZWVlboq+8+iVhPJabRn62tjE/pGEBdlZmZHLtKS4TohOna2W/vI/U/no6Og27w9w0003cdNNNwHw4IMPkpiYCNQuFpSTk0N8fDw5OTmnLStaVlbGggULWLZsGTNmzODDDz/krbfe4s033+SWW26pL5ecnMzWrVtxOp1eM5Ggd0SpVBPSMvKZ0K91Q18bCwqovft6ZbouRKT+q64Z6eDBg7z33ntce+21AFx55ZUsWLAAgAULFjB79uxT9nv66af55S9/SWBgIOXl5YgINpvttCuJgQMHMn78eB5++OH6n7usrCw+/PDDzq5au2mSUF7pSFE5mbltG/ra2NThceSVVLJd775WLldffTVJSUnMmjWLF198kR49egBw//33s3z5cgYPHszy5cu5//776/c5cuQIGzdurE8c99xzD+eccw4LFizguuuuO+0Yf/vb3zh69CiDBg1i5MiR3HLLLfVLkpaVlZGYmFj/+OMf/+iGWp+ZNjcpr1Q39DV1aOuHvjY2ZWjt3dcr0vMYldi9o0JTXuzLL79s8vWoqChWrlzZ5HsJCQmnrFU9d+5c5s6d2+wxunbtyquvvtrke06nsw3RuodeSSiv1J6hr41FhQeT0rcHK3bpFB1KNUeThPI6rZ31tTWmDY9jV04xh4v07mulmqJJQnmdsxn62ljdyKaVOuGfUk3SJKG8ztkMfW1sYEw4A6LDWK5NTh5BR5p1vLP9P9UkobzO2Qx9bcq0pDjW7ivUta8tFhISQmFhoSaKDmSMobCwkJCQkHZ/ho5uUl6lbujrgynDOuwzpw2P45Uv9vHF7gIuHxXfYZ+r2iYxMZHs7Gzy8/M79TgVFRVn9aXpSVpTl5CQkPqbAtvDrUlCRGYCzwN24G/GmKcavd8HWAB0d5W53xiz1J0xKs+Wluma9fUshr42Nq5Pd3qEBrIiPVeThIUCAwPp379/px8nLS2NsWPHdvpx3MEddXFbc5OI2IEXgUuBJOBaEUlqVOw3wDvGmLHAPOAld8WnvENaZh69undh0FkMfW0swG5jyrBYVmXk4ajxvHHqSlnJnX0SE4E9xph9xpgqYCEwu1EZA3R1Pe8GHHFjfMrDVbkm5EsdGnPWQ18bmz48jhPl1WzYf7xDP1cpb+fOJNELONRgO9v1WkOPAD8SkWxgKfBz94SmvMHG/cc4WVVzVndZN+fCITEEBdj4bNfRDv9spbyZuGskgYjMBS4xxtzs2r4emGiM+XmDMne7YnpWRCYDfwdGGGOcjT7rVuBWgLi4uJQzzd9+JqWlpR65XGB7+ENdFmZUsuKAgxemhhIS0LFXEgB/2lRBdomTZy7q0iFXKv5wTryR1qXWlClTNhljxrdY0BjjlgcwGVjWYPsB4IFGZXYCvRts7wNiz/S5KSkppr1Wr17d7n09jT/UZdqzaeaHr67ttOMuXH/A9L1vidlxuKhDPs8fzok30rrUAjaaVnx3u7O5aQMwWET6i0gQtR3TixuVOQhMBRCR4UAI0Lnj4ZRXyD5eRlZeaYfcZd2cqcPjsAl8tlNvrFOqjtuShDHGAdwJLAPSqR3FtFNEHhORK13F7gFuEZGtwL+B+a6Mp/xc3dDXzuiPqBMdHsz4vpF8pndfK1XPrfdJmNp7HpY2eu2hBs93Aee5MyblHdIy8+gd2YWBMWGdepwZyXE8/nE6h46V0TsytFOPpZQ30Gk5lMerdNTw9Z5CUoec/ayvLZnumvBv2U4d5aQUaJJQXmDdvmOUV9dw8bDOa2qq0zcqjGE9I7TJSSkXTRLK463OzCM4wMY5A6LccrwZSXFs3H+MwtJKtxxPKU+mSUJ5vLTMfCYPjKJLkN0tx5uR3BOngZXpeW45nlKeTJOE8mjfFZzku4KTHTqhX0uSE7qS2KMLn2q/hFKaJJRnS8us/WvenUlCRJiZ3JOvsgp0jQnl9zRJKI+2KiOPATFh9Ily73DUmSN6UlXjZFWGNjkp/6ZJQnmssioH6/Ydc+tVRJ1xfXoQExGsQ2GV39MkoTzWN3sKqapxumXoa2M2m3BJchyrM/Ipr6px+/GV8hSaJJTHWp2ZR1iQnfH9elhy/JnJ8ZRX1/BFlk4fpvyXJgnlkYwxpGXmc+6gaIID3DP0tbFJAyLpHhrIpzu0yUn5L00SyiNl5pZwuKjckqamOoF2G9OGx7EiPZcqhy5rqvyTJgnlkepGFVnRad3QpSN6UlLh4Ju9BZbGoZRVNEkoj7Q6I4/khK707BZiaRznDYomPDiAT7Zrk5PyT5oklMcprTJsOnDc0qamOiGBdqYNj2XZrqNU12iTk/I/miSUx9leUIPT4BFJAuCykfEUlVXzzd5Cq0NRyu00SSiPszXfQVRYEKMTu1sdCgAXDokhPDiApdtyrA5FKbfTJKE8iqPGyfaCGi4aGoPN1rkLDLWWNjkpf6ZJQnmUbw8VcbIapg6LszqUU9Q1Oa3RJiflZzRJKI+yKiMPu8AFQ6KtDuUUdU1OH2uTk/IzmiSUR1mdkcfgHja6hgRaHcopQgLtTNUmJ+WHNEkoj5F9vIyMoyWMjgmwOpQmXa6jnJQf0iShPEbdXdZjYqyZq6kl/21yOmJ1KEq5jSYJ5TFWpufRPzqM+HDP/LEMCbQzIzmOT3ccpdKh04cr/+CZv43K75ysdLBmbyFTPeQGuubMGp1AcYWDL3brXE7KP2iSUB7hy6wCqmqcTB3uWUNfGzt/UDQ9QgP5aKs2OSn/oElCeYRVGblEhARYtsBQawXabVw6Mp7lu3Ipq3JYHY5SnU6ThLKc02lYlZFP6tBYAu2e/yM5a1QC5dU19R3tSvkyz/+NVD5va3YRBaWVTBvu2f0RdSb2jyQ2IpjFW7TJSfk+TRLKcqsy8rDbhIuGxFgdSqvYbcIVoxJIy8ynuKLa6nCU6lSaJJTlVqTnkdK3B91Dg6wOpdVmjY6nqsbJMl3/Wvk4tyYJEZkpIpkiskdE7m+mzA9EZJeI7BSRt9wZn3K/w0XlpOcUe/zQ18bG9O5On8hQFusoJ+Xj3JYkRMQOvAhcCiQB14pIUqMyg4EHgPOMMcnAXe6KT1ljxa5cAKYnefbQ18ZEhKvGJPD1ngJyiyusDkepTuPOK4mJwB5jzD5jTBWwEJjdqMwtwIvGmOMAxhgdPuLjVqTnMiAmjAEx4VaH0mazx/bCadB7JpRPc2eS6AUcarCd7XqtoSHAEBH5WkTWishMt0Wn3K64opq1+wq97iqizsCYcEYnduP9bw9bHYpSncad0202tcyYabQdAAwGUoFE4EsRGWGMKTrlg0RuBW4FiIuLIy0trV0BlZaWtntfT+ONdVmb46C6xhBdcYS0tNz6172pLiPCq3kzo4o3l6yiV6M5p7ypHi3Rungmt9TFGOOWBzAZWNZg+wHggUZl/grMb7C9Ephwps9NSUkx7bV69ep27+tpvLEuP39rsxn32GfGUeM85XVvqktecYUZ8MDH5ulP0097z5vq0RKti2c6m7oAG00rvrvd2dy0ARgsIv1FJAiYByxuVOYDYAqAiERT2/y0z40xKjeprnGyOjOPqcNjsXvIWtbtERMRzAWDo/ng2yM4nY0vjJXyfm5LEsYYB3AnsAxIB94xxuwUkcdE5EpXsWVAoYjsAlYD/2OM0RVefNC6fccoqXAwzcMn9GuNq8b04nBRORsPHLc6FKU6XKv7JEREgB8CA4wxj4lIH6CnMWZ9az/DGLMUWNrotYcaPDfA3a6H8mEr0nMJCbRxwWDvuMv6TGYkxxEaZOf9b7OZ2D/S6nCU6lBtuZJ4idp+hWtd2yXU3vegVJsYY1i+K5fzB8XQJcgzV6Fri9CgAC4dEc+SrTlUVOtiRMq3tCVJTDLG3AFUAJjaexm8Zx4F5TF25RRzuKic6UnedZf1mcxJSaSk0sGynTpNh/ItbUkS1a67pg2AiMQAzk6JSvm0ZTtzsQkev8BQW0zqH0lijy4s2pRtdShKdai2JIk/A+8DsSLyBPAV8LtOiUr5tM92HmV8v0iiw4OtDqXD2GzC1eMS+WpPATknyq0OR6kO0+okYYx5E/h/wJNADnCVMeY/nRWY8k0HCk+ScbSES5J7Wh1Kh7t6XCLGwHub9Q5s5TvadMe1MSYDyOikWJQfqGuzn+GlU3GcSZ+oUCb2j+TdTdncnjrQ6nCU6hBtGQL7UFOvG2Me67hwlK9btjOX5ISu9I4MtTqUTjEnJZH/t2gbmw8WtVxYKS/Qlj6Jkw0eNdRO+d2vE2JSPiqvpILNB4/7ZFNTnctGxtMl0M6iTYdaLqyUF2j1lYQx5tmG2yLyDKdPq6FUs5bvysUYfDpJhAcHcPmoeBZvOcJFF/pOx7zyX2czLUcoMKCjAlG+b9nOXPpFhTIkzvvWjmiLeRN6c7KqhvVHHVaHotRZa0ufxHb+O7W3HYgBftsZQSnfU1xRzZq9BfzkvP7UzvDiu1L69mBQbDhfZJdZHYpSZ60to5uu4L9JwgHkGmN0DgLVKivTc6muMczw4aamOiLCvAm9efzjdHbnljAkLsLqkJRqtxabm0SkRESKge3ADtcjAzjuel2pFi3dfpT4biGM7d3d6lDc4vvjErELLFyvHdjKu7WYJIwxEcaYrs093BGk8m6llQ4+353PzBE9sXnx2hFtERkWREqcnfe+zdZJ/5RXc+eiQ8pPrcrIo8rh5LKR8VaH4lYXJQZSVFbNZ7tyWy6slIdqU5IQkR4iMlFELqx7dFZgynd8sj2H2IhgUvr0sDoUtxoeZaN3ZBf+ve6g1aEo1W6tThIicjPwBbWrxz3q+veRzglL+YqyKgerM/P8qqmpjk2EeRP6sGZfIXvySqwOR6l2acuVxC+BCcABY8wUYCyQ3ylRKZ+RlplPRbWTS0f4V1NTnWsm9CbQLryxVq8mlHdqS5KoMMZUAIhIsGuyv6GdE5byFUu35xAVFuS3y3pGhwdz6Yh43t2cTVmV3lynvE9bkkS2iHQHPgCWi8iHwJHOCUv5gorqGlZl5HHJiJ7Y/aypqaHrJ/elpMLB4i3666K8T1vmbvqe6+kjIrIa6AZ82ilRKZ+QlplPWVUNl47w/RvozmR83x4MjYvg9bUHuGZCb5+/41z5lrZ0XP9KRBIBjDGfG2MWG2OqOi805e2WbDtCVFgQkwdEWR2KpUSEH03uy84jxWw5pFOIK+/SluamrsAyEflSRO4QEd9bNUZ1mLIqByvT87h0ZE8C7Ho7zvfG9iIsyM7raw9YHYpSbdKW5UsfNcYkA3cACcDnIrKi0yJTXm1leh7l1TVcMSrB6lA8QnhwAN8fl8iSrTkUlFZaHY5SrdaeP/HygKNAIRDbseEoX/HR1iPEdQ1mQj//HNXUlBvP7UdVjZO39OY65UXa0ifxMxFJA1YC0cAtxphRnRWY8l7FFdWk7c7nspHxfj2qqbFBseFcNCSG19ceoMrhtDocpVqlLVcSfYG7XE1OfzHG7OqkmJSXW74zlyqHk1mjtampsZ+c35/8kko+3q7DYZV3aEufxP3GmC2uzaWdFI/yAUu2HaFX9y5+My14W1w4OJqBMWH88+v9GGNa3kEpi7V32Im2IagmHT9ZxZdZBVwxOl7vB2iCiDD/vP5syz7B5oPHrQ5HqRa1N0m82qFRKJ/x6c6jOJyGWTqqqVlXj+tF15AA/vHVfqtDUapFbem4vlRE1olIJpAqIpM7MS7lpd7/9jADY8JITtD1qJoTGhTAtZP68OnOoxw6putgK8/WliuJl4C7gXOAl4E/iMi1bTmYiMwUkUwR2SMi95+h3BwRMSIyvi2fr6yVfbyM9d8d43tje2lTUwt+fG5/bAJ//+o7q0NR6ozakiRyjTFfG2OOG2NWApcAv27tziJiB14ELgWSgGtFJKmJchHAL4B1bYhNeYAPXRPYzR7Ty+JIPF/PbiFcNaYXCzcc5NhJnd1Gea62JIn9IvK4iAS5tquBtqykMhHYY4zZ55rzaSEwu4lyvwWeBira8NnKYsYYPvj2MBP69aB3ZKjV4XiFWy8cQEW1k9fX6FQdynO1JUkY4PvAIRH5CtgDpInI4Fbu3ws41GA72/VaPREZC/Q2xixpQ1zKA+zKKSYrr1SvItpgcFwEU4fFsmDNfiqqa6wOR6kmSVvHaotICDACGN3gMcAY07uF/eYClxhjbnZtXw9MNMb83LVtA1YB840x+113d99rjNnYxGfdCtwKEBcXl7Jw4cI21aFOaWkp4eHh7drX01hdl4UZlSw/4OD5KaGEB51df4TVdekoralH5rEanlxfwQ1JQVzcJ9BNkbWdr5wT0LrUmTJlyiZjTMv9vsYYtzyAycCyBtsPAA802O4GFAD7XY8Kahc1Gn+mz01JSTHttXr16nbv62msrIujxmkmPL7c3LxgQ4d8nq+cl9bUw+l0mtkvfGUu+P0qU+2o6fyg2slXzokxWpc6wEbTiu9ud87hvAEYLCL9Xf0a84DFdW8aY04YY6KNMf2MMf2AtcCVpokrCeVZ1uwtJK+kku+N1aamthIRfpY6kIPHyvhom07VoTyP25KEMcYB3AksA9KBd4wxO0XkMRG50l1xqI737uZsIkICuHiYTgrcHtOHxzGsZwQvrNqD06lTdSjP4tbVYIwxS40xQ4wxA40xT7hee8gYs7iJsql6FeH5iiuq+WRHDleOTiAk0G51OF7JZhPumDKIvfkn+WTHUavDUeoUumSYOisfb8uhotrJ3PFnHLegWnDZyHgGxITxf6uy9GpCeRRNEuqs/GfjIQbFhjM6sZvVoXg1u024c8ogMo6WsCI91+pwlKqnSUK12568UjYfLGJuSqJOw9EBrhydQJ/IUP5v1R6dRlx5DE0Sqt0WbcrGbhO+N05HNXWEALuNO6YMZPvhE6xMz7M6HKUATRKqnRw1Tt7bnE3qkBhiI0KsDsdnfH9cIv2iQnnms0ztm1AeQZOEapcvswrIK6lk7vhEq0PxKYF2G3dNG0LG0RI+3p5jdThKaZJQ7bNww0Eiw4K4eFic1aH4nFmjExgSF86flu/GUeO0Ohzl5zRJqDbLLa5gRXoec1MSCQrQH6GOZrcJd08fyr6Ck7z/7WGrw1F+Tn/DVZu9s+EQNU7DtRP7WB2Kz7okOY6Rvbrx/MosKh06Q6yyjiYJ1SY1TsPCDYc4b1AU/aLDrA7HZ4kI914ylOzj5byx9qDV4Sg/pklCtckXWfkcLirnuol9rQ7F5104OJoLBkfzf6uyOFFWbXU4yk9pklBt8ta6g0SHBzE9STusO5uI8MClwzlRXs2LaXusDkf5KU0SqtWOnqhgVUYec8f31g5rN0lK6MrV4xJ57ev9HDpWZnU4yg/pb7pqtbeIEIKHAAAWL0lEQVTrOqwnaIe1O90zYwg2G/xhWabVoSg/pElCtUqVw8mb6w5w0ZAY+kSFWh2OX4nv1oWbzx/A4q1H2HzwuNXhKD+jSUK1yic7csgrqWT+ef2sDsUv/TR1ILERwTyyeKdO16HcSpOEapXXvtnPgOgwLhocY3Uofik8OIAHLxvOtuwTvLPxkNXhKD+iSUK1aMuhIr49WMSN5/bDZtMpwa0ye0wCE/r14OllmTokVrmNJgnVogXf7Cc8OICrU3QyPyuJCI9eOYKisir+uFw7sZV7aJJQZ5RXUsGSbUeYOz6R8OAAq8Pxe0kJXfnROX15fe0Bdh45YXU4yg9oklBn9Obagzichhsn97M6FOVyz/ShRIYF8eB726nRTmzVyTRJqGaVV9Xw+toDXDw0Vudp8iDdQgN5aFYyW7NPsOCb/VaHo3ycJgnVrP9sOsSxk1X8NHWg1aGoRmaNiid1aAzPfJbJ4aJyq8NRPkyThGqSo8bJK1/sY1yf7ozv28PqcFQjIsJvZ4/AGPjfD3ZgjDY7qc6hSUI1aemOo2QfL+enFw1ERIe9eqLekaHcM2MIqzLyWLz1iNXhKB+lSUKdxhjDX9P2MjAmjGnDdbZXTzb/3H6M6d2dhz7cSV5xhdXhKB+kSUKd5qs9BezKKea2CwfqzXMeLsBu49kfjKbSUcN9727TZifV4TRJqNO8tHovcV2DmT02wepQVCsMjAnnvpnDWJ2Zr1N2qA6nSUKdYt2+QtbsK+SWCwYQHGC3OhzVSjdO7sfkAVE89tEuXXdCdShNEuoUf1qxm5iIYH50ji5P6k1sNuEPc0dhE+GXC7+lusZpdUjKR2iSUPXW7C1k7b5j/OyigYQE6lWEt0nsEcrvvj+SzQeLeG7FbqvDUT7CrUlCRGaKSKaI7BGR+5t4/24R2SUi20RkpYjon7Nu9NyK3cRGBHPdJF15zlvNGp3AvAm9eSltL19lFVgdjvIBbksSImIHXgQuBZKAa0UkqVGxb4HxxphRwCLgaXfF5+++2VvAuu+O8bNUvYrwdg/PSmZgTDh3vb2F/JJKq8NRXs6dVxITgT3GmH3GmCpgITC7YQFjzGpjTF2v21pA56Z2A2MMf1q+m7iuwVw7Ua8ivF2XIDsvXjeOkopqfv7vzTi0f0KdBXcmiV5Aw/F52a7XmnMT8EmnRqQAWL4rlw37j/PziwfrVYSPGNozgie/P5K1+47x1CcZVoejvJi46+YbEZkLXGKMudm1fT0w0Rjz8ybK/gi4E7jIGHPa9bKI3ArcChAXF5eycOHCdsVUWlpKeHh4u/b1NO2tS43T8OuvayeIe+K8Ltg94OY5XzkvnlCPN3ZVsuKgg9tGBTM5of3rgXhCXTqK1qXWlClTNhljxrdY0BjjlgcwGVjWYPsB4IEmyk0D0oHY1nxuSkqKaa/Vq1e3e19P0966vL5mv+l73xKzbEdOxwZ0FnzlvHhCPaocNWbuX74xQ3+z1Ow4XNTuz/GEunQUrUstYKNpxXesO5ubNgCDRaS/iAQB84DFDQuIyFjgZeBKY0yeG2PzS6WVDp5bsZuJ/SKZnqRzNPmiQLuNF384ju5dgrh5wUZydX4n1UZuSxLGGAe1TUjLqL1SeMcYs1NEHhORK13F/gCEA/8RkS0isriZj1Md4JUv9lFQWsWDlw/XmV59WExEMH+fP54T5dXctGADJysdVoekvIhbFy02xiwFljZ67aEGz6e5Mx5/ln28jFe+2MsVo+IZ07u71eGoTpac0I0XrhvLzQs28suF3/Ly9eM9ov9JeT6949pPPfrRLgThwcuGWx2KcpOLh8XxyJXJrEjP49GPduqMsapV3HoloTzDqoxclu/K5f5Lh5HQvYvV4Sg3umFyP7KPl/PKF/vo3iWQu2cMtTok5eE0SfiZiuoaHl68k0Gx4fzkvP5Wh6Ms8MClwygqq+LPq/bQtUsgN18wwOqQlAfTJOFnXlq9h0PHyvn3LecQFKCtjf5IRHjy+6MoqXDw+MfpRIQEcM0EvdNeNU2/JfxI5tES/vr5PmaPSWDywCirw1EWstuE5+aN4cIhMdz/3nbe3nDQ6pCUh9Ik4Seqa5zc+5+tRIQE8NAVjedVVP4oOMDOK9encOHgGO57dztvrjtgdUjKA2mS8BN/SdvL9sMnePyqEUSFB1sdjvIQIYF2Xr4+hYuHxfLr93ew4Jv9VoekPIwmCT+w88gJ/rwyi1mjE7h0ZLzV4SgPExJo568/SmFGUhwPL97JHz/L1OGxqp4mCR9X6ajh3v9so3toEI9dmWx1OMpDBQXYeOmH47hmfG/+vGoPD7y3XacYV4CObvJ5T3ycTnpOMX+7YTw9woKsDkd5sAC7jaeuHkls12D+b9Ue8ksqef7asYQH69eEP9MrCR+2ZNsR/rXmADef359pOoGfagUR4Z4ZQ/ntVSNI253P91/6moOFZS3vqHyWJgkf9V3BSe5/dztj+3TnvkuHWR2O8jLXn9OXBT+eSG5xJbNf/Ir0whqrQ1IW0SThg8qrarj9zc0E2IUXrhtHoF1Ps2q78wdH88Ed5xEZFsQfNlbwl7S9OJ3aoe1v9NvDxzidhrvf2ULG0WL++IPR9NK5mdRZ6B8dxgd3nEdKnJ3ff5rBTQs2cPxkldVhKTfy6x6pu+66i+7dzzxN9hVXXMG9994LQGpqKvPnz2f+/PkUFBQwZ86cFo/RuPw999zDrFmzyMzM5Lbbbmtx/8blf/e733HuuefyzTff8OCDD9aXKyoqonv37hw8VsaRonL6RoXx2JoQQhqVf/nllxk6dCgfffQRzz77bIvHb1x+0aJFREdH89prr/Haa6+1uH/j8mlpaQA888wzLFmypMl96uoCnFJ+zZo1vPvuuwA88MADrFmz5ozHjoqKOqV8YWEhr7zyCgC33noru3fvPuP+Q4YMOaV8VFQUTz75JABXX301hYWFZ9y/V69epKam1pefPHnyKT9LLfGUn70jB/aR9dr93HD7/7JwTyEX3PNXAr99m25dAs+4f+OfVU/52Xv77bd55JFHWtzfG372WlOPs6VXEj4kr6SSI0XlxHUNIb5biNXhKB8zc0Q8791+LiGBNtJzitlfcBKn3k/h88Tbb5oZP3682bhxY7v2TUtLa9VfdN7guXdW8MKWKiYPjOKf8ycQ4MX9EL5yXnylHnBqXSqqa3j600z+8fV39IsK5XffH8m5A6OtDbANfPW8tJWIbDLGjG+pnPd+k6h6n+/O54VvK0lK6MqLPxzn1QlCeb6QQDsPzUrirVsm4TRw3avruPvtLRSWVlodmuoE+m3i5b7ZU8Ct/9pIQriNf/1kIl1DztxOrFRHOXdgNJ/96kLunDKIj7Yd4eJnP+efX39Htd6p7VM0SXixL7PyuWnBRvpFhfE/E0LoHqp3VCv3Cgm0c+8lQ1n6iwsY0asrj360i0v+9AUrduXq/E8+QpOEl3pvczY//ucG+kaF8sbNk4gI0kXtlXUGx0Xwxk2T+PuN40Hg5n9tZM5f1/D1ngJNFl7Or4fAeiNjDH/5fC9Pf5rJuQOj+Ov1KdrEpDyCiDB1eBwXDonh7Q2HeGHVHn74t3VM6h/JnRcP4vxB0YjoHzPeRpOEFymrcvCb93fw3reHuXJ0As/MHa1LkCqPE2i38aNz+jInJZGF6w/yUtperv/7epLiu3LbRQO4bGS8zgLgRfRMeYk9eaVc9eLXvL/lML+aNoTnrhmjCUJ5tJBAO/PP68+X903h91ePpMJRwy8XbuH836/iuRW7yS2usDpE1Qp6JeHhjDG8s/EQj320i+BAO//6yUQuGBxjdVhKtVpwgJ1rJvRhbkpvVmXk8fraAzy3IosXVu1hyrBY5qQkMmVorP7R46E0SXiwg4Vl3P/eNr7ZW8ik/pE8N28M8d10LiblnWw2YVpSHNOS4thfcJJ/rz/Ie98eZvmuXCLDgrh8ZDyXj4pnQr9I7Dbtu/AUmiQ8UHlVDX//ah8vrt6L3SY88b0RXDuhDzb9xVE+ol90GA9cNpz/uWQoX2YVsGhTNv/ZdIjX1x4gJiKYGUlxTBsex+SBUYQE2q0O169pkvAgjhon727O5o/Ld5NbXMmMpDgenZ2sVw/KZwXYbUwZFsuUYbGcrHSwOjOPj7fl8P63h3lz3UG6BNo5d2AU5w+O5oLB0QyMCdcRUm6mScIDlFU5+M/GbP7+1XccPFbG2D7deeG6cUzoF2l1aEq5TVhwAFeMSuCKUQlUVNewdl8hK9Pz+CIrn5UZeQDERgQzsX8kE/tHMr5vJEPiwnUamk6mScJCGUeLeXdTNu9szOZEeTXj+nTn15cPZ0ZSnP61pPxaSKCd1KGxpA6NBeDQsTK+zCpg7b5CNuw/xpJtOQB0CbQzMrEboxO7kZzQjeSErgyICdc+jQ6kScKNjDHszT/JyvRcFm89ws4jxQTYhOlJcdx8QX9S+uqVg1JN6R0ZynWT+nDdpD4YY8g+Xs7mg8f59mARWw4VseCbA1S55owKDrAxKDacIXERDIoNZ2BMGANiwukTGar9G+3g1iQhIjOB5wE78DdjzFON3g8G/gWkAIXANcaY/e6MsaPlFVewfv8xNnx3jM9357Pftaj8yF7deHhWEleOTiAqPNjiKJXyHiJC78hQekeGMntMLwCqa5zszS9l5+Fi0nOK2Z1Xypq9hbz/7eEG+0F81xAibFV8nL+VXj26kNC9CwndutCzWzBxXUMIDw7Qq/hG3JYkRMQOvAhMB7KBDSKy2Bizq0Gxm4DjxphBIjIP+D1wjbtiPBsV1TUcKCzju4KTZB4tYVfOCXYeKSb7eDlQe1k8sX8kN53fn4uHx+myokp1oEC7jWE9uzKsZ9dTXi+pqGZ/QRn7CkrZl3+SQ8fL2L4vh89355NfWknjaaVCg+zERAQTEx5MdHgwUeFBRIUF0SMsiMiwILp1Cax/dO0SSNeQQJ+/v8OdVxITgT3GmH0AIrIQmA00TBKzgUdczxcBL4iIGDfPEOZ0GiocNZRV1VBWWUNxRTUlFQ6KK6o5frKKwpNVFJZWkVtcQc6JcnJOVHC0uKL+B04E+keFMbp3d26c3I8J/SNJTuiqUxEo5WYRIYGMTOzGyMRu9a+lpRWRmppKpaOG3BOVHC4qJ7e4gtzi2t/jgtIq8ksq2JNfyvr9VRwvqzotmTQUHGAjIiSA8OAAwuoeQXZCgwLoEmSnS6CdLkF2QgLthATaCAmofR4cYCMowFb/b1CAjSB77b8BNhtBAUKAzUaAXQi027DbhECbDbtdCLDVPtzx1ejOJNELONRgOxuY1FwZY4xDRE4AUUBBRwfzzoZD/OnLMgLXr8ZR46SqxlDpqKHK4aTS0fJ8+KFBdnp2DaFntxAmD4yib2QY/aJD6R8dxsCYcMKCtbtHKU8WHGCnT1QofaJCz1iuxmkoKquiqLyaorJqisqqKK6oprjcQXF5NaWVDkoqHZRUOCirdHCyykHhySoOHiujvKqG8uraPzhb873SVjckBTGlwz/1VO78Jmuqoa9xGmxNGUTkVuBWgLi4uPoFy9vicJ6Dnl2chARWYrcJAQKBdgi02Qmw2Qm2Q7BdCLZDaKDQJUAIDYCIICEiSAiy14Va4XoUQREcK4Jje9oczlkrLS1t1/+DJ/KVuvhKPUDr0pAd6OF6YAdCXY9mS9d2ljuNweGEqhqodhqqnVDteu5wQrUTakzD57UJqsaA04DD9ZrT1JZzGugZVNnp58WdSSIb6N1gOxE40kyZbBEJALoBxxp/kDHmFeAVqF3juj1rvKYCY3WtW4/kK3XxlXqA1sVTuaMu7mwk3wAMFpH+IhIEzAMWNyqzGLjR9XwOsMrd/RFKKaX+y21XEq4+hjuBZdRef/3DGLNTRB4DNhpjFgN/B14XkT3UXkHMc1d8SimlTufW3lVjzFJgaaPXHmrwvAKY686YlFJKNU/HZCqllGqWJgmllFLN0iShlFKqWZoklFJKNUuThFJKqWaJt9+GICL5wIF27h5NJ0z5YRGti+fxlXqA1sVTnU1d+hpjYloq5PVJ4myIyEZjzHir4+gIWhfP4yv1AK2Lp3JHXbS5SSmlVLM0SSillGqWvyeJV6wOoANpXTyPr9QDtC6eqtPr4td9Ekoppc7M368klFJKnYHfJgkR2S8i20Vki4hstDqethCRf4hInojsaPBapIgsF5Es1789rIyxNZqpxyMicth1XraIyGVWxthaItJbRFaLSLqI7BSRX7pe98bz0lxdvOrciEiIiKwXka2uejzqer2/iKxznZO3XUsXeLQz1OU1EfmuwTkZ0+HH9tfmJhHZD4w3xnjdeGkRuRAoBf5ljBnheu1p4Jgx5ikRuR/oYYy5z8o4W9JMPR4BSo0xz1gZW1uJSDwQb4zZLCIRwCbgKmA+3ndemqvLD/CicyMiAoQZY0pFJBD4CvglcDfwnjFmoYj8FdhqjPmLlbG25Ax1+SmwxBizqLOO7bdXEt7MGPMFp6/YNxtY4Hq+gNpfao/WTD28kjEmxxiz2fW8BEinds12bzwvzdXFq5hapa7NQNfDABcDdV+q3nJOmqtLp/PnJGGAz0Rkk2vNbG8XZ4zJgdpfciDW4njOxp0iss3VHOXxzTONiUg/YCywDi8/L43qAl52bkTELiJbgDxgObAXKDLGOFxFsvGSBNi4LsaYunPyhOuc/ElEgjv6uP6cJM4zxowDLgXucDV9KOv9BRgIjAFygGetDadtRCQceBe4yxhTbHU8Z6OJunjduTHG1BhjxgCJwERgeFPF3BtV+zSui4iMAB4AhgETgEigw5sy/TZJGGOOuP7NA96n9gfIm+W62pLr2pTzLI6nXYwxua5fBifwKl50Xlxtxe8Cbxpj3nO97JXnpam6ePO5McYUAWnAOUB3EalblTMROGJVXO3RoC4zXU2DxhhTCfyTTjgnfpkkRCTM1SGHiIQBM4AdZ97L4y0GbnQ9vxH40MJY2q3uC9Xle3jJeXF1LP4dSDfG/LHBW153Xpqri7edGxGJEZHuruddgGnU9q+sBua4innLOWmqLhkN/gARavtWOvyc+OXoJhEZQO3VA9Su8/2WMeYJC0NqExH5N5BK7QyQucDDwAfAO0Af4CAw1xjj0Z3CzdQjldrmDAPsB26ra9P3ZCJyPvAlsB1wul5+kNq2fG87L83V5Vq86NyIyChqO6bt1P5B/I4x5jHX7/9CaptnvgV+5PpL3GOdoS6rgBhAgC3ATxt0cHfMsf0xSSillGodv2xuUkop1TqaJJRSSjVLk4RSSqlmaZJQSinVLE0SSimlmqVJQimlVLM0SSillGqWJgmlOpiIpInIUNfzqIbrZSjlbTRJKNXxBgFZruejqL1zWSmvpElCqQ4kIn2Bw65J8KA2SWyzMCSlzoomCaU61hhOTQopaJJQXkyThFIdazQQAiAig6ldmU6bm5TX0iShVMcaA9hEZCvwELVTU9945l2U8lw6C6xSHUhE9gBjXWtDK+X19EpCqQ7iWsjKqQlC+RK9klBKKdUsvZJQSinVLE0SSimlmqVJQimlVLM0SSillGqWJgmllFLN0iShlFKqWZoklFJKNUuThFJKqWb9f1DaMij+YGBvAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mus   = np.linspace(max(0., nug - 15), nug + 15, 100)\n",
    "tmus  = [ht.tmu(x0, mui, parbest = parbest) for mui  in mus]\n",
    "ptmus = [htcore.tmu_pvalue(tmui)            for tmui in tmus]\n",
    "htplot.plt_wilks_pvalues(beta, mus, ptmus);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: Consider now a uniform backgrund $1.3 \\times 10^{-4}$ c/(kg y keV) and and two gaussian peaks with sigma 1%, one signal peak, $^{136}$Xe at 4458 keV with unknown number of events and other background peak ($^{214}$Bi) at 4447 keV with 0.06 c/(keV y). Consider the range 2400-2500 keV. \n",
    "\n",
    "Compute the upper limit in number of signal events at 90% CL vs exposure in the case there is no signal. \n",
    "\n",
    "Compute the observation potential in the case the number of signal is 0.12 c/(keV y)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}