Repository: cornellius-gp/gpytorch
Pull request #918 (2019-11-01 16:54:02)

Add fully Bayesian modelling via NUTS
Alright, the time came to take matters in to our own hands on this :-)

This is mostly done in my mind, two things left to do that I'll finish hopefully today or tomorrow:
- [x] Add `expand` methods to the rest of the priors (100% necessary because of how extending torch Distributions work).
- [x] Add a few more custom `pyro_load_from_samples` override to do additional cleanup. These will all be two liners that do one thing and call super. For example, we need to update `batch_shape` in the kernels.
- [x] Unit tests + docs + text explaining stuff in the notebook.

The interface is pretty clean in my opinion (see the example notebook). Here's sampling:
```python
model = ExactGPModel(train_x, train_y, likelihood)

# Register normal GPyTorch priors
model.mean_module.register_prior("mean_prior", UniformPrior(-1, 1), "constant")
model.covar_module.base_kernel.register_prior("lengthscale_prior", UniformPrior(0.01, 0.2), "lengthscale")
model.covar_module.register_prior("outputscale_prior", UniformPrior(1, 2), "outputscale")
likelihood.register_prior("noise_prior", LogNormalPrior(-1.5, 0.1), "noise")

mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)

def pyro_model(x, y):
    model.pyro_sample_from_prior()
    output = model(x)
    loss = mll.pyro_factor(output, y)
    return y

nuts_kernel = NUTS(pyro_model, adapt_step_size=True)
mcmc_run = MCMC(nuts_kernel, num_samples=100, warmup_steps=200)
mcmc_run.run(train_x, train_y)
```

And here's evaluation:
```python
model.pyro_load_from_samples(mcmc_run.get_samples())
model.eval()
expanded_train_x = train_x.unsqueeze(-1).repeat(100, 1, 1)
print(expanded_train_x.shape)
with gpytorch.settings.debug(False):
    output = model(expanded_train_x)
```

@Balandat @gpleiss @dme65 

Balandat: cc @qingfeng10
jacobrgardner: Test failures are just the incomplete Distribution/Prior refactor.
jacobrgardner: @Balandat @gpleiss @dme65 This is ready for review.

Patch:
diff --git a/docs/source/conf.py b/docs/source/conf.py
index e35bfcae..acfc9d85 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -74,7 +74,7 @@ class Mock(MagicMock):
             return ModuleMock
         elif "Distribution" in name:
             return _Distribution
-        elif "Normal" in name or "Gamma" in name or "Wishart" in name:
+        elif "Normal" in name or "Gamma" in name or "Wishart" in name or "Uniform" in name:
             return _SubDistribution
         elif "Kernel" in name or "Parallel" in name:
             return _Kernel
diff --git a/examples/01_Simple_GP_Regression/Simple_GP_Regression_Fully_Bayesian.ipynb b/examples/01_Simple_GP_Regression/Simple_GP_Regression_Fully_Bayesian.ipynb
new file mode 100644
index 00000000..a4c6dbd9
--- /dev/null
+++ b/examples/01_Simple_GP_Regression/Simple_GP_Regression_Fully_Bayesian.ipynb
@@ -0,0 +1,256 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Sampling Hyperparamters with GPyTorch + NUTS\n",
+    "\n",
+    "In this notebook, we'll demonstrate how to integrate GPyTorch and NUTS to sample GP hyperparameters and perform GP inference in a fully Bayesian way.\n",
+    "\n",
+    "The high level overview of sampling in GPyTorch is as follows:\n",
+    "\n",
+    "1. Define your model as normal, extending ExactGP and defining a forward method.\n",
+    "2. For each parameter your model defines, you'll need to register a GPyTorch prior with that parameter, or some function of the parameter. If you use something other than a default closure (e.g., by specifying a parameter or transformed parameter name), you'll need to also specify a setting_closure: see the docs for `gpytorch.Module.register_prior`.\n",
+    "3. Define a pyro model that has a sample site for each GP parameter, and then computes a loss. For your convenience, we define a `pyro_sample_from_prior` method on `gpytorch.Module` that does the former operation. For the latter operation, just call `mll.pyro_factor(output, y)` instead of `mll(output, y)` to get your loss.\n",
+    "4. Run NUTS (or HMC etc) on the pyro model you just defined to generate samples. Note this can take quite a while or no time at all depending on the priors you've defined.\n",
+    "5. Load the samples in to the model, converting the model from a simple GP to a batch GP (see our example notebook on simple batch GPs), where each GP in the batch corresponds to a different hyperparameter sample.\n",
+    "6. Pass test data through the batch GP to get predictions for each hyperparameter sample."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 39,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "The autoreload extension is already loaded. To reload it, use:\n",
+      "  %reload_ext autoreload\n"
+     ]
+    }
+   ],
+   "source": [
+    "import math\n",
+    "import torch\n",
+    "import gpytorch\n",
+    "import pyro\n",
+    "from pyro.infer.mcmc import NUTS, MCMC\n",
+    "from matplotlib import pyplot as plt\n",
+    "\n",
+    "%matplotlib inline\n",
+    "%load_ext autoreload\n",
+    "%autoreload 2"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 40,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Training data is 11 points in [0,1] inclusive regularly spaced\n",
+    "train_x = torch.linspace(0, 1, 6)\n",
+    "# True function is sin(2*pi*x) with Gaussian noise\n",
+    "train_y = torch.sin(train_x * (2 * math.pi)) + torch.randn(train_x.size()) * 0.2"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 41,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# We will use the simplest form of GP model, exact inference\n",
+    "class ExactGPModel(gpytorch.models.ExactGP):\n",
+    "    def __init__(self, train_x, train_y, likelihood):\n",
+    "        super(ExactGPModel, self).__init__(train_x, train_y, likelihood)\n",
+    "        self.mean_module = gpytorch.means.ConstantMean()\n",
+    "        self.covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.PeriodicKernel())\n",
+    "    \n",
+    "    def forward(self, x):\n",
+    "        mean_x = self.mean_module(x)\n",
+    "        covar_x = self.covar_module(x)\n",
+    "        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Running Sampling\n",
+    "\n",
+    "The next cell is the first piece of code that differs substantially from other work flows. In it, we create the model and likelihood as normal, and then register priors to each of the parameters of the model. Note that we directly can register priors to transformed parameters (e.g., \"lengthscale\") rather than raw ones (e.g., \"raw_lengthscale\"). This is useful, **however** you'll need to specify a prior whose support is fully contained in the domain of the parameter. For example, a lengthscale prior must have support only over the positive reals or a subset thereof."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 42,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "sample: 100%|██████████| 300/300 [00:15, 19.16it/s, step size=6.34e-01, acc. prob=0.910]\n"
+     ]
+    }
+   ],
+   "source": [
+    "from gpytorch.priors import LogNormalPrior, NormalPrior, UniformPrior\n",
+    "# Use a positive constraint instead of usual GreaterThan(1e-4) so that LogNormal has support over full range.\n",
+    "likelihood = gpytorch.likelihoods.GaussianLikelihood(noise_constraint=gpytorch.constraints.Positive())\n",
+    "model = ExactGPModel(train_x, train_y, likelihood)\n",
+    "\n",
+    "model.mean_module.register_prior(\"mean_prior\", UniformPrior(-1, 1), \"constant\")\n",
+    "model.covar_module.base_kernel.register_prior(\"lengthscale_prior\", UniformPrior(0.01, 0.5), \"lengthscale\")\n",
+    "model.covar_module.base_kernel.register_prior(\"period_length_prior\", UniformPrior(0.05, 2.5), \"period_length\")\n",
+    "model.covar_module.register_prior(\"outputscale_prior\", UniformPrior(1, 2), \"outputscale\")\n",
+    "likelihood.register_prior(\"noise_prior\", UniformPrior(0.05, 0.3), \"noise\")\n",
+    "\n",
+    "mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)\n",
+    "\n",
+    "def pyro_model(x, y):\n",
+    "    model.pyro_sample_from_prior()\n",
+    "    output = model(x)\n",
+    "    loss = mll.pyro_factor(output, y)\n",
+    "    return y\n",
+    "\n",
+    "nuts_kernel = NUTS(pyro_model, adapt_step_size=True)\n",
+    "mcmc_run = MCMC(nuts_kernel, num_samples=100, warmup_steps=200)\n",
+    "mcmc_run.run(train_x, train_y)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Loading Samples\n",
+    "\n",
+    "In the next cell, we load the samples generated by NUTS in to the model. This converts `model` from a single GP to a batch of `num_samples` GPs, in this case 100."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 43,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "model.pyro_load_from_samples(mcmc_run.get_samples())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 44,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "model.eval()\n",
+    "test_x = torch.linspace(0, 1, 101).unsqueeze(-1)\n",
+    "test_y = torch.sin(test_x * (2 * math.pi))\n",
+    "expanded_test_x = test_x.unsqueeze(0).repeat(100, 1, 1)\n",
+    "output = model(expanded_test_x)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 45,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "torch.Size([100, 101])"
+      ]
+     },
+     "execution_count": 45,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "output.mean.shape"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Plot Mean Functions\n",
+    "\n",
+    "In the next cell, we plot the first 25 mean functions on the samep lot. This particular example has a fairly large amount of data for only 1 dimension, so the hyperparameter posterior is quite tight and there is relatively little variance."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 54,
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQcAAADFCAYAAABZ/DesAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO2dd3hUVf7/3+fe9J5IAoSahCo1ZDCA1CQgINKWbl1R3FV3XX+sCOouj66soqJfEXVBcBcVEKWD1IQWOgMJAUIJTAIEUkiZFNLvPb8/PpmWuQkJJBDwvJ6Hh2RmcufW9zmfehjnHAKBQFAV6X7vgEAgaJwIcRAIBJoIcRAIBJoIcRAIBJoIcRAIBJoIcRAIBJo43O8dAIDZs2eLeKpAcB/5+OOPWdXXGoU4AMD7779fq89lZmYiICCggfdGUBfENWl81OWazJ07V/N1YVYIBAJNhDgIBAJNGo1ZIfh9k5+fj9zcXJSXl9/vXXkoUFUV+fn5Nq8xxuDt7Y0mTZrUahtCHASNgps3b6JFixZwdnYGY3a+MUEdKS8vh6Ojo81riqIgJSVFiIPgwYJzDhcXl/u9Gw81siyjLoWWwucgEFgRFxeHuLi4Bv8eo9GIdevWNfj33A1CHAQPHGlpaYiKikJ6evodbyMuLg5Lly5FTEwMli5dCoPBAADw9vbGmjVr6mtXq8XHx0fze+Li4tC5c2esW7cO69atw4IFC8z7pkVN790twqwQPHB89NFHOHToEP79739j4cKFdf57o9GITz/9FCtXrjS/Nm3aNKxcuRJ+fn71uas14uvra/daaGgogoKCMH78ePNrI0eOxNatW+0+azAYsGzZMsybN69B9k+Ig+CBwcfHByUlJebflyxZgiVLlsDFxQVGo7HW21mzZg0iIiJsXvP19UVMTAzCwsIQFxeHmJgYxMfHY/r06Thx4gQAID4+HuPGjcPu3bvh5+eHoKAgJCcnY82aNQgKCkLHjh2xfft2rFy5Eq+99hpmzpwJADafDwoKwrJly9CzZ0+cPHmy1sdtmiHs3r0bABAREYH4+HgkJycjLi4O3t7eNu+1atWq1uejOoRZIXhgOHfuHCZPngxXV1cAgKurK6ZMmYLz58/XeVt5eXnVvhcaGorIyEj07NkTy5YtQ3x8PHbv3o2IiAi89957CAsLMwtDREQEfH19MW/ePDz33HPmbYwfPx7BwcF2n3/33Xcxbtw4REZGIigoqE77HBwcDD8/P/j5+WH9+vWIiIhAUFAQQkND7d6rD4Q4CB4YmjdvDi8vL5SWlsLFxQWlpaXw8vJCs2bN6rSdiIgI82zARHJyMiIjI21eM5kY48aNw/Tp0/Hpp5+irKwM3t7eCA0NNc8+fHx8bLa9YMEChIWFmV+r+vm6YjQaERwcjAULFsDb2xs9e/Y0vw6QeVHde3eDMCsEDxSZmZl4+eWXMX36dCxbtuyOnJLBwcF46623sHTpUgQFBSE+Ph6LFi0yv280Gm3MCpMZEBERgYiICCxbtsw8apum9UajET4+PpgwYQLeffdds2B8+OGHNp+fOXMm1q9fj549e5r/NjQ01PzdcXFxSE5ONkcykpOTzftm+r68vDwYDAbk5ubCaDQiOTnZ7r2UlBT4+/vfzakGawwNZmfPns1F4dWDS31ck0uXLqFdu3b1tEcCrSQoQPs8z507V7MqU5gVAoFAEyEOAoFAEyEOAoFAEyEOAoFAEyEOAoFAEyEOAkE9cLuCrQeh0KoqQhwEvzuMRqO56Op2hU21JSgoCEuXLq32/eoKrQASlr59+9qIi8FgsHvtXtNgSVA6nW5W5Y+9ARzX6/WfNNR3CQR1wVRbERwcXG/btM6SrA6tQiuA0rVNSVlff/01AErvNqVG3y8aRBx0Ot1ivV7/itXvJ3Q6HYRACG5HWRlw5crdd4Jq04bDyUn7vYiICLz++uuIjIzEuHHjEBwcDIPBYFO4lJycjAULFmDmzJnYvXs3wsLCkJOTg7CwMHOx1YQJExATE4N58+aZxcG0nboWWnl7e1f7XtVtVi2yqrqv8+bNw549eyDLMuLj480FYHWl3s0KnU7nA6BqYvdiAHPq+7sEgjshODgYW7duRc+ePfH666/DYDDYFS5FRkbCx8cHkZGRMBqN5vRpUwGWqax6/PjxWLZsmXnbd1NoNX78eCxdutQupbrqNm+3rwBsisXulIaYOfgBmFU5e7A25m4/7xL87nFyAtq3b9iU/qVLl+Kll15CZGQk/Pz8kJycbK53MD2ARqPRxgwwzS5qg7e3N4KDgxEUFISYmJha71dkZCSmTZtmU7Sltc0FCxbUuK8AMHbsWDg4OOC9997DN998Uyuzpyr1PnOoFISwKsIwFEB0fX+XQHCnmDot7d692zyqJycnIz4+HgaDAVu3bkVycjIMBgOSk5PNhVimak5Tz4fdu3dj+vTp5oKpt99+G8uWLUNMTAxycnLMhVYxMTHmQquqxMXFYdmyZTAYDOYZgin6ERcXZy7eMm3zdvsaFxeHDRs2mEvK70QYgHtQeFVpZiTDXjDMzJ49m7/22mu12p6p+k3QeKiPa2I0Guvc3+B+YTQasWDBAvzrX/+637tSLYqiQJZlu9eTk5PtrtXXX39935bD+xVAZHXCYKIuVX2iKrPxcbfXJD8/X7OKsDGSkJCAhISERr+/WvsnSVKtr1WDioNOp5sPYL5er69dPyyB4AEgMjLSrjHMw0iDJUHpdLoJAHbp9froyt97NdR3CR58GGN1WlNBUHc453VaMKih8hyiQFGL6Eqfgx+AyQDEDEKgibOzs9l3IVa8ahhKS0vh4FD7R77exaFSDHZV/rrY6q2GXwxA8MDSrFkzpKenIycnR8wg6gFVVSFJtoYBY6xOrePqXRz0er0RgJB+QZ1wcHBAy5Yt7/duPDTUR+s+UXglEAg0EeIgEAg0EeIgEAg0EeIgEAg0EeIgEAg0EeIgEAg0EeIgEAg0EeIgEAg0EeIgEAg0EeIgEAg0EeIgEAg0EeIgEAg0EeIgEAg0EeIgEAg0EeIgEAg0EeIgEAg0EeIgEAg0EeIgEAg0EeIgEAg0EeIgEAg0EeIgEAg0uRfL4Qms4Bw4fpwhPZ3B1DmcMUBV6Wd3d46BAznqsLyAQNAgiFuwHklLS8Ozzz6Ln376Cc2aNbN7PzGR4exZhsceU/HYY9prMxiNwLZtEsrLgcGDVfj5NfReCwTaCLOiHvnoo49w6NAh/Pvf/7Z5nXNgyxYJRUXAxIkq2rSpfhs+PsBTT6kYO1bFyZMMa9dKKCxs4B0XCDQQM4d6wMfHByUlJebflyxZgiVLlsDJyQXr1+dhwwYShlu3JJSUAOXlQEUFoCiAJAGenoCXF4e3N+DmBjg5kaDIMqDTqdi9W4KjI/DEEyokIeeCe4QQh3rg3Llz+OtfZ2PHjk0oKyuGLLvC3X0cvLw+wfPPO8DVlcPPj4TA5FsAyNegKEBmJgAwuLsDHh6Alxfg6cnRqRPHsWMMaWkMnp7AsmUOGDRIxeuvqxDLSQoaGiEOd0FKCnDwoISYmFY4e9YbZWWlYMwFilIKV1dPeHg0Q1AQCYOjI+Djw+HqSqIgywBjHJJEAlFUxHDrFpCXx5CbC9y4wXD8OIOrKxAUxFFQwBAYyHHqFDBihAM6duQYM0ZFRIRYV1LQMAhxqCOKAuzdy7BvH8O1awy5uQyXLwPXr2fCzW0GfHxeBmPfobQ0HWFhHC1bkhhwDjg5cfTqxSudjNw8i3B0ZCgq4rh0iSE1lUwOACguBq5fJ8EoKACSkxkcHBhatOA4f57h7FkZq1ZxBAUBf/2rAg+P+3VWGhcXLzJcvUozsOocv4LbI8ShlpSXAzt3SoiOZigpARITgStXGFQV6NgRCAxcg4oKDmdnICtrERwdaWZx+TLMJgDnDGvXknnh4QEEB9PnS0oAZ2fA15f+lyQSE09PoEcPjqwswMGBPmcwMGRmMpSWklAlJDAkJgLHjjmgQweOv/5Vwe95PdrsbBKHqCgVZ84wXLzI4OHBceMGXYSwMC5MsloixOE2qCoQHc2wZw9DdjbD6dMMKSnkA+jSRUVgIBAbS4LRujXlLhQWMjz+uIqFCxW4utpvs7CQZgFbtzIUFDC4utIMIzMT8PcHHnmEPufsTOKRnAxcvQo4ODD07q3i6lXyT1y9ypCXBxQVAefPA6mpDAcOOECn43jrLQWtWt3bc9UYiI6WMGkS+WTCwjjeeUfGiBEqOnTgKCkBVqyQMHq0Ci+v+72njR8hDjVw/DjDgQOUm3DyJENqKoOHB/DYYypKSxny8hiys4GWLYG//11BeTnDwYMML72k4pdfGPR6hgED7Ke1Hh5At24c3brRe0YjcOCABDc3oLTU8jmjkSE9HWjaFAgNBTjncHDgCAkBkpJISFJSyJ9x7RqZHzk5wLZtDDExDujXj+PTTxX4+NyjE3afOXuW4dFHudmMW71awgsvKEhJkdCsGZ3r1q1VrFghYdw4VZhht0GIgwbp6cCuXRKOHmWIiSG/gpsbR8eOHGVlAMAwdaqCsjKGGzcYysuBRx/lSE0FAgM5nJwY2rTh+PBDGW++qWD48JrtXh8fYNQockBcvQrEx0vgHAgI4Bg0iOPcOZphAGTKKArQoweQk8Ph6Qnk5wNt2wJnzwKtWtHUOi2NITqaoXt3R4wYoeLbbxVI0u0TtR4UOAd+/VWCuzvNxCZNUnH+PMMf/kDncedOCYMHq2jWDDh3jmaAGRlA8+bAtGkqfvlFwtSp6m2+5feNEAcrVBXYvl3CqVMcy5fLyM0FfH05PD05fHyA0FCOBQsUODqSqbF2LUOXLhyvvKLCaAS++IKmsKGhKoYMAV59lSE6WgJjKp54onaOsdataXQDaEZx7BhDaSkJgyQBffpwyDJHaqqERx6B2fdQUQE0a0ZRkKZNgSZNOHJyGPLzOXbsYGje3BFPP62gvNySqLVw4cIGO5cNze7dDIMHqwgIIHNs7VqGpk3pHCcmMvj5cZi0r08fFd99JyEwkM7TuHEqevfmOHyYoW9f4bCsDiEOlVy+DOzbJ+GrryRcvcoQEMDh5QVIEsObb1ZgxgyLIyslBVi7VsLkySpycxk4B+bOlfGHP6iYMkXF9u0SOnZU8eabCr75RkZyMsP8+RK6d+dQFKqfaNWKHJI1JTX5+ADDhnEAdANzDqSmAikp5AiVJDJpAgKAy5cZZJn+T0sjU6N5cw5fX8DHhyEmxh3ffGOfqOXi4gKj0dhwJ7YBoIQyukYAHf+FCwzh4SoqKoAzZxgmTbLMCi5dIp/QU0+pSEsD1q8nsyIhQUJ5OYej4/06ksbN714cVJUe9OXLGQ4ckODhAXTvruLqVQnTpyt4+20Vsmz5POfAggUy+vYl8QgOVvGf/0jo21fFoEF0sw4frmLdOgleXhzXrzOkpADDhytITZUwbZqK4mLg2jWGbdskKAptt0sXFSEhNe8rY2Q2tGplEQwtSkpoZP3vfyVcu0bi1atXEpKSZqGgYANUtRiMueLRR8fiL3+Zj7Q0mm4/KGzbRk5Fa7p0AQ4fluDiAjz5pOU9zoGsLIZOnSg03Lw54O/PkZFBGafR0RJGjBDmhRa/a3FISAAWLZKxdi2DqjIMHKgiO5vB2Rk4frxcs+hpwQIZ7dtTAtK+fRJCQjgyMhhatKCRGgAOHmQ4dYrhlVdUTJyoIi+PITNTwpNPqvj1V5pxdOpEGZAA3cCnTzNs3EgPcnAwOSvvNOTm4gKMHMkxcqSCsjLg668lbN3aHIriBVUtBeACzktx8aI3TpygubePD+DqCgwapMLF5c6+915QUUHhXuvR/vp18vVkZpIAuLtb3jt8mKFPHxXe3sD+/RKeeEJFv34cmzZJGDtWRUkJmWXWA0BDc+sWRVVkmdLmw8I49uwhP5Ppmru7k7/pfqbL/y7FISsL+OwzGT/+yJCfz9CrF9mge/bQAz1jhvZIsmcPgyxztG4NeHvTaz/8QOZC9+4ceXnAb79J6NdPxTvvKIiJobDanDkyOnXiyMxkGDRIxebNEp56yvIdjMG8DYBMnM2b6a6QJKBNG3KGOjnZ7xPnZHNfvUr5D6YbzPS/ogC5uQzNmnE4OGTAz28GnJxmIDd3CSQpHStXyti+neP//T8VQ4eq2LOH6kAee0xtlKHQQ4cY+va1vT56Pc0kVq+W4OwMWM+qMjMZ+vWj303lL6Y8Es6BIUNU7N7NMHTovfM9bNggYcoUqpP5738lHD7MMHOmalOmn51Nn/P355oRr3vB70oc0tMpvPX11xIyMuiBefFFFZcv08i9a1eF5mzh9GmGQ4cYjh1jaNKEpvUbN0rIz6ciKRcXehCjoy0XHQDKykhEAgM5cnMZEhKAiRM5Hn1URUwMQ2Sk9kUPCQFCQugB4By4coWyMsvKaFgxPfymn/39OVq1ohHIeqQpKqLj7dKFw80NeP/9n5GYyLBzp4TY2K9w5YqEgACKwLz7royjR4GxY4ExY1To9QxHj5JwBgfX2yW4a3JzGR55xP68GQxA+/YcV69aBDInB/D2tnzW25vDaKRZUni4iiNHyCGZl8cQEwPk5zM4OpKpWV5OwnG3JfOpqcDJkxICAzl0Oo7du2mGKknAL7/QTCYjg/JnQkMt+/rII8D48SquXAFWrpQwZoxqMyO6F/wuxCE5mezUH36QcOkSg5sbPQAdO3LExkoYMIDjnXcq7KbxmZnAnj0S2rWjkOHQoRweHhR2/PBDCYrC0L07R1ISg68vN4fRTHTqRGHInj2Bo0cpzBYfz9CzJ0dBAeVR9O5d86jAGIUp27at2c9QleJiuqnOn2dwcqKHoaSEMisDAoDp0zl27eI4ehRwcqKU7M2bZej1HAYDMGIER58+Ko4fp3yNoUNV+PrW+usbBFMVqzWUOMZx8qSEiRNVBAZSFKJfP44DByQb/0O/fhwxMeRjCAwEjhxhOHUKuHmTHsZx4yyfVVVgxw4JzZtz9Ox5ZyN3dDQly40apeLsWYbVqyX4+VHC2+rVEqKiVDRpArRowfHrrxJ69KCUeusZRJs2QMuWKtauldCnj4rWre9oV+6Ih1oc4uPJ6XfiBEdsrAxXV46OHYGxYxVcvMgQH8/w/vsKwsPtL35sLENxMcXPV66U4O1N8fQnnlChqpSgFBWlYssWCeHh1L2pKp06cWzeLGHAABUnTsi4fp2hqAjo2ZNuuP37GS5cYOjYsX6mjWVlwL59DMXFDD/8wCDLDMOHq0hKYmCMISuLHpjwcI4hQ8gf8sYbDjh7ljz2bm5ATg7D/Pkyrl1T8PjjDGPHqggL49i1i57KYcPuX0XokSMUkbDm+HEJrVtb/DfNm9PnAG4uezfh5ITKPBX6/8QJhqAgjj//WcX69RKsxVeSgBEjVOzde2fXKDsbUBTL/nbrRhW23t4MmzeT2JoyYQEgIkLFu+/K6NWLm/NbTMgy3Ye7dlEdT2GhZQZpQlUBDw8yQa5dAxITnTBqVJ122Y6HThwUBdixg+HgQQm5uRxbtjA4OjK0bk1FT/7+HAkJDG3aAF9+qcDNzfbvi4vJ1uvTR0VQEBATQw6tY8ck+PqSI+z77yUMG6bi8GEJTZpwtG9/+yQnilzQw3X0KDMLyo4dEpydOdq2vfNjzsig6lBnZ4qefP89jaJhYTTd/sMfSCC+/VZCRARHRIRaOWNSsWpVBV56ScblywwdO6o4eVJCbi7w448y4uLIqz9ihIonnlBx8ybNRiIi1PsS3cjOZnj8cdtzrSgUxpwwwSIabm7kpPT1tb8ubm4cZ88Cp09LeO01BUlJ9IR5eHDk58MurXrwYI4tW2jE9/ev/b7u3EkmponycpqxGQxA167cRhgqKsjHFBrK8eSTlDOzcSOZEibi4igRLjsb6NDBVjxMGI3Ad99JyMlheOaZstrvbDU8NK1DsrOBb7+V8N57MmJiGNatY9ixQ0bHjvTgjRlDDp+sLMqQ+/BDe2E4d45hyxYJEyaQMKSkkGpfuECe5KFDVRw4QBmKZWUMzZtzhIZS2nJGhvZ+mcJmHTrQTZCXx8zhRYBmIgkJEm7cqPsxp6cDa9ZIuHCBYdw4Fa6uHIcOWR7edu04iopoanv2LMPEiRxTp6pISCDbevduCYwBL7ygon9/FW5uFIbt04duysREhs8/l/HzzxJ27WLw9weeflrF6dMMhw/f2+kD19Df0lI6B507277Zv7+K5ctl9Olj+zrnVBpvenADA+l6ACQC+/ZpPw5PPqmaZ061IS6OWgFaj+x79jB0705dwK5ftz13e/YwjB+vYtw4ypEJDAT8/GgwOXGCYc0aCouPH6/i5Zcp8WvrVvv9uX6doVMnjuHDVWze7Gb3fl1pMHHQ6XSzdDrdBJ1ON0On081oiO/gnBx1//ynjEWLJMTFMezcybBli4zOnVX4+1Nac7duFG4MDAQ+/VTBiBH2o4+pW9PEiSocHWnaefiwhMcf58jOpof84kW6qFevMvj5qZWOPo5hw1Ts3at9Knv35tDrKYKRm8tgMFBR1oEDlhtk9GiamSQn1+64s7JIFJKSyI8RFsaxciVNrz08yOQZOJD8HZ06cezbx8zixhgQFUWzh2vXGDZskNC/P5lbQUEcbdtyhIcDzz1H4pmbC3z1lYwVKyT89790joYNo7LzdeskzYe2IUhKYujQwfbLDh1icHEBunSxfd3dnZyx1tGdmzep6KpvXxL+qjg6Wkrlq8IYMGCAiv37ayeIV64wu5yVW7cYDh2SMHKkisceI18OQPdwQQGDlxftg4cHR24uObIXLJDh48MxYYJtDkznzhw9e6pYtcqSJ5OQwJCVRaZsRgbD5MlFtdrXmmgQs0Kn080HcFyv168x/a7T6SaYfr9bkpKoojEjg0KLiYkM16/TA9ChA9lcnNNDWF5OeQvPPaeia1f7OzkpieHECYZRo2wLcUxZdNHR9AC0aAHcuGFKqFGRnS3BxYWjb1/KR9DpVBw7xuz6B8gyiY+vL5kW2dmAnx+F2CoqLF2mx44lm/LWLWjuJwDk5VFExM+PnJ+MAadOkeBMmaJiwwa6WcaPp5H/wgXyGZw5Qz4Ta2eehwfw1lsKPv9cwpo1Enx9OW7doszLkBAVqsqgKBTJiY2V8NtvDOfPU1+JkSMptOrvz/HjjxIGD76ry1krLlxg5voTE+fP25sZAJ1vFxfKJ3B3p4Swigqa9TBGpkVREZkfwcEcly9ThKhbN47Tp5m5IM6aVq3oAczP5zVWdBYXwy5PJCGBQVE4wsPpXgkJAU6eZNDpSLgHDLAcV8uWwLx5Ml5+mWYJRUXkP6lKYCDNhletktC0KTUR6tSJY/9+Mimpu9jd0VAzhxlVhGAXgFfudqM7dzJ88IE3PvtMwr59DElJFF4sLiaPr7s71fL7+gIlJVRa3bw5OWlu3KACJRMVFSQAubnAlCm2wnDgAINOR2aIwUAjUFoa2byMcYweTc6u4mJm/ruQEMp6LC+3328HB7I5g4Lohty/nxyF27fbnv6hQzlKS2m/rEex0lKa2Rw9Sg/7kCE0uqxeLcHBgbzsKSlAQQGVKbu40Pc5OFCUxMWFo1cv7dyNMWNUtG9P9vSFCwxGI0NCAoXeXn9dQceOVHHKOcOFC8DChTL+9z8qYQcAnY7jm288sXMnmRqnT5OJVd8zCusEIdPvV64wc26INXo9wzPPKNi8WcKKFRI6deIYNsySVDZgAD1EAM06zpyhn9u1o4Y71TF8uIqdO2t+ZPbtkzBokL2IOTgwc0j45EmGwkLgs88k7Nwp4dgxmpXNmSMjI4OiGz4+HJ07c7NPRAs3NyAkhCMxUUJmJnDkCAlDfVHv4qDT6XppvJwDIOput71zZzrWrRuJwsJ0dOxII6mPD43OFy+SL6BVKzqp77yjYO3aCsyZo+Kpp1QMGMCh19M0+vvvJaxcSTHmqiP9pUum0ZMqMy9epOlsWhp5mF1cGBwcUJlsY8vIkaqmLRgeTk7Ivn1VXLnCUFhI/RicnOghtyYsjMyUzZslbNkiYds2hq1bKfwWGUmmwvr1Ek6epPoB05R6714JAQEwe9X37aN4+o0bDBUVrNoQWEgI9ZYID+d44w0FFRUUOlu6VMKePRJ69aIHZvJkFZ6eVLfx7bcy/vEPGYsWUShg6tQipKYyeHpSkVpKCu3zpk0SNm68M3+KNaWlsKt/2LiRCq+0yMhguHyZBo9p08i3YI2rqyUhqmrkxSTkWsgyJaRdvlz9vhYXw6aHR0UFze6efFI1+w98fDj++EcVsszw9NMqSkspR8PkBzMaGRYskGE00vdVZ25GR5NZ9eijlOlJaflUeTtx4kSkp6dXv6O1oCFmDn4gMbDGCAA6ne6uOgskJs5DTs5hxMV9iIsXqQlrUhKFmtq0AZ54gmPxYgWffKLYZZW5ugJDhlDILiSEo2lTevi3b5eQlUWfyc8nZ1LXrpReu2KFhKlTFZw5Qw4jU1bk4cMMrVqpaNnS/js8Pbmdc/KRR0zJO9RHsqiIIgxRUdxu9gDQVPjxx1VkZ8OcaBUdLWHHDqodGDdORVSUZSQ8eJAhPx82NQKFhQyMAa6u/Lapwc2bk/OrbVu6Sb286KblnATvyy8VhIbSzTtmjApFIQfnzz9TWM7VVcWLL6q4dEmCJAHh4eR1Hz2a/iUn00ORmnrbS6wJmWu2QnDkiKRZCp+fT8LYuTNHv34qCgq0t+ntbRFmNzdubv8/cGDNvgWTD0mLa9dgd0/89hslkq1fT/06JkywJJUVFHBs2EDC0b8/XafwcI6xY1W88oqC5cslpKYyLF0q4/p1y2yMc4oapaaSCIaGckybxjFliorUVIann/4Yx44dt1sioa4wXs/zP51ONwHAd3q93tfqNR8AuQBC9Hq9oerfzJ49m7/22mvVbrN9+/Yote6CYsYFnToZ8dln2ejeXalxvwoLGbZtc8XQoSXw8bHcaKWlQEKCE27elHD0qAv69aP3r1+XkZUlo2VLBSNGFFfa907w91eQkOAESeKIiiqxS8rhHNi0yQ1jxtg6hHbudMGwYSXYv98ZWVky3Nw4hg8vxtWrMtLSHBAebjm+EyeckJ8vYfDgktvmFKgq8MEH3pg5Mx+enqbmMRLOnXPErX1NEpsAAB8wSURBVFsMTk4c3bqVw9e3+ukm58DWra548slilJYCH3/sjb//PR87d7pW5noUobSUYetWVzRrpuDAASesXOlh7m3RqlUpevdW0L17GUpKJAwcWIyQEPvrceiQM0pKGIYMuf1xaZ07EykpMnbscMUrr9gu6HHhggNOnHBG587lCA0tQ1kZsH+/C6KiSqpuEuXlwJ49tN2iIga93gkDB9I12LbNFSNGFFe7P2lpMlJSHNC3r+09uWOHK4YNK7Y5tn/8wxtdu1Zg7Nhb5tlmRQV9R2amjIAABaNGFWuej5gYF4SFleHkSSe4uqooLJRw6xZDfLwT+vQpxeDBJXBxAdLSJOzb54o5cwKhKPbPibOzM5KSkqo9nq+//hoff/yx3R40hENSq/7XlIRadUZhJiAgoNoNnj9/HrNnz8avv26CqhZDklwREjIWUVGfYMgQhsJCfyQmcgwerC10BgOlsP75zyokyb79T2Ag8NNPEr78UoWTkxNKS4GJEx3wz38q0Ok4AE8AtO5EVJSKK1coBNismfYDN3AgQ2qqJ3r1suxP164MxcVeGD0a+OADuTLnwhMBAZSJxzlFV9avl9C1K680D27fy+zzzyVMn84REmKxc44epczALVtoPzt2vL0dGhgowdWV9ucf/wC+/LIZPvhAwY8/SoiN9cLkySpefJHOU7duwJgxwK5dlNKdkeGEy5cZevRwgYcHsGWLJ4KCaHYWEWFJ+x07lqIG27b5YNKk2hd4eXlJCAiwHMOqVTImT1YREGAJ1+3dS52627Uj888k2u7utn9rjfV7CQmWn9u2ZZAkTzRpor0/AQHAlSsS3NxsfVVeXhKaNvU0/37gAENZmYRXX1Ugy3QSsrOpZ0jbthyjRpFDOi7OC1FR9vfulCnAzz9L6N9fhV4vISiIqnwnTuQoKHDFt9/SzCEnhxL0vL0vobBwFsrKNgAohqurK8aMGYOPP/64xuerOhrCrMgBUNV88AEAvV5/R40DmjdvDi8vqiiUZaoo9Pf3whdf+GPgQF6ZsUhTrarhqLg4huRkSpIx3TDXr9PnAZo5/PQT5TY4OtJN9uKLDnj//YpKYSBMqbsXLjC0b19ztVzHjrzSt2B5rUcPjtOnJTg5kV/E1ZXj5EkS62HD6EFetkzC0KFqrbPxDh5kUBSGsDD70GxWFjV8qW1VX2QkRWYAahrTpQt5wlu0oA5XixZJ2LaN/BqcA6++qmLaNA4vL46KCoZz54CFCyXk5FBTHEdHoH17Oq7//c8S8vT3B6ZOVbFmjVQrj3pqqqXaFSDfEufcHNbknAS1VSuKBpj6XJh45BGOmze1t92+PceFC5YB07SP/ftzHDxY84l78kkVv/1m+QyFHy37WVwMfPqpjG+/Vcxm3YUL1Bbg2WdVFBQwNG1KJmdeHjOHJK1hjM7VunUSPv9cxm+/McTFUY7DwYMM2dl0bS5eBC5dApydm8PFxQtAKWTZGaWlpfDy8rrjjl/1Lg56vf4k7GcPfgCi72a7mZmZGDDgZQwatBvdu8+ALKdjzRoJiYmUANS1Kz0UP/xgif3u20chLFOB0+XLlB9w9SqlTv/wg4SFC2UMHaoiJkbC+vUSbt4E+vVTERpq+/2masALF8i+r877b2L0aBWbNllOr3Wx1IABHAYDw9WrdGNev07hL1nWdnRqcfQog15Po5I1KSnkxDp6lFbJut1+mnBwIO+3KaIzbRpFPTp2pAzQ/HwK+40Zo6JnT45ffpHw1lsKNmxQ0KdPCRwdGW7dYli8WMayZeTH+eorGYxxBAZyvPmmjPh42rajI4UVDxyQcPVqzfsVFyfZiF98PDNnkyoK5S4MGEB5AKWlsKtc7duX48gR7du8a1eOs2fpGnTpoiIxkX42hZ9rsrhlmSppKVWbcmJM1Z+KAnz5pYzBgy2zo0OHGG7epOhQWZntfg4dah+5unUL+OorCW+/TY7JP/1JQV4ew+DBFA7fvJkc0wYD3TO9epGfx9MzA488MgN/+ct2vPzyy8ioLjuvFjRU+vSSKnkNQwEsvpsNrl69GooCPP10Bf785y/h6kon8NIl8kq3aEEZZCtW0EjVtCmNDB06UCn1oUOUFxEczJGdzXD2LCWmAJSLP3o0RRSOHKG02qpkZzP4+5MTMCOD2WXfVUWWyam4fbuE4cPpAXV2pg7Ijz7K8dNPlLR0+DDNMJ5+mm6a1aup92F15dKKQqm1rVtTwlLVSr34eEq7vXKFVvLWqhupjqgoFRs3Shg/nmZZnTpRhuXo0RwdOlBTnOvXqVho3jwFv/5KyV3/+U8O/u//mmP9eoaSEobTp4HERBnPPquguJghLEzFF18oeOstGe7uCtq3J7EcP55mFqWl1aegc26ZCcTHMzg6Ul1KURH1kJwwwWK2HDnCzNmd1tehuuQmgHJP8vIoarNpEzNHf3r1orVKq87KrOncmSM2lhzi5eWWB556W3L06UMzmU2bJHTowPHoo5ZIUlgYdaXKy6MO4gkJQG6uBC8v4MwZCsW3a0eRlh49KAr07bcO2LjRAYpC2+nWjbqTe3lRK8GMDIZu3dZg9GgVspyH55//8naXvEYaJM9Br9e/DSC4MkNyFoDL9ZEAdeAAw7/+lYs9e6h4aMwYFVOnqggKol4Js2bJ0OsZVq6kpB1KlqKwZW4uhfYKCymMNWWKig8+UPDBBwoqKuhvKCzE7WzhsjIa7fLzqQlHbWnThmLnphlEeDg3px336UNJT3Fxll4CTk7As8+qMBgY1q2TcO2aZVvFxWSrrl9PIdj0dIYnnrB9EEwjncFAo3xdC6RkmdJ2Td/bvTslBXFOQhEeruKHH6g3BWNUDJSQICElxQELFij44AMVjo4c5eUM3t4cP/4oY/FiiuMnJTF88omCRYtkm+MaNUrFtWtk/lXl1i2KtpiOjUwAqqpdv17C00/bljHn5THNTtum6lgtBg/m5uxW65lC27Ywz+xqYsAAWmDINEDv2kUl5ZJE5tB331FBl8FAqflbtkiIjZVw5gw1J5ZlulZ//ztd95ISyr3o0YNjxAiO3r1V/P3vMkJDHWEwAGVlVJfRrRsgy9RkqF8/FcuXV+Cpp1T8+msFXnxRrdGhWlsarPBKr9d/Ut/bHDSIIzOT47XXFCxYIGPuXAXFxajsDk0P/MGD1KPg2DGG9esZ2rdXK0OMwNmzlD7bpAndBapK2XOM0c9JSUwzicSUM3DgADXfCAqqfaJJu3YcHh4ca9dSYdSpUwwFBTSiODqSfV81wWfQIA5V5YiLI/NHkgBHR45Bg+hYiioDIVVF7PBhVplTIaFHj5pX866OwYPJZPD2prUdnniCwnDt2tEN/913FXjnHcrg69GDY+RIsr3376fai8WLK/DOO47Iy6P1ODIzgVmzZEyfrmDSJJpxfPihjFmzFHOvhIgIykHZtcu26cqhQ5ZGLbt2kVm3ZYuE48ep3V7VqsTq6NyZY/16ya4GA7CUR5eVWepgmjal16orxqqKlxdHkybAJ59IcHWlZLhr1yib9LHHKIvWJHL5+YCLi20vD4OB6lxatVJx7pyEP/1JRVERxxtvOGDvXtOxUduA1q2BsDBKAZdlOnetWnH8/LNkTv2vLx6oqsyTJ8kzPmIEtW2fM0fGkCF0g+7ZQ1Vrb7+tgDHy3r7+ugyjkaaKphvJZH+aHsjHH1fh60sPhIuLfU4A59QExNOTshevX2c2jsra0KwZVUaqKlBcTGaGkxNw/boET0/Ks6+ajCVJlBClNa3dudO2T4GJjAxqXkJt52y7TdWFCRNU/PgjRU3atuVIS6PEnBdfpO0tWKDgk08oQYy6aJVhwwaKXPz5z0CnTuV44w0HXLliGs05/vMfGXv3qpg9W8WIESp++knCSy+p5uI3nY6yWFesoJlRkyaUgeruTit+FRYC//mPjIkTVZumKCYOHbI3KaxxcbFPUDIxbJiKHTssER5Tf8ohQzh++822OlKL/HxabMhoZEhOJiGfM0dB9+72nz1wgCJeAM1gf/tNgpMTXbPHHwe8vVX85S8OOHiQzF7GaKbw/POWoqu//lXGxYvAq68qKClh2LSJ2uxrdQq7Gx4ocejVi+PECQVz58pQFOr0s38/Q1ERpaz6+VlmBDt3UuhyyxYJ3bqpNZZEFxeblpujHHhrgTh6lBJw0tNJkG7evPNqREmiGzE2lkaOP/1JxezZ1Ky2tms65uaSaVN1hMjLowSsU6cYevSgeoo77bsgSVSLkpZGqb8vvEBLy+3dSw4xZ2fgb3+jGUVmJsNnn/kgMpKKuVaupIjLSy+RubBhA4mInx/H5csMb74p47HHOCIjuXntCJMTNjCQKmb37qWWd+fPM0iShN27GUaMIEeoljAAVPPSv3/153DIEGqBp9VM1tWVZnEFBbYzENOsorzcPkPThKqST+Onn2Rz4ZtOB01hAGB2Rh49ynDlCqvsI0kO0Q8+cMDOnSTEsszRpQs3J/QxRo1t/vY3GS+8oKJzZ8pEZQz1mjJtzQMlDtevA8eOOeMvf1HQpAk1ZPHzoxDh4cPkpEtJYSgrIw+wlxfQurWCL76Q8c472is/VVRQu65p0yiNdetW2xH3xg1yPm7cSHFmk0lyp3h7o3KBGhKh1q3Jy56cDM1qwapER0t2HacAap46fDiJ4cCBXLOXQV1gjB7WwEDaTng4hc1WrZLg50e9ElxcUBm1yEOTJv44cECCpyeN/l27cvj7M3z8sYLFi2UcOkQLBhcXk4f/5EkJY8YoWLWKzr1p1GOMRuzNmxnmzFGwZYuEOXMUnD3L7GoWTJjWGrWmooJyKioqaPodEECDQFUTzsTQoRQyDA7muHIFZpMsMpIiWSanclWOHWP47TeGV1+lbMymTcn/pUV6OvlQVq2S4ObGcesWXfuVK2WcPAmkp1PS2ogRCl57jRznDg4UXr98mf599JFi9rNUrUatbx4oceAc+MMfitCkCWWeDBhAttaUKTRlT0mhikHrkb95c+CNNxTMmSPjo49sBcIUCps0iWw1R0dLU5YWLWjUbNfOcgEuXND2SdQVk9/AzY1CVO+8I0OSpNv6Ms6fZwgJsc9d4JxGpIoKmj4fOSLZOSvrg5YtKe5eXg5zPsNTTwHff++KZ56hVGuA3tu8WUJMDEVe5s+vwPffU56EKU+hqIjjhx9k9OjBcf488PrrqnkBYEWhf9HREsLCqI9jVpZl0ZqqxMYytGxJ9SimB9/BgfI8TKXYCQkUFZo/X0LXrtTzwfpeIKGmisu4OAlt2tCxeHiQY7Rq/gRAofFPP5WxcGEFvL3Jb5CXxzTNOc6BJUtktGtHVcSSBFy8KMFg4DhzhmYKzz+v4KuvrP0GZFKdOkXNil566d620H+gxMG02Kw1ffuqZsdVdY1QW7QAPv9cwezZMiIiOFq0oEKqsjKaklnboYMGkbJPmkT9HqdMUXHuHIWqapszcDuoNwCNRo6OFKrKy6OkKOusSmsUBeb9qYrJYbpnDzV62blTsulDWN84Olqcdo6OwIQJRYiO9jJndvr5kY38/PMq3nhDxvnz1HhlzhwFv/wiISaGWtiZptc3bkg4c0ZCeLiKvn054uPJuTl4MDlVt2+nzltVycmhcvKDBxmefVbFqFE1tbCj87pmDQlnbCxVoHbvzs0DQHg4DTaOjrYzjKgoKqizLhmPjiazJypKRYsWNIMaO5acs1X9VqmpJHS0CLPF+R0fT9GXV19V8Omn9vteVgbExFCeh/Ugda94oMRBizZt6KEpKCBvbnW4ugIffaSYG3n26FH9SkeDBql47z0yRQDqPg3gtovO1BY3NxqNTDfg+PEqFi6UcfQotafXerC3bpXMI3NVsrMZmjTh5jUYTJ7xewVj5HA9fpyqXp96yrIQ0HvvKTh+XDLf6OHh1Hznf/+TUVBAzraMDOpvkZMj4exZSrueOZMcy2lptH3riMH587S4cXk5FU+1a8dx7ZolqcwUkvTyoo7P1h2/+vWjfJaICA6ABPmXXyRERlJPx1GjaB1NU5NagExBb28yNwICKIzq4cERFETNYuPjqXL3wAHb3gwVFXTd4uPJ8ejjQ8dz6hRljP7xjwoWLdK+pklJtHjz2LFqrRPj6psHXhwAWuZsxQoJzzxTc/NTNzdg8mRa38CULq2F0WhZX6KggBxS9e306ddPxcGD5ERzc6OcgitXyKavehx799LCvFohtZMnGUJDqRw4LIy2eb/WOejdm8rlN26khjSDBlGfCD8/mk089RTldhw5ImHsWAXHjzNcukRdkMrLOa5eBZKS0qEo03Dx4iq0bdsUzs5Uvbp5s4SzZ8lf07IlrR3SoQMQH0+5IVoYjZRiTg1TKLckMJBeU1Uyz3r1IicnzbbIUdq9Oz3oJnEAyIT94gsJzZuT+eHsTJErJyeO9HRg4ECODRsks0/q6FHqSrZhA+XP5OaSWVpRwTBmjIrlyxXNQUBRSFD8/TkmT76/K3E9FOIgSRSP37799kububhQNuKmTbRaVY8etg/SwYOU9DN8OMetW5Sf8Mor9X+Rmjen7zJNd4cPp/BeeTnw44/Uhaq4mIqoQkIs2XXWcE4jzOTJHCdOUIg1Lo7d85mDNR4eNBPKzKRIhSyTEO7ZQ3Uaw4dzDB+u4Px56mnh7k7+AFPiVF7eh8jKOogjR/6F5OSv4e9PK4b37UuhPOvGrIcOVe+kBKjXB+VNUF7JoUMMBw/SMgLbtlnCwYzR/ZOWRj6oUaMoO3L3bpphcE7Rr549OaKjqfR6yBAyNahNIJXeDx1KJu7ChTIuX6aZQ3Ex9cNgjCEiQsXChYrZt1KVxERqljNypFrjLPheUe8l23fC7Nmz+fvvv1+rz2ZmZlZbYXbmDENqKqvWs6z1eYPBsoK1olAtganN+c6dDI8+yqu9mHfLlSuwSXE+dYohPd3UcITyAzp14ppRFoA6Y3fpQkVeAQGWTlg1hfQagpquiaJQWvPNmwyxsbROR1gYR5s2tIIX5/TQDhvmDUWxL60GXNC7dyGcnOghbt2aHHddulDOwNixdRfu1FTqIB4SQkJm7XNSVWDLFuqE9cMPFCmJiZHQuzd15u7fX8XNmwxHjjAMGaKic2cgKQmYOdMBqkoFb5QPQ+Lg7k6LJy1apL0EAkB/ExNDSVpana3uhJquSVXmzp17z0q27xtdu1IIb8UKsiFvV4zWtStH165cM7wVHc0qFxRpuP1t04ZamnXuTCYDee4ltGzJ0bQpr7aXJEBT5tJS8jUcPEgCs2GD/QKz9xuqMaHRe8wYqgs5fpxWA6+oYObzvnjxBfzvf2/j6NGNKC+nhX6dncfB0/MTpKWRiDg40MP70Ucyrl6lhsGxsTT9d3YGXFzI3GjZktt1FremZUvgH/9QsXw5dbsqK6PIRng4+aFGj6ZFaMrKgFdeccALL1CtypQpZO75+nIsXy7hxg0ZiYlU39OyJTlSS0up/iMwkNKie/fm+NOfFM0FgbKyqIOXnx/HpEn3bz2Q6nioxAGgyMTUqeR0OnaMkkyqrh9p/b/ppgsIoBwGg4EhL48e2PpabKYmxo0jc+LZZ+nmGDeOVjdq1ozW19AaSUzFPM88Q36RESNUsxPufi68ejsYo0a66enAsWOSzb76+DSHm5sXKipK4ejogoqKUjRt6glv72bmfIWyMnIOl5dT6XdpKdXbxMYysxi4uVGqvLc3ZXeGhMCm5wJj5AMJCOAYOFBFYiLNPrKyKKJAXa4os7N/f46YGGDFCoZevRi2bKGQ6s2blvU2ExLoeqSn01KJzs7kcH3xRUovP3xYshEGzqmOJCWFzBtTkVtj5KETB4AuXN++NFrVhooK8phnZVHLLR+fezctlyTykK9aRfkaTk5ko1++TJ2oY2NtHYylpdQAZPx4am8fEkKj5O3ShxsTzZpBc4bz448ZmDHjZUyfPh3Lli1Deno6Vq6swE8/0SJFeXkkEDk59ICmptLDzhjljTRtSslQOTlAaiqzKcF2dKTMRAcH+rwsU8fmggKG776jDuUZGVTi7+dHnzt1in7PyqJ+CU2aUD/HgoI0FBVNg7v7Kri5NYOzM3UXHzKEmvK2aEHZlhs3UnFYRQWtP5GRQf6sHj3UakPWjYmHUhzqioMDzThatLg/F8zPjwRi5Uoyh9q0oXwOWabKvrVrybFHKd40MzpyhMHb29LGPjXV1rv+ILJ69Wrzz19+aSk3fv558jMkJTHs2MFgMNAIbDRS+7+SEhKKjAzq40hLw6HSqcchSfTwu7tzcwl3fj7DtWskNqaGvwEBJC63btF5zsuzzFgKC6nbtbs7UFLyIRTlIFT1Xxg16is89xylv3NOUYn9+ylKYVoty8EB6NlTrVP5fGNAiEMjwcuLoiixsbS6tYcH9SxwdWXmdOnychrNNm2SEBpqWehk716G/v0fjFnD3dC+PfV9UFUKFV64QO3yi4uBLl0oXVpV6VxmZNDDXFpKvgNVBVSVmU1JR0dL7YSHB1VLFhWRA1VVLZ9xdjaZKcCFCx4oKbE4TQsLF2Pp0sVYvtwFP/+cD6ORRLpVK47585V6rZC8HwhxaEQwhsoFeakb8pkzNALt2kXZn46OQOfOqk1VaGEhTY1NNRC/B0xmo6kCNTGRok75+RQ5yMtjCAoiMSkoIFMMoAdeVel3RbEsScg5mRyZmUCrVtSC3zT7AEhEevVS0azZObz77hwcO7YRxcXFcHJyRXj4WDz//HxwTl3NTan4DwNCHBopHh5Anz4cffooOHGCISuL2a1wXV5OXYeee+7hnzVUB2NUgGRdhHTjBq17WlxM0ZDiYjJBcnOpDNzUH6Eq7dpR0VTLlhwzZqg2EY+sLGD79hZo394TsbGlcHFxQVlZKTp39sQzzwQAePiugRCHB4CwMGpy8+uvFMUIDCSH18WLDFOnqna5/L93LNWkdzabuniROoi1acMrO2Mx5ObSWh6bN2fi5ZdtnaYPK0IcHhACAqgtW1YW1VJ4eXFMnvz7MSXuJR06UO/Rq1fJWdm9Ozf3CqnOafowIsThAaNJE9x1TwlB7aAlBH+/57qRpl8IBIL7jRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgSb2veKXT6WZV/tgbwHG9Xv9JfX+HQCBoeOpVHHQ63WK9Xv+K1e8ndDodhEAIBA8e9WZW6HQ6HwDGKi8vBjCnvr5DIBDcO+pz5uAHYFbl7MFg9bpPbf44MzOzVl9iNFbVH8H9RlyTxkd9XJN6Ewe9Xm/Q6XRhVYRhKIDo2vx9QEBArb+rLp8V3BvENWl83O01qVefg16vP2n6udLMiAIQVp/fIRAI7g23FYfKh7xa9Hp9dfOXXwFEVplJCASCB4QaxUGn000AmQY1fcao1+vfrvLafADzrWcSAoHgwaJGcdDr9WsArKnLBisFZZder4+u/L2XEAmB4MGjvvMcokBRi+hKc8QPwGQAQhwEggeMehOHSjHYVfnrYqu36jTzEAgEjYP6DGUaAbD62p5AILi/iMIrgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgiRAHgUCgSb12n74b5s6de793QSAQWME45/d7HwQCQSNEmBUCgUATIQ4CgUATIQ4CgUCTRuOQFAgEtadyTdreVRexruazswAYQMtTQq/XL6nNdzQ6cajrgdzpgQtqT13OceVNGwxaBjEHwAwAa/R6veEe7OpDT+V6tL0ADAVdk9t9fj6A45WLYkOn083X6XQTTL/XRKMyKyoPxKDX69dU3oAhlTdbvXxeUHfu4Bz7AZgP4DKA5Mq/FcJQT+j1+mi9Xv8Jar849YwqQrALwCu1+cNGJQ6o+4Hc8YELas2dnGNfACF6vd63NiOUoGHQ6XS9NF7OARBVm79vNGZFXQ/kbg9ccHvu9BxXLqpsbJCdEtQFP9D1ssYIADqdzqfyOlVLY5o51Hgg9fB5Qd25o3Os0+lm6HS6CZX/z2rIHRTUiA8q/URWmK5n1dftaEziUNcDuasDF9SKOznH0QB+qeKjmNFQOyioEa2Zgem6VRV9OxqNWYG6H8hdHbigVtT5HGs4H3eBHJQiinTvyQEJvDU+gNn0q5HGNHOo64Hc1YELakWdzrFOp/PR6XS8islhBIU2BfcYvV5/EvYC7wea3d2WRiMOdT2Quz1wwe25w3P8SRXhCEYt4vGC+kGn0wVXCTUvqfL7UACLa7OtRiMOldR4IPV54IJaU+trUikK2VX+fiKA22bxCWqHTqfrVenknQBgkk6nm1UlqhQFq1BzZQZlcKWDeBaAy7UNLze6km2rbLxgAEbrbLxKx9ZEvV4/tDafF9QPdbkmlSbFDNCMIwRW2XmCB4tGJw4CgaBx0NjMCoFA0EgQ4iAQCDQR4iAQCDQR4iAQCDQR4iAQCDQR4iAQCDQR4iAQCDQR4iAQCDQR4iAQCDT5/9DQZtbijvGlAAAAAElFTkSuQmCC\n",
+      "text/plain": [
+       "<Figure size 288x216 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "with torch.no_grad():\n",
+    "    # Initialize plot\n",
+    "    f, ax = plt.subplots(1, 1, figsize=(4, 3))\n",
+    "    \n",
+    "    # Plot training data as black stars\n",
+    "    ax.plot(train_x.numpy(), train_y.numpy(), 'k*', zorder=10)\n",
+    "    \n",
+    "    for i in range(25):\n",
+    "        # Plot predictive means as blue line\n",
+    "        ax.plot(test_x.numpy(), output.mean[i].detach().numpy(), 'b', linewidth=0.3)\n",
+    "        \n",
+    "    # Shade between the lower and upper confidence bounds\n",
+    "    # ax.fill_between(test_x.numpy(), lower.numpy(), upper.numpy(), alpha=0.5)\n",
+    "    ax.set_ylim([-3, 3])\n",
+    "    ax.legend(['Observed Data', 'Sampled Means'])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "anaconda-cloud": {},
+  "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.1"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 1
+}
diff --git a/examples/01_Simple_GP_Regression/index.rst b/examples/01_Simple_GP_Regression/index.rst
index f9f40201..734e2146 100644
--- a/examples/01_Simple_GP_Regression/index.rst
+++ b/examples/01_Simple_GP_Regression/index.rst
@@ -5,6 +5,7 @@
    :hidden:
 
    Simple_GP_Regression.ipynb
+   Simple_GP_Regression_Fully_Bayesian.ipynb
    Simple_GP_Regression_CUDA.ipynb
    Simple_MultiGPU_GP_Regression.ipynb
    Simple_GP_Regression_With_LOVE_Fast_Variances_CUDA.ipynb
diff --git a/gpytorch/distributions/distribution.py b/gpytorch/distributions/distribution.py
index 7aed809d..a6d675be 100644
--- a/gpytorch/distributions/distribution.py
+++ b/gpytorch/distributions/distribution.py
@@ -3,7 +3,7 @@
 from torch.distributions import Distribution as TDistribution
 
 
-class Distribution(TDistribution):
+class _DistributionBase(TDistribution):
     """
     The base class of Distributions. (Same as torch.distribution.Distribution
     or pyro.distribution.Distribution).
@@ -20,3 +20,17 @@ class Distribution(TDistribution):
 
     def __mul__(self, other):
         raise NotImplementedError()
+
+
+try:
+    # If pyro is installed, add the TorchDistributionMixin
+    from pyro.distributions.torch_distribution import TorchDistributionMixin
+
+    class Distribution(_DistributionBase, TorchDistributionMixin):
+        pass
+
+
+except ImportError:
+
+    class Distribution(_DistributionBase):
+        pass
diff --git a/gpytorch/distributions/multivariate_normal.py b/gpytorch/distributions/multivariate_normal.py
index f24e4b00..c7df69e9 100644
--- a/gpytorch/distributions/multivariate_normal.py
+++ b/gpytorch/distributions/multivariate_normal.py
@@ -13,7 +13,7 @@ from .distribution import Distribution
 from ..utils.broadcasting import _mul_broadcast_shape
 
 
-class _MultivariateNormalBase(TMultivariateNormal, Distribution):
+class MultivariateNormal(TMultivariateNormal, Distribution):
     """
     Constructs a multivariate normal random variable, based on mean and covariance.
     Can be multivariate, or a batch of multivariate normals
@@ -185,7 +185,7 @@ class _MultivariateNormalBase(TMultivariateNormal, Distribution):
             return super().variance
 
     def __add__(self, other):
-        if isinstance(other, _MultivariateNormalBase):
+        if isinstance(other, MultivariateNormal):
             return self.__class__(
                 mean=self.mean + other.mean,
                 covariance_matrix=(self.lazy_covariance_matrix + other.lazy_covariance_matrix),
@@ -211,52 +211,6 @@ class _MultivariateNormalBase(TMultivariateNormal, Distribution):
         return self.__mul__(1.0 / other)
 
 
-try:
-    # If pyro is installed, add the TorchDistributionMixin
-    from pyro.distributions.torch_distribution import TorchDistributionMixin
-
-    class MultivariateNormal(_MultivariateNormalBase, TorchDistributionMixin):
-        """
-        Constructs a multivariate normal random variable, based on mean and covariance.
-        Can be multivariate, or a batch of multivariate normals
-
-        Passing a vector mean corresponds to a multivariate normal.
-        Passing a matrix mean corresponds to a batch of multivariate normals.
-
-        :param torch.tensor mean: Vector n or matrix b x n mean of mvn distribution.
-        :param ~gpytorch.lazy.LazyTensor covar: Matrix n x n or batch matrix b x n x n covariance of
-            mvn distribution.
-        """
-        def confidence_region(self):
-            """
-            Returns 2 standard deviations above and below the mean.
-
-            :rtype: (torch.Tensor, torch.Tensor)
-            :return: pair of tensors of size (b x d) or (d), where b is the
-              batch size and d is the dimensionality of the random
-              variable. The first (second) Tensor is the lower (upper) end of
-              the confidence region.
-            """
-            return super().confidence_region()
-
-        def get_base_samples(self, sample_shape=torch.Size()):
-            """Get i.i.d. standard Normal samples (to be used with rsample(base_samples=base_samples))"""
-            return super().get_base_samples(sample_shape=sample_shape)
-
-        @lazy_property
-        def lazy_covariance_matrix(self):
-            """
-            The covariance_matrix, represented as a LazyTensor
-            """
-            return super().lazy_covariance_matrix
-
-
-except ImportError:
-
-    class MultivariateNormal(_MultivariateNormalBase):
-        pass
-
-
 @register_kl(MultivariateNormal, MultivariateNormal)
 def kl_mvn_mvn(p_dist, q_dist):
     output_shape = _mul_broadcast_shape(p_dist.batch_shape, q_dist.batch_shape)
diff --git a/gpytorch/kernels/kernel.py b/gpytorch/kernels/kernel.py
index 0b1e1be5..87ba18a3 100644
--- a/gpytorch/kernels/kernel.py
+++ b/gpytorch/kernels/kernel.py
@@ -241,6 +241,11 @@ class Kernel(Module, _ClassWithDeprecatedBatchSize):
 
         self.initialize(raw_lengthscale=self.raw_lengthscale_constraint.inverse_transform(value))
 
+    def local_load_samples(self, samples_dict, memo, prefix):
+        num_samples = next(iter(samples_dict.values())).size(0)
+        self.batch_shape = torch.Size([num_samples]) + self.batch_shape
+        super().local_load_samples(samples_dict, memo, prefix)
+
     def covar_dist(
         self,
         x1,
diff --git a/gpytorch/mlls/exact_marginal_log_likelihood.py b/gpytorch/mlls/exact_marginal_log_likelihood.py
index bce6eed7..877a8ded 100644
--- a/gpytorch/mlls/exact_marginal_log_likelihood.py
+++ b/gpytorch/mlls/exact_marginal_log_likelihood.py
@@ -60,3 +60,9 @@ class ExactMarginalLogLikelihood(MarginalLogLikelihood):
         # Scale by the amount of data we have
         num_data = target.size(-1)
         return res.div_(num_data)
+
+    def pyro_factor(self, output, target, *params):
+        import pyro
+        mll = self(output, target, *params)
+        pyro.factor("gp_mll", mll)
+        return mll
diff --git a/gpytorch/mlls/marginal_log_likelihood.py b/gpytorch/mlls/marginal_log_likelihood.py
index 37ea1ab6..0fd3bc43 100644
--- a/gpytorch/mlls/marginal_log_likelihood.py
+++ b/gpytorch/mlls/marginal_log_likelihood.py
@@ -50,3 +50,9 @@ class MarginalLogLikelihood(Module):
                 Additional arguments to pass to the likelihood's :attr:`forward` function.
         """
         raise NotImplementedError
+
+    def pyro_factor(self, output, target):
+        """
+        As forward, but register the MLL with pyro using the pyro.factor primitive.
+        """
+        raise NotImplementedError
diff --git a/gpytorch/models/exact_gp.py b/gpytorch/models/exact_gp.py
index 261d3a89..a6d573b7 100644
--- a/gpytorch/models/exact_gp.py
+++ b/gpytorch/models/exact_gp.py
@@ -45,6 +45,15 @@ class ExactGP(GP):
             self.train_targets = fn(self.train_targets)
         return super(ExactGP, self)._apply(fn)
 
+    def local_load_samples(self, samples_dict, memo, prefix):
+        # Pyro always puts the samples in the first batch dimension
+        num_samples = next(iter(samples_dict.values())).size(0)
+        self.train_inputs = tuple(
+            tri.unsqueeze(0).expand(num_samples, *tri.shape) for tri in self.train_inputs
+        )
+        self.train_targets = self.train_targets.unsqueeze(0).expand(num_samples, *self.train_targets.shape)
+        super().local_load_samples(samples_dict, memo, prefix)
+
     def set_train_data(self, inputs=None, targets=None, strict=True):
         """
         Set training data (does not re-fit model hyper-parameters).
diff --git a/gpytorch/module.py b/gpytorch/module.py
index 1c6fd774..14ce34dc 100644
--- a/gpytorch/module.py
+++ b/gpytorch/module.py
@@ -18,6 +18,8 @@ class Module(nn.Module):
         self._priors = OrderedDict()
         self._constraints = OrderedDict()
 
+        self._strict_init = True
+
     def __call__(self, *inputs, **kwargs):
         outputs = self.forward(*inputs, **kwargs)
         if isinstance(outputs, list):
@@ -34,6 +36,9 @@ class Module(nn.Module):
                 "Invalid parameter name {}. {} has no module {}".format(parameter_name, type(self).__name__, module)
             )
 
+    def _strict(self, value):
+        _set_strict(self, value)
+
     def added_loss_terms(self):
         for _, strategy in self.named_added_loss_terms():
             yield strategy
@@ -88,7 +93,10 @@ class Module(nn.Module):
                 try:
                     self.__getattr__(name).data.copy_(val.expand_as(self.__getattr__(name)))
                 except RuntimeError:
-                    self.__getattr__(name).data.copy_(val.view_as(self.__getattr__(name)))
+                    if not self._strict_init:
+                        self.__getattr__(name).data = val
+                    else:
+                        self.__getattr__(name).data.copy_(val.view_as(self.__getattr__(name)))
 
             elif isinstance(val, float):
                 constraint = self.constraint_for_parameter_name(name)
@@ -198,16 +206,16 @@ class Module(nn.Module):
 
         """
         if isinstance(param_or_closure, str):
-            if param_or_closure not in self._parameters:
+            if param_or_closure not in self._parameters and not hasattr(self, param_or_closure):
                 raise AttributeError(
                     "Unknown parameter {name} for {module}".format(
                         name=param_or_closure, module=self.__class__.__name__
                     )
-                    + "Make sure the parameter is registered before registering a prior."
+                    + " Make sure the parameter is registered before registering a prior."
                 )
 
             def closure():
-                return self._parameters[param_or_closure]
+                return getattr(self, param_or_closure)
 
             if setting_closure is not None:
                 raise RuntimeError("Must specify a closure instead of a parameter name when providing setting_closure")
@@ -270,6 +278,46 @@ class Module(nn.Module):
             raise RuntimeError("Must provide inverse transform to be able to sample from prior.")
         setting_closure(prior.sample())
 
+    def pyro_sample_from_prior(self):
+        """
+        For each parameter in this Module and submodule that have defined priors, sample a value for that parameter
+        from its corresponding prior with a pyro.sample primitive and load the resulting value in to the parameter.
+
+        This method can be used in a Pyro model to conveniently define pyro sample sites for all
+        parameters of the model that have GPyTorch priors registered to them.
+        """
+        return _pyro_sample_from_prior(module=self, memo=None, prefix="")
+
+    def local_load_samples(self, samples_dict, memo, prefix):
+        """
+        Defines local behavior of this Module when loading parameters from a samples_dict generated by a Pyro
+        sampling mechanism.
+
+        The default behavior here should almost always be called from any overriding class. However, a class may
+        want to add additional functionality, such as reshaping things to account for the fact that parameters will
+        acquire an extra batch dimension corresponding to the number of samples drawn.
+        """
+        self._strict(False)
+        for name, (prior, closure, setting_closure) in self._priors.items():
+            if prior is not None and prior not in memo:
+                memo.add(prior)
+                setting_closure(samples_dict[prefix + ("." if prefix else "") + name])
+        self._strict(True)
+
+    def pyro_load_from_samples(self, samples_dict):
+        """
+        Convert this Module in to a batch Module by loading parameters from the given `samples_dict`. `samples_dict`
+        is typically produced by a Pyro sampling mechanism.
+
+        Note that the keys of the samples_dict should correspond to prior names (covar_module.outputscale_prior) rather
+        than parameter names (covar_module.raw_outputscale), because we will use the setting_closure associated with
+        the prior to properly set the unconstrained parameter.
+
+        Args:
+            :attr:`samples_dict` (dict): Dictionary mapping *prior names* to sample values.
+        """
+        return _pyro_load_from_samples(module=self, samples_dict=samples_dict, memo=None, prefix="")
+
     def update_added_loss_term(self, name, added_loss_term):
         from .mlls import AddedLossTerm
 
@@ -314,6 +362,53 @@ def _validate_module_outputs(outputs):
         )
 
 
+def _set_strict(module, value, memo=None):
+    if memo is None:
+        memo = set()
+
+    if hasattr(module, "_strict_init"):
+        module._strict_init = value
+
+    for mname, module_ in module.named_children():
+        _set_strict(module_, value)
+
+
+def _pyro_sample_from_prior(module, memo=None, prefix=""):
+    try:
+        import pyro
+    except ImportError:
+        raise RuntimeError("Cannot call pyro_sample_from_prior without pyro installed!")
+    if memo is None:
+        memo = set()
+    if hasattr(module, "_priors"):
+        for prior_name, (prior, closure, setting_closure) in module._priors.items():
+            if prior is not None and prior not in memo:
+                if setting_closure is None:
+                    raise RuntimeError(
+                        "Cannot use Pyro for sampling without a setting_closure for each prior,"
+                        f" but the following prior had none: {prior_name}, {prior}."
+                    )
+                memo.add(prior)
+                prior = prior.expand(closure().shape)
+                value = pyro.sample(prefix + ("." if prefix else "") + prior_name, prior)
+                setting_closure(value)
+
+    for mname, module_ in module.named_children():
+        submodule_prefix = prefix + ("." if prefix else "") + mname
+        _pyro_sample_from_prior(module=module_, memo=memo, prefix=submodule_prefix)
+
+
+def _pyro_load_from_samples(module, samples_dict, memo=None, prefix=""):
+    if memo is None:
+        memo = set()
+    if hasattr(module, "_priors"):
+        module.local_load_samples(samples_dict, memo, prefix)
+
+    for mname, module_ in module.named_children():
+        submodule_prefix = prefix + ("." if prefix else "") + mname
+        _pyro_load_from_samples(module_, samples_dict, memo=memo, prefix=submodule_prefix)
+
+
 def _extract_named_added_loss_terms(module, memo=None, prefix=""):
     if memo is None:
         memo = set()
diff --git a/gpytorch/priors/__init__.py b/gpytorch/priors/__init__.py
index 0ece1e34..7e026d85 100644
--- a/gpytorch/priors/__init__.py
+++ b/gpytorch/priors/__init__.py
@@ -4,7 +4,7 @@ from .prior import Prior
 from .horseshoe_prior import HorseshoePrior
 from .lkj_prior import LKJCholeskyFactorPrior, LKJCovariancePrior, LKJPrior
 from .smoothed_box_prior import SmoothedBoxPrior
-from .torch_priors import GammaPrior, MultivariateNormalPrior, NormalPrior
+from .torch_priors import GammaPrior, MultivariateNormalPrior, NormalPrior, LogNormalPrior, UniformPrior
 
 
 # from .wishart_prior import InverseWishartPrior, WishartPrior
@@ -17,9 +17,11 @@ __all__ = [
     "LKJPrior",
     "LKJCholeskyFactorPrior",
     "LKJCovariancePrior",
+    "LogNormalPrior",
     "MultivariateNormalPrior",
     "NormalPrior",
     "SmoothedBoxPrior",
+    "UniformPrior",
     # "InverseWishartPrior",
     # "WishartPrior",
 ]
diff --git a/gpytorch/priors/horseshoe_prior.py b/gpytorch/priors/horseshoe_prior.py
index dbc263a7..bc8c9587 100644
--- a/gpytorch/priors/horseshoe_prior.py
+++ b/gpytorch/priors/horseshoe_prior.py
@@ -5,7 +5,8 @@ from numbers import Number
 
 import torch
 from gpytorch.priors.prior import Prior
-from torch.distributions import constraints
+from gpytorch.distributions import Distribution
+from torch.distributions import constraints, HalfCauchy, Normal
 from torch.nn import Module as TModule
 
 
@@ -47,3 +48,16 @@ class HorseshoePrior(Prior):
         lb = self.K / 2 * torch.log(1 + 4 * A)
         ub = self.K * torch.log(1 + 2 * A)
         return torch.log((lb + ub) / 2)
+
+    def rsample(self, sample_shape=torch.Size([])):
+        local_shrinkage = HalfCauchy(1).rsample(self.scale.shape)
+        param_sample = Normal(0, local_shrinkage * self.scale).rsample(sample_shape)
+        return param_sample
+
+    def expand(self, expand_shape, _instance=None):
+        new = self._get_checked_instance(HorseshoePrior)
+        batch_shape = torch.Size(expand_shape)
+        new.scale = self.scale.expand(batch_shape)
+        super(Distribution, new).__init__(batch_shape)
+        new._validate_args = self._validate_args
+        return new
diff --git a/gpytorch/priors/prior.py b/gpytorch/priors/prior.py
index e70fce9d..2286f32f 100644
--- a/gpytorch/priors/prior.py
+++ b/gpytorch/priors/prior.py
@@ -2,7 +2,7 @@
 
 from abc import ABC
 
-from torch.distributions import Distribution
+from ..distributions import Distribution
 from torch.nn import Module
 
 
diff --git a/gpytorch/priors/torch_priors.py b/gpytorch/priors/torch_priors.py
index ea105f70..8e35fec5 100644
--- a/gpytorch/priors/torch_priors.py
+++ b/gpytorch/priors/torch_priors.py
@@ -1,6 +1,6 @@
 #!/usr/bin/env python3
 
-from torch.distributions import Gamma, MultivariateNormal, Normal
+from torch.distributions import Gamma, MultivariateNormal, Normal, LogNormal, Uniform
 
 from .prior import Prior
 from .utils import _bufferize_attributes, _del_attributes
@@ -10,7 +10,8 @@ MVN_LAZY_PROPERTIES = ("covariance_matrix", "scale_tril", "precision_matrix")
 
 
 class NormalPrior(Prior, Normal):
-    """Normal (Gaussian) Prior
+    """
+    Normal (Gaussian) Prior
 
     pdf(x) = (2 * pi * sigma^2)^-0.5 * exp(-(x - mu)^2 / (2 * sigma^2))
 
@@ -23,6 +24,35 @@ class NormalPrior(Prior, Normal):
         _bufferize_attributes(self, ("loc", "scale"))
         self._transform = transform
 
+    def expand(self, batch_shape):
+        return Normal.expand(self, batch_shape, _instance=self)
+
+
+class LogNormalPrior(Prior, LogNormal):
+    """
+    Log Normal prior.
+    """
+    def __init__(self, loc, scale, validate_args=None, transform=None):
+        TModule.__init__(self)
+        LogNormal.__init__(self, loc=loc, scale=scale, validate_args=validate_args)
+        self._transform = transform
+
+    def expand(self, batch_shape):
+        return LogNormal.expand(self, batch_shape, _instance=self)
+
+
+class UniformPrior(Prior, Uniform):
+    """
+    Uniform prior.
+    """
+    def __init__(self, a, b, validate_args=None, transform=None):
+        TModule.__init__(self)
+        Uniform.__init__(self, a, b, validate_args=validate_args)
+        self._transform = transform
+
+    def expand(self, batch_shape):
+        return Uniform.expand(self, batch_shape, _instance=self)
+
 
 class GammaPrior(Prior, Gamma):
     """Gamma Prior parameterized by concentration and rate
@@ -38,6 +68,9 @@ class GammaPrior(Prior, Gamma):
         _bufferize_attributes(self, ("concentration", "rate"))
         self._transform = transform
 
+    def expand(self, batch_shape):
+        return Gamma.expand(self, batch_shape, _instance=self)
+
 
 class MultivariateNormalPrior(Prior, MultivariateNormal):
     """Multivariate Normal prior
@@ -73,3 +106,6 @@ class MultivariateNormalPrior(Prior, MultivariateNormal):
         module = self._apply(lambda t: t.cpu())
         _del_attributes(module, MVN_LAZY_PROPERTIES)
         return module
+
+    def expand(self, batch_shape):
+        return MultivariateNormal.expand(self, batch_shape, _instance=self)


Test patch:
diff --git a/test/examples/test_simple_gp_regression.py b/test/examples/test_simple_gp_regression.py
index 7e2d96b9..b0d71fd2 100644
--- a/test/examples/test_simple_gp_regression.py
+++ b/test/examples/test_simple_gp_regression.py
@@ -9,7 +9,7 @@ from gpytorch.distributions import MultivariateNormal
 from gpytorch.kernels import RBFKernel, ScaleKernel
 from gpytorch.likelihoods import GaussianLikelihood
 from gpytorch.means import ConstantMean
-from gpytorch.priors import SmoothedBoxPrior
+from gpytorch.priors import SmoothedBoxPrior, LogNormalPrior, NormalPrior, UniformPrior
 from gpytorch.constraints import Positive
 from gpytorch.test.base_test_case import BaseTestCase
 from gpytorch.test.utils import least_used_cuda_device
@@ -19,8 +19,8 @@ from torch import optim
 class ExactGPModel(gpytorch.models.ExactGP):
     def __init__(self, train_inputs, train_targets, likelihood):
         super(ExactGPModel, self).__init__(train_inputs, train_targets, likelihood)
-        self.mean_module = ConstantMean(prior=SmoothedBoxPrior(-1, 1))
-        self.covar_module = ScaleKernel(RBFKernel(lengthscale_prior=SmoothedBoxPrior(exp(-3), exp(3), sigma=0.1)))
+        self.mean_module = ConstantMean()
+        self.covar_module = ScaleKernel(RBFKernel())
 
     def forward(self, x):
         mean_x = self.mean_module(x)
@@ -453,6 +453,47 @@ class TestSimpleGPRegression(BaseTestCase, unittest.TestCase):
 
             self.assertLess(torch.max(var_diff / noise), 0.05)
 
+    def test_pyro_sampling(self):
+        try:
+            import pyro
+            from pyro.infer.mcmc import NUTS, MCMC
+        except:
+            return
+        train_x, test_x, train_y, test_y = self._get_data(cuda=False)
+        likelihood = GaussianLikelihood(noise_constraint=gpytorch.constraints.Positive())
+        gp_model = ExactGPModel(train_x, train_y, likelihood)
+
+        # Register normal GPyTorch priors
+        gp_model.mean_module.register_prior("mean_prior", UniformPrior(-1, 1), "constant")
+        gp_model.covar_module.base_kernel.register_prior("lengthscale_prior", UniformPrior(0.01, 0.2), "lengthscale")
+        gp_model.covar_module.register_prior("outputscale_prior", UniformPrior(1, 2), "outputscale")
+        likelihood.register_prior("noise_prior", LogNormalPrior(-1.5, 0.1), "noise")
+
+        mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, gp_model)
+
+        def pyro_model(x, y):
+            gp_model.pyro_sample_from_prior()
+            output = gp_model(x)
+            loss = mll.pyro_factor(output, y)
+            return y
+
+        nuts_kernel = NUTS(pyro_model, adapt_step_size=True)
+        mcmc_run = MCMC(nuts_kernel, num_samples=3, warmup_steps=20)
+        mcmc_run.run(train_x, train_y)
+
+        gp_model.pyro_load_from_samples(mcmc_run.get_samples())
+
+        gp_model.eval()
+        expanded_test_x = test_x.unsqueeze(-1).repeat(3, 1, 1)
+        output = gp_model(expanded_test_x)
+
+        self.assertEqual(output.mean.size(0), 3)
+
+        # All 3 samples should do reasonably well on a noiseless dataset.
+        self.assertLess(torch.norm(output.mean[0] - test_y) / test_y.norm(), 0.2)
+        self.assertLess(torch.norm(output.mean[1] - test_y) / test_y.norm(), 0.2)
+        self.assertLess(torch.norm(output.mean[2] - test_y) / test_y.norm(), 0.2)
+
     def test_posterior_latent_gp_and_likelihood_fast_pred_var_cuda(self):
         if torch.cuda.is_available():
             with least_used_cuda_device():