RejectionSampling_CC_1.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "RejectionSampling_CC_1.ipynb",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oKCcyTumsloK",
        "colab_type": "text"
      },
      "source": [
        "# Rejection Sampling\n",
        "\n",
        "The idea behind rejection sampling is to sample from a distribution of interest by sampling uniformally under its desnity function. This is useful if your distribution of interest is computationally intractable. \n",
        "\n",
        "More formally, suppose sampling from distribution f(x) is difficult. There must exist another distribution g(x) from which we can easily sample. To do this, we can reject samples that are outside of f(x) and accept those that are within its boundaries.\n",
        "\n",
        "The algorithm is as follows:\n",
        "\n",
        "1. Draw from $$U \\sim Uniform(0, 1)$$\n",
        "2. Draw from $$X \\sim g(x)$$\n",
        "3. If $$ U \\leq \\frac{f(x)}{c * g(x)}$$\n",
        "then accept the sample, otherwise reject it. \n",
        "\n",
        "The \"c\" constant here helps adjust the candidate distribution g(x) so it can enevelop f(x)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "H_TmpmZouKdA",
        "colab_type": "text"
      },
      "source": [
        "## Imports\n",
        "Our imports include the scipy library that contains the functions that define the density functions of both the uniform and beta distributions."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Jx8TD6mGvYc9",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 71
        },
        "outputId": "ab50aa31-7b04-42fb-9c82-bc2013ff0ce2"
      },
      "source": [
        "import numpy as np\n",
        "from scipy.stats import beta as beta_dist\n",
        "from scipy.stats import uniform as uniform_dist\n",
        "from scipy.stats import uniform\n",
        "import matplotlib.pyplot as plt\n",
        "import seaborn as sns"
      ],
      "execution_count": 15,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "/usr/local/lib/python3.6/dist-packages/statsmodels/tools/_testing.py:19: FutureWarning: pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.\n",
            "  import pandas.util.testing as tm\n"
          ],
          "name": "stderr"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-uMHGXswuYSm",
        "colab_type": "text"
      },
      "source": [
        "## Code\n",
        "\n",
        "The code below uses the uniform distribution as our candidate g(x) to sample from the distribution of interest f(x). This code uses a variant of rejection sampling known as \"Empirical Supremum Rejection Sampling\", which interatively updates the constant C based on the ratio of f(x) and g(x)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kixLuEEC3Vzj",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "C = 2\n",
        "a = 0\n",
        "b = 1\n",
        "alpha = 2\n",
        "beta = 2\n",
        "\n",
        "samples = []\n",
        "\n",
        "for i in range(10000):\n",
        "    U = np.random.uniform(0, 1)\n",
        "    X = np.random.uniform(a, b)\n",
        "    interest_dist = beta_dist.pdf(X, alpha, beta)\n",
        "    candidate_dist = uniform_dist.pdf(X, a, b)\n",
        "    if U <= (interest_dist / (C * candidate_dist)):\n",
        "      samples.append(X)\n",
        "    C_next = max(C, interest_dist / candidate_dist)\n",
        "    C = C_next"
      ],
      "execution_count": 20,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "g2fP5To-88-D",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 283
        },
        "outputId": "c5a0fb9a-db08-47c5-f0fd-fcfc50244122"
      },
      "source": [
        "sns.distplot(samples)"
      ],
      "execution_count": 21,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.axes._subplots.AxesSubplot at 0x7fbe5b0b97f0>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 21
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD5CAYAAAA3Os7hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3xcV5338c9vZtS7rGJ1ucmxLHfFLSR2uhOCTSCA05OFeAkJ7MIuCyxLeQLsPsCyCzyEDQmkmnSS4CQODinGibvcuy3bqpbVe/FIM+f5Q3JWkWVrZI10p/zer5dembn3aOZ7I/mnO+eee44YY1BKKeX/bFYHUEop5R1a0JVSKkBoQVdKqQChBV0ppQKEFnSllAoQWtCVUipAOIZqICKPAzcBNcaYgvO0WQr8EggB6owxS4Z63aSkJJObmzussEopFex27NhRZ4xJHmzfkAUdeBL4DfD0YDtFJB74LbDMGFMmIimehMrNzaWoqMiTpkoppfqISOn59g3Z5WKM2QA0XKDJbcArxpiyvvY1w06olFJqxLzRh54HJIjIehHZISJ3eeE1lVJKDZMnXS6evMY84GogAtgsIluMMUcHNhSRVcAqgOzsbC+8tVJKqbO8cYZeAawzxrQbY+qADcCswRoaYx41xhQaYwqTkwft01dKKXWRvFHQ/wx8QkQcIhIJLAAOeeF1lVJKDYMnwxafA5YCSSJSAfyA3uGJGGMeMcYcEpG/AHsBN/B7Y8z+0YuslFJqMEMWdGPMrR60+Tnwc68kUkopdVH0TlGllAoQWtCVUipAeGPYolJB69mtZR61u22BDtNVo0/P0JVSKkBoQVdKqQChBV0ppQKEFnSllAoQWtCVUipAaEFXSqkAocMWlRohtzFUt3TR1NFN+5keABKjQkmKCSM2PMTidCqYaEFXahBDjS83xlBS30FRSQNHa9o+KuQDZcRHMCc7nhsKxpMQFToaUZX6iBZ0pYbBGMPh062sP1JDeWMn4SE28lJjyEuJISU2jKhQBwaobz9DVVMXeyuaeGNvFX87Wss/XpPHXYtyCLFrT6caHVrQlfJQdUsXb+6tori2jcSoUJbPSmdudgKhjnMLdGJUKFNSYrgiL5lTTZ3srWzmR28c5LltZfx65Rzy02MtOAIV6LSgKzWEbpeb9w7X8MGxWkIdNm6amcaCCeOw28Sj70+Pj+CfrsvjvcM1/Our+7j5txv58acL+Fxh1ignV8FGC7pSF3Cyrp1Xd1VQ1+ZkbnY8ywrSiA4b/j8bEeHqaam8mRXP157bxTdf3svBqha+98l8bB7+YVBqKFrQlU+zavKrrm4X6w6cZuvJBhIiQ7j3slympMRc9Ov1P44bZ6QB8MTGEvaUN3HznMyPzvZ1Ei81ElrQlerHGMOaPaf45TtHae3q4bJJ47g2f/yg/eQXyybCJ2ekERFq591DNZzpcbPy0myPu3CUOh8t6ErRW8g3H6/nF389yo7SRjLiI7h9QQ5ZiZGj8n4iwtWXpBLmsLN2XxWv7Kzgs/MyR+W9VPDwZE3Rx4GbgBpjTMEF2l0KbAZWGmNe9l5EpYZ2sV0zrV3dvHOomic3lrCnopnkmDB+dstMnD1ubDL6Z8yfmJyEs8fFO4dqCA+xc/uCbGQM3lcFJk/O0J8EfgM8fb4GImIHfgq87Z1YSl3YmW4XDR1OWjq76XC6cBuD24DDJoTYbYQ6bL3/tdsIcfRuc7kMu8ubON3cycGqVvaUN7H5eD1Ol5uccZH85OYCPjs3k/AQu8d/ILzhyqkpdHW7+bC4jkf+doL7l04as/dWgcWTRaI3iEjuEM2+CvwJuNQLmZQ6h8ttKK5p42BVC2UN7dS0nMFcxOv81ztHAbAJTEyO5q5FOdwwI405WfGWjTYREW4oGE9LVzc/W3eYiclRXD99vCVZlH8bcR+6iGQANwNXogVdeVlzZzcbi+vYWdZIh9NFmMNGzrhICjLiSIkJJy7cQVSYA5tNEHoLv9PlprvHjdNl6Ha5P3rusAvLCsaTFB3GlJQYIkLtVh/eR0SEz87NRET4+gu7eenLi5ieHmd1LOVnvHFR9JfAt4wx7qH6/kRkFbAKIDtbh2ep86tp6eK/3znKi9srMBjy0+OYkxXPlJRoHCO4df6qS1K9mNK7Quw2HrtzHise3sh9TxXx2oOXkRITbnUs5Ue8UdALgef7inkScKOI9BhjXhvY0BjzKPAoQGFh4cV8YlYBztnj5slNJ/n1u8U4e9wU5iZwxZTkoJnYKiU2nMfuKuRzj2xm1dM7eH7VQsJDfOeThPJtIx5ca4yZYIzJNcbkAi8DXxmsmCs1lGPVrdz82438+9rDzJ+QyNtfv4IVszOCppifVZARxy9XzmZ3eRP/8vJejNFzH+UZT4YtPgcsBZJEpAL4ARACYIx5ZFTTqaBgjOGZLaX85M1DRIU5eOSOeSwr6L0ouOl4vcXprHH99PH8y7Kp/OwvR5g6PoYHrpxsdSTlBzwZ5XKrpy9mjLlnRGlU0Ol0uvj2K3v58+5TXDk1mZ/eMlP7jfvcv2QSR0638p9vH2FqagzX5Ptu/7/yDXqnqLJMZVMnX3qqiMOnW/jm9VO5f8mkoJ+oauD497nZCWwvaeCBZ3dy/5JJpMT2/rHTOV/UYHSmfWWJA6eaufnhjVQ0dvD4PZfywJWTg76YDybEbuOOBb2LYjyzpZROp8vqSMqH6Rm68hpP767MTozk758pIi4ihD/dv5i81IufxTAYxEeGcvuCbH7/wUme317GXYtyrY6kfJQWdDWmDlW18L0/7yc5Oow7F+VSVNJIUUnjmL3/WN7S700546JYPjudV3dVsu7Aae5clGN1JOWDtMtFjZl9lc38cWspaXHh3Hf5ROIiQqyO5FcuzU1k4cRxfFhcx592VFgdR/kgLehqTOyvbOaF7WVkJUbyd5dN8Knb7v3JJ2ekMSEpiu+8uo/d5U1Wx1E+Rgu6GnVHTrfwwvZyMhMiuWdxrt75OAJ2m3Db/GxSYsL4+2eKqG7psjqS8iFa0NWoOlHbxh+3ljE+Lpx7FucS5tBiPlJRYQ4eu6uQtq4e7n1iO21neqyOpHyEFnQ1aqpbuli9tZTEqFDu1TNzr5qWFstvbp/LkepWvvrsTnpcbqsjKR+go1zUqGjt6ubpzSWE2GzcvTiXyDD9VfOms6N1PjUzndd2V3LbY1tZMTv9nNWO9Aak4KL/ypTXOXvcPLOllLYzPdx3+UQSIoNrcq2xNH9CIg3tTjYcqyUhKpQleclWR1IW0oKuvMptDC8WlVPZ2MkdC3PITBidRZbV/7pueipNnU7WHThNQmQIMzPjrY6kLKIFXXnVW/uqOFjVwk0z05iWFmt1nKBg61vtqLmzm5d2VBAd7mBiUrTVsZQF9KKo8prtJxvYeLyexZPGsXhSktVxgkqI3cadC3NIjAxl9ZZSTutwxqCkBV15xZ7yJtbsPcWUlGhunJFmdZygFBnq4J7Lcgm123hy40maOpxWR1JjTLtc1Ig1tDu5f/UOYsIdfKEwC9sQa8uq0ZMQGcrdi3N5dMMJntxUQpjD7tFduToaJjDoGboaEZfb8LXndlHX7uT2+Tk6PNEHpMVFcMfCHOrbnDyzpZRuHaMeNLSgqxH5xdtH+LC4jh+vKCAjIcLqOKrPpORobinMpKS+nZeKynHruqRBYciCLiKPi0iNiOw/z/7bRWSviOwTkU0iMsv7MZUvWnfgNL9df5xb52fx+UuzrI6jBpiVGc+NBePZf6qFN/dW6WLTQcCTM/QngWUX2H8SWGKMmQH8CHjUC7mUjyur7+CfX9zDzMw4fvCp6VbHUefxiSnJfGJyEptP1PPBsTqr46hR5ski0RtEJPcC+zf1e7oFyBx5LOVLBi4K4XIbfrfhON1uN9fnj+eVnZUWJVOeWFYwnpaubv5y4DTxeuNRQPN2H/oXgbfOt1NEVolIkYgU1dbWevmt1Vh551A1FY2d3Dwnk4Qova3f19lEuGVuJjnjInl5RwUVjR1WR1KjxGsFXUSupLegf+t8bYwxjxpjCo0xhcnJOueEPzpe28aGo7UU5iQwIyPO6jjKQw67jdsX5BAd7mD1llJaOrutjqRGgVcKuojMBH4PrDDG1HvjNZXvaT/Tw4tF5SRFh3HTzHSr46hhig5zcOfCHLq63azeqsMZA9GIC7qIZAOvAHcaY46OPJLyRcYY/rSzgg6ni5Xzswh16IhXf5QWF8HnCzOpaOzk1V2VOvIlwAx5UVREngOWAkkiUgH8AAgBMMY8AnwfGAf8tm8u5h5jTOFoBVbW2HKinsOnW7lpZhppcTre3J/lp8dxbX4qfz1YTWpMGEumplgdSXmJJ6Ncbh1i/5eAL3ktkfI5ta1neGv/aaamxrBo4jir4ygvWJqXTHVLF28frCY9Xv9ABwr93KwuyOXu7WoJsdu4eW7GOSviKP8kInxmTibJMWG8WFSui00HCC3o6oKe2HiSsoYObpqZRmx4iNVxlBeFOmzcNj+bbpfhq8/u0nVJA4AWdHVeJ2rb+Pm6I0wbH8PsLL0ZJRClxIazYnY620oa+NW7x6yOo0ZIp8ZTg3K5Dd98eS9hDhsr5mhXSyCbk52AAR5+v5gleckU5iZaHUldJD1DV4N6YuNJdpQ28sPl07WrJQj84FP5ZCRE8I8v7Ka1S2868lda0NU5zna1XDMthZvnZFgdR42BmPAQfvmF2Zxq6uSHaw5aHUddJC3o6mOMMXz7lX2EOWz85OYZ2tUSROblJPLAlZP5084K3jtcbXUcdRG0oKuPeXVXJdtONvCdG6eRGhtudRw1xh68ajJ5qdH86yv7adGuF7+jBV19pLmjm39fe4jZWfF8oVAXrAhGYQ47P7tlFjWtXfzH2sNWx1HDpAVdfeQ/3z5CQ7uTH3+6AJtNu1qC1eyseL50+USe21bGlhM6154/0WGLQaz/whWnmjpZvaWUhRPHsbeimb0VzRYmU1b7+jV5rN1Xxfde28+bX7tcJ2PzE/pTUhhjWLuviohQO9dMS7U6jvIBEaF2/s/y6RyraePxjSetjqM8pAVdcaiqlRN17Vw9LZWIULvVcZSPuHpaKtfmp/Krd45R2dRpdRzlAS3oQa7H7eat/VUkx4QxX+8QVAP84FP5GAw/fkPHpvsDLehBbtvJBurbndxYkIZdL4SqATITIvnK0sm8tf80W/UCqc/Tgh7EzvS4eP9wDROTo5g6PsbqOMpH3Xf5RNLjwnnojYO43LrCkS/TUS5BbNPxetqdLq7LH291FOXDIkLtfOuGS/iH53fzLy/vYV7O0F1zty3IHoNkaiBPlqB7HLgJqDHGFAyyX4BfATcCHcA9xpid3g6qvKu5o5sPjtUybXwM2YmRVsdRFus/hHUwxhiyEiJ4+0A1BRlxhDn04rkv8qTL5Ulg2QX23wBM6ftaBfzPyGOp0fa7Dcc50+3mmnwdpqiGJiJ8ckYarWd62Fisfem+asiCbozZADRcoMkK4GnTawsQLyJp3gqovK+x3cmTm0qYkRmnCz4rj2WPiyI/LZYPjtXSdqbH6jhqEN64KJoBlPd7XtG3TfmoJzaepMPp4kpd7V0N03X5qTh73PztSI3VUdQgxnSUi4isEpEiESmqra0dy7dWfVq6unliUwnLpo/X2RTVsKXEhjM3J4EtJxto7HBaHUcN4I2CXgn0n5ovs2/bOYwxjxpjCo0xhcnJyV54azVcz2wupbWrhweunGx1FOWnrr4kBQHeP6xn6b7GGwV9DXCX9FoINBtjqrzwusrLOp0uHv/wJEunJjMjM87qOMpPxUeGUpibyM6yRhra9Szdl3gybPE5YCmQJCIVwA+AEABjzCPAWnqHLBbTO2zx3tEKqzxzviFom0/UU9/uJC8lZshhakpdyJK8ZLaXNLD+SA2fmZtpdRzVZ8iCboy5dYj9BnjAa4nUqHAbw6biOrISIsgZp+PO1cjERYRwaW4C2042cOXUFBKiQq2OpNBb/4PG4aoW6tudfGJKsq4TqrxiSV4KIsL6ozrAwVdoQQ8SHxbXER8ZQn5arNVRVICIiwihMCeBnaWNOuLFR2hBDwIVjR2U1HeweFKSzqiovGpJXu9otb8d0bN0X6AFPQh8WFxHmMNGYU6C1VFUgImPDGVebgI7Shtp0rN0y2lBD3CtXd0cqGxhbk4C4SE6oZLyvqVnz9K1L91yWtAD3PaSRlzGsHDCOKujqAAVHxnKvJwEikobae7stjpOUNOCHsBcbsP2kgYmp0STHBNmdRwVwJZMTcYYw8biOqujBDUt6AHsUFULzZ3denauRl1CZCgzM+PZVtJAp9NldZygpQU9gG05WU98RAiXpOnycmr0XT4lCWePm60ndb50q2hBD1C1rWc4UdvO/AmJ2PRGIjUG0uIimJISzabj9XR161m6FbSgB6jtJQ3YBObpUEU1hq7IS6btTA+v7hp0wlU1yrSgB6Ael5udZY1MS4slJjzE6jgqiExMiiIjPoLHNpzA5TZWxwk6WtAD0MGqFjqcLi7NHXp1dqW8SUS4Ii+ZE3Xt/PVgtdVxgo4W9AC0vaSB+MgQJqdEWx1FBaHp6bFkJ0byyN+O0zsZqxorWtADTGl9O8dr2ynM0Yuhyho2Ee67YiK7y5vYXtJodZygogU9wLywvRxBL4Yqa31uXibjokL53d+OWx0lqGhBDyAut+GVnZXkpcYQF6EXQ5V1wkPs3L04l3cP11Bc02Z1nKChBT2AfFhcx+mWLubq2bnyAbctyCbUbuPpzSVWRwkaHhV0EVkmIkdEpFhEvj3I/mwReV9EdonIXhG50ftR1VBeKionPjKEaeP1zlBlvaToMD41K52Xd1TQ0qWTdo2FIQu6iNiBh4EbgHzgVhHJH9Ds34AXjTFzgJXAb70dVF1Yc0c3bx+s5tOzM3DY9YOX8g33LM6lw+ni5aIKq6MEBU/+5c8Hio0xJ4wxTuB5YMWANgY4u7ZZHHDKexGVJ9bsqcTZ4+aWeboCu/IdMzLjKMxJ4KnNJbj1RqNR50lBzwDK+z2v6NvW3w+BO0SkAlgLfNUr6ZTHXt5RwSXjY5iermuGKt9y9+JcSus7WH+0xuooAc/hpde5FXjSGPMLEVkEPCMiBcYYd/9GIrIKWAWQnZ3tpbdWR6tb2VPRzPduykd07LnyAc9uLfvoscttiA138B9rD3O6+czH2t22QOuAN3lyhl4JZPV7ntm3rb8vAi8CGGM2A+FA0sAXMsY8aowpNMYUJicnX1xidY6Xispx2IRPz063OopS57DbhAUTx3Gspo2ali6r4wQ0Twr6dmCKiEwQkVB6L3quGdCmDLgaQESm0VvQdYHBMdDtcvPqrkqunpbCuGhdlUj5pktzE3HYhM0ndK700TRkQTfG9AAPAuuAQ/SOZjkgIg+JyPK+Zv8E3Ccie4DngHuMTuIwJtYfqaWuzcnn5mUN3Vgpi0SHOZiZGceusiZd0WgUedSHboxZS+/Fzv7bvt/v8UHgMu9GUwP175c8a/WWUqLCHFQ1dw26XylfsWhiEjvLmthR1sgnJp/TI6u8QAcs+7G2Mz0cPt3CnKx47Da9GKp8W0ZCBDmJkWw5UY9bP8CPCi3ofmxPeRNuoxNxKf8xf0IiDe1OTta1Wx0lIGlB91PGGHaUNpKZEEFqbLjVcZTySEFGHOEhNraXNFgdJSBpQfdTp5q7eifiytazc+U/Quw2ZmfFc+BUCx1neqyOE3C0oPupnaWNOGzCrMx4q6MoNSyX5ibicht2lTdZHSXgaEH3Qz0uN7vLm5iWFktEqN3qOEoNS1pcBJkJERSVNugSdV6mBd0PHTrdSme3Sy+GKr9VmJNIdcsZPUv3Mi3ofmhnaSOx4Q5dBFr5rVmZcYTabTy/Te+d8CYt6H6mubObo9WtzM1J0EWgld8KC7EzIzOO1/dU0aqLX3iNFnQ/s6usEQPM09Etys9dmptIZ7eL1/dUWR0lYGhB9yPGGIpKG5mQFKUTcSm/l5UQQV5qNC9s124Xb9GC7kdK6jtoaHfqxVAVEESElZdms6eimYOnWqyOExC0oPuRHaUNhDlsFKTHWR1FKa+4eU4GoXabnqV7iRZ0P9Ha1c2+ymZmZsYT6tAfmwoMCVGhXJufyut7q+h2uYf+BnVBWhn8xBt7q+h2GQq1u0UFmJvnZNDQ7mTDUV0TZ6S0oPuJF4vKSYkJIzMhwuooSnnVFXnJJESG8OqugStbquHSgu4HjlW3squsiXk5CboItAo4oQ4bN81M568Hq2nRMekjogXdD7zYtwj0HB17rgLUzXMzONPj5i/7T1sdxa95VNBFZJmIHBGRYhH59nnafF5EDorIARF51rsxg1e3y80rO3sXgY4O82jFQKX8zpyseHLHRfKadruMyJAVQkTswMPAtUAFsF1E1vStI3q2zRTgO8BlxphGEUkZrcDB5t1DNdS3O/l8YRbVLWesjqOUV/VfB3dScjTvHa7hf9YfJy4i5GPtbluQPdbR/JInZ+jzgWJjzAljjBN4HlgxoM19wMPGmEYAY0yNd2MGr+e2lZEaG8aSvGSroyg1qmZnxWPoXVpRXRxPCnoGUN7veUXftv7ygDwR2SgiW0RkmbcCBrOy+g42HKtl5aXZOOx6uUMFtnHRYWQnRrKrvFHnSb9I3qoSDmAKsBS4FXhMRM5ZSkdEVolIkYgU1dbqmNOhPLe9DAFWzs+yOopSY2J2VjzVLWc43dJldRS/5ElBrwT6V5TMvm39VQBrjDHdxpiTwFF6C/zHGGMeNcYUGmMKk5O1C+FCnD1uXioq56pLUkmL07HnKjjMzIjDJrC7TLtdLoYnBX07MEVEJohIKLASWDOgzWv0np0jIkn0dsGc8GLOoPP2wdPUtTm5faFeDFLBIzLMwdTUGPZUNOHWbpdhG7KgG2N6gAeBdcAh4EVjzAEReUhElvc1WwfUi8hB4H3gm8aY+tEKHQye3VpGRnwEV0zRTzIquMzOTqClq4fjtW1WR/E7Hg1sNsasBdYO2Pb9fo8N8I2+LzVCx2vb2HS8nm9ePxW7Te8MVcHlkvExhDls7C1vZkpKjNVx/IoOnfBBz20tw2ETPleYaXUUpcZciN3G9PRYDlQ16wyMw6QF3cd0dbt4eWcF108fT0pMuNVxlLLEzMx4urrdHK1utTqKX9GC7mPe2l9FU0e33hmngtqk5GiiQu3sqWi2Oopf0YLuY57dWsaEpCgWTRxndRSlLGO3CTMy4zhc1cKZbpfVcfyGzvbkA87OZ3G6pYvtJY3cUDCe57eXD/FdSgW2WZnxbDnRwMEqXW/UU3qG7kM2H6/HYRPm6TS5SpGVGEl8RAh7tdvFY1rQfUSn08Xu8kZmZ8UTqdPkKoVNhJmZcRyraaWx3Wl1HL+gBd1HFJU20O0yLJqkfedKnTUzMx63gbX7q6yO4he0oPsAtzFsOVHPhKQonbdFqX7S4sJJjg7jz7tPWR3FL2hB9wGHq1pp7OjWkS1KDSAizMqKY3tJA1XNnVbH8Xla0H3A5hN1xEWEMC0t1uooSvmcWZnxGANv7NFul6FoQbfY0epWjte2s3BCos7botQgxkWHMTMzjjV7tNtlKFrQLfbkphIcNqEwN9HqKEr5rOWz0tlX2czJunaro/g0LegWau7o5tWdlczKiidKhyoqdV43zUxHBNboxdEL0oJuoReLyunsdunFUKWGMD4unPm5ifx5T6WuN3oBWtAt4nIbnt5SwvzcRNLjdaiiUkNZPjudE7XtHDilUwGcjxZ0i7x3uIbyhk7uuSzX6ihK+YUbC9Jw2ITX9eLoeWlBt8hTm0pIiwvnuvxUq6Mo5RcSokK5Ii+Z1/ecwu3WbpfBeFTQRWSZiBwRkWIR+fYF2n1WRIyIFHovYuA5Vt3Kh8V13LEwB4dd/6Yq5anls9I51dzFjrJGq6P4pCGriYjYgYeBG4B84FYRyR+kXQzwD8BWb4cMNE9tLiHUYePW+bqIhVLDcW1+KuEhNh3tch6enB7OB4qNMSeMMU7geWDFIO1+BPwU6PJivoDT3NnNn3ZUsmJWOolRoVbHUcqvRIU5uHpaKmv3VdGj642ew5OCngH0X22hom/bR0RkLpBljHnzQi8kIqtEpEhEimpra4cdNhC81DdU8e7FuVZHUcovLZ+VTn27k43H662O4nNG3IErIjbgv4B/GqqtMeZRY0yhMaYwOTl5pG/td1xuw9ObS7k0N4GCjDir4yjll5ZOTSYm3KHdLoPwpKBXAln9nmf2bTsrBigA1otICbAQWKMXRs/1/uEayho6uGfxBKujKOW3whx2bigYz7oDp+nS9UY/xpP7zbcDU0RkAr2FfCVw29mdxphmIOnscxFZD/yzMabIu1H9z9m1Qs96fONJYsMdNLQ7z9mnlPLc8lkZvFhUwfuHa7hhRprVcXzGkGfoxpge4EFgHXAIeNEYc0BEHhKR5aMdMFDUtHRRXNPGwonjdFZFpUZo0aRxJEWH6QyMA3g0I5QxZi2wdsC275+n7dKRxwo8m0/U66yKSl2kwT7RTkmJ5q8Hq/nDByeJCLUDcNuC4B4KrHe1jIFOp4tdZU3MzIwnWmdVVMor5mTH0+M27D/VbHUUn6EFfQzsKGvE6XLrAtBKeVFGfATJMWHs1LtGP6IFfZSdXQA6JzGSDJ1VUSmvERHmZsVTWt9BfdsZq+P4BC3oo+zo6VYa2p16dq7UKJidnYAAu8qbrI7iE7Sgj7LNJ+qJDXcwPV1vJFLK2+IiQpiUHM2uskZd+AIt6KOqprWLYzVtzJ+gQxWVGi1zsuNp7OimpL7D6iiW04I+iracqMduE+ZP0KGKSo2W6elxhDps7CzVi6Na0EdJS1c3O0ubmJUZp0MVlRpFoQ4bszLj2FvZRGtXt9VxLKUFfZS8XFTRO1RxYtLQjZVSI1KYk0i3y/D6niqro1hKC/oocLsNT20uITsxkowEHaqo1GjLTIggNTaMF4rKh24cwLSgj4L1R2sore/QoYpKjRERoTAnkT3lTRw+3WJ1HMtoQR8FT2wsISUmjAIdqqjUmJmTFU+o3cbz24L3LF0Lupcdqmrhg2N13L04V4cqKjWGIsMcXDc9lVd2VtDpDM550rWge9nvPzhJRIid24N81pPDfxYAAAvFSURBVDelrHDnwhxaunp4PUin1dWC7kU1LV2s2VPJ5wsziY/UBaCVGmvzJyQyNTWGp7eUBOWdo1rQveipzSX0uA33XqZLzCllBRHhjkU57K9sYXcQzu+iBd1LOpw9rN5SxnX5qeQmRVkdR6mgdfOcDKLDHDyzudTqKGPOo4IuIstE5IiIFIvItwfZ/w0ROSgie0XkXRHJ8X5U3/anHRU0d3Zz3+UTrY6iVFCLDnPwmbkZvLG3Kuim1R2yoIuIHXgYuAHIB24VkfwBzXYBhcaYmcDLwM+8HdSXudyGP3x4kllZ8czLSbA6jlJB765FuThdblZvCa7F2D05Q58PFBtjThhjnMDzwIr+DYwx7xtjzk51tgXI9G5M3/buoWpK6ju47/IJiOhQRaWsNjklmqsuSeHpzSV0dQfPEEZPCnoG0H+kfkXftvP5IvDWSEL5m99/cJKM+AiWTR9vdRSlVJ/7Lp9IfbuTV3dVWh1lzHj1oqiI3AEUAj8/z/5VIlIkIkW1tbXefGvL7ClvYltJA/delovDrteYlfIVCycmUpARy2MfnMDtDo4hjJ7M61oJZPV7ntm37WNE5Brgu8ASY8ygVyKMMY8CjwIUFhb67f/hZ7f+b7/c6i2lhIfYsIl8bLtSyloiwn2XT+Qfnt/Ne4druCY/1epIo86TU8rtwBQRmSAiocBKYE3/BiIyB/gdsNwYU+P9mL6pqrmTg1UtLJ6URHiI3eo4SqkBbpyRRkZ8BL9dXxwUNxoNeYZujOkRkQeBdYAdeNwYc0BEHgKKjDFr6O1iiQZe6rsoWGaMWT6KuX3C+iO1hDlsXDZJ5zxXyhcM9im5MDeBP+8+xUOvH2RKagwAtwXo1BweLaVjjFkLrB2w7fv9Hl/j5Vw+r6ali/2VzVyRl0xEqJ6dK+Wr5uUksP5ILe8ermFySnRAj0TTq3gXaf3RWhx24bLJenaulC9z2GwsnZpMWUMHxTVtVscZVVrQL0JVcyd7yptYNHGcrheqlB+Yl51AXEQI7x6uCei+dC3oF2HdgdOEhdhYkpdidRSllAccdhtXTk2hrKGDA6cCd0UjLejDtOl4HUer21ial6J950r5kXk5CaTEhPGXA6c50xOYd49qQR8GYww/feswcREhul6oUn7GbhM+OSONhnYnT28KzJkYtaAPw6u7KtlT0cw101II0btClfI7U1JjyEuN5tfvHQvImRi1KnmoqcPJT948xJzseOZk64yKSvmrGwvS6HS6+MnaQ1ZH8Tot6B766V+O0NTZzU8+PQNbAI9jVSrQpcSG85Wlk3hlZyXvHa62Oo5XaUH3wI7SRp7bVsa9i3PJT4+1Oo5SaoQeuGoyeanR/Osr+2np6rY6jtdoQR9C25ke/vmlPaTHhfOP1+ZZHUcp5QVhDjs/u2UWNa1d/Oj1g1bH8Rot6EP4/mv7Ka1v57+/MFtvIlIqgMzOiuf+pZN4aUcFLxaVD/0NfkAL+gW8srOCV3ZV8rWrp7Bgog5TVCrQfP2aPC6bPI5/e20/+yqarY4zYlrQz2N3eRPffXU/CyYk8tWrplgdRyk1Chx2G79eOYekqFC+vHoHNa1dVkcaES3ogzhe28a9T2wjOSaM/3fbHOw2HdWiVKAaFx3GI3fOo7HDye2PbfXr8ela0AeobOrkrj9sw24Tnv67+aTEhFsdSSk1ymZmxvOHuy+lrKGDO/6wjaYOp9WRLooW9H72VTRz88Mbaens5ol75pObFGV1JKXUGFk0aRyP3VXI8Zo2PvM/m/xyql0t6H3W7qvi87/bTIjdxsv3L2ZGZpzVkZRSY+yKvGRWf2kBLZ3dfPrhjfxl/2mrIw1L0Bf0mpYu7l+9g6/8cSd542N49YHFTB0fY3UspZRF5k9IZM2Dn2BichRfXr2D+1fv4FRTp9WxPOLRwGoRWQb8it41RX9vjPm/A/aHAU8D84B64AvGmBLvRvWu6pYuntxUwurNpZxxufnm9VNZdcVEnXRLKUV6fAQvfXkRj204wW/eL2b9kVpunZ/NPYtzyR4XaXW88xqyoIuIHXgYuBaoALaLyBpjTP/bq74INBpjJovISuCnwBdGI/BI1LWd4f3DNfz1YDXvH6nB5TYsKxjPN6+/hAnaX66U6ifMYefBq6bw6TkZ/Oe6Izy9uYQnNp3ksklJXDMthasuSSUrMcKn1iiVoZZjEpFFwA+NMdf3Pf8OgDHmP/q1WdfXZrOIOIDTQLK5wIsXFhaaoqKiiwptjMHlNnS7DE6XG2ePm25X79eZHjfNnd00tjupbj1DRWMHJXXt7K9sobLvY1NaXDg3FKRx9+IccsYNv5APtrK4Usp/3LYge9jfU93SxR+3lPLGvipO1LYDkBgVyvT0WCYlR5OZEEFqbDixESHEhDuIDXcQEx5CeIidELtgtwkOmw2bMKI/AiKywxhTONg+T7pcMoD+98VWAAvO18YY0yMizcA4oG74cS/szb1VPPjcTjxdFjDUbiMzIYK5OQncvTiHRROTKMiI9am/qkop35caG843rpvKN66bysm6dj4srmN/RTP7TzWzq6yCtjM9Hr/W/Usn8a1ll3g945hOTiIiq4BVfU/bROTIWLzvMeB977xUEqPwR2qM6TH4Bj0GC93+vw8tOYZv/xS+ffHfnnO+HZ4U9Eogq9/zzL5tg7Wp6OtyiaP34ujHGGMeBR714D19kogUne+jjr/QY/ANegy+IRCOoT9PhnRsB6aIyAQRCQVWAmsGtFkD3N33+BbgvQv1nyullPK+Ic/Q+/rEHwTW0Tts8XFjzAEReQgoMsasAf4APCMixUADvUVfKaXUGPKoD90YsxZYO2Db9/s97gI+591oPslvu4v60WPwDXoMviEQjuEjQw5bVEop5R/0tkillAoQWtAHISLLROSIiBSLyDmji0QkTERe6Nu/VURyxz7lhXlwDN8QkYMisldE3hWR8w6FsspQx9Cv3WdFxIiIz41W8OQYROTzfT+LAyLy7FhnHIoHv0vZIvK+iOzq+3260Yqc5yMij4tIjYjsP89+EZFf9x3fXhGZO9YZvcYYo1/9vui98HscmAiEAnuA/AFtvgI80vd4JfCC1bkv4hiuBCL7Ht/vj8fQ1y4G2ABsAQqtzn0RP4cpwC4goe95itW5L+IYHgXu73ucD5RYnXtAviuAucD+8+y/EXgLEGAhsNXqzBf7pWfo55oPFBtjThhjnMDzwIoBbVYAT/U9fhm4Wnzr1tMhj8EY874xpqPv6RZ67y/wJZ78HAB+RO/cQb64dpgnx3Af8LAxphHAGFMzxhmH4skxGCC273EccGoM8w3JGLOB3tF357MCeNr02gLEi0ja2KTzLi3o5xpsqoOM87UxxvQAZ6c68BWeHEN/X6T3DMWXDHkMfR+Ns4wxb45lsGHw5OeQB+SJyEYR2dI3s6kv8eQYfgjcISIV9I6G++rYRPOa4f578Vljeuu/8j0icgdQCCyxOstwiIgN+C/gHoujjJSD3m6XpfR+StogIjOMMU2WphqeW4EnjTG/6JvM7xkRKTDGuK0OFmz0DP1cw5nqgAtNdWAhT44BEbkG+C6w3BjjayvjDnUMMUABsF5ESujt+1zjYxdGPfk5VABrjDHdxpiTwFF6C7yv8OQYvgi8CGCM2QyE0ztHir/w6N+LP9CCfq5AmOpgyGMQkTnA7+gt5r7WbwtDHIMxptkYk2SMyTXG5NJ7HWC5Mebi5mQeHZ78Lr1G79k5IpJEbxfMibEMOQRPjqEMuBpARKbRW9BrxzTlyKwB7uob7bIQaDbGVFkd6qJYfVXWF7/ovep9lN6r+9/t2/YQvQUDen9hXwKKgW3ARKszX8QxvANUA7v7vtZYnXm4xzCg7Xp8bJSLhz8Hobfr6CCwD1hpdeaLOIZ8YCO9I2B2A9dZnXlA/ueAKqCb3k9EXwS+DHy538/g4b7j2+eLv0eefumdokopFSC0y0UppQKEFnSllAoQWtCVUipAaEFXSqkAoQVdKaUChBZ0pZQKEFrQlVIqQGhBV0qpAPH/AS2+uftpOKIkAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    }
  ]
}