01_CloudBlog_extended.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Statistical models and PPL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Bayesian Statistics in Climate and Weather\n",
    "\n",
    "In Bayesian statistics, probability is viewed as a degree of belief of an event occurring (in contrast to frequentist statistics where probability is the frequency of an event occurring). The Bayesian view of probability is valuable in climate and weather forecasting. It allows us to quantify how likely we expect events to be, such as the probability of rain on a certain day.\n",
    "\n",
    "In Bayesian statistics, we can incorporate a prior probability: our degree of belief of event occurring before observing any evidence. On an average day, we might expect rain to occur with some probability. We write the prior probability of rain as $P({\\sf event})$. Then we observe some evidence and update our degree of belief about the event occurring. Let's say we look out the window and observe some data: we see there are many dark clouds in the sky.  So we can update our degree of belief of rain occurring this afternoon. This is called conditioning on observations. \n",
    "\n",
    "We write this updated probability as $P({\\sf event}|{\\sf obs})$. This is a conditional probability, or a probability of one event (rain) assuming that another event has occurred (cloud). In particular, we call this the posterior probability of rain, because we have observed some data (clouds) and updated our prior probability. We can use Bayes' theorem to calculate it:\n",
    "\n",
    "$$P({\\sf event}|{\\sf obs}) \\propto P({\\sf obs}|{\\sf event}) \\times P({\\sf event}) $$\n",
    "\n",
    "There are two additional terms here. $P({\\sf obs})$. is the evidence or the overall probability of seeing the cloud, regardless of whether it rains or not. $P({\\sf obs}|{\\sf event})$. is the likelihood, the probability of observing the data if you know the outcome, i.e. the probability of seeing the cloud if it ends up raining later. \n",
    "\n",
    "We are usually concerned with learning the posterior probability $P({\\sf event}|{\\sf obs})$, which is why Bayes' theorem is useful. But it can also be used more generally when we want to learn anything given some data, e.g. we might want to learn parameters of a model. This is called Bayesian inference.\n",
    "\n",
    "This happens to be really useful and has many applications in climate and weather science. It becomes even more powerful when we think about probability distributions, which describe the probabilities associated with all the possible outcomes. One important use of Bayesian statistics is in data assimilation: incorporating observed data from satellites about the current state of the atmosphere into the weather model for prediction of the future state of the atmosphere."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What is a probabilistic programming language (PPL)?\n",
    "\n",
    "The traditional approach to constructing statistical models can be time-consuming and requires expertise. Typically it involves statisticians writing down the problem on paper and carrying out a mathematical procedure to design a bespoke statistical model. This would then be translated into code to construct a model ready for use. Since the model is specifically designed for the task at hand, this method must be repeated for each new problem.\n",
    "\n",
    "Probabilistic programming languages are designed to make this process simpler, quicker and less technical. The idea is to remove the need for the expert statistician to hand-design the statistical model, by embedding this step within the programming language itself. Statistical tasks such as encoding observations (conditioning) and learning latent variables (inference) are automated and some of the complicated mathematical and statistical steps are hidden from the user. This makes statistical modelling more accessible to non-experts and tasks that traditionally require a lot of thought and hand-engineering can be automated."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What does Probabalistic Programming involve? \n",
    "\n",
    "PPLs are just ordinary programming languages, equipped with tools that make statistical modelling easier for the user, reducing the need for hand-designing programs. The main components of a PPL are:\n",
    "\n",
    "1) Sampling: Our model is probabilistic- it requires drawing values at random from probability distributions.\n",
    "\n",
    "2) Conditioning: We have some observed data that can be used to update probabilities\n",
    "\n",
    "3) Inference: We can learn something from the known data and model. This could be a \"latent variable\", some underlying factor that you can't directly observe but might influence the data. This is usually the hard part.\n",
    "\n",
    "Probabilistic programming is a suitable choice when we have a probabalistic model, that relies on sampling from distributions in order to make predictions. PPL makes it easier to do conditioning on observed data and to learn something about the model. Note that probabilistic programming languages are not new. TThey have been around since Simula in 1966. However, with the rise of deep learning, they are growing in popularity, with many new PPLs designed with inference for machine learning in mind (e.g. Pyro, Edward, Infer.Net, webppl).\n",
    "\n",
    "Let's look at a simple example to show how probabilistic programming works. We'll be using Pyro, a probabilistic programming language built on top of pytorch in python. Most of this example comes from the Introduction to Inference in Pyro tutorial."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A simple model for a moving cloud\n",
    "\n",
    "Let's say we have a cloud moving through the atmosphere with some speed in 1 dimension. We observe it after 10 seconds but with some measurement error. We want to use a probabilistic approach to find the distribution of the wind speed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1) Statistical Model\n",
    "The wind speed is something we can't directly observe- its a hidden variable. But we have some prior knowledge for what kind of wind speed to expect. We use this to sample from a prior distribution for the particle speed with a guess wind speed and variance. Then, the particle is moved a distance of speed x time and its final position is observed with some measurement error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "import torch\n",
    "import pyro\n",
    "import pyro.distributions as dist\n",
    "np.random.seed(123)\n",
    "def propagate_cloud(guess, guess_sd, time):\n",
    "    \"\"\" Returns speed and observed position of cloud.\n",
    "    guess = our initial guess for the speed of particle\n",
    "    guess_sd = our guess for the standard deviation of the particle speed\n",
    "    time = length of time before observing \n",
    "    \"\"\"\n",
    "    # sample from prior distribution\n",
    "    speed = pyro.sample(\"speed\", dist.Normal(guess, guess_sd))  \n",
    "    position_final = speed * time\n",
    "    # final position is observed with a bit of noise from our measurement tools\n",
    "    position_observed = pyro.sample(\"position\", dist.Normal(position_final,\n",
    "                                                           torch.tensor([10.])))\n",
    "    return(speed, position_observed)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "guess = torch.tensor([5.])\n",
    "guess_sd = torch.tensor([2.])\n",
    "time = torch.tensor([10.])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(tensor([4.6267]), tensor([44.4646]))"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "propagate_cloud(guess, guess_sd, time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The  `pyro.sample` statements are the first component of PPL, allowing us to sample from distributions. We provide pyro with unique names for each variable, so they can be tracked. We can run this model many times to find the distribution of particle speeds and observed positions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "def plot_hist(hist1, hist2=None, title1=\"\", title2=\"\", true_dist=None):\n",
    "    fig = plt.figure(figsize=(10,4))\n",
    "    ax1 = plt.subplot2grid((1, 2), (0, 0));\n",
    "    ax1.hist(hist1, 25, density=True, label=\"Guide\", alpha=0.7);\n",
    "    ax1.axvline(np.mean(hist1), c=\"k\", ls=\"--\");\n",
    "    plt.title(title1);\n",
    "    if hist2 is not None:\n",
    "        ax2 = plt.subplot2grid((1, 2), (0, 1));\n",
    "        ax2.hist(hist2, 25, density=True, label=\"Guide\", alpha=0.7);\n",
    "        ax2.axvline(np.mean(hist2), c=\"k\", ls=\"--\");\n",
    "        plt.title(title2);\n",
    "    if true_dist is not None:\n",
    "        # Exact posterior distribution of speed: provide tuple for mu and sigma\n",
    "        exact_mu,exact_sigma =  true_dist\n",
    "        x = np.linspace(-5., 15., 100);\n",
    "        f = 1./(exact_sigma*np.sqrt(2*np.pi))*np.exp(-(x-exact_mu)**2/(2*exact_sigma**2))\n",
    "        for ax in fig.axes:\n",
    "            ax.plot(x, f, lw=2, label=\"Exact\", color=\"red\");\n",
    "            ax.axvline(exact_mu, ls=\"--\", color=\"red\");\n",
    "            ax.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl0AAAEICAYAAABlKUHyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3X28VWWd9/HPNwicUlERHeVBdMAmtKIJsfvOHkbSwDFxJrw79gQTDTUjdzk2rzuYSsuytCm1RnvAYCRL0TCnU4MZik1RiaCRCkocUXkcRTF8KDTod/+xrkPb7d6HfeCcdZ19zvf9ep0Xe1/rWtf+rcPhx+9c61prKSIwMzMzs+71ktwBmJmZmfUFLrrMzMzMSuCiy8zMzKwELrrMzMzMSuCiy8zMzKwELrrMzMzMSuCiy3oNST+R9IHccZhZ3yBphKRnJPXroM8zko4pMy7ruVx0WZeQdJKkX0jaLmmbpJ9LOiF3XGZmlSQ9LOn3qRh6VNJ/SNp/b8aKiPURsX9E7Epjv+gXv7R9XVfEbs3PRZftM0kHAj8E/h04BBgKfBp4LmdcZmZ1vD0i9gf+CjgB+ETmeKyPcNFlXeFYgIi4LiJ2RcTvI+LHEXGPpGlp1uvf0yzYA5ImtO8oaZCkuZK2SNok6bOVU/WS3i/pfklPSrpF0lEV205J422XdAWgUo/azJpaRGwCbgaOl3SkpNY0U98m6R/a+0kaL2mFpKfS7NilqX2kpJDUX9JFwBuBK9Is2hWpT0galV4PkvQtSVslPSLpE5JekrZNk7RU0hdTvntI0qSyvyfWvVx0WVf4DbBL0nxJkyQdXLX9RGAdcChwAfA9SYekbfOBncAo4LXAqcAHACSdCfwr8HfAEOBnwHVp26HAjRS/oR4KPAi8obsO0Mx6H0nDgdOAX1Hklo3AkcAU4HMVvyB+GfhyRBwI/AVwQ/VYEfFxihw1M51SnFnjI/8dGAQcA7wZeB/w9xXbTwTWUOS0LwBzJfmXyV7ERZfts4h4CjgJCOAqYGv6jfHw1OUx4PKI+ENEXE+RVP4mbZ8EnBsRz0bEY8BlQEva74PA5yPi/ojYCXwOGJtmu04DVkfEwoj4A3A58D/lHLGZNbn/lPRbYCnw38Acihz2sYjYERErgW8C7039/wCMknRoRDwTEXd09gPTDP47gdkR8XREPAx8qeIzAB6JiKvSGrH5wBHA4S8azJqWiy7rEqkwmhYRw4DjKX5bvDxt3hQvfLL6I2n7UcBLgS2SfpuS4DeAw1K/o4AvV2zbRnEKcWjaf0PF50flezOzDpwZEQdFxFER8U8U+WRbRDxd0ecRilwDMJ1iGcUDkpZLOn0vPvNQYEAat9ZnQMUvjhHxu/Ryrxb5W8/kosu6XEQ8AFxNUXwBDK2aIh8BbKYokp4DDk0J8KCIODAijkv9NgAfrNh2UET8WUT8AtgCDG8fMI0/HDOzztsMHCLpgIq2EcAmgIhYGxFnU/xCeAmwUNLLa4wTNdraPU4xY3ZURdvuz7C+wUWX7TNJfynpo5KGpffDgbOB9in4w4APS3qppLOAVwKLImIL8GPgS5IOlPQSSX8h6c1pv68DsyUdl8YdlPYH+C/gOEl/J6k/8GHgz8s4XjPrXSJiA/AL4POS9pP0aorZre8ASHqPpCER8Ufgt2m3XTWGepRivVatz9hFsRbsIkkHpGUS5wHf7tqjsZ7MRZd1hacpFoAuk/QsRbF1H/DRtH0ZMJriN72LgCkR8UTa9j6KKffVwJPAQop1DETETRS/VS6Q9FQac1La9jhwFnAx8EQa/+fdepRm1pudDYykmPW6CbggIhanbROBVZKeoVhU3xIRO2qM8WVgSrr68Cs1tv9f4FmKC4uWAtcC87r0KKxH0wuX2ph1LUnTgA9ExEm5YzEzM8vJM11mZmZmJXDRZWZmZlYCn140MzMzK4FnuszMzMxK0D93ANUOPfTQGDlyZO4wrARr1qwB4BWveEXmSCy3u+666/GIGJI7jq7gHNacnI9sb3Umf/W4omvkyJGsWLEidxhWgre85S0A/OQnP8kah+Un6ZE992oOzmHNyfnI9lZn8pdPL5qZmZmVoMfNdFnf8YlPfCJ3CGZmgPORlcNFl2Xz1re+NXcIZmaA85GVw6cXLZuVK1eycuXK3GGYmTkfWSk802XZnHvuuYAXrppZfs5HVgbPdJmZmZmVwEWXmZmZWQlcdJmZmZmVwEWXmZmZWQm8kN6y2f8N7wVg+tXL6/aZO+2EssIxsz6gXr6pzEfOO9ZdXHRZNoePfnXuEMzMAOcjK4dPL1o2j669h0fX3pM7DDMz5yMrhYsuy+auG7/KXTd+NXcY1mQkTZS0RlKbpFk1tg+UdH3avkzSyNR+iqS7JN2b/jy5Yp/XpfY2SV+RpNR+iKTFktamPw8u6zitXM5HVgYXXWbWNCT1A64EJgFjgLMljanqNh14MiJGAZcBl6T2x4G3R8SrgKnANRX7fA2YAYxOXxNT+yzgtogYDdyW3puZ7RUXXWbWTMYDbRGxLiKeBxYAk6v6TAbmp9cLgQmSFBG/iojNqX0VsF+aFTsCODAifhkRAXwLOLPGWPMr2s3MOs1Fl5k1k6HAhor3G1NbzT4RsRPYDgyu6vMO4FcR8Vzqv7HOmIdHxJY01hbgsFpBSZohaYWkFVu3bu30QZlZ3+Ciy8yaiWq0RWf6SDqO4pTjBzsxZociYk5EjIuIcUOGDOnMrmbWh/iWEZbNie86L3cI1nw2AsMr3g8DNtfps1FSf2AQsA1A0jDgJuB9EfFgRf9hdcZ8VNIREbElnYZ8rCsPxnoO5yMrg2e6LJvBI45l8Ihjc4dhzWU5MFrS0ZIGAC1Aa1WfVoqF8gBTgCUREZIOAv4LmB0RP2/vnE4bPi3p9emqxfcB368x1tSKdutlnI+sDC66LJtNq+5k06o7c4dhTSSt0ZoJ3ALcD9wQEaskXSjpjNRtLjBYUhtwHn+64nAmMAr4pKSV6at9jdY/At8E2oAHgZtT+8XAKZLWAqek99YLOR9ZGXx60bL59Q/mATD0uPGZI7FmEhGLgEVVbedXvN4BnFVjv88Cn60z5grg+BrtTwAT9jFkawLOR1YGF13W5Tp6lqKZmVlf5dOLZmZmZiVw0WVmZmZWAhddZmZmZiXwmi7L5g1TZ+cOwcwMeGE+amRd6txpJ3RnONZLueiybAYdcVTuEMzMAOcjK4dPL1o261f+jPUrf5Y7DDMz5yMrhWe6LJv7fvQdAEaMfWPmSMysr3M+sjJ4psvMzMysBC66zMzMzErQUNElaaKkNZLaJM2qsf08Sasl3SPpNklHVWybKmlt+ppava+ZmZlZX7DHoktSP+BKYBIwBjhb0piqbr8CxkXEq4GFwBfSvocAFwAnAuOBCyQd3HXhm5mZmTWHRhbSjwfaImIdgKQFwGRgdXuHiLi9ov8dwHvS67cBiyNiW9p3MTARuG7fQ7dm96Z/+HTuEMzMAOcjK0cjRddQYEPF+40UM1f1TAdu7mDfoZ0J0Hqv/QcfnjsEMzPA+cjK0UjRpRptUbOj9B5gHPDmzuwraQYwA2DEiBENhGS9wbpliwE45sRTMkdiZn2d85GVoZGF9BuB4RXvhwGbqztJeivwceCMiHiuM/tGxJyIGBcR44YMGdJo7NbkHrj9Rh64/cbcYZiZOR9ZKRopupYDoyUdLWkA0AK0VnaQ9FrgGxQF12MVm24BTpV0cFpAf2pqMzPbKw1cTT1Q0vVp+zJJI1P7YEm3S3pG0hUV/Q+QtLLi63FJl6dt0yRtrdj2gbKO08x6nz2eXoyInZJmUhRL/YB5EbFK0oXAiohoBf4N2B/4riSA9RFxRkRsk/QZisIN4ML2RfVmZp1VcTX1KRQz6csltUbE6opu04EnI2KUpBbgEuCdwA7gk8Dx6QuAiHgaGFvxGXcB36sY7/qImNlNh2RmfUhDjwGKiEXAoqq28ytev7WDfecB8/Y2QDOzCnu8mjq9/1R6vRC4QpIi4llgqaRR9QaXNBo4DPBD+Mysy/mO9GbWTBq5Inp3n4jYCWwHBjc4/tkUM1uVF/y8I934eaGk4bV2kjRD0gpJK7Zu3drgR5lZX+MHXls2J59zce4QrPk0ckV0w1dc19ACvLfi/Q+A6yLiOUkfAuYDJ79o8Ig5wByAcePGNfpZ1oM4H1kZPNNl2ex3wEHsd8BBucOw5tLIFdG7+0jqDwwC9riWVNJrgP4RcVd7W0Q8UXE19lXA6/Y+dOvJnI+sDJ7psmzWLv0hAKNPOn2fxpl+9fI99pk77YR9+gzrMXZfTQ1sopiZeldVn1ZgKvBLYAqwpOp0YT1nU/W0DElHRMSW9PYM4P59iN16sK7KR2YdcdFl2TjJWWc1eDX1XOAaSW0UM1wt7ftLehg4EBgg6Uzg1IorH/8PcFrVR35Y0hnAzjTWtG47OMvK+cjK4KLLzJpKA1dT7wDOqrPvyA7GPaZG22xg9t7GamZWyWu6zMzMzErgosvMzMysBD69aGZmTa+RC2rMcnPRZdmc+s+X5w7BzAxwPrJyuOiybPoP3C93CGZmgPORlcNruiyb+5cs5P4lC3OHYWbmfGSlcNFl2Tx05608dOetucMwM3M+slK46DIzMzMrgdd0Waf4CiEzs8b4EWVWzTNdZmZmZiVw0WVmZmZWAp9etGxOm/X13CGYmQHOR1YOz3SZmZmZlcBFl2Vz783f5t6bv507DDMz5yMrhYsuy2bDr5ey4ddLc4dhZuZ8ZKVw0WVmZmZWAhddZtZUJE2UtEZSm6RZNbYPlHR92r5M0sjUPljS7ZKekXRF1T4/SWOuTF+HdTSWmdnecNFlZk1DUj/gSmASMAY4W9KYqm7TgScjYhRwGXBJat8BfBL4lzrDvzsixqavx/YwlplZp7nosmz6DRhIvwEDc4dhzWU80BYR6yLieWABMLmqz2Rgfnq9EJggSRHxbEQspSi+GlVzrL0P33oq5yMrg+/TZdm87bwv5w7Bms9QYEPF+43AifX6RMROSduBwcDjexj7PyTtAm4EPhsR0ehYkmYAMwBGjBixF4dluTkfWRk802VmzaTWLFPsRZ9q746IVwFvTF/v7cxYETEnIsZFxLghQ4bs4aPMrK9y0WXZrGydy8rWubnDsOayERhe8X4YsLleH0n9gUHAto4GjYhN6c+ngWspTmPu1VjWnJyPrAwuuiybzauXs3n18txhWHNZDoyWdLSkAUAL0FrVpxWYml5PAZakU4U1Seov6dD0+qXA6cB9ezOWNS/nIyuD13SZWdNI66pmArcA/YB5EbFK0oXAiohoBeYC10hqo5iVamnfX9LDwIHAAElnAqcCjwC3pIKrH3ArcFXape5YZmad5aLLdpt+tX/Ls54vIhYBi6razq94vQM4q86+I+sM+7o6/euOZX2b86XtDZ9eNDMzMyuBZ7osm4H7D8odgpkZ4Hxk5XDRZdlMmOmbe5tZz+B8ZGVw0WU9Wletm2h0nLnTTuiSzzMzM6vmNV2WzYrvXsmK716ZOwwzM+cjK4Vnuiybxx68N3cIZmaA85GVo6GZLkkTJa2R1CZpVo3tb5J0t6SdkqZUbdslaWX6qr6JoZmZmVmfsMeZLkn9gCuBUygeibFcUmtErK7oth6YBvxLjSF+HxFjuyBWMzMzs6bVyOnF8UBbRKwDkLQAmAzsLroi4uG07Y/dEKOZmZlZ02vk9OJQYEPF+42prVH7SVoh6Y702I0XkTQj9VmxdevWTgxtzezlBx/Gyw8+LHcYZmbOR1aKRma6VKOtMw98HRERmyUdAyyRdG9EPPiCwSLmAHMAxo0b54fJ9hFv/uCFuUMwMwOcj6wcjcx0bQSGV7wfBmxu9AMiYnP6cx3wE+C1nYjPzMzMrFdopOhaDoyWdLSkAUAL0NBViJIOljQwvT4UeAMVa8Gsb7vj2ku549pLc4dhZuZ8ZKXY4+nFiNgpaSZwC9APmBcRqyRdCKyIiFZJJwA3AQcDb5f06Yg4Dngl8I20wP4lwMVVVz1aH7Zt/W9yh2BmTaCrnkzREecjK0NDN0eNiEXAoqq28yteL6c47Vi93y+AV+1jjGZmZmZNz48BMjMzMyuBiy4zayoNPCFjoKTr0/Zlkkam9sGSbpf0jKQrKvq/TNJ/SXpA0ipJF1dsmyZpa8VTNT5QxjGaWe/kZy9aNgf++YjcIViTafAJGdOBJyNilKQW4BLgncAO4JPA8emr0hcj4vZ0sdBtkiZFxM1p2/URMbMbD6tPK2O9ViOcj6wMLrosm5Om/WvuEKz57PEJGen9p9LrhcAVkhQRzwJLJY2qHDAifgfcnl4/L+luaqxRtd7N+cjK4NOLZtZMGnlCxu4+EbET2A4MbmRwSQcBbwduq2h+h6R7JC2UNLzOfn6qhpntkYsuy2bp1Z9j6dWfyx2GNZdGnpCxV0/RkNQfuA74SvtMGvADYGREvBq4FZhfa9+ImBMR4yJi3JAhQ/b0UdYDOR9ZGXx60bJ56n/W5w7Bmk8jT8ho77MxFVKDgG0NjD0HWBsRl7c3RMQTFduvolgfZr2Q85GVwTNdZtZMGnlCRiswNb2eAiyJiA5nuiR9lqI4O7eq/YiKt2cA9+9D7GbWx3mmy8yaRiNPyADmAtdIaqOY4Wpp31/Sw8CBwABJZwKnAk8BHwceAO6WBHBFRHwT+LCkM4CdaaxppRyomfVKLrrMrKk08ISMHcBZdfYdWWfYWuvAiIjZwOy9CtTMrIqLLsvmkBHH5g7BzAxwPrJyuOiybF7/rvNyh2BmBjgfWTm8kN7MzMysBC66LJv//sb5/Pc3zt9zRzOzbuZ8ZGXw6UXL5tknH8sdgpkZ4Hxk5fBMl5mZmVkJXHSZmZmZlcBFl5mZmVkJvKbLsjnsL16VOwQzM8D5yMrhosuyGXfWOblDMDMDnI+sHD69aGZmZlYCF12WzW1XfIzbrvhY7jDMzJyPrBQ+vWjZPPfM9twhmJkB+fLR9KuXN9Rv7rQTujkSK4OLrj6i0X/YZmZm1j18etHMzMysBC66zMzMzErgosuyOXLMCRw5xusUrHMkTZS0RlKbpFk1tg+UdH3avkzSyNQ+WNLtkp6RdEXVPq+TdG/a5yuSlNoPkbRY0tr058FlHKOVz/nIyuCiy7IZe8Z0xp4xPXcY1kQk9QOuBCYBY4CzJY2p6jYdeDIiRgGXAZek9h3AJ4F/qTH014AZwOj0NTG1zwJui4jRwG3pvfVCzkdWBhddZtZMxgNtEbEuIp4HFgCTq/pMBuan1wuBCZIUEc9GxFKK4ms3SUcAB0bELyMigG8BZ9YYa35Fu5lZp7nosmxuufQj3HLpR3KHYc1lKLCh4v3G1FazT0TsBLYDg/cw5sY6Yx4eEVvSWFuAw2oNIGmGpBWSVmzdurXBQ7GexPnIyuCiy7LZ9fxz7Hr+udxhWHNRjbbYiz770v/FnSPmRMS4iBg3ZMiQzuxqPYTzkZXBRZeZNZONwPCK98OAzfX6SOoPDAK27WHMYXXGfDSdfmw/DfnYXkduZn2eiy4zaybLgdGSjpY0AGgBWqv6tAJT0+spwJK0VqumdNrwaUmvT1ctvg/4fo2xpla0m5l1mu9Ib2ZNIyJ2SpoJ3AL0A+ZFxCpJFwIrIqIVmAtcI6mNYoarpX1/SQ8DBwIDJJ0JnBoRq4F/BK4G/gy4OX0BXAzcIGk6sB44q/uP0sx6Kxddls3w15yUO4QXaeRxSX4GWl4RsQhYVNV2fsXrHdQpjiJiZJ32FcDxNdqfACbsQ7h9WjM9fqwn5iPrfVx0WTavmvSe3CGYmQHOR1aOhtZ0NXAH6DdJulvSTklTqrZNTXdzXitpavW+ZmZmZn3BHouuBu8AvR6YBlxbte8hwAXAiRQ3NbzAj9Gwdosu/hCLLv5Q7jDMzJyPrBSNzHTt8Q7QEfFwRNwD/LFq37cBiyNiW0Q8CSzmT4/XMDMzM+szGim6GrkD9D7t67s5m5mZWW/XSNG1L3drbmhf383ZzMzMertGiq5G7gDdHfuamZmZ9RqN3DJi9x2ggU0UNxp8V4Pj3wJ8rmLx/KnA7E5Hab3S0ePfmjsEMzPA+cjKsceiq5E7QEs6AbgJOBh4u6RPR8RxEbFN0mcoCjeACyOio2egWR/yypOn7LmTmVkJnI+sDA3dHLWBO0Av54UPjK3sNw+Ytw8x2h40012fK+18bgcA/QfulzkSM+vrnI+sDH7gtWXz48vO5ceXnZs7DDMz5yMrhYsuMzMzsxL42YtmZmY9XCPLSOZOO6GESGxfeKbLzMzMrAQuuszMzMxK4NOLls3ok07PHYKZGeB8ZOXwTJdlM/qk053orNMkTZS0RlKbpFk1tg+UdH3avkzSyIpts1P7GklvS22vkLSy4uspSeembZ+StKli22llHaeVy/nIyuCZLstmx9O/BWC/Aw7KHIk1C0n9gCuBUygeM7ZcUmtErK7oNh14MiJGSWoBLgHeKWkMxRM1jgOOBG6VdGxErAHGVoy/ieJmz+0ui4gvdvexNZNmvTdgR5yPrAye6bJsllw5iyVXvmiiwqwj44G2iFgXEc8DC4DJVX0mA/PT64XABElK7Qsi4rmIeAhoS+NVmgA8GBGPdNsRWI/kfGRlcNFlZs1kKLCh4v3G1FazT0TsBLYDgxvctwW4rqptpqR7JM2reI6smVmn+fRiD9Ybp/DN9pFqtEWDfTrcV9IA4AxgdsX2rwGfSf0+A3wJeP+LgpJmADMARowYUT96M+vTPNNlZs1kIzC84v0wYHO9PpL6A4OAbQ3sOwm4OyIebW+IiEcjYldE/BG4ihefjmzvNycixkXEuCFDhuzVgZlZ7+eZLrNO8p2hs1oOjJZ0NMWC9xbgXVV9WoGpwC+BKcCSiAhJrcC1ki6lWEg/GrizYr+zqTq1KOmIiNiS3v4tcF8XH4+Z9SEuuiybv/zrd+QOwZpMROyUNBO4BegHzIuIVZIuBFZERCswF7hGUhvFDFdL2neVpBuA1cBO4JyI2AUg6WUUV0R+sOojvyBpLMXpxYdrbLdewvnIyuCiy7I55sRTcodgTSgiFgGLqtrOr3i9Azirzr4XARfVaP8dxWL76vb37mu81hycj6wMLrosm2eeKJbO7D/48MyRmFlf1xvykZc+9HxeSG/Z/PSqC/jpVRfkDsPMzPnISuGiy8zMzKwELrrMzMzMSuCiy8zMzKwELrrMzMzMSuCrFy2b4ye+O3cIZmaA85GVw0WXZTNi7Btzh2BmBjgfWTl8etGy2b7lEbZveSR3GGZmzkdWChddls3P53+en8//fO4wzMycj6wULrrMzMzMSuCiy8zMzKwELrrMzMzMSuCiy8zMzKwEvmWEZfOat78/dwhmZoDzkZXDRZdlM/S48blDMDMDnI+sHD69aNk8sf43PLH+N7nDsCYjaaKkNZLaJM2qsX2gpOvT9mWSRlZsm53a10h6W0X7w5LulbRS0oqK9kMkLZa0Nv15cHcfn+XhfGRlcNFl2Sy79lKWXXtp7jCsiUjqB1wJTALGAGdLGlPVbTrwZESMAi4DLkn7jgFagOOAicBX03jt/joixkbEuIq2WcBtETEauC29t17I+cjK4KLLzJrJeKAtItZFxPPAAmByVZ/JwPz0eiEwQZJS+4KIeC4iHgLa0ngdqRxrPnBmFxyDmfVRXtNlZs1kKLCh4v1G4MR6fSJip6TtwODUfkfVvkPT6wB+LCmAb0TEnNR+eERsSWNtkXRYVx5MTzT96uW5QzDrtVx0mVkzUY22aLBPR/u+ISI2p6JqsaQHIuKnDQclzQBmAIwYMaLR3cysj2no9OLeLlyVNFLS79Pi1JWSvt614ZtZH7MRGF7xfhiwuV4fSf2BQcC2jvaNiPY/HwNu4k+nHR+VdEQa6wjgsVpBRcSciBgXEeOGDBmy1wdnZr3bHme6KhaunkKRtJZLao2I1RXddi9cldRCsXD1nWnbgxExtovjtl7gde/4p9whWPNZDoyWdDSwiWJh/Luq+rQCU4FfAlOAJRERklqBayVdChwJjAbulPRy4CUR8XR6fSpwYdVYF6c/v9+tR2fZOB9ZGRo5vbh74SqApPaFq5VF12TgU+n1QuCKtHDVrK7DR786dwjWZNIarZnALUA/YF5ErJJ0IbAiIlqBucA1ktooZrha0r6rJN1Akbt2AudExC5JhwM3pZTVH7g2In6UPvJi4AZJ04H1wFmlHWw38Hqt+pyPrAyNFF37snAV4GhJvwKeAj4RET+r/oC+uB7CyQ8eXXsP4GRnnRMRi4BFVW3nV7zeQZ3iKCIuAi6qalsHvKZO/yeACfsYsjUB5yMrQyNruvZl4eoWYEREvBY4j2Jq/8AXdfR6iD7prhu/yl03fjV3GGZmzkdWikaKrr1euJruh/MEQETcBTwIHLuvQZuZmZk1m0aKrt0LVyUNoFgf0VrVp32xKbxw4eqQ9js+SzqGYuHquq4J3czMzKx57HFN174sXAXeBFwoaSewC/hQRGzrjgMxMzMz68kaujnq3i5cjYgbgRv3MUYzMzOzpuc70ls2J77rvNwhmJkBzkdWDhddls3gEb6mwsx6BucjK0NDjwEy6w6bVt3JplV35g7DzMz5yErhmS7L5tc/mAfA0OPG76Fn79TIDXLnTjuhhEjMrK/nIyuHiy4zMzPbrdEnpviXws5z0dUN/IgfMzMzq+aiy8zMrI/wpEBeXkhvZmZmVgLPdFk2b5g6O3cIZmaA85GVw0WXZTPoiKNyh9BtPIVv1lx6cz6ynsOnFy2b9St/xvqVP8sdhpmZ85GVwjNdls19P/oOACPGvjFzJGbW1zkfWRk802VmTUXSRElrJLVJmlVj+0BJ16ftyySNrNg2O7WvkfS21DZc0u2S7pe0StJHKvp/StImSSvT12llHKOZ9U6e6TKzpiGpH3AlcAqwEVguqTUiVld0mw48GRGjJLUAlwDvlDQGaAGOA44EbpV0LLAT+GhE3C3pAOAuSYsrxrwsIr5YzhGaWW/mmS4zaybjgbaIWBcRzwMLgMlVfSYD89PrhcAESUrtCyLiuYh4CGgDxkfEloi4GyAingbuB4aWcCxm1se46DKzZjIU2FDxfiMvLpB294mIncB2YHAj+6bktuRQAAAJjklEQVRTka8FllU0z5R0j6R5kg6uFZSkGZJWSFqxdevWzh6TmfURPr1o2bzpHz6dOwRrPqrRFg326XBfSfsDNwLnRsRTqflrwGdSv88AXwLe/6JBIuYAcwDGjRtXHY81AecjK4OLLstm/8GH5w7Bms9GYHjF+2HA5jp9NkrqDwwCtnW0r6SXUhRc34mI77V3iIhH219Lugr4YZcdifUozkdWBp9etGzWLVvMumWLc4dhzWU5MFrS0ZIGUCyMb63q0wpMTa+nAEsiIlJ7S7q68WhgNHBnWu81F7g/Ii6tHEjSERVv/xa4r8uPyHoE5yMrg2e6LJsHbr8RgGNOPCVzJNYsImKnpJnALUA/YF5ErJJ0IbAiIlopCqhrJLVRzHC1pH1XSboBWE1xxeI5EbFL0knAe4F7Ja1MH/WvEbEI+IKksRSnFx8GPljawVqpnI+sDC66zKyppGJoUVXb+RWvdwBn1dn3IuCiqral1F7vRUS8d1/jLYsfPWXW8/n0opmZmVkJPNNlZmZmndbI7OrcaSeUEEnz8EyXmZmZWQk802XZnHzOxblDMDMDnI+sHC66LJv9Djgodwh9ik8FmNXnfGRl8OlFy2bt0h+ydqnvNWlm+TkfWRk802XZtCe40SednjmSnqvR2wB4hqp38+0gup/zkZXBM11mZmZmJfBMVyf5N04zMzPbG57pMjMzMyuBZ7oSz2BZM/PPr5lZz+eiy7I59Z8vzx2CmRngfGTlcNFl2fQfuF/uEMzMAOcjK4eLLsvm/iULAXjlyVMyR2JmfZ3zUffwTZlfyAvpLZuH7ryVh+68NXcYZmbOR1aKpp/p8gJis67j30rL5fxl1piuzE0581xDM12SJkpaI6lN0qwa2wdKuj5tXyZpZMW22al9jaS3dV3oZtYXdUc+qjempKPTGGvTmAO6+/jMrPfa40yXpH7AlcApwEZguaTWiFhd0W068GREjJLUAlwCvFPSGKAFOA44ErhV0rERsaurD8TMypHz0UTdkY/SPvXGvAS4LCIWSPp6GvtrXX5gZtYnNDLTNR5oi4h1EfE8sACYXNVnMjA/vV4ITJCk1L4gIp6LiIeAtjSemdne6I58VHPMtM/JaQzSmGd247GZWS/XyJquocCGivcbgRPr9YmInZK2A4NT+x1V+w6t/gBJM4AZ6e0zktY0FH3XOxR4PNNn19PrY5r3911Sh/f671MXKS2meX/fcNdDgaMa7Ntd+ajWmIOB30bEzhr9XyBjDusJP2O9KoZ9zEe96ntRVgydyBWdGWefvg+djKnR/NVQ0aUabdFgn0b2JSLmAHMaiKVbSVoREeNyx1HJMTXGMTWmB8c0stHuNdr2NR/VmvFvOH9BvhzWE/4+HUPPisMx9JwYamnk9OJGYHjF+2HA5np9JPUHBgHbGtzXzKxR3ZGP6rU/DhyUxqj3WWZmDWuk6FoOjE5X8QygWIjaWtWnFZiaXk8BlkREpPaWdDXR0cBo4M6uCd3M+qDuyEc1x0z73J7GII35/W48NjPr5fZ4ejGtiZgJ3AL0A+ZFxCpJFwIrIqIVmAtcI6mN4jfKlrTvKkk3AKuBncA5PfzKxeynOGtwTI1xTI1p6pi6Kx/VGjN95MeABZI+C/wqjd2T9IS/T8fwJz0hDsdQ6AkxvIiKX+bMzMzMrDv5MUBmZmZmJXDRZWZmZlYCF13Jnh4tkiGe4ZJul3S/pFWSPpI7JijuCC7pV5J+mDuWdpIOkrRQ0gPp+/W/Msfzz+nv7D5J10naL1Mc8yQ9Jum+irZDJC1Oj7VZLOngHhDTv6W/u3sk3STpoDJjalY5cla9vJTj56o6F5X9yKZaeafs70OtXNPd34fO5BUVvpJ+Ru+R9FfdHEfdXKIe8khCF1284NEik4AxwNkqHhmS007goxHxSuD1wDk9ICaAjwD35w6iypeBH0XEXwKvIWN8koYCHwbGRcTxFAuzWzKFczUwsaptFnBbRIwGbkvvc8e0GDg+Il4N/AaYXXJMTSdjzqqXl3L8XFXnovZHNo0GnqR4ZFN3qpV3Svs+dJBruvv7cDWN55VJFFcJj6a4eXBXPkKrVhw1c4le+AiwicBX07+h0rnoKjTyaJFSRcSWiLg7vX6a4h90zbthl0XSMOBvgG/mjKOSpAOBN5GuKouI5yPit3mjoj/wZyru7/QyMt3bKSJ+SnH1XqXKR+SU/libWjFFxI8r7vp+B8X9sKxjWXJWB3mp1J+r6lwklfvIpg7yTtn/vqpzzRa6+fvQybwyGfhWFO6guO/dEd0VRwe5pMc8ktBFV6HWo0WyFjiVJI0EXgssyxsJlwP/D/hj5jgqHQNsBf4jnWr4pqSX5womIjYBXwTWUyTA7RHx41zx1HB4RGyB4j9Q4LDM8VR7P3Bz7iCaQPacVZWXyv65qs5FDT+yqYvUyzulfR9q5RrgLsr9PrSrd9w5f04rc0n2fy/tXHQVGn7cR9kk7Q/cCJwbEU9ljON04LGIuCtXDHX0B/4K+FpEvBZ4lvJPme2W1jJMBo4GjgReLuk9ueJpJpI+TnH66ju5Y2kCWXNWzrxUJxeV/f3Inndq5RqK03nVcv5fluXntEYu6TH/x7voKvTIxxVJeilFYvtORHwvczhvAM6Q9DDFqYyTJX07b0hA8Xe3MSLaZwEXUiTDXN4KPBQRWyPiD8D3gP+dMZ5qj7ZP76c/H8scDwCSpgKnA+8O3zywEdlyVp28VObP1YtyEcXMV5mPbKqXd8r8PtTLNTkeXVXvuEv/Oa2TS3rM//EuugqNPFqkVGmNwlzg/oi4NGcsABExOyKGpQcTt1A8WiX7DE5E/A+wQdIrUtMEijuO57IeeL2kl6W/wwn0rAsPKh+R0yMeayNpIsWd38+IiN/ljqdJZMlZHeSl0n6u6uSid1PiI5s6yDtl/vuqlWtWk+fRVfWOuxV4X7qK8fUUyy22dFcQHeSSnvNIwojwV1EMn0ZxtcODwMd7QDwnUUx/3gOsTF+n5Y4rxfYW4Ie546iIZyywIn2v/hM4OHM8nwYeAO4DrgEGZorjOoq1Hn+g+E1vOsXal9uAtenPQ3pATG0U6y3af86/nvtnqhm+cuSsenkp189VZS6iWGd1Z/p5+m53/7urlXfK/j7UyjXd/X3oTF6hOK13ZfoZvZfiSsvujKNuLgE+nuJYA0wq4+ez1pcfA2RmZmZWAp9eNDMzMyuBiy4zMzOzErjoMjMzMyuBiy4zMzOzErjoMjMzMyuBiy4zMzOzErjoMjMzMyvB/wf4ZJZqwJdszQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "uncond = [propagate_cloud(guess=torch.tensor([5.]),\n",
    "                          guess_sd=torch.tensor([2.]),\n",
    "                          time=torch.tensor([10.])) for i in range(1000)]\n",
    "speed_uncond, pos_uncond = zip(*uncond)\n",
    "plot_hist(speed_uncond, pos_uncond, \"Speed\", \"Position\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2) Conditioning\n",
    "\n",
    "What if we observe the particle at a certain point, x=30 m? We want to use this information in our model to learn more about the latent variable, speed. This is when we condition on our observations.\n",
    "In pyro, this is done with the pyro.condition statement. This statement returns a new function almost identical to the original, but overwrites the   `pyro.sample(\"position\",...)` by fixing it at the observed value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "position_observed = torch.tensor([30.])\n",
    "conditioned_propagate = pyro.condition(propagate_cloud,\n",
    "                                       data={\"position\": position_observed})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(tensor([4.5267]), tensor([30.]))"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "conditioned_propagate(guess, guess_sd, time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl0AAAEICAYAAABlKUHyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XuYXXV97/H3RxCwcpdAEQgXjVbwgjWA54BK5SJYC7bCI3hUolH0VE6l2ucI6iOVnipaRTgVq2hSsIpWQdvUYjFyEdEDJtTILSAXMQQiICAXrWjwe/7YK7odJ2SFTNbaM/N+Pc88s/e6fmaS9Z3v/q2190pVIUmSpPXrcX0HkCRJmg5suiRJkjpg0yVJktQBmy5JkqQO2HRJkiR1wKZLkiSpAzZdmjKSXJLkDX3nkDQ9JJmZ5KEkGzzKMg8l2a3LXBpdNl2aEEn2S/LtJPcnuTfJt5Ls1XcuSRqW5NYk/9U0Q3cm+cckmz6WbVXVsqratKoeabb9Oy/8mvm3TER2TX42XVpnSTYHvgL8PbA1sAPwXuDhPnNJ0mr8SVVtCvwhsBfw7p7zaJqw6dJEeBpAVX2uqh6pqv+qqq9V1VVJ5jSjXn/fjIJdn+SAVSsm2SLJvCQrktye5P8MD9UneX2SpUnuS3JBkp2H5h3UbO/+JB8F0ulPLWlSq6rbga8Cz0zy5CQLmpH6m5K8cdVySfZOsjjJA83o2KnN9F2SVJINk/wt8ALgo80o2kebZSrJU5vHWyT5dJK7k/wwybuTPK6ZNyfJZUk+1NS7HyQ5tOvfidYvmy5NhO8DjyQ5O8mhSbYaM38f4BZgG+Ak4EtJtm7mnQ2sBJ4KPBc4GHgDQJKXA+8E/gyYAXwT+FwzbxvgPAavULcBbgb2XV8/oKSpJ8lOwEuB7zKoLcuBJwNHAO8beoF4OnB6VW0OPAX4wthtVdW7GNSo45pTiseNs8u/B7YAdgNeBLwWeN3Q/H2AGxjUtA8C85L4YnIKsenSOquqB4D9gAI+CdzdvGLcrlnkLuC0qvplVf0zg6Lyx838Q4Hjq+qnVXUX8BHgqGa9NwHvr6qlVbUSeB+wZzPa9VLguqo6t6p+CZwG/Kibn1jSJPcvSX4CXAZ8AziTQQ17R1X9vKqWAJ8CXtMs/0vgqUm2qaqHqurytd1hM4L/SuDEqnqwqm4FPjy0D4AfVtUnm2vEzga2B7b7nY1p0rLp0oRoGqM5VbUj8EwGrxZPa2bfXr99Z/UfNvN3Bh4PrEjyk6YIfgLYtlluZ+D0oXn3MjiFuEOz/m1D+6/h55L0KF5eVVtW1c5V9ecM6sm9VfXg0DI/ZFBrAOYyuIzi+iSLkrzsMexzG2CjZrvj7QOGXjhW1c+ah4/pIn+NJpsuTbiquh44i0HzBbDDmCHymcAdDJqkh4FtmgK4ZVVtXlV7NMvdBrxpaN6WVfWEqvo2sALYadUGm+3vhCStvTuArZNsNjRtJnA7QFXdWFVHM3hB+AHg3CRPHGc7Nc60VX7MYMRs56Fpv96HpgebLq2zJH+Q5O1Jdmye7wQcDawagt8W+Iskj09yJPAM4PyqWgF8Dfhwks2TPC7JU5K8qFnv48CJSfZotrtFsz7AvwN7JPmzJBsCfwH8fhc/r6SppapuA74NvD/JJkmezWB067MASV6dZEZV/Qr4SbPaI+Ns6k4G12uNt49HGFwL9rdJNmsuk3gb8JmJ/Wk0ymy6NBEeZHAB6BVJfsqg2boGeHsz/wpgFoNXen8LHFFV9zTzXstgyP064D7gXAbXMVBVX2bwqvLzSR5otnloM+/HwJHAKcA9zfa/tV5/SklT2dHALgxGvb4MnFRVC5t5hwDXJnmIwUX1R1XVz8fZxunAEc27D//vOPP/F/BTBm8sugw4B5g/oT+FRlp++1IbaWIlmQO8oar26zuLJEl9cqRLkiSpAzZdkiRJHfD0oiRJUgcc6ZIkSerAhn0HGGubbbapXXbZpe8Y6tkNN9wAwNOf/vSek6gLV1555Y+rakbfOSaCNWx6smZNX2tTv0au6dpll11YvHhx3zHUs/333x+ASy65pNcc6kaSH655qcnBGjY9WbOmr7WpX55elCRJ6sDIjXRJAO9+97v7jiBJrVmz1IZNl0bSgQce2HcESWrNmqU2PL2okbRkyRKWLFnSdwxJasWapTYc6dJIOv744wEvSpU0OViz1IYjXZIkSR2w6ZI07SXZKcnFSZYmuTbJW5vpf53k9iRLmq+X9p1V0uTl6UVJgpXA26vqP5NsBlyZZGEz7yNV9aEes0maImy6JE17VbUCWNE8fjDJUmCHflNJmmpsutSbuWctWu28Tfd9TYdJpN9IsgvwXOAKYF/guCSvBRYzGA27b5x1jgWOBZg5c2ZnWdWfsfVrVc0anj5vzl6dZtLo85oujaTtZj2b7WY9u+8YmmaSbAqcBxxfVQ8A/wA8BdiTwUjYh8dbr6rOrKrZVTV7xowpcQtJrSVrltpwpEsj6c4br2oe+UpR3UjyeAYN12er6ksAVXXn0PxPAl/pKZ5G3KqaZeOlR+NIl0bSled9jCvP+1jfMTRNJAkwD1haVacOTd9+aLE/Ba7pOpsmB2uW2nCkS5IG1269Brg6yaqPFX8ncHSSPYECbgXe1E88SVOBTZekaa+qLgMyzqzzu84iaery9KIkSVIHbLokSZI64OlFjaR9XvW2viNIUmvWLLVh06WR9KSZT+s7giS1Zs1SGzZdGkm3X/ud5pGf0yVp9K2qWTvssXfPSTTKbLo0kr73b/MHD/7uLf0GkaQWVtUsmy49Gi+klyRJ6oBNlyRJUgdsuiRJkjpg0yVJktQBL6TXSNr3mBP7jiBJrVmz1IZNl0bSFtvv3HcESWrNmqU2bLo0kpYt+WbzyM/pkjT6VtWsmXu+oOckGmU2XRpJ1/zHZ5tH3lpD0uhbVbNsuvRovJBekiSpAzZdkiRJHbDpkiRJ6oBNlyRJUge8kF4j6YVvfG/fESSpNWuW2rDp0kja9Enb9R1BklqzZqkNmy6NpFuuWDh4MMfP6ZI0+lbVrN32OajnJBplra7pSnJIkhuS3JTkhHHmvy3JdUmuSnJhkp2H5h2T5Mbm65iJDK+p6/qLz+P6i8/rO4YktWLNUhtrbLqSbACcARwK7A4cnWT3MYt9F5hdVc8GzgU+2Ky7NXASsA+wN3BSkq0mLr4kSdLk0Gaka2/gpqq6pap+AXweOHx4gaq6uKp+1jy9HNixefwSYGFV3VtV9wELgUMmJrokSdLk0abp2gG4bej58mba6swFvro26yY5NsniJIvvvvvuFpEkSZImlzZNV8aZVuMumLwamA383dqsW1VnVtXsqpo9Y8aMFpEkSZImlzbvXlwO7DT0fEfgjrELJTkQeBfwoqp6eGjd/cese8ljCarp5cVvOaXvCJLUmjVLbbQZ6VoEzEqya5KNgKOABcMLJHku8AngsKq6a2jWBcDBSbZqLqA/uJkmPapNNtuSTTbbsu8YktSKNUttrHGkq6pWJjmOQbO0ATC/qq5NcjKwuKoWMDiduCnwxSQAy6rqsKq6N8nfMGjcAE6uqnvXy0+iKeXGy74yeODndEmaBFbVrFn7vaznJBplrT4ctarOB84fM+09Q48PfJR15wPzH2tATU+/brrw1hqSRp9Nl9rwhteSJEkdsOmSJEnqgE2XpGkvyU5JLk6yNMm1Sd7aTN86ycLmNmYLvaOGpHVh0yVJsBJ4e1U9A3g+8JbmdmcnABdW1Szgwua5JD0mrS6kl7p28F+e1ncETSNVtQJY0Tx+MMlSBnfPOJzffNbg2Qw+Z/AdPUTUiLNmqQ2bLo2kDTfepO8ImqaS7AI8F7gC2K5pyKiqFUm2Xc06xwLHAsycObOboBop1iy14elFjaSlF53L0ovO7TuGppkkmwLnAcdX1QNt1/NWZrJmqQ1HujSSfvCdr/cdQdNMksczaLg+W1VfaibfmWT7ZpRre+Cu1W9B09mqmvWMFx/RcxKNMke6JE17GdxKYx6wtKpOHZq1ADimeXwM8K9dZ5M0dTjSJUmwL/Aa4OokS5pp7wROAb6QZC6wDDiyp3ySpgCbLo20uWctWuMy87w/o9ZRVV0GZDWzD+gyi6Spy9OLkiRJHXCkSyPppSd8vO8IktSaNUttONIlSZLUAZsujaSrv/oZrv7qZ/qOIUmtWLPUhqcXNZJu+95lADzr0FevcVkvtpfUt7WpWZq+HOmSJEnqgE2XJElSB2y6JEmSOuA1XRpJG2y0cd8RJKk1a5basOnSSHrJ207vO4IktWbNUhueXpQkSeqATZdG0pIF81iyYF7fMSSpFWuW2rDp0ki647pF3HHdmj9/S5JGgTVLbdh0SZIkdcCmS5IkqQM2XZIkSR3wIyM0kjbedIu+I0hSa9YstWHTpZF0wHEf6DuCJLVmzVIbnl6UJEnqgE2XRtLiL57B4i+e0XcMSWrFmqU2PL2okXTXzVf3HUGSWrNmqQ1HuiRJkjpg0yVJktQBmy5JkqQOtLqmK8khwOnABsCnquqUMfNfCJwGPBs4qqrOHZr3CLDqZPeyqjpsIoJranviVttO6PbmntXunmjz5uw1ofuVND1MdM3S1LTGpivJBsAZwEHAcmBRkgVVdd3QYsuAOcBfjbOJ/6qqPScgq6aRF73p5L4jSFJr1iy10Waka2/gpqq6BSDJ54HDgV83XVV1azPvV+shoyRJ0qTX5pquHYDbhp4vb6a1tUmSxUkuT/Ly8RZIcmyzzOK77757LTatqeryc07l8nNO7TuGJLVizVIbbUa6Ms60Wot9zKyqO5LsBlyU5Oqquvm3NlZ1JnAmwOzZs9dm25qi7l32/b4jSFJr1iy10Wakazmw09DzHYE72u6gqu5ovt8CXAI8dy3ySZIkTQltRroWAbOS7ArcDhwFvKrNxpNsBfysqh5Osg2wL/DBxxpWk0PbdwpKkjSdrHGkq6pWAscBFwBLgS9U1bVJTk5yGECSvZIsB44EPpHk2mb1ZwCLk3wPuBg4Zcy7HiVJkqaFVp/TVVXnA+ePmfaeoceLGJx2HLvet4FnrWNGTUOb//7MviNoGkkyH3gZcFdVPbOZ9tfAG4FV7+55Z1MLpd9hzVIb3vBaa6WrU4f7zXlnJ/uRGmcBHwU+PWb6R6rqQ93H0WRjzVIb3gZI0rRXVZcC9/adQ9LUZtOlkXTZWe/jsrPe13cM6bgkVyWZ37wxSBqXNUtt2HRpJD3wo2U88KNlfcfQ9PYPwFOAPYEVwIdXt6Af8Cxrltqw6ZKkcVTVnVX1SFX9Cvgkg1uirW7ZM6tqdlXNnjFjRnchJU0qNl2SNI4k2w89/VPgmr6ySJoafPeipGkvyeeA/YFtms8cPAnYP8meDG57divwpt4CSpoSbLo0krae+bS+I2gaqaqjx5k8r/MgmrSsWWrDpksj6fmvelvfESSpNWuW2vCaLkmSpA7YdGkkfeMT7+Ebn3jPmheUpBFgzVIbnl7USPrpfXf1HUGSWrNmqQ1HuiRJkjpg0yVJktQBmy5JkqQOeE2XRtK2T3lW3xEkqTVrltqw6dKvzT1rUd8Rfm32kW/pO4IktWbNUhueXpQkSeqATZdG0oUffQcXfvQdfceQpFasWWrD04saSQ8/dH/fESSpNWuW2nCkS5IkqQM2XZIkSR3w9KI0pM07OOfN2auDJJKkqcamSyPpybvb2EiaPKxZasOmSyNpz8Pm9h1BklqzZqkNr+mSJEnqgE2XRtIFp76VC059a98xJKkVa5ba8PSiRtIjv3i47wiS1Jo1S2040iVJktQBmy5JkqQO2HRJkiR1wGu6NJJ2es5+fUeQpNasWWrDpksj6VmHvrrvCJLUmjVLbXh6UZIkqQM2XRpJ55/yZs4/5c19x5CkVqxZaqNV05XkkCQ3JLkpyQnjzH9hkv9MsjLJEWPmHZPkxubrmIkKLkmSNJmsselKsgFwBnAosDtwdJLdxyy2DJgDnDNm3a2Bk4B9gL2Bk5Jste6xJUmSJpc2I117AzdV1S1V9Qvg88DhwwtU1a1VdRXwqzHrvgRYWFX3VtV9wELgkAnILUmSNKm0abp2AG4ber68mdZGq3WTHJtkcZLFd999d8tNS5IkTR5tPjIi40yrlttvtW5VnQmcCTB79uy229ZamHvWor4jrJVd9z6w7wiS1Jo1S220abqWAzsNPd8RuKPl9pcD+49Z95KW62oae8aLj1jzQpI0IqxZaqPN6cVFwKwkuybZCDgKWNBy+xcAByfZqrmA/uBmmvSoVj78c1Y+/PO+Y2iaSDI/yV1JrhmatnWShc07rxf6JiA9GmuW2lhj01VVK4HjGDRLS4EvVNW1SU5OchhAkr2SLAeOBD6R5Npm3XuBv2HQuC0CTm6mSY/qax85nq995Pi+Y2j6OIvffZPPCcCFVTULuLB5Lo3LmqU2Wt0GqKrOB84fM+09Q48XMTh1ON6684H565BRktarqro0yS5jJh/Oby6POJvBpRHv6CyUpCnHT6SXpPFtV1UrAJrv265uQd+BLakNmy5JWkdVdWZVza6q2TNmzOg7jqQRZdMlSeO7M8n2AM33u3rOI2mSa3VNl9S1Wfu9rO8I0gLgGOCU5vu/9htHo8yapTZsujSSLGDqUpLPMbhofpvmndgnMWi2vpBkLoP7yx7ZX0KNOmuW2rDp0kj6+YM/AWCTzbbsOYmmg6o6ejWzDug0iCYta5ba8JoujaSLzjiBi87wY5EkTQ7WLLVh0yVJktQBmy5JkqQO2HRJkiR1wKZLkiSpA757USPpD/7oFX1HkKTWrFlqw6ZLI2m3fQ7qO4IktWbNUhueXtRIeuieO3nonjv7jiFJrViz1IZNl0bSpZ88iUs/eVLfMSSpFWuW2vD0orSW5p61aI3LzJuzVwdJJEmTiSNdkiRJHbDpkiRJ6oBNlyRJUge8pksj6ZmH/I++I0hSa9YstWHTpZE0c88X9B1BklqzZqkNTy9qJN2/4ofcv+KHfceQpFasWWrDpksj6Vtnv59vnf3+vmNIUivWLLVh0yVJktQBmy5JkqQO2HRJkiR1wKZLkiSpA35khEbSc/7k9X1HkKTWrFlqw6ZLI2mHPfbuO4IktWbNUhueXtRIumfZ97ln2ff7jiFJrViz1IZNl0bSFeecyhXnnNp3DElqxZqlNmy6JEmSOuA1XVPA3LMW9R1BkiStgU2XtB60aYTnzdmrgySSpFHh6UVJkqQOtBrpSnIIcDqwAfCpqjplzPyNgU8DzwPuAV5ZVbcm2QVYCtzQLHp5Vb15YqJrKnveK/687wiS1Jo1S22sselKsgFwBnAQsBxYlGRBVV03tNhc4L6qemqSo4APAK9s5t1cVXtOcG5NcdvNenbfESSpNWuW2mhzenFv4KaquqWqfgF8Hjh8zDKHA2c3j88FDkiSiYup6ebOG6/izhuv6juGJLVizVIbbZquHYDbhp4vb6aNu0xVrQTuB57UzNs1yXeTfCPJC9Yxr6aJK8/7GFee97G+Y0gkuTXJ1UmWJFncdx6NJmuW2mhzTdd4I1bVcpkVwMyquifJ84B/SbJHVT3wWysnxwLHAsycObNFJEnq1B9V1Y/7DiFpcmsz0rUc2Gno+Y7AHatbJsmGwBbAvVX1cFXdA1BVVwI3A08bu4OqOrOqZlfV7BkzZqz9TyFJkjTi2jRdi4BZSXZNshFwFLBgzDILgGOax0cAF1VVJZnRXIhPkt2AWcAtExNdkjpRwNeSXNmMyv+OJMcmWZxk8d13391xPEmTxRpPL1bVyiTHARcw+MiI+VV1bZKTgcVVtQCYB/xTkpuAexk0ZgAvBE5OshJ4BHhzVd27Pn4QSVpP9q2qO5JsCyxMcn1VXTq8QFWdCZwJMHv27LGXX0gS0PJzuqrqfOD8MdPeM/T458CR46x3HnDeOmbUNLTPq97WdwQJgKq6o/l+V5IvM3hH96WPvpamG2uW2vA2QBpJT5r5O5f+SZ1L8kTgcVX1YPP4YODknmNpBFmz1IZNl0bS7dd+B4Ad9ti75ySa5rYDvtx87OCGwDlV9R/9RtIosmapDZsujaTv/dt8wAKmflXVLcBz+s6h0WfNUhve8FqSJKkDNl2SJEkd8PTiCJt71qK+I2g9avPvO2/OXh0kkSR1wZEuSZKkDjjSpZG07zEn9h1BklqzZqkNmy6NpC2237nvCJLUmjVLbXh6USNp2ZJvsmzJN/uOIUmtWLPUhiNdGknX/MdnAZi55wt6TiJJa2bNUhuOdEmSJHXApkuSJKkDNl2SJEkdsOmSJEnqgBfSayS98I3v7TuCJLVmzVIbNl0aSZs+abu+I0hSa9YsteHpRY2kW65YyC1XLOw7hiS1Ys1SG450aSRdf/F5AOy2z0E9J5GkNbNmqQ1HuiRJkjrgSFdP5p61qO8IkiSpQ450SZIkdcCRLmmEtR0RnTdnr/WcRJK0rmy6NJJe/JZT+o4gSa1Zs9SGTZdG0iabbdl3BElqzZqlNrymSyPpxsu+wo2XfaXvGJLUijVLbTjSpZG0qnjN2u9lPSeZOtpcH+a1YdJjY81SG450SZIkdcCRrvXAz+CSJEljOdIlSZLUAUe6pCnA0VVJGn02XWvJP27dOPgvT+s7giS1Zs1SGzZdGkkbbrxJ3xEkqTVrltrwmi6NpKUXncvSi87tO4YktWLNUhuOdDU8bThafvCdrwPwjBcf0XMSSVoza5bacKRLkiSpA61GupIcApwObAB8qqpOGTN/Y+DTwPOAe4BXVtWtzbwTgbnAI8BfVNUFE5YeP2Vb0vq1pvonSW2tselKsgFwBnAQsBxYlGRBVV03tNhc4L6qemqSo4APAK9MsjtwFLAH8GTg60meVlWPTPQP8mg8dSi1M5HHylR4sdOy/klSK21OL+4N3FRVt1TVL4DPA4ePWeZw4Ozm8bnAAUnSTP98VT1cVT8Abmq2J0mTQZv6J0mttDm9uANw29Dz5cA+q1umqlYmuR94UjP98jHr7jB2B0mOBY5tnj6U5IZW6VdvG+DH67gN9z8CGea/7jH36H3/Dqb9/ue/bq32v/N6S7Ju2tS/9VHDJkLf/wfWxWTNvs381+3969zzX9dnlLUyaX/fjEbu1vWrTdOVcaZVy2XarEtVnQmc2SJLK0kWV9Xsidqe+598Gdz/9N7/BOqlhk2EyfxvMFmzm7tbkzF3m9OLy4Gdhp7vCNyxumWSbAhsAdzbcl1JGlXWMEkTpk3TtQiYlWTXJBsxuDB+wZhlFgDHNI+PAC6qqmqmH5Vk4yS7ArOA70xMdEla79rUP0lqZY2nF5trtI4DLmDwlun5VXVtkpOBxVW1AJgH/FOSmxiMcB3VrHttki8A1wErgbd09M7Fvof5p/v+of8M7n96739CrK7+9Ryrrcn8bzBZs5u7W5MudwYDUpIkSVqf/ER6SZKkDth0SZIkdWDKNl1J/i7J9UmuSvLlJFt2tN9DktyQ5KYkJ3Sxz6F975Tk4iRLk1yb5K1d7n8oxwZJvpvkKz3se8sk5zb/9kuT/LeO9/+Xze/+miSfS7JJB/ucn+SuJNcMTds6ycIkNzbft+p4/70cf9PB6o7zJM9J8v+SXJ3k35Jsvpr1ezlGJiB358dWs99Nknwnyfea/b+3mb5rkiuaY+yfmzdajLf+ic3fgxuSvKSLzOuaO8lBSa5s/k2uTPLirnKva/ahbcxM8lCSv+oueQtVNSW/gIOBDZvHHwA+0ME+NwBuBnYDNgK+B+ze4c+8PfCHzePNgO93uf+hHG8DzgG+0sO+zwbe0DzeCNiyw33vAPwAeELz/AvAnA72+0LgD4FrhqZ9EDiheXzC+vz/v5r9d378TZev1R3nDN5p+aJm+uuBv1nN+r0cI+uSu69jq9lXgE2bx48HrgCe32Q4qpn+ceB/jrPu7s3fgY2BXZu/DxtMgtzPBZ7cPH4mcHsXmSci+9A2zgO+CPxVl9nX9DVlR7qq6mtVtbJ5ejmDz9dZ33q9ZUhVraiq/2wePwgsZZw7AKxPSXYE/hj4VJf7bfa9OYMGYB5AVf2iqn7ScYwNgSdk8Hl1v0cHn+lUVZcyeNfwsOFbc50NvLzL/fd0/E0Lj3KcPx24tFlsIfCKsev2eYysS+5G58cWQA081Dx9fPNVwIsZ3PYOVn+M9XYrvHXJXVXfrapVv99rgU2SbLyeIw/vf11+5yR5OXALg+wjZco2XWO8HvhqB/sZ75YhnTY9qyTZhcGrlSs63vVpwP8GftXxfmEwwng38I/N6c1PJXliVzuvqtuBDwHLgBXA/VX1ta72P8Z2VbWiybUC2LanHNDd8TftjDnOrwEOa2YdyW9/qOsqvR4jq6xt7r6PreaSiSXAXQwaw5uBnwy9sFhdre/1b8I65B72CuC7VfXw+kv6ux5r9ub/8zuA93aVdW1M6qYrydeb8/tjvw4fWuZdDD4j7LNdRBpnWuefyZFkUwZDq8dX1QMd7vdlwF1VdWVX+xxjQwanuf6hqp4L/JTBqbVONNdNHc7gNMKTgScmeXVX+x9FHR9/08o4x/nrgbckuZLB6btfjLNar8cIPLbcfR9bVfVIVe3JYMR2b+AZ4y02zrRe/yasQ24AkuzB4PKAN62fhKu3DtnfC3xkaKRspLS59+LIqqoDH21+kmOAlwEHVHOSdz3r/ZYhSR7PoKB9tqq+1OW+gX2Bw5K8FNgE2DzJZ6qqq+K4HFheVatG986l2z8oBwI/qKq7AZJ8CfjvwGc6zLDKnUm2r6oVSbZn8GqxUz0cf9PGeMd5VV3P4Fo6kjyNwWn+sXo9RtYh90gcW1X1kySXMLi+aMskGzYjL6ur9b3/TYDHlHvVpSJfBl5bVTd3FnaMx5B9H+CIJB8EtgR+leTnVfXRzkI/ikk90vVokhzCYIjxsKr6WUe77fWWIUnC4FqNpVV1alf7XaWqTqyqHatqFwY/+0UdNlxU1Y+A25I8vZl0AIO7IXRlGfD8JL/X/FscwOC6lT4M35rrGOBfu9x5T8fftLC64zzJts33xwHvZnCh8W/p8xhZl9z0eGwlmZHm3bdJnsCgAVwKXMzgtnew+mOst1vhrUvuZr1/B06sqm91kXfM/h9z9qp6QVXt0vwdOg1436g0XMCUfveMgWnjAAAA+ElEQVTiTQzOpS9pvj7e0X5fyuBdOTcD7+r4Z96PwXDrVUM/90t7+v3vTz/vXtwTWNz8Dv4F2Krj/b8XuJ7BdSr/BGzcwT4/x+A6l18yeGU9F3gScCFwY/N9647338vxNx2+VnecA29tas/3gVP4zR1HngycP7R+L8fIBOTu/Nhq9vts4LtN7muA9zTTd2PQQN3E4F1yGzfTDwNOHlr/Xc3fgxuAQzv8f/KYczNofn869O+0BNh2MmQfs52/ZsTevehtgCRJkjowZU8vSpIkjRKbLkmSpA7YdEmSJHXApkuSJKkDNl2SJEkdsOmSJEnqgE2XJElSB/4/6UdI/kr67FUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cond = [conditioned_propagate(guess=torch.tensor([5.]), \n",
    "                              guess_sd=torch.tensor([2.]), \n",
    "                              time=torch.tensor([10.])) for i in range(1000)]\n",
    "speed_cond, pos_cond = zip(*cond)\n",
    "plot_hist(speed_cond, pos_cond, \"Speed\", \"Position\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3) Inference\n",
    "The next step is inference- the difficult part. The aim is to find the posterior distribution of the speed, given our model and observation above.\n",
    "$$P({\\sf speed}|{\\sf position}=30) \\sim \\, ?$$\n",
    "Being able to calculate this exactly is rare, so we will estimate the distribution, we can use the stochastic variational inference algorithm. This requires a distribution called a guide (or variational distribution) that is an approximation the posterior distribution we want to know: $P({\\sf speed}|{\\sf position}=30)$. The guide relies on additional parameters that will be tuned to \"guide\" it towards the posterior distribution.\n",
    "\n",
    "For this example, we will sample the wind speed from a normal distribution with mean and variance given by the new parameters, a and b. These are first defined with the statement pyro.param(...) and are stored in the \"Pyro Parameter Store\". The values we provide here for a and b aren't too important, as we will optimise these in the next step."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this example, we will sample the wind speed from a normal distribution with mean and variance given by the new parameters, a and b. These are first defined with the statement `pyro.param(...)` and are stored in the \"Pyro Parameter Store\". The values we provide here for a and b aren't too important, as we will optimise these in the next step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from torch.distributions import constraints\n",
    "def guide(guess, guess_sd, time):\n",
    "    \"\"\" The guide (or variational distribution) approximates \n",
    "    the posterior distribution, using extra variational parameters,\n",
    "    a and b, to sample from the speed distribution. \"\"\"\n",
    "    a = pyro.param(\"a\", torch.tensor(guess))\n",
    "    b = pyro.param(\"b\", torch.tensor(guess_sd), constraint=constraints.positive)\n",
    "    speed = pyro.sample(\"speed\", dist.Normal(a, b))\n",
    "    return speed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check the distribution of the guide before optimising a and b. In our simple example we can actually calculate the exact posterior distribution using rules of normal distributions, to find that posterior distribution of the cloud speed has a mean of 3.4 m/s with a standard deviation of 0.89 m/s, shown in red above. Our initial guess is far from the ground truth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact Mean = 3.40, Exact Standard Deviation = 0.89\n"
     ]
    }
   ],
   "source": [
    "speed_obs_mean = 30./10.\n",
    "speed_obs_sd = 10./10.\n",
    "guess = 5.\n",
    "guess_sd = 2.\n",
    "exact_mean = (speed_obs_sd**2 * guess + guess_sd**2 * speed_obs_mean)/(speed_obs_sd**2 + guess_sd**2) \n",
    "exact_sd = np.sqrt(guess_sd**2/((guess_sd**2 + speed_obs_sd**2)))\n",
    "\n",
    "print(\"Exact Mean = {:.2f}, Exact Standard Deviation = {:.2f}\".format(\n",
    "    exact_mean, exact_sd))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda3/envs/infolab/lib/python3.7/site-packages/ipykernel_launcher.py:6: UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).\n",
      "  \n",
      "/anaconda3/envs/infolab/lib/python3.7/site-packages/ipykernel_launcher.py:7: UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).\n",
      "  import sys\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASYAAAEICAYAAADyYlmcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8VPW9//HXJwkhLIFAEgRZoyKbNcimFcGlaHHfK/baqw/1Ui22F9ve/rRatS6t9br1ca+2oLXYXhVRry1yccValRYEarBssskSlgRCCGsgy+f3x5kJQ5hkzsycySzn83w85nEyM2c+8x1J3n6/3/mec0RVMcaYVJKV7AYYY0xzFkzGmJRjwWSMSTkWTMaYlGPBZIxJORZMxpiUY8Fk0o6IfCQitya7HSZxLJhMVETkLBH5m4jUiMguEZkvIqOT3S6TWXKS3QCTPkSkCzAHuB2YBeQC44BDyWyXyTzWYzLROBlAVV9R1QZVPaiq76nqFyJyU6D39F+B3tQqEflG8IUi0lVEfici20Rki4g8LCLZIc/fLCIrRaRaRN4Vkf4hz50fqFcjIv8NSJt+atPmLJhMNFYDDSLyoohcKCLdmj1/OrAeKALuB/5XRLoHnnsRqAdOAk4DLgBuBRCRK4CfAlcBxcAnwCuB54qAN4B7A3XXAWMT9QFNarBgMq6p6h7gLECB54AdIjJbRI4L7FIJPK2qdar6KvAlcHHg+QuBqaq6X1UrgaeASYHXfRf4paquVNV64BfA8ECv6SJghaq+rqp1wNPA9rb5xCZZLJhMVALhcZOq9gFOAY7HCQuALXr0UeEbA8/3B9oB20Rkt4jsBqYBPQL79Qd+HfLcLpzhWu/A6zeHvL+G3jeZyYLJxExVVwEzcAIKoLeIhM7/9AO24gTJIaBIVQsCty6qOiyw32bguyHPFahqB1X9G7AN6BssGKjfF5PRLJiMayIyWER+JCJ9Avf7AtcDCwK79AB+ICLtRORaYAgwV1W3Ae8BT4hIFxHJEpETReTswOt+C9wtIsMCdbsGXg/wf8AwEblKRHKAHwA92+LzmuSxYDLR2Iszwb1QRPbjBNIy4EeB5xcCA4GdwCPANapaFXjuX3GWF6wAqoHXgV4Aqvom8CtgpojsCdS8MPDcTuBa4FGgKlB/fkI/pUk6sRPFGS+IyE3Arap6VrLbYtKf9ZiMMSnHgskYk3JsKGeMSTnWYzLGpJykHcRbVFSkAwYMSNbbm2h8+aWzHTTIg1JfBkrFX8uknyVLluxU1eJI+yUtmAYMGMDixYuT9fYmGuec42w/+siDUucESsVfy6QfEdnoZj8byhljUo6dj8lEdu+9HpbyrpbJXBZMJrIJEzws5V0tk7ksmExkZWXOdvhwD0qVBUrFXytV1dXVUV5eTm1tbbKbkjR5eXn06dOHdu3axfR6CyYT2dSpztaDCeupgVqZPPldXl5Ofn4+AwYM4OiTLfiDqlJVVUV5eTklJSUx1bDJb2M8VltbS2FhoS9DCUBEKCwsjKvHaMFkTAL4NZSC4v38FkzGnYYG+PxzsEOYTBuwYDKtW7wYVq6E+fNhxAh49dVkt8i4kJ2dzfDhw5tujz76qGe1y8rKmDt3rmf1wrHJb9MyVbj8cqisPPLYW2/BpEktvyaCX/ziFx40zETSoUOHpm9AvVZWVsbixYu56KKLElIfrMdkWrNxI2zdCt27w5w5zmN//Wtcw7kzzzyTM88806MGmmjU1NQwaNCgpuMVr7/+ep577jkAbr/9dkaNGsWwYcO4//77m16zaNEizjzzTEpLSxkzZgw1NTXcd999vPrqqwwfPpxXE9WDVtWk3EaOHKkmxc2apQqqZ5yh+sknqt27O/fXr4+55Pz583X+/PkeNjL1rFix4sgdJ8a9v0WQlZWlpaWlTbeZM2eqqup7772nZ5xxhr7yyiv6zW9+s2n/qqoqVVWtr6/Xs88+W5cuXaqHDh3SkpIS/eyzz1RVtaamRuvq6vT3v/+9TpkyJbr/Dk3/OVisLvLBhnKmZYsWOdtt25zDUsaNgz//2ek1xbg+5ac//SmQ2euYUkFLQ7nzzz+f1157jSlTprB06dKmx2fNmsX06dOpr69n27ZtrFixAhGhV69ejB49GoAuXbq0WfttKGda9tlnzjY/39meHbioyccfJ6c96ShRfaYYNTY2snLlSjp06MCuXbsA+Oqrr3j88ceZN28eX3zxBRdffDG1tbWoatKWPVgwmfAaGmDJEufnYDCNH+9sLZjS1lNPPcWQIUN45ZVXuPnmm6mrq2PPnj106tSJrl27UlFRwdtvvw3A4MGD2bp1K4sCPee9e/dSX19Pfn4+e/fuTWg7LZhMeKtWwb590K8f5OY6j5WWOiG1bh1s2ZLc9plWHTx48KjlAnfddRerV6/m+eef54knnmDcuHGMHz+ehx9+mNLSUk477TSGDRvGzTffzNixYwHIzc3l1Vdf5fvf/z6lpaWcf/751NbWcu6557JixYqETn7bHJMJLzi/NGYM7Njh/JyTA2PHwjvvwCefxLVswCRWQ0ND2MdXrlzZ9POTTz7Z9POMGTPC7j969GgWLFhwzOPBXlSiWDCZ8IK/eKNHwwUXHHl8/HgnmD7+OKZgevrppz1qoMlkFkwmvODE95gxR5/uJLgGKTj/FKVMPt2J8Y7NMZljHToES5eCCIwcCR984NwATjrJ2W7YEFPpDz74gA+CtYxpgfWYzLGWLYO6OhgyxJnsfvhh5/EJE6BXL2jXzjlM5cAB6NgxqtIPB2rZmSxNa6zHZI61fr2zHTz42OeysqB/f+fnja4ueGFM1CyYzLE2bXK2/fqFfz54PcAYh3PGRGJDOXMsCyZP3TLD26/Wf3fTaFf7VVRUcOedd7JgwQK6detGbm4uP/nJT7jyyitbfM1FF13Eyy+/TEFBwVGPP/DAA3Tu3Jkf//jHcbXdLesxmWMFh2gWTGlLVbniiisYP34869evZ8mSJcycOZPy8vJWXzd37txjQikZLJjMsZr3mKZNc25BcQTTtGnTmBZayyTEhx9+SG5uLrfddlvTY/379+f73/8+M2bM4I477mh6/JJLLmk6qHrAgAHs3LkTgEceeYRBgwYxYcKEplOlAKxbt46JEycycuRIxo0bx6pVqzxvvw3lzLGaB9OgQUc/HzyzQAzBNKh5LZMQy5cvZ8SIETG/PtjD+vzzz6mvr2fEiBGMHDkSgMmTJ/Pb3/6WgQMHsnDhQr73ve/x4YcfetV0wGUwichE4NdANvC8qoY9T6eIXAO8BoxW1cWetdK0nf37oarKOT6uRw/nsbfecraXXups4+gxvRWodWmwlmkTU6ZM4dNPPyU3N5cpU6ZE3P+TTz7hyiuvpGNgOchll10GwL59+/jb3/7Gtdde27TvoUOHPG9vxGASkWzgGeB8oBxYJCKzVXVFs/3ygR8ACz1vpWk7mzc72759naUBAE884WyDYdKzpxNclZVOkHXq5Lr8E4FaFkyJNWzYMN54442m+8888ww7d+5k1KhR5OTk0NjY2PRcS5dZCnfKk8bGRgoKChJ22t4gN3NMY4C1qrpeVQ8DM4HLw+z3EPAY4N/Lj2aCSN/Iga1lSgPnnXcetbW1/OY3v2l67MCBA4Azj1RWVkZjYyObN2/ms+DhRyHGjx/Pm2++ycGDB9m7d29TT7dLly6UlJTw2muvAc4ke+gJ57ziZijXG9gccr8cOD10BxE5DeirqnNEpMXvE0VkMjAZoF9rv/gmedwEEzjDuTVrnOHc0KGJblVac/v1vpdEhD/96U/ceeedPPbYYxQXF9OpUyd+9atfMXbsWEpKSvja177GKaecEnYuasSIEVx33XUMHz6c/v37M27cuKbnXnrpJW6//XYefvhh6urqmDRpEqWlpZ62300whTuFXdMp9EQkC3gKuClSIVWdDkwHGDVqlF2gLBVFE0xgSwZSWK9evZg5c2bY51566aWwj28I+fe85557uOeee47Zp6SkhHfeeceTNrbEzVCuHOgbcr8PsDXkfj5wCvCRiGwAzgBmi8gorxpp2pAFk0kBbnpMi4CBIlICbAEmAd8OPqmqNUBR8L6IfAT82L6VS1PhgumPfzx2vxiD6Y/hahnTTMRgUtV6EbkDeBdnucALqrpcRB7EuRTL7EQ30rShcMHUt++x+8UYTH3D1cpAyTyRfyrQOC8l72odk6rOBeY2e+y+FvY9J64WmeRpbDyyXCA0mILndb7uuiOPxRhMwXNEXxdaK8Pk5eVRVVVFYWGhL8NJVamqqiIvLy/mGrby2xxRUQGHD0NR0dHnWQp+5RwaJsG1TDt2RLWWKfj1dSYHU58+fSgvL2dH8FzpPpSXl0efPn1ifr0FkznC7cQ3OGuZeveGr75yLogZPLOloV27dpTEeEFQ47CDeM0R0QQTOL0mgO3bE9Me41sWTOYICyaTIiyYzBGhx8m5YcFkEsTmmMwRwYDp1evox19/Pfz+MQTT6y3VMiaEBZM5oqLC2R533NGPFxUdu2/oflEEU1FLtYwJYUM5c0QwYJoH04wZzq25GHpMM2bMaPFy1MYEWTCZI1rqMUUKpuDrXLBgMm5YMBnH4cNQXe2sTyosdPcam/w2CWLBZByVlc62Rw/Iznb3mmDPqqLCOZzFGI9YMBlHS8O41uTlQUGBcznx6urEtMv4kgWTccQSTGDDOZMQtlzAOFr6Rg5g7txjHwvq2RNWrXJeP2xYxLeZ21otYwIsmIyjtR5T6JkGmotyLVPH1moZE2BDOeMIBlNwaBbq2WedWzhRDuWeffZZnm2pljEBFkzG0VqPadYs5xZOlGuZZs2axayWahkTYMFkHDb5bVKIBZNxWDCZFGLBZBwWTCaFWDAZZ4FkVZVzOEq0R/9bMJkEsOUC5sjhKMXF4Q9H+eijll9bVAQisHOnE3Dt2rX6Vh+1VsuYAOsxmdiHcQA5OU6gqTpXTDHGAxZMJnIwPf64c2tJFMO5xx9/nMdbq2UMFkwGIgfTnDnOrSVRrGWaM2cOc1qrZQwWTAbiG8rBkWDats2b9hjfs2AyR4Zg4Q5HcaNHD2cbnEQ3Jk4WTCb+HlNxsbO1yW/jEVsuYCIHU4cOrb8+imDqEKmWMVgwGYgcTG+/3frrg0M5F8H0dqRaxmBDOQNH5oZsKGdShAWT3zU0OKu2oeWrozz0kHNrSRTB9NBDD/FQa7WMwYLJ7NrlrNru1q3lw0nmzXNuLQkGU2WlU6sV8+bNY15rtYzBgskEeznBcIlFp07OFVNqa2H/fm/aZXzNgsnvvAgmkagmwI2JxILJ77wIptDXWzAZD9hyAb9zE0xuLhnuMpgK3V5+3Piaq2ASkYnAr4Fs4HlVfbTZ87cBU4AGYB8wWVVXeNxWkwhugumNNyLXCZ0Ab8UbbmoZ34s4lBORbOAZ4EJgKHC9iAxtttvLqvo1VR0OPAY86XlLTWLYUM6kIDdzTGOAtaq6XlUPAzOBy0N3UNU9IXc7Aa1/Z2xSh5tguvtu59Yal8F09913c3ekWsb33AzlegObQ+6XA6c330lEpgA/BHKB8zxpnUk8N8H0979HruPyW7m/u6llfM9Nj0nCPHZMj0hVn1HVE4H/B9wbtpDIZBFZLCKLd1iXPzXYUM6kIDfBVA70DbnfB9jayv4zgSvCPaGq01V1lKqOKo73D8F4w4LJpCA3wbQIGCgiJSKSC0wCZofuICIDQ+5eDKzxrokmYVSPHCfnVTDZyeKMByLOMalqvYjcAbyLs1zgBVVdLiIPAotVdTZwh4hMAOqAauDGRDbaeGT3bqivh/x8aN++5f369Ilcy2WPqY+bWsb3RCMcdJkoo0aN0sWLFyflvU3A6tUwaBCccAKsWxdfLVXnhHKHDjnHy3Xs6E0bTUYRkSWqOirSfnZIip+FXugyXiI2z2Q8Y8HkZ24nvqdOdW6RuAimqVOnMtVNLeNrdqycn7kNprIyd/VcBFOZ21rG16zH5GfBAAkujoyXfTNnPGLB5GderWEKsjkm4xELJj/zOpjsZHHGIzbH5Gdug+nkk93Vc9FjOtltLeNrFkx+5jaYpk93V89FME13W8v4mg3l/MzmmEyKsmDyK1X3wTR5snOLxMW3cpMnT2aym1rG12wo51d798Lhw86hI5EOH1m92l1NFz2m1W5rGV+zHpNfeT2MA+ja1blo5r59zjXmjImRBZNfJSKY7Hg54xELJr9KRDCF1rNgMnGwOSa/Ck5QuzkcZfhw93UjBNPwaGoZ37Jg8qtoekxPP+2+boRv5p6OppbxLRvK+ZWX52IKZUM54wELJr+K5swCN9zg3NyIcLzcDTfcwA1uaxnfsqGcX0UzlCsvd183Qo+pPJpaxresx+RX0Ux+R8OGcsYDFkx+ZcsFTAqzYPKjaI6Ti5adxdJ4wOaY/GjvXucySx07QqdOkff/+tfd147QY/p6NLWMb9l15fxo3To46STo3x82bPC2dmOjc/HM+nrneLnWLqRpfMeuK2da5vVFCEJlZUFRkfNz8PLjxkTJgsmPol1cefXVzs2tVoZzV199NVdHU8v4ks0x+VG0E99VVdHVbyWYqqKtZXzJekx+lMihHNg3cyZuFkx+lKjj5IJsLZOJkwWTHyVqDVOQBZOJk80x+VG0h6N84xvR1W/lQN5vRFvL+JIFkx9F22P62c+iqx/SY7plxqKjn+s7EYBbZizidzeNjq6u8Q0byvlRW01+21DOxMiCyW9Uo5/8vvBC5+ZWK9/Kvfvkv/Puk//uvpbxJRvK+U0015MLOngwuvdopcfUcPhQdLWML1mPyW8SPYwD6N7dOTRl926y6+sT9z4mY1kw+U2i1zDBUcfL5e+tTtz7mIxlweQ3bdFjAjjuOAC67NmV2PcxGcnVHJOITAR+DWQDz6vqo82e/yFwK1AP7ABuVtWNHrfVeCGWxZWXXBL9+/TsCf/85zHB1Lf0rOhrGd+JGEwikg08A5wPlAOLRGS2qq4I2e1zYJSqHhCR24HHgOsS0WATp1iGcj/+cfTvE+gxda05+qDdr11oV0gxkbkZyo0B1qrqelU9DMwELg/dQVX/oqoHAncXAH28babxTKIPRwkKDuVq7GwCJnpugqk3sDnkfnngsZbcArwd7gkRmSwii0Vk8Q5bfJcc27c725493b/mnHOcWzRamGOa++htzH30tuhqGd9xM8ckYR4Lez5eEbkBGAWcHe55VZ0OTAfn1Lou22i8FEswxSJQv7XJ72MOVwlhh6v4m5tgKgf6htzvA2xtvpOITADuAc5WVVtFl6oqKpxtooOphTkmY9xwM5RbBAwUkRIRyQUmAbNDdxCR04BpwGWqamcHS2Vt1WOy5QImDhGDSVXrgTuAd4GVwCxVXS4iD4rIZYHd/hPoDLwmImUiMruFciaZDh2CXbsgOxsKCxP7Xk2T3xZMJnqu1jGp6lxgbrPH7gv5eYLH7TKJEHoepqwo1tZ+61vRv1dxMWRl0Wl/Ddn19TTkOL9qJWPsV8VEZgfx+kms80vf+17075WdDUVFZFVWkr+3mt3dnOUJQ867JvpaxnfskBQ/iXV+6cAB5xatMPNM9YdqqT9UG30t4ysWTH4SDKZAYLh20UXOLVphgum9p6by3lNTo69lfMWCyU/a6hu5oOBaJlsyYKJkweQnbR1MtmTAxMiCyU/aanFlkC2yNDGyYPKTWOeYYmU9JhMjWy7gJ7EO5W66Kbb3C3O83MCzYji3k/EdCyY/iTOYoj7oNszqbwsm44YFk1/s3w/79kH79tC1a3Sv3bkz4i7hQqtLzU6e4ug5ptq9uwHIyy+Irg3GV2yOyS+CE9/HHQcS7kw2rbjmGucWpX2dC2gUodP+GrIanKulfPjMXXz4zF1R1zL+YsHkF229VABozM5hb343slTJD/SUjHHDgskvkhBMAHu6dAfsmzkTHZtj8ou2XsMUEAymrjVVR52fOVZ21kt/sGDyi7ZewxSwp6tz3qdoD0tpLYBM5rNg8ot4hnK33+5sD0b/0t0FTjAV7Ha+2Rt87tXRFzG+Y8HkF/EE03WBSwTG0Iup7ub00LpVOyepO+H086N/f+M7Fkx+EU8wbY59dih4grhgMO2rcua6Ohe27ZDSpBf7Vs4v4plj+s53nFsMdnXrARwJpo+fu5+Pn7s/plrGPyyY/KChAbYGrrjVu7VrlXqvunsgmHbZxXOMexZMflBRAfX1UFQEeXlt+tZ7unSnUbLI31tNdn1dm763SV8WTH5QXu5s+/Ztfb8EaMzOYXdBEVmqTd/MGROJBZMfBCev+/RJyts3nwA3JhL7Vs4P4u0x/ehHzjbGE1FWd+sBLKegupJTJv5LbEWMr1gw+UGwxxRrMF16qbONcTV2deCbue7VlfT7pgWTicyCyQ/iHcp9+WVcb18dHMrtqqRm20YAuvbqH1dNk9lsjskP4h3Kffe7zi1GwR5TQfUO5r/4S+a/+MuYaxl/sGDyg3iHcnGqbrbI0phILJgyXRIXVwYFF1l2t2AyLlkwZbrt251w6tHDOd93EgTnmLru3oGoJqUNJr1YMGW64PxSktYwAdS3a8/ezl3JaWigKHDub2NaY9/KZTov5pfuvdfZlsdeorpbD/L31XD21yeypqd9I2daZz2mTOdFME2Y4NzisDswAX5ytx70HjYmrlom81mPKdN5MZQrK4u7GcHTn+SuW0ZVUS8K+50cd02TuazHlOm86DFNnerc4hD8Zu7wgndY+PKTcdUymc96TJkuyh5TuIsA/Mf2vXE3Y3eB881cTzv1iXHBekyZLsmLK4OCPabedYeT2g6THlwFk4hMFJEvRWStiBxzfWcRGS8i/xCRehGJ/lrSJjHq65O+uDKosofTY+tfdyip7TDpIWIwiUg28AxwITAUuF5EhjbbbRNwE/Cy1w00cdi+HRobnfN85+YmtSm7uvekISubXvV1tG9sTGpbTOpzM8c0BlirqusBRGQmcDmwIriDqm4IPGe/calk0yZnG+cw7n+v/l7cTWnIyaGqsCc9dmxh4rlXYTNNpjVugqk3HHV153Lg9FjeTEQmA5MB+vXrF0sJE401a5ztSSfFVWbdwFM9aIwznOuxYwtDO3RiqScVj2aXD88cbuaYJMxjMR3wpKrTVXWUqo4qLi6OpYSJRjCYBg6Mq8yJa77gxDVfxN2cHYF5przln8Vdy2Q2N8FUDoSOBfoAWxPTHOMpj4Lpqjee5ao3no27OcEJ8JzFH8Zdy2Q2N8G0CBgoIiUikgtMAmYntlnGEx4Fk1eCwTTAvpkzEUQMJlWtB+4A3gVWArNUdbmIPCgilwGIyGgRKQeuBaaJyPJENtq4oJpywRQcyvU/bMFkWudq5beqzgXmNnvsvpCfF+EM8UyqqKiAffugWzcoLEx2awDYUeyspepXd4ishnoas+3AAxOerfzOVCnWWwI43D6PbTntaAd031WR7OaYFGb/y8pUHgbTzG//MO4aQbt6n0CvjV9SXLmFncVttxq9paUEtowgNVmPKVN5GEyb+53MZo9OU1Ld12lPj8rNEfY0fmY9pkwVIZhaW4zY3JDAuqOVHpzgbS1wFtCjckvctUzmsmDKVB72mC556wXAm2BasG4ZNwHFlXGcp9dkPBvKZSJVWLvW+TmFJr8BNrRzrtTSw4LJtMKCKRNt3QoHDkBRERQUJLs1R9mQ6wRT8Y4tToAaE4YFUyZKwaUCQTXZOezKzibv0EFbMmBaZMGUiVI4mACWt+8IQN/Na5LcEpOqLJgy0apVztajYPrDjXfzhxvv9qTW2BvvpnrE2QD03bTak5om81gwZaJ//MPZDh/uSbmKXv2p6OXNRSq79upPxRBnUWO/jV96UtNkHgumTNPYeCSYRo70pGRp2SeUln3iSa1NZZ/w2QHnqiv9NluPyYRn65gyzbp1sGcP9Orl3DxwwTsvAbB0+Li4ay175yVWqVKXk0vxjq10OLCXgx3z465rMosFU6ZZssTZjhwZ1erutlQvwpY+JzBgwyr6bl7D6kEjkt0kk2JsKJdpQoIplW3qNwiAfhttOGeOZT2mTBMaTFXJbUprNgUOCu6b5Hkmu4BBarIeUyZR9XziO1GCZyvoZ0sGTBjWY8ok69ZBTY0z6X388YA3R/A//28/96QOwPhArfKOnWkU4fgt68mur6Mhp51n72HSnwVTmgsdioxe+D63AWU9TuS/PJz4ri48zrNanQO1anEuTtCzYjPHb/3Ks/M9mcxgQ7kM0n/jSgA2Dhjsad3RC99n9ML3Pam1fuH7rA/UCoZRyVcrWnuJ8SELpgwyYINzKMrGAUM8rXvOX97gnL+84UmtVX95g1WBWqtPPg2AIStSc1mDSR4LpgyRU3eIE9YtA+CrEm+DKVGWnXIGAEOXL0QaG5LcGpNKLJgyxOBV/6D94Vo29hvEnq5FyW6OK5XH9WVHUS86799D/w123Jw5woIpQ5y69FMAvigdm+SWREGE5YFe0ynL/p7kxphUYsGUCVQ59Yv5QJoFE7B8mBNMw5YvTHJLTCqx5QIZoNe2DRTv2MrezgV8VTLU8/q/mfKoZ7XOa1Zr5dBRNGRlc8K6f5J3cB+1HTp79l4mfVmPKQMEh3H/PPXraFa25/X35RewL9+bc4fn5ReQF1LrYMd81p8wjJyGBgavXOLJe5j0Z8GUAU5d6gzjlpaelZD6Yz+dw9hP53hSa82nc1jTrNaKYacDNs9kjrBgSnMd9+9h4JqlNGRlN00ke+3MT+dwZgKDqSxwnqcxn31A7qFaT97HpDcLpjQ37uPZZDc2sGrwyLQ94dqmAYNZf8IwOu3fwxkL3kl2c0wKsMnvNNDSqTmy6+t49P2ZALx/wfVt2STPfTDhOiZPv4/zPpjFx+MvB5FkN8lOiZJE1mNKY6cveJfu1ZWU9z6Rf556ZrKbE5fFo79BTZfu9C1fy8lffp7s5pgks2BKU9LYyMR3/geAdy68ISV6GPFoyGnHX8+5CoBvzJuV5NaYZLOhXJoqLfuE3lvWs6tbDz47/YKEvtev73zas1oXtFLro3Ov4qL/+z0jlnzESWuWsnZgqWfva9KLBVMa6lKzk3998ZcAvPfNbyf8JGuH2+d5ViunlVo1BUW8+80buHjui9z63P2CZ1SlAAAFbklEQVQ88PP/SdkFly3NP9nckzcsmFKIm6uaZDXU893f/oyue3axavBI5k34VsLbde6HrwPwl/OuibvWykCtIS3U+vOVkzll2QL6b/qS619+it/f8rO437Mt2YS5NyyY2lg8l1SSxgaum/k0g1ctYXfXQqbd9hCN2Yn/Jxz12QeAN8H0VaBWS8HUkNOO5yb/nPt+fiNnffoWO4qPZ86lN6f9HJqJjk1+p4mC6h386D/vYMIHs2jIymbabY+kzelNorWt9wm8dMN/0CjClW9O47bf/JT2tQeS3SzThlz971ZEJgK/BrKB51X10WbPtwf+AIzEuWjQdaq6wdum+lPPbRsY/9c/cdanc+i0fw81Xbrz3OQHWT04sy8S+en4y9jTpRuTp93H6EXzGLxyCR+dexUfn30Fuwp7Jrt5JsEiBpOIZAPPAOcD5cAiEZmtqqEnar4FqFbVk0RkEvAr4DpPWrhjB2za5EmpWDz41vIWn7vv0mHhn1AN/5gqJeuWIyiiSlZjA9n19eQ01JN7qJb2hw+Sv7earrurOG77JvpvXEX36sqmEsuGnc7v/u3+jO0pNffF8HE88rMXuPl3D3LC+uVc+tYLXPrWC1QV9uSrkqHsKD6e3QU92N8pn9q8jhzOzaMhO4eGnBwaJZvGrCwQQXGGgdp8NNjGw8MHH1jZ4nMt/S7F9PvX1oYMgY4dPS3ppsc0BlirqusBRGQmcDkQGkyXAw8Efn4d+G8REdVwf6FR+tOfYPLkuMvE6r7Wnozhqkb3Rrn/gQ6dWTRmAn89+wo2JuCUJqlu2/ElPHLvC5y0Zinnvz+TocsXUli1ncKq7clumrda+F3y+vcvIcrKoNTbpR0SKTtE5BpgoqreGrj/HeB0Vb0jZJ9lgX3KA/fXBfbZ2azWZCCYMoOAaM6nWgTsjLhXevPDZwT7nJkk2s/YX1WLI+3kpscUrr/bPM3c7IOqTgemu3jPYxshslhVR8Xy2nThh88I9jkzSaI+o5tv5cqBviH3+wBbW9pHRHKArsAuLxpojPEfN8G0CBgoIiUikgtMAmY322c2cGPg52uADz2ZXzLG+FLEoZyq1ovIHcC7OMsFXlDV5SLyILBYVWcDvwP+KCJrcXpKkxLQ1piGgGnGD58R7HNmkoR8xoiT38YY09Zs5bcxJuVYMBljUk7aBJOIPCAiW0SkLHC7KNlt8pKITBSRL0VkrYjclez2JIqIbBCRfwb+DRcnuz1eEJEXRKQysJ4v+Fh3EXlfRNYEtt2S2UYvtPA5E/J3mTbBFPCUqg4P3OYmuzFeCTns50JgKHC9iGTyMu9zA/+GmbLGZwYwsdljdwHzVHUgMC9wP93N4NjPCQn4u0y3YMpUTYf9qOphIHjYj0kDqvoxx67buxx4MfDzi8AVbdqoBGjhcyZEugXTHSLyRaBLmfZd4xC9gc0h98sDj2UiBd4TkSWBQ5Qy1XGqug0gsO2R5PYkkud/lykVTCLygYgsC3O7HPgNcCIwHNgGPJHUxnrL1SE9GWKsqo7AGbZOEZHxyW6QiUtC/i5T6gyWqjrBzX4i8hzgzaVhU4Obw34ygqpuDWwrReRNnGHsx8ltVUJUiEgvVd0mIr2AyoivSEOqWhH82cu/y5TqMbUm8I8bdCWwrKV905Cbw37Snoh0EpH84M/ABWTWv2Oo0MO0bgT+nMS2JEyi/i5TqscUwWMiMhxniLMB+G5ym+Odlg77SXKzEuE44E1xTtCWA7ysqml/TXAReQU4BygSkXLgfuBRYJaI3AJsAq5NXgu90cLnPCcRf5d2SIoxJuWkzVDOGOMfFkzGmJRjwWSMSTkWTMaYlGPBZIxJORZMxpiUY8FkjEk5/x/Wj0X4Diq1egAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Distribution of guide before any optimisation\n",
    "speed_guide = [guide(guess=torch.tensor([5.]), \n",
    "               guess_sd=torch.tensor([2.]), \n",
    "               time=torch.tensor([10.])).item() for i in range(1000)]\n",
    "plot_hist(speed_guide,title1=\"Speed\",true_dist=(3.4,0.89))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stochastic Variational Inference\n",
    "Next we want to optimise the guide function to get a the best approximation to the posterior distribution, according to some loss function- we'll use the Evidence Lower Bound Operator or \"ELBO\" (link to why). We use stochastic variational inference , `pyro.infer.SVI` :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 29.6 s, sys: 468 ms, total: 30.1 s\n",
      "Wall time: 30.3 s\n"
     ]
    }
   ],
   "source": [
    "%%time \n",
    "pyro.clear_param_store()\n",
    "import pyro.infer\n",
    "import pyro.optim \n",
    "adam_params = {\"lr\": 0.001, \"betas\": (0.90, 0.999)}\n",
    "optimizer = pyro.optim.Adam(adam_params)\n",
    "svi = pyro.infer.SVI(model=conditioned_propagate, \n",
    "                     guide=guide,\n",
    "                     optim=optimizer,\n",
    "                     loss=pyro.infer.Trace_ELBO())\n",
    "\n",
    "losses, a, b  = [], [], []\n",
    "num_steps = 20000\n",
    "for t in range(num_steps):\n",
    "    losses.append(svi.step(guess=5., guess_sd=2., time=10.))\n",
    "    a.append(pyro.param(\"a\").item())\n",
    "    b.append(pyro.param(\"b\").item())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(tensor(3.4276, requires_grad=True), tensor(0.9108, grad_fn=<AddBackward0>))"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pyro.param(\"a\"), pyro.param(\"b\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAGoCAYAAABbtxOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xl8XXWd//HXJ2vTNm3SNt0LaVkKFGiBsskiCAICgozb6KigOIyOjvobZxSXUdBRUXHXgQFhABdwQWRTKELLUii1hZbSfUv3JUnbNM2e3M/vj3NSTkNucpPeJTf3/Xw87iPnnvM953zON/ee+7nf+z3fY+6OiIiIiIgE8jIdgIiIiIjIQKIEWUREREQkQgmyiIiIiEiEEmQRERERkQglyCIiIiIiEUqQRUREREQilCCLiIiIiEQoQRYRGWDMrMrMmszsQOTxczO7zsxeiLPOPDNrDsvWmdlzZnZSlzInmNkj4fJ6M5trZm9Jz1GJiGQPJcgiIgPTO919eOTx6QTW+bS7DwdGA/OAX3UuMLOjgPnAMmAqMBF4CJhjZmcnPXoRkSymBFlEZJBx93bgAeCEyOybgJfc/Svuvsfd6939pwRJ9HczEKaIyIClBFlEZJAxsyLgn4AFkdlvB/7QTfHfA+eY2dB0xCYikg0KMh2AiIh0689m1h55/p9AWy/r/NTMbgWGAk3AP0SWjQF2dLPODoLGknKgsf/hiogMHmpBFhEZmN7l7mWRx50JrPMZdy8DhgBXAn80s5PDZTXAhG7WmQDEgL1JiVpEZBBQgiwiMsi4e8zdnwfWAZeEs/8GvLeb4u8j6Jus1mMRkZC6WIiIZBczsyHRGe7e3E2hswku0lsezroZ+LuZfQv4AUF3jeuAj/BGEi0iIqgFWURkoHq0yzjID4Xz30LQv/jgw8w6Gzt+3lmeYHSKr7r7XwHcfS1wLjATqCLoe/xu4FJ3n5+2oxIRyQLm7pmOQURERERkwFALsoiIiIhIhBJkEREREZEIJcgiIiIiIhFKkEVEREREIrJumLcxY8Z4ZWVlpsMQERERkSyyePHiGnevSKRs1iXIlZWVLFq0KNNhiIiIiEgWMbNNiZZVFwsRERERkQglyCIiIiIiEUqQRUREREQilCCLiIiIiEQoQe6HVTv3s31fU6bDEBEREZEUyLpRLAaCy378PABVt1yR4UhEREREJNnUgiwiIiIiEqEEWUREREQkIi0JspkNMbOFZrbUzJab2c3h/Klm9rKZrTWz35lZUTriERERERGJJ10tyC3A29x9JjALuMzMzgK+C/zI3Y8B9gLXpykeEREREZFupSVB9sCB8Glh+HDgbcAfw/n3Au9KRzwiIiIiIvGkrQ+ymeWb2RJgN/AUsB7Y5+7tYZGtwKQ4695gZovMbFF1dXV6AhYRERGRnJS2BNndO9x9FjAZOAM4vrticda9w91nu/vsioqKVIYpIiIiIjku7aNYuPs+YB5wFlBmZp1jMU8Gtqc7HhERERGRqHSNYlFhZmXhdAlwMbASmAu8Jyx2LfBwOuIREREREYknXXfSmwDca2b5BEn57939MTNbATxgZv8NvArclaZ4RERERES6lZYE2d1fA07pZv4Ggv7IIiIiIiIDgu6kJyIiIiISoQRZRERERCRCCbKIiIiISIQSZBERERGRCCXIh2l99QHO/e4zVNe3ZDoUEREREUkCJciH6ZfPb2Tr3ibmrNiZ6VBEREREJAmUIIuIiIiIRChBFhERERGJUIIsIiIiIhKhBFlEREREJEIJsoiIiIhIhBJkEREREZEIJcgiIiIiIhFKkEVEREREIpQgi4iIiIhEKEEWEREREYlIS4JsZlPMbK6ZrTSz5Wb22XD+TWa2zcyWhI/L0xGPiIiIiEg8BWnaTzvweXd/xcxKgcVm9lS47Efufmua4hARERER6VFaEmR33wHsCKfrzWwlMCkd+xYRERER6Yu090E2s0rgFODlcNanzew1M7vbzMrjrHODmS0ys0XV1dVpirRv3DMdgYiIiIgkQ1oTZDMbDjwIfM7d9wO3AUcBswhamH/Q3Xrufoe7z3b32RUVFWmLV0RERERyT9oSZDMrJEiOf+PufwJw913u3uHuMeBO4Ix0xZNsZpmOQERERESSIV2jWBhwF7DS3X8YmT8hUuwa4PV0xCMiIiIiEk+6RrE4B/gwsMzMloTzvgx8wMxmAQ5UAf+SpnhERERERLqVrlEsXgC664Twl3TsPxliMWdjbQNHVQzPdCgiIiIikkK6k16CfvbMOi76wbOs2VWf6VBEREREJIWUICdo0aY9AOyoa85wJCIiIiKSSkqQRUREREQilCCLiIiIiEQoQRYRERERiVCCLCIiIiISoQRZRERERCRCCXISPfTqVrbsacx0GCIiIiJyGJQgJ4k7/L/fLeWa/5mf6VBERERE5DAoQT5M1uX+gDUHWjMTiIiIiIgkhRJkEREREZEIJch91KXBGPeMhCEiIiIiKaIEWUREREQkQgmyiIiIiEiEEuQ+itejQj0tRERERAYHJch99PTKXT0u7zqqhYiIiIhkl7QkyGY2xczmmtlKM1tuZp8N548ys6fMbG34tzwd8RyO+17alOkQRERERCSF0tWC3A583t2PB84CPmVmJwA3Ak+7+zHA0+HzAae9I8aBlvaeC2k4CxEREZFBIS0JsrvvcPdXwul6YCUwCbgauDcsdi/wrnTE01efeeBVXt28L6Gy6mEhIiIikt3S3gfZzCqBU4CXgXHuvgOCJBoYG2edG8xskZktqq6uTleoB/1l2c6071NEREREMiOtCbKZDQceBD7n7vsTXc/d73D32e4+u6KiInUBHobODhYx9bQQERERyWppS5DNrJAgOf6Nu/8pnL3LzCaEyycAu9MVT7KpC7KIiIjI4JCuUSwMuAtY6e4/jCx6BLg2nL4WeDgd8SSLu3P/ws0Hp0VEREQk+xWkaT/nAB8GlpnZknDel4FbgN+b2fXAZuC9aYonKZZvf6OXiNJjERERkcEhLQmyu79A/AEeLkpHDKnQoQ7HIiIiIoOO7qQnIiIiIhLR5wTZzD5rZiMscJeZvWJml6QiuGyiLsgiIiIig0N/WpA/Fg7RdglQAXyUoC+xiIiIiEjW60+C3NmX+HLg/9x9KTl6Azmz7qdFREREJHv1J0FebGZzCBLkJ82sFIglN6zsoy4WIiIiIoNDf0axuB6YBWxw90YzG0XQzUJEREREJOv1pwX5bGC1u+8zsw8BXwXqkhtWdmho6ch0CCIiIiKSZP1JkG8DGs1sJvAFYBNwX1KjyhK/X7Sl2/lLt+xLcyQiIiIikiz9SZDbPbiv8tXAT9z9J0BpcsPKPtEuyFf/Yn7G4hARERGRw9OfPsj1ZvYlgltHn2dm+UBhcsPKPq6r9EREREQGhf60IL8faCEYD3knMAn4flKjEhERERHJkD4nyGFS/BtgpJldCTS7e072QW7tyPnR7UREREQGnf7cavp9wELgvcD7gJfN7D3JDiwbPP7ajoPTifaw+N9n11N54+PqkiEiIiIyQPWni8VXgNPd/Vp3/whwBvBfyQ1r8Pi/+RupvPFx2sLW5u8+sQqAmENHzGlp11BxIiIiIgNJfxLkPHffHXle28/tDCpO9y3CP3pqDQCN3YyZ/MlfL2b6V59IaVwiIiIi0jf9SWyfMLMnzew6M7sOeBz4S08rmNndZrbbzF6PzLvJzLaZ2ZLwcXk/Yslqc1bsynQIIiIiItJFn4d5c/f/NLN3A+cABtzh7g/1sto9wM958w1FfuTut/Y1hsFAfZBFREREBqb+jIOMuz8IPNiH8s+ZWWV/9pUtest3HWdjTQOxsFxbR/cr1B5owcwYNawoyRGKiIiISCISTpDNrB667WhrgLv7iH7s/9Nm9hFgEfB5d9/bj20MCE1t3V9sZ2YHpy+8dd7B6W88tqLb8qf9998AqLrliuQFJyIiIiIJS7gPsruXuvuIbh6l/UyObwOOAmYBO4AfxCtoZjeY2SIzW1RdXd2PXaVerI89Jl7dnLXfBUREREQGtYyNPuHuu9y9w91jwJ0Ew8XFK3uHu89299kVFRXpC7Iv1KdYREREZFDIWIJsZhMiT68BXo9XdjBatbP+4HSb7sgnIiIiMmCkJUE2s/uBl4DpZrbVzK4Hvmdmy8zsNeBC4P+lI5ZMqa5vibvs1idXpzESEREREelJv0ax6Ct3/0A3s+9Kx77TpbcOFm//0XNxl63ZVR93mYiIiIikV87fAS9Z1AVZREREZHBQgpwkkdHcRERERCSLKUFOklicJuS6prZe1527OvGh637/9y3srGtOuLyIiIiI9I0S5CT5xdz1Kd/H3oZWvvDga1x798KU70tEREQkVylBHqCeWrGLe1+sOmRee3g3ktqG+CNiAPxi7jo21TakKjQRERGRQU0Jcgr1NLRbb/75vkV8/ZHlfV5v9/5mvv/karUyi4iIiPSTEuQU+t9nk9vtwnsdTO6NW143tXUkdd8iIiIiuUIJclbqfcgMDTsnIiIi0j9KkLOS8/CSbazauf9NSzTcnIiIiMjhScud9HJRW0esT+V3728mPy+x7LbmQCuffWAJAFW3XNHn2EREREQkPiXIKdLXPsBnfPvppO5fPSxERERE+kddLFKkrrGNX76wMdNhiIiIiEgfqQU5Rb728OtJ3+b/za+Ku+zVzXuJOUwuL0n6fkVERERyiRLkFDnQ0p7U7bk7t82LP2zcNf/zIgAPf+qcpO5XREREJNeoi8Ugc/Uv5mc6BBEREZGspgQ5RZI9DnFft9fXUTREREREJJC2BNnM7jaz3Wb2emTeKDN7yszWhn/L0xVPNqlrbKMtFj/h3bKn8U3zmlp1Jz0RERGR/khnC/I9wGVd5t0IPO3uxwBPh88HhWQ1ILs7M78xh+lffaLb5S9vqOW8781N2v7nrd5NXWNbP9cWERERyX5pS5Dd/TlgT5fZVwP3htP3Au9KVzyp5knqY3HyTXN6XL5mV32cAPq+r70NrVz3f3/nhl8t6vvKIiIiIoNEpvsgj3P3HQDh37EZjmfAqe/vaBj9uOV0S3vQjWNjTUP/9ikiIiIyCGQ6QU6Imd1gZovMbFF1dXWmw0mIWT8y1H745mMru18QaUH+0C9f5puPreh1Wx6ulKbQRURERAakTCfIu8xsAkD4d3d3hdz9Dnef7e6zKyoq0hpgfyWri0VvWhMYreKFdTXc1cNd/ZpaO7jvpaqDI2VYf5qfRURERAaJTN8o5BHgWuCW8O/DmQ0neTpi6UmQ44m5U9/cRumQwl7LfveJVdzzYhWxmFqQRURERNI5zNv9wEvAdDPbambXEyTGbzeztcDbw+eDQmd/3kxpjzkn3TSHrXvfPARcV3saWgFoCIeGi+bHdY1tbNvXlIoQRURERAaktLUgu/sH4iy6KF0xpFO6+iD3Zuve3pPbmB/achyN/cIfzGNPQytVt1yRkvhEREREBppM90EetKrrWzIdQo+WbNnHzY8u77GvdM2BloOtyyIiIiK5ItN9kAetmgMDI0H+y7Id3c5/z20v0h5z3nbc2IP9pQ9epBc2IH/nL6vSEaKIiIjIgKIEeZC776VNPS7/8F0L3zSvM0H2pN0PUERERCR7qItFDuqIOe3djLJxsA9yeJlee4cSZBEREck9SpBzSOWNj1N54+N878nuu0507WKxaU/vI2DsqGui9TBH7Hh1817ed/tLtLR3HNZ2RERERJJBCXIOmreq+7sRbqgObjF9cAyLXm528quXqjj7O8/wpT8tO2T+xpoGtiSQXHf60p+WsbBqD+t36xbXIiIiknlKkOWgB1/ZCkBVbSOxmLN0a90hy92dtvDOfdX1LfzXw8sBeHrVLto7YmyuDZLiC2+dx3nfm9vn/b/vf186nPBFREREkkIJcg460NLea5lfvrDhTfNunbOaY77yV5rbOg65U6AB3/nrKs7//lx21jWnNC4RERGRVFOCnIMSuTPet7sM8bZ1byO/mLseCBLZ6H1QzIz562oANG6yiIiIZD0lyJKQc7/7RpeJ2f/9t4TWWbxpb6rCEREREUkZJcjSL11vpN15e+pFm/YcnPfrBW+MwfyLueu44b5FvW63rrEtKfGJiIiI9JduFCL90hIZ2m1PQ+vBrhVfCy/cA2iI9Cn+/pOrAXhxfQ1vOWpM3O1WH2hh5NDCZIcrIiIikjC1IEu/PLl8Z69l5qzY9aZ5H7zzZWr7cBvu1Tvr1a9ZRERE0kotyNIvCzfu6b1QWK7r8G2tHT3dWMQ50NLOiV9/kiNHD2VTbSPjRhTz8pcvPoxoRURERBKnBFn6pbvW4e50N7Zxnhnuzsod9cxZsZNVO+sPLnOHzz3wKgCbwnGVd+1PrMX5yp89z4ghhfz2n89KqLyIiIhId5QgS9pV17dw5ref7naZA39buftN8z9+79/528rdLP36JazcsZ+zpo0+ZHlHzHl92/6Dz9s6Yuza38zk8qH9ivFnT69lUnkJ/3Dq5H6tLyIiItlrQCTIZlYF1AMdQLu7z85sRJJK33xsRdxll/zouW7ndybNM2+eA8CKb1zK0KLg5XvTI8u558Wqg2V/vWATa3bVc99Lm3jlv97OqGFFQDCW85jhxQwpzO81xh88tQZACbKIiEgOGhAJcuhCd6/JdBCSei8n2H+5J39YtJW/rdzFceNLD0mOAb7659cPTu9vamPUsCLW7Ko/mHx/7JypbN7TwKSyEq4/dxq3PLGSb73rJMrDRLqr3768Gcf5pzOP7DGm9o4Y9720iQ+ddSRFBbr+VUREJFsNpARZJGFffyQYTu75tT1/p1q4cQ+X/vi5Q4alu3v+xoPTf1u5m237mvjLsp1U3XIFAJU3Pn7INr780DIArpo5kdIhhazbXU/l6GEU5B+aBN+/cDPfeGwFTW0dfOrCo/t/cCIiIpJRA6WZy4E5ZrbYzG7outDMbjCzRWa2qLq6OgPhSbb6woOvHZIcdxW97fZ//mHpm5LjqLf94Fk21zZy8Q+f4923vcjqnfXc/Ohy3J2Lf/gs33x8JQDrdh9g9/5mGlraWbChlsobH2dTbUPyDkpERJKqoaWdA5Gx+0XM3TMdA2Y20d23m9lY4Cng39y9286os2fP9kWLer8jWzL1lDTJ4FZaXEB9kk6al5wwjitOnsBVMyfS0NrB8OI3fsCJxZzwZoSs3X2AEUMKGT9ySFL2K8nz4voazpw6muXb61hUtZePnTu113X+7f5XmTWljOt7KNvWESPfjLy8rveofENHzGnriDGkMJ+2jhgFeXbwDpYy8DS3dbCnoZXVu+r57cub+fH7ZzGsWD/aZtKL62rYWNtwsLtcR8zZVNvAnc9v4P6FWwC47Z9O5axpo9le18TNj6xgYdUejh03nHs+egYTy0qAQ3OCj55TybCiAhZv2su9HzuDfY2tYDBiSGGP17vUHGjhQHM7lWOGdbvc3Yk55PdwTsgGHTHnpfW1nHvMGFbt3M/EshIaWtppaGnn6LGlaY/HzBYnep3bgEiQo8zsJuCAu9/a3XIlyCKD3+cuPoaYw4HmdkYPL6Ij5sTc+cnTaxlgp6yU+Y9LjmVIYT7/Hf4y0ZsrTprA48t2pDgqSaYxw4up6cONk1JpUlnJIb+oZaN3nzqZ1bv2HzKikaTOv5w/jefW1rByR//q+3/+6VQuP2lCkqPqWVYlyGY2DMhz9/pw+ingG+7+RHfllSCLiIiIZLe5/3EBU+O0oKdKXxLkgfB7zzjgofCnwgLgt/GSY5FsV5BntMeCL6UfP3cqz62t5vKTJjC0KJ/HX9vB0q11zJw8ktIhhby6eS8NrR1ccfIEHn/t0JbBC6ZXcPmJE1i8aS9/XrKNgjzj9KmjWLChlua2GG85ajTlQ4sOtiiOHlbEeceM4elVu5k6ZhgVw4t5etWh402PLS1md33vrVmjhhUdvP13dPrLlx/H/82vYvzIITS0tNMecz7x1qNo64ixs66ZYcUFLNy4h2dW7eaqmRN5fm0108eXcuH0sVwyYzz7m9pYuHEPpx5ZzkmTRuI47mAG+Wbk5xlrdx/gh3PWMLm8hHefNhkzeH3bfl7fVsf5x45hc20jy7fv5w+Lt3Yb+/tnT2H5jjryzHhtax1TRpVwzSmT+cOiLeyoaz6k7LtPnczVsyZSUVrM8OICGls7cIL/3dIt+8jPy+PlDbXsb24DYH9TO9/5h5No64hRc6CVwnwj5hDz4DiGFxdQOqSAoUX5tHbEONDSzp4DrZQOKWRIYR5DiwsYW1pMQ0s7Le0xRpYU0tIeoyN8vTS3ddDY2s7EshIMw3Hy84ym1g46Ys6IkkI21TaybV8T+xpb2b6vmZMmjWRaxTBqDrRgGCNLCnl1y1627WuiubWD7XXN/DFOXfXkH06ZxCUzxjN/XQ2/WrCpz+t3lWcQcxg/YghjSot4fdt+hhcXJNQn9BcfPJXq+mbeOn0sjy3dzkXHj+Pynz7/pnJ3XTub7fua+Nkz69hd38LFx4+lpT3G6p31XDC9gh11zawLuzfNmDSCP72y7ZD1Lz5+7CFjtP/npdM5bnwpMYeZU0ayv6kNd2hpj7G3sZWykiLqW9qoqmmksbWdtxw1BrPgPTO8uID2mFOQZzS1dXDjg8tYsb2OfU1tNLZ2HNzH1bMm8vCS7XGP/S1HjebF9bWJVDGTy0uoKC3m1c37+OQFR1GUn8fIkkJGDy/iyNHDmDhyCM1tMeat2c33nliNuzNl1FBuOH8af6/aQ11TG0dXDGfn/qCeOmJOdX0L2+uaKR1SwL+cP401uw5wVMVwNtU2cP6xFSzYUMv08aU0tXXw82fWccP507j3xSouO3E8VTWNFBfmMWPiCPY2tnGguZ0rT57AM6t2c8oRZdQ3t7N9X/MhF1UDfO3KE/jJ02upawred5+56BiOqhjGtDHDOWHiCNpjMTbXNrK3sY3CfGPznkaOGz+CoUX5rNixn2FFBTS1ddDeEaP6QAtrdtUzYWQJVTUNzDqijMrRw1hffYCNNQ2s232AfY1tTB9fSkVpMVNHD+N/5q1jxqSRbNvbxDWnTIq8rwvIM9iyt5FJZUOZVF5Ca3uMUcOK2NfYSlFBHvsa22hobae4IJ/9TW2UFOUzrKiAjbUNlA8tZMSQQiDoVlF9oIVYzBk3YgiPLN1OaXEBtz27nrKSQmYdUcYn33o0NQ0tLN9Wx61z1sT9v3/gjCncv3ALpUMKeM9pkynKz2NyeQlVtY0s3LiH848dw7HjSrlt3npOOaKM+xdu4fTKcj55wVFsrm3kpkdXkGdw0qSRfPjsSo4bX0p9czsVpUVMLh+KOzhOW4fT2Rtk1/4WOmLO8CEFbA9/lRg9rIi1uw9QXJDHqp31TB9XmvbkuK8y3oLcV2pBFoAjRg3l4uPHcey44ZwxdRT7mto4euxwSiN9/NQ/s+86zwcvrKvBHWZOLqOuqY0NNQe4YPrYDEcn6XCgpZ18M9piMVraggtcK0qL+7ydvQ2tjCwppLUjxuY9jRw7rpSW9g6mfzVo/6i65Qo21zZSkG8H+3YORMu311FckJeR/pIATa0dFBfkxe2fvqH6AGNHDKEw3zBs0A4xWdfUxsiSwkyHkVWaWjsoKep93P9cklVdLPpKCfLAddz4Ut5x4gRe3ljL244by9tPGEdzW4xJ5SWHXJAmIrmrs1VY5wQRSbds62IhaXDfx86gtT3Gx+9bxNeuPIGd+5s5YcIIag60cPyEEZw8eST14QVRxQWH+43zmKTELCKDjxJjEckGOlNlgfecNvlgP8FvX3MSHe5MGzOM4yeMoKW9g7KSIoYU5iXUpaDzZhjdKR2in69ERERElCAnyVUzJ9LaHqOxrYPS4gLeM3syCzfuYVxpMVW1jUwsG8K7Zk0i5jBuRDHu9DjmaVe3vndmCqMXERERkU5KkBPQU6trTy7s4aImXT8mIiIiMjANzstdRURERET6SQmyiIiIiEiEEmQRERERkQglyCIiIiIiEUqQRUREREQisu5OemZWDWzKwK7HADUZ2G82Ul0lRvWUONVV4lRXiVE9JU51lTjVVeIyUVdHuntFIgWzLkHOFDNblOjtCXOd6ioxqqfEqa4Sp7pKjOopcaqrxKmuEjfQ60pdLEREREREIpQgi4iIiIhEKEFO3B2ZDiCLqK4So3pKnOoqcaqrxKieEqe6SpzqKnEDuq7UB1lEREREJEItyCIiIiIiEUqQRUREREQilCD3wswuM7PVZrbOzG7MdDyZYGZTzGyuma00s+Vm9tlw/k1mts3MloSPyyPrfCmss9Vmdmlk/qCvTzOrMrNlYZ0sCueNMrOnzGxt+Lc8nG9m9tOwPl4zs1Mj27k2LL/WzK7N1PGkgplNj7xulpjZfjP7nF5TATO728x2m9nrkXlJew2Z2Wnha3RduK6l9wiTJ05dfd/MVoX18ZCZlYXzK82sKfL6uj2yTrd1Eq/es1Gcukrae87MpprZy2Fd/c7MitJ3dMkTp55+F6mjKjNbEs7P9ddUvPwg+89X7q5HnAeQD6wHpgFFwFLghEzHlYF6mACcGk6XAmuAE4CbgP/opvwJYV0VA1PDOszPlfoEqoAxXeZ9D7gxnL4R+G44fTnwV8CAs4CXw/mjgA3h3/JwujzTx5ai+soHdgJH6jV18HjPB04FXk/FawhYCJwdrvNX4B2ZPuYk19UlQEE4/d1IXVVGy3XZTrd1Eq/es/ERp66S9p4Dfg/8Yzh9O/DJTB9zsuqpy/IfAF/Ta6rH/CDrz1dqQe7ZGcA6d9/g7q3AA8DVGY4p7dx9h7u/Ek7XAyuBST2scjXwgLu3uPtGYB1BXeZyfV4N3BtO3wu8KzL/Pg8sAMrMbAJwKfCUu+9x973AU8Bl6Q46TS4C1rt7T3fIzKnXlLs/B+zpMjspr6Fw2Qh3f8mDT5/7ItvKOt3VlbvPcff28OkCYHJP2+ilTuLVe9aJ87qKp0/vubBV723AH8P1s7aueqqn8DjfB9zf0zZy6DUVLz/I+vOVEuSeTQK2RJ5vpefEcNAzs0rgFODlcNanw59J7o78TBSv3nKlPh2YY2aLzeyGcN44d98BwQkFGBvOz/W6AvhHDv2w0Wuqe8l6DU0Kp7telKo+AAAgAElEQVTOH6w+RtDq1Gmqmb1qZs+a2XnhvJ7qJF69DybJeM+NBvZFvpgM1tfVecAud18bmafXFG/KD7L+fKUEuWfd9XPJ2XHxzGw48CDwOXffD9wGHAXMAnYQ/OwE8estV+rzHHc/FXgH8CkzO7+HsjldV2EfxauAP4Sz9Jrqu77WTc7UmZl9BWgHfhPO2gEc4e6nAP8O/NbMRpBDddKNZL3ncqUOP8ChX+j1mqLb/CBu0W7mDcjXlRLknm0FpkSeTwa2ZyiWjDKzQoIX/2/c/U8A7r7L3TvcPQbcSfDTG8Svt5yoT3ffHv7dDTxEUC+7wp+KOn962x0Wz+m6IvgS8Yq77wK9pnqRrNfQVg7tcjAo6yy8yOdK4J/Cn2YJuwvUhtOLCfrSHkvPdRKv3geFJL7nagh+Li/oMn/QCI/tH4Dfdc7Ta6r7/IBBcL5SgtyzvwPHhFfmFhH8FPxIhmNKu7DP1V3ASnf/YWT+hEixa4DOK34fAf7RzIrNbCpwDEEn+0Ffn2Y2zMxKO6cJLhZ6neA4O6/KvRZ4OJx+BPhIeGXvWUBd+HPUk8AlZlYe/uR5SThvsDmkNUavqR4l5TUULqs3s7PC9/ZHItsaFMzsMuCLwFXu3hiZX2Fm+eH0NILX0YZe6iRevQ8KyXrPhV9C5gLvCdcfdHUFXAyscveDP/nn+msqXn7AYDhfHe5VfoP9QXDF5RqCb4VfyXQ8GaqDcwl+0ngNWBI+Lgd+BSwL5z8CTIis85WwzlYTueJ0sNcnwZXdS8PH8s5jJOif9zSwNvw7KpxvwC/C+lgGzI5s62MEF8asAz6a6WNLQV0NBWqBkZF5ek0Fx3Q/wU+3bQQtKNcn8zUEzCZIhNYDPye8q2o2PuLU1TqC/oyd56vbw7LvDt+XS4FXgHf2Vifx6j0bH3HqKmnvufD8tzCs/z8AxZk+5mTVUzj/HuATXcrm+msqXn6Q9ecr3WpaRERERCRCXSxERERERCKUIIuIiIiIRChBFhERERGJUIIsIiIiIhKhBFlEREREJEIJsohIljKzz5nZ0EzHISIy2GiYNxGRLGVmVQTjiNZkOhYRkcGkoPciIiKSaeGdGX9PcKvVfIIbMUwE5ppZjbtfaGaXADcDxQSD6n/U3Q+EifTvgAvDzX3Q3del+xhERLKFuliIiGSHy4Dt7j7T3U8EfgxsBy4Mk+MxwFeBi939VGAR8O+R9fe7+xkEd6L6cZpjFxHJKkqQRUSywzLgYjP7rpmd5+51XZafBZwAzDezJcC1wJGR5fdH/p6d8mhFRLKYuliIiGQBd19jZqcBlwPfMbM5XYoY8JS7fyDeJuJMi4hIF2pBFhHJAmY2EWh0918DtwKnAvVAaVhkAXCOmR0dlh9qZsdGNvH+yN+X0hO1iEh2UguyiEh2OAn4vpnFgDbgkwRdJf5qZjvCfsjXAfebWXG4zleBNeF0sZm9TNAwEq+VWURE0DBvIiKDnoaDExHpG3WxEBERERGJUAuyiIiIiEiEWpBFRERERCKUIIuIiIiIRChBFhERERGJUIIsIiIiIhKhBFlEREREJEIJsoiIiIhIhBJkEREREZEIJcgiIiIiIhFKkEVEREREIpQgi4iIiIhEKEEWEREZwMzsAjPbmsTtVZqZm1lBsrbZh30fYWYHzCy/n+t/2cx+mey4RLpSgixJYWZVZtZqZmO6zF8Snogr0xzPBWYWC0/E9Wa22sw+ms4Y+iusr6MzHYeIvMHMzjWzF82szsz2mNl8Mzs9XHadmb2Q6RiTJTyfN4Xnzn3hcX/CzA47Z3D3ze4+3N07EojjTV8M3P3b7v7xw42jm31dZ2Yd4WfG/vCz68pk7yfZMvllZ7BTgizJtBH4QOcTMzsJKMlcOGx39+HACOCLwJ1mdkJfN9Lflo5M0ElSJPnMbATwGPAzYBQwCbgZaMlkXIk4jHPCO929FDgSuIXgHHpX0gIbmF4KPzPKCI7192Y2qi8bsEDW5Fb6zIgva/6JkhV+BXwk8vxa4L5oATMrNrNbzWyzme0ys9vNrCRcVm5mj5lZtZntDacnR9adZ2bfDFtu6s1sTtcW6+544M/AXuCEcFt/MLOdYWvQc2Y2I7Kfe8zsNjP7i5k1ABea2RVm9mrYsrDFzG6KlO/8Bv/RcNnesLXldDN7LWyB+XmXeviYma0Myz5pZkeG858LiywNWzLeH86/MmzR6GzNOTmyrSoz+6KZvQY06IQnknTHArj7/e7e4e5N7j7H3V8zs+OB24Gzw/fsPoAEzxnXhufCGjP7SmR5SXge2mtmK4DTo8GY2Y1mtj48D64ws2siy64Lz5E/MrM9wE1mlh+ed2vMbANwRaIH7u517v4I8H7gWjM7MdxPT+fyldHWVzMrCPd9atcWz/C8uTI8lg1m9i/h/GHAX4GJYb0eMLOJZnaTmf06su2rzGx5eG6cF/4/OpdVmdl/hOfhOjP7nZkNSeCYY8DdBA080xL8bPqWmc0HGsN1uj2usPwFZrbVzL5gZrvNbIeZvcvMLjezNRb8QvHlSPm8yP+81syiiXvnZ8a+sI7ODtfp9jMmXOZm9ikzWwus7a0+cpa766HHYT+AKuBiYDVwPJAPbCFofXCgMiz3Y+ARglaYUuBR4DvhstHAu4Gh4bI/AH+O7GMesJ7gw6okfH5LnHguALaG03nANUAbMD2c97FwH8VhTEsi694D1AHnhOsOCbd3Uvj8ZGAX8K6wfGV4jLeHZS8BmoE/A2MJWpt2A28Ny78LWBfWUwHwVeDFyP4dODry/NRw/TPDer02rO/iSN0vAaYAJZl+Leihx2B7EPwKVQvcC7wDKO+y/DrghS7zEjln3Bmey2YStEYfHy6/BXg+PE9OAV7vPJ+Fy98LTAy3/X6gAZgQiaUd+Lfw/FICfAJYFW5rFDA33H9BnOOtAi7uZv5m4JPhdE/n8q8Bv4msdwWwqsuxF0SWHQUY8FaCBPPUSB1u7RLDTcCvw+ljw2N/O1AIfCE8txZFjmNhWFejgJXAJ+Ic88H/YVhvnwXqgZEk9tm0GZgRrluYwHG1h/VUCPwzUA38Ntz+DILPkGlh+c8BC4DJBJ9Z/wvc3119hvMS+Yx5KqwTfWbEe99nOgA9BseDNxLkrwLfAS4L34AF4ZuxMjxRNABHRdY7G9gYZ5uzgL2R5/OAr0ae/yvwRJx1LwBiwD5gD0EC+Y9xypaFMY4Mn98D3NfL8f4Y+FE43XmCmhRZXgu8P/L8QeBz4fRfgesjy/LCk+eR4fOuCfJtwDe77H81byTcVcDHMv0a0EOPwfwIk417gK1hcvMIMC5cdh1dEuRu1u/unDE5snxh5zkK2ABcFll2A10SxS7bXgJcHYllc5flzxBJDAm+xPcnQV4AfKW3czlwNEFyOTR8/hvga12OPd6+/wx8Npy+oOtxc2iC/F/A7yPL8oBtwAWR4/hQZPn3gNvj7Pe68P+6D6gJj/VNdRCW7e6z6Ru9/P+7HlcTkB8+Lw3r5MxI+cW88YVqJXBRZNkEggafgu7qk8Q+Y96W6ffUQH/op1hJtl8R/OQzlS7dK4AKgm/gi82sc54RtIpiZkOBHxEk1+Xh8lIzy/c3LujYGdleIzC8h1i2u/vkrjMt6FP8LYJWmAqCRBpgDEHLMQSt39F1ziRo1TkRKCL4Fv+HLpveFZlu6uZ5Z6xHAj8xsx9Ed0HQ0rypm+M4kuCnzX+LzCsiaBXptAURSRl3X0mQRGFmxwG/Jkh6P9Bd+QTPGfHOZxM59D19yHnBzD4C/DtBckS4XrS7WdfzQY/b64NJBA0OPZ7L3X2dma0E3mlmjwJXAad0t0EzewfwdYLW4Lxwu8sSjGcikWNx95iZbQnj7NS1jqPnza4WuPu53cSYyGdT18+M3o6rNrJuU/i3p8+Mh8wsFlneAYyLcxyJfMboM6MX6oMsSeXumwgu1rsc+FOXxTUEb/oZ7l4WPkZ6cFEEwOeB6QTfokcA54fzjeT6IHA1QYv3SN74kInux7us81uCFqMp7j6SoDtFf+PaAvxLpA7K3L3E3V/sofy3upQf6u739xCviKSIu68iaE0+sXNWN8UO55yxg6A7RKcjOifCvqR3Ap8GRrt7GUEXjJ7OX3G3lygLRuyYBLxA7+dygPsJvjxcDaxw93XdbLOY4Ne1Wwla48uAv0SOpbfz2naCZLBze0ZwnNv6eny9SOSz6WCsCRxXX20B3tHlM2CIu2+j+zpK5DNGnxm9UIIsqXA9wc83DdGZHlz4cCfwIzMbC2Bmk8zs0rBIKcFJd194AcLXUxRfKUF/v1qCb/XfTnCdPe7ebGZnECTZ/XU78CULLww0s5Fm9t7I8l3AtMjzO4FPmNmZFhhmwQVApYcRg4gkyMyOM7PPd16YZWZTCJK/BWGRXcBkMyuKrHY454zfE5wjysN9Rn89GkaQ3FSHsXyUNxL1nrb3GTObbGblwI2JBmJmI8IL7h4g6NqwLIFzOWH5S4BPEnxZ6E5ny3o10B62ul4SWb4LGG1mI3s4rivM7CIzKyRIZFuAeI0N/dXXz6bejquvbge+ZW9czF1hZleHy6oJfgWd1qV8T58xkgAlyJJ07r7e3RfFWfxFgosHFpjZfuBvBN/MIfi5soQ3+n89kaIQ7yP4mWkbsII3PuR68q/AN8ysnuDCit/3d+fu/hDwXeCBsA5eJ7jwp9NNwL3hVdnvC+vyn4GfE4zEsY7wp14RSYt6gotkX7ZgZJsFBO/bz4fLnwGWAzvNrCacdzjnjJsJzlEbgTkEXdcAcPcVwA+AlwgSyJOA+b1s707gSWAp8Apv/nWvO4+GsW8h6Hf8QyA6lnxP53LcfUcY41uA33W3A3evBz5DUDd7Cb5EPBJZvoqgJXpDeD6c2GX91cCHCIbfqwHeSTA8XWsCx9cXffps6u24+uEn4fpzwv/JAoLXI+7eSNBlcH5YR2cl8BkjCTB3tbKLiIiIiHRSC7KIiIiISIQSZBERERGRCCXIIiIiIiIRSpBFRERERCKy7kYhY8aM8crKykyHISLSJ4sXL65x94pMx5EsOheLSLbpy3k46xLkyspKFi2KN4KYiMjAZGb9vXvZgKRzsYhkm76ch9XFQkREREQkQgmyiIiIiEhEShNkM6sys2VmtsTM3vRbXHjb3J+a2Toze83MTk1lPCIiucbMppjZXDNbaWbLzeyz3ZTRuVhEJCIdfZAvdPeaOMveARwTPs4Ebgv/iohIcrQDn3f3V8ysFFhsZk+FtyzupHOxiEhEprtYXA3c54EFQJmZTUjmDjpizgfvXMBt89Ync7MiIlnB3Xe4+yvhdD2wEpjUpVjKz8Wvbt7L5T95nuXb65K5WRGRlEh1guzAHDNbbGY3dLN8ErAl8nwrbz5xY2Y3mNkiM1tUXV3dpwDy84x9jW3MXxevEVtEJDeYWSVwCvByl0UpPxfnmbFix36272vu03oiIpmQ6gT5HHc/leDnu0+Z2fldlls36/ibZrjf4e6z3X12RUXfhxE9vbKcVzbvpb0j1ud1RUQGAzMbDjwIfM7d93dd3M0qST0Xjx1RDEB1fUuf1hMRyYSUJsjuvj38uxt4CDijS5GtwJTI88nA9mTHMbtyFI2tHazY0fUzQURk8DOzQoLk+Dfu/qduiqT8XDx6WJAg765XC7KIDHwpS5DNbFh4QQhmNgy4BHi9S7FHgI+EV1CfBdS5+45kxzK7shyARVV7k71pEZEBzcwMuAtY6e4/jFMs5efiooI8xpYWs21vUzI3KyKSEqkcxWIc8FBwbqYA+K27P2FmnwBw99uBvwCXA+uARuCjqQhkwsgSJpWVsHDjHj527tRU7EJEZKA6B/gwsMzMloTzvgwcAek9F1eOHkZVbUMqNi0iklQpS5DdfQMws5v5t0emHfhUqmKIOvfoMfz19R20d8QoyM/04B0iIunh7i/QfR/jaJm0nIunjBrKi+t1wbSIDHw5kymed+wY9je389o2DTEkIpIJk8qGsGt/M226YFpEBricSZDPOWoMZvDCWrVeiIhkwqTyEmIOO+t0oZ6IDGw5kyCXDyvipEkjeX5t38buFBGR5JhUNhSALXsbMxyJiEjPciZBBjjvmDG8snkfdY1tmQ5FRCTnHDEqSJC37tFIFiIysOVUgvy248bSEXPm6yIREZG0m1A2hPw8Y/MetSCLyMCWUwnyzMllDC8u4AXddlpEJO0K8/OYWDZEXSxEZMDLqQS5ID+Ps6aNZr4SZBGRjDhi1FC1IIvIgJdTCTLAuUePZlNtI1t0ghYRSbsjRg1lc63OvyIysOVegnzMGAC1IouIZMCUUUOpbWiloaU906GIiMSVcwnyURXDGTeimHmrNdybiEi6dY5koW4WIjKQ5VyCbGZcNXMic1bspKqmIdPhiIjklCNHDQNgk7pZiMgAlnMJMsD1507DgUeXbs90KCIiOeWI0Z0tyGqgEJGBKycT5PEjh3DaEeU89tqOTIciIpJTRpYUUja0kCq1IIvIAJaTCTLA1bMmsnpXPSt37M90KCIiOeX48SNYXLU302GIiMSVswnyFSdPJD/PeETdLERE0uotR41m9a56jWQhIgNWzibIo4YVce7RY3j41W20d8QyHY6ISEqY2d1mttvMXo+zfKSZPWpmS81suZl9NNUxHTOuFID11QdSvSsRkX7J2QQZ4ANnHMH2umaeXrU706GIiKTKPcBlPSz/FLDC3WcCFwA/MLOiVAZ0zLjhAKzbrQRZRAamnE6QLz5+LCOGFPDMSiXIIjI4uftzwJ6eigClZmbA8LBsSvs+TC4vwUxDvYnIwJXTCXJBfh5nThvNC+tqcPdMhyMikgk/B44HtgPLgM+6e7f9zszsBjNbZGaLqqv7f7Ol4oJ8Jo4sYYtuFiIiA1ROJ8gAbz22gm37mtQXTkRy1aXAEmAiMAv4uZmN6K6gu9/h7rPdfXZFRcVh7XTKqBI2KUEWkQFKCfKxwUn+2TU1GY5ERCQjPgr8yQPrgI3Acane6TFjS1m1Y78ukhaRASnnE+Qpo4ZyVMUwnl3T/58LRUSy2GbgIgAzGwdMBzakeqdnThtFQ2sHy7bVpXpXIiJ9lvMJMsBbjx3Lgg21NLZqTE4RGVzM7H7gJWC6mW01s+vN7BNm9omwyDeBt5jZMuBp4IvunvKf1M6aNhozeEajCInIAFSQ6QAGgrcdN5a752/kxXW1XHzCuEyHIyKSNO7+gV6WbwcuSVM4B40ZXszJk8tYvEl31BORgUctyMDpU8sZWpTP3NVqyRARSZepo4dqqDcRGZCUIBMMOXTeMWOYu2q3hnsTEUmTqWOGs72uifrmtkyHIiJyCCXIofOPrWB7XTMbaxoyHYqISE6YXVmOOyzY0NN9TERE0k8Jcuico8YAMH+dhnsTEUmH0ytHUVSQx0vrazMdiojIIZQgh44cPZRJZSXMX6cTtYhIOhQV5HH+MRX87u+bdbMmERlQlCCHzIxzjh7NSxtq6YipH7KISDrcfPUM2mPOL5/fmOlQREQOUoIccc7RY6hrauN1DVwvIpIWk8pKuOj4sczTKEIiMoAoQY54S9gP+QX1QxYRSZtTppSzo66ZmgMtmQ5FRARQgnyIitJiZkwcwVMrdmU6FBGRnDFj0ggAlm/fn+FIREQCKU+QzSzfzF41s8e6WXadmVWb2ZLw8fFUx9Oba06ZxJIt+1i7qz7ToYiI5IQZE0YCqHubiAwY6WhB/iywsoflv3P3WeHjl2mIp0dXzZpInsGjS7dnOhQRkZwwcmghx44bzoINGkVIRAaGlCbIZjYZuALIeOKbqLGlQ5hdOYo56mYhIpI2Z08bzd+r9uiueiIyIKS6BfnHwBeAWA9l3m1mr5nZH81sSorjScglJ4xj1c56Ntc2ZjoUEZGc8PYTxtPcFuOVzfsyHYqISOoSZDO7Etjt7ot7KPYoUOnuJwN/A+6Ns60bzGyRmS2qrq5OQbSHunTGeAAeW6ZuFiIi6TBjYnCh3pqduv5DRDIvlS3I5wBXmVkV8ADwNjP7dbSAu9e6e+e4PncCp3W3IXe/w91nu/vsioqKFIYcmDJqKKcdWc6jS3ekfF8iIgLlw4oYW1rMal0gLSIDQMoSZHf/krtPdvdK4B+BZ9z9Q9EyZjYh8vQqer6YL60uOn4sK3fsp7pe43KKSPYys7vNbLeZvd5DmQvCkYSWm9mz6Ywvavr4UtYoQRaRASDt4yCb2TfM7Krw6WfCE/JS4DPAdemOJ57zjg5aqueu0t2dRCSr3QNcFm+hmZUB/wNc5e4zgPemKa43mT6ulNU762lp78hUCCIiQJoSZHef5+5XhtNfc/dHwukvufsMd5/p7he6+6p0xJOIEyeNYNqYYTz06rZMhyIi0m/u/hywp4ciHwT+5O6bw/IZaxU47chyWtpjumGIiGSc7qQXh5lx2YnjWVi1h7omDTskIoPWsUC5mc0zs8Vm9pF4BVN9wfSsI8oAWLpFI1mISGYpQe7BRcePoyPmPLcm9SNniIhkSAHBBdJXAJcC/2Vmx3ZXMNUXTE8YWcK4EcUsUYIsIhmmBLkHs6aUUTa0kGfUD1lEBq+twBPu3uDuNcBzwMxMBTNrSplakEUk45Qg9yA/z7j0hPHMWb6Txtb2TIcjIpIKDwPnmVmBmQ0FziSDIwrNmlJOVW0jextaMxWCiIgS5N5cPWsiDa0d6mYhIlnJzO4HXgKmm9lWM7vezD5hZp8AcPeVwBPAa8BC4JfuHndIuFSbOWUkAEu2qhVZRDKnINMBDHRnTB1F+dBC/rJsJ5edOKH3FUREBhB3/0ACZb4PfD8N4fTq5Mll5Bks2byPC6ePzXQ4IpKj1ILci4L8PC47cQJ/W7mLplaNzSkikkrDiws4ZmypLtQTkYxSgpyAd548gcbWDuau1sV6IiKpdvyEUtbqjnoikkFKkBNw5rTRjBlezGOvbc90KCIig97RY4ezva6ZhhZdHC0imaEEOQH5ecblJ43n6ZW72deoK6tFRFLpqIrhAGyobshwJCKSq5QgJ+h9s6fQ0h7j8WU7Mh2KiMigdvTYIEFeV61uFiKSGUqQEzRj4gimjhnGX5ftzHQoIiKD2pGjh5GfZ6zfrRZkEckMJcgJMjMuO3E8CzbUqpuFiEgKFRXkMXXMMJZtq8t0KCKSo5Qg98GlM8bTHnOeXqnRLEREUunco8ewYEMtB3ShnohkgBLkPjh50kgmjBzCE8vVzUJEJJUuP2kCLe0x5ml4TRHJACXIfZCXZ1w6YzzPrammsVWtGiIiqXLakeWMHlbEUyt2ZToUEclBSpD76NIZ42lpj/Hs6upMhyIiMmjl5xmzppSxeqdGshCR9FOC3EenV5ZTPrRQ3SxERFKscswwqmobiMU806GISI5RgtxHBfl5XHT8OOatrqa9I5bpcEREBq1pFcNobouxbV9TpkMRkRyjBLkfLpheQV1TG0u3aggiEZFUmTFxJADLt+tcKyLplXCCbGblZnaGmZ3f+UhlYAPZuUePIc/g2TXqhywi6WNmQ8zs383sT2b2oJn9PzMbkum4UuW48aUU5BmvqTFCRNIsoQTZzD4OPAc8Cdwc/r0pdWENbGVDizj1iHL+umwH7uobJyJpcx8wA/gZ8HPgeOBXPa1gZneb2W4ze72XcqebWYeZvSdp0R6mIYX5HDOulOXb92c6FBHJMYm2IH8WOB3Y5O4XAqcAOd18+p7TJrN29wGWbNmX6VBEJHdMd/fr3X1u+LgBOLaXde4BLuupgJnlA98laPwYUGZNGcmiqj3UNbVlOhQRySGJJsjN7t4MYGbF7r4KmJ66sAa+K06ewJDCPP64eGumQxGR3PGqmZ3V+cTMzgTm97SCuz8H7Ollu/8GPAgMuLtyvOe0KTS0duiGISKSVokmyFvNrAz4M/CUmT0MbE9dWANf6ZBC3nHiBB5Zup3mto5MhyMig5iZLTOz14AzgRfNrMrMNgIvAYd1PYiZTQKuAW5PoOwNZrbIzBZVV6fnR8RZU8ooLS7g71W95fgiIslTkEghd78mnLzJzOYCI4EnUhZVlnjPaZN56NVtPLViF++cOTHT4YjI4HVlCrf9Y+CL7t5hZj0WdPc7gDsAZs+enZYLMPLzjOMnjFA/ZBFJq4QS5Ch3fzYVgWSjs6eNZlJZCfcv3KwEWURSxt03pXDzs4EHwuR4DHC5mbW7+59TuM8+OWPqKG57dj3b9zUxsawk0+GISA7QOMiHIS/P+PDZR/Li+lqWaRgiEclC7j7V3SvdvRL4I/CvAyk5Bnj/6VPIzzO+/+TqTIciIjlCCfJh+uCZR1BaXMBtz67LdCgiIm/y/9u77/CoyvT/4+97Jj2ENBIgBQiELj2ACixSFEQBGyr2FQuu+nPX3a9l13Utq66iq2tbFhvK2l0rYgMUkB56JwECJARSIKRA2szz+2OGMCEFAjOZSbhf1zVXJmfOzPnMmZln7nnOc84RkQ9xjFXuKiKZIjJFRKaKyFRvZztViVEh3DC4PbPX7+NAYam34yilzgJaIJ+hlkH+XH9ue77buJ/d+SXejqOUUtUYYyYbY9oaY/yNMQnGmLeMMdONMTV2yjPG3GKM+cwbOU/mpvPaU2EzzFyS4e0oSqmzgBbIbjB5UCJWEa7891IOH9FjdSqllLt1aBXKiK4xvLskg/JKu7fjKKWaOS2Q3aB9dCgzbhpAXnEZ/1utx0VWSilPuKxfPEfKbaTnFHs7ilKqmdMC2U1GdmtN38QIPlixR08/rZRSHtCjbUsANmbpTtFKKc/SAtmNrhvcjvScYlZmHPJ2FKWUanY6xbQgMSqYD1dqR4RSyrM8XiCLiFVE1ojI7FpuCxSRj0UkXUSWi0gHT+fxpPG94wgL8uPdpRnejqKUUs2OxSLcMawja/YUsHZvgbfjKKWascboQb4P2FLHbVOAQ8aYZOBF4NlGyOMxwQFWbji3PXM2ZLMzV8fIKaWUu03sF0+A1aIYHCQAACAASURBVMI367K9HUUp1Yx5tEAWkQTgEuDNOmaZCLzrvP4ZMEpOdq5TH3frkCQCrBZe+3mHt6MopVSz0zLInwu6xjB7/T5sdh1moZTyDE/3IL8EPADUdUyeeGAvgDGmEjgMRJ84k4jcISKpIpKam5vrqaxuERMWyE3nteeLNZlkFRz1dhyllGp2JvSNI6eojBW7Dno7ilKqmfJYgSwilwI5xphV9c1Wy7QaXQLGmBnGmBRjTEpMTIzbMnrKdYPbYzfww8b93o6ilFLNzqhurQkJsPL1un3ejqKUaqY82YM8BJggIhnAR8BIEfnvCfNkAokAIuIHhANNvksgqVUo3dqE8ZU23kop5XbBAVYu7NGa7zZmU2HTk4YopdzPYwWyMeZh52lNOwDXAvONMTecMNvXwM3O61c552kWg8quGZjIur0FLEnP83YUpZRqdsb3jqPgSAWLtY1VSnlAox8HWUSeEJEJzn/fAqJFJB24H3iosfN4ynWD2xETFshrv6R7O4pSSjU7Qzu3ws8iLNvZ5Dc6KqV8UKMUyMaYX4wxlzqvP2qM+dp5vdQYM8kYk2yMGWSM2dkYeRpDoJ+V24clsTg9n9V79MQhSinlTkH+VgZ3jOKrtVmUV+owC6WUe+mZ9Dzo+sHtiQoN4B/fbdWzPimllJvdNqwj2YdLmbVst7ejKKWaGS2QPSg00I/7L+zCil0H+XlbjrfjKKVUs3JBlxj6JITz/A/byCsu83YcpVQzogWyh10zMJG24UHMWqo9HEqpxicib4tIjohsrOP260VkvfOyRET6NHbG0yUivHB1X8oqbUx8dTFFpRXejqSUaia0QPYwf6uFCX3jWJiWpz0cSilvmAmMref2XcBwY0xv4ElgRmOEcpfk2BY8d1UfsgqO8vK8NG/HUUo1E1ogN4JJAxIR4JEvNupYZKVUozLGLKSe48sbY5YYY47tSbwMSGiUYG501YAEJg9qxxuLdrFkhx72TSl15rRAbgTJsS14cGw3vt+0n09S93o7jlJK1WUK8F1dN4rIHSKSKiKpubm5jRjr5P42vgftokL4yxcbKa2weTuOUqqJ0wK5kUwZmsTADpE8/+N2jpRXejuOUkpVIyIjcBTID9Y1jzFmhjEmxRiTEhMT03jhTkGQv5WnLj+HXXkleuQgpdQZ0wK5kVgswoNju5FbVMa/f9nh7ThKKVVFRHoDbwITjTH53s5zuoZ1juHagYnMXJLB3C165CCl1OnTArkRpXSIYmLfOKYv2MHW/YXejqOUUohIO+Bz4EZjzHZv5zlTT152DnHhQTz8+Xr2Hjzi7ThKqSZKC+RG9rfxPQkN9OOpb7d4O4pS6iwgIh8CS4GuIpIpIlNEZKqITHXO8igQDbwuImtFJNVrYd3A32rhvSmDKCqt5Ma3lrMyQ09FrZRqOC2QG1lUaAB3/qYTi9Ly2JKtvchKKc8yxkw2xrQ1xvgbYxKMMW8ZY6YbY6Y7b7/NGBNpjOnrvKR4O/OZSo4N462bB5JXXM7UWauw2XU8slKqYbRA9oJrByYS5G/hjYU7vR1FKaWapaGdW/HU5eeQX1LO2r0F3o6jlGpitED2gsjQAG4Y3J4v12ZpL7JSSnnIiG6xBPpZeGfxLm9HUUo1MVoge8k9I5MJD/ZnysyVLErzreOJKqVUc9AyyJ87h3di9vps3vpVi2Sl1KnTAtlLIkICmDVlMHYDt85cyc7cYm9HUkqpZufekcmM6BrDc99v1c4IpdQp0wLZi86JD+ebe4ciCG8s0vHISinlbv5WC09f0YvwYH9uensFqXpUC6XUKdAC2ctiwgK5vF88H67Yy/ytB7wdRymlmp224cHM/9MFtAj048W52/WoFkqpk9IC2Qc8cVlPYsMCmfbDdsor7d6Oo5RSzU6LQD/uG9WZxen5/O79VXoqaqVUvbRA9gGBflaevOwctmQX8vuP13CkvNLbkZRSqtm5dUgStw9L4odNB/h2Q7a34yilfJgWyD5iTM823DsymTkb9nPnrFWUVdq8HUkppZoVi0V4cGw3eieE88Bn68k8pKeiVkrVTgtkH/LHi7ry3JW9WZSWx+PfbPZ2HKWUanb8rBZev74/ANf8Zxm780u8nEgp5Yu0QPYxVw9M5LahSXywfA8rdune1kop5W4JkSE8MKYrWQVHGT7tF6Yv2MGspRl8tTaLnMJSb8dTSvkAP28HUDXdf1EXvlq3jwc+W8enU88nJizQ25GUUqpZuWVIEvGRIdz9/mr+8d3WqukBVgsbHr+IQD+rF9MppbxNe5B9UEiAH69O7se+glImTV+iO+0ppZQHXNijNZufGMOsKYN455aBXNY3jnKbnV/T8rwdTSnlZVog+6jBHaP5z40DyMg/wvM/bPd2HKWUapb8rBaGdY5hRLdYnruqDy2D/Ph8TZa3YymlvEwLZB82olss1w1ux7tLM1i3t8DbcZRSqlkL8LMweVA7vl2fzWs/p3s7jlLKi7RA9nEPjOlKqxYB/HbmSg4fqfB2HKVUEyMib4tIjohsrON2EZGXRSRdRNaLSP/GzuhLHhzbjbE92zDth20sTs/TE4oodZbSAtnHRYQE8MTEczhYUs6d/02ltEKPj6yUapCZwNh6br8Y6Oy83AH8uxEy+SyLRXh8Yk9CA6xc/+Zyuj/6Pa//or3JSp1ttEBuAsb0bMP/jenKsp0H+fPnG6i06emolVKnxhizEKjvmJETgfeMwzIgQkTaNk4639S6ZRCf/24I945MxmY3/PPH7Tw9ZwuPfb2JVbsPUlhaoT3LSjVzepi3JuLuEcnkF5fz9uJdzN6QzQNjujJlaBIi4u1oSqmmLR7Y6/J/pnNajXMxi8gdOHqZadeuXaOE85aubcLo2qYrl/Ruy6TpS5mxcCcAM5dkADCscysevrg7PeJaejGlUspTtAe5CXl0fA+mXdWbQD8Lf/92C9+sr/H9pZRSDVXbr+xau0eNMTOMMSnGmJSYmBgPx/IN3dq0ZNUjF5L+1MWs+PMorugXD8DaPQVc+soiftp8wMsJlVKeoAVyEzMpJZG1j15EtzZh/PXLjWQfPurtSEqppi0TSHT5PwHY56UsPinAz4Kf1UJsyyD+eU1fMv5xCQseGEG7qBBufy+VV+aleTuiUsrNtEBugqwW4d83DKC0wsaTszfrWDil1Jn4GrjJeTSLc4HDxhjdPHUSUaEB/H50FwBe+Gk7N729guIyPamTUs2FxwpkEQkSkRUisk5ENonI47XMc4uI5IrIWuflNk/laW6SWoVyy5AOzNmwn+kLdno7jlLKR4nIh8BSoKuIZIrIFBGZKiJTnbPMAXYC6cAbwO+8FLXJuaxfPOsevYiU9pEs3J7LtTOWsigt19uxlFJu4Mmd9MqAkcaYYhHxB34Vke+ce0m7+tgYc48HczRbD47pxqasQl6au52rBiQQExbo7UhKKR9jjJl8ktsNcHcjxWl2wkP8+XTqeQx6eh4bswq58a0V/N+YrmzMOsyzV/WmZZC/tyMqpU6Dx3qQnYcMKnb+6++86FgAN7JYhMcm9KSs0s7Ap+YyY+EOb0dSSqmzjojw5d1DeHxCT7q3bcm0H7bx3cb93Pz2CnbkFp/8AZRSPsejY5BFxCoia4Ec4CdjzPJaZrvSefamz0QksZbbEZE7RCRVRFJzc3Xzlavk2BZM7BsHwNNztjJX96hWSqlGFx8RzM3nd+D5Sb0Z3b01AGv2FDDqhQVMnbWK7zfu93JCpU7dfxbs4IPleyg4Us57SzNYvjPf25EanTTGDl4iEgF8AdxrjNnoMj0aKDbGlDnHw11tjBlZ32OlpKSY1NRUzwZugopKK5jw6mJ25ZVw7cBEHrq4G/kl5bzw4zYiQwL47ZAkkmNbeDumUmctEVlljEnxdg530ba4fpU2Owu25zLl3ePr6LHxPQgN9OP85FbERwR7MZ3yhAqbHT+LNNnzE5RW2Phx8wGiQwO4/k1Hf6YIHCsT37llICO6xXox4ZlrSDvcKAUygIj8DSgxxjxfx+1W4KAxJry+x9FGuW4HS8qZ9sNWPlq5F9eXNcBqwd8qfHH3EJJaheJv1YOXKNXYtEA+O6UdKMLPauGO91JJyzk+3OLXB0eQEBnixWTKnXbkFjP6nwu4f3QX7h3V2dtxGqy80k6XR76r9bbxfeL4Zt0++iZG8OXdQ6qmV9jsvPZzOpv2FZJbVMaEPnHcOjTJozlX7T7IS3PTKKu0896tgwjytzbo/g1phz22k56IxAAVxpgCEQkGRgPPnjBPW5fDCU0Atngqz9kgKjSAZ67ozTUD23HZa4vpFR/Of24cAMDF/1rERS8uBODz351P/3aRblnm9gNFzNmQzfAuMfRz02MqpVRz0bl1GAAf33keT87eTFFpBXO35PDk7M1Mv2FAk+1tdHWopJzr3lzOtv2FvHhNXyb2ja8xT1mljaLSSkID/LAZQ4tAR/mxavdBcovKGN29NX5NtPPm5605/HbmSsBxyL/L+sWTGNU0fvyUVtjYur+IBz9bX236/Rd24YKuMeQWlTGqe2uSokN4eX465/ztBz676zy2Zhfx+4/XVrvP2r0F5BaX8X8XdcVicd/7ekduMZe9upi+7SJYlJZXNX3Zznwu6Oq5Hm2P9SCLSG/gXcCKY6zzJ8aYJ0TkCSDVGPO1iDyDozCuBA4Cdxljttb3uNprcWoqbPZqPcWvzEvjhZ+2A45NJm/dnMIFXWL5Zv0+Fm7P4+4RnYiLCG7Qr7GDJeUMe3Y+JeU2AK7oF8/fxvckPMQ9e22XVtg4Wm4jv6SM5FjHl0zBkXI+Sd3LpAGJRIYGuGU5p+qHTfsprbCREBlMv8RItzYAqvnTHmR1zOPfbOKdxRm8eE0fWrcM4uHPN7A7/wjjerXhgTHd6NAq1NsRAUeb17FVaFWRX5ulO/KZ/Mbxg1NFhviz+q8XsnRnPhU2w+PfbGJnbkmN+638y2henpfGrGW7q6bdeG57Hrm0O4F+DesV9KZv12dz9werAcc49KwCx8m7Lu8XzzNX9GpwD6erSpu96kdDfnEZ0S3ce6Sqw0crGD7tZwqOVAAQFuTHkodGsnZvAUOTW1X78ZZTWMqgp+fVeIzxfeIY2S2GikrDA/9zFNkX9mjNGzedeVNXWmHjidmb+WD5nmrTR3ePZd7WHF6Y1Icr+ic06DF9coiFu2ijfHpsdsOmfYcJCfBjzEsLsdlrf907RIdwRf8E7hmRXGcBePhIBbe9t5It2UUUl1Xy/m2DeXPRTn7e5tiB8uGLu3HtoHaEB1cvlI0xp9xbMnXWKr7fVH2nlkt7t2W28/Tal/RuyyvX9qvKWFphI9DPQqXdeGQIyYlfAh1jQplxY0qDx3UfKiln9vp9PPv9Nu4b1Znbf9OxwVlO/PFTl9IKG3bn5/uXbbmkZhyiV0JLEiND6Na2JS0C/ThQWMqq3YcYlBRFK5fGN6vgKEvS8+gR15KurcNYs7eAhz/fwItX96VXQr2joLymtMJ2Rl9Gp2JP/hGiWgRU9X41hBbI6hhjDOc+M48DhWW13j6iawz/mtzP7YeIs9sNO3KLqxW8O3KLmf7LDj5dlclFPVrz1OW9CLBaeH/Fbp77fhv+VuHHPwyndctAJk1fSl5xGcH+VoZ2bsXdI5K59OVfyS8pZ9KABIL8rdUK3tMRHuzPExN7snlfISXllTxySQ+Pf66LyyrZmHWYlPaRNXqxcwpLSc8tJjTAjw7RodU6gJ79fiv//mUH8RHBfH3PECJDAvhgxR4e+bJqVyuWPTyKNuFBp5xlS3Yhd85aRUSIP7vySrhjWEf2HjrCJ6mZ3DsymaRWobSPDuHdJbtJzymmVVggg5OiuHtEcrXH2Zh1mHWZBUwakEiA3/HndKS8ksFPz6N7m5YEBVhZuN3xvX3Due3487juhATU3bYdKilnZcZB7pi1itAAK1/fO5ROMce/Ayttdsa9vIjtB4p5fEJPxvVqy+Gj5Szdkc+N53U45XVwzCNfbuC/yxzF8eRBifRJiCA5tgUD2kdiN46TpjWUFsiqXqt2H2L5rnye+34bYYF+3H9RF95YuJN9h0urzdc3MYIhydEMTopmWOdWfLEmi4Xbc5mzYT/lNjutWwZyz4hkbjyvA8YYnvluKzMWHj9pyaCkKJKiQ7l6YCKrdx/iqTlbEIE/XdS1xoc57UAR8ZHBvLVoF/+cu71qDHV8RDDlNjt5xWVV0+LCg9h3uJRbzu/AlKFJPPDZepbuzCfY30qAn4V7RiQzZ2M2mYeOcm7HaIYmR1NuM3RtHcagpKh6101WwVEiQ/yrNRLllXau/s9SNmQdZlyvtnyzznEW3lYtAvn5T8MJDfCr9mPi8NEK/rtsN+v2FvDXS3uwdEc+LYL8uPfDNTV+mNw+LIk/j+uOiHC03EZwgJWyShuvzEtnzoZsIkMDSI5pQZ/ECK5OSeCB/61nzoZsfjskiQXbctmcXUhybAvuHtGJS3vH4W+1kJ5TzKer9vJrWh6b9hXiZxEqT1huSICVkAArecXl1aa9d+sgNmcX8uhXm+pcRxEh/hQcqWB4lxhuPLc9o3u0rnedNob3lmbw99lbGJIczYpdB3l5cj9GdW9NblEZqRkHGdq5FRYRDhSW0iE6tEG9//O3HmDaD9sZnBTFzCUZ9Gjbkm//39AGbxrXAlm5Wp9ZwPVvLqeotJIPbh/M+Z1asWRHHte9cfxgT30SwhmUFMWQ5Fb0bx/Jk99s5rZhHfGzCvO35HDlgASiQgOw2U2txYIxhjtmrWLNnkPERwSzLvMwAFenJHBl/wSumXHiaQlgSHI0G7MKOXy0otr0AKuFcpu91ufy4x9+Q5fWYeQUlTLoqeO9jOHB/rxxUwqJUcFEhwZytMJGgNXCNTOWsj7zMDFhgcz/43BaBPrx0cq9PPz5hhqP7W8Vhia34pkrehMZ6l+td9kYw8GScqYv2MGA9lGMPafNSda640dCcXklFZV2olsEsvfgEYY993PV7bueGVf12bbbDRNe+5WNWYVVt2947CLCgvzZkVvMqBcWALD4oZHVdro0xvDE7M28szgDgFcm92NMzzYE+FnYur+QuZsPcN3g9kS5bAW12w33frSGb9ef3kksv75nCL0TIjhYUs7HK/fy7PfHN8ZvfHwMK3blMzQ5hpfnpfHqz+lVt/3ugk48MLZbg5b1a1oe3duG1dqjnVdcxuWvL2bvwaPVpq/+64XVnq8xBmPgaIWND5bv4coBCVTa7eQWldG1dRh/+GQd36zbxyW92vLqdf3cNhRJC2R1Sg4frcAYQ0SI4027r+AoVotwxetLqjYT1eX2YUn85ZIe1aYZY8g8dJQt2YX8Z+FOVu0+VOf9/a3CExPPoXNsCw4UllVtojqmV3w47946iKjQAOx2g80YftmWy/rMAu6/sAv3f7KOL9ZkVbtPx5jQWjfluUppH0mfxAje+nUX/lZheJdYnp/Um09TM5n2w7aqL4CXJ/dj3pYDrN1bQEyLQFJ3H+L5SX24akACdrvhrV938dQcx5D5YH8r4/u0pbisknlbciirrP1LBBxnQLxreCcu6BrDVdOXsufgEQa0j2RfwVGyD5dW20QH0CLQr0Gnr/3mnqGMf/XXGtOvH9yOmLBAXpqbVm16oJ+F0d1bM3fLgRq548KDiG0ZxNq9BcRHBPPkZT2ZvzWH3flHqo0DO/bl6E52uzmlInZ9ZgGPf7O5zvdaoJ+l1tfjqgEJTB3eCT+LEBcRzK/pufRLjKwxbGd3fgnDp/1S9X/rloE8fXkvRnVv+I8CLZDViUorbJRW2Kra4GNe/Gk7/5qXVut9/K1Chc3xvR0SYOXqlERmLsngozvO5dyO0YBji+Fbv+7k6Tn1jlisMvveoXRu3YKnvt3Ce0sdPcBJrUL500Vd+d/qTOZvzQFg0oAEHpvQkwA/C1+syeKBz9bXKLBmLt5Fhc1w0/nt6xwqYbcb5m/NoXdCOLEtj/eufrU2izkbsmkbHszKjINs2ldY7X4XdI3h0t5x7M4vYcmO/Bqf+6UPj6RteO1HB9l/uJRb3lnB1v1F9a6LqcM7cbCkjE9SMzmvYzRLazm8Wf92EazeU0Cwv5VFD46otvXtGJvdcOGLC6p9J7l+R8WFBzHz1kH8b1UmpRU2jpTb+HRVJuD4/skvLuOqAQn8ffYWlu3K56bzOvBp6l78rRbaRYVwbscoducfIaVDFFP/u4p2USGM7BbLzCUZVcsL9rdytMJWI1t0aAABfhauHdiO+0a7f4fCTfsO89cvN7J6T0HVtMv6xvHStf2otNmZvT6bBz5bX+cPrsFJUSzfdZCxPdswbVJvwty4JUULZHVGHDuRHGBocgxvLtrJhyv2UFhayZ2/6cgN57bnx80HuOm89vVu5jfGsSnv8NEKbns3lYTIEN68OYXcojIufaVmAXeMv1VY+ZfRNb4wanv81XsK+GnzAdpFhTB5UCIiwq68EnbmFjMwKYqDxeW8PD+NQyXlTEpJ5L6P1lR9sTTU3SM68X9jqv/Knrl4F499s7nGvEOSoykus5Ec04LVew4xvk8cGMOUoR2rbZ4zxjBj4U6+WJNFaYWNjPwjVbddN7gdT1/eC3D8kLn3wzUs3J7LPSOSueuCTvzju62MPacN53eKZmFaHuv3FlSNMQd47srejDmnDXnFZQT5W6t6N7buL8QYR8+OwdFQBvlbsdsNRaWV/HbmCopKK3n/9sHEhjm+uI71bLvKPnyUxen5PPb1JvyswqhurVmYlss7twzknPjqQzBsdsNzP2xl9e5DjOgWy6W94mgXXfcOLIvT8/jDx2u5ckACfxjdpdrmQdd1tyO3mJvfXklWwVEu7xdPtzZh/Lwth7+N78lT327h1/TjRXy/dhHkF5eTHNui6su+Ns9d2Zt//rSd/YWltG4ZyKGSCsptdm4bmkSAn4W7RyQTehrDK0ALZHXqjDEs33WQvokRLNiey8cr99Z439bWo3tJr7bkFJWyMuN44TiwQySPXNKDbQeKGNOzDblFZdzzwWpyisp49bp+nN+pVdW8x3pFr0lJ5NmreldN/zUtj+W78rlvVOdG3ZFuT/4RKu12bnsvtc7Oj0t6tSUmLLCqMHzq8nO4fnB7wDGcINjfiojw1y831jn841/X9iWlQxRD/jG/xm3XDW7H3yeeg80Yps5axTzn6xDkb+HZK3vXukPiMcYYHvt6E+8urb7cY0eFOFF8RDAf3D6Y9tHVx6CfbHji/R+v5XOXDqPR3Vvz4NiudG4dxgOfreOT1Myq27q3bck7twwkNizQ4/vR2O2G0kobo15YQPbhUq5JSeSLtVmU19OJdMzkQe145opebs+kBbLyacYY9heWMnNJBrPXZZNVcJQ/XdSFiX3jKbfZq41pcqdt+4vYe/AIC7bncvP5HegUE8obi3ayKC2P8ko794xMZljnGHbmFvPS3DTahgfRNjyI+MgQRnWLrbUxMcaQX1LOV2v3kdQqhBFdY89oU5Axho1ZhfSIa1ltk2mFzc7O3BK6tqm7p/bpOVuYsXAnr1/fn3G92p7W8m12g0U45eewK6+ES19eVLWjZoCfhWlX9WbulhyC/S20bhnEF2uyyDx0vFfcInD7sI5c0DWWAe0jCfCzkJFXwsvz08jIK6nW6wDQKSaU168fwIcr9nCgsJSsgqOsd24qBup8vja7YV1mAb3jw6t9qR8qKWdReh6Ltufy6apMurUJq/YDxc8iJLUKZf/hUtpGBPHYhJ7ViojTpQWyOlPGGNbsLSAhMpjYsCDmbMjGIkJcRBATXl1cbd4Aq4Uv7j6fnnEN22fAdccwX3KgsJTb30tlQPtIyivt9E2M4Ir+CVXt5KTpS6p+GPSKD+fO4R2554M1/Pv6/iTHtmDcy4sY3zuOxyf2JCTAr9YhKT9u2s99H63l3I5R3Hx+B7buL+K3QzpU9YTnFJZy8b8WcdcFnbht2KnvP7In/wiRof7MXJzBxb3akBwbxqylGfz1q02EB/sztmcb5m/L4ZM7zyPpNHbQPFJeydwtOXSIDiHAz0K3Ni2r3Z5TVEqrUM8XxHXZk3+E30z7udq05X8eReuWQVXF//cbs0mMCmFHbglzNx/guat6e2TsuRbIqklpyM57qm52u6GorLLGzpGedqikHItF2LyvkIc+X89ul55wcBTEfx7XnSlDk1i+6yD3fLCGvOLjOyddP7gd75+wl/Ivf7qAP3yyljUnFMvg2DHDZjfERwRz1wWduH5wO7e8fw4UlvL+st2M7xNX7177p0sLZOVppRU2th8oold8+FnZpm7ad5hLXq6+hTImLJDcIkd7s+ShkcT50AlazmSn36bmYEk5a/YcosJmTmmsuKdogayU8opN+w4zZ0M2l/eLZ1feEYL9rZzXKbpab83avQW8NHc7v2yrftr4W4ckMXV4RyJDA/C3Wqiw2bHZDd+uz+aPn64DYP4fh5PUKpTiskq3jktrDFogK+V5+cVlvLskgz0Hj7Azr6Rqa9OZbFlTzYcWyEqpJsFuNyxIy6WkrJJLe8fVOd/SHfnERQTVGJvXlGiBrFTj2ldwlBd/2s59ozvrWQsV4CNn0lNKqZOxWIQRp3AmpPM6RTdCGqVUcxIXEcy0SX28HUM1Ub43El8ppZRbichYEdkmIuki8lAtt7cTkZ9FZI2IrBeRcd7IqZRSvkILZKWUasZExAq8BlwM9AAmi0iPE2Z7BPjEGNMPuBZ4vXFTKqWUb9ECWSmlmrdBQLoxZqcxphz4CJh4wjwGOHZsqHCg5kFalVLqLKIFslJKNW/xwF6X/zOd01w9BtwgIpnAHODe2h5IRO4QkVQRSc3Nza1tFqWUaha0QFZKqeattgPinnj4osnATGNMAjAOmCUiNb4fjDEzjDEpxpiUmJgYD0RVSinf0OSOYrFq1ao8Ean9fJH1awXknXQuz/OVHOA7WTRHTb6SxVdygO9kOd0c7d0d5BRld2DuygAACSBJREFUAoku/ydQcwjFFGAsgDFmqYgE4XiedZ6X+zTbYl95DcF3svhKDvCdLJqjJl/J4is54PSynHI73OQKZGPMaXVbiEiqLxyD1FdygO9k0Rw1+UoWX8kBvpPFV3I0wEqgs4gkAVk4dsK77oR59gCjgJki0h0IAuodQ3E6bbEvrTtfyeIrOcB3smiOmnwli6/kAM9n0SEWSinVjBljKoF7gB+ALTiOVrFJRJ4QkQnO2f4I3C4i64APgVtMUzuLlFJKuVGT60FWSinVMMaYOTh2vnOd9qjL9c3AkMbOpZRSvups6kGe4e0ATr6SA3wni+aoyVey+EoO8J0svpKjKfKldecrWXwlB/hOFs1Rk69k8ZUc4OEsolvRlFJKKaWUOu5s6kFWSimllFLqpLRAVkoppZRSykWzL5BFZKyIbBORdBF5yAOPnygiP4vIFhHZJCL3Oac/JiJZIrLWeRnncp+HnXm2icgYd2YVkQwR2eBcZqpzWpSI/CQiac6/kc7pIiIvO5e3XkT6uzzOzc7500Tk5gZm6OryvNeKSKGI/L6x1omIvC0iOSKy0WWa29aBiAxwruN0531rOxFDXTmmichW57K+EJEI5/QOInLUZd1MP9ny6npODcjittdDRJJEZLkzy8ciEtCAHB+7ZMgQkbWeXidS9+e20d8nZ4vT+Sw38PG1La6ZQdviunM0eltcR45Gb4fryaJtsStjTLO9AFZgB9ARCADWAT3cvIy2QH/n9TBgO9ADx6lb/1TL/D2cOQKBJGc+q7uyAhlAqxOmPQc85Lz+EPCs8/o44DscZ9o6F1junB4F7HT+jXRejzyD12A/joNzN8o6AX4D9Ac2emIdACuA85z3+Q64uAE5LgL8nNefdcnRwXW+Ex6n1uXV9ZwakMVtrwfwCXCt8/p04K5TzXHC7S8Aj3p6nVD357bR3ydnw6W+944bl6Ft8clfA22Lj09r9La4jhxuey04xXa4riwn3H7Wt8XNvQd5EJBujNlpjCkHPgImunMBxphsY8xq5/UiHMcZja/nLhOBj4wxZcaYXUC6M6cns04E3nVefxe4zGX6e8ZhGRAhIm2BMcBPxpiDxphDwE84z7J1GkYBO4wx9Z1xy63rxBizEDhYyzLOeB04b2tpjFlqHJ+891we66Q5jDE/GsdxaQGW4TirWZ1Osry6ntMpZalHg14P56/xkcBnJ8tSXw7n41yN4zi8dXLHOqnnc9vo75OzhLbFx5epbfFZ2hb7Sjt8sizaFjs09wI5Htjr8n8m9TeYZ0REOgD9gOXOSfc4NwG87bJ5oa5M7spqgB9FZJWI3OGc1toYkw2ONyMQ20hZwHHWLtcPmTfWCbhvHcQ7r7sj0604fs0ekyQia0RkgYgMc8lX1/Lqek4N4Y7XIxoocPmyOd11Mgw4YIxJc5nm8XVywufWF98nzYG2xQ7aFvvmZ8zbbbEvtcOgbTHQ/Avk2saZeOS4diLSAvgf8HtjTCHwb6AT0BfIxrG5or5M7so6xBjTH7gYuFtEflNfbE9mcY5/mgB86pzkrXVSb8wGLttd6+YvQCXwvnNSNtDOGNMPuB/4QERaumt5dXDX6+GujJOp/gXu8XVSy+e2zlnrWKY337tNibbF2hafNGYDl91c2mJfa4dB22Kg+RfImUCiy/8JwD53L0RE/HG8sO8bYz4HMMYcMMbYjDF24A0cm0Xqy+SWrMaYfc6/OcAXzuUecG5mOLZJJKcxsuD4YlhtjDngzOSVdeLkrnWQSfVNcQ3O5Nx54FLgeucmH5yb0fKd11fhGGPW5STLq+s5nRI3vh55ODZz+Z0w/ZQ573sF8LFLPo+uk9o+t/Xcv9HfJ82MtsXaFh/jM58xX2iLfakddubXtvgYcxqD/ZvKBceptHfiGOB+bDB7TzcvQ3CMaXnphOltXa7/AcdYIoCeVB94vxPHoPszzgqEAmEu15fgGK82jeqD3Z9zXr+E6oPdV5jjg9134RjoHum8HnUa6+Yj4LfeWCecsFOBO9cBsNI577EB/+MakGMssBmIOWG+GMDqvN4RyDrZ8up6Tg3I4rbXA0fPlOvOIb871Rwu62VBY60T6v7ceuV90twv9b133LgMbYvrzqNtsY+0xbXk8Eo7XFsWl/WibbExzbtAdq6YcTj2itwB/MUDjz8UR3f9emCt8zIOmAVscE7/+oQPwV+cebbhsjflmWZ1vnHXOS+bjj0GjrFJ84A0599jbxoBXnMubwOQ4vJYt+LYKSAdl4a1AVlCgHwg3GVao6wTHJuGsoEKHL8ep7hzHQApwEbnfV4FxxkpTzFHOo5xUsfeK9Od817pfM3WAauB8SdbXl3PqQFZ3PZ6ON97K5zP71Mg8FRzOKfPBKaeMK/H1gl1f24b/X1ytlzqeu+48fG1La49i7bFPtIW15Gj0dvhurI4p89E22KMMXqqaaWUUkoppVw19zHISimllFJKNYgWyEoppZRSSrnQAlkppZRSSikXWiArpZRSSinlQgtkpZRSSimlXGiBrM5KIvJ7EQnxdg6llDqbaVusfJUe5k2dlUQkA8fxE/O8nUUppc5W2hYrX+V38lmUatpEJBT4BMcpJq04DqAeB/wsInnGmBEichHwOI6zFu3AcZDxYmfj/TEwwvlw1xlj0hv7OSilVFOnbbFqSnSIhTobjAX2GWP6GGPOAV7CcS72Ec4GuRXwCDDaGNMfSAXud7l/oTFmEI4z8LzUyNmVUqq50LZYNRlaIKuzwQZgtIg8KyLDjDGHT7j9XKAHsFhE1gI3A+1dbv/Q5e95Hk+rlFLNk7bFqsnQIRaq2TPGbBeRATjO7/6MiPx4wiwC/GSMmVzXQ9RxXSml1CnStlg1JdqDrJo9EYkDjhhj/gs8D/QHioAw5yzLgCEikuycP0REurg8xDUuf5c2TmqllGpetC1WTYn2IKuzQS9gmojYgQrgLhyb574TkWzn2LdbgA9FJNB5n0eA7c7rgSKyHMcPyrp6NpRSStVP22LVZOhh3pSqhx6CSCmlvE/bYtXYdIiFUkoppZRSLrQHWSmllFJKKRfag6yUUkoppZQLLZCVUkoppZRyoQWyUkoppZRSLrRAVkoppZRSyoUWyEoppZRSSrn4/wD6lNN9H1ajAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(10,6));\n",
    "ax1 = plt.subplot2grid((2, 2), (0, 0), colspan=2);\n",
    "ax1.plot(losses);\n",
    "ax1.set(xlabel=\"step\", ylabel=\"loss\", title=\"ELBO\");\n",
    "ax2 = plt.subplot2grid((2, 2), (1, 0));\n",
    "ax2.plot(a);\n",
    "ax2.set(xlabel=\"step\", ylabel=\"a\", title=\"Mean Parameter\");\n",
    "ax3 = plt.subplot2grid((2, 2), (1, 1));\n",
    "ax3.plot(b);\n",
    "ax3.set(xlabel=\"step\", ylabel=\"b\", title=\"Standard Deviation Parameter\");\n",
    "plt.tight_layout();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now this has converged, we can use the optimised values of `a` and `b` to sample from the guide. We can also compare this to the analytic speed, which can be calculated for this simple example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean = 3.41, Standard Deviation = 0.89\n"
     ]
    }
   ],
   "source": [
    "print(\"Mean = {:.2f}, Standard Deviation = {:.2f}\".format(np.mean(a[10000:20000]),np.mean(b[10000:20000])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlYAAAEICAYAAACdyboFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8VNX9//HXSQKEJRB22UFFRFSoglatFi1a11qtVm21Wu2Pun6/WrvYWnfaqnXrt1/XVsVvXVCrtUqpWvddQQWVTQRZIhBCIBsQsp3fH+feZAgzyZ2ZO5nMzPv5eMxjMjN3PnMmZD585pxzzzHWWkREREQkeXnpboCIiIhItlBhJSIiIhISFVYiIiIiIVFhJSIiIhISFVYiIiIiIVFhJSIiIhISFVZZxBhzgTGm1BhTY4zpn+72xMMYs9AYMzXB5/7bGHN2yE1KWVwRcYwx3Y0xzxljKo0xT6a7PUEkkxeMMfcYY65KQZtSElcSY7SOVedhjFkJDAYagXrgHeB8a+2aAM/tAlQBX7fWLkhlO5NljJkJlFhrf5vutviMMdcCu1trz0x3W0SyjTHmNWAisIu1dnvE/WcBlwAHW2sbjDHnAD+x1n4jLQ1tpTPmhc72O5Kdqceq8znBWtsLGAKUAn8O+LzBQCGwMN4XNI7+FkQkdMaY0cChgAW+0+rhUcDn1tqGkF6rIIw4Ikmx1urSSS7ASmBaxO1jcUnHv90NuAVYjSu67gG6A3sAW3CJqwZ4xTv+YGAuUOldHxwR6zXgd8DbwDZgd6APcD+wDvgKmAHkx2hrN+AOYK13uQPo5j02FSgBfgNs9N7XD73HpuN64+q8tj7X+r0D1wJPAg8D1cCn3nv8NbABWAMc1eq9/MT7eXfgde89bwQejzjuT95zq4APgUO9+4/22lPvtWlBlLh5wG+BVV4b/g/o4z022vvdn+3922wErkz335MuunSGC3C1l2duA2ZH3H9dq8/dRUAtrse+Bqjwjoua97zH/FzzK2A98Lcorx/kszvdy2PrgMu9x4LkhXO893Y7UAGswOXdc7xcswE4O6ItM4EZ3s8DgNne8zYBbwJ53mNXAMtx+W8RcJJ3//gYv6PmuN7t/wd84cV9Fhga8ZgFzgeWAZuBO/FGr3QJ56Jeik7KGNMDOA14L+Lum3AFxiRcATEMuNpa+zkwwTum2Fp7hDGmH/Av4H+A/rik9q9Wc6/OwiWUIlzSeQho8GJ/DTgK+EmMJl4JfN1ry0TgAFzy8u2CSxzDcAXHfcaYcdba+4BHgJuttb2stSfEiH8C8DegL/Ax8AIuQQ4DrgfujfG8G4AXvecNZ8cev7lee/sBjwJPGmMKrbXPA7/HFWG9rLUTo8Q9x7scDuwK9AL+t9Ux3wDGAd8CrjbGjI/RRpFc8iPcZ/4R4NvGmMEA1tpr2PFzdyfuP/x3vdvF3vOj5r2I+LvgPtOjcPmstXNo/7N7ODAWl/OuMMZMC5gXAA4EPsHl2UeBWcAUr61nAv9rjOkV5XmX44rCgbgRh9/gih5wRdWhuC+71wEPG2OGWGsXx/gdNTPGHAH8Afg+buRjldemSMd7bZzoHfftGO9NEqDCqvN5xhhTgetVORL4I7jhOty3kMustZustdW4D/3pMeIcByyz1v7NWttgrX0MWIIrWHwzrbULreuG7wccA1xqrd1ird2A+xYWK/4PgeuttRustWW4D/9ZrY65ylq73Vr7Oq7I+34cv4c3rbUveG17Epd8brTW1uOSxGhjzE5JBfftchTuG1qttfYt/wFr7cPW2nLv93Er7pvwuIDt+SFwm7V2hbW2Btd7dnqroYfrrLXbrJvjtgCXtERyljHmG7jP4xPW2g9xBcMP4nh+kLzXBFzj5ZptUcIE/exusdZ+CjwInBHH2/zSWvugtbYReBwYgcuN2621L+J6vXaP8rx6XOEzylpbb61901qvS8naJ621a621Tdbax3G9SwcEbM8PgQestR9ZN5/t18BB3pCs70ZrbYW1djXwKq5olZCosOp8vut9C+kGXAy8bozZBVdY9AA+NMZUeMXX89790QzFfVOJtAr3bc8XOSl+FNAFWBcR/15gUMD4q7z7fJuttVvaeLw9pRE/bwM2eonLvw3um2drvwQM8IF3puG5/gPGmMuNMYu9M5AqcN8GBwRsT7T3W4D7pulbH/Hz1hjtE8klZwMvWms3ercf9e4LKkjeK7PW1rYRI8hnd02rx5PJVVhrW98XLRf8ETdc96IxZoUx5gr/AWPMj4wx8yPe894kmKu8YrKcHXO/clUKaaJfJ+UVEU8bY+7FDTE9jfuATrDWfhUgxFpcsRRpJC4pNb9MxM9rgO3AABtsIqkf358sP9K7z9fXGNMzorgaCXwW5XVDZa1dj/uG639bfskY8wbum+GvcMN0C621TcaYzbgiLEibWv8+R+KGTUtxQ44iEsEY0x3XS51vjPH/I+8GFBtjJtroZy+3/hxupP28F8ZndwSuR99/3M9lqcxV1bjhwMuNMROAV40xc3HF1l9wuepda22jMWY+CeYqY0xP3DBlkP83JATqseqkvDP1TsTNFVpsrW3CfdhuN8YM8o4ZZoyJNTY+B9jDGPMDY0yBMeY0YC/cZMmdWGvX4eYm3WqM6W2MyTPG7GaM+WaM+I8BvzXGDDTGDMDNeXi41THXGWO6GmMOxY3p++vUlOLmOoTOGHOqMcZPlptxSagRN4+sASgDCowxVwO9I55aihtejPWZeAy4zBgzxpsv4c+9COVsJpEs9F3cZ28v3FDTJNzk6zdx866iKQWGG2O6AiSQ96IJ8tm9yhjTwytwfowb0vPb01ZeSJgx5nhjzO7ecGcV7nfVCPTE5a0y77gf43qsfDv8jqJ4FPixMWaSMaYb7v2+b61dGfZ7kOhUWHU+zxljanAftN/hzijxe4V+hfs2854xpgp4iRhzhKy15bhi5nJcN/AvgeMjuuSj+RHQFXcWymbg77ienmhmAPNwkzY/BT7y7vOt92KsxU1aPd9a638jvB/Yy+vmfqaN9iRiCvC+9zt8Fvhva+2XuMnv/wY+x3WT17Jj979f9JUbYz6KEvcB3GT6N4AvvedfEnLbRbLJ2cCD1trV1tr1/gU3cfyHMZZGeAXXC77eGOPnqsB5L4Ygn93Xvdd4GbjFmxsF7eeFZIzFvZca4F3gLmvta9baRcCt3n2lwD64Mw990X5Hzay1LwNXAU/hznLcjdhzZSUFtECohM64FdQfttZqiExEOi1vQveXQBf1PktY1GMlIiIiEhIVViIiIiIh0VCgiIiISEjUYyUiIiISkrStYzVgwAA7evTodL28xGPpUnc9Lp4TcWKFWuqFSj6WZJ4PP/xwo7U21qK2GUP5K7NsnT8fgB6TtMC4JC5o/kpbYTV69GjmzZuXrpeXeEyd6q5fey2EUFO9UMnHksxjjGm9G0BGUv7KLPOL3e5Xk/RvJkkImr80FCgiIiISEm1pI+377W9DDBVeLBGRIJp+85t0N0FyiAorad+0aSGGCi+WiEgQ+/3yl+luguQQFVbSPm/iJyFM/JzvxZqUxZNI6+vrKSkpoba2Nt1NSZvCwkKGDx9Oly5d0t0UEZY+7rb+G3faaWluSWbI9RyWbP5SYSXtu/RSdx3ChPNLvVjZPHm9pKSEoqIiRo8ejdtfNbdYaykvL6ekpIQxY8akuzkibPvpT90PKqwCyeUcFkb+0uR1kZDV1tbSv3//nEtIPmMM/fv3z9lvuyKZLpdzWBj5S4WVSArkYkKKlOvvXyTT5fJnONn3rsJKgmlshI8/Bm2BJCIZqLCxESor090MyQEqrKRt8+bB4sXw9tuw337gTQKVzi0/P59JkyY1X2688cbQYs+fP585c+aEFk8kZSoq4Kqr2Kuqij1rauAHP0h3iySgTM5hmrwusVkLJ54IGza03Pfcc3D66QmH/P3vfx9Cw6Q93bt3bz4DM2zz589n3rx5HHvssSmJLxKaG2+Em26iq3/7P/+Bbduge/d0tkoCyOQcph4riW3VKli7Fvr1g9mz3X2vv57UcODBBx/MwQcfHFIDJR6VlZWMGzeueb/GM844g7/85S8AXHDBBUyePJkJEyZwzTXXND9n7ty5HHzwwUycOJEDDjiAyspKrr76ah5//HEmTZrE4+rBlE7s86efdz/cfTdMnAj19fD+++ltlCQsY3KYtTYtl/33399KJ/fEE9aCtV//urVvvmltv37u9ooVCYd8++237dtvvx1iIzufRYsWtdxwZWj4l3bk5eXZiRMnNl9mzZplrbX2xRdftF//+tftY489Zr/97W83H19eXm6ttbahocF+85vftAsWLLDbt2+3Y8aMsR988IG11trKykpbX19vH3zwQXvRRRfF93to/nUwz6Yp54R5Uf7KAPX1trZrobVgF/7xj3bD4Ye7z85116W7ZZ2eclhy+UtDgRLb3Lnuet06t63NoYfCP//peq0SXN/jN97WEtm8jlVnEKsb/cgjj+TJJ5/koosuYsGCBc33P/HEE9x33300NDSwbt06Fi1ahDGGIUOGMGXKFAB69+7dYe0XSdrixXSrq6Vs4FDqZsygpr6egeDyl3R6mZzDNBQosX3wgbsuKnLX3/ymu37jjfS0JxOl6vtegpqamli8eDHdu3dn06ZNAHz55ZfccsstvPzyy3zyySccd9xx1NbWYq3N6VOuJcN5+evL0XsBUFPg9SO8+y7U1aWrVZlHOSxuKqwkusZG+PBD97NfWB12mLtWYZWxbr/9dsaPH89jjz3GueeeS319PVVVVfTs2ZM+ffpQWlrKv//9bwD23HNP1q5dy1yv57K6upqGhgaKioqorq5O59sQaZ/3d/vlrhMAaDQGxo93k9f93CYZJxNymAoriW7JEqipgZEjoat3Ts3Eia7IWr4cvvoqve2TNm3btm2HU5WvuOIKPv/8c/76179y6623cuihh3LYYYcxY8YMJk6cyNe+9jUmTJjAueeeyyGHHAJA165defzxx7nkkkuYOHEiRx55JLW1tRx++OEsWrRIk9elc/N6rFaOGd9yn74cZoxMzmGaYyXR+fOrDjgAysrczwUFcMgh8Pzz8OabSS27IKnV2NgY9f7Fixc3/3zbbbc1/zxz5syox0+ZMoX33ntvp/v9b4AinVJtLXz6KU0mj1Wj9my5/7DD4N57XWH1q1+lr33SrkzOYSqsJDr/j27KFDjqqJb7DzvMFVZvvJFQYXXHHXeE1EARkRjmz4eGBtYO243thT3ofu+97n6vJ4O33nLTHfLz09dGyVoqrCQ6f+L6AQfApEkt9/trUCU4R2FSZCwRkVTwJ67v6iaujzvttJbHRoyANWvclIY99khH6yTLaY6V7Gz7dliwAIyB/feHl15yF4Ddd3fXK1cmFPqll17iJT+WiEgq+BPXx7jC6qObb+ajm292jyWZw0Taox4r2dlnn7kVisePd5PVZ8xw90+bBkOGQJcubpubrVuhR4+4Qs/wYk2bNi3sVouIOF6P+kqvsMrzt9L65S9h9Gj3sworSRH1WMnOVqxw13vuufNjeXkwapT7edWqjmuTiEgQ1jbnsPW7jNr5cRVWkmIqrGRnq1e765Ejoz+uxCQinVVZmZvO0K8f2wuj9Kgrf0mKaShQdqbCKlTnzQz3tN77z5kS6LjS0lIuu+wy3nvvPfr27UvXrl355S9/yUknnRTzOcceeyyPPvooxcXFO9x/7bXX0qtXL37+858n1XaRlFP+Cl06clgm5y/1WMnO/CE+JaaMZa3lu9/9LocddhgrVqzgww8/ZNasWZSUlLT5vDlz5uyUlEQySpTCamtdI1vrGl2BoPzV6WV6/lJhJTtrnZjuvdddfEkkpnvvvZd7I2NJSrzyyit07dqV888/v/m+UaNGcckllzBz5kwuvvji5vuPP/745k2xR48ezcaNGwH43e9+x7hx45g2bRpLly5tPn758uUcffTR7L///hx66KEsWbKkY96USBBRCqtZ589g1vneSTjDhrnFjtetcwuJSqeT6flLQ4Gys9aJady4HR8fM8ZdJ1BYjWsdS1Ji4cKF7Lfffgk/3/+G+PHHH9PQ0MB+++3H/vvvD8D06dO55557GDt2LO+//z4XXnghr7zySlhNF0lOlMJqy74Htzyen+8eW7HCHau1rDqdTM9fgXqsjDFHG2OWGmO+MMZc0cZxpxhjrDFmcnhNlA61ZQuUl7v9AQcNcvc995y7+JLosXruued4LjKWdIiLLrqIiRMnMmVKsPlZb775JieddBI9evSgd+/efOc73wGgpqaGd955h1NPPZVJkybx05/+lHXr1qWy6UlT/soxUQqroU/fw9Cn72k5xstht933QujzhyR8mZa/2u2xMsbkA3cCRwIlwFxjzLPW2kWtjisC/gt4P/RWSsdZs8ZdjxjhllYAuPVWd33CCe56l11c4bVhgyvEevYMHP5WL9YJfixJiQkTJvDUU081377zzjvZuHEjkydPpqCggKampubHamMMhxhjdrqvqamJ4uJi5s+fH36jU0D5KwdFFlbeCNAxLz0BwP0ne0NLXmE1oGxtBzdOgsj0/BWkx+oA4Atr7QprbR0wCzgxynE3ADcDGrTOZO2dUQNayyoDHHHEEdTW1nL33Xc337d161bAzUOYP38+TU1NrFmzhg/87YsiHHbYYfzjH/9g27ZtVFdXN/cy9u7dmzFjxvDkk08CbpLpggULOuAdJUz5K9e0d/INNBdW/Td27t7WXJXp+SvIHKthwJqI2yXAgZEHGGO+Boyw1s42xsQ8n9EYMx2YDjCyrT96SZ8ghRW4xLRsmRsO3GuvVLcqowVdHiFMxhieeeYZLrvsMm6++WYGDhxIz549uemmmzjkkEMYM2YM++yzD3vvvXfUuQz77bcfp512GpMmTWLUqFEceuihzY898sgjXHDBBcyYMYP6+npOP/10Jk6c2JFvLx7KX7lk2zbXk15Q4HrWidEj5fdYlauwCqKjc1im568ghdXO/Wlgmx80Jg+4HTinvUDW2vuA+wAmT55s2zlc0iGewgp0ynInNmTIEGbNmhX1sUceeSTq/Ssj/j2vvPJKrrzyyp2OGTNmDM8//3wobewAyl+5xD8df/hwN0k9Fr+w2qihwM4qk/NXkKHAEmBExO3h7Pg1oAjYG3jNGLMS+DrwrCaAZigVVpJdlL9ySZz5S0OBkgpBeqzmAmONMWOAr4DTgR/4D1prK4EB/m1jzGvAz62188JtqnSIaInpb3/b+bgEC6u/RYslkjrKX7kkRmH18OX/s+NxQ4dCQQHFleV0qdO0OglXu4WVtbbBGHMx8AKQDzxgrV1ojLkemGetfTbVjZQOFC0xjRix83EJFlYjosXKQtbaqGel5AprO8dImfJXjolRWG3fbe8dj4tYy6pfeWkHNS6z5HIOSzZ/BVog1Fo7B5jT6r6rYxw7NakWSfo0NbUstxCZmB5/3F2fdlrLfQkWVo97sU6LjJVlCgsLKS8vp3///jmZmKy1lJeXU1hYmO6mAMpfOcUvrPyzlj2jHnHLvKz64eUtd44eDStWaJ5VFLmcw8LIX1p5XVqUlkJdHQwYAD0idoX3T3mNLIb8tazKyuJay8o/fTabC6vhw4dTUlJCWVlZupuSNoWFhQwfPjzdzZBcE6PHatrb/wLg/taFFZrAHk2u57Bk85cKK2kRdOInuLWshg2DL790e27tvntq25ZBunTpwhh/2x8R6Tjx5DDvP87iivIUNigzKYclR5swS4t4khJ468QA69enpj0iIkFZ25LDgszl9PJX70oVVhIuFVbSQoWViGSqjRuhthb69oWiovaP9/JXHxVWEjIVVtIicp/AIFRYiUhnkWD+6lOlwkrCpTlW0sIvkIYM2fH+v/89+vEJFFZ/jxVLRCQZXh76tKkHd8ycu8NDD131YPPP53mPDSgr5SbUYyXhU2ElLUq99VwGD97x/gEDdj428rg4CqsBsWKJiCTDy19Vvfvt9FDDkFE73VfVuz/gzbGyFnJsWQFJHQ0FSgu/QGpdWM2c6S6tJdBjNXPmTGZGiyUikgwvD1X12bmw2u2BGez2wIwd7qvrVsi2wp50aaiHiooOaaLkBhVW0iJWj1V7hVVp8JWLVViJSEo091j13+mhqfNeYeq8V3a6v9IvwuLIYSLtUWElTl0dbN7s1qfqv3NiikqT10Wks2hjKDCWqj5erlMOkxCpsBJnwwZ3PWiQ20crCL9nq7TUbYcjIpIuXmFV2SfgF8PIY1VYSYhUWIkTaxiwLYWFUFwM9fWut0tEJF0S6LFSYSWpoMJKnEQKK9BwoIh0Dm1MXo9FhZWkgpZbECfWGYEAc+bEft4uu8CSJe75Eya0+zJz2oolIpKI+nrYtAny8qjp1Wenhx+86emoT2ue6K7CSkKkwkqctnqsevSI/bw417Lq0VYsEZFE+HNEBw7E5u08R7SpqDjq05p7t1RYSYg0FCiOX1j5Q3uR7rrLXaKJcyjwrrvu4q5YsUREEtFW/gLG330l4+++cqf7NRQoqaAeK3Ha6rF64gl3feGFOz8W51pWT3ixLowWS0QkEe3MET3403cBWNzq/ubCSutYSYjUYyWOJq+LSKZqa45oG6qLvKHADRugsTHkRkmuUmEljgorEclUCeavxoICqnsVu3X4Nm5MQcMkF6mwEkeFlYhkqkTzF5pnJeFTYSXuVOXycredzYAB8T1XhZWIpFs7k9fbosJKwqbJ67LDqcpRt7N57bXYzx0wAIxx3ej19dClS5sv9VpbsUREEhHZY/XVzg/ff9fOGzD7mldqV2ElIVGPlSTVjU5BgSvIrIWysnDbJSISRIKT10EbMUv4VFhJ+4XVLbe4SyxxDAfecsst3NJWLBGReLWTw/a9/Wfse/vPoj6moUAJmworab+wmj3bXWKJYy2r2bNnM7utWCIi8QgwR3TKsvlMWTY/6mNay0rCpsJKkhsKhJbCat26cNojIhKUNwWhslcx5/3to7if3lxYKX9JSFRYSUsXeAJn1AAwaJC79ifBi4h0FO+LYfO+f3Gq7t3X/aD8JSFRYSXJ91gNHOiuNXldRDqa98Wwqnf/hJ5eXeQVVspfEhIttyDtF1bdu7f9/DgKq+7txRIRiYffY9U7do/V9vzY/9XV9Cp2P5SXuxXY89TfIMlRYSXtF1b//nfbz/eHAgMUVv9uL5aISDwCFFYP//nFmI81FhRA376weTNs2hT/Iskirag0l5a5BRoKFJFM4+WvygTnWAHKYRIqFVa5rrGxZfPR/jHmKNxwg7vEEkdSuuGGG7ihrVgiIvHwCqsaf65UFJNvvpDJN18YO4YKKwmRCqtct2mTWzW9b9/Y29G8/LK7xOInpQ0bXKw2vPzyy7zcViwRkXh4xVB1UXHMQyauXMLElUtix4jMYSJJUmGV6/xvaH5iSUTPnlBYCLW1sGVLOO0SEQmiubCK3WPVLvVYSYhUWOW6MAorY+KawC4iEpoAPVbtUv6SEKmwynVhFFaRz1diEpGOErH5u3qspLPQcgu5LkhhFWtSe6SAial/kFgiIkFs2QK1tdR16cb2brHXyKts4zFAhZWEKlBhZYw5GvgTkA/81Vp7Y6vHzwcuAhqBGmC6tXZRyG2VVAhSWD31VPtxAk7+fCpILJEQKX9lschhQGNiHvb32//VdhxNXpcQtTsUaIzJB+4EjgH2As4wxuzV6rBHrbX7WGsnATcDt4XeUkkNDQVKFlP+ynJhDAOC8peEKsgcqwOAL6y1K6y1dcAs4MTIA6y1VRE3ewJtn3MvnUeQwurXv3aXtgRMTL/+9a/5dXuxRMKj/JXNAk5cP/iGczn4hnNjH6DCSkIUZChwGLAm4nYJcGDrg4wxFwE/A7oCR4TSOkm9IIXVu++2HyfgWTXvBoklEh7lr2zm5Zu2FgcFGL9uJQDvxDrAz38bN2q/QElakL+eaAPXO32js9beaa3dDfgV8NuogYyZboyZZ4yZV6ZvBp2DhgIluyl/ZbMwlloA6NoV+vRxO1FUVITQMMllQQqrEmBExO3hwNo2jp8FfDfaA9ba+6y1k621kwcm+x+5hEOFlWQ35a9sFtYcK1AOk9AEKazmAmONMWOMMV2B04FnIw8wxoyNuHkcsCy8JkrKWNuyT2BYhZXOqpHORfkrm4XVYwXKYRKadudYWWsbjDEXAy/gTld+wFq70BhzPTDPWvsscLExZhpQD2wGzk5loyUkFRXQ0ABFRdCtW+zjhg9vP1bAb3vDg8QSCYnyV5YL2GNV1rOo/VheDrvz8Xf5aHkhAPefMyW59klOCrSOlbV2DjCn1X1XR/z83yG3SzpC0GHAhx9uP5ZfnG3d6i49esQIFSCWSIiUv7JYwMLq2T/+s/1Y3gk4RdWbk26W5Dad+pDL/C7vMOaLGKM5CiLSsVIwFKjCSpKlwiqXBe2xuvRSd2lPgMLq0ksv5dIgsURE2uMXVr3b7rGaes2ZTL3mzLZjNRdWOitQkqO9AnNZ0MJq/vxg8QIUVvODxhIRaUttLdTUQJcubOveq81DdytzJ4K+1tZBXv7qpR4rSZJ6rHKZXwD5i3smS2fViEhHifxi2MY+gYF5+au3CitJkgqrXBbWGlY+zbESkY6SovylOVaSLBVWuSzsxBRwWxsRkaSlKH9pjpUkS3OsclnQxLTHHsHiBeix2iNoLBGRtsRRWK0p7t9+vOY5VhVu8eQwhhclJ6mwymVBE9N99wWLF6Cwui9oLBGRtsRRWL34+yfbj9etGxQVUVBdTfdtNWzrEWBRUZEoNBSYyzTHSkQyVdj5KyJWUZXmWUniVFjlKmuDJ6bp092lPQHOCpw+fTrTg8QSEWlLHIXVUb85laN+c2r7MXVmoIRAQ4G5qroa6urc1jMxtp9p9vnnwWIG6LH6PGgsEZG2RBZWVW0fOqKiPFhMnRkoIVCPVa5KRTd6nz7QpYtbtK+2Nry4IiKtpSKH6cxACYF6rHJVKpKSv1/g2rUu/ogR4cUWEYkUmcOWVycV6ryZcwH43voGjkWrr0ty1GOVq1JRWEXG0wR2EUmlFOSw6iK356CGAiUZ6rHKVf4E8yDb2UyaFDxuO4XVpHhiiYhEU1enUJVmAAAgAElEQVQHFRWQnw99296AGWD5wKGBwlYXFQNQVKWhQEmcCqtcFc+3vTvuCB63nTMD74gnlohINBs3uusBAyCv/YGX1657OFBYv8dKZwVKMjQUmKv8wkdDgSKSaVKUv5p7rFRYSRJUWOUqv/AJMhR45pnuEkQ7+wWeeeaZnBk0lohINPHkL+A7vziR7/zixHaPq+7dD/C2tRFJkIYCc1U8Q4ElJcHjttNjVRJPLBGRaLz88kFNPvd6Z/S1ZeCWYGcN7tBjZW3i7ZOcph6rXBXP5PV4aChQRFLNy19VvdufuB6P+q6F1HbrTpeGegprt4QaW3KHCqtcpeUWRCRTefnFn2weppZeKw0HSmJUWOWiePYJjFeA/QJFRJLi5Re/CApT85mBVZtCjy25QXOsclF1NWzf7vYI7Nmz/eMPOih47HZ6rA6KJ5aISCvnzZzLxR9+zteA6qJ+gZ6zeMjowPFrvMJKE9glUSqsclG8vVV/+EPw2MXFUFAAVVWueOvWrVWoOGKJiEThFz3VvYP1WL1z1QOBY/vztrTkgiRKQ4G5KM5TleOSl+cW7YOWRfxERELUu8oVPVUBe6zioTlWkiwVVrko3sX1vvc9dwmqjeHA733ve3wvnlgiIq34vUk1AedYnXLZcZxy2XGBjtV+gZIsDQXmoniHAsvL44vfRmFVHm8sEZEIBfV19NhWQ2NePlt7FAV6Tp/t2wLHr1FhJUlSj1UuSuVQIOjMQBFJmV413vyqomJsgH0C49WyEbMKK0mMCqtclKp9An1ay0pEUsSf+5SKNawi42qOlSRKhVUuStUaVj4VViKSIv76UkHnV8VLGzFLsjTHKhfFu53Nt74VX/w2NmL+VryxREQi+EstVPUOfkbggtF7Bj7Wj1tUXeEWUzYmvgZKzlNhlYvi7bG66qr44kf0WJ3XeoPUEUcDbpG/+8+ZEl9cEcl5fo9VPKuuz/vlXYGPrevWne1du9Gtbjts2QK9esXdRsltGgrMRR01eV1DgSISslTPsYKWMwOVwyQRKqxyjbXxT14/5hh3CaqNswJfuO2/eeG2/w4eS0Qkgj/3KZ7C6sxLjuLMS44KfHxzb5jObJYEaCgw11RXQ12d2yewR49gz9kWfA0YoM0eq8a67fHFEhGJ4C+DEE9h1a2xIa7XqFaPlSRBPVa5JtXDgAD9+rmtbSoqyG+IL6GJiLTF77Hy9/RLheYeKxVWkgAVVrkm1WtYwQ77BeqUZREJUyJDgfGq9vcg1FCgJECFVa7piB4rgMGDgZYzeEREwtBcWKWwx6qyj1dYlZam7DUkewUqrIwxRxtjlhpjvjDGXBHl8Z8ZYxYZYz4xxrxsjBkVflMlFIksDnr88e4Sj112AXYurEZM/AYjJn4jvlgiSVD+yiLbt9Nj2xYa8oPvEwgwd+wk5o6dFPj4qj793Q8qrCQB7U5eN8bkA3cCRwIlwFxjzLPW2kURh30MTLbWbjXGXADcDJyWigZLkhIZCvz5z+N/Ha/Hqk/ljpsu73PMmfHHEkmQ8leW8b4Y1vQqjmvhzk8uuy2ul2lefHT9+rieJwLBeqwOAL6w1q6w1tYBs4ATIw+w1r5qrd3q3XwPGB5uMyU0qd7OxucPBbYqrEQ6mPJXNvHyVyrnVwFU9tZQoCQuSGE1DFgTcbvEuy+W84B/R3vAGDPdGDPPGDOvTGdbpIf/Dcwbqgtk6lR3iUeMOVZzbjyfOTeeH18skcQpf2UTL381D9UFdN6FR3DehUcEPr5Kc6wkCUHWsYrW32qjHmjMmcBk4JvRHrfW3gfcBzB58uSoMSTFEimsEhFjjlWknba7iaDtbiQkyl/ZxMtfzZPLU6S6qC9NJo+88nKor4cuXVL6epJdgvRYlQAjIm4PB9a2PsgYMw24EviOtVarQHZW/jewVBdWMeZYiXQw5a9s4uWvyjh7rOJl8/LdWlbWai0riVuQwmouMNYYM8YY0xU4HXg28gBjzNeAe3FJSQt/dGYd1WOl5Rakc1D+yiYJDgUmokrzrCRB7RZW1toG4GLgBWAx8IS1dqEx5npjzHe8w/4I9AKeNMbMN8Y8GyOcpNP27bBpE+TnQ/8UJ6bmyesqrCR9lL+yjD8U2FuFlXRegfYKtNbOAea0uu/qiJ+nhdwuSQV/qYVBg9zq6EF9//vxv9bAgZCXR88tleQ3NNBY4P7UxhygPxXpWMpfWcQrcqrinGP1zj4Hxf1SWstKEqVNmHNJovOrLrww/tfKz4cBA8jbsIGi6s1U9HXLO4w/4pT4Y4mIQMTk9fh6rBZf8Lu4X0prWUmitKVNLkl0ftXWre4SryjzrBq219KwvTb+WCIiCRZWedUV5FVXxPUcbWsjiVJhlUv8wsoreAI79lh3iVeUwurF2y/lxdsvjT+WiOS22lqoqPC2s+kd11N//KuT+fGvTo7rOZpjJYlSYZVLOuqMQJ+/lpWWXBCRZHlzRKt698fGM0c0QZpjJYlSYZVLOrqw0pILIhKW5qUWUrs4qE/b2kiiVFjlko5aHNSnRUJFJCwJzq9KVHMBp8nrEicVVrkk0TlWiVKPlYiExe+x6t0xPVY1vYrBGCgvh4aGDnlNyQ5abiGXJDoUeM45ib1elP0Cx37j+MRiiUhuS2I7m9cmB9+A2deUX+DW49uwwW1rM2RI3DEkN6mwyiVJFlZxb5ocZfV1FVYikpAkhgKXn/vbxF5z8GBXWJWWqrCSwFRY5YotW6CmBrp1gz594nvuxo3tHhKt6OpduZHb2XGOVa23lkxhUXF8bRCR3JbEPoEF61YB0DBkVHxPHDwYPv1U86wkLppjlSv8ieuDB7t5A/E45RR3iVNNr2KajKHnlkryGt0chVfuvIJX7rwi7lgikrvOmzmXZZ98ASQ2x+rsG37M2Tf8OP4X9uejlpZy3sy5bfbai/hUWOWKjl5qATdHobqoL3nWUhTnqsciIpH8uZoddVYg0JIvteSCxEGFVa5IQ2EFLd8udWagiCTDn1LQoYVVRI+VSFCaY5UrOnoNK49fWPWpLGdNCPHinkAvIhmvW+1WCrdvY3vXbtQW9uy4F/YLq/XrYZ+Oe1nJbCqsckVHr2Hl8SeaxrutjeYyiIjP7/Gu6t0v/jmiydBQoCRAhVWuSGYo8IIL3PW2+J9aUewKq+IKd2bhnod/L/4gIpLT/GHARM4IBHjpkOMSe2F/iYW1axN7vuQkFVa5IpnC6rTT3HUCvUib+7oesr6b3Qaqux54ZPyvLyI5Ldn5Vat+eHliLzx8uLsuKUns+ZKTVFjlimQKqzWJz46q6DsQaCmsaspdl3qv/h07JCkimat3koVVt+WfAbB9t73je2JxMfToAdXVFG6robZ7r4ReX3KLzgrMFcnMsTrrLHdJwKa+g4CWwuqNv1zDG3+5JqFYIpKbmocCE9wn8Mxb/4szb/2v+J9oDAwbBkDfzWUJvbbkHhVWuaCxsWWOgJckOsrmfl5htWlDh76uiGQP/4vZZu+LWofyhgP9Noi0R4VVLigtdbuzDxgAhYUd+tJVvfvRZPIoqt5MfkN9h762iGSHfpvcFILN/dIwhaC5sFKPlQSjwioX+BMvR4zo8Jduyi+gongAedY2nxkoIhIPv8d7U7809Fj5Q4GbtOSCBKPCKhf4k8/9M1w6WOsJ7CIigVnb0mPVt2N7rM6bOZeHSxoB5S8JTmcF5oJke6wu905Vjm+Nz2ZuXsRCijdvYO+jf5hYEBHJTZs3062ulm2FPdnWI7Gz8v497fsJv7z/xbBYQ4ESkAqrXOD3WCVaWJ1wgrtOcDV0f8Jpv80bGPltFVYiEgfvi2Eyw4BrTz4/4eduishfIkGosMoFyQ4FLl2a1Mtv9ocCN22gct0qAPoMGZVUTBHJEV7+SuaMwJ6fvAPAln0Pjvu5OrNZ4qXCKhckOxT405+663P+mNDT/YRYvLmMtx/6AwDHXnFPYm0RkdziFVabkjgj8PR7fgvA/Xe9Evdzq4v60ZCfT1FNBQX12xNug+QOTV7PBckOBSZpc6tFQkVEAvO+GKZlqQXA5uVRWezPs9KZzdI+FVbZLo2Lg/r8rnTNURCRuDX3WKVhqQVPyzwrLbkg7VNhle3Wr3fF1aBB0K1bWprgz7HqU1GGsTYtbRCRDJXmHivQmYESHxVW2c6fX5WmNawAGrp0o7pXHwoaGxnQ2JC2dohIBvJ7rNKxnY1H0xkkHpq8nu3CmF/1Wzfxk5LEQ2zuO4iimkq+edDRLNtFZwSKSADWtpwVmESP1XPHn5NUM5qnM+jMQAlAhVW2C6OwmjbNXSe4jhVARd9BjFyzjD36DmLrhAMSb4uI5I5Nm6C2lq3de1HbvWfCYTYc+6OkmrG52DuzuUJDgdI+FVbZLoyhwPnzk26G343fdflnlA8YQv+ReyQdU0SyXHNvVXLDgL3f/w8AVQcemdDztZaVxENzrLJdGD1Wl17qLknwE1Pde8/z/qO3JRVLRHJECGtYAZz60B841VtDLxGaYyXxUI9Vtouzx+q8KMN9v1hfnXQzKrx1YHZpqE86lojkCP+MwDROXAeoKB4AQJ/KcmhogAL91ymxqccq26V5cVCf32M1rL4ure0QkQwSUo9VshoLulDZux/5TY2wbl1a2yKdX6DCyhhztDFmqTHmC2PMFVEeP8wY85ExpsEYc0r4zZSENDSkfXFQ34ZBrsdslLaEkA6m/JXBmjdgTm9hBS05jOXL09sQ6fTaLayMMfnAncAxwF7AGcaYvVodtho4B3g07AZKEtavh6YmGDwYunZNa1M29duFxrx8hjTU062pKa1tkdyh/JXhQpq8HoayQd6XUxVW0o4gA8UHAF9Ya1cAGGNmAScCi/wDrLUrvcf0P2Znsnq1u05yGPDp712YdFMaCwoo778Lg8q+4ujDT0YzraSDKH9lMi+HJdtj9fSpFyfdlLKB6rGSYIIUVsOANRG3S4ADE3kxY8x0YDrAyJEjEwkh8Vi2zF3vvntSYZaP3TeExriu9EFlX7FX954sCCXijqJNvPfdf86UFLyiZADlr0xVVwcrV0JeHhsHDE0q1ObDT066Oc1DgV98kXQsyW5B5liZKPcltOGbtfY+a+1ka+3kgQMHJhJC4uEXVmPHJhVmt2WfsNuyT5JuTpmXmAoXfpB0LJGAlL8y1ZdfuqkMI0fS0CW5qQx9X32avq8+nVSMDYPVYyXBBCmsSoDIsaThwNrUNEdCFVJhdfJTd3HyU3cl3Rz/G1/BvFeSjiUSkPJXhvrT3f8CYGHP5Ceun/zk/3Lyk/+bVIwNAyN6rLSZvLQhSGE1FxhrjBljjOkKnA48m9pmSShCKqzC4hdWo3VmoHQc5a8MNbjUjeCWDk7vUjG+mqJithX2hKoqKC9Pd3OkE2u3sLLWNgAXAy8Ai4EnrLULjTHXG2O+A2CMmWKMKQFOBe41xixMZaMlAGs7XWHlDwWOqlNhJR1D+StzdbbCCmO05IIEEmj5WGvtHGBOq/uujvh5Lq6LXTqL0lKoqYG+faF//3S3BoCyge505ZH128lrbKApX6sXS+opf2WmQV5htaGzFFa4JRdGrV7qhgMPTOgcCMkBWnk9W3Wy3iqAum6FrCvoQheg36bSdDdHRDqxTtdjhRYJlWDUZZCtQiysZv3gZ0nH8G0atitDVi1l4Iav2Diw41aDj7UUg5ZhEOmEamvpt2k9jXn5bByQfJ548uxfh9ColukMKqykLSqsslWIhdWakXskHcO3ecRYWLWUQRvWsHjCAaHFFZEssmIFedZSOnAIjSFseFx14JEhNApKtZaVBKDCKlu1U1i1tZhma+O9dafCKIS+AL4BDNrwVdKxRCRLefkrrPlVg+b8n4t37I+SiqMeKwlCc6yyVYg9Vsc/9wDHP/dA0nEA3lv+GQADN5SEEk9EspCXv8KaX3XC7JmcMHtm0nE29x3o9l31Tw4SiUKFVTaytqWruhNNXgdY2aUbAINUWIlILCEXVmGxefms7e+217nmtn+muTXSWamwykZr18LWrTBgABQXp7s1O1jZ1RVWA8u+0urFIhJd81Bg59uTsWyQm0yv6QwSiwqrbNQJl1rwVeYXsCk/n8Lt27TkgohE10l7rKClTUPWfZnmlkhnpcIqG3XiwgpgYbceAIxYsyzNLRGRTmfrVigpoSE/n/L+u6S7NTspGeHOklb+klh0VmA2WrLEXYdUWP1fSGvAABxy9q/ZPOcheGs2I1Z/zoJJh4YWW0SywOefA7Bx4LDQdmeYdf6MUOIArPaWnxmx+vPQYkp2UWGVjT76yF1PmhRKuNIho0KJA9BnyChKx0+Bt2YzctXS0OKKSJbw8tfqEeGtn7dl34NDi7V26Bga8gvYpXQNVFdDUVFosSU7aCgw2zQ1tRRW++8fSsiJ899k4vw3Q4m1ev6bfLC1GoCRa/SNT0Ra+fBDAFaN3jO0kEOfvoehT98TSqzGgi6sHbaru/HJJ6HElOyiwirbLF8OVVUwZIi7hOCo5x/hqOcfCSXWZ88/wvNzX6a+oCsDy9bS3SuyRESAlBRWx7z0BMe89ERo8fzhQObPDy2mZA8NBWYbLymx//5xra7ekRqM4avhuzJ65RJGrFnG5+P2S3eTRKQzaGiABQsAWDUqvMIqbGtGqLCS2NRjlW0iCqvObPXIcQCMXKXhQBHxLFoEtbWw665s7dk73a2JafUor7D6+OP0NkQ6JfVYZZvIwqo8vU1pS/OZNWmeZ9VWr97950zpwJaISKZ8MWzusfrsM6ivhy5d0tsg6VTUY5VNrA194nqqrPEKq5E6ZVlEfBlSWG3r0YuygUNh+3ZYqrObZUfqscomy5dDZaWbtD50KBDOlgt//X/XhRIH4DAvVkmPXjQZw9CvVpDfUE9jgb7xieS65c+/zm7ALRt7weDw4j58+f+EF8yzeuQ4BpatdcOBe+8denzJXCqsMlzkUNaU9//D+cD8Qbvx5xAnrm/uH16G6+XFqgU2DBrOLqVrGLr2y+YeLBHJUQ0NDPdWM189alyoobfvFn7hs3rkHuz/4atuAvtZZ4UeXzKXhgKzyKhVi4FwT1MGV7BNef8/ocRa8f5/WOHF8oupMV8uCiW2iGSwxYvpVredsoFD2dKrT6ihRz1yK6MeuTXUmM1fBud2zrOvJX1UWGWR0SvdVjarRo8PNe7UV59i6qtPhRJryatPscSL9fkeXwNg/CIlJpGc569flYJlFqa9/S+mvf2vUGN+sfu+kJcH770HNTWhxpbMpsIqSxTUb2fX5Z8B8OWYcAurVPls768DsNfC9zFNjWlujYik1euvA7Bi1wlpbkgwW3r1gSlT3FmBr72W7uZIJ6LCKkvsueQjutXVsmrkOKr6DEh3cwLZMHgEZQOG0GtLFaNW6swakZzV1ARz5gDwaYj7+qXcUUe56xdfTG87pFNRYZUl9l3wFgCfTDwkzS2JgzEs9Hqt9v7s3TQ3RkTSZt482LCBjf2HsHboruluTWB/aBwJwPon/5nmlkhnosIqG1jLvp+8DWRYYQUsnOAKqwkL309zS0Qkbf7l5j99MvEQMCbNjQluxa4T2Nq9J7usXw0rV6a7OdJJaLmFLDBk3UoGlq2lulcxX47ZK/T4d190Y2ixjmgVa/Fek2nMy2fX5Z9SuK2G2u69QnstEckQfmG1b2q+GD501YMpiduUX8Di8VPY/6PX3HDg9OkpeR3JLOqxygL+MOCn+x6EzcsPPX5NUTE1RcWhxCosKqYwIta2HkWs2HUCBY2N7Ln4w1BeQ0QyyLp17ozA7t1ZMj41K643DBlFw5BRKYntT2fQPCvxqbDKAvsucMOACyZ+IyXxD3lrNoe8NTuUWMvems2yVrEWTTgQ0DwrkZzkTVrnW9+ivmthSl5itwdmsNsDM1ISe+HeLn/x0kvuDEHJeSqsMlyPLVWMXbaAxrz8lm9OITv4rdkcnMLCav6kQwE44IOX6Lq9NpTXEZEMMdvLB8cdl7KXmDrvFabOeyUlsTcOHMZXw3Z124k980xKXkMyiwqrDHfoG8+S39TIkj33Z1uPonQ3JyGrR+/Jil0n0HNLFV9/7/l0N0dEOsqaNa6wysuD449Pd2sS9trUk90Pf/5zehsinYImr2eA82Ls+5ffUM+N/5kFwH+OOqMjmxS6l6adxvT7ruaIl57gjcNO7BRnBsX6vQPcf86UDmyJSJb605+goQFOOw2GDwfWpbtFCXnnkOP44bP3wptvur0DJ01Kd5MkjdRjlcEOfO8F+m3eQMmw3TJrUb0o5k35FpW9+zGi5Av2WPpxupsjIil28V2vUHvn3e7GL36R3sYkqbZ7T/jxj90N9VrlPBVWGco0NXH08w8D8PwxZ3aKHp5kNBZ04XWvO/1bLz+R5taISKpNffVpCmu3smivKbB/as4G7FAXXeSuH30UNm5Mb1skrTQUmKEmzn+TYV+tYFPfQXxw4FEpfa0/XXZHaLGOaiPWa4efzLH/epD9PnyN3Zct4IuxE0N7XRHpRLZsYZo3jeH5Y87i1jaG3cPw4E1PpzQ+wHnvVHLpPgexz6fvwjXXwJ13pvw1pXNSYZWBeldu5EcP/QGAF7/9AxoLuqT09eq6hXcKdEEbsSqLB/DCt8/kuDkP8ZO/XMO11z3caRcMjTX/SnOvRNphLVxwAcWV5awaOY6F3nIrqdQU0jp87XnqlIsYv2guBXfd5c5yPPbYDnld6VxUWHUibU2W9uU1NvDTe66iT9Umluy5Py9P+37K23X4K38H4NUjTkk61mIv1vgYsf550nT2/uw9Rq1eyhmP3s6D512V9Gt2JE14F4ntvJlzOez1Zzj7b39je9dC/jr92g6ZxjD+7isBWHzB71L6OmtG7sEzJ/2UU/5+J5VnnMXVNzxGTe++gD7/uUSFVQcLUjzFYpoaOW3WHey55EMq+vTn3vNvoCk/9f+Ekz94CQinsPrSixWrsGos6MJfpl/H1dedzTfeeo6ygUOZfcK5GT+HTERg3JIP+cHDtwDw0Dm/Ye2w3TrkdQ/+1C0+vLgDXuv5Y85kn0/fYdzSj7n4z7/gzktuprp3vw54ZeksNHk9QxRvLuPyP17MtJeeoDEvn3vP/x1VfQaku1kpsW7Yrjxy5i9oMoaT/nEv59/9G7rVbk13s0QkUY2NMGMGP7/5Iro01PHa1JN4/6Cj092qlLB5+dz/k2vZ1HcQY7/4hKuuP4eRq5amu1nSgQIVVsaYo40xS40xXxhjrojyeDdjzOPe4+8bY0aH3dBctcu6lXx/1h1cf9UZjF/yIZW9+3H7z/7E53vul+6mpdRbh32HP//XLWwr7MmUuS9z0y++y3efvod+5evT3TTJMMpfaVRVBffc4876u+oq8mwT/zruHB45M7OXV2hP+YAhzLh6Jst325v+5eu56rqz4ZRT4I033LpdktXaHUcyxuQDdwJHAiXAXGPMs9baRRGHnQdsttbubow5HbgJOC2UFpaVwerVoYRKxPXPLYz52NUnTIj+gLXR77OWMcsXYrAYa8lraiS/oYGCxga6bq+lW902iqo306einMHrVzNq1RL6bd7QHOKzCQdy//+7Jmt7qlr7ZNKh/O6qBzj3/uvZdcVCTnjuAU547gHK++/Cl2P2omzgUCqKB7GlZxG1hT2o61pIY34BjQUFNJl8mvLywBgsbhjRth5N7ODhxeuvjT0QEetvKaG/v442fjz06JHuVkSV9vzV1AQfZ+C6bK1zmH/by2M0NblLfb27bNsGW7ZARYXbVHn1avjoI1i0yPVWARV9+vPAeVezcJ+DOvjNpEdl8QBu/tXdfH/W//DN1/9B3lNPwVNPUVvYg8JDDoI993SLog4aBEVF0LMndO3qLvn57mLMjpfWNEUieYMHe4vThifIBJ0DgC+stSsAjDGzgBOByMR0InCt9/Pfgf81xhhro1UYcXrmGZg+Pekwibq6rQeviz/eb+M8fmv3Xsw9YBqvf/O7rBqzV/wvmOHWDR3D7377ALsvW8CR/5nFXgvfp3/5evpnW89VjL+lsP/+UmL+fJjYaZfGSG/+qquDyZOTDpOx8vNh6lTuHfctPtp/Kg1duqW7RR2qoUs3Hj3rF8w5/hwOf+XvHPj+CwwsWwsvv+wukn6XXw633BJqSNNe7jDGnAIcba39iXf7LOBAa+3FEcd85h1T4t1e7h2zsVWs6YBfJY0D4hl4HgBk+6prufAeQe8zm8T7HkdZawemqjGtKX91uFx4n7nwHkHvM5pA+StIj1W0vsbW1ViQY7DW3gfcF+A1d26EMfOstVn91S8X3iPofWaTDHiPyl8dKBfeZy68R9D7TEaQyeslwIiI28OBtbGOMcYUAH2ATWE0UEQkCcpfItKhghRWc4GxxpgxxpiuwOnAs62OeRY42/v5FOCVUOYniIgkR/lLRDpUu0OB1toGY8zFwAtAPvCAtXahMeZ6YJ619lngfuBvxpgvcN/0Tk9BWxPqgs8wufAeQe8zm3Tq96j81eFy4X3mwnsEvc+EtTt5XURERESC0crrIiIiIiFRYSUiIiISkowprIwx1xpjvjLGzPcux6a7TWFqb9uNbGGMWWmM+dT7N5yX7vaEwRjzgDFmg7cekn9fP2PMf4wxy7zrvulsYxhivM+s/lyGKZt/V8pfmS0XclhH5q+MKaw8t1trJ3mXOeluTFgitt04BtgLOMMYk83LrB/u/RtmyxopM4HWO8peAbxsrR0LvOzdznQz2fl9QpZ+LlMk635Xyl9ZYSbZn8Nm0kH5K9MKq2zVvO2GtbYO8LfdkAxgrX2Dndc9OhF4yPv5IeC7HdqoFIjxPkWUvzJcLuSwjsxfmVZYXWyM+cTr0svobslWhgFrIm6XePdlIwu8aIz50NsiJFsNttauA/CuB6W5PamUrZ/LVMjG35XyV3bKlRwW+meyUxVWxpiXjDGfRbmcCNwN7AZMAtYBt6a1seEKtKVGljjEWrsfbtjgImPMYRIKNgsAAAFGSURBVOlukCQlmz+XccvRHKb8JZkqJZ/JIHsFdhhr7bQgxxlj/gLMTnFzOlKQbTeygrV2rXe9wRjzD9wwwhvpbVVKlBpjhlhr1xljhgAb0t2gVLDWlvo/Z+HnMm45msOUv7JT1uewVOWvTtVj1RbvH9Z3EvBZrGMzUJBtNzKeMaanMabI/xk4iuz6d4wUuU3K2cA/09iWlMnyz2Wosvh3pfyVnbI+h6XqM9mpeqzacbMxZhKui3kl8NP0Nic8sbbdSHOzUmEw8A9jDLi/vUettc+nt0nJM8Y8BkwFBhhjSoBrgBuBJ4wx5wGrgVPT18JwxHifU7P1c5kCWZnDlL8yXy7ksI7MX9rSRkRERCQkGTMUKCIiItLZqbASERERCYkKKxEREZGQqLASERERCYkKKxEREZGQqLASERERCYkKKxEREZGQ/H9BJOp/QgsMFQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Distribution of guide after optimisation\n",
    "new_speed_guide = [guide(guess=5., guess_sd=2., time=10.).item() for i in range(1000)]\n",
    "plot_hist(speed_guide, new_speed_guide, \"Before optimisation\", \"After optimisation\",\n",
    "          true_dist=(3.4,0.89))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is that after optimisation, the distribution of the guide looks almost identical to the exact posterior. So we have a good approximation to the posterior distribution which we can easily draw samples from.\n",
    "\n",
    "This example used stochastic variational inference, which gives us a way to approximate a posterior distribution when we don't know the exact form of the distribution. In this case, we used normal distributions, so we could have used Markov Chain Monte Carlo to learn a posterior distribution exactly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, we can specify the number of samples and number of steps we want to run this optimisation over, and the number of samples for the Monte Carlo posterior approximation. Then we can obtain the marginal distribution for the latent variable in fewer lines of code, e.g."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "svi = pyro.infer.SVI(model=conditioned_propagate, \n",
    "                     guide=guide,\n",
    "                     optim=optimizer,\n",
    "                     loss=pyro.infer.Trace_ELBO(),\n",
    "                     num_samples=1000, \n",
    "                     num_steps=1000)\n",
    "svi.run(guess=5., guess_sd=2., time=10.)\n",
    "marginal = pyro.infer.EmpiricalMarginal(svi, sites=\"speed\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Draw samples from marginal\n",
    "marginal()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Distribution of guide after optimisation\n",
    "new_speed_guide = [marginal().item() for i in range(1000)]\n",
    "plot_hist(speed_guide, new_speed_guide, \"Before optimisation\", \"After optimisation\",\n",
    "          true_dist=(3.4,0.89))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MCMC\n",
    "Here we will carry out MCMC to show that we can do inference exactly. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We set up the mcmc object below, with a Hamilton Monte Carlo (HMC) kernel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "import pyro.infer.mcmc\n",
    "hmc_kernel = pyro.infer.mcmc.HMC(conditioned_propagate,\n",
    "                                   adapt_step_size=True)\n",
    "mcmc = pyro.infer.mcmc.MCMC(hmc_kernel,\n",
    "                            num_samples=1000,\n",
    "                            warmup_steps=500)\n",
    "mcmc_run = mcmc.run(guess=torch.tensor([5.]),\n",
    "                    guess_sd=torch.tensor([2.]), \n",
    "                    time=torch.tensor([10.]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now `mcmc_run` is an object from which we can get the marginal distribution for any latent variables, `p(z)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "marginal = pyro.infer.EmpiricalMarginal(mcmc_run, sites=\"speed\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mcmc_speed = [marginal().item() for i in range(10000)]\n",
    "plot_hist(new_speed_guide, mcmc_speed, \"SVI\", \"HMC MCMC\",\n",
    "          true_dist=(3.4,0.89))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:infolab] *",
   "language": "python",
   "name": "conda-env-infolab-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}