model_physics.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import torch\n",
    "import pyro\n",
    "from pyro import distributions as dist\n",
    "from pyro import poutine\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import tqdm\n",
    "#sns.set()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Toy model\n",
    "\n",
    "Input X - corresponds to any kind of phase space input - some **distribution** with prior.\n",
    "Input $\\mu$ - corresponds to parameters of surrogate, ie - size of detector, magnetic field, etc. - a **number**.\n",
    "\n",
    "Define some functions f, g - those can be varied. Probably, bettern to stcik with some types of them, ie:\n",
    "continous, pereodic, polynomial, ...\n",
    "\n",
    "$$X \\sim U(0,1), \\ \\ \\mu \\sim \\delta(\\mu - \\mu_0) $$\n",
    "\n",
    "$$latent\\_x \\sim  \\mathcal{N} (latent\\_x | X, 1)$$\n",
    "\n",
    "$$laten\\_mu \\sim \\mathcal{N}(laten\\_mu |\\mu, 1)$$\n",
    "\n",
    "$$ f(\\alpha) = a + b \\alpha + sin(\\alpha) $$\n",
    "$$ g(\\alpha) = a * \\alpha $$\n",
    "\n",
    "$$y \\sim \\mathcal{N}(y|f(latent_x), abs(g(latent_mu))$$\n",
    "\n",
    "### Problem statement\n",
    "\n",
    "1. Fit GAN/GP/NV to match p(y|X,mu)\n",
    "\n",
    "2. Solve optimisaiton problem of the choise:\n",
    " - $argmin_{mu} (y - Y_0) ^2 2$\n",
    " - $argmax_{mu} p(y|mu) = \\int p(y|X,mu) p(X) dx$\n",
    " - $etc$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make model with `pyro`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu_dist = dist.Delta(torch.tensor(-50., requires_grad=True))\n",
    "p_dist = dist.Exponential(rate=0.01)\n",
    "theta_dist = dist.Uniform(0, 0.01)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample = p_dist.sample([1000]) + 20\n",
    "pt = 1 / sample * torch.sin(theta_dist.sample(sample.shape)) * 10000 # normalising kostyl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFu9JREFUeJzt3X+s5XV95/HnywF/rTq0MLXID8cK2YDG8mOCGDeEpbZLWbiTppjidv3RaKYx2qpx01YNqJDdrMkGaqPRZcUUXaN2Ry13CMZSf/VHInpnABVG06GrYRYiCIJlbXWHvveP80Vv73y/33PuuefOufc7z0dycs/3+/mc7/l8mct7vvP+/EpVIUkalifNuwGSpNkzuEvSABncJWmADO6SNEAGd0kaIIO7JA2QwV2SBsjgLkkDZHCXpAE6Zl5fnCefUDx9e2vZuacd2bZI0maxd+/e71fVtnH15hbcefp2uGCptWhp8cg2RZI2iyTfnaSeaRlJGiCDuyQNkMFdkgZobjn3c08zty5J68Und0kaIIO7JA2QwV2SBijz2mYveU7B77aWVb3zCLdGkjaHJHurase4ej65S9IAGdwlaYAM7pI0QAZ3SRogg7skDZDBXZIGyOAuSQM0cXBPsiXJ7Ulubil7SpJPJjmQ5LYk22fZSEnS6qxm4bA3AfuBZ7WUvRb4QVWdluQK4D3Ab/Vebetz4AInK0nSepjoyT3JycC/Bz7UUWUncGPzfjfwK0my9uZJkqYxaVrmj4E/AP65o/wk4F6AqjoEPAocv+bWSZKmMja4J7kUeKCq9vZVazl32KI1SXYlWUqyxE8eXEUzJUmrMcmT+0uBhSTfAT4BXJTkf66ocxA4BSDJMcBW4OGVF6qq66tqR1Xt4MljN++WJE1pbHCvqrdV1clVtR24AvhCVf3HFdUWgVc37y9v6sxnuUlJ0vTb7CW5GliqqkXgBuCjSQ4wemK/YkbtkyRNYVXBvaq+BHypeX/VsvP/BLx8lg2TJE3PGaqSNEBTp2XW6tzTYGlxXt8uScPmk7skDZDBXZIGaG5pmb0HIAvtZWW6RpLWxCd3SRogg7skDZDBXZIGyOAuSQNkcJekAZrbaBkevQ/2vLuj0B2aJGktfHKXpAEyuEvSABncJWmADO6SNEAGd0kaoEk2yH5qkq8muTPJXUkOG+KS5DVJHkxyR/N63fo0V5I0iUmGQv4YuKiqHktyLPA3ST5bVV9ZUe+TVfXG2TdRkrRaY4N7s9H1Y83hsc1rBptfPxt4y9ovI0k6zEQ59yRbktwBPADcWlW3tVT7zSRfT7I7ySkzbaUkaVUmCu5V9XhVnQWcDJyX5IUrquwBtlfVi4C/BG5su06SXUmWkizBQ2tptySpR0ZZl1V8IHkn8H+r6r91lG8BHq6qrf3XObvgy61lVc9aVZsk6WiRZG9V7RhXb5LRMtuSHNe8fxrwMuBbK+qcuOxwAdi/uuZKkmZpktEyJwI3Nk/kTwL+rKpuTnI1sFRVi8DvJ1kADgEPA68Ze9WtW+ACn9AlaT2sOi0zsy8+bkdxwVJrmXuoSlK7maVlJEmbj8FdkgZojpt1PA57fthRaC5ektbCJ3dJGiCDuyQNkMFdkgZofjl3vgdc11HmBtmStBY+uUvSABncJWmADO6SNEAGd0kaoDl2qLoTkyStF5/cJWmADO6SNEAGd0kaICcxSdIATbLN3lOTfDXJnUnuSvLuljpPSfLJJAeS3JZk+3o0VpI0mUnSMj8GLqqqXwbOAi5Ocv6KOq8FflBVpzF6HH/PbJspSVqNsWmZGu3D91hzeGzzWrk3307gXc373cD7kqSm3MOv5R8Hy9pjykaSxpko595sjr0XOA14f1XdtqLKScC9AFV1KMmjwPHA91dcZxewa3R0Cu3j3Lvy8JKkSU00WqaqHq+qs4CTgfOSvHBFlbR9rOU611fVjtHmrsevvrWSpImsaihkVT0CfAm4eEXRQUaP4iQ5BtgKPDyD9kmSpjDJaJltSY5r3j8NeBnwrRXVFoFXN+8vB74wbb5dkrR2k+TcTwRubPLuTwL+rKpuTnI1sFRVi8ANwEeTHGD0xH7FurVYkjRW5vWAneN2FBcstZbV4hFujCRtEkn2jvot+7n8gCQNkMFdkgbI4C5JAzS3hcOeetqPeN7i7R2lZx/RtkjS0PjkLkkDZHCXpAGa43ru3bLQXeYwSUkab27B/QU8naWW3HpfYJckTca0jCQNkMFdkgbI4C5JAzS3nPveRyB2jkrSuth4o2X2dG+xN+I2e5I0jmkZSRogg7skDdAkOzGdkuSLSfYnuSvJm1rqXJjk0SR3NK+r1qe5kqRJTJJzPwS8tar2JXkmsDfJrVV194p6f11Vl078zffcBzvH5dclSdMY++ReVfdX1b7m/T8A+4GT1rthkqTprSrnnmQ7o/V4b2spfkmSO5N8NskLZtA2SdKUJh4KmeQZwKeAN1fVD1cU7wOeW1WPJbkE+HPg9JZr7AJ2jY62TtlkSdI4Ez25JzmWUWD/WFV9emV5Vf2wqh5r3t8CHJvkhJZ611fVjtHmrk9fY9MlSV3GPrknCXADsL+qru2o84vA96qqkpzH6C+Nh/qv/GzgLattryRpApOkZV4KvBL4RpI7mnNvB04FqKoPApcDr09yCPhH4IqqqnVoryRpAplXDE7OLvjyqj9X9ax1aI0kbQ5J9o5S2/2coSpJAzTHhcO+B1x3+OnLXBhMktbKJ3dJGiCDuyQN0MZbz32MpHs9mipTOpIE8wzuz38OXNsejGuh/SN9gV2S9DOmZSRpgAzukjRABndJGqAN2aGarFx0UpK0GnML7uceB0stHae9gd0JTpI0EdMykjRABndJGqC5pWX27n3c3LokrZMN2aHaac+4vwxcDliSwLSMJA3S2OCe5JQkX0yyP8ldSd7UUidJ/iTJgSRfT3LO+jRXkjSJSdIyh4C3VtW+JM8E9ia5taruXlbn14HTm9eLgQ80PyVJczA2uFfV/cD9zft/SLIfOAlYHtx3Ah9p9k39SpLjkpzYfLbVueduYWnp8Bx5/zh3c+qSNIlV5dyTbAfOBm5bUXQScO+y44PNuZWf35VkKcnSgw8+uLqWSpImNnFwT/IM4FPAm6tq5eN1Wj5y2M7bVXV9Ve2oqh3btm1bXUslSRObKLgnOZZRYP9YVX26pcpB4JRlxycD9629eZKkaUwyWibADcD+qrq2o9oi8Kpm1Mz5wKN9+Xb42SSmlS9J0tpNMlrmpcArgW8kuaM593bgVICq+iBwC3AJcAD4EfA70zfpuu6iPW8Z81k7XCUJJhst8ze059SX1yngDbNqlCRpbZyhKkkDZHCXpAGa38Jhz98C17bkyHf25dV78vEAuJmHJIFP7pI0SAZ3SRqgzbWe+xhd4+SrHCIp6eiy8TbI7lsc7HU9OfWdToCSpCeYlpGkATK4S9IAGdwlaYDmlnPfewDSknOXJK3d5hot09tp2jfByclNko4upmUkaYAM7pI0QJsrLTOl5N2dZVWmbCQNz9jgnuTDwKXAA1X1wpbyC4GbgP/dnPp0VV099psfvQ/2dAfdVpf1BeKOstV+hyQNwCRP7n8KvA/4SE+dv66qS2fSIknSmo3NuVfVXwEPH4G2SJJmZFYdqi9JcmeSzyZ5wYyuKUma0iw6VPcBz62qx5JcAvw5cHpbxSS7gF0Ap556Kt/97uF58r7OzysX39FZds3Cf15VoyVpyNb85F5VP6yqx5r3twDHJjmho+71VbWjqnZs27ZtrV8tSeqw5uCe5BeTpHl/XnPNh9Z6XUnS9CYZCvlx4ELghCQHGY05PBagqj4IXA68Pskh4B+BK6qq1q3FkqSxxgb3qnrFmPL3MRoquSp7H4EsrvZTkqRJbLwZqj0Tla7pWUXyjMXbO0pcelLS0ce1ZSRpgAzukjRABndJGqD55dzvebx9842bntX5kTMWuvLqcPffntN6Pg+NGbhjSl7SAPnkLkkDZHCXpAHaeEMhj7C+sfZlykbSJjW/4L51C1zQkl/v2QR7P8/vLMtl7bn17vHvsH/x7O72SdImZlpGkgbI4C5JA2Rwl6QBOuo7VPvY2Spps9pwwf3Kek9nWd9uS10dp/t7IvSVtbv7uxbd2UnS5mVaRpIGyOAuSQM0yU5MHwYuBR6oqhe2lAd4L3AJ8CPgNVW1b9YN3Wi6sj3m4iVtBBm3I16SC4DHgI90BPdLgN9jFNxfDLy3ql489otP21Fcu3R4Qc8kJi7rXlSs056e6/W6rrvopu4NRQzuktZTkr1VtWNcvbFpmar6K+Dhnio7GQX+qqqvAMclOXHypkqSZm0WOfeTgHuXHR9szkmS5mQWQyHTcq4115NkF7ALgG2nzuCrN57k3Z1lVd3pHEmapVkE94PAKcuOTwbua6tYVdcD1wPs2LGjllry0+nZrKPPlQvvaD3fN15938KZnWWf4Tc6y3az+jH1knQkzSItswi8KiPnA49W1f0zuK4kaUqTDIX8OHAhcEKSg8A7gWMBquqDwC2MRsocYDQU8nfWq7GSpMmMDe5V9Yox5QW8YWYtkiSt2dzWltn7SP/CXENkZ6ukI2XsJKZ1++LOSUzdARDe0l3U1RE75aSoOqltENDImR9Y/QTccZ2tBndJk5jZJCZJ0uZjcJekATK4S9IAzW+zjnseb8+HX9ade+7akANg/+LZref31fmdnzln8e7OsrRPsh3pWBysr3199yVJs+aTuyQNkMFdkgbI4C5JA7ThNsg+WiVt4/F7NgzBsfGSus0xuH+PtuB15eJPOj9xTZ7cWfYuzmk9f85l3R2j9YfdE5XyUE+HasfuTpezu/Mju3vmMO3P87sLJWkKpmUkaYAM7pI0QObcN7Gr6NigJH/Y+Zmq6TZDkbS5zG/hsON2FBe0LBwGVEd+um9VRW7q6Fz8UPdH9i1OtxPTNQsduzu9rvu7zljomeDU4e7Xt/cjAFz1gbd3lhncpeFy4TBJOopNFNyTXJzk20kOJPmjlvLXJHkwyR3Nq+cZVpK03ibZZm8L8H7gVxlthv21JItVtXJhlk9W1Rtn0aje9IvW5MzOzb37hmP2jLfvSocB1bEGj6T1N0mH6nnAgar6e4AknwB2At2rbq3Fniny6tCZW+8dy74wXX9D1wJhH+O3Oz/Tl8PvcvsHzugs283l3R/s2ISkd2EzSYMySVrmJODeZccHm3Mr/WaSryfZneSUtgsl2ZVkKckSP3lwiuZKkiYxSXBve/Rd+ci7B9heVS8C/hK4se1CVXV9Ve2oqh08edvqWipJmtgkwf0gsPxJ/GTgvuUVquqhqvpxc/g/gHNn0zxJ0jQmybl/DTg9yfOA/wNcAfyH5RWSnFhV9zeHC8D+mbZSm1L7YmhwRt3T+ZnOjcTHbHbSNTdCOlqNDe5VdSjJG4HPAVuAD1fVXUmuBpaqahH4/SQLwCHgYeA10zfpLZ0lVy60z8gEuGZn+8Sd7Hm0+6tu6i7qm3TU1XF6zkJ3H/OViz1t75h0dA3dk5G6Ok2hv+P0btp3rArtgXik+8+kU9suW5KOmImWH6iqW4BbVpy7atn7twFvm23TJEnTcoaqJA2QC4cdZboniE2RetlAMsWEqc41jMZca7X5/a6+h59ez/V+tA7mt3BYnlPwuy0lPUHmpp48c0+OvDPPPOVM2CurfUOR3+AznZ85Z3H1c776+hj6dC5s1jdBjO6dnbpmtQLsX2j/b9u3iBqsfvZqVz8r0Ls4XG8bDO7ahFw4TJKOYgZ3SRqgoyLnPk0+9mjUPS59Hb6r68/kCK4n2pku6RlmCh1t7015Hbn+jP5UY1c7NsZG7NOmSd0ovt38gvvW58AFh/+hTDMeHGA/Hasa9vyP2pU7B7hmsSNvDVzTEZiuofszffc1TRv6+hjOWLy9tZ8hmS7I7F/syKtPa5oced+4+Z6+GGjP7/fmwfv+gpkqv98fPMHgpNkzLSNJA2Rwl6QBMrhL0gAdFR2qR6P2sel9uy1pXvo6Ets6C4/0TmXTfN9G6OQc1+6N0Mb1NMdJTGcXfHlVn+ldTbBjMk1vB21Ph+W+hTM7y87JV1rP97Wvz91/e07r+bxntjtFTb2VXt+KjHs6Oib7Ojl3TrHbVl9H5hSTs8ZNLOq/3uH3tl4B98gE976O9nGdwe2mCZz9fybd7ej6rqEGdycxSdJRzOAuSQNkzl2Dt5YUzJG43kb5riNhuklWa7nm6j83yzTPtKmhWaTeJgruSS4G3stos44PVdV/XVH+FOAjjLbXewj4rar6Tu9Ft26BC1rysl05XMZMpunIu3YuogXUSW3bw47ksG1if+aM6shpdy2iBb0TY7LzXav6HoDL2d1Z9p+e9l9az2+9qSeH/6Ep849dufVZb9bRk1fvnYzWM/GtU29/Qdd99QSmvhmvU81s7Q+Cq12IbC2BZLZ56+u6g91Cx/dMuRjeke6UnoexaZkkW4D3A78OnAm8IsnK3sbXAj+oqtMY9Xy8Z9YNlSRNbpKc+3nAgar6+6r6CfAJYOeKOjuBG5v3u4FfSdL9WCxJWleTBPeTgHuXHR9szrXWqapDwKPA8bNooCRp9caOc0/ycuDfVdXrmuNXAudV1e8tq3NXU+dgc3xPU+ehFdfaBexqDv818O1Z3cgcnAB8f96NOIK832HzfjeP51bVtnGVJulQPQicsuz4ZOC+jjoHkxwDbAUeXnmhqroeuH6C79zwkixNMpFgKLzfYfN+h2eStMzXgNOTPC/Jk4ErgJUbjS0Cr27eXw58oeY19VWSNP7JvaoOJXkj8DlGQyE/XFV3JbkaWKqqReAG4KNJDjB6Yr9iPRstSeo30Tj3qroFuGXFuauWvf8n4OWzbdqGN4j00ip4v8Pm/Q7M3BYOkyStH9eWkaQBMrh3SPLhJA8k+eaycz+f5NYkf9f8/LnmfJL8SZIDSb6epH0N3w0sySlJvphkf5K7krypOT/Ie07y1CRfTXJnc7/vbs4/L8ltzf1+shlEQJKnNMcHmvLt82z/NJJsSXJ7kpub48HeK0CS7yT5RpI7kiw15wb5+9zG4N7tT4GLV5z7I+DzVXU68PnmGEZLM5zevHYBHzhCbZylQ8Bbq+oM4HzgDc0yE0O95x8DF1XVLwNnARcnOZ/R0hnXNff7A0ZLa8Awlth4E7B/2fGQ7/UJ/7aqzlo27HGov8+HqypfHS9gO/DNZcffBk5s3p8IfLt5/9+BV7TV26wv4CbgV4+GewaeDuwDXsxoYssxzfmXAJ9r3n8OeEnz/pimXubd9lXc48mMgtlFwM1Ahnqvy+75O8AJK84N/vf5iZdP7qvz7Kq6H6D5+QvN+UmWaNg0mn+Gnw3cxoDvuUlT3AE8ANwK3AM8UqMlNOBf3tNmX2Ljj4E/AP65OT6e4d7rEwr4iyR7m9nxMODf55Vcz3022hZJ25TDkJI8A/gU8Oaq+mHP+m+b/p6r6nHgrCTHAZ8Bzmir1vzctPeb5FLggaram+TCJ063VN3097rCS6vqviS/ANya5Fs9dYdyzz/lk/vqfC/JiQDNzwea85Ms0bDhJTmWUWD/WFV9ujk96HsGqKpHgC8x6ms4rllCA/7lPf30fvuW2NigXgosJPkOo1VdL2L0JD/Ee/2pqrqv+fkAo7+8z+Mo+H1+gsF9dZYvs/BqRnnpJ86/qulxPx949Il/+m0WGT2i3wDsr6prlxUN8p6TbGue2EnyNOBljDobv8hoCQ04/H435RIbVfW2qjq5qrYzmj3+har6bQZ4r09I8q+SPPOJ98CvAd9koL/Pread9N+oL+DjwP3A/2P0t/prGeUdPw/8XfPz55u6YbShyT3AN4Ad827/FPf7bxj9M/TrwB3N65Kh3jPwIuD25n6/CVzVnP8l4KvAAeB/AU9pzj+1OT7QlP/SvO9hyvu+ELh56Pfa3Nudzesu4B3N+UH+Pre9nKEqSQNkWkaSBsjgLkkDZHCXpAEyuEvSABncJWmADO6SNEAGd0kaIIO7JA3Q/wdEwyuGyAJWbgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "my_cmap = plt.cm.jet\n",
    "my_cmap.set_under('white')\n",
    "plt.hist2d(sample.numpy(), pt.numpy(), bins=50, cmap=my_cmap, cmin=0.0000001);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [],
   "source": [
    "def y_sample(offset=20):\n",
    "    mu = pyro.sample('mu', mu_dist)\n",
    "    x = pyro.sample(\"x\", dist.Normal(mu, 1))\n",
    "    \n",
    "    p = pyro.sample('p', p_dist) + offset\n",
    "    theta = torch.sin(pyro.sample('theta', theta_dist))\n",
    "    #pt = 1 / p * torch.sin(theta) * 10000\n",
    "    \n",
    "    theta_s = 1 / p\n",
    "    y_var = theta_s * abs(x) ** 3 / 6\n",
    "    y_mean = x ** 2 / p\n",
    "    return pyro.sample(\"y\", dist.Normal(y_mean, torch.sqrt(y_var)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {},
   "outputs": [],
   "source": [
    "conditioned_y_sample = poutine.condition(y_sample, data={\"mu\": torch.tensor(24.)})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAEJ9JREFUeJzt3W+M5VV9x/H3p6Bg1br8Gehmd+1g3FhsU9BMEEMfKFgFNEIbaWiMbs02+wQNRhO71KSNSR/gExGT1oSIdW2sSlHKBolKF4jpA9FZQQRXykqpTJayawVsYzRFv31wz7aXZZa5M3tn7p2z71dyc3+/c8/M/e7c4XMP53fumVQVkqR+/dqkC5AkrS6DXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktS5EyddAMDpp59es7Ozky5DktaVvXv3/riqZpbqNxVBPzs7y/z8/KTLkKR1Jcm/j9JvpKmbJI8m+V6S+5LMt7ZTk9yR5OF2f0prT5JPJNmf5P4kr135P0OSdKyWM0f/xqo6t6rm2vlOYE9VbQX2tHOAS4Ct7bYD+OS4ipUkLd+xXIy9DNjVjncBlw+1f7YGvglsSLLxGJ5HknQMRg36Ar6eZG+SHa3tzKp6HKDdn9HaNwGPDX3tQmuTJE3AqBdjL6iqA0nOAO5I8oPn6ZtF2p6z6X17w9gB8PKXv3zEMiRJyzXSiL6qDrT7g8AtwHnAE4enZNr9wdZ9Adgy9OWbgQOLfM8bqmququZmZpZcHSRJWqElgz7Ji5O89PAx8GbgAWA3sK112wbc2o53A+9uq2/OB54+PMUjSVp7o0zdnAnckuRw/3+oqq8m+TZwU5LtwI+AK1r/24FLgf3Az4D3jL1qSdLIlgz6qnoEOGeR9v8ELlqkvYCrxlKdJOmYTcUnYzUeszu/smj7o9e+dY0rkTRN3NRMkjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUufcAuE4MLw1gtshSMcfR/SS1DmDXpI6Z9BLUueco1/njrY18Sj9na+Xjg+O6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnXF65Di13SaWk45sjeknqnEEvSZ0z6CWpc87RrxPOy0taKUf0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdc9XNccw/QiIdHxzRS1LnRg76JCckuTfJbe38rCT3JHk4yReTvLC1n9TO97fHZ1endEnSKJYzor8a2Dd0/lHguqraCjwJbG/t24Enq+qVwHWtnyRpQkYK+iSbgbcCn2rnAS4Ebm5ddgGXt+PL2jnt8Ytaf0nSBIw6ov848CHgV+38NOCpqnqmnS8Am9rxJuAxgPb4063/syTZkWQ+yfyhQ4dWWL4kaSlLBn2StwEHq2rvcPMiXWuEx/6/oeqGqpqrqrmZmZmRipUkLd8oyysvAN6e5FLgZOA3GIzwNyQ5sY3aNwMHWv8FYAuwkORE4GXAT8ZeuSRpJEuO6KvqmqraXFWzwJXAnVX1TuAu4B2t2zbg1na8u53THr+zqp4zopckrY1jWUf/58AHkuxnMAd/Y2u/ETittX8A2HlsJUqSjsWyPhlbVXcDd7fjR4DzFunzc+CKMdQmSRoDPxkrSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjq3rG2K1a/ZnV951vmj1751QpVIGjeDfoodGb6StBIGvZY0/IbjSF9af5yjl6TOGfSS1DmDXpI65xy9FuWFYKkfjuglqXMGvSR1zqCXpM4Z9JLUOS/Galn88JS0/jiil6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpc0sGfZKTk3wryXeTPJjkI639rCT3JHk4yReTvLC1n9TO97fHZ1f3nyBJej6jjOh/AVxYVecA5wIXJzkf+ChwXVVtBZ4Etrf+24Enq+qVwHWtnyRpQpYM+hr473b6gnYr4ELg5ta+C7i8HV/WzmmPX5QkY6tYkrQsI83RJzkhyX3AQeAO4IfAU1X1TOuyAGxqx5uAxwDa408Dp42zaEnS6EYK+qr6ZVWdC2wGzgPOXqxbu19s9F5HNiTZkWQ+yfyhQ4dGrVeStEzLWnVTVU8BdwPnAxuSHN4UbTNwoB0vAFsA2uMvA36yyPe6oarmqmpuZmZmZdVLkpY0yqqbmSQb2vGLgDcB+4C7gHe0btuAW9vx7nZOe/zOqnrOiF6StDZG2aZ4I7AryQkM3hhuqqrbknwf+EKSvwbuBW5s/W8E/j7JfgYj+StXoW5J0oiWDPqquh94zSLtjzCYrz+y/efAFWOpTlPNveml9cFPxkpS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1LlRPjClNTS8Nl2SxsGg11j44Slpejl1I0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzboEwBdzfRtJqckQvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5V91oVR1tRZF/nERaO47oJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUudcXqmxc5M2abosOaJPsiXJXUn2JXkwydWt/dQkdyR5uN2f0tqT5BNJ9ie5P8lrV/sfIUk6ulGmbp4BPlhVZwPnA1cleTWwE9hTVVuBPe0c4BJga7vtAD459qolSSNbMuir6vGq+k47/i9gH7AJuAzY1brtAi5vx5cBn62BbwIbkmwce+WSpJEs62JsklngNcA9wJlV9TgM3gyAM1q3TcBjQ1+20NokSRMwctAneQnwJeD9VfXT5+u6SFst8v12JJlPMn/o0KFRy5AkLdNIQZ/kBQxC/nNV9eXW/MThKZl2f7C1LwBbhr58M3DgyO9ZVTdU1VxVzc3MzKy0fknSEkZZdRPgRmBfVX1s6KHdwLZ2vA24daj93W31zfnA04eneCRJa2+UdfQXAO8Cvpfkvtb2F8C1wE1JtgM/Aq5oj90OXArsB34GvGesFXfCteaS1sqSQV9V/8Li8+4AFy3Sv4CrjrEuSdKYuAWCJHXOoJekzhn0ktQ5NzXTRAxfjPbvx0qryxG9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pzLKzVVXHYpjZ9Br4lzgzdpdTl1I0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzboGwhvyov6RJcEQvSZ0z6CWpcwa9JHXOOXpNLfeml8bDEb0kdc4RvdY1R/3S0gx6rQsuTZVWzqkbSeqcQS9JnTPoJalzBr0kdW7Ji7FJPg28DThYVb/b2k4FvgjMAo8Cf1xVTyYJcD1wKfAz4E+r6jurU7r0bK7AkRY3yoj+M8DFR7TtBPZU1VZgTzsHuATY2m47gE+Op0xJ0kotGfRV9Q3gJ0c0Xwbsase7gMuH2j9bA98ENiTZOK5iJUnLt9I5+jOr6nGAdn9Ga98EPDbUb6G1PUeSHUnmk8wfOnRohWVIkpYy7ouxWaStFutYVTdU1VxVzc3MzIy5DEnSYSsN+icOT8m0+4OtfQHYMtRvM3Bg5eVJko7VSrdA2A1sA65t97cOtb83yReA1wFPH57ikSbF1Tg63o2yvPLzwBuA05MsAH/FIOBvSrId+BFwRet+O4OllfsZLK98zyrULElahlQtOoW+pubm5mp+fn7SZaw6N+aaLo7utd4l2VtVc0v1c/fKVWa4S5o0g17HLefudbxwrxtJ6pxBL0mdM+glqXPO0UtHcO5evXFEL0mdc0Qv4TJY9c0RvSR1zhH9KnB02I8jX0vn7LUeOaKXpM4Z9JLUOadupBVyGabWC0f0ktQ5R/TSGDi61zRzRC9JnXNELy3DKEtnHd1r2jiil6TOGfSS1DmDXpI65xz9mLjtgaRpZdBLq2i5F2a9kKvVYNBLa+RoIX60/xs09DUuztFLUucMeknqnFM3x8ALsFopf3e0lgx6aR0YZX7feXwdjUG/TI7EtB74BqBhBr20zjjY0HIZ9FInfAPQ0Rj0UuecxpFBPwJHSpLWM9fRS1LnVmVEn+Ri4HrgBOBTVXXtajyPpOVZ7v+dOtXTh7EHfZITgL8B/gBYAL6dZHdVfX/cz7WanK6RRp/fX+3rAF5nODarMaI/D9hfVY8AJPkCcBmwroJe0rONOvjxw13TZzWCfhPw2ND5AvC6VXgeYPk7AkpaO8vdmfNY/ibvcr/ncvusxLS86aWqxvsNkyuAt1TVn7XzdwHnVdX7jui3A9jRTl8FPDTWQp7f6cCP1/D5Vmo91LkeaoT1Uac1js96qHMcNf5WVc0s1Wk1RvQLwJah883AgSM7VdUNwA2r8PxLSjJfVXOTeO7lWA91rocaYX3UaY3jsx7qXMsaV2N55beBrUnOSvJC4Epg9yo8jyRpBGMf0VfVM0neC3yNwfLKT1fVg+N+HknSaFZlHX1V3Q7cvhrfe0wmMmW0AuuhzvVQI6yPOq1xfNZDnWtW49gvxkqSpotbIEhS546roE9yRZIHk/wqydwRj12TZH+Sh5K8ZVI1tloubnXsT7JzkrUMS/LpJAeTPDDUdmqSO5I83O5PmXCNW5LclWRfe62vnrY6k5yc5FtJvttq/EhrPyvJPa3GL7bFDBOX5IQk9ya5rZ1PVZ1JHk3yvST3JZlvbVPzeg/VuSHJzUl+0H4/X79WdR5XQQ88APwR8I3hxiSvZrA66HeAi4G/bVs5rLmhLSQuAV4N/Emrbxp8hsHPZ9hOYE9VbQX2tPNJegb4YFWdDZwPXNV+ftNU5y+AC6vqHOBc4OIk5wMfBa5rNT4JbJ9gjcOuBvYNnU9jnW+sqnOHlitO0+t92PXAV6vqt4FzGPxM16bOqjrubsDdwNzQ+TXANUPnXwNeP6HaXg987Wi1TfoGzAIPDJ0/BGxsxxuBhyZd4xH13spg36WprBP4deA7DD49/mPgxMV+DyZY32YGAXQhcBuQaasTeBQ4/Yi2qXq9gd8A/o12XXSt6zzeRvRHs9i2DZusZSRnVtXjAO3+jAnX83+SzAKvAe5hyups0yH3AQeBO4AfAk9V1TOty7S87h8HPgT8qp2fxvTVWcDXk+xtn7iHKXu9gVcAh4C/a9Ngn0ryYtaozu7+8EiSfwZ+c5GHPlxVtx7tyxZpm9RypGmqZd1K8hLgS8D7q+qnyWI/1smpql8C5ybZANwCnL1Yt7Wt6tmSvA04WFV7k7zhcPMiXSf9+3lBVR1IcgZwR5IfTLiexZwIvBZ4X1Xdk+R61nA6qbugr6o3reDLRtq2YY1MUy2jeCLJxqp6PMlGBiPUiUryAgYh/7mq+nJrnro6AarqqSR3M7iesCHJiW20PA2v+wXA25NcCpzMYPrh40xZnVV1oN0fTHILgx10p+31XgAWquqedn4zg6BfkzqduhnYDVyZ5KQkZwFbgW9NqJb1toXEbmBbO97GYE58YjIYut8I7Kuqjw09NDV1JplpI3mSvAh4E4MLc3cB72jdJv6zrKprqmpzVc0y+D28s6reyRTVmeTFSV56+Bh4M4NFF1PzegNU1X8AjyV5VWu6iMHW7WtT5yQvUEzggsgfMnhn/QXwBM++6PlhBvOkDwGXTLjOS4F/bfV8eNI/t6G6Pg88DvxP+zluZzBnuwd4uN2fOuEaf5/BVML9wH3tduk01Qn8HnBvq/EB4C9b+ysYDDD2A/8InDTp13yo5jcAt01bna2W77bbg4f/e5mm13uo1nOB+fa6/xNwylrV6SdjJalzTt1IUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOve/4BcyoWZOtOQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist([conditioned_y_sample().item() for _ in range(10000)], bins=100);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return x.sin() + x.tanh()\n",
    "\n",
    "#@pyro.condition(data={\"x\": torch.tensor(1.)})\n",
    "def y_sample():\n",
    "    mu = pyro.sample('mu', mu_dist)\n",
    "    x = pyro.sample(\"x\", dist.Normal(mu, 1))\n",
    "    x = f(x)\n",
    "    return pyro.sample(\"y\", dist.Normal(x, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Statistic of interest\n",
    "\n",
    "$$r = \\frac{1}{N} \\sum\\limits_{i=0}^{N} || y_i - Y ||_2^2$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def R(ys: torch.tensor, Y_0=-5):\n",
    "    return (ys - Y_0).pow(2).mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Approaches to solve the problem\n",
    "\n",
    "### First one\n",
    "\n",
    "  1. Make approximation of this model:\n",
    "\n",
    "  $$p(y | \\mu)$$\n",
    "\n",
    "  With GANs or some other techniques.\n",
    "  \n",
    "\n",
    "  2. Sample $\\{ y_i \\}_{i=0}^{B} \\sim p(y | \\mu)$\n",
    "  \n",
    "  \n",
    "  3. Calc $r = R(\\{ y_i \\}_{i=0}^{B})$\n",
    "  \n",
    "  \n",
    "  4. Find $\\frac{\\partial r}{\\partial \\mu}$\n",
    "  \n",
    "### Second one\n",
    "\n",
    "\n",
    "  1. Make approximation of this model:\n",
    "\n",
    "  $$p(r | \\mu)$$\n",
    "  \n",
    "  \n",
    "  2. Sample $r \\sim p(r | \\mu)$\n",
    "  \n",
    "  \n",
    "  3. Find $\\frac{\\partial r}{\\partial \\mu}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Investigation of the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "samples = torch.stack([y_sample() for _ in range(1000)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEBCAYAAABysL6vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAEc1JREFUeJzt3X+MZXdZx/H3zO7O7EJnsW4G2/LDmkCfICI/C4a2aELTBMU0Db8iJIhSVmLVmhQIYkkpaiBRGpIiVsFitQoYmkLoIgG7kABFLCo0UnxCtMRAN2YzEuk2dGc6M/5x57Znp3Nnzjn33rn3fOf9Sgh37pxzz/O9c86np9/73G9n1tfXkSSVZ3bSBUiSxsOAl6RCGfCSVCgDXpIKZcBLUqEMeEkqlAEvSYUy4CWpUAa8JBXKgJekQhnwklSo/RM45jxwIXACWJ3A8SWpi/YB5wJ3A6fr7DCJgL8Q+NIEjitJJbgE+HKdDScR8CcAfvCDB1lba76S5ZEjZ7G0dGrkRe2mro+h6/VD98dg/ZO322OYnZ3h7LMfDxsZWsckAn4VYG1tvVXA9/ftuq6Poev1Q/fHYP2TN6Ex1J7a9kNWSSqUAS9JhTLgJalQBrwkFcqAl6RCGfCSVCgDXpIKNYk+eGlqLRw+xMH53mXx0OmHeeCHPxrLPtJuMOClioPz+/nlaz4FwKffdzkPjGkfaTc4RSNJhTLgJalQBrwkFarWHHxEXAe8euPHY5n5toj4CHAx8ODG89dn5u1jqFGS1MKOAR8RlwKXAc8F1oHPRsQVwAuAl2Rm7aUrJUm7p84d/AngmsxcBoiIbwNP3fjfzRHxJOB2enfwa2OrVJLUyI4Bn5nf6j+OiKfTm6q5BPgF4DeB/wPuAN4IfKjugY8cOathqY9aXFxove+06PoYul4/1BtDm3Hu1nvT9b9B1+uH6R9D7T74iHgmcAx4a2YmcEXldzcCr6dBwC8tnWq1WP7i4gInT3a707jrY+h6/TB4DJsv2DrjbLPPsLr+N+h6/bD7Y5idnWl8Y1yriyYiLgLuBN6embdExLMi4hWVTWaAlUZHliSNVZ0PWZ8CfBJ4TWYe33h6Bnh/RBwHTgFHgVvGVqUkqbE6UzRvAQ4CN0RE/7mbgPcAXwEOALdl5kfHUqEkqZU6H7JeDVw94NcfHG05kqRR8ZusklQoA16SCuVywSrCONZkX15ZfaQF8vTyKvNz+0b6+tK4GfAqwjjWZJ87sO+M13TNd3WNUzSSVCgDXpIKZcBLUqEMeEkqlAEvSYUy4CWpULZJauqNo8dd2gsMeE29cfS4S3uBUzSSVCgDXpIKZcBLUqEMeEkqlAEvSYUy4CWpULZJqrZqP/rC4UON+tHH0ctefc0621TXdF9eWR36+NK0M+BV2zD96OPoZd/8mnW22Wl7qSRO0UhSoQx4SSqUAS9JhTLgJalQBrwkFcqAl6RC2SYpdZTr5GsnBrzUUa6Tr504RSNJhTLgJalQtaZoIuI64NUbPx7LzLdFxKXADcAh4OOZee2YapQktbDjHfxGkF8GPBd4DvD8iPgV4GbgcuAZwIUR8bJxFipJaqbOFM0J4JrMXM7MFeDbwAXAdzLzvsx8GLgVeNUY65QkNbTjFE1mfqv/OCKeTm+q5kZ6wd93AnjyyKuTJLVWu00yIp4JHAPeCjxM7y6+bwZYa3LgI0fOarL5GRYXF1rvOy32+hja7ru8slpr31Ft02bfnV53eWWVuQP7HvN42O27dk51rd6tTPsY6n7IehFwG/C7mfmxiPh54NzKJucA9zc58NLSKdbW1pvsAvTe0JMnu93x29UxbD6Zm4xhVPvOHdhXa033/utvdwFuVUPdC7a6b9OxLS4unDGGNttvVWeXzqmuXgNVuz2G2dmZxjfGOwZ8RDwF+CTwmsw8vvH013q/iqcB9wGvpfehqyRpStS5g38LcBC4ISL6z90EvIHeXf1B4DPAJ8ZQnySppTofsl4NXD3g188ebTmSpFHxm6ySVCgDXpIKZcBLUqFcLlgaoWqf/jjWaK/7PQAJDHhppDb36Y+6S7ru9wAkcIpGkoplwEtSoQx4SSqUAS9JhTLgJalQBrwkFco2Sakhe9HVFQa81FC1Fx3sR9f0copGkgplwEtSoQx4SSqUAS9JhTLgJalQBrwkFco2SbFw+BAH53unwunlVebn9gHt1jOvvtY41kMflWove3XMJejK30DjZ8CLg/P7z1hjfJj1zDe/1qjXQx+Vzeuql7TGelf+Bho/p2gkqVAGvCQVyoCXpEIZ8JJUKANekgplwEtSoWyT1ETtlZ7t6jirqv34uzn+vfK+73UGvCZqr/Rsbx5n3+Z+/N0a/1553/c6p2gkqVAGvCQVqvYUTUQcBu4CXp6Z342IjwAXAw9ubHJ9Zt4+hholSS3UCviIeBHwIeCCytMvAF6SmSfGUZgkaTh1p2jeBFwF3A8QEY8DngrcHBH3RMT1EeF0jyRNkVqhnJlXZuaXKk+dAxwHfh34OeAS4I2jL0+S1FarNsnM/C/giv7PEXEj8Hp60zi1HDlyVptDAzzSN9xlXRnDdnX2f7e8ssrcga3XUx+0/1bPn7FG+8oq8wNes45peX+b1lHnPd1K9b1rU8M43q/qGLYaz7T8jYYx7WNoFfAR8Szggsy8beOpGWClyWssLZ1ibW298bEXFxc4ebLbXbvTNobtTtJqnZu36/9ucXFh4Hrq1W12en6Ua7QPOu5ua1pHnfd0K5vfu61es2rQ32OUNo9h87k0TddAG7s9htnZmcY3xm2/6DQDvD8ijgOngKPALS1fS5I0Bq0+GM3Me4D3AF8B7gW+kZkfHWVhkqThNLqDz8zzK48/CHxw1AVJkkbD1kZJKpQBL0mFMuAlqVAGvCQVyoCXpEIZ8JJUKANekgplwEtSoQx4SSqUAS9JhTLgJalQbVeT1B6w3RrjKtPC4UMcnO/FwkOnH+aBH/5owhVpGAa8BqquMQ7t1mVXtxyc33/GGu7dXrFdTtFIUqEMeEkqlAEvSYUy4CWpUAa8JBXKgJekQtkm2THT0qdcp0d+0Dbj7q+3f3/06px31W00HfxrdMy09ClXe+QH9ccP2qbOvuOuTc3UOe+q2/S302Q5RSNJhTLgJalQBrwkFcqAl6RCGfCSVCgDXpIKZZukNCZN+/G70r9frdM146ebAS+NSdN+/K7072+u0zXjp5dTNJJUKANekgpVa4omIg4DdwEvz8zvRsSlwA3AIeDjmXntGGuUJLWw4x18RLwI+DJwwcbPh4CbgcuBZwAXRsTLxlmkJKm5OlM0bwKuAu7f+PmFwHcy877MfBi4FXjVmOqTJLW04xRNZl4JEBH9p84DTlQ2OQE8eeSVSZKG0qZNchZYr/w8A6w1fZEjR85qceieLvQK72RUY+i/zvLKKnMH9j3m8SB1tlEZ6pxrw6zbX/dc3ryd1/H4tQn47wHnVn4+h0enb2pbWjrF2tr6zhtusri4wMmT3e68HWYMm0+o/ussLi6c0Zu80+tv3l7l2upc2O486qvbl7/VvjvVsdev4zZmZ2ca3xi3CfivARERTwPuA15L70NXSdIUadwHn5kPAW8AbgPuBf4D+MRoy5IkDav2HXxmnl95fCfw7HEUJEkaDb/JKkmFMuAlqVAGvCQVyuWCO6wr64drcuqs3T7MeeQ5ON0M+A7ryvrhmpw6a7cPcx55Dk43p2gkqVAGvCQVyoCXpEIZ8JJUKANekgplwEtSoWyTLFC1N/n08irzc/se81h7zyR71jf347e1cPgQB+f3P/I6W/X161EGfIE29yYPeqy9ZZI966M69sH5/Tv29etRTtFIUqEMeEkqlAEvSYUy4CWpUAa8JBXKgJekQtkmKWmqVXvf1Yx38JKmWr/3vd//rvoMeEkqlAEvSYUy4CWpUAa8JBXKgJekQhnwklQom0s7wD5gSW14B98B9gFLasOAl6RCGfCSVKihJnYj4gvAE4GVjad+IzO/NnRVkqShtQ74iJgBLgB+MjPb/1d0JUljMcwUTWz8/+ci4psR8VujKEiSNBrDTNGcDdwJ/DZwAPhiRGRmfr7OzkeOnNX6wIuLC633HbfllVXmDux7zOPNpnkM0jj1r4vFxYWB18t2107VpK+jSR9/J60DPjO/Cny1/3NE/CXwi0CtgF9aOsXa2nrj4y4uLnDy5AON99sti4sLj7Qzfvp9l29Za9MxTPtJJDUxd2DfltfIoGtnu/N/klmw21k0OzvT+Ma49RRNRFwcES+tPDXDox+2SpImbJgpmh8D3h0RL6Y3RfOrwJtHUpUkaWit7+Az8w7gGPBvwL8AN29M20iSpsBQffCZ+U7gnSOqRZI0Qn6TVZIKZcBLUqFcg3aMlldWH2nxeuj0wzzwwx9NuCKpfNXltU8vrzI/1+un34vXoAE/Rpv7fae3e18qR395behdd3v5GnSKRpIKZcBLUqEMeEkqlAEvSYUy4CWpUAa8JBXKNskdjLun1p5d6bGq3yGpo3od6VG+IzsYd0+tPbvSY23+DslONl9H6nGKRpIKZcBLUqEMeEkqlAEvSYUy4CWpUAa8JBWq022S4+ohr9NT23St9+r2yyurrWuTStK0331Ux6pzzVZzoLp99fnDT3jcI7kzKIM258lufsel0wE/rh7yOj21Tdd6b9rXK+0Fu3ldNL1mN+fAAwOe3ymDqtvXPfaoOEUjSYUy4CWpUAa8JBXKgJekQhnwklQoA16SCtW5Nsk6fbOD+l0H9c1XH4+rpqbb72Z/sNRFw1wjTXvih3n9SepcwNfpmx3U71qnf3W71x2mpmHHIOlMw1wjTXvid7O2UXKKRpIKZcBLUqEMeEkq1FBz8BHxWuBa4ADw/sz805FUJUkaWus7+Ih4EvBHwMXAc4CjEfHToypMkjScYe7gLwWOZ+b/AkTEJ4BXAu/eYb99ALOzM60P/MSzDzV6XD1Wne3bHGMvPJ6WOqbt8bTUMW2Pp6WOphkxqP6mObLd+9Im/yr71O7pnllfX298IICI+D3g8Zl57cbPVwIvzMyjO+x6MfClVgeVJF0CfLnOhsPcwc8C1X86zABrNfa7m16BJwD/yxeSVM8+4Fx6GVrLMAH/PXpB3XcOcH+N/U5T858+kqQz/GeTjYcJ+H8E3hURi8CDwCuAnaZnJEm7pHUXTWZ+H/h94AvAN4C/y8x/HlVhkqThtP6QVZI03fwmqyQVyoCXpEIZ8JJUKANekgrVuf/gR1VE/BLwgcz8qUnX0lREXAFcT+/LC3cDRzNzebJV1RcRl9Orfwa4D/i1zPzBZKtqJyL+mt6yG3816Vp2UsoCfxFxGLgLeHlmfnfC5TQSEdcBr9748Vhmvm2S9WxnagM+Io4C5216ei4z37Hx+58A/oRewEylQWOgt0jbB4DnZeb/RMTHgDcAf7G7FW5vm/rfC/wZcGFmfj8i3g28C7h6dyvc2XbnUUScB/w58FLg+K4X11Blgb/n0/vC4F0R8YXMvHeylTUTES8CPgRcMOlamoqIS4HLgOfS+yb/ZyPiisy8fbKVbW1qAx64KzP/vfpERPxM5ccP07uDfO+uVtXMlmPIzAcj4vzMXImIxwFPBKbx7nfQ3+AAcNXGdyEA7gFet9vF1bTdefQ64FPA0q5X1U7bBf6mzZuAq4C/mXQhLZwArun/23ZEfBt46mRLGmyaA36giPgd4F+Bf5p0LW1thPvLgFuB7wOfm3BJtWXmEnA7QEQcAt4O3DjRolrIzD8GiIiLJ11LTefRC5i+E8ALJ1RLa5l5JUBETLqUxjLzW/3HEfF0elM1F02uou117kPWjbuvVwB/MOlahpWZ/5CZR4A76E15dEpEPAE4BnwzM2+ZdD17QNsF/jRiEfFM4PPAWzPzO5OuZ5DOBTzwKnorqn0d+AxwXkR0avnhiPjxiLis8tTfAj87qXraiIhz6S37fA9w5YTL2Su+R+/c76u7wJ9GKCIuAu4E3j7tNzadm6LJzOuA6wAi4nzgi5l5ybY7TZ8Z4NaIeEFm/je9f2h1ZoXNiNgHfBr4+8z8w0nXs4e4wN+ERcRTgE8Cr8nMqf9gfpoD/sUR8cpNz80B75hEMS1tOYaNDo6jwB0RsQ7cC7x598vb0aC/wd3A84D9ld9/vT+3OmVKOI+A3gJ/EdFf4G8O+LAL/O26twAHgRsqnyHclJk3Ta6kwVxsTJIK1cU5eElSDQa8JBXKgJekQhnwklQoA16SCmXAS1KhDHhJKpQBL0mF+n8Zi2QLzrjV8wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(samples.detach().cpu().numpy(), bins=100);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "18.617748260498047"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R(samples).item()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Simple way to change `mu`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "conditioned_y_sample = poutine.condition(y_sample, data={\"mu\": torch.tensor(-2.)})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "samples = torch.stack([conditioned_y_sample() for _ in range(1000)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEBCAYAAABysL6vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAEuBJREFUeJzt3X+MZXV5x/H3zLK/kEHoZCigWGqUR2sVfwAmAlqVktDabAiiKU0sjUiNtKUJaKzFwNI22lg3JliLVdfS0ootiAaWGi1Yo2KtWpVU9AlpMY0yaTZToyzq7uzM9I97L5wd5u7c3+fMd96vhHjumXPv+Tj33GfOfs9zv2dqZWUFSVJ5pusOIEkaDwu8JBXKAi9JhbLAS1KhLPCSVCgLvCQVygIvSYWywEtSoSzwklQoC7wkFcoCL0mFOqaGfW4HzgbmgaUa9i9JG9EW4BTgq8DBXp5QR4E/G/hCDfuVpBKcD3yxlw3rKPDzAD/84WMsLzdrJsvZ2eNYWDhQd4wnaWouaG62puaC5mYzV/8mmW16eooTT3wKtGtoL+oo8EsAy8srjSvwQCMzQXNzQXOzNTUXNDebufpXQ7aeh7a9yCpJhbLAS1KhLPCSVCgLvCQVygIvSYWywEtSoSzwklSoOvrgpXXNHL+THdtbh+fPDh7m0R//tOZE0sZjgVcj7dh+DL9xzacAuOu9u3i05jzSRuQQjSQVygIvSYWywEtSoSzwklQoL7JqQ7G7RuqdBV4bit01Uu8copGkQlngJalQFnhJKpQFXpIKZYGXpEJZ4CWpUBZ4SSpUT33wEXEj8FpgBfhIZu6JiAuAPcBO4OOZed34YkqS+rXuGXxEvAJ4FfAC4Czg9yPiTGAvsAt4LnB2RFw0zqCSpP6sW+Az8/PAKzPzMHASrbP+E4CHMvPh9vpbgUvHmlSS1JeexuAzczEidgMPAvcCpwLzlU3mgaePPp4kaVA9z0WTmddHxJ8DdwFn0BqP75gClvvZ8ezscf1sPjFzczN1R1hTU3PBZLJ120dn/aHFJbZt3fL48qRyDaqp2czVvyZnW7fAR8RzgB2Z+c3M/ElEfILWBdelymYnA4/0s+OFhQMsL6+sv+EEzc3NsH9/86avamouGF+21R+azj6Otr46CVn1Z03T1PfTXP2bZLbp6am+T4x7OYN/JrA7Is6jdda+C/gg8J6IeBbwMHAZrYuukqSG6OUi6z3APuAbwNeB+zPzNuBy4A5a4/LfBW4fX0xJUr96GoPPzBuAG1atuxc4c/SRJEmj4DdZJalQ3tFJE78NXnV/Bw8tsX3blontW9pMLPCa+G3wVu/PW/BJ4+EQjSQVygIvSYWywEtSoSzwklQoC7wkFcouGk1EtTVS0mR4Bq+J6LRGdloiJY2fBV6SCmWBl6RCWeAlqVAWeEkqlG0NOsKhxaXH75pUnQgMmjcZWDWrpCezwOsI27ZuWXMisM7jJk0GtjqrpCM5RCNJhbLAS1KhLPCSVCgLvCQVyousm0i3W+U1hV0x0mhZ4DeRo90qrwnsipFGyyEaSSpUT2fwEXE98Lr2w32Z+baI+ChwHvBYe/3uzLxzDBklSQNYt8BHxAXAhcCLgBXg0xFxMXAW8PLMnB9vREnSIHo5g58HrsnMQwAR8R3gGe3/9kbE04A7aZ3BL48tqSSpL+sW+Mz8dmc5Ip5Na6jmfOBXgLcAPwLuBt4IfKjXHc/OHtdn1MloahdHU3KtleOpJxzLtq2tjpxDi0uPL4/KoN01TfmdraWp2czVvyZn67mLJiKeB+wD3pqZCVxc+dlNwBvoo8AvLBxgeXmlj6jjNzc3w/79TZptpWVUuUZxIK7OMTc386Tul7WyDrPvQbtrmvheQvnH2ag1NRdMNtv09FTfJ8Y9ddFExLnAvcDbM/OWiHh+RFxS2WQKWOxrz5KkserlIutpwCeB12fmfe3VU8D7IuI+4ABwJXDL2FJKkvrWyxDNtcAOYE9EdNbdDLwL+BKwFbgjMz82loSSpIH0cpH1auDqLj/+wGjjSJJGxW+ySlKhnItGA6lOXFZVbWls2i3+qqr5m5xTGoYFXgNZPXFZx+qWxmY2tz05f1NzSsNwiEaSCmWBl6RCWeAlqVAWeEkqlBdZ1bONcku9zoRnc3MzdshoU7PAq2cb5ZZ6G6WTRxo3h2gkqVAWeEkqlAVekgplgZekQnmRVUXrd24c56hRSSzwKlq/HTXOUaOSOEQjSYWywEtSoSzwklQoC7wkFcoCL0mFsotG6oHtk9qILPBSD2yf1EbUU4GPiOuB17Uf7svMt0XEBcAeYCfw8cy8bkwZJUkDWHcMvl3ILwReBLwQeElE/CawF9gFPBc4OyIuGmdQSVJ/ernIOg9ck5mHMnMR+A5wBvBQZj6cmYeBW4FLx5hTktSndYdoMvPbneWIeDatoZqbaBX+jnng6SNPJ0kaWM8XWSPiecA+4K3AYVpn8R1TwHI/O56dPa6fzSemqbeka2quo2niLf7WynO0nL2s79wicPVyvzmawFz9a3K2Xi+yngvcAfxhZt4WEa8ATqlscjLwSD87Xlg4wPLySj9PGbu5uRn2729ef8Sock36QGziLf46v8fq7+JoOdfavrq+87Pq89d7r0o/zkatqblgstmmp6f6PjFet8BHxGnAJ4HXZ+Z97dVfaf0ongU8DFxG66KrJKkhejmDvxbYAeyJiM66m4HLaZ3V7wDuAW4fQz5J0oB6uch6NXB1lx+fOdo4kqRRcS4aSSqUUxVIXTSxC0jqh2fwUhed7ppOh4y00VjgJalQFnhJKpQFXpIKZYGXpEJZ4CWpUBZ4SSqUBV6SCmWBl6RCWeAlqVAWeEkqlAVekgrlZGPaNJw8TJuNZ/DaNJw8TJuNBV6SCmWBl6RCWeAlqVAWeEkqlAVekgplgZekQlngJalQPX/RKSKOB+4HXpOZ34uIjwLnAY+1N9mdmXeOIaMkaQA9FfiIeCnwIeCMyuqzgJdn5vw4gkmShtPrEM2bgKuARwAi4ljgGcDeiHggInZHhMM9ktQgPZ3BZ+YVABHRWXUycB/wFuBHwN3AG2md5fdkdva4fnJOTFPnKlkr16HFJbZt3QLAwcUltreXq+s1Hkc7Tno5hjbScdYETc0Fzc420GRjmfnfwMWdxxFxE/AG+ijwCwsHWF5eGWT3YzM3N8P+/Y/WHeNJuuWam5t5fF6Vu96764jlbttrNKq/39W/1/WOoY12nNWtqblgstmmp6f6PjEeaFglIp4fEZdUVk0Bi4O8liRpPAadLngKeF9E3AccAK4EbhlZKknS0AY6g8/MB4B3AV8CHgS+mZkfG2UwSdJw+jqDz8zTK8sfAD4w6kCSpNHwjk4bzMzxO9mx/ehvW/XORT87eJhHf/zTSUST1DAW+A1mx/ZjjuiWWUvnzkWdbZrZfyBp3PxykiQVygIvSYWywEtSoSzwklQoC7wkFcoumsJVWybVPNW2V1taNWoW+MKtbplUs6xue7WlVaPkEI0kFcoCL0mFssBLUqEs8JJUKC+yTojdEpImzQI/IXZLSJo0h2gkqVAWeEkqlAVekgplgZekQnmRtaGqXTeHFpdqTqOqXuf3sXNKdbPAN1Qvt+ZTParz+0D398fOKdXNIRpJKlRPZ/ARcTxwP/CazPxeRFwA7AF2Ah/PzOvGmFGSNIB1z+Aj4qXAF4Ez2o93AnuBXcBzgbMj4qJxhpQk9a+XIZo3AVcBj7QfnwM8lJkPZ+Zh4Fbg0jHlkyQNaN0hmsy8AiAiOqtOBeYrm8wDTx95MknSUAbpopkGViqPp4Dlfl9kdva4AXY9fpO6vV1nP4cWl9i2dcuTluvKpeF1e6+q69fbptfjYtSaepw1NRc0O9sgBf77wCmVxyfzxPBNzxYWDrC8vLL+hhM0NzfD/v3jaWZbfRB09jM3N3NEK111fbft1Wy9vIe9bLPWcTFO4zz+h9HUXDDZbNPTU32fGA9S4L8CREQ8C3gYuIzWRVdJUoP03QefmT8DLgfuAB4EvgvcPtpYkqRh9XwGn5mnV5bvBc4cRyBJ0mj4TVZJKpRz0Ugj1OtEZIO+ZnXSMicz03os8NIIVSciG9Ukcatfs9Oz4WRmWo9DNJJUKAu8JBXKAi9JhbLAS1KhvMi6AYyjM0NS+TyD3wA6XRTV28RJ0nos8JJUKAu8JBXKAi9JhbLAS1Kh7KJpELtlylV9bw8tLtWcRpuFBb5BxjGPiZrB91Z1cIhGkgplgZekQlngJalQFnhJKpQFXpIKZRdNDWyH3Ny6vf/DHBfV5x48tMT2bVsAb+W32Vnga2DL3ObW7f0f5rhY/Vxv5ScYssBHxOeAk4DF9qrfzcyvDJ1KkjS0gQt8REwBZwC/kJmHRxdJkjQKw1xkjfb/fiYivhURvzeKQJKk0RimwJ8I3AtcDLwaeHNE/OpIUkmShjbwEE1mfhn4cudxRHwE+DXgs708f3b2uEF3PVZH62I4tLjEtq2t7oSDi0tsby9X10vjNEinTbeOnbWO5eNPOLaRx3WTu86anG2YMfjzgO2ZeW971RRPXGxd18LCAZaXVwbd/VjMzc2wf3/3noO5uZmunQpHe17nudKwBum0WevYPNqx3M9xPQnrfS7rNMls09NTfZ8YD9NFcwJwY0S8DNgK/Dbw5iFeT5I0QgOPwWfm3cA+4BvA14G97WEbSVIDDNUHn5nvBN45oiySpBFyLhpJKpRTFYzYzPE72bG99Wt1HhDVrdscNdocLPAjtmP7Mc4DosY42hw1Kp9DNJJUKAu8JBXKAi9JhbLAS1KhvMi6jmpXjKTB2F1WDyvXOlZ3xUjqn91l9XCIRpIKZYGXpEJZ4CWpUBZ4SSqUBV6SCrWpumh6bdXqtzWy223UBrm9mtRE1c9EddKy6nIT2x83e3vmpirwvbZq9dsa2e02aoPcXk1qotWfiW7LTWt/3OztmQ7RSFKhLPCSVCgLvCQVygIvSYXacBdZu10V7+Uqf1W3W5nZ9SK1jKsLrPq63Tpbqp/nQ4tLI88wSau78ibZzbPhCny3q+K9XuXv6HYrs9XbSZvVuLrAVr/uWp0tJU3yV/3/ApPt5nGIRpIKNdQZfERcBlwHbAXel5l/OZJUkqShDXwGHxFPA/4MOA94IXBlRPzSqIJJkoYzzBn8BcB9mfl/ABFxO/Ba4MZ1nrcFYHp6auAdn3TizseXq69TXT/M8ihfq6TlpuRo2nJTcjRh+Wif636f0+vrrmdUr9NNL6+5+ngZJEflOU/uGuliamVlpe8dAUTEHwFPyczr2o+vAM7JzCvXeep5wBcG2qkk6Xzgi71sOMwZ/DRQ/eswBSz38Lyv0go4D2zs/idJmpwtwCm0amhPhinw36dVqDtOBh7p4XkH6fGvjyTpCP/Vz8bDFPh/AW6IiDngMeASYL3hGUnShAzcRZOZPwD+GPgc8E3gHzLz30cVTJI0nIEvskqSms1vskpSoSzwklQoC7wkFcoCL0mF2nDTBU9KRPw68P7M/MW6swBExMXAblpfdvgqcGVmHqo3VUtE7KKVbQp4GPidzPxhvameEBF/S2tajb+pOUdjJ+eLiOOB+4HXZOb3ao7zuIi4Hnhd++G+zHxbnXk6IuJGWlOzrAAfycw9NUda06Ys8BFxJXDqqtXbMvMd7Z//PPAXtApW7bloTer2fuDFmfm/EXEbcDnw1w3I9m7gr4CzM/MH7QP/BuDqOnNl5jsi4lTgg8CrgfsmkaebyuR8L6H1Zb/7I+JzmflgnbkAIuKlwIeAM+rOUhURFwAXAi+iVUg/HREXZ+adNed6BfAq4AW0/lg/GBH7MjPrzLWWTVnggfsz8z+rKyLilysPP0zrjPTdE03VJVdmPhYRp2fmYkQcC5wETPoMudvvbCtwVft7EQAPAL/VgFy0c3wKWJhgnm4GnZxvEt4EXAX8Xd1BVpkHrun8SzUivgM8o95IkJmfj4hXZubh9h/uY2h92bNxNmuB7yoi/gD4D+Df6s5S1S7uFwG3Aj8APlNzJAAycwG4EyAidgJvB26qNVRbZr4HICLOqzsLrX9lzFcezwPn1JTlCJl5BUBE1B3lCJn57c5yRDyb1lDNufUlekL787gbuBb4J1qfycbxImtF+8zvEuBP6s6ylsz858ycBe6mNSzSGBHxVGAf8K3MvKXuPA006OR8m15EPA/4LPDWzHyo7jwdmXk9MAecRutfQY1jgT/SpbRma/sacA9wakTUPrVxRPxcRFxYWfX3tMb/GiEiTqE1BfQDwBU1x2mq79M6tjp6nZxvU4uIc4F7gbc35cQhIp4TES8EyMyfAJ+gQZ/HKodoKtp/ka8HiIjTgX/NzPOP+qTJmAJujYizMvN/aP0hasSMnBGxBbgL+MfM/NO68zSYk/P1KSJOAz4JvD4za71Ivsozgd3tob8VYBewt95Ia9usBf5lEfHaVeu2Ae+oI0zFmrnaHSFXAndHxArwIPDmJmSj1bL5YuCYys+/1hnXrTFX3e/lEdodRp3J+bYBH3ZyvnVdC+wA9lSuD9ycmTfXFwky856IOAf4Bq17WtyRmbfVmakbJxuTpEI5Bi9JhbLAS1KhLPCSVCgLvCQVygIvSYWywEtSoSzwklQoC7wkFer/AdiAdyP0zbDvAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(samples.detach().cpu().numpy(), bins=100);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 100/100 [00:15<00:00,  6.34it/s]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "rs = []\n",
    "mus = np.linspace(-100, 100, 100)\n",
    "for mu in tqdm.tqdm(mus):\n",
    "    conditioned_y_sample = poutine.condition(y_sample, data={\"mu\": torch.tensor(mu)})\n",
    "    samples = torch.stack([conditioned_y_sample() for _ in range(1000)])\n",
    "    rs.append(R(samples).item())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fde191f2630>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEBCAYAAACdctWRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X+UXOV93/H3rtAvW5KRlk0FRrJlGb7Q9c8InKrI1sYmPlFCMRAbJ7gWBPPDB5yDGwN1MS7g2LEdAnVjkO3KBonQAw6WUAxCnJyCUfghJynlx6kQ37hUtKotNctKtRAWi9Bu/5g7y+xoZue5s/feuffO53UOh907d2a+88zqe5/7fZ773J6xsTFERKScejsdgIiIpEdJXkSkxJTkRURKTEleRKTElORFREpMSV5EpMSU5EVESkxJXkSkxJTkRURKTEleRKTElORFRErsqA6850zgVGA3cLgD7y8iUkTTgGOBfwBGQp/UiSR/KvBoB95XRKQMPgg8FrpzJ5L8boB9+15hdDT+Cph9fXMYHj6QeFBTpbjiUVzx5TU2xRVPu3H19vYwf/6bIcqhoTqR5A8DjI6OtZXkq8/NI8UVj+KKL6+xKa54phhXrDK3Bl5FREpMSV5EpMSU5EVESkxJXkSkxDox8CoiJfLIk7tYd/92hveP0DdvJuesXMrygYWdDksisZK8mf05cIy7X2Bm1wEXAvuih9e6+61JBygiydu2fQ8bt74w5cS8bfse7njQGTlUmfAxvH+E9VueB1Ciz4ngJG9mHwHOBzZHm04Bft/dt6URmIikY9v2Pazf8jyvvT4KTC0xb9z6wniCr3rt9VE2bn1BST4ngmryZrYA+BrwpzWbTwGuMbNnzewWM5uVRoAikqyNW18YT/BV1cQc1/D+xlfXN9su2QvtyX8P+BKwCMDM5gBPAVcB/wNYB3w52idIX9+cOHFO0N8/t+3npklxxaO44guJ7ZEnd3HHlh28tO8gx8yfzepVJzO4bNH443ubJOC9+0eOeP1Wr9U/fzZD+w4eGef82bloxzzE0EiWcbVM8mZ2EbDL3R8yswsA3P0A8Ds1+9wE3EaMJD88fKCtq776++cyNPRy7OelTXHFo7jiC4mtvhQztO8g3/6rp9n/8qvj5ZMF82Y27GkvmDdzwuuHvNZZK5ZMqMkDzDiql7NWLOl4O+b1u2w3rt7enrY6xyE9+U8Cx5rZ08ACYI6Z3Q486u63Rfv0AIdiv7uIJGqyUkw1MZ+zcumE5A2VxHzOyqWxX2v5wELmzZ1V+Nk1SQ5EJ/E6SWqZ5N39t6o/Rz35QeBqYIeZ/QR4EbgcuDeVCEUkWEiNvJp0WiWj0Hr74LJFDCw+eiphd1RSA9FJDmgnqa158u4+ZGaXAvcBM6gse3lTkoGJSHx9TUoxffNmTvh9+cDClokn9LWKLuSMJcvXSVqsJO/u66gMsuLuG4ANyYckIu0KLcVk/Vp5ltQMobzONNIVryI5kUQ9N7QUk/Vr5VlSZyx5PfNRkhfJgSTruSGlmE68Vl4ldcaS1zMfJXmRHMhrPbcbhJ6xtDrTyuuZj5K8SA6E1nPzOEUvRN7jbnXGEnqmlcczHyV5kRwIqec+8uSuXE7RayU0Qeb5QFDkMy2tJy+SA+esXMqMoyb+c6yv596xZUdia85kKWStnOqBoHqgqx4Itm3fk2mszeR15kwIJXmRHFg+sJDzV5003nPvmzeT81edNKGX+FKDNWIg/4kmJEEmuWhaGprNkOn0zJkQKteIxJRWWaFVPfeYJouB5T3RhJSi0uwpJ/F95XXmTAj15EVi6GRZYfWqk1uWdPIopBSVVk85qe8r5Ewrr9STF4mhkwNwg8sWsf/lV3M7ONlMyNTCtHrKSX5feZw5E0JJXiSGTg/AFTXRtIo7rTnmnf6+8kBJXiSGvF66XgZpHMD0fakmLxJLSH1Z8kPfl3ryIrG0W1bI84U+ZZbXpQaypCQvElNIWaE2qb951jRGDo3y+uHK7S6LcqVqEYQcPIs6jpEUJXmRhNVfxv/Kq4eP2Kcol8TnUW1ir6WDZ2OqyYskrNG0vUa6aYZHUurnvdfL01WyeaEkL5Kw0OTdTTM8khJyANXBcyIleZGEhSTvbpvhkZSQBK6D50RK8iIJazRtb1oPzJldGQIr0iXxedMqgevgeSQNvIokTNP20tNo+YMqtXNjwUnezP4cOMbdLzCz9wHfB+YBfwt81t1fTylGkcLp9ml7adEBNL6gJG9mHwHOBzZHm+4ELnL3n5rZD4CLge+kE6JIeekiqfh0AI2nZZI3swXA14A/Bd5rZm8DZrv7T6Nd1gE3oCQvEkvobfEkvm3b97DpsW0M7TvY9QfPkIHX7wFfAvZFvx8H7K55fDdwfMJxiZRe3u+GVFTVg2f1Bit5u5Vg1ibtyZvZRcAud3/IzC6INvcCYzW79QCtr/yo09c3J+5TxvX3z237uWlSXPF0e1x7m0wH3Lt/pGkM3d5mITY9tq3hwXPTYzs5c/CEDkU1UZbt1apc80ngWDN7GlgAzKGS4I+t2Wch8Iu4bzw8fIDR0bHWO9bp75/L0NDLsZ+XNsUVT1nimkpNfUGTZXAXzJvZMIaytFnaGt0isbo9D3G22169vT1tdY4nLde4+2+5+7vc/X3Avwd+7O5/CLxqZqdFu30a2BL7nUUKbqq3ltMyuOko8k2309DuxVCfAv6DmT1PpXf/F8mFJFIMU62pF/m+oXmmg+dEwfPk3X0dlZk0uPszwAfSCUmkGJK4tZymAzY2lTJYdb9Nj+3U7Bp0xatI23RruXQkMbV0+cBCzhw8IRc1+E5TkhdpU6NL7NMoC1R7tXv3j7CgC3qlk5XBkvzc3XIhmpK8SJuyuMS+Gy+YSqIM1ko3tauSvEiNR57cxbr7twcn7bRr6ln1avMkizJYN7WrlhoWiWzbvodb7nmm7SmRaciiV5s3WcyO6aZ2VZIXiWzc+gIjhybej7XTywx045zvLKaWdlO7qlwjEslj7y6rwd0shQx4pl0GK2O7NqMkLxLJ45TI2sHdMsyuycuAZzetS68kLxI5Z+VS7njQjyjZDO8f4ao1j3csCVR7tXlbI6YdeRrw7JYL0VSTF4ksH1jI5z7x3oY99zwMwpZBHktiZackL1JjcNkibrzstIaJvtODsGXQTQOeeaFyjUgD6nGmo2gDnmW4KlZJXqSBPA7ClkGRBjzzMkg8VUryIg0UrcdZJEUZ8MzTIPFUKMmLNFCkHqekoywlOyV5kSaK0uOE/NSO8xJHEspSstPsGpGCm+ptCMsWR1LKcocp9eRFCi4vteO8xJGU0JJd3s9elORFCi4vteO8xJGkViW7IszAUblGpODycoFRXuLI0lRv5p4FJXmRgstL7TgvcWSpCGcvQeUaM/sK8HFgDPiBu99sZrcDK4BXot1ucPd70wlTRJrJy3TPvMSRpSLMwGmZ5M1sJfBh4D3AdOA5M9sMnAJ8yN13pxuiiLSSxXTPPKwDnzdFuGiuZbnG3bcCv+nurwO/RuXAcBBYDNxmZs+a2Q1mptKPSEmVbXpkUrK4i9VUBZVr3P2Qmd0AXAncQ6VH/zBwGfBL4H7gM8DalOIUkQ6abIDxzMETOhRVPuT97KVnbGwseGczexNwH/BDd/9PNdvPBla7+9kBL/N2YGfMOEWkg878wl/TKFP0AD++6WNZh9PtlgAvhu4cUpM/CZjl7k+7+6/MbCPwSTMbdvcN0W49wKE4UQ4PH2B0NPwAU5XXu+MorngUV3ydjG1BkwHGBVGZIo9tltfvst24ent76OubE/t5IeWadwA3mNkKKrNrPgZsBb5lZg8DB4BLgPWx311EMhN6ZWaj/YowwCiNhQy8PgBsBp4CngSecPevAF8HHgeeA55297vSDFRE2hc6cNpsPyD3A4zSWOjA6/XA9XXb1gBrkg9JRJIWuq7MZPvdeNlpSupTUD1D2rt/hAUZXkOgtWtEukDolZlFuIKziDq5xo2SvEgXCL0yswhXcBZJ7fhGvaxW6NQFTCJdIHRdmW5cfyYt9eMbjWRxhqSevEgXCF1XphvXn0lLo/GNelmcISnJi3SJ0Csz834FZ1G06qVndYakJC/S5fJ+Z6Oiaja+UX1Ms2tEJHXtzvqoPTD0z5/NWSuW6MBQp9kFZH907vsYWHx0ZnFo4FWki7VzZ6P6AcWhfQe1ImUDzVaoHFy2KNM41JMXKbFWpZh25sWX7YbdacrD+IZ68iIlFbKUQTv3ZdUFU8WinrxISYX0uNtZeEwXTBVrsFpJXqSkQnrc7cyLDz0wFCkRxtHJJQraoSQvUlKhPe64deP6A0Oj2TVFS4RxFG1MQklepKTSXAO+9sDQ6CYYRUuEcRRtTEJJXqSkOrlEQdESYRxFG5NQkhcpsU5N4StaIoyjaHfJUpIXCVDWQcS0FC0RxlG0RdyU5EVa6MQgYqODypmDc1N5rzSEJsI8HjxDYsrDRU6hlORFWsh6ELHZQWXe3FkT1jzJY4Ks1SoR5nEGTh5jmipd8SrSQtaDiM0OKnds2TH+e+iNufOsnXVz0pbHmKZKPXmRFkIHEZPqWTc7eLy07+D4z2WYohh68Axt1/r9LjhjIPZqj2WcFaSevEgLIbfES7Jn3WwGyjHzZ4//XIZkFLJuTmi7Ntrvlnueid3+7azlk3dBSd7MvmJmz5nZdjP742jb6Wb2rJn9zMy+mm6YIp3TbMnY2t5knNP8bdv3cNWax7nwGw9z1ZrHj0hEzQ4qq1edPP57GZJRyMEztF0b7Tdy6HDsMksZ73HbslxjZiuBDwPvAaYDz5nZQ8BtwEpgF7DZzFa5+5Y0gxXplFaDiHFKD60G9prNTBlctmj8ytIyTFEMmYET2q5JlX6KNj0yRMsk7+5bzew33f11M3tr9JyjgZ+5+04AM7sT+ASgJC9dKbRuH1pLb3VQKUsyavU5Q9s1ZL/QmTNFmh4ZImjg1d0PmdkNwJXAPcBxwO6aXXYDx8d5476+OXF2n6C/P5/zhRVXPGWK64IzBrjlnmcYOXR4fNvM6dO44IyBCa+3t0mPc+/+kaD3rd3nzMG5nDl4QuxY05DWdxnariH7bXpsW8MD7KbHdmbejln+7QfPrnH368zsm8B9wInAWM3DPcBowyc2MTx8gNHRsdY71mm0GFIeKK54yhbXwOKjWf3bdkTPemDx0RNeb0GTHueCeTNbvm/Z2ixEaLs22q86u6a631DN7KRaQ/sOZtqu7bZXb29PW53jkJr8ScAsd3/a3X9lZhuBjwOHa3ZbCPwi9ruLlEjIaX4ZaulZCy2f1O9Xn0zLvJ7OZEJm17wDWGtmM81sBvAx4HuAmdk7zWwacB6qx4u0FDJTR9JRxpkzIUIGXh8wsw8AT1HpvW9w97vNbAjYAMwCHgB+lGqkIiVRtoG9oijLYHVcoQOv1wPX1217CHhv8iGJiKSjGw+wuuJVRKTElORFREpMSV5EpMSU5EVESkxJXkSkxJTkRURKTEleRKTElORFREpMSV5EpMSU5EVESkxJXkSkxJTkRURKTEleRKTElORFREpMSV5EpMSU5EVESkxJXkSkxJTkRURKTEleRKTElORFREpMSV5EpMSOCtnJzK4Dzo1+3ezuV5vZ7cAK4JVo+w3ufm8KMYqISJtaJnkzOx34KPB+YAx40MzOBk4BPuTuu9MNUURE2hXSk98NfMHdXwMwsx3A4ui/28zsrcC9VHryo6lFKiIisbVM8u6+vfqzmZ1ApWzzQWAQuAz4JXA/8BlgbSpRiohIW3rGxsaCdjSzAWAzcJ27r6977GxgtbufHfBSbwd2xoxTREQqlgAvhu4cOvB6GrAB+Ly7321m7wZOdPcN0S49wKE4UQ4PH2B0NOwAU6u/fy5DQy/Hfl7aFFc8iiu+vMamuOJpN67e3h76+ubEfl7IwOsiYBPwSXd/ONrcA3zLzB4GDgCXAOubvISIiHRISE/+SmAWcLOZVbd9F/g68DgwHdjg7nelEqGIiLQtZOD1CuCKJg+vSTYcERFJkq54FREpMSV5EZESU5IXESkxJXkRkRJTkhcRKTEleRGRElOSFxEpMSV5EZESC1q7RrK1bfseNm59geH9I/TNm8k5K5eyfGBhp8NKXMjnDG2LbmkzkbiU5BPyyJO7WHf/9iknrG3b97B+y/O89nplaf7h/SOs3/I8QKmS1iNP7mr5OUPbIs5+OhBIt1G5JgHbtu/hlnueYXj/CPBGktm2fc+EfdZveX7SfQA2bn1hPFlVvfb6KBu3vpDyp8jWHVt2tPycoW0Rsl9o+4uUjXryAVr1ADdufYGRQ4cnPKeaZKr7TZaIal+rmoTqNdueV63a7KV9Bxs+r/ZzhrZFyH6h7S9SNkryLYSUAkKSTGjC6ps3s+G+ffNmHhFXXksPIW12zPzZDDVI9LWfM7QtQvYry8FTJC6Va1oIKQXUJ51G20P2AThn5VJmHDXxa5lxVC/nrFw6/nveSw8hbbZ61cktP2dIW4TuF9r+ImWjJN9CSA/wnJVLmTl92oTH201YywcWcv6qk8aTT9+8mZy/6qQjykOdrNtv276Hq9Y8zoXfeJir1jx+xMElpM0Gly1q+TlD2iJ0v9D2FymbUpdrkpiiF1IKWD6wkHlzZ006u6a2Nt+qxLJ8YOGkpZdOlh5CSjGhZZZWnzN0n5D94rS/SJmUNsmHJKOQfc5ZuXTCPtC4Bzi4bBEDi4+eNKbQhNVKnLr9pse2MbTvYGJJLWQAM7TNspZU+4sUSWnLNSEljZB9QksGWYpTt68ObiZVtw85i8hjm4l0q9L25JOc8ZK3HmBI6SGtKYNJlmJEJH2lTfIhySg0YeVRp+r2eS3FiEhjpS3XhJQ0yjzjIq0pgyrFiBRLUE/ezK4Dzo1+3ezuV5vZ6cDNwGzgh+5+bUoxHiFk1kxISaPMMy7a6XGHXmClUoxIcbRM8lEy/yjwfmAMeNDM/gD4JrAS2AVsNrNV7r4lzWAhbGGrqiSn6BVN9TNtemznpLNrahN7rbIujCbSbUJ68ruBL7j7awBmtgM4EfiZu++Mtt0JfAJIPclPtrCVktFEywcWcubgCQwNvdzw8foppPXUriLF1zLJu/v26s9mdgKVss23qST/qt3A8YlH10DIwlYSptEMnHpqV5FiC55dY2YDwGbgKuB1Kr35qh5g8mxRp69vTpzdeeTJXdyxZQdjTR7vnz+b/v65sV4zaZ1+/2aaxbU3IIGn2a5Fa688yGtsiiueLOMKHXg9DdgAfN7d7zazlcCxNbssBH4R542Hhw8wOtosZU/Uqqww46hezlqxpGlZIgv9/XM7+v7NTBbXgiZTSKtq2zXpVS+L2F6dltfYFFc87cbV29sTu3MMAVMozWwRsAk4z93vjjb/XeUhe6eZTQPOI8V6/GRlBU3ha1+jKaRVte2a91UvRaS5kJ78lcAs4GYzq277LnABld79LOAB4EcpxAdMXhe+8bLT0nrb0gudQqobbogUV8jA6xXAFU0efm+y4TQW58rUPN9MI49CppDqhhsixVWIK15Dr0xVWSEduuGGSHEVYu2a2rLC3v0jLFBZIbapnOFovRqR4ipEkoc3ygqTjUx3a1mhVQIPWTd/MmVe/kGk7AqT5EMUeVXJdoUk8CTOcMq6/INI2RWiJh+qzKtKNhNy45NuPcMRkZL15LuxrBCSwLvxDEdEKkqV5CGsrFCmaZYhCVwDpyLdq1TlmhBlm2YZUqLSjT5EulfpevKtlG2aZWiJSgOnIt2p65J8GQch85LAy1QGEymLrkvyGoRMx1Tn4ndS7cGpf/5szlqxJPcxi4Tqupp8N06zzELIVM48qh+jGdp3sNBjNCL1uq4nX7RplkUpgYSWwaqfZ7LlKbLQ7N62UOwxGpF6XZfkoTjTLItUAgkpg+Xl87S6CQ0Ue4xGpFbXlWtC5GWaZZFKICFlsLx8npB722qMRsqiK3vyrYROs0x7wK5IM4FCymB5+Tyt3k9jNFImSvINhCSj+lP+6oAdhJceWpWEijYTqFUZLC+fp1kcgGbXSOmoXNNAyE0yplp6CCkJlW0mUF4+T7M4Lv5X/5zbrv2oEryUipJ8AyHJaKqlh5CDRNmWI6j9PD107vOUrV1FJqNyTQMh9eWplh5CDxJ5uZo1KSE3f8liZlPZ2lWkGSX5Jlolgamu7JiX+nTeJDHNMg/TX0XyQuWaNtWf8vfPnx3rlD8v9em8yWKsQ6SbBPfkzWwe8ARwhru/aGa3AyuAV6JdbnD3e1OIMbdqe/uTlR+aPReKc+VtVtIc6+j2tpXuFJTkzew3gLXAiTWbTwE+5O670wisyELLBaoLHymrsQ6RbhHak78YuBz4SwAzexOwGLjNzN4K3EulJz/5ZYRdoFFNee19z7H2vufUWw8QOtbR7ECqsQ6RiXrGxsaCdzazF4FBKrX8m4DLgF8C9wN3ufvagJd5O7AzXpjFceFX/4ahfQebPj5z+jQ+94n3MrhsUYZRFcsjT+7iji07eGnfQY6ZP5vVq05mcNmi8e2N2rfargC33PMMI4cOH/GY2lxKYgnwYujObSV5d3+xbvvZwGp3PzvgZd4O7BwePsDoaPh7V8WtfWelGteF33i45b5982Zy42WnZRBV/tsrVMiiYtV2ncrsmry2F+Q3NsUVT7tx9fb20Nc3B2Im+bamUJrZu4ET3X1DtKkHONTOaxVNu0sR1FJ9OL6QRcWq7Vo71lH9vlQuk27V7jz5HuBbZvYwcAC4BFifWFQ5Ndkc7jMH5wKNa8r1VB8+UquDZ8iBsb5d87K0sUgntZXk3f1ZM/s68DgwHdjg7nclGlkOTTY978zBE4Ajp0bWazYXvpsv4AlJxq3OkBq1q6ZTisRM8u7+9pqf1wBrkg6ok9rtTU62FEFI8u72HmdIMp7sDKlZu2o6pUgXLWvQKtlOpTc5WfklZC58t/c4Q5JxOxePaTqlSJck+ZAE3m5vMomlCLq9xxmajONePJbW9yVSJF2R5EMSeFq9yRBF7nGGlqPq96kOVEN6yVhLR4iUIMmHJJmQBJ5WbzJEkkmuvj0uOGOAgcVHJxnuhPdqdYbUbJ95c2eNx5VmMtbSEdLtCp3kQwcsQxJ4J0/tk0pyjdrjlnueYfVvWyqJLuQMqdk+d2zZwTcvXT6+TclYJB2FTvKhA5YhCbzTp/YhSa7VWUuj9hg5dHjSG5BP9jmTmG3UbJ+XJln6QUSSU+gkH2dKI7RO4HnuTYactbRzA/JmZz9JzTZqts8x82eHfGwRmaJC3zQk5IbbVcsHFnLjZadx2xc/zI2XnZbbZN5MyM00krwBech+ITc+abbP6lUnN4xVRJJV6CTfTXdXCumlN2qPmdOntXUD8tDZRq1uiN1sH60IKZKNQpdrOl1Hz1JIaaRRe9TPrgmdRZTkbKM8l8FEyq7QSR66J4GEzv6pb4/6ZU1DX0cXEomUQ+GTfLdI6qwlziB0Eu8nIp2lJF8gSZ21hL5Ot5wliZRZoQdeRURkckryIiIlpiQvIlJiSvIiIiXWiYHXaVC583i7pvLcNCmueBRXfHmNTXHF005cNc+ZFud5PWNjY7HfbIpWAI9m/aYiIiXxQeCx0J07keRnAqcCu4HDWb+5iEhBTQOOBf4BCL5tXCeSvIiIZEQDryIiJaYkLyJSYkryIiIlpiQvIlJiSvIiIiWmJC8iUmJK8iIiJVaI9eTN7A7gYXdfF/2+GLgT+DXAgU+5+wEzOxr4z8A7gCHgXHffk3Js7wPW1WzqB/a5+7vM7HzgG8D/jR7b7O5fSjOeutiuAy4E9kWb1rr7rVHM3wfmAX8LfNbdX88wro8BNwA9wE7gD919Xw7a6zzgWmA68C13vzWr924Qy3XAudGvm939ajO7ncoV469E229w93szjusnVP7dHYo2XQospYPtZmYXAZ+r2bQE+EvgzXSovcxsHvAEcIa7v2hmpwM3A7OBH7r7tdF+qf9b7HiSN7NLgOPqNs9w92vM7Djge8BHgIdrHl8DrHH3u83sy8CXgX8LfBV41N1/18w+DfxH4JNpxwm8L9rnTcDfA5+NHj8F+GN3vyuJGNqI6xTg9919W93jdwIXuftPzewHwMXAd7KIi0oS/w5wqrv/3My+AlwPXEEG7dWMmb0V+BqwjMrVhE+Y2U/c/bkOxHI68FHg/cAY8KCZnU2lfT7k7ruzjimKqwc4EXhbNRFF7XY3HWw3d/8+lUSJmQ0Am6j8Tf2EDrSXmf0GsJZKW2Fms4HbgJXALmCzma1y9y1k8G+x40keeMLd/3vtBjN7V/Tjp4C/BoZrHpsOfAg4K9q0DthKJcn/bvQYwF3ArWY23d2rvY604qz6d8BWd6+uK3EqcIKZXQM8A/yRu+8jWZPFdQpwjZm9jUov4UrgnwGz3f2n0T7rqPSqE/3DmiSu6cDl7v7zaPOzVL5nyKa9mjmdytni3ijWHwEfB76S0fvX2g18wd1fi2LZASyO/rstSqz3UumZjjZ/mcRZ9P+/MbM+KonsZfLTblD5O74G+BWda6+LgcupnE0AfAD4mbvvBDCzO4FPmNlzZPBvMdc1eXe/MTpK1zoG2F9zSrMbOD76+bjod6LH91Mpn6TOzN4CXELlS6raDfwJ8B4qR/BbsoglimcO8BRwFfDrwNFUznjG26gmxuOPeIGUuPtw9ZQ56uF8kUrPqxpLR9qLDrdLLXffXv2Hb2YnUCnbPEjlbPZC4F9QWaTqMxmHNh94CDibytn1Z6kk0ly0W3QGNNvd7wEW0qH2cveL3L12EcZmf1uZ/M3loScfVy+VU9ha1aNz/fqdPTWPpe1fA5vc/Z+qG9z97OrPZvZnwAsZxYK7HwB+p+b9b6JyyvgAE9svyzYaFx0U7wWecff10Nn24si/q460S62o9LAZuMrdnUpyrT72bWA1ld50JqKy33jpLyov3EylTFrVyXa7lEo8uPv/pMPtVaPZ31Ymf3O57sk38U/AW8ysuqbyscAvop9/TuUIjpkdBcylptSTsrOo1CaJ3v8tZvZvah7vAbIc3FxsZhfWvf8h4P9QabOqhbzRflnFdiyV5aafBS6KtnW0vchBu9Qys9Oo9Jq/6O7rzezdZvZ7NbtUv88sY1phZh+Qe6gRAAAByklEQVSpi+FFctBuZjaDSs37x9HvHW+vGs3+tjL5mytcko/q64/yxoDqamBL9PMD0e9Ejz+aUD1+UtGA1DJqejnAAeDqaBAGKqP/Wc6EOAj8mZktieK7HLjX3f8X8GqURAA+zRvtl7ro4Hwf8Ffu/nl3r/ZkOt1e/wX4iJn1RwPov0elRJI5M1tEpYR1nrtXOw49wLfMbH40LnUJ2bYPVEp+N5rZLDObC5xP5Qw2D+32HuAf3b06kyYP7VX1d4CZ2Tujv//zgC1Z/VvMQ7nmX5rZx+u2zaAyeNLMZcB6M7sW+N/AH0TbvwysM7PtwP/jjQG9tOPsB15z91erD7j7YTM7F/hOVHv+R944ACWpYVzR7KRLqSTUGVRuMnBT9PingLXRNK//BvxFVnFRWQv714Gjah7/r+5+UUbt1VA00+dLVGZkzAC+7+5/n9X717kSmAXcbFYd6+S7wNeBx6kMXm/IehaSu98fHYSforK2+a3u/nhO2u0dVHrG1VifNbOOtldNLK+a2QXABirf6wPAj6KHU/+3qPXkRURKrHDlGhERCackLyJSYkryIiIlpiQvIlJiSvIiIiWmJC8iUmJK8iIiJaYkLyJSYv8fEOp+UnBbbBEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(mus, rs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}