Exact and approximate gamma samplers

Reparameterization trick through a gamma rejection sampler.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from autograd import grad\n",
    "import autograd.numpy as np\n",
    "import autograd.numpy.random as npr\n",
    "import autograd.scipy.special as sp\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "npr.seed(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def exact_sample_gamma(shape, scale, zs, us,\n",
    "                       error_on_failure=False):\n",
    "    \"\"\"\n",
    "    This takes in an (ideally) infinite stream of\n",
    "    random normals and uniforms and then generates\n",
    "    a rejection sampler for the gamma distribution\n",
    "\n",
    "    :param shape: shape parameter >= 1 (talk to me for details on shape < 1)\n",
    "    :param scale: scale of gamma distribution\n",
    "    :param zs:    infinite stream of normal random variates\n",
    "    :param us:    infinite stream of uniform random variates\n",
    "    \"\"\"\n",
    "    assert shape >= 1.0\n",
    "    \n",
    "    b = shape - 1/3.\n",
    "    c = 1. / np.sqrt(9*b) \n",
    "\n",
    "    # Rejection sample\n",
    "    for z, u in zip(zs, us):\n",
    "        v = 1 + c*z\n",
    "        if v < 0:\n",
    "            continue\n",
    "        v = v ** 3\n",
    "        \n",
    "        # First test the \"squeeze\" -- a simple function to determine acceptance\n",
    "        if (u < 1.0 - 0.0331 * z**4) or \\\n",
    "           (np.log(u) < 0.5 * z**2 + b * (1 - v + np.log(v))):\n",
    "            return b * v * scale\n",
    "\n",
    "    if error_on_failure:\n",
    "        raise Exception(\"Failed to accept!\")\n",
    "    else:\n",
    "        return b * v * scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Approximate this without the rejection step\n",
    "def approx_sample_gamma(shape, scale, z, u):\n",
    "    return exact_sample_gamma(shape, scale, [z], [u])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## First show that the truncated rejection sampler is very close to the exact sampler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape: 1.0\t shape 1.0\n",
      "Shape: 2.0\t shape 2.0\n",
      "Shape: 5.0\t shape 5.0\n",
      "Shape: 10.0\t shape 10.0\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaoAAAI7CAYAAAC9Xd+NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3X28VWWd///XGxHIGziKCQocCM1udNR0Ih1Dj5YKNopT\namllqDM1ahmald0J0q05jWY6Y/5SE/2RaVOgqWEzejTSjEKwRFNDONyr3N8pd5/vH2uBm+O+O3DO\nXmuf834+Hvtx9l7Xtdb+7Evcn31d61rXUkRgZmaWV92yDsDMzKwcJyozM8s1JyozM8s1JyozM8s1\nJyozM8s1JyozM8s1JyqrS5LGSroj6zg6E7ep5ZUTldWz3F0EKOlgSb+R9IqkzVXUP1zSnyStlTRN\n0mG1iLOMdmlTSRenn+c1SbdWUf9SSYskrZD0E0m7tkcc1jk4UZm1r43Az4HzK1VMv4wnAROAhvTv\nZEndOzTC2lgAfBO4pVJFSScDXwKOBwYDBwBXdWh0VlecqCzXJH1Z0nxJqyQ9K+n4guKekm5Py/4i\n6YhW+72Ylv1V0ukFZZ+SNFXSj9Jf8LMknVBQ3jv9Vb9Q0jxJ35SkauKNiOcj4jZgVhXVm4BdIuL6\niNgYET8CBJxQfrftPuOb2kZSN0lfTT//yrRnMyAtu05SS8H295c5/lGSfi9puaSnJB1XTVwAETEp\nIu4FllVR/Vzgloh4LiJWAuOB86p9L+v8nKgstyQdBFwMHBkRvYGTgTkFVU4FJgJ9gPuAGwvKXgSO\nSfe7CrhTUr+C8vcBLwB9gXHALyU1pGW3AxuAocB7gBOBf01jGiRpmaSB7fARDwaebrVtZrq9rApt\n8wXgo8CIiOhD0rtbl5b9ETgU2Iuk7e6R1KPI8QcAvwbGR8RewOXA/0jqm5Z/WdK91X/Usg4m+dxb\nzQT2lbRXOx3f6pwTleXZZqAHcIik7hHREhEvFZRPjYgpkSxYeQfJFzAAEfE/EbEkfX4PSVIaVrDv\nkrQnszki7gb+BnxI0r7ASODSiHgtIl4FrgPOTo81LyL2joj57fD59gBWttq2Ctizin3Ltc0FwNci\n4sU05r9ExPL0+cSIWBERWyLiWqAn8I4ix/84cH9ETEn3+z/gT8Ap6eurI+K0tnzYMlq3wyqSnmU1\n7WBdgBOV5VZE/B0YQ9LjWSJpoqT+BVUWFzxfB/SS1A1A0rnpcNVySctJfrXvU1B/Qau3mwvsT3KO\nZFdgUdpzWg7c1Grf9rIG6N1qWx9gdaUdK7TNIGB2sf0kXZ4OdW5tl94U/2yDgbPSNtjaDscA+1Xx\nudqqdTv0IZnUUbEdrGtworJci4i7ImI4yRcnwNWV9pHUCNwMXBQRe6VDV8+Q/ErfakCr3RqBhcA8\n4DWgb9pz2isiGiLiUNrfMxT0AlOHptsrKtM280gmJGwnPR/1ReCMgnbZ2ntpbR4wIW2Dre2wZ0R8\nv5rY2ugZoHC24+EkPd7lHfBeVoecqCy3JB0k6fj0HMoGYD2wpdwu6d/d03qvphMLzgMOaVV3X0mf\nk9Rd0pnAO4EHImIx8BBwraQ9lRgq6dg2xN2TZEhNknoWOweUagY2p3H0kHRJGvfD6XGOk1T081Zo\nm58A35R0YFr3HyTtTTKUthFYmr7flZQeXrsTOFXSSWkb9krj2b/KNthFUi9gF6B72g67lKg+AbhA\n0rvS81JfB26r5n2sa3CisjzrCXwPeIWkt/NW4Ctl6gdARDwL/AD4A8nw4MHA1FZ1nwTeDrxKMo36\nIwW/4M8lOf8zi2TW2j1Af9g2mWJVqckUkgaTJI2/pPGsB54rKH9A0hVpnBuB04FPAcvT9x0VEZvS\n6oOA35f4rOXa5j+Bu4GHJK0kSVy9gCnp43ngJZLh0nnFDp6egxsFfDV9j7kkEyq2Dq1+RdL9JWKD\nJNmsA75Mcr5rHfC1dN/t2jA9D/Z94JE0rr+TDGmaAaAsbpwo6Rbgn0m6928aUkmnwU7mjXH2X0bE\nt2oYonVikj4FXBARVfeSsiDpZuCeiPht1rGYZSmrCwtvA35E0uUv5bF2nFVkVnci4tNZx2CWB5kM\n/UXEVJKhjnKqusDSzMw6tzyfozpa0gxJ90t6d9bBWOcREbfnfdjPzN6Q1zXF/gw0RsQ6SSNJ1kM7\nqHUlSblblNTMzKoXERVHz3LZo4qINRGxLn3+ILBrOr22WF0/KjzGjh2beQz18HA7uY3cTrV9VCvL\nRCVKnIcqXJNN0jCS2YnVLG5pZmadTCZDf5Imkqwc3VdSCzCW5LqViIibgTMkXUhyceJ6kgU2zcys\nC8okUUXEORXKb2T7lbBtJzQ1NWUdQl1wO1XmNqqO26l9ZXLBb3uRFPUcv5lZVyaJqGIyRV5n/ZmZ\n1dyQIUOYO3du1mF0OoMHD2bOnDk7vL97VGZmqfQXftZhdDql2rXaHlUup6ebmZlt5URlZma55kRl\nZma55kRlZma55kRlZma55unpZmZlXHfllaxoaemw4zc0NjJm/PiK9YYMGcLLL79M9+7diQgkMXr0\naK6//vp2j+nRRx/lE5/4BPPmFb0BdM05UZmZlbGipYVxQ4Z02PHHVXl9kSTuv/9+jj/++A6LZaut\niTAvPPRnZlYnil2LdNFFF3HGGWdse/3lL3+ZE088EYAVK1Zw6qmnsu+++9K3b19OPfVUFi5cuK3u\n8uXLOf/88xkwYAB9+/blwx/+MOvWreOUU05h4cKF7LnnnvTu3ZvFixd3/Icrw4nKzKyO/eAHP+Cv\nf/0rEyZM4He/+x233XYbEyZMAGDLli2cf/75zJs3j5aWFnbbbTcuvvjibft+4hOfYP369Tz77LO8\n/PLLXHrppey22248+OCD7L///qxevZpVq1bRv3//rD4e4KE/M7O6cfrpp293juqaa67hggsuYMKE\nCYwcOZLevXtzww03sN9++wGw99578y//8i8A9OzZk6985St84AMfAGDRokVMmTKFZcuW0bt3bwCG\nDx+ezQeroO4T1QsvvFCybPDgwfTo0aOG0ZiZdZzJkycXPUc1bNgwhg4dyiuvvMKZZ565bfv69esZ\nM2YMU6ZMYcWKFUQEa9asISKYP38+e++997YklWd1n6iuvvqPRbevWbOUiy46nGOPPbbGEZmZdYxS\n6xDeeOONbNiwgf3335+rr76aK664AkiGBV944QWmTZvGW9/6VmbOnMkRRxxBRDBo0CCWLVvGqlWr\n3pSs8jSRAjpBoho48ONFt8+Z08yWLVtqHI2ZWW09//zzfOMb3+Cxxx6jV69eDBs2jFNOOYVDDz2U\n1atX85a3vIXevXuzbNkyxo0bt22//v37M3LkSC666CJuuOEG9thjD5544gmGDx9Ov379WLp0adEk\nloW6T1RmZh2pobGx6inkO3r8ap166qnssssu285RnXjiiSxYsICvfOUrHHLIIQB85zvf4ZOf/CR/\n+tOfGDNmDOeccw777LMPAwYM4Atf+AL33nvvtuPdcccdjBkzhne+851s3LiR448/nuHDh/OOd7yD\ns88+m6FDh7JlyxZmzZqV6YSKur/Nx9ixxeOfM6eZ0aN9p00zq55v89ExfJsPMzPr1JyozMws15yo\nzMws1zJJVJJukbRE0tNl6lwv6QVJMyQdXsv4zMwsP7LqUd0GnFyqUNJI4ICIeDvwGeCmWgVmZmb5\nkkmiioipwPIyVUYBE9K6TwJ9JPWrRWxmZpYveT1HNQAovBHKgnSbmZl1MXV/wW9z87htz4cMaWLI\nkKbMYjEzs9Kam5tpbm5u8355TVQLgEEFrwem296kqWlcLeIxM7Od1NTUtN0iDFdddVVV+2WZqJQ+\nirkXuBj4uaSjgBURsaRYxZdffrnoAZYvX86GDbu3R5xm1oVdeeV1tLSs6LDjNzY2MH78mKrrNzU1\n8fTTT7NkyRJ23XXXDosrTzJJVJImAk1AX0ktwFigBxARcXNEPCDpFEkvAmuB80od65Unnii6ffHy\nmcya1cBJJ53U7vGbWdfR0rKCIUPGddjx58yp/thz585l6tSpNDQ0cO+99/KRj3yk3ePZvHkzu+yy\nS7sfd2dkNevvnIjYPyJ6RkRjRNwWET+OiJsL6nw2Ig6MiMMiYnqpYx3cp0/Rx17d8zqqaWa2YyZM\nmMDRRx/N6NGj+elPf7pt+3nnnceFF17ISSedRO/evTn++ONpaWnZVt6tWzd+9KMfccABB7Dvvvvy\npS99aVvZ7bffzvvf/34uu+wy9tlnH6666ioigm9961sMGTKE/v37M3r0aFavXg3A3XffzdChQ1mz\nZg0ADz74IPvttx9Lly7tsM+d11l/ZmbWyoQJE/jEJz7BOeecw5QpU3jllVe2lU2cOJGxY8eydOlS\nDjvsMD7+8e1vgTRp0iSmT5/O9OnTmTx5Mrfeeuu2sieffJIDDzyQl19+ma997Wvbbmf/6KOPMnv2\nbFavXr3tFvZnnXUWxxxzDJdccgnLli3jX//1X7n11lvp27dvh31uJyozszowdepUWlpaOOusszji\niCM48MADmThx4rbyD33oQxxzzDHsuuuufPvb3+aJJ55gwYI35qBdccUV9OnTh4EDBzJmzBh+9rOf\nbSsbMGAAF110Ed26daNnz55MnDiRyy67jMGDB7Pbbrvx3e9+l7vuumvbPf5uuOEG/u///o+mpiZG\njRrFyJEjO/SzO1GZmdWBCRMmcNJJJ7HXXnsBcPbZZ3P77bdvKx806I2J0rvvvjt77703Cxcu3LZt\n4MCB254PHjx4u7LCfQEWLlzI4MGDt6u/adMmlixJ5rT16dOHM888k2eeeYbLLrusnT5haVWdyJHU\nHTgTODrdtDuwGVgHPA1MjIjXOiRCM7Mu7rXXXuPuu+9my5Yt7LfffgC8/vrrrFy5kqefTpZMnTfv\njTUS1qxZw7Jlyxgw4I11EubNm8e73vUuAFpaWth///23lbW+9fz+++/P3Llzt72eO3cuu+66K/36\nJQsEzZgxg1tvvZWzzz6bz33uczz44IPt/Im3VzFRSXovMBz4bUT8rEj5AcCnJc2MiEc7IMYd9uPr\n7+SnN/y6ZHm//fdgymOTahiRmVnb/epXv6J79+7MnDlzuynpZ511FhMmTADggQce4PHHH+cf//Ef\n+cY3vsHRRx+9XTK65pprGDZsGKtXr+aHP/whl19+ecn3O/vss/n+97/PiBEj2Gefffja177Gxz72\nMbp168Zrr73GJz/5Sb73ve8xevRo3vve9/Lf//3fXHjhhR32+avpUb0WEf9ZqjAi/g5cL2mopB4R\nsaH9wts5a1bBRw/5esnySfO/VcNozKweNTY2tGkK+Y4cv5IJEyZw/vnnb9dDArj44ov5/Oc/zwc/\n+EHOOeccxo0bxxNPPMGRRx7JnXfeuV3dUaNGceSRR7Jq1SrOO+88zj///JLvd/7557No0SKOPfZY\nXn/9dUaMGMH1118PwFe/+lUGDx7Mpz/9aSC5nf0JJ5zASSedxAEHHNDWj1+VNt2KXtLbgEV5GeaT\nFGOPe6Ro2R/mT+WZFZO44JD/KLn/pPnfYsaL/9tR4ZlZnanXW9Gfd955DBo0iPHjxxct79atGy++\n+CJDhw6tcWSJWt+K/nLgqPQNhkt6fxv3NzMza5O2XhX7R2CIpLdFxO8knd4RQdXKsqVLGTd6dNGy\nhsZGxpT4dWJmlietJ0O0tTzv2pqoBgGzgcskHQw8DtTvbIRNmxg3ZEjRonFz5tQ0FDOzHVV48W4x\nmzdvrlEkHaOtiWo28IuImCipL/DhDoipZla+vpHRk2YULXspVjKutuGYmVkRbU1UPwcOA6YDQ4H+\n7R5RDW3a0oshDcVXLZ7hGYFmZrlQNlFJ6gnsERFLASJiM0mSIiKmAdMK6g6KiHlFD2RmZraDyiaq\niHhd0omS9gQmRcT61nUkNQBnAbPY/vbxda3cRAvwZAuzzmjw4MF1P/EgjwqXY9oRFYf+IuLXkvoD\nl0raF+gF7EqyhNJaYD7wk4hYuVOR5E2ZiRbgyRZmndEc/3+dS1Wdo4qIxcB3OjiWXCk30QI82cLM\nrFbaNJlC0r8BZ5P0qO4ovNFhZ1NuogV4soWZWa20ddbf0og4QdJewGmSroiI73VEYHnni4XNzGqj\nrYmql6Qj0lvD3y7p1I4Iqh6sXLueOTNWFC176am5jHGeMjNrF21NVIcCR0j6NhDABkmrgUERcUe7\nR5djvgbLzKw22pqoJpOsuH65pB7APwL/BJwDdKlEVY6ntpuZtZ82JaqIeKLg+QaStf4el/TLthxH\n0gjgOpLV22+JiKtblR9HkhRnp5t+GRH1003x1HYzs3bT1h5VURExu3KthKRuwA3AB4CFwDRJkyPi\nuVZVH4uI09ojvlrz1HYzs/bTLomqjYYBL0TEXABJdwGjgNaJqm4vD/fUdjOz9pNFohrA9kstzSdJ\nXq0dLWkGsAD4YkTMqkVwteCp7WZm1csiUVXjz0BjRKyTNJLknlcHFavYPOen254PaTicIQ2H1yTA\nneGp7WbWFTU3N9Pc3Nzm/bJIVAuAxoLXA9Nt20TEmoLnD0r6L0l7R8Sy1gdrGjK6o+LsMJ7abmZd\nUVNTE01NTdteX3XVVVXtl0WimgYcKGkwsAj4GMmyTNtI6hcRS9Lnw0imxL8pSXVGntpuZra9mieq\niNgs6bPAQ7wxPf1ZSZ9JiuNm4AxJFwIbgfXAR2sdZ1bKDQuChwbNrOvJ5BxVRPwGeEerbT8ueH4j\ncGOt48qDSjMGH555qSdimFmXktfJFFZKmYuJfSGxmXVGTlR1ptzFxE8s+hv4/JaZdTJOVHWm3NDg\nw3//nJduMrNOx4mqE/HSTWbWGTlRdSKeiGFmnZETVVfiiRhmVoecqLqQckODjyyay5zR40ru29jY\nwPjxpXtrZmYdxYmqCyk3NPjqrMuYM6P0OomznrrNicrMMuFEZYktW2hqaChZfOvMFp/fMrNMOFFZ\nVbziu5llxYnKqlL2+q0yswnBPS4z2zlOVLbTKi2k+8SDT7KipaVomZOYmVXiRGU7rdL1W7+a9e8e\nNjSzHeZEZR2uXCL71R8u5PADP1hy337778GUxyZ1VGhmVgecqCxTmzb15PSBXy9Zfm2ZROYkZtY1\nOFFZrpVLZOWSGDiRmXUWTlRWt9wbM+sanKis09qZ3tiCVxYw4K0DipY5yZnVlhOVdUmVemPXzBnD\n6e/xkKNZHjhRmbXRzgw5luupgZOcWTFOVGbtrFwiK9dTAyc5s2IySVSSRgDXAd2AWyLi6iJ1rgdG\nAmuB0RFR+ta1VtacFTMY0lB6ZXRL5KGdskpyGze8ypknFP/shauHNDc309TUVPI4lnA7ta+aJypJ\n3YAbgA8AC4FpkiZHxHMFdUYCB0TE2yW9D7gJOKrWsXYWefgCrgf13k47leQeL716yH33PM5PJz4O\nwOJls+m/99DtyislwdUb1jL8hJOLlnXW+5w5UbWvLHpUw4AXImIugKS7gFHAcwV1RgETACLiSUl9\nJPWLiCU1j9asCyi3eshrr4/ZlgCbN/2UpoGjtyuvlASvebz0vc4m3fMN7p3465L7lkuClRLkzuxb\nLrnOnv03hg59R8l9Gxsb6NatZLHtgCwS1QBgXsHr+STJq1ydBem2NyWqqQsfKPomy9a/vFNBmlk7\nKXOvs2mv77LDsy8rJsid2bdMcp3+3K/ptuqjJfeddM836NVzLZPunPqmso5MrvWY1KuliNjpg7Tp\nDaWPACdHxKfT158AhkXEJQV17gO+GxGPp6//F/hSRExvdazaBm9mZu0qIlSpThY9qgVAY8Hrgem2\n1nUGVahT1Qc0M7P6lsVI6jTgQEmDJfUAPgbc26rOvcC5AJKOAlb4/JSZWddU8x5VRGyW9FngId6Y\nnv6spM8kxXFzRDwg6RRJL5JMTz+v1nGamVk+1PwclZmZWVvU7SRKSSMkPSfpeUlfzjqePJJ0i6Ql\nkp7OOpa8kjRQ0sOSnpH0F0mXVN6r65HUU9KTkp5K22ls1jHllaRukqZLan1Kw1KS5kiamf57+mPF\n+vXYo0ovGn6egouGgY8VXjRsIOn9wBpgQkQcmnU8eSSpP9A/ImZI2gP4MzDK/5beTNJuEbFO0i7A\n74FLIqLil0xXI+lS4Eigd0SclnU8eSRpNnBkRCyvpn699qi2XTQcERuBrRcNW4GImApU9Q+hq4qI\nxVuX54qINcCzJNfsWSsRsS592pPk/Hb9/crtYJIGAqcAP8k6lpwTbcg/9Zqoil007C8X2ymShgCH\nA09mG0k+pUNaTwGLgd9GxLSsY8qha4Ev4iReSQC/lTRN0r9VqlyvicqsXaXDfr8APp/2rKyViNgS\nEe8hua7xfZLenXVMeSLpQ8CStIeu9GHFHRMRR5D0Pi9OT1OUVK+JqpqLhs2qIqk7SZK6IyImZx1P\n3kXEKuARYETWseTMMcBp6fmXnwHHS5qQcUy5FBGL0r+vAL/izcvobadeE1U1Fw1bwr/sKrsVmBUR\nP8w6kLyStI+kPunztwAnsv1C0l1eRHw1IhojYijJd9LDEXFu1nHljaTd0hEMJO0OnAT8tdw+dZmo\nImIzsPWi4WeAuyLi2Wyjyh9JE4HHgYMktUjyhdOtSDoG+DhwQjpVdnp6vzTb3n7AI5JmkJzDmxIR\nxVeENiuvHzA1Pd/5B+C+iHio3A51OT3dzMy6jrrsUZmZWdfhRGVmZrnmRGVmZrnmRGVmZrnmRGVm\nZrnmRGVmZrnmRGVmZrnmRGVmZrnmRGVmZrnWPesAzGx76Y0JPwoMJbmdzTDgPyLipUwDM8uIe1Rm\n+XMYyWrus0kWFL4HWJRpRGYZcqIyy5mImB4RG4CjgUcjojkiXss6LrOsOFGZ5Yyk90rqCxwcES9J\nGp51TGZZ8jkqs/wZQXK798clnQ68mnE8ZpnybT7MzCzXPPRnZma55kRlZma55kRlZma55kRlZma5\n5kRlZma55kRlZma55kRlZma55kRlZma55kRlZma55kRldUnSWEl3ZB1HZ+I2tbxyorJ6lrv1vySd\nK+lPklZKapF0taSS/59JOjytv1bSNEmH1TLeIna6TSX1kPQTSXPSdpguaUSFfS6VtEjSinTfXXc2\nDus8nKjM2tdbgM8DfYH3AR8ALi9WMf0yngRMABrSv5Ml1fti0d2BFmB4RPQBvgHcLamxWGVJJwNf\nAo4HBgMHAFfVKFarA05UlmuSvixpvqRVkp6VdHxBcU9Jt6dlf5F0RKv9XkzL/pquQr617FOSpkr6\nUfoLfpakEwrKe6e/6hdKmifpm5JUTbwR8eOI+H1EbIqIRcD/DxxTonoTsEtEXB8RGyPiRyQ3Sjyh\nRP2q2kZSN0lfTT//yrSnNiAtuy7t6W3d/v4yxz9K0u8lLZf0lKTjqmyDdRExPiLmpa/vB14Cjiyx\ny7nALRHxXESsBMYD51XzXtY1OFFZbkk6CLgYODIiegMnA3MKqpwKTAT6APcBNxaUvQgck+53FXCn\npH4F5e8DXiDp+YwDfimpIS27HdhAciv49wAnAv+axjRI0jJJA6v8GMcCz5QoOxh4utW2men2siq0\nzRdIbmU/Iu3RnA+sS8v+CBwK7EXSdvdI6lHk+AOAXwPjI2Ivkl7h/6T3ydqaJO+tFGdatx/wdsq3\nw8yC1zOBfSXtVc3xrfNzorI82wz0AA6R1D0iWiLipYLyqRExJZJ71dxB8gUMQET8T0QsSZ/fQ5KU\nhhXsuyTtyWyOiLuBvwEfkrQvMBK4NCJei4hXgeuAs9NjzYuIvSNifqXgJZ1P0ov4jxJV9gBWttq2\nCtiz0rEp3zYXAF+LiBfTmP8SEcvT5xMjYkVEbImIa4GewDuKHP/jwP0RMSXd7/+APwGnpK+vjojT\nKgWZDmPeCfw0Ip4vUa11O6wi6VlW0w7WBThRWW5FxN+BMSQ9niWSJkrqX1BlccHzdUCvrRMX0kkN\nT6XDVstJfrXvU1B/Qau3mwvsT3KOZFdgUdpzWg7c1GrfitKhxm+T9GqWlai2BujdalsfYHWl41do\nm0HA7BJxXZ4OdW5tl94U/2yDgbPSNtjaDscA+1WKreC9RJKkXgc+V6Zq63boQzKpo2I7WNfgRGW5\nFhF3RcRwki9OgKsr7ZOetL8ZuCgi9kqHrp4h+ZW+1YBWuzUCC4F5wGtA37TntFdENETEoVQpneH2\nY+CfI2JWmarPUNALTB1K6SGy7ZRpm3kkExJax/V+4IvAGQXtsrX30to8YELaBlvbYc+I+H41saVu\nIUmCH46IzWXqPQMUznY8nKTHu7wN72WdmBOV5ZakgyQdn55D2QCsB7aU2yX9u3ta79V0YsF5wCGt\n6u4r6XOSuks6E3gn8EBELAYeAq6VtKcSQyUdW2XMJ5D0Ij4SEX+uUL0Z2JzG0UPSJWncD6fHOk5S\n0c9boW1+AnxT0oFp3X+QtDfJUNpGYGn6fldSenjtTuBUSSelbdgrjWf/KtvhJpI2PS0iNlSoPgG4\nQNK70vNSXwduq+Z9rGvIJFFJukXSEkmtTyRvLT8unY01PX18vdYxWi70BL4HvELS23kr8JUy9QMg\nIp4FfgD8gWR48GBgaqu6T5Kc4H8V+CZJYtn6C/5ckvM/s4BlwD1Af9g2mWJVmckUXycZxnpA0uq0\n7v1bCyU9IOmKNM6NwOnAp4Dl6fuOiohNafVBwO9LvE+5tvlP4G7gIUkrSRJXL2BK+nieZBbeOpKe\n05uk5+BGAV9N32MuyYSKrUOrXyn8XIXSHu2nSXtGBe1wdlq+XRum58G+DzySxvV3kiFNMwCUnIeu\n8ZsmQxBrSIYW3jSkkk6D/UI1J2vN2krSp4ALIqKqXlJWJN0M3BMRv806FrMsZXJhYURMlTS4QrWq\nrlsx66wi4tNZx2CWB3k+R3W0pBmS7pf07qyDMTOzbGQy9AeQ9qjuKzH0twewJSLWSRoJ/DAiDqp5\nkGZmlrlcrikWEWsKnj8o6b8k7d36ehRJuVuU1MzMqhcRFU/zZDn0J0qchypc6kbSMJKeX9GLJiPC\njwqPsWPHZh5DPTzcTm4jt1NtH9XKpEclaSLJgpx9JbUAY0mmA0dE3AycIelCkms+1pOsW2ZmZl1Q\nVrP+zqm/ghsFAAAgAElEQVRQfiPbLzBqZmZdVJ5n/Vk7aWpqyjqEuuB2qsxtVB23U/vKbNZfe5AU\n9Ry/mVlXJonI+WQKMzOzipyozMws15yozMws15yozMws15yozMws15yozMws15yozMws15yozMws\n15yozMws15yozMws15yozMws15yozMws15yozMws15yozMws15yozMws15yozMws15yozMws15yo\nzMws15yozMws1zJJVJJukbRE0tNl6lwv6QVJMyQdXsv4zMwsP7LqUd0GnFyqUNJI4ICIeDvwGeCm\nWgVmZmb5kkmiioipwPIyVUYBE9K6TwJ9JPWrRWxmZpYveT1HNQCYV/B6QbrNzMy6mLwmKjMzMwC6\nZx1ACQuAQQWvB6bb3mTcuHHbnjc1NdHU1NSRcZmZ2Q5qbm6mubm5zfspIto/mmreWBoC3BcR/1Ck\n7BTg4oj4kKSjgOsi4qgi9SKr+M3MbOdIIiJUqV4mPSpJE4EmoK+kFmAs0AOIiLg5Ih6QdIqkF4G1\nwHlZxGlmZtnLrEfVHtyjMjOrX9X2qDyZwszMcs2JyszMcs2JyszMcq2qyRSSugNnAkenm3YHNgPr\ngKeBiRHxWodEaGZmXVrFyRSS3gsMB34bEX8pUn4A8CFgZkQ82iFRlo7NkynMzOpUtZMpqklU/1As\nQRWpNxSYHxEbqg9z5zhRmZnVr3ZLVK0O+jZgUV6G+ZyozMzqV0dNT78cOCp9g+GS3r8jwZmZmVWr\nrYnqj8AQSW+LiN8B+3RATGZmZtu0NVENAjYAl0l6GPjH9g/JzMzsDW1d62828IuImCipL/DhDojJ\nzMxsm7Ymqp8DhwHTgaFA/3aPqE5cd+WVrGhpKVrW0NjImPHjaxyRmVnnVHbWn6SewB4RsbTigaRB\nETGvUr32lOWsv+MOO563qU/RsicW/Y2zR76v5L5OZGZm7XSbj4h4XdKJkvYEJkXE+iJv1ACcBcxi\n+9vHd2or1+7CkIFjipb9ata/M2fGipL7PvLgTGa0FD892NjYwPjxxY9rZtYVVRz6i4hfS+oPXCpp\nX6AXsCvJEkprgfnATyJiZYdGWmMnH3s6SxauKVn+0vxXkvsOF7FpSy+GNJRONq/Ouow5Mw4vWjbr\nqducqMzMClR1jioiFgPf6eBYcmXJwjWcPvDrJcuvmbMTyWTLFpoaGooWTZq/dsePa2bWCbVpMoWk\nfwPOJulR3RERN3dIVF3YsqVLGTd6dMlyn98ys66mrbP+lkbECZL2Ak6TdEVEfK8jAuuqVq5dX/b8\n1ktPzWWM85SZdSFtTVS9JB0REdOB2yWd2hFBdWWVzm89PPPSkj0u97bMrDNqa6I6FDhC0reBADZI\nWg0Miog72j26DlZuwkS5yRJZKtfjcm/LzDqjtiaqySTXXl0uqQfJEkr/BJwD1F2iKjdhYqcmS3Sg\ncj2uGfO/VeNozMw6XpsSVUQ8UfB8A/A48LikX7blOJJGANeRrDV4S0Rc3ar8OJKkODvd9MuI8Ldw\nBZ6IYWadUVt7VEVFxOzKtRKSugE3AB8AFgLTJE2OiOdaVX0sIk5rj/i6jE2bGDdkSMnicXPm1CwU\nM7P20i6Jqo2GAS9ExFwASXcBo4DWiarishq2vZWvb2T0pBkly1+KlYyrXThmZu0ii0Q1gO2XWppP\nkrxaO1rSDGAB8MWImFWL4OpZpRmDPodlZvUoi0RVjT8DjRGxTtJIYBJwULGK48aN2/a8qamJpqam\nWsRXlxYvXc7o0eOKlnmNQTPraM3NzTQ3N7d5vywS1QKgseD1wHTbNhGxpuD5g5L+S9LeEbGs9cEK\nE5WVt3pteI1BM8tM687EVVddVdV+WSSqacCBkgYDi4CPkSzLtI2kfhGxJH0+jGRK/JuSVCVXXnkd\nLS2lV3lYtHRNLq+V6jBeY9DM6lDNE1VEbJb0WeAh3pie/qykzyTFcTNwhqQLgY3AeuCjO/Jev5n8\nCLvpvJLla9feuyOHNTOzGsrkHFVE/AZ4R6ttPy54fiNw486+z4a1azllYPEeBMC0LVt29i06jZfm\nL+DwAz9Ysrzf/nsw5bFJNYzIzCyR18kUVmObNvUse1uTSZ4xaGYZcaKyqnjGoJllxYnKquIZg2aW\nFScqq45nDJpZRpyobKd5IoaZdSQnKttpnohhZh2pW9YBmJmZleMelXW4ckODHhY0s0qcqKzDlRsa\n9LCgmVVS94lq7drSM84iooaR2I7wRAwzq6TuE9XnPvdfRbdv2rSBda+/XuNorK08EcPMKqn7RNXY\n+MWi2+fMaSbigRpHY+3N57fMrO4TlXVuPr9lZk5UVreWLV3KuNGjS5Y3NDYyZvz42gVkZh2i7hPV\nX//0p6LbX1n6NzZt2lTjaKyWVq5dz5wZpW+M+dJTcxnjPGVW9+o+UfVbsqTo9ldeXsjmzZtrHI3V\n0qYtvRjSUHox3F/94UKf3zLrBOo+Ub11992Lbn/LrrvWOBLLm3Lnt26deamHDc3qRN0nKrMd4WFD\ns/rhRGVd0s4MG27c8CpnnlD83lzg3phZe3OiMiui3LDhtY//e9ne2BMPPsmKlpaiZU5iZm3nRGXW\nRhV7Y7NKJ7JHHpzJjJbSNy1obGzw3ZLNWskkUUkaAVxHcpuRWyLi6iJ1rgdGAmuB0RExo7ZRdh5z\nVsxgSEPpoSpLtFc7lUtkr866jDkzSr/HpHu+wb0Tf120LA8zFZubm2lqaso0hnrgdmpfNU9UkroB\nNwAfABYC0yRNjojnCuqMBA6IiLdLeh9wE3BUrWPtLJyoqlOTdtqyhaaGhpLF017fpfSQY5nzZgAL\nXlnAgLcOKFrWXknOX8DVcTu1ryx6VMOAFyJiLoCku4BRwHMFdUYBEwAi4klJfST1i4jiF02ZdQGV\nFvC9Zs4YTn9P+ye5wrLFy2Yz6c6p25XnoadnnVsWiWoAMK/g9XyS5FWuzoJ025sS1QtLny/6Jqte\nX75TQZp1JjuT5ArLmjf9lKaBo7crb68k2Nby1RvWMvyEk0vu6/N9nYdqfc8mSR8BTo6IT6evPwEM\ni4hLCurcB3w3Ih5PX/8v8KWImN7qWL7hlJlZHYsIVaqTRY9qAdBY8Hpguq11nUEV6lT1Ac3MrL6V\nnifbcaYBB0oaLKkH8DHg3lZ17gXOBZB0FLDC56fMzLqmmveoImKzpM8CD/HG9PRnJX0mKY6bI+IB\nSadIepFkevp5tY7TzMzyoebnqMzMzNoii6G/diFphKTnJD0v6ctZx5NHkm6RtETS01nHkleSBkp6\nWNIzkv4i6ZLKe3U9knpKelLSU2k7jc06pryS1E3SdEmtT2lYStIcSTPTf09/rFi/HntU6UXDz1Nw\n0TDwscKLhg0kvR9YA0yIiEOzjiePJPUH+kfEDEl7AH8GRvnf0ptJ2i0i1knaBfg9cElEVPyS6Wok\nXQocCfSOiNOyjiePJM0GjoyIqq4jqtce1baLhiNiI7D1omErEBFTAV9QVkZELN66PFdErAGeJblm\nz1qJiHXp054k57fr71duB5M0EDgF+EnWseScaEP+qddEVeyiYX+52E6RNAQ4HHgy20jyKR3SegpY\nDPw2IqZlHVMOXQt8ESfxSgL4raRpkv6tUuV6TVRm7Sod9vsF8Pm0Z2WtRMSWiHgPyXWN75P07qxj\nyhNJHwKWpD10pQ8r7piIOIKk93lxepqipHpNVNVcNGxWFUndSZLUHRExOet48i4iVgGPACOyjiVn\njgFOS8+//Aw4XtKEjGPKpYhYlP59BfgVb15Gbzv1mqiquWjYEv5lV9mtwKyI+GHWgeSVpH0k9Umf\nvwU4ke0Xku7yIuKrEdEYEUNJvpMejohzs44rbyTtlo5gIGl34CTgr+X2qctEFRGbga0XDT8D3BUR\nz2YbVf5Imgg8DhwkqUWSL5xuRdIxwMeBE9KpstPT+6XZ9vYDHpE0g+Qc3pSIeCDjmKw+9QOmpuc7\n/wDcFxEPlduhLqenm5lZ11GXPSozM+s6nKjMzCzXnKjMzCzXnKjMzCzXnKjMzCzXnKjMzCzXnKjM\nzCzXnKjMzCzXnKjMzCzXumcdgJltL70x4UeBoSS3sxkG/EdEvJRpYGYZcY/KLH8OI1nNfTbJgsL3\nAIsyjcgsQ05UZjkTEdMjYgNwNPBoRDRHxGtZx2WWFScqs5yR9F5JfYGDI+IlScOzjsksSz5HZZY/\nI0hu9/64pNOBVzOOxyxTvs2HmZnlmof+zMws15yozMws15yozMws15yozMws15yozMws15yozMws\n15yozMws15yozMws15yozMws15yorC5JGivpjqzj6EzcppZXTlRWz3K3/pekT0naJGmVpNXp32PL\n1D9c0p8krZU0TdJhtYy3iHZpU0nNktYXtMOzFepfKmmRpBWSfiJp1/aIwzoHJyqz9vd4RPSOiD3T\nv48Vq5R+GU8CJgAN6d/JkjrDYtEBXFTQDu8qVVHSycCXgOOBwcABwFW1CdPqgROV5ZqkL0uan/4y\nf1bS8QXFPSXdnpb9RdIRrfZ7MS37a7oK+dayT0maKulH6S/4WZJOKCjvnf6qXyhpnqRvSlIHfLwm\nYJeIuD4iNkbEj0hulHhC+d22xVm0bSR1k/TV9POvTHtqA9Ky6yS1FGx/f5njHyXp95KWS3pK0nFt\n/HzVttm5wC0R8VxErATGA+e18b2sE3OistySdBBwMXBkRPQGTgbmFFQ5FZgI9AHuA24sKHsROCbd\n7yrgTkn9CsrfB7wA9AXGAb+U1JCW3Q5sILkV/HuAE4F/TWMaJGmZpIFlQn+PpJclPSfp65JK/X92\nMPB0q20z0+1lVWibL5Dcyn5ERPQBzgfWpWV/BA4F9iJpu3sk9Shy/AHAr4HxEbEXcDnwP+l9srYm\nyXsrhPndtB1+VyHJHUzyubeaCewraa8Kx7cuwonK8mwz0AM4RFL3iGiJiJcKyqdGxJRI7lVzB8kX\nMAAR8T8RsSR9fg9JUhpWsO+StCezOSLuBv4GfEjSvsBI4NKIeC0iXgWuA85OjzUvIvaOiPklYn4U\nOCQi9gU+ku73xRJ19wBWttq2CtizbKskyrXNBcDXIuLFNOa/RMTy9PnEiFgREVsi4lqgJ/COIsf/\nOHB/RExJ9/s/4E/AKenrqyPitDLxfYkk0Q8A/j/gPklvK1G3dTusIumNVdMO1gU4UVluRcTfgTEk\nPZ4lkiZK6l9QZXHB83VAr629F0nnpsNVyyUtJ/nVvk9B/QWt3m4usD/JOZJdgUVpz2k5cFOrfcvF\nPCci5qbPnyEZxjqjRPU1QO9W2/oAq6t4n3JtMwiYXWw/SZenQ51b26U3xT/bYOCstA22tsMxwH6V\nYkvjmxYRa9MhzQnA70mTXBGt26EPyTmuiu1gXYMTleVaRNwVEcNJvjgBrq60j6RG4GaSk/l7pUNX\nz7D9OZMBrXZrBBYC84DXgL5pz2mviGiIiEPZcaXO1TxDQS8wdWi6vaIybTOPZELC9kEk56O+CJxR\n0C5bey+tzQMmpG2wtR32jIjvVxNbsXBLvA8kn7dwtuPhJD3e5Tv4XtbJOFFZbkk6SNLx6TmUDcB6\nYEu5XdK/u6f1Xk0nFpwHHNKq7r6SPiepu6QzgXcCD0TEYuAh4FpJeyoxtNwU81Yxj0iHD5H0TuDr\nJDP7imkGNqdx9JB0SRr3w+n+x0kq+nkrtM1PgG9KOjCt+w+S9iYZStsILE3f70pKD6/dCZwq6aS0\nDXul8exfRRv0SffrKWkXSR8HhgO/KbHLBOACSe9Kz0t9Hbit0vtY1+FEZXnWE/ge8ApJb+etwFfK\n1A+AiHgW+AHwB5LhwYOBqa3qPgm8HXgV+CbwkYJf8OeSnP+ZBSwD7gH6w7bJFKvKTKb4APC0pNUk\nkxF+AXx3a6GkByRdkca5ETgd+BSwPH3fURGxKa0+iGTIrJhybfOfwN3AQ5JWkiSuXsCU9PE88BLJ\ncOm8YgdPz8GNAr6avsdckgkVW4dWvyLp/hKx7Qp8C3g53ffi9HO9mO67XRum58G+DzySxvV3kiFN\nMwCUnIeu8ZtKtwD/TNK9f9OQSjpDaDJvjLP/MiK+VcMQrROT9CnggoioqpeUFUk3A/dExG+zjsUs\nS1ldWHgb8COSLn8pj1WYVWTWqUXEp7OOwSwPMhn6i4ipJEMd5XTEBZZmZlZn8nyO6mhJMyTdL+nd\nWQdjnUdE3J73YT8ze0Ne1xT7M9AYEeskjSSZNXVQ60qScrcoqZmZVS8iKo6e5bJHFRFrImJd+vxB\nYNd0em2xun5UeIwdOzbzGOrh4XZyG7mdavuoVpaJSpQ4D1W4JpukYSSzE5fVKjAzM8uPTIb+JE0k\nWTm6r6QWYCzJdSsRETcDZ0i6kOTixPUkC2yamVkXlEmiiohzKpTfyPYrYdtOaGpqyjqEuuB2qsxt\nVB23U/vK5ILf9iIp6jl+M7OuTBJRr5MpzMzMtnKiMjOzXHOiMjOzXHOiMjOzXHOiMjOzXHOiMjOz\nXHOiMjOzXHOiMjOzXHOiMjOzXHOiMjOzXHOiMjOzXHOiMjOzXHOiMjOzXHOiMjOzXHOiMjOzXHOi\nMjOzXHOiMjOzXHOiMjOzXHOiMjOzXHOiMjOzXMskUUm6RdISSU+XqXO9pBckzZB0eC3jMzOz/Miq\nR3UbcHKpQkkjgQMi4u3AZ4CbahWYmZnlSyaJKiKmAsvLVBkFTEjrPgn0kdSvFrGZmVm+5PUc1QBg\nXsHrBek2MzPrYrpnHcDOGjdu3LbnTU1NNDU1ZRaLmZmV1tzcTHNzc5v3U0S0fzTVvLE0GLgvIg4t\nUnYT8EhE/Dx9/RxwXEQsaVUvsorfzMx2jiQiQpXqZdmjUvoo5l7gYuDnko4CVrROUp3VlVdeR0vL\nipLlv3t4Cnv22L1oWb/992DKY5M6KjQzs0xkkqgkTQSagL6SWoCxQA8gIuLmiHhA0imSXgTWAudl\nEWcWfjP5EXZT6Y+7ZPG9fPKfvlu0bNL8b3VUWGZmmckkUUXEOVXU+WwtYsmbDWvXcsrAhpLl07Zs\nqWE0ZmbZq/vJFPXo5GNPZ8nCNUXLXpr/CgzcseMuXrqc0aPHlSxvbGxg/PgxO3ZwM7OMOFFlYMnC\nNZw+8OtFy66Zs+OJZPXaYM6M0ot4zHrqNicqM6s7VSUqSd2BM4Gj0027A5uBdcDTwMSIeK1DIrTq\nbdlCU0PpYcNJ89fWMBgzs/ZRMVFJei8wHPhtRPysSPkBwKclzYyIRzsgRjMz68Kq6VG9FhH/Waow\nIv4OXC9pqKQeEbGh/cIzM7OuruISShHxl63PJb1NUq8S9WY7SZmZWXtr62SKy4F7gGZJw0mue5ra\n/mFZR1i2dCnjRo8uWtbQ2MiY8eNrG5CZWRXamqj+CAyR9LaI+J2k0zsiKOsYK9euZ86M4qtevPTU\nXMY4T5lZDrU1UQ0CZgOXSToYeBzwmj2tlLtOCnbuWqmdsWlLL4Y0FJ+ePsOrWphZTrU1Uc0GfhER\nEyX1BT7cATHVvXLXScHOXStlZtbVtPV+VD8HDkmfDwX6t284ZmZm2yvbo5LUE9gjIpYCRMRmYHr6\nfBowraDuoIiYV/RAZmZmO6hsjyoiXgeOlnS2pLcUqyOpQdKngcEdEaCZmXVtFc9RRcSvJfUHLpW0\nL9AL2JVkCaW1wHzgJxGxskMjtQ710vwFHH7gB0uW+15XZpaVqiZTRMRi4DsdHItlaNOmnmUngPhe\nV2aWlTZNppD0b5IelvS7dLjPzMysQ7V11t/SiDgBOA14XdIVHRCTmZnZNm1NVL0kHRERyyPiduCZ\njgjKzMxsq7Ze8HsocISkbwMBbJC0GhgUEXe0e3RmZtbltTVRTQYUEZdL6gH8I/BPwDmAE1Un5gVt\nzSwrbUpUEfFEwfMNJGv9PS7pl205jqQRwHUkQ4+3RMTVrcqPI0mKs9NNv4yIXE07K7eeX1Zr+XUk\nL2hrZllpa4+qqIiYXblWQlI34AbgA8BCYJqkyRHxXKuqj0XEae0RX0cot55fZ1zLzwvamllW2jqZ\noj0MA16IiLkRsRG4CxhVpJ5qG5aZmeVRFolqAFC4JuD8dFtrR0uaIel+Se+uTWhmZpY37TL01wH+\nDDRGxDpJI0nueXVQxjGZmVkGskhUC4DGgtcD023bRMSagucPSvovSXtHxLLWBxs3bty2501NTTQ1\nNbV3vGZm1g6am5tpbm5u835ZJKppwIGSBgOLgI8BZxdWkNQvIpakz4eRTIl/U5KC7ROVZcML2ppZ\nNVp3Jq666qqq9qt5ooqIzZI+CzzEG9PTn5X0maQ4bgbOkHQhsBFYD3y01nFa9bygrZl1pEzOUUXE\nb4B3tNr244LnNwI31jouMzPLnyxm/ZmZmVXNicrMzHLNicrMzHItr9dRWSdSblagZwSaWSVOVCVc\neeV1tLQUX4QVYNHSNZ1u4dmOUm5WoGcEmlklTlQl/GbyI+ym80qWr117bw2jMTPrupyoStiwdi2n\nDGwoWT5ty5YaRmNm1nU5UVmmyt2QEXxTRjNzorKMlbshI/imjGbmRGUZK3dDRvBNGc3M11GZmVnO\nOVGZmVmueejPcs0XC5uZE5Xlmi8WNjMP/ZmZWa65R2V1y3cWNusaunSiOvnY01mycE3Rspfmv+K1\n/HLOdxY26xq6dKJasnBNyS+6a+aUvrbH6sPipcsZPXpc0bLGxgbGj/d/Y7N60KUTlXVuq9cGc2Yc\nXrRs1lO3OVGZ1QknKuu8tmyhqaH4wsKT5q+tcTBmtqOcqKxL8kQMs/qRSaKSNAK4jmR6/C0RcXWR\nOtcDI4G1wOiImFHbKK0zqzQR49o/XOgLjc1youaJSlI34AbgA8BCYJqkyRHxXEGdkcABEfF2Se8D\nbgKOqnWsncWcFTMY0lD8XI29obCdyiWyckkMOncia25upqmpKeswcs/t1L6y6FENA16IiLkAku4C\nRgHPFdQZBUwAiIgnJfWR1C8ilrTljcpNP4euMwXdiao61bZTpd7YTTO/0GlnG/oLuDpup/aVRaIa\nAMwreD2fJHmVq7Mg3famRDVr1qySb7Ro/io+3HhlyXJPQbeOUG624aR7vsG9E39dct8FryxgwFsH\nFC3rzD01s3LqfjLF+Z+6pGTZmvXraxiJWarMbMNpr+9Stjd2zZwxnP6eHRtyLJfkypW1Zd/Fy2Yz\n6c6pHXLsPO27ccOrnHlC8R8bf5s9m3cMHVpy3ydmr2DJqhXMmTPuTWW/e3gKe/bYveS+/jFSnCKi\ntm8oHQWMi4gR6esrgCicUCHpJuCRiPh5+vo54LjWQ3+Sahu8mZm1q4hQpTpZ9KimAQdKGgwsAj4G\nnN2qzr3AxcDP08S2otj5qWo+oJmZ1beaJ6qI2Czps8BDvDE9/VlJn0mK4+aIeEDSKZJeJJmefl6t\n4zQzs3yo+dCfmZlZW9Tt/agkjZD0nKTnJX0563jySNItkpZIejrrWPJK0kBJD0t6RtJfJJWendOF\nSeop6UlJT6XtNDbrmPJKUjdJ0yXdm3UseSVpjqSZ6b+nP1asX489qvSi4ecpuGgY+FjhRcMGkt4P\nrAEmRMShWceTR5L6A/0jYoakPYA/A6P8b+nNJO0WEesk7QL8HrgkIip+yXQ1ki4FjgR6R8RpWceT\nR5JmA0dGxPJq6tdrj2rbRcMRsRHYetGwFYiIqUBV/xC6qohYvHV5rohYAzxLcs2etRIR69KnPUnO\nb9ffr9wOJmkgcArwk6xjyTnRhvxTr4mq2EXD/nKxnSJpCHA48GS2keRTOqT1FLAY+G1ETMs6phy6\nFvgiTuKVBPBbSdMk/VulyvWaqMzaVTrs9wvg82nPylqJiC0R8R6ShcfeJ+ndWceUJ5I+BCxJe+hK\nH1bcMRFxBEnv8+L0NEVJ9ZqoFgCNBa8HptvM2kxSd5IkdUdETM46nryLiFXAI8CIrGPJmWOA09Lz\nLz8Djpc0IeOYcikiFqV/XwF+xZuX0dtOvSaqbRcNS+pBctGwZ9gU5192ld0KzIqIH2YdSF5J2kdS\nn/T5W4AT2X4h6S4vIr4aEY0RMZTkO+nhiDg367jyRtJu6QgGknYHTgL+Wm6fukxUEbEZ2HrR8DPA\nXRHxbLZR5Y+kicDjwEGSWiT5wulWJB0DfBw4IZ0qOz29X5ptbz/gEUkzSM7hTYmIBzKOyepTP2Bq\ner7zD8B9EfFQuR3qcnq6mZl1HXXZozIzs67DicrMzHLNicrMzHLNicrMzHLNicrMzHLNicrMzHLN\nicrMzHLNicrMzHLNicrMzHKte9YBmNn20hsTfhQYSnI7m2HAf0TES5kGZpYR96jM8ucwktXcZ5Ms\nKHwPsCjTiMwy5ERlljMRMT0iNgBHA49GRHNEvJZ1XGZZcaIyyxlJ75XUFzg4Il6SNDzrmMyy5HNU\nZvkzguR2749LOh14NeN4zDLl23yYmVmueejPzMxyzYnKzMxyzYnKzMxyzYnKzMxyzYnKzMxyzYnK\nzMxyzYnKzMxyzYnKzMxyzYnKzMxyzYnKckvSWEl3ZB1HVyLpJUknZB2HWSEnKvt/7d17mB11fcfx\n9wdDkkKAQIQg2SyRW9VYiKBRRHDx0hCsQL3ipRhqhVosDdRblQeSaKlULYjGIhUpyBNjEQhgEVBh\nwZCCeSCLSEAuuZHbGnKD3QRz+/aPmQ0nZ8+cPRt2z8zZ/byeZ5+cM7/fzPnOsJzv/m4zRVe4e3xJ\nGi/pTklrJG2vUL6/pFskdaRf/B/r4XgXSFolaYOkH0ras/+irw9Je0q6MT3/HZJOqlDnMknPp9fx\nGz0c792Snkiv6a8lNfdf9FY0TlRmvbcV+Cnwtxnl3wdeAg4EPgn8p6TXV6ooaRLwReBk4FDgcGB6\nXweck98An6DCs7QknQucBvwFcDTwfknnVDpIeif5m4CvAgcAD5NcfxsknKgsd5K+JGm5pBfSv5pP\nLikeJum6tOwxSceW7fdMWvb79E7jXWWfkjRX0nfTlsrC0i4tSfumrZeVkp6T9DVJqiXeiHgqIq4F\nFkhaG4YAABGISURBVFY4l72ADwAXRcTmiHgAuBX4m4zDnQVcExFPRsRGYAZwdi1xSBol6XZJ6yWt\nlXRfSVmTpJsk/TFtsVyZbj8sbZE8n5bdIGnfjONL0pfTa7xG0mxJI2uJLSK2RsSVETEP2JFx3t+O\niFURsQr4FjAl43AfAH4fETenz+maBhwj6ahaYrHG50RluUq/bM4DjouIfYFJwJKSKu8HZgH7AbcD\nM0vKngFOSPebDtwgaXRJ+VuBp4FRJF9uN5d80V4HbCF53PubgPcCf5fGNFbSOklNu3FKRwFbI+LZ\nkm2PAuMz6o9Py0vrHiRp/xo+659JHlU/CjgI+AqApD2AnwOLgWZgDDA73UfApcDBwOuBJpJrU8n5\nJK2eE4FDgPUkrUXSz3lU0pk1xFlJpfOu6RpFxCaS//ZZ9W2AcaKyvG0HhgJvlDQkIpZFxOKS8rkR\ncVckz6P5MUk3EQARcVNEtKevbyRJShNL9m1P/6rfHhH/A/wBeJ+kg4DJwAUR8VJEPA9cAXwsPdZz\nEXFARCzfjfMZAbxQtu0FYJ8q9TeW1VWV+qW2Aq8BXpue4wPp9onp9i+m57clbdkQEc9GxK8jYltE\nrAUuB96Zcfxzga+mrZ6tJK29D6WJkIg4JiJmZ+zbk0rnPaLGul31a7lGNgA4UVmu0pbHVJK/6tsl\nzZJ0cEmV1SWvNwHDu74oJZ0laUHa9bWe5C/sV5fUX1H2cUtJWgaHAnsCq9KW03rgqrJ9d1cHUN6V\nth/wYo319yOZQJJVv9S/A88Cd6fdc19Kt48FlkZEty43SQdJ+kna1boBuIHs8z4UuCW9RutIujq3\nAqMz6vdGpfPuqLFuV/1arpENAE5UlruImB0RJ5J8MQJc1tM+6ayvq4F/iIj9I2J/4HGS1kiXMWW7\nNQMrSbrLXgJGpS2n/SNiZEQczSv3FDBE0uEl245JY6vk8bS8ywSSluD6nj4oIjoj4vMRcThJF92F\n6fjec0BzV0IvcynJmNH4iBhJMtkja2xuGTA5vUZd12nvdEzplap03tWu0YSuN5L2Jpl0klXfBhgn\nKsuVpKMknSxpKMmY0WYqD77v3CX9d++03vOS9pB0NvDGsroHSfpHSUMkfRh4HXBHRKwG7gYul7RP\nOmngsEpTqKvEPQwYlrzUsDT+rvGTm4EZkvaS9A6Scbas9WDXA5+W9Pp0XOoi4NqSz7lW0o8yYnhf\nSUJ8EdiWXpPfksy0+0YawzBJb0/r7UPSQnlR0hjgC1VO8wfApV1TwSUdKOm0atelLL6hkoanb4el\n16z0vC+UdEgax4Wl513mFmC8pL9Oj3EJ0BYRT9UaizU2JyrL2zDgG8AaktbOgcC/VKkfABHxBPBt\n4EGS7sHxwNyyug8BRwLPA18DPljSUjmLZGxsIbAOuJFkgkHXZIoXsiZTSDqUJKE+lsazGXiypMp5\nwF7AH0m61v4+jbfbsSPiLpIuvHtJJj88y66TG8ZWOK8uRwK/kvQi8AAwMyLuS7v83p+WLyNpYX0k\n3Wc6cBywgWRyyk1lxyxdt/YdkhmLd0vaCMyjZAwwnWlZbY3YH4BOku7WO4FNXUkvIn6Qfv5jJBMl\nbouI/6p07HQM8YMkrcF1wJuB3Z3EYQ1IyRh1nT9Uugb4K5Iujm7dLZLeSfI/yKJ0080R8fU6hmgN\nTtKngE9HRM2tpKJRsvC3DTg6IrotLDYbLIbk9LnXAt8laf5nuT8iau5mMBto0pl2noJtg14uXX8R\nMZdkTUY1NS2+NDOzga3IY1THS2qT9L+S3pB3MNZYIuK6Ru72M7OX5dX115OHgeaI2CRpMjCHZMX/\nLiQV7oalZmZWu4josfeskC2qiOhIp/kSEb8A9pR0QEZd//Twc8kll+QeQyP8+Dr5Gvk61fenVnkm\nKpExDlV6vzZJE0lmJ66rV2BmZlYcuXT9SZoFtACjJC0jWcA3FIiIuJrkfmKfJbldy2bgo3nEaWZm\n+cslUUXEx3son8mud8m2V6ClpSXvEBqCr1PPfI1q4+vUt3JZ8NtXJEUjx29mNphJIhp1MoWZmVkX\nJyozMys0JyozMys0JyozMys0JyozMys0JyozMyu0ot7rzzJMOukM2ld2VCwbfcgI7rp/Tp0jMjPr\nX05UDaZ9ZQdnNF1UsWzOcj9b0swGHieqAWTx8hVMOOI9meVucZlZI3KiGkC2bRuW2doCt7jMrDF5\nMoWZmRWaE5WZmRWaE5WZmRWax6gK5oqLL2bDsmWZ5RvWroWmOgZkZpYzJ6qC2bBsGdPGjcss/9G2\nh+sXjJlZAThRFcy9C5aypG1DZvmGP9UxGDOzAnCiKpiNna9iXNPUzPLtO7LLelJtnZXXWJlZUTlR\nDSLV1ll5jZWZFZVn/ZmZWaE5UZmZWaHlkqgkXSOpXdLvqtS5UtLTktokTahnfGZmVhx5taiuBSZl\nFUqaDBweEUcC5wJX1SswMzMrllwSVUTMBdZXqXI6cH1a9yFgP0mj6xGbmZkVS1HHqMYAz5W8X5Fu\nMzOzQabhp6dPmzZt5+uWlhZaWlpyi8XMzLK1trbS2tra6/2KmqhWAGNL3jel27opTVRmZlZc5Y2J\n6dOn17Rfnl1/Sn8quQ04C0DS24ANEdFer8DMzKw4cmlRSZoFtACjJC0DLgGGAhERV0fEHZJOlfQM\n0AmcnUecg8m6tWuZNmVKZvnI5mamzphRv4DMzFK5JKqI+HgNdT5Xj1gssbFzc9Wb4S5esJSpzlNm\nloOijlFZnW3bMZxxI7NveNvmewGaWU6KOj3dzMwMcIsqF5NOOoP2lR0VyxYvX+Mn+JqZlagpUUka\nAnwYOD7dtDewHdgE/A6YFREv9UuEA1D7yo7Mx218c8nuP2/KzGwg6jFRSXoLcCLwy4j4SYXyw4Fz\nJD0aEff1Q4xmZjaI1dKieiki/iOrMCKeBa6UdJikoRGxpe/CMzOzwa7HyRQR8VjXa0mvlTQ8o94i\nJykzM+trvZ3193ngbQCSTpT0jr4PyczM7GW9TVS/BcZJem1E/AZ4dT/EZGZmtlNvE9VYYAtwoaR7\ngDf3fUhmZmYv6+06qkXAzyJilqRRwAf6ISYzM7Odetui+inwxvT1YcDBfRuOmZnZrqq2qCQNA0ZE\nxFqAiNgOPJK+ng/ML6k7NiKeq3ggMzOz3VQ1UUXEnyS9V9I+wJyI2FxeR9JI4CPAQnZ9fLwNIIuX\nr2DCEe+pWDb6kBHcdf+cOkdkZoNFj2NUEfFzSQcDF0g6CBgO7ElyC6VOYDnww4jY2K+RWq62bRuW\nedunOb6zupn1o5omU0TEauDSfo7FzMysm15NppD0GUn3SPqNpHP6KygzM7MuvZ31tzYi3gWcBvxJ\n0pf7ISYzM7Odepuohks6NiLWR8R1wOP9EZSZmVmX3i74PRo4VtK/AgFskfQiMDYiftzn0ZmZ2aDX\n20R1K6CI+LykoSS3UHo78HHAicrMzPpcrxJVRPxfyestwDxgnqSbe3McSacAV5B0PV4TEZeVlb+T\nJCkuSjfdHBENMwf64ouvYNmyDZnlq9Z2+HHzZmY16m2LqqKIWNRzrYSkPYDvAe8GVgLzJd0aEU+W\nVb0/Ik7ri/jq7c5b72UvnZ1Z3tl5Wx2jMTNrbH2SqHppIvB0RCwFkDQbOB0oT1Sqd2B9ZUtnJ6c2\njcwsn79jRx2jMTNrbHkkqjHsequl5STJq9zxktqAFcAXImJhPYKz3lu3di3TpkzJLB/Z3MzUGTPq\nF5CZDSh5JKpaPAw0R8QmSZOBOcBRlSpOmzZt5+uWlhZaWlrqEZ+V2Ni5mSVt2WNyixcsZarzlNmg\n19raSmtra6/3yyNRrQCaS943pdt2ioiOkte/kPR9SQdExLryg5UmKsvHth3DGTdyamZ5m+8FaGZ0\nb0xMnz69pv16u+C3L8wHjpB0aDrF/Uxgl9kFkkaXvJ5IMiW+W5IyM7OBr+4tqojYLulzwN28PD39\nCUnnJsVxNfAhSZ8FtgKbgY/WO04zMyuGXMaoIuJO4M/Ltv2g5PVMYGa94zIzs+LJo+vPzMysZk5U\nZmZWaE5UZmZWaE5UZmZWaEVd8GsDyOLlK5hwxHsqlo0+ZAR33T+nzhGZWSNxorJ+t23bMM5ouqhi\n2RwvBjazHrjrz8zMCs2JyszMCs1df7up2sMR/WBEM7O+40S1m6o9HNEPRjQz6ztOVLup2sMR/WDE\n2q1eu54pU6Zlljc3j2TGjOw7s5vZwOdEZbl6sTNY0jYhs3zhgmudqMwGOScqy9eOHbSMrNwyBZiz\nvLOOwZhZEXnWn5mZFZoTlZmZFZoTlZmZFZrHqKzQfJ9AM3OiskLzfQLNzF1/ZmZWaG5RWcPyYmGz\nwcGJKsOkk86gfWVHZvni5Wt8P7+cebGw2eCQS6KSdApwBUnX4zURcVmFOlcCk4FOYEpEtNUzxvaV\nHZljIwDfXNI4X4BLNrQxbmT2F3rD6mGx8FWPLs9scVVqbbW2ttLS0tKHAQ48vka18XXqW3Ufo5K0\nB/A9YBIwHviYpNeV1ZkMHB4RRwLnAlfVO86BZMmGuub4wuhqcVX6ufPWe7vVb21trX+QDcbXqDa+\nTn0rjxbVRODpiFgKIGk2cDrwZEmd04HrASLiIUn7SRodEe11j9YaV5UW1+UPPtVt2vvqdYuYc8Nc\nwFPfzYokj0Q1Bniu5P1ykuRVrc6KdFu3RHXeeedlftDMmTOrBlJtHMpjUANbpWnvrdv+m5amKQBc\n/uBnM9dvrVizgjEHjsk8drVyJ0Cz3lNE1PcDpQ8CkyLinPT9J4GJEXF+SZ3bgX+LiHnp+18BX4yI\nR8qOVd/gzcysT0WEeqqTR4tqBdBc8r4p3VZeZ2wPdWo6QTMza2x5LPidDxwh6VBJQ4EzgfJH4t4G\nnAUg6W3ABo9PmZkNTnVvUUXEdkmfA+7m5enpT0g6NymOqyPiDkmnSnqGZHp65We+m5nZgFf3MSoz\nM7PeaNh7/Uk6RdKTkp6S9KW84ykiSddIapf0u7xjKSpJTZLukfS4pMcknd/zXoOPpGGSHpK0IL1O\nl+QdU1FJ2kPSI5LKhzQsJWmJpEfT36ff9li/EVtU6aLhp4B3AytJxr3OjIgnq+44yEh6B9ABXB8R\nR+cdTxFJOhg4OCLaJI0AHgZO9+9Sd5L2iohNkl4FPACcHxE9fskMNpIuAI4D9o2I0/KOp4gkLQKO\ni4j1tdRv1BbVzkXDEbEV6Fo0bCUiYi5Q0y/CYBURq7tuzxURHcATJGv2rExEbEpfDiMZ3268v3L7\nmaQm4FTgh3nHUnCiF/mnURNVpUXD/nKxV0TSOGAC8FC+kRRT2qW1AFgN/DIi5ucdUwFdDnwBJ/Ge\nBPBLSfMlfaanyo2aqMz6VNrt9zPgn9KWlZWJiB0R8SaSdY1vlfSGvGMqEknvA9rTFrrSH6vshIg4\nlqT1eV46TJGpURNVLYuGzWoiaQhJkvpxRNyadzxFFxEvAPcCp+QdS8GcAJyWjr/8BDhZ0vU5x1RI\nEbEq/XcNcAvdb6O3i0ZNVLUsGraE/7Lr2Y+AhRHxnbwDKSpJr5a0X/r6z4D3suuNpAe9iPhKRDRH\nxGEk30n3RMRZecdVNJL2SnswkLQ38JfA76vt05CJKiK2A12Lhh8HZkfEE/lGVTySZgHzgKMkLZPk\nhdNlJJ0AfAJ4VzpV9pH0eWm2q9cA90pqIxnDuysi7sg5JmtMo4G56Xjng8DtEXF3tR0acnq6mZkN\nHg3ZojIzs8HDicrMzArNicrMzArNicrMzArNicrMzArNicrMzArNicrMzArNicrMzArNicrMzApt\nSN4BmNmu0gcTfhQ4jORxNhOBb0XE4lwDM8uJW1RmxXMMyd3cF5HcUPhGYFWuEZnlyInKrGAi4pGI\n2AIcD9wXEa0R8VLecZnlxYnKrGAkvUXSKGB8RCyWdGLeMZnlyWNUZsVzCsnj3udJOgN4Pud4zHLl\nx3yYmVmhuevPzMwKzYnKzMwKzYnKzMwKzYnKzMwKzYnKzMwKzYnKzMwKzYnKzMwK7f8Bw3GWxmPs\ncPcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10eb22e48>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "N_samples = 100000\n",
    "N_reject = 10\n",
    "shapes = [1.0, 2.0, 5.0, 10.0]\n",
    "N_shapes = len(shapes)\n",
    "bins = np.linspace(0, 5.0, 50)\n",
    "\n",
    "plt.figure(figsize=(6,8))\n",
    "for i,shape in enumerate(shapes):\n",
    "    scale = shape\n",
    "    print(\"Shape: {}\\t shape {}\".format(shape, scale))\n",
    "    exact_samples = \\\n",
    "        np.array(\n",
    "            [exact_sample_gamma(shape, 1./scale,\n",
    "                                npr.randn(N_reject), npr.rand(N_reject))\n",
    "             for _ in range(N_samples)])\n",
    "\n",
    "    approx_samples = \\\n",
    "        np.array(\n",
    "            [approx_sample_gamma(shape, 1./scale,\n",
    "                                 npr.randn(), npr.rand())\n",
    "             for _ in range(N_samples)])\n",
    "\n",
    "    plt.subplot(N_shapes, 1, i+1)\n",
    "    plt.hist(exact_samples, bins, normed=True, color='r', alpha=0.5, label=\"Exact\")\n",
    "    plt.hist(approx_samples, bins, normed=True, color='b', alpha=0.5, label=\"Approx\")\n",
    "    plt.title(\"shape: %.1f, scale: %.1f\" % (shape, scale))\n",
    "    plt.yticks(np.linspace(0,1.5,4))\n",
    "    plt.ylabel(\"$p(x)$\")\n",
    "    plt.xlabel(\"$x$\")\n",
    "\n",
    "    if i == 0:\n",
    "        plt.legend(loc=\"upper right\")\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##     Simple example of black box variational inference using reparameterization gradients for a gamma-Poisson model. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Simulate data from a gamma-Poisson model\n",
    "N = 10\n",
    "alpha0 = 1.0\n",
    "beta0 = 3.0\n",
    "lmbda_true = npr.gamma(alpha0, 1./beta0)\n",
    "y = npr.poisson(lmbda_true, size=N)\n",
    "assert np.sum(y) > 0, \"Trivial dataset, try larger beta0.\"\n",
    "\n",
    "# Compute the true posterior\n",
    "alpha_star = alpha0 + np.sum(y)\n",
    "beta_star = beta0 + N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Define the log probabilitiy, \n",
    "def logprob(lmbda):\n",
    "    \"\"\"\n",
    "    log p(y, lmbda)\n",
    "    \"\"\"\n",
    "    poisson_lp = - np.sum(sp.gammaln(y)) \\\n",
    "                 + np.sum(y * np.log(lmbda)) \\\n",
    "                 - N * lmbda\n",
    "            \n",
    "    gamma_lp = alpha0 * np.log(beta0) - sp.gammaln(alpha0) + \\\n",
    "               (alpha0-1) * np.log(lmbda) - beta0 * lmbda\n",
    "        \n",
    "    return poisson_lp\n",
    "\n",
    "def gamma_entropy(alpha, beta):\n",
    "    \"\"\"\n",
    "    E[-log q(lmbda; alpha)]\n",
    "      = alpha - log(beta) + log(Gamma(alpha)) + (1-alpha) digamma(alpha)\n",
    "    \"\"\"\n",
    "    return alpha - np.log(beta) + sp.gammaln(alpha) + (1-alpha) * sp.digamma(alpha)\n",
    "\n",
    "# Assume the variational factor is given beta_star\n",
    "def objective(alpha, beta, Zs, Us):\n",
    "    \"\"\"\n",
    "    E[log p(y | lmbda)] where\n",
    "    lmbda ~ q(lmbda; alpha, beta)\n",
    "\n",
    "    We approximate this with samples of uniform and normal streams\n",
    "    which are then converted into a stream of gammas with the given\n",
    "    parameters.\n",
    "    \"\"\"\n",
    "    lmbdas = np.array([exact_sample_gamma(alpha, 1./beta_star, zs, us)\n",
    "                       for zs, us in zip(Zs, Us)])\n",
    "    lps = [logprob(lmbda) for lmbda in lmbdas]\n",
    "    return gamma_entropy(alpha, beta) + np.mean(lps)\n",
    "\n",
    "# Compute the gradient wrt alpha using autograd\n",
    "gradient = grad(objective, argnum=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x111041cf8>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEPCAYAAACzwehFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xm8HGWd7/HPN2RPgGwQyEJCwLAkLAYNRECPOIjDPgwK\nw1xF43WbURQULyoOR2cBR0UdZ9QrYIbxJeAdRLYBBAaPc1mcsGUhbBGiSYiEfQkhh4T85o+nm5wc\nzjm9VXd1d77v16te3VVdXf2rl5HveZ6n6ilFBGZmZgMZlHcBZmbW/BwWZmZWksPCzMxKcliYmVlJ\nDgszMyvJYWFmZiXlGhaSLpG0VtKSPj77vKTNksblUZuZmW2Rd8tiAXBU742SpgBHAn9oeEVmZvYm\nuYZFRNwOPN/HR98Bzm5wOWZm1o+8WxZvIul4YFVELM27FjMzSwbnXUBPkkYAXyZ1Qb2xOadyzMys\noKnCAtgDmA4sliRgCnCvpLkR8VTvnSV5YiszsypEREV/iDdDN5QKCxHxQETsEhEzImJ3YDXw1r6C\noigi6r6cd955Dfkdn0tzno//nfl82ulcIqr7GzvvS2cvA+4EZkpaKekjvXYJ3A1lZpa7XLuhIuK0\nEp/PaFQtZmbWv2bohmp6HR0deZeQmXY6F2iv82mnc4H2Op92Opdqqdr+q2YgKVq5fmsNkqru5zVr\nRoV/0y03wG1mZk3OYWFmZiU5LMzMrCSHhZmZldTyYeFxRzOz+mv5sHj55bwrMDNrfy0fFs89l3cF\nZmbtz2FhZmYltXxYPN/Xo5PMzCxTLR8WblmYmdWfw8LMzEpyWJiZWUktHxYeszAzq7+WDwu3LMzM\n6s9hYWZmJbV8WLgbysys/lo+LNyyMDOrv1zDQtIlktZKWtJj2z9KekjSIkm/kLTDQMdwWJiZ1V/e\nLYsFwFG9tt0MzIqIA4HlwJcGOoDDwsys/nINi4i4HXi+17ZbI2JzYfW3wJSBjvHaa9DdXacCzcwM\nyL9lUcp84MaBdhg3zoPcZmb11rRhIekrwMaIuGyg/caNc1eUmVm9Dc67gL5I+jBwNHBEqX3Xrevk\nW9+C3XaDjo4OOjo66l2emVlL6erqoqurq6ZjKHJ+Lqmk6cB1EbFfYf19wLeBd0bEsyW+G8ccE3zi\nE3DccXUv1bZRksj7/ydmWSr8m1Yl38n70tnLgDuBmZJWSvoI8H1gNHCLpPsk/WCgY7gbysys/nLt\nhoqI0/rYvKCSY4wd67AwM6u3ph3gLpevhjIzq7+2CAu3LMzM6sthYWZmJbV8WIwd624oM7N6a/mw\ncMvCzKz+HBZmZlaSw8LMzErK/Q7uWkiKjRuD4cPT7LODWj76rBn5Dm5rNy13B3cWBg+GUaPgpZfy\nrsTMrH2VvINb0ulAOX9WCSAiLq21qEoVu6LGjGn0L5uZbRvKme7jhYi4ppyDSTqhxnqq4stnzczq\nq2Q3VLlBUem+WfIgt5lZfZU9kaCkwcD7gXmFTaOA14H1wBLgsojYkHmFZXBYmJnVV1lhIentwOHA\nLRFxeR+f7wF8XNLiiPhNxjWW5LAwM6uvclsWGyLiwv4+jIjHgH+SNEPS0Ih4LZvyyuMxCzOz+irr\n0tmIWFrmfo83OijALQszs3qr+D4LSdMkHS/psML6SdmXVRmHhZlZfVXzpLyRwMnAXEnrSI9FvSrT\nqirkbigzs/qqJixOBD4XEc9JGgEckXFNFdtpJ1i7Nu8qzMzaVzXTfayOiOcAIuLVWn5c0iWS1kpa\n0mPbWEk3S3pE0q8k7VjqOFOnwqpVtVRiZmYDqSYsnpZ0uaTjJO0P7FPD7y8Ajuq17Rzg1ojYC7gN\n+FKpg0yenFoWmzbVUImZmfWrqllnJc0ETgeGAhdFxKNVFyBNA66LiP0L6w8D74qItZJ2AboiYu9+\nvhvF+idPht/+NrUyzLLkWWet3VQz62w1YxYUwuEr1Xy3DDtHxNrC7zwpaedyvlTsinJYmJllr6qw\naLAB/6Tr7OwE0hTlN93UwTve0dGAkszMWkdXVxddXV01HaPibihJdwEnR8QTkiYBdwEfI80V9R+V\n3pTXRzfUQ0BHj26oX0dEn+MiPbuhzjoLdt0Vzj67otMxK8ndUNZuGtUNdUyPq6HWSNonItZXcZwi\nFZaia4EPA98gjYuUNZPt1KmwYkUNVZiZWb8qvhqqGBQAkoYAR1b745IuI93UN1PSSkkfAS4AjpT0\nCPCewnpJvnzWzKx+Km5ZSNoeOAb4M2AqcD1l/vXfW0Sc1s9Hf1LpsXbbzWFhZlYv5U5RPpb0LIsT\nSGMTLwNfjYhFdaytIm5ZmJnVT7ndUP9CetDR/IjoIM0NdaqkWfUqrFITJ8ILL8CGXB6/ZGbW3sq6\nGkrSmIh4ode2QcCFwL0R8dM61VeqruhZ/+67wy23wJ575lGNtStfDWXtppqrocp9nsULfWzbHBGf\nA2qaHypL7ooyM6uPkmEhaZik8f19HhFX9tg31/unPchtZlYfJcMiIrqBeZL+ojAl+ZtIGiPp48C0\nrAushFsWZmb1UdbVUBFxfeFu6jMl7QQMJ00i+DqwHlgFXBwRL9at0jJMnQqLF+dZgZlZeyr7PouI\neBL4B0n/G1gOLKz1eRZZmzoVrr8+7yrMzNpPNdN9DAYOAz4jaSQpOO4EpkbEt7IsrlLuhjIzq49q\nHn60LiL+PiJOJt2k9wSpW+qkTCurgsPCzKw+qmlZbJK0ALgOeBjYGBGXSsr9bu5x42DjxjRd+Q47\n5F2NmVn7qGYiwSuA84EDgU8Ctxe25z60LLl1YWZWD+XcZzGn97aIeDQi/iYizoiIuwfat9EcFmZm\n2SunG+poSbPZ+pkTfc19IGAKcF8WhVXLYWFmlr2SYRERfyfpxIi4uhEF1cphYWaWvXLHLM6X1Cnp\naEkTen4gqamGkj3lh5lZ9sq9Guo6YCXpmRYXS3oVuBtYCOwJ/FV9yqvc1KmwcmXeVZiZtZdyw+Kr\nhTmifiKp+FzstwJzgLn1Kq4a7oYyM8teuWHxAaD4zIooTFn+a+DXkpbVpbIqTZ0Kq1dDRLqU1szM\nalfumMU3JP1fSR8Edu712ZKMawJA0pmSHpC0RNLPJA0t53ujR8OwYfDMM/Woysxs21RuWHwbuBzY\nFThE0iJJN0v6RuGzTEmaBHwGmBMR+5NaQKeW+/2ZM+GRR7Kuysxs21XuFOXFQOgqbpM0ETiYdCd3\nPWwHjJK0GRgJrCn3i7NmwbJlcNhhdarMzGwbU83cUABExFrgWkmZd/hExBpJ3yZdgbUeuDkibi33\n+8WwMDOzbFQdFkURcWcWhfQkaQxpRttpwIvAlZJOi4jLeu/b2dn5xvuOjg46OjqYNQtuuCHrqszM\nWlNXVxddXV01HUMRfc3ckS9JJwNHRcTHCusfBA6OiE/32i/6qn/1ajjoIFi7tiHlWpuTRDP+/8Ss\nWoV/0xVdL1rN8yyKPzah9F5VW0kaSB8uScB7gIfK/fLkydDd7SuizMyyUnVYAB/KrIpeImIhcCVw\nP7CYNEnhj8v9vgT77utxCzOzrNQSFnW95S0ivhYR+0TE/hFxekRsrOT7HuQ2M8tO2QPckoaQ7rOA\nFBRjJO3WY5fVEbE5y+Jq4bAwM8tOJVdDjQB277E+HphOCo4AngZezayyGs2eDb/8Zd5VmJm1h7LD\nIiJeAn5TXJf0toj4r7pUlYFiy8JzRJmZ1a6WMYumtssusHkzPPVU3pWYmbW+WsLimsyqqAPJ4xZm\nZlmpOiwi4ndZFlIPs2bBAw/kXYWZWetr224oSIPcblmYmdWu5AB34cl45cx1IICIuLTWorIyaxZc\ncUXeVZiZtb5yroZ6ISLKGp+QdEKN9WTKV0SZmWWjKScSLFd/Ewn2tNNOsGhRmi/KrBqeSNDaTTUT\nCVZyB/dg4P3AvMKmUcDrpOdNLAEui4gNlfx4IxRbFw4LM7PqlRUWkt4OHA7cEhGX9/H5HsDHJS2O\niN+86QA5mj07XRH13vfmXYmZWesq92qoDRFxYUQs7evDiHgsIv4JWCVpaHbl1e5tb4OFC/Ouwsys\ntZUVFv2FhKQhPQe1I+LxiHgtq+KycNhhcPvtaZDbzMyqU/F9FpK2l3SqpJ+T5oqalX1Z2dljD9i0\nCVauzLsSM7PWVe6YxVjS4PYJpIHtl4GvRsSiOtaWCQkOPTS1LqZNy7saM7PWVG7L4l9IVz7Nj4gO\n4GTgVElN3aooKnZFmZlZdcq6z0LSmIh4ode2QcCFwL0R8dM61VeqrpL3WQDcfTfMnw9L+xx5MRuY\n77OwdlPNfRY135Qn6eSIuLKmg1T/22WFxcaNMG5cGrcYO7YBhVlbcVhYu6kmLGqeSLBeQSFpR0n/\nLukhScskHVztsYYMgblz4a67sqzQzGzb0cyzzn4PuCEi9gEOAB6q5WDFQW4zM6tcU4aFpB2AwyNi\nAUBEbCo81rVqhx0Gd9yRSXlmZtucmsJC0jRJj0rqkHSMpGEZ1bU78IykBZLuk/RjSSNqOeAhh8C9\n90J3d0YVmpltQ8qeSLAvEfEHSR0RsSarggoGA3OAv46IeyR9FzgHOK/3jp2dnW+87+jooKOjo88D\n7rADzJwJ990H8+b1uYuZWVvq6uqiq6urpmNUfDWUpCuAdcCdwB0R8UhNFfT9GxOBuyJiRmH9MOD/\nRMRxvfYr62qoojPOgKlT4eyzMy3X2pyvhrJ205CroSLiVOAC4DXgDEkPS/pGlhMIRsRa0qSEMwub\n3gM8WOtxDz3U4xZmZtWopmVxcOF7vy2svx9YDBwbERdmVph0AHAxMAR4HPhIRLzYa5+KWhZr1qQp\ny9euTZfTmpXDLQtrN3V9+FEPRwIbJZ0JvAKsBJ4GlldxrH5FxGLg7Vkec9IkmDEjtS76GdowM7M+\nlOyGkjRM0vgem34J/CYiTomI+RHRCewJvCJpap3qzMwJJ8A1ZT1R3MzMisqdG+pYYHvg6oh4tY/P\ni7PSPhgRDbv1rdJuKIAlS+DEE+Gxx9KMtGaluBvK2k1d54aStAswH9gZGE4aS3id1BW1Gri495hC\nvVUTFhGpK+raa2G//epUmLUVh4W1m1wmEsxTNWEB8LnPwYQJcO65dSjK2o7DwtpNQy6dlfQxSbdJ\n+v+SPl7p95uBxy3MzCpTzXQfz0bEEcDxQLekczKuqe4OPxwefxyeeCLvSszMWkM1YTFc0pyIeD4i\nLgWWZV1UvQ0eDEcfncYtzMystGrCYn/gNEk3SroB+GhhIsEPZlxbXR1/vMPCzKxc1dzBPa/wvTsL\nU3y8DXgHcFpEzKlDjQPVUtUAN8DLL8PkyakravvtMy7M2ooHuK3dNOpJeYdHxJ0AEfEasIrU2rio\nimPlZvvt4Z3vhKuvzrsSM7PmV01YjJd0vaR9C+tnAX8HvOlmvWY3fz5cfHHeVZiZNb9qwmJhRBxL\nmuIDYBrwO9L8UC3luOPgkUfSYmZm/asmLOZI+itgH0l7kMJiBDAq08oaYMgQOP10uOSSvCsxM2tu\n1Qxw7wwcAiwBZpAunZ0P3BMRt2Re4cC1VD3AXbR8eXo+96pVMDSzJ3JYO/EAt7WbRg1wPw1MAL4I\nzIyItRFxfqODIitveQvss48vozUzG0g1YfF50mD2L4BXJZ2VbUmN97GPeaDbzGwg1XRDnRAR1/RY\n/0BE/L/MKyuvlpq7oQBefRWmTIF774Xp02uvy9qLu6Gs3TSqG2pPSXMlzZD0LmCPKo7RVEaMgNNO\n80C3mVl/qmlZjAK+AMwFlgLnRUR3HWorp5ZMWhYADz+cbtJ77DHf0W1bc8vC2k0uz7OQdFRE/Kqm\ng/R/7EHAPcDqiDi+j88zCwuAU06Bgw6CL34xs0NaG3BYWLupW1gU7qs4nfRUvK0+AvaOiF0r+dFy\nSToTOAjYoRFhsXQpHHlkmr585MjMDmstzmFh7aaeYxaPAodFxBG9lncDH6m40jJImgIcDTTsOqX9\n9oN58+Cilprlysys/kqGReHZFbdGxMa+Po+Im3rum2Ft3wHOBhr6J92558I3vwkbNjTyV83Mmtvg\nMvY5WtJsUpdTX4r/MRcwBbiv1qIkHQOsjYhFkjoG+G06OzvfeN/R0UFHR0dNv33QQXDAAbBgAXzq\nUzUdysysKXR1ddHV1VXTMWoe4K4HSf8A/C9gE2neqe2BqyLiQ732y3TMouiuu+DUU+HRR2HYsMwP\nby3GYxbWbnK5GqreCvdyfL4RA9w9nXACHHIIfOlLdTm8tRCHhbWbasKinG6ovn5oIukqpdHAcGAj\n8BLwq4jYVM0xm813vgNz58IHP5ju7jYz25ZV1bKQtEdEPNbH9rdExPJMKiuvjrq1LAD+5m9SV9QV\nV9TtJ6wFuGVh7abh3VCSJkTEM5KG9He1VD3VOyzWr4d9902D3e9+d91+xpqcw8LaTaPmhuppZ0nf\nIs1E23ZGjoRvfxvOOAM2NjwKzcyaR61hsR/wr8B/115KczrpJNhllzSGYWa2rao1LER6Ut6LGdTS\nlKR0R/e3vgWLFuVdjZlZPmoNizuAfwIOyKCWpjV9Olx4YZrGfP36vKsxM2u8pr/PYiD1HuDuKSKF\nxfjx8M//3JCftCbhAW5rNw0f4Ja0X+H1QEn71nKsZifBD38I110HN9yQdzVmZo1VazfUqMLrEGC5\npB1qPF5TGzMG/u3fYP58WLEi72rMzBqn1rBYKOl9wEkRsTEiXsqiqGb2rnfBV74Cxx0HL7X92ZqZ\nJdXewX0g8H5gF6ALuCwiXs+2tLLqaNiYRU8RaUba1avhmmtgu+0aXoI1kMcsrN00csxiNvC3EfFR\n4FVgnyqP05Ik+P7305VR55yTdzVmZvWXydVQkk6OiCszqKfS382lZVH07LPpyXqf/nS6y9vak1sW\n1m4aOevs8Ih441lyxaCQNCwiuqs5ZisaPx5uuSWNYwwbBp/4RN4VmZnVR1VhAUySNAtYBDwB7ADM\nAV4G7s6otpYwbRrcemuaaHD4cDj99LwrMjPLXlVhERGPS1oJHAwcCjwP3B4Rr2RZXKvYc8/Uwjji\nCBg8GP7yL/OuyMwsW9W2LIiITYXnYwfp8ac7Az/NqK6Ws/fecPPN8Kd/Ck89BWeemXdFZmbZqTos\nCh4Dfk5qXcyqvZzWNns23HEHHHUU/PGPcMEFMKjWO1nMzJpArQ8/mg58FlgeET/IqKZKfj/Xq6H6\n8+yz6aa9GTPg4ovTWIa1Ll8NZe2mbvdZSBojaU9Jh0jatbBtBmnM4so8gqKZjR+fBr1few0OPxxW\nrsy7IjOz2pTbSbIY2Bd4PCL+CGmQG7gO+Pusi5I0RdJtkpZJWiqp5e5iGDkSfv5zOOUUOPhguO22\nvCsyM6teuWHxw4i4NiKeknSapAslnRYR64F63Iy3CTgrImYB84C/lrR3HX6nriT4whfgZz9LV0h9\n/euwaVPeVZmZVa7csHiq+CYiLgO6C68AmT8OKCKejIhFhffrgIeAyVn/TqMccQTccw/cfjscdhgs\nX553RWZmlSk3LA6U9M7iAgzp8b6uT8krDKIfSIs/53vyZLjppvQApXnz4Ec/gs2b867KzKw8ZV0N\nJWkJcA/pmdu9HRQR+2ddWOF3R5Nmtf3biLimj8/jvPPOe2O9o6ODjo6OepSSqQcfhI9+NF1W+6Mf\nwX775V2RDcRXQ1mr6+rqoqur6431r33taxVfDVVuWPxJRNxa6We1kDQYuB64MSK+188+TXnpbDk2\nb4Yf/xi++tX0MKVzz4Xtt8+7KuuLw8LaTd0unR0oDOoRFAU/AR7sLyha3aBB8MlPwtKlsGYN7LUX\nXHQRvN7wp4KYmZWWyRTlWZN0KPBfwFLSdCIBfDkibuq1X8u2LHq75x74/Ofhuefg/PPhmGPS1VSW\nP7csrN1U07JoyrAoVzuFBaQn8F1zTeqaGjECOjvTXFMOjXw5LKzdNDwsJE2IiGckDYmIjVUfqPrf\nb6uwKNq8GX7xC/ja19LNfV/8IvzZn/nxrXlxWFi7ySMs9gXmA89ExAVVH6j632/LsCjavBmuvhq+\n+U14+mk466z0vIxRo/KubNvisLB208hncBftB/wrLX4PRLMaNAhOOgnuugsuvTTNN7XbbvDZz8Ij\nj+RdnZltS2oNCwHLgBczqMUGcOihcNVVcP/9MHp0epTrEUfAT38K6zO/h97MbGu1dkNNBb4I3BcR\nCzKrqvzfb+tuqIF0d8N118GCBanl8ed/nuafOvxwj21kzd1Q1m7qOmYh6V+B8yOiaTpAtuWw6OmJ\nJ9JkhZdfnp7Sd8op8IEPwNy5fvhSFhwW1m7qHRY3AucDs4GbI+J3lZeYLYfFmz30EFxxBVx5Jbz4\nYrqK6sQTU4tj6NC8q2tNDgtrN/UOixsi4ujC+/cCbwH+IyJ+X2mhWXFYDOzhh+GXv0zL8uVw5JFw\n7LHpsa8TJ+ZdXetwWFi7aVhY9Nj298D4iPhkJT+aFYdF+Z58Em68Ea6/Pj2Iado0eO97U4Acemi6\nn8P65rCwdlP3bqiI+NPC+wNJT8ibCHw1Im6stNgsOCyqs2kTLFwIv/oV/Od/wqJFMGcOvPvdqbtq\n3jzfy9GTw8LaTd1bFsAXgL8FpgNf72va8EZyWGRj3Tq44w749a/TA5oWLYJ9900tjnnz0jJ1at5V\n5sdhYe2m3mHxMrAc+LuIuKqK+jLnsKiPV1+Fu++GO+9Ml+XedRcMGQJvf3ta5s5NLZHx4/OutDEc\nFtZu6h0Wfx4Rv6iqsjpxWDRGBKxYkQKkuNx/P4wbl0LjwAPhgAPSMm1a+0186LCwduNZZ61hNm+G\nxx6De++FxYu3LOvWwezZ6el/s2fDrFmpS2vnnVs3RBwW1m4cFpa7Z56BZcvSQ50eeCA9QnbZstQ6\n2Wcf2HvvtOy1F8ycCTNmNP/9Hw4LazcOC2tKEenO8ocfThMgPvxwWpYvh1WrYNIk2HPPLcsee6QQ\n2X33NA9W3hwW1m4cFtZyNm5M4yG/+13q1iq+rliRltGjYfr0Lcu0aWnZbbe07Lhj/bu3HBbWbhwW\n1lYiYO1a+P3v07JiBaxcCX/4w5ZXSJf1FpfJk2HKlPQ6aVJ6HT++tjmyHBbWbhwWtk2JSPNfrVqV\nlieegNWrtyx//GPatm5dmt5k0iTYdde07LJLep04cetlxIg3/47DwtpNNWExuF7F1ErS+4Dvkp65\ncUlEfCPnkqzJSDBmTFr226///TZsSNOdrFmTlrVr0/rChSlQ1q7dsgwblq7cmjgRdtopvQf47nfT\n+oQJWy8jR7buVV5mlWjKloWkQcCjwHuANcDdwKkR8XCv/dyysMxEwEsvbQmOp59Oyyc/KT7zmeCZ\nZ9LVXk8/Dc8+m14jUjdX72XcuPQ6dmx6P3bs1svo0Q4Zy0/bdENJOgQ4r8dcVOcA0bt14bCwRhio\nG2r9+hQgzz67ZXn+eXjuua3fF5cXXkjbXnttS6tozJgUIMX3O+448LLDDmlp9kuOrXm1UzfUZGBV\nj/XVwNycajHr18iRW67MqkR3dxpvKYbHCy9sWS++f/LJ9FpcXnpp69ftttsSHNtvv+W1r2X06K3f\n915GjXL42MCaNSzK1tHRQUdHBwC///3vmT59Op2dnQB+9Wsmr0VZHnfYMPjBD0rvN2FC39sj4Nxz\nO3ntNfjEJzp58UX4/vc76e6GY4/t5OWX4aqr0vpBB3Wybh3ceWfaf9KkTl55BVasSOvbbZf2j+hk\n6FCYMKGT0aPhxRc7GTIEZs3qZNQoWL48rXd0pPU77kjrH/hAJyNHwtVXp/VPfSqtX3RRWj/33LR+\n/vmdSPn/77ktvnZ1db2xXvzvZaWauRuqMyLeV1h3N5TlZlu4GioitXZeeWXLsm7d1uvFZf36La/9\nvX/11S3biu83bkxXm40YkVpkxff9LcOHv/m19/v+lmHDtrwfOtSPF+6tncYstgMeIQ1w/xFYCPxF\nRDzUaz+HhdXdthAWjfD66yk4iuFRfF9c37AhLT23F5fu7q3fF/crfqe4FD/r+b67OwXGsGFbQqT4\nvr+luE/xez2/33t7celre8/3Q4Zsvf/QoTB4cD4XOrTNmEVEvC7p08DNbLl09qESXzOzJrbddlvG\nSBpp8+Z0QUF399YBUu7S87uvvJIuVOjuTi2l4md9vRbf99xv48at119/ve8QGTJky/a+Xnu/r2Sp\ndmyqKcMCICJuAvbKuw4za22DBm3pkmo2mzdvCY++QqX3Z8X3PbeXWjZsePO2ajRlN1S53A1ljeBu\nKGs31XRDedjHzMxKcliYmVlJDgszMyvJYWFmZiU5LMzMrCSHhZmZleSwMDOzkhwWZmZWksPCzMxK\ncliYmVlJDgszMyvJYWFmZiU5LMzMrCSHhZmZleSwMDOzkhwWZmZWksPCzMxKarqwkPSPkh6StEjS\nLyTtkHdNZmbbuqYLC+BmYFZEHAgsB76Ucz10dXXlXUJm2ulcoL3Op53OBdrrfNrpXKrVdGEREbdG\nxObC6m+BKXnWA+31D6WdzgXa63za6Vygvc6nnc6lWk0XFr3MB27Muwgzs23d4Dx+VNItwMSem4AA\nvhIR1xX2+QqwMSIuy6FEMzPrQRGRdw1vIunDwMeAIyKie4D9mq94M7MWEBGqZP9cWhYDkfQ+4Gzg\nnQMFBVR+smZmVp2ma1lIWg4MBZ4tbPptRPxVjiWZmW3zmi4szMys+TT71VC5kTRF0m2SlklaKumM\nvGuqlaRBku6TdG3etdRK0o6S/r1wA+cySQfnXVMtJJ0p6QFJSyT9TNLQvGsql6RLJK2VtKTHtrGS\nbpb0iKRfSdoxzxor0c/5tOTNwn2dS4/PPi9ps6Rx5RzLYdG/TcBZETELmAf8taS9c66pVp8FHsy7\niIx8D7ghIvYBDgAeyrmeqkmaBHwGmBMR+5PGEk/Nt6qKLACO6rXtHODWiNgLuI0muLm2An2dT9Pd\nLFymvs4FSVOAI4E/lHsgh0U/IuLJiFhUeL+O9B+jyflWVb3CP46jgYvzrqVWhb/qDo+IBQARsSki\nXsq5rFptB4ySNBgYCazJuZ6yRcTtwPO9Np8AXFp4fylwYkOLqkFf59OMNwuXo5//bQC+Q7qQqGwO\nizJImg6lsvyaAAACqElEQVQcCPx3vpXUpPiPox0GqXYHnpG0oNCt9mNJI/IuqloRsQb4NrASeAJ4\nISJuzbeqmu0cEWsh/eEF7JxzPVlq6ZuFJR0PrIqIpZV8z2FRgqTRwJXAZwstjJYj6RhgbaGlpMLS\nygYDc4B/iYg5wHpSt0dLkjSG9Jf4NGASMFrSaflWlbl2+COl5W8WLvxR9WXgvJ6by/muw2IAhS6B\nK4GfRsQ1eddTg0OB4yU9DlwOvFvSv+VcUy1Wk/4yuqewfiUpPFrVnwCPR8RzEfE6cBXwjpxrqtVa\nSRMBJO0CPJVzPTUr3Cx8NNDKQb4HMB1YLGkFqTvtXkklW34Oi4H9BHgwIr6XdyG1iIgvR8RuETGD\nNHB6W0R8KO+6qlXo3lglaWZh03to7YH7lcAhkoZLEul8Wm3AvneL9Vrgw4X3pwOt9sfWVufT42bh\n40vdLNyE3jiXiHggInaJiBkRsTvpD6+3RkTJMHdY9EPSocBfAkdIur/QN/6+vOuyN5wB/EzSItLV\nUP+Qcz1Vi4iFpNbR/cBi0v+xf5xrURWQdBlwJzBT0kpJHwEuAI6U9Agp/C7Is8ZK9HM+3wdGA7cU\n/lvwg1yLLFM/59JTUGY3lG/KMzOzktyyMDOzkhwWZmZWksPCzMxKcliYmVlJDgszMyvJYWFmZiU5\nLMzMrCSHhZmZleSwMDOzkgbnXYBZOylMCbMTaYK2q4H1EVH2A2bMmpXDwiwjhYkNT4+Ivyg8qvK7\npDmfHBbW8twNZZad04HLACLiOWAu8FyuFZllxGFhlp2hFFoRhYfMrCs81tKs5XnWWbOMFLqhjmfL\nsyiOAroi4qr8qjLLhsPCzMxKcjeUmZmV5LAwM7OSHBZmZlaSw8LMzEpyWJiZWUkOCzMzK8lhYWZm\nJTkszMyspP8B/0RHsmjNa/AAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1110c2b70>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Sample randomness from which to evaluate the gradient\n",
    "N_samples = 10\n",
    "N_reject = 10\n",
    "Zs = npr.randn(N_samples, N_reject)\n",
    "Us = npr.rand(N_samples, N_reject)\n",
    "\n",
    "# Lay down a grid of parameters at which to evaluate the gradient\n",
    "alpha_max = 2 * np.ceil(alpha_star)\n",
    "alphas = np.linspace(1.0, alpha_max, 100)\n",
    "\n",
    "# Plot the gradient as a function of alpha\n",
    "# assuming that the model is given beta_star\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot(alphas, [gradient(a, beta_star, Zs, Us) for a in alphas])\n",
    "ax.plot(alphas, np.zeros_like(alphas), ':k')\n",
    "\n",
    "# Plot the true alpha* (since it is known in this model)\n",
    "ylim = ax.get_ylim()\n",
    "ax.plot(alpha_star * np.ones(2), ylim, '-k')\n",
    "\n",
    "# Limits etc\n",
    "ax.set_ylim(ylim)\n",
    "ax.set_xlim(1., alpha_max)\n",
    "ax.set_xlabel(\"$\\\\alpha$\")\n",
    "ax.set_ylabel(\"$\\\\nabla_{\\\\alpha} \\\\; \\\\mathrm{E}_{q(\\\\lambda; \\\\alpha)}\"\n",
    "              \"[\\\\log p(y, \\\\lambda)] + H[q(\\\\lambda; \\\\alpha)]$\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [Root]",
   "language": "python",
   "name": "Python [Root]"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}