CBayesTutorial_Example2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# An Introduction to \"Consistent Bayesian\" Inversion \n",
    "\n",
    "## Example 2: A higher dimensional problem and template for working with data sets\n",
    "\n",
    "This tutorial shows how to perform consistent Bayesian inversion directly on any data set. \n",
    "\n",
    "In particular, we use a data set for the 100-dimensional example shown in  \n",
    "\n",
    "_[***Combining Push-Forward Measures and Bayes' Rule to Construct Consistent Solutions to Stochastic Inverse Problems\n",
    "***, T. Butler, J. Jakeman, T. Wildey, SIAM J. Sci. Comput., 40(2), A984–A1011](https://epubs.siam.org/doi/abs/10.1137/16M1087229)_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Libraries\n",
    "\n",
    " * We make use of `numpy` to handle the sample sets defined as arrays, and we also use the `random` subpackage within `numpy` to generate random samples of various distributions.\n",
    " \n",
    " \n",
    " * We use a standard kernel density estimator in `scipy.stats` to approximate the densities on sample sets.\n",
    " \n",
    " \n",
    " * We use `matplotlib.pyplot` to visualize results.\n",
    " \n",
    " \n",
    " * We use `scipy.io` to load a data set. The `scipy.io` can read/write many data formats and is an overall useful way to connect Python scripts to computational codes written in almost any language. Read more about `scip.io` here https://docs.scipy.org/doc/scipy/reference/io.html."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The libraries we will use\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import chi2, gaussian_kde as GKDE\n",
    "from scipy.stats import norm \n",
    "import scipy.io as sio"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the data set containing both QoI and parameter samples\n",
    "\n",
    "Loading the data set `elliptic_kde100_10k.mat` creates a dictionary where\n",
    "\n",
    "\n",
    " * QoI samples are \"tagged\" with the key `qq` and are a $10000\\times 3$ array, but we only consider the first QoI samples defined by the first column of this array.\n",
    "\n",
    "\n",
    " * Parameter samples are \"tagged\" with the key `pp` and are a $10000\\times 100$ array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_set = sio.loadmat('elliptic_kde100_10K.mat')\n",
    "qvals = data_set['qq'] # QoI samples\n",
    "qvals = qvals[:,0] # Only using first QoI here\n",
    "lam = data_set['pp'] # parameter samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What are the observations?\n",
    "\n",
    "While it is possible/reasonable for the observed data or a description of the observed density to be included with the data set, we instead consider the situation where observations are collected independently/separately from the experiments/simulations. \n",
    "\n",
    "For the sake of argument, we assume that the observations lead to a $N(0.7,0.01^2)$ density being defined on the QoI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define an observed density\n",
    "obs_vals = norm.pdf(qvals, loc=0.7, scale=1.0e-2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What do we need to do with the data set to construct the consistent posterior?\n",
    "\n",
    "The fundamental assumption is that the data set we loaded comes from some ***prior*** knowledge, i.e., the parameter samples come from a prior and the QoI data are associated with these prior samples. \n",
    "\n",
    "All we need to do with these samples is perform density estimation on the QoI samples.\n",
    "\n",
    "***We do not need to know anything about this prior to define the re-weighting of these samples, which subsequently means we can use this existing sample set to generate i.i.d. samples from the posterior.*** "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Compute the pushforward of the prior\n",
    "q_kde = GKDE( qvals, 'silverman' )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now we perform rejection sampling (see previous lecture for a review)\n",
    "\n",
    "Notice how ***fast*** this is because all the computations are taking place in a 1-dimensional data space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rejection_sampling(r):\n",
    "    # Perform accept/reject sampling on a set of proposal samples using\n",
    "    # the weights r associated with the set of samples and return\n",
    "    # the indices idx of the proposal sample set that are accepted.\n",
    "    N = r.size # size of proposal sample set\n",
    "    check = np.random.uniform(low=0,high=1,size=N) # create random uniform weights to check r against\n",
    "    M = np.max(r)\n",
    "    new_r = r/M # normalize weights \n",
    "    idx = np.where(new_r>=check)[0] # rejection criterion\n",
    "    return idx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "27.58\n"
     ]
    }
   ],
   "source": [
    "# Compute the posterior\n",
    "r = np.divide(obs_vals,q_kde(qvals))\n",
    "    \n",
    "samples_to_keep = rejection_sampling(r)\n",
    "\n",
    "post_q = qvals[samples_to_keep]\n",
    "post_lam = lam[:,samples_to_keep]\n",
    "\n",
    "accept_rate = samples_to_keep.size/lam.shape[0]\n",
    "print(accept_rate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Don't forget to do some diagnostics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7001278189914161\n",
      "0.010121988757525296\n",
      "0.9967766306208634\n",
      "0.7371862445858198\n"
     ]
    }
   ],
   "source": [
    "print(np.mean(post_q))\n",
    "print(np.sqrt(np.var(post_q)))\n",
    "print(np.mean(r))\n",
    "print(np.mean(r*np.log(r)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing results\n",
    "\n",
    " * We almost always start by checking the push-forwards against the observed density.  \n",
    " \n",
    " \n",
    " * If the data space is higher-dimensional, we may only visualize the marginals of the densities, or we may instead opt to visualize certain statistics of the distribution or the marginals as we show below for the 100-dimensional parameter space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEKCAYAAAARnO4WAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8FdX5+PHPkwUSCBDWAAlrWEQQIqsbqxsqRUDrhltFqVot/hSK1lqxdWurpdaqrXWvVKkoatXaL2VxRYUoKMi+J2xhCUgSsp7fHzNJ7kxuyM29c5Ob3Of9et1X7jmZOXMmF55Mzpx5jhhjUEop1fjF1HcHlFJK1Q0N+EopFSU04CulVJTQgK+UUlFCA75SSkUJDfhKKRUlNOArpVSU0ICvlFJRQgO+UkpFibi6PFi7du1M9+7d6/KQSinV4GVmZh4wxrQPtZ06Dfjdu3dn5cqVdXlIpZRq8ERkhxft6JCOUkpFCQ34SikVJTTgK6VUlKjTMXylVN0rLi4mKyuL48eP13dXVA0SEhJIS0sjPj4+LO0HHPBFJBZYCWQbYyaISA/gdaAN8DVwjTGmKCy9VEoFLSsrixYtWtC9e3dEpL67o6phjOHgwYNkZWXRo0ePsByjNkM6M4B1PuXfAXONMb2Bw8A0LzumlPLG8ePHadu2rQb7CCcitG3bNqx/iQUU8EUkDbgIeM4uCzAOWGBv8jIwKRwdVBFg3jyYMgWefrq+e6KCpMG+YQj35xTokM6fgF8ALexyWyDXGFNil7OAVI/7piLB4sVw9dXW+4ULITkZrrqqfvuklApKjVf4IjIB2G+MyfSt9rOp38VxRWS6iKwUkZU5OTlBdlPVmz//2Vl+4on66Ydq0LKysrj44ovp3bs36enpzJgxg6KiIl566SVuu+22+u5eFUlJSfXdhbAIZEjnTGCiiGzHukk7DuuKP1lEyv9CSAN2+9vZGPOsMWaoMWZo+/YhPxms6tLevfD++866r76CNWvqpz+qQTLGMGXKFCZNmsSmTZvYuHEjx44d49577w3L8UpKSmreKErVGPCNMfcYY9KMMd2BK4AlxpipwFLgUnuz64B3wtZLVT9efRVKS6vWv/hi3fdFhU4kvK9qLFmyhISEBH7yk58AEBsby9y5c3nhhRfIz89n165djB8/nr59+/LAAw8AkJeXx0UXXcSgQYMYMGAA8+fPByAzM5PRo0czZMgQzj//fPbs2QPAmDFj+OUvf8no0aN56KGH6N69O2VlZQDk5+fTpUsXiouL2bJlC+PHj2fIkCGMHDmS9evXA7Bt2zZOP/10hg0bxn333Re2j6C+hfLg1WzgThHZjDWm/7w3XVIRwRh44QX/3/vHP6C4uG77oxqstWvXMmTIEEddy5Yt6dq1KyUlJXz11VfMmzePVatW8cYbb7By5Uo+/PBDOnfuzOrVq1mzZg3jx4+nuLiY22+/nQULFpCZmckNN9zg+CshNzeXjz76iPvvv59Bgwbx0UcfAfDvf/+b888/n/j4eKZPn86TTz5JZmYmjz32GLfeeisAM2bM4JZbbmHFihV07Nix7n44daxWD14ZY5YBy+z3W4Hh3ndJRYSvvoJ16/x/LyfHGuqZpBOzVM2MMX5nn5TXn3vuubRt2xaAKVOm8Omnn3LhhRcyc+ZMZs+ezYQJExg5ciRr1qxhzZo1nHvuuQCUlpbSqVOnivYuv/xyx/v58+czduxYXn/9dW699VaOHTvG559/zo9//OOK7QoLCwH47LPPePPNNwG45pprmD17tvc/iAigT9oq/6q7ui/34osa8FVA+vfvXxFMyx09epRdu3YRGxtb5ZeBiNCnTx8yMzP54IMPuOeeezjvvPOYPHky/fv3Z/ny5X6P07x584r3EydO5J577uHQoUNkZmYybtw48vLySE5OZtWqVX73j4apq5pLR1WVnw+vv+6s+9WvnOX337du6qqGw5jwvqpx9tlnk5+fzyuvvAJYV+Z33XUX119/Pc2aNWPRokUcOnSIgoIC3n77bc4880x2795Ns2bNuPrqq5k5cyZff/01ffv2JScnpyLgFxcXs3btWr/HTEpKYvjw4cyYMYMJEyYQGxtLy5Yt6dGjB2+88Yb94zCsXr0agDPPPJPX7X/z8+bN8+xHHmk04Kuq3noLjh6tLKekwK9/Df37V9aVllpj+UrVQERYuHAhb7zxBr1796ZPnz4kJCTw8MMPA3DWWWdxzTXXkJGRwSWXXMLQoUP57rvvGD58OBkZGTz00EP86le/okmTJixYsIDZs2czaNAgMjIy+Pzzz6s97uWXX86rr77qGOqZN28ezz//PIMGDaJ///6884411+SJJ57gqaeeYtiwYRw5ciS8P5B6JOYEv5m9NnToUKMLoDQAF1wAH34IQClw5+DBvH3gAKPbt+eZzEwq/nAeOBDsKyQVudatW0e/fv3quxsqQP4+LxHJNMYMDbVtHcNXVX3zTcXbPwF//vprAP6xcye9gF+Xf3PNGigqgiZN6rqHSqkg6JCOcvrhB9i3D7Ay4j1kVzdp0oTf/va3zPadslZWBjs8WXlNKVUHNOArp61bK94+ghX0AYqKinjsscf4wb0I/ZYtddUzpVSINOArp82bAdgJuLLoMHv2bNqddJLf7ZVSkU8DvnKyA/ivgUKf6tTUVGbMmAG9evndXikV+TTgK6ctW/gWeMVV/cADD9CsWbOqAV+HdJRqMDTgK6fNm3kYZ67rfv36cd1111mF9PSK+jygcOPGuuydaoBiY2PJyMhgwIAB/PjHPyY/P99RX/7avn17wG3OmjWL/v37M2vWrKD7dcYZZwS9b0Ol0zKV0+bNfOSqeuSRR4iLs/6pLNq5k38DnwOrgHe2buWi0lKIja3jjqqGIjExsSKdwdSpU/nrX//KnXfe6aivrb/97W/k5OTQtGnTWu9bWlpKbGzsCR/aqm6fhk6v8FWl48fJ3bUL34QJ8fHxXHTRRRXld5Ys4UkgE+uhrM9KSiA7u447qoJRT9mRHUaOHMnmAO/7GGOYNWsWAwYM4JRTTqlIkTxx4kTy8vIYMWJERV25OXPmcM011zBu3Dh69+7N3//+dwCWLVvG2LFjueqqqzjllFOAykVOqjuOv30aOr3CV5W2bcOdH7N3794VV/dg/Rn81FNPVZQ/B+vGbdeuddJF1XCVlJTwn//8h/HjxwNQUFBARkYGAD169GDhwoWO7d966y1WrVrF6tWrOXDgAMOGDWPUqFG8++67JCUlVfvXwbfffssXX3xBXl4ep556asUFy1dffcWaNWvo0aNHQMc50T4NlQZ8VWnz5ioB/yTXNEz3uOdXQPGGDcSPGxfevqkGyzewjxw5kmnTpgHUOKTz6aefcuWVVxIbG0tKSgqjR49mxYoVTJw48YTHu/jii0lMTCQxMZGxY8fy1VdfkZyczPDhw/0G7uqO07Jly2r3aag04KtKW7ZUCfjunB7dunWjU1ISe44dA6AAWL18OUNvuaVu+qganGDH6oPN8+Uv3TI40ycHepzq9mmoAlnEPEFEvhKR1SKyVkQesOtfEpFtIrLKfmWEv7sqrPxc4ftJ4sQZrqv+zzWBWoNQT9mRgzZq1Cjmz59PaWkpOTk5fPzxxwwfXvOaS++88w7Hjx/n4MGDLFu2jGHDhoXlOA1RIDdtC4FxxphBQAYwXkROs783yxiTYb+Cu92uIseWLXQFelP5D8NflsUzTj/dUf58586wd01Fn8mTJzNw4EAGDRrEuHHj+P3vfx/Q8oPDhw/noosu4rTTTuO+++6jc+fOYTlOQ1Sr9Mgi0gz4FLjFfr1njFkQ6P6aHjnC9e5d8eRsIbB5wQJ6TZhQZerbFx98wOk+M3e6iLCztDTwqRqqTkVTeuQ5c+aQlJTEzJkz67srQQtneuSApmWKSKyIrAL2A4uMMV/a33pIRL4Vkbki4ndCrIhMF5GVIrIyJycn1P6qcCkpAZ8HX5oC/c8/3+8851PHjcO3dpcx7LJTKCulIldAAd8YU2qMyQDSgOEiMgC4BzgJGAa0Afyu+muMedYYM9QYM7R9+/YedVt5budOK+iXS0kBe56yW9OEBIa6bmZ9/v774eydUgGZM2dOg766D7daPXhljMkFlgHjjTF7jKUQeBFonHc5ooX7YRh3zhyXM1JTHeXPP/3U6x4ppTwWyCyd9iKSbL9PBM4B1otIJ7tOgEnAmnB2VIWZO+D75Mzx5wzf9W2Bz9fox69UpAvkCr8TsFREvgVWYI3hvwfME5HvgO+AdsCD4eumCrstWyjzLddwhX/6WWc5yt/s3UteXp73/VJKeabGB6+MMd8Cp/qp10crG5HSTZtoA3TBujHT7+uvub+kxJFWwVfK4MGkA+XJkUuNITMzs+KRdKVU5NHkaQqA7evXcxRYC7wJ/HXZsmqDPQDp6ZyL9WDGhcCNTZvSsmXLuuiqamAiNT3yiTz88MNB7XfjjTfy/fffe9wb79RqHn6odB5+hCor472EBH5UXFxRddZpp/HJ8uUn3IdmzaDQZ12sgwehTZswdlQFo77n4SclJXHMTsUxdepUhgwZwp133umor62WLVsGnR45EMH0LZgUyv72qfd5+KqR272bdT7BHqBfTelgY2Kq3tjV1a8aBBEJ6jVkyJCQj11f6ZGra2vPnj2MGjWq4i+QTz75hLvvvrsi4dvUqVMBePXVVxk+fDgZGRn89Kc/pbS0FLB+Mfz6179mxIgRLF++nDFjxlB+Ufvaa69xyimnMGDAAGbPrpy17t6nThlj6uw1ZMgQoyLQJ5+YG6xFripef/zjH2ve70c/cqZTmT8//H1Vtfb99987yrg+60BfgwcPDur4zZs3N8YYU1xcbCZOnGiefvppY4wxMTExZtCgQWbQoEFm0qRJVfZbsGCBOeecc0xJSYnZu3ev6dKli9m9e7ejTbf777/fDBw40OTn55ucnByTlpZmsrOzq23rscceMw8++KAxxpiSkhJz9OjRKu1///33ZsKECaaoqMgYY8wtt9xiXn755Yqf5Xyff/ejR482K1asMNnZ2aZLly5m//79pri42IwdO9YsXLjQ7z5u7s/L3mel8SAGa7ZMZV3hu6oCGgJwzcVnzx7PuqQaj0hIj1xdW8OGDeOGG26guLiYSZMmVfTT1+LFi8nMzKxIwlZQUECHDh0A6z7EJZdcUmWfFStWMGbMGMofNp06dSoff/wxkyZNqnafuqABX2GCDfjupFQa8JUfkZAeubq2Ro0axccff8z777/PNddcw6xZs7j22mur9OO6667jkUceqbJ/QkKC33H7E/W9un3qgo7hK/Zt2kSuT7lZfDxdunSpeUefgG+A3O3bg/5PqupOsMMBmZmZddpPL9MjV9fWjh076NChAzfddBPTpk3jazsnVHx8PMX2fa2zzz6bBQsWsH//fgAOHTrEjh07TtiHESNG8NFHH3HgwAFKS0t57bXXGD16dIg/kdDpFb6qchOtb6dOxMQEcC3QqRMXARuA3UDB/Pnk/OUvtGvXLhzdVFFm8uTJLF++nEGDBiEitU6PvHPnzor0yNW19fLLL/OHP/yB+Ph4kpKSeOWVVwCYPn06AwcOZPDgwcybN48HH3yQ8847j7KyMuLj43nqqafo1q1btX3o1KkTjzzyCGPHjsUYw4UXXsjFF1/s2c8mWDotU/H6gAFcuXZtRXnyGWfw1mef1bzjqlX0O/VU1vtUrV69moEDB3rfSRW0+p6WWZc0PfKJ6ZCOItuVtjo10AXJO3XCvbTE7t27vemUUspzOqSjyMrNdZTTasijU6F9ezqLONa3212LpyWV8tqcOXPquwsRTa/wo11+PllFRY6qtD59Ats3JobOrrz4uzdu9KpnykN6M71hCPfnpAE/2u3ZQ5arKi3QIR2gc+vWjvLubds86JTyUkJCAgcPHtSgH+GMMRw8eJCEhISwHUOHdKLd7t10wVq7MhtrLdu0tLSAd++ckgK7dlU2l+X+9aHqW1paGllZWegSo5EvISGhVv//aksDfrTbs4d/2W8NcPCCC2jdvXvAu3dOSwOfmVfZ+/Z52j0Vuvj4eHr06FHf3VARQAN+tPOZVSNAux49oBZPAXZ2BZLdhw971TOllMcCWeIwQUS+EpHVIrJWRB6w63uIyJcisklE5otIk/B3V3nOPY2yU6da7d7JdYN3b15eRSZBpVRkCeSmbSEwzhgzCGu9i/EichrwO2CuMaY3cBiYFr5uqrBx579x58epQUK3bvhmwC8zpuIRdKVUZKkx4NvZOctXAoi3XwYYByyw61/GWshcNTTuK/xaBnx9+EqphiOgaZkiEisiq7AmcyzCWso01xhTYm+SBaRWs+90EVkpIit1lkAECnFIh86dNeAr1UAEFPCNMaXGmAwgDRgO+EvM4XeSrzHmWWPMUGPM0PLc0Cpy3Lp1KxnABOBmYPWhQ7VroF0762lbH/q0rVKRqVazdIwxuSKyDDgNSBaROPsqPw0rYaJqSPLyWFtUxGpgtV11ieup2xrFxJCalES7H36gM9AZSImP97afSilPBDJLp72IJNvvE4FzgHXAUuBSe7PrgHfC1UkVJiE+ZVvutyedRA7WL43/AJP69/egc0oprwVyhd8JeFlEYrF+QfzLGPOeiHwPvC4iDwLfAM+HsZ8qDMzu3VUDfhBP+YmufKVUg1BjwDfGfAuc6qd+K9Z4vmqgDmzYgO8ATsu4OFq0aFH7htwBX2/aKhWRNHlaFMtav95RTgsm2IMGfKUaCA34USzLldkyrW3b4BpyT+XUIR2lIpIG/CiW5cpsmRrAeqF+ua7wTXa2puJVKgJpwI9i2a4UCGldugTXUOfO/BIYDfQGkpYtY8WKFaF2TynlMc2WGcWyXJktA17a0K1TJ1YBH5eXjdGnbZWKQHqFH8Wy8vIc5bR+/h6gDoC/p2115SulIo4G/GiVl0eWK41x2sknB9dWTAydk5IcVbs3bQq2Z0qpMNGAH6X8PnQV7Bg+0Nk1w2f3jh1Bt6WUCg8N+FHqyMaN+A7oJMbE0Nq1IHltdE5JcZR3Z2cH3ZZSKjw04Eepw1u2kA40tcupiYmIaxy+NlJSndmxcw4eDL5zSqmw0Fk6UapHaSmbsRcuB3InhbZ+TQfXwuf7jxwJqT2llPc04Ecre9qkAO2AdsHesLV1cE3p3J+XhzEmpL8alFLe0iGdaOVOfxDsU7a25t2708ynXFRWxhG9ylcqomjAj1Z79zrLtV3a0K1jRzq4qnQxc6Uiiwb8aOUO+CFe4dOxIymuKg34SkWWQFa86iIiS0VknYisFZEZdv0cEckWkVX268Lwd1d5xuuA3759lSv8fVnumf5KqfoUyE3bEuAuY8zXItICyBSRRfb35hpjHgtf91RYFBUx6eBBDJBiv2Y1bUrLUNqMi6NDQgIcP15RtX/r1tD6qZTyVCArXu0B9tjvfxCRdUDqifdSkczs28d/geM+dbPiQp+w1aFlS+T4cdoCHYDEwsKQ21RKeadWY/gi0h1rucMv7arbRORbEXlBRIJ/TFPVqaNbtjiCfaJIcEsbutw3cCBFQA6wFrh+6NCQ21RKeSfggC8iScCbwB3GmKPAM0A6kIH1F8Dj1ew3XURWisjKnJwcD7qsQrVvwwZHOSUhwZP58ompqc4/Gd33CZRS9SqggC8i8VjBfp4x5i0AY8w+Y0ypMaYM+DvVLGhujHnWGDPUGDO0ffv2XvVbhWDfli2Ockrz5t407L7xqwFfqYgSyCwdAZ4H1hlj/uhT7ztxezKwxvvuqXDYu3Ono9wxhKRpzoY04CsVyQK5U3cmcA3wnYissut+CVwpIhlY6Vi2Az8NSw+V5/a5VqNKadfOm4Y14CsV0QKZpfMpVsoVtw+8746qC3tdD0SlhPqUbTkN+EpFNH3SNgrtc61l2zGEhU+cDXXkPuAKYCzQf8UK9rhz9iil6o1my4xC+44dc5RTevb0puGOHVmINSUTgOJi9u3dSyev/oJQSoVEr/CjjTHsPX7cUZXSu7c3bbdqRYpreue+7du9aVspFTIN+NHm2DH2lZU5qjq6ctkHTYQOzZo5qvbrYuZKRQwN+FHG7NnDPlddSqiJ03x0cD2xu18XM1cqYmjAjzJHNm/GN8NNYkwMSUlJnrXfwTWnf79mzFQqYmjAjzIHtmwh1qfc0aO0CuVSOjiTJO/TqZlKRQydpRNlegFFwCFgL1Awfryn7XdwzcjZf+CAp+0rpYKnAT/a7N1LDPbC5QADB3rafIeuXR3l/bqurVIRQwN+tPF48XK3Dq45/fvz8k64/f79sGCB9fWqq6BPH0+7o5TyoWP40cbrpQ1dUvr2dZT3FxZijHHUGQPvvw9TpkBqKvzsZ/DAA9YfG3/4A5SWetolpZRNA360CXPAb96jB74z8QuN4ejRoxVlY+CGG2DCBFi4EEpKfLYthF/8AkaOhI0bPe2WUgoN+NEnzAGflJQqi5nv9xlGmjcPXnrpxE0sXw6DB8MXX3jbNaWinQb8aFJayi/27uVW4AHgr8Ch+Hhvj5GQQIfYWEfV/s2bATh6FGbNqrrLSSdBmzbOurw8uPRSa2xfKeUNDfjR5OBBXjeGZ4A5wC3A4YICzw+TkphY8b4VcGzXLsAap/f9A6NJE1iyBL7/HtasgR/9yNlOdrZ1I1fH9JXyhgb8KOI3rUJKiufHefKUU9gJHAdygfPT01m7Fp54wrndL34BY8eCCHTqBO+8A3fd5dxm8WKYM8fzLioVlQJZ4rCLiCwVkXUislZEZtj1bURkkYhssr96tE6eCpcjW7ZQ5FNu5nFahXLdunenC9DULps9e7n9dueVeteucM89zv1E4NFHYdQoZ/2DD1qzepRSoQnkCr8EuMsY0w84DfiZiJwM3A0sNsb0BhbbZRXB9rqmvnR0Zbb0jOtG8DuLk1i61LnJ3Lng7/BxcTB/ftV7yTfdBGEYfVIqqtQY8I0xe4wxX9vvfwDWAanAxcDL9mYvA5PC1UnlDXdu+hRXZkvPuKL144szHOXzzoPJk0+8+/z54Hvvd88eeOYZLzupVPSp1Ri+iHQHTgW+BFKMMXvA+qUAVWbjqQizz755Wi7FPTXGKz4BfxWD+HS38+nbxx+3hm9OZNQo+H//z1n3yCPgWqxLKVULAQd8EUkC3gTuMMYcrWl7n/2mi8hKEVmZk5MTTB+VR/a65uB37BCm39E+Af8pfub41rhxMGBAYM3Mng2+txgOHIA//9mLDioVnQIK+CISjxXs5xlj3rKr94lIJ/v7nQC/M6aNMc8aY4YaY4a2b9/eiz6rIO1zZa5MSU0Ny3EK27Thd8DPaMqLLAEuqfjebbcF3k67dnDHHc66P/wBcnM96aZSUSeQWToCPA+sM8b80edb7wLX2e+vA97xvnvKS/tcmStTunULy3FiO3fmHuBpCinldeAtoJiuXavOta/JXXdBcnJlOTfXuuGrlKq9QK7wzwSuAcaJyCr7dSHwKHCuiGwCzrXLKoLtdWWu9GwtW5e4lBQr9bLDHm65xZqFUxvJyVXn5s+dCwcPhtBBpaJUILN0PjXGiDFmoDEmw359YIw5aIw52xjT2/56qC46rIJUUMA+30xlQEp6eniOFRtLq9i2jqr4uG1MmxZcczNmQFuf5n74AV5/PYT+KRWl9EnbaJGdXfUp286dw3a4o+JcWGXEyWsI9hZOixZVx/LffDPIjikVxTTgRwmTlVUnaRUAVq+G/SWnOuoGdFgRUptXXOEsf/SRJlZTqrY04EeJA+vXO9IqtIiLC0taBbBm0kB3R13x0XUhtdmrFwwaVFkuK4O33w6pSaWijgb8KNHu6FEOAKuBD4DnzjknLMfZvr18fN05A2iHOw9/EC691FnWYR2lakcDfpSQ3btpCwwELgAuGzMmLMeZO7c8SVp3R/32Q6Hf07/kEmd5yRLwoFmlooYG/GiRne0sh+Ghq4MH4bnnykvOK/ydeXmUlZWF1H6/fnDyyZXlkhJ4992QmlQqqmjAjxZ1EPCffhry88tLrRBaVXyvyJgqqR2C4R7WWbAg5CaVihoa8KNFmAN+fn7VPDedKjLiW3a4snUGwz2ss2gRuB4gVkpVQwN+NCgrg927nXUeB/wXXrCSm5Vr2RIyYp059ravWRPycU45BXr3riwXFcF774XcrFJRQQN+NMjJ4YGSEn6DlRTpw2bNKIjx7qMvLi6filnp5puhd8tER50XAV+k6rDOW2/531Yp5aQBPxpkZ/MEcD9wI3BBfj65Hqac/Oc/YefOynLTplYu+24++RBigFz3XxlBmjLFWV6yRBc6VyoQGvCjQMHWrRz2KceK0MGjXPhlZdY6tL5uuMFKiT/llFNYAmzFWtD8d+ef78kxBw925tbJzYVvvvGkaaUaNQ34USD7++8d5c7NmhHru35gCN5+G9avryzHxsKsWdb7bn37MhboAcRD1RvHQYqJgbFjnXWLF3vStFKNmgb8KJC9ebOjnOqbYD4ExsDDDzvrrroKevQoP5DrxrBHAR/g7LOd5f/9z7OmlWq0NOBHgewdOxzlVI+GcxYtgsxMZ93dd/seqO4C/qefwvHjnjWvVKOkAT8KZO/Z4yinejQl84knnOVJk5xPwoYz4PfqBV26VJaPH4flyz1rXqlGKZAlDl8Qkf0issanbo6IZLtWwFIRKsuVcCa1YswleNnZ8OGHzjrH1T2ENeCLVL3K13F8pU4skCv8l4Dxfurn+q6A5W23lJeyjzofgEo76aSQ23zpJWuGTrmMDBg+3LVRSgrfifAkMBO49OBBXn3xxZCPXU4DvlK1E8gShx8DmpOwoSooILu42FGV2q9fSE2WlcHzzzvrpk2zrrod4uL4b4sW/Bx4HHgT+GLp0pCO7WvcOGd5xQpw/W5TSvkIZQz/NhH51h7yae1Zj5S3srNxD6Sk+g5+B2HZMti2rbLctClMnep/2+6udQ23u2YMhaJzZyuDZrnSUmslLKWUf8EG/GeAdCAD2IN1AeeXiEwXkZUisjInJyfIw6lgle7cyR5XXag3bd1X95dcAq2r+ZXfLS3NUd6RlRXSsd10WEepwAUV8I0x+4wxpcaYMuDvgHv01nfbZ40xQ40xQ9sHu4q1Ctr+devjVgpAAAAdn0lEQVQo8Sm3jo8nMTGx2u1rcvhw1ZWmpk2rfvvu6emO8vb9+zHGBH18N52Pr1Tgggr4ItLJpzgZCD0rlgqL7A0bHOXUli1Dam/ePCgsrCz36AEnWjyrXXo6vr9ejhUWcvjw4Wq3r60xY6wnb8utXQv73Ku1K6WAwKZlvgYsB/qKSJaITAN+LyLfici3wFjg/4W5nypI2b6D7UCabxKaILiHc264wRlw3SQtzbXYIWz3IC9+ueRkGDLEWbdsmWfNK9WoxNW0gTHmSj/Vz/upUxHonLg4VgNZQDaQctFFQbe1ejWsWlVZjomB66+vYafUVLoB63yqtm/fzuDBg4Puh9uYMdYMnXJLl8Lll3vWvFKNhj5p28g137ePgcCFwE3AxIkTg27r5Zed5fPOA9c92apSU8N6hQ9VE6l5OPNTqUZFA35j59HShsXF1vi9rxqv7u3j9XRVbdq4Mag+VOess6wsneU2bqy6wJdSSgN+4+bh0ob//S/s319ZbtUKAvpjoUUL+rpmBW10pWsOVYsWMGyYs07H8ZWqSgN+Y5aTAyU+kzKTk6FZs6Cacg/nXHYZBDq7s0/Hjo7yRo+v8EGHdZQKhAb8xsyj4ZxDh+Ddd511110X+P49u3fHd7mVrH37yMvLC6ov1dGAr1TNNOA3YnlbtlQuXA5836pVUO3Mnw9FRZXlXr3gjDMC379Jly6483Nu2rQpqL5U58wzIT6+srxlC+za5ekhlGrwNOA3YrtWr65cuBy46Ntvg2rHPZxz7bV+EqWdSFoafXyKiXFx7NnjTvgQmmbNYMQIZ51e5SvlpAG/Ectet85RTm3TptZtbNgAX37prLvmmlo20rMn9wCLgB3AsSlTuOCCC2rdl5q4n/jVG7dKOWnAb8R2btniKKd27lzrNp55xlkeMwa6d69lI+npnAWcA3QFYrZurXU/AqHj+EqdmAb8RmyTKzNlr1oufHL0KLzwgrPuhhuC6IgrgRphCvinnw5NmlSWt2+3Xkopiwb8xsoYNrqSlPV2J52pwUsvwQ8/VJZTUqzpmLWWmmolzS936BDk5gbR0IklJlpB35de5StVSQN+Y7VnDxt91yAE+px6asC7l5XBk086626+2Rm3AxYTY6XV9OUabvKKe1hHF0RRqpIG/EaqbNMm3GtL9enbN+D9P/gAfBenio+3An7QeroSLIQp4I8e7SxrwFeqkgb8Rip75UoKfMqt4+NpW4vUyE884SxfcQW4HpitnfR0vgWeBmYAF9x3H2+6V1LxwIgRVcfxd+70/DBKNUga8Bupjd984yj3ad8eCXDy/Nq1VVeOmjEjxA6lp/MG8DPgz8CHGzeywjensUcSE2G4a/21jz/2/DBKNUga8Bspd76aPl27Bryv++r+zDOrLjJSa+npjoevIDw5dUCHdZSqTiArXr0gIvtFZI1PXRsRWSQim+yv1SxhrerLRldegT79+gW0344d1uwcXyFf3UO9Bny9wlfKEsgV/kvAeFfd3cBiY0xvYLFdVhFk48GDjnKfAC/RH3nEyn1frkcPmDzZgw51714l4G/evJnS0lIPGnc6/fSq+fE9zuSgVINUY8A3xnwMHHJVXwyUZ1h5GZjkcb9UKI4eZaNv1AZ6uxPN+LFjR9UHrX71K4ircSHMACQm0jo1lfY+VYWFhewMwx3VpCQYOtRZp1f5SgU/hp9ijNkDYH/tUN2GIjJdRFaKyMqcnJwgD6dqo3j9era56noH8JTtQw85r+579gwib86J6Di+UvUq7DdtjTHPGmOGGmOGtm/fvuYdVMjiduxgO9ZY2zPAr9LTSUpKOuE+27fDiy866+67z5lyOGR1GPBHjXKW9QpfKQj2j/V9ItLJGLNHRDoB+2vcQ9UZ2baNNCANGAcw3n0LpqqHHnIujpWeDldf7XHH6jDgn3WW9YBv+cPGa9fCgQPQrl1YDqdUgxDsFf67QPmaR9cB73jTHeUJ91Os7uRlfjZ3z8y57z6Pxu599exZZwG/VSvIyHDW6VW+inaBTMt8DVgO9BWRLBGZBjwKnCsim4Bz7bKKFLUM+HPmOK/ue/WCqVO97xbp6biTO6xfvz4MB7LosI5STjVewxljrqzmW2d73BflFXfAd+ex8fHddzBvnrPu/vvDcHUPkJ5OOiCAsat27txJfn4+zYJcXP1ERo+GP/2psqwLoqhop0/aNjZFRVWTx5wg4N93HxhTWR4wAK6s7ld8qNq0IaFVq7Cvb1tu5EhnefVqcD2eoFRU0YDfyOStX8+MsjKewl5SsH17a8FXP774At5x3X156CHnQ0ueEqnTYZ22bauO42t+fBXNNOA3Mps+/ZQ/A7cB5wFn+65g4sMY+OUvnXWnnQY/+lGYO5ieju8TAUkJCRw5ciRshxs3zllesiRsh1Iq4mnAb2Q2ff21o9ynmpTI//tf1avdhx+2LsLDqmdPbsR6RiAbODpzJtOnTw/b4c523WlavDhsh1Iq4mnAb2Q2uoZH+nTpUmWb0lKYOdNZd845VVeLCov0dE7Gej6gMyBhWt+23MiRzhvQGzeCa6lfpaKGBvxGZqPrhq2/Va5efBG+/dZZ9/DD4eyVD/cU0c3udbm81aJF1fz4OqyjopUG/Ebm6717HeU+rixiP/xgJUTzNXUqDBsW7p6Vd8j16NWaNdafHGGkwzpKWTTgNyI5Gzawxif7WQww9JJLHNs8+ijs21dZTkiow6t7gNRUZ36D/HzYsCGsh/R349Z3KqpS0UIDfiPysesJqlMTE0nu1KmivGMHPP64c5+ZM6EWi2GFTqTq8lmZmWE95OmnW7/YymVlQZim/isV0TTgNyLLXAvRjnE9cDV7NhQWVpY7drTq6tyQIWwFnsSaPnrOvffywAMPhO1wTZtaydR86Ti+ikYa8BuRZWvXOspjfKLce+/B/PnO7R980FospM4NGcIq4OfAU8DiXbv47LPPwnpIHcdXSgN+o5GTk8Oao0cryjHAWfb4/ZEjcPPNzu0zMuD66+uufw5DhtRpEjWoGvCXLq1MnaxUtNCA30h8/MEHjnIGkGxf4c+aBdnZld+LjYXnngtjCoWadO1KrzZtHP/4du3aRV5eXtgOOXiwlTK53MGDsGpV2A6nVETSgN9ILFu40FEe064dJCayZAn8/e/ObWfNqnrftE6J0HTo0CpJ1MKVGx+sX27uB8veey9sh1MqImnAbySWffmlozxm0CCOHYObbnJu17evlf643vkZ1tkQ5umZEyY4y+++G9bDKRVxQgr4IrJdRL4TkVUistKrTqnaKSoqIrmkhPLlZwUYee653HEH+GYuEIHnn3dOUaw3gwfjXlY93OP4EyY4cwVlZmqaBRVdvLjCH2uMyTDGDK15UxUOTZo04ZNWrTiMlRL5CWBR8SSef9653e23w5ln1kMH/fF3hR/mgJ+SAiNGOOt0WEdFEx3SaQxyc2HLFpoD5wATpTs3PeZMYXDSSfDII/XRuWp0785Jrjmh690JfsJg4kRnWYd1VDQJNeAb4P9EJFNEwpfjVp2YT0rkUmK4OnEBR45Ujl00aQKvvVbtOij1Q4S+p57qqNq4dStlYZ4refHFzvLixVZ+IaWiQagB/0xjzGDgAuBnIjLKvYGITBeRlSKyMicnJ8TDKb98UhM8xL18mu+cgvO731Vd+SkSdDj9dHxmSpJfWEhWmAfV+/VzJuwsKoL/+7+wHlKpiBFSwDfG7La/7gcWAsP9bPOsMWaoMWZo+/btQzmcqo4d8D/hLB7AOQVn/Hj4+c/ro1M1k6FD6eeqywxzXh0RHdZR0SvogC8izUWkRfl7rBX11njVMRWYdevWUZqZySFaM5V5lFH5NFWHDvDSSxATqXdqhgzhNFfVxx99FPbDugP+++9DSUnYD6tUvQslFKQAn4rIauAr4H1jzIfedEsFoqCggCFDhtBx82ZOJYNdZAKVueVfecWamRKxevRgVPPmjqqP6mB85ayzoHXryvLBg7B8edgPq1S9CzrgG2O2GmMG2a/+xpiHvOyYqtn//vc/CgoKOADsZClWOjLrI505E84/vz57FwARRvqsvJIK9E9ICPuN27g4uOgiZ91bb4X1kEpFhEj9Y18F4N0qg88TAWHYMHiogfz6bXfJJSwAtgK7gH/ExRFTB2NQ7mGdefOsG7hKNWYa8BuosrIy3n3bHfAvpkXzUl57zZqK2SBMmcIlQA+sJ4RZscJaqSXMJkxwJlPLyYF//zvsh1WqXmnAb6C+/HIF+w/s96lpAYzm2b/HVFknPKJ17gxnnOGsq4PxlcREuPpqZ91zz4X9sErVKw34DdS997qv7sfzs1NXcsWV4nf7iHbppc7ym2/WyWFvvNFZ/u9/YefOOjm0UvVCA34D9K9/wdKlzoDfk5N5/M8NZRzHZcoUZ/mzz2D37rAfNiPDmSbaGHjxxbAfVql6owG/gXnxRbjiiq04H3mI5a2U92l6ZgPNX9etG/jM1skHttfR+Ir7Kv+FF6C01P+2SjV0GvAbkKeeghtuAGNeddSfQksGXXmWM/dvA5N1zjncA5wBJAM3zJ1bJ8e98kprPL/czp3gWgteqUZDA34DkJUF110Ht90GcBj4o+P7N3AY7PVrG6ric8/lUWA5UAwsz82lcNeusB+3VSu47DJnnXuFMKUaCw34EezIEbjvPujTx3pq1vIYcKRim9bATzp0qDrTpYHpPmYMafHxFeXjwMo6usp3D+ssXAhrNEmIaoQ04EegzExracLOneHBB6GgoPw7+7GWN6n0C6DVpZdGcMKcwIgIo/o5U6m9/be/Wb/1wuzMM+HkkyvLZWVw553WTVylGpOGHSUakZISa/bNiBEwdKg1Jzw/373Vo0BeRakDcDvAtGl11s9wGnPVVY7yn/LzybzllrAfV6Tqk8mLFsEHH4T90ErVKQ349SwvD5580hq2ufxy+Oor/9ulpMCt1w7hZKn8yH4JNL/xRhg8uG46G2aX3XwzKT7J1EqA6157jcIwp0wGa2GUMWOcdXfdBcXFYT+0UnVGA349OXQIfvMba0biz38O27b53y41FX77W9i0CZ4qeo9vTRmvAmOAnyYnR9i6haFp1aoVz7omwq8F7p88OezjKyIwd65zotOGDfDMM2E9rFJ1yxhTZ68hQ4aYaHb8uDHvv2/MjTca07y5MVYU8/8aP96Yd94xprjY3nnRoqobPf10vZ5PuFw/bpzBWj7TACYGzEf33lsnx542zfkjbt3amB076uTQSlULWGk8iMFi6vDO1NChQ83KlSvr7Hj1qbjYukJcs8Z6ffcdLF164vVTmzSBa6+1bhg67l+++ab1Dd9B/cGDrfGf2Ngq7TR0R3JzGdCpE1nHj1fUNQF+M3kyd/3rX8TFxYXt2Hv3Qu/ecOxYZV3PntZn17Vr2A6r1AmJSKYxJvQnK734rRHoq7Fe4ZeVGbNxozHPPGPMtdcaM2iQMfHxJ76C9321bGnM3Xcbs2dPZZv79+83H/7nP8Y88ID/nT7/vP5OuA7894UXHFf55a8RKSlm/Zo1YT32o49W/XH37KlX+qr+4NEVfmg7w3hgA7AZuLum7RtLwC8oMOaLL4x58kkrwHfpEnhw93117GjM735nTG6uMQUFBWbt2rVm4cKF5urLLjNN4uKMgHnW34733VffP4I6cc/55/sN+u9162bM3LnG7N0bluMWFxtz6aVVf+w9ehjz5pvGFBWF5bBKVcurgB/0kI6IxAIbgXOBLGAFcKUx5vvq9hlw0gDzxvMLanUcYyA+Lp5WLVtV/tdDMMbKebI/5zAFx8soKhby8iEvP4a8fCgoEIpLhKIiKCm1bsbFxlROVzdGaNasLSUlwvFCa/vjhULBcThy9Afy8osoKhIKi6G4SMgvEA4eFg7lxnA4VzCmujQGrQF/Qw7fAytp3mw7XTpuISlxBzHmAEd/OMqRY8fYm5tLdZ/F77Dm2xMXB08/bU3SjwLGGP41cya3/vGPHLLrOmD9Y4sHazhrxAg46STo04enN24kt7SU9J496dmrFy2Sk2narBlNmzUjrkkT68MXARHE5+5sXHw8rZOTKw8sQnExXHVzSxb8Ow/fZSMB2rct47KLj9M3vZSOKWWktC8joakhLg7iYiEmxtjbtStvzrdpfjj2A4VFhUH9TJJbJvsd0io4XkBefp6fPWqW1DyJhKYJVeqLi4s58kNwz0EkNE0gqXlSlXpjDAcPHwyqzbjYOJJbJfv93uHcw5SWBZcEqV2bdn7rI+Vzio+Lp/dpvep3SAc4HfivT/ke4J4a9vF7xVbza/QJrpR7B9lm0gnavDLINjGwoaKdVHaZ8XxgZvJ7cyH9g24zGczm1q2N+eijcFw8RLw9ixaZSYmJBjB3nOBPpr5B/nxHVNNeEXGmFV2D/NziTvDv6/oQ/n2trqbNv4TQ5gvVtLk8hDZvrabNohDaHHGCn+kpjfhzGm3w6Ao/lLtfqVir0pXLAkaE0F6jcCePcyGbGcRq2lF5JfMIUNvneNKBW4DrR4+m7YsvQo8eHva04eh4zjm8tW0b702eTJ8vvrD+37gUAJs8Pm48JXRhJ+F/1lepuhFKwPc3plHlf6KITAemh3CcBuWnPEsfP/X+/xCtJEAXoDfQF5jYqRPnTp9OzLXXWtNEopykpPCjzz+H7Gz45z/h5Zdh7dqK768DwrH0ecPNP6pUVaEE/CysGFUuDaiyaoUx5lngWQCROCM1hr6q4kmkleQAIBjE/hojZRwsa0YJrX2+5/uCqr+DBDDE0ZRrE14gllISpYBEjpMghSTKcd4sXMOqkgQQq83ytmMos78aP+1aYlu3tsbaKw4nEBND36IifpyfT9f4eLo1bUqXpk1p36IFrVq1Ijk5mXapqSScfLL1yG3fvtCrV4NOdxw2qakwaxbMnGmtfbt+PWzcSMrXX/PnVavYcvgwW44eZVdBAQVlZRQaQ2FZWeUovDFVPrnkuDgoH8N3/fXQOjeXdj4J8g1CkWlCCXGUEWO/Kv+1lW8jxNJO9lepBzhm4iikbVCn34qjxMmBKvXHTTF5QbbZnCIS/LRZbI5xNMg2E4iluZ82jSnmUJBtxtGMVn7aBMg1SZQG1W4sbatpM1I+p3gSKQqqF1WFctM2Duum7dlANtZN26uMMWur2yea5uErpZRXvJqHH/QVvjGmRERuA/4LxAIvnCjYK6WUql8hPbJojPmA2t+LVEopVQ80eZpSSkUJDfhKKRUlNOArpVSU0ICvlFJRQgO+UkpFiTrNhy8iP2Bl12ys2gH+n+JoHBrz+TXmcwM9v4aurzGmRaiNhG8lCf82ePHwQKQSkZV6fg1TYz430PNr6ETEkydWdUhHKaWihAZ8pZSKEnUd8J+t4+PVNT2/hqsxnxvo+TV0npxfnd60VUopVX90SEcppaKEZwFfRMaLyAYR2Swid1ezzWUi8r2IrBWRf/rUXycim+zXdV71ySshnlupiKyyX+/WXa8DV9P5ichcn3PYKCK5Pt+L6M8OQj6/xvD5dRWRpSLyjYh8KyIX+nzvHnu/DSJyft32PDDBnp+IdBeRAp/P76913/sTC+DcuonIYvu8lolIms/3av9/z4t1ErHSI28BegJNgNXAya5tegPfAK3tcgf7axtgq/21tf2+tRf9qu9zs98fq+9zCPX8XNvfjpUKO+I/u1DPr7F8fljjv7fY708Gtvu8Xw00BXrY7cTW9zl5eH7dgTX1fQ4hntsbwHX2+3HAP+z3Qf3f8+oKfziw2Riz1RhTBLwOXOza5ibgKWPMYQBjzH67/nxgkTHmkP29RcB4j/rlhVDOrSEI5Px8XQm8Zr+P9M8OQju/hiCQ8zNAS/t9KypXprsYeN0YU2iM2QZsttuLJKGcX6QL5NxOBhbb75f6fD+o/3teBXx/C5qnurbpA/QRkc9E5AsRGV+LfetTKOcGkCAiK+36SeHubBAC/vmLSDesK8Eltd23HoVyftA4Pr85wNUikoW1fsXttdi3voVyfgA97KGej0RkZFh7WnuBnNtq4BL7/WSghYi0DXDfKrx60jaQBc3jsIY+xmCtf/uJiAwIcN/6FPS5GWNyga7GmN0i0hNYIiLfGWO2hLXHtVObn/8VwAJjTPkir5H+2UFo5weN4/O7EnjJGPO4iJwO/KOB/N+D0M5vD9bnd1BEhgBvi0h/Y8zRMPc5UIGc20zgLyJyPfAx1nKyJQHuW4VXV/iBLGieBbxjjCm2/3zcgBUkA1oMvR6Fcm4YY3bbX7cCy4BTw93hWqrNz/8KnMMdkf7ZQWjn11g+v2nAvwCMMcuBBKzcM43l8/N7fvZQ1UG7PhNrvLxP2HscuBrPzRiz2xgzxRhzKnCvXXckkH398ujmQxzWTYMeVN586O/aZjzwsv2+HdafI22xbjpsw7rx0Np+36a+b6h4dG6tgaY+9Zs4wQ3DSD0/e7u+wHbsZzdM5Y2jiP3sPDi/RvH5Af8Brrff97MDgwD9cd603Urk3bQN5fzal58P1o3R7Ej69xngubUDYuz3DwG/sd8H9X/Py85fCGzE+i16r133G2Ci/V6APwLfA98BV/jsewPWDaPNwE/q+4Pw6tyAM+zyavvrtPo+l2DOzy7PAR71s29Ef3ahnF9j+fywbvx9Zp/HKuA8n33vtffbAFxQ3+fi5flhjX2vteu/Bn5U3+cSxLldinWhsRF4DvsCxP5erf/v6ZO2SikVJfRJW6WUihIa8JVSKkpowFdKqSihAV8ppaKEBnyllIoSGvCVZ0QkTUTesbP3bRWRv4hI0zAcZ4yInOFTvllErrXfXy8inT06zkg7++kqEUl0fe/zAPa/Q0SaedGXExwjwzf7pVInogFfeUJEBHgLeNsY0xvrSeNE4PdhONwYrDnyABhj/mqMecUuXg94EvCBqcBjxpgMY0yB7zeMMWdUs4+vO4BaBXwRia3N9kAG1lxupWpW3w8e6KtxvICzgY9ddS2Bw0ASViD+i8/33gPG2O+fAVZiPSTzgM8224EHsB6a+Q44CSvl7V6spyZXASOxHpqaifWQyjGsh4hWARcBC33aOxd4q5q+f2Mf4wWsJ09vBA5hPcE4z88+x+yvY7BSLiwA1gPzsB7E+zlQZLe51N72PGC5fT5vAEk+5/lr4FOs9A7pwIdAJvAJcJK93Y+BNVgPEn2M9XTmTiDHPt/L6/vfgb4i+1XvHdBX43jZAW6un/pvsK5CTxTw29hfY+3gOdAubwdut9/fCjxnv58DzPRpq6Js7z/Ufi92EG5vl/+J62lLrLwru4A+dvkV4A77/UvApdWcr2/AP4KVyyTGDuhn+fS/nf2+nR2km9vl2cCvfbb7hU/bi4He9vsRwBL7/XdAqv0+2f7q+LnqS18neumQjvKK4D9bn7+sfm6XicjXWL8c+mM9Kl/uLftrJtbVfcCMMQb4B1bq3GTgdKy8K776AtuMMRvt8svAqNocB/jKGJNljCnDutL218/TsFMAiMgq4Dqgm8/35wOISBLWcNUb9nZ/AzrZ23wGvCQiN2H9clSqVrxKj6zUWirzdgMgIi2BFKwhlgE47xkl2Nv0wBqOGWaMOSwiL5V/z1Zofy0luH+vLwL/Bo4DbxhjSlzfD+QXUk0Kfd5X10/BWrDiymrayLO/xgC5xpgM9wbGmJtFZATWUNUqEamyjVInolf4yiuLgWY+s2VigcexhhsKsIYtMkQkRkS6ULmyUkusYHdERFKACwI41g9Ai0C+Z6z0xruBX2EN0bitB7qLSC+7fA3wUQB9CIRvX74Aziw/jog0E5EqqXqNlat9m4j82N5ORGSQ/T7dGPOlMebXwAGs9Lgn+lko5aABX3nCHj6ZDFwqIpuAg0CZMeYhe5PPsG6Afgc8hnXjEmPMaqyhnLVYN0w/C+Bw/wYm29Ml3asYvQT81TWVch6wyxjzvZ9+Hwd+gjWE8h1QBni12PWzwH9EZKkxJgdrvP01EfkW6xfASdXsNxWYJiKrsX4u5cva/UFEvhORNVj3A1ZjLXt3sn2+l3vUb9VIabZMFRb2PPnXgCnGWnyiPvvyF+AbY8zz9dkPpeqbBnzVqIlIJtaQ0bnGmMKatleqMdOAr5RSUULH8JVSKkpowFdKqSihAV8ppaKEBnyllIoSGvCVUipKaMBXSqko8f8BWnGz+WqFKS8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Compare the observed and the pushforwards of prior and posterior\n",
    "qplot = np.linspace(0.6, 0.9, num=100)\n",
    "plt.clf()\n",
    "\n",
    "q_kde_plot = q_kde(qplot)\n",
    "obs_vals_plot = norm.pdf(qplot, loc=0.7, scale=1.0e-2)\n",
    "postq_kde = GKDE( post_q, 'silverman' )\n",
    "postq_kde_plot = postq_kde(qplot)\n",
    "\n",
    "oplot = plt.plot(qplot,obs_vals_plot, 'r-', linewidth=4, label=\"Observed\")\n",
    "prplot = plt.plot(qplot,q_kde_plot,'b-', linewidth=4, label=\"PF of prior\")\n",
    "poplot = plt.plot(qplot,postq_kde_plot,'k--', linewidth=4, label=\"PF of posterior\")\n",
    "plt.xlim([0.6,0.9])\n",
    "plt.xlabel(\"Quantity of interest\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEKCAYAAAAB0GKPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGdNJREFUeJzt3XuULXV55vHvo4hcFZWoCByPGtRlXEaxw5ioEYMaBUWdYAQvAwwuJmpUdJxIbupoMoMZb8zSyJyoEa8gqAMKEhFBvATC4SKIiCAiHDkjx3hHvKDv/FHVy22zu7u6e9/7+1mrV+9du2rX05tDvfXWb1dVqgpJkpZzh3EHkCRNBwuGJKkTC4YkqRMLhiSpEwuGJKkTC4YkqRMLhiSpEwuGJKkTC4YkqZPtxh1gLXbffffauHHjuGNI0lS5+OKLv1NVv7XS5aa6YGzcuJHNmzePO4YkTZUk31zNch6SkiR1YsGQJHViwZAkdTK0gpHk3UluTvLlnml3T3J2kmva33drpyfJ/05ybZLLk+w7rFySpNUZZofxHuDJC6YdC5xTVfsA57TPAZ4C7NP+HA28Y4i5JEmrMLSCUVXnA99dMPnpwInt4xOBZ/RMf281LgB2S7LHsLJJklZu1GMY96qqrQDt73u20/cEbuyZb0s7TZI0ISZl0Dt9pvW9d2ySo5NsTrJ527ZtQ44lSZo36hP3vp1kj6ra2h5yurmdvgXYu2e+vYCb+r1BVW0CNgHMzc15Q3KNxcZjz1jy9euPO2hESaTRGXWHcTpwePv4cOC0nun/qf221KOAH8wfupIkTYahdRhJPgTsD+yeZAvwGuA44MNJjgJuAJ7Vzn4mcCBwLfAT4Mhh5ZK6Wq6LkNaboRWMqjpskZcO6DNvAS8eVhZJ0tpN9cUHpWFayzjEsLsTx0g0DpPyLSlJ0oSzw5CGyE5As8QOQ5LUiR2GtAJdxybsLDSL7DAkSZ3YYUgrsLBz8FwNrSd2GJKkTuwwpDXo7TjsNjTr7DAkSZ3YYUgDspJuw29RaRrZYUiSOrFgSGOw8dgzHPPQ1LFgSJI6cQxDGqD5rsExCs0iOwxJUicWDGmA7Cw0yywYkqROHMOQBsBvPGk9sMOQJHVihyENwCCuYuv4hyadHYYkqRM7DGmI7Bo0S+wwJEmd2GFIQ2BnoVlkhyFJ6sSCIQ2Y52RoVlkwJEmdOIYhDdj1xx20aJfh2IammR2GJKkTOwxplZYaq7CT0CwaS8FI8nLgBUABVwBHAnsAJwF3By4Bnl9VPx9HPqkLi4LWm5EXjCR7Ai8FHlJVtyb5MHAocCDwlqo6KckJwFHAO0adTxoExzA0i8Z1SGo7YMckvwB2ArYCfwQ8p339ROC1WDA0pSwMmkUjLxhV9a0kbwRuAG4FPgVcDHy/qm5rZ9sC7DnqbNK8tZ5LYcHQLBrHIam7AU8H7gd8HzgFeEqfWWuR5Y8GjgbYsGHDkFJqvXODL93eOA5JPQH4RlVtA0jyUeAPgN2SbNd2GXsBN/VbuKo2AZsA5ubm+hYVaVBW22lYcDSLxlEwbgAelWQnmkNSBwCbgXOBQ2i+KXU4cNoYskm/wQ2/9GvjGMO4MMmpNF+dvQ24lKZjOAM4KcnftdPeNeps0lr0diMWGs2iVE3vUZ25ubnavHnzuGNI0lRJcnFVza10Oc/0lobAs8A1iywY0hBYFDSLvPigNGDeD0OzyoIhDZjdhWaVh6SkZaymY7BoaBZZMKRluPGXGhYMaQhW2pVYlDQNLBjSEFgANIsc9JaGwG9KaRZZMKQhsMPQLPKQlDQECzsMC4hmgQVDGgILhGaRh6SkAdp47BmOX2hm2WFIA2RnoVlmwZCGYFBdhgVIk8SCIQ2BG3rNIscwpCFwHEOzyIIhDYEdhmaRBUOS1IkFQ5LUiYPe0gp0HZvwkJRmkQVDWgELgdYzC4a0Bot1HBYWzSILhrQGFgatJxYMaQhWex6GBUiTzIIhDYEbfs0iC4Y0IL1dhQVDs8iCIQ2IRUKzzoIhDYHdhmaRBUMaAouEZpEFQ+pgpd96smBoFo2lYCTZDXgn8FCggP8MXA2cDGwErgf+tKq+N4580kIWAGl8HcbxwFlVdUiS7YGdgL8Czqmq45IcCxwLvGpM+aQVWdiBWGA0i1JVo11hchfgS8D9q2flSa4G9q+qrUn2AM6rqgct9V5zc3O1efPm4QaWVsET9zTJklxcVXMrXW4clze/P7AN+OcklyZ5Z5KdgXtV1VaA9vc9+y2c5Ogkm5Ns3rZt2+hSS9I6N45DUtsB+wIvqaoLkxxPc/ipk6raBGyCpsMYTkRpcOwaNCuW7TCS7JTkb5P8U/t8nyRPXcM6twBbqurC9vmpNAXk2+2hKNrfN69hHZKkAevSYfwzcDHw++3zLcApwCdWs8Kq+n9JbkzyoKq6GjgA+Er7czhwXPv7tNW8vzRplhrPsPvQNOlSMB5QVc9OchhAVd2aJGtc70uAD7TfkLoOOJKm2/lwkqOAG4BnrXEdkqQB6lIwfp5kR5rzJUjyAOBna1lpVV0G9BuhP2At7yt1tdpvMXVl56BZ1KVgvAY4C9g7yQeARwNHDDOUJGnyLFswqursJJcAjwICvKyqvjP0ZNIQLdcBDLsDkabRogUjyb4LJm1tf29IsqGqLhleLEnSpFmqw3hT+3sHmvGGL9F0GA8DLgQeM9xo0visZQzC7kSzatHzMKrq8VX1eOCbwL5VNVdVjwQeAVw7qoCSpMnQZdD7wVV1xfyTqvpykocPMZM0VewotF50KRhXJXkn8H6ar9Y+D7hqqKkkSROnS8E4Engh8LL2+fnAO4aWSJowdhBSo8vXan+a5O3Ap2k6jKur6hdDTyZJmijLFowk+wMn0twFLzQn8B1eVecPN5o0XnYW0m/qckjqTcCT2gsFkuSBwIeARw4zmCRpsnQpGHeaLxYAVfW1JHcaYiZpInQ5F8MuROtJl4KxOcm7gPe1z59Lc7lzSdI60qVgvBB4MfBSmjGM84F/HGYoaVr060LsOjSrunxL6mfAm9sfSdI6tdTFBy9fasGqetjg40iD5d6+NDhLdRi/ojnv4oPAx4FbR5JIkjSRlrr44MOBw4BdaIrG3wO/A3yrqr45mnjS2lx/3EHe/U4akEULBkBVfbWqXlNV+9J0Ge8FXj6SZJKkibLkoHeSPYFDgWcC36MpFh8bQS5poOwypLVbatD7s8CuwIdp7uH93fal7ZPcvaq+u9iykqTZs1SHcV+aQe//AhzdMz3t9PsPMZckacIsWjCqauMIc0iSJlyXM70l9ehybodjJppFS35LSpKkeXYYWvc8G1zqZlUdRpIbBh1EkjTZVtthZKAppDGws5BWZrVjGDXQFJKkibfUiXuvWOwlmutLSVPNO+pJK7NUh7HrIj+7AMevdcVJ7pjk0iSfaJ/fL8mFSa5JcnKS7de6DknS4Cw1hvGuqtrS74UkTxvAul8GXAXcpX3+BuAtVXVSkhOAo4B3DGA90qqt9nyKUXYmnvOhUVmqwzgnycaFE5McCbx1LStNshdwEPDO9nmAPwJObWc5EXjGWtYhSRqspTqMlwNnJzmwqq4BSPKXwHOAx61xvW8F/oLmEBfAPYDvV9Vt7fMtwJ5rXIfUyTC6Aff6NYuWupbUmUl+BnwyyTOAFwC/B/xhVX1vtStM8lTg5qq6OMn+85P7RVhk+aNpL4a4YcOG1caQJK3QkudhVNU5SY4AzgO+CBxQVT9d4zofDRyc5EBgB5oxjLcCuyXZru0y9gJuWiTTJmATwNzcnF/v1ZoNoxtYTddiV6JJt+gYRpIfJfkh8EmajfoBwM0901elqv6yqvZqr4Z7KPCZqnoucC5wSDvb4cBpq12HJGnwljoktetirw3Jq4CTkvwdcCnwrhGvX+vcIMcy7BY0i8Z68cGqOo/mcBdVdR2w3zjzSJIW59VqNZM8Q1saPO+HIUnqxA5DM2kcYwh2NZp1dhiSpE7sMLSu2RVI3dlhSJI6scPQujassY6FnYvnZWgW2GFIkjqxw5CW4XWhpIYdhiSpEzsMaRnLdQt+00rrhR2GJKkTOwxpjfp1ILPUdTgeo3l2GJKkTuwwNDHcK5cmmwVDE8ONrDTZLBgaiVnqHrqw+GkWWTA0Em5ApelnwZBWqEu3ZIHULLJgSCtkMdB6ZcHQujDqMRSLimaRBUPrghtwae0sGFIfa+1ILFCaRRYMqQ83+NLtWTCkIRjFmIlFTaNmwZCGwI25ZpEFQ+veMLoBC4ZmkQVD654bd6kbC8Y6t96u8TQqFiHNIgvGOueGTVJXFgyN1ax2OBZizaKRF4wkewPvBe4N/ArYVFXHJ7k7cDKwEbge+NOq+t6o82m03LBK02McHcZtwH+tqkuS7ApcnORs4AjgnKo6LsmxwLHAq8aQT1NuEroWC6Fm0cgLRlVtBba2j3+U5CpgT+DpwP7tbCcC52HB0Cq4sZaGY6xjGEk2Ao8ALgTu1RYTqmprknuOMZrWaBL28sfJoqVZNLaCkWQX4CPAMVX1wyRdlzsaOBpgw4YNwwuoNXGDKc2esRSMJHeiKRYfqKqPtpO/nWSPtrvYA7i537JVtQnYBDA3N1cjCSyt0FIdlsVU02oc35IK8C7gqqp6c89LpwOHA8e1v08bdTZpUCwKmkXj6DAeDTwfuCLJZe20v6IpFB9OchRwA/CsMWSTBmKlYzgWGE2DcXxL6vPAYgMWB4wyizQsFgDNIs/01royqm9vWTA0iywYWlfckEurZ8GQhmC9n4cySBb5yWHBkIbAjZxmkQVDGgI7jOlhce/OgiENgRshzaI7jDuANGvsLjSrLBjSgNldaFZ5SEoagmnuMix4WowFQxoCN7qaRRYMaQgG1WFYeDRJLBjSELih1yyyYEhDMKljGBYyrYUFQxoCN8yaRRYMaQgcw9AssmBIQ+CGXrPIgiENwaSOYWjyTNPOhQVDGqNp2lhIXhpEktSJHYa0iLUcVrJz0Cyyw5AkdWKHoYnhQLE02ewwJEmd2GFoYszScX+7pcGapX8b08wOQ5LUiR2GpIGzI5hNdhiSpE7sMKQRcs9b08wOQ5LUiR3GFPGbN9PDTkKzaKI6jCRPTnJ1kmuTHDvuPJKkX5uYDiPJHYG3A08EtgAXJTm9qr4y3mSTw73W6TGObtB/Hxq2Seow9gOurarrqurnwEnA08ecSZLUmpgOA9gTuLHn+RbgP4wpizR11uMYl13VaE1Sh5E+0+p2MyVHJ9mcZPO2bdtGEEuSBJPVYWwB9u55vhdw08KZqmoTsAlgbm7udgVFmgTu+WoWTVKHcRGwT5L7JdkeOBQ4fcyZJEmtiekwquq2JH8O/AtwR+DdVXXlmGNJkloTUzAAqupM4Mxx55Ak3d4kHZKSJE0wC4YkqRMLhiSpEwuGJKkTC4YkqRMLhiSpk1RN78nSSX4EXD3uHB3sDnxn3CE6MOfgTENGMOegTUvOB1XVritdaKLOw1iFq6tqbtwhlpNkszkHZxpyTkNGMOegTVPO1SznISlJUicWDElSJ9NeMDaNO0BH5hysacg5DRnBnIM20zmnetBbkjQ6095hSJJGZKoKRpK7Jzk7yTXt77stMe9dknwrydtGmbFd97I5k9w3ycVJLktyZZI/m9CcD0/yr23Gy5M8exJztvOdleT7ST4xwmxPTnJ1kmuTHNvn9TsnObl9/cIkG0eVbUGO5XL+YZJLktyW5JBxZGxzLJfzFUm+0v5bPCfJfSc0558luaL9//vzSR4yiTl75jskSSVZ+hteVTU1P8A/AMe2j48F3rDEvMcDHwTeNok5ge2BO7ePdwGuB+4zgTkfCOzTPr4PsBXYbdJytq8dADwN+MSIct0R+Dpw//a/55eAhyyY50XACe3jQ4GTR/nZrSDnRuBhwHuBQ0adcQU5Hw/s1D5+4QR/nnfpeXwwcNYk5mzn2xU4H7gAmFvqPaeqwwCeDpzYPj4ReEa/mZI8ErgX8KkR5Vpo2ZxV9fOq+ln79M6Mp9vrkvNrVXVN+/gm4Gbgt0aWsNHpv3tVnQP8aFShgP2Aa6vquqr6OXASTdZevdlPBQ5I0u/+9cO0bM6qur6qLgd+NeJsvbrkPLeqftI+vYDmVs6j1iXnD3ue7gyMY7C4y79PgNfT7JT9dLk3nLaCca+q2grQ/r7nwhmS3AF4E/DfRpyt17I5AZLsneRy4Eaavebb3cN8yDrlnJdkP5o9la+PIFuvFeUcoT1p/tvN29JO6ztPVd0G/AC4x0jS9cnQ6pdzEqw051HAJ4eaqL9OOZO8OMnXaTbGLx1Rtl7L5kzyCGDvqup0GHfizvRO8mng3n1e+uuOb/Ei4MyqunGYO3IDyElV3Qg8LMl9gP+b5NSq+vagMsJgcrbvswfwPuDwqhr4Xuigco5Yv39gC/cku8wzbJOQoYvOOZM8D5gDHjfURP11yllVbwfenuQ5wN8Ahw872AJL5mx3rt8CHNH1DSeuYFTVExZ7Lcm3k+xRVVvbDdjNfWb7feCxSV5EMzawfZIfV9WiAz5jytn7XjcluRJ4LM1hi4nKmeQuwBnA31TVBYPMN8icY7AF2Lvn+V7Awi5xfp4tSbYD7gp8dzTxbpdhXr+ck6BTziRPoNmReFzPYd1RWunneRLwjqEm6m+5nLsCDwXOa3eu7w2cnuTgqup76ZBpOyR1Or+u0ocDpy2coaqeW1Ubqmoj8ErgvYMuFh0smzPJXkl2bB/fDXg0o7+QYpec2wMfo/kcTxlhtl7L5hyTi4B9ktyv/ZwOpcnaqzf7IcBnqh1pHKEuOSfBsjnbQyj/Bzi4qsa149Al5z49Tw8CrhlhvnlL5qyqH1TV7lW1sd1eXkDzuS5+nalRj9yvcdT/HsA5NB/+OcDd2+lzwDv7zH8E4/mW1LI5gScCl9N8c+Fy4OgJzfk84BfAZT0/D5+0nO3zzwHbgFtp9q7+eATZDgS+RjOu89fttNe1/+MB7ACcAlwL/Btw/1H/d+6Y8/faz+wW4N+BKyc056eBb/f8Wzx9QnMeD1zZZjwX+J1JzLlg3vNY5ltSnuktSepk2g5JSZLGxIIhSerEgiFJ6sSCIUnqxIIhSerEgqHO2qtZvq/n+XZJto3y6rAL8jy4vRropUkeMKYMr2tPJFvNsvsP4rNL8tIkVyX5wFrfa5H3PybJTqtYbtWfjSbTxJ3prYl2C/DQJDtW1a0055J8a4x5ngGcVlWvGdYKktyxqn652OtV9ephrXsFXgQ8paq+MaT3PwZ4P/CT5Wac135uK/pslvusNX52GFqpT9KcuQpwGPCh+ReS7Jzk3Ukuavf6n95O35jkc+39Fi5J8gft9P2TnJfk1CRfTfKBfldyTXNPjgvaeyB8LMndkhxIsyF7QZJz+yzz4yRvSHPPkU8n2a9d13VJDu6Q69wkHwSuaKf9bZvx7CQfSvLKdvp70t4/Isn1Sf57+15XJHlwO32/JF9sP5MvJnlQn7yPa7ul+Y5p1z7zvCLJl9ufY9ppJ9Bcvvr0JC9fMP8RSU5Lc5+Qq5O8Zpn32jnJGUm+1E5/dpKX0lzW/tz5zznJk9LcI+WSJKck2aXn7391ks8Dz1rw2RzQ/l1XtP9G7txvmYV/sybMOM4+9Gc6f4Af09wz4VSaM5gvA/anvf8E8D+A57WPd6M5w3RnYCdgh3b6PsDm9vH+NFdv3Ytm5+Vfgcf0We/lNNcNguYs1be2j18LvHKRrEWz1w3NpU0+BdwJ+F3gsnb6UrluAe7XPp9r/9Ydaa6/c838eoH30N4/guaeJi9pH7+IX58tfxdgu/bxE4CP9Kxn/rP7OPDo9vEu8/P3/D2PpCleO7evXwk8ome9u/f5DI6guX/JPdrsX27/lr7vBfwJ8E89y9914fsDu9PcO2Hn9vmrgFf3zPcXPcu/h+ZyKDvQXDX1ge309wLH9FvGn8n+8ZCUVqSqLk9z17jDgDMXvPwk4OD5vW+aDcUGmguevS3Jw4Ff0tyUad6/VdUWgCSX0dzI5/PzLya5K80Nmz7bTjqR5lIby/k5cFb7+ArgZ1X1iyRXtOuApoAslWv+EM9jaA593dpm+vgS6/1o+/ti4D+2j+8KnJjm+kLVrnehLwBvbschPjr/mfR4DPCxqrqlzfBRmotVXrpEFoCzq+rfe5Z5TJuh33udBbwxyRtoCtnn+rzfo4CHAF9om8HtaQr9vJP7LPMg4BtV9bX2+YnAi4G3LrGMJpAFQ6txOvBGmj3k3ns7BPiTqvqNiygmeS3N9X9+l6aT6L1RS+/VRn/J4P5N/qLaXViamwL9DKCqfpXmqrEAL18i1y29f8IK1jv/9/T+La8Hzq2qZ7bF9ryFC1XVcUnOoLn2zwVJnlBVX11lht946z7P+75XVX0tzc3HDgT+Z5JPVdXrFswWmiJ02CLru6XPtOWy91tGE8gxDK3Gu4HXVdUVC6b/C/CS+XGINFcWhWYPe2s199F4Ps2tIzupqh8A30vy2HbS84HPLrHISnTN9XngaUl2aI/XH7TIfEutZ/7LAUf0myHJA6rqiqp6A7AZePCCWc4HnpFkpyQ7A8+kudjicp6Y5p7oO9J8SeALi71Xmvuy/KSq3k+zQ7Bv+x4/ojkUB80VTR+d5Lfb3Dsl6e3M+vkqsHF+GQb731AjZIehFWsPlxzf56XX0xxmuLwtGtcDTwX+EfhIkmfRXLlzpXuUhwMnpPlq53XAkauMvlCnXFV1UZLTaa4s/E2aDfoPVrCef6A5JPUK4DOLzHNMksfTdCZfYcGd5KrqkiTvobniLTTjI8sdjoKm2L0P+G3gg9VeurrfeyX5Y+B/JfkVzRWKX9i+vgn4ZJKtVfX4JEcAH5ofuKa5OdD84abbqaqfJjkSOKXt7i4CTuiQXRPGq9VKHSTZpap+3Bat82kuR3/JuHMtpd2wz1XVn487i2aDHYbUzaYkD6EZyD9x0ouFNAx2GJKkThz0liR1YsGQJHViwZAkdWLBkCR1YsGQJHViwZAkdfL/AcW6fGBBR70sAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.clf()\n",
    "param = np.arange(100)\n",
    "plt.barh(param,np.mean(post_lam,1))\n",
    "plt.xlim([-0.4,0.4])\n",
    "plt.xlabel(\"Mean of marginals of posterior\")\n",
    "plt.ylabel(\"KL Mode\")\n",
    "plt.show()"
   ]
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}