Lecture08_02 - Rejection sampling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rejection Sampling\n",
    "\n",
    "##### Keywords: sampling, rejection sampling, rejection sampling on steroids"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Setup Finished\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "print(\"Setup Finished\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic Rejection Sampling\n",
    "\n",
    "The basic idea, come up with by von Neumann is:\n",
    "\n",
    "If you have a function you are trying to sample from, whose functional form is well known, basically accept the sample by generating a uniform random number at any $x$ and accepting it if the value is below the value of the function at that $x$.\n",
    "\n",
    "This is illustrated in the diagram below:\n",
    "\n",
    "![](../images/Rejection.png)\n",
    "\n",
    "\n",
    "### The process\n",
    "\n",
    "1. Draw $x$ uniformly from $[x_{min},\\, x_{max}]$\n",
    "2. Draw $y$ uniformly from [0,$y_{max}$]\n",
    "3. if $y$ < f($x$), accept the sample\n",
    "4. otherwise reject it\n",
    "5. repeat\n",
    "\n",
    "This works  as more samples will be accepted in the regions of $x$-space where the function $f$ is higher: indeed they will be accepted in the ratio of the height of the function at any given $x$ to $y_{max}$.\n",
    "\n",
    "The reason this all works is the frequentist interpretation of probability in each $x$ sliver\n",
    "As we have more samples the accept-to-total ratio reflects the probablity mass in that sliver better.\n",
    "\n",
    "### Example\n",
    "\n",
    "The following code produces samples that follow the distribution $P(x)=e^{-x}$ \n",
    "  for $x=[0,10]$ and generates a histogram of the sampled distribution. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Count 99359 Accepted 10000\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAFJCAYAAABU5W56AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4VOXB/vHvbNknG4TIFvawI5uIIiAo4loVUUz8xQXr\ngihClUJRAYtVqK+8Wiri3jYtIIqtWHm1FQVEEDEskR3CvgcSSCaQyWRmfn8EIiiELDM5M5n7c11c\nycycmbnPU6/eec48c47J6/V6ERERkYBkNjqAiIiIXJiKWkREJICpqEVERAKYilpERCSAqahFREQC\nmIpaREQkgFmNDnA+ubmFPn29hIQo8vNP+vQ1Q43GsOY0hjWnMfQNjWPN+XoMk5LsF3wsJGbUVqvF\n6AhBT2NYcxrDmtMY+obGseZqcwxDoqhFRESClYpaREQkgKmoRUREApiKWkREJICpqEVERAKYilpE\nRCSAqahFREQCWECe8ERERKSqVq/+gYkTf0fz5i0wmUw4nU6uu+56hg69m//85/8IDw+nf/+Bv3he\nXt4x/vKXd/jNb8YZkPriVNQiIlJn9OjRk+effwmAkpIS0tPvoH//gXz++WdMn/7n8z4nMbEeUVHR\nrFmTRbduPWozbqWoqEVExKeiJz9L+Kf/8ulrOm+5jaLJL1TpOSdPnsRsNrN48VdcdllvAObPn0d2\n9lqef/5FXnhhEh06dGLIkDsZNOh63n33TRW1iIiIP2Vl/cDjjz+M2WzGarUyZsxYvvji/7jpplsA\nuOOOu/jhh5X84Q+TcblcDBlyJwDNm7cgO3utkdEvqM4XddS0P0BcNDw62ugoIiIhoWjyC1We/frK\n2Ye+z5g3bw4JCfXKb99zz/08+ugDvPvu38vvs1gsWK1WPB4PZnNgrbOu80UdvvBTHDt2MTy/M5hM\nlX7ee+N/ueBARESCT0JCIg5H2VUZXS4Xf/rTK4wdO4FXXpnK66+/jc1mw+v1YrFYAq6kIQS+nuVu\n3pIY50niTp0wOoqIiBigW7cebNy4HoA33vgTV155FbfeOoTLL7+CWbNmAJCTs51OnTobGfOC6vyM\n2t2qNQCN8w9wIire4DQiIuIv3bv3pHv3nr+4/5prBvG73z1Nevq9jBr1VPn9Dz74SPnv//nPQoYM\nGVYrOauq7s+oTxd1o/wDBicREREjREVFc/31N7F48aLzPn7s2FGKioq49NKutZysckJqRi0iIqHp\nhhtuvuBj9erVZ+zYCbWYpmrq/Iy6tKVm1CIiErzqfFF7k5IoCovSjFpERIJSnS9qTCYOJDSi4fGD\nmD1uo9OIiIhUSd0vamB/QiPC3C7qFx41OoqIiEiV1PnFZFBW1FC2oOxIXLLBaURE6r7hU7/y6etV\n9iRUmZl/4YcfvsftLsVkMjFy5GjatWvvsxyrV//AJ5/MZ+bM81/gwx9CoqgPnC7qRvkHWNO8m8Fp\nRETEH3bu3MG33y7ljTfexWQysW3bFl54YTJ//esco6PVSEgU9dkzahERqZtiYmI4fPgQn332CZdf\nfiVt2rTl7bf/ypo1Wbz//tt4PB5OnTrFpEkvYLPZmDjxdyQnJ3Pw4EGuueY6du7MYevWLVx55VU8\n8shIHn/8YZo1a87u3bsAeP75F895v6+++pIPPvgHZrOZLl26MmLEE2Rnr+XPf34Vq9VKREQEL7ww\njaio6BrtV0gU9YH4n2bUIiJSNyUlNWDq1OnMn/8B7733NhERETz88GPk5eUxceIU6tdP4m9/e4+v\nv/6S6667gYMH9/O///s6Tmcxd955K//610LCwyMYOvQWHnlkJACdOnVh7NgJfPzxh2Rmvk+/fgMA\nOH78OO+99ybvvJNJREQEU6Y8x6pV3/H99ysZOPBa7rornWXLllJQUKiiroxT4VHkRSfQ+LiKWkSk\nrtq3by/R0dFMmDAJgM2bN/L006MYOXI0r776MpGRUeTmHqFz50sBaNiwMTExMdhsNhITE4mNjQPA\ndNYFnHr0uAyAzp27sGzZkvL79+zZw/Hj+Tz99Cig7NrX+/fvIyPjAf72t/d48skRJCU1oEOHTjXe\nr5BY9Q1ls+oGJ45gLXUZHUVERPwgJ2cb06f/EZer7P/nmzZNISbGzmuvvcKECZN45pnJ1K+fVL69\nqRJXVNyyZRMA2dnraNGiZfn9TZo0oUGDZF59dSZ//vNbDB06jI4dO/Of/yzkxhtvZsaMN2nRoiUL\nFnxc4/0KiRk1lH1O3Wn/BhqeOMTeek2NjiMiIj7Wv/9Adu3aya9/fS9RUZF4PF4ee+xJ1q1bzWOP\nPURkZAQJCfU4ejS30q+5cOG/+eCD2URERPDcc78nJ2c7AImJiQwbdg+PP/4wbrebhg0bMXDgIEpK\nXEyd+gKRkZGYTCZ++9tnarxfJq/X663xq/hYbm6hT19v+NSvuH3VPxn+zV/5w6/G813r3hd9jq5H\nfa6kJLvP/3cJNRrDmtMY+obGsXIef/xhxo6dQLNmzX/xmK/HMCnJfsHHQufQd4IWlImISPAJqUPf\noK9oiYhI5fz5z28ZHQEIoRn1obhLcJvMmlGLiEhQCZmiLrXaOBKbpBm1iIgElZApaij7ilbCyeNE\nOYuMjiIiIlIpIVXU+xMbA9A4b7/BSURERConpIp6b2ITAJrm7TM4iYiISOWEVlGfPtGJilpERIJF\nhV/PcrlcTJgwgf3791NSUsKIESNo2LAhjzzyCM2bNwcgLS2NG2+8kXnz5jF37lysVisjRoxgwIAB\nFBcXM3bsWI4dO0Z0dDTTpk0jMTGxNvbrvMpn1Mf2GpZBRESkKios6gULFhAfH8/LL7/M8ePHue22\n2xg5ciQPPPAAw4cPL98uNzeXzMxM5s+fj9PpJD09nT59+jBnzhxSU1N54okn+Oyzz5g5cybPPvus\n33fqQgqi4iiIsGtGLSIiQaPCQ9/XX389Tz75JABerxeLxcL69etZvHgx99xzDxMmTMDhcJCdnU23\nbt0ICwvDbreTkpLC5s2bycrKom/fvgD069ePFStW+H+PLmJvvSYknziMrbTE6CgiIiIXVeGMOjq6\n7BqaDoeDUaNGMXr0aEpKSrjzzjvp1KkTb7zxBq+//jrt2rXDbref8zyHw4HD4Si/Pzo6msLCyp0X\nNSEhCqvVUt19qtDexCZ03L+JRvkH2J3U/ILbVXTe1VClMak5jWHNaQx9Q+NYc7U1hhc9hejBgwcZ\nOXIk6enp3HLLLRQUFBAbGwvAoEGDmDJlCj179qSo6KfvJhcVFWG324mJiSm/v6ioqPx5F5Off7I6\n+1Ip+xJ/WlBWUVHrhPXn0kn8a05jWHMaQ9/QONZcwFyU4+jRowwfPpyxY8cydOhQAB588EGys7MB\nWLFiBR07dqRLly5kZWXhdDopLCwkJyeH1NRUunfvzpIlZRfaXrp0KT169PDVPlXbntMrv5voc2oR\nEQkCFc6oZ82aRUFBATNnzmTmzJkAjB8/nhdffBGbzUb9+vWZMmUKMTExZGRkkJ6ejtfrZcyYMYSH\nh5OWlsa4ceNIS0vDZrPxyiuv1MpOVWSfvkstIiJBJGSuR13O62Xen9M4FHcJo+599YLP0fWoz6VD\nZTWnMaw5jaFvaBxrLmAOfddJJhP7EhvTOH8/Zo/b6DQiIiIVCr2ipuzwd5jbRfKJI0ZHERERqVBI\nFvVPC8p0hjIREQlsIVnUWlAmIiLBIiSLem+iLs4hIiLBISSL+mD8JbjMVpocU1GLiEhgC8mi9pgt\nHExoWDajDrxvp4mIiJQLyaKGsnN+R5ecJLEo3+goIiIiFxTCRa2V3yIiEvhCtqjPfEWr2VEVtYiI\nBK6QLerd9VMASDm22+AkIiIiFxayRX0goRGlZgvNj+4xOoqIiMgFhWxRl1ps7EtoTMrR3Zi8HqPj\niIiInFfIFjXA7vrNiHIVk1SQa3QUERGR8wrxoi77nLrZUX1OLSIigSmki3pX/WYANDumz6lFRCQw\nhXRR7zld1M1zNaMWEZHAFNJFfSQ2iZO2CB36FhGRgBXSRe01mdlTL4Um+fuxul1GxxEREfmFkC5q\nKFv5bfW4aZR/wOgoIiIiv6CiLl/5rQVlIiISeFTUZxaU6XNqEREJQCrq8nN+a0YtIiKBJ+SL+kRU\nPPlRcZpRi4hIQAr5ogbYUy+FS04cJrLklNFRREREzqGi5qczlDU9pmtTi4hIYFFRc9YZynT4W0RE\nAoyKmrPO+a2iFhGRAKOipmzltwcTLXJ3Gh1FRETkHCpqwGmL4EBCo7Ki9nqNjiMiIlJORX3ajgYt\niHGepEHBEaOjiIiIlFNRn7YzqQUALXX4W0REAoiK+rScBi0BaHlERS0iIoFDRX3amRm1FpSJiEgg\nUVGfdjw6nrzoBM2oRUQkoKioz7IzqQUNCnOJOVVodBQRERFARX2OnAZaUCYiIoFFRX0WfU4tIiKB\nRkV9lh1nZtT6nFpERAKEivosB+MbcsoWoRm1iIgEDBX1WbwmMzuTmpNybC8UFxsdR0REREX9czuT\nWmDxerBu2WR0FBERERX1z+04fYYy6/ofDU4iIiKiov6FHadXflvXZxucREREBKwVPehyuZgwYQL7\n9++npKSEESNG0Lp1a8aPH4/JZKJNmzZMmjQJs9nMvHnzmDt3LlarlREjRjBgwACKi4sZO3Ysx44d\nIzo6mmnTppGYmFhb+1Yte+o1xW0yY/1RRS0iIsarcEa9YMEC4uPjmT17Nu+88w5TpkzhpZdeYvTo\n0cyePRuv18uiRYvIzc0lMzOTuXPn8u677zJ9+nRKSkqYM2cOqampzJ49m9tuu42ZM2fW1n5VW4kt\nnL2JTcoOfbvdRscREZEQV2FRX3/99Tz55JMAeL1eLBYLGzZsoFevXgD069eP5cuXk52dTbdu3QgL\nC8Nut5OSksLmzZvJysqib9++5duuWLHCz7vjG9suaYPpZBGWbVuNjiIiIiGuwkPf0dHRADgcDkaN\nGsXo0aOZNm0aJpOp/PHCwkIcDgd2u/2c5zkcjnPuP7NtZSQkRGG1Wqq1Q76wPbkVgzYsInHHJujb\ny7AcgSYpyX7xjaRCGsOa0xj6hsax5mprDCssaoCDBw8ycuRI0tPTueWWW3j55ZfLHysqKiI2NpaY\nmBiKiorOud9ut59z/5ltKyM//2RV98OntiW3BuDUN8tx3DjE0CyBIinJTm6uLlZSExrDmtMY+obG\nseZ8PYYVlX6Fh76PHj3K8OHDGTt2LEOHDgWgQ4cOrFy5EoClS5fSs2dPunTpQlZWFk6nk8LCQnJy\nckhNTaV79+4sWbKkfNsePXr4ap/8aldSc7xWK9a1q42OIiIiIa7CGfWsWbMoKChg5syZ5QvBnnnm\nGV544QWmT59Oy5YtGTx4MBaLhYyMDNLT0/F6vYwZM4bw8HDS0tIYN24caWlp2Gw2XnnllVrZqZpy\nWcMo7dCpbEGZywU2m9GRREQkRJm8Xq/X6BA/5+tDMsOnflXl58zL/YTIzPfJW7QMd+cuPs0TjHSo\nrOY0hjWnMfQNjWPNBcyh71BW2rUbADYd/hYREQOpqC/gTFFb164xOImIiISyi676DlUPLDjMPIuN\nA18sYXRS5Q6dvzd+oJ9TiYhIqNGM+gLcFis7klrQ7OhubKUlRscREZEQpaKuwPbk1lg9bprn7jI6\nioiIhCgVdQW2X9IKgDaHtxucREREQpWKugLbktsAKmoRETGOiroC+xIbU2wNp/UhFbWIiBhDRV0B\nj9lCTnJLmubtI6LklNFxREQkBKmoL2LrJW2weD20PpxjdBQREQlBKuqL2NKwLQDtDm4xOImIiIQi\nFfVFbD5d1G1V1CIiYgAV9UUcs9cnN6ZeWVEH3vVLRESkjlNRV8KWhm1JOHmC5BOHjY4iIiIhRkVd\nCZsb6XNqERExhoq6EvQ5tYiIGEVFXQk7GrTEZbZqRi0iIrVORV0JLmsYOQ1a0iJ3F+Eup9FxREQk\nhKioK2lLo7ZYPW5a6bzfIiJSi1TUlbRZJz4REREDqKgraUv5grKtBicREZFQoqKupFx7fY5FJ5bN\nqHXiExERqSUq6soymdjcqC2JRfk0KDhidBoREQkRKuoq0AU6RESktqmoq2Bjo/YAdNy/0eAkIiIS\nKlTUVZCT3BKnNYwO+1TUIiJSO1TUVVBqsbG5YVuaH9tDzKlCo+OIiEgIUFFX0YbGHQDocGCTwUlE\nRCQUqKiraEOTsqLuqMPfIiJSC1TUVbSlYVtKzRY67t9gdBQREQkBKuoqctoi2N6gFa2O7CDcVWx0\nHBERqeNU1NWwoUkHrB63rk8tIiJ+p6KuhjMLyjrpc2oREfEzFXU1bDp94pMOOvGJiIj4mYq6GhyR\ndnbWb0a7g1uwul1GxxERkTpMRV1NGxt3ILy0hFaHc4yOIiIidZiKuprOfE6t836LiIg/qairaUOT\njgB03rve4CQiIlKXqairKS8mkb2JTei4fyMWd6nRcUREpI5SUddAdtPORLqKST20zegoIiJSR6mo\na2BdShcAuuzNNjiJiIjUVSrqGvixSSc8mLh0j4paRET8Q0VdA45IOzsatKDdwS2Eu5xGxxERkTqo\nUkW9bt06MjIyANi4cSN9+/YlIyODjIwMFi5cCMC8efMYMmQId911F19//TUAxcXFPPHEE6Snp/PQ\nQw+Rl5fnp90wTnbTLtjcpbTX9alFRMQPrBfb4O2332bBggVERkYCsGHDBh544AGGDx9evk1ubi6Z\nmZnMnz8fp9NJeno6ffr0Yc6cOaSmpvLEE0/w2WefMXPmTJ599ln/7Y0B1qV0YUjWv+iy50ejo4iI\nSB100Rl1SkoKM2bMKL+9fv16Fi9ezD333MOECRNwOBxkZ2fTrVs3wsLCsNvtpKSksHnzZrKysujb\nty8A/fr1Y8WKFf7bE4NsbNyeUrNFC8pERMQvLjqjHjx4MPv27Su/3aVLF+688046derEG2+8weuv\nv067du2w2+3l20RHR+NwOHA4HOX3R0dHU1hYWKlQCQlRWK2Wqu6LIYrDItnSMJV2B7ZgsbkhPt7o\nSH6TlGS/+EZSIY1hzWkMfUPjWHO1NYYXLeqfGzRoELGxseW/T5kyhZ49e1JUVFS+TVFREXa7nZiY\nmPL7i4qKyp93Mfn5J6say1DZTbvQcf8mTnz6BSXX32h0HL9ISrKTm1u5P7Tk/DSGNacx9A2NY835\negwrKv0qr/p+8MEHyc4uO8y7YsUKOnbsSJcuXcjKysLpdFJYWEhOTg6pqal0796dJUuWALB06VJ6\n9OhRzV0IbGe+T237ZrGxQUREpM6p8ox68uTJTJkyBZvNRv369ZkyZQoxMTFkZGSQnp6O1+tlzJgx\nhIeHk5aWxrhx40hLS8Nms/HKK6/4Yx8Mt+WSVJzWMMK+WULRxTcXERGpNJPX6/UaHeLnfH1IZvjU\nr3z6eucz6ePf03PXau5/6B2O2etX+nnvjR/ox1S+o0NlNacxrDmNoW9oHGsuoA99y/mtbt4dgO67\n1hicRERE6hIVtY9ktegGqKhFRMS3VNQ+ciC+EYfikum6Zy1mj9voOCIiUkeoqH3FZGJ1827EOE/S\n9uBWo9OIiEgdoaL2oazyz6lXG5xERETqChW1D/3YtDMus5Ue+pxaRER8REXtQ6fCItnYuD1tDm8n\n7uRxo+OIiEgdoKL2sdXNy1Z/d929zuAkIiJSF6iofezM96l76HNqERHxARW1j+2q34xj0Ql027UW\nk9djdBwREQlyKmpfM5lY3bw78adO0ObQdqPTiIhIkFNR+8H3rS4DoNeOVQYnERGRYKei9oM1zbpS\nYrFx+faVRkcREZEgp6L2A6ctgnUpl9L82B6Sjx8yOo6IiAQxFbWfrNThbxER8QEVtZ+sallW1Jfn\nfG9wEhERCWYqaj/Ji0lkyyVt6LRvA9HFDqPjiIhIkFJR+9H3rXph8Xp08hMREak2FbUfrTx9+LtX\njj6nFhGR6lFR+9Hu+s04HNuAHruysLpdRscREZEgpKL2J5OJla16EeM8Scd9G41OIyIiQUhF7Wff\ntbocgCu3rTA4iYiIBCMVtZ9taNKB45FxXLF9BWaP2+g4IiISZFTUfuYxW1jRpjcJJ0/QYb8Of4uI\nSNWoqGvBt22uBKDPVh3+FhGRqlFR14Ifm3aiIMLOlTr8LSIiVaSirgUes4UVrXuTWJRP+wObjY4j\nIiJBREVdS75NPXP4e7nBSUREJJioqGtJdtPOZYe/t63A5PUYHUdERIKEirqWuC1WvmvVi3pFebQ7\nsMXoOCIiEiRU1LVo+enD31dt/dbgJCIiEixU1LVoXUoXCiLs9N2yTKu/RUSkUlTUtajUYmNZ2z4k\nnDzOpXuyjY4jIiJBQEVdy75ufzUAV29aYmwQEREJCirqWra5YVsOxSVzxfbvCHcVGx1HREQCnNXo\nACHHZGJxu37cvfJDem9fyfCpEVV+iffGD/RDMBERCUSaURtgSfv+APTfvNTgJCIiEuhU1AbYl9iE\nbcmt6L5rDXEnjxsdR0REApiK2iCL2/fH4vVw1RZ9p1pERC5MRW2QpW374jaZGbBpsdFRREQkgKmo\nDXI8OoE1zbrS9tA2Uo7uMTqOiIgEKBW1gf7b6VoABq3/0uAkIiISqFTUBvq+1WUcj4xjwKbFWEtd\nRscREZEAVKmiXrduHRkZGQDs3r2btLQ00tPTmTRpEh5P2SUb582bx5AhQ7jrrrv4+uuvASguLuaJ\nJ54gPT2dhx56iLy8PD/tRnAqtdj4ukN/4k4VcPmO742OIyIiAeiiRf3222/z7LPP4nQ6AXjppZcY\nPXo0s2fPxuv1smjRInJzc8nMzGTu3Lm8++67TJ8+nZKSEubMmUNqaiqzZ8/mtttuY+bMmX7foWDz\n306DABj0ow5/i4jIL120qFNSUpgxY0b57Q0bNtCrVy8A+vXrx/Lly8nOzqZbt26EhYVht9tJSUlh\n8+bNZGVl0bdv3/JtV6xY4afdCF576zVlU8O2dNu9lqSCXKPjiIhIgLnoKUQHDx7Mvn37ym97vV5M\nJhMA0dHRFBYW4nA4sNvt5dtER0fjcDjOuf/MtpWRkBCF1Wqp0o4Es/92upb2B7dwzYZFzL3i7otu\nn5Rkv+g2/mDU+9YlGsOa0xj6hsax5mprDKt8rm+z+adJeFFREbGxscTExFBUVHTO/Xa7/Zz7z2xb\nGfn5J6saK6gta3sVDy1+l2s3LOKD3nfhNVV8oCM3t3J/8PhSUpLdkPetSzSGNacx9A2NY835egwr\nKv0qr/ru0KEDK1euBGDp0qX07NmTLl26kJWVhdPppLCwkJycHFJTU+nevTtLliwp37ZHjx7V3IW6\n7VRYJN+0vYrkgly671xtdBwREQkgVS7qcePGMWPGDIYNG4bL5WLw4MEkJSWRkZFBeno69913H2PG\njCE8PJy0tDS2bdtGWloaH3zwAY8//rg/9qFOWHjpDQDctO7/DE4iIiKBxOT1er1Gh/g5Xx+SGT71\nK5++nr/8cc442h7cyiPDZ3IovuEFtzPiMpc6VFZzGsOa0xj6hsax5gL60Lf4z2ddb8SMlxvXfW50\nFBERCRAq6gDybZsryY+K49r1iwh3OY2OIyIiAUBFHUBKrTa+6HwddqeDfpuXGh1HREQCgIo6wHze\nZTBuk5mb1i6EwFs+ICIitUxFHWCO2evzXevLaZW7k/YHNhsdR0REDKaiDkD/7noTALdlfWJwEhER\nMZqKOgCtb9KRbcmt6b19JQ3zDxodR0REDKSiDkQmEx/3vA0zXm5brVm1iEgoU1EHqOVtruBwbAOu\nWf8VsSdPGB1HREQMoqIOUB6zhU+6/4pwdwk36rSiIiIhS0UdwP7b6RoKw2O4ee1CwnQCFBGRkKSi\nDmDFYZH836XXE3eqgIEbvzY6joiIGEBFHeA+7XYTJRYbQ1d9jMVdanQcERGpZSrqAHc8OoEvOl9H\ncsERBmxabHQcERGpZSrqIDD/sttxWazctfJDzB630XFERKQWqaiDwDF7ff7b6VoanjhMf12sQ0Qk\npKiog8RHl92By2xl2HcfgluzahGRUKGiDhK5sUks6jiQxscPEP7Pj4yOIyIitURFHUQ+7HUHpWYL\nUdP/CKVaAS4iEgpU1EHkSFwy/+10Ldbt24iY+w+j44iISC1QUQeZub2H4Y2MJOrll+DUKaPjiIiI\nn1mNDiBVkxeTyEedbuTOVfP59P/9ln9ednulnvfe+IF+TiYiIv6gGXUQmn/Z7TjCo7nz+/lEFzuM\njiMiIn6kog5CRRExfNjrDuxOB3es+qfRcURExI9U1EHq311v4lh0Ir9a8yn1C3ONjiMiIn6iog5S\nJbZwMvvcQ3hpCfd9k2l0HBER8RMVdRD7quMAtiW34urNS2m/f5PRcURExA9U1EHMazLz1oCHAHho\n8TuYvB6DE4mIiK+pqIPc5kbtWNyuH20O5zBww9dGxxERER9TUdcBf+l7L8XWcO5blkmk86TRcURE\nxIdU1HXAMXt9Puo1hISTx7ln+Wyj44iIiA+pqOuIj3vezv74Rty8diGtD203Oo6IiPiIirqOcFnD\nmHnto1i8HkZ+OROzR9esFhGpC1TUdUh2ShcWdRhA6yM7uGXNZ0bHERERH1BR1zHv9X+Aggg79yyf\nTVKBzlgmIhLsVNR1TEFkLO/1v59IVzEjv3wDvF6jI4mISA2oqOugRR0GsrpZV3rsWs116780Oo6I\niNSAirouMpmYcd1IHOFRPLj4XRqcOGx0IhERqSYVdR111J7EWwMeIspVzJNfzACPTi8qIhKMVNR1\n2Nftr+a7Vr3osm89ke/MMjqOiIhUg4q6LjOZeP3axzgRGUv0lElYNqw3OpGIiFSRirqOOx4dz2vX\nPY7J6ST24fuhqMjoSCIiUgUq6hCwqlUvTj7yGNZtW4l55rdGxxERkSqwVveJt99+OzExMQA0adKE\nRx99lPHjx2MymWjTpg2TJk3CbDYzb9485s6di9VqZcSIEQwYMMBn4aXyip59HtuK5UTOzsTVtz/O\nO+4yOpKIiFRCtYra6XTi9XrJzMwsv+/RRx9l9OjRXH755UycOJFFixbRtWtXMjMzmT9/Pk6nk/T0\ndPr06UMhgBwIAAAQRElEQVRYWJjPdkAqZ/j/fkvDHo/w6qYxWEc9zriVxexLbFLhc94bP7CW0omI\nyIVU69D35s2bOXXqFMOHD+fee+9l7dq1bNiwgV69egHQr18/li9fTnZ2Nt26dSMsLAy73U5KSgqb\nN2/26Q5I5R1MaMiM6x4nylXMM5+8pGtXi4gEgWrNqCMiInjwwQe588472bVrFw899BBerxeTyQRA\ndHQ0hYWFOBwO7HZ7+fOio6NxOBy+SS7VsqztVbQ5tJ0hWf/iN5+/you/Go/XpKUKIiKBqlpF3aJF\nC5o1a4bJZKJFixbEx8ezYcOG8seLioqIjY0lJiaGorNWGRcVFZ1T3BeSkBCF1WqpTjSphL/2zaBl\n7g5653zPXSs/5IPew867XVKSvcLbUnUaw5rTGPqGxrHmamsMq1XUH330EVu3bmXy5MkcPnwYh8NB\nnz59WLlyJZdffjlLly6ld+/edOnShVdffRWn00lJSQk5OTmkpqZe9PXz83VI1p88Zgsv3/Q00//x\nFOnL57K7fjO+a937F9vl5haW/56UZD/ntlSdxrDmNIa+oXGsOV+PYUWlX61jnkOHDqWwsJC0tDTG\njBnDiy++yDPPPMOMGTMYNmwYLpeLwYMHk5SUREZGBunp6dx3332MGTOG8PDwau+I+E5BZCwv/up3\nlFjDeHrhdNoc2mZ0JBEROQ+T1xt410H09V96w6d+5dPXq0suy/meZxZMpSDSzti0P3I4Lrn8sbNX\nfesv8JrTGNacxtA3NI41F/Azaqk7VrXqxVsDfk3CyRNM/OcUoou12E9EJJCoqIWFXW/knz1uJSVv\nH88seAlbaYnRkURE5DQVtQDwfr/7WNbmSjrv28D4f/8Ri7vU6EgiIoKKWk7zmsy8csMYVjfrSq8d\nP/Cbz18Ft9voWCIiIU9FLeVKrTZe/NXv2NC4Pf22LCPmqVHg8RgdS0QkpKmo5RxOWzi/v+1ZtiW3\nInJ2JjG/e1plLSJioGpfPUvqrpPh0UwaMok/fPgcLd5/h29W7eL1a0fgMV/8bHG6kIeIiG9pRi3n\nVRgZyzN3TmFbciuuW/8lYz5/DbNHn1mLiNQ2FbVcUGFkLM8O/T0bG7Xj6s1LGffvl7GWuoyOJSIS\nUlTUUqEzh8Gzm3Tiyu3fMfmfv9dJUUREapGKWi6qOCyS529/juWte3Pp3h+Z9sEE6hfmGh1LRCQk\nqKilUkps4Uy7eSwLut1Ms2N7+J/Z42hxZIfRsURE6jwVtVSax2zh7QG/5p3+w6lXlMfUDyZwxbYV\nRscSEanTVNRSZZ/0+BUv3fxbzF4vEz6dxj3f/gOTV9+1FhHxBxW1VMvy1CsZmzaNQ3HJ3L3yQ579\n5EWinEVGxxIRqXNU1FJtu5KaMyb9f1jT7FJ67fiBV//+G6xrsoyOJSJSp6iopUYckXYm3z6Reb3u\nIPnEEeJvGkTk63/SaUdFRHxERS015jFbyLwqg+eGTsaTWI+Y558lLu0OTEeOGB1NRCToqajFZ7JT\nLiX/6+U4rxlE2NeLSOzXi/B/fgRer9HRRESCli7KIT71wLs/Yuoygpvdzbhv2d+IfWQ4K/74JjOv\nfZTj0Qm/2F4X8RARqZhm1OJzXpOZT7vfzBMZr7G+cUeuyFnJzL88wTUbFulrXCIiVaSiFr85mNCQ\nCXdN4Y2BD2P1lDL6ixlM/WCCzmgmIlIFKmrxK6/JzMKuNzLyvhl82+YKOhzYzP/+42keWfSWLu4h\nIlIJ+oxaakVubBJTbxlHt11rePjrt7l53UL6bfkGc/IeTj3wa4iIMDqiiEhAMnm9gbckNze30Kev\nN3zqVz59PakZa6mLW1d/ytBVHxHjPMkRexJ/75POknb98JgtF31+sC5AS0qy+/y/7VCjMfQNjWPN\n+XoMk5LsF3xMh76l1pVabczvNYSHh8/i4x63EX/yOL/5/DVeyxzDlVuXa8GZiMhZVNRimMLIWN7v\nfz+PPjCTLzsOpGnePn737z/y+l9HMWDj15g9bqMjiogYTkUthsuNTeK1waN47P4Z/LfjNTQ8fpDf\nfP4ab77/GDeuXUhEySmjI4qIGEZFLQHjQEJj/jT4CR4e/gafXXoDiY48Rnz1Fu+//SDDl7xH8vFD\nRkcUEal1WvUtASc3tgGzrnmEub3v4vrsL7hx3efcnrWAW7M+5ftWl2Hr5cbVfyBYLr7wTEQk2Kmo\nJWAdj05g7hV381GvO+izdTm3rPk3vXO+h7vv4GhMPb7qMIAvO17DwYSGF32tYF0pLiKiopaAV2qx\nsaR9f5a070+bQ9u4dv0i+m9eyl3ff8Rd33/E+sYdWNK+P8tb96YgKs7ouCIiPqWilqCy7ZI2bLuk\nDe/1f4De279j0PovuXTvj3Tav5FHF71JdkpnlqVexYrWl1MYGWt0XBGRGlNRS1By2sLLZ9lJBblc\nuW05fbcso9vudXTbvY7HvnyDH5t24ocWPfm+ZU+j44qIVJvOTCZ1SoMTh+mzdTl9t35Lm8Pby+8v\nbdmKkkGDKbl2MK7LrzDklKU6G1TNaQx9Q+NYc7V5ZjIVtdRZiY48eu7MoueOH+i6Zx2RrmIAnJYw\nNjdqS3bTzvzYtDPbLmlNqcX2i+f7egGa/s+x5jSGvqFxrLnaLGod+pY6Ky8mkf90HsR/Og/CWuqi\n0/4N9NyZRee9P3Lp6X8AxdZwNjZuz8bG7dnSsC3bkltTFBFjcHoRkTIqagkJpVYba5t1ZW2zrgDE\nniqg0971dN67ni77fqT77rV03722fPs9iU2IOXg1pT0uw9WtB+627SAszKj4IhLCVNQSkgoiY1me\neiXLU68EIO7kcVIPbqXdwa20PbiFNoe2ETnn7zDn7wB4rVbcqe0o7diJ0o6dy39669c3cjdEJASo\nqEWAE1HxrGrVi1WtegFg9rh5//bG2LJWYV2zGuvGH7Fu2oh143r4cG7589wNknG3ScXdqg3uNm1w\nt25DaetUPE2a6sxpIuITKmqR8/CYLdz3ySGgKSQ1hf63Yu7r5pITh2h5ZCctcnfRIncXKcf2kPTt\nMsK+/eac53vDw3G3bI27eQvcKc3wpKRAp3ZY4hrgSUnBG3PhhSMiImdTUYtUksds4UBCYw4kNGZZ\n26vK7w9zOWl0/CBN8vbROH8/TfL2l/3cnkPkpg3nvEbi6Z8FEXYOxzUg155E16svxXNJQzzJl5z1\n8xK8cfFgMtXiHopIIFJRi9RQiS2cXUnN2ZXU/NwHvF5iTxWQfOIwyQVHaFBwhOSCIySfKPs95ehe\n2hzOge3fnfd1nZYw8mISyItJ5ERkHAWRsZyIiqMg0s6wO3vjqVcfb716eBLL/hEd7f+dFZFap6IW\n8ReTiYKoOAqi4tjWMPWXj3u9xJ06QaIjj0RHHvWK8st+L8o762c+7Q5sweL1nPvcxe/98uUiI/Ek\nJOKNjcUbG4cnNrb8959ux5XdFxeHx37696govFHReKOjy04Eo1m8SEDxe1F7PB4mT57Mli1bCAsL\n44UXXqBZs2b+fluRwGcycSIqnhNR8exs0PKCm5k9bmKKHcSeKiD2VAFxp3+W/Ssk7tQJYk8Vlt12\nFBB1bC9RJVuw/bzcK8GDiWJbOE5bBMW2cJIa1T+nyL1RUeU/CY/AGx6ONzwCwsPwnr5NeDjesHCI\nKPvpDY846/fTj595ji0MbDYwm/UHgsgF+L2ov/zyS0pKSvjggw9Yu3YtU6dO5Y033vD324rUGR6z\npXxmXmleL5GuYqKcRUQ7TxLtLCLKWUSMs4iokpNn3XeSCFcxEa5iwl1OIkqLiSgpJsLlJMJVTPGO\nnUS4nNg8bv/t4GkusxW32UKpxYLHZKHUcvq22YL79GM5p38vtVho3aweWG1gs+K12cBqw2u1lq22\nN5vBYsF75nezBcymn+4zmX+2nfmn+ywWvOXPMYPFfJ77zryO6ac/Ms78POuf98zvmC6wDefZ/pev\n88vXNV/wNcr/nX4vLz+7HyAxGkte0bl/HP38D6WKblf4GBd+DMry+Po9/bSt9+c7c+Yxkwnq195J\nkfxe1FlZWfTt2xeArl27sn79en+/pYiYTJwKi+RUWCTHfLDA3Op2lRW566cSt7ldhJWWYHO7sJW6\nym67S376/fRPm7vk3NulZduFlbqweEqxeNxYPG6sHjcWTylWtxuL143F7cbqKcXmLiXC5Ty9TWn5\n/ZYDVT9iID9JvPgmUpEJE2D0+Fp5K78XtcPhICbmp788LBYLpaWlWK0XfuuKznlaHZ++cqtPX09E\nRCSplt7H7O83iImJoaioqPy2x+OpsKRFRETkJ34v6u7du7N06VIA1q5dS2rqeVa/ioiIyHn5/TKX\nZ1Z9b926Fa/Xy4svvkirVq38+ZYiIiJ1RkBej1pERETK+P3Qt4iIiFSfilpERCSA1emi9ng8TJw4\nkWHDhpGRkcHu3buNjhR0XC4XY8eOJT09naFDh7Jo0SKjIwWtY8eO0b9/f3JycoyOEpTefPNNhg0b\nxpAhQ/jwww+NjhN0XC4XTz31FHfffTfp6en677CK1q1bR0ZGBgC7d+8mLS2N9PR0Jk2ahMfj3+/0\n1+miPvusaE899RRTp041OlLQWbBgAfHx8cyePZt33nmHKVOmGB0pKLlcLiZOnEhERITRUYLSypUr\nWbNmDXPmzCEzM5NDhw4ZHSnoLFmyhNLSUubOncvIkSN59dVXjY4UNN5++22effZZnE4nAC+99BKj\nR49m9uzZeL1ev09g6nRR66xoNXf99dfz5JNPAuD1erFYLAYnCk7Tpk3j7rvvpkGDBkZHCUrLli0j\nNTWVkSNH8uijj3L11VcbHSnotGjRArfbjcfjweFw6HwWVZCSksKMGTPKb2/YsIFevXoB0K9fP5Yv\nX+7X96/T/0tV56xocq7o05dOdDgcjBo1itGjRxucKPh8/PHHJCYm0rdvX9566y2j4wSl/Px8Dhw4\nwKxZs9i3bx8jRozg888/x6QLeVRaVFQU+/fv54YbbiA/P59Zs2YZHSloDB48mH379pXf9nq95f/t\nRUdHU1hY6Nf3r9Mzap0VzTcOHjzIvffey6233sott9xidJygM3/+fJYvX05GRgabNm1i3Lhx5Obm\nGh0rqMTHx3PVVVcRFhZGy5YtCQ8PJy8vz+hYQeUvf/kLV111FV988QWffPIJ48ePLz+UK1VjNv9U\nnUVFRcTGxvr3/fz66gbTWdFq7ujRowwfPpyxY8cydOhQo+MEpX/84x/8/e9/JzMzk/bt2zNt2jSS\nkmrrLMF1Q48ePfjmm2/wer0cPnyYU6dOER8fb3SsoBIbG4vdXnYdhbi4OEpLS3G7/X9VtLqoQ4cO\nrFy5EoClS5fSs2dPv75fnZ5eDho0iG+//Za77767/KxoUjWzZs2ioKCAmTNnMnPmTKBsYYUWRUlt\nGjBgAKtWrWLo0KF4vV4mTpyo9RJVdP/99zNhwgTS09NxuVyMGTOGqKgoo2MFpXHjxvHcc88xffp0\nWrZsyeDBg/36fjozmYiISACr04e+RUREgp2KWkREJICpqEVERAKYilpERCSAqahFREQCmIpaREQk\ngKmoRUREApiKWkREJID9fx57eZnmmKzSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10f8ddcf8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "P = lambda x: np.exp(-x)\n",
    "\n",
    "# domain limits\n",
    "xmin = 0 # the lower limit of our domain\n",
    "xmax = 10 # the upper limit of our domain\n",
    "\n",
    "# range limit (supremum) for y\n",
    "ymax = 1\n",
    "#you might have to do an optimization to find this.\n",
    "\n",
    "N = 10000 # the total of samples we wish to generate\n",
    "accepted = 0 # the number of accepted samples\n",
    "samples = np.zeros(N)\n",
    "count = 0 # the total count of proposals\n",
    "\n",
    "# generation loop\n",
    "while (accepted < N):\n",
    "    \n",
    "    # pick a uniform number on [xmin, xmax) (e.g. 0...10)\n",
    "    x = np.random.uniform(xmin, xmax)\n",
    "    \n",
    "    # pick a uniform number on [0, ymax)\n",
    "    y = np.random.uniform(0,ymax)\n",
    "    \n",
    "    # Do the accept/reject comparison\n",
    "    if y < P(x):\n",
    "        samples[accepted] = x\n",
    "        accepted += 1\n",
    "    \n",
    "    count +=1\n",
    "    \n",
    "print(\"Count\",count, \"Accepted\", accepted)\n",
    "\n",
    "# get the histogram info\n",
    "hinfo = np.histogram(samples,30)\n",
    "\n",
    "# plot the histogram\n",
    "plt.hist(samples,bins=30, label=u'Samples');\n",
    "\n",
    "# plot our (normalized) function\n",
    "xvals=np.linspace(xmin, xmax, 1000)\n",
    "plt.plot(xvals, hinfo[0][0]*P(xvals), 'r', label=u'P(x)')\n",
    "\n",
    "# turn on the legend\n",
    "plt.legend();\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that $y_{max}$ was just assumed here. In general we might have to do a maximization. This has a cost. We want to keep this cost low, or we might be spending some time there. If the optimization is complex, it might be cheaper to  do something else..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rejection Sampling with Steroids\n",
    "\n",
    "The simple rejection sampling method has fundamental problems.  \n",
    "\n",
    "For our simple example, it's quite easy to determine the supremum.  In practice, while you may know how to quickly (i.e. constant time) evaluate your function everywhere on the domain of interest, finding a bound very close to the  supremum may not be a feasible calculation.  In addition, even if you find a tight bound for the supremum, basic rejection sampling will still be very inefficient as you will reject many samples (especially in low density regions). \n",
    "\n",
    "Furthermore, if you support is infinitely large, you are not going to be able to reject from an infinitely long box in finite time!\n",
    "\n",
    "This is a hard problem to solve and we will need other techniques to address this problem of **low acceptance probability**.\n",
    "\n",
    "However, it is possible to do a more efficient job while still taking advantage of the simplicity of rejection sampling.  Our modified technique will introduce a **proposal density** $g(x)$. This notion of a proposal density is one used in many sampling techniques, in different ways, but its importance will always lie in figuring ways to increase the acceptance rate.\n",
    "\n",
    "The proposal density will have the following characteristics:\n",
    "\n",
    "- $g(x)$ is easy to sample from and (calculate the pdf)\n",
    "- Some $M$ between 1 and $\\infty$ exists so that $M \\, g(x) > f(x)$ in your entire domain of interest\n",
    "- ideally $g(x)$ will be somewhat close to $f$ so that you'll sample more in high density regions  and much less in low density regions\n",
    "\n",
    "Its probably obvious that an optimal value for M is the supremum over your domain of interest of $f/g$. At that $x$ you will accept stuff with probability 1. In general you want $M$ as close to 1 as possible, since the probability of acceptance is $1/M$.\n",
    "\n",
    "You can see that this is the case by finding the proportion of samples from $g(x)$ that are accepted at each $x$ and then averaging over $x$:\n",
    "\n",
    "$$\\int dx g(x) prop(x) =  \\int dx g(x) \\frac{f(x)}{Mg(x)} = \\frac{1}{M}\\int dx f(x) = \\frac{1}{M}$$\n",
    "\n",
    "Once you've picked a proposal distribution g, your modified rejection sampling technique is as follows:\n",
    "\n",
    "1. Draw $x$ from your proposal distribution $g(x)$\n",
    "2. Draw $y$ uniformly from [0,1]\n",
    "3. if $y$ < f($x$)/$M g(x)$, accept the sample\n",
    "4. otherwise reject it\n",
    "5. repeat\n",
    "\n",
    "The entire process is illustrated in the diagram below:\n",
    "\n",
    "![](../images/rejsteroid.png)\n",
    "\n",
    "**Example**\n",
    "\n",
    "The following code produces samples that follow the distribution $P(x)=e^{-x}$ \n",
    "  for $x=[0,10]$ and generates a histogram of the sampled distribution. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Count 23692 Accepted 10000\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAFJCAYAAABU5W56AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XlcVPX+x/HXwLAJg4AMuIGKijuKey6YltcW65qWBv3M\nsuxmZmrZ1crSSiszl7LboqkVZabpzSxtc981VEwUFTT3BQVkUVnn94fK1dyQGZgB3s/HgwfDmTPn\nfM4n8s3ZvsdgsVgsiIiIiENysncBIiIicn0KahEREQemoBYREXFgCmoREREHpqAWERFxYApqERER\nB2a0dwHXkpSUbtPl+fpWICXlrE2XWd6oh9ZTD62nHtqG+mg9W/fQbDZd971ysUdtNDrbu4RSTz20\nnnpoPfXQNtRH65VkD8tFUIuIiJRWCmoREREHpqAWERFxYApqERERB6agFhERcWAKahEREQemoBYR\nEXFgCmoRESnToqNnER+/85rvJSYmMHPmtBKu6NYoqEVEpMw6ceI4CQl7qV+/4TXfr127DkeOHOLI\nkcMlXFnhOeQQoiIiUnp5jhmF26LvbbrMrPt6kDlm7A3nWbx4EatXr+Ds2bOkpqby+ONPsnt3PJ07\n3wHAf/7zPs7Ozjz11DMMGzaIPn0eoV27DnTu3JUFC+YyePDzNq3ZVrRHLSIiZca5c+eYPPk/TJ78\nIVOnTmbr1hhq164LwL/+NYgtW/5g3LjRNGjQiHbtOgBQp05dtm6NsWfZN1Tm96gn3X8XFSq48/Qc\n2/51JyIi15Y5ZuxN936LS7NmzXFycsLPrxImkzepqSn4+fkBYDQa6d07krFjR7NgwU8Fn6lUyZ8z\nZ87Ypd7CKPNBvWhrDIeystj41u8YnK4+gDBzZBc7VCUiIsVh9+54AJKTT5OZmYmfnx/p6Rl4enqR\nlpZGdPQsBg8exvjxYxk/fjIA6elp+Pr62bPsGyrzh77r+PiSBrglH7J3KSIiUsySk08zZMhAXnxx\nKC+8MIIWLVqxc+cOAN55502ioh6lV68+eHtXZN68OQDs3LmDFi1a2bPsGyrze9S1q1WHE8dxO7yD\nbP8a9i5HRESKUbNmzRk4cHDBzzVr1uLDD6fQpcudvPXWhILpr7wypuD10qW/MWDAwJIs85aU+T3q\n2qH1AMg/nmDnSkREpKRVrlyFOnXqXvc+6oSEvVSrVp2qVauVcGWFV+b3qEOaNYc5X3Mu2XHvkRMR\nEevdc89915z+2GNPXvczderUpU6dusVVkk2U+T3qmhcvv09JP2XnSkRERG5dmQ9qv9B6+ABHz6fb\nuxQREZFbVqigjo2NpW/fvgCcPn2agQMH8sgjj/Dwww9z8OBBAObOnUvPnj3p3bs3y5cvB+D8+fMM\nHjyYqKgoBgwYQHJycjFtxvUZnJwINrpxMC8XcrNLfP0iIiLWuOk56unTp/PDDz/g4eEBwIQJE7jv\nvvu455572LBhA/v27cPDw4Po6Gjmz59PVlYWUVFRtG/fnm+++YbQ0FAGDx7MTz/9xEcffcSoUaOK\nfaP+LrBCRbanncT92B7OBzUu8fWLiIgU1U33qIODg5k6dWrBz1u2bOHEiRM89thjLFq0iNatW7N9\n+3bCw8NxdXXFZDIRHBxMfHw8MTExdOzYEYCIiAjWr19ffFtyAxUrBgJgPHLtq/5ERKRs+vXXJaxc\nueya7yUnn2bSpPElXNGtu+kedbdu3Th8+H9XTB85cgRvb28+//xzPvzwQ6ZPn07NmjUxmUwF83h6\nepKRkUFGRkbBdE9PT9LTC3ee2Ne3Akaj861uy3W5+QfDoT/JPbnvqvfMZtM1PiHXol5ZTz20nnpo\nG+Whj2fPnmXZsl+YMWPGNd83m034+/uyf/8uWrdufcvLL6ke3vLtWT4+PnTpcmHYzS5dujB58mQa\nN25MZmZmwTyZmZmYTCa8vLwKpmdmZuLt7V2odaSknL3Vsm7IUqUebP2JzJSj/D3+k5J0kVlhmM0m\n9cpK6qH11EPbKO4+jhkzikU2fnrWfff1YMwNxg/PyjrPm2+O5vTpJAICAtm2bSv9+z9F06YtSUpK\nZ/78uWzfvo3XX3+LsWNH07BhY3r2fIj27bswffqn1KrV4JbqsXUPbxT6t3zVd4sWLVi5ciUAmzdv\npk6dOoSFhRETE0NWVhbp6ekkJiYSGhpK8+bNC+ZdtWoVLVq0KOImWOdc9QvnpU9llvzFbCIiUvwW\nLvwvVatW5eOPZ9K//79ISUlm69YY6tSpA0CvXr3Jzs5i3Lgx5OTk0LPnQ8CFkcu2b99mz9Jv6pb3\nqEeMGMGoUaOYM2cOXl5eTJw4kYoVK9K3b1+ioqKwWCwMGzYMNzc3IiMjGTFiBJGRkbi4uDBx4sTi\n2IabyvHyo6rBwOGsswTapQIRkfJjzJixN9z7LQ4HDuynTZt2ANSoURMfH1/OnEnF17dSwTyPPPIY\nTz/9ODNmfFUwzdnZGaPRSH5+Pk7XeHCTIyhUUFevXp25c+cCUK1aNWbNmnXVPL1796Z3795XTPPw\n8OCDDz6wQZnWq+FagfVZmVgyUzF4+ti7HBERsaGQkNrs2LGdiIjbOXLk8MWQ9iMj48Lh6ZycHD74\nYCIvvvgyEye+w3/+Mx0XFxcsFgvOzs4OG9JQDgY8uaSy14W/qlwPbbdzJSIiYmvdu/+T48ePMWjQ\nAGbO/BRXV1fCw1sUPDnr448/oF27Dvzznz1p0+Y2Pvnkwt1MiYkJNG7cxJ6l31SZH+v7Em+/anD6\nIJYju6B+hL3LERERG9qzZzfdu/+T1q3bcujQQf78czt33NGVl14aTlTUozz33AsF8z7xxL8KXv/6\n62J69uxjj5ILrdzsUbsGXrig4OypA3auREREbK1q1WpER89i4MD+vPHGKJ5/fgQVKnhy1133smLF\n0mt+5vTpU2RmZtK0abMSrvbWlJs96tyLI5KdPnOCKnauRUREbKtSJX+mTv30qul33939hp958cWX\ni7Msmyg3e9Qp5hCqAkfOnbF3KSIiIoVWboI61+hCHaMrx/Nyycmy7YAqIiIixaXcBDVANY+KAOQe\n32vnSkRERAqnXAW1n09lAJwO7bBzJSIiIoVTbi4mA3ANCIFDf5JzMgE3excjIlKG9X/n2k+sKqqZ\nI7sUar7o6M/5449N5OXlYjAYGDRoKPXr39o43jeyZcsfLFw4n48++tBmy7yZchXU+UFNIGYhqSlH\n8bJ3MSIiYlP79+9j7dpVfPzxDAwGA3v37mbs2DF88cU39i7NKuUqqJMq16UqcDwzher2LkZERGzK\ny8uLEyeO89NPC2nTph1169Zj+vQv2Lo1hlmzppOfn8+5c+cYPXosLi4uvPbaSwQGBnLs2DHuuOMf\n7N+fyJ49u2nXrgP/+tcgnn32KWrUqMmBA38B8Prrb12xvmXLfufbb7/GycmJsLBmDBw4mO3bt/Hh\nh1MwGo24u7szdux4KlTwtGq7ylVQp1WoSKiTMytys8nNPofR1cPeJYmIiI2YzQG8884k5s//lpkz\np+Pu7s5TTz1DcnIyr732Jv7+Zr78cibLl//OP/5xN8eOHWHy5P+QlXWehx76J99/vxg3N3cefPA+\n/vWvQQA0bhzGiy++zIIF84iOnkVERGcAUlNTmTnzUz77LBp3d3fefPNVNm/ewKZNG+nS5U56945i\nzZpVpKWlK6hvVXUPb8hM4dzJ/ZiqN7R3OSIiYiOHDx/C09OTl18eDUB8/E6GD3+OQYOGMmXKBDw8\nKpCUdJImTZoCUKVKNby8vHBxccHPzw9v7wt3BhkMhoJltmjRCoAmTcJYs2ZlwfSDBw+SmprC8OHP\nAXD27FmOHDlM376P8+WXMxkyZCBmcwANGza2ervK1VXfcNmV3wf1cA4RkbIkMXEvkya9S05ODgBB\nQcF4eZl4//2JvPzyaF55ZQz+/uaC+S8P5OvZvXsXANu3x1KrVkjB9OrVqxMQEMiUKR/x4YfTePDB\nPjRq1IRff13MPfd0Z+rUT6lVK4Qfflhg9XaVuz1qr4A6cGQXWcd3Y93BCBERcSSdOnXhr7/28+ST\nj1Khggf5+RaeeWYIsbFbeOaZAXh4uOPrW4lTp5IKvczFi3/k229n4+7uzquvvkFiYgIAfn5+9Onz\nCM8++xR5eXlUqVKVLl26kp2dwzvvjMXDwwODwcC///2K1dtlsFgsFquXYmNJSek2Xd7ltwk0TtjE\nuz+8RajJTJ0B0wt9yX95ZzabbP7fpbxRD62nHtqG+lg4zz77FC+++DI1atS86j1b99BsNl33vXJ3\n6PtEQC1CgcOZyTjg3ygiIiJXKHeHvpNM/jQ2OBGfn8e5tJP2LkdERBzUhx9Os3cJQDnco8ZgoIan\nDwDpJ/fZuRgREZEbK39BDfj7VgUgT2N+i4iIgyuXQe1eORSAcyf0FC0REXFs5TKoMyrXxQwkpRy1\ndykiIiI3VC6D+nClYJoCp85nkJ6eZu9yRERErqtcBvUxn8qEcWFEmri4ODtXIyIicn3lMqjznZwJ\n8vIFIG6HhhIVERHHVS6DGsDP78KDLuP+2GTnSkRERK6vUEEdGxtL3759r5i2aNEi+vTpU/Dz3Llz\n6dmzJ71792b58uUAnD9/nsGDBxMVFcWAAQNITk62YenWcQoMxQ3YsW2rvUsRERG5rpsG9fTp0xk1\nahRZWVkF03bu3Ml3331XMARnUlIS0dHRzJkzhxkzZjBp0iSys7P55ptvCA0NZfbs2fTo0YOPPvqo\n+LbkFh0116ApEHdg/xXbJiIi4khuGtTBwcFMnTq14OeUlBQmTZrEyy+/XDBt+/bthIeH4+rqislk\nIjg4mPj4eGJiYujYsSMAERERrF+/vhg2oWgO+AfTEsjJy2PXLl1QJiIijummQd2tWzeMxgtDgufl\n5fHKK6/w0ksv4en5v4dEZmRkYDL978kfnp6eZGRkXDHd09OT9HTHeVrLUd+qNL/4LNKtW7fYuRoR\nEZFru6WHcsTFxXHgwAHGjBlDVlYWCQkJjBs3jrZt25KZmVkwX2ZmJiaTCS8vr4LpmZmZeHt7F2o9\nvr4VMBqdb6W0W5br7EJ170A4c5zd8X/e8BFjcoF6ZD310HrqoW2oj9YrqR7eUlCHhYXx008/AXD4\n8GGef/55XnnlFZKSkpgyZQpZWVlkZ2eTmJhIaGgozZs3Z+XKlYSFhbFq1SpatGhRqPWkpJy99S0p\nApfAOlQ4c5wNq9fo2aw3oefXWk89tJ56aBvqo/VK8nnUNnnMpdlspm/fvkRFRWGxWBg2bBhubm5E\nRkYyYsQIIiMjcXFxYeLEibZYnc0cNteg+Z41rEtMIDMz84rD+SIiIo7AYLl06bYDsfVfev3fWXbN\n6a0SN5G68C0mA4sW/UqbNm1tut6yRH+BW089tJ56aBvqo/VKco+63A54AnDQvwYtL77eti3GrrWI\niIhcS7kO6pPeZlp4VABgmwY+ERERB1Sug9picKJWgwZUBLZt1R61iIg4nnId1AB5DRvTEkjcl0ha\n2hl7lyMiInIFBXX9BrS6+FoDn4iIiKMp90Gd26AR7S6+/kNP0hIREQejoK7fkNsuvt68eaNdaxER\nEfm7ch/UFrMZP38zdY1G/vhjM/n5+fYuSUREpEC5D2qA3AYNaZebS1raGfbs2W3vckRERAooqLkY\n1Bdf6/C3iIg4EgU1kFdfQS0iIo7JJg/lKM36v7OM0GPnmABUcDKy6JflZFwcG3zmyC72LU5ERMo9\n7VEDB/yDAQPNXNzJTDlK9rk0e5ckIiICKKgByHJx56hvVTrnZgGQckwXlImIiGNQUF+0L6AWt+fl\nAJB8JN7O1YiIiFygoL5ov7kWbQAng4HkIzvtXY6IiAigoC6QGBCCCahZwZfU43vJy8myd0kiIiIK\n6kv2m2sBcJuLK5b8XFKO6fC3iIjYn4L6olRPH5I9fbk76ywApw/F2bkiERERBfUV9ptr0f1cGgaD\ngdOHd9i7HBEREQX15RIDalERqOxThdTjezh37py9SxIRkXJOQX2ZS+epm3j6kp+XS0zMZjtXJCIi\n5Z2C+jL7Ai4EdQQGANauXW3PckRERBTUlzvmU4VzLu50P3sGMLB+/Vp7lyQiIuWcgvoyFoMT+801\naZxyBB9zTf74Y5POU4uIiF0pqP9mv7kWzpZ8QvyDyc7OZtOmDfYuSUREyjEF9d/sCwgBoI2HNwAr\nViyzZzkiIlLOFSqoY2Nj6du3LwC7du0iKiqKvn378sQTT3Dq1CkA5s6dS8+ePenduzfLly8H4Pz5\n8wwePJioqCgGDBhAcnJyMW2G7ey7eOV3l9xs3NzcWLlyuZ0rEhGR8uymQT19+nRGjRpFVtaFsa/H\njRvHq6++SnR0NF27dmX69OkkJSURHR3NnDlzmDFjBpMmTSI7O5tvvvmG0NBQZs+eTY8ePfjoo4+K\nfYOsdbBSEHkGJxqeOkibNu3YsWM7J0+etHdZIiJSTt00qIODg5k6dWrBz5MmTaJBgwYA5OXl4ebm\nxvbt2wkPD8fV1RWTyURwcDDx8fHExMTQsWNHACIiIli/fn0xbYbtZLu4ccivOrWS9tMpohMAq1ev\nsG9RIiJSbhlvNkO3bt04fPhwwc8BAQEAbNmyha+++oqvv/6a1atXYzKZCubx9PQkIyODjIyMgume\nnp6kp6cXqihf3woYjc63tCG2tLdyXWqePkjPZo15E9iwYTVPP/2E3epxFGaz6eYzyQ2ph9ZTD21D\nfbReSfXwpkF9LYsXL+bjjz9m2rRp+Pn54eXlRWZmZsH7mZmZmEymK6ZnZmbi7e1dqOWnpJwtSlk2\nkxBYm65xS6l5+AT+/v788suvnDx5YQzw8spsNpGUVLg/tOTa1EPrqYe2oT5az9Y9vFHo3/JV3wsX\nLuSrr74iOjqaoKAgAMLCwoiJiSErK4v09HQSExMJDQ2lefPmrFy5EoBVq1bRokWLIm5CydobWAcA\ntz9jiYi4nePHj7F7tx57KSIiJe+WgjovL49x48aRmZnJ4MGD6du3Lx988AFms5m+ffsSFRVFv379\nGDZsGG5ubkRGRrJ3714iIyP59ttvefbZZ4trO2zqL3NNcp2cMW7bwu233wHA8uVL7VyViIiURwaL\nxWKxdxF/Z+tDMv3fufV7oSd/9Ty1zxwhbv0WmjRvREREZ777bqFN6ypNdKjMeuqh9dRD21AfrefQ\nh77Li4TAOhiysqiakkLTpuGsX7+G9PQ0e5clIiLljIL6Oi6dp3bZtoWuXbuRk5OjUcpERKTEFemq\n7/IgofKFoF77+Y9sCPsHAKMnzuK/cf+7cn3myC52qU1ERMoP7VFfx4FKwWQ7u1D3RAIVA0Jw8/Tl\n5P4tWPLz7F2aiIiUIwrq68hzNrLPXIsapw7gmpdLQK2WZJ87Q+rxBHuXJiIi5YiC+gYSAutgzM+j\nZtJfBIa0BODE/j/sXJWIiJQnCuobSKhcG4C6JxLwr9EUJ2cXTiRutnNVIiJSniiob2BvYF3gQlAb\nXdypFNSE9FN/cS4tyc6ViYhIeaGgvoHDftU4b3SjzsXz0pcOfx9P3GTPskREpBxRUN9AvpMziYEh\nBCUfxj37HJXrtAXg+F7Hf1yniIiUDQrqm9hTuS7OlnzqnEjE3csP36r1OX1kJ1lnU+1dmoiIlAMK\n6pvYXaUeAPWP7QagSt12YMnneMJGe5YlIiLlhIL6JuIvBnW9gqC+DYBje9fZrSYRESk/FNQ3cdrk\nT5JXpQtBbbHg4W2mYmBdTh/8k+Tk0/YuT0REyjgFdSHsrlIP37NnCDxzAoAqobdhseTz88+L7VyZ\niIiUdQrqQoiv+vfz1BcOf//4Y/l9PrWIiJQMBXUh/P08tadPFbzNtVi5cjlnzujqbxERKT4K6kLY\nFxBCjpOxYI8aoEpoO3Jycvjpp0V2rExERMo6BXUh5BhdSQwIoVbSX7jlZAFQrX4EAN999609SxMR\nkTJOQV1Iu6vWw5ifR+0TF4YTrVAxkDZtbmPt2tUcPXrEztWJiEhZpaAupPi/DXwC0KtXbywWCwsW\nfGevskREpIxTUBfS7oILyvYUTLv//h64uLjo8LeIiBQbBXUhJZn8Oe3pd2GP2mIBwM+vEnfc8Q92\n7tzBzp1xdq5QRETKIgV1YRkMxFeth19mCgFpJwsmP/RQHwDmz59rr8pERKQMU1Dfgr8/oAOga9e7\nMJm8mT9/Lvn5+fYqTUREyigF9S3YWbUBAI2O7CyY5u7uzv339+Do0SOsWrXCTpWJiEhZVaigjo2N\npW/fvgAcOHCAyMhIoqKiGD16dMFe5Ny5c+nZsye9e/dm+fLlAJw/f57BgwcTFRXFgAEDSE5OLqbN\nKBmJgSFkGV1peHjnFdMfeeRRAL7++kt7lCUiImXYTYN6+vTpjBo1iqysCwN9vP322wwdOpTZs2dj\nsVhYunQpSUlJREdHM2fOHGbMmMGkSZPIzs7mm2++ITQ0lNmzZ9OjRw8++uijYt+g4pTr7EJ8lXrU\nPH0Qr3PpBdNbtGhF/foNWLx4EadP64laIiJiOzcN6uDgYKZOnVrwc1xcHK1btwYgIiKCdevWsX37\ndsLDw3F1dcVkMhEcHEx8fDwxMTF07NixYN7169cX02aUnLhqDQFoeHRXwTSDwcAjjzxKTk4O8+Z9\nY6/SRESkDDLebIZu3bpx+PDhgp8tFgsGgwEAT09P0tPTycjIwGQyFczj6elJRkbGFdMvzVsYvr4V\nMBqdb2lDSkpc9QtB3ejwTszm/23z008/yZtvjmbOnK8YNWpkQY/Kksu3V4pGPbSeemgb6qP1SqqH\nNw3qv3Ny+t9OeGZmJt7e3nh5eZGZmXnFdJPJdMX0S/MWRkrK2Vstq8TsrlKPXCdnGh2JIynp8j88\n3Ljnnu58//0ClixZSqtWbexWY3Ewm01/2165Veqh9dRD21AfrWfrHt4o9G/5qu+GDRuyceNGAFat\nWkXLli0JCwsjJiaGrKws0tPTSUxMJDQ0lObNm7Ny5cqCeVu0aFHETXAcWS7uJATUpvbJfXDZHycA\njzzSD9BFZSIiYju3HNQjRoxg6tSp9OnTh5ycHLp164bZbKZv375ERUXRr18/hg0bhpubG5GRkezd\nu5fIyEi+/fZbnn322eLYhhIXV70hxvw8XGI2XzG9Y8dOBAfX5L///Y7U1BQ7VSciImWJwWK5OB6m\nA7H1IZn+7yyz6fJaJW7itYVvkfnCCM6OeOWK9z788H3eeONVXn/9LQYOLBt/mIAOldmCemg99dA2\n1EfrOfShb4FdFwc+cdl49VXsUVH/h7u7OzNnTiMvL6+kSxMRkTJGQV0EGR4m9vvXuHDoOzv7ivf8\n/CrRq1dvDhz4i2XLfrNThSIiUlYoqItoZ7WGGM6dwxi79ar3+vd/CoDPPvu0pMsSEZEyRkFdRJcG\nPnHZcPXh7yZNwmjT5jaWL19KYuLeki5NRETKEAV1EcVVbwSAy7rV13z/iScu7FXPnDm9xGoSEZGy\nR0FdRMlefuTWDcV1/TrIybnq/XvvvZ8qVary9dfRulVLRESKTEFthV8q1MZwNpMJz0+j/zvLrvhy\ncXHhqaee4ezZTD7/fIa9SxURkVJKQW2F2OAwAMIObb/m+48++hgmkzfTp3/C+fPnS7I0EREpIxTU\nVvizemPyMdD04LWD2mTy5rHHniAp6STz5s0p4epERKQsUFBbIcPDxL6AWtQ/thu3nKxrzjNgwNO4\nuLjw0UcfkJ+fX8IViohIaaegttL2oDBc8nJpcNnzqS9XuXIVHnroYRITE1iy5KcSrk5EREo7BbWV\nCs5TH/zzuvM888xzGAwGpkx5DwccWl1ERByYgtpKO6s1INfJ+boXlAGEhtbjvvt6EBu7ld9++7kE\nqxMRkdJOQW2l864e7K4SSp0TiXiez7jufC+8MAKA9957R3vVIiJSaApqG9geFIazJZ9GR3Zed54G\nDRpy//0PsG3bVn7//ZcSrE5EREozBbUNXDpPfb3btC55/vl/AzBhwtvaqxYRkUJRUNvA7sqhZBld\nbxrUDRs24r77erBt21aWLv21hKoTEZHSTEFtA7lGF/6s3pgapw9SKf3UDee9dK56/Pi3dF+1iIjc\nlILaRrbUbA5A87+ufj715Ro2bESPHj2Jjd3KokXfl0RpIiJSiimobSSmVjhw86AGGDnyVYxGI+PG\nvU52dnZxlyYiIqWYgtpGjvpU5XjFQJod3IZTft4N5w0JqU2/fv3566/9REfPKqEKRUSkNFJQ24rB\nwJaa4XhlnaXesT03nf3550fg6enFxInjychIL4ECRUSkNFJQ21BMwXnqLTed12w28+yzQzh16hT/\n+c8HxV2aiIiUUgpqG/ozqAk5TkZaFOI8NcDTTz9LQEAgH388lWPHjhZzdSIiUhopqG3onKsHO6s1\noO6JBAxJSTed39PTk5deepWzZ8/yxhuvlUCFIiJS2iiobWxLzQtXf7uuXFao+SMj/49mzcKZP38u\nGzduKM7SRESkFCpSUOfk5PDCCy/w8MMPExUVRWJiIgcOHCAyMpKoqChGjx5dMJjH3Llz6dmzJ717\n92b58uU2Ld4RXbqf2nXZ74Wa38nJibfemgDASy8NJy/vxleMi4hI+VKkoF65ciW5ubnMmTOHQYMG\nMWXKFN5++22GDh3K7NmzsVgsLF26lKSkJKKjo5kzZw4zZsxg0qRJZf6+4b/8a3Da0xfXFUuhkCOP\ntWzZmj59otixYztfffVF8RYoIiKlSpGCulatWuTl5ZGfn09GRgZGo5G4uDhat24NQEREBOvWrWP7\n9u2Eh4fj6uqKyWQiODiY+Ph4m26AwzEY2FKzOU6nTmHcGlPoj40a9TpeXibefvsNUlKSi7FAEREp\nTYxF+VCFChU4cuQId999NykpKXzyySds3rwZg8EAXLhIKj09nYyMDEwmU8HnPD09yci4/jObL/H1\nrYDR6FyU0hzCptqt6Bq3FN81y+CuLoX6jNlsYsyY0QwfPpx3332Tzz77rJirvHVms+nmM8kNqYfW\nUw9tQ320Xkn1sEhB/fnnn9OhQwdeeOEFjh07Rr9+/cjJySl4PzMzE29vb7y8vMjMzLxi+uXBfT0p\nKWeLUpZSUyGoAAAgAElEQVTD2FqjGRY3N/IW/JeUISMK/bnIyMeZNesLZsyYQffuPWnfvmMxVnlr\nzGYTSUkamMUa6qH11EPbUB+tZ+se3ij0i3To29vbuyBwK1asSG5uLg0bNmTjxo0ArFq1ipYtWxIW\nFkZMTAxZWVmkp6eTmJhIaGhoUVZZqmS5uJMdcTvGXTtx+mt/oT/n4uLCpEkf4OTkxAsvPMf58+eL\nsUoRESkNirRH/dhjj/Hyyy8TFRVFTk4Ow4YNo3Hjxrz66qtMmjSJkJAQunXrhrOzM3379iUqKgqL\nxcKwYcNwc3Oz9TY4pOmWEJ4F5v/7AxY1v++K92aOvP7h8PDwFgwY8DSffvoRkye/y0sv6f5qEZHy\nzGCxWCz2LuLvbH1Ipv87hbun2Zb8MpL5Ylp/YoOaMOqhN69470ZBDZCRkUFERBuOHz/G77+vpmHD\nRsVZaqHoUJn11EPrqYe2oT5az+EPfcvNJXv5sbtyXRofjsPz/M0voLucl5cX7747idzcXIYMeeaK\n8/8iIlK+KKiL0abarXG25NOiEA/p+Ls77+zGww8/QmzsViZPnlAM1YmISGmgoC5GG0NaAdA6cXOR\nPj927DtUrx7E5MkT2HoL92SLiEjZoaAuRgf8a3DCO4AWf8VgzLv1w9fe3hX54IOPycvLY9Cgpzh3\n7lwxVCkiIo5MQV2cDAY21m6NV9ZZGh3eWaRFdOgQwVNPDSQhYS/jxo2xbX0iIuLwFNTFbEPtNgC0\n27u+yMt45ZUx1K0byrRpH7N8+VJblSYiIqWAgrqYxVVvSKpHRW5LWI9TftGejOXh4cFHH03HxcWF\nQYMGcPz4MRtXKSIijkpBXczynZxZX7ctvmfP0PBI0Q5/AzRtGs6YMWM5deoUAwc+qcdhioiUEwrq\nErC2bjsA2u8p+uFvgCeffJq77+7O2rWrmThxvC1KExERB6egLgF/BjUmzd1EOysOfwMYDAbef/8/\nBAUFM3HieFavXmnDKkVExBEpqEtAvpMz6+u0xS8zhQZHrXset4+PL9OmzcLZ2Zmnn36Co0eP2KhK\nERFxRArqErI29NLh73VWL6tFi1aMGTOWpKSTPP74I3rKlohIGaagLiHbg5pcOPy9dz3k51u9vAED\nBtK7dyRbt27hxReH4oDPVhERERtQUJeQPGcjG2q3plJmMsbNm6xensFgYMKEKTRrFs63385mxoxP\nbVCliIg4GgV1CVp38fC32w8LbLI8Dw8PZs36Gn9/M6+++hJr1qyyyXJFRMRxKKhLUGxwGGnuJty/\nXwC5uTZZZrVq1Zk5MxonJycef/z/2LNnt02WKyIijkFBXYJynV1YU689TkkncVm1wmbLbdu2HRMn\nfsCZM6lERT3IyZMnbbZsERGxL6O9Cyhvlje4nXtif2bL6x8wedPVfyfNHNmlSMt9+OFHOHjwAO+9\n9w6PPtqHBQt+okKFCtaWKyIidqY96hIWX6UexysGclvCBtxybHtb1YsvvsRDDz3Mli0xPPPMAA0z\nKiJSBiioS5rBwIr6EXjknKdtwkYbL9rA5Mkf0qFDBIsXL+Lf/35et22JiJRyCmo7WNmgEwCd4m1/\nlbarqyuff/41jRuHER09i7Fjx9h8HSIiUnIU1HZw2K86ewNr0/yvrVQ8m2rz5Xt7V+Tbb/9L7dp1\nmDp1Mh98MMnm6xARkZKhoLaTFQ064WzJp8PutcWyfLPZzLx5C6lWrTpjx47h889nFMt6RESkeCmo\n7WRVvY7kGZzovGtFsa2jevUgvvtuIf7+/owY8Txz5nxdbOsSEZHioaC2k1RPX7bWaEa943sJPnWw\n2NZTu3Zd5s5diI+PD0OGPMPs2dHFti4REbE9BbUd/db4TgC67vi9WNfTuHET5s//EV9fX4YOHcRX\nX31RrOsTERHbUVDb0abarUj1qEjnXSsw5uYU67ouhXWlSpV4/vnBfPnlrGJdn4iI2EaRg/rTTz+l\nT58+9OzZk3nz5nHgwAEiIyOJiopi9OjR5F98lOPcuXPp2bMnvXv3Zvny5TYrvCzIdXZhecNOVDyX\nRpt91j9R62YaNWrM/Pk/4u/vz/DhQ/j00/8U+zpFRMQ6RQrqjRs3snXrVr755huio6M5fvw4b7/9\nNkOHDmX27NlYLBaWLl1KUlIS0dHRzJkzhxkzZjBp0iSys7NtvQ2l2m+NuwLQ9c/iPfx9ScOGjViw\n4CcCAyvz6qsv8dZbb2hQFBERB1akoF6zZg2hoaEMGjSIp59+mttvv524uDhat24NQEREBOvWrWP7\n9u2Eh4fj6uqKyWQiODiY+Ph4m25AaXeoUhC7qtQj/MA2zGlJJbLO+vUb8OOPv1KrVghTprzH8OFD\nNNyoiIiDKtJDOVJSUjh69CiffPIJhw8fZuDAgVgsFgwGAwCenp6kp6eTkZGByWQq+JynpycZGRk3\nXb6vbwWMRueilFYq/db4Thoc280dcUsxm58skXWazU1Yv34dd999N9HRn5OZmcbXX3+Nu7v7DT5j\nuu57UjjqofXUQ9tQH61XUj0sUlD7+PgQEhKCq6srISEhuLm5cfz48YL3MzMz8fb2xsvLi8zMzCum\nXx7c15OScrYoZZVaa+p1YMCKGdwZt5SkE2fAqWSu8XNyqsC8eT/w6KORLFiwgDvu6MqsWV/h6+t3\n1bxms4mkpPQSqausUg+tpx7ahvpoPVv38EahX6REaNGiBatXr8ZisXDixAnOnTvHbbfdxsaNFx4y\nsWrVKlq2bElYWBgxMTFkZWWRnp5OYmIioaGhRduKMuycqwer63UgMC0J12W/lei6vb0rMmfOAu69\n937WrVvD3Xffwb59CSVag4iIXF+Rgrpz5840aNCABx98kIEDB/Laa68xYsQIpk6dSp8+fcjJyaFb\nt26YzWb69u1LVFQU/fr1Y9iwYbi5udl6G8qExU3vBsB95vQSX7e7uzszZnzJ4MHD2Lcvkbvu6sLa\ntatLvA4REbmaweKAl/za+pBM/3eW2XR5xeXdb0ZQ//gekjdsJb9WiF1qmD07muHDh2AwGHjvvfeJ\njPw/QIfKbEE9tJ56aBvqo/Uc/tC3FI+fmt2DwWLBw44P0IiK6su8eQvx9PRkyJBneOWVf5OTU7yD\nsYiIyPUpqB3I2rrtyPc34/5NNJy13wV17dt35Oefl1GvXn2mT/+Enj27c+zYMbvVIyJSnimoHUiu\n0YVzjz6GU2oq7v/9zq61hITUYcmSZdx//wNs3LieFi1asHHjBrvWJCJSHhXp9iwpPoPS6zHD4MTJ\ntyYy9GA1uHhvOsDMkV1KtBYvLy+mT/+c5s1b8uabr/HAA/fw+uvjePLJpwvumRcRkeKlPWoHc9rk\nz4Y6baidtJ8GR+0/ipvBYOCZZwbz+++/4+PjwyuvjOCxxx4hJSXZ3qWJiJQLCmoH9GOzewHoEbPQ\nzpX8z+23387SpWto374jS5b8SOfO7Vm3bo29yxIRKfMU1A5oR/VG7A2sQ9uEjVRJcZyLuKpUqcp3\n3/3ASy+9yokTx+nZszvvvDOW3Nxce5cmIlJmKagdkcHAgpY9cMJCjy2Os1cN4OzszLBhL/LDDz9T\nvXoQkya9y/3336XRzEREiomC2kGtq3sbJ7wDuGPHMrzPnrF3OVdp1aoNy5atoWfPB/njj0107tye\nadM+KngOuYiI2IaC2kHlOzmzsPn9uOVlc0/sEnuXc03e3hX55JOZfPbZF1SoUIFRo0bSo8c97NuX\naO/SRETKDAW1A/ut8R2ku3nRfdtiXHOy7F3Odd1//wOsWrWJ7t3/yYYN6+jcuR3Tp3+svWsRERtQ\nUDuw864eLGl6FxXPpdFl53J7l3NDZrOZGTO+ZNq0WXh4ePDKKyO499472bHjT3uXJiJSqimoHdyi\n8HvJdnbhwc0LwMHH3DYYDPTo0YtVqzbxwAO9iIn5g65dI3jttZfJyMiwd3kiIqWSgtrBpXr68kuT\nfxCYdhL3eXPsXU6hBAQE8Omns/j22/8SFBTMJ598SIcOrVi8+Ecc8GFtIiIOTUFdCsxv9QA5zkYq\nTJ4Apeie5c6d72Dlyg08//y/SUo6yWOPRREZ2Ys9e3bbuzQRkVJDY32XAqdN/vzW+E7uif2Z6P6v\ns7xh56vmKelxwAvLw8ODkSNH0atXb0aOHM6yZb+zcuVy+vcfwPDhI/H19bN3iSIiDk171KXEd616\nkeNkpM+GeTjl59m7nFtWt24o3323kC+++IagoGCmT/+Etm3DmTHjU41sJiJyAwrqUiLJ28zSRl2o\nlnqUjrtL5xjbBoOBu+++l9WrNzFmzDhyc/N46aUX6dy5HUuW/KTz1yIi16CgLkXmte5FrpMzD2+Y\nWyr3qi9xc3PjmWcGs2HDVvr2fZy9e/fQr18k99xzJ2vWrLJ3eSIiDkVBXYqcrBjIb43vpHrKEe6M\nW2bvcqxmNpuZOPF9Vq3aSPfu/yQmZjM9e3bnoYf+ybZtW+xdnoiIQ1BQlzJz2vYhy+hK5Po5Dj1a\n2a0IDa3HzJnR/PrrCjp16szKlcv5xz9up3//vuzcGWfv8kRE7EpBXcoke/nxQ/h9+Gec5t5ti+1d\njk01a9acefMWMn/+Ipo3b8GPPy7k9ttvo1+/KGJjt9q7PBERu1BQl0LzWz1AhpsnD22aj+f5sjfi\nV8eOnViyZBlffz2XFi1asmTJj3Tt2omoqAfZvHmjvcsTESlRCupSKNPdi3mte2HKyqDX5v/au5xi\nYTAY6Nr1LhYvXsq8eQu57bb2/P77r9x7b1d69bqPFSuW6SpxESkXFNSl1I/N7uW0px/3b12Ef3qS\nvcspNgaDgU6dOrNw4RIWLlxCp06dWb16Jb1796Bz5/bMmfM12dnZ9i5TRKTYKKhLqWwXN6LbP4Jb\nbjb9Vkfbu5wScdtt7Zk3byG//rqCBx7oxe7du3juuYG0aNGY99+fSEpKsr1LFBGxOauC+vTp03Tq\n1InExEQOHDhAZGQkUVFRjB49uuBZxHPnzqVnz5707t2b5csd+1GNpc2yRp3ZG1ib2+NXYdy4wd7l\nlJhmzZrz6aez2LQplqeffpbMzEzGjXud8PCGjBz5Art3x9u7RBERmylyUOfk5PDaa6/h7u4OwNtv\nv83QoUOZPXs2FouFpUuXkpSURHR0NHPmzGHGjBlMmjRJhyltyGJwYlrnAQB4jRoBF/84Ki+CgoJ5\n44232LZtJ2PGjMPX14+ZM6fTsWNrHnjgXn744b/kOPijQUVEbqbIQT1+/HgefvhhAgICAIiLi6N1\n69YAREREsG7dOrZv3054eDiurq6YTCaCg4OJj9feji3FV63PivoRuMRuxe3b2fYuxy68vSvyzDOD\n2bx5O7NmfU1ERGfWrl3Nk0/2Izy8Ie+8M5ajR4/Yu0wRkSIp0tOzFixYgJ+fHx07dmTatGkAWCwW\nDAYDAJ6enqSnp5ORkYHJZCr4nKenJxkZN7+dyNe3Akajc1FKK5c+7/gotx/8A++3XofHHgFv72JZ\nj9lsuvlMdvbYY1E89lgUu3fv5pNPPmHWrFlMmvQu778/ke7du9O/f3/uvvtuXFxc7FJfaeiho1MP\nbUN9tF5J9bBIQT1//nwMBgPr169n165djBgxguTk/13Ik5mZibe3N15eXmRmZl4x/fLgvp6UlLNF\nKavcOm3yJ3PwMDzHj+PsiyPJHDve5uswm00kJaXbfLnFxc+vKi+//AZDhozg++/nM2vWZyxcuJCF\nCxdiNgfw0EMPExn5f9SrV7/EaiptPXRE6qFtqI/Ws3UPbxT6RTr0/fXXX/PVV18RHR1NgwYNGD9+\nPBEREWzceGEwilWrVtGyZUvCwsKIiYkhKyuL9PR0EhMTCQ0NLdpWyA2dHTSE3JDaeHz2KUaNk13A\n09OTRx55lN9/X8XSpWsYMOBpcnNz+OijD+jYsTV3392FL7+cRVraGXuXKiJyTTa7PWvEiBFMnTqV\nPn36kJOTQ7du3TCbzfTt25eoqCj69evHsGHDcHNzs9Uq5XLu7mRMmIIhPx+v4UNBz3i+SpMmYYwb\n9y7bt+/hs8++4I47urJ16xaGDx9C48Z1eeKJR/npp0WcP3/e3qWKiBQwWBxweCdbH5Lp/07pf9JU\nYQ39+X3u2Lmczzr1Z2GL+wumzxzZxarlltVDZceOHeXbb2czd+43JCTsBcBk8qZ79/t54IEH6dAh\nAqOxSGeIrlJWe1iS1EPbUB+t5/CHvsVxzez0OGnuJh5ZNxtzWtkdscxWqlSpytChw1m79g+WLl3N\noEFDMJlMfPPNV/Tu3YOmTevz8ssvsmnTxoKxAURESpKCuoxJ8/BmZqfH8Mg5z6DfPwbHO2DikAwG\nA02aNGX06DfZsiWOH374mX79niAvL5fPPvuU7t270rRpfUaMeJ5Vq1aQq1MLIlJCFNRl0NKGXdhS\noxkt/trCP3b8bu9ySh0nJyfatm3HhAmT+fPPvcyePY/IyP8jOzuLWbM+48EH76dRo9oMGfIMv/66\nhKyssvFccBFxTDpHXUb5pycx9cshOOXnM/jR9zlZMfC68xbm/LXOaUFubi7r16/lxx8Xsnjxj5w4\ncRwALy8Td97Zla5d76JLl65UqlTpmp9XD62nHtqG+mg9naMWq50ymZnWeQAVcs4z5JepGCw6v2ot\no9FIx46dGD9+ErGx8fz0028888xz+PlV4vvvFzBo0FM0bBjCPffcyZQp77Fjx596FKeIWE1BXYYt\nb3A7G2q3JuzwDrpv/cne5ZQpTk5OtGrVhjFjxrJ5cywrV25g1KjXad26LVu2/MFbb71Bly7tCQ9v\nyPDhQ/nllyVXDP4jIlJYOvRdxvlkpvLhl8/hkX2OF6Im8Je55lXz6NC3baWkJLN8+VJ+/fVnli//\nnZSUFABcXV1p3botnTp1plOnzjRp0hRnZw2Veyv0e2gb6qP1SvLQt4K6HGiVuInXFr7FIb/qDHvk\nPbJc3K94X0FdfHJzc4mJ+YPffvuZ1auXs3Xr1oL3fH196dChU0Fw16hR036FlhL6PbQN9dF6JRnU\nthnJQRza5tqtWdj8Pv65ZRFPLZvO1G6D7V1SuWE0GmnTpi1t2rTFbJ7Irl37Wb16BatWrWDlyuUs\nWvQ9ixZ9D0DNmrWIiOhM+/YdaNeuA4GBle1cvYg4AgV1OfF5h0dpdDiOf8QtZXtwGCsbdLJ3SeWS\nv78/DzzwIA888CAWi4V9+xJYsWI5K1cuZ+3a1Xz55Uy+/HImACEhtWnXrgNt27ajXbsOVK8eZOfq\nRcQedOi7HKmScowpXw0D4IVH3uOwX3VAh75Lys16mJuby7ZtW1i3bi3r169h48YNZGT8b/7g4Brc\ndlv7gq+aNWsVPFq2vNDvoW2oj9bTOWoFdbHpsHsNI356j8O+1Xg+agLn3CooqEvIrfYwNzeXuLg/\nC4J7w4Z1pKamXra8AFq2bE2rVm1o2bI1TZs2w8PDozhKdxj6PbQN9dF6OkctxWZNvQ7UPZ5Az5jv\nef7nKbx1/0h7lyTXYTQaado0nKZNwxk48Fny8/PZtWsn69evYf36dWzevJElS35kyZIfAXBxcaFx\n4yYFwd2yZWuqVate7va6RcoaBXU59EXHvoQk7aNt4iZ6b5wH3GnvkqQQnJycaNSoMY0aNebJJ5/G\nYrFw5Mhh/vhjU8HXn39uZ+vWLUyb9jEAlStXoUWLVoSHN78Y+s3w8fG185aIyK3Qoe9yyvtcGpO+\nfgFz2inSP/+a7Hu633B+HSqzXkn08Ny5c8TGbrsivE+ePHHFPDVr1qJZs3CaNm1Os2bhhIU1xWTy\nLta6bEW/h7ahPlpP56gV1CUi5OQ+xs95CTcXJ1K/X0xueIvrzqv/sa1njx5aLBaOHj3Ctm1biY3d\nyrZtW4iN3VowCMslderUpWnTcJo1C6dx4zAaNWrskHve+j20DfXRegpqBXWJaZW4iVcXvYPFrxIp\nS5aSf51BN/Q/tvUcpYcWi4WDBw+wbduWggCPjd1GenraFfNVrx508VB7k4tfjalZsxZOTvYbedhR\neljaqY/WU1ArqEvUbHMCppeGkxtaj9Qff8VyjT0p/Y9tPUfuYX5+Pvv3J7Jt21bi4nYQF/cncXE7\nrjps7unpRcOGjQoCvHHjJtSr1wAvL68SqdORe1iaqI/W01XfUqKikurQv8U/eSBmIce63MPonqPJ\nMboChbvHWko/JycnateuS+3adenVq3fB9JMnTxaE9oXvf7Jlyx9s3rzxis8HBQVTr1596tVrcPF7\nferWrVdiAS5SlimoBYBZEf0wpyXRYe86Rv74Lm/dN5I8Z/16lHcBAQEEBNxB5853FEw7f/48e/bE\ns2PHheCOj49n9+5d/P77r/z++69XfD44uAb16tUnNPRCeNev34C6devh6elZ0psiUmrpX2IBwGJw\nYuLdw6iQfZbW+/7g+Z+nMPHuYfYuSxyQu7s7YWHNCAtrdsX0lJRkdu+OJz5+F3v2xBe8/u23X/jt\nt1+umLd69SBq165DnTp1qV27DrVr16VOnbpUq1bdrufARRyRgloK5BpdeOv+l3h9wRgidq/hvNEN\nRt4B+odTCsHX14+2bdvRtm27K6YnJ58uCO3du3exe3c8iYkJrFx5YYzzy3l4eFCrVu2LIf6/AK9T\npy7e3hVLcnNEHIYuJpOrVMjKZOx3r1H3RCLnHn+SjLffwxxYURefWEkX8FwpIyOdxMQEEhL2kpiY\nQGLiXhISEkhMTODs2cyr5vf3N1OvXijVqgVTs2aty75CqFSpkkZguwX6XbServpWUNud6Vwa4+a9\nSq1TB/i18Z38586B5Ds5XzGPLjS7NfrHsXAsFgvHjx8jIWEvCQl72bcvoeD14cOHyMvLu+ozXl4m\natSo+bcAv/BVrVp1jEYdPLycfhetp6u+xe7SPbx55aE3eX3B6/xjx++45mYz+a4hV4W1iK0ZDAaq\nVKlKlSpV6djxysex+vi4s2VLHH/9tZ8DB/7ir7/2F3zt359IXNyfVy3PaDQSFBRMjRo1CQqqQVBQ\nEEFBwVSvHkxQUBCBgZVxdtbvtTguBbVcV7qHN6MefIPR/32T2+NX4ZqbzYR7XiDX6GLv0qSccnFx\nISSkNiEhta96z2KxcPLkyYvBva8gwC8F+ooV1z6y5uLiQtWq1QgKCr4Y4EFXvK5atRouLvqdF/sp\nUlDn5OTw8ssvc+TIEbKzsxk4cCB16tRh5MiRGAwG6taty+jRo3FycmLu3LnMmTMHo9HIwIED6dy5\ns623QYrRWTdPRvcczavfj6NdwgbG/PcN3r5vBJnuuj9WHIvBYCAwMJDAwEDatGl71fsZGRkcOXKY\nQ4cOcOjQIQ4dOsjhwwcLXq9Zs+qay3VycqJKlapUrx5EtWrVqVatOlWrVqVKlWoF381ms65Wl2JT\npKD+4Ycf8PHxYcKECaSmptKjRw/q16/P0KFDadOmDa+99hpLly6lWbNmREdHM3/+fLKysoiKiqJ9\n+/a4urraejukGJ139eD1B17lhSWTaZewgfHfvsyYnq/auyyRW+Ll5VUwGMu1nD9/niNHDnHo0CEO\nHz70t0A/xObNG9m4cf01P+vi4kLlylWoUqUqVatWpWrVK8O8atVqBAQE6hC7FEmRgvquu+6iW7du\nwIXDTc7OzsTFxdG6dWsAIiIiWLt2LU5OToSHh+Pq6oqrqyvBwcHEx8cTFhZmuy2QEpHt4sb47i/y\nxMpZ3L/1R96bPQLDww3Ia6L/llI2uLu7F4zOdi05OTmcOHGco0ePcvToYY4ePcqxY0cu/nyEY8eO\n8scfm8jPz7/m552dnQkMrEyVKlWpXLkKlStXJjDw0lcgAQEXXleqVEl753KFIgX1pVGFMjIyeO65\n5xg6dCjjx48vuD3C09OT9PR0MjIyMJlMV3wuIyPjpsv39a2A0ai/PB1NvpMz0zs/yUnvAJ5cORP+\neRd88QX07Gnv0kqNG13ZKYVjzx5WrepHeHjD676fm5vLiRMnLu6VH77mV2zsVmJiNl93GUajkcqV\nK1OlSpVrfl3YU69CQECAVVez63fReiXVwyL/Vz527BiDBg0iKiqK++67jwkTJhS8l5mZibe3N15e\nXmRmZl4x/fLgvp6UlLNFLUtKwMIW95Nk8mfksg8x9OpF5vMvcvbfr2hglJvQLTHWKw09dHX1pnbt\nRtSu3eia7+fl5XH69GlOnDjOyZPHOXHiBCdOHOfEieMcP3784vQTxMbGsnnz9QPdYDDg728mMLAy\nZrMZf/8LX2ZzAP7+/gQEBBT8XKmS/xWnHEtDHx2dw9+ederUKfr3789rr73GbbfdBkDDhg3ZuHEj\nbdq0YdWqVbRt25awsDCmTJlCVlYW2dnZJCYmEhoaWrStEIeyLrQdKS/0oGK/KDwnTcC440/SP5qO\nRaNHidyQs7PzxTHUA4DrnzqyWCykpqZcEeQXXh+7Ytq+fYns2LH9puutWNGnINCrV6+KyeRzWbBf\n+G42++Pvb8Zk8tYAMg6kSAOejB07liVLlhASElIw7ZVXXmHs2LHk5OQQEhLC2LFjcXZ2Zu7cuXz7\n7bdYLBb+9a9/FZzbvhENeFJ6eJ1L59+L3yP8QCzHKgYy4d7h7K1cV4OhXIP2YqynHl7b2bNnOXUq\niVOnkkhKuvT95DWmJXH69Clu9s++m5sbfn6VCr4qVfK77Ge/y6b/bx4PD48S2lrHoJHJFNSlilN+\nHo+sm82DmxaQ7+TElx36ct/cKToU/jcKGeuph9bLy8vDySmbXbv2XRbkJzl16tQV4X769GmSk5NJ\nT08r1HIrVKhwVZhXqlQJX1+/q0Ld19cXHx9fPDw8Su2eu8Mf+ha5XL6TM9Ed+hIbHMYLS6bQf9Xn\nZEceIm3qp1gCAuxdnohcxtnZGbM5AIOhcHvA2dnZpKQkXwzuC1+nT58mJSW54PWF6Rd+TkjYe82x\n2q/F1dUVHx9ffH19qVjRpyDAfXx8rvj+9/crVvQpV7e6KajFZrYHN+W5vpMZ+vMHtFy+FL+I1mS8\n/SZOz7wAAAzjSURBVB5ZPXpBKf2rWaS8c3V1LbiNrLDOnTt3VbhfHuqpqSmkpqaSmppCSkoKSUkn\n2bt3z3VvbbsWb++KNwx1Hx8fvL298fauSMWKFS++vjDNzc2tKK2wGwW12NSZCj688cAoZvvvxXPs\nGLz/1Z+s7xeQ/u5kLIGB9i5PREqAh4cHHh7VqFq1WqE/k5+fT0ZGOikpKVcF+ZkzqVd8v/z9hIQ9\nnD17a3cKubu7/y3AL7w2ma6edingL59+6RblkqJz1FJsqqQc47lfP6TxkTjS3byYcfvjLGvYmRkv\n3Wnv0uxC51etpx7aRlnrY1ZWFqmpqVcE+ZkzqaSlnSEtLY0zZ8787XUqaWlppKWd4cyZM+Tk5NzS\n+oxGI++++y7/939P2mwbdI5a7OKYbxVe7v0md8f+zGOrv2ToL1P5x5+/4dw9QCOaiYjNuLm5FYzz\nfqssFgvnz58vCO0zZ1JJT0+7+PrMFYGenn7he0ZGBkFBQcWwJdemoJZiZTE4sbjZPWwOacUTK2fS\nfu968u6MYEnYXXzVPuqaD/fQrV0iUlIMBsPFQ/Uet3QeviSPSuj+GSkRSd5m3rlvBK/1HM0xn8p0\nj13MtJkD+WfMQlxys+1dnojI/7d398FRlHcAx793e2+5xLvkIAnIezDhRVCENPUlIVJlIm0pLc1I\nzDSMomNlGBVKGSij0RmtQmfqaJ1CBEeYQV4sQkdqHbSAJkKmiJQ4BChBC2hCEl4uJCTkZW93+8eF\nCyBxDOGyd+fvM3PzPPvsy/PbJcNvd2f32YgliVr0qQPD7+DJotdYmz0bq6HxWOkaStbMY8rhj7Hq\nmtnhCSFExJFELfpcwGZnS9ZMHp9TwtZJvyTx4nl+t/01Xlu3gLuryqEHr2gIIUSsk0QtTHMhzsOa\n3Id54pEV7Lj1JwzxV/OH9/9EUk4Wzr9thEDA7BCFEMJ08nqWiBg3N9SQ/9lWphz5BJuuUedN5e+T\nZrBr7BTaHMFRlKL5QbNYeyXGDHIMbww5jr3Xl0OIyhW1iBinkgbxl7wneXzOSv55+zR8zX7m7lrF\nmtWPMqf0LVLP15kdohBC9DlJ1CLinPGkUHLfb3n0sVWsv6sAVXHwq/3bWPXWXDyzC7Dv+hdo8uCZ\nEOKHQd6jFhHrfHwSm+4q4N2sX3NPVTnTD7zPqO0f4Nz+AWcT+rFr7BR23HoftUkDgei+LS6EEN2R\nRC0iXkCxUzoml9IxuaTXHeP+yp3k/reMBz97lwc/e5fKQWMpHZOL5extGP37mx2uEELcUJKoRVQ5\nNiCdYwPSeSv3Ee788t9MrdzB7d8cZFzNYYzxq1CzJ9M+YybtP/05hq+f2eEKIUSvSaIWUand7gxd\nZSc3neHuY+XkHN3NqNKPcZR+jHvh0xwcMo7PR2TyWVomf1z+G7NDFkKI6yKvZ4mYktJYzz1V5eRU\n7SG9/stQeyBtJB1T8+i4Pw/1x3eBy9XnsckrMb0nx/DGkOPYe335epYkahGzfM1+Mo/vJ/N/n3Nn\nbSWWiy0AGC4XamYWavZkOu6ZTOCOieBwhD0e+c+x9+QY3hhyHHuvLxO13PoWMcuf4OOj8VP5aPxU\nbAGVcTWHyDy+n/HfHCRtdxmO3WXEA4bbjZp1J+qP70Kd9CMCd0zE8CaaHb4QQgCSqMUPRMBmp2LY\nBCqGTQDA09rEuG8qGf9NJbdVH2ToJ7twfNJ15+Vr32CODhzFnY/8AvWOSWijRvfJVbcQQlxNErX4\nQWqK81CecTflGXcD4L14nozaKkbXVjGq9ijpdccY6q+G3+8EwLDZ0DJGE7h1HIFbx4dKeR1MCBFu\nkqiFABrdiewbmcW+kVkAWHWNwf5qlk2wYjvwH2yHD2I7chjb4UrYvCm0npaSipaegTYyHS09He2W\ndAK3ZKAPHgKKYtbuCCFiiCRqIa5Btyp83X8YhdVA8hDInYE1R2NAYx1pp48z4swJRpw5wdBzX5O8\nZzeOPZ9esX6HYqcm6WbqvAOo96Zy2pNCvTeFek8Kp72p/LX4Z+bsmBAi6kiiFuJ70q0Kp5IGcSpp\nELtHZYfaHWo7N5+vZbC/mkENNQz214TKEWdPXntb65PQhg5HHzwEbeBA9AED0VMHXFYOCD7QZrH0\n1e4JISKUJGoheqnD7uRE8nBOJA+/coZh4GltIrWxntSm06Q0nSa16TSpjcF6SuUhnF8c6Ha77YoD\nf0IS/gQfjXFemuI85Nx/G4avH7qvH3q//hj9Ouu+fhAfH94dFUKYQhK1EOFisdDk9tLk9nJsYMa3\n5xsG3tZGfM1+fM1++rU0BOst/svKBkafOopi6MF1Dn7UbXdGXBx6kg/D48HweNE9nlC9a9obbPN6\n0W/qrLvdGO54jPj44EAwchUvREQJe6LWdZ3nn3+eo0eP4nA4ePHFFxk2bFi4uxUi8lksNLoTaXQn\ncjwlrdvFrLpGQlszntYmPK1NeDvL4O8C3tZGPK0XuC3RwOr3Y62rxVJ1FIuu9zgkw2LBcMeD292V\nwC9L5IbbHSpxujCcTgynC5wOjM5pnE4MhxNcwdJwui6rd86/tI7dAXY7WK1ygiBEN8KeqHfs2EFH\nRwfvvPMOFRUVLFu2jJUrV4a7WyFihm5VQlfm35thEKe24W5vIb79IvHtLbjbW0hob8HdcfGytou4\n1DZcahtOtR1XoA1XRxsutR3X2SacgdO41Hbsevi//23Y7WC3Y9jsYFOCpd0ONltwns0GLieJFivY\n7Bg2G9jsYL80v7NNUYKJX1EwLtWtClgtXW0W61XLWbvaFAUjtI4VFOs12i5tx9J1knGpvOxnXKpj\n6WYZrrH8t7fz7e1au91G6NfZl8FV7QC+eBR/y5UnR1efKH3X9HfOo/t5EIznRvcZpmWNq3fm0jyL\nBfon0FfCnqj3799PTk4OABMmTKCysjLcXQohLBZaHXG0OuI41/3IhN+bTVODiVztTOJqG3ZNxRHo\nwK6p2ANqcFrr6Kp3lnat48rpQHA5R0BF0QMouoaia9h0DUUPYNM0FEND6dCwtXWg6K2XLRNA0TpL\no+d3DEQXn9kBRLulS2H+kj7pKuyJurm5mYSErjMPRVEIBALYbN13/V1jnl6Pf/x5xg3dnhBCCJHc\nR/1Yw91BQkICLS0toWld178zSQshhBCiS9gT9cSJEykrKwOgoqKCjIxrPP0qhBBCiGsK+2cuLz31\nXVVVhWEYvPTSS4wcOTKcXQohhBAxIyK/Ry2EEEKIoLDf+hZCCCHE9ZNELYQQQkSwmE7Uuq5TXFzM\nrFmzKCoq4uTJa38gQXRPVVUWLVpEYWEh+fn57Ny50+yQota5c+fIzc3lq6++MjuUqPTGG28wa9Ys\nZs6cyebNm80OJ+qoqsrChQspKCigsLBQ/g576IsvvqCoqAiAkydP8tBDD1FYWMhzzz2Hfh2jAPZE\nTCfqy0dFW7hwIcuWLTM7pKizbds2EhMT2bBhA2+++SYvvPCC2SFFJVVVKS4uxuVymR1KVNq7dy8H\nDhxg48aNrFu3jrq6OrNDijqlpaUEAgE2bdrEvHnzePXVV80OKWqsXr2aZ555hvb2dgBefvll5s+f\nz4YNGzAMI+wXMDGdqGVUtN574IEHePrppwEwDANFUUyOKDotX76cgoICUlJSzA4lKu3evZuMjAzm\nzZvHE088wb333mt2SFFnxIgRaJqGrus0NzfLeBY9MHToUF5//fXQ9KFDh8jKygJg8uTJlJeXh7X/\nmP6Xup5R0cSV4js/ndjc3MxTTz3F/PnzTY4o+mzduhWfz0dOTg6rVq0yO5yo1NDQwKlTpygpKaG6\nupq5c+eyfft2LPIhj+/N7XZTU1PDtGnTaGhooKSkxOyQokZeXh7V1dWhacMwQn978fHxXLhwIaz9\nx/QVtYyKdmPU1tYye/ZsZsyYwfTp080OJ+ps2bKF8vJyioqKOHLkCIsXL+bMmTNmhxVVEhMTyc7O\nxuFwkJaWhtPpxO/3mx1WVFm7di3Z2dl8+OGHvPfeeyxZsiR0K1f0jNXalTpbWlrweDzh7S+sWzeZ\njIrWe2fPnmXOnDksWrSI/Px8s8OJSuvXr+ftt99m3bp1jBkzhuXLl5Oc3FejBMeGSZMm8emnn2IY\nBvX19bS2tpKYmGh2WFHF4/Fw003B7yh4vV4CgQCaFv6vosWisWPHsnfvXgDKysrIzMwMa38xfXk5\ndepU9uzZQ0FBQWhUNNEzJSUlNDU1sWLFClasWAEEH6yQh6JEX5oyZQr79u0jPz8fwzAoLi6W5yV6\n6OGHH2bp0qUUFhaiqioLFizA7XabHVZUWrx4Mc8++yyvvPIKaWlp5OXlhbU/GZlMCCGEiGAxfetb\nCCGEiHaSqIUQQogIJolaCCGEiGCSqIUQQogIJolaCCGEiGCSqIUQQogIJolaCCGEiGCSqIUQQogI\n9n+9g+OEHs8kCgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10f9a4da0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "p = lambda x: np.exp(-x)  # our distribution\n",
    "g = lambda x: 1/(x+1)  # our proposal pdf (we're thus choosing M to be 1)\n",
    "invCDFg = lambda x: np.log(x +1) # generates our proposal using inverse sampling\n",
    "\n",
    "# domain limits\n",
    "xmin = 0 # the lower limit of our domain\n",
    "xmax = 10 # the upper limit of our domain\n",
    "\n",
    "# range limits for inverse sampling\n",
    "umin = invCDFg(xmin)\n",
    "umax = invCDFg(xmax)\n",
    "\n",
    "N = 10000 # the total of samples we wish to generate\n",
    "accepted = 0 # the number of accepted samples\n",
    "samples = np.zeros(N)\n",
    "count = 0 # the total count of proposals\n",
    "\n",
    "# generation loop\n",
    "while (accepted < N):\n",
    "    \n",
    "    # Sample from g using inverse sampling\n",
    "    u = np.random.uniform(umin, umax)\n",
    "    xproposal = np.exp(u) - 1\n",
    "    \n",
    "    # pick a uniform number on [0, 1)\n",
    "    y = np.random.uniform(0,1)\n",
    "    \n",
    "    # Do the accept/reject comparison\n",
    "    if y < p(xproposal)/g(xproposal):\n",
    "        samples[accepted] = xproposal\n",
    "        accepted += 1\n",
    "    \n",
    "    count +=1\n",
    "    \n",
    "print(\"Count\", count, \"Accepted\", accepted)\n",
    "\n",
    "# get the histogram info\n",
    "hinfo = np.histogram(samples,50)\n",
    "\n",
    "# plot the histogram\n",
    "plt.hist(samples,bins=50, label=u'Samples');\n",
    "\n",
    "# plot our (normalized) function\n",
    "xvals=np.linspace(xmin, xmax, 1000)\n",
    "plt.plot(xvals, hinfo[0][0]*p(xvals), 'r', label=u'p(x)')\n",
    "plt.plot(xvals, hinfo[0][0]*g(xvals), 'k', label=u'g(x)')\n",
    "\n",
    "\n",
    "\n",
    "# turn on the legend\n",
    "plt.legend();\n"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Raw Cell Format",
  "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": 1
}