VAEs.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "VAEs.ipynb",
      "version": "0.3.2",
      "provenance": [],
      "collapsed_sections": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/robertmoni/modelbasedrl/blob/master/VAEs.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "qOt-sj9JSJEU",
        "colab_type": "code",
        "outputId": "db374f3d-20a6-4aa0-8497-42ded7eb0621",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 187
        }
      },
      "source": [
        "!pip install pyro-ppl pystan"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Requirement already satisfied: pyro-ppl in /usr/local/lib/python3.6/dist-packages (0.3.3)\n",
            "Requirement already satisfied: pystan in /usr/local/lib/python3.6/dist-packages (2.19.0.0)\n",
            "Requirement already satisfied: graphviz>=0.8 in /usr/local/lib/python3.6/dist-packages (from pyro-ppl) (0.10.1)\n",
            "Requirement already satisfied: six>=1.10.0 in /usr/local/lib/python3.6/dist-packages (from pyro-ppl) (1.12.0)\n",
            "Requirement already satisfied: torch>=1.1.0 in /usr/local/lib/python3.6/dist-packages (from pyro-ppl) (1.1.0)\n",
            "Requirement already satisfied: contextlib2 in /usr/local/lib/python3.6/dist-packages (from pyro-ppl) (0.5.5)\n",
            "Requirement already satisfied: tqdm>=4.31 in /usr/local/lib/python3.6/dist-packages (from pyro-ppl) (4.32.1)\n",
            "Requirement already satisfied: opt-einsum>=2.3.2 in /usr/local/lib/python3.6/dist-packages (from pyro-ppl) (2.3.2)\n",
            "Requirement already satisfied: numpy>=1.7 in /usr/local/lib/python3.6/dist-packages (from pyro-ppl) (1.16.3)\n",
            "Requirement already satisfied: Cython!=0.25.1,>=0.22 in /usr/local/lib/python3.6/dist-packages (from pystan) (0.29.7)\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "k6mZkupmSSRx",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import pyro\n",
        "import pyro.distributions as dist\n",
        "from pyro.infer.mcmc import MCMC, HMC, NUTS\n",
        "from pyro.optim import Adam\n",
        "from pyro.infer import SVI, Trace_ELBO\n",
        "import torch.distributions.constraints as constraints\n",
        "\n",
        "import pystan\n",
        "\n",
        "from statsmodels.distributions.empirical_distribution import ECDF\n",
        "from sklearn.neighbors import KernelDensity\n",
        "from sklearn.model_selection import GridSearchCV\n",
        "\n",
        "\n",
        "import numpy as np\n",
        "import scipy.stats as stats\n",
        "import pandas as pd\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "import pandas as pd\n",
        "\n",
        "import torch\n",
        "from torch import Tensor\n",
        "import torch.nn as nn\n",
        "import torch.nn.functional as F\n",
        "import torch.optim as optim\n",
        "from torch.utils.data import DataLoader, TensorDataset"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "iuXefrGrXSpa",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Setup some data\n",
        "theta = 0.6\n",
        "n = 1000\n",
        "y = stats.bernoulli.rvs(theta, size=n)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Xs1w4UJwYj49",
        "colab_type": "text"
      },
      "source": [
        "## Get MCMC samples for this model using Stan"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vJ5RL6khBrMg",
        "colab_type": "code",
        "outputId": "22be7502-4167-450d-efe5-fa2c45c6419e",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        }
      },
      "source": [
        "#Compile model\n",
        "\n",
        "model_code = \"\"\"\n",
        "data {\n",
        "    int<lower=0> n;\n",
        "    int y[n];\n",
        "}\n",
        "parameters {\n",
        "    real<lower=0, upper=1> theta;\n",
        "}\n",
        "model {\n",
        "    // likehood:\n",
        "    y ~ bernoulli(theta);\n",
        "    \n",
        "    // prior:\n",
        "    theta ~ beta(0.5, 0.5);\n",
        "}\n",
        "\"\"\"\n",
        "\n",
        "sm = pystan.StanModel(model_code=model_code)"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_cbff015f4398d041dc3dc90aad40531c NOW.\n"
          ],
          "name": "stderr"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "c_BJJVHhGJ_h",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Sample model\n",
        "data_dict = {'y': y, 'n': n}\n",
        "fit = sm.sampling(data=data_dict, iter=1000, chains=4)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "AxBYSv1QDwqB",
        "colab_type": "code",
        "outputId": "2aded813-3440-4a15-a8f3-cac6ce23108f",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 354
        }
      },
      "source": [
        "# Extract samples\n",
        "theta = fit.extract(permuted=True)['theta']\n",
        "\n",
        "# Print some statistics\n",
        "print(\"Some samples:\", theta[:10])\n",
        "print(\"Mean:\", np.mean(theta, axis=0))\n",
        "print(\"Standard deviation:\", np.std(theta, axis=0))\n",
        "\n",
        "# Prepare plots\n",
        "_, ax = plt.subplots(2, 2)\n",
        "\n",
        "# histograms\n",
        "# warning: for a caveat about using histograms see\n",
        "# https://stats.stackexchange.com/a/51753\n",
        "ax[0, 0].hist(theta, 15)\n",
        "ax[0, 1].hist(theta, 30)\n",
        "\n",
        "# Empirical cumulative distribution\n",
        "ecdf = ECDF(theta)\n",
        "ax[1, 0].plot(ecdf.x, ecdf.y)\n",
        "\n",
        "# Density estimation using KDE (with tuning parameter chosen by 3 fold CV)\n",
        "params_for_kde_cv = {'bandwidth': np.logspace(-2, 3, 10)}\n",
        "grid = GridSearchCV(KernelDensity(), params_for_kde_cv, cv=3)\n",
        "grid.fit(theta.reshape(-1, 1))\n",
        "x_kde = np.linspace(0, 1, 10000).reshape(-1, 1)\n",
        "y_kde = np.exp(grid.best_estimator_.score_samples(x_kde))\n",
        "ax[1, 1].plot(x_kde, y_kde)"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Some samples: [0.61295504 0.62972802 0.61487219 0.62385754 0.57044825 0.61062003\n",
            " 0.5997759  0.60374794 0.584482   0.56985944]\n",
            "Mean: 0.6007145269718976\n",
            "Standard deviation: 0.015428704817597625\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "[<matplotlib.lines.Line2D at 0x7f08cf6bcbe0>]"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 6
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAD8CAYAAACW/ATfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XmcXFWd9/HPr7vT3Ul3h87SBAhJ\nOkhCBJcAzSaMI5sGxsfghsEHCRjNo8LMMOozwjivWRhnBp8ZF3jwwYmABEUWUYcIOIqIgyIBQiQJ\nJEJ2s5F963R6qe7f88c91VQ61d1VXVVdXX2/79erX33r3Ft1zq069atzzz33XHN3RERkeCsrdgFE\nRKTwFOxFRGJAwV5EJAYU7EVEYkDBXkQkBhTsRURiQMFeYsvMJpnZ02a20sxeNbO/DOljzexJM1sd\n/o8J6WZmt5vZGjNbbmZnFHcPRDKnYC9xlgC+4O6nAucC15vZqcBNwFPuPg14KjwGuAyYFv7mA3cO\nfpFFBkbBXmLL3be5+9KwfBBYBUwEZgMLw2YLgSvC8mzgPo8sBurN7PhBLrbIgFQUuwAA48eP98bG\nxmIXQ4apl156aZe7N/S1jZk1AqcDzwMT3H1bWPUGMCEsTwQ2pTxtc0jbRi9Ut6WQMqnbSUMi2Dc2\nNrJkyZJiF0OGKTPb2M/6WuBHwI3ufsDMute5u5tZVnOKmNl8om4eJk+erLotBdNf3U6lbhyJNTMb\nQRTo73f3H4fk7cnumfB/R0jfAkxKefqJIe0I7r7A3ZvcvamhIaNGl0jBKdhLbFnUhL8bWOXuX09Z\ntQiYG5bnAo+mpF8TRuWcC+xP6e4RGdKGRDeOSJGcD3wCWGFmL4e0vwFuBR42s3nARuDKsO4J4HJg\nDdACXDe4xRUZOAX7EtB40+P9brPh1j8bhJIML+7+W8B6WX1xmu0duL6ghYqRdPVa9bhw1I0jIhID\natkXUSYtdhGRfFDLXkQkBtSyF5FBoSPZ4lKwF5G8G2hg7/k8nbDNH3XjiIjEgIK9iEgMKNiLiMSA\ngr2ISAwo2IuIxICCvYhIDCjYi4jEgIK9iEgMKNiLiMRAv8HezKrN7AUzW2Zmr5rZP4b0qWb2vJmt\nMbOHzKwypFeFx2vC+sbC7oKIiPQnk+kS2oCL3L053MLtt2b2M+DzwDfc/UEz+zYwD7gz/N/r7ieb\n2Rzgq8DHClR+ERnGNOd9/vTbsvdIc3g4Ivw5cBHwSEhfCFwRlmeHx4T1F1vqHZxFRGTQZdRnb2bl\n4bZtO4AngbXAPndPhE02AxPD8kRgE0BYvx8Yl+Y155vZEjNbsnPnztz2QkRE+pTRrJfu3gnMNLN6\n4CfAjFwzdvcFwAKApqYmz/X14k63LhSRvmQ1Gsfd9wFPA+cB9WaW/LE4EdgSlrcAkwDC+mOA3Xkp\nrYiIDEgmo3EaQoseMxsJXAqsIgr6HwmbzQUeDcuLwmPC+l+FGzWLiEiRZNKNczyw0MzKiX4cHnb3\nx8xsJfCgmX0F+D1wd9j+buB7ZrYG2APMKUC5RUQkC/0Ge3dfDpyeJn0dcHaa9Fbgo3kpnYiI5IWu\noBURiQEFexGRGNANx0WkpOiq2oFRsBeRnGRyjYcUn4J9gegLICJDifrsRURiQMFeRCQGFOxFRGJA\nffYikhWdjypNatmLiMSAgr2ISAyoG0dESl7PriVdZHU0texFRGJAwV5EJAYU7EVEYkDBXkQkBjK5\nLeEkM3vazFaa2atm9pchfayZPWlmq8P/MSHdzOx2M1tjZsvN7IxC74SIiPQtk9E4CeAL7r7UzOqA\nl8zsSeBa4Cl3v9XMbgJuAr4EXAZMC3/nAHeG/1JkmVwMo1EMIsNTJrcl3AZsC8sHzWwVMBGYDbwn\nbLYQ+DVRsJ8N3BduMr7YzOrN7PjwOiJDhpndA7wf2OHubwtpY4GHgEZgA3Clu+81MwNuAy4HWoBr\n3X1pMco9mHS17PCRVZ+9mTUS3Y/2eWBCSgB/A5gQlicCm1KetjmkiQw19wKzeqTdRHTEOg14KjyG\nI49Y5xMdsYqUjIyDvZnVAj8CbnT3A6nrQives8nYzOab2RIzW7Jz585sniqSF+7+DLCnR/JsoiNV\nwv8rUtLv88hioN7Mjh+ckorkLqMraM1sBFGgv9/dfxyStye7Z0Kl3xHStwCTUp5+Ykg7grsvABYA\nNDU1ZfVDUWw6tB3Wsj1iPap70szmE7X+mTx5cuFKKpKFTEbjGHA3sMrdv56yahEwNyzPBR5NSb8m\njMo5F9iv/nopRQM5Yg3PW+DuTe7e1NDQUICSiWQvk5b9+cAngBVm9nJI+xvgVuBhM5sHbASuDOue\nIDqJtYboRNZ1eS2xSGHldMQqMlRlMhrnt4D1svriNNs7cH2O5RIpluQR660cfcR6g5k9SDSUWEes\nUlI066XElpk9QDR8eLyZbQb+Hh2xyjClYC+x5e5X9bJKR6wlLt0girhfMKi5cUREYkAtexHppmHF\nw5da9iIiMaBgLyISAwr2IiIxoGAvIhIDCvYiIjGgYC8iEgMK9iIiMaBgLyISAwr2IiIxoGAvIhID\nCvYiIjGgYC8iEgP9ToRmZvcA7wd2uPvbQtpY4CGgEdgAXOnue8MtDG8jmve7BbjW3ZcWpugiIpnr\nOclb3KY8zqRlfy8wq0faTcBT7j4NeCo8BrgMmBb+5gN35qeYIiKSi0xuS/iMmTX2SJ5NdIcfgIXA\nr4EvhfT7wo0eFptZffJ+nvkqsIjkh6YzjpeBzmc/ISWAvwFMCMsTgU0p220OaQr2JSKTABC3w1+R\n4SDnE7ShFe/ZPs/M5pvZEjNbsnPnzlyLISIifRhosN9uZscDhP87QvoWYFLKdieGtKO4+wJ3b3L3\npoaGhgEWQ0REMjHQYL8ImBuW5wKPpqRfY5Fzgf3qrxcRKb5Mhl4+QHQydryZbQb+HrgVeNjM5gEb\ngSvD5k8QDbtcQzT08roClLmgdNJKRIajTEbjXNXLqovTbOvA9bkWSkSk0NI17Ibz4ANdQSsiEgMK\n9iIiMaBgLyISAwr2IiIxoGAvIhIDA50uQURKiIYUi1r2IiIxoGAvIhID6sYRGYbUbTMww/lCK7Xs\nRURiQMFeRCQGYtONo8NaEYmz2AR7yZ9MfziHS1+nyHCgbhwRkRhQsBcRiQF144iUOJ2PKqye72+p\ndk8q2IuIZKFUx+IXJNib2SzgNqAcuMvdby1EPklq2QxNmXwupfAlSaW6LaUq7332ZlYOfAu4DDgV\nuMrMTs13PiKDTXVbSlkhWvZnA2vcfR2AmT0IzAZWDuTF1LIZ3kqs9a+6LWmVQj0uRLCfCGxKebwZ\nOKfnRmY2H5gfHjab2Wv9vO54YFdeSpi9uOZd7PzHA7vsqzm/zpTciwKobivvHKTU43zmnXHdLtoJ\nWndfACzIdHszW+LuTQUskvIeYvkXe98HSnVbeQ/FvAsxzn4LMCnl8YkhTaTUqW5LySpEsH8RmGZm\nU82sEpgDLCpAPiKDTXVbSlbeu3HcPWFmNwA/Jxqedo+7v5qHl874sLgA4pp3sfMv9r4fQXVbeZdy\n3ubuxchXREQGkebGERGJAQV7EZEYKEqwN7NZZvaama0xs5vSrL/WzHaa2cvh71Mp6zpT0helpE81\ns+fDaz4UTqDlLW8zuzAl7WUzazWzK8K6e81sfcq6mQPd97DNlWa20sxeNbMfpKTPNbPV4W9uSvqZ\nZrYivObtZmb5zNvMZprZcyFtuZl9LGX7jPY9x/3O6TMfTHGt23Gt13nY98Gr2+4+qH9EJ7bWAicB\nlcAy4NQe21wL3NHL85t7SX8YmBOWvw18Nt95p2wzFtgDjAqP7wU+kqd9nwb8HhgTHh+bkue68H9M\nWE5u8wJwLmDAz4DL8pz3dGBaWD4B2AbUZ7rvueSd62euul34uh3Xel1qdbsYLfvuS87dvR1IXnI+\nYOEX/yLgkZC0ELiigHl/BPiZu7dk+bxM8v808C133wvg7jtC+vuAJ919T1j3JDDLzI4HRrv7Yo9q\nxn0MfN/T5u3ur7v76rC8FdgBNAzSfqeVxWc+mOJat+Nar3Pd97QKVbeLEezTXXI+Mc12Hw6HVo+Y\nWeqFLNVmtsTMFicPNYFxwD53T/TzmrnmnTQHeKBH2j+H53zDzKrSPCfT/KcD083s2bCPs/p57sSw\n3N8+5ZJ3NzM7m6gFszYlub99zzXvXD7zwRTXuh3Xep2P/Aetbg/VE7Q/BRrd/R1Ev/QLU9ZN8ehS\n448D3zSztwxi3oQWx9uJxlon3QzMAM4iOhz9Ug75VxAd9r0HuAr4jpnV5/B6ecs77Pv3gOvcvSsk\n52vf+8q70J/5YIpr3Y5rve4v/0Gr28UI9v1ecu7uu929LTy8CzgzZd2W8H8d8GvgdGA3UG9myYvE\neruMPae8gyuBn7h7R8pztnmkDfgu0aFdOplcbr8ZWOTuHe6+HnidqKL09twtYbmv18w1b8xsNPA4\n8GV3X5zlvueUd46f+WCKa92Oa73OOf9Brds9O/EL/Uf0K7cOmMqbJzRO67HN8SnLHwQWh+UxQFVY\nHg+sJpwMAX7IkSc0PpfPvFPSFgMXpnsO0YmkbwK35rDvs4CFKfu4ieiwbiywPrwHY8LyWE9/Iuvy\nPOddCTwF3Jjmdfvd9xzzzukzV90ufN2Oa70utbpdrC/F5US/bmuJflEBbgE+EJb/FXg1vHFPAzNC\n+ruAFSF9BTAv5TVPCpVjTXijqvKZd1jXSPQLW9bjNX8VyvMK8H2gNod9N+DrRHOkr0h+4GHdJ8P+\nrSE65EymN4W81wJ3EK6MzlfewNVAB/Byyt/MbPY9h7xz/sxVtwtftwf6+ZZ6vS6luq3pEkREYmCo\nnqAVEZE8UrAXEYkBBXsRkRgo2m0JU40fP94bGxuLXQwZpl566aVd7p7tlZF5obothZRN3c4q2JvZ\nPcD7gR3u/rY06w24jejsdAtwrbsv7e91GxsbWbJkSTZFEcmYmW0sVt6q21JI2dTtbLtx7iUaM9qb\ny4guFpgGzAfuzPL1RUSkALJq2bv7M2bW2Mcms4H7PBrPudjM6s3seHfflkMZZZhwdxJdTluii/bw\n15bopKW9k8MdnXR2OYnOMC64+zngePj/5us4QMo6gAumjaeqonzwd0yKas2Ogxw7uprR1SOKXZQh\nLd999r1NCnRUsDez+UStfyZPnpznYshg2tfSzrb9rexqbmPnwZS/5iOX9x/uoJCXdbz45UtoqFOw\nj5M39rdy2W2/4ZK3TuDOq3vO/iCpinaC1t0XEG6829TUpCu7Ssj6XYd4cuUbrNp2kLU7m1m+ef9R\n24wcUc6xo6toqK3i5GNrOe8t4zhm5AiqKsqorCijqqKcyooyKsvLGFlZzsjKckaUlVFWBmVmGGDh\nXhVmhMcQLaWmWfe6+lFq2cXN79buoqPT+dkrbxS7KENevoN9JpMCSQlqT3Tx/Prd/L+n1/Lcut0A\nnHBMNZPHjeLGS6YxfUIdDXVRcG+oq6KmakgM9JJhbsOuQ93L+w93cMxI/eD3Jt/fyEXADWb2IHAO\nsF/99aXtcHsn3/7vtdzz2/UcbEswvraSL1w6nQ+deSIT60cWu3gFEeZ5vw+YQHSqYIG732ZmY4GH\niOaR2QBc6eGGFFIcW/e3di9v2XtYwb4P2Q69fIBoTubxZrYZ+HtgBIC7fxt4gmjY5RqioZfX5bOw\nMnj2tbRzz2/Xc+/vNnCgNcEFJ4/nmvOmcP7J4+PQak8AX3D3pWZWB7xkZk8S3dbvKXe/Ndxr9CZy\nm+dccrRt/2HMohP5W/cd5tQTRhe7SENWtqNxrupnvQPX51QiKaqW9gR//+irPLJ0M+5w2duO4xPn\nTuHck8ZRVpb2fs/DTjga3RaWD5rZKqKBBrOJGjsQ3fjj1yjYF9W2/a28feIxLN+8nz2H2otdnCFt\n2DfRJHMrtx7g+h8sZcPuQ1z3rql8YOYJzJw0WDcTGprCUOPTgeeBCSndkm8QdfNIEe1r6WDmpPoo\n2Lco2PdFwV4AWPDMWr76X69RXmb8x9Vn8t7Tjit2kYrOzGqBHxHd3OJAcnQQREexZpZ2FJmGFQ+e\n5tYEx9ZVU1VRxl617PukidBizt256zfr+Jcn/sBFM47ldzddpEAPmNkIokB/v7v/OCRvD/crTd63\ndEe657r7AndvcvemhoaiTMkTC22JTto7u6irrmBsTaW6cfqhYB9jic4ubnjg93zl8VU0TRnD/73q\ndMbXVhW7WEUX5ni6G1jl7l9PWbUImBuW5wKPDnbZ5E3NrQkAaqsqGDOqkr3qxumTunFiyt35lyf+\nwOPLt/GXF0/jhotOZkS5fvuD84FPACvM7OWQ9jfArcDDZjYP2Eh0g24pkua2N4N9/agR7G3p6OcZ\n8aZgH1MPL9nEPc+u54OnT+TGS6aR2h8dd+7+W5KX6h7t4sEsi/TuYLJlX11BXXUFu3e1FLlEQ5ua\ncjG0evtBvvLYKmZOquffP/pOBXopScmWfV1VBXXVIzjYqpZ9XxTsY2bvoXY+cfcLVFaU8Y2PzaQ8\nJmPnZfhp7tGyT7b0JT0F+xhp7ejkU/ctYVdzG3defSZTx9cUu0giA3ao/c0++7qqCprbE3R1aU7F\n3ijYx0Rnl/PXjyznpY17+dqV7+TsqWOLXSSRnBzZZz8C9zd/AORoCvYx8aUfLWfRsq18/tLpzJ45\nsdjFEclZ6mic2uqKI9LkaAr2MfDdZ9fzyEub+eT5U/nzi04udnFE8qK5NUGZRfdOqAvBXv32vVOw\nH+ZeWL+Hf3liFRfNOJabL5+hkTcybDS3JaitqsDMqAu3JFSw752C/TDW3Jbgc/e/xHHHVPO1j75T\nF03JsHKwNdEd5Gurki17Db/sjS6qGsYefOGP7Gpu54efOY8xNZXFLo5IXjW3dXQH+dHqxumXmnrD\n1LJN+/j3X7zGn05voGnKmGIXRyTvmtsS3SdmaxXs+6VgPwy9smU/8xYuYXxtFV+7UlfIyvDU3NbZ\n3bJP/j+k0Ti9UrAfZtoTXfzFg7/nQGsH3/r4GZrFUoat5taO7hZ9TaWGXvZHwX6Y+fHSzazbeYjb\n55zOO2N+lykZ3prbEtSGIF9WZtRUlivY90HBfhjZtKeFWx5bySkT6rjkrccWuzgiBdXc+mafPUT9\n9urG6Z2C/TBy04+XU15m3H1tExUaZinDWGeXc6j9zT57gJqqCg4q2PdKEWGYeHjJJp5ds5u/umQ6\nJ44ZVeziiBRUcg6cupSWfV2VWvZ9UbAfBjbvbeEfF73K2Y1juea8KcUujkjBpd6SMKmmqqI7XY6m\nYF/iusJslgD/9tF3qPtGYiHZgj+iz76qQido+6DIUOK+t3gjv1u7m799/6lMGaf56SUeDrYd3bJX\nsO+bgn0Je337Qf7psZVceEoDc86aVOziiAyaZHdNnUbjZEzBvoR955l1jCgv42tXztRVshIrb85l\nP6I7rUYt+z4p2JeovYfaeXTZVj585kTGapIziZlky76mqrw7rbaqgo5Opy3RWaxiDWkK9iXI3fnf\njyyjq8u5+lyNvpH4SfbZ16W07JP99xqRk56CfQlatGwrv1y1g7+6dDozjhtd7OKIDLreWvYAh9rU\nsk9Hwb7EdHU5d/56LY3jRvHpPzmp2MURKYrmtg5Gjig/YqhxTfIGJm26gUk6CvYl5j9f3sIf3jjI\n5997CpUV+vgknprbOo8YYw9vjsxRyz49RYsS0trRydd+8Tpvn3gM73/78cUujkjRNLclqKs6Mtgn\nW/bNatmnpWBfQr733Ea27DvMTZfNoKxMQy0lvlLnsk/qPkGrln1aCvYlorWjkzueXsO7pzdw/snj\ni10ckaJqbksccfUsaDROfxTsS8Rjy7ex/3AHn3m3TsqKHGxNdHfbJNVW69aEfck62JvZLDN7zczW\nmNlNadZfa2Y7zezl8Pep/BQ1vtydHzy/kcZxozjvLeOKXRyRokvXZz9qRDQMU3Pap1fR/yZvMrNy\n4FvApcBm4EUzW+TuK3ts+pC735CnMsbegy9uYukf93HL7NM0LYIIoRunR599WZlRqznte5Vty/5s\nYI27r3P3duBBYHb+iyVJh9oSfPOXr3PG5Ho+oatlRXD36JaEVUe3VWuqytVn34tsg/1EYFPK480h\nracPm9lyM3vEzNJOx2hm881siZkt2blzZ5bFiI+Fz21g+4E2br78rWrVDxIzu8fMdpjZKylpY83s\nSTNbHf6PKWYZ46wt0UWiy49q2UOY5rhdwT6dQpyg/SnQ6O7vAJ4EFqbbyN0XuHuTuzc1NDQUoBil\n70BrB3f9Zj1/Or2BsxrHFrs4cXIvMKtH2k3AU+4+DXgqPJYiaO6eF6eXYK+WfVrZBvstQGpL/cSQ\n1s3dd7t7W3h4F3DmwIsXb/c+u4E9h9r5/KXTi12UWHH3Z4A9PZJn82bDZSFwxaAWSrp135IwXcte\nc9r3Kttg/yIwzcymmlklMAdYlLqBmaVe2vkBYFVuRYynHQdb+c5v1vGeUxp456T6YhdHYIK7bwvL\nbwATilmYOEs3l31STaXmtO9NVqNx3D1hZjcAPwfKgXvc/VUzuwVY4u6LgL8wsw8ACaLW0bV5LnMs\n/PvPX6O1o5ObL3trsYsiPbi7m5n3tt7M5gPzASZPnjxo5YqLg2luNp5UW61g35usgj2Auz8BPNEj\n7e9Slm8Gbs69aPH1hzcO8KOlW5hz1iROOa6u2MWRyHYzO97dt4Wj1x29bejuC4AFAE1NTb3+KMjA\nNKe5/2yS7kPbO11BOwR9/RevU11Rpr76oWURMDcszwUeLWJZYi050Vlvo3EOtSVw129sTwr2Q8wz\nr+/kFyu389n3vIVxtVXFLk4smdkDwHPAKWa22czmAbcCl5rZauCS8FiKIN3NxpNqum9N2DXYxRry\nsu7GkcLp6nL+6bGVTKwfyac1B07RuPtVvay6eFALImkd6KPPvi5lfpzqEeVHrY8zteyHkEeWbmb1\njmb+6tLpVFWoooqkc7A1QWV5WdpgXlOZnOZY/fY9KdgPEbua27jlpys5c8oYPnh6uouSRQSiPvt0\n/fXwZj++gv3RFOyHiCdWbKO5LcG/fujtlOvGJCK9OtiaSNtfD5rTvi8K9kPAvpZ2/uO/1/GWhhqm\nHVtb7OKIDGmZBPtDmh/nKDpBOwTc//wf2bLvMD/41Dma7EykH73NeAlv3of2oFr2R1HLvsjaE10s\n/N0G/mTaeN6l2w2K9OtAawd11UdPlQCpo3F0H9qeFOyL7J5n17PjYBvzLpha7KKIlIR0d6lKqum+\n6XjHYBapJCjYF9GmPS3c/tRqzphcz7unaZpnkUz01Wc/akQ5ZtCslv1RFOyLpLktwafvW0J5mXHb\nnNMp0wgckX65e9Sy76Ubp6zMopkv1Wd/FJ2gLYKuLufGB3/P6h3N3HvdWUwaO6rYRRIpCYc7Ouns\n5S5VSboPbXpq2RfBE69s45erdvDly9/Kn6j7RiRjB/uYFyeppqpcF1WloWA/yA61JfjaL17nxDEj\nueY83UBcJBt9zWWfVFs9QsE+DXXjDLL/819/YP2uQ3x/3jlUlOu3ViQb+w9Ho2xGj0zfZw9Qq5Z9\nWoo2g+jRl7dw3+KNzD1vChdM05h6kWzta2kHYMyoyl630U3H01OwHyS/fm0Hn394GedMHctfz5pR\n7OKIlKS9LVHLfsyo3lv2NbpbVVoK9oNg+eZ9XH//Uk6ZUMd3rmnqvvBDRLKTbNnX99Gyr1OwT0vB\nvsB2Hmxj3sIl1I+q5O5rm3odHywi/dtzqJ3yMmN0n6NxdGvCdBTsC+yLP1zGgcMdLLjmTI4/ZmSx\niyNS0va2dFA/ckSfEwbWVleQ6NKtCXtSsC+gRcu28t+v7+SL7z2F0044ptjFESl5+1raqe+jvx5S\n5rRXV84RFOwL5KlV2/niD5cxfUItn9B4epG82NvSztia3vvrIeXWhBqRcwQF+wJY+se9fPb7S5k8\ndhT3f+pc3fhYJE/2tXT0eXIWYExN1PLfG07mSkTBPs/W7mzmUwuXcHx9NT/8X+fRUFdV7CKJDBu7\nmtsYX9v3d6qhthqAHQfbBqNIJUNjAPNo9faDXHPPCxhw73VnM6afw00RyVxHZxe7mtuZMLrvYH9s\nWL9Twf4ICvZ5snXfYT505+9whx98+hymjq8pdpFEhpVk8J4wurrP7cbVVGKmln1PCvZ5sGlPC5+9\n/yXaE108esP5zDhudLGLJDLsbD/QCtBvy76ivIxxNZVq2fegYJ+jNTsO8pFvP0dnl3PHx89QoBcp\nkO0HouB9bF3fLXuAhrpqdh5sLXSRSoqCfQ6+v3gjtzy2kuqKMn7yufM5+djaYhdJZNh6Y/9hoP9u\nHICGuqruHweJaDTOAN31m3X87X++wjsmHsOTn/9TBXqRAtuwu4WaynLG1/Y/8GHy2JFs3H1IUyak\nUMs+Sy3tCb70oxX8dNlWLnnrsdx59ZmM0Lz0IgW3YfchGsfX9DlVQlLjuBoOtCbY29LR70VYcaFg\nn4U39rfyFw/+nhc37OGGC0/mzy8+WYFeZJCs33WIt03MbNqRkxpqwnOaGVsztpDFKhmKVBno6Ozi\nseVbufz237Bi836+ceVMvvi+U6iq0JWxIoOhuS3BH/e0MC3D7tKp46Pt1u44VMhilRS17Puw/3AH\nP1yyiXt+u56t+1uZcVwdd3z8DPXPiwyyFZv34w7vnFSf0faTx46irqqClzfv48qzJhW4dKVBwb6H\nlvYEz67ZzdOv7eAnS7dwuKOT0yfX83f/4zQunNGg1rxIETy/fjdmMPPEzIJ9eZlxZuMYXly/p8Al\nKx1ZB3szmwXcBpQDd7n7rT3WVwH3AWcCu4GPufuG3ItaGInOLhav28Oyzft4+g87WLZ5Hx2dTmVF\nGe9/x/H8z3Mmc8bkMRmdFJLhr7/6L/nn7vxsxRs0TRmT1RQkF5w8nq88voo1O5p1NE6Wwd7MyoFv\nAZcCm4EXzWyRu69M2WwesNfdTzazOcBXgY/lq8ADkejsYuOeFrYfaGXHgTa2H2hl094W1u86xEsb\n99LaEd3k4LQTRjPvgpM4/+RxnNU4VrNVyhEyrP+SZ4+v2MZr2w/y1Q+/PavnzZ45ka/+1x+441er\n+cbHZsa+wZZty/5sYI27rwO+HKyKAAAImklEQVQwsweB2UBqZZ8N/ENYfgS4w8zMBzDg9afLtrJ1\n32ESXU6i00l0ddHR6SQ6u6K0ri4SnU57ZxeH2ztpbkvQ0t7JobYEh9oTtLR1cqg90R3MU9VVV3BS\nQy2zTjuO894yjveddly/U6dK7GVS//uVPBcEkPqtcDxNWsqyH7ldz23fTEu/3tOmHb3tES+Z+loD\neH66Mqdb7K3Muw+189jyrbzjxGP40Bknko2Guio++56Tuf2p1Wzd18rpU+qpq6pgRHkZ5WWlE/hP\nHDOKWW87LufXyTbYTwQ2pTzeDJzT2zbunjCz/cA4YFfqRmY2H5gPMHny5LSZfe+5jbywYU/Kc2BE\nWRkV5UZFmVFRXkZFmTGivIxRleXUVFVQU1XOuJpR1FZVMDKkjRxRzgn11UwaO4oJo6uZMLq6+242\nIlnIpP73W7f3HmrnK4+vKlAR8y+1QWzdaXZUWuq2RponHfH81LSU17IjtxtZWc6HzziRL82aMaBh\nzjdePI36kSN44IU/cvdv1pPoKr2LrC48paEowT5v3H0BsACgqakp7Sfw3evOAqCi3BhRVkZZCf0a\nS3z1V7cnjR3Fin94b/fjdIGzv2CYKt226QJ0al7pAnTPsgwHZWXGJy+YyicvmIq709HpdHR20VVC\nV9bm6ygk22C/BUgdx3RiSEu3zWYzqwCOITpRm7Uatb5laMmk/vervMyoq+77PqqSf2ZGZYVRWRHP\ny4uy3esXgWlmNtXMKoE5wKIe2ywC5obljwC/Gkh/vcgQlEn9FxmSsmo6hz74G4CfEw09u8fdXzWz\nW4Al7r4IuBv4npmtAfYQfSFESl5v9b/IxRLJiA2FRreZ7QQ29rPZeHqc5B1Ecc272PnnK+8p7t6Q\nh9fJWj91u9ifbdJQKQcMnbIMlXJA32XJuG4PiWCfCTNb4u5Nyjs++Rd73wttqOzfUCkHDJ2yDJVy\nQP7KEs8zFSIiMaNgLyISA6UU7Bco79jlX+x9L7Shsn9DpRwwdMoyVMoBeSpLyfTZi4jIwJVSy15E\nRAaoKMHezGaZ2WtmtsbMbkqz/loz22lmL4e/T6Ws60xJX5SSPtXMng+v+VC46CVveZvZhSlpL5tZ\nq5ldEdbda2brU9bNHOi+h22uNLOVZvaqmf0gJX2uma0Of3NT0s80sxXhNW+3Xq55H2jeZjbTzJ4L\nacvN7GMp22e07znud06f+WDLoI5VhfKuCeVvTFl3c0h/zczeNwhl+Xx4z5eb2VNmNiVlXdr3vUDl\n6Os7n7beF7As30gpx+tmti9lXT7fk3vMbIeZvdLLegvf5zXh8zkjZV3274m7D+of0cUoa4GTgEpg\nGXBqj22uBe7o5fnNvaQ/DMwJy98GPpvvvFO2GUt0wdio8Phe4CN52vdpwO+BMeHxsSl5rgv/x4Tl\n5DYvAOcSTXnyM+CyPOc9HZgWlk8AtgH1me57Lnnn+pkP0fr9OeDbYXkO8FBYPjVsXwVMDa9TXuCy\nXJhSjz+bLEtf73uBypH2e9dXvS9UWXps/+dEF8/l9T0Jr/Vu4AzglV7WXx6+zxa+38/n8p4Uo2Xf\nPU2su7cDyWliByy0ZC8imlIZYCFwRQHz/gjwM3dvyfJ5meT/aeBb7r4XwN13hPT3AU+6+56w7klg\nlpkdD4x298Ue1YT7GPi+p83b3V9399VheSuwA8jmIqVc9jutLD7zwZbJvs4mKi9E5b847M9s4EF3\nb3P39cCa8HoFK4u7P51SjxcTzfeTb7l879LW+0Esy1XAAznk1yt3f4ao0dib2cB9HlkM1Ifv+4De\nk2IE+3TTxE5Ms92Hw6HLI2aWOvlUtZktMbPFFrpRiKZQ3ufuiX5eM9e8k+ZwdAX45/Ccb1h0t650\nMsl/OjDdzJ4N+zirn+dODMv97VMueXczs7OJWkRrU5L72/dc887lMx9smezrEdOAA8lpwDOtn/ks\nS6p5RC3JpHTveyHLke57V7T3JHRpTQV+lZKcr/ckE31957N+T4bqCdqfAo3u/g6iX62FKeumeHQ1\n2ceBb5rZWwYxb8Iv69uJ5kdJuhmYAZxFdGj1pRzyryDq0ngPUaviO2aW2Y03c9dn3mHfvwdc5+7J\nO8Lka9/7yrvQn3nsmdnVQBPwbynJg/m+9/m9K5I5wCPu3pmSVrJ1sRjBvt9pYt19t7u3hYd3Ed3P\nNrluS/i/Dvg1cDrRFMr1Fk2pnPY185F3cCXwE3fvSHnOtnCo1QZ8l94PvTOZInczsMjdO8Kh/OtE\nQbC3527hyEPvAe97H3ljZqOBx4Evh0PKbPY9p7xz/MwHWzbTgGNHTgOelymUsywLZnYJ8GXgAyl1\nv7f3vSDl6ON7V5T3JDjqCD6P70km+vrOZ/+e5HKCYSB/RC24dUSHR8kTJKf12Ob4lOUPAovD8hig\nKiyPB1YTTq4AP+TIk3Wfy2feKWmLgQvTPYfoRMo3gVtz2PdZwMKUfdxEdIg/Flgf3oMxYXls2K7n\nCdrL85x3JfAUcGOa1+1333PMO6fPfIjW7+s58gTtw2H5NI48QbuO3E7QZlKW04m65Kb1SO/1fS9Q\nOXr7zvda7wtVlrDdDGAD4VqkfL8nKa/ZSO8naP+MI0/QvpDLe1KsL8TlRC23tUQtRYBbiFoWAP8K\nvBo+iKeBGSH9XcCKkL4CmJfymicRBb01IQhU5TPvlA9mC1DW4zV/FcrzCvB9oDaHfTfg60T3NV1B\nCGZh3SfD/q0h6kpJpjeFvNcCd6RW0HzkDVwNdAAvp/zNzGbfc8g75898CNbv6lDeNaH8J6U898vh\nea+RZlRVAcryS2B7yue6qL/3fTC/833V+0KVJTz+B3o0XArwnjxANLKtg+jIdh7wGeAzKd+Jb4Vy\nrgCacnlPdAWtiEgMDNUTtCIikkcK9iIiMaBgLyISAwr2IiIxoGAvIhIDCvYiIjGgYC8iEgMK9iIi\nMfD/AT+2hBnhzuJYAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 4 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gHrplkFWY3ld",
        "colab_type": "text"
      },
      "source": [
        "## Get MCMC samples for this model using Pyro"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "CZZUnlYYYRuj",
        "colab_type": "code",
        "outputId": "d8c42569-ac1b-4960-b8cf-30984e882e14",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        }
      },
      "source": [
        "y_tensor = torch.as_tensor(y, dtype=torch.float32)\n",
        "\n",
        "def model(y_tensor):\n",
        "    prior_dist = dist.Beta(torch.Tensor([.5]), torch.Tensor([.5]))\n",
        "    theta = pyro.sample('theta', prior_dist)\n",
        "    with pyro.plate('observe_data'):\n",
        "        pyro.sample('obs', dist.Bernoulli(theta), obs=y_tensor)\n",
        "\n",
        "nuts_kernel = NUTS(model, adapt_step_size=True)\n",
        "mcmc_run = MCMC(nuts_kernel, num_samples=500, warmup_steps=300).run(y_tensor)\n",
        "posterior = mcmc_run.marginal('theta').empirical['theta']"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Sample: 100%|██████████| 800/800 [00:06<00:00, 120.08it/s, step size=1.07e+00, acc. rate=0.897]\n"
          ],
          "name": "stderr"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "liU25QagUqms",
        "colab_type": "code",
        "outputId": "d7e9d217-f7ef-4e94-95b5-73ba1a354d4f",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 85
        }
      },
      "source": [
        "print(\"Some samples:\", posterior.sample((10,)).reshape(-1))\n",
        "print(\"Mean:\", posterior.mean)\n",
        "print(\"Standard deviation:\", posterior.stddev)"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Some samples: tensor([0.5688, 0.5738, 0.5995, 0.5835, 0.5992, 0.5903, 0.5898, 0.6041, 0.5919,\n",
            "        0.6167])\n",
            "Mean: tensor([0.5998])\n",
            "Standard deviation: tensor([0.0151])\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OMfd4QjNnATb",
        "colab_type": "text"
      },
      "source": [
        "## Get replications (new instances of similar to data) from MCMC samples"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "o5BKWMvHnJK6",
        "colab_type": "code",
        "outputId": "2ae0c157-a00f-49be-e641-2726f9d7b34b",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 320
        }
      },
      "source": [
        "no_replications = 10000\n",
        "\n",
        "posterior_samples = posterior.sample((no_replications,))\n",
        "posterior_samples = posterior_samples.numpy().reshape(-1)\n",
        "\n",
        "replications = stats.bernoulli.rvs(posterior_samples)\n",
        "bins = np.arange(0, replications.max() + 1.5) - 0.5\n",
        "_, ax = plt.subplots()\n",
        "ax.hist(replications, bins)\n",
        "ax.set_xticks(bins + 0.5)"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "[<matplotlib.axis.XTick at 0x7f08cf736c88>,\n",
              " <matplotlib.axis.XTick at 0x7f08cf7cf550>,\n",
              " <matplotlib.axis.XTick at 0x7f08cf7cf940>]"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 9
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAD0tJREFUeJzt3W+snnV9x/H3RypuUWNBzpqmrSmJ\njQaXKOSkYFzMBrEUXFYeKMEtekKa9EldNFuywZ50A1nwiUyySdKMbsU4sUENjRLxpGLMEkEOwlCo\nrGcoaRugR09BHVFT992D8yu7V8/h3AdOzy39vV/Jnft3fa/vdd2/KyX9cP2576aqkCT15zWjnoAk\naTQMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnVo16Ai/lvPPOq40bN456GpL0\nqvLQQw/9uKrGFuv7rQ6AjRs3MjU1NeppSNKrSpKnhunzEpAkdcoAkKROGQCS1KmhAiDJ6iR3JflB\nkoNJ3p3k3CSTSQ6193Nab5LcmmQ6yaNJLhrYz0TrP5Rk4nQdlCRpccOeAXwa+FpVvR14J3AQuA44\nUFWbgANtGeAKYFN77QBuA0hyLrALuBjYDOw6GRqSpJW3aAAkeRPwXuB2gKr6VVU9B2wD9ra2vcBV\nbbwNuKPm3A+sTrIWuByYrKrZqjoOTAJbl/VoJElDG+YM4HxgBviXJA8n+eckrwfWVNXTrecZYE0b\nrwMOD2x/pNUWqkuSRmCYAFgFXATcVlUXAv/N/13uAaDm/l3JZfm3JZPsSDKVZGpmZmY5dilJmscw\nAXAEOFJVD7Tlu5gLhGfbpR3a+7G2/iiwYWD79a22UP3/qardVTVeVeNjY4t+kU2S9DIt+k3gqnom\nyeEkb6uqJ4DLgMfbawK4ub3f3TbZD3w0yZ3M3fB9vqqeTnIv8PcDN363ANcv7+FIsPG6r456Chqx\nH938/lFP4VVh2J+C+HPgc0nOBp4ErmXu7GFfku3AU8DVrfce4EpgGnih9VJVs0luBB5sfTdU1eyy\nHIUkacmGCoCqegQYn2fVZfP0FrBzgf3sAfYsZYKSpNPDbwJLUqcMAEnqlAEgSZ0yACSpUwaAJHXK\nAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwA\nSeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnq1FABkORHSb6X5JEkU612bpLJ\nJIfa+zmtniS3JplO8miSiwb2M9H6DyWZOD2HJEkaxlLOAP6oqt5VVeNt+TrgQFVtAg60ZYArgE3t\ntQO4DeYCA9gFXAxsBnadDA1J0sp7JZeAtgF723gvcNVA/Y6acz+wOsla4HJgsqpmq+o4MAlsfQWf\nL0l6BYYNgAK+nuShJDtabU1VPd3GzwBr2ngdcHhg2yOttlBdkjQCq4bs+4OqOprk94DJJD8YXFlV\nlaSWY0ItYHYAvOUtb1mOXUqS5jHUGUBVHW3vx4AvM3cN/9l2aYf2fqy1HwU2DGy+vtUWqp/6Wbur\naryqxsfGxpZ2NJKkoS0aAElen+SNJ8fAFuD7wH7g5JM8E8Ddbbwf+Eh7GugS4Pl2qeheYEuSc9rN\n3y2tJkkagWEuAa0BvpzkZP+/VdXXkjwI7EuyHXgKuLr13wNcCUwDLwDXAlTVbJIbgQdb3w1VNbts\nRyJJWpJFA6CqngTeOU/9J8Bl89QL2LnAvvYAe5Y+TUnScvObwJLUKQNAkjplAEhSpwwASeqUASBJ\nnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQp\nA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjo1dAAkOSvJw0m+0pbPT/JA\nkukkX0hydqu/ri1Pt/UbB/Zxfas/keTy5T4YSdLwlnIG8DHg4MDyJ4FbquqtwHFge6tvB463+i2t\njyQXANcA7wC2Ap9JctYrm74k6eVaNUxTkvXA+4GbgL9IEuBS4E9by17gb4HbgG1tDHAX8I+tfxtw\nZ1X9EvhhkmlgM/DtZTmSZuN1X13O3UnSGWvYM4B/AP4K+J+2/Gbguao60ZaPAOvaeB1wGKCtf771\nv1ifZxtJ0gpbNACS/DFwrKoeWoH5kGRHkqkkUzMzMyvxkZLUpWHOAN4D/EmSHwF3Mnfp59PA6iQn\nLyGtB4628VFgA0Bb/ybgJ4P1ebZ5UVXtrqrxqhofGxtb8gFJkoazaABU1fVVtb6qNjJ3E/cbVfVn\nwH3AB1rbBHB3G+9vy7T136iqavVr2lNC5wObgO8s25FIkpZkqJvAC/hr4M4knwAeBm5v9duBz7ab\nvLPMhQZV9ViSfcDjwAlgZ1X9+hV8viTpFVhSAFTVN4FvtvGTzD3Fc2rPL4APLrD9Tcw9SSRJGjG/\nCSxJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkA\nktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJ\nnTIAJKlTiwZAkt9J8p0k/5HksSR/1+rnJ3kgyXSSLyQ5u9Vf15an2/qNA/u6vtWfSHL56TooSdLi\nhjkD+CVwaVW9E3gXsDXJJcAngVuq6q3AcWB7698OHG/1W1ofSS4ArgHeAWwFPpPkrOU8GEnS8BYN\ngJrz87b42vYq4FLgrlbfC1zVxtvaMm39ZUnS6ndW1S+r6ofANLB5WY5CkrRkQ90DSHJWkkeAY8Ak\n8F/Ac1V1orUcAda18TrgMEBb/zzw5sH6PNsMftaOJFNJpmZmZpZ+RJKkoQwVAFX166p6F7Ceuf9r\nf/vpmlBV7a6q8aoaHxsbO10fI0ndW9JTQFX1HHAf8G5gdZJVbdV64GgbHwU2ALT1bwJ+MlifZxtJ\n0gob5imgsSSr2/h3gfcBB5kLgg+0tgng7jbe35Zp679RVdXq17SnhM4HNgHfWa4DkSQtzarFW1gL\n7G1P7LwG2FdVX0nyOHBnkk8ADwO3t/7bgc8mmQZmmXvyh6p6LMk+4HHgBLCzqn69vIcjSRrWogFQ\nVY8CF85Tf5J5nuKpql8AH1xgXzcBNy19mpKk5eY3gSWpUwaAJHXKAJCkThkAktQpA0CSOmUASFKn\nDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdcoA\nkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHVq0QBIsiHJfUkeT/JYko+1+rlJJpMc\nau/ntHqS3JpkOsmjSS4a2NdE6z+UZOL0HZYkaTHDnAGcAP6yqi4ALgF2JrkAuA44UFWbgANtGeAK\nYFN77QBug7nAAHYBFwObgV0nQ0OStPIWDYCqerqqvtvGPwMOAuuAbcDe1rYXuKqNtwF31Jz7gdVJ\n1gKXA5NVNVtVx4FJYOuyHo0kaWhLugeQZCNwIfAAsKaqnm6rngHWtPE64PDAZkdabaH6qZ+xI8lU\nkqmZmZmlTE+StARDB0CSNwBfBD5eVT8dXFdVBdRyTKiqdlfVeFWNj42NLccuJUnzGCoAkryWub/8\nP1dVX2rlZ9ulHdr7sVY/CmwY2Hx9qy1UlySNwDBPAQW4HThYVZ8aWLUfOPkkzwRw90D9I+1poEuA\n59ulonuBLUnOaTd/t7SaJGkEVg3R8x7gw8D3kjzSan8D3AzsS7IdeAq4uq27B7gSmAZeAK4FqKrZ\nJDcCD7a+G6pqdlmOQpK0ZIsGQFX9O5AFVl82T38BOxfY1x5gz1ImKEk6PfwmsCR1ygCQpE4ZAJLU\nKQNAkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0y\nACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1atEASLInybEk3x+o\nnZtkMsmh9n5OqyfJrUmmkzya5KKBbSZa/6EkE6fncCRJwxrmDOBfga2n1K4DDlTVJuBAWwa4AtjU\nXjuA22AuMIBdwMXAZmDXydCQJI3GogFQVd8CZk8pbwP2tvFe4KqB+h01535gdZK1wOXAZFXNVtVx\nYJLfDBVJ0gp6ufcA1lTV0238DLCmjdcBhwf6jrTaQnVJ0oi84pvAVVVALcNcAEiyI8lUkqmZmZnl\n2q0k6RQvNwCebZd2aO/HWv0osGGgb32rLVT/DVW1u6rGq2p8bGzsZU5PkrSYlxsA+4GTT/JMAHcP\n1D/Snga6BHi+XSq6F9iS5Jx283dLq0mSRmTVYg1JPg/8IXBekiPMPc1zM7AvyXbgKeDq1n4PcCUw\nDbwAXAtQVbNJbgQebH03VNWpN5YlSSto0QCoqg8tsOqyeXoL2LnAfvYAe5Y0O0nSaeM3gSWpUwaA\nJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhS\npwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHVq\nxQMgydYkTySZTnLdSn++JGnOigZAkrOAfwKuAC4APpTkgpWcgyRpzkqfAWwGpqvqyar6FXAnsG2F\n5yBJYuUDYB1weGD5SKtJklbYqlFP4FRJdgA72uLPkzwxyvm8Sp0H/HjUk9DIdP/nn0+OegYj97Zh\nmlY6AI4CGwaW17fai6pqN7B7JSd1pkkyVVXjo56HRsM/fyWZGqZvpS8BPQhsSnJ+krOBa4D9KzwH\nSRIrfAZQVSeSfBS4FzgL2FNVj63kHCRJc1b8HkBV3QPcs9Kf2xkvofXNP38N9d9Aqup0T0SS9FvI\nn4KQpE4ZAGcQf2ajb0n2JDmW5PujnotGI8mGJPcleTzJY0k+9pL9XgI6M7Sf2fhP4H3MfcHuQeBD\nVfX4SCemFZPkvcDPgTuq6vdHPR+tvCRrgbVV9d0kbwQeAq5a6O8BzwDOHP7MRueq6lvA7KjnodGp\nqqer6rtt/DPgIC/xawsGwJnDn9mQ9KIkG4ELgQcW6jEAJOkMk+QNwBeBj1fVTxfqMwDOHIv+zIak\nM1+S1zL3l//nqupLL9VrAJw5/JkNqXNJAtwOHKyqTy3WbwCcIarqBHDyZzYOAvv8mY2+JPk88G3g\nbUmOJNk+6jlpxb0H+DBwaZJH2uvKhZp9DFSSOuUZgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqU\nASBJnTIAJKlT/ws+BxGfDzPA0AAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pSFfXqFzY6rt",
        "colab_type": "text"
      },
      "source": [
        "## Get approximate Bayesian inference for Pyro and stochatisc variational inference"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "3OEK2MFJUvEx",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def model(y_tensor):\n",
        "    prior_dist = dist.Beta(torch.Tensor([.5]), torch.Tensor([.5]))\n",
        "    theta = pyro.sample('theta', prior_dist)\n",
        "    with pyro.plate('observe_data'):\n",
        "        pyro.sample('obs', dist.Bernoulli(theta), obs=y_tensor)\n",
        "\n",
        "def guide(y_tensor):\n",
        "    alpha = pyro.param(\"alpha\", torch.Tensor([1.0]),\n",
        "                       constraint=constraints.positive)\n",
        "    beta = pyro.param(\"beta\", torch.Tensor([1.0]),\n",
        "                       constraint=constraints.positive)\n",
        "    theta = pyro.sample('theta', dist.Beta(alpha, beta))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "yHyBUAE9ax1R",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# set up the optimizer\n",
        "pyro.clear_param_store()\n",
        "adam_params = {\"lr\": 0.2, \"betas\": (0.90, 0.999)}\n",
        "optimizer = Adam(adam_params)\n",
        "\n",
        "# setup the inference algorithm\n",
        "svi = SVI(model, guide, optimizer, loss=Trace_ELBO())\n",
        "\n",
        "n_steps = 100\n",
        "# do gradient steps\n",
        "for step in range(n_steps):\n",
        "    svi.step(y_tensor)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "2il2SgT5aumi",
        "colab_type": "code",
        "outputId": "96e78127-5320-41f6-97c7-31be5667d1ca",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 388
        }
      },
      "source": [
        "alpha = pyro.param(\"alpha\")\n",
        "beta = pyro.param(\"beta\")\n",
        "    \n",
        "inf_distribution = stats.beta(alpha.data.numpy(), beta.data.numpy())\n",
        "print(\"Some samples:\", inf_distribution.rvs(10))\n",
        "print(\"Mean:\", inf_distribution.mean())\n",
        "print(\"Standard deviation:\", inf_distribution.std())\n",
        "\n",
        "_, axes = plt.subplots(2)\n",
        "\n",
        "# Plot the posterior\n",
        "x_svi = np.linspace(0, 1, 10000)\n",
        "y_svi = inf_distribution.pdf(x_svi)\n",
        "axes[0].plot(x_svi, y_svi)\n",
        "\n",
        "# Plot replications\n",
        "posterior_samples = posterior.sample((no_replications,))\n",
        "posterior_samples = posterior_samples.numpy().reshape(-1)\n",
        "\n",
        "replications = stats.bernoulli.rvs(posterior_samples)\n",
        "bins = np.arange(0, replications.max() + 1.5) - 0.5\n",
        "axes[1].hist(replications, bins)\n",
        "axes[1].set_xticks(bins + 0.5)"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Some samples: [0.74527141 0.51609851 0.36864775 0.58797763 0.41954004 0.59100899\n",
            " 0.731861   0.52022881 0.56102338 0.64184099]\n",
            "Mean: [0.55838448]\n",
            "Standard deviation: [0.10022922]\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "[<matplotlib.axis.XTick at 0x7f08d05548d0>,\n",
              " <matplotlib.axis.XTick at 0x7f08d054d438>,\n",
              " <matplotlib.axis.XTick at 0x7f08d054dc50>]"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 12
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl4VeW5xuHfm5khjAmDQAgIMggi\nGEWl0gpWEaygttah6lErPVZtbR2qttYBPdZzqhXbastxrrYOaCvOVYTaVqEGEWQmzHOY54QM7/lj\nLziRErKBZK/srOe+Lq6s4dt7vR8Znr2mb5m7IyIi0ZMSdgEiIhIOBYCISEQpAEREIkoBICISUQoA\nEZGIUgCIiESUAkBEJKJCCQAzSzWz6Wb2ZhjbFxGR8PYAfgjMDWnbIiICpCV6g2bWERgB3A/8+GBt\nc3JyPD8/PxFliYg0GNOmTdvg7rk1tUt4AACPALcC2TU1zM/Pp7CwsO4rEhFpQMxsWTztEnoIyMzO\nAYrdfdpB2ow2s0IzK1y/fn0CqxMRiZZEnwMYBJxrZkuBF4EhZvZ81QbuPs7dC9y9IDe3xj0YERE5\nTAk9BOTutwO3A5jZ14Cb3f07iaxBJGxbd5cxeX4xUxZvZOG6HazbXgJAdmY6XXObUNC5JUN7taVT\nq8YhVyoNXRjnAEQiacmGnTw2qYg3Zq6mpKyS7Kw0erdvxgl5LUkxY9OuPUxfvoU3Z67h7jfmcFr3\nHL43+Gi+0j0n7NKlgQotANx9MjA5rO2LJMr2kjL+5735vDB1ORmpKZzXvyPfKuhIv44tSE2xf2u/\nZMNO3pyxmj9MWcZ3npzKkJ5tuOsbvencukkI1UtDZvX5gTAFBQWuq4AkmX2yaCM3vzKDNVt3c+nA\nzvxgaHdyszPjem1peQXPfryUX08sosKde0f24YIBHTD799AQqcrMprl7QU3tNBSESB1wd574+2Iu\neWIKGWkpjL/2VMaM6hP3H3+AzLRURg8+mvd+NJi+HZpz8yszuOPPX1BWUVmHlUuU6ByASC0rr6jk\np3+exUuFKzi7TzseurAfjTMO/1ftqBaN+OM1J/PQX+fz2ORFrNy8m8e/cwJNM/XrK0dGewAitais\nopIb/jSdlwpXcMOQbvz2kgFH9Md/r9QU49ZhPXnwgr58vGgjVz79L3aWltdCxRJlCgCRWrKnvJLv\nv/AZ78xay53n9OamM3uQcoCTvEfi2yfm8ehF/fls+RaufPpTdu1RCMjhUwCI1ILKSueW8TN4f846\n7h15LFd/pUudbWvEce155NvHU7hsEz/403QqKuvvhRxSvykARGrBg+/O4/XPV3PrsB5cfkp+nW/v\nG/2O4p5zj+WDucXc+8Zs6vPVfFJ/6SySyBF65p9L+P1Hi7ns5M5c+9WjE7bdy07JZ8Xm3Yz7aDFd\nc5tyxan5Cdu2NAzaAxA5Ah8v2sCYt+ZyRq+23H3usQm/Rv+2YT05o1dbxrw5h2nLNid025L8FAAi\nh2n1lt3c8Mfp5LduzK++3e+Ad/XWtZQU46EL+9GhZSOu/+NnbNxRmvAaJHkpAEQOQ0lZBdc+P43S\n8kp+f1kB2VnpodXSvFE6j106gI0793DjS59TqZPCEicFgMhheODtucxYuZVffqsf3do0Dbscjj2q\nOfeceyx/X7iBpz9eGnY5kiQUACKHaNK8Yp79ZBlXDerCsD7twi5nn4tO7MQZvdrw4LvzWLBue9jl\nSBJQAIgcgvXbS7ll/Ax6tsvm1mE9wi7nS8yMB84/juzMNH700ufsKdeYQXJwCgCROLk7t46fwfaS\nch69uD9Z6alhl/RvcrMzeeD8vsxevY1HJy4Muxyp5xQAInF6fsoyJs1fz09H9OKYttlhl1OtM49t\nx/kDOvC7vy1i3tptYZcj9ZgCQCQOKzbt4oF35jH4mFwuO7lz2OXU6GcjetOsUTq3vfqFhoqQaikA\nRGrg7tz+2hcY8MD5fZPigSytmmTw83N68/mKLTw/ZVnY5Ug9pQAQqcHLhSv4R9EGbh/eiw4tGoVd\nTtxGHn8Up3XP4b/fncfqLbvDLkfqIQWAyEGs21bCfW/NZWCXVlxyUl7Y5RwSM+P+UX2pcOeeN2aH\nXY7UQwoAkWq4Oz/98yzKKip58ILjan1s/0TIa92YG4Z0573Z6/howfqwy5F6RgEgUo23vljDB3PX\ncdPXe5Cf0yTscg7bd0/rQn7rxtz9xmzdGyBfogAQOYDtJWXc+8Yc+nZozlV1+HCXRMhMS+Xn3+jN\n4vU7eebjJWGXI/WIAkDkAH71/kLW7yjlvlF9Qhnls7YN6dmWIT3bMPaDhRRvKwm7HKknFAAi+5mz\nehvPfLyES07Ko1+nFmGXU2t+fk5vyiqcX7wzL+xSpJ5QAIhUUVnp3Pn6LFo0zuCWs+rXWD9HKj+n\nCd89rQuvTV/FtGWbwi5H6gEFgEgV46etZNqyzdx+dk9aNM4Iu5xad93p3WiTncmYN+fqOcKiABDZ\na/POPTzwzlwKOrfkggEdwy6nTjTJTOPms3rw+YotvDFzTdjlSMgUACKB/35vPttKyhkzqk9SXvMf\nrwsGdKRX+2Y8+M48Ssoqwi5HQqQAEAGmL9/Mi58u58pT8+nVvlnY5dSp1BTjZyN6sWrLbp7+59Kw\ny5EQKQAk8ioqnZ/9ZRZtsjO58evHhF1OQgzqlsPQnm14bFIRG/Qg+chSAEjkPT9lGbNXb+POc3rT\nNDMt7HIS5vbhvdhVVsEjHywIuxQJiQJAIq14ewm//Ot8vtIthxF924ddTkJ1a9OUSwfm8ad/raCo\nWM8QjiIFgETaA2/Po7SskntHHpsU4/zXth8O7U7j9FT+623dHBZFCgCJrCmLN/Ln6asYPbgrXXOb\nhl1OKFo3zeT6Id34cF4x/1i4IexyJMEUABJJZRWV3PmXWXRo0YjrTu8WdjmhuuLUfDq2bMR9b83R\n4yMjRgEgkfT0P5ewsHgHd597LI0yUsMuJ1RZ6ancdnZP5q3dzvhpK8IuRxJIASCRs2brbh75YCFD\ne7bh673bhl1OvTCib3sG5LXgl39dwI7S8rDLkQRJaACYWSczm2Rmc8xstpn9MJHbFwG47825VFQ6\nd597bNil1Btmxs/O6c367aX8/m+Lwi5HEiTRewDlwE3u3hs4GbjOzHonuAaJsI8WrOetL9Zw3end\n6NSqcdjl1CsD8lpybr+jGPfRYj1EPiISGgDuvsbdPwumtwNzgQ6JrEGiq7S8grsmzCa/dWNGD+4a\ndjn10q3DeuDA/7w3P+xSJAFCOwdgZvlAf2DqfstHm1mhmRWuX6+HWEvt+d+PFrNkw07uGdmHrPRo\nn/itTseWjfnuV7rw5+mrmLFiS9jlSB0LJQDMrCnwKnCju2+rus7dx7l7gbsX5ObmhlGeNEArNu3i\nN5OKOLtPO756jH6uDub7p3cjp2kG9701R88MaOASHgBmlk7sj/8L7v5aorcv0XTPG3NIMePOc3TK\nqSZNM9O46cwefLp0M+/OWht2OVKHEn0VkAFPAnPd/eFEblui64M56/hg7jp+MLQ7R7VoFHY5SeHC\ngk70bJfNA+/Mo7RczwxoqBK9BzAIuAwYYmafB/+GJ7gGiZDtJWXc+foserTN5qpBXcIuJ2mkphg/\nHdGL5Zt28dzHy8IuR+pIQse+dfd/ANEbcUtC89BfF7B2WwmPXTqAjDTd93goTuuey+k9cnn0w4Vc\ncEJHWjVpeM9Ijjr9RkiD9dnyzTz7yVKuOCWf/nktwy4nKd0xvBe79lQwVs8MaJAUANIg7Smv5PZX\nv6BdsyxuPqtH2OUkre5ts7nkpDyen7qcouIdYZcjtUwBIA3SuI8WMX/ddsaM7BOpp3zVhRvPiD0z\n4J43Zuuy0AZGASANzuL1O3j0wyJG9G3PGRrs7Yi1bprJj888hr8v3MB7s3VZaEOiAJAGpaLSuXX8\nTLLSUrjrG7rmv7ZcdnJnerbL5t435rBrj0YLbSgUANKgPPmPxRQu28zd5x5Lm2ZZYZfTYKSlpjBm\nVB9Wby3h1x8WhV2O1BIFgDQYC9dt55d/XcCZvdtyXn+NMVjbTsxvxQUDOvLE3xfrhHADoQCQBqG8\nopKbXplBk4xU7j+vbyQf8J4Itw/vSVZ6KndNmKUTwg2AAkAahMcnL2Lmyq3cf15fcrMzwy6nwcpp\nmsktZ/Xgn0UbeXPmmrDLkSOkAJCkN3PlFsZOXMg3+h3F8L7twy6nwbt0YGf6dmjOPW/MZvPOPWGX\nI0dAASBJbXtJGTf8aTptsjMZM1KPeEyE1BTjwQuOY8uuMsa8OSfscuQIKAAkabk7P/vLLFZs2sXY\ni/vTorHGqkmU3kc149qvHc1r01cxaX5x2OXIYVIASNIaP20lr3++mhvPOIYT81uFXU7kXD+kG93a\nNOWnr33B9pKysMuRw6AAkKRUVLyDuybM5uSurbju9G5hlxNJmWmpPHjBcazZVsKD784Luxw5DAoA\nSTrbS8r43h8KaZSeyiPf7k9qii75DMsJnVty1aAuPD9lOZN1KCjpKAAkqVRWOje9PIOlG3fxm0sG\n0K657vYN2y1n9aBH22xufmUmG3eUhl2OHAIFgCSVxyYX8dc567hjeC9OObp12OUIkJWeyiMXHc+2\n3WX85NWZukEsiSgAJGlMnLuOh95fwKjjj+KqQflhlyNV9GrfjJ+c3ZMP5hbzwtTlYZcjcVIASFKY\nuXIL1/9xOn2Oas4D5x+noR7qoStPzee07jmMeXMOs1dvDbsciYMCQOq9FZt2cdUzhbRqksGT/1FA\no4zUsEuSA0hJMR6+8HhaNE7n2uc/Y+suXRpa3ykApF7buquMK5/5lD3lFTxz5Ym0ydZJ3/osNzuT\nxy49gTVbd3PjS9OprNT5gPpMASD11vaSMi5/+l8s37iLcZcX0L1tdtglSRxO6NySn5/Tm0nz1zN2\n4sKwy5GD0MNSpV7aWVrOlU9/yuxVW3ns0gGc3FVX/CST75zcmc9XbGXsxIV0yWnCKD2foV5SAEi9\ns3tPBd99tpDPlm/m1xcP4Mxj24VdkhwiM+O/zu/Dys27uGX8DNo2y9Jlu/WQDgFJvbJ1dxmXPzWV\nKUs28tCF/RhxnIZ3TlaZaamMu6yAzq2b8L0/FLJg3fawS5L9KACk3ijeXsK3f/8Jn6/Ywm8uHsB5\n/TuGXZIcoeaN03nmyhPJTE/l0iemsni9HiVZnygApF4oKt7Ot373Ccs27uLJK07UJ/8GpGPLxvzx\nuwOprHQu+d+pLNu4M+ySJKAAkNBNmlfMeb/9mJ2l5bxwzUAGH5MbdklSy7q3zeb57w6kpLyCi8dN\n0Z5APaEAkNBUVjqPTS7iqmc/Ja91Y16//isMyGsZdllSR3q1b8bzVw+kpLySb/4udqhPwqUAkFAU\nbyvhiqf/xX+/O58Rfdsz/j9PpUOLRmGXJXWsT4fmvHrtqTTJTOXicVP4cN66sEuKNAWAJNy7s9Yy\nbOzf+XTpJh44vy+/vri/hneIkC45TXj12lM5uk0Trn62kEcnLtQdwyFRAEjCrN6ym2ueK+Q/n59G\nu2ZZvHnDV7j4pDwN7BZBbbKzeOV7pzLq+A48/P4CrnmukE0794RdVuToRjCpc7v2lPPk35fw+N8W\n4Q53DO/JlYO6kJ6qzx9R1igjlYcv7Ef/vBaMeXMOZ/7qb9x/Xl/O0o1/CaMAkDpTUlbBK4UrGDux\niA07Sjmzd1vuPKc3nVo1Drs0qSfMjMtPyefE/Fbc/MoMvveHaYw4rj13DO+lc0IJYPX56T0FBQVe\nWFgYdhlyiDbt3MPzU5bx3CdL2bBjDyd1acVPhvXkhM66wkeqV1ZRyeOTF/HbSUUAjB7clWsGd6VZ\nVnrIlSUfM5vm7gU1tlMASG2oqHT+UbSBV6et5L3Zayktr+RrPXK55rSunHp0ax3nl7it2rKbB9+Z\nx4QZq8nOTOOyUzpz5aAu5GZnhl1a0lAASJ0rKavg40Ub+HBeMe/PWce6baU0b5TON/q15/JT8jlG\nwzfLEZi1aiuPT17E27PWkJZifL13W75V0InB3XNJTdEHioOptwFgZsOAsUAq8IS7/6K6tgqA+mVH\naTnTl2+mcOlmpi3bTOGyTZSUVdI4I5XTuucw6vgODOnVhsw0XdIptWfx+h28MHU5f56+ik0799Cq\nSQZf65HL0J5tGdStNS0aZ4RdYr1TLwPAzFKBBcDXgZXAp8DF7j7nQO0VAIlXVlHJ2q0lrNlawsrN\nu1hYvIOF67azYN0OVmzehTuYQc92zRjYpRVDerZhYNdW+qMvdW5PeSUfzivmvdlrmTS/mC3BIyeP\nzm3CCZ1b0rdDc47ObUq3Nk3Jzc6M9GHHeAMg0VcBnQQUuftiADN7ERgJHDAADldFpVNSVkHVaNsb\ndL5vfu/E3i/+5eVfanvg1/r/v/gArznw++7/nvtvM57XHrQ/OKXllZSWV7In+FpaVrFvWWl5BTtK\nytm6u4ytu8vYsruMbbvL2LKrjHXbSli/o/RL9aSnGl1ymtC3Y3O+eUJHju/Ugv55LcjWiTlJsIy0\nFIb1acewPu0or6hk+oot/GvJJqYv38z7c9bxcuHKfW2bZqbRvnkWbZtl0aZZJm2ys2jWKI3srHSy\nM9NomplGk8w0MtJSyEhNIS3VSE9NIT34mpZqpAQBYrAvTGLTBNMWWxAs2xs3ZkZtRE9qipGVXrcf\nrBIdAB2AFVXmVwIDa3sjs1ZtZeRv/1nbb9ugpKUYzRul07xxOs0bpdO6aQa92mfTvnkjjmqRRfvm\njejQshF5rRrren2pd9JSUzgxvxUn5rcCYh+I1m0rZdH6HRQV72DJhp2s3VrCuu0lTF28k+LtJZRV\n1N/znQdyznHt+c0lA+p0G/XuPgAzGw2MBsjLyzus92jfPIs7hveMvV+VLN5/j7Bqqlddbwdqs/+6\nal5bdZv7v+ZLnxz221DVTw8HfM3+rz1IfzJSU8hMTyEzLYXMtFQy01LISo9NZ6Sl0DQzjcYZqZHe\nRZaGxcxo1zyLds2zGNQt59/Wu8f2jLeVlLGjpJwdpeXsKClnT0UlZRVOeUUleyoqKa9wyioqKat0\n8Ng++L49cPcv7XFXtzfv1E7QdM1pWivvczCJDoBVQKcq8x2DZfu4+zhgHMTOARzORto0y2L04KMP\nt0YRaWDMYodTstJTaaOL0/ZJ9L79p0B3M+tiZhnARcCEBNcgIiIkeA/A3cvN7HrgPWKXgT7l7rMT\nWYOIiMTU6xvBzGw9sOwI3iIH2FBL5SSDqPUX1OeoUJ8PTWd3r/HRevU6AI6UmRXGcy1sQxG1/oL6\nHBXqc93Q9X0iIhGlABARiaiGHgDjwi4gwaLWX1Cfo0J9rgMN+hyAiIhUr6HvAYiISDUUACIiEZX0\nAWBmw8xsvpkVmdltB1ifaWYvBeunmll+4qusXXH0+cdmNsfMZprZRDPrHEadtammPldpd4GZuZkl\n/SWD8fTZzC4MvtezzeyPia6xtsXxs51nZpPMbHrw8z08jDpri5k9ZWbFZjarmvVmZo8G/x8zzax2\nR4dz96T9R+xu4kVAVyADmAH03q/N94HfBdMXAS+FXXcC+nw60DiYvjYKfQ7aZQMfAVOAgrDrTsD3\nuTswHWgZzLcJu+4E9HkccG0w3RtYGnbdR9jnwcAAYFY164cD7xAbH/JkYGptbj/Z9wD2PV/A3fcA\ne58vUNVI4Nlgejww1JJ7GMwa++zuk9x9VzA7hdige8ksnu8zwBjgQaAkkcXVkXj6fA3wW3ffDODu\nxQmusbbF02cHmgXTzYHVCayv1rn7R8CmgzQZCTznMVOAFmbWvra2n+wBcKDnC3Soro27lwNbgdYJ\nqa5uxNPnqq4m9gkimdXY52DXuJO7v5XIwupQPN/nY4BjzOyfZjYleNxqMounz3cD3zGzlcDbwA2J\nKS00h/r7fkjq3fMApPaY2XeAAuCrYddSl8wsBXgY+I+QS0m0NGKHgb5GbC/vIzPr6+5bQq2qbl0M\nPOPuD5nZKcAfzKyPu1eGXVgySvY9gBqfL1C1jZmlEdtt3JiQ6upGPH3GzM4Afgqc6+6lCaqtrtTU\n52ygDzDZzJYSO1Y6IclPBMfzfV4JTHD3MndfQux5290TVF9diKfPVwMvA7j7J0AWsUHTGqq4ft8P\nV7IHQDzPF5gAXBFMfxP40IOzK0mqxj6bWX/g98T++Cf7cWGooc/uvtXdc9w9393ziZ33ONfdC8Mp\nt1bE87P9F2Kf/jGzHGKHhBYnsshaFk+flwNDAcysF7EAWJ/QKhNrAnB5cDXQycBWd19TW2+e1IeA\nvJrnC5jZvUChu08AniS2m1hE7GTLReFVfOTi7PP/AE2BV4Lz3cvd/dzQij5Ccfa5QYmzz+8BZ5rZ\nHKACuMXdk3bvNs4+3wT8r5n9iNgJ4f9I5g90ZvYnYiGeE5zXuAtIB3D33xE7zzEcKAJ2AVfW6vaT\n+P9ORESOQLIfAhIRkcOkABARiSgFgIhIRNXrk8A5OTmen58fdhkiIkll2rRpGzyOZwLHFQBm1gJ4\ngti11g5cBcwHXgLygaXAhe6+ORhmYSyxM9e7iJ2l/yx4nyuAnwVve5+7P8tB5OfnU1iYzFfyiYgk\nnpkti6ddvIeAxgLvuntPoB8wF7gNmOju3YGJwTzA2cRuRukOjAYeDwpqRewSp4HExvy4y8xaxrl9\nERGpZTUGgJk1JzZi3ZMA7r4nuNW86iBrzwKjgunqBi86C3jf3TcFg1e9DyT72CUiIkkrnj2ALsTu\ntHs6GIP7CTNrArStckfaWqBtMF3d4EV1OqiRiIgcmnjOAaQRG6/6BnefamZj+f/DPQC4u5tZrdxR\nZmajiR06Ii8vrzbeUiIm/7aGMiCoHK6lvxgRdglJIZ49gJXASnefGsyPJxYI6/aOSx183TvmTHWD\nF8U1qJG7j3P3AncvyM2t8SS2iIgcphoDwN3XAivMrEewaCgwhy8PsnYF8HowXd3gRXvHLWkZnPw9\nM1gmIiIhiPc+gBuAF4IR+hYTG5AoBXjZzK4GlgEXBm0POHiRu28yszHERvwDuNfdD/YkHBERqUNx\nBYC7f07swSL7G3qAtg5cV837PAU8dSgFiohI3dBQECIiEaUAEBGJKAWAiEhEKQBERCJKASAiElEK\nABGRiFIAiIhElAJARCSiFAAiIhGlABARiSgFgIhIRCkAREQiSgEgIhJRCgARkYhSAIiIRJQCQEQk\nohQAIiIRpQAQEYkoBYCISEQpAEREIkoBICISUQoAEZGIijsAzCzVzKab2ZvBfBczm2pmRWb2kpll\nBMszg/miYH1+lfe4PVg+38zOqu3OiIhI/A5lD+CHwNwq8w8Cv3L3bsBm4Opg+dXA5mD5r4J2mFlv\n4CLgWGAY8JiZpR5Z+SIicrjiCgAz6wiMAJ4I5g0YAowPmjwLjAqmRwbzBOuHBu1HAi+6e6m7LwGK\ngJNqoxMiInLo4t0DeAS4FagM5lsDW9y9PJhfCXQIpjsAKwCC9VuD9vuWH+A1IiKSYGk1NTCzc4Bi\nd59mZl+r64LMbDQwGiAvL++QX59/21u1XZKISIMUzx7AIOBcM1sKvEjs0M9YoIWZ7Q2QjsCqYHoV\n0AkgWN8c2Fh1+QFes4+7j3P3AncvyM3NPeQOiYhIfGoMAHe/3d07uns+sZO4H7r7pcAk4JtBsyuA\n14PpCcE8wfoP3d2D5RcFVwl1AboD/6q1noiIyCGp8RDQQfwEeNHM7gOmA08Gy58E/mBmRcAmYqGB\nu882s5eBOUA5cJ27VxzB9kVE5AgcUgC4+2RgcjC9mANcxePuJcC3qnn9/cD9h1qkiIjUPt0JLCIS\nUQoAEZGIUgCIiESUAkBEJKIUACIiEaUAEBGJKAWAiEhEKQBERCJKASAiElEKABGRiFIAiIhElAJA\nRCSiFAAiIhGlABARiSgFgIhIRCkAREQiSgEgIhJRCgARkYhSAIiIRJQCQEQkohQAIiIRpQAQEYko\nBYCISETVGABm1snMJpnZHDObbWY/DJa3MrP3zWxh8LVlsNzM7FEzKzKzmWY2oMp7XRG0X2hmV9Rd\nt0REpCbx7AGUAze5e2/gZOA6M+sN3AZMdPfuwMRgHuBsoHvwbzTwOMQCA7gLGAicBNy1NzRERCTx\nagwAd1/j7p8F09uBuUAHYCTwbNDsWWBUMD0SeM5jpgAtzKw9cBbwvrtvcvfNwPvAsFrtjYiIxO2Q\nzgGYWT7QH5gKtHX3NcGqtUDbYLoDsKLKy1YGy6pbvv82RptZoZkVrl+//lDKExGRQxB3AJhZU+BV\n4EZ331Z1nbs74LVRkLuPc/cCdy/Izc2tjbcUEZEDiCsAzCyd2B//F9z9tWDxuuDQDsHX4mD5KqBT\nlZd3DJZVt1xEREIQz1VABjwJzHX3h6usmgDsvZLnCuD1KssvD64GOhnYGhwqeg8408xaBid/zwyW\niYhICNLiaDMIuAz4wsw+D5bdAfwCeNnMrgaWARcG694GhgNFwC7gSgB332RmY4BPg3b3uvumWumF\niIgcshoDwN3/AVg1q4ceoL0D11XzXk8BTx1KgSIiUjd0J7CISEQpAEREIkoBICISUQoAEZGIUgCI\niESUAkBEJKIUACIiEaUAEBGJKAWAiEhEKQBERCJKASAiElEKABGRiFIAiIhElAJARCSiFAAiIhGl\nABARiSgFgIhIRCkAREQiSgEgIhJRCgARkYhSAIiIRJQCQEQkohIeAGY2zMzmm1mRmd2W6O2LiEhM\nQgPAzFKB3wJnA72Bi82sdyJrEBGRmETvAZwEFLn7YnffA7wIjExwDSIiQuIDoAOwosr8ymCZiIgk\nWFrYBezPzEYDo4PZHWY2P8x6klQOsCHsIiQ0kf/+24NhVxC6HvE0SnQArAI6VZnvGCzbx93HAeMS\nWVRDY2aF7l4Qdh0SDn3/xcwK42mX6ENAnwLdzayLmWUAFwETElyDiIiQ4D0Ady83s+uB94BU4Cl3\nn53IGkREJCbh5wDc/W3g7URvN2J0CC3a9P2XuH4GzN3ruhAREamHNBSEiEhEKQAaEA2zEW1m9pSZ\nFZvZrLBrkXCYWSczm2Rmc8xstpn98KDtdQioYQiG2VgAfJ3YDXafAhe7+5xQC5OEMbPBwA7gOXfv\nE3Y9knhm1h5o7+6fmVk2MA24KqlvAAAA2ElEQVQYVd3fAe0BNBwaZiPi3P0jYFPYdUh43H2Nu38W\nTG8H5nKQ0RYUAA2HhtkQkX3MLB/oD0ytro0CQESkgTGzpsCrwI3uvq26dgqAhqPGYTZEpOEzs3Ri\nf/xfcPfXDtZWAdBwaJgNkYgzMwOeBOa6+8M1tVcANBDuXg7sHWZjLvCyhtmIFjP7E/AJ0MPMVprZ\n1WHXJAk3CLgMGGJmnwf/hlfXWJeBiohElPYAREQiSgEgIhJRCgARkYhSAIiIRJQCQEQkohQAIiIR\npQAQEYkoBYCISET9Hx0XrFdiuR+oAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 2 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Skv0GLRriN5m",
        "colab_type": "text"
      },
      "source": [
        "## Using GPU and data subsampling with Pyro "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "eHr1j4ZeiSXa",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Setup some data for another model\n",
        "mu = -0.6\n",
        "sigma = 1.8\n",
        "\n",
        "n2 = 10000\n",
        "y2 = stats.norm.rvs(mu, sigma, size=n2)\n",
        "y2_tensor = torch.as_tensor(y2, dtype=torch.float32).cuda()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "eT_xpEVMiffT",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def model(y2_tensor):\n",
        "    # Priors:\n",
        "    prior_dist_mu = dist.Normal(torch.Tensor([0.]).cuda(),\n",
        "                                torch.Tensor([1.]).cuda())\n",
        "    mu = pyro.sample('mu', prior_dist_mu)\n",
        "    \n",
        "    prior_dist_sigma = dist.Gamma(torch.Tensor([1.]).cuda(),\n",
        "                                  torch.Tensor([1.]).cuda())\n",
        "    sigma = pyro.sample('sigma', prior_dist_sigma)\n",
        "    \n",
        "    # Likelihood:\n",
        "    with pyro.plate('observe_data', size=len(y2_tensor),\n",
        "        subsample_size=5000, use_cuda=True) as ind:\n",
        "        pyro.sample('obs', dist.Normal(mu, sigma),\n",
        "            obs=y2_tensor.index_select(0, ind))\n",
        "        \n",
        "\n",
        "def guide(y2_tensor):\n",
        "    alpha_mu = pyro.param(\"alpha_mu\", torch.Tensor([0.0]).cuda())\n",
        "    beta_mu = pyro.param(\"beta_mu\", torch.Tensor([3.0]).cuda(),\n",
        "        constraint=constraints.positive)\n",
        "    mu = pyro.sample('mu', dist.Normal(alpha_mu, beta_mu))\n",
        "\n",
        "    alpha_sigma = pyro.param(\"alpha_sigma\", torch.Tensor([1.0]).cuda(),\n",
        "        constraint=constraints.positive)\n",
        "    beta_sigma = pyro.param(\"beta_sigma\", torch.Tensor([1.0]).cuda(),\n",
        "        constraint=constraints.positive)\n",
        "    sigma = pyro.sample('sigma', dist.Gamma(alpha_sigma, beta_sigma))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "HDEiFzsFkyPO",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# set up the optimizer\n",
        "pyro.clear_param_store()\n",
        "adam_params = {\"lr\": 0.2, \"betas\": (0.90, 0.999)}\n",
        "optimizer = Adam(adam_params)\n",
        "\n",
        "# setup the inference algorithm\n",
        "svi = SVI(model, guide, optimizer, loss=Trace_ELBO())\n",
        "\n",
        "n_steps = 10\n",
        "# do gradient steps\n",
        "for step in range(n_steps):\n",
        "    svi.step(y2_tensor)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "fiq-YmJdmOU4",
        "colab_type": "code",
        "outputId": "18540fd2-4e24-40f7-f081-a3cf9584539b",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 286
        }
      },
      "source": [
        "# Generate replications\n",
        "\n",
        "alpha_mu = pyro.param(\"alpha_mu\")\n",
        "beta_mu = pyro.param(\"beta_mu\")\n",
        "alpha_sigma = pyro.param(\"alpha_sigma\")\n",
        "beta_sigma = pyro.param(\"beta_sigma\")\n",
        "    \n",
        "mu_distribution = stats.norm(alpha_mu.item(), beta.data.item())\n",
        "sigma_distribution = stats.gamma(alpha_sigma.item(), beta_sigma.item())\n",
        "\n",
        "mu_samples = mu_distribution.rvs(no_replications)\n",
        "sigma_samples = sigma_distribution.rvs(no_replications)\n",
        "\n",
        "data_replications = stats.norm(mu_samples, sigma_samples).rvs()\n",
        "\n",
        "# Density estimation using KDE (with tuning parameter chosen by 3 fold CV)\n",
        "params_for_kde_cv = {'bandwidth': np.logspace(-2, 3, 10)}\n",
        "grid = GridSearchCV(KernelDensity(), params_for_kde_cv, cv=3)\n",
        "grid.fit(data_replications.reshape(-1, 1))\n",
        "x_kde = np.linspace(-20, 20, 10000).reshape(-1, 1)\n",
        "y_kde = np.exp(grid.best_estimator_.score_samples(x_kde))\n",
        "plt.plot(x_kde, y_kde)"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "[<matplotlib.lines.Line2D at 0x7f08cfc7fb00>]"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 16
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd4VGXe//H3N51QAoSEkgABQu8Q\nQBGpFkAE7GABKyii7trL6toeV10fdVEXRSyIIrooiqKggmJBgYTeCSGQhBYChJBe7t8fOfjLkw1k\ngCT3lO/runJx5sw9w2eOEz6eLsYYlFJKKT/bAZRSSrkHLQSllFKAFoJSSimHFoJSSilAC0EppZRD\nC0EppRSghaCUUsqhhaCUUgrQQlBKKeUIsB3gdDRq1MjExMTYjqGUUh4lISHhkDEmorJxHlUIMTEx\nxMfH246hlFIeRUR2uzJONxkppZQCtBCUUko5tBCUUkoBWghKKaUcWghKKaUALQSllFIOLQSllFKA\nh52HoJQnKSkxpB3NZcfBLA4cyyczt5D8whICA4Ragf40rhdCs/q1aNWoNmG1Am3HVUoLQamqtDsj\nm+83H+C3xEOs3HWY7IJil17XqlFtukaFcV5sOIPaRdIkLKSakyr137QQlDpLx/OLmL8mjfmrU1m9\n5ygArRvVZmzPKLpEhdE2sg7N6tcirFYgIYH+FBaXkFNQzP7MPFKP5LDj4HHWpRzlj6QMFqzbC0Cn\npvUY3aMZY3tEaTmoGiPGGNsZXBYXF2f00hXKXezLzOX935KZs3IPWXlFtG9cl7E9oxjVrSnNG4ae\n9vsZY9i6P4uftqXz3eb9rNlzFBEYENuIm86LYXC7SPz8pBo+ifJ2IpJgjImrdJwWglKn58CxPKYt\n2cEnq1IoMYYRXZty64BW9GzRoEr/nl2Hspm/Jo3/xKewLzOPNhG1mTSwNVf0iibAX48HUa6r0kIQ\nkeHAvwB/YKYx5vlyzwcDHwC9gQzgGmNMsoj0BWacGAY8aYyZ77wmGcgCioEiV8JqISibjuYUMH3Z\nTmYtT6ao2DCub3MmD2xzRmsDp6OwuISF6/fx9i9JbNp7jNYRtXnw4vZc3LkJIrrGoCpXZYUgIv7A\nduBCIBVYBYw3xmwuM2YK0M0Yc7uIjAMuM8ZcIyKhQIExpkhEmgLrgGbO42QgzhhzyNUPpYWgbCgq\nLmHOyj28tHgbWflFjO0RxV8vaEeL8OotgvKMMfyw5SAvLtrKjoPHiWvZgOcu70q7xnVrNIfyPK4W\ngis7lfsCicaYJOeN5wJjgM1lxowBnnSm5wGvi4gYY3LKjAkBPGf7lFLAquTDPPHlJrbsO8aA2Eb8\nbVRHOjSpZyWLiHBhp8YMaR/BvIRUnl+0lUum/cIdg9pw59BYggP8reRS3sOVDZFRQEqZx6nOvArH\nGGOKgEwgHEBE+onIJmADcLvzPJSWw3cikiAik072l4vIJBGJF5H49PR0Vz6TUmct43g+f/1kLVe9\n+TuZOQVMv64Xs2/pa60Mygrw92Nc3xYsuXcQo7o1Y9rSRC57Yzk704/bjqY8XLXvmTLGrDDGdAb6\nAI+IyIlj6AYYY3oBI4A7RWTgSV4/wxgTZ4yJi4io9IY/Sp0VYwxfrdvLha/8zNfr9zJ1SCxL7hvM\niK5N3W57fXidYF65pgfvTIxjX2Yuo6b9yqfxKZW/UKmTcKUQ0oDmZR5HO/MqHCMiAUAYpTuX/2SM\n2QIcB7o4j9OcPw8C8yndNKWUNQez8rj9wwTu+ngNzRvUYuHd53P/xe2pFeTem2KGdWzMt/cMpHvz\nMB6ct56/fbGBwuIS27GUB3KlEFYBbUWklYgEAeOABeXGLAAmOtNXAkuNMcZ5TQCAiLQEOgDJIlJb\nROo682sDFwEbz/7jKHVmftp2kOGv/sKP29J5ZEQHPrujv0ftrG0SFsJHt57D5IGt+fCPPdz43kqO\n5hTYjqU8TKU7lZ0jgqYCiyk97PRdY8wmEXkaiDfGLADeAWaLSCJwmNLSABgAPCwihUAJMMUYc0hE\nWgPznVXwAGCOMWZRVX84pSpTWFzCS4u38dbPSXRoUpdPJ59DbKTnFEFZ/n7CIyM7EhtZh8fmb+Ty\nfy9n9q39iKpfy3Y05SH0xDTlsw5m5XH77ARW7znK9ee04G+XdCIk0L03D7lqVfJhbn5/FXWCA5h9\nSz9iI+vYjqQscvWwUz3dUfmkjWmZjHn9N7bsy+L1a3vy7NiuXlMGAH1iGvLJpHMpLC7h6rd+Z0Nq\npu1IygNoISifs2jjfq5683cEmHfHuYzq1sx2pGrRqVk9/nN7f2oF+nPdzD/YmKaloE5NC0H5lE9X\npTDlowTaN6nLF1PPo3OzMNuRqlWrRrX5ZPI51A0J5IZ3VrBtf5btSMqNaSEonzHzlyQe/Gw9A9pG\n8PFt5xBZ1zcuKx3dIJSPbu1HoL8f181cQZKewKZOQgtB+YTXl+7g2YVbGNm1CTMnxLn9uQVVLaZR\nbebc1g9jDNfPXMGBY3m2Iyk3pIWgvN7MX5J46bvtXN4zitfG9yIowDe/9rGRdZl1c1+O5hZy8/ur\nyM4vqvxFyqf45m+G8hkfr9zz55rBi1d2w9/HbzDTJSqMN67rxdb9Wdw5ZzVFekazKkMLQXmtr9fv\n5dH5GxjSPoJXr+mpN5VxDGkfyTNjuvDTtnT+vmCT7TjKjeg9lZVXSth9mHs/XUdcywZMv763z24m\nOplr+7Vgz+Ec3ly2ky5RYYzv28J2JOUG9LdEeZ09GTnc9kECzcJCmHFDnFedcFaVHri4PQPbRfD3\nLzexZs8R23GUG9BCUF4lM6eQG99fSYkxvHtjHxrUDrIdyW35+wnTxvWgcVgwd3y4mvSsfNuRlGVa\nCMprFBWXcOec1aQczuGt63vTOkKv31OZ+qFBvHV9HEdzC7j74zUUl3jOtc1U1dNCUF7jpe+282vi\nIf5nbFf6tQ63HcdjdGpWj6fHdOH3pAzeXLbTdhxlkRaC8gqLNu7jzWU7ubZfC67u07zyF6j/46re\n0VzSrSkvf79d9yf4MC0E5fESDx7nvk/X0b15ff5+aSfbcTySiPDcZV1pUi+Ee+au5bietOaTtBCU\nRzueX8Tk2fGEBPoz/bpeBAfoEUVnKqxWIK+O60HqkRye/krPT/BFWgjKYxljeHDeOnYdyua18T1p\npncGO2t9YhoyeVAbPo1P5Zcd6bbjqBqmhaA81tu/JPHNhv08NLwD/WMb2Y7jNe4Z1pbWEbV5+LMN\nuunIx2ghKI/0R1IGLyzaxoguTZg0sLXtOF4lJNCff17Zjb2Zuby4aKvtOKoGaSEoj3PwWB53fbyG\nlg1DefHKboj49gXrqkPvlg25qX8rPvh9Nyt3HbYdR9UQLQTlUYqKS7jr4zVk5RXy7+t7UTck0HYk\nr3X/xe2Iql+Lx7/YSKFeFdUnaCEoj/LSd9tZseswz13WlQ5N6tmO49VCgwJ44tJObDuQxazlybbj\nqBqghaA8xvebD/Dmsp2M79uCy3tF247jEy7q1JjB7SN49YcdHNS7rHk9LQTlEfZk5HDvp2vpElVP\nTz6rQSLCk5d2pqCohOe+2WI7jqpmLhWCiAwXkW0ikigiD1fwfLCIfOI8v0JEYpz5fUVkrfOzTkQu\nc/U9lTohr7CYOz5KQIDp1/XWy1nXsJhGtZk8qDVfrN2rO5i9XKWFICL+wBvACKATMF5Eyv8v2i3A\nEWNMLPAK8IIzfyMQZ4zpAQwH3hKRABffUymMMTy5YBOb9h7j5at70LxhqO1IPmnK4Fia1AvhfxZu\npkSviOq1XFlD6AskGmOSjDEFwFxgTLkxY4BZzvQ8YJiIiDEmxxhz4syWEODEN8mV91SKD1fsYe6q\nFKYMbsMFnRrbjuOzagX5c99F7ViXmsnXG/bZjqOqiSuFEAWklHmc6syrcIxTAJlAOICI9BORTcAG\n4HbneVfeU/m4FUkZPLVgE0M7RHLfRe1tx/F5l/eKpmPTery4aCv5RcW246hqUO07lY0xK4wxnYE+\nwCMiEnI6rxeRSSISLyLx6el6bRVfkXokhykfraZFeCivjuuBv5+efGabv5/w6MgOpB7J5YPlu23H\nUdXAlUJIA8peYD7amVfhGBEJAMKAjLIDjDFbgONAFxff88TrZhhj4owxcRERES7EVZ4ut6CYybMT\nKCgq4e0JcdTTk8/cxvltIxjYLoLXlu7gaE6B7TiqirlSCKuAtiLSSkSCgHHAgnJjFgATnekrgaXG\nGOO8JgBARFoCHYBkF99T+aCSEsMD89axed8xpo3vSRu9DabbeXRkB7Lyi5j+k95dzdtUWgjONv+p\nwGJgC/CpMWaTiDwtIqOdYe8A4SKSCNwLnDiMdACwTkTWAvOBKcaYQyd7z6r8YMozPb9oK1+v38dD\nwzswpEOk7TiqAh2a1GNsjyhm/Z7MwSw9Wc2biDGecwhZXFyciY+Ptx1DVZN3ft3FM19vZsK5LXlq\ndGe9aJ0bSz6UzbCXl3HDOS15cnRn23FUJUQkwRgTV9k4PVNZuYUF6/byzNebGdGlCX+/VMvA3cU0\nqs2VvaKZs2IPe4/m2o6jqogWgrJueeIh7vt0LX1bNeSVa/SIIk9x17BYDIbXlibajqKqiBaCsmrz\n3mNMmp1Aq0a1efuGOL0shQeJbhDK+L4t+E98Crszsm3HUVVAC0FZk3I4hxvfW0ndkABm3dyXsFA9\nvNTT3DkkFn8/4V9LdtiOoqqAFoKy4kh2ARPfW0leYTGzbu5L07BatiOpM9C4Xgg3nNOSL9akkZR+\n3HYcdZa0EFSNyy0o5pZZq0g9ksvMiX1o17iu7UjqLEwa1JpAfz/e+FHPS/B0WgiqRp24BeaalKNM\nG9eDvq0a2o6kzlJk3RDG923BF2vT2JORYzuOOgtaCKrGGGN4/MtN/LDlAE+N7szwLk1tR1JV5PZB\nbfAXYfoyXUvwZFoIqsZMW5LIxyv3MGVwGyacG2M7jqpCTcJCuLpPNPMSUvS8BA+mhaBqxNyVe3jl\nh+1c0SuaBy7WS1l7o9sHtcEYeEvXEjyWFoKqdku3HuCxLzYyqF0Ez1/RVc9C9lLRDUK5olc0H69K\n4eAxvcaRJ9JCUNVqfepR7vxoDR2b1uXf1/Ui0F+/ct5sypA2FJcYZvycZDuKOgP626mqzZ6MHG5+\nfxXhdYJ498Y+1A4OsB1JVbOW4bUZ070ZH63YQ8bxfNtx1GnSQlDV4kh2ATe+t5LCYsP7N/Ulsu5p\n3ShPebApQ2LJKypm5q+7bEdRp0kLQVW5vMJibv0gntSjucycGEdspN7kxpfERtbhkq5N+WB5st5V\nzcNoIagqZYzh0c83kLD7CK9e04M+MXrimS+aOjSW7IJi3l+ebDuKOg1aCKpKzfxlF5+vSePeC9sx\nsqueeOarOjSpx0WdGvPeb8lk5RXajqNcpIWgqsyP2w7yj2+3cEnXptw1NNZ2HGXZ1KGxZOYW8uEf\ne2xHUS7SQlBVIvHgce6es4YOTerxz6u66bkGim7R9RnULoKZvySRW1BsO45ygRaCOmvZ+UVMnh1P\nUIAfb0+MIzRIDy9Vpe4aGktGdgEfr9S1BE+ghaDOijGGR+dvYNehbF67tidR9fW+Bur/i4tpyDmt\nG/LWzzvJL9K1BHenhaDOyscrU/hy7V7+ekE7+rdpZDuOckN3DW3LgWP5zEtItR1FVUILQZ2xjWmZ\nPPnVJga2i+DOIboTWVWsf5tweraoz/SfdlJYXGI7jjoFLQR1RrLyCpk6ZzUNQ4N45eru+PnpTmRV\nMRHhrqGxpB7J5cu1e23HUaeghaDOyFNfbWbP4Rxeu7Yn4XWCbcdRbm5I+0g6Na3Hv39MpLjE2I6j\nTsKlQhCR4SKyTUQSReThCp4PFpFPnOdXiEiMM/9CEUkQkQ3On0PLvOYn5z3XOj+RVfWhVPVatHEf\n8xJSuXNIrJ6JrFxyYi0h6VA2Czfssx1HnUSlhSAi/sAbwAigEzBeRDqVG3YLcMQYEwu8ArzgzD8E\nXGqM6QpMBGaXe911xpgezs/Bs/gcqoYcPJbHI59voFt0GHcPa2s7jvIgF3duQmxkHd5YmkiJriW4\nJVfWEPoCicaYJGNMATAXGFNuzBhgljM9DxgmImKMWWOMObHRcBNQS0R0+4KHMsbwwLz15BYW88o1\nPfTeBuq0+PkJU4fEsu1AFt9vOWA7jqqAK7/RUUBKmcepzrwKxxhjioBMILzcmCuA1caYshdJf8/Z\nXPS4nOTUVhGZJCLxIhKfnp7uQlxVXT78YzfLtqfz2MiOtInQK5iq0zeqW1Nahofy+tJEjNG1BHdT\nI/+LJyKdKd2MNLnM7OucTUnnOz83VPRaY8wMY0ycMSYuIiKi+sOqCqUeyeEf325lYLsIrj+npe04\nykMF+PsxZXAbNqRlsmy7/g+eu3GlENKA5mUeRzvzKhwjIgFAGJDhPI4G5gMTjDF/3n3bGJPm/JkF\nzKF005RyQ6VnI29EgOcu66LXKVJn5bKe0TQLC+E1XUtwO64UwiqgrYi0EpEgYBywoNyYBZTuNAa4\nElhqjDEiUh9YCDxsjPntxGARCRCRRs50IDAK2Hh2H0VVl89Xp/Hz9nQeHN6B6AahtuMoDxcU4Mft\ng9uQsPsIfyQdth1HlVFpITj7BKYCi4EtwKfGmE0i8rSIjHaGvQOEi0gicC9w4tDUqUAs8ES5w0uD\ngcUish5YS+kaxttV+cFU1UjPyueZhZuJa9mAG3RTkaoiV8c1J6JuMK//uMN2FFWGS5elNMZ8A3xT\nbt4TZabzgKsqeN2zwLMnedversdUtjz51SZy8ot5/opuejayqjIhgf5MOr81//PNFhJ2H6F3ywa2\nIyn0TGV1Cku3HmDh+n3cPSxW74usqty1/VrQIDSQN35MtB1FObQQVIXyCov5+4JNxEbWYdLANrbj\nKC9UOziAWwa0YunWg2xMy7QdR6GFoE7izWU7STmcy9OjOxMUoF8TVT0m9I+hbkiAriW4Cf1NV/9l\nT0YO03/ayahuTekfq/c4UNWnXkggN/aP4duN+9l+IMt2HJ+nhaD+y9Nfb8LfT/jbJeUvWaVU1bvp\nvFaEBvnzb11LsE4LQf0fS7Yc4IctB7lnWFuahIXYjqN8QMPaQVx/TksWrNtL8qFs23F8mhaC+lNe\nYTFPfbWZ2Mg63HReK9txlA+59fxWBPj7Mf2nnZUPVtVGC0H96d3fdrHncA5P6Y5kVcMi64Ywvk9z\nPludSsrhHNtxfJb+1isAMo7nM/3HnVzQMZLzdEeysuD2wW3w8xNeX6r7EmzRQlAATFuyg5zCYh4e\n0cF2FOWjmobV4tq+LZi3OpXdGbovwQYtBMXO9ON8tGIP4/s2Jzayru04yodNGdyGAD/hX0v0Gkc2\naCEoXvh2KyGB/vzlgna2oygfF1kvhAnntuSLNWnsTD9uO47P0ULwcSuSMvhu8wHuGNyGRnX07qbK\nvsmD2hAc4M80XUuocVoIPqykxPDcN1toGhbCzXqYqXITjeoEM7F/DAvW7dWzl2uYFoIP+2r9Xtal\nZnL/Re2pFeRvO45Sf5o8sDWhgf786wddS6hJWgg+Kq+wmBcXbaNT03pc1jPKdhyl/o8GtYO4eUAr\nFm7Yx5Z9x2zH8RlaCD5q1vJk0o7m8rdLOuqNb5RbunVAa+qGBPDqD9ttR/EZWgg+6Eh2Aa//mMjQ\nDpF6NVPltsJCA7l1QGsWbzqg90uoIVoIPmja0h1k5xfxiJ6EptzcTQNiCKsVyCvf61pCTdBC8DG7\nDmUz+/fdjOvbgraN9SQ05d7qhQQyaWBrlmw9yNqUo7bjeD0tBB/z4qKtBAf48ZcL2tqOopRLbuwf\nQ8PaQby0eJvtKF5PC8GHxCcf5tuN+7l9UBsi6+q9DpRnqB0cwJTBbfg18RDLEw/ZjuPVtBB8hDGG\nZxduoXG9YG49v7XtOEqdluvPaUmzsBBeWLwNY4ztOF5LC8FHLNywj7UpR7lPT0JTHujEtbbWpRzl\nu80HbMfxWi4VgogMF5FtIpIoIg9X8HywiHziPL9CRGKc+ReKSIKIbHD+HFrmNb2d+YkiMk1E9GD4\napJfVMwLi7bSoUldrugVbTuOUmfk8l5RtImozUuLt1FcomsJ1aHSQhARf+ANYATQCRgvIuXvvn4L\ncMQYEwu8ArzgzD8EXGqM6QpMBGaXec104DagrfMz/Cw+hzqF2b/vJuVwLo9d0hF/PQlNeagAfz/u\nu6g9Ow4eZ/6aNNtxvJIrawh9gURjTJIxpgCYC4wpN2YMMMuZngcMExExxqwxxux15m8CajlrE02B\nesaYP0zpBsEPgLFn/WnUfzmSXcC0JTsY1C6C89tG2I6j1FkZ0aUJXaPCeOX77eQXFduO43VcKYQo\nIKXM41RnXoVjjDFFQCYQXm7MFcBqY0y+Mz61kvdUVeBfS3ZwPL+IR0d2tB1FqbMmIjw4vD1pR3OZ\ns2KP7Thep0Z2KotIZ0o3I00+g9dOEpF4EYlPT0+v+nBeLPFgFrP/2M21/VrQvomehKa8w4DYRpzb\nOpzXlyaSnV9kO45XcaUQ0oDmZR5HO/MqHCMiAUAYkOE8jgbmAxOMMTvLjC+7d7Oi9wTAGDPDGBNn\njImLiNBNHqfj2YVbCA3y5696JzTlRUSEB4a3JyO7gHd/3WU7jldxpRBWAW1FpJWIBAHjgAXlxiyg\ndKcxwJXAUmOMEZH6wELgYWPMbycGG2P2AcdE5Bzn6KIJwJdn+VlUGT9tO8hP29K5Z1hbwvVOaMrL\n9GrRgIs6NWbGz0kcyS6wHcdrVFoIzj6BqcBiYAvwqTFmk4g8LSKjnWHvAOEikgjcC5w4NHUqEAs8\nISJrnZ9I57kpwEwgEdgJfFtVH8rXFRWX8OzCLcSEhzLh3BjbcZSqFvdf3J7jBUVMX7az8sHKJQGu\nDDLGfAN8U27eE2Wm84CrKnjds8CzJ3nPeKDL6YRVrpmzcg+JB4/z9oQ4ggL03EPlndo1rstlPaOY\ntTyZm86LoWlYLduRPJ7+a+FlMnMKefn77fRvE84FHSMrf4FSHuyvF7TDGPTy2FVEC8HL/GvJDo7l\nFvL4qE7oyd/K2zVvGMqEc1syLyGVrfv1VptnSwvBiySlH+eD35O5pk8LOjatZzuOUjVi6tBY6gQH\n8Py3W21H8XhaCF7kuW+2EBLoz30X6WGmynfUDw1i6tBYftqWzm96eeyzooXgJX7Zkc4PWw4ydWgs\njfQwU+VjJpwbQ1T9Wjz3zRZK9MJ3Z0wLwQsUFZfw7NdbaNEwlJvOi7EdR6kaFxLozwMXt2fT3mN8\nuU4vfHemtBC8wNxVKWw7kMWjIzsQHKD3OlC+aXT3ZnSJqsdLi7eTV6gXvjsTWggeLjO39DDTfq0a\ncnHnJrbjKGWNn5/w6IiOpB3NZdbyZNtxPJIWgod748dEjuQU6GGmSgH9YxsxpH0Er/+YqJe0OANa\nCB4s+VA27/22i6t6R9MlKsx2HKXcwsMjOpKdX8TrPybajuJxtBA82HPfbCHI34/7L2pvO4pSbqN9\nk7pc1bs5H/yeTMrhHNtxPIoWgodanniI7zYfYMqQWCLrhdiOo5Rb+euF7fD3E15cvM12FI+iheCB\niksMT3+9maj6tbhlQCvbcZRyO03CQrjt/NZ8tW4vq/ccsR3HY2gheKD/xKewdX8Wj47sSEigHmaq\nVEVuH9SGyLrBPPXVZj1ZzUVaCB4mK6+Ql77bRp+YBozsqoeZKnUytYMDeHB4B9alHOWLtXqymiu0\nEDzMGz/u5NBxPcxUKVdc3jOK7tFhvLBoq95/2QVaCB5kT0YO7/66iyt6RdMtur7tOEq5PT8/4YlL\nO3HgWD5v6Z3VKqWF4EH+8e0W/P2EB4frYaZKuap3y4aM7t6Mt35OIvWIHoZ6KloIHuKPpAy+3bif\nOwa3obEeZqrUaXloRAdE0HsmVEILwQMUlxie+XozzZxD6ZRSpyeqfi0mDWzD1+v3sSr5sO04bksL\nwQN8tjqVTXuP8dCIDtQK0sNMlToTtw9qTZN6ITyth6GelBaCmzueX8Q/F2+jZ4v6jO7ezHYcpTxW\naFAAD4/owIa0TOYlpNqO45a0ENzc9J8SSc/K5wk9zFSpszamRzN6t2zAC4u2cjRHr4ZanhaCG0s5\nnMPbv+xibI9m9GzRwHYcpTyeiPDMmC4cySngn3qdo/+iheDGXli0FT+BB4d3sB1FKa/RqVk9JvaP\nYc7KPaxLOWo7jltxqRBEZLiIbBORRBF5uILng0XkE+f5FSIS48wPF5EfReS4iLxe7jU/Oe+51vmJ\nrIoP5C3ikw/z9fp9TB7Yhmb1a9mOo5RXuffCdkTUCebxLzdSrDuY/1RpIYiIP/AGMALoBIwXkU7l\nht0CHDHGxAKvAC848/OAx4H7T/L21xljejg/B8/kA3ijEudqpk3qhTB5kB5mqlRVqxsSyGOXdGR9\naiZzVu6xHcdtuLKG0BdINMYkGWMKgLnAmHJjxgCznOl5wDAREWNMtjHmV0qLQblo/po01qdm8tCI\n9oQGBdiOo5RXGt29Gee2Duefi7Zy6Hi+7ThuwZVCiAJSyjxOdeZVOMYYUwRkAuEuvPd7zuaix+Uk\nh9CIyCQRiReR+PT0dBfe0rNl5xfx4uKtdI8OY0z38otZKVVVRIRnxnYmt7BYz2B22NypfJ0xpitw\nvvNzQ0WDjDEzjDFxxpi4iIiIGg1ow1vLdnLgWD5PXNoJPz89zFSp6hQbWZdbz2/NvIRUlu88ZDuO\nda4UQhrQvMzjaGdehWNEJAAIAzJO9abGmDTnzyxgDqWbpnxa2tFc3vo5iUu7N6N3y4a24yjlE+4e\n2paW4aE8+vkGcguKbcexypVCWAW0FZFWIhIEjAMWlBuzAJjoTF8JLDXGnHTXvYgEiEgjZzoQGAVs\nPN3w3ubFRaWrrQ/p1UyVqjG1gvz5x+VdSc7I4dUfttuOY1WlheDsE5gKLAa2AJ8aYzaJyNMiMtoZ\n9g4QLiKJwL3An4emikgy8DJwo4ikOkcoBQOLRWQ9sJbSNYy3q+5jeZ6E3Uf4cu1ebju/NdENQm3H\nUcqn9G/TiPF9m/P2L0msT/XdcxPkFP8j73bi4uJMfHy87RhVrqTEMPbfv3HgWB5L7xtM7WA9skip\nmpaZW8hFryyjQWgQX901gEAdBRWgAAAOPUlEQVR/7zlvV0QSjDFxlY3znk/sweYlpLI+NZOHR3TQ\nMlDKkrBagTw7titb92f57N3VtBAsO5ZXyIuLt9KrRX3G9tDDTJWy6cJOjbmkW1OmLUlk+4Es23Fq\nnBaCZa8t2UFGdgFPju6sVzNVyg08NbozdUMC+MvctRQUldiOU6O0ECzamX6c935L5qre0XSLrm87\njlIKaFQnmH9c3pXN+47xryW+ddSRFoIlxhie/moztQL9eeBivZqpUu7kos5NuDoumuk/7STeh265\nqYVgydKtB1m2PZ27h7Ulom6w7ThKqXKeuLQzUQ1qce+n6zieX2Q7To3QQrAgv6iYZ77eTOuI2kzs\nH2M7jlKqAnWCA3j56h6kHMnh2a83245TI7QQLJixLInkjByeGNWJoAD9T6CUu+oT05DJA9swd1UK\n32zYZztOtdN/jWrY7oxsXv8xkZFdmzC4vd4TSCl3d++F7ejevD4PzVvPnowc23GqlRZCDTLG8MSX\nmwj09+OJUZ1tx1FKuSAowI/Xx/dEBKZ+vNqrD0XVQqhB32zYz7Lt6dx3UTuahIXYjqOUclHzhqG8\neGV31qdmevW9E7QQakhWXiFPfbWJLlH1uOGclrbjKKVO0/AuTbixfwzv/raL7zbttx2nWmgh1JCX\nFm8j/Xg+/zO2KwFedNEspXzJIyM70CWqHvf/Zx27M7Jtx6ly+i9TDViRlMGs33cz8dwYujfXM5KV\n8lTBAf78+9reiAiTZyeQU+Bd5ydoIVSznIIiHvxsPS0ahvKg3vhGKY/XIjyU18b3ZPuBLB6Ytx5P\nuoVAZbQQqtk/F29jd0YOL1zRjdAgvbS1Ut5gYLsIHri4AwvX72PGz0m241QZLYRqtHLXYd5fnszE\nc1tybptw23GUUlXo9kGtuaRrU15YtJVfdqTbjlMltBCqSVZeIff/Zx3NG4Ty0Ai9eJ1S3kZEePHK\nbrSNrMvUOWu8YiezFkI1eeLLTaQeyeHlq7vrpiKlvFTt4ABmTOiNCNz0/ioycwptRzorWgjV4PPV\nqcxfk8Y9w9oRF9PQdhylVDVqGV6bt67vTcrhHO74KIHCYs89k1kLoYrtOpTN419spG+rhkwdGms7\njlKqBvRrHc7zl3dj+c4M/jZ/o8ceeaTbMqpQXmExd3+8hgB/P169pgf+fnpLTKV8xRW9o9l1qPTi\nla0jajN5UBvbkU6bFkIVMcbwty82siEtk7cnxNGsfi3bkZRSNezeC9uxKyOb5xdtpWV4bYZ3aWI7\n0mnRTUZV5MM/djMvIZW7h8ZyYafGtuMopSzw8xP+96rudI+uzz1z15Cw27Nuv+lSIYjIcBHZJiKJ\nIvJwBc8Hi8gnzvMrRCTGmR8uIj+KyHEReb3ca3qLyAbnNdNExGO3r8QnH+aprzYzpH0Ef7mgne04\nSimLQgL9mTkxjqZhIdwyK57Eg1m2I7ms0kIQEX/gDWAE0AkYLyKdyg27BThijIkFXgFecObnAY8D\n91fw1tOB24C2zs/wM/kAtu09mssdH60mukEtXh3XEz/db6CUz2tUJ5gPbu5HgJ8fE95Zyb7MXNuR\nXOLKGkJfINEYk2SMKQDmAmPKjRkDzHKm5wHDRESMMdnGmF8pLYY/iUhToJ4x5g9Tujv+A2Ds2XwQ\nG47lFXLTe6vIKyhmxoQ4wmoF2o6klHITLcJDef+mPhzLK2Liuys94hwFVwohCkgp8zjVmVfhGGNM\nEZAJnOpaDVHO+5zqPd1aYXEJd360mp3px5l+fW/aNa5rO5JSys10iQpjxg292XUom9s+iCevsNh2\npFNy+53KIjJJROJFJD493T2uF2KM4bH5G/hlxyGeu7wrA9o2sh1JKeWm+sc24uWre7Bq92Hu/ngN\nRW584porhZAGNC/zONqZV+EYEQkAwoCMSt4zupL3BMAYM8MYE2eMiYuIiHAhbvV7fWkin8aXHlF0\ndVzzyl+glPJpl3ZvxhOjOvHd5gM89NkGSkrc88Q1VwphFdBWRFqJSBAwDlhQbswCYKIzfSWw1Jzi\nVD1jzD7gmIic4xxdNAH48rTTW/DhH7v53++3c1nPKP56oR5RpJRyzU3nteKvF7Tjs9WpPPnVJrc8\nm7nSE9OMMUUiMhVYDPgD7xpjNonI00C8MWYB8A4wW0QSgcOUlgYAIpIM1AOCRGQscJExZjMwBXgf\nqAV86/y4tS/XpvH4lxsZ2iGSF6/shgcfKauUsuDuYbFkFxQx4+ckQoMCeGh4e7f6d8SlM5WNMd8A\n35Sb90SZ6TzgqpO8NuYk8+OBLq4Gte2HzQe499N19GvVkH9f14tAvS+yUuo0iQiPjOhAdn4Rby7b\nSZ1gf6YObWs71p/00hUuWJ54iClzVtOlWT1mTuxDSKC/7UhKKQ8lIjwzpgu5BcW89N12agUFcMuA\nVrZjAVoIlVqz5wi3fhBPTHgo79/UlzrBusiUUmfHz6/05jo5BcU88/Vmagf5M65vC9ux3P+wU5u2\n7j/Gje+tolGdYD68pR8NagfZjqSU8hIB/n5MG9+Twe0jePjzDXyyao/tSFoIJ5N8KJsb3llJSKAf\nH93aj8h6IbYjKaW8TFCAH29e35tB7SJ46LMNzF1ptxS0ECqwLzOX62auoKi4hA9v6UfzhqG2Iyml\nvFRIoD9v3dD7zzWFjy2WghZCORnH87l+5gqO5Rbywc39aKuXpFBKVbOQQH/evL43Q9pH8MjnG5iz\nwk4paCGUkZlbyIR3V5J6JJd3buxD1+gw25GUUj4iJNCfN28oLYVH52/goxW7azyDFoLjeH7pFQm3\nH8jizRt607dVQ9uRlFI+JjigtBSGdojksfkbmf1HzZaCFgKQU1DEze+tYmNaJm9c24sh7SNtR1JK\n+ajgAH+mX9+LCzpG8vgXG5n5S1KN/d0+Xwh5hcXcOiue+N2HeXVcDy7q7Fn3QFVKeZ/SUujNJV2b\n8uzCLUxbsqNGrn3k04WQX1TM7R8m8HtSBv97dXdGdWtmO5JSSgEQ6O/Hv8b14PJeUby5bCdpR6v/\nrms+e9ptYXEJd81Zw0/b0nn+8q5c1jO68hcppVQNCvD346UruzNlcCzRDar/8HefXEMoKi7hL5+s\n5bvNB3hqdGe3OGVcKaUq4ucnxEbWqZm/q0b+FjdSVFzC/f9Zx8L1+3hsZEcm9o+xHUkppdyCT20y\nKiwu4d5P1/HVur08cHF7bhvY2nYkpZRyGz5TCAVFJdwzdw3fbtzPoyM7MGlgG9uRlFLKrfhEIeQX\nFTN1zhq+33yAx0d1cptrjyullDvx+kIoLC7hjg9Xs3TrQZ4e05kJ58bYjqSUUm7J6wshwE9o1ag2\nz13WlWv76dFESil1Ml5fCCLC46M62Y6hlFJuz+cOO1VKKVUxLQSllFKAFoJSSimHFoJSSilAC0Ep\npZRDC0EppRSghaCUUsqhhaCUUgoAqYnbslUVEUkHzvSu042AQ1UYp6portOjuU6P5jo93pqrpTEm\norJBHlUIZ0NE4o0xcbZzlKe5To/mOj2a6/T4ei7dZKSUUgrQQlBKKeXwpUKYYTvASWiu06O5To/m\nOj0+nctn9iEopZQ6NV9aQ1BKKXUKXl8IIvJPEdkqIutFZL6I1C/z3CMikigi20Tk4hrOdZWIbBKR\nEhGJKzM/RkRyRWSt8/OmO+RynrO2vMrleFJE0soso5G2sjh5hjvLJFFEHraZpSwRSRaRDc4yireY\n410ROSgiG8vMaygi34vIDufPBm6Sy/p3S0Sai8iPIrLZ+V28x5lf/cvMGOPVP8BFQIAz/QLwgjPd\nCVgHBAOtgJ2Afw3m6gi0B34C4srMjwE2WlxeJ8tldXmVy/gkcL/t75aTxd9ZFq2BIGcZdbKdy8mW\nDDRygxwDgV5lv9fAi8DDzvTDJ34v3SCX9e8W0BTo5UzXBbY7v3/Vvsy8fg3BGPOdMabIefgHEO1M\njwHmGmPyjTG7gESgbw3m2mKM2VZTf5+rTpHL6vJyY32BRGNMkjGmAJhL6bJSDmPMz8DhcrPHALOc\n6VnA2BoNxUlzWWeM2WeMWe1MZwFbgChqYJl5fSGUczPwrTMdBaSUeS7VmecOWonIGhFZJiLn2w7j\ncLflNdXZDPiujc0NZbjbcinLAN+JSIKITLIdppzGxph9zvR+oLHNMOW4y3cLEYkBegIrqIFl5hX3\nVBaRH4AmFTz1mDHmS2fMY0AR8JE75arAPqCFMSZDRHoDX4hIZ2PMMcu5atSpMgLTgWco/QfvGeB/\nKS179X8NMMakiUgk8L2IbHX+r9itGGOMiLjL4Y5u890SkTrAZ8BfjDHHROTP56prmXlFIRhjLjjV\n8yJyIzAKGGacDXBAGtC8zLBoZ16N5TrJa/KBfGc6QUR2Au2AKtspeCa5qIHlVZarGUXkbeDr6srh\nghpdLqfDGJPm/HlQROZTunnLXQrhgIg0NcbsE5GmwEHbgQCMMQdOTNv8bolIIKVl8JEx5nNndrUv\nM6/fZCQiw4EHgdHGmJwyTy0AxolIsIi0AtoCK21kLEtEIkTE35luTWmuJLupADdaXs4vwwmXARtP\nNrYGrALaikgrEQkCxlG6rKwSkdoiUvfENKUHV9hcTuUtACY60xMBd1kztf7dktJVgXeALcaYl8s8\nVf3LzObe9BraY59I6Tbetc7Pm2Wee4zSI0S2ASNqONdllG5vzgcOAIud+VcAm5ysq4FL3SGX7eVV\nLuNsYAOw3vklaWr5OzaS0iNBdlK62c1aljKZWlN6xNM65/tkLRfwMaWbQgud79YtQDiwBNgB/AA0\ndJNc1r9bwABKN1mtL/Pv1siaWGZ6prJSSinABzYZKaWUco0WglJKKUALQSmllEMLQSmlFKCFoJRS\nyqGFoJRSCtBCUEop5dBCUEopBcD/AzBQFPfQcH74AAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lk1-rqM-HiCu",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# define the PyTorch module that parameterizes the\n",
        "# diagonal gaussian distribution q(z|x)\n",
        "class Encoder(nn.Module):\n",
        "    def __init__(self, z_dim, hidden_dim, input_dim):\n",
        "        super(Encoder, self).__init__()\n",
        "        # setup the three linear transformations used\n",
        "        self.fc1 = nn.Linear(input_dim, hidden_dim)\n",
        "        self.fc21 = nn.Linear(hidden_dim, z_dim)\n",
        "        self.fc22 = nn.Linear(hidden_dim, z_dim)\n",
        "        # setup the non-linearities\n",
        "        self.softplus = nn.Softplus()\n",
        "\n",
        "    def forward(self, x):\n",
        "        # then compute the hidden units\n",
        "        hidden = self.softplus(self.fc1(x))\n",
        "        # then return a mean vector and a (positive) square root covariance\n",
        "        # each of size batch_size x z_dim\n",
        "        z_loc = self.fc21(hidden)\n",
        "        z_scale = torch.exp(self.fc22(hidden))\n",
        "        return z_loc, z_scale\n",
        "\n",
        "\n",
        "# define the PyTorch module that parameterizes the\n",
        "# observation likelihood p(x|z)\n",
        "class Decoder(nn.Module):\n",
        "    def __init__(self, z_dim, hidden_dim, input_dim):\n",
        "        super(Decoder, self).__init__()\n",
        "        # setup the two linear transformations used\n",
        "        self.fc1 = nn.Linear(z_dim, hidden_dim)\n",
        "        self.fc21 = nn.Linear(hidden_dim, input_dim)\n",
        "        self.fc22 = nn.Linear(hidden_dim, input_dim)\n",
        "        # setup the non-linearities\n",
        "        self.softplus = nn.Softplus()\n",
        "\n",
        "    def forward(self, z):\n",
        "        # define the forward computation on the latent z\n",
        "        # first compute the hidden units\n",
        "        hidden = self.softplus(self.fc1(z))\n",
        "        \n",
        "        mu = self.fc21(hidden)\n",
        "        sigma = torch.exp(self.fc22(hidden))\n",
        "        return mu, sigma\n",
        "\n",
        "\n",
        "# define a PyTorch module for the VAE\n",
        "class VAE(nn.Module):\n",
        "    # by default our latent space is 50-dimensional\n",
        "    # and we use 400 hidden units\n",
        "    def __init__(self, input_dim,\n",
        "        z_dim=50, hidden_dim=400, use_cuda=False):\n",
        "        super(VAE, self).__init__()\n",
        "        # create the encoder and decoder networks\n",
        "        self.encoder = Encoder(z_dim, hidden_dim, input_dim=input_dim)\n",
        "        self.decoder = Decoder(z_dim, hidden_dim, input_dim=input_dim)\n",
        "\n",
        "        if use_cuda:\n",
        "            # calling cuda() here will put all the parameters of\n",
        "            # the encoder and decoder networks into gpu memory\n",
        "            self.cuda()\n",
        "        self.use_cuda = use_cuda\n",
        "        self.z_dim = z_dim\n",
        "\n",
        "    # define the model p(x|z)p(z)\n",
        "    def model(self, x):\n",
        "        # register PyTorch module `decoder` with Pyro\n",
        "        pyro.module(\"decoder\", self.decoder)\n",
        "        with pyro.plate(\"data\", x.shape[0]):\n",
        "            # setup hyperparameters for prior p(z)\n",
        "            z_loc = torch.zeros(x.shape[0], self.z_dim, dtype=x.dtype, device=x.device)\n",
        "            z_scale = torch.ones(x.shape[0], self.z_dim, dtype=x.dtype, device=x.device)\n",
        "            # sample from prior (value will be sampled by guide when computing the ELBO)\n",
        "            z = pyro.sample(\"latent\", dist.Normal(z_loc, z_scale).to_event(1))\n",
        "            # decode the latent code z\n",
        "            mu, sigma = self.decoder.forward(z)\n",
        "            # score against actual images\n",
        "            pyro.sample(\"obs\", dist.Normal(mu, sigma).to_event(1), obs=x)\n",
        "            # return the loc so we can visualize it later\n",
        "            #return loc_img\n",
        "\n",
        "    # define the guide (i.e. variational distribution) q(z|x)\n",
        "    def guide(self, x):\n",
        "        # register PyTorch module `encoder` with Pyro\n",
        "        pyro.module(\"encoder\", self.encoder)\n",
        "        with pyro.plate(\"data\", x.shape[0]):\n",
        "            # use the encoder to get the parameters used to define q(z|x)\n",
        "            z_loc, z_scale = self.encoder.forward(x)\n",
        "            # sample the latent code z\n",
        "            pyro.sample(\"latent\", dist.Normal(z_loc, z_scale).to_event(1))\n",
        "\n",
        "    # define a helper function for reconstructing images\n",
        "    def reconstruct_img(self, x):\n",
        "        # encode image x\n",
        "        z_loc, z_scale = self.encoder(x)\n",
        "        # sample in latent space\n",
        "        z = dist.Normal(z_loc, z_scale).sample()\n",
        "        # decode the image (note we don't sample in image space)\n",
        "        loc_img = self.decoder(z)\n",
        "        return loc_img\n",
        "      \n",
        "    def new_instances(self, size=1):\n",
        "         z = stats.norm.rvs(size=(size, self.z_dim))\n",
        "         mu, sigma = self.decoder.forward(torch.as_tensor(z,\n",
        "             device=torch.device('cuda'), dtype=torch.float32))\n",
        "         return stats.norm.rvs(mu.data.cpu().numpy(), sigma.data.cpu().numpy())"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rA1jNJOJ92fu",
        "colab_type": "text"
      },
      "source": [
        "## Variational autoencoders"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6wQWH15L91jW",
        "colab_type": "code",
        "outputId": "647789c1-5973-42a9-daab-77246a6af4e4",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 187
        }
      },
      "source": [
        "# clear param store\n",
        "pyro.clear_param_store()\n",
        "\n",
        "no_instances = 20000\n",
        "input_dim = 3\n",
        "mu = stats.norm.rvs(size=input_dim)\n",
        "\n",
        "# Generate a positive definite matrix\n",
        "sigma = stats.norm.rvs(size=(input_dim, input_dim))\n",
        "sigma[np.triu_indices(input_dim)] = 0\n",
        "sigma += np.diag(np.abs(stats.norm.rvs(size=input_dim)))\n",
        "sigma = np.matmul(sigma.transpose(), sigma) # inverse cholesky decomposition\n",
        "\n",
        "dataset = stats.multivariate_normal.rvs(mu, sigma, size=no_instances)\n",
        "dataset = torch.as_tensor(dataset, dtype=torch.float32)\n",
        "dataset = TensorDataset(dataset)\n",
        "train_loader = DataLoader(dataset, batch_size=1000, shuffle=True,\n",
        "     num_workers=1, pin_memory=True, drop_last=False)\n",
        "\n",
        "# setup the VAE\n",
        "vae = VAE(use_cuda=True, input_dim=input_dim)\n",
        "\n",
        "adam_args = {\"lr\": 0.001}\n",
        "optimizer = Adam(adam_args)\n",
        "\n",
        "# setup the inference algorithm\n",
        "svi = SVI(vae.model, vae.guide, optimizer, loss=Trace_ELBO())\n",
        "\n",
        "train_elbo = []\n",
        "for epoch in range(100):\n",
        "    # initialize loss accumulator\n",
        "    epoch_loss = 0.\n",
        "    # do a training epoch over each mini-batch x returned\n",
        "    # by the data loader\n",
        "    for x, in train_loader:\n",
        "        x = x.cuda()\n",
        "        epoch_loss += svi.step(x)\n",
        "\n",
        "    # report training diagnostics\n",
        "    if not epoch % 10:\n",
        "        normalizer_train = len(train_loader.dataset)\n",
        "        total_epoch_loss_train = epoch_loss / normalizer_train\n",
        "        train_elbo.append(total_epoch_loss_train)\n",
        "        print(\"[epoch %03d]  average training loss: %.4f\" %\n",
        "             (epoch, total_epoch_loss_train))"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "[epoch 000]  average training loss: 12.0625\n",
            "[epoch 010]  average training loss: 5.7273\n",
            "[epoch 020]  average training loss: 5.3826\n",
            "[epoch 030]  average training loss: 5.3734\n",
            "[epoch 040]  average training loss: 5.3898\n",
            "[epoch 050]  average training loss: 5.3576\n",
            "[epoch 060]  average training loss: 5.3630\n",
            "[epoch 070]  average training loss: 5.3428\n",
            "[epoch 080]  average training loss: 5.3604\n",
            "[epoch 090]  average training loss: 5.3348\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "RgPXUrLyQSZv",
        "colab_type": "code",
        "outputId": "d6adf493-c6c3-483e-9197-195805526f41",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 238
        }
      },
      "source": [
        "# Generating new instances (replications) from the trained VAE\n",
        "new_instances = vae.new_instances(100000)\n",
        "\n",
        "print(\"True means\")\n",
        "print(mu)\n",
        "print(\"Empirical means of replications:\")\n",
        "print(new_instances.mean(0))\n",
        "\n",
        "print(\"----------------------------------------\")\n",
        "\n",
        "print(\"True covariance matrix\")\n",
        "print(sigma)\n",
        "print(\"Empirical covariance matrix of replications:\")\n",
        "print(np.cov(new_instances, rowvar=False))"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "True means\n",
            "[ 0.96114028 -0.20429133 -0.83652706]\n",
            "Empirical means of replications:\n",
            "[ 0.96612779 -0.23635043 -0.91886219]\n",
            "----------------------------------------\n",
            "True covariance matrix\n",
            "[[ 4.74371511 -2.46075695 -1.97981744]\n",
            " [-2.46075695  2.67486053  1.06406042]\n",
            " [-1.97981744  1.06406042  1.70799343]]\n",
            "Empirical covariance matrix of replications:\n",
            "[[ 4.73416853 -2.35541449 -1.90316439]\n",
            " [-2.35541449  2.74554436  1.07713679]\n",
            " [-1.90316439  1.07713679  1.70246573]]\n"
          ],
          "name": "stdout"
        }
      ]
    }
  ]
}