CS146 Session 3.1 PCW

CS146 Session 3.1 PCW.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Normal likelihoods and normal-inverse-gamma priors\n",
    "\n",
    "Today we explore how samples from a prior distribution can be interpreted as instances of the likelihood function. Specifically, we look at how samples from a normal-inverse-gamma (NIG) distribution can be interpreted as normal distributions.\n",
    "\n",
    "**In short:** Each sample from the NIG distribution is a pair $(x, \\sigma^2)$. These values specify the mean and variance of a normal distribution and so we can think of the sample (the pair of values) as an instance of the normal distribution (which will be our likelihood function). More below.\n",
    "\n",
    "## Normal-inverse-gamma in SciPy\n",
    "\n",
    "Even though SciPy does have classes defined for the normal distribution (`scipy.stats.norm`) and the inverse-gamma distribution (`scipy.stats.invgamma`), it does not have one defined for the normal-inverse-gamma distribution. To help you, the functions below implement the probability density function and a sampler for the normal-inverse-gamma distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as sts\n",
    "import matplotlib.pyplot as plt\n",
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Function definitions for the normal-inverse-gamma distribution. The parameters\n",
    "of the distribution, namely mu (μ), either lambda (λ) or nu (ν), alpha (α),\n",
    "beta (β), are used as defined here:\n",
    "\n",
    "  https://en.wikipedia.org/wiki/Normal-inverse-gamma_distribution\n",
    "\n",
    "Note that we use the symbol nu (ν) rather than lambda (λ) for the third\n",
    "parameter. This is to match the notation used in the conjugate priors table on\n",
    "Wikipedia:\n",
    "\n",
    "  https://en.wikipedia.org/wiki/Conjugate_prior#Table_of_conjugate_distributions\n",
    "'''\n",
    "\n",
    "def norminvgamma_pdf(x, sigma2, mu, nu, alpha, beta):\n",
    "    '''\n",
    "    The probability density function of the normal-inverse-gamma distribution at\n",
    "    x (mean) and sigma2 (variance).\n",
    "    '''\n",
    "    return (\n",
    "        sts.norm.pdf(x, loc=mu, scale=np.sqrt(sigma2 / nu)) *\n",
    "        sts.invgamma.pdf(sigma2, a=alpha, scale=beta))\n",
    "\n",
    "def norminvgamma_rvs(mu, nu, alpha, beta, size=1):\n",
    "    '''\n",
    "    Generate n samples from the normal-inverse-gamma distribution. This function\n",
    "    returns a (size x 2) matrix where each row contains a sample, (x, sigma2).\n",
    "    '''\n",
    "    # Sample sigma^2 from the inverse-gamma distribution\n",
    "    sigma2 = sts.invgamma.rvs(a=alpha, scale=beta, size=size)\n",
    "    # Sample x from the normal distribution\n",
    "    x = sts.norm.rvs(loc=mu, scale=np.sqrt(sigma2 / nu), size=size)\n",
    "    return np.vstack((x, sigma2)).transpose()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 1\n",
    "\n",
    "1. Generate 10 samples from the normal-inverse-gamma (NIG) distribution with parameters as\n",
    "   provided below.\n",
    "   \n",
    "   Each sample corresponds to the mean and variance of a normal\n",
    "   distribution.\n",
    "   \n",
    "   With these NIG parameters, the prior 95% confidence interval for\n",
    "   the mean is about [-10, 10] and for the variance [0.1, 10].\n",
    "   \n",
    "   In practice you would\n",
    "   work the other way around: use confidence intervals (or other information) to determine values for the\n",
    "   prior hyperparameters.\n",
    "\n",
    "\n",
    "2. Plot the 10 normal distributions corresponding to your 10 samples. To see the functions\n",
    "   clearly, plot your graphs on the domain [-15, 15].\n",
    "   \n",
    "   You should see that the 10 samples\n",
    "   (normal distributions) are all quite different. This means the prior is quite broad\n",
    "   (uncertain) over the mean and variance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Normal-inverse-gamma prior hyperparameters\n",
    "mu_0 = 0        # The prior mean is centered around 0.\n",
    "nu_0 = 0.054    # The smaller ν₀ is, the more uncertain we are about the prior mean.\n",
    "alpha_0 = 1.12  # α₀ and β₀ control the marginal prior over the variance.\n",
    "beta_0 = 0.4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOy9eZAc2X3f+XmVdZ99N4DuBhoNNIABMABmBnNpSJNDjnhYK1LrK8Sww6FYrxXhDUnrsHZX2tgNhVZrx27Yf9hmrLy2vKZsaS3L8sqUSIni0BxqKA45BwAOjgFmcHYDfaDP6uq6r8y3f2Rl9VFX1tFdPY33iUAAnZmV9QB0f/NX3/c7hJQShUKhUHz8cXR7AQqFQqHoDErQFQqFYp+gBF2hUCj2CUrQFQqFYp+gBF2hUCj2Cc5uvfHAwIAcHx/v1tsrFArFx5IrV66sSCkHq53rmqCPj49z+fLlbr29QqFQfCwRQjysdU5ZLgqFQrFPUIKuUCgU+wQl6AqFQrFPUIKuUCgU+wQl6AqFQrFPUIKuUCgU+wQl6AqFQrFPUIKuaArd0PnDO3/IWnat20tRKBTbUIKuaIqvvv9Vfv3tX+eP7/1xt5eiUCi2oQRdYZvXp1/nax98DYA7a3e6vBqFQrGdrpX+Kz5+/MmDP2EkOMJoaJS7sbvdXo5CodiGitAVtplPznO85zin+05zP3afolHs9pIUCsUmlKArbCGlZD45z0hwhMneSQpGgYfxmj2CFApFF1CCrrBFPB8nWUhyKHiIyd5JAO6uKdtFodhLKEFX2GIuOQfAaHCUicgEmtDUxqhCscdQgq6wxXxyHoBDwUO4NTfj4XEVoSsUewwl6ApbWBH6oeAhACZ7J1Wmi0Kxx1CCrrDFXHKOkCtExBMBYDQ0ymJqESlll1emUCgslKArbDGXnCtH5wB93j6Kskg8H+/iqhQKxWaUoCtsYaUsWvR5+wCIZqPdWpJCodiGEnRFQ6SUVSN0UIKuUOwllKArGrKWWyNTzKgIXaHY4yhBVzRkKb0EwIHAgfKxfl8/ANGMEnSFYq+gBF3RkFguBkCPp6d8zPqzitAVir2DEnRFQ2LZSkF3Opz0eHpYza52a1kKhWIbStAVDSlH6N6eLcf7vH0qQlco9hBK0BUNsQQ94o5sOa4EXaHYWyhBVzRkPbdOwBXApbm2HFeCrlDsLZSgKxoSy8W2+OcWStAVir2FEnRFQ2oKuq+P9dw6BaPQhVUpFIrtKEFXNCSWrS7o/d7+8nmFQtF9bAm6EOILQojbQoh7QohfrXL+sBDiz4UQ7wshrgsh/nLnl6roFrFcrNxlcTOqWlSh2Fs0FHQhhAb8JvBF4DTwFSHE6W2X/a/AH0gpnwF+FvgXnV6oonus59ZreuiAykVXKPYIdiL0F4B7UsoHUso88PvAl7ddI4Fw6c8RYL5zS1R0k6JRJFFI1BV0FaErFHsDO4I+Asxs+nq2dGwzvw78LSHELPAt4Ber3UgI8fNCiMtCiMvLy8stLFex26zn1gGqWi693l5A9XNRKPYKdgRdVDm2fUzNV4B/K6UcBf4y8LtCiIp7Syl/S0p5UUp5cXBwsPnVKnYdS9At8d5MyB1CINSQC4Vij2BH0GeBsU1fj1Jpqfwd4A8ApJRvA15goBMLVHSXcpVolQjdIRwE3UES+cRuL0uhUFTBjqBfAiaFEEeFEG7MTc9vbLvmEfBZACHEU5iCrjyVfcBabg2gqocOEHaHlaA/YdxOZbm8nur2MhRVaCjoUsoi8AvA68CHmNksN4UQvyGE+FLpsl8G/q4Q4hrwH4Cfk2p68L7AslzqCbqyXJ4c/mhxjc9fvs3PXrtP3jC6vRzFNpx2LpJSfgtzs3PzsV/b9OdbwCudXZpiL1CtF/pmQu6QitCfEKYzOf7erYcc8riYyxV4N5bik32hbi9LsQlVKaqoSywXw+Vw4XP6qp4PuUMqQn9CuLSeQgL/+uw4Hofgv6yq//e9hhJ0RV3Wc+v0enoRolqyk4rQnySuJdL4HA7Oh/z8RE+Q7ypB33MoQVfUJZaNEfaEa55Xgv7kcC2e4VzIhyYEP9kf5kEmx/10ttvLUmxCCbqiLslCkrC7vqCni2mKRnEXV6XYbYqG5INkmvMhPwCf7Te/J36wluzmshTbUIKuqEsinyDoDtY8b4l9Mq9+sPczd9NZMobkXMjcSznsdeNzCKbSuS6vTLEZJeiKuqQKKQKuQM3zIbeZ5aBsl/3N1UQaoByhCyE44vMwnVWCvpdQgq6oS7KQJOSqnZpmnYsX1AbZfuZ6IkNAc3DM7ykfG/e5mUrnu7gqxXaUoCvq0shysSL0eE4J+n7mbirLqYAXx6Zsp3Gfh0fZHIaqIdwzKEFX1CSv5ykYBYKuxoKuLJf9zVwuz6jXveXYuM9D1pAs5NQIwr2CrUpRxZOJJdJ2NkWVoLdHfjbB+renMdIFQn9pFP+FoW4vqYwhJfO5Al8c2CroR32m/TKdyXNom9gruoOK0BU1SRXMBkz1InQrR10Jensk335M/mEcPZ4n8cO9NR9mJV8kZ0hGvK4tx8d9pohPZ9TG6F5BCbqiJolCKUKvI+h+px+HcKjy/zaQhiR7O4r3dD/Blw9RmEmgJ/bOZuNszlzL2LYofMTjximUoO8llKAramLlltezXIQQqlq0TQrzSYxkAe/JXrxPmWP9sh/tnSlQc1nTIx/ZJuhOh2DM62Yqs3cePk86StAVNUkWSoJeJ0IHM3XRiuYVzZO9vQYCvCd6cR0MoEU8ZD7cS4JuCvaIx1Vxbtzn4aGK0PcMStAVNbEToYPq59Iu2dtRXKMhtKAbIQTep/rI3V1DFvZGv/G5XJ6A5iDi1CrOHfF5eJhVEfpeQQm6oiZ2I3Q1tah1ZMEgP5PAe3yj37znWARZMCgspbu4sg3msgVGPO6qHTcPul2sF3Uy+t54+DzpKEFX1KQcoTeyXNwhVVjUIsXVDEhwHfCXj7mGzD8Xl/eGoM/m8ox6K+0WgCGPmfm8lFe56HsBJeiKmiQLSTyaB5dW/YfZQlkurVMoibZzYEPQnf0+EFBYznRrWVuYzVYWFVkMu83vjUVVXLQnUIKuqEmykGwYnYOaWtQOxSVTtJ2DGxOhhNOBs8+7JyL0tG4QLeiMeGoIemmjdDGv2ifvBZSgK2qSzCfLpf31CLqDZPUsBUNFac1SXMmgRTw43Fs3HJ2D/rLYd5P5Ug769qIiiyG3abksKstlT6AEXVGTZCFZt3WuhdVxMV3ofkT5caOwnN4SnVs4h3wUVjJIo7uNr6w+LQeqpCwC9LucOAUsKctlT6AEXVGTZD7ZMGURKIu+lRWjsIeUkuJypqqguwb9UDTQ17o74i1a0AFTuKvhEIJBt0tZLnsEJeiKmtj10C3RV1OLmsNIFJA53RTvbVgi3+2N0dWCKdS1BB1M20VluewNlKAramJX0FWE3hrlDJdqlsvg3khdXC1F3r11BH3Y7VJZLnsEJeiKmtjdFLU8dKs7o8IexeXKDBcLLeDCEXCWr+kWq4UiPU4Nl6OyqMhi2KMsl72CEnRFVQxpNJwnahFwlyJ0Zbk0RTGaAadAC3uqntd6vRS77qEX69otYFouq4UihS5v4CqUoCtqkC6kkUh7aYslW0ZZLs2hx3I4Ix5EjejXGfGgr3e38dVqvkhfA0G3iouWlY/edZSgK6piibOdCF0Jemvo63m0SPXoHEDr8aDH8sguzuxcLRTpd1c25dqMKi7aOyhBV1TFbqdFAJ/Th0M4lOXSJPp6rr6gRzzIvI7M6ru4qq3Ys1xMQVeZLt1HCbqiKnY7LYI55CLgCqhN0SaQhkSP59B66gs60DXbRUpJtKDbsFxK1aIq06XrKEFXVKUZQbeuU5aLfYxEHgzQIrWHK1tiX4x1R9DjRZ2ClA0j9MGyh64sl26jBF1RFcs+sbMpCqbXriwX+xRLUXcjywW6F6GvWlWi7vqC7nIIwk4H0YIS9G6jBF1RFWuk3OZN0aWl15mZ+bdVrw+5Q8pyaQLdjqCH3CDMbJhuELVRJWrR53IqQd8DNP6fUjyRpPKmOIfcIaSUTE19lanprwLgcHgYGfnKlusDrgBr2bVdX+fHFT1mdjGsJ+hCE2hhdxcjdFOgG3noYIr+qhL0rmMrQhdCfEEIcVsIcU8I8as1rvkbQohbQoibQojf6+wyFbtNopBAIPA5fayvX2Fq+qscPPBX6Ov7JLfv/G8kEje3XB90BVWE3gT6eg7hcuDw1xdLLeLpWoRulf03slzAitC7l42jMGko6EIIDfhN4IvAaeArQojT266ZBP5n4BUp5Rng7+/AWhW7SKqQIugK4hAOFha/gcPh5cSJX+fsmX+GEE7m5v9gy/UBV0BNLWoCK2Wx2pzOzWg9HvR4d4Yw22nMZaEsl72BnQj9BeCelPKBlDIP/D7w5W3X/F3gN6WUawBSyqXOLlOx2yTyCYLuIIZRZGnpzxgY+AxOZwCXq4f+vk+wsvLdLQUvezlCX5yK8xf/8Q7f+Tc3Sa51t/LSwhT02hkuFlrEQzGW60px0WqhiM8h8GuNZaLPpREtFLtaBKWwJ+gjwMymr2dLxzZzAjghhPihEOIdIcQXqt1ICPHzQojLQojLy8vLra1YsSsk8+Zwi7W1H1EoRDkw/NPlcwMDr5HLLZBIbtgue3VqUSGn86f/4hofvjXP/StLvP31e91eEtC4qMhCi3igaGCkdz/6XS00Lvu36Hc5yRqStG7s8KoU9bAj6NU+E25/DDuBSeDTwFeA/0cI0VPxIil/S0p5UUp5cXBwsNm1KnaRVCFFyB1icelbaFqQvr5Plc8NDLwKCFaW3ygfs/LV99rUops/mCOTKPCl//4CFz53mDvvLbI43d35p1KX6PF83aIii26mLq7mi7b8c4C+0nVqY7S72BH0WWBs09ejwHyVa/5YSlmQUk4BtzEFXvExJVFIEHAFiMUu0dv7Epq2IT5udz+RyLMsr3y3fMxKb9xLPnoxr/P+dx4xcrKXg8d7eO7zR/CFXLz7x/e7ui4jVQBZSktsgBY2rzESu++jRwu6Lf8cNnx2tTHaXewI+iVgUghxVAjhBn4W+Ma2a/4IeBVACDGAacE86ORCFbtLqpCiz+khk5kmEr5Qcb6v75Mkkx9SLG7t+bKXfPTpG6uk43me+8IRANw+J0+9cojZ2zFyme5FknqylLJoR9CDZhWmnth9K2vVRh8Xi76yoKsIvZs0FHQpZRH4BeB14EPgD6SUN4UQvyGE+FLpsteBVSHELeDPgf9RSrm6U4tW7DyJfIIhpzlcIRw+V3E+HDoLSBLJD4G92XFx5qMobq/GyIkN9+/ImT6kIZn7qHs581a07bAh6NY1ehci9GYE3bpOWS7dxdb/lpTyW8C3th37tU1/lsA/KP1S7AOS+ST9IgGIqoIeCp0FIJH4gN6e58uCvpci9NmP1jh0ohfHpiyN4YkILq/Gw1urTDzTnX0cPWlG21b0XQ+HW0N4tF23XDK6QVo37HvoLrPFrorQu4sq/VdUkNfz5I08YbmK338Mp7Oyn4vHM4jbPUQi/gGwMbVor3jo8ZUM8eUMo6d6txzXNAejJ3uZuRntWoqdFW07go0jdDCtmd2O0Jsp+wcIOzU0oTz0bqMEXVGBaZtIvMXHRKpE5xbh0Nly6uJei9Bnb5uWynZBBzh8pp9ENEtssTsZOUaygHA7cHjqD46wcHRB0DfK/m2uUQh6naq4qNsoQVdUkMqn6NUkmpEiHD5f87pQ6Cyp1H10Pb3nPPTZj9bwh930HaycuDT2lCnyc3diu70swNwUtRudA2ghF0ZydzdFy2X/NiN0MFsErKoWul1FCbqigkQhwSGXWSASCp2ueV0ofBYwSCRu7bmpRQv31zk02VO1tD484MPjd7Iy0x17yEjkbWW4WGgh966X/5ctF5seOmxUiyq6hxJ0RQXJfJIDLtNfDgRqlxOEyxujN8tTi/ZChJ5NFUhEswwert7LXQjBwFiQ5UfdEXQ9WcBhY0PUwhFyI/M6Rm73/OlmOi1a9KmOi11HCbqigmQhybDLwOHsq7ohauF2D+F0hkmlzUKdvdLPZXXWfKgMjNaetjQ4FmJ1LoXehVJ1I5m3leFiYUXzRnL3ovTVgo4mIOK056GDac+oTdHuogRdUUGykGTYKfH4jtS9TgiB33+MdKok6O7gnrBcVixBH6v9MBo8HEIvGsQWdndjVOoGRqrYtOUCu5uLHi31cXE06Aa5mX6Xk7VCEUM16OoaStAVFSRzCYZdBoHA8YbXBvxHSaengL0Toa/MJPCH3fjDtUXTEvvdtl2MlLm52cymaDeKi1bz9ouKLPpcTgwgVlRRerdQgq6oIJNbwOuASPBkw2v9/mPk8osUi4k946EvzyYZGKs/3Lpn2I/T7WB5lzdGrRJ+LdSM5WJea+xi+X8znRYtVHFR91GCrqhAz80CEA6danhtIDABQDo9RdAV7Lqg60WDtccpBkbrD7d2OAQDo6Fdj9CtPi5NReh+Fzh2OUJvouzfwsqIiarUxa6hBF1RgSgsABDwN7Zc/H5T0FPpB2aE3mUPPfo4haHLhhE6mJumq3Opna8Y1YuQM/9drCi7mU1R4RA4grtbXBQt2G+da9Gn+rl0HTUkWlGBs7hC1nDgdg80vNbnO4wQGunUfULuUNc99NU5Uzj7RxoLes8BP/lMkUyiUNdvb4uVu/Cffg6iU/D8f4Pu+DnAXmOuzWgh9671c9GlZK2g264StehTLXS7jorQFRV4jRjr0tdw3iWAw+HG6x0jnZ4i4Ap0fWpRbCGNwyGIDPkaXtt3wKwiXVvYoYdQdAp+61WIz8Pka/Cj/wvj5l8g3BoOd3NiqQVd5aZeO020UETSXJUoqBa6ewEl6IoK/CRJi8YRrkUgcIxU+v5GP5d896L02GKa8KAPzcYczJ4DfgDWdip18fv/GIwC/Pyfw9/4HXjpv0NfnEfzNZ/77gjuXoTezHDozfg1Bz6HQ1kuXUQJumILhpEjQJ6cI2L7NX7/UTKZaQJOMyru5sbo2mKanmG/rWuDPR6cHm1nIvTlO3D99+H5/xZ6x81jn/4VDMcQjtyjpm+nBV3oqcKudIiM5k3LZKBJDx3MTBcl6N1DCbpiC5nMHEKArvXZfo3PdwTDyBPSzMizWz66YUjWlzL02hR04RD0Dvt3prjorX8KTh+88vc3jnkj6L6jaPlHsPRhU7dzBN2gS+QuTFpqpezfot/lLD8QFLuPEnTFFjKZhwAIl/3hDz7fYQC80ozMuxWhJ1az6EWjbKXYoWfY33nLpZiDD78JZ/8KBLf+OxpFPw4Rgxv/qalbWrnou+Gjt2q5gPkQUB5691CCrthCOmPaAZrnoO3X+LzmDHG3vg7QtdRFq7+5XcsFoPeAn0Q0S6GTUeX9P4d8As78zJbDsmhgZAy0vl5T0JuwT6y89d1IXYy2E6G7laB3EyXoijIz0TTfu3mJrAGpfNj267zeQ4ADTY8C3YvQLS+8t4kIvbeU6dLRYRe3/hi8PXD0U1sO61bZ/5FTEHsEs5ds37JcLboLDbpW80UiTg2Xw34fFwvloXcXJegKAK48jPLl3/wha/EHrBYFv/OjJe4u2quidDhcpqgXloHueeixxTSegBNfE1WYveVMlw6tuZiH238KJ/8yaFuLh6wsFe34s+BwmbaMTTYi9N2xXJrNQbfoczlJ6gY5Y/e7WCqUoCuAeLbAz//OFcJeJ+cOplgpOnAJHz/325eIpuxFhD7fYYp5s8K0WxF6bDFte0PUIjJoZuasL2U6s4hHb0N2HZ766YpTlv/t6IvA2Asw9X3bt3X4nOBgVyYXtVL2b2HZNGuquKgrKEFX8C/fvM9qKs9Xf/YCFBdZLQr+/mfO8Xg9w7/+wQNb9/D5xshl57o6tSi2mKZnqDlBd7o1AhE38ZUOCfr0WyA0GP9ExalyhB50m3bM4+uQjtq6rXAIHAF3uRfMTtJK2b9Fvyou6ipK0J9wHq9n+DdvTfEzFw5xYjALssBKUXDu0AG++PRBfvfth6xnGkeFPu9hCoVVet2+rkToxbxOaj1vq0J0O+FBH+vLHRT0QxfAW7kHYUXoWsgFE58CJEz/wPatteDuzBZdzestbYjCpn4uqkFXV1CC/oTze+8+oqAb/PLnTpLJzACwUnQQcAX4e586RjJX5P9952HD+/h8ZqbLqMfbFQ89vpoFzHmhzRIZ8BFfyba/iHwa5i5Xjc7B3NAUHg3h0mDkOXAH4UETtkto5yN0KWV7lkuppYHaGO0OStCfYAxD8p9/PMcrxwcY6/OTzZqCHtUFIXeIsyMRPjk5wO++/RDDqJ9iZwn6sFvriuViWSatCHp40EcqlqPYru87ewn0PBypLuj65uHQmguOvNKUj64FXTveEz2pGxSkbFnQleXSXZSgP8G8OxVlLpbhrz03CkAmO48E1ormwGeAv/rsKAvxLFcerdW9l1VcNOjqTpZLW4Jeek3bUfr0WyAccPilqqeN7cOhxz8Bq/cguWzr9laEvpPl/+WiohY99F6naqHbTZSgP8H84Y9nCXqcfO70AQCy2VmKIoDOhqC/dnoYt9PBn15/XPdeTmcEpzNEr6aTKFRPd8w/esTKv/yXRP/dv0OPxTr6d4kvZ3F6NHxNTAKysDJd2t4YffgjOHi+qn8O2yJ0gNHnzd/nLtu6vRZ07Xj5v+V9t+qhOx2CHqemsly6hBL0J5SCbvD6zQW+cPYAvpLvmc3MkhNBAq4ADmF+awQ9Tl49Oci3bjxGr2O7CCHwekcIi0LVCH3t9/8j9z//BZb/+VdZ/D/+T+5+5rOk3nuvY3+f9ZUMkQGvrZa/29mI0NsQdEOH+fdh9IWal+jbI/SD582MmFm7gu4u32enaKfs30KV/3cPJehPKD9+uEYiW+S1p4bLx7LZeZJ4y21wLX7q3CGWEjkuT9dPsfN6RwmIbIWHnvz+91n4jd8g8IlPcPx7b3D0j76O68AB5n7xl8g/ar7zYDXiK5mW7BYAX8iF06O1l+myfBsKKXOzswqyaCAzxbIoA+D2w4GztiN062Gwk9WiG4LeWmERmNWiStC7gxL0J5Q/v72M0yF45Xg/AFLqZHPzJAw3IffWeZyfOTWE0yF48059r9frHcErU6Q2pS3q8Tjz/9Ov4Dl1ktF/9k9xHTyI99Qpxv7vfwHA3C//D217wlJK4sutC7oQgsiAtz0Pfe6K+fvIs1VPl4uKtltCIxdh7sdgo7LSsmt2NELPdypCV5ZLN1CC/oTy5u0lnh/vI+Q1BSaXW0TKImu6o+yfWwQ9Tp493Mtbd1fq3tPnHUWjgENuTC1a/drX0NfXOfQP/yGOwMZ93UeOMPQrv0L2xg0S3/1uW3+XdDxPsWC0LOhg2i5tWS5zV8ATgb5jVU9bUbW2vS3B6EXIxWHlTsO3KEfoO9iga7VQxOsQ+G0MCKmFsly6hxL0J5D5WIaPFhK8emqjtWs2Ow/ASkESdFdOK/rE5AAfzK/XbQXg9Y4A0KtJkvkkxeVlov/udwj/1E/hPX264vrIl34a99GjLP/zf47UW4/orMja2txshfCgj/hypvVPC/M/NguKHNV/pKwuidr2WaJNbIw6/C5w7GyEHi3o9LucLe1FWFiCvhvDOBRbUYL+BPIXJevk0yeHyscy2VkAlgt6hYcO8MnJAaSEH96rHaV7faag9zkNkvkk0X//75G5HIO/9ItVrxdOJ4O/9Ivk790n8cYbLf99NlIWvS3fIzLgo1gwSMdbiH4LGVi8WdM/h40eLFs2RcGM6D0R03ZpgFn+79rRFrrtFBVlMhl+9KMfMX3jGllDcvfRTIdXp2iEEvQnkHenogwEPUwObQh3NmMK+nw2W1XQz432EPY669ouPq+Zz96nSeKZNdb/89cJfvKTuI8cqfma0Oc+h3N4mNgf/mGrfx1T0AWE+lsX9LCVutjKxujCDTCKNf1z2BShb7dcHA4YPgOLH9h6Ky3o3tHy/9V8saWUxUQiwde+9jW+853v4EiZeyj/6vf+A1euXOn0EhV1sCXoQogvCCFuCyHuCSF+tc51f00IIYUQFzu3REUnkVLy7oNVXjzat+VjdTY7h9s9wHohU7EpCqA5BD9xbIC36kToTmcEHF76nJLMWz+iuLREz1//a3XXIzSNyH/9M6R+8BaFhYWW/k7x5Yw5H7SNzIxIyX9fb8VHf3zN/P3QMzUvMZIFhFdDuKr8yB04a0b4NjZGHUHXjpb/r7bQmCuXy/Hbv/3bxGIx/vbf/tv89S98DoDBY8f55je/ybVr13ZiqYoqNBR0IYQG/CbwReA08BUhRIUhKoQIAb8EvNvpRSo6x+xahvn1LC8c3TozNJudw+MZIatnKzZFLV6c6GMulmE+Vl30hBA43cP0aRLHN99AGxwg+KlPVb12Mz1/9a+CYbD+9a83/xfCFOF2NkShFN2LNiJ0Xy+ER2peUlFUtJkDT0M+CbHphm+10xF6tAXL5Y033iAajfI3/+bfZGJiohzhX3z1Mxw5coRvfetbJBL2eusr2sNOhP4CcE9K+UBKmQd+H/hylev+d+AfAx3ocqTYKd6bMnPJX5zYKuiZ7CxOt+mpV4vQAZ4fN19z+WHtNgBe7wgDDgPvpVv0fPnLCFfjyk332Bj+F19k/Y/+uKWNNDNlsXW7BUBzOgj2elpLXVy4YYpynY1EPZmv9M8ths9u3KcBjpBrx8r/c4ZBUjeaGm4xMzPDe++9xwsvvMD4+DhA+fXruuSnf/qnKRaLfPvb3+74ehWV2BH0EWDz7sZs6VgZIcQzwJiU8k/q3UgI8fNCiMtCiMvLy/b6Vyg6y7tTq0R8Lk4MbYi2lAbZ7GNwDQDUjNBPHQgRcGt1C4z8vsP0axKhG4Q+/wXb6wp/8QvkHz4kf++e7dfARtvcdiN0MG2XpouL9CIs3YID5+peZiQKtSP0oafMHjALjX10LeiGokTmOp/nXc5Bb8Jy+e53v0swGOSzn/1s+VhvuUGXzsDAAJ/85Ce5efMms7OznV2wogI7gl4t7CiHB0IIB/BPgV9udCMp5W9JKS9KKS8ODtqfKq/oHO9NRXl+vA/HpnmR+fwyUuYxtF4AQq7qEbpTczXeFIIAACAASURBVPDM4V4uTdeO0EP+I7hdkB7x4z17xva6gp/5DAhB/L/8F9uvgfba5m4nPNhCLvrqPShmzQi9DnoyX7khauHyQf+krY1RR2jnhkVHmyz7n5mZ4eHDh7zyyit4PJ7y8YhTw7Hpfi+//DJer5cf/vCHHV+zYit2BH0WGNv09Sgwv+nrEHAWeFMIMQ28BHxDbYzuPVaTOaZX01wc791y3EpZzDtMIQ+4q0foABfHe/loIU48W93H9QrzQb3w8lBTucyuoSF8Fy40XWRkCXA7OegW4QEf6XieQr6J6NeySeoIuiwYyKxeWSW6mQNP27JctHL5f+d99NVSdaddQX/rrbfw+Xw8++zW7B5NCHo2DYv2eDw8//zzfPjhh6ys1C9OU7SHHUG/BEwKIY4KIdzAzwLfsE5KKdellANSynEp5TjwDvAlKaW9BhWKXePqjNnh8NnDWwU9m5kDICPM8W21InQwfXQpzV4w1ZAfLgGweKZ5Tzv02mvkbn1IfnbO9mvaaZu7HSvTpamN0YXroLlh4ETNS/RaVaKbOXAW1mcgU79N8caw6M5H6JYA20lbXF1d5fbt27zwwgtbonOL/m3Voi+++CKapvHOO+90bsGKChoKupSyCPwC8DrwIfAHUsqbQojfEEJ8aacXqOgc7z+KoTkET49EthzPZk0BTRqmWITd1du/AlwY68EhzHtVo/DWbQDig81HkKHXTB82+eabtl/TTtvc7YRbaaO7cMP0wLXa718uKqrloQMMlyL8xZt1325HI/QmPPSrV68ihOC556oXU23v5xIMBjl79iw3btwgn9/5uahPKrby0KWU35JSnpBSHpNS/qPSsV+TUn6jyrWfVtH53uT9mTWeOhgqt8u1yGRncbn6iBdzAIQ9tQU94HFyfCjI9dlKQZdSkn3zErIg0GS86fW5jxzBdfgwqbfesv2adtrmbqeci95MhL74wYYY12CjqKiB5QINbRdHwAWCHclFjxaKaAJ6nPWzXAzD4OrVqxw/fpxwuPr3SrV+LhcuXCCXy/HRRx91bM2KrahK0ScE3ZBcm1nnmbHeinPZ7Bxe7wjxvCnC1SpFN3NutIfrs+sVqXO5O3fQl1fIFzx4ZWtj6IKfeIXUe+9h2Izi2mmbux1PwInbq5U3WhuSXIbUMgxX9qnZjCW+dSP00DAEBhtmuljl/zvjoRfpdTpxNHg43r9/n0QiwYULF2peU62F7pEjR+jp6eHq1asdWa+iEiXoTwj3lpIkc0UujPVUnMtm5/B5R0nkEwRdQTRH/Qjt/GiE1VSeuW0FRqkfmBPs054IgRbLEQKf+AQynSbz48a9TaSUpqD3d0bQhRDNZbosleyRofqCbs0B1QINbKHhs7BoZ2PUvWMeuh3//Nq1a/h8Pk6ePFnzmmoNuhwOBxcuXODBgwesr693ZM2KrShBf0J4vzQT9JnDWwVdSlmK0A8Rz8Xr+ucW50sPhWszW38ok2/9EM+JE+S9/YQcrQmO/4UXwem0Zbuk43mKeaPsfXeC8IDP/qbo0ofm7w0EXU/mET4nwtngx+3AWfOeev3o2xHaoQg9X6TfXf9hXigUuHPnDk899RROZ23x73M5KUqIF7dmDD39tGktffjhh+0vWFGBEvQnhOtz64S8To4ObE1JzBdWMYwcXp8Zodfzzy1OHQjj1hxbfHQjmyVz5QqBT3wCnH34HZJ8ofkoTAsG8D/zDMm3GucsJ8o56O1ViW4m3O8lvpq1V4m5eBP8/RAcqnuZkcij2dm0PXAO9Dys3K17mRZ074iHbqfT4v3798nn85yu0g55MwNua1j0VkHv7+9neHiYW7dutbdYRVWUoD8h3Jhd59xopGLz0Oqy6POOEs/Ha5b9b8btdPDUwRDXNgl65tp1ZKGA/4XncbrNsXZryfstrdX/8kvkbt9uOEja2rzsRA66RXjAh263je7SLTM6b+A568lC/ZRFC6sFQIMCI0fQhZ4odLz8347lcuvWLXw+H0ePHq17nfVgWK0y6OL06dM8evSIeLz5jXNFfZSgPwHkijofLcQ5uy1dESBbKiqyNkXtWC5gbox+MBfHKA2OTl++BELgf/ZZPB6zM0SsRUEPPP88SEm6QetVy+tup23udmy30TUMWPrIbH3bACORr78hajEwCZrHzG2vg1n+b3S0/F+XkrXScItaFItFbt++zalTp9C0+taMFaGv5CutoaeeegpAZbvsAErQnwDuLCQp6JJzI5UbopmMJej2I3SAc6MRkrkiD1bMbJbMlSt4Tp1CC4cJ+MzC4mT6YUvr9Z47h3C7Sb93qe518ZUMgYi7rba52wmXHg4NM11i0+ZQ6Ab+OVgRug3LRXPB4AlYrG9HWE2+Ojm5KFbQkdTPQZ+amiKXy5UFuR4DpQfDSpUIfWhoiIGBAeWj7wBK0J8Ars+Z1sW50coIPZOdweXqw+kMmB66zQh988aoLBRIv38V/0Wz20PAd5C8AelMaxNrHB4PvvPnSV9qJOjZjm6IwqY2uo0yXSzRbRChG3kdmdPtRehg2i5L9QXdavJldNBHt6yRgToR+p07d3C5XA3tFth4MFjFSts5ceIEDx8+JJfLtbBaRS2UoD8B3JhdJ+JzMdpbKX7ZzCxe7wgFvUCmmLEt6McGg/jdGtdnY2Rv3UJmMvhLVYNhT5g1XZDPPW55zf7nnyf70UfodXzWTuagWzhdGsEeT2PLxcpwGTxV9zIrG8VWhA5mxJ94DOnaHS3LEXqicxF6o7J/KSV37txhYmICl42WyB6Hg5DmqBqhA0xOTmIYBlNTU60vWlGBEvQngBtz1TdEwYzQfb6xclGRXcvFaiFwdXad9GXT6/ZfLAm6K0y0KNALiy2v2f/CC2AYNX10vWCQjOU6LuhgRukNLZelm9A7Dp76RVhWvrj9CL1k4dRpAbAjEXqDsv+lpSXW19c5caJ2z5rtDLidrNSI0MfGxnC73dy9Wz+jR9EcStD3OdmCzu2FREX/FgApdbLZ+XKGC9Qv+9/O+bEePpyPk7x0CffRozgHzH7qQXeQqO5AFGtHmY3wXTiPcLlIX6reRSIRzYKESAdTFi0iAzaKixZv2fLPDTuNuTZjZbrUsV0cfqv8v3MR+koDy+XOnTuAGVnbpd9VW9CdTicTExPcvXt3R4Z1PKkoQd/n3F5IUDRkVUHP5ZaQsoDXN0Yib44Is2u5gOnJF4pFUpevlKNzKAl6UeAw0hSLqZbW7fB68Z47V9NHX+9gl8XthAd9JGM5ioUaWSTFnNkH3c6GqFUlard5WHAYfH11UxeFJnD4XRgdrBZdzhcQ1G6de/fuXQ4ePFizd0s1BtzOqmmLFpOTk8TjcdSwm86hBH2fc33OLO55uuqG6NYcdGhO0M+P9nAkvoBIJsobogAuh4ukNFuqWp0cW8H//EWyt26hJysfCpbHvSOC3u8FuVG4VMHybZB6wx4usBGhOxqV/VsIYW602sh06WSEvpwv0uvScDoqbblsNsvMzAzHjh1r6p4DLldNDx3g+PHjAMp26SBK0Pc5N2Zj9AXcjPRU2xA1s1B8vjHiueYFfbTXxwuJRwBbBB02hmW0J+jPg66Teb+yr0t8JYPmcuAP27QymsB6SNT00S07ZKhxDrqeLODwOxFaEz9qQ6fNTVfDqHmJFnJ31ENfyRcZdFd/6ExPTyOlbF7Q3WY/F6OGpRKJRBgaGlKC3kGUoO9zbsyZBUXVN0RnAYHXe6hsudjdFAWzmdWLiUdEg324RrZOvJfOfqBNQX/mGXA6q+ajx1eyhPu9iCoRZbs0LC5aumUOtehvLHC63aKizQyfMXPcY9M1L9F2IEKv5Z/fv38fl8vF2NhY1fO16Hc50SXEirULoCYnJ3n06BHZrJot3wmUoO9jsgWdO4sJzlXxz8GM0D2eYRwOT0ubolJKJh7f5VrvONltfrPHPYAuRdnWaQWH34/vzJmqPnp8NdPxHHQLf9iN5nLU3hhdvAUDJ+sOtbAw7BYVbcbKba9juziC7o566CuFAoM1Mlzu37/P+Ph43WZc1dioFq3vo6v0xc6hBH0fc+txHN2QVUv+wYzQvd5RABL5BB7Ng0erHCdWi8KjR3jja9zon+DW46354hFPD+uGq60IHcD/wvNkPvgAI7MhrlJK4sudz0G3EEKUm3RVZemWOaXIBi1F6FZue51MFy3kQhYMjA6V/y/ni1UFPRqNEo1Gm7ZbYFO1aB1BHxsbw+PxKNulQyhB38d8UNoQrVYhCpDJzODzmYLeTNm/RfqymVJ4o3+C6zNbG2n1eHpY00Xbgu577jkoFslc2+hvkksVyWf18oShnaBmX/R0FOJzZqtbGxjJvP2URQtPEHqP1s10sWaLdsJHz+gGSd2o6qE/ePAAoCVBL1eL1tkY1TSNiYkJ7t27p9IXO4AS9H3M9dl1BoJuDkYqc7UNI08ut4DPa/qizTTmskhfuozW10fu4Gg5m8Yi4omwlDfKvWJaxf/MMyAE6Ssb+egbKYudz0G3sPqiV4jMkr2SfwAjpyPzhv2Uxc00yHTROtjPZbnUQGugSoR+//59wuEwA6Uag2ao189lMxMTE8TjcVZXV5t+D8VWlKDvY27MrtfcEM1m5wGJd1OE3rSgX76M/7nnODfaw43ZrYLe4+lhtQiFwiq63sSMzm1okQieyUkymypG4zuYg24R7veSz+rk0tvEyKrgHG4coZdTFpuN0MHMdIneh0L1fzvLxtHttPltgGWJbN8U1XWdBw8ecOzYsZZmtva5andc3IwV/d+/31p3TsUGStD3KZm8zt2l2huiGznopQg915zlUlhYoDA7i//5izw9GuHecpJUbkP8Ip4IUV1sea9W8V98jvTVa8iief+daJu7nXLq4nbbZeFGaajFcMN76OulKtFWUiuHz4A0zJz3KmgRc69Dj7ff3Gq5FEFvt1zm5+fJ5XIt2S0AToegz6WxXMdDB+jr66Onp6ds7yhaRwn6PuXW43UMCU+PVrbMha056ACxXIxeb+UA6VpY/Vt8zz3HudEIUsLN+Y2NUTNCLwl65lFLfwcL33PPIdNpsh+a/bPjK1l8IRdub3NZF81gDc1Y3566uHjTFFsbEasltpb4NkU506V6TxeH3wlOUX5otIMluNs3Ra2IeWJiouV7D7pdDQUdzCh9amoKXe9cj/cnESXo+xTLAqlW8g9m1CyEC4/HjDRjuRg9nuriX4305Us4AgG8p06Vs2g2j6SLeCKsFM1vr3YF3eriaPnoO9FlcTtW9L8lQjd0s+DHht0CbUbofRPg9NbMdBFCoIU9nYnQLQ/dVSnohw4dwu/3t3zvYbeTxQaWC5gPjXw+z9xce5voTzpK0Pcp1+fWGQx5GA5Xjw4zmRm83oMIoZEpZsgUM01F6JkrV/A99yxC0xgKeTkY8XJj08Zoj6eHlAFSeNsWdNeBA7hGRshcMStGd0PQ3V4nvpBra4QenYJixtaGKJgRunBrOFr5JOHQYPBk/a6LYXdHIvSVfJGw04F3UzVrNptldna2ZbvFYsjtYjHXWNCPHj2KEEL56G2iBH2fcmN2nXM1NkTBHD1n+eexrBlZ93rsCXpxbY3c3Xv4Lz5fPvb0SITrs1sFHQRFradtQYeSj37lCsWiTiKa29EMF4ueIT/rS5sE3UojtBuhx/NokTZaEwyfrS/oEQ9Ghzz0gW09zq1y/3bsFoADHhdL+WLDlES/38+hQ4eUj94mStD3IalckfvLyZoFRVCK0EsZLtGc2ebWboRuZZxs7rB4bjTC1EqK9YwZjYXcIQSCjAh1RNB9zz2HHo0SvXofaUh6hlu3AewSGfIRW0pvHFi8CcLRcKiFhb6ea80/txg6DaklSK1UPa2F3RTX823nby/nK6tEWy33386w20VBSqK1OlduYmJigtnZWdUGoA2UoO9DbsyZG6IXxqp74sViikIhWhGh93n7bN0/fekywuPBe3YjUj1X2ny9WbJdNIdG2BMmKb1kMrNI2d5ml+WjL10yPeWeod0QdD/p9Tz5bGlTb/Em9E+Cy96nA30935p/btFg2IUW9pjDojONNx3rsZIvVuSgt1ruv50hj/n6JRs++rFjx5BSMj093dZ7PskoQd+HXJupPUMUNhpmlSP0rBmh290UTV+5gu/8eRzuDbGyNl+vb/PR13QnUubJ5VqfXgTgnphA6+0letsca7cbgm69R9lHX/zAtn8uDYmeaDNCt6ydWoIe6Uwu+kKuwPCmlMVYLEY0Gm3bbgHK97WzMTo6OorL5VI+ehsoQd+HXJ2JcbjPT3+wuphkt+Wgx3IlD92G5aInU2Rv3dpitwD0BtyM9fm2FBhFPBGWS8Fju7aLEALfc88SW0zh8TvxNtvwqgUiQ6XUxaUMZOMQe2hb0I1kAYwWM1wsgkPgHzDH3VWhnIu+3rqPnirqJHSDg56Nf09LUNvdEIVNgp5r/CnC6XQyPj6uBL0NlKDvQ67NxDhfw24BSGceApT7uKxl19CEZquwKPP++2AYFf3PAc6N9HB9biN1scfTw3wpw6EjG6PPXSRp+In07byYw0Yuemwpvank3+6GaBs56JsZPl2zBYD1sGgnQn9cipwPbBL0Bw8eEAqFGBwcbPm+Fs1YLmD66NFolFgs1vhiRQVK0PcZS/Es8+vZmv45QCb9EKczhMtl9ixfy63R4+nBIRp/O6SvXAanE9+FCxXnnh6NMBPNsJYyBSbijjCXzSCEk3RHBP1Z0r4hAiTbvpcd3F4n/oib9aX0pgwXmymLpai5rQgdzCEaSx+aOfDbsIZFtxOhL5QeuFaEbhgGDx48YGJioqVy/+0ENI2g5ii/TyNUG4D2UIK+z7ha8s8vjNXOcEmnp/D7jpZ/YNeya7YzXNKXL+M9cxpHlWITq82AlY8e8URYy63j9Y6QKX0qaAft2Aly3j5867tXfFJOXVy8CZ4IREZtvc6KmtuP0M+Yue/Ryn7hwukwR9G1E6HntkboCwsLZDKZjvjnFgc8LlseOsDg4CChUEilL7aIEvR9xrXZGE6H4MyhOoKemcLvP1r+2q6gG7kc2WvX8T9XabcAnNlWMdrj6SFdTOPzjZNOtz/AIBEzfVjXTP15m50kMuQjtpxpquQfSlWimrA/S7QWB542f1+4XvW0FvF0JkIved2dKPffzpDbzEW3gxCCiYkJHjx4gFFnBJ+iOkrQ9xlXZ2KcOhjC69Kqntf1HNnsPL7Ngl6yXBqRvX4dWShU9c8BIj4XRwcC5QIj6yHhcB8knZ5CyvZ+QGOLZk64+/5V9Hi8wdWdoWfITyaeJ/f4gW27BUo56CF3+yPyhp4ChwseX6t6ut1q0ce5AiHNQcBpfr/cv3+f4eFhQqHmeuPXY9jttFUtanHs2DEymQwLCwsdW8OTgi1BF0J8QQhxWwhxTwjxq1XO/wMhxC0hxHUhxBtCiCOdX6qiEYYhuT6zzvkaDbmAkvUh8fs2/oti2ZitHPT0lSsgBP7nnq15zbnRSNly6feaHn3R2Y9hZMnl2vsBtYp8fOllMlevtnUvu/QeMK2ltXQPHKrcN6hFMdZmyqKF02OKeg1Bd/Z6KcZaL8RZyBXKdks+n2dmZqaj0TnAkMfFUr5guwDKen/lozdPQ0EXQmjAbwJfBE4DXxFCnN522fvARSnlOeD/A/5xpxeqaMyDlRSJXLHuhmg6Y1ofluWiG7rtxlypd9/Fc+IEWqS2nfP0SITH61mW4ln6faagJwma792m7bK2kCYQduEUxapzRneC3gMB872Lo3DwvO3X6WtZnH0dak9w8Lwp6FUEUev1ILM6Rrq1QRePc4XyhujDhw/Rdb0j6YqbGXa7yBiShG7vE1owGGR4eFj56C1gJ0J/AbgnpXwgpcwDvw98efMFUso/l1JaNdLvAPZ2jhQdZWNDtI6gp6cB8PvHAXOwhUQ29NCNXI7Mj98n8NJLda979oh5nx8/WmPQb6a9rRlmNkYq3d4P6NrjFL2HgvjOnSP1zrtt3csu4QEvDofBmjFuu+Rf6oZpufR2IEIHU9AzUViv7Cvv7DUfGsW11nz0hfxGhP7gwQM0TePw4cOtr7UKh0r3n8/Zt4YmJiZ49OgR+XznBmE/CdgR9BFgZtPXs6Vjtfg7wJ9VOyGE+HkhxGUhxOXl5WX7q1TY4tpMjKDHybHBYM1r0ukp3O5BnE7TI13LrgGNG3Nl3n8fmcvhf7m+oJ85FMatOfjxoxgDPnNs2XIuj6YFSLch6NKQRBfS9B0MEHjpJbI3b6Kvrzd+YZs4NAe93lXWtFOg2dvg1NfzIDfEtm0OlqyeKraLVnoPfa1520WXkqV8gYMe84F7//59Dh8+jNvdZqrlNka85v1ms8356Lqu8+hR++muTxJ2BL3ark5VM0wI8beAi8A/qXZeSvlbUsqLUsqLnShaUGzl6kyMc6MRHHU24tLpbRkuuZKgN4jQU2+/A5pWc0PUwuPUeHo0wpWHa3g0DyF3iJXsKn5/e5kuibUsxZxO78EAgZdfAsMg9d57Ld/PNoZBLw+IFurFMFspRk1x1Tol6MNnzKZgVQTdWfoU0EqEvpwvokszrTCRSLC0tNRx/xxg1Gs+CGez9qPtw4cPo2ma8tGbxI6gzwKbW66NAvPbLxJCvAb8L8CXpJTt9/RUNEW2oPPh43hduwWsHPTx8tflCL2RoL/zNr6nn0YL1o7+LZ493MONuXXyRYMB3wDLmWX8/om2BH3tseno9R0M4Dt/HuHzkX77nZbvZ/+Np+gVU8TTfop5ew3GrGi5Yx662w8DJ+Fx5Uaw8DkRHq2lCP3xppTFTpb7b2fY7cIlRFOC7na7OXz4sBL0JrEj6JeASSHEUSGEG/hZ4BubLxBCPAP8K0wxX+r8MhWNuD67TtGQPHO4tjDn81EKhVUCgePlY0tp87/LskeqoScSZG980NBusXjuSC/5osHN+XUGfAOsZlbx+yfIZufQ9dYyMtYWUoAp6MLtxn/xIql3dkHQH1+l1zkDiK2tdOtQXMuCoL1e6NsZeRbmrlRsjAohzEyXFgR9oeRpH/C4uHPnDsFgkAMHDnRkuZtxCMEhj4u5JgQdYHJykqWlJdUGoAkaCrqUsgj8AvA68CHwB1LKm0KI3xBCfKl02T8BgsB/EkJcFUJ8o8btFDvEpWmzY+LFI7UFPZW6C0AgcKJ8bCm9hNPhrJu2mL50CQyDwEsv21rLs6WHypWHawz4BljJrJRsHkk6M23rHtuJPk7hC7nKTbkCL71E/sEDCovtdXFsyNyP6fOa+z3Wp4RG6GtmyqLQOljmMfIcpFdhbbrilNbraSlCnytF6ENOwb179zhx4gQOx86Upox63U156AAnTpjfp3fu3NmJJe1LbP3vSSm/JaU8IaU8JqX8R6Vjvyal/Ebpz69JKYellBdKv75U/46KTnNpOsrkUJDeQO2oMJkyfzACwcnysaX0EoO+wbp9XFJvv4PwevE9Yy8PeyjsZazPx+XptbLlEiw9RFLJ1n441x6n6DsYKH8d+ImXS2t7u6X72Wb2Ej2jQwhhPlTsUFzLds4/txgt7V3MXak4ZUbouaYHXcxk8vgcgvT8HPl8npMnT3ZipVUZ8bqYbSLLBWBgYID+/n5u3769Q6vaf6hK0X2AbkiuTK/x/NH6xUGp1F2czhAe93D52FJ6iWH/cJ1XQfqdt/E/++yW/ueNeGG8n0vTUQa8A2SKGXAdQAgnyVTzP5xSSqKP0/RuEnTPyZNovb0766MX8zB/Fe3wM0SG/KzO2WsKpkc7mINuMXQGnL6qgq71epE5velBFw+zOQ77PNy5cwen08nRo0cbv6hFRr1uFnIF8k2W8584cYLp6WlyObUtZwcl6PuAjxbiJHJFXhhvIOjJOwQCJ7Z00VtMLzLkH6r5muLysjk/1KZ/bvHiRB+rqTx6wdxEXc3F8PsnSCabF/TkWo58prglQhcOB/4XXyT1zjttj2CryeIN0HMw+jwDY0Fbgi6LBnoiX84+6Ria06xUnb1cccrZV8p0iTZnuzzM5DnidXPnzh0mJiY6nq64mVGvG8nGRqxdTp48ia7ranPUJkrQ9wGXpkr++Xht/1xKSTJ1l0BgcsuxRoJuFfDY9c8tXix9WlhYM8VmJbNCMHiSVKp5y2VlJgHA4OGt/UUCL79McXGR/NR00/e0hSWeo88zMBokvpIl16AisxjLgexgyuJmRp4zUxeLW60LrVxcZF/QpZQ8yuYZlEVisVjZr94pxjxWLnpztsvY2Bher1fZLjZRgr4PeG86yqGIl9He2mPZ8vllisVY2csGSBaSZIqZupZL6p23cYTDeE8/1dSaDvf5ORD28mDB/BYzffSTZLNzFIuJpu61PJMEAf0jW1MmA6VPDam3f9TU/WwzewlChyAywsCo+TBpFKUXV81xdc7+HRJ0PWd+ctiEs98cxFFcydi+1WpBJ6UbuGJmMLDTgj7aQnERgKZpTE5OcvfuXdV90QZK0D/mGIbkR/dXeflY7bRD2LQhGti6IQowHKgu6FJKUn/xAwIvv4zQqndvrIUQgheO9vHBI/OHcDWzSiB4Ysta7LIyk6BnyI/Ls3UNrrExXKOjpH7wVlP3s83Me+XNyIEx82GyPNNA0JdKgj64AzNPD5dsr0db9w0cHg0t4im/tx0eZUxPujg3w8GDBwmHwx1bZjUOtVBcZHHy5EnS6TSzs5WtDxRbUYL+MefW4zixdIFXjvfXvc7KLrFEFWAxZab81bJcsrduUVxeJvjqp1ta24sTfSyvaziFk+W0GaEDTfvoyzMJBscqC5qEEAQ/9SlS77yDkW2942BV4o/NGaJjLwLgD7vxhVyszjYQ9OU0joATrd0+6NUIH4Leo/Cw8hOJc9BHYdleWiXAo5Kw5mcf7Wh2i4XH4WDI7Ww60wXg+PHjOBwOZbvYQAn6x5wf3lsB4JXj9SP0ROImHvcwHvfGdYvp+oKefPNNEILgX/pLLa3tFKB3KgAAIABJREFUk8cHAYFP62E5s4zXO4KmBZvy0bPJAslojoGx6v25g5/+NDKbJf1uh5t1Pfyh+fv4K4D58BgYDbLSQNALS+mdic4tjvyEKejbNoKdgz6KyxnbG8QPM6awhrLpHbdbLCZ8HqbSzWereL1ejhw5ogTdBkrQP+b88P4qx4eCDIfre7bxxA1C4a0Dji3Lpaagf/8v8J07h7Ovca/0ahzu93Ok3w/FPuaT82ZEHTxJImF/4tDKbGlDtIag+194HuH3k3jzzZbWWJPpH4AnDAfOlQ/1j4ZYnU+i12kDW1zO4BraYUHPRGF5q7i5hvzInI6RsBcBP8zmCOkFBsMhDh48uBMrreCY38PdFgQd4NSpU6ysrLC0pArR66EE/WNMvmhwaSrKK8fq2y3FYpJ0+gGh0NNbji+ll+jx9ODRKlPsiisrZK9fb9lusfjk5ADxeIiZhNmwMxw+RyJxE8OwlzNtedYDVSwXAIfHQ+AnXib55vc7m744/RYcfhkcG7794OEgRlESnateYKSnChipAs5BX+fWsZ0jP2H+/mir7WK9Z8Gmjz6VyuBPJThz5kxHhkHb4bjfy2qhyFqhuXx5gNOnTyOE4IMPPtiBle0flKB/jLk8HSVT0G3YLbcASTi0NUJfTC/WzHBJvPE9AIKvvtrWGj9xfJB8rpel9BI5PUc4dA7DyJbbEDRi6WGcYK8HX6h2jnTo1VcpPn5M9laHZo0mFmD1Hox/YsvhA0fNwR4LD6q37S2WPOwdtVx6j0LoIEz/cMthV+k9izZ99AeJNOFMirNnzza+uEMc95uBw/0WovRQKMT4+DgffPDBztUd7AOUoH+MeeOjJdxOhw1BN9PcQuHKCL2W3ZJ4/XXcR47gadNffflYPxT7kUjmk/OEw6aFEU9UH3q8nYX76xw4VntCEkDwM58BTSPx+nfaWmuZ6VLWTMk/twj1e/GF3SxM1RJ0Mzp27WSELoT5oJn6C9iUxucIuxFurbyGeqR0nSUDDqDvSDOuWhz3m7bg3XRrG9hnzpwhGo2qWaN1UIL+MeZ7Hy3x8kQ/AY+z7nXxxA08ngNbNkSllMwl5zgYqPRPi2trpN59l9DnP9/2x/GIz8WJvnEAZhOz+HxHcDojxOPVZ2RuJhHNklzLcWCivqA7e3sJvPgi8de/3Znober7Jf9868g5IQQHjoZZfFB9QHVhKQ1Ox84UFW3m2GchtbQlH10IgXPIXqbL1aUoUgguDvXvmt0CMOZ14xaipQgd4KmnnsLhcHD9ur1g4ElECfrHlAfLSaZWUrz2VO0qT4tE4gPC2/zztdwa8Xyc8ch4xfXJ730PdJ3Q5z/XkbW+dtwsSrqx+AAhBOHwOeLxxj+UlrVxsEGEDhD6/OcpPHxErt1MCCnh7nfh2Ktmuf02DkxEWF/OkKmy+VhcSuMa8CLqDBjpCMc+Y/5+740th12DfgqLjQX9jbtmGf1PnppscGVncToE4z4P91qM0AOBAJOTk1y/fh1dt9eb/klDCfrHlO99ZO72v3qqvqDn81HS6amy1WExvT4NwHh4vOI18T/7Nq7RUbynt88Cb40vP30Sabh4d8b0zcPhc6RSd9D1+vbAwv11nG4H/aONh2qEfvI10DTif/bt9ha7eBMS83D8J6uetj4tLExVRun5+RSug43X2jahYRh+Gu5/b8th16EgRjyPXifTRUrJlcVVXIbB+YP1m7LtBJMBD/dajNABnnnmGVKpFPfu3evgqvYPStA/pnzn5iKnDoTqlvsDxGKXAOjpeWHL8en4NFAp6IXFJVI/+hHhn/qpjn0cPzoYxGUMcDf6EIBw+DxS6sTjN+q+7vH9dYbHw2g2+oo7+/oIvPQS8W9+E9lOifjdkg9//LWqpwePhHA4RMXGqB7PYSTyuGw8fDrC8c+YFaO5jTYK7tJ75+u0J5ibm2PO4eSIS6Dtot1iccznYTqTo2C0Zo1NTk7i9/u5erVyepNCCfrHksfrGS49jPJTTzfOH16LvYvD4a2M0OPTuBwuDgUPbTke/5NvgmEQ+fKXO7rmA4ER1guLrKXy9EQuAoK1WO3Wt4WczspssqF/vpnIz/wMhfl50pcqOxLa5t53zdzzcPV/W5dbY/BIiPk7a1uO50sFR+7R6vnyHef4a2AU4MGbG2s7FAQBhdnavXKuXLnCWjDC+T77/66d5HjAS1HCVKa1KF3TNM6dO8ft27dJJu21M36SUIL+MeRPrz9GSvivzh9qeG0s9h6RyLM4HFvT/qbXpzkcOoy2Kc9aSsn6H/0RvvPn8Ux0tjf208NHEa5V/uT6PC5XhFDoDGtrtYdTzN1ZQxqSQyfqz0jdTOi1z+IIBln/+tdbW2Q6aka9k9XtFouxp/pYnIpv6byYn02AANemFr87yuGXwdcLtzaGgzk8Gs5BX80IPZ1O896HH5FyezkT3qVPEtt4OmhmAF1P2G9TsJ2LFy9iGAZXrlT2hn/SUYL+MeSb1x9zdiTM0YH64lEoxP7/9s49PKrqXPi/NfckJJOZZBJIINwhgEK4CFJAEFABRVFPT9Far9VjW489fWoftZ6v1lOPPR61n/a0R7+2Wlvv2lYEi3KrLRULEggCIZCEkCthJvfJZTLX9f2xJzAkM8lMMrmA+/c8SXb2evfa7157zbvXvGvt96Wt7TiW1EU9ysqd5T0mRDuPFuIuKcV8443xVBeA2ZmTEFoPfzykTFpaLItpaSmI6EevPNaITq8ha2r0Bl2TkEDK2jU4t23D3xZddqHzKNoC0g8zek+4NW6GFSmh5sS5XJee6jb0mYloDLEFMes3Wj3kXgvFH4Pv3GjXkJ189ttCdwoKCqgzKi66GUmDvBInAtOSTCRqNRQ4+2/Q09PTmTx5Mvn5+erkaDdUg36BUdHQzhdVzayfHc3ofD8gSbWcb9B9AR9VrVWMTxl/3v6mN99EJCaSsm5tPFUGOPvwOFJ3gsqGDqyWryCll+bm8O6RysIGsqdb0OljM5CpN9+MdLlo2fxB7Eoe/SNYJ8GYOb2KZU5KQW/UUlmkhJ6VUuKtaUWfPUTuli5mbgC383y3S/YoAq0e/M7zXRp+v5/9+/cjx04AYPowGXStEMwelUDBAEboAAsXLqS1tZWioqI4aXZxoBr0C4y391ehEXB9Xt8GvbFxDxqNCXM3/3lNWw2+gO+8CVFfYyPOP/8Z8w3Xox2EUKq51lwANKbTvF9QQ2rqAoTQhXW7tNR10OJwkTMr9hgypjlzMF1yCU2vvxHb5GibQ4nfMusm5eWdXtBqNWRPt1AVNOj+ZjeBdh+GCOEJBo2Jy8FohsJNZ3ednRjtNko/duwYzc3NtI8dz2iDnjHGQYgGGSVzUxIpbHPFnI4ulKlTp2K1WtmzZ4/65mgIqkG/gPD4AryXX8XK3EzGmHt/G1FKSV39DtKsy9Bozo/VUuFUVptMNJ/zkze/9wekx4P161+Pv+JAekI6GYkZjLbV825+FRIT5pS5NDTu7iFbWagYypyZvceoCYcQAus3bsNTVkb7ZzEkkD72AcgAXHJzVOLjZlhx1rlotnfgqVCWMA7ZhGgXOgPMuE5xFbkVA67PGgVagbv83CqcQCDA3//+d9JtNoqklkWpSUP6QlF35qYk4Q5Iitr7H/JYo9GwdOlSamtr1SWMIagG/QJixzE79W0evn55Tp+yra1HcLvPYLP1nOA70aj4sbsMesDtpun110lcfDnGKVPiq3QIM9NmYkg8TU2zi51FDmy2q2lrK6Kjo/w8ufLD9aSkmzBn9O8V+uS1a9Gmp9P4299Gd4CUcPB3SiLmzOjW3k+co7x1W3rQQWdxE5pEnWJMh5p5t4OnFQqViWCNQYtxQgru4nP+/eLiYhwOB1MXL6HW7WWReYgmbiOQl6zc14H40QFmz55NSkoKu3fvVkfpQVSDfgHxu3+Uk52awBVTbX3K1tVtRwgt6ekre5QVOAqYZJ6E2agsXWt+9z18dXWk/8v98Vb5PGZaZ1LXWU2WRcOrn50iI2MNAA7HuZeBOpweqooambogs9+jSI3BQNqdd9C+Zw8dBQV9H1CdD2eOwGV3R32OZKuJMZPNlO4/Q2dJM8YpqYP/hmg4xi2C9Olw4NWzu4xTLXjPtONv9eD3+9m1axdWq5Xm0WMBuDx1eFa4dDHOZCBNrxuwQdfpdCxdupSqqipKSqIL9naxoxr0C4T95Y18fqqRu5dORBuF4XDU7SDVfBl6/fmJowMywKG6Q8zNmKv873bT8Otfk7BgPomLFoarKm7MTJuJRLI6z8feskZKG5JISZmDo+6jszIl+XakhGkLBxY0ynLLLWgtFup/8cu+hfNfBsMomP21mM4xZUEm3jMdBFo9mKZFTtA9qAgB8++Amnw4o4SWNU1VdOksbebQoUPU1dWxevVq9jtdmHVacodpQrQLIQSXpyaxu6l1wCPrefPmYbVa2b59u7riBdWgXzD8z19KSUsycOvCvt0tztajdHSUYguOgEMpay6j1dNKXkYeAE2vv4HP4cD2ne8Mul91RpoS0yVndBMpJh3P7ywhI2Mtra1HcbkqASj+3E76uFFYswbmFtAkJZH2zXto37OH9n2fRxZsc8DRP8GcjWCMzQc+ZX4GGXrlI2ScOkwGHWDOLaBPhM9+Dihr4TVJepxFdj755BPGjRvHjBkz2NvSxmXmJDTD6D/v4qq0FGrdXgrbos+DGg6dTsdVV11FfX29ui4d1aBfEBysbGJ3cR3fXDaJhCjWOdfUvIlGY2J0Zs+3PQvqFBfE3Iy5eB0O6n/5S0atWEHS4sVx17s7tgQbaaY0ypzHue+KSewsslPnXQpoqDn9Do2n23GUOwc8Ou/Ccuut6LOysD/5JNIbIdv8p88rb1wu+lbM9SemGMhJNdAmQQxGDtGoFbHCgrvhyHvQcBKhERinpPK3kn20t7ezZs0ayl0eSjvcLB1md0sXq9JSEMD2hvCRK2MhNzeXiRMnsnPnTlpawoc2/rKgGvQRTiAgeWJzIRnJRr6xeHyf8j5fK2fObCYzcz16fc/lh4cch7CarOQk5+B45lmk10vmo48Mhuo9EEKwYPQC9pzewzcW52BNMvDMzmZstquoqXmbgp0n0Ok15C6Oj0HXJCSQ8egjuEtKaHzjjZ4CzlrF3TLnFkiPfTLY1+xmlNtPjdtPab49DhoPgK88CFoDfPozAOoyOzkuq7ls+lyys7P5wKGEKlifEf2LWoOJzaBnbkoiO+oHbtCFEKxfvx4pJR9++OGXeoJUNegjnPcOVPFFdQs/XDeDUX3EPQeorf0TgYCLsdm39iiTUnLQfpA8Wx6t27bj3LKFtHvvxTC+7wdFvFg5biWNnY2cai3i+1dPY9+pRk62r8Pna8Zu38KMJVkkjIqcnShWklevJmn5FdQ9/wLu7svb/vY0BHxwxQ/6VXdH/hmEgGaziYKdVcNrSJIzlVH6oTdpLd3Lhwd3YiaRee4JALzvaGaROYlsU/zadqBclZZCQWsHDneEb08xYLVaWbVqFSUlJezdGzlG0MWOatBHMDXNLp7aepzLJli4IYoXifz+DsorXsJsnt8jGBfAsYZjVLdVc6XxUmoffxzT7Nmkf2twV7Z0Z9nYZeg0Ov5S+RduuSyHBeMtPLFdj/ROwjJtG7NXxTeDjhCCMT/5CZqEBGq+/xABd/ANyorP4MBvYeF9YI09bo0MSNrz7RgnpzLj6hwaqtsoP1wfV91jZvnD+BNs/PHdt+ns7OS6WavwnXBS6HByor2TDZnD6OcPw7U25dvCW7WNcalv0aJF5ObmsmPHDsrLy+NS54WGatBHKF5/gAffKsAfkDz71TlRTVhWVL6Mx+NgypSHw5ZvKt2E2Wdg1tMfgM9H9jP/jdAPre832ZDMotGL2FW5CyHgv26+lCSXpGrvOgzJdpyuP8T9nPqMDMb89CncJ05Q++ijyM422PwgpObAyn/vV52dxU34m90kXZbJ9MtHYxmTxKfvleDzDN9Ki4AxhU1p91PuSWX9dCMTrsiFgOS1ozVoBVxnG54Ii5GYlmRiuSWZV2rqBvTWaBdCCDZs2IDFYuHtt9/+UqaqUw36CERKyWPvH+FARRNP3XQp49P6XvHhclVRWfkrbLY1pJrn9yj3+D3sKt7Kjz9MxHfyFNkvvDCkrpZQVuaspLK1kqLGIianj+IOQwr1Zy6l0TOXslMv4HbXxf2cyStWkPGDh3Bu/Qj7/dch60th/c/BEPtqGumXtHx0Cq3VRMKsdLRaDVdsnIazvpMDH1fEXfdo8Pv9bNmyhSNVTlbZ6plz7KfoOw7SfImFN2UnN6eZsRmGceI2AveNs2H3+NjsaO5bOApMJhO33XYber2e11577Utn1FWDPsIIBCRPbDnGu/nVPLhyCtdHESLX73dz5OgDCKFj6pRHw8r8tXAL//r7RrJONDDmP59k1NIlYeWGgmsmXEOyPpkXD73I3s1leO2dkGfh2c/X4fW5OVr4IIHAwP2q3bHedRfWFVNp2muntnoFMmdpv+pp31+Lz95B6rqJCJ3yERo73cK0RZkc+KiciqMN8VS7T1wuF2+99RYFBQUsX76cZfc8BWmT4d3beXWKlwBwX1nkLEbDyZXWZKYmGvm/5XZc/oGP0gEsFgu33347Go2GV155heLi4rjUeyGgGvQRRJvbx7feOMCrn5Vz95KJfO+qaX0eEwh4KTr+CK2tR5k54xkSEsb2kGnYt4ekex9nSq1gzHPPkrphw2CoHzVmo5k7L7kT+wE3Bz+uYObSLL5/71zW5i3i5aMbaW7+nMKix5EyPh9wAPw+xNaHyMj8G+lXTaZlzwkqbvsGnqqqmKrxnmmn5aNyjJPMmGadH2tmxa25WLNHsf03R3FUDHz1RjSUlpby4osvUlZWxvr167nyyivBlAK3vsMey3ze6gjwzy4vqXsduI4Os48/DBoheGrqWE663DxVdjpu9dpsNu69916sVitvvvkmW7duxe3uf+q7CwUxXDPzCxYskPn5A8gscxEhpWRboZ0nthRid3by79fO5K4lE/r0m3u9LRwr+gH19buYPOkhJkw4fy21t6aGuv/5BS2bNnEmFaxPP8mc5dEFnxpMAv4An24q5siO0zTYKvjuD2/CnGBGSslLfyuj9OQzrJu4A2m8khWLXkCnG2DsEUcRbPo2nD4IS74Lq5/A+fHH1P7ocaTXS9rdd2G9884+o0z6Gjup+3+HkQFJxnfy0KUae8g4G1xseq4AV5uH1XfOZPK8vpN494fq6mp2795NcXExaWlp3HTTTWRnZ58tL3e5Wbf/OGnt1XyY/wBu3S/wdphJv+dSjBNGli8d4LHial6uqeeF3By+Nib2KJuR8Hq97Ny5k3379pGUlMTy5cvJy8vDYBg5q31iRQhxQEq5IGxZNAZdCLEGeAHQAr+RUv5Xt3Ij8HtgPtAAfE1KWd5bnapBh8Z2Dx8dreX1vZUU1TqZnpnMT2++lHk5va9G8Ps7sds3c7LsOTyeRqZPe5yxY29TytraaN/zGc4//5nWXbsICMmWBRLtnRv5wYofDcFVRSbgD1B+pIF9m8toPN2OdZ6GZw0PcUnmLJ6/8nksJuW695U18MGe57gy613avBZIvp9ls79GenIML8UEAlD5D8h/RYlznmCBa5+DS246K+KtrcXxzLM4t25FJCZivu46UtauIWH+fDQhH3gZkLgO19G06SRIie2+2RiyI+vS3uJm6/8exlHRyoTZ6Vx27QQyxg8sJLGUkoaGBkpKSjhy5AinT5/GaDSybNkyFi1ahD44uS2lZFu9k+8dr0QCH12SxcRPHsX/xTYc3p/hlzZSVthIXjkNEWOs+cGkwx/g9sNlfNrcxvfGZ/Kv4zNJjCKXbLRUV1ezbds2qqqqMJlMzJo1i9zcXMaNG4fJNLyhEGJlQAZdCKEFioGrgGpgP3CLlPJYiMy3gdlSyvuFEBuBG6WUvQbG+DIYdH9A0tbpw9nppcXlxdHaSVldOyfr2jhS00LhaacStyRzFPddMZkb8rLQh3RiKSWBQCdu9xk6XBV0dFTQ0pRPY9On+PxOkrSTGd/5TxhOG3CXltB5+DDuUsXoeM2J/CPPxFuzWli9cCOPLHwEnabvdezxQEqJzxugvcmNs96Fs96F/ZSTymONdDg9pKSbWHLzVCbmpbO9YjsP736YBF0CG3M3sjR7KeNTxjNKZ+bD/O24G55ldGIlrZ4kKtvnIPWXkmGZRKZlIrbk0WSYNKRoOjG6m9E0l0FDGTgKoXwPtDuUGC2X3aO8eJOUHlbfzqIiGn//Gs6PP0a6XGBMxDRzLsZp89Ak5xBoTyHgEugyjKRuyME41oroY4Tn9wX4YlcV+VvL8br9WLOSGJdrxZqdhNmWgNmWgDFJjy4YOsDn8+HxePB6vbS3t9PW1kZraysNDQ3Y7Xbsdjvt7UoWptGjR5OXl8fcuXMRej0Oj4/TnR7ynR28b2/iSJuLS0cl8NKs8UxODBqris8I7HiWprKFuAJL0GjaScg8g3GcAV32aLSjs9BYrYhEK2iHpp90xx0I8P3jVfzB3kSGQceNGRaWWZOZnmTCptdhGqCBl1JSVVXF559/zokTJ/AG3x7OzMwkKysLi8WCxWIhJSUFk8mEyWTCaDRiMBjQaEaOd3qgBn0x8GMp5TXB/x8FkFL+NERmW1DmH0IIHXAGsMleKu+vQf/e757mr2Mvjfm4UCS9uzKic0INYh0icg0SQcCvI+DXEwho46RH78g+l0zK+OsRUp0QfjRaH0ITQISt6dxeKfvSte/ln0PSP6JYhtp1HhHhhH40uDTnu6PG+iq50r2TJe6/o8cX9jit0GHUGNEKXQ+3XtdHVsal50TL+ecq0k7mA8M1HNFNxyfOrczRSy+JsgMdfgQSgUSDBCQaKc/2gvB9JNJZg1c6xOFtrj51iKe/+cN+HdubQY/mUZwNhM4cVQPdk1SelZFS+oQQLUAacN4sjBDiPuA+gJycvoNMhSOp08OYKJa19ev+hPSDnp0iXI29y0TsWNH0NymQUqP87TLgPgMyoAX8wLkJnoGbsMhGI0Shvuvo5djo74cM+d0TjcaL0HnRaL1otD4QEiEC5/4Spt27n1yG6NPt4SlDfkkCPQzb2eMijFW6mcc+JGTkZ7eUCCk5v+3ObUlAyADJgTZSfC2Y/S2Mc1dg9iuTsX70+Im8TLGdAAIPeo0BrdChEzo0QosQmuBZBALR44rOajKIBjCXJnJ5m05h4JQ+i9P6dFo1ibRrTHSIBPxCEzTngkBQv4AQZ/f1F0mIkRfn9sXyeItWzuQO/7AdKNEY9NgtWXgZpJS/An4Fygg9inP34Ml/+T/9OUxFRUXloicax1A1MC7k/7FA9/VFZ2WCLhczEJ/3eVVUVFRUoiIag74fmCqEmCiEMAAbgc3dZDYDdwS3/wn4S2/+cxUVFRWV+NOnyyXoE38A2IaybPEVKWWhEOI/gHwp5WbgZeA1IUQpysh842AqraKioqLSk6jWJ0kptwJbu+37Uch2J/DV+KqmoqKiohILI2dxpYqKiorKgFANuoqKispFgmrQVVRUVC4SVIOuoqKicpEwbNEWhRB1QH+zAaTT7S3UEYKqV2yoesXOSNVN1Ss2BqLXeCmlLVzBsBn0gSCEyI8Uy2A4UfWKDVWv2Bmpuql6xcZg6aW6XFRUVFQuElSDrqKionKRcKEa9F8NtwIRUPWKDVWv2Bmpuql6xcag6HVB+tBVVFRUVHpyoY7QVVRUVFS6oRp0FRUVlYuEEWvQhRBfFUIUCiECQogF3coeFUKUCiFOCCGuiXD8RCHEPiFEiRDinWDo33jr+I4Q4lDwp1wIcSiCXLkQ4khQbtATqQohfiyEqAnRbV0EuTXBNiwVQjwyBHo9I4Q4LoQ4LIR4XwiRGkFuSNqrr+sXQhiD97g02JcmDJYuIeccJ4T4RAhRFOz/3w0js0II0RJyf4ck+3df90Uo/DzYXoeFEPOGQKfpIe1wSAjhFEL8WzeZIWsvIcQrQgiHEOJoyD6rEGJH0BbtEEKEzQIvhLgjKFMihLgjnEyfSClH5A8wA5gO/BVYELJ/JvAFYAQmAicBbZjj3wU2BrdfAr41yPo+B/woQlk5kD6Ebfdj4KE+ZLTBtpsEGIJtOnOQ9boa0AW3nwaeHq72iub6gW8DLwW3NwLvDMG9GwPMC24noyRo767XCuDDoepP0d4XYB3wEUoGs8uBfUOsnxYln/H44Wov4ApgHnA0ZN9/A48Etx8J1+8BK1AW/GsJbltiPf+IHaFLKYuklCfCFN0AvC2ldEspTwGlwMJQAaFkvl0J/CG463fAhsHSNXi+fwbeGqxzDAILgVIpZZmU0gO8jdK2g4aUcruUsiuZ4l6U7FfDRTTXfwNK3wGlL60S3bMqxxkpZa2U8mBwuxUoQsnZeyFwA/B7qbAXSBVCjBnC868CTkop+/sG+oCRUu6mZ7a20H4UyRZdA+yQUjZKKZuAHcCaWM8/Yg16L4RLWt29w6cBzSHGI5xMPFkG2KWUJRHKJbBdCHEgmCh7KHgg+LX3lQhf8aJpx8HkbpTRXDiGor2iuf7zkp8DXcnPh4Sgi2cusC9M8WIhxBdCiI+EELOGSKW+7stw96mNRB5UDUd7dZEppawF5YENZISRiUvbRZXgYrAQQuwERocpekxK+UGkw8Ls61fS6miIUsdb6H10vkRKeVoIkQHsEEIcDz7J+01vegEvAj9BueafoLiD7u5eRZhjB7yGNZr2EkI8BviANyJUE/f2CqdqmH2D1o9iRQgxCvgj8G9SSme34oMoboW24PzIJmDqEKjV130ZzvYyANcDj4YpHq72ioW4tN2wGnQp5ep+HBZN0up6lK97uuDIKpxMXHQUSlLsm4D5vdRxOvjXIYR4H+Xr/oAMVLRtJ4T4NfBhmKJo2jHuegUne64DVsmg8zBMHXFvrzDEkvy8Wgxh8nMhhB63Tmf3AAAB90lEQVTFmL8hpfxT9/JQAy+l3CqE+F8hRLqUclCDUEVxXwalT0XJWuCglNLevWC42isEuxBijJSyNuiCcoSRqUbx9XcxFmX+MCYuRJfLZmBjcAXCRJQn7eehAkFD8QlKwmpQElhHGvEPlNXAcSlldbhCIUSSECK5axtlYvBoONl40c1veWOE80WT/Dveeq0BHgaul1J2RJAZqvYakcnPgz76l4EiKeXPIsiM7vLlCyEWonyOGwZZr2juy2bg9uBql8uBli5XwxAQ8VvycLRXN0L7USRbtA24WghhCbpIrw7ui42hmPnt52zxjShPLTdgB7aFlD2GskLhBLA2ZP9WICu4PQnF0JcC7wHGQdLzVeD+bvuygK0henwR/ClEcT0Mdtu9BhwBDgc705juegX/X4eyiuLkEOlViuInPBT8eam7XkPZXuGuH/gPlAcOgCnYd0qDfWnSELTRUpSv2odD2mkdcH9XPwMeCLbNFyiTy18ZAr3C3pduegngl8H2PELI6rRB1i0RxUCbQ/YNS3uhPFRqAW/Qft2DMu+yCygJ/rUGZRcAvwk59u5gXysF7urP+dVX/1VUVFQuEi5El4uKioqKShhUg66ioqJykaAadBUVFZWLBNWgq6ioqFwkqAZdRUVF5SJBNegqKioqFwmqQVdRUVG5SPj/pwsI5DliXe4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "samples = norminvgamma_rvs(mu_0, nu_0, alpha_0, beta_0, size = 10)  # YOU HAVE TO COMPLETE THIS\n",
    "\n",
    "# YOU HAVE TO PLOT THE NORMAL PDF CORRESPONDING TO EACH SAMPLE ABOVE\n",
    "# distribution = np.linspace(-15, 15, 100)\n",
    "\n",
    "# distribution = (samples[i,0], samples[i,1], mu_0, nu_0, alpha_0, beta_0)\n",
    "\n",
    "#for i in range(10):\n",
    " #   print(norminvgamma_pdf(samples[0,0], samples[0,1], mu_0, nu_0, alpha_0, beta_0))\n",
    "\n",
    "x = np.linspace(-10, 10, 200)\n",
    "\n",
    "for i in range(10):\n",
    "    sts.norm(samples[i,0], math.sqrt(samples[i,1]))\n",
    "    sts.norm(samples[i,0], math.sqrt(samples[i,1])).pdf(x)\n",
    "    plt.plot(x, sts.norm.pdf(x, samples[i,0], math.sqrt(samples[i,1])))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 2\n",
    "\n",
    "Generate 1,000,000 samples from the normal-inverse-gamma prior above and calculate\n",
    "approximate 95% confidence intervals over the mean and the variance using the\n",
    "samples. You can use the `numpy.percentile` function for this.\n",
    "\n",
    "Your confidence intervals should approximately match the intervals [-10, 10] and [0.1, 10]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-9.96781016 10.02077367]\n",
      "[ 0.10177289 10.05310737]\n"
     ]
    }
   ],
   "source": [
    "moresamples = norminvgamma_rvs(mu_0, nu_0, alpha_0, beta_0, size = 1000000) \n",
    "\n",
    "print(np.percentile(moresamples[:,0], [2.5, 97.5]))\n",
    "print(np.percentile(moresamples[:,1], [2.5, 97.5]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 3\n",
    "Code the equations for calculating the posterior normal-inverse-gamma hyperparameters\n",
    "from the prior hyperparameters and data.\n",
    "\n",
    "The equations are found on found [Wikipedia](https://en.wikipedia.org/wiki/Conjugate_prior#Table_of_conjugate_distributions) and reproduced below, using $d_i$ for a datum.\n",
    "Note that $n$ is the number of data, and $\\overline{d}$ is the mean of the data.\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "  \\mu_{\\text{post}} &= \\frac{\\nu_0\\mu_0 + n\\overline{d}}{\\nu_0 + n} \\\\\n",
    "  \\nu_{\\text{post}} &= \\nu_0 + n \\\\\n",
    "  \\alpha_{\\text{post}} &= \\alpha_0 + \\frac{n}{2} \\\\\n",
    "  \\beta_{\\text{post}} &=\n",
    "     \\beta_0 +\n",
    "     \\frac{1}{2}\\sum_{i=1}^n (d_i-\\overline{d})^2 +\n",
    "     \\frac{n\\nu_0}{\\nu_0+n}\\frac{(\\overline{d}-\\mu_0)^2}{2}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Once you have the update equations implemented, Bayesian inference is done!\n",
    "\n",
    " * The likelihood function is the normal distribution with unknown mean and variance.\n",
    " * The posterior distribution is of the same type as the prior – normal-inverse-gamma.\n",
    " * The posterior parameters below give you the exact posterior normal-inverse-gamma distribution.\n",
    " * No approximation or numerical integration is needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.array([1, 2, 3, 4])  # In class you will get a larger data set.\n",
    "                               # This is just to get you started.\n",
    "mu_post = ((nu_0*mu_0)+len(data)*np.mean(data))/(nu_0+len(data))\n",
    "nu_post = nu_0+len(data)\n",
    "alpha_post = alpha_0+(len(data)/2)\n",
    "beta_post = beta_0 + 0.5*(2.25+0.25+0.25+2.25)+(4*nu_0*(np.mean(data)-mu_0)**2)/((nu_0+len(data))*2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 4 (optional)\n",
    "\n",
    "You are told that the prior information we used above is incorrect. Actually, the prior 95%\n",
    "confidence interval on the mean should be [-15, 15] and on the variance [0.5, 2]. So, the prior\n",
    "over the mean is less certain (broader) than we had before, but the prior over the variance is\n",
    "more certain (narrower).\n",
    "\n",
    "Determine prior hyperparameters for the normal-inverse-gamma distribution that match the\n",
    "prior information above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}