Codes/08_LindleysParadox.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lindley's Paradox"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To ilustrate the paradox we are going to consider the Laplace's problem of estimating the probability of a male birth.\n",
    "\n",
    "Suppose that a total of 48,870 girls and 49,581 boys were born. Let be $\\theta$ the **probability of a male birth** and $Y|\\theta$ the observed male births. As Laplace, we are going to consider a uniform prior for $\\theta$. The model is then $$\\theta\\sim \\textsf{Uniform}(0,1)$$ and $$Y|\\theta\\sim\\textsf{Binomial}(n,\\theta),$$ where $n$ is the total number of births."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import norm, binom\n",
    "\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "fem_births = 48870\n",
    "mal_births = 49581\n",
    "births = fem_births + mal_births"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We get a simple estimate for $\\theta$ as the observed proportion of male births."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5036109333577109"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "theta_hat = mal_births/births\n",
    "theta_hat "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To ilustrate Lindley's Paradox we are going to perform the next hypothesis testing $H_0:\\theta=0.5$."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we are going to estimate the frecuentist p-value. Since we have large data, we can approximate the distribution of $Y|\\theta$ by a normal distribution $$Y\\overset{\\cdot}{\\sim}\\textsf{N}(\\mu,\\sigma^2),$$ with $\\mu=n\\theta$ and $\\sigma^2=n\\theta(1-\\theta)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta0 = 0.5\n",
    "mu = births * theta0\n",
    "sigma_2 = births * theta0 * (1-theta0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.023451501665971836"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# p-value for a double side hypothesis test\n",
    "2 * norm.sf(mal_births, mu, np.sqrt(sigma_2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.023256586093285063"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# p-value for a double side hypothesis test considering the exact distribution of Y\n",
    "2 * binom.sf(mal_births, births, theta0)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That is, from a frecuentist perspective **we reject the hypothesis $H_0: \\theta=0.5$.**\n",
    "\n",
    "To get a better insight of why, we are rejecting $H_0$, let's calculate the best estimate for $\\theta$ under $H_0$, and plot the density distribution for $Y$ under $H_0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "49225.5"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "births*theta0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f0a712f5fd0>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj0AAAGdCAYAAAD5ZcJyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABhg0lEQVR4nO3de1xUdf4/8NfMwMxwH5A7gqKgeL+gEqZpRWLZhWrNdd0s86vuZmVrbWVb2u7Wz9Jqu1l2+abWt9Jst5u5FmqpJaLi/YagKHgZLgIz3AdmPr8/hjk6isogcObyej4e84A553PmvD+gM28+V4UQQoCIiIjIzSnlDoCIiIioMzDpISIiIo/ApIeIiIg8ApMeIiIi8ghMeoiIiMgjMOkhIiIij8Ckh4iIiDwCkx4iIiLyCF5yB+BMLBYLzpw5g4CAACgUCrnDISIiolYQQqCqqgrR0dFQKi/fnsOk5wJnzpxBbGys3GEQERFRGxQVFaFr166XPc+k5wIBAQEArD+0wMBAmaMhInJDNTVAdLT1+zNnAD8/eeMht2A0GhEbGyt9jl8Ok54L2Lq0AgMDmfQQEXUEler894GBTHqoXV1taAoHMhMREZFHYNJDREREHoFJDxEREXkEjukhIiK3JYRAU1MTzGaz3KHQNVCpVPDy8rrm5WSY9BARkVsymUw4e/Ysamtr5Q6F2oGvry+ioqKgVqvb/BpMeoiIyO1YLBYUFBRApVIhOjoaarWai866KCEETCYTSktLUVBQgMTExCsuQHglTHqIiMjtmEwmWCwWxMbGwtfXV+5w6Br5+PjA29sbJ0+ehMlkglarbdPrcCAzERG5rba2CJDzaY/fJf81EBERkUdg0kNEREQegUkPEREReQQmPURERE7kwQcfREZGxiXHf/nlFygUClRWVgIA9u3bh9GjR0Or1SI2NhaLFi3q3EBdEGdvEZHLE0LgaHE1fssvg9kicOfgaEQEtm12B5ErMBqNGDduHNLS0rB06VLs378fDz30EHQ6HWbOnCl3eE6LSQ8RuaRGswXf7TmDzXml+C3/HMqqG6Rzr6w7gvH9I/HAyO4Y1i2Y67MQAGtyXNcoz8rMPt6qdv13+Nlnn8FkMuHjjz+GWq1Gv379sGfPHrz++utMeq6ASQ8RuZw6kxl/+r8cbDpaKh3TeisxvHsI6hvN2HGiAmv2ncWafWfRNyoQz93eByN7hsoYMTmDukYz+s7/UZZ7H/pHOnzV7feRm5WVhRtuuMFudeL09HS88sorqKioQHBwcLvdy520aUzPkiVL0L17d2i1WqSkpGD79u1XLL969WokJSVBq9ViwIABWLt2rd15IQTmz5+PqKgo+Pj4IC0tDXl5edL5EydOYPr06YiPj4ePjw969uyJBQsWwGQy2ZVRKBSXPLZt29aWKhKRkzLWN+KBj7dj09FS+HirMPvGnvhixnXYu2AcPp2egtV/GokfHhuFScNiofFS4tBZI6Yt24GdJ8rlDp2o1dasWQN/f3+7x6233iqd1+v1iIiIsLvG9lyv13dqrK7E4bRz1apVmDt3LpYuXYqUlBS88cYbSE9PR25uLsLDwy8pv3XrVkyePBkLFy7E7bffjs8//xwZGRnYtWsX+vfvDwBYtGgR3nrrLaxYsQLx8fF4/vnnkZ6ejkOHDkGr1eLIkSOwWCx4//33kZCQgAMHDmDGjBmoqanBq6++ane/9evXo1+/ftLzLl26OFpFInJS5TUmTP04GwdOGxGg8cKyacMxrHvIJeX6RQfhld8NxLzbkvCXVXvwc24pHlq+A6v/NBK9IwNkiJycgY+3Cof+kS7bvR1x44034r333rM7lp2djT/+8Y/tGZbnEQ4aMWKEmD17tvTcbDaL6OhosXDhwhbL33fffWLChAl2x1JSUsSsWbOEEEJYLBYRGRkpFi9eLJ2vrKwUGo1GfPHFF5eNY9GiRSI+Pl56XlBQIACI3bt3O1olicFgEACEwWBo82sQUcc4W1kn0l77RXR7eo0Y8o+fxP5Tla26rrahSdy95FfR7ek1YsRLmaKovKaDI6Urqq4WArA+qqs77DZ1dXXi0KFDoq6ursPu0VEeeOABcdddd11y/OeffxYAREVFhbj//vsvKbNx40YBQJSXl3dOoJ3sSr/T1n5+O9S9ZTKZkJOTg7S0NOmYUqlEWloasrKyWrwmKyvLrjxg7Xe0lS8oKIBer7crExQUhJSUlMu+JgAYDAaEhFz6F96dd96J8PBwjBo1Ct99990V69PQ0ACj0Wj3ICLnU2tqwh8+3Ia8kmpEBmrx5axU9I8JatW1PmoVPn5wOBLD/VFsbMDUj7ejvMZ09QuJnFhqaio2b96MxsZG6VhmZiZ69+7N8TxX4FDSU1ZWBrPZ3GI/4uX6EC/X72grb/vqyGvm5+fj7bffxqxZs6Rj/v7+eO2117B69Wr88MMPGDVqFDIyMq6Y+CxcuBBBQUHSIzY29rJliUg+b27Iw/GyGkQGarH6T6lICPd36HqdrxqfTB+B6CAtjpfWYNqy7agzyTOLh6g9/OEPf4Barcb06dNx8OBBrFq1Cm+++Sbmzp0rd2hOzeUWJzx9+jTGjx+PiRMnYsaMGdLx0NBQzJ07FykpKRg+fDhefvll/PGPf8TixYsv+1rz5s2DwWCQHkVFRZ1RBSJywOGzRny0pQAA8GJGf8SGtG3H7KggH3wyPQXBvt7Ye8qA9zYda88wiTpVUFAQfvrpJxQUFCA5ORlPPPEE5s+fz+nqV+HQQObQ0FCoVCoUFxfbHS8uLkZkZGSL10RGRl6xvO1rcXExoqKi7MoMHjzY7rozZ87gxhtvxMiRI/HBBx9cNd6UlBRkZmZe9rxGo4FGo7nq6xCRPMwWgXn/2Q+zReDW/pFI6xtx9YuuICHcHy/dPQAPf7YL7286hknDYxGj82mnaInax/Lly1s8PnbsWAghpOcDBw7Eli1bOikq9+BQS49arUZycjI2bNggHbNYLNiwYQNSU1NbvCY1NdWuPGDtd7SVj4+PR2RkpF0Zo9GI7Oxsu9c8ffo0xo4di+TkZCxbtqxVW8zv2bPHLpEiItfyefZJ7CmqhL/GCy/c2e/qF7TCrf0jMSI+BA1NFrz83yPt8ppE5BocnrI+d+5cPPDAAxg2bBhGjBiBN954AzU1NZg2bRoAYOrUqYiJicHChQsBAHPmzMGYMWPw2muvYcKECVi5ciV27twptdQoFAo8/vjjePHFF5GYmChNWY+Ojpb2HrElPN26dcOrr76K0tLzC5LZWopWrFgBtVqNIUOGAAD+85//4OOPP8ZHH33U9p8OEcmm2FiPRetyAQBPje/dbttKKBQKzL+9L+5451d8v/cMHkjt1uK0dyJyPw4nPZMmTUJpaSnmz58PvV6PwYMHY926ddJA5MLCQrtWmJEjR+Lzzz/Hc889h2effRaJiYn45ptvpDV6AOCpp55CTU0NZs6cicrKSowaNQrr1q2DVmt9k8vMzER+fj7y8/PRtWtXu3gubOr75z//iZMnT8LLywtJSUlYtWoVfve73zlaRSJyAn///iCqGpowOFaHKSnd2vW1+8cEYdKwWKzcUYS/f38I386+Hkolt6ogcncKcWHW4OGMRiOCgoJgMBgQGBgodzhEHmvjkWI8tHwnVEoF1jw6Cn2i2v//Y2lVA2589RdUNzRh8e8GYuIwzt7sFDU1gH/z7LvqasDPr0NuU19fj4KCAsTHx0t/QJNru9LvtLWf3y43e4uI3JsQAv/KtG5DM31UfIckPAAQFqDBozclAAAW/ZiLmoamDrkPETkPJj1E5FSyC8qx/7QBGi8lZt3Qo0Pv9eD13dGtiy9Kqxrw7i/5HXovIpIfkx4iciofbj4OAPhdcld08e/YJSU0Xio8e1sfAMDy306gqr7xKlcQkStj0kNETiO/pAobjpRAobB2bXWGcX0jkBDujxqTGV/lnOqUexKRPJj0EJHTsK28fEufCPQIc2yribZSKBR4YGR3AMCKrSdgsXBuB7mG7t2744033pA7jHbTGfVh0kNETqGkqh7/2XUaADBrTMeO5bnYPUNiEKD1wolztdiUV3r1C4g6WFFRER566CFER0dDrVajW7dumDNnDs6dOyd3aC6NSQ8ROYVPs07CZLZgaJwOyd06d7FAP40X7muesr78txOdem+iix0/fhzDhg1DXl4evvjiC+Tn52Pp0qXS7gfl5eWyxGU2m2GxWGS5d3th0kNEsqs1NeHTbScBADM7eMbW5UxN7QaFAth0tBTHS6tliYEIAGbPng21Wo2ffvoJY8aMQVxcHG699VasX78ep0+fxt/+9jepbFVVFSZPngw/Pz/ExMRgyZIl0jkhBF544QXExcVBo9EgOjoajz32mHS+oaEBTz75JGJiYuDn54eUlBT88ssv0vnly5dDp9Phu+++Q9++faHRaPDRRx9Bq9WisrLSLuY5c+bgpptukp7/+uuvGD16NHx8fBAbG4vHHnsMNTU10vmSkhLccccd8PHxQXx8PD777LN2/AleHpMeIpLdVzmnUFnbiG5dfHFL35Y3L+5o3br44abe4QCAT7JOyhIDdTAhrIsjyvFo5TrA5eXl+PHHH/Hwww/Dx8d+M9zIyEhMmTIFq1atknYjWLx4MQYNGoTdu3fjmWeewZw5c6SNtv/973/jX//6F95//33k5eXhm2++wYABA6TXe+SRR5CVlYWVK1di3759mDhxIsaPH4+8vDypTG1tLV555RV89NFHOHjwIKZMmQKdTod///vfUhmz2YxVq1ZhypQpAIBjx45h/PjxuPfee7Fv3z6sWrUKv/76Kx555BHpmgcffBBFRUX4+eef8dVXX+Hdd99FSUmJg7/QNhAkMRgMAoAwGAxyh0LkMZrMFjH6lY2i29NrxCdbC2SNZfPREtHt6TWi3/x1wlhnkjUWt1VdLYQ1BbB+30Hq6urEoUOHRF1dXcv37uxHK+u6bds2AUB8/fXXLZ5//fXXBQBRXFwsunXrJsaPH293ftKkSeLWW28VQgjx2muviV69egmT6dJ/yydPnhQqlUqcPn3a7vjNN98s5s2bJ4QQYtmyZQKA2LNnj12ZOXPmiJtuukl6/uOPPwqNRiMqKiqEEEJMnz5dzJw50+6aLVu2CKVSKerq6kRubq4AILZv3y6dP3z4sAAg/vWvf132Z9Pi77RZaz+/2dJDRLL6Lb8MheW1CPLxxu+S5d0KYlRCKHqG+aG6oQn/5vR1kpFoZctQamrqJc8PHz4MAJg4cSLq6urQo0cPzJgxA19//TWamqwrj+/fvx9msxm9evWCv7+/9Ni0aROOHTsmvZ5arcbAgQPt7jFlyhT88ssvOHPmDADgs88+w4QJE6DT6QAAe/fuxfLly+1eNz09HRaLBQUFBTh8+DC8vLyQnJwsvWZSUpJ0fUdyeMNRIqL29PVu64ytOwZFwUetkjUW2/T1+d8exCdZJzE1tTs3InUnvr7W/b7kuncrJCQkQKFQ4PDhw7j77rsvOX/48GEEBwcjLCzsqq8VGxuL3NxcrF+/HpmZmXj44YexePFibNq0CdXV1VCpVMjJyYFKZf//zt///HIRPj4+UCjs/w8MHz4cPXv2xMqVK/HnP/8ZX3/9NZYvXy6dr66uxqxZs+zGD9nExcXh6NGjV429ozDpISLZ1DQ0Yd0BPQDg7iFdZY7G6p6hXbF4XS6Ol9VgS34ZxvS6+ocLuQiFosM2OG0vXbp0wS233IJ3330Xf/nLX+zG9ej1enz22WeYOnWqlIhs27bN7vpt27ahT58+0nMfHx/ccccduOOOOzB79mwkJSVh//79GDJkCMxmM0pKSjB69GiH45wyZQo+++wzdO3aFUqlEhMmTJDODR06FIcOHUJCQkKL1yYlJaGpqQk5OTkYPnw4ACA3N/eSwdEdgd1bRCSbHw/qUddoRvcuvhgap5M7HACAv8YL9wyNAQB8vYtdXNT53nnnHTQ0NCA9PR2bN29GUVER1q1bh1tuuQUxMTF46aWXpLK//fYbFi1ahKNHj2LJkiVYvXo15syZA8A6++p///d/ceDAARw/fhz/93//Bx8fH3Tr1g29evXClClTMHXqVPznP/9BQUEBtm/fjoULF+KHH364aoxTpkzBrl278NJLL+F3v/sdNJrzW8Y8/fTT2Lp1Kx555BHs2bMHeXl5+Pbbb6WBzL1798b48eMxa9YsZGdnIycnB//zP/9zycDtjsCkh4hkY+vayhgSc0kTupzuHGxNen46VIw6k1nmaMjTJCYmYufOnejRowfuu+8+9OzZEzNnzsSNN96IrKwshIScX8fqiSeewM6dOzFkyBC8+OKLeP3115Geng4A0Ol0+PDDD3H99ddj4MCBWL9+Pb7//nt06dIFALBs2TJMnToVTzzxBHr37o2MjAzs2LEDcXFxV40xISEBI0aMwL59+6RZWzYDBw7Epk2bcPToUYwePRpDhgzB/PnzER0dLZVZtmwZoqOjMWbMGNxzzz2YOXMmwsPD2+PHd0UK0drRUh7AaDQiKCgIBoMBgYGBcodD5NZKjPW4buEGWASw6a9j0a2L83Q7CCEwetHPOFVRh3f+MAS3D4y++kXUOjU1gG3MSHV1h3U31dfXo6CgAPHx8dBqtR1yD+pcV/qdtvbzmy09RCSLb/ecgUUAyd2CnSrhAawDmu8YZE10vttzRuZoiKi9MOkhIln8p7lr6+4hMTJH0rI7m5OeX3JLYahrlDkaImoPTHqIqNMd0Rtx+KwRapUStw+MkjucFiVFBiAx3B8mswU/HtTLHQ4RtQMmPUTU6b5u3k39xqQw6HzVMkfTMoVCIbX2fL+XXVxE7oBJDxF1KrNF4Js9tq4t51ib53Js43p+yy9DaVWDzNEQ0bVi0kNEnSrr2DkUGxsQ5OONG5Oce+G/7qF+GNQ1CBYBrN1/Vu5wqA04Qdl9tMfvkkkPEXUqWyvPhIFR0HjJu+1Ea0izuNjF5VK8vb0BWHcJJ/dg+13afrdtwW0oiKjTNJktWH+4GABwh4usfXPHoGi8tPYwck5WoKi8FrEhrdtDieSlUqmg0+lQUlICAPD19XWqBTCp9YQQqK2tRUlJCXQ63SV7hTmCSQ8RdZqdJytQWdsIna83hncPljucVokI1OK6+C7IOn4O3+87g4fHtryfEDmfyMhIAJASH3JtOp1O+p22FZMeIuo0mYesrTw3JYXDS+U6vet3Do5G1vFz+G4Pkx5XolAoEBUVhfDwcDQ2cq0lV+bt7X1NLTw2THqIqFMIIfDTIet6N+P6RsgcjWNu7R+J5745gCP6KnZxuSCVStUuH5jk+lznTy0icmm5xVUoKq+D2kuJ0YnOPWvrYjpfNZK7WbvjNjSPSSIi18Okh4g6ReZBa7IwOiEUfhrXa2RO62PdAXrDEY4PIXJVTHqIqFNkNreQ3OJiXVs2N/exxr3t+DlU1XN8CJErYtJDRB3urKEO+04ZoFCcTx5cTc8wf8SH+qHRLLAlr0zucIioDZj0EFGHW988a2tIrA5hARqZo2m7m5Oau7gOs4uLyBUx6SGiDvfTIVvX1rWtsSG3m5rH9fycWwKzhdsbELkaJj1E1KGM9Y3YdvwcAGBcP9fs2rIZ3j0EAVovlNeYsKeoQu5wiMhBTHqIqENtyi1Fo1mgR5gfeob5yx3ONfFWKTG2t7W1Zz27uIhcDpMeIupQmYdce9bWxWxT1zcy6SFyOUx6iKjDmJos+Ll5XRtXW4X5csb0CoNKqWhebJE7eBO5EiY9RNRhdp4oR1VDE0L91Rgc6xobjF4NV2cmcl1Meoiow2zKKwUA3NDcOuIuuDozkWti0kNEHWbzUesifmN6udZeW1fD1ZmJXBOTHiLqECVV9Th81ggAuD4hVOZo2hdXZyZyTUx6iKhD/NqcDPSPCUSov+uuwnw5NzWvzrwpt1TmSIiotZj0EFGHsLWAjE50r64tm9GJ1tarX/PLIARXZyZyBUx6iKjdWSwCW2yDmN006UmJ7wK1SonTlXUoKKuROxwiagUmPUTU7g7rjSirNsFXrZKmd7sbnwvq9ms+x/UQuQImPUTU7myztlJ7dIHay33fZkY1d3FxMDORa3DfdyMiks3mo9auLdu4F3dlq9+2Y+fQZLbIHA0RXQ2THiJqV7WmJuw8WQ7AuiihO+sXHQSdrzeqGpqw91Sl3OEQ0VUw6SGidrXt+Dk0mgW6BvsgPtRP7nA6lEqpwPU9ra09ti49InJeTHqIqF3ZPvxHJ4ZBoXCfrScuZ9QFU9eJyLkx6SGidrW5ear6mF7uPZ7HZlTzatN7iiph5JYURE6NSQ8RtZtTFbU4XloDlVKB1J6ekfTEhviiexdfmC0C246dkzscIroCJj1E1G5sU7cHx+oQ5OMtczSdx7bqNLu4iJwbkx4iaje2VZjdfar6xaRxPVyvh8ipMekhonZhsQhsbe7e8bSkJ7VnF6iUChwvq8Hpyjq5wyGiy2hT0rNkyRJ0794dWq0WKSkp2L59+xXLr169GklJSdBqtRgwYADWrl1rd14Igfnz5yMqKgo+Pj5IS0tDXl6edP7EiROYPn064uPj4ePjg549e2LBggUwmUx2r7Nv3z6MHj0aWq0WsbGxWLRoUVuqR0RtcERfhcraRviqVRjYVSd3OJ0qUOuNQV2DAAC/5nHXdSJn5XDSs2rVKsydOxcLFizArl27MGjQIKSnp6OkpKTF8lu3bsXkyZMxffp07N69GxkZGcjIyMCBAwekMosWLcJbb72FpUuXIjs7G35+fkhPT0d9fT0A4MiRI7BYLHj//fdx8OBB/Otf/8LSpUvx7LPPSq9hNBoxbtw4dOvWDTk5OVi8eDFeeOEFfPDBB45WkYjaYOsxa9fO8O4h8FZ5XiPyqOZxPdySgsiJCQeNGDFCzJ49W3puNptFdHS0WLhwYYvl77vvPjFhwgS7YykpKWLWrFlCCCEsFouIjIwUixcvls5XVlYKjUYjvvjii8vGsWjRIhEfHy89f/fdd0VwcLBoaGiQjj399NOid+/era6bwWAQAITBYGj1NURkNX35dtHt6TXivV/y5Q5FFtsLzoluT68Rg//+ozCbLXKH47yqq4UArI/qarmjITfR2s9vh/4cM5lMyMnJQVpamnRMqVQiLS0NWVlZLV6TlZVlVx4A0tPTpfIFBQXQ6/V2ZYKCgpCSknLZ1wQAg8GAkJAQu/vccMMNUKvVdvfJzc1FRUVFi6/R0NAAo9Fo9yAix5ktAtkF1q0nUnt0kTkaeQyO1cFPrUJFbSOO6KvkDoeIWuBQ0lNWVgaz2YyIiAi74xEREdDr9S1eo9frr1je9tWR18zPz8fbb7+NWbNmXfU+F97jYgsXLkRQUJD0iI2NbbEcEV3ZwTMGVNU3IUDjhX7RgXKHIwtvlRLDulv/ENt2nOv1EDkjl+t4P336NMaPH4+JEydixowZ1/Ra8+bNg8FgkB5FRUXtFCWRZ8lqnrU1Ij4EXh44nsfmuuZWLiY9RM7JoXen0NBQqFQqFBcX2x0vLi5GZGRki9dERkZesbzta2te88yZM7jxxhsxcuTISwYoX+4+F97jYhqNBoGBgXYPInJcVvOHfGpPz+zasrmuh7WlJ7ugHBaLkDkaIrqYQ0mPWq1GcnIyNmzYIB2zWCzYsGEDUlNTW7wmNTXVrjwAZGZmSuXj4+MRGRlpV8ZoNCI7O9vuNU+fPo2xY8ciOTkZy5Ytg1JpH3pqaio2b96Mxsbze99kZmaid+/eCA4OdqSaROSARrMFO5rH81znoeN5bPrHBMFPrYKhjuN6iJyRw+3Qc+fOxYcffogVK1bg8OHD+POf/4yamhpMmzYNADB16lTMmzdPKj9nzhysW7cOr732Go4cOYIXXngBO3fuxCOPPAIAUCgUePzxx/Hiiy/iu+++w/79+zF16lRER0cjIyMDwPmEJy4uDq+++ipKS0uh1+vtxur84Q9/gFqtxvTp03Hw4EGsWrUKb775JubOnXstPx8iuor9pw2oMZkR5OONvlGe3VrKcT1Ezs3L0QsmTZqE0tJSzJ8/H3q9HoMHD8a6deukQcOFhYV2rTAjR47E559/jueeew7PPvssEhMT8c0336B///5Smaeeego1NTWYOXMmKisrMWrUKKxbtw5arRaAtcUmPz8f+fn56Nq1q108QlibkIOCgvDTTz9h9uzZSE5ORmhoKObPn4+ZM2c6/lMholazjedJiQ+BUqmQORr5XdejCzYdLcW24+fw0Kh4ucMhogsohC1rIBiNRgQFBcFgMHB8D1Er3f+/2diSV4YFd/TFtOv5Ib+7sAJ3v7sVQT7e2P38LUwEL1ZTA/j7W7+vrgb8/OSNh9xCaz+/PXeaBRFdM1OTBTtPWNfB8vRBzDYDOK6HyGkx6SGiNtt7qhJ1jWaE+KnRKzxA7nCcgpdKieHxHNdD5IyY9BBRm9nG81zXg+N5LsT1eoicE5MeImozW9LjqVtPXI4t6eF6PUTOhUkPEbVJfaMZOYUcz9OS/tGB0riew3ru6UfkLJj0EFGb7C6shKnJgrAADXqG+csdjlOxH9dTLnM0RGTDpIeI2mR78yrMKfEhUCg4nudiHNdD5HyY9BBRm2w/0bwoIcfztMiW9GznuB4ip8Gkh4gc1mi2YNfJSgDAiOZtF8gex/UQOR8mPUTksAOnDahrNEPn643EcI7naQnH9RA5HyY9ROQw23ie4d25Ps+VSFPXOa6HyCkw6SEih+04YU162LV1ZcObfz47T1aA2xwSyY9JDxE5xGIR2NG839aIeCY9VzIgJghabyXKa0w4VlotdzhEHo9JDxE55GhJFQx1jfBVq9Av+vK7GROg9lJicKwOALC9oELeYIiISQ8ROcY2nie5WzC8VHwLuRpbF+D2Ao7rIZIb37GIyCHZBRzP4wjbDC5blyARyYdJDxG1mhACO2xJD8fztMrQuGColAqcrqzD6co6ucMh8mhMeoio1U6eq0VJVQPUKiUGNY9VoSvz03ihf/PYJ1vCSETyYNJDRK22vXmq+qDYIGi9VTJH4zpsU9dtPz8ikgeTHiJqte3s2moTaVwPW3qIZMWkh4ha7cKVmKn1bD+vvJJqVNSYZI6GyHMx6SGiVtEb6lFYXgulwjpdnVovxE+NhOY9ynawi4tINkx6iKhVbONR+kUHIUDrLXM0rsfW2sOkh0g+THqIqFV2sGvrmoyIt7aObed6PUSyYdJDRK3CQczXxpYsHjxtQK2pSeZoiDwTkx4iuipDbSNyi6sAAMO7czxPW3QN9kV0kBZNFoHdhZVyh0PkkZj0ENFV5RRaW3l6hPmhi79G5mhcl23q+nZOXSeSBZMeIroq275Rw7uxa+tacDAzkbyY9BDRVe1s/pAexq6ta2IbD7WrsAKmJovM0RB5HiY9RHRFDU1m7D1lAAAM48yta5IQ5o9gX2/UN1pw8IxB7nCIPA6THiK6ogOnDTA1WRDqr0b3Lr5yh+PSlEoFkpu7CHNOcuo6UWdj0kNEV2QbzzOsWwgUCoXM0bg+Wxchx/UQdT4mPUR0RRzP075sU/53nqiAEELmaIg8C5MeIrosi0VgZ3M3DFdibh/9Y4Kg9lLiXI0JJ87Vyh0OkUdh0kNEl3WstBqVtY3w8Vahb3Sg3OG4BY2XCoO6BgFgFxdRZ2PSQ0SXZRvPMzhWB28V3y7ai20W3E4mPUSdiu9iRHRZtg9lbj3RvoZ1ax7XwxlcRJ2KSQ8RXdaOk7ZBzBzP056Sm5Oe46U1OFfdIHM0RJ6DSQ8RtajYWI+i8jooFcCQOJ3c4bgVna8aieH+ANjaQ9SZmPQQUYt2No/nSYoMRIDWW+Zo3I+t9YyLFBJ1HiY9RNSiHRzP06GGc5FCok7HpIeIWrST43k61LDm7SgOnDagvtEsczREnoFJDxFdorqhCYfOGAFwJeaOEhvig/AADRrNAnuLKuUOh8gjMOkhokvsLqyARQAxOh9EBfnIHY5bUigU0irXHMxM1DmY9BDRJaRNRtnK06FsU9e5SCFR52DSQ0SXyOF4nk5xYUuPxcLNR4k6GpMeIrLTZLZgd2ElAM7c6mh9ogLgq1ahqr4JR0uq5A6HyO0x6SEiO4fPVqHWZEaA1gu9wgPkDseteamU0sKPtnWRiKjjMOkhIju2qepD44KhVCpkjsb92aauc1wPUcdj0kNEdmwziWybYlLHsg0W5wwuoo7HpIeIJEIIqcWBg5g7x5C4YCgVwKmKOugN9XKHQ+TWmPQQkeRURR2KjQ3wUiowOFYndzgewV/jhT5RgQDOdy0SUcdg0kNEEtvml/2iA+GjVskcjecYJq3Xwy4uoo7EpIeIJNxvSx7J0no9bOkh6khMeohIYmtp4CDmzmVbD+nw2SrUNDTJHA2R+2pT0rNkyRJ0794dWq0WKSkp2L59+xXLr169GklJSdBqtRgwYADWrl1rd14Igfnz5yMqKgo+Pj5IS0tDXl6eXZmXXnoJI0eOhK+vL3Q6XYv3USgUlzxWrlzZlioSeRxDXSNyi60L5CVzUcJOFRXkgxidD8wWgT3cfJSowzic9KxatQpz587FggULsGvXLgwaNAjp6ekoKSlpsfzWrVsxefJkTJ8+Hbt370ZGRgYyMjJw4MABqcyiRYvw1ltvYenSpcjOzoafnx/S09NRX39+JoPJZMLEiRPx5z//+YrxLVu2DGfPnpUeGRkZjlaRyCPtLqyAEEBciC/CA7Ryh+Nxkjmuh6jDOZz0vP7665gxYwamTZuGvn37YunSpfD19cXHH3/cYvk333wT48ePx1//+lf06dMH//znPzF06FC88847AKytPG+88Qaee+453HXXXRg4cCA++eQTnDlzBt988430On//+9/xl7/8BQMGDLhifDqdDpGRkdJDq+WbN1Fr2AYxc5NReZxfr4fjeog6ikNJj8lkQk5ODtLS0s6/gFKJtLQ0ZGVltXhNVlaWXXkASE9Pl8oXFBRAr9fblQkKCkJKSsplX/NKZs+ejdDQUIwYMQIff/wxhLj8Jn4NDQ0wGo12DyJPtcO2Pk83DmKWg+3nvruwEmZuPkrUIRxKesrKymA2mxEREWF3PCIiAnq9vsVr9Hr9Fcvbvjrympfzj3/8A19++SUyMzNx77334uGHH8bbb7992fILFy5EUFCQ9IiNjXXofkTuotFskcaSsKVHHr0jAxCg8UJ1QxOO6PkHGFFH8JI7gPb0/PPPS98PGTIENTU1WLx4MR577LEWy8+bNw9z586VnhuNRiY+5JEOnTGivtGCIB9vJIT5yx2OR1IpFRgcp8OWvDLknKxAv+gguUMicjsOtfSEhoZCpVKhuLjY7nhxcTEiIyNbvCYyMvKK5W1fHXnN1kpJScGpU6fQ0NDQ4nmNRoPAwEC7B5EnsnVtJXfjJqNyGt68Xs8ODmYm6hAOJT1qtRrJycnYsGGDdMxisWDDhg1ITU1t8ZrU1FS78gCQmZkplY+Pj0dkZKRdGaPRiOzs7Mu+Zmvt2bMHwcHB0Gg01/Q6RO7ONog5mevzyMq2PlIOd1wn6hAOd2/NnTsXDzzwAIYNG4YRI0bgjTfeQE1NDaZNmwYAmDp1KmJiYrBw4UIAwJw5czBmzBi89tprmDBhAlauXImdO3figw8+AGBdW+fxxx/Hiy++iMTERMTHx+P5559HdHS03XTzwsJClJeXo7CwEGazGXv27AEAJCQkwN/fH99//z2Ki4tx3XXXQavVIjMzE//v//0/PPnkk9f4IyJyb0II7qzuJAbH6aBSKnDGUI/TlXWI0fnIHRKRW3E46Zk0aRJKS0sxf/586PV6DB48GOvWrZMGIhcWFkKpPN+ANHLkSHz++ed47rnn8OyzzyIxMRHffPMN+vfvL5V56qmnUFNTg5kzZ6KyshKjRo3CunXr7Kabz58/HytWrJCeDxkyBADw888/Y+zYsfD29saSJUvwl7/8BUIIJCQkSNPriejyCstrUVrVAG+VAoO4yaisfNVe6BcdiH2nDNh5ohwxg2PkDonIrSjEleZ0exij0YigoCAYDAaO7yGP8e+cU3hi9V4MjdPhPw9fL3c4Hu/v3x/Est9OYGpqN/zjrv5Xv8DV1NQA/s2D5aurAT8/eeMht9Daz2/uvUXk4WyL4Q3nJqNOwbZeD1dmJmp/THqIPJxtphB3VncOtnWSjuiNqKpvlDkaIvfCpIfIg5XXmJBfUg2AM7ecRUSgFrEhPrAI6+rMRNR+mPQQeTDbVPWeYX4I8VPLHA3ZDJe6uDh1nag9Mekh8mAcz+OchnGRQqIOwaSHyIPt5HgepzS8eVzP7qIKNJotMkdD5D6Y9BB5qPpGM/adqgRw/kOWnEPPMH/ofL1R32jBwTPcfJSovTDpIfJQ+04Z0GgWCAvQIC7EV+5w6AJKpeKCqesc10PUXpj0EHko2yajw7sHQ6HgJqPOxtb6toNJD1G7YdJD5KF2SjurczyPM7KNs9p5ogJcOJ+ofTDpIfJAFouQpqtzPI9z6h8TCI2XEudqTDheViN3OERugUkPkQfKK6mGsb4JvmoV+kZxnzlnpPFSSRvAclwPUftg0kPkgWzjRIbE6eCl4tuAszo/rofr9RC1B77bEXkgW8vBMI7ncWrnx/WwpYeoPTDpIfJAtpYDrsTs3IbGBUOhAE6cq0VJVb3c4RC5PCY9RB7mrKEOpyvroFQAg+N0codDVxDk443eEQEAgBx2cRFdMyY9RB7GtvVE3+hA+Gu8ZI6GrsbWGredXVxE14xJD5GH2cHxPC5lePz59XqI6Now6SHyMNsLrEnPiHgmPa7ANoPr4BkDqhuaZI6GyLUx6SHyIIbaRuQWVwHgIGZXERXkgxidDywC2FNYKXc4RC6NSQ+RB9l5shxCAD1C/RAWoJE7HGol7sNF1D6Y9BB5kO3SJqNs5XEltvV6mPQQXRsmPUQehON5XJPt97WrsAKmJovM0RC5LiY9RB6izmTG/lMGAEx6XE1CmD+Cfb1R32jBgTMGucMhcllMeog8xO7CCjRZBKKCtOga7CN3OOQApVJxfr2eAnZxEbUVkx4iD3HheB6FQiFzNOQoW+sckx6itmPSQ+QhOJ7HtaXEdwFgHcxstgiZoyFyTUx6iDyAqcmCXYXWFX2Z9LimPlEB8Nd4oaq+CUf0RrnDIXJJTHqIPMCBMwbUN1oQ7OuNhDB/ucOhNvBSKZHczbpeD7u4iNqGSQ+RB9jR/CE5rHsIlEqO53FVHNdDdG2Y9BB5ANuHZAq7tlxaygVJjxAc10PkKCY9RG7OYhHSSr5cidm1DegaBI2XEudqTDhWWiN3OEQuh0kPkZvLLa6Csb4JvmoV+kUHyh0OXQONlwpD4nQA2MVF1BZMeojcnK2VJ7lbMLxU/C/v6kY0T13fXnBO5kiIXA/fAYncXLZtfR52bbkF27iebI7rIXIYkx4iNyaEkGZuDecgZrcwJE4HL6UCZw31OFVRJ3c4RC6FSQ+RGzt5rhYlVQ1Qq5QYHKuTOxxqB75qLwzoGgSA43qIHMWkh8iNbTtuHfcxOFYHrbdK5miovXC9HqK2YdJD5MZsSc91Pdi15U6k9XpOMOkhcgSTHiI3JYTAtuPWD8XrenaRORpqT8ndQqBQAAVlNSgx1ssdDpHLYNJD5KZOnKuF3lgPtUqJoXHBcodD7SjIxxt9Iq1rLmWzi4uo1Zj0ELkpaTxPHMfzuKOU5i5L2++ZiK6OSQ+Rmzo/noddW+4otfn3msWkh6jVmPQQuSHreB4OYnZnKfFdoFAAx0trUMxxPUStwqSHyA2dOFeLYmMDx/O4sSBfb/SNso7rYRcXUesw6SFyQxzP4xmkLq5jTHqIWoNJD5EbsiU9qRzP49ZSm5ciYEsPUesw6SFyM/bjeZj0uLPh8SFQKqzdmWcN3IeL6GqY9BC5GWk8j5cSQ+J0codDHShQ640BMdZ9uNjFRXR1THqI3Iztw28I99vyCLbVtpn0EF0dkx4iN8OuLc/C9XqIWo9JD5Eb4XgezzO8ewi8lAqcqqhDUXmt3OEQOTUmPURupKCsBiVVHM/jSfw0XhjYtXlcD1t7iK6ISQ+RG7Htqs7xPJ7F1qq3jeN6iK6ISQ+RG8li15ZHsq3Xk3X8HIQQMkdD5LzalPQsWbIE3bt3h1arRUpKCrZv337F8qtXr0ZSUhK0Wi0GDBiAtWvX2p0XQmD+/PmIioqCj48P0tLSkJeXZ1fmpZdewsiRI+Hr6wudTtfifQoLCzFhwgT4+voiPDwcf/3rX9HU1NSWKhK5HItFYGt+GQDg+oRQmaOhzjSsWwi8VQqcNdSjkON6iC7L4aRn1apVmDt3LhYsWIBdu3Zh0KBBSE9PR0lJSYvlt27dismTJ2P69OnYvXs3MjIykJGRgQMHDkhlFi1ahLfeegtLly5FdnY2/Pz8kJ6ejvr685vomUwmTJw4EX/+859bvI/ZbMaECRNgMpmwdetWrFixAsuXL8f8+fMdrSKRS8otrsK5GhN81SoMjtXJHQ51Ip8Lfuecuk50BcJBI0aMELNnz5aem81mER0dLRYuXNhi+fvuu09MmDDB7lhKSoqYNWuWEEIIi8UiIiMjxeLFi6XzlZWVQqPRiC+++OKS11u2bJkICgq65PjatWuFUqkUer1eOvbee++JwMBA0dDQ0Kq6GQwGAUAYDIZWlSdyJh9uPia6Pb1GPPhxttyhkAxe+/GI6Pb0GvHYF7vkDuXKqquFAKyP6mq5oyE30drPb4daekwmE3JycpCWliYdUyqVSEtLQ1ZWVovXZGVl2ZUHgPT0dKl8QUEB9Hq9XZmgoCCkpKRc9jUvd58BAwYgIiLC7j5GoxEHDx5s8ZqGhgYYjUa7B5Gr+pVdWx7twkUKBcf1ELXIoaSnrKwMZrPZLrEAgIiICOj1+hav0ev1Vyxv++rIazpynwvvcbGFCxciKChIesTGxrb6fkTOxNRkQXbzzC0mPZ5paFwwNF5KlFQ1IK+kWu5wiJySR8/emjdvHgwGg/QoKiqSOySiNtldWIG6RjNC/dXoHREgdzgkA623CiPiQwAAW/LKZI6GyDk5lPSEhoZCpVKhuLjY7nhxcTEiIyNbvCYyMvKK5W1fHXlNR+5z4T0uptFoEBgYaPcgckW/NXdtjewZCqVSIXM0JJdRza18v+aVyhwJkXNyKOlRq9VITk7Ghg0bpGMWiwUbNmxAampqi9ekpqbalQeAzMxMqXx8fDwiIyPtyhiNRmRnZ1/2NS93n/3799vNIsvMzERgYCD69u3b6tchckW28Tyj2LXl0UYlWn//2QXlMDVZZI6GyPl4OXrB3Llz8cADD2DYsGEYMWIE3njjDdTU1GDatGkAgKlTpyImJgYLFy4EAMyZMwdjxozBa6+9hgkTJmDlypXYuXMnPvjgAwCAQqHA448/jhdffBGJiYmIj4/H888/j+joaGRkZEj3LSwsRHl5OQoLC2E2m7Fnzx4AQEJCAvz9/TFu3Dj07dsX999/PxYtWgS9Xo/nnnsOs2fPhkajucYfE5HzqqpvxN5TBgDA9YlMejxZn8hAdPFT41yNCbsKK7hIJdFFHE56Jk2ahNLSUsyfPx96vR6DBw/GunXrpEHDhYWFUCrPNyCNHDkSn3/+OZ577jk8++yzSExMxDfffIP+/ftLZZ566inU1NRg5syZqKysxKhRo7Bu3TpotVqpzPz587FixQrp+ZAhQwAAP//8M8aOHQuVSoU1a9bgz3/+M1JTU+Hn54cHHngA//jHPxz/qRC5kOzj5TBbBOJD/RCj85E7HJKRUqnA9Qmh+G7vGfyaV8akh+giCsG5jRKj0YigoCAYDAaO7yGX8cJ3B7F86wn88bo4vJgxQO5wSGZf7izCU1/tw6BYHb6dfb3c4Vyqpgbw97d+X10N+PnJGw+5hdZ+fnv07C0id/Abx/PQBUY3d3HuP1UJQ22jzNEQORcmPUQurNhYj7ySaigUQGoPJj0ERAX5oGeYHywC2HqMU9eJLsSkh8iF2Vp5BsQEIcjXW+ZoyFmMTgwDAGzJZ9JDdCEmPUQujFtPUEvOr9fDpIfoQkx6iFyUEILjeahF1/XsAi+lAoXltSg8Vyt3OEROg0kPkYs6VlqNYmMDNF5KJHcLljscciL+Gi8MidMBALbkc3VmIhsmPUQu6pdc64fZiPgQaL1VMkdDzmZUgnVcD7u4iM5j0kPkojYdtSY9Y3qFyRwJOSPblhRbj52D2cLl2IgAJj1ELqnW1ITs4+UAgLG9w2WOhpzRoK5BCNB6wVDXiP2nDXKHQ+QUmPQQuaCsY+dgMlvQNdi6JgvRxbxUSqQ2b0Ox5SjH9RABTHqIXJJtPM/Y3mFQKBQyR0POakxva9fnL0x6iAAw6SFyOUII/HK0BAAwthe7tujybF2fuwsrUFFjkjkaIvkx6SFyMcfLalBUXge1SomRCdxFmy4vRueD3hEBsAhgcx5be4iY9BC5mAunqvuqvWSOhpzdjUnW1p6fj5TIHAmR/Jj0ELmYX3Kbu7Z6c6o6Xd2Nzf9ONh0t5dR18nhMeohcSJ3JjOwC21R1Jj10dUO7BSNA64WK2kbsPVUpdzhEsmLSQ+RCso6XwdRkQYzOBz3D/OUOh1yAt0qJG5oXsPyFXVzk4Zj0ELmQTc3jecZwqjo54MbmWVwbc5n0kGdj0kPkQmzrrYzl1hPkANtWJQdOG1FirJc5GiL5MOkhchEFZTU4ea4W3ioFRiaEyh0OuZCwAA0GdQ0CwIUKybMx6SFyEbZZW8O7h8Bfw6nq5BjbQoWcuk6ejEkPkYu4cOsJIkfZ1uvZkleGRrNF5miI5MGkh8gFVDc0IevYOQDATUnceoIcNzAmCF381KhuaMLOExVyh0MkCyY9RC5gU24pTGYLeoT6cao6tYlSqZA2IP2Zs7jIQzHpIXIBPx3SAwBu6RvBqerUZjdyXA95OCY9RE6u0WzBxuYPqVv6RsgcDbmyG3qFQaVUIK+kGifP1cgdDlGnY9JD5OSyj5ejqr4Jof5qDIkLljsccmFBPt5IiQ8BAPx0sFjmaIg6H5MeIieX2dy1dXNSBFRKdm3RtUnvFwkA+PGgXuZIiDofkx4iJyaEQOYh61/k4/qxa4uune3fUU5hBUqquDozeRYmPURO7OAZI84Y6uHjrcL1XIWZ2kFUkA8GxeogBKSEmshTMOkhcmI/NX8o3dArFFpvlczRkLtIb27t+ZHjesjDMOkhcmI/NY+7GNc3UuZIyJ3YxvVkHSuDsb5R5miIOg+THiInVVReiyP6KqiUCq7CTO2qZ5g/EsL90WgWXLOHPAqTHiInZevaGt49GMF+apmjIXdzvouLs7jIczDpIXJSmdIqzOzaovZn6+L6+Ugp6hvNMkdD1DmY9BA5oYoaE3Y0bwo5jqswUwcYEBOE6CAt6hrN2JJXJnc4RJ2CSQ+RE1p/uBhmi0BSZABiQ3zlDofckEKhwDguVEgehkkPkRP6ft9ZAMBtA6JkjoTcma2La8PhYjSZLTJHQ9TxmPQQOZnyGhN+y7d2N9w+kEkPdZzh3YMR7OuNitpGbD9RLnc4RB2OSQ+Rk/nvgbMwWwT6xwSiR5i/3OGQG/NSKZHWxzpmbN0BdnGR+2PSQ+Rkvt97BgBw+8BomSMhT3Bbc2vi2v1n2cVFbo9JD5ETKTHWI7vA2s0wgeN5qBOMSghFiJ8aZdUmbD12Tu5wiDoUkx4iJ/LD/rMQAhgap+OsLeoU3iolbhtgHdD87Z4zMkdD1LGY9BA5kTXNs7bYtUWd6a7BMQCsU9e5UCG5MyY9RE7iVEUtck5WQKEAJnDWFnWi5LhgxOh8UN3QhI3ci4vcGJMeIifxQ3Mrz4juIYgI1MocDXkSpVKBOwZZWxe/YxcXuTEmPUROwta1ZfvwIepMdw22/rvbmFsCQ12jzNEQdQwmPUROoKCsBvtPG6BSKnBrf24wSp0vKTIAieH+MDVZuC0FuS0mPUROYE3z2jwje3ZBF3+NzNGQJ1IoFFJrD7u4yF0x6SGSmRAC3zUnPezaIjndOcg6i2vrsTKUVNXLHA1R+2PSQySzfacMyCuphsZLKW0ASSSHuC6+GBKng0WcH1hP5E6Y9BDJ7MudRQCAW/tHIsjHW+ZoyNPd1dzayIUKyR0x6SGSUZ3JLI2fmDgsVuZoiIAJA6OhVAB7iipRUFYjdzhE7YpJD5GMfjyoR1VDE2J0Pkjt0UXucIgQFqDBDb3CAACrdhTJHA1R+2LSQySj1TnWD5WJw7pCqVTIHA2R1e+HxwEAvso5hUbuvE5upE1Jz5IlS9C9e3dotVqkpKRg+/btVyy/evVqJCUlQavVYsCAAVi7dq3deSEE5s+fj6ioKPj4+CAtLQ15eXl2ZcrLyzFlyhQEBgZCp9Nh+vTpqK6uls6fOHECCoXikse2bdvaUkWiDldUXovf8s9BoQB+l9xV7nCIJDf3CUeovwZl1Q3YcJjbUpD7cDjpWbVqFebOnYsFCxZg165dGDRoENLT01FS0vJ/jK1bt2Ly5MmYPn06du/ejYyMDGRkZODAgQNSmUWLFuGtt97C0qVLkZ2dDT8/P6Snp6O+/vyUySlTpuDgwYPIzMzEmjVrsHnzZsycOfOS+61fvx5nz56VHsnJyY5WkahTfJVzCgBwfc9QdA3mjurkPLxVSikRX7mjUOZoiNqRcNCIESPE7Nmzpedms1lER0eLhQsXtlj+vvvuExMmTLA7lpKSImbNmiWEEMJisYjIyEixePFi6XxlZaXQaDTiiy++EEIIcejQIQFA7NixQyrz3//+VygUCnH69GkhhBAFBQUCgNi9e7ejVZIYDAYBQBgMhja/BlFrmM0WMXLhBtHt6TXim92n5A6H6BLHS6tFt6fXiO7PrBGnK2rb74Wrq4UArI/q6vZ7XfJorf38dqilx2QyIScnB2lpadIxpVKJtLQ0ZGVltXhNVlaWXXkASE9Pl8oXFBRAr9fblQkKCkJKSopUJisrCzqdDsOGDZPKpKWlQalUIjs72+6177zzToSHh2PUqFH47rvvrlifhoYGGI1GuwdRZ9h67BxOV9YhUOvFtXnIKcWH+uG6HiEQ4vyyCkSuzqGkp6ysDGazGREREXbHIyIioNe3vFeLXq+/Ynnb16uVCQ8Ptzvv5eWFkJAQqYy/vz9ee+01rF69Gj/88ANGjRqFjIyMKyY+CxcuRFBQkPSIjeWUYeoctg+RuwbHQOutkjkaopZNHmEd0Lx65ymYLULmaIiunZfcAbSX0NBQzJ07V3o+fPhwnDlzBosXL8add97Z4jXz5s2zu8ZoNDLxoQ5nqG3EuuYNHe/j2jzkxNL7WRfMPF1Zhy15pRjbO/zqFxE5MYdaekJDQ6FSqVBcXGx3vLi4GJGRLTfRR0ZGXrG87evVylw8ULqpqQnl5eWXvS8ApKSkID8//7LnNRoNAgMD7R5EHe2bPadharIgKTIA/WP4b46cl9ZbhbuHWPfjWrmdXVzk+hxKetRqNZKTk7FhwwbpmMViwYYNG5CamtriNampqXblASAzM1MqHx8fj8jISLsyRqMR2dnZUpnU1FRUVlYiJydHKrNx40ZYLBakpKRcNt49e/YgKirKkSoSdSiLRWDF1hMArF0HCgXX5iHn9vsR1tbI9YeLUVrVIHM0RNfG4e6tuXPn4oEHHsCwYcMwYsQIvPHGG6ipqcG0adMAAFOnTkVMTAwWLlwIAJgzZw7GjBmD1157DRMmTMDKlSuxc+dOfPDBBwAAhUKBxx9/HC+++CISExMRHx+P559/HtHR0cjIyAAA9OnTB+PHj8eMGTOwdOlSNDY24pFHHsHvf/97REdb94lZsWIF1Go1hgwZAgD4z3/+g48//hgfffTRNf+QiNrLprxSHC+rQYDGC/dybR5yAUmRgRgcq8Oeokr8e9cp/GlMT7lDImozh5OeSZMmobS0FPPnz4der8fgwYOxbt06aSByYWEhlMrzDUgjR47E559/jueeew7PPvssEhMT8c0336B///5Smaeeego1NTWYOXMmKisrMWrUKKxbtw5arVYq89lnn+GRRx7BzTffDKVSiXvvvRdvvfWWXWz//Oc/cfLkSXh5eSEpKQmrVq3C7373O4d/KEQdZflvJwAA9w2Phb/GbYbUkZubPCIWe4oq8X/bTmLG6B5QcfVwclEKIQSH5DczGo0ICgqCwWDg+B5qd/kl1Uh7fRMUCmDTkzcirgsXJCTXUN9oRurCDaiobcR7U4bi1gHXMGygpgbw97d+X10N+Pm1T5Dk0Vr7+c29t4g6ySdZJwAANydFMOEhl6L1VuH+67oBAD7cclzmaIjajkkPUScw1DVK2048dH13eYMhaoM/pnaDWqXErsJK5JyskDscojZh0kPUCVbvLEKtyYzeEQFI7dlF7nCIHBYeoMVdg60TRz7+tUDmaIjahkkPUQczWwRWNHdtPXh9d05TJ5c1fXQ8AOC/B86iqLxW5miIHMekh6iDbThcjKLyOuh8vZExOEbucIjaLCkyEKMTQ2ERwPLm9aaIXAmTHqIOtqx5mvrvh8fBR819tsi1TR9lbe1ZtaMIxvpGmaMhcgyTHqIOtLuwAlnHz0GlVOD+1G5yh0N0zcb0CkNiuD+qG5rw5Q5uTUGuhUkPUQd6a0MeAOCeITGI0fnIHA3RtVMoFFJrz7LfTqDJbJE5IqLWY9JD1EH2FlXi59xSqJQKzL4xQe5wiNpNxpAYdPFT43RlHb7be0bucIhajUkPUQd5e6O1leeuwdHoHspVZ8l9aL1VeKi5tefNDXls7SGXwaSHqAMcOG3A+sMlUCqAR9jKQ27owZHd0cVPjZPnavGfXaflDoeoVZj0EHWAN5vH8tw5KBo9wvxljoao/flpvKQd19/ckAdTE1t7yPkx6SFqZwfPGJB5qBgKBfDITYlyh0PUYf54XTeEBWhwurIOX+7kTC5yfkx6iNrZ2xvyAQC3D4xGQjhbech9+ahVmD3W2trzzsZ81DeaZY6I6MqY9BC1o8NnjVh3UA+FAnj0Jo7lIff3+xFxiArSQm+sxxfbC+UOh+iKmPQQtaPFP+YCAG7rH4VeEQEyR0PU8bTeKjzSnOAv+fkY6kxs7SHnxaSHqJ1sOlqKjUdK4KVUYO64XnKHQ9RpJibHomuwD8qqG/DpthNyh0N0WUx6iNpBk9mCF9ccAgA8MLI7enLGFnkQtZcSj91sHbS/5OdjKK8xyRwRUcuY9BC1g8+yC5FXUo1gX288xhlb5IHuGRKDPlGBMNQ1YvGPR+QOh6hFTHqIrlFlrQn/Wn8UADB3XG8E+XrLHBFR5/NSKfGPu/oBAFbuKMLeokp5AyJqAZMeomv0xvo8VNY2ondEACYPj5U7HCLZDO8egruHxEAIYP53B2GxCLlDIrLDpIfoGuSXVOHTbScBAM/f3hdeKv6XIs8279Yk+KlV2FtUia9yTskdDpEdvkMTtZEQAi/+cBhmi0BanwiMSgyVOyQi2YUHavF4mnX24ivrjsBQ2yhzRETnMekhaqM1+87il9xSeKsU+NuEPnKHQ+Q0Hry+OxLC/XGu5vx4NyJnwKSHqA1Kqurx/LcHAAAPj01AfKifzBEROQ9vlRIv3GEd1PxJ1gkcOG2QOSIiKyY9RA4SQuBvXx9AZW0j+kYFSqvREtF5oxJDMWFAFCwCeOLLvdyXi5wCkx4iB3275wwyDxXDW6XAqxMHwZuDl4la9Pe7+iHUX43c4iq82rxFC5Gc+G5N5IASYz0WfHcQAPDYTYnoGx0oc0REzivUX4NX7h0IAPjo1wJszS+TOSLydEx6iFpJCIFnv94PQ10jBsQE4U9je8odEpHTu7lPBCaPiAMAPLF6Lwx1nM1F8mHSQ9RKX+WcwvrDJVCrlOzWInLAcxP6oHsXX5w11Et71BHJge/aRK1w8IxBmq01Jy0RvSMDZI6IyHX4abzw+qTBUCkVWLPvrNzhkAdj0kN0FRU1Jsz6NAf1jRaM6RWGP41htxaRo4bGBWP2jZzpSPJi0kN0BWaLwGMrd+NURR3iQnzx1u+HQKVUyB0WkUt69KYEDOsWLD2vqDHJGA15IiY9RFew+MdcbMkrg4+3Ch9MTeYO6kTXwFulxBu/Hyw9n/vlHjSaLfIFRB6HSQ/RZazZdwZLNx0DACyeOBBJkZyeTnStuvhrpO+3HS/HP77nwGbqPEx6iFqwp6gST321DwAw64YeuH1gtMwREbkfhQL4dNtJfJp1Qu5QyEMw6SG6yKEzRkz932zUmswYnRiKv6b3ljskIrf0l+bd2F/4/hB+zi2RORryBEx6iC6QX1KN+/83G8b6JiR3C8bSPybDi+vxEHWI/xkdj3uGxMBsEZj1aQ42HS2VOyRyc3w3J2pWVF6LP36UjXM1JvSPCcTHDw6Hn8ZL7rCI3JZCocDL9w7ELX0jYGqyYMYnO7GZiQ91ICY9RAD0hnr84aNt0BvrkRjuj08eSkGQD2dqEXU0tZcSS/4w1C7x2ZLHxIc6BpMe8nj5JdWY+P5WFJXXoVsXX/zf/6QgxE8td1hEHsOW+KT1iUBDkwX/s2Infs3j5qTU/pj0kEfbdvwc7nn3NxSVWxcf/L/pKYgI1ModFpHHUXsp8e6UoUjrE46GJgseWr4D/845JXdY5GaY9JDH+nr3KWnQ8pA4Hb5+eCRiQ3zlDovIY6m9lFgyZShu7R8Jk9mCJ1bvxcL/HobZIuQOjdwEkx7yOEIIvLUhD39ZtReNZoHbBkTiixnX2S2aRkTy0HipsOQPQ/HoTdZ9ut7fdBwzPtmJqvpGmSMjd8CkhzxKSVU9pi3fgdczjwKwLjz4zuSh0HqrZI6MiGyUSgWeGNcbb00eAo2XEhuPlOCed7fieGm13KGRi2PSQx5j/aFi3PrGFvySWwq1lxL/7+4BmHdbHyi5gSiRU7pzUDS+nJWK8AAN8kqqcdtbW7DstwJY2N1FbcSkh9xerakJz369H//zyU6cqzEhKTIAax4dhT+kxMkdGhFdxaBYHb5/dBRGJYSivtGCv39/CJM/3IbCc7Vyh0YuiEkPuS0hBP67/yzS39iMz7MLAQD/Myoe3z5yPXpFBMgcHRG1VkSgFp9OH4EXM/rDV61CdkE5xr+5GZ9knUATd2knByiEEGwnbGY0GhEUFASDwYDAQO6o7coOnDbgn2sOIbugHAAQGajFqxMHYVRiqMyREXm4mhrA39/6fXU14Ofn0OWF52rx16/2Sv+3E8P98fT4JNzcJxwKBbuqPVVrP7+Z9FyASY/rKyqvxTsb8/FlThGEADReSsy6oQf+NLYnfNXcUoJIdteY9ACAxSLw6baT+Nf6o6istc7qGtE9BM/cloShccHtGS25CCY9bcCkx3UdOG3A+5uPY+3+s9KaHncMisYztyYhRucjc3REJGmHpMfGUNeIpZuO4eNfC9DQZO3mGtMrDNOu744bEsM4ScGDMOlpAyY9rqXRbMGm3FJ8/FsBth47Jx0fnRiKOTcnYlj3EBmjI6IWtWPSY3PWUId/ZR7FVzmnYJvY1SPMDw+O7I57h3blxsEegElPGzDpcX5CCOwpqsQ3u0/j+31nUV5jAgColArcMTAKM27ogX7RQTJHSUSX1QFJj83JczVYsfUkVu8sQlVDEwDAX+OFcX0jcPugKIxKCIPai/N33BGTnjZg0uOcGprM2HmiApuOliLzUDEKymqkc6H+Gtw9JBoPXh/PbiwiV9CBSY9NdUMTvtpZhBVZJ+3eL4J8vDG+XyTG9YtASo8u8GcLkNtg0tMGTHqcg6nJgsNnjdhVWIHf8suw9dg51JrM0nkfbxXS+0UgY0gMRiWEwkvFv9yIXEYnJD02FovArsIKrNl3Fj/sP4vSqgbpnJdSgcGxOlyfEIrrE0IxICYIPmquzO6qOjTpWbJkCRYvXgy9Xo9Bgwbh7bffxogRIy5bfvXq1Xj++edx4sQJJCYm4pVXXsFtt90mnRdCYMGCBfjwww9RWVmJ66+/Hu+99x4SExOlMuXl5Xj00Ufx/fffQ6lU4t5778Wbb74Jf9t/HgD79u3D7NmzsWPHDoSFheHRRx/FU0891ep6MenpfHUmM46VVuNocRUOnzVid2El9p82SIMSbcICNBjTKwxjeoXhxqRw/oVG5Ko6Mem5kNkikF1wDmv3n8WWvDKcvGhxQ5VSgcRwfwzqqsPA2CD0iQpEQrg/ArXenRIfXZvWfn47/MmxatUqzJ07F0uXLkVKSgreeOMNpKenIzc3F+Hh4ZeU37p1KyZPnoyFCxfi9ttvx+eff46MjAzs2rUL/fv3BwAsWrQIb731FlasWIH4+Hg8//zzSE9Px6FDh6DVagEAU6ZMwdmzZ5GZmYnGxkZMmzYNM2fOxOeffy5VeNy4cUhLS8PSpUuxf/9+PPTQQ9DpdJg5c6aj1aR2YrEIVNSacNZQj1MVdThVUYtTFXUoKq9Ffmk1Cstr0VLaHeTjjSFxOgzvHoKxvcPQJzKQMzGIqM1USgVG9gzFyJ7WtbqKymvxa34Zfs0vQ/bxcpRVN+CIvgpH9FVYtbNIui4sQIOEMH/0CPND12BfxAT7IEanRYzOF2EBGqj4vuRSHG7pSUlJwfDhw/HOO+8AACwWC2JjY/Hoo4/imWeeuaT8pEmTUFNTgzVr1kjHrrvuOgwePBhLly6FEALR0dF44okn8OSTTwIADAYDIiIisHz5cvz+97/H4cOH0bdvX+zYsQPDhg0DAKxbtw633XYbTp06hejoaLz33nv429/+Br1eD7VaDQB45pln8M033+DIkSOtqhtbei4lhEBDkwV1JjPqGq2PmoYmVNc3obrB+qiqb0JlbSMq60ww1DaiotaEsmoTSqsaUFbdgKar7JMT4qdGYrg/ekcGYFBXHYbE6RAf6seFxojckUwtPVcihIDeWI+9RQbsP12JfacMOFpchWJjwxWvUyqs71+h/hp08bd+1fl4I8jHG4HNXwO03vDXeMFPo4Kfxgt+Gi9ovZTQequg9VYxaWonHdLSYzKZkJOTg3nz5knHlEol0tLSkJWV1eI1WVlZmDt3rt2x9PR0fPPNNwCAgoIC6PV6pKWlSeeDgoKQkpKCrKws/P73v0dWVhZ0Op2U8ABAWloalEolsrOzcffddyMrKws33HCDlPDY7vPKK6+goqICwcGXLljV0NCAhobz/6iNRqMjP45WW3+oGFvySi85bksFbGmnaD5yYRoqpOcCQli/F7bvYf/cIuy/mi0CFmF7AE0WAbPFArNFwGwRaDTbvlrQZBFoMltgarLAZLagocn6/cXdTG0V6q9GTLAvYoN9EBvii67BPogP9UOviACE+mva5R5ERG2hUCgQFeSDqCAfjO8fKR2vqm/E8dIaHCutxvHSGpyprMOpyjqcqayD3lCPJotAWbX1j7y28lYpoPVSwdtLCbVKCbWXEt4qBbxVSnipFPBSKuGlVMBLpYBKqYBSoYCX8vz3SoUCSqW1DtbngAKAUqGAQqGAovm59WvzcwUAnD8H2I5Zy1z43O7n1MLPrS3S+kTItjq+Q0lPWVkZzGYzIiIi7I5HRERctjVFr9e3WF6v10vnbceuVObirjMvLy+EhITYlYmPj7/kNWznWkp6Fi5ciL///e+Xr3A72VVYgRVZJzv8Ph1NrVJC462Ev8ar+S8XL+n7YD9v6HzV0Pl4Q+frjVB/DcICNAgP0KKLvxreHGxMRC4mQOuNQbE6DIrVXXLObBE4V9OAsioTyqqtrdrnqk0w1DW3etc1wVDXiKr6RtQ2mFFjakJNQxNqGswwXbBfWKNZoNHcBFy5UcmthAdqXCPpcTfz5s2za4UyGo2IjY1t9/uk9uxi14RplxtflClfmHVfmHHbHW++RnlBFm/73vZVpbRm+armzF+psP5l4KVSnP9eaf/XhLdKAbWX9S8N218cGi8VfNQqaL2UnCVFRNRMpVQgPECL8ACtw9eaLQINTWbUN1pQ32hGQ5MFjebzreumJmuLfKPFArNZoMliQaPZ2mpva6m3tuTbWvabjwlrVx2A5rLnewIA67nzPQuX72WQOhsu6Ha4eJBCawbGiEuuspJzqxCHkp7Q0FCoVCoUFxfbHS8uLkZkZGSL10RGRl6xvO1rcXExoqKi7MoMHjxYKlNSUmL3Gk1NTSgvL7d7nZbuc+E9LqbRaKDRdHzXyujEMIxODOvw+xARkfNTKRXwVXvBV331stS+HPrTXa1WIzk5GRs2bJCOWSwWbNiwAampqS1ek5qaalceADIzM6Xy8fHxiIyMtCtjNBqRnZ0tlUlNTUVlZSVycnKkMhs3boTFYkFKSopUZvPmzWhsbLS7T+/evVvs2iIiIiIPIxy0cuVKodFoxPLly8WhQ4fEzJkzhU6nE3q9XgghxP333y+eeeYZqfxvv/0mvLy8xKuvvioOHz4sFixYILy9vcX+/fulMi+//LLQ6XTi22+/Ffv27RN33XWXiI+PF3V1dVKZ8ePHiyFDhojs7Gzx66+/isTERDF58mTpfGVlpYiIiBD333+/OHDggFi5cqXw9fUV77//fqvrZjAYBABhMBgc/bEQEVFrVFdLczJEdbXc0ZCbaO3nt8NJjxBCvP322yIuLk6o1WoxYsQIsW3bNuncmDFjxAMPPGBX/ssvvxS9evUSarVa9OvXT/zwww925y0Wi3j++edFRESE0Gg04uabbxa5ubl2Zc6dOycmT54s/P39RWBgoJg2bZqoqqqyK7N3714xatQoodFoRExMjHj55ZcdqheTHiKiDsakhzpAaz+/uQ3FBbhODxFRB3PCdXrI9bX285vTcYiIiMgjMOkhIiIij8Ckh4iIiDwCkx4iIiLyCEx6iIiIyCMw6SEiIiKPwKSHiIiIPAKTHiIiIvIITHqIiIjIIzi0y7q7sy1ObTQaZY6EiMhN1dSc/95oBMxm+WIht2H73L7aJhNMei5QVVUFAIiNjZU5EiIiDxAdLXcE5GaqqqoQFBR02fPce+sCFosFZ86cQUBAABQKhdzhtAuj0YjY2FgUFRV5xH5irK97Y33dG+vr3jqyvkIIVFVVITo6Gkrl5UfusKXnAkqlEl27dpU7jA4RGBjoEf+pbFhf98b6ujfW1711VH2v1MJjw4HMRERE5BGY9BAREZFHYNLj5jQaDRYsWACNRiN3KJ2C9XVvrK97Y33dmzPUlwOZiYiIyCOwpYeIiIg8ApMeIiIi8ghMeoiIiMgjMOkhIiIij8Ckx4m9/PLLUCgUePzxx6Vjer0e999/PyIjI+Hn54ehQ4fi3//+9yXX/vDDD0hJSYGPjw+Cg4ORkZFhd76wsBATJkyAr68vwsPD8de//hVNTU12ZX755RcMHToUGo0GCQkJWL58eQfU8ry21vfo0aO46667EBoaisDAQIwaNQo///yzXRlnrC/Qcp2PHTuGu+++G2FhYQgMDMR9992H4uJiu+vKy8sxZcoUBAYGQqfTYfr06aiurrYrs2/fPowePRparRaxsbFYtGjRJfdfvXo1kpKSoNVqMWDAAKxdu7ZD6mnTlvqeOHEC06dPR3x8PHx8fNCzZ08sWLAAJpPJ7rXdpb4XamhowODBg6FQKLBnzx67c+5WX3d5z2pNfV3pPeuFF16AQqGweyQlJUnn6+vrMXv2bHTp0gX+/v649957L6lve9VlyZIl6N69O7RaLVJSUrB9+3bHKyTIKW3fvl10795dDBw4UMyZM0c6fsstt4jhw4eL7OxscezYMfHPf/5TKJVKsWvXLqnMV199JYKDg8V7770ncnNzxcGDB8WqVauk801NTaJ///4iLS1N7N69W6xdu1aEhoaKefPmSWWOHz8ufH19xdy5c8WhQ4fE22+/LVQqlVi3bp3T1TcxMVHcdtttYu/eveLo0aPi4YcfFr6+vuLs2bNOW9/L1bm6ulr06NFD3H333WLfvn1i37594q677hLDhw8XZrNZunb8+PFi0KBBYtu2bWLLli0iISFBTJ48WTpvMBhERESEmDJlijhw4ID44osvhI+Pj3j//felMr/99ptQqVRi0aJF4tChQ+K5554T3t7eYv/+/U5V3//+97/iwQcfFD/++KM4duyY+Pbbb0V4eLh44okn3LK+F3rsscfErbfeKgCI3bt3u2193eU9q7X1daX3rAULFoh+/fqJs2fPSo/S0lLp/J/+9CcRGxsrNmzYIHbu3Cmuu+46MXLkSOl8e9Vl5cqVQq1Wi48//lgcPHhQzJgxQ+h0OlFcXOxQfZj0OKGqqiqRmJgoMjMzxZgxY+ySAD8/P/HJJ5/YlQ8JCREffvihEEKIxsZGERMTIz766KPLvv7atWuFUqkUer1eOvbee++JwMBA0dDQIIQQ4qmnnhL9+vWzu27SpEkiPT39Wqt3iWupb2lpqQAgNm/eLJ03Go0CgMjMzBRCOF99hbh8nX/88UehVCqFwWCQylZWVgqFQiHV59ChQwKA2LFjh1Tmv//9r1AoFOL06dNCCCHeffddERwcLNVPCCGefvpp0bt3b+n5fffdJyZMmGAXV0pKipg1a5ZT1bclixYtEvHx8dJzd6zv2rVrRVJSkjh48OAlSY871ded3rNaU19Xe89asGCBGDRoUIvnKisrhbe3t1i9erV07PDhwwKAyMrKate6jBgxQsyePVt6bjabRXR0tFi4cKFD9WH3lhOaPXs2JkyYgLS0tEvOjRw5EqtWrUJ5eTksFgtWrlyJ+vp6jB07FgCwa9cunD59GkqlEkOGDEFUVBRuvfVWHDhwQHqNrKwsDBgwABEREdKx9PR0GI1GHDx4UCpz8f3T09ORlZXlVPXt0qULevfujU8++QQ1NTVoamrC+++/j/DwcCQnJztlfa9U54aGBigUCrvFu7RaLZRKJX799VcpVp1Oh2HDhkll0tLSoFQqkZ2dLZW54YYboFar7eqTm5uLiooKqYzcv+PW1LclBoMBISEh0nN3q29xcTFmzJiBTz/9FL6+vpe8vjvV153es1pTX1d8z8rLy0N0dDR69OiBKVOmoLCwEACQk5ODxsZGuziSkpIQFxcnxdEedTGZTMjJybEro1QqkZaW5nB9mfQ4mZUrV2LXrl1YuHBhi+e//PJLNDY2okuXLtBoNJg1axa+/vprJCQkAACOHz8OwNoP+9xzz2HNmjUIDg7G2LFjUV5eDsA6TubCf4AApOd6vf6KZYxGI+rq6pymvgqFAuvXr8fu3bsREBAArVaL119/HevWrUNwcLDT1fdqdb7uuuvg5+eHp59+GrW1taipqcGTTz4Js9mMs2fPSrGGh4fbXefl5YWQkJCr1sd27kplbOfby7XW92L5+fl4++23MWvWLOmYO9VXCIEHH3wQf/rTn+wS2wu5U33d6T2rNfV1tfeslJQULF++HOvWrcN7772HgoICjB49GlVVVdDr9VCr1dDpdJfE0R7vRba6lJWVwWw2t8u/ZyY9TqSoqAhz5szBZ599Bq1W22KZ559/HpWVlVi/fj127tyJuXPn4r777sP+/fsBABaLBQDwt7/9Dffeey+Sk5OxbNkyKBQKrF69utPq0hrtUV8hBGbPno3w8HBs2bIF27dvR0ZGBu64447LfmjK6Wp1DgsLw+rVq/H999/D398fQUFBqKysxNChQ6FUut5/1/au7+nTpzF+/HhMnDgRM2bM6IwqOKQ96vv222+jqqoK8+bN6+zwHdYe9XWn96zW1NfV3rNuvfVWTJw4EQMHDkR6ejrWrl2LyspKfPnll3KH1iZecgdA5+Xk5KCkpARDhw6VjpnNZmzevBnvvPMOcnNz8c477+DAgQPo168fAGDQoEHYsmULlixZgqVLlyIqKgoA0LdvX+k1NBoNevToITVJRkZGXjLq3TbaPjIyUvp68Qj84uJiBAYGwsfHx2nqu3HjRqxZswYVFRUIDAwEALz77rvIzMzEihUr8MwzzzhNfVtT54aGBowbNw7Hjh1DWVkZvLy8oNPpEBkZiR49ekixlpSU2L1uU1MTysvLr1qf1tTZdt5Z6mtz5swZ3HjjjRg5ciQ++OADu3PuVN+NGzciKyvrkv2Jhg0bhilTpmDFihVuVV93es9q7e/Xld6zLqbT6dCrVy/k5+fjlltugclkQmVlpV1rz4X/ztqjLiqVCiqVql3+Pbven45u7Oabb8b+/fuxZ88e6WF7o9uzZw9qa2sB4JK/gFUqlfTXUnJyMjQaDXJzc6XzjY2NOHHiBLp16wYASE1Nxf79++0+ODMzMxEYGCi98aSmpmLDhg1298nMzERqaqpT1fdyZZRKpVTGWeoLXL3OKpVKKhsaGgqdToeNGzeipKQEd955pxRrZWUlcnJypLIbN26ExWJBSkqKVGbz5s1obGy0q0/v3r2lJnRn+B23pr6AtYVn7NixUivAxb9vd6rvW2+9hb1790rX26aZr1q1Ci+99JLb1ded3rNaU19Xe8+6WHV1NY4dO4aoqCgkJyfD29vbLo7c3FwUFhZKcbRHXdRqNZKTk+3KWCwWbNiwwfH6OjTsmTrdhTMDTCaTSEhIEKNHjxbZ2dkiPz9fvPrqq0KhUIgffvhBumbOnDkiJiZG/Pjjj+LIkSNi+vTpIjw8XJSXlwshzk8hHDdunNizZ49Yt26dCAsLa3EK4V//+ldx+PBhsWTJkg6fwt2W+paWloouXbqIe+65R+zZs0fk5uaKJ598Unh7e4s9e/Y4fX0vrrMQQnz88cciKytL5Ofni08//VSEhISIuXPn2l0zfvx4MWTIEJGdnS1+/fVXkZiYaDdlvbKyUkRERIj7779fHDhwQKxcuVL4+vpeMqXZy8tLvPrqq+Lw4cNiwYIFHTqlua31PXXqlEhISBA333yzOHXqlN3UWXes78UKCgoumb3lbvV1l/es1tTX1d6znnjiCfHLL7+IgoIC8dtvv4m0tDQRGhoqSkpKhBDWKetxcXFi48aNYufOnSI1NVWkpqZK17dXXVauXCk0Go1Yvny5OHTokJg5c6bQ6XR2s8Jag0mPk7v4P9TRo0fFPffcI8LDw4Wvr68YOHDgJVO6TSaTeOKJJ0R4eLgICAgQaWlp4sCBA3ZlTpw4IW699Vbh4+MjQkNDxRNPPCEaGxvtyvz8889i8ODBQq1Wix49eohly5Z1VDUlbanvjh07xLhx40RISIgICAgQ1113nVi7dq1dGWetrxCX1vnpp58WERERwtvbWyQmJorXXntNWCwWu2vOnTsnJk+eLPz9/UVgYKCYNm2aqKqqsiuzd+9eMWrUKKHRaERMTIx4+eWXL7n3l19+KXr16iXUarXo16+fXfLcURyt77JlywSAFh8Xcpf6XqylpEcI96qvO71ntaa+rvSeNWnSJBEVFSXUarWIiYkRkyZNEvn5+dL5uro68fDDD4vg4GDh6+sr7r77brs/SNqzLm+//baIi4sTarVajBgxQmzbts3h+iiEEMKxtiEiIiIi18MxPUREROQRmPQQERGRR2DSQ0RERB6BSQ8RERF5BCY9RERE5BGY9BAREZFHYNJDREREHoFJDxEREXkEJj1ERETkEZj0EBERkUdg0kNEREQegUkPEREReYT/D90iTxTyruXTAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "birth_vect = np.linspace(48500, 50000, 100)\n",
    "\n",
    "plt.plot(birth_vect, norm.pdf(birth_vect, mu, np.sqrt(sigma_2)), label='H0')\n",
    "plt.axvline(mal_births, color='r', label='Observed')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's clear that the distribution of the data under $H_0$ cannot explain the observed number of male births, thus we reject $H_0$."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let be $H_1: \\theta\\neq 0.5$. From a Bayesian point of view, we need to calculate the posterior probability of $H_0$, which (using Bayes' theorem and total probability law) can be calculated as\n",
    "\n",
    "$$P(H_0|Y)=\\frac{P(Y|H_0)P(H_0)}{P(Y|H_0)P(H_0)+P(Y|H_1)P(H_1)}.$$\n",
    "\n",
    "We can consider prior probabilities $P(H_0)=P(H_1)=0.5$.\n",
    "\n",
    "Furthermore, $P(Y|H_0)$ and $P(Y|H_1)$ can be calculated as\n",
    "\n",
    "$\\begin{aligned}\n",
    "P(Y|H_0)&=\\binom nY \\theta_0^Y (1-\\theta_0)^{n-Y} \\\\\n",
    "P(Y|H_1)&=\\int\\limits_0^1 \\binom nY \\theta^Y (1-\\theta)^{n-Y}d\\theta=\\binom nY B(k+1, n-k+1)=1/(n+1),\n",
    "\\end{aligned}$\n",
    "\n",
    "where $\\theta_0=0.5$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9505229569718407"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_y_H0 = binom.pmf(mal_births, births, theta0)\n",
    "p_y_H1 = 1/(births+1)\n",
    "\n",
    "p_H0_y = p_y_H0 / (p_y_H0 + p_y_H1)\n",
    "\n",
    "p_H0_y"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus, from a Bayesian perspective **we do not reject the hypothesis $H_0$.** Furthermore, the posterior probability of $H_0$ is huge!"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What has happened!?"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next explanation was taken from Wikipedia."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The apparent disagreement between the two approaches is caused by a combination of factors. First, **the frequentist approach above tests $H _0$ without reference to $H_1$.** The Bayesian approach evaluates $H_0$ as an **alternative** to $H_1$, and finds the first to be in better agreement with the observations. This is because the latter hypothesis is much more diffuse, as $\\theta$ can be anywhere in $[0,1]$ (except 0.5 [which has measured zero]), which results in it having a very low posterior probability. To understand why, it is helpful to consider the two hypotheses as generators of the observations:\n",
    "\n",
    "- Under $H_0$, we choose $\\theta=0.5$, and ask how likely it is to see 49,581 boys in 98,451 births.\n",
    "- Under $H_1$, we choose $\\theta$ randomly from anywhere within 0 to 1, and ask the same question.\n",
    "\n",
    "Most of the possible values for $\\theta$ under $H_1$ are very poorly supported by the observations. In essence, the apparent disagreement between the methods is not a disagreement at all, but rather two different statements about how the hypotheses relate to the data:\n",
    "\n",
    "- The frequentist finds that $H_0$ is a poor explanation for the observation.\n",
    "- The Bayesian finds that $H_0$ is a far better explanation for the observation than $H_1$.\n",
    "\n",
    "The ratio of the sex of newborns is improbably 50/50 male/female, according to the frequentist test. Yet 50/50 is a better approximation than most, but not all, other ratios. The hypothesis $\\theta=0.504$ would have fit the observation much better than almost all other ratios, including $\\theta=0.5$.\n",
    "\n",
    "Thus, if we replace $H_1$  with $H_2:\\theta=\\hat\\theta$, i.e., the maximum likelihood estimate for $\\theta$, the posterior probability of $H_0$ would be only 0.07 compared to 0.93 for $H_2$ (of course, one cannot actually use the MLE as part of a prior distribution)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.07126646849013898"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_y_H0 = binom.pmf(mal_births, births, theta0)\n",
    "p_y_H2 = binom.pmf(mal_births, births, theta_hat)\n",
    "\n",
    "p_H0_y = p_y_H0 / (p_y_H0+ p_y_H2)\n",
    "\n",
    "p_H0_y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.928733531509861"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_H2_y = p_y_H2 / (p_y_H0+ p_y_H2)\n",
    "\n",
    "p_H2_y"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Another example"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose we observe $y=1$ from the model $y\\sim\\textsf{N}(\\theta,1)$ with a uniform prior density on $\\theta$, thus $\\theta|y\\sim\\textsf{N}(y,1)$. And consider the hypothesis $H_0:\\theta=0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.31731050786291415\n",
      "0.8413447460685429\n"
     ]
    }
   ],
   "source": [
    "# p-value for a double side hypothesis test\n",
    "print(2*norm.sf(1, 0, 1))\n",
    "\n",
    "# posterior probability of the hypothesis\n",
    "print(norm.sf(0, 1, 1))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order for a Bayesian analysis to yield a nonzero probability for a point null hypothesis, it must begin with a nonzero prior probability for that hypothesis; in the case of a continuous parameter, such a prior distribution (comprising a discrete mass, of say 0.5, at $\\theta_0$ mixed with a continuous density elsewhere) usually seems contrived. In fact, most of the difficulties in interpreting hypothesis tests arise from the artificial dichotomy that is required between $\\theta=\\theta_0$ and $\\theta\\neq\\theta_0$. Difficulties related to this dichotomy are widely acknowledged from **all perspectives** on statistical inference. In problems involving a continuous parameter $\\theta$ (say the difference between two means), the hypothesis that $\\theta$ is exactly zero is rarely reasonable, and it is of more interest to estimate a posterior distribution or a corresponding interval estimate of $\\theta$ or the posterior probability that $\\theta>\\theta_0$, which is a more satisfactory and informative conclusion than the dichotomous verdict ‘reject’ or ‘do not reject.’"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "e7370f93d1d0cde622a1f8e1c04877d8463912d04d973331ad4851f04de6915a"
  },
  "kernelspec": {
   "display_name": "Python 3.10.2 64-bit",
   "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.11.8"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}