{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# (DEPRECATED) An Introduction to Inference in Pyro\n", "\n", "## WARNING\n", "\n", "***This tutorial has been deprecated*** in favor of the updated [Introduction to Pyro](https://pyro.ai/examples/intro_long.html). It may be removed in the future.\n", "\n", "Much of modern machine learning can be cast as approximate inference and expressed succinctly in a language like Pyro. To motivate the rest of this tutorial, let's build a generative model for a simple physical problem so that we can use Pyro's inference machinery to solve it. However, we will first import the required modules for this tutorial:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import torch\n", "\n", "import pyro\n", "import pyro.infer\n", "import pyro.optim\n", "import pyro.distributions as dist\n", "\n", "pyro.set_rng_seed(101)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## A Simple Example\n", "Suppose we are trying to figure out how much something weighs, but the scale we're using is unreliable and gives slightly different answers every time we weigh the same object. We could try to compensate for this variability by integrating the noisy measurement information with a guess based on some prior knowledge about the object, like its density or material properties. The following model encodes this process:\n", "\n", "$${\\sf weight} \\, | \\, {\\sf guess} \\sim \\cal {\\sf Normal}({\\sf guess}, 1) $$\n", "$${\\sf measurement} \\, | \\, {\\sf guess}, {\\sf weight} \\sim {\\sf Normal}({\\sf weight}, 0.75)$$\n", "\n", "Note that this is a model not only for our belief over weight, but also for the result of taking a measurement of it. The model corresponds to the following stochastic function: " ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "def scale(guess):\n", " weight = pyro.sample(\"weight\", dist.Normal(guess, 1.0))\n", " return pyro.sample(\"measurement\", dist.Normal(weight, 0.75))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Conditioning\n", "\n", "The real utility of probabilistic programming is in the ability to condition generative models on observed data and infer the latent factors that might have produced that data. In Pyro, we separate the expression of conditioning from its evaluation via inference, making it possible to write a model once and condition it on many different observations. Pyro supports constraining a model's internal `sample` statements to be equal to a given set of observations.\n", "\n", "Consider `scale` once again. Suppose we want to sample from the distribution of `weight` given input `guess = 8.5`, but now we have observed that `measurement == 9.5`. That is, we wish to *infer* the distribution:\n", "$$({\\sf weight} \\, | \\, {\\sf guess}, {\\sf measurement} = 9.5) \\sim \\, ? $$\n", "\n", "Pyro provides the function `pyro.condition` to allow us to constrain the values of sample statements. `pyro.condition` is a higher-order function that takes a model and a dictionary of observations and returns a new model that has the same input and output signatures but always uses the given values at observed `sample` statements:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "conditioned_scale = pyro.condition(scale, data={\"measurement\": torch.tensor(9.5)})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Because it behaves just like an ordinary Python function, conditioning can be deferred or parametrized with Python's `lambda` or `def`:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def deferred_conditioned_scale(measurement, guess):\n", " return pyro.condition(scale, data={\"measurement\": measurement})(guess)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In some cases it might be more convenient to pass observations directly to individual `pyro.sample` statements instead of using `pyro.condition`. The optional `obs` keyword argument is reserved by `pyro.sample` for that purpose:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def scale_obs(guess): # equivalent to conditioned_scale above\n", " weight = pyro.sample(\"weight\", dist.Normal(guess, 1.))\n", " # here we condition on measurement == 9.5\n", " return pyro.sample(\"measurement\", dist.Normal(weight, 0.75), obs=torch.tensor(9.5))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, in addition to `pyro.condition` for incorporating observations, Pyro also contains `pyro.do`, an implementation of Pearl's `do`-operator used for causal inference with an identical interface to `pyro.condition`. `condition` and `do` can be mixed and composed freely, making Pyro a powerful tool for model-based causal inference." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Flexible Approximate Inference With Guide Functions\n", "\n", "Let's return to `conditioned_scale`. Now that we have conditioned on an observation of `measurement`, we can use Pyro's approximate inference algorithms to estimate the distribution over `weight` given `guess` and `measurement == data`. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "Inference algorithms in Pyro, such as `pyro.infer.SVI`, allow us to use arbitrary stochastic functions, which we will call *guide functions* or *guides*, as approximate posterior distributions. Guide functions must satisfy these two criteria to be valid approximations for a particular model: \n", "1. all unobserved (i.e., not conditioned) sample statements that appear in the model appear in the guide.\n", "2. the guide has the same input signature as the model (i.e., takes the same arguments)\n", "\n", "Guide functions can serve as programmable, data-dependent proposal distributions for importance sampling, rejection sampling, sequential Monte Carlo, MCMC, and independent Metropolis-Hastings, and as variational distributions or inference networks for stochastic variational inference. Currently, importance sampling, MCMC, and stochastic variational inference are implemented in Pyro, and we plan to add other algorithms in the future.\n", "\n", "Although the precise meaning of the guide is different across different inference algorithms, the guide function should generally be chosen so that, in principle, it is flexible enough to closely approximate the distribution over all unobserved `sample` statements in the model. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the case of `scale`, it turns out that the true posterior distribution over `weight` given `guess` and `measurement` is actually ${\\sf Normal}(9.14, 0.6)$. As the model is quite simple, we are able to determine our posterior distribution of interest analytically (for derivation, see for example Section 3.4 of [these notes](http://www.stat.cmu.edu/~brian/463-663/week09/Chapter%2003.pdf)).\n", "\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def perfect_guide(guess):\n", " loc = (0.75**2 * guess + 9.5) / (1 + 0.75**2) # 9.14\n", " scale = np.sqrt(0.75**2 / (1 + 0.75**2)) # 0.6\n", " return pyro.sample(\"weight\", dist.Normal(loc, scale))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Parametrized Stochastic Functions and Variational Inference\n", "\n", "Although we could write out the exact posterior distribution for `scale`, in general it is intractable to specify a guide that is a good approximation to the posterior distribution of an arbitrary conditioned stochastic function. In fact, stochastic functions for which we can determine the true posterior exactly are the exception rather than the rule. For example, even a version of our `scale` example with a nonlinear function in the middle may be intractable:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "def intractable_scale(guess):\n", " weight = pyro.sample(\"weight\", dist.Normal(guess, 1.0))\n", " return pyro.sample(\"measurement\", dist.Normal(some_nonlinear_function(weight), 0.75))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What we can do instead is use the top-level function `pyro.param` to specify a *family* of guides indexed by named parameters, and search for the member of that family that is the best approximation according to some loss function. This approach to approximate posterior inference is called *variational inference*.\n", "\n", "`pyro.param` is a frontend for Pyro's key-value *parameter store*, which is described in more detail in the documentation. Like `pyro.sample`, `pyro.param` is always called with a name as its first argument. The first time `pyro.param` is called with a particular name, it stores its argument in the parameter store and then returns that value. After that, when it is called with that name, it returns the value from the parameter store regardless of any other arguments. It is similar to `simple_param_store.setdefault` here, but with some additional tracking and management functionality.\n", "\n", "```python\n", "simple_param_store = {}\n", "a = simple_param_store.setdefault(\"a\", torch.randn(1))\n", "```\n", "\n", "For example, we can parametrize `a` and `b` in `scale_posterior_guide` instead of specifying them by hand:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "def scale_parametrized_guide(guess):\n", " a = pyro.param(\"a\", torch.tensor(guess))\n", " b = pyro.param(\"b\", torch.tensor(1.))\n", " return pyro.sample(\"weight\", dist.Normal(a, torch.abs(b)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As an aside, note that in `scale_parametrized_guide`, we had to apply `torch.abs` to parameter `b` because the standard deviation of a normal distribution has to be positive; similar restrictions also apply to parameters of many other distributions. The PyTorch distributions library, which Pyro is built on, includes a [constraints module](https://pytorch.org/docs/master/distributions.html#module-torch.distributions.constraints) for enforcing such restrictions, and applying constraints to Pyro parameters is as easy as passing the relevant `constraint` object to `pyro.param`:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "from torch.distributions import constraints\n", "\n", "def scale_parametrized_guide_constrained(guess):\n", " a = pyro.param(\"a\", torch.tensor(guess))\n", " b = pyro.param(\"b\", torch.tensor(1.), constraint=constraints.positive)\n", " return pyro.sample(\"weight\", dist.Normal(a, b)) # no more torch.abs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pyro is built to enable *stochastic variational inference*, a powerful and widely applicable class of variational inference algorithms with three key characteristics: \n", "\n", "1. Parameters are always real-valued tensors\n", "2. We compute Monte Carlo estimates of a loss function from samples of execution histories of the model and guide\n", "3. We use stochastic gradient descent to search for the optimal parameters. \n", "\n", "Combining stochastic gradient descent with PyTorch's GPU-accelerated tensor math and automatic differentiation allows us to scale variational inference to very high-dimensional parameter spaces and massive datasets. \n", "\n", "Pyro's SVI functionality is described in detail in the [SVI tutorial](svi_part_i.ipynb). Here is a very simple example applying it to `scale`:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "a = 9.11483097076416\n", "b = 0.6279532313346863\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAkZklEQVR4nO3dd3hc1Z3/8fdXxb3bwrgAwgYMBkwTxEAophpwYpYQyiZZssuuQ8qG7C7FoSRkQ1gSSiBA2HUSEiCsSZaaXzDdBlPtyNi44G4M7pKbLBfV+f7+mJGskUbSaKSZke79vJ7Hj2bu3HvPuRr5M2fOPfdcc3dERCQ8crJdARERySwFv4hIyCj4RURCRsEvIhIyCn4RkZBR8IuIhIyCX0QkZBT8IoCZrTWzfWa2u8G/h83sm2b2bjPbvGVmFbF1y8xstpkd22idSWY218z2mNk2M3vKzEZm5qhEElPwi+z3JXfv0+Df95LY5nvu3gcYDLwFPFn3gpldDvwv8CAwBDgaqATeNbOBHV57kSQp+EU6gLvXAE8DYwHMzID7gDvd/Sl33+fum4F/BnYD/5a1ykroKfhFOoCZdQO+BnwYWzQGOBj4v4bruXsEeBY4P6MVFGkgL9sVEOlEXjCzmgbPbwSqW9nmV2Z2L9AL2AdcFls+JPZzU4JtNjV4XSTj1OIX2e9Sdx/Q4N9vktjm++4+AOgBTAKeMbNxwNbY68MSbDOswesiGafgF+kA7h5x93eAVcAFwHJgPfDVhuuZWQ7wFeDNjFdSJEZdPSKtMzPr0XCBu1ckWOlUoid3l7i7m9kNwG/MbD3wPNAfuAvoB/wy/dUWScw0H79IdBw/MBSobbD4deBF4PcJNskH3gDGA3XnBTYDj7h7faib2WTgNqIfCJXAq8BN7r6ugw9BJGkKfhGRkFEfv4hIyCj4RURCJm3Bb2aPmVmJmS1utPxfzWy5mS0xs1+kq3wREUksnS3+PwATGy4wswnAZGCcux8N3JvG8kVEJIG0Ded099lmVtho8beBu929MrZOSTL7GjJkiBcWNt6ViIi0ZN68eVvdvaDx8kyP4z8COMPMfgZUADe4+99a26iwsJDi4uK0V05EJEjM7LNEyzMd/HnAQKJjn08G/mxmozzBmFIzmwJMATj44IMzWkkRkSDL9Kie9cBzHjUXiNDMZFXuPs3di9y9qKCgyTcVERFJUaaD/wXgHAAzOwLohiarEhHJqLR19ZjZdOBsYEhsrpIfA48Bj8WGeFYB1yTq5hERkfRJ56ieq5t56evpKlNERFqnK3dFREJGwS8iEjKBD/6S8gpeW7I529UQEek0Ah/8V0/7kClPzqOqJpLtqoiIdAqBD/512/cB4GjwkIgIhCD4RUQknoJfRCRkFPwiIiET+OBX376ISLzAB7+IiMQLfPAblu0qiIh0KoEPfhERiafgFxEJmcAHv07uiojEC3zw11Ffv4hIVGiCXy1/EZGowAe/WvoiIvECH/wiIhIvbcFvZo+ZWUns/rqNX7vBzNzMhqSrfBERSSydLf4/ABMbLzSzg4Dzgc/TWHY99e2LiMRLW/C7+2xge4KXfgncBJlNZPX1i4hEZbSP38y+DGxw94+TWHeKmRWbWXFpaWkGaiciEg4ZC34z6wXcCvwomfXdfZq7F7l7UUFBQbvLV5ePiEhUJlv8o4FDgY/NbC0wEvjIzA5MZ6Hq4hERiZeXqYLcfRFwQN3zWPgXufvWtJarlr6ISJx0DuecDnwAjDGz9WZ2bbrKEhGR5KWtxe/uV7fyemG6yhYRkebpyl0RkZAJTfC7uvpFRIAQBb+IiEQp+EVEQkbBLyISMgp+EZGQUfCLiISMgl9EJGQU/CIiIaPgFxEJGQW/iEjIKPhFREJGwS8iEjKBD37N0SMiEi/wwV9HHwAiIlGBD37TnRdFROIEPvhFRCSegl9EJGTSec/dx8ysxMwWN1h2j5ktM7OFZva8mQ1IV/l11LcvIhIvnS3+PwATGy17HTjG3ccBK4AfprF8ERFJIG3B7+6zge2Nlr3m7jWxpx8CI9NVfpP6oKa/iAhkt4//n4CXm3vRzKaYWbGZFZeWlqZciEb1iIjEy0rwm9mtQA3wVHPruPs0dy9y96KCgoLMVU5EJODyMl2gmV0DTALOddepVxGRTMto8JvZROBm4Cx335uJMvXRIiISL53DOacDHwBjzGy9mV0LPAz0BV43swVm9t/pKl9ERBJLW4vf3a9OsPh36SqvNWr5i4hEBf7KXY3qERGJF/jgFxGReIEPfnXxiIjEC3zwi4hIPAW/iEjIKPhFREImNMGvrn4RkajQBL+IiEQp+EVEQkbBLyISMgp+EZGQUfCLiIRMaIJfU/+LiESFJvhFRCRKwS8iEjIKfhGRkFHwi4iEjIJfRCRk0nnP3cfMrMTMFjdYNsjMXjezlbGfA9NVfmMa0yMiEpXOFv8fgImNlk0F3nT3w4E3Y89FRCSD0hb87j4b2N5o8WTg8djjx4FL01W+iIgkluk+/qHuvgkg9vOA5lY0sylmVmxmxaWlpRmroIhI0HXak7vuPs3di9y9qKCgINvVEREJjEwH/xYzGwYQ+1mS7gJ1UldEJF6mg/8vwDWxx9cAL2aqYE3VIyISlc7hnNOBD4AxZrbezK4F7gbON7OVwPmx52ll6S5ARKSLyUvXjt396mZeOjddZYqISOs67cldERFJDwW/iEjIBD74dU5XRCRe4INfRETihSf41fQXEQFCEPwazikiEi/wwS8iIvECH/zq4RERiRf44BcRkXgKfhGRkAlN8Ls6fUREgBAEv0b1iIjEC3zwq50vIhIv+MEfm4i/vKImyzUREekcAh/8kViTf+pzC7NbERGRTiKp4Dez682sn0X9zsw+MrML0l25jrS1vCrbVRAR6RSSbfH/k7vvAi4ACoB/JAN3z+pIGtUjIhKVbPDXDY65GPi9u3+MBsyIiHRJyQb/PDN7jWjwv2pmfYFIqoWa2b+Z2RIzW2xm082sR6r7EhGRtkk2+K8FpgInu/teIJ9od0+bmdkI4PtAkbsfA+QCV6WyLxERabtkg/9UYLm77zSzrwO3AWXtKDcP6GlmeUAvYGM79iUiIm2QbPA/Cuw1s+OAm4DPgCdSKdDdNwD3Ap8Dm4Ayd3+t8XpmNsXMis2suLS0NJWiREQkgWSDv8ajV0JNBh509weBvqkUaGYDY/s5FBgO9I59i4jj7tPcvcjdiwoKClIpqtH+2r0LEZFASDb4y83sh8A3gJfMLJdoP38qzgM+dfdSd68GngNOS3FfIiLSRskG/5VAJdHx/JuBEcA9KZb5OTDezHqZmQHnAktT3JeIiLRRUsEfC/ungP5mNgmocPdU+/jnAM8AHwGLYnWYlsq+2lRuugsQEekikp2y4QpgLvBV4Apgjpldnmqh7v5jdz/S3Y9x92+4e2Wq+xIRkbbJS3K9W4mO4S8BMLMC4A2iLfcuQZcZi4hEJdvHn1MX+jHb2rBtp6CuHhGRqGRb/K+Y2avA9NjzK4EZ6alSerjGc4qIAEkGv7vfaGZfAU4n2msyzd2fT2vNREQkLZJt8ePuzwLPprEuIiKSAS0Gv5mVk7h73AB3935pqVUH2bBzX7arICLS6bQY/O6e0rQMncWOPbrrlohIY11qZI6IiLRfaIJfY3pERKJCE/wdbV9VLeUV1dmuhohImyn4U3T6z2dy7B1NbiMgItLpKfhTtF0njkWki1Lwi4iEjIJfRCRkFPwiIiGj4BcRCZnwBL8G8ouIAGEKfhERAbIU/GY2wMyeMbNlZrbUzE7NRj1ERMIo6WmZO9iDwCvufrmZdQN6ZakeIiKhk/HgN7N+wJnANwHcvQrQ1VAiIhmSja6eUUAp8Hszm29mvzWz3o1XMrMpZlZsZsWlpaXtLlTndkVEorIR/HnAicCj7n4CsAeY2ngld5/m7kXuXlRQUJDpOoqIBFY2gn89sN7d58SeP0P0gyDjtuyqoHDqS3z0+Y5sFC8ikhUZD3533wysM7MxsUXnAp+ko6wP12xrWC5le6t5c+mW+mXvrdoKwBPvr01H8SIinVK2xvH/K/CUmS0EjgfuSkchMxZtint+3R/nce3jxZSWVwJgFl2u/n8RCZOsDOd09wVAUbrLsbpkjz1eu20PANW1kegyLOF2IiJBFugrdxvGunvz7foWXhIRCZxAB39OgxZ/omz/37mfA7C7siZDNRIRyb5AB3/DJn+iTp25n24HYHNZRWbqIyLSCQQ6+OO6epp5nOi5iEiQBTr4G3b1tKSl/n8RkaAJdPBbK109dSItBP/W3ZUcdssM5n2mi7xEJBgCHfytndytf63BixXVtfXj/AHmrNlOTcT53btr0lBDEZHMC3TwN2zxt9Sbs7JkN4ffOoMde6q45rG5nPyzN9JfORGRLAl08DfWUndPda0zf90O5sRG+oiIBFWggz/Zk7siImES6OBvmPufb99LRIN3RESCHfyNW/ybd0Uv1Dr97pls2dX0oi2N6hSRMAh48Df/2sL1ZW3alz4URCQoAh38LUn0maBTAiISBoEOfmtjkqtVLyJhEOjgb6mrJ6eVI9c0DiISVIEO/pZutNLaTVga5766gUQkKIId/C2FdStB3tL8PSIiXVnWgt/Mcs1svpn9NX1ltPBaK9tqzL+IBFU2W/zXA0vTWUBLJ3cXb2h5OGfjFr++AIhIUGQl+M1sJHAJ8Nu0ltPCa/e+tqLJMoW7iIRBtlr8DwA3AZHmVjCzKWZWbGbFpaWlKRXS2nDOb501qtnX1McvIkGV8eA3s0lAibvPa2k9d5/m7kXuXlRQUJBSWS0N5wTo2z2vUd32P062j//tFaUUTn2JnXur2lg7EZHsyEaL/3Tgy2a2FngaOMfM/piOglqbnbPxN4KGjfzmWvwPvLGCwqkv1T//9axVAHyyaVeKtRQRyayMB7+7/9DdR7p7IXAVMNPdv56Ostoz9N6b6YR64I2V8evVl6WB/iLSNQR8HH/qYewt3qxxv7mxG7foAi8R6SryWl8lfdz9LeCtdO2/rWGc7K0aE9lc1nSaZxGRzijYLf42rt8w7Ns6pucHf1rQxi1ERLIj0MHfnlsv1k3S1lFdODW1EX704mJ9MxCRrAt08LcW2i3NwOn163RMXd5bvY0nPviMqc8tbPe+fvvOGmYtK+mAWolIGGW1jz/dWgv+FxZsbPa1jrx+67Nte1hVshvomDmA7nwpOtPF2rsvaf/ORCR0Ah78LSd/XRin21n3vJWW/dbURti6u4oD+/dIy/5FJJiC3dXTxvU/XLOt/nGywzkb+uFzi4hkcFrPn/71E8b/15u6alhE2iTYwd/G5P/tu5/uf5JCfk+f+zlrtu5Jev0nP/yMNaWpf+uYuTzaz79rX03K+xCR8Al08LdnVM/dLy9Ledttuyv5tJUPAHfn9hcWM/mR91Iup35fqXxKiUhoBTr42zMS87n5G6L7aPNOnLPueYsJ977V4lp1PULlFam31us+2DSRqIi0RbCDv52D8EvKK/jOUx+1ebvdlfvDfFdFdcJ1OuJm7nVHp9wXkbYIePC3b/slG/fPuFm2r5oH3mh685bG9lTWxj0fd8drCdfryHPAHfEhIiLhEezhnO2dMbNBnr6/ehvvr97W/Lox3396fpK7bhrWSzaWMaBXN0YM6JnUPuq7epJaW0QkKtAt/tZuxNKa+et2tnmbz7btTWq9RI30S371LqffPbP++cL1O3lm3vrmd2LN70tEpDmBbvHntDP5f/XmytZXSlEyYf3lh6Mjfi4/aWRre2t/hUQkNALd4u/MOuKevnUfa6XluoBLRJIX6ODvzDdH6Yg2et2opat/82EH7E1EwiLYXT1ZTv512xP39x/741c5aFCvuPV6d2/7W5GpuYZEJFgCHfzZdsYvZiVcXl5ZE3dz9qumfcj5Y4e2ur9Xl2xmYK9uCV+rjTi1EadbXqC/xIlIB8h48JvZQcATwIFABJjm7g+mo6z2jupJh0Rj7jeV7eMP769tdps3PtnCuUcdwLeenNfsOqNvmQFoqmYRaV02Wvw1wH+4+0dm1heYZ2avu/snHV1Qtrt6Enln5dYmyxpfzFW2t5rNu/bfqeufnyjmtNGDk9q/u7f7imURCbaM9wu4+yZ3/yj2uBxYCoxIR1ldNf6u+J8PuPCB2XHLkrl4DGDppnIANu7cx49fXExNbSTpches28m+qtrWVxS27q7k5mcWUlGt35d0PVntEDazQuAEYE6C16aYWbGZFZeWlqZaQLvqly3Lt5SnvG1txJn4wGxOu3smj3/wGdP/tq7JOuUV1eytqmHxhrL6ZSXlFVz6yHvc+MzHlFdUs7TBOYiwqqiupbqZD867X17Gn4rX8deFmzJcK0nGii3l/G3t9mxXo9PKWvCbWR/gWeAH7t4kZdx9mrsXuXtRQUFB5ivYRW3dU8myzfs/OG5/YXHc66tKyjn2jtcY+6NXmfTQu/XhX1EVDbgF63ZyzWNzuejBdwI3B9BzH62ncOpL9a30b/5+Lsfe8Wqz6x95+ytc/t8fJHytI67DkOTM+2wHo2+ZQWl5ZVLr//lv67jgl7P5ajPvnWQp+M0sn2joP+Xuz6WrnM54cjfdfvKXJU2Wbd8TvcBrycYyftfwZjPApIfeZfnmcnbE7uJVG3E++nwnABXVkTZ1FXVW23ZXsm13Jfe/Hp1kr2RXNEDeWl7a6rTYH7dx2o7FG8oonPoS769uei5HUvPbd9ZQG3HmfJpcd+dNzy5Mc426vowHv0XPPP4OWOru96e1rC7by5+6tQnmCirbF50a+pJfvcv0uU27fi58YHb9DWE2le0/qXzUj17hgl/ObrI+RLuLtjQ4Ad3Y7soa5n22g9eWbGbxhjJWbilnV0U1c9Yk9583Gdc9OY83l25psrymNsI5973Fq0s2A3DSnW9w0p1vdFi5DVXW1MZ9M6oL/DeXliRcf/GGMp78YG1a6tKcsn3V3P/a8i77Id5Fe2w7tWyM6jkd+AawyMwWxJbd4u4zOrog/cFEVVTXcuJPX09p2zVb93DD/33MwYN6UbavmilnjmJovx6Mv+tN9lTVMvvGCQzsnc/DM1cxadxwjh3Zn/mf7+C7T33ExrL4D4YvHDqIOZ9uZ9EdFxBx6N8zP+Vjqo04ryzZzCtLNjcZwlq2r5o1pXuY+uxCLjz6wPrlVTX7g2/r7v3dBpGIk5NjlO2tZk9VDcMbzY66duseSndXUji4NwV9u+PuVMb2devzi7n1+cX1dahrbCT601u4fmf9/EvfOLUw5WNvq7tfXsb0uZ8z5sB+XDJuWMbKraiuZdayEi46tn1l1v1O6z5f5366nYMG9WRY//3vk7vz0MxVXHZiy+NE3lpeQnlFDV86bni76pQuNbUR3lm1leK127nhgjFpG6GX8eB393fJ0IAb5X7URQ++067tG84Q2rir6Mx79l+k9j+z13DsiP4sanDSuKE5n0ZPtl39mw9ZvGEXP5o0lvPHDq2/irmyppbZK7Zy6ujB5OUYPfJzm+zjxQUb6JmfS37u/i+rhVNf4ttnj+aUwkGMG9m/fjqMHXurufOv+0cJl8T6iGsiEYru3F/v7z89nwevOoHT7o5+mD33ndM4oG/3+tfPbnA3tbV3X8Lj76/lpUYndR+ZtYrvTjis/vne6lr+5Ylibr34KFaX7ubMIwrqQx+gujbCsk3lDO3fnZJdlRw9vB9jbn+Ff51wGN875zCWbNzFnE+3M2tZCV8tGslTcz5n2aZdLLzjwoS/2zruzrTZa/jSccMZPqAnG3buY/7nO+p/v7UR58E3VvC18Yewr6qWkvJKTi4cmJaAuezX7/PJpl08++3TOOmQga2u//onW+jVLZfTDxtSvywScfZUxXfHXfE/0b77+684jstOjE5guLJkN/e/voLXP4n/BjjliWL+67Jj6dsjn255OXzz938D4AujBnFA3x5x624uq2DnvipGDenDjr1VFPTpTk6O1Q+RrqypJT8n+ndXuruSt1eUUra3mleWbGbCmAKOGtaPd1ZuJeLOl44bTt8eeRx+QF9yLP6mUFU1ETbu3MeIgT3JMeOhmSv55mmF9MjP5dz73mbDzn0AfPm4EYw5sG9Sv+u2sq5wAq+oqMiLi4vbvN3DM1dy72ut3zxFuoYcS+4GNj+99JgmJ7U7ypPXnsIdf1nC6tKW76nc2LVfPLTJh2YqTjpkIH9/ysH07JbL1t2VPDNvPQvXRz9ob554JEcM7cO1j7f9/0pL8nON6lrn1FGDqa6Nhlbdt7lp3ziJWctLmD53HUce2Jce+bmcPaaAB97YP7Nt7265nHPUUAb37kZlTYS8HGNwn25s2VXB9Lnr+MKhgzjvqKH8bMbS+m0uGTesyYcrwMvXn5FyQ6buG2dXctgBfXjj389KeXszm+fuRU2WBzn4H5m1inteXZ6GGomIZMbquy4mN8WRKs0FvyZ2EelkzjnygBZfP3ZE/wzVJGpIn+7NvtaWHqIjhvapf9yrW9NuvKBo7v0rHNwr4fLWvLSo468VCfQkbTq5mxmHDO7F/Vccx1cejfa93jzxSGYtK2Fg73wM44SDB7B22x6+eFgB40cN4r3V2zi5cGD9ybm/+/V7zI8NIQX4+viDGdirGw/NXAVE+9Wfn7+eNaV7uPykkeTl5rBkQxnjRw9mx54quuflMqBXPt1jE9SVlleyeVcFB/brQb+e+WzcuY991bVUVEc4Ymgf+vbIrx+JYxbtw925t5rKmgjb9lRSWRPhs217GNCrG0cM7UtVTYQNO/Yx5sC+mMEL8zcwf91OppwxiuMOGpD076k24pzx85lNTno39tt/KGLOp9s5dfRg9lbV0Ktb+v6bfrhmG3265zHpoXeB/f3mkYgz6pYZDOrdjeLbzktqX/uqasnPNfJy29eefPSt1Sxcv5MLjh7KyIG9GD6gZ/3tSPdW1bCnspYde6vo3zOfPZU1RDx6zuSTjbso6Nudx99fy+QTRvDl44bz6dY9HDSwJ9/93494dcn+/v+Z/3EW59z3dov1+JczDuU377S9e+7kwkHMXNZ0VNfk40fwYOzmTv/99RNZtrmcTzbu4udfGUfPbrmUlldStq+aMQf2JT83h4rqWt5btbXVhkAqAt3V8+hbq/n5K8vSUKOurW+PvCbj1x+86niuf3pBUtvfPPFILjtxBD/5f0uYsWgz79w0gYMG9aqfJvqwA/q0soemCqe+VP/47RvP5pDBvTnpp69zybhh/OfkY9q8v87qhfkb+MGfFvDAlcfz2ba9jCrozfABPeiZn8fFv4r2XWdjor27X17Gq0s2M+uGs+uX7a6sIcdI6wdPpn24Zhv9euQzdng/qmoifL59L+fdH/0A6Ns9jy+MGkxlTS3vrNzKjyaN5T9jgwN+9nfHcOvzizlqWD+uP/dwrvtj8xMm/nTy0dz+YvR6mh+cdziL1pcx+YQRXDB2KI+99yn/csaouMEJ6dRcV09w3tEEwtziH9a/R9yY/Be+ezqzlpXw4Jsref47p3PNY3PrRw9MGFPA5ONHMPn4EVRU1/LrWav4zoTDOPL2VwB4b+o5vLp4M//51084f+xQvn32aAAeuPIEvjdhd/2onFQCv86KOy/iiNtejls27/bzU95fZ3XpCSO49IS0TE3VLlMvOpKbJ46JW9YnhXtEdHbjR+2f7LBbXg498vcH8Ns3TWBQ7248MmsV76zcyqEFvetfmzRuOLc+v5hrTj2EQ4f0prE//OPJ7Kms5b3VW7n8pIMaBP8Rcet95+zDmmybDcF7Zxvoyrn/3QmjeWTW6lbXW3jHBXTLzaFkV2Xc0Mr3p57Dii2767tbCvp25/iDBvBv50f/EF/6/he577UVHDK4F1/7wiH12/XIz+XfL4gGwJ+mjKdsXzUjBvTkzCOiQ+zOPGL/9Bnd8nIYO7xfhxyv7iMQlZfFy83DOKtrpME1bXXnHa47azQnHjyQ8aMGMWncMK48+SD698yP+yY295ZzOeWuN+ufnz0m2h2Tyesk2iPYwZ+Gv+Przz28vp+uPb511iiWbNjF+FGDuO6s0Zx971tMPPpAbrhwDN3zcjAzbrzwSIC4/miIXv056aF3ue+rx9GvR/QiqIMH92LqRUdy98vL6tdtaQzwgF7d+OmlLXehfKFB6+iwA/oy77bzGNQ78Y1gOsItFx/JXTOWtXgyMcjevvHslO7EJqkb2Dv6/+emiWPqrxvJzTFOjU2D/vDfn5hwuwP69Ui4vKsI9F9ZslM2XHTMgby8eHNS627b0/pEUQf07U5txHn8n07hmBH96y8Aqa6NUBuJXvXZ+KrVd28+p9n9NW6JHTOif8J+4OvOGs2kccPSNrXy4DQH8pQzRzPlzNFpLaMzO2Rw0y4ESa++PfJZc9fFKTUSLz72QGYsSpwb137x0IyPvmqLYAd/km/mo18/Ke7kYkN1F96MHdaPy04cwT+cWkjfHvls2VVBJOIUFUYvPumWl9Nsa7guuPNzc8jPJeEVqR1l5MDUhoyJhFVOit1rv/5aNDfOOHxIk9dunzS2vdVKq0AHfzKmnDkKgD9/69T6S8EBfvGVcZx22GBGDuzV5E28eeKRGa2jiHROc245t11zTmVLoIO/Z+xkzeUnjYybb6bOZSeMYGosxE85dJDuVysibTK0i/b1Bzr4ryg6iJJdlVx31mi+deYofvjcIob06c4vrzye3BzTSBIRCaVAX8AlIhJmmqtHREQABb+ISOgo+EVEQkbBLyISMlkJfjObaGbLzWyVmU3NRh1ERMIq48FvZrnAI8BFwFjgajPr3Je5iYgESDZa/KcAq9x9jbtXAU8Dk7NQDxGRUMpG8I8A1jV4vj62LI6ZTTGzYjMrLi0tzVjlRESCLhtX7iaaEanJVWTuPg2YBmBmpWb2WYrlDQG2prhtV6VjDgcdczi055gPSbQwG8G/HjiowfORwMaWNnD3gpZeb4mZFSe6ci3IdMzhoGMOh3Qccza6ev4GHG5mh5pZN+Aq4C9ZqIeISChlvMXv7jVm9j3gVSAXeMzdl2S6HiIiYZWV2TndfQYwI0PFTctQOZ2JjjkcdMzh0OHH3CVm5xQRkY6jKRtEREJGwS8iEjKBDv6gzglkZmvNbJGZLTCz4tiyQWb2upmtjP0c2GD9H8Z+B8vN7MLs1Tx5ZvaYmZWY2eIGy9p8jGZ2Uux3tcrMfmVmqd1ZOwOaOeY7zGxD7L1eYGYXN3gtCMd8kJnNMrOlZrbEzK6PLQ/se93CMWfuvXb3QP4jOmJoNTAK6AZ8DIzNdr066NjWAkMaLfsFMDX2eCrw89jjsbFj7w4cGvud5Gb7GJI4xjOBE4HF7TlGYC5wKtELB18GLsr2sbXxmO8AbkiwblCOeRhwYuxxX2BF7NgC+163cMwZe6+D3OIP25xAk4HHY48fBy5tsPxpd69090+BVUR/N52au88Gtjda3KZjNLNhQD93/8Cj/0ueaLBNp9PMMTcnKMe8yd0/ij0uB5YSncIlsO91C8fcnA4/5iAHf1JzAnVRDrxmZvPMbEps2VB33wTRPyzggNjyIP0e2nqMI2KPGy/var5nZgtjXUF1XR6BO2YzKwROAOYQkve60TFDht7rIAd/UnMCdVGnu/uJRKe2/q6ZndnCukH+PdRp7hiDcOyPAqOB44FNwH2x5YE6ZjPrAzwL/MDdd7W0aoJlXfK4Exxzxt7rIAd/m+cE6ircfWPsZwnwPNGumy2xr37EfpbEVg/S76Gtx7g+9rjx8i7D3be4e627R4DfsL+bLjDHbGb5RAPwKXd/LrY40O91omPO5Hsd5OAP5JxAZtbbzPrWPQYuABYTPbZrYqtdA7wYe/wX4Coz625mhwKHEz0h1BW16RhjXQTlZjY+NtrhHxps0yXUhV/M3xF9ryEgxxyr4++Ape5+f4OXAvteN3fMGX2vs32GO81nzy8mesZ8NXBrtuvTQcc0iugZ/o+BJXXHBQwG3gRWxn4OarDNrbHfwXI66UiHBMc5nejX3WqiLZtrUzlGoCj2H2g18DCxq9U7479mjvlJYBGwMBYAwwJ2zF8k2j2xEFgQ+3dxkN/rFo45Y++1pmwQEQmZIHf1iIhIAgp+EZGQUfCLiISMgl9EJGQU/CIiIaPgF0mSmf3AzHplux4i7aXhnCJJMrO1QJG7b812XUTaIyv33BXp7GJXRf+Z6GXwucD/AcOBWWa21d0nmNkFwE+ITpe7GvhHd98d+4D4EzAhtru/d/dVmT4Gkeaoq0cksYnARnc/zt2PAR4gOg/KhFjoDwFuA87z6IR5xcC/N9h+l7ufQvRqygcyWnORVij4RRJbBJxnZj83szPcvazR6+OJ3iDjPTNbQHQ+mUMavD69wc9T011ZkbZQV49IAu6+wsxOIjqHyn+Z2WuNVjHgdXe/urldNPNYJOvU4hdJwMyGA3vd/Y/AvURviVhO9FZ5AB8Cp5vZYbH1e5nZEQ12cWWDnx9kptYiyVGLXySxY4F7zCxCdLbMbxPtsnnZzDbF+vm/CUw3s+6xbW4jOhssQHczm0O0cdXctwKRrNBwTpEOpmGf0tmpq0dEJGTU4hcRCRm1+EVEQkbBLyISMgp+EZGQUfCLiISMgl9EJGT+PzOKRwsxAHziAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "guess = 8.5\n", "\n", "pyro.clear_param_store()\n", "svi = pyro.infer.SVI(model=conditioned_scale, \n", " guide=scale_parametrized_guide,\n", " optim=pyro.optim.Adam({\"lr\": 0.003}),\n", " loss=pyro.infer.Trace_ELBO())\n", "\n", "\n", "losses, a, b = [], [], []\n", "num_steps = 2500\n", "for t in range(num_steps):\n", " losses.append(svi.step(guess))\n", " a.append(pyro.param(\"a\").item())\n", " b.append(pyro.param(\"b\").item())\n", " \n", "plt.plot(losses)\n", "plt.title(\"ELBO\")\n", "plt.xlabel(\"step\")\n", "plt.ylabel(\"loss\");\n", "print('a = ',pyro.param(\"a\").item())\n", "print('b = ', pyro.param(\"b\").item())" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXycZbnw8d+Vfd/3pU2XdN8opWVnZC0oIKDHoogH9SCKSjkv5wDqC0frUpXjy1FULIqoIOCRqpSlUBCoUui+pUvadEvS7Pu+3+8fz2SapkmbZZ6Zycz1/XzyIfPMnck1JTPX3Nt1izEGpZRSytcEeTsApZRSaiiaoJRSSvkkTVBKKaV8kiYopZRSPkkTlFJKKZ8U4u0ARislJcXk5eV5Owyl2L59e40xJtXbcYyVvpaUrxjutTThElReXh7btm3zdhhKISInvB3DeOhrSfmK4V5LOsSnlFLKJ2mCUkop5ZM0QSmllPJJmqCUUkr5JE1QSimlfJImKKV8nIg8LSJVIlIwzP0iIj8VkSIR2SMiiz0do1J20ASllO97Blh+lvuvB/KdX3cDv/RATErZThOUUj7OGLMRqDtLk5uB3xvLh0CCiGR6Jjql7GNrghKR+0SkQET2icjKIe7/jHNIYo+IbBKRhXbGo0ZHzwqbMLKBkgG3S53XziAid4vINhHZVl1dPeSDGWPo6O51f5RKjZJtCUpE5gH/BiwFFgIfE5H8Qc2OAVcYYxYAq4A1dsWjRqarp4+vPb+T2365ieWP/4Pmjm5vh6TOTYa4NuSnC2PMGmPMEmPMktTUoas0rXrlAEu/95Y741NqTOzsQc0GPjTGtBljeoD3gFsGNjDGbDLG1Dtvfgjk2BiPGoHnNp9g3e4ytp+op7CymZd3l3k7JHVupUDugNs5wJj/x8VEhNDc2UNfn/aglXfZmaAKgMtFJFlEooAbOP1FNNgXgNeHumMkwxJq/Lp7+/jBawcBeOrOJSRFh7HpSK2Xo1Ij8DJwp3M134VAozGmfKwPFhcRgjHQ3NnjvgiVGgPbisUaYw6IyA+BDUALsBsY8i9eRD6ClaAuHeax1uAc/luyZIl+rLPJxkPVdPX28ZvPLeGq2emsL0hjw/4Kunv7CA3W9TTeIiLPAw4gRURKgUeBUABjzJPAa1gfAIuANuCu8fy+uMhQAJrau4l3fq+UN9hazdwY8xvgNwAi8n2soYjTiMgC4NfA9cYY/bjuRW8dqCI2IoTLZ1hzE9fMSeOlHaXsLG5g6ZQkL0cXuIwxt5/jfgPc667fFxfhTFA6/6i8zO5VfGnO/04CbgWeH3T/JGAt8FljzCE7Y1HndqiymdmZca7e0sXTUwgOEt47VOXlyJQnxUVan1ub2nWIT3mX3eM2L4nIfmAdcK8xpl5E7hGRe5z3PwIkA78QkV0ioofTeIkxhkOVzcxIj3Fdi4sIZfGkBDYeqvFiZMrT+of1Gtq6vByJCnR2D/FdNsS1Jwd8/0Xgi3bGoM7U22f41cYjCMKXLp9KUJBQ2dRJc0cPM9JjT2t7eX4q/73hEDUtnaTEhHspYuVJabERAFS3dHo5EhXodOY7AL1bWMWP1hfyw/UHeXLjEQAOVzUDMD0t5rS2/fNR7xdpLypQJEeHERwkVDZ1eDsUFeA0QQWgv+06tUXmR+sL6ezp5VBlC8AZPah52fGEhQSxr6zJozEq7wkKElJjwqlq0h6U8i5NUAFoZ0k9H1uQyeOfWgTAzG+t55U9ZSRGhZIcHXZa2+AgITcxkuLaNm+EqrwkPS6cymZNUMq7NEEFmO7ePsoaOshLjub6+Rmu6zuLG7hkegoiZ1bNmZISzfHaVk+GqbwsLS6CKh3iU15m6yIJ5XvKGtrp7TNMSo4iPCSY9Ssv4/nNxXT3GR69cc6QP5OXHM0/DtfQ22cIDhqq7JvyN2mx4Ww7frYC6krZTxNUgCmus4bqJidFATArI45v3zzvrD8zJyuOzp4+1z4p5f/S4yKob+ums6eX8JBgb4ejApQO8QWYYzXWUN3k5OgR/8ySyVYViW0n6s/RUvmL9DhrS0G1zkMpL9IEFWCOVrcSHRbsegMaidykSFJjw9l+vI7Onl42HanRs6L8XFqctReqUlfyKS/SBBVgDpQ3MT09dsjFEMMREZZMTmTbiXp+suEQn35qMxv2V9oYpfK2tFjrA4wulFDepAkqgBhj2F/exLys0c8jLZuSRGl9Oy9tt+r9vntIjz3xZ+nOHlSVDvEpL9IEFUCO17bR3NHD3Kz4Uf/sVbPTAahpseqzbTxUrcN8fiwpKoyQIKG8UXtQyns0QfmowopmSurcuzn22Q9PIAKX5aeM+mdzEiNd31841epNbddFE34rKEiYlBzF8Rrd/6a8x+8SlMPh4JlnngGgu7sbh8PBs88+C0BbWxsOh4MXX3wRgMbGRhwOB2vXrgWgpqYGh8PBunXrAKioqMDhcLB+/XoASkpKcDgcvPXWWwAcPXoUh8PBe++9B0BhYSEOh4NNmzYBUFBQgMPhYOvWrQDs2rULh8PBrl27ANi6dSsOh4OCggIANm3ahMPh4ODBg1z3+EYu+9E7XLL8Fo4ePQrAW2+9hcPhoKSkBID169fjcDioqKgAYN26dTgcDmpqrLp5a9euxeFw0NjYSEldG7/951GiK3aTbI3e8Oyzz+JwOOjuts79eeaZZ3A4HK5/y6eeeoqrr74asOahUkKt3tM3b7D2Sz342w2cd89P+MQvN9HbZ3jssce47bbbXD+/evVqVqxY4bq9atUq7rjjDtftRx55hLvuOnW23sMPP8zdd9/tuv3AAw9w772njjlauXIlK1eudN2+9957eeCBB1y37777bh5++GHX7bvuuotHHnnEdfuOO+5g1apVrtsrVqxg9erVrtu33XYbjz32mOv2TTfdRCDLT4tx1WhUyhv8LkH5g+rWU+fwtCVMc8tj/vb944hAwslNY36Mj6dUk33sNebnxLPiglyOdMZQnzCTbSfq2VGsvSl/Mz0thuO1bXT19Hk7FBWgZKLNIyxZssRs2+bfx0a9U1jFXb+1el23Ls7mJ/+yaNyPufzxjaTGhvOHLywb92MBvLS9lP/zv7tdt7/smMaDy2e55bEnChHZboxZ4u04xupcr6W/7jzJyhd3seH+y8kfVERYKXca7rWkPSgfdKjCGlY5b1IC+06eu4p4VVMH3b3Wp9zq5k66evro6T31qbeutYuDFc0szXPfse3LplqPtXhSAkunJI34OI72rl7WF1TwTqGe0uvr+o9eOVzV4uVIVKDSBOWDCsqayE6I5JJpKRRVt9DR3Tts24rGDpb94G3++03rUMELvvcWM771OtO/+bprgnurs6baRdOS3RZjTmIU61dexotfuohFuQkcLG92Jcmz+drzO7jn2e3c9dut1LeeOrHVGMPP3yliT2kDLZ161LgvmJYagwgcrtQEpbxDa/H5oL2lDczPjmdedhy9fYbCimYW5iYM2fb9ohqMgSffO3JGL8bx2Lt875Z57CtrIiI0iPk5o19efjazMqz9VHOz4ujq7aOoqmXIWn17Sxt5/0gNty7O5q0Dp3pOO0vqccxIo6Kpg6v++z3au3v58RuFAPzg1vmsuCB3VBuKlXtFhgWTkxipCyWU12gPysc0tnVzvLaN+TnxrqS06UjtkG3LG9tPmwfae7IRx8xUjv3gBh6+3poP+uZfCvjj5mKWTUm2rejnXOfG3+EONVy9/gCrXz/ITT97H4BXvnYpUWHBbNhfybfX7ePi1X+nfVAv8eG1e1m3p9yWeNXI5afFUqRDfMpLNEH5mE1HrF7QotwEMuMjmZ8dz+sF5UNuiv3yszsAuHFhluva/3zqPESEL10xjV999nzXdXcO7w02JSWGmPAQPjw6dCItrLDe4CqaOrhyVhrzsuNZPjeD57eU8LsPTrjaffDwldwwP4M1zrjfPajzVN6WnxbD0ZrW0+Y0lfIUTVA+5plNx8lOiOTCqVZCWbE0lz2ljXzjLwWs3VHqSlRlDe3sKmngU0tyefxTi1i/8jLevP9y4qNCXY913dwMshOsDbbzxlA9YqSCg4QbF2by6p5yKho7Tpsza2zrpqalk0+en8PHFmTyw9sWAHD3FVNdba6encZLX76IzPhIfvGZ87l2bgbXzU336tL1dwqrWPnCzoB/Y85Pj6Wrp4+iau1FKc/TBOVDGtq62Hysjk8uyXEdDPjppZOICA3i+S3F/PufdruG+/552Opp/esleQQHCbMy4pgxxFLgb310Nkvzklg0aeg5LHf5zLLJtHf3cuEP3uYrz+2gpqWToqpmDlRYw343zM/kiU8vJtVZhHRWRhy/+/xSVlyQy+MrzuP8yaevMJyVEceJurazLhAZrL61i5Uv7Bx39QNjDHf9dit/3VXGewFec/Dy/BRE4PuvHfR2KCoAaYLyIf09hmVTTg3HiQgb7r+Cjy3IBOCeP2ynsKKZX208Ql5yFDPPsT/l+vmZ/Omei4gJt3c9zLzseOY4F0j8/WAVS777Flf/ZCMr1nwInFqyPNAVM1JZfduCIWPLT4/BGOt4kOG8treczz+zlQ+cSfvp94/x111lPP3+sVHFfu9zO/gP51xeSV0bJXXtAExLjWbRMItTAkVaXAQ3Lsjin4erA743qTzP1gQlIveJSIGI7BORlUPcP0tEPhCRThF5YKjH8HdlDe3sK7NKEX3+mW2IcMabYm5SFE98ejGrb51Pc2cP1z2+kSPVrXzFMZ0gHzqC/Y//tozff37pkPf1DzWO1Kk9OMOvIPv2un38/WAVtz/1IcYYV23Av49i7qqutYtX95bzv9tLyXvoVS770Ttc/uN3APjp7eeRHDPyc7P81dIpSfQZqB2wLUApT7DtY7WIzAP+DVgKdAHrReRVY8zhAc3qgK8DH7crDl/2flENn/n15tOu/deNc4kMG3q13Yqlk8hOjOSzv9nCxdOS+cT5OZ4Ic8QSosK4fEYqB1ctp7q5k5zESM5btYFZGbGjTqRTUqIJEiiqamFXSQMHypu4femk09oEDViCXlrf7jotuLS+nca27tPm44az92TjsPf1L6MPdKmus6E6XcdwKOUJdo77zAY+NMa0AYjIe8AtwI/6GxhjqoAqEfmojXH4rOc2nzjt9jVz0vncxXln/ZnL8lPZ/I2rSIsN99k9QhGhweQmRQHwzwevJDR49HGGhwSTlxzN4coWfvb3IgAuyEtkepo1pNnR3UtFUweOmam8W1jNzpIGKpo6WJSbwK6SBg5WNLHMudCkpbOHk/XtzMw4czh0x4l6ggTef+hKokJD6OjpZdn33yYvOco1Dxjo+g8vrG7pAOxbbKPUYHYmqALgeyKSDLQDNwBjKqInIncDdwNMmjTpHK0njoPlzVw3N53v3TKf/7fhEJ+/dMqIfm4ifYodz9zXrMxY/nHo1ObjX7x7hMc+sZCgIKG0vg1jrJWK7x2qZsP+SoyBG+ZnsKukgW/9tYALpyZzuKqZD49alTS2f+vq04bstp+o43/ePkxIkJAZbw1BxhPKX75yMZOTo8cct7/p70FV6+GFysNsS1DGmAMi8kNgA9AC7AbGVMPGGLMGWANWgUu3BelF7V29HKtt5aZFWaTEhPO9W+Z7OySfc9HUZF7bW+G6vXbHSf5+sIrZGXGUN1oLGeZkxjElOZo3Cqx2/asBD1e1nFFDbvuJeq6dm+G6/d1XDwBw9+VTT2t33qRE9z+ZCWzgEJ9SnmTrIgljzG+MMYuNMZdjzTcdPtfPBIrDVc0YwzlX4QWyjy04tQH5l59ZzJSUaBrauvngaC3Ha63DHOdmxTHbWWoJYFJSFI9/6lT194SoUP71Ymsp/p7SU/NNLZ097C5p4EtXTOU/A6wK+2iFhwQTHxmqx78rj7N17bGIpBljqkRkEnArcJGdv28i2e18sxxqXkRZEqPD+OMXl1Hf1s318zO5fn4mhyub2Xq8nqrmDj5xfg4hwUHMyYzj1T3lhIUEkRITxsfPy+ay/BQqmjpcx9tvPlbHngELIrafqKfPwCXTRn+6cCDKjI+gokmPf1eeZXex2Jecc1DdwL3GmHoRuQfAGPOkiGRgzUvFAX3OpehzjDHnPmNiAqtr7eKpjUeZkR7DlBSd6zibi6efnkDy02PPOJuof/9VaJC4Fo4kx4SfNt+0IDueN/dXYIxBRNhUVENosLAkT4fzRiIjPoKKRk1QyrNsTVDGmMuGuPbkgO8rAN9aK22zvj7Dmo1HKa5r49kvLPPZlXgTSX+V9rMtMpmfE8+L20o42dBOTmIU7x+pYfGkRKLCtKD/SGTGR1AwgrPJlHInfXV62ENr9/CnbaVcNDWZS/N1eMkdUmLC2f3ItWfd97TAmcR2FDfw9oEqCk428X+umeGpECe8jLhIalo66ezpta0qvlKDaYLyoIKTjfxpWymhwcID1+mbozuda1PuzIxYQoOFrz+/03Xt6jnpdoflN3KTrGX4pfXtTEs9s2yVUnbQWnwe9Ma+CoIENn/j6jOKoyp7hYcEn7YgZfWt84c8XFENrX9f2HgL8So1GpqgPKSvz7B2x0mW5CWRFB3m7XACUlqstcH5oetnsWKp/2z49oT+xTzHNEEpD9IhPg95+v1jnGxo554rpp67sbLFFy+dQnNHN7dfoMlptJKiw0iMCuWoJijlQZqgPOSlHScJCwnijgsnezuUgHXx9JQzlq2rkZuaGsOxsxx/opS76RCfBxScbORAeRMPLp+ly8rVhJUeF05Vs+6FUp6jCcpmxhj+8897SI4O47bF2d4OR6kxS4kJp6ZFz4RSnqMJymY/fqOQ/eVNPHDdTBKidHGEGj0RWS4ihSJSJCIPDXF/ooj8RUT2iMgW51lsbpcSE05jezedPb12PLxSZ9AEZbMPjtYSGRrMrdp7UmMgIsHAz4HrgTnA7SIyZ1CzbwC7jDELgDuB/7Ejlv6q5rXai1IeognKRr19hoPlzXzqglzdfa/GailQZIw5aozpAl4Abh7UZg7wNoAx5iCQJyJu34Wc4qxtWNOiVc2VZ2iCstGxmlbau3uZm6UbQtWYZQMlA26XOq8NtBvrtABEZCkwmWFqXIrI3SKyTUS2VVdXjyqQlBhriFoTlPIUTVA22l9uFdecowlKjd1Qyz4HH9q5GkgUkV3A14CdDHM4qDFmjTFmiTFmSWpq6qgCcfWgmnWIT3mG7oOy0b6yRkKDhfw0PfNJjVkpkDvgdg5QNrCB83iauwDE2sdwzPnlVq6j37UHpTxEe1A22l/WRH5aLGEh+s+sxmwrkC8iU0QkDFgBvDywgYgkOO8D+CKw0Y4z1SJCg4kND6FaT9ZVHqI9KJsYY9hf1sSVs9K8HYqawIwxPSLyVeANIBh42hizb+DBn8Bs4Pci0gvsB75gVzwpseE6B6U8RhOUTcobO6ht7dL5JzVuxpjXgNcGXRt48OcHQL4nYkmJCdMEpTxGx55ssvV4HQAX5OmxGsp/aDUJ5UmaoGzy5r5KYsND9Mwh5VdSdYhPeZAmKBt09fSx4UAlNy3KIjhIi8Mq/5ESE05DWzddPX3eDkUFAE1QNthf3kRXTx+X6NEOys/074WqbdVelLKfJigb/OztwwCcPznRy5Eo5V6uahK6WVd5gCYoN6tt6eTtg1UszIknPS7C2+Eo5VYpsVqPT3mOJig3e72gAoD/+7HBBaeVmvhSnUN8ullXeYKtCUpE7hORAhHZJyIrh7hfROSnznNu9ojIYjvjsVtlUwf/9fI+YiNCWJCT4O1wlHK7zPgIwkOCOFzV7O1QVACwLUE5D037N6zjAhYCHxORwZsJr8faYJgP3A380q54POG5D0/Q02f42e3naXkj5ZdCgoOYnx3Pe4eqMWZwzVql3MvOd9HZwIfGmDZjTA/wHnDLoDY3A783lg+BBBHJtDEm2+wqaeDn7x7hyllpOGZqeSPlv/7lglwOVbaw5Vidt0NRfs7OBFUAXC4iySISBdzA6VWZYWRn3YzrDBtP+dV7R+jtM/zn8pneDkUpW10z2zoLce/JRi9HovydbQnKGHMA+CGwAViPdaja4DNqRnLWzbjOsPEEYwybj9Vx2+IcZmVo5Qjl3xKiQokND6Gkrs3boSg/Z+tEiTHmN8aYxcaYy4E64PCgJuc862YiOFHbRl1rF4ty470dilK2ExEmJUdxqLLF26EoP2f3Kr40538nYR1J/fygJi8DdzpX810INBpjyu2MyQ5/22Xl1I/o0RoqQFw+I5Utx+tobOv2dijKj9m91OwlEdkPrAPuNcbUi8g9/WfZYB0hcBQoAp4CvmJzPLb4666TXDo9hZzEKG+HopRHXDUrjd4+w8LvvElVc4e3w1F+ytbzoIwxlw1xbeA5Nga4184Y7Ha8ppVjNa18Ztkkb4eilMcsyk0gLDiIrt4+3i+q4ZbzcrwdkvJDullnnH73wXFCgoTr50/I1fFKjUlIcBAHVy0nKiyY3SW6mk/ZQxPUOH1wpJZL81PIToj0dihKeVRQkJARF0G11uVTNtEENU5lDe1MStK5JxWYUmLDqW7SBKXsoQlqHFo7e2jq6CEjXquWq8CUFhuuPShlG01Q41DeaK1eyorX4T0VmFJjw7WyubKNJqhxqHAmKO1BqUCVGhtOS2cPbV2Di8QoNX6aoMahvLEdsI4gUCoQ6flQyk6aoMahf4hPT85VgSop2joCvl4rSigbaIIah+O1raTFhhMRGuztUJTyioQoZ4Jq7fJyJMofaYIahwPlzczK1OrlKnCd6kFpglLupwlqjLp7+yiqamZ2Zqy3Q1HKaxKjQgEd4lP20AQ1RsV1bXT3GmakaYJSgSsuIpQg0SE+ZQ9NUGNUXGsd1paXolUkVOAKChKSY3QvlLKHJqgxKnaeJpqrZY5UgMtKiKTMueVCKXfSBDVGx2paiQ4Ldu0DUSpQZSdEcLJBE5RyP01QY3SkuoVpaTGIiLdDUcqrsuIjKWtoxzreTSn30QQ1RkVVLUxLjfF2GEp5XVZCJB3dfbqST7mdJqgxaOnsobyxg+lpmqCUynKehVamw3zKzTRBjcELW4oBdA+UUuA6rFPnoZS7aYIapaqmDr776gEALpme4uVolPK+rASrFqX2oJS7aYIapV9tPEpIkPDOAw7CQ7QGn1JJ0WGEhwRpglJupwlqFAormvnNP4+xeHIiU1KivR2OUj5BRMhOiKSsocPboSg/owlqhN4vquG6xzcC8GXHNC9Ho5RvyUqI1Dko5Xa2JigRuV9E9olIgYg8LyIRg+5PFJG/iMgeEdkiIvPsjGc8HnuzEIBnv7CMj8xM83I0SvmWrIQIHeJTbmdbghKRbODrwBJjzDwgGFgxqNk3gF3GmAXAncD/2BXPeJTUtbGzuIH/XD6TS/N1YYRSg2UnRFHV3ElnT6+3Q1F+xO4hvhAgUkRCgCigbND9c4C3AYwxB4E8EUm3OaZRe7ewCoDr52V6ORKlfFP/Sr7KRi0aq9zHtgRljDkJPAYUA+VAozHmzUHNdgO3AojIUmAykDP4sUTkbhHZJiLbqqur7Qp5WB8eqyMzPoK8ZC0MqzxPRJaLSKGIFInIQ0PcHy8i60Rkt3NI/S5Px6h7oZQd7BziSwRuBqYAWUC0iNwxqNlqIFFEdgFfA3YCPYMfyxizxhizxBizJDU11a6Qh2SMYcuxOpZOSdK6e8rjRCQY+DlwPdaIw+0iMmdQs3uB/caYhYAD+G8RCfNknFpNQtkhxMbHvho4ZoypBhCRtcDFwLP9DYwxTcBdzvsFOOb88hlFVS1UN3eybEqyt0NRgWkpUGSMOQogIi9gffDbP6CNAWKdr6EYoI4hPujZKSPeGuLTHpRyJzvnoIqBC0UkyvnCuQo4MLCBiCQM+KT3RWCjM2n5jA0HKgG4cpau3FNekQ2UDLhd6rw20BPAbKw53r3AfcaYPs+EZ4kIDSY7IZLCymZP/lrl5+ycg9oM/BnYgfWiCQLWiMg9InKPs9lsYJ+IHMQawrjPrnjGasuxOmakx7g+ISrlYUONKw8+1+I6YBfWUPoi4AkRiRvywWycz108OZFtx+v02A3lNrau4jPGPGqMmWWMmWeM+awxptMY86Qx5knn/R8YY/KdbW41xtTbGc9o1bd28W5hNfOy470digpcpUDugNs5nLka9i5grbEUYQ2Tzxrqweycz71oajKVTZ3sK/OpQRA1gWklibP46vM7AJiTOeSHUaU8YSuQLyJTnMPhK4CXB7UpxhpCx7lNYyZw1KNRAtfPyyAkSPjbrpOe/tXKT9m5SGJC6+nt4/2iWsJDgvjcxXneDkcFKGNMj4h8FXgDa7P708aYff3D5M7RiFXAMyKyF2tI8EFjTI2nY02MDsMxM42n/nGM8sYOvvvxeSREeXQxofIzmqCG0T/Z+6NPLCA0WDuaynuMMa8Brw269uSA78uAaz0d11A+cX4Obx2o5JU95czPjudLV2jdSjV2+s47jD2ljQAszEnwciRKTRzXzU3nmbsuYFJSFH/ZqUN9anw0QQ1jf1kTseEhTErS6hFKjZSI4JiZxvXzMzha3Upfn67oU2OnCWoY+8ubmJUZS1CQVo9QarSy4iPp6u2jrq3L26GoCUwT1BB6+wz7y5qYm6XLy5Uai0znvsGKRj3EUI2dJqghHKtpob27l/m6/0mpMcmM19p8avxGvIrPWfw1H3CVVDDGbLQjKG/bUdwAoBt0lds5D+38CnApVkWIfwK/NMb4VVejv/JKRZNfPS3lYSNKUCLyRawyRDlYJVUuBD4ArrQvNO9obOvm8Q2HyE2KJD8txtvhKP/ze6AZ+Jnz9u3AH4BPei0iGyRHhxEaLJTrEJ8ah5H2oO4DLgA+NMZ8RERmAd+2Lyzv6Ozp5eLVb9Pa1ctvPrdEF0goO8x0HovR7x0R2e21aGwSFCRkxEdQWq9DfGrsRjoH1dE/BCEi4c7Tb2faF5Z3vLqnnNauXq6clcZVs33uYF/lH3aKyIX9N0RkGfC+F+OxzeyMOApONtLepcfAq7EZaYIqFZEE4K/ABhH5G2cWrJzwNh2pBWDNZ8/3ciTK34jIXhHZAywDNonIcRE5hjVUfrl3o7PH4smJHKtpZfYj66nSuSg1BiMa4jPG3OL89ouiL9oAACAASURBVL9E5B0gHlhvW1Resqe0gY/MTCVESxsp9/uYtwPwtBsXZrH69YMA7Cxp4Lq5GV6OSE00o34nNsa8Z4x52RjjVzvwWjp7OFzVwsJcLW2k3M8Yc+JsX96Ozw7ZCZH84z8/AsCukgYvR6MmIu0qOBWcbMQYrb2nlDvlJkWxMCee7cd96qg3NUFognLaU2p9wluQo3uflHKny2eksuV4HQcr9CBDNTqaoJx2nGggNymS5Jhwb4eilF/5lyXWgcDbT2gvSo2OJijAGMO2E/UsmZzk7VCU8js5iZHEhodwsLzZ26GoCUYTFFBc10ZNSyfnT070dihK+R0RYUpqNCfq2rwdippgNEEBW50TuEvyNEEpZYfcpChKNEGpUdIEBewuaSAmPIQZabHeDkUpvzQpKYrS+jZ69QBDNQoBn6CMMWw6UsP87HitvaeUTSYnRdHda/T4DTUqtiYoEblfRPaJSIGIPO88amDg/fEisk5Edjvb3WVnPEM5VNnCkepWblyY5elfrVTAmJISDcDRmlYvR6ImEtsSlIhkA18Hlhhj5gHBwIpBze4F9jurOzuA/xaRMLtiGkr/3gydf1LKPlNTraNrCk428rtNx+no1gKy6txGfGDhOB4/UkS6gSjOLDBrgFgRESAGqAN6bI7pNCdqrYnbSUlRnvy1SgWUlJgwwkKC+PEbhQCEBgfx6WWTvByV8nW29aCMMSeBx4BioBxoNMa8OajZE8BsrMS1F7jPGNM3+LFE5G4R2SYi26qrq90a5/HaVjLjI4gIDXbr4yqlThERLs9Pcd1+p7DKi9GoicLOIb5E4GZgCpAFRIvIHYOaXYd1Qm8WsAh4QkTiBj+WMWaNMWaJMWZJamqqW+M8UdvG5GTtPSllt19/7gKeuesCLpqazLbjdd4OR00Adi6SuBo4ZoypNsZ0A2uBiwe1uQtYayxFwDFglo0xneFEbSt5ydGe/JVKBSzHzDQcM1Opb+umsb3b2+EoH2dngioGLhSRKOcc01XAgSHaXAUgIulYp/QetTGm07R09lDT0sVkTVBKeUyec0XfMV3Rp87BzjmozcCfgR1Y80tBwBoRuUdE7nE2WwVcLCJ7gbeBB40xNXbFNNiJWusFkqdDfEp5zOwMaxR/X1mjlyNRvs7WVXzGmEeBRwddfnLA/WXAtXbGcDb9K/i0B6WU5+QmRRIfGcre0kZY5u1olC8L6EoSx509KF0koZTniAgLcuLZXao9KHV2AZ2gTtS0kRITTnS43dvBlFIDLcpNoLCiidZOj257VBNMQCeo47WtOv+klBcszEmgz8CBcj1lVw0voBOUtQdK55+U8rSZGdbJAYcqW7wcifJlAZugOrp7qWjq0B6UUl6QnRBJdFgwhyr1lF01vIBNUEVV1ie3Kanag1LK04KChPz0WAorNEGp4QVsgurfgzE3K97LkSgVmGZlxGoPSp1VwCaogxXNRIYGM1mrmCvlFfnpsdS2dlHb0untUJSPCtgEdbiyhfz0GD1FVykvyU6wzi+taOrwciTKVwVsgjpU2Ux+Wqy3w1AqYKXEhANQ3aw9KDW0gExQRVXNVDV3MitDE5RS3tKfoGpaurwcifJVAZmgHn15HwCXDjhATSnlWSmx/QlKe1BqaAGZoI5UtTIzPZbZmWecjaiU8pDosGAiQoOo0SE+NYyAS1CtnT1UNnfw0QWZ3g5FqYAmIqTGhlOtPSg1jIBLUIcqmzEGnX9SE4aILBeRQhEpEpGHhrj/P0Rkl/OrQER6RSTJG7GOVnpsBJW6ik8NI+AS1M7iBgAd3lMTgogEAz8HrgfmALeLyJyBbYwxPzbGLDLGLAIeBt4zxtR5PtrRS4+PoLJJe1BqaAGXoP68vZSZ6bHk6gZdNTEsBYqMMUeNMV3AC8DNZ2l/O/C8RyJzg8y4CMob2zHGeDsU5YMCKkE1d3Szv7xJ55/URJINlAy4Xeq8dgYRiQKWAy8N92AicreIbBORbdXV1W4NdCwy4iPo6O6jqV3PhVJnCqgEVVrfDsC01BgvR6LUiA1V6mS47saNwPtnG94zxqwxxiwxxixJTU11S4DjkRFvVZMob2r3ciTKFwVUgiqpawMgJzHSy5EoNWKlQO6A2zlA2TBtVzCBhvcAMvsTVKMulFBnCqgE1d+D0vknNYFsBfJFZIqIhGEloZcHNxKReOAK4G8ejm9c0mKtBFWtCyXUEEK8HYAnldS3ERUWTGJUqLdDUWpEjDE9IvJV4A0gGHjaGLNPRO5x3v+ks+ktwJvGmFYvhTomqc5qElXN2oNSZwqoBFVa305uYhQiWsFcTRzGmNeA1wZde3LQ7WeAZzwXlXtEhAYTFxFClVaTUEMIqCG+kro2nX9SysekxUVQpUN8agi2JigRuV9E9jl3tz8vIhGD7vfYDnhjDCfr23X+SSkfkx4XrkN8aki2JSgRyQa+DiwxxszDGj9fMbCNJ3fA17d109zZoz0opXxMemwEJxt0s646k91DfCFApIiEAFEMvzwWbN4BX1TVAsD0NN0DpZQvuWBKEpVNnRwob/Z2KMrH2JagjDEngceAYqAcaDTGvDlU23PtgHfH7vfCiiZAE5RSvuaCvEQADjpfo0r1s3OILxGrZtgUIAuIFpE7hml+1h3w7tj9vvlYHRlxEWQn6BCfUr5kcnI0ocHCocoWb4eifIydQ3xXA8eMMdXGmG5gLXDxMG1t3wG//UQ9F0xJ0iXmSvmY0OAgshMiKa1v83YoysfYmaCKgQtFJEqsrHAVcGBwI0/sgG9s66a8sYN5WXrEhlK+KF2Xmqsh2DkHtRn4M7AD2Ov8XWtE5J7+XfBOtu+AP9mgJY6U8mXpcRFU6lJzNYitlSSMMY8Cjw667PEd8GXOBNVfmFIp5VvS48KpaOzAGKPD8MolICpJlDdaCUoXSCjlm9LjIujsOf1cqO7ePteHSxWYAiJBnWzoIDRYSIkJ93YoSqkhpMdZoxv9w3x/P1hJ/jdf5+LVf6da6/QFrIBIUAfKm5iaEkNQkA4dKOWL+hNU/7lQz35Y7LrveO2EKtCu3MjvE1Rfn2FXSQOLJyd4OxSl1DD6N9DvLW0A4FhNK9NSowF0+XkA8/sEdaS6hcb2bs7LTfR2KEqpYSRFhzErI5bH3jzEn7aVcKymlZsWZgNwsl7noQKV3yeofxbVAHBJfoqXI1FKnc1Ni7IA+M8/7wHgkunJpMSEu07CVoHH7xNUeWMHEaFBZOkSc6V82p0X5REfeeq063nZ8WQnRrr2MarA4/cJqqalk+TocN1boZSPiwkP4Z0HHCRHh/FlxzQiQoPJSYzUHlQA8/sj32tbukiOCfN2GEqpEUiKDmPbt652faDMSYhkw75K+vqMrsINQH7fg6pt7SQ5WhOUUhPFwNGOnMRIunr7qGnxvb1Qm4/W0tbVc+6Gasz8PkHVtXSRrBt0lZqQchKt+pklPjbMt+VYHZ9a8yFPvnfU26H4Nb9OUMYYalp1iE+piSo70SpP5mt7ofaebATghG4itpVfJ6iWzh66evp0iE+pCaq/fqavreQ7UG6d/tvd2+flSPybXyeoutYuAJKjdYhPqYkoOjyElJhwfrS+kE1Harway7rdZfxt10kA9pVZCaqsQY8IsZNfJ6iaFmeC0iE+pSasy5yb7D/91GavxdDR3cvXnt/JfS/s4ou/2+bqQfWflKDs4dcJqta58kermCs1cf37NTNc3zd1dHslhpUv7HJ9/9aBSgAW5iZQ1dypw3w28usE1T/El6RzUEpNWLlJUbzytUsBeH5z8Tlaj11LZw83//x9vvLcdjYeqnZdb2zrZv2+CkTgS5dPBSAuIoRPnJ+DMVDROLGH+dYXVPDOwSpvhzEkv96oW6sJSim/MC87nrlZcfyzqIYvXTHNlt/x/OZidpc0sLsEXttbwe5HriU+KpRtJ+oAeO6Ly7h4WgpfvXI6fQYKnCv5SurbyE2KsiUmux2vaeWeZ7eTEhPOtm9d7e1wzuDXPajali6iw4KJCA32dihKqXFakBPP3pONGGNsefzntxYTG37qM3uJc2n7luN1hAYLiydZJyLERoQSHxnKJGdSKq2buPNQW45bybempZPePnv+XcfDrxNUXWsnSbpAQim/kJ8WS0Nbt2vo3p0qGjs4Wt3Kymtm8NKXLwag2jmHvfVYHfOz48/4oJsWZ81tVzRN3CG+/c7ViOCbz8OvE1RtaxdJusRcKb8wJcU6wPBYjfs3x+52HpR43qQE0mKt94x3D1bxq/eOsPdkIxdMSTrjZ8JDgkmICp3QR9IfrDiVoMps2mt2oLyJHcX1Y+r5+nWCqmvt0k26SvmJ/gR11IYEtf1EPcFBwpzMOFKdCep3H5zgB68fpLvXsDTvzAQFkBoTfkaCOl7TSp8PDpcNZoyhsKLZ9dzsOhjyqY1HuecP28d0ooRfJ6jali5dIKGUn8hJjCQ8JIjDlc1ufdyNh6pZs/EoydFhRIQOPWd94dTkIX82LS78tKGxqqYOHI+9y/dfO+DWGO3wwZFa6tu6udF5UKRdQ3x7TjayICd+TD9ra4ISkftFZJ+IFIjI8yJyxqmBIuIQkV3Odu+563cbY7QHpZQfCQkOIj89hoMVwyeovaWNpy0RP5ealk7ufHoLACuvPrXfavncDNf387PjiQ4fesHz5ORojjvr8RVVtXDAGdsLW0tGHIO3vH+khuAg4ROLcwgJEpra3b/HrKWzhyPVLczPThjTz9u2zFxEsoGvA3OMMe0i8idgBfDMgDYJwC+A5caYYhFJc9fvb+nsoau3T3tQSvmRGemxbCqqPe1aX5/hRF0b0WHB3PjEPwH43i3zuOW8bKLChn+Le7+ohs/82qpOcdPCLD69bJLrvh9/cgHf+thsIkODz/oY01NjaGjrJu+hVwGrlwfQ1eP7m3cPVbYwJSWayLBg4iJDbdkEvbO4HmPwzR4UVgKMFJEQIAooG3T/p4G1xphiAGOM23aL9dfISo/To96V8heTk6KpaOqgo7vXde3b6/bxkcfe5aG1e13XvvmXAu5/cddQD+Ey8P7/WbHotPtiI0LJSYwiOSacyLDht6lcPTv9tNv9p//29PWNOkl9/pmt5D306mkLF+x0uLKZGekxgLXxuKnd/Wdb/XVnGbHhIcMOkZ6LbQnKGHMSeAwoBsqBRmPMm4OazQASReRdEdkuIncO9VgicreIbBORbdXVI+u+7yyuB8aeuZVSvmdysrX36PiAYy7e2GeVHvr7wSoW5MTz09vPc10/W6+gP/E8/qlFY5rAB5iUHOXqefUvrggNFvoMozpkcU9pA393VnP4267Bn+PHrqSujZ4hSjG1d/Vyoq6N/LRYAFt6UH19hncLq7h6TvpZk/zZ2JagRCQRuBmYAmQB0SJyx6BmIcD5wEeB64D/KyIzBrXBGLPGGLPEGLMkNTV1RL9/R3E9iVGhrpU/SqmJb1629YFz+eP/4IUtxXzjL3tPm9y/bXEONy3MYs1nzwfgR+sPDvk4nT29lNa387Urp/Px87LHFdP3b5nP8dUfZes3r+bQd6/nsU8uBKB9QC/vbHp6+1j1yn6ix/gmPpyGti4u+9E7fOeV/a5rpfVtLPnuWzy3+QTGWEOmAHERoW6fg9pZ0kBtaxeXz0gZ82PYOcR3NXDMGFNtjOkG1gIXD2pTCqw3xrQaY2qAjcBCd/zy3SWNLMpNGPMnI6WU75mWGk26c4PsQ2v38kdnbb6Xv3oJT96xmM84ezPXzEknISqUv+0sG3L/zdHqVnr7DNPTYtwaX1hIEJHOVYDtXSNLUNtP1LP1eD33XzODSUlRY1rubYw5oxJEoXPBxusFFa5r/zhcQ01LJ9991Vpl6Briiwyhwc0J6i87SwkLDuKqQcOgo2FngioGLhSRKLGyxFXA4LWXfwMuE5EQEYkClg3RZtT6+gzHa1vd/senlPIuEeHP91zsKjME8F83zmFBTgLL52USEhzkavflK6bR3NnDL987csbjHK5qAU71INypfzirY4Q9qP6jO25cmEV2QiQv7y7j+S2jK4r7pT9sZ/njG0+7Vuhcjh8aZH1Ir2vtYsuxOtf9ocFCnnOEKScxitL69iGHA8eiu7ePP20r5WMLM4mLCB3z49g5B7UZ+DOwA9jr/F1rROQeEbnH2eYAsB7YA2wBfm2MKRjv765q7qSzp49JyTq8p5S/yU2K4vr51jLw2xbn8K+XTBmyXf/c0NodJ8+473BlM8FBwtRU979H9Peg2kbYgyqsbCYhKpS02HAy4q1FXQ+v3Utx7ciOue/u7ePN/ZUcrmo5bYPw7hKrmG1lcyfNHd0sXrWBv+w89W9x4dRkQp0JfVZGLF09fW7bBH2ito2unj4umTb24T2weRWfMeZRY8wsY8w8Y8xnjTGdxpgnjTFPDmjzY2PMHGebx93xe4vrrP+xkyZohWGl1NndsWwyNy3M4qtXTh+2TWxEKHdfPpXiurYzhr8OVTYzOTmK8BD3F5LuX5be2jmyVXH7y5qYlRGLiJy26vh/t599L1VnTy/PfniCzUdP9YqqByzM6C/f1NtnXMN9ABdOtSpHXDtgr1d/T/KQmzZB95ejGu8HAL+sJHHCucJnsiYopfxSblIUP739vHMugpqRbvUMjtW0nHb9cGULM9LcP7wHEBthJajmjnMnqLrWLvaXN7Ewx9rI+uUrprlWIf7s70VnrZrx4tYSvvXXAr7y3HbXtVJnBfamjm6OVLe4lsG/dcBaIXjnRZN56s4lrLp5LrctPrU4ZHpaDEFi7Y1yh/734LxxjmL5ZYIqqWsjSCArIdLboSg1biKyXEQKRaRIRB4apo0tFVkmujmZcQAcKD/1Rt/Z08vx2lby0+2Zo+6fc2ke1IP64EgteQ+9yjf/Yu3XKmtoZ/GqDXT3Gi51HmsfHxXKTQuzuO+qfAD2lQ2/J+rDo9aG5aYBibDEefRHQWkjxsDyeVYv6fWCcgC+euV0YiNC+exFeadtQI4IDWZSUhRFVe7pQRXXtREbEUJC1Njnn8BPE9SJujayEiIJC/HLp6cCiIgEAz8HrgfmALeLyJxBbforstxkjJkLfNLjgfqovBRrFKV/2B+sFXx9BvJtWCABEOPqQZ2+Ku4X7xYB8Jxz5WH/vqfL8lO4dPrpczUrluYCZ5/HGjhst2RyIiKw9XgdJ2pb2eUc3rtqVhphwUGcqG0jOyGStNjhCxfkp8fy2t4KtxxncrCimakp0eNeRe2X7+DFdW06/6T8xVKgyBhz1BjTBbyAtb9wILdXZHE4HDzzzDMAdHd343A4ePbZZwFoa2vD4XDw4osvAtDY2IjD4WDt2rUA1NTU4HA4WLduHQAVFRU4HA7Wr18PQElJCQ6Hg7feeguAo0eP4nA4eO89q+NXWFiIw+Fg06ZNABQUFOBwONi6dSsAu3btwuFwsGuXVQli69atOBwOCgqs9VWbNm3C4XBQWFhIVFgI8eHCL//wvxw9ehSAP77xAQDJQVbSWr9+PQ6Hg4oKazn2unXrcDgc1NTUALB27VocDgeNjdaigxdffBGHw0Fbm/Xzzz77LA6Hg+5uKyH94fe/Q3q7XEN8Tz31FFdec52rxwNw3Q0fZV9ZE/GRoZzf/CE333zqf+ljjz3GV77wOevftr2b1atXs2LFCtf9q1at4tN33MnxAYsoKvZvITshkuc2F3PFj9/lmZffYUpKNInRYYT2WMN2iydbBy6uXLmSlStXun723nvv5YEHHmBultXb/Oi3nubhhx923X/XXXfxyCOPuG7fcccdrFq1ynV7xYoVrF692nX7459cwbZjtVziTLo33XQTY+WfCaq2zbXjXKkJLhsYOFte6rw20IgqssDYqrJMdGlRQnfEqWKlx5r6COpuIzfevjqdQb2dp/WguiMS6O41XDzNKvnTGxLJ/vIm5mTGMVQnQ0wvYSFBNLQP3ZtpDkumt89w2+IcBENSVyWJUaeeT1VYFvOdm5p7sRaCXD377KVOv3T5NABag8c39FmaspQ+hOvnZY7rcQBrg9dE+jr//PPN2TS2d5nJD75ifvFO0VnbKTVewDZj89871nDdrwfc/izws0FtngA+BKKBFOAwMONcj32u15K/uO/5HeaS1W+7bt/6i/fNJ5/cZOvvvPFn/zB3/PpD1+31BeVm8oOvmMfeOGgmP/iK2Vlcb2Z88zXznXX7hn2MJd/dYB56afcZ149Wt5jJD75iJj/4ijlR02paOrqNMcZ88XdbXdcnP/iKeWrjEWOMMf9vQ6FZ+O03TLOz3dnc/fut5pqfvDvap+vS3dNrZnzztVH/+w73WvK7HlSRcwOebtJVfqIUyB1wO4cziy7bVpHFHyTHhLvmVYwxFFW12P7+MDM9lr0nG+np7eMfh6td70vnTbJ6ch//+ft09vQx27mIYygJkaE0tJ1Z3WHrcWtZeUiQkJsU6ToK5L6r8vnIzFOl4D7v3B9231X5fPjwVcQMc2TIQKmxZx7AOBqHq1ro7OlzVfQYL79LUEc0QSn/shXIF5EpIhKGdWTNy4Pa2FKRxV8kRYfR1tVLR3cvta1dNLZ3Mz3V3veHK2el0dDWzaMv7+Ozv9nCj98oJC4ixLWcvN+sjOEXasRHhtI4RPmhwgprk/H+7yw/bRHCvOx4nv7XC1gyOZF7PzKNIGcFCREZ8hDGoaTEhFPf1j3m40L6Vx3OzXJPkW7bzoPylqLqFsKCg8hN1CXmauIzxvSIyFeBN4Bg4GljzL4B1VieNMYcEJH+iix9uKkii79IibHmZmpbuyhxruaz+wPsNXPSXYsW+k1LiyE5Jpz7rsrnl+8eIS8l6qyllhKiQl3HBg10sKKJednxQ65SFhH+/OXBJU9Hrn9fWVFVC3Oyhu/dDedodQshQeK2NQB+l6COVLWSlxLlqsml1ERnjHkNeG3QtScH3f4x8GNPxjVRJEVbxWVrWzpdQ23TbE5QIcFBfOOG2dz7xx2ua1NTrN95/zUzuP+aMw5tOENcZOhp+7fAGqI8UN7MtXPGXoD1bPoXVuw92TDGBNXKpKQoVwml8fK7d/FjNS2uPwSllOo/Vbu2tYuiqhaiw4LJ9MBBph9dkMnR799ArHPuZ1ra6KoqpMSEU9PS6aqvV9vSSXljB3WtXWcdGhyPvORoYiNC2F3aOKafP17b6tb6hn7Xg6po7ODyGSM7M0op5f/6h/jqWro4VNnM9PRY1/yM3YKCxFVRYtoo571yk6Lo7OmjsrmDD47U8u9/2u26b9ZZFleMR1CQsCg3gZ3FDaP+2ZK6Ng5WNHPFTPe9//pVD6qls4fWrl495l0p5XKqB9XJocpmZtpU4mg4n7toMnBqBd9IzXTOT720vfSMiuxjGX4bqcz4COpaR7+Sb81GayP0jQuy3BaLX/Wgqpwna/YfaKaUUjHhIYSFBLHteD01LV0szB1dohivh2+YzZeumHbWMkNDuSDPqvxQWt9ORVMHWfERXDs3g2lpMeM6Y+lchls9eC6FFc3Mz453nXrsDn6VoCqbrKyfPso/BKWU/xIRkqPDXLXvLshL8ujvjwgNHlPhahFhZnosta1dFNe18a8X5/GNG2bbEOHpEqLC6Ojuo7Ond8THkXT39rHnZAO3L3XP/qd+fjXEV9Vs9aDSdIhPKTVAckwYPX2GiNAgpp7jiA5fkhobzv6yJrp6+jxWXzQu0uqdDdeL2lPawBv7Kk679vaBKjq6+1g2xb3J378SVH8PSof4lFID9C81n50ZN6G2oKTEhHGywTpCI9dDCSremaCaBiSoN/dV0N7VS1VzBzc98T5f+sP2/jJbAPxpWwkZcRGu86fcxc+G+DqIDA0eUUkPpVTg6F/JN7iSg69LjT31YduupeWDxQ/qQe0ra+TuP2wnMz6Cf7tsqqtdaX07uUlRNHV0s/FQNZ+/dIrbk//E+SgxAnVtXSRFh437DBKllH+55bxs5mXH8Ynzc7wdyqgsGTBflhbrmZGhROchg7UtVv3C7SfqAShv7OA7r+x3tdtVYi1Ff6Oggp4+c8aZVu7gV12NhrZuEqPtW92ilJqYLstP5bL8ibc/8qpZadwwP4OleUke++Cdm2gNJf7i3SNcOzeDrcfrSYwKpaO7j/buXq6dk86b+yspqW9j05Ea/uPPewA433nelDv5VYKqa+067UwUpZSayEKCg/jFZ8736O/sP6Z9V0kDRVXNrNtdxvXzMrh2bjr3v7ibWxfnsKO4gSNVrUQ6i9CK4Kqq7k5+NcTX0NZFgiYopZQaMxHh1vOsMzFf2GKdlfnZCydzy3k57HrkGpbPyyA5OoyXdpTy4lbr/s0PX2VLLH6VoOrbukmK0iE+pZQaj89caFW/WLenjISoUC5yngTc3wH4ykes03cPVjQzPS3Gtq09tiYoEblfRPaJSIGIPC8iEYPud4hIo4jscn49MtxjnUtPbx+N7d3ag1JKqXHq3ytW2dTJ3Ky4M+a/bl6UzdevygcgOdq+91zbEpSIZANfB5YYY+ZhnWWzYoim/zDGLHJ+fWesv69/SWSi9qCUUmpcEqPDiAqz5pfmDXP44GX51qq97DFUyRgpu4f4QoBIEQkBojjzqGq3qXcejZxoYzZXSqlA8bEFmQDDVic/f1Ii3/34PB69ca5tMdi2is8Yc1JEHgOKgXbgTWPMm0M0vUhEdmMlrweMMfsGNxCRu4G7ASZNGrrWU05iJH+79xJy9CRdpZQat+/cPI9/v2YmGfFDzy8FBQl3OOeq7GLnEF8icDMwBcgCokXkjkHNdgCTjTELgZ8Bfx3qsYwxa4wxS4wxS1JTh87mEaHBLMxNIDlGyxwppdR4RYQGD5ucPMXOIb6rgWPGmGpjTDewFrh4YANjTJMxpsX5/WtAqIi4fzuyUkqpCcfOBFUMXCgiUWItAbkKODCwgYhkOO9DRJY646m1MSallFIThJ1zUJtF5M9Yw3g9wE5gjYjc47z/SeATwJdFpAdrnmqFGAFMVAAAA2lJREFUGVgiVymlVMCytdSRMeZR4NFBl58ccP8TwBN2xqCUUmpi8qtKEkoppfyHJiillFI+SROUUkopn6QJSimllE/SBKWUUsonyURb1S0i1cCJszRJAWo8FI636XP1rsnGmIl3TKvTOV5LvvjvbZdAeq7gm893yNfShEtQ5yIi24wxS7wdhyfoc1V2CaR/70B6rjCxnq8O8SmllPJJmqCUUkr5JH9MUGu8HYAH6XNVdgmkf+9Aeq4wgZ6v381BKaWU8g/+2INSSinlBzRBKaWU8kl+k6BEZLmIFIpIkYg85O14xkJEnhaRKhEpGHAtSUQ2iMhh538TB9z3sPP5ForIdQOuny8ie533/bT/zC1fIiK5IvKOiBwQkX0icp/zul8+34lEX0sT62/Lr19LxpgJ/wUEA0eAqUAYsBuY4+24xvA8LgcWAwUDrv0IeMj5/UPAD53fz3E+z3BgivP5Bzvv2wJcBAjwOnC9t5/bEM81E1js/D4WOOR8Tn75fCfKl76WJt7flj+/lvylB7UUKDLGHDXGdAEvADd7OaZRM8ZsBOoGXb4Z+J3z+98BHx9w/QVjTKcx5hhQBCwVkUwgzhjzgbH+4n4/4Gd8hjGm3Bizw/l9M9Zpy9n46fOdQPS1NMH+tvz5teQvCSobKBlwu9R5zR+kG2PKwfpDBNKc14d7ztnO7wdf91kikgecB2wmAJ6vj9PX0gT+2/K315K/JKihxkn9ff38cM95Qv1biEgM8BKw0hjTdLamQ1ybcM93AgjEf0+/+Nvyx9eSvySoUiB3wO0coMxLsbhbpbPrjfO/Vc7rwz3nUuf3g6/7HBEJxXpBPWeMWeu87LfPd4LQ19IE/Nvy19eSvySorUC+iEwRkTBgBfCyl2Nyl5eBzzm//xzwtwHXV4hIuIhMAfKBLc6ufLOIXOhcgXPngJ/xGc7YfgMcMMb8ZMBdfvl8JxB9LU2wvy2/fi15ewWKu76AG7BWrxwBvunteMb4HJ4HyoFurE8zXwCSgbeBw87/Jg1o/03n8y1kwGobYAlQ4LzvCZwVQ3zpC7gUa/hgD7DL+XWDvz7fifSlr6WJ9bflz68lLXWklFLKJ/nLEJ9SSik/owlKKaWUT9IEpZRSyidpglJKKeWTNEEppZTySZqglFJK+SRNUEoppXzS/wcyLrgcqSUo6wAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.subplot(1,2,1)\n", "plt.plot([0,num_steps],[9.14,9.14], 'k:')\n", "plt.plot(a)\n", "plt.ylabel('a')\n", "\n", "plt.subplot(1,2,2)\n", "plt.ylabel('b')\n", "plt.plot([0,num_steps],[0.6,0.6], 'k:')\n", "plt.plot(b)\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Note that SVI obtains parameters very close to the true parameters of the desired conditional distribution. This is to be expected as our guide is from the same family.**\n", "\n", "Note that optimization will update the values of the guide parameters in the parameter store, so that once we find good parameter values, we can use samples from the guide as posterior samples for downstream tasks.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Next Steps\n", "\n", "In the [Variational Autoencoder tutorial](vae.ipynb), we'll see how models like `scale` can be augmented with deep neural networks and use stochastic variational inference to build a generative model of images." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.6" } }, "nbformat": 4, "nbformat_minor": 4 }