00_getting_started_flexible.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Getting started with `sbi`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note, you can find the original version of this notebook at [https://github.com/sbi-dev/sbi/blob/main/tutorials/00_getting_started_flexible.ipynb](https://github.com/sbi-dev/sbi/blob/main/tutorials/00_getting_started_flexible.ipynb) in the `sbi` repository."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`sbi` provides a simple interface to run state-of-the-art algorithms for simulation-based inference.\n",
    "\n",
    "**The overall goal of simulation-based inference is to algorithmically identify model parameters which are consistent with data.**\n",
    "\n",
    "In this tutorial we demonstrate how to get started with the `sbi` toolbox and how to perform parameter inference on a simple model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "\n",
    "from sbi.analysis import pairplot\n",
    "from sbi.utils import BoxUniform\n",
    "from sbi.inference import SNPE, simulate_for_sbi\n",
    "from sbi.utils.user_input_checks import (\n",
    "    check_sbi_inputs,\n",
    "    process_prior,\n",
    "    process_simulator,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameter inference in a linear Gaussian example\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each of the implemented inference methods takes three inputs:  \n",
    "1. observational data (or summary statistics thereof) - _the observations_  \n",
    "2. a candidate (mechanistic) model - _the simulator_  \n",
    "3. prior knowledge or constraints on model parameters - _the prior_\n",
    "\n",
    "If you are new to simulation-based inference, please first read the information on the [homepage of the website](https://sbi-dev.github.io/sbi/) to familiarise with the motivation and relevant terms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this illustrative example we consider a model _simulator_ that takes in 3 parameters ($\\theta$). For simplicity, the _simulator_ outputs simulations of the same dimensionality and adds 1.0 and some Gaussian noise to the parameter set. \n",
    "\n",
    "> Note: This is where you instead would use your specific  _simulator_ with its parameters.\n",
    "\n",
    "For the 3-dimensional parameter space we consider a uniform _prior_ between [-2,2].\n",
    "\n",
    "> Note: This is where you would incorporate prior knowlegde about the parameters you want to infer, e.g., ranges known from literature. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_dim = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "from Matlab2Py.offline_matlab import OfflineEngine\n",
    "import numpy as np\n",
    "\n",
    "def init_simulator(SIMULATION_PATH, model_name,):\n",
    "    #Here you setup your simulation configuration & initalzie necessary variables \n",
    "\n",
    "\n",
    "    MODEL_PATH = f\"{SIMULATION_PATH}{model_name}.slx\"\n",
    "    \n",
    "    eng = OfflineEngine(model_path = MODEL_PATH, sim_path = SIMULATION_PATH, model_name = model_name)\n",
    "    eng.load_engine()\n",
    "    #Run the simple_sim/init.m to initalize theta\n",
    "    eng.run_engine_script('init')\n",
    "    eng.set_simulation_mode(s_mode='normal')\n",
    "\n",
    "    return eng\n",
    "\n",
    "SIMULATION_PATH = 'simple_sim/'\n",
    "model_name = 'test'\n",
    "eng = init_simulator(SIMULATION_PATH=SIMULATION_PATH, model_name=model_name)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def run_simulation(eng, theta : list[3]):\n",
    "    if type(theta) == torch.Tensor:\n",
    "        theta = theta.tolist()\n",
    "       \n",
    "    seeds = np.random.randint(0, np.iinfo(np.int32).max, size = num_dim)\n",
    "\n",
    "    o = eng.run_full_sim(blocks_inputs={'gau' : (seeds, 'seed'), 'theta' : (theta, 'Value')})\n",
    "    return eng.get_simoutput(o, 'simout')\n",
    "\n",
    "\n",
    "def simulator(eng):\n",
    "    def _simulator(inputs):\n",
    "        output = (torch.Tensor(run_simulation(eng, theta)).squeeze() for theta in inputs)\n",
    "        return torch.stack((*output,) )\n",
    "    return _simulator\n",
    "\n",
    "\n",
    "simul = simulator(eng)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have to ensure that your _simulator_ and _prior_ adhere to the requirements of `sbi` such as returning `torch.Tensor`s in a standardised shape. \n",
    "\n",
    "You can do so with the `process_simulator()` and `process_prior()` functions, which prepare them appropriately. Finally, you can call `check_sbi_input()` to make sure they are consistent which each other."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check prior, return PyTorch prior.\n",
    "prior = BoxUniform(low=-2 * torch.ones(num_dim), high=2 * torch.ones(num_dim))\n",
    "\n",
    "prior, num_parameters, prior_returns_numpy = process_prior(prior)\n",
    "\n",
    "# Check simulator, returns PyTorch simulator able to simulate batches.\n",
    "simul = process_simulator(simul, prior, prior_returns_numpy)\n",
    "\n",
    "# Consistency check after making ready for sbi.\n",
    "check_sbi_inputs(simulator, prior)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we instantiate the inference object. In this example, we will use neural perform posterior estimation (NPE):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note: In the `sbi` toolbox, NPE is run by using the `SNPE` (Sequential NPE) class for only one iteration of simulation and training. \n",
    "\n",
    "> Note: This is where you could specify an alternative inference object such as (S)NRE for ratio estimation or (S)NLE for likelihood estimation. Here, you can see [all implemented methods.](16_implemented_methods.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inference = SNPE(prior=prior)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we need simulations, or more specifically, pairs of parameters $\\theta$ which we sample from the _prior_ and corresponding simulations $x = \\mathrm{simulator} (\\theta)$. The `sbi` helper function called `simulate_for_sbi` allows to parallelize your code with `joblib`.\n",
    "\n",
    " > Note: If you already have your own parameters, simulation pairs which were generated elsewhere (e.g., on a compute cluster), you would add them here.  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "28cefebf703a4b2c8abbaeb0c5c32034",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Running 2000 simulations.:   0%|          | 0/2000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "theta.shape torch.Size([2000, 3])\n",
      "x.shape torch.Size([2000, 3])\n"
     ]
    }
   ],
   "source": [
    "theta, x = simulate_for_sbi(simulator, proposal=prior, num_simulations=2000)\n",
    "print(\"theta.shape\", theta.shape)\n",
    "print(\"x.shape\", x.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then pass the simulated data to the inference object. Both `theta` and `x` should be a `torch.Tensor` of type `float32`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inference = inference.append_simulations(theta, x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we train the _neural density estimator_ to learn the association between the simulated data (or data features) and the underlying parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Neural network successfully converged after 81 epochs."
     ]
    }
   ],
   "source": [
    "density_estimator = inference.train()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we use this _density estimator_ to build the posterior distribution $p(\\theta|x)$, i.e., the distributions over paramters $\\theta$ given observation $x$. \n",
    "\n",
    "The `posterior` can then be used to (among other features which go beyond the scope of this introductory tutorial) sample parameters $\\theta$ from the posterior via `.sample()`, i.e., parameters that are likely given the observation $x$. \n",
    "\n",
    "We can also get log-probabilities under the posterior via `.log_prob()`, i.e., we can evaluate the likelihood of parameters $\\theta$ given the observation $x$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Posterior conditional density p(θ|x) of type DirectPosterior. It samples the posterior network and rejects samples that\n",
      "            lie outside of the prior bounds.\n"
     ]
    }
   ],
   "source": [
    "posterior = inference.build_posterior(density_estimator)\n",
    "\n",
    "print(posterior) # prints how the posterior was trained"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualisations of the inferred posterior for a new observation\n",
    "\n",
    "Let's say we have made some observation $x_{obs}$ for which we now want to infer the posterior:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note: this is where your experimental observation would come in. For real observations, of course, you would not have access to the ground truth $\\theta$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta_true = prior.sample((1,))\n",
    "# generate our observation\n",
    "x_obs = simulator(theta_true)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Given this observation, we can sample from the posterior $p(\\theta|x_{obs})$ and visualise the univariate and pairwise marginals for the three parameters via `analysis.pairplot()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5e7e4586f08d434ba653da3bcb266b38",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Drawing 10000 posterior samples:   0%|          | 0/10000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAIjCAYAAAD87kn6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAtmUlEQVR4nO3df2xd5Z3n8c9zzr3Xvr52bMcmpjhhyUyW2YVQIkySTdNlGciARLtRKBSGP8IviZmsElSUioZ/UNBsEUgwCNGwTnekhUwEQ4QEgwiVRpCgELoEGHupyOwopFMI1Gx+GGwnse+Pc8959o9r3zgkkNi+9rnP9fslWTEntvvNUdR3nvPTWGutAACAU7y4BwAAABNHwAEAcBABBwDAQQQcAAAHEXAAABxEwAEAcBABBwDAQQQcAAAHJeIeAKi0v/B+GvcIqHJvRi/P2P8Wfx9xLpP9+8gKHAAABxFwAAAcRMABAHAQAQcAwEEEHAAABxFwnKZvMKu+wWzcYwAAzoGAo6xvMKtVf7tHq/52DxEHgCpHwFE2MFxQNgiVDUINDBfiHgeobcaUPoBJIuAAKuNsMSJSZ8c+QQXwJDYAU2eMZDxJ0Te2j60RRrdbO5NTVYez/sPGk/GMbGQlRbNzv2DKCDiAyRuNk0kkJW80VNFojDwj4/tSFMlaK4WhJMmG4awKlldXV/7cjv65zdh+k2SLxVLI7Sz+Rw4mhYADmBrjnYq1JJnREPm+5HmS78uEoaw0GvcwpkFjkkyWPzXRuCMUnidFo6vvqFjajzY6tWIn5DgHAg5JpSvQ9/cNxT0GXGJGo2280iozmSjFx/fLK8yxGNkgkAkjKQwVZcdW4bNjxellGkqfGFOKtlQOtw0jmXxeduxIRRhJNjp1aF2q+f2DySPgUN9gVisf3x33GHDJaLxNIiElkzJzmmTTpYhHDSlZz8hYlUIVSV4+kAqBTDGUN3RCNp+XDUPZoKhaPwccXdgmSbJJX2F9QvIkLx/KyxWlYiRvOCuTy8tGkUw2J1ssSmP7xkbn+OmYzQg4Trtl7KYrLtRvPj4c4zSoeuXVtSd5nozvSamkbF1Ktj6hYmNK1jcykZUp2lKcfSPjeVJQlFJJqVgsnf8Niqd+Zo1GPGwoHUIP6xMK056sMUokPPmekSlGMmEoI8kERdlCIFkrE0Wqzb2BSiLgOM3cTCruEVDNRq82N56RSSZkUkmZdFphW5OCOSkVG3xl23yFKY2uwCUvlOoHQ6UGi/LyoZKFoHRBWyGQ8YuyoWp6pXlyQVqSVKw3ChqNrCclcr78XEp+wSrdn1LiREFeLijd15vLlY5wjF74ZyMrRbPsugGcFwKO0yy+qDnuEVDNjFc67+0Zmbo6mbqUbFNGIxelNdLuq9BiNNwZKWoYO38rmcAo80VCmf/nKTli5eUa5UeRrJ+XKRRGz/l6KpW89gxe6skaqZixKrYEkm9lcr68nJGfM8r8MaX0VwmlTkRKFyMZ3yvtmzAs/YOpWKzlf99gCgg4ynbe/8O4R4AjjDGl28YSCSmZULHeU5g2CjJS1FJUsmH0tIw1Kga+ggFfxXojE0o2OXpluufJerX/LKmgqXQwvDgnVKolL9+PlM8lFWYTiuo8FQd8FUc8ecXSeXLjeaf2jbWlUxU1fIoBk0fAcVZf8ShVfBuvdDW1qa+XzaQVNtUp32yUb5XyHaEu/5M+/Uljv/JRQtkwqcFCgz4euVjJ4YSiROmcsF+XPLXCrHFmwYgk6cKWk/p+25dK+4H+ONKiw8NzNDiSVnZ4jkxoZD1P6UxKJghlokjGGNlZsH8weQQcp2nNpJRO+lq3vUdv/fy/qLMlHfdIqDJm9FYx1aUUNdWr0JxSfq5Rrj1SY+dx/bf5b2t53Vcaiqy+iur0ZbFV//34TTo+2KbINypmfCVTCSlIlG6rMp5q+d7w/3TJp5Kkq+Z8rv/auF8NRvpdoU0f5+br9yPz9OaJ/6hcsU4yRsXGpLx8UV5QLO3jsHb3C6au9o9fYUI6W9LauraLF5rgnKxnTn0YSb5V0g/V4o2o2atXs2fU4hXU4o2oPlmU9SXrqbSqnEUry8ZEQXMSeTV5OTV5Rk1eQk1eVg1eQY1+Xl7CyvpS5EuRb2R9b1btH0weK3CcoY0r0fEd7OhtTmb0kalm7NRsJBWKCf1bYZ5avC8U2IRyNqGvwkYFof/NH3LmD67R87xBVPqzHwma9VG+RfVeoD8U5qk/aNJg0KCoaORZSUaKUkZRXUJ+MlE+VVF+RG2N7h9MHgEHMDGRLYfEWJXuW7aSrFEQ+PpjoU0t/oj80UeqnojqVQy98tebse+NrGTM6ENfTK1ehK6iLf3ZB4oN+r/5TiVNqP6gSV8FGQ0F9bKhV376bJQwCus8JZLjnmYHfAsCDmBiPHNqNRhFMqGVn5f8rFFhJKVPs+1KekVFo+HqDxp1crhefk7yCio93CUs/SPAWnvq5Sc16vOTrZKkE3V1ykcJJU2or4OMhgr16s82ymR9+XkjryB5RSsvtOOekc79Y/h2BBzAuY29PWvspSWjL+Iw+VCJ4aIajiTlFTzlh+v0lv0P2pv5ExWLvqKCL1vwlP4iqYbDVsmTkVIDeXkjOSmXP/WGsnGr+lrzeW+nJOnTTKT/01yQ8ayiwJMt+DJ5T5nPfdX3WyWHrer7C/KHC/JO5hSNPfCmxv+Bg8kj4OCWMZy/0XuSjTFSGMkERflZT3UDKXlFyS8YRYk6hfUpJUZX5V5Rqu+3Sg+E8rOR/OGCVAhki+HoSrO2V5lNfyj9GjT5yrfWy3pSMpRM0cjPSfXHSvsmMRLJP5mXGclL+ULpKWzW1vz+weQR8FmubzCrddt7lE76av3GxWuEHWewpTeKWUmmWJQpBJJnlBwpSkaynq+6r42iOiMvkLx86bBw3YlIiZOh/FzpEaoqFqUoLL8fu5bVHS/9Gb1QkjWyvuQVJVOU/MLovhmJ5I8UZXJB6ZnoxeLoi2Bqf/9g8gj4LDcwXFA2CLXt3mXle765FxxnGA2tHX02t/GKijyvFPFcUskwUiKVVH19QumjdYqSpbeRmdDKRFaJ43l5J3KlOJ0cKb+NrPys7xpeZTYfHJZUeplJscGX9Yy8sPSiFy+M5J8slN5MVghkTo7IBqU3t5VfLRqGNXt6AVNDwCHp9FvHxu4Fv+t/faCB4QIBxynWSjaUtUYmCEpXkYdR6W1aqaRsLqG6bCDrj17oZowURfJO5qSRbClK+fypFeYs4A2cLP2aTChRl5Q1pvQGstEL+VQorbpVDGVHrwuwxeK494IDZ0fAcVbcC45zsWFUOhRurRSMXtRWDKViWHp1qEoPezFhJOXysvnRF5cUi6VVZWRrfvUtSWY4W7robzTSRirtk7ErzQuBbBSdCncYlo901Pq+wdQQcAATZ+3ood3RwASB5Pulc+O+P+7LbOm91kFw6pax8SvLWfCazGjoeOkT35dJlB7QYsOw9A+gcbeL2fGvDx3brxw6x3cg4AAmx0alh68YTwqjU7EZuzVs/D3e3zwcPItWljYolj4JR6+697zyKlvf2Cd27NnnhBvngYADmBxrRx/oUgq5saV7xa1UivpopO24iJ/x/bPAWJSNNSr/icetss84zz1L9gumjoADmLyx2Ixe2HYKK8myaOxhNfr2t4uxnzAJBBxAZRChcxu/j8aedc5+wyQRcACIA+HGFPE+cAAAHETAAQBwEAEHAMBBBBwAAAcRcAAAHETAAQBwEAEHAMBBBBwAAAcRcAAAHETAAQBwEAEHAMBBBBwAAAcRcAAAHETAAQBwEAEHAMBBBBzf6avhQtwjAADOgoDjrFozKaWTvtZt71HfYDbucQAA30DAZ7lvW2F3tqS1dW2XskGoAVbhAFB1CPgs1jeY1brtPUonfbVmUmf8fttZtgEAqkMi7gEQn4HhgrJBqG33LlNnSzrucQAAE8AKHKy0AcBBBBwAAAcRcAAAHETAAQBwEAEHAMBBBBwAAAcRcAAAHETAAQBwEAEHAMBBBBwAAAcRcAAAHETAAQBwEAEHAMBBBBwAAAcRcAAAHETAAQBwEAEHAMBBBBwAAAcRcAAAHETAAQBwEAEHAMBBBBwAAAcR8Fnsq+FC3CMAACaJgM9SfYNZrdveo3TSV2sm9Z1fS+gBoPoQ8FlqYLigbBBq69oudbakz/o1rZmU0klf67b3qG8wO8MTAgC+CwGf5dq+Y/Xd2ZLW1rVdygahBliFA0BVIeD4Tt8VeABAfAg4AAAOIuAAADiIgAMA4CACDgCAgwg4AAAOIuAAADiIgAMA4CACDgCAgwg4AAAOIuAAADiIgAMA4CACDgCAgwj4LMU7vgHAbQR8FuobzGrd9h6lk75aedsYADgpEfcAmHkDwwVlg1Db7l2mzpZ03OMAACaBFfgsxru+AcBdBBwAAAcRcJwXLnoDgOpCwPGdWjMppZO+1m3vUd9gNu5xAACjCDi+U2dLWlvXdikbhBpgFQ4AVYOA45zGLnbjMDoAVA8CjnPiMDoAVB8CjnPiMDoAVB8CjvPCPeMAUF0IOAAADiLgAAA4iIADAOAgAg4AgIMIOAAADiLgAAA4iIADAOAgAg4AgIMIOAAADiLgsxAvJQEA9xHwWaZvMKt123uUTvpq5fGoAOCsRNwDYGYNDBeUDUJtu3eZOlvScY8DAJgkVuCzFC8nAQC3EXAAABxEwAEAcBABBwDAQQR8luEWMgCoDQR8FuEWMgCoHdxGNotU4hYyVvAAUB1Ygc9Ck7mFrDWTUjrpa932HvUNZqdhKgDARBBwnJfOlrS2ru1SNgg1wCocAGJHwHHeePgLAFQPAj6LcP4aAGoHAZ8luAIdAGoLV6HPEpV8iQkreQCIHyvwWWYq57G5Eh0AqgcBx3njSnQAqB4EfJao1GFvrkQHgOpAwGeB6biAjfPgABAvAj4L/P7oSWWDUFvXdk35Arbx58H39w1VaEIAwEQR8BrUN5gtf+z55Fh59b1oXuOUf3ZnS1ovr1shSfrp1vf0T/9ymAvaACAG3EZWY/oGs1r1t3uUDcLytnTS18vrVkx59T1mcWezXl63Qj/+1bv66+09kqTfPnRdxX4+AODcjLXWxj0EAACYGA6hAwDgIAIOAICDCDgAAA4i4AAAOIiAAwDgoCnfRmat1YkTJyoxC2pcU1OTjDFxjwEANWHKAe/v79e8efMqMQtq3NGjR3XBBRfEPQYA1IQpBzyVKj1b+4svvtCcOXOmPJDrjh8/rgULFrA/xhnbJ2N/VwAAUzflgI8dEp0zZw7BGof9cSYOnwNA5XARGwAADiLgAAA4aMoBr6ur0+bNm1VXV1eJeZzH/jgT+wQAKo+XmQAA4CAOoQMA4CACDgCAgwg4AAAOIuAAADhoSgEPgkCbNm3SFVdcoUwmo4suukh33nmnvvzyy0rN55xnn31Wl1xyierr67V8+XJ98MEHcY8Ui8cee0xLly5VU1OT5s2bpzVr1ujAgQNxjwUANWNKAR8ZGVFvb68efvhh9fb26pVXXtGBAwe0evXqSs3nlB07dmjjxo3avHmzent7deWVV+rGG2/U0aNH4x5txu3Zs0fr16/Xvn379OabbyoIAt1www0aHh6OezQAqAkVv43sww8/1LJly3To0CFdfPHFlfzRVW/58uVaunSptmzZIkmKokgLFizQ/fffr4ceeijm6eJ17NgxzZs3T3v27NE111wT9zgA4LyKnwMfGhqSMUYtLS2V/tFVrVAoqKenR6tWrSpv8zxPq1at0nvvvRfjZNVhaGhIkjR37tyYJwGA2lDRgOdyOW3atEl33HHHrHuRR39/v8IwVEdHx2nbOzo6dPjw4Zimqg5RFOmBBx7QypUrtXjx4rjHAYCaMKGAv/DCC2psbCx/7N27t/x7QRDotttuk7VW3d3dFR8U7lq/fr3279+vl156Ke5RAKBmTOh1oqtXr9by5cvL/93Z2SnpVLwPHTqk3bt3z7rVtyS1t7fL930dOXLktO1HjhzRhRdeGNNU8duwYYN27typd955R/Pnz497HACoGRNagTc1NWnRokXlj3Q6XY73wYMH9dZbb6mtrW26Zq1qqVRKXV1d2rVrV3lbFEXatWuXVqxYEeNk8bDWasOGDXr11Ve1e/duLVy4MO6RAKCmTGgF/k1BEOjWW29Vb2+vdu7cqTAMy+d7586dq1QqVZEhXbFx40bddddduvrqq7Vs2TI9/fTTGh4e1j333BP3aDNu/fr1evHFF/Xaa6+pqamp/PeiublZ6XQ65ukAwH1Tuo3ss88++9aV1dtvv61rr712sj/aWVu2bNETTzyhw4cPa8mSJXrmmWdOO+0wWxhjzrr9ueee09133z2zwwBADeJ1ogAwjf7C+2ncI6DKvRm9PKnv41noAAA4iIADAOAgAg4AgIMIOAAADiLgACqubzCrvsFs3GMANW1K94EDwDf1DWa16m/3KBuEkqSd9/9QizubY54KqD2swAFU1MBwoRxvSdrfNxTjNEDtIuAAADiIgAOYVq2Z2fVIZWCmEPApstbqqaee0sKFC9XQ0KA1a9ZoaIhDhsCYzhaefQ9MBwI+RQ8++KC6u7u1bds27d27Vz09PXrkkUfiHgsAUOMI+BS8//77euqpp7Rjxw5dc8016urq0n333aff/OY3cY8GAKhxBHwKnnzySV1//fW66qqryts6OjrU398f41QAgNmAgE9SPp/XG2+8oZtvvvm07blcTs3N3PMKAJheBHySent7lc1m9fOf/1yNjY3lj1/84he69NJLJUk333yzWltbdeutt8Y8LQCg1hDwSfrkk0+UyWT08ccf66OPPip/XHLJJVq5cqUk6Wc/+5n+/u//PuZJgZn11XBBkvT4T66IeRKgthHwSTp+/Lja29u1aNGi8kcymdTBgwd1yy23SJKuvfZaNTU1xTwpMHP6BrNat71H6aTP/d/ANONZ6JPU3t6uoaEhWWtljJEkPfroo7rpppt02WWXxTwdEI+xx6huu3eZ2kYDPrYiB1BZrMAn6brrrlMul9Pjjz+uTz/9VL/85S/1+uuvq7u7O+7RgNi1ZVJqzaSUTvpat72HN5MB04CAT1JHR4eef/55dXd36/LLL9e+ffv07rvvasGCBXGPBlSFzpa0tq7tUjYINcAqHKg4DqFPwe23367bb7897jGAqtXGeXBg2hDwabRq1Sr97ne/0/DwsObPn6+XX35ZK1asiHssAEANIODT6K233op7BABAjeIcOABMJ2NKH0CFsQIHgJkwFnFrTw+6tfHMA+cRcACYTuYbBzrPWIxHpV8IOSaIgAOYdrP5YS7GM2OfyPilmNvxsQ5D2ciWwh6FMz8gnMU5cADThoe5qLQCN57kGcnzJN+X8X0ZY0pPcTTeqcgDE0DAAUwbHuYieY2Z0secOfKa58hraS7/aprnlH6voUFeul4mmZJJJCTPj3tsOIBD6ACm1Wx/mItpbZYk2WRCUTola4xMFMkUI8lamWxeKoYyQSA7kpUtFkuH1QsR58XxnQg4AEwjW5cs/xo2JCVjZMLRgIdWXhiVX4gkv1CKehTJGk8SEce3I+AAMI1O/lmrJClIeyo0GVlP8vNSIm/lBVZ1A/VKjBTlZYvyPSPl8lLekwkj2VAi4vg2BBwAptHAvy/932zQZFVojWR9q8RJX4lhT15BajjiqW4ooeTJUOlCUcbzZCWZXF6SShG3XJ2OMxFwAJhGQePor01Wag7k+1ZF38r6vry8UfGEUSJn5BU92aQvJfxSxD1PCkfvETeGVTjOQMABYBpFf3ZSxkjfaz6pq9q/UNoP1Jdt0bFso/pHGjRUN1fFtKe6QaPkibSSCU+eJJPNlsJdkGwxivuPgSpEwAFgGv3nS/4gSbpqziHdlPlXNXlGfywm1Bc260DuIv3P7ErlwkZFCaN0f0ImskoWiqVbyiIr6xVLP4hVOL6BgAOomNn8xLVvk0mUzmVnvLwynlGDSarBC1QfBarzAvl+pKJXej669VR6ItvYVele6WEvZBtnQ8ABVETfYFbrtvconfTVOsvv/R6vLTks30Rq8UdUZzwlja+kCkqaUEkTyvcj2YRkv/nsFs+MPsHNO/U8dS5mwzgEHEBFDAwXlA1Cbbt3mTpb0nGPUzWaEyOSpAaTV1K+PBn5RkopVMoU5Rsr69lTq2+pdK+4MbJe6VWkxjOjz0vnMDpOIeAAKmq2P3ntm3qPXyzfWJ0M6yX9mzJeXp8F/05/LMzVoVybhoYalDjpKTEi+TkrrxDJBKFsGElhKFlbirdEvHEaAg4A02jvv1wqGem9pkv09txLlTCR+kcadGK4XkE2qfQf6lR/zCp1IlLdQEH+iZzMSE42KMgGRdmgKFke5oIzEXAAmEaJ/tKjVAs5X58XfXl+pMLJlMzJhBI5o9SgVHfcKjkcycsFMvmiFBTHrcCJN86OgAPANGo4XDqxXRzyVRxMy3pSesTIz6n0JLZjkeoGQyWyRXkjBZl8QTYISi80GX/4HPgGAg4A02juv5ZurbOekfWNZCRTtPKCSF5olRjMyxvJywRF2eER2WJRthCUDp+PrcCBsyDgADCNkseDM7aZMJKJrEwxkhnOyeQKUrEoFUorbxFunAcCDgDTKPHVyTM3RqWAy1oplx89ZB7JFkqvE7VBsXTonPPf+A4EHACm0+FjZ90cjYU5DKUokrW29LlUinfEQ1vw3Qg4AEwjWxh9vKznnf4bUekQubVWGl1tn7rfm8PnODcCDgDTyBaLp/6j/EjUU4E+7Srzse0cNsd5IOAAMI1OD/jos1IJNCqAgAPATCHcqCDv3F8CAACqDQEHAMBBBBwAAAcRcAAAHETAAQBwEAEHAMBBBBzAjPhquBD3CEBNIeAAplVrJqV00te67T3qG8zGPQ5QMwg4gGnV2ZLW1rVdygahBliFAxVDwAFMu7ZMKu4RgJpDwAEAcBABBwDAQQQcAAAHEXAAABxEwAEAcBABBwDAQQQcAAAHEXAAABxEwAEAcBABBwDAQQQcAAAHEXAAFcHrQoGZRcABTFnfYFbrtvconfTVyotLgBmRiHsAAO4bGC4oG4Tadu8ydbak4x4HmBVYgQOoGF4bCswcAg4AgIMIOIApO98L2LjQDagcAg5gSs7nArbWTErppK9123vUN5id4QmB2kTAAUzJ2AVsW9d2fesFbJ0taW1d26VsEGqAVThQEQQcQEWc6wI2LnADKouAAwDgIAIOAICDCDgAAA4i4AAAOIiAAwDgIAIOAICDCDgAAA4i4AAAOIiAAwDgIAIOAICDCDiAGcUbyYDKIOAAZgRvJAMqi4ADmBG8kQyoLAIOYMbwRjKgcgg4AAAOIuAAADiIgAMA4CACDgCAgwg4AAAOIuAAADiIgAOYEp6sBsSDgAOYtL7BrNZt71E66auVe7yBGZWIewAA7hoYLigbhNp27zJ1tqTjHgeYVViBA5gynrAGzDwCDgCAgwg4AAAOIuAAADiIgAOYcdx6BkwdAQcwY1ozKaWTvtZt71HfYDbucQCnEXAAM6azJa2ta7uUDUINsAoHpoSAA5hR3HIGVAYBBwDAQQQcAAAHEXAAABxEwAEAcBABBwDAQQQcwKTxQBYgPgQcwKTwLnAgXrwPHMCk8C5wIF6swAFMCQ9mAeJBwAHEgvPnwNQQcAAziheaAJVBwAHMKF5oAlQGAQcw4zhvDkwdAQcwKZzDBuJFwAFMGPeAA/HjPnAAE8Y94ED8WIEDmLSpnsvmMDwweQQcwIzjVjJg6gg4gBnHrWTA1BFwALHgVjJgagg4gAnj3DUQPwIOYEK4hQyoDtxGBmBCuIUMqA6swAFMyNjh80qdw/7xr97V/r6hivwsYDZhBQ7gnMZu9RoYLlTs8Pn479/zyTFlg1AXtaRZ1QPnyVhrbdxDAKhe//r/jusn/+N/KxuEkqR00tfL61ZocWfzlH/2h599rZ9ufe+0bb996LqaivhfeD+NewRUuTejlyf1fQQcAAAHcQ4cAAAHEXAAABxEwAEAcBABBwDAQQQcAAAHcR84aoq1VidOnIh7DDigqalJxpi4xwAmjYCjpvT392vevHlxjwEHHD16VBdccEHcYwCTRsBRU1Kp0tO9vvjiC82ZMyfmaeJ3/PhxLViwgP0xztg+Gfu7AriKgKOmjB0SnTNnDsEah/1xJg6fw3VcxAYAgIMIOAAADiLgqCl1dXXavHmz6urq4h6lKrA/zsQ+Qa3gZSYAADiIFTgAAA4i4AAAOIiAAwDgIAIOAICDCDhqRhAE2rRpk6644gplMhlddNFFuvPOO/Xll1/GPVpsnn32WV1yySWqr6/X8uXL9cEHH8Q9Uiwee+wxLV26VE1NTZo3b57WrFmjAwcOxD0WMCUEHDVjZGREvb29evjhh9Xb26tXXnlFBw4c0OrVq+MeLRY7duzQxo0btXnzZvX29urKK6/UjTfeqKNHj8Y92ozbs2eP1q9fr3379unNN99UEAS64YYbNDw8HPdowKRxGxlq2ocffqhly5bp0KFDuvjii+MeZ0YtX75cS5cu1ZYtWyRJURRpwYIFuv/++/XQQw/FPF28jh07pnnz5mnPnj265ppr4h4HmBRW4KhpQ0NDMsaopaUl7lFmVKFQUE9Pj1atWlXe5nmeVq1apffeey/GyarD0NCQJGnu3LkxTwJMHgFHzcrlctq0aZPuuOOOWfcij/7+foVhqI6OjtO2d3R06PDhwzFNVR2iKNIDDzyglStXavHixXGPA0waAYezXnjhBTU2NpY/9u7dW/69IAh02223yVqr7u7uGKdEtVm/fr3279+vl156Ke5RgCnhdaJw1urVq7V8+fLyf3d2dko6Fe9Dhw5p9+7ds271LUnt7e3yfV9Hjhw5bfuRI0d04YUXxjRV/DZs2KCdO3fqnXfe0fz58+MeB5gSVuBwVlNTkxYtWlT+SKfT5XgfPHhQb731ltra2uIeMxapVEpdXV3atWtXeVsURdq1a5dWrFgR42TxsNZqw4YNevXVV7V7924tXLgw7pGAKWMFjpoRBIFuvfVW9fb2aufOnQrDsHy+d+7cuUqlUjFPOLM2btyou+66S1dffbWWLVump59+WsPDw7rnnnviHm3GrV+/Xi+++KJee+01NTU1lf9eNDc3K51OxzwdMDncRoaa8dlnn33ryurtt9/WtddeO7MDVYEtW7boiSee0OHDh7VkyRI988wzp512mC2MMWfd/txzz+nuu++e2WGACiHgAAA4iHPgAAA4iIADAOAgAg4AgIMIOAAADiLgAAA4iIADAOAgAg4AM2R/35AueegN7e8binsU1AACDgAzZCzcez45FvMkqAUEHABmSGum9DjfJ/7pgPoGszFPA9cRcKCGWGv11FNPaeHChWpoaNCaNWs0NMTh2mpRn/TLnw8MF2KcBLWAgAM15MEHH1R3d7e2bdumvXv3qqenR4888kjcY0FS32BW67b3xD0GaggBB2rE+++/r6eeeko7duzQNddco66uLt133336zW9+E/doUGnFnQ1CPXjjn0mSfvyrd7mYDVNCwIEa8eSTT+r666/XVVddVd7W0dGh/v7+GKfCN7VlTr3WloBjKgg4UAPy+bzeeOMN3Xzzzadtz+Vyam5ujmkqANOJgAM1oLe3V9lsVj//+c/V2NhY/vjFL36hSy+9VF988YWuvfZaXXbZZfr+97+vl19+Oe6RAUxRIu4BAEzdJ598okwmo48++ui07T/60Y+0cuVKJRIJPf3001qyZIkOHz6srq4u3XTTTcpkMvEMDGDKCDhQA44fP6729nYtWrSovO3QoUM6ePCgbrnlFn3ve9/T9773PUnShRdeqPb2dn399dcEHHAYh9CBGtDe3q6hoSFZa8vbHn30Ud1000267LLLTvvanp4ehWGoBQsWzPSYACqIFThQA6677jrlcjk9/vjj+su//Eu98MILev311/XBBx+c9nVff/217rzzTv3d3/1dTJMCqBRW4EAN6Ojo0PPPP6/u7m5dfvnl2rdvn959993TVtn5fF5r1qzRQw89pB/84AcxTju7tY67jWz858BEEXCgRtx+++36/PPPNTIyop07d+pP//RPy79nrdXdd9+t6667TmvXro1xSnS2pPXrtV3lz4HJIuDALPDb3/5WO3bs0D/+4z9qyZIlWrJkiT7++OO4x5pVvhr37HPCjUrgHDgwC/zwhz9UFEVxjzFrjT0HPZ301ZpJlV9k8hUvNMEUsAIHgGk29hz0rWu71NmSVmsmpXTS17rtPbxWFJNGwAFghow9B72zJa2ta7uUDUJeK4pJI+AAEIM2rkDHFBFwAAAcRMABAHAQAQcAwEEEHAAABxFwAAAcRMABAHAQAQcAwEEEHAAABxFwAAAcRMABAHAQAQcAwEEEHACm0f6+If34V+/GPQZqEAEHgGm0v28o7hFQowg4AMToK14nikki4AAQg9ZMSumkr3Xbe9Q3mI17HDiIgANADDpb0tq6tkvZINQAq3BMAgEHgJi0ZVJxjwCHEXAAABxEwAEAcBABBwDAQQQcAKZR67jz3K2c80YFEXAAmEadLWlJ0q/XdpU/ByqBgAPADCDeqDQCDgDTiCetYboQcACYRuu29yid9Dn/jYpLxD0AANSybBBq273LOISOimMFDgDTjCeuYToQcAAAHETAAQBwEAEHAMBBBBwAAAcRcAAAHETAASBmP/7Vu9rfNxT3GHAMAQeAmIx/uEvfYDbGSeAiAg4AMelsSevXa7vKnwMTQcABIEaEG5NFwAEAcBABB4AqwFvLMFEEHABi1JpJKZ30tW57DxeyYUIIOADEqLMlra1ru5QNQg2wCscEEHAAiBlvK8NkEHAAABxEwAEAcBABBwDAQQQcAAAHEXAAABxEwAEAcBABBwDAQQQcAAAHEXAAABxEwAEAcBABBwDAQQQcAAAHEXAAABxEwAEAcBABBwDAQQQcAAAHEXAAABxEwAEAcBABBwDAQQQcAAAHEXAAABxEwAEAcBABBwDAQQQcAAAHEXAAABxEwAEAcBABBwDAQQQcAAAHEXAAABxEwAGgSnw1XIh7BDiEgANAzFozKaWTvtZt71HfYDbuceAIAg4AMetsSWvr2i5lg1C/P3oy7nHgCAIOAFVg0bxGVuGYEAIOAFVg/Cp8gHPhOA8EHACqRFsmFfcIcAgBBwDAQQQcAAAHEXAAABxEwAEAcBABBwDAQQQcAAAHEXAAABxEwAEAcBABBwDAQQQcAAAHEXAAABxEwAEAcBABBwDAQQQcAAAHEXAAABxEwAEAcBABBwDAQQQcAKrMV8OFuEeAAwg4AFSJ1kxK6aSvddt71DeYjXscVDkCDgBVorMlra1ru5QNQg2wCsc5EHAAqCJtmVTcI8ARBBwAAAcRcAAAHETAAQBwEAEHgCr041+9q/19Q3GPgSpGwAGgirSOu4iNW8nwXQg4AFSRzpa0fr22q/w58G0IOABUGcKN80HAAaBK8UhVfBcCDgBVhkeq4nwQcACoMjxSFeeDgANAFRp7pCorcHwbAg4AVWgs3H/NYXR8CwIOAFVo/KHz3x89GeMkqFYEHACq0OLO5vLnXMyGsyHgAFCFFnc267PHf6Rt9y5TNgj14adfa3/fkI4ez8U9GqqEsdbauIcAAAATwwocAAAHEXAAABxEwAEAcBABBwDAQQQcAAAHJeIeAABqlbVWJ06ciHsMOKCpqUnGmAl9DwEHgGnS39+vefPmxT0GHHD06FFdcMEFE/oeAg4A0ySVKr2Q5IsvvtCcOXNiniZ+x48f14IFC9gf44ztk7G/KxNBwAFgmowdEp0zZw7BGof9caaJHj6XuIgNAAAnEXAAABxEwAFgmtTV1Wnz5s2qq6uLe5SqwP4401T2CS8zAQDAQazAAQBwEAEHAMBBBBwAAAcRcAAAHETAAQBwEAEHgGkQBIE2bdqkK664QplMRhdddJHuvPNOffnll3GPFptnn31Wl1xyierr67V8+XJ98MEHcY8Ui8cee0xLly5VU1OT5s2bpzVr1ujAgQMT/jkEHACmwcjIiHp7e/Xwww+rt7dXr7zyig4cOKDVq1fHPVosduzYoY0bN2rz5s3q7e3VlVdeqRtvvFFHjx6Ne7QZt2fPHq1fv1779u3Tm2++qSAIdMMNN2h4eHhCP4f7wAFghnz44YdatmyZDh06pIsvvjjucWbU8uXLtXTpUm3ZskWSFEWRFixYoPvvv18PPfRQzNPF69ixY5o3b5727Nmja6655ry/jxU4AMyQoaEhGWPU0tIS9ygzqlAoqKenR6tWrSpv8zxPq1at0nvvvRfjZNVhaGhIkjR37twJfR8BB4AZkMvltGnTJt1xxx2z7k1c/f39CsNQHR0dp23v6OjQ4cOHY5qqOkRRpAceeEArV67U4sWLJ/S9BBwAKuCFF15QY2Nj+WPv3r3l3wuCQLfddpusteru7o5xSlSb9evXa//+/XrppZcm/L28DxwAKmD16tVavnx5+b87OzslnYr3oUOHtHv37lm3+pak9vZ2+b6vI0eOnLb9yJEjuvDCC2OaKn4bNmzQzp079c4772j+/PkT/n5W4ABQAU1NTVq0aFH5I51Ol+N98OBBvfXWW2pra4t7zFikUil1dXVp165d5W1RFGnXrl1asWJFjJPFw1qrDRs26NVXX9Xu3bu1cOHCSf0cVuAAMA2CINCtt96q3t5e7dy5U2EYls/3zp07V6lUKuYJZ9bGjRt111136eqrr9ayZcv09NNPa3h4WPfcc0/co8249evX68UXX9Rrr72mpqam8t+L5uZmpdPp8/453EYGANPgs88++9aV1dtvv61rr712ZgeqAlu2bNETTzyhw4cPa8mSJXrmmWdOO+0wWxhjzrr9ueee0913333+P4eAAwDgHs6BAwDgIAIOAICDCDgAAA4i4AAAOIiAAwDgIAIOAICDCDgAAA4i4AAAOIiAAwDgIAIOAJgW1lo99dRTWrhwoRoaGrRmzRoNDQ3FPVbNIOAAgGnx4IMPqru7W9u2bdPevXvV09OjRx55JO6xagbPQgcAVNz777+vFStW6J//+Z911VVXSZL+5m/+Ri+88IIOHDgQ83S1gRU4AKDinnzySV1//fXleEtSR0eH+vv7Y5yqthBwAEBF5fN5vfHGG7r55ptP257L5dTc3BzTVLWHQ+gAgIp677339IMf/ED19fXyfb+8PQgC/fmf/7leeuklrVq1SsViUcViUT/72c903333xTixmxJxDwAAqC2ffPKJMpmMPvroo9O2/+hHP9LKlSvV1NSkd955Rw0NDRoeHtbixYv1k5/8RG1tbfEM7CgCDgCoqOPHj6u9vV2LFi0qbzt06JAOHjyoW265Rb7vq6GhQVLpcLu1VhwMnjjOgQMAKqq9vV1DQ0OnRfnRRx/VTTfdpMsuu0ySNDg4qCuvvFLz58/Xgw8+qPb29rjGdRYBBwBU1HXXXadcLqfHH39cn376qX75y1/q9ddfV3d3d/lrWlpa9Lvf/U6ffvqpXnzxRR05ciTGid1EwAEAFdXR0aHnn39e3d3duvzyy7Vv3z69++67WrBgwVm/9sorr9TevXtjmNRtXIUOAJhRR44cUUNDg5qamjQ0NKSVK1fqH/7hH3TFFVfEPZpTuIgNADCjDh06pL/6q78qX7x2//33E+9JYAUOAICDOAcOAICDCDgAAA4i4AAAOIiAAwDgIAIOAICDCDgAAA4i4AAAOIiAAwDgIAIOAICDCDgAAA4i4AAAOOj/A4Ap7W8VhxahAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 600x600 with 9 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "samples = posterior.sample((10000,), x=x_obs)\n",
    "_ = pairplot(samples, limits=[[-2, 2], [-2, 2], [-2, 2]], figsize=(6, 6),labels=[r\"$\\theta_1$\", r\"$\\theta_2$\", r\"$\\theta_3$\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Assessing the posterior for the known $\\theta, x$ - pair "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this special case, we have access to the ground-truth parameters that generated the observation. We can thus assess if the inferred distributions over the parameters match the parameters $\\theta_{true}$ we used to generate our test observation $x_{obs}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a3d1719fd43441e3809bae68c87c6429",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Drawing 10000 posterior samples:   0%|          | 0/10000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAIjCAYAAAD87kn6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAsvElEQVR4nO3df3DUdZ7n8de3v92ddH6QBEKihDiyxzq3igNjBAqZcl1l9dbZo3B0nHGq/FnnDFdgaeE5WFfl4k6to7V6nutiwdZWrTIUKueWriPO7a2ChTgH4iSjNe5uIXOrWQYGMJIESPrHt7/fz/3RSRN+xCGdTn/70/18VGUNDcm891voM59vf77fr2OMMQIAAFaJhD0AAACYOAIOAICFCDgAABYi4AAAWIiAAwBgIQIOAICFCDgAABYi4AAAWCga9gBAsf1x5Nthj4Ay91bwSsn+t/j7iN+l0L+PrMABALAQAQcAwEIEHAAACxFwAAAsxCY2nGKM5A3nPo/VSY4T7jwAgHGxAscp3rD041m5j9GQAwDKEgEHAMBCBBwAAAvxHjiAqTF2D4Ux4c0RttHjUM3HAFOCgAOYvLGxdiJjPnVkAiMpqN6AjR6PM/eEmmDM51V6bDApBBxA4UbC7bhuLlQRZ+RzJx8lJwhkjJECI5lAxverKliReGzkk0juuEQiku/njoMk42UlBWd/YRUdIxSGgAMozOiq24mcFm8nHsv92gRSYGSMkeP7khPIGGck5H64s5eS60oa+SEnNvKfXM859UNOYM4+HOYcQQfOQMCRd2gwqVljP59ZH+o8KGNj4u24rpxYVE40KsWichIJKerm4uQHcoyRyXhSNivH9xVIMp5ykaqCVabTkPv3aPT4yHEkLyt5nhQYOem0TDYrBYGMH+TOUgSR6vohBwUh4MjrH/LyAe8f8jRrZqjjoMw5riu5rpxoVE5dnZx4TKY+IW9Gg0w8IvlGTmAkI7lDaTnJjBwvq8jgCZlkUsYPZLJexUfctE+XJAU1Mfl1URnHUSTtK5LK/UATOZGUk0rLZH1peFgmm5UjP3e2osKPDSaHgAMozOjqe/S0eTwmUxtTtj4qv9aVExg5fu5Dklwp9/5vNJoLvyTjR1TpG9yCurgkya+NKlvvykQcuVFHruvI8aNysoEcSU7Gk9Ju7v1xJyKJFTi+HAEHMDGOk4t3PCantkZOba38thb5DXFlmmI6OcuVX+vICSQFUsQ3qu2PKT4YU3TYV8z3cxuyvYycbDZ3ptip3NXm0OyEJClTH1GmyZGJSG7KlZuW3IxR3ecxxQbSiqS83I05Uikp40nJILeDv0reasDEEXAA52803q6bO23eUKdgWp1O/l6DkjMiSrc4GvpKVkr4ujjTr1tPfKTZ3qAOt0/Xm18s0OHB6ZqWnaZYYKSkmzt1PBqpCnXsq67kSB0NffpWrFud/oAOOC16Jb5A/56ZqcZP46o/6ip2MlCtpMjJuJzhpALPk/wg9wMO74fjHAg4gII4bkQmFpWJR+XVReQ1OPIajdwmT99LdevPD26XUe7yZyPpv0z/uf48ulw/P/AfZWJROV5UikRGrhWv3E1bXqPRzdlf6vHhN047Hiu1Wz9s+ab+seFKeScicnzJxKO5YzPyNoOMkWMc3g/HOXErVQAT4kQcOW5Eio2+5x1Tcqaj5AVGmY6M/mTGP+tHvW/LlVFUJv/PiGO0btpP1d58XKZmZEf26PXjFWzOrN/o8cwbZx8PGf1l/5tqb/lCqRmOUs0R+fUxmURcisdO7ezP3wiGpwPidJX9bw6A4nIiuei6rkxNTEFtTF5DVKk2I6dzWF/9ymGtOb7znK1xJBlH+k+zPlJQG5OpiclxHClS2WH6b9p59l3YNPKSI91S84FSM43S0x15DVH5dTGZ2rgUi+ci7kYq/occFIa/FQAKF3FkXOU2ZkUD1bpZNQ6lx/3jjqQL44Onfb2UW9VXqrbkSTnjnP52jNFsfzB3DF3JREZu8MJqG+eBgAM4b07EGVk15/7T4fhGkaxRxJMy6aiOpep0ONE47tcbSYfTTbnrw/3cndpOffPKjNZn8RaZcf5/M46j37jNcnzJ8ZU7QPnfPOP4AGcg4AAmJhKR4zgjN2kxuUvFMo6ClKsTqRq92/l7kjm9RVLu146R3vx8geSbc2/KqsCI//30r33J8TB6uW5+LuCBTq3UTe74mJEP4FwIOIDCGCP5RhEvUCQjOUlXQ8M1+mV0tn565SUykgIn95GVo0CO/iyzXIeGW+QEgRw/qIo4fejP0p9dtExGTv54+I4jI+mJ3/9D/Xt2ptyUkzuGWSPHD06dnQgq9/I6TB6XkSHv2HAm7BFgCWOMnKyvSNpTdCiq+kNxuRlX6YE6bc98VbsSv6dn5w7o2/2/0qz0Cf325Az970Nf15HjLWroH5aTzOTuB+77uTuPBeOsyCuA9y/T9L+cq/Te9Ev0vcgHmp0dVK9a9HL8Ch3oa1X9wYhqvzCKnwxUcyytyMmUnOGUgkwmd49036/o6+RROAIOSdLBgaTuf+mX6nHDngRlLwhyp7qzWcnLyh3OKHGsRk7gyk05GnZqlK2pUW+2Qf8z3SnHkxoOGtX/1lMinVHkZO6e6LmHeQQVfyOX+t/k/tk3bab+cvqfyLhSNCm5SUf1GSnxuVHNcV/Rk74iQ2k5wymZVDoX7sCMPE9dFfsDDgpHwCFJ6h/KKOn5IzesBs5tdKXsGCOT9XMhzriKDvmqcSXHRJRNRBTEc5uy3IwU8aT4yUDRYV+RdFZOOrf6NtlsVZxCj58YeS66yW1aMxEpmpLclDl1bE76ig5nc7dQzfpS4J9+bKrgOGHiCDiAiRk95Z1KSYEvJVOq8bKqqYnLr4+rtr9Wfk0kv0PdyZrcqeHBYTleVubksEw6LTNy+rySV9+SNO3TpCQpqHGVrXNzTyPL5vYOOL5RdDCtyHA6d2xODOV+sMlmc19c4ccGk0PAAZy/kWdVO5FA8jyZIJDjZnPPtY66ctMJ1UQcBTFXjh/kn0gWGRyWM5TMhSmdlslkcjusq+D0cLR/WJJkYq6iJ0aeBx7k4q1soMhQUsp4uWPjZXLPBB85fZ77QiKOcyPgACZmJOLyc4/BNMaM3M/clROJyD0eUyQakQLJGQmRM5ySSadzD+fIZnORMkFVxMkZTuX+6boyUffUk9dGduIrlZYZfXDJyMY+Y8zIca7cH2wweQQcwPnLX6fs565rHgmw8bJSxJGTSstJpeRE3Nzv+bnfz++oDsaEqUoekxkcG5Ck3A1w3JFNJqPHLTAjP9CMPMhl5J8mMLm3J4AvQcABFGZ0JS5J8nNPzRpdMUYiud3qIyvN0RX3aTvOqyDekmQyI5dnjrmDnaT8Nd75sxHSmLcUKv/MBCaPgENS7jIyYEKMkXQqNMZXbqt1dszd1EaCZM7crFYl8ZZ0alUtnf1QkrGnyav1+KBgBBw6OJDUDzZ3KxH2ILDPmaExvkyGEI2V31EOFBkBh/qHuAMbiohoAyXBvdBxTtxWFQDKGwFH3t+vvCr/+f0v/ZL3xQGgjBFwnFPS8zm1DgBljIADAGAhAg4AgIUIOAAAFiLgAABYiIADAGAhAo5xfcEudAAoWwQc55SIuVq5uZtrwQGgTBFwnNNf3fZ1rgUHgDJGwHFO0+viYY8AAPgSBBwAAAsRcLBZDQAsRMCr3MGBpFZu7lYi5qqlPhb2OACA88TzwKtc/1BGSc/XpnsWaVZTIuxxAADniRU4JEkz6tm0BgA2IeAAAFiIgAMAYCECDgCAhQg4AAAWIuAAAFiIgAMAYCECDgCAhQg4AAAWIuAAAFiIgFc5HmQCAHYi4FXs9AeZcCtVALAJAa9iow8y2Xh7lzqaeZAJANiEgIMHmQCAhQg4AAAWIuAAAFiIgAMAYCECDgCAhQg4AAAWIuAAAFiIgONLcac2AChPBBzn1FIfUyLmauXmbh0cSIY9DgDgDAQc5zSrKaGNt3cp6fnqZxUOAGWHgGNc3KENAMoXAQcAwEIEHAAACxFwAAAsRMDxO3EpGQCUHwKOcbXUx7mUDADKFAHHuDqauZQMAMoVAceX4lIyAChPBLyK8d42ANiLgFepgwNJrdzcrUTMVQurbACwTjTsARCO/qGMkp6vTfcsUkdzIuxxAAATxAq8yvEeNwDYiYADAGAhAg4AgIUIOAAAFiLgAABYiIADAGAhAg4AgIUIOAAAFiLgAABYiIADAGAhAg4AgIUIOM4LTy4DgPJCwPGlWurjSsRcrdzcrYMDybDHAQCMIOD4Uh3NCW28vUtJz1c/q3AAKBsEHL8TTywDgPJDwAEAsBABBwDAQgQcAAALEfAqxWVhAGA3Al6FDg4ktXJztxIxVy1sUAMAK0XDHgCl1z+UUdLztemeRepoToQ9DgCgAKzAqxiXhwGAvQg4AAAWIuAAAFiIgAMAYCECDgCAhQg4AAAWIuAAAFiIgAMAYCECDgCAhQh4FSr0PujcPx0AygcBrzKF3Ae9pT6uRMzVys3dOjiQnOIJAQDng4BXmdH7oG+8veu874Pe0ZzQxtu7lPR89bMKB4CyQMCr1ETvg8590wGgvBBwAAAsRMABALAQAQcAwEIEHAAACxFwAAAsRMABALAQAa8yk72bGndjA4DyQMCrSCF3YRvF3dgAoLwQ8CpSyF3YRnE3NgAoLwS8ChV6VzXuxgYA5YOAAwBgIQIOAICFCHgVYQc5AFQOAl4lJrMDHQBQfqJhD4DSGN2BvumeRRPegX4mVvIAED5W4FViNLqT2UnOteAAUD4IeBUo1unzsdeC//roySJOCACYKAJeBX599GTBN3A509y2hvwq/P/882FW4gAQEgJe4T4+OJhffc9ta5j09+toTuiVlUuU9Hz9YHO3lj6xg4gDQAgIeAX7198e17c37pYkvbJyyaRX36PmdTRp0z2L8r/+4NNjOno8VZTvDQA4P44xxoQ9BAAAmBhW4AAAWIiAAwBgIQIOAICFCDgAABYi4AAAWGjS90I3xujEiRPFmAUVrrGxUY7jhD0GAFSESQe8r69PbW1txZgFFe7o0aOaOXNm2GMAQEWYdMDj8dy9tQ8cOKBp06ZNeiDbHT9+XJ2dnRyPMUaPyejfFQDA5E064KOnRKdNm0awxuB4nI3T5wBQPGxiAwDAQgQcAAALTTrgNTU1WrdunWpqaooxj/U4HmfjmABA8fEwEwAALMQpdAAALETAAQCwEAEHAMBCBBwAAAtNKuCe52nt2rW6/PLLVV9fr1mzZumOO+7QoUOHijWfdZ577jldfPHFqq2t1eLFi7V3796wRwrF448/roULF6qxsVFtbW1asWKF9u3bF/ZYAFAxJhXw4eFh9fT06JFHHlFPT49effVV7du3T8uXLy/WfFbZunWr1qxZo3Xr1qmnp0fz58/XDTfcoKNHj4Y9Wsnt3LlTq1at0p49e/TWW2/J8zxdf/31GhoaCns0AKgIRb+M7IMPPtCiRYvU29uriy66qJjfuuwtXrxYCxcu1Pr16yVJQRCos7NT9913nx5++OGQpwvX559/rra2Nu3cuVNXX3112OMAgPWK/h744OCgHMdRc3Nzsb91WctkMuru7tayZcvyr0UiES1btky7d+8OcbLyMDg4KEmaPn16yJMAQGUoasBTqZTWrl2r2267reoe5NHX1yff99Xe3n7a6+3t7Tp8+HBIU5WHIAj0wAMPaOnSpZo3b17Y4wBARZhQwLds2aKGhob8x65du/K/53mebr31VhljtGHDhqIPCnutWrVKH3/8sV5++eWwRwGAijGhx4kuX75cixcvzv+6o6ND0ql49/b2aseOHVW3+pak1tZWua6rI0eOnPb6kSNHdMEFF4Q0VfhWr16tbdu26d1339Xs2bPDHgcAKsaEVuCNjY2aO3du/iORSOTjvX//fr399tuaMWPGVM1a1uLxuLq6urR9+/b8a0EQaPv27VqyZEmIk4XDGKPVq1frtdde044dOzRnzpywRwKAijKhFfiZPM/TLbfcop6eHm3btk2+7+ff750+fbri8XhRhrTFmjVrdOedd+rKK6/UokWL9Mwzz2hoaEh333132KOV3KpVq/Tiiy/q9ddfV2NjY/7vRVNTkxKJRMjTAYD9JnUZ2WeffTbuyuqdd97RNddcU+i3ttb69ev15JNP6vDhw1qwYIGeffbZ0952qBaO45zz9eeff1533XVXaYcBgArE40QBYAr9ceTbYY+AMvdW8EpBX8e90AEAsBABBwDAQgQcAAALTWoXOgCMyxjJG859HquTxtnYCKAwrMABTA1vWPrxrNzHaMgBFA0BBwDAQgQcAAALEXAAACxEwAEAsBABnyRjjJ5++mnNmTNHdXV1WrFihQYHB8MeCwBQ4Qj4JD300EPasGGDNm3apF27dqm7u1uPPvpo2GMBACocAZ+E999/X08//bS2bt2qq6++Wl1dXbr33nv1s5/9LOzRAAAVjoBPwlNPPaXrrrtOV1xxRf619vZ29fX1hTgVUB52/Zp/D4CpRMALlE6n9eabb+qmm2467fVUKqWmpqaQpgLKw8GBpL7/k1+EPQZQ0Qh4gXp6epRMJvXggw+qoaEh//HDH/5Ql1xyiSTppptuUktLi2655ZaQpwVKq38oE/YIQMUj4AX65JNPVF9fr1/96lf68MMP8x8XX3yxli5dKkm6//779ZOf/CTkSQEAlYiAF+j48eNqbW3V3Llz8x+xWEz79+/XzTffLEm65ppr1NjYGPKkAIBKRMAL1NraqsHBQRlj8q899thjuvHGG3XppZeGOBkAoBrwONECXXvttUqlUnriiSf03e9+V1u2bNEbb7yhvXv3hj0aAKAKsAIvUHt7u1544QVt2LBBl112mfbs2aP33ntPnZ2dYY8GAKgCrMAn4Tvf+Y6+853vhD0GAKAKEfAptGzZMn300UcaGhrS7Nmz9corr2jJkiVhjwVMuS+4jAyYcgR8Cr399tthjwCU3MGBpFZu7g57DKDi8R44gKLqH8oo6fl6YNklYY8CVDRW4ACmxPS6WNgjlB/HOfX5mEtQgUKwAgcwJVrq4/nPDw0mQ5wkZI5z6kOSnMjZr4/9feA8sQIHMCUubErkP+8f8jRrZojDhMiJnuNMhHFkgtwK3Inkwm0CIxm/lKPBcgQcAKZS5Bwr6yAiJxKc/brjcGod542AA8AUijTUj3ziyom6kiST9aXAlwIjmUAKjIzvy6TTIyvxgJDjdyLgADCVWppy/4zH5NfFJceRk/bkpDzJGDleVsr6kufl/lw2K/m+jD9yOp2QYxwEHACmUjz3HniQiMlPxKRI7qx6xHEk3+SC7jhyJMl1c/F2IpIzshIHxkHAAWAKDX8ltwL3GiJKN0dkIpKbkmJJo0jWqKbfkzvkKZLyFDFGSqflZDwFydSpzW4Bm9twNgIOAFPoxEVRGUe6sPmYljf2qMMM6FDQrFezV+g36Rmq/21Etf1RRYfiqvF8OcOujJuWk81KQSBlszKGzW04GwEHgCmUaXS0wv2lfhT/qUxaciQZSffq/+q/1y3XPyUWKJt05GQjMjFXjuvKiURkXDcXcGAcBBwAplDr7x/Wjw7+VK5OX0EbST/2f6pfdF6kvtgMxRMRxYYTcl1XkYgjJ5XOfYUxks8pdJyNO7EBwBS6L7ort+w+gzPyf76V+IVSrVJqhqPMtJj8aXGZuhopGpXjRuRwhzaMg4ADKCoeJXq6C9Mnxn3/2jFGnf6ATMTIOJLhv8iYAP66ACia0UeJJmKuWup5mIkk9SXqxr3PuXGk38SaFERz8XZ8yckaOdng1A1eRuPPShxnIOAAimb0UaIbb+/SrDH3Qq9m73Z+RY4xOnMNbpS71PuVlstz/yV2JCcwcvyRu7CZMfEGzoGAAyi6GWOeRFbtPgi+oid+/w9lJAXOqQ8j6dk/uEr/lmpX9KSj6LAUyQS51XfWz21cC3KrcOBc2IUOAFOo9/+16W/Upn/6yld1R/oX6vCOq9dt0Uv1X9e/nZyput6o6g4b1ZwIFDueUWQ4IyeVUeBlZbJZ7saGcRFwAJhCsX5XcqTe2gv05w1/ItUaKRORczyiaNJRzUAu3tGhQJGMn1t9B8HICtzkHzsKnImAAyia8XagHxuu3p3p9Ydym8/8GkfZuoiMI0W83IebkRJ9vmr6s3JTvpzhtJx0JvdUMjPmXui8F45zIOAAiuL0HehxSadO/d7/0i/1xoNt6miuvo1trR8lJUlBLCK/JpLbrOYbRTyjiB8oOpiWk8zI8bIyQ8MyGU8mk8k91CTggSYYHwEHUBSjO9A33bMoF+rMUP73kp6v/qFMVQY8OpgLuHFdRWtyzwN3siOb1YIgF++Mlzt1nvFy73sTb5wHAg6gqNiBfrrI4MgPMpGITHQk4H4u3jJGSmdGoh3IpNNSEMj4Y8LN6XOMg4ADwBQKjg3kPx+9LWpgzKmAB0Hu/e7AyIze85yVN84DAQeAqeR5+U9NZOTWG6PRlnLXeZuAU+aYMAIOAFPIZLPnfn3s5WHsNkcBCDgATKHxAg5MFrdSBQDAQgQcAAALEXAAACxEwAEAsBABBwDAQgQcAAALEXAAACxEwAEAsBABB1AS4z0rHEBhCDiAKZeIuVq5uVsHB5JhjwJUDAIOYMr91W1fzz8THEBxEHAAU256Hc8IB4qNgAMAYCECDgCAhQg4AAAWIuAAioLLxIDSIuAAJu3gQFIrN3crEXPVUs+GNaAUomEPAMB+/UMZJT1fm+5ZpI7mRNjjAFWBFTiAopnB6hsoGQIOAICFCDgAABYi4AAAWIiAAwBgIQIOAICFCDgAABYi4AAAWIiAAwBgIQIOAICFCDgAABYi4AAAWIiAAwBgIQIOoGR4ZjhQPAQcwJRrqY8pEXO1cnO3Dg4kwx4HqAgEHMCUm9WU0Mbbu5T0fPWzCgeKgoADKAmeFQ4UFwEHAMBCBBwAAAsRcACTxu5yoPQIOIBJOTiQ1MrN3UrEXLXwPjdQMtGwBwBgt/6hjJKer033LFJHcyLscYCqwQocQFGwyxwoLQIOAICFCDgAABYi4AAAWIiAAwBgIQIOAICFCDgAABYi4AAAWIiAAwBgIQIOAICFCDgAABYi4AAAWIiAAwBgIQIOAICFCDiAkvpiKBP2CEBFIOAASqKlPq5EzNXKzd06OJAMexzAegQcQEl0NCe08fYuJT1f/azCgUkj4ABKZkZ9POwRgIpBwAEAsBABB1BybGQDJo+AA5iUicSYjWxA8RBwAAU7OJDUys3dSsRctZzH+9tsZAOKJxr2AADs1T+UUdLztemeRepoTpzX17CRDSgOVuAAJo0oA6VHwAEAsBABB1AwdpMD4SHgAAoy0Q1sAIqLTWwAClLIBjYAxcMKHMCksIENCAcBBwDAQgQcAAALEXAAoWAHOzA5BBxASXE/dKA4CDiAkuJ+6EBxEHAAJcfOdWDyCDgAABYi4AAKwiY0IFwEHMCEcRtVIHzcShXAhHEbVSB8rMABFIzNaEB4CDgAABYi4AAmrFgb2NgIBxSOgAOYkGJsYONubMDkEXAAEzK6gW3j7V0Fb2DjbmzA5BFwAAWZ7AY2NsABk0PAAYSK98GBwhBwAOft4EBSvz56sijfa+z74Ds/+Zz3woEJ4kYuAH6no8dT+tfDJ7Ryc7eSnl+UO7B1NCf0ysol+vbG3brz7/YqEXO18fYuzaiPq62xRm3Taos0PVCZHGOMCXsIAAAwMZxCBwDAQgQcAAALEXAAACxEwAEAsBABBwDAQlxGhopijNGJEyfCHgMWaGxslOM4YY8BFIyAo6L09fWpra0t7DFggaNHj2rmzJlhjwEUjICjosTjuZuLHDhwQNOmTQt5mvAdP35cnZ2dHI8xRo/J6N8VwFYEHBVl9JTotGnTCNYYHI+zcfoctmMTGwAAFiLgAABYiICjotTU1GjdunWqqakJe5SywPE4G8cElYKHmQAAYCFW4AAAWIiAAwBgIQIOAICFCDgAABYi4KgYnudp7dq1uvzyy1VfX69Zs2bpjjvu0KFDh8IeLTTPPfecLr74YtXW1mrx4sXau3dv2COF4vHHH9fChQvV2NiotrY2rVixQvv27Qt7LGBSCDgqxvDwsHp6evTII4+op6dHr776qvbt26fly5eHPVootm7dqjVr1mjdunXq6enR/PnzdcMNN+jo0aNhj1ZyO3fu1KpVq7Rnzx699dZb8jxP119/vYaGhsIeDSgYl5Ghon3wwQdatGiRent7ddFFF4U9TkktXrxYCxcu1Pr16yVJQRCos7NT9913nx5++OGQpwvX559/rra2Nu3cuVNXX3112OMABWEFjoo2ODgox3HU3Nwc9igllclk1N3drWXLluVfi0QiWrZsmXbv3h3iZOVhcHBQkjR9+vSQJwEKR8BRsVKplNauXavbbrut6h7k0dfXJ9/31d7eftrr7e3tOnz4cEhTlYcgCPTAAw9o6dKlmjdvXtjjAAUj4LDWli1b1NDQkP/YtWtX/vc8z9Ott94qY4w2bNgQ4pQoN6tWrdLHH3+sl19+OexRgEnhcaKw1vLly7V48eL8rzs6OiSdindvb6927NhRdatvSWptbZXrujpy5Mhprx85ckQXXHBBSFOFb/Xq1dq2bZveffddzZ49O+xxgElhBQ5rNTY2au7cufmPRCKRj/f+/fv19ttva8aMGWGPGYp4PK6uri5t3749/1oQBNq+fbuWLFkS4mThMMZo9erVeu2117Rjxw7NmTMn7JGASWMFjorheZ5uueUW9fT0aNu2bfJ9P/9+7/Tp0xWPx0OesLTWrFmjO++8U1deeaUWLVqkZ555RkNDQ7r77rvDHq3kVq1apRdffFGvv/66Ghsb838vmpqalEgkQp4OKAyXkaFifPbZZ+OurN555x1dc801pR2oDKxfv15PPvmkDh8+rAULFujZZ5897W2HauE4zjlff/7553XXXXeVdhigSAg4AAAW4j1wAAAsRMABALAQAQcAwEIEHAAAC3EZGQCUgjGSN5z7PFYnjbMzHjhfrMABoBS8YenHs3IfoyEHJoGAAwBgIQIOAICFCDgAABYi4EAFMcbo6aef1pw5c1RXV6cVK1ZocHAw7LEATAECDlSQhx56SBs2bNCmTZu0a9cudXd369FHHw17LABTgIADFeL999/X008/ra1bt+rqq69WV1eX7r33Xv3sZz8LezQAU4CAAxXiqaee0nXXXacrrrgi/1p7e7v6+vpCnArn8gd/9o/6+CBvbWByCDhQAdLptN58803ddNNNp72eSqXU1NQU0lT4MgQck0XAgQrQ09OjZDKpBx98UA0NDfmPH/7wh7rkkkt04MABXXPNNbr00kv1ta99Ta+88krYIwOYJG6lClSATz75RPX19frwww9Pe/2b3/ymli5dqmg0qmeeeUYLFizQ4cOH1dXVpRtvvFH19fXhDAxg0gg4UAGOHz+u1tZWzZ07N/9ab2+v9u/fr5tvvlkXXnihLrzwQknSBRdcoNbWVh07doyAAxbjFDpQAVpbWzU4OChjTP61xx57TDfeeKMuvfTS0/5sd3e3fN9XZ2dnqccEUESswIEKcO211yqVSumJJ57Qd7/7XW3ZskVvvPGG9u7de9qfO3bsmO644w797d/+bUiTAigWVuBABWhvb9cLL7ygDRs26LLLLtOePXv03nvvnbbKTqfTWrFihR5++GFdddVVIU4LSfpiKBP2CLCcY8aecwNQkYwx+t73vqevfvWr3JktLJmh3KNEJf1B6u+UVK1+/vC16mhOhDwYbMUKHKgCP//5z7V161b9wz/8gxYsWKAFCxboV7/6VdhjVb1+VuGYBN4DB6rAN77xDQVBEPYYAIqIFTgAlNh//cP/EPYIqAAEHABK7MKm2rBHQAUg4AAAWIiAA0AJHBpM5j9vqY+HOAkqBQEHgBLoH/LGfJ7bfc4TyTAZBBwASuxffntckvTxIQKOwhFwACixSy+cJkmaN4tntaNwBBwASuzSkXDP6yDgKBwBBwDAQgQcAAALEXAAKIFjw9z3HMVFwAFgih0cSOr+l34Z9hioMAQcAKZY/1BGSc8PewxUGAIOAICFCDgAABYi4AAAWIiAAwBgIQIOACH5079+jweaoGAEHABCRMBRKAIOAICFCDgAlFhLfWzM5/EQJ4HNCDgAlNispoT+5vYuSVJHcyLkaWArAg4AISDcmCwCDgCAhQg4AAAWIuAAAFiIgAMAYCECDgCAhQg4AAAWIuAAMIU+PjioP/3r98b9/S+GMiWcBpWEgAPAFBrvXuct9XElYq5Wbu7WwYFkiadCJSDgABCCjuaENt7epaTnq59VOApAwAEgJDO4DzomgYADAGAhAg4AgIUIOABMIR4XiqlCwAFgCo0+deyvb/t6yJOg0hBwACiBC5t4fCiKi4ADQMi4mQsKQcABICTczAWTQcABICTczAWTQcABYAr9rtPj3MwFhSLgADCFVm7uViLmqqU+FvYoqDDRsAcAgEqW9HxtumeRZrELHUXGChwAphinyTEVCDgAABYi4AAAWIiAAwBgIQIOAICFCDgAABYi4AAAWIiAAwBgIQIOAICFCDgAABYi4AAAWIiAAwBgIQIOAGXgdz12FDgTAQeAELXUx5WIuVq5uVsHB5JhjwOLEHAACFFHc0Ibb+9S0vPVzyocE0DAASBkPG4UhSDgAABYiIADAGAhAg4AgIUIOAAAFiLgAABYiIADAGAhAg4AgIUIOAAAFiLgAABYiIADAGAhAg4AgIUIOAAAFiLgAABYiIADAGAhAg4AZeILngeOCSDgABCylvq4EjFXKzd36+BAMuxxYAkCDgAh62hOaOPtXUp6vvpZheM8EXAAKAMz6uNhjwDLEHAAACxEwAEAsBABBwDAQgQcAAALEXAAACxEwAGgjHAzF5wvAg4AZYCbuWCiCDgAlAFu5oKJIuAAUCa4mQsmgoADAGAhAg4AgIUIOAAAFiLgAABYiIADAGAhAg4AgIUIOAAAFiLgAABYiIADAGAhAg4AgIUIOAAAFiLgAFBmeKQozgcBB4AyMfpI0Tv/bq8++OxY2OOgzBFwACgTHc0JrfvPl0qSvr1xN88Fx5ci4ABQpnguOL4MAQeAMjKvoynsEWAJAg4AZWReR5O23feNsMeABQg4AAAWIuAAAFiIgAMAYCECDgBlisvI8GUIOACUmdFw/2BzNxHHuAg4AJSZjuZE/nOuBcd4CDgAlBkuJcP5IOAAUMYODiQ5jY5zIuAAUMZ+sLlby/7HTiKOsxBwAChDo08mk6Sk5/NeOM4SDXsAAMDZOpoTevvBP9THBwf1g83d+vXRk/nfa2usUdu02hCnQzlwjDEm7CEAAMDEcAodAAALEXAAACxEwAEAsBABBwDAQgQcAAALcRkZAEwRY4xOnDgR9hiwQGNjoxzHmdDXEHAAmCJ9fX1qa2sLewxY4OjRo5o5c+aEvoaAA8AUicfjkqQDBw5o2rRpIU8TvuPHj6uzs5PjMcboMRn9uzIRBBwApsjoKdFp06YRrDE4Hmeb6OlziU1sAABYiYADAGAhAg4AU6Smpkbr1q1TTU1N2KOUBY7H2SZzTHiYCQAAFmIFDgCAhQg4AAAWIuAAAFiIgAMAYCECDgCAhQg4AEwBz/O0du1aXX755aqvr9esWbN0xx136NChQ2GPFprnnntOF198sWpra7V48WLt3bs37JFC8fjjj2vhwoVqbGxUW1ubVqxYoX379k34+xBwAJgCw8PD6unp0SOPPKKenh69+uqr2rdvn5YvXx72aKHYunWr1qxZo3Xr1qmnp0fz58/XDTfcoKNHj4Y9Wsnt3LlTq1at0p49e/TWW2/J8zxdf/31GhoamtD34TpwACiRDz74QIsWLVJvb68uuuiisMcpqcWLF2vhwoVav369JCkIAnV2duq+++7Tww8/HPJ04fr888/V1tamnTt36uqrrz7vr2MFDgAlMjg4KMdx1NzcHPYoJZXJZNTd3a1ly5blX4tEIlq2bJl2794d4mTlYXBwUJI0ffr0CX0dAQeAEkilUlq7dq1uu+22qnsSV19fn3zfV3t7+2mvt7e36/DhwyFNVR6CINADDzygpUuXat68eRP6WgIOAEWwZcsWNTQ05D927dqV/z3P83TrrbfKGKMNGzaEOCXKzapVq/Txxx/r5ZdfnvDX8jxwACiC5cuXa/Hixflfd3R0SDoV797eXu3YsaPqVt+S1NraKtd1deTIkdNeP3LkiC644IKQpgrf6tWrtW3bNr377ruaPXv2hL+eFTgAFEFjY6Pmzp2b/0gkEvl479+/X2+//bZmzJgR9pihiMfj6urq0vbt2/OvBUGg7du3a8mSJSFOFg5jjFavXq3XXntNO3bs0Jw5cwr6PqzAAWAKeJ6nW265RT09Pdq2bZt838+/3zt9+nTF4/GQJyytNWvW6M4779SVV16pRYsW6ZlnntHQ0JDuvvvusEcruVWrVunFF1/U66+/rsbGxvzfi6amJiUSifP+PlxGBgBT4LPPPht3ZfXOO+/ommuuKe1AZWD9+vV68skndfjwYS1YsEDPPvvsaW87VAvHcc75+vPPP6+77rrr/L8PAQcAwD68Bw4AgIUIOAAAFiLgAABYiIADAGAhAg4AgIUIOAAAFiLgAABYiIADAGAhAg4AgIUIOABgShhj9PTTT2vOnDmqq6vTihUrNDg4GPZYFYOAAwCmxEMPPaQNGzZo06ZN2rVrl7q7u/Xoo4+GPVbF4F7oAICie//997VkyRL94he/0BVXXCFJ+tGPfqQtW7Zo3759IU9XGViBAwCK7qmnntJ1112Xj7cktbe3q6+vL8SpKgsBBwAUVTqd1ptvvqmbbrrptNdTqZSamppCmqrycAodAFBUu3fv1lVXXaXa2lq5rpt/3fM8/dEf/ZFefvllLVu2TNlsVtlsVvfff7/uvffeECe2UzTsAQAAleWTTz5RfX29Pvzww9Ne/+Y3v6mlS5eqsbFR7777rurq6jQ0NKR58+bpW9/6lmbMmBHOwJYi4ACAojp+/LhaW1s1d+7c/Gu9vb3av3+/br75Zrmuq7q6Okm50+3GGHEyeOJ4DxwAUFStra0aHBw8LcqPPfaYbrzxRl166aWSpIGBAc2fP1+zZ8/WQw89pNbW1rDGtRYBBwAU1bXXXqtUKqUnnnhCn376qf7iL/5Cb7zxhjZs2JD/M83Nzfroo4/06aef6sUXX9SRI0dCnNhOBBwAUFTt7e164YUXtGHDBl122WXas2eP3nvvPXV2dp7zz86fP1+7du0KYVK7sQsdAFBSR44cUV1dnRobGzU4OKilS5fqpZde0uWXXx72aFZhExsAoKR6e3v1/e9/P7957b777iPeBWAFDgCAhXgPHAAACxFwAAAsRMABALAQAQcAwEIEHAAACxFwAAAsRMABALAQAQcAwEIEHAAACxFwAAAsRMABALDQ/wc/ipEkQKIhtgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 600x600 with 9 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "samples = posterior.sample((10000,), x=x_obs)\n",
    "_ = pairplot(samples, points=theta_true, limits=[[-2, 2], [-2, 2], [-2, 2]], figsize=(6, 6), labels=[r\"$\\theta_1$\", r\"$\\theta_2$\", r\"$\\theta_3$\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The log-probability should ideally indicate that the true parameters, given the corresponding observation, are more likely than a different set of randomly chosen parameters from the prior distribution. \n",
    "\n",
    "Relative to the obtained log-probabilities, we can investigate the range of log-probabilities of the parameters sampled from the posterior."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first sample an alternative parameter set from the prior\n",
    "theta_diff = prior.sample((1,))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "high for true theta : tensor([3.7909])\n",
      "low for different theta : tensor([-326.5020])\n",
      "range of posterior samples: min: tensor(-8.1930)  max : tensor(4.0474)\n"
     ]
    }
   ],
   "source": [
    "log_probability_true_theta = posterior.log_prob(theta_true, x=x_obs)\n",
    "log_probability_diff_theta = posterior.log_prob(theta_diff, x=x_obs)\n",
    "log_probability_samples = posterior.log_prob(samples, x=x_obs)\n",
    "\n",
    "print( r'high for true theta :', log_probability_true_theta)\n",
    "print( r'low for different theta :', log_probability_diff_theta)\n",
    "print( r'range of posterior samples: min:', torch.min(log_probability_samples),' max :', torch.max(log_probability_samples))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Next steps\n",
    "\n",
    "To learn more about the capabilities of `sbi`, you can head over to the tutorial on [inferring parameters for multiple observations ](https://sbi-dev.github.io/sbi/tutorial/01_gaussian_amortized/) which introduces the concept of amortization. \n",
    "\n",
    "Alternatively, for an example with an __actual__ simulator, you can read our [example for a scientific simulator from neuroscience](https://sbi-dev.github.io/sbi/examples/00_HH_simulator/)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}