taku-y/advi.ipynb

## advi.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Warning: patsy not found, not importing glm submodule.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/taku-y/anaconda/anaconda/envs/py35con/lib/python3.5/site-packages/IPython/html.py:14: ShimWarning: The `IPython.html` package has been deprecated. You should import from `notebook` instead. `IPython.html.widgets` has moved to `ipywidgets`.\n",
      "  \"`IPython.html.widgets` has moved to `ipywidgets`.\", ShimWarning)\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import sys, os\n",
    "sys.path.insert(0, os.path.expanduser('~/git/github/taku-y/pymc3'))\n",
    "\n",
    "import numpy as np\n",
    "from scipy import stats\n",
    "import pymc3 as pm\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data = np.random.randn(100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Applied log-transform to sd and added transformed sd_log to model.\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as model: \n",
    "    mu = pm.Normal('mu', mu=0, sd=1, testval=0)\n",
    "    sd = pm.HalfNormal('sd', sd=1)\n",
    "    n = pm.Normal('n', mu=mu, sd=sd, observed=data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sampling function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import theano\n",
    "import theano.tensor as tt\n",
    "from theano.sandbox.rng_mrg import MRG_RandomStreams\n",
    "from pymc3.variational.advi import ADVIFit\n",
    "from pymc3.sampling import NDArray, _choose_backend\n",
    "from pymc3.backends.base import MultiTrace\n",
    "\n",
    "def sample_vp(draws, model, vparams, seed=1):\n",
    "    \"\"\"Draw samples from variational posterior. \n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    draws : int\n",
    "        Number of random samples. \n",
    "    model : pymc3.Model\n",
    "        Probabilistic model. \n",
    "    vparams : dict or pymc3.variational.ADVIFit\n",
    "        Estimated variational parameters of the model. \n",
    "    seed : int\n",
    "        Seed of random number generator. \n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    trace : pymc3.backends.base.MultiTrace\n",
    "        Samples drawn from the variational posterior. \n",
    "    \"\"\"\n",
    "    if type(vparams) is ADVIFit:\n",
    "        vparams = {\n",
    "            'means': vparams.means, \n",
    "            'stds': vparams.stds\n",
    "        }\n",
    "\n",
    "    # Make dict for replacements of random variables\n",
    "    r = MRG_RandomStreams(seed=seed)\n",
    "    updates = {}\n",
    "    for var in model.free_RVs:\n",
    "        u = theano.shared(vparams['means'][str(var)]).ravel()\n",
    "        w = theano.shared(vparams['stds'][str(var)]).ravel()\n",
    "        n = r.normal(size=u.tag.test_value.shape)\n",
    "        updates.update({var: (n * w + u).reshape(var.tag.test_value.shape)})\n",
    "    vars = model.free_RVs\n",
    "        \n",
    "    # Replace some nodes of the graph with variational distributions\n",
    "    samples = theano.clone(vars, updates)\n",
    "    f = theano.function([], samples)\n",
    "\n",
    "    varnames = [str(var) for var in model.unobserved_RVs]\n",
    "    trace = NDArray(model=model, vars=model.unobserved_RVs)\n",
    "    trace.setup(draws=draws, chain=0)\n",
    "\n",
    "    for i in range(draws):\n",
    "        # 'point' is like {'var1': np.array(0.1), 'var2': np.array(0.2), ...}\n",
    "        point = {varname: value for varname, value in zip(varnames, f())}\n",
    "        trace.record(point)\n",
    "\n",
    "    return MultiTrace([trace])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inference with NUTS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\r",
      " [-----------------100%-----------------] 500 of 500 complete in 0.2 sec"
     ]
    }
   ],
   "source": [
    "with model:\n",
    "    step = pm.NUTS()\n",
    "    trace_nuts = pm.sample(500, step)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inference with ADVI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 0 [0%]: ELBO = -215.01\n",
      "Iteration 2000 [10%]: ELBO = -343.54\n",
      "Iteration 4000 [20%]: ELBO = -159.81\n",
      "Iteration 6000 [30%]: ELBO = -156.86\n",
      "Iteration 8000 [40%]: ELBO = -151.9\n",
      "Iteration 10000 [50%]: ELBO = -148.11\n",
      "Iteration 12000 [60%]: ELBO = -148.88\n",
      "Iteration 14000 [70%]: ELBO = -153.17\n",
      "Iteration 16000 [80%]: ELBO = -149.75\n",
      "Iteration 18000 [90%]: ELBO = -150.09\n",
      "Finished [100%]: ELBO = -149.35\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(-1000, 0)"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAg0AAAFkCAYAAACjCwibAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4U2X+9/FP0jTd0p0ChZZ9X1qhKEoBFXUUAdcCwjhu\nuIDWcQHFBREcwPmNyvwerc4zM27MKANSxVHnmUUdBQE30AKKZRGQsogFCjQptGlznj8KgdjtAG2S\npu/XdXFdzVnu8z05LfnknPvcx2IYhiEAAIAGWANdAAAAaB4IDQAAwBRCAwAAMIXQAAAATCE0AAAA\nUwgNAADAFFugCzgdhmFo1qxZ2rhxo+x2u+bOnav09PRAlwUAQEhrlmcaPvjgA1VUVGjRokWaOnWq\nnnzyyUCXBABAyGuWoWHNmjUaNmyYJCkzM1PffPNNgCsCACD0NcvQ4HQ6FRsb631ts9nk8XgCWBEA\nAKGvWYYGh8Mhl8vlfe3xeGS11r8rjJYNAMCZaZYdIQcOHKiPPvpIl112mQoKCtSjR48G17FYLCou\nLvVDdWhsKSmxHLtmjOPXfHHsmreUlNiGFzpFzTI0XHLJJVq5cqWuu+46SaIjJAAAftAsQ4PFYtHs\n2bMDXQYAAC1Ks+zTAAAA/I/QAAAATCE0AAAAUwgNAADAFEIDAAAwhdAAAABMITQAAABTCA0AAMAU\nQgMAADCF0AAAAEwhNAAAAFMIDQAAwBRCAwAAMIXQAAAATCE0AAAAUwgNAADAFEIDAAAwhdAAAABM\nITQAAABTCA0AAMAUQgMAADCF0AAAAEwhNAAAAFMIDQAAwBRCAwAAMIXQAAAATCE0AAAAUwgNAADA\nFEIDAAAwhdAAAABMITQAAABTCA0AAMAUQgMAADCF0AAAAEwhNAAAAFMIDQAAwBRCAwAAMIXQAAAA\nTCE0AAAAU1pMaDhSXqn1W/fLMAxJ0tufbNXU51eqwl2l3ftcevWf3+ntT7Zqw/YDKi2rUP7H36u0\nrEKHyyp0y2//q3mvrfG2VVnlUWWVRx99tVNvLvteVR6Pz7zDrgrv64POch1yVWhT0UE5j7hVUlqu\nZ/PXad/BI6qs8ujb7Qe0cUeJNhUdlOdYbZL0+Ya92ltS5q23ssqjzzfs1fe7D2lT0UFt23NYm3ce\nVMGWfdqxt1TuyuoaNu4oUdnRSp/9PnD4qL7eXCzDMLTv4BEV/eT0ziuvqNKWXYckSSWl5Tpw+Kj2\n7Hfpo693qcJdpSPlldqz36Vtew5r6+7D+vFAmTweQwWb96nwhxLtPVDmrW/H3lJVeTz66eAR7d7n\nkmEY2rPfpXJ3lX4qKZO7sspbV9lRt5xH3D77d/yfJO97YRiGKqs82nfwiKo8Hnk8J94jj8dQhbvK\nu55hGDpSXinDMOQxDO+6xxknvb+GYcjjMXy2c/xfbcsCACSL0UL+R7z5iX9r36GjNaZ3S4vXlp2H\nfKZl9UzRmo3FNZYNt1m9H84/FxcdrrgYu3YWu7zTuqfFa/PP2jbjnpwM/Z/8dae8XnprhzcQSFJq\ncrT27C875XYANK0ObRzasdfZ8IIBkBQX4f3i0zY5Rjt/cio1OVoHSsvldnu8QTstJUY7i13q0i5O\nPTsk6L9f7ZLdZlV5RZUsVovKK6rUOTVOew+Uqay8+otMfIxdkfYwJcVFKjY6XFt2HdKR8ipVuKuU\nGBuh2Gi7tu05LEnq2DbWu52V639Um6RoWVT93pWWueWu9Oigs1ztW8Vo7ff7lRgbodYJUdq936XS\nMrdaxUdq36GjGp6Zqm+2HZDzSPU6WT1bKzoiTMUHjyrSHqYqj6Gyo5Wyh1v1w4+lSoiNkMdjqLTM\nraMVVeqUGqtIe5jKK6rkPOJWemuHvtq0T+mtY7T/cLlcR9zq3SlRhT8cVL/OSWqTFKXig0e1cUeJ\nYqPtKndXKatniiqrDG3fc1hhVosSYiMUZrXIXenR/sNH1b5VjEpKy9U+xaGNOw7KapWqPIZcR9yK\nd0Tohx9LNahnisLCrFr3/X6VV1TJYxjqmZ6gIxWVKiktV0pClKo8hpxlFerQJlZVHkNz7xza6L8f\nLSY0jJn690CXAACA37z7zJWN3maLuTwBAADODKEBAACYYgt0AQAQ6m4b3Ud/fm/DKa83qGeK0tvE\naunyrXrkV1kqdVXo2+0HdMhZoTWbfPtdjTy3g/752Y4628rqmaJLz+7g7dQ96ryO6pwapw9WF6lw\nx0GfZYdmpGrFuj2SpKemDNHytbt10FmuXh0TZbeF6WhFpb747iddPChN67/frw/W7NRdV/dXVs8U\nbdxRov9Z+LV3m307J8lZ5tbIczsozFr9PbW0rEKvv79JX3z3kyTp6uFdNGZIJ1VWefTjgTIlx0Vq\n1itf6KqhXdQtLV5Ll2/VV5uK5TEMVVYZmnxlX2V2a6Wtuw6p9IhbtjCr8t5ar/vGZSolIUqHXRVa\n9c2PuuHSnrr1dx9Jkn435Ty9/ck2XZHdSa0To2UYhvLeWq8+nZJ0UVaaz/6XlJbroLNcaSkxOlJe\npe0/lmrr7kM6UFqumy7rpY8Ldunjr3drZ7FTD04YoFYJkTpaXqUf9paqe3qCVq3fo3dWbldaSowG\n92mjj7/epVtH99Gygt36bMNe/T43W7JYFGa16J0V21S446ByLuii/12yTrYwq/70wAWSqjuI28Ks\nOuSq0Hurtuvycztqw/YDeukf33lrffrOIUqKi9TmnQf15GtfVb/vPVJ06+g+5n/RTgF9GgC0GI6o\ncDmPuL2v7x2boR9+LNWqb37U3pIjkqRbLu+t5Wt3e+8qOm76xAH6ZtsBFf3k1Lrv93unX3pOulon\nRis+xq4ER4SWr92l5Wv3eOfPuXWw2rWK0Ydrdur19zdJku6+pr9WrN+jvp2TlN7aoSdf+0rn9G6t\nrJ6t9ad3vlXVsbuEsvu31aRRtf/n/8V3e/V///6tJOn5+4YrKsKmf372g5Z8/L2uyO6kof1T9fr7\nm5TexqHvdx3WpFG9lRQXWe/7M/mZj2W3henZe4ZJklJSYlVcXFrvOoZh6JCrQgmOCO+0w2UVio0K\nl8ViqXddqfpOKauJ5U7Xog83a//ho7rr6v5Nto3GUnbUrbAwqyLCw+pdbv3W/fr9G2uV1TPFZ7+e\nf2u91mwq1g2X9dQFZ7VXSkpso9dIaAAQMFdkd9I7K7crLjpch8uqP8yvG9FNi/67xdT6XdvH6YHr\nBqii0qOoiDB5PIbueHpZrcsmx0XoqTuzdaS8Ui++t0GDerXWeX3bSqr+Zvn4y18oJSFSj94wSH97\nf7M+/GqnJKlzaqy27SnVnx+8wPtN+c1l36tr+3id1a1Vrdt6Z+U2vf3JNknSyw+NaHA/jpRXKtIe\n5v2Qnfr8SpWUlmto/1TdMqp3retUVnn04nsbNGJgmnqkJ0iqvg156+7D6twu1lvrqais8shikXdd\nM6EBgbFjb6lSk2MUbjtxnI+UV+rbbQc0sEeKrFZLk4QGLk8A8LsJF3XXRYPSZLVYdMXQzrJImjJ/\nmSrcHv3inA6yWCxqkxSttsnRWrtln9747xZVeQzFRNp01bAu3m/sj/5qkCTJfuybmcdT+y3Rd1zR\n1/vBGhVh093XZvjMT4yN8H67liSPqr9LOaLC9egNg+R2e3w+hK89v2vjvBHHREXU8V9xPV/AbWFW\nTb6yn880q9Wibmnxp12HLYxubs1FhzY1A0FUhE2DerVu0u3yGwK0QD3TEzTqvI6nvF6XdnGSpCuH\ndtbl5/quf+vomt+IWydE1dpOvy5J3lPSVotFFotFefcO1x+nnS9JuuTsdGV0TVbrhChdMihdF5zV\n3rvuhQPb19pmfQb3aaPE2IiGFzzm+PlXq6W6vgh7/aeLazZwaovX3H51A0130h44PZxpAFqgq4Z1\nVs8OifrHpz9IklrFR6ptcrS+2XpAibERKiktlyTdPqaPUpNjNPvVLyVJU67sJ6vVogSHXZVVHv2/\nz6rXT0uJ0ZB+qeqelqAwq0XTXliltJQYVbhr/+Zf20XR+r7ltm8dI0nq1TFRVotF40d08xlB9bhw\nm1XXXdRdiz7cbP7NqMUlg9L0+Ya9dV4aaMiZXvO9/hc9lffW+hod9IBAIzQAIezCge310Ve7JElx\nMXbvSH/HT9Wf17etPv32R0nS/ePOksdj6A9vf6M1pdU989PbxKp9qxhvexaLvN/Yw21h6tg2Vj/8\nWOrtYJdy7MzC73OzFRMVroLN+/TC29/UqOtUu1INz2in+Gi7enVMlCRdek6HOpf9xdnp3tBw37jM\nBjv/1SY1OUbP3zf8lNc7Li2l+j3L6Jp8WusP7JFiqi8E4G9cngBCWJT9xPeCIcc6/dlt1jp7tVut\nJ6ZH2MN8AoOkmuvV8dkf74iQLcyqQb1a69mpF9SYf6rfxK1Wiwb0SKn72v/P9OuSpJGDO6h/l+Qa\n++APA3ukaOr4szT5yr5+3zbQlAgNQID976+H6n9/7TtGfGpy9Cm30yo+Uuef1a7WeWFW81fHRxw7\nJX7zyF4NLjssM1WSdG6fNqbbl050XGwq9487S2Mv7Nak26iPxWJR385JirRzMhehhd9oIEA6p8bp\nsRure/+f/GRUSRqW0U4d2jj03qrt3oF35t42WHtLjujZYw8z+93k85QcH6lJ/1M9eE16a4fGXtBN\nywp2S6oneJyUH2o74dC7Y6JefPBCn7MOdS0/YmCazu7VWrHR9gb3V5Jm33KOtuw8WGcHSQDBLShC\ng9Pp1LRp0+RyueR2u/Xwww8rMzNTBQUFmjdvnmw2m4YMGaLc3FxJUl5enpYtWyabzaaHH35YGRkZ\nDWwBCLzhmak+g/6c/AF88r3WknTxoDTZwqzq0ylJO/aWyhEVrqS4SKUmnzjV3urYB+/xpw1KUnSk\nTTNuGCRHlE2t4qO09JOt3uWH9G+rf32xQ7+8pId32nl922rVNz/qssG+fQRqCwxS7b35GwoMrROr\nw0tG12Slt3YovbWj3uUBBK+gCA2vvPKKhgwZohtuuEHbtm3T1KlT9dZbb2nWrFnKy8tTWlqabr/9\ndhUWFsrj8Wj16tVasmSJ9uzZo7vvvlv5+fmB3gWgVgO6t9JFWdWD79jCrD6h4WRRETbdOzZTh5zl\n6pwa53MnQW33Y5/seP+A4/0Njt8W+XNpKQ69NP1Cn34JfTsn6Q9Tz29wBDqv0xi5LyYqXC/cP9z8\nNgAEraAIDTfffLPs9upvK5WVlYqIiJDT6ZTb7VZaWvX11aFDh2rlypWy2+3Kzs6WJKWmpsrj8aik\npESJiYkBqx+hb9Ko3j7jvZsVGx2uPp2SvK8fu3GQfrNgda3Lnm5P++OpwczHeW0dIE/lw/x0xw3g\n2j4QGvz+l5yfn68FCxb4THvyySfVr18/FRcX68EHH9Sjjz4ql8slh+PEacyYmBgVFRUpMjJSCQkJ\n3unR0dFyOp2EBjSpy4d1Pa3QEBlp9xnKtfrn6tBgs1nPaJjX4+uGHbu0ERFpq9Fe9EmXDhpjSNlW\nrRyKd5gfJKkxt43A4NjhZH4PDTk5OcrJyakxfePGjZo2bZqmT5+uQYMGyel0yul0eue7XC7Fx8cr\nPDxcLpfLZ3psLL/UaFr79p3e+PtHj7rrHLvf7fac1rj+943LVJTd5l23qrJ6AKXy8soa7ZWVnehg\n2RjPENi/36mKIxUNL3gSnl/QfHHsmremCHxBccvlli1bdO+99+rpp5/W0KHVt545HA7Z7XYVFRXJ\nMAytWLFCWVlZGjBggFasWCHDMLR7924ZhuFz5gEItME+tx82/vPg+ndJ9nm+gD+fOGfmqYUAQldQ\nXGicP3++KioqNHfuXBmGobi4OD3//POaNWuWpk2bJo/Ho+zsbO9dEllZWRo/frwMw9DMmTMDXD1w\ngiMqXDnnd9XnG/b6bZv+eE6B3WZVRaVHdltQfM8AECBBERpeeOGFWqdnZmZq8eLFNabn5uZ6b78E\nmqurh3du3Aab8CzA03dl67CroskHZQIQ3IIiNAChJCbK5FDHnU/zbokAcESFyxEVHugyAAQY5xqB\nYzo2MB5CXVrFn3ggksXy89sLm74PgHEKt1wCwJkgNADHPH7z2fXMrWOERIs05ap+3tdtk6pHP4x3\nmBtWuTGcGNyp5rzjQ0n36kBnYQBnjssTQAOu/0WPGtNaJ0SpTVK0xmR3UufUOP3fqefr/dVFyu6f\n6vf6hvRrq6XLt+qsbq1qzDu3b1uF28LUtxPjmAA4c4QGoAEjBqbJfWwshOOiIm26b1ym97U9PEyj\nzut0YgE/3gc5+ryOOq9vG7WKr/kQKKvForN7tfZfMQBCGpcn0CINz6z9EdKNzR/DGlgslloDAwA0\nNkIDWqQz/TA3u7rhz5GXAKCJERrQIjX5CQBuZQAQgujTAJxkfm62qqoMVXk8euiPnzXJNi7KSlN4\nGHkdQPNDaECLM/OmQVpesLvWeQkmn+DY4OWNei5L/PKSmndjAEBzwNcdtDid2sb5p4ei/LYZAPAL\nQgNwWkgDAFoeQgNahHP7tml4IQBAvQgNaBEmjert89rMeYJu7eOV1TOlenlOLAAAHSHRQpkIAY/8\nKuu0m2d4BgChiDMNCCnjLuxmarl+nZJOqd2fD9LEmQcALRGhASHlssEdTC03oEfKGW2nocxApgAQ\niggNaBEsfIwDwBkjNCDkTLvurKbfSAMZhD4NAEIRoQEhp0+nJL00/UJdPaxzoEvh/AaAkEJoQEiy\nWCyyNerzHXzPHURHhDdi2wDQPBAaEPR+Mbjjaa3XVJcIsvu11Y2X9Wyi1gEgeBEaEPTuzMnUvNvP\nPbNGGvE6waTRfZQUF9l4DQJAM0FoQNALs1rUNin6lNczfj64whmhdwIAEBoAU7gfAgAIDQAAwBRC\nAwAAMIXQAL/o3TGxSdqNiqj7mWsnd2nwe4+ERu1PAQDBgdAAv7h3bGaTtHvb6D5N0u7PnXYG4MlW\nAEIIoQF+EW7z/68a3/UBoHERGtCsGUQDAPAbQgMAADCF0ICQMefWwYEuwWtgjxRJUtd2cQGuBAAa\nT91dz4FG0DM9QRuLDvplW+1axdQ5z1JLh8TU5Gjt2V9mqu1wm1Wt4iPVv2uyqeUnXtJD2Rmp6pJK\naAAQOggNaFIPThxwSncePPTLgfrHpz9o/db9Z77xBjY888az9cZHW/TR17sabMpiseh3U4aY3rQt\nzKqu7eJNLw8AzQGXJ9CkLBaLrFbztx0mxkac2pgOZ9APMsIedlrPtACAlorQgKAzLDO1UdrhvgoA\naFyEBgSdmMhwPXpDlrmFGTsJAPyGPg0IuE5tY7X9x1KfaY6o8DNv+GenGh75VZbiY+xn3i4AtFCE\nBgSlNomN39egW3s6JgLAmeDyBJq3ejou0KcBABoXoQEB17+LubEPAACBRWhAwF0xtFOgSwAAmEBo\nQMCFWU/8GtrDwwJYCQCgPnSERFD4za2DVfRTqf/vbuCWTQAwjdCAoNC+VYza1/PsiLqEhdV9ssww\nM341vSUBwDRCA5q1vp0TNbR/qrL7tw10KQAQ8ggNaNbCrFbdMqp3oMsAgBaBjpAAAMAUQgNC1qk8\nkhsA0DBCA0IeN0gAQOMgNAAAAFMIDQhZXJ0AgMYVVKHh+++/16BBg1RRUSFJKigo0Lhx4zRx4kTl\n5eV5l8vLy9PYsWM1YcIErVu3LlDlorng+gQANIqgCQ1Op1O/+93vFBER4Z02a9YszZ8/XwsXLtS6\ndetUWFioDRs2aPXq1VqyZInmz5+vJ554IoBVIxCuPb+LqeX6d0mSJF02uENTlgMALUbQhIaZM2fq\n/vvvV2RkpKTqEOF2u5WWliZJGjp0qFauXKk1a9YoOztbkpSamiqPx6OSkpKA1Y3g1T0tQc/eM0w5\n53cNdCkAEBL8PrhTfn6+FixY4DOtXbt2GjVqlHr27Okd+tflcsnhcHiXiYmJUVFRkSIjI5WQkOCd\nHh0dLafTqcTERP/sAJoVR1R4oEsAgJDh99CQk5OjnJwcn2mXXnqp8vPztWTJEu3bt0+TJk3SH/7w\nBzmdTu8yLpdL8fHxCg8Pl8vl8pkeGxvrt/pxalJSGj42p7pMTExErdNPh8PReG2hbry3zRfHDicL\nimGk//3vf3t/HjFihF5++WWFh4fLbrerqKhIaWlpWrFihXJzcxUWFqann35at9xyi/bs2SPDMHzO\nPCC4FBeXNvoyLlf5Ka1bH6ez8dpC7VJSYnlvmymOXfPWFIEvKELDySwWi/cSxezZszVt2jR5PB5l\nZ2crIyNDkpSVlaXx48fLMAzNnDkzkOUiABp1pEfurAAA04IuNHz44YfenzMyMrR48eIay+Tm5io3\nN9efZSFUMZgDAJgWNHdPAACA4EZoAAAAphAa0OxY6IcAAAFBaAAAAKYQGgAAgCmEBgAAYAqhAQAA\nmEJoAAAAphAaAACAKUE3IiRajjm3DtYhV0WgywAAmERoQMC0axWjdq1iAl0GAMAkLk8AAABTCA0A\nAMAUQgNatMTYCElS26ToAFcCAMGPPg1o0Qb2TNFNI3spo2tyoEsBgKDXYGioqKjQv/71L61fv16S\n1L9/f1122WWy2+1NXhzQ1KwWi4Zntgt0GQDQLNR7eaKkpETXXnut/vrXv8pms8kwDP3lL3/Rtdde\nq5KSEn/VCEiS7r6mvy4b3EExUeGBLgUAWqR6zzQ89dRTGjNmjG6//Xaf6S+88IKeeuopzZs3r0mL\nA042oEeKBvRI0ccFuwJdCgC0SPWeaVi/fn2NwCBJd955p9asWdNkRQEAgOBTb2hwu911zgsLC2v0\nYgAAQPCqNzS0adNGn332WY3pn376qVJTU5usKAAAEHzq7dMwdepU3XnnnbruuuuUkZGhqqoqff31\n13rrrbf04osv+qtGhIjBfdro8w17z7wh48ybAACcunpDQ0ZGhl599VW99NJL+te//iWLxaKMjAwt\nXLhQHTp08FeNCBG3j+mj8SO66f68lYEuBQBwGhocp6Fbt2568skn/VELQpzFYpGD2yUBoNmqt09D\nRUWF/va3v+mDDz6Q0+nUpEmTNHDgQP3qV7/Stm3b/FUjAAAIAvWGhkceeUSrVq3SokWLdP3116t3\n795auHChRowYoZkzZ/qrRgAAEATqvTxRWFio9957TxUVFRo+fLimTZsmSerVq5fefPNNvxSI4DH2\nwq5a8tH3gS4DABAg9Z5psNmqM4Xdblfbtm1rnYeWYfrEAXXOG9QzxY+VAAACpd7QYLFYav25ttcI\nbT07JNY5b+S5Hf1YCQAgUOo9XfDdd9+pd+/ekiTDMHx+JjQAANCyNNinAQAAQGrg8kR9xowZ05h1\nAACAIHfaoWHnzp2NWQdgGqNIA0BgnHZooE8DAAAty2mHBgAA0LLU2xGyV69etZ5R4O4JAABannpD\nw2OPPaZf/vKXkqRNmzapR48e3nlz5sxp2soAAEBQqffyRH5+vvfn6dOn+8xbs2ZN01QEAACCUr2h\nwTCMWn+u7TUgSZNG9Q50CQCAJmK6IyTDSMOM7P6pgS4BANBETD97AmgMnKACgOar3o6Qmzdv1kUX\nXSRJ2rt3r/dnwzBUXFzc9NXhjMXF2HXYVRHoMgAAIaDe0PDvf//bX3UApnH+CwACo97Q0L59e3/V\ngSA297bBkiRLkHxcc4UDAAKDESHRoNTkmEZsjY98AGiuCA0AAMAUQgNMMzhLAAAtGqEBAACYQmgA\nAACmEBrgVwzuBADNF6EBAACYQmgAAACmEBoAAIApQREaPB6P5s6dq4kTJyonJ0fLli2TJBUUFGjc\nuHGaOHGi8vLyvMvn5eVp7NixmjBhgtatWxeosgEAaFHqHUbaX/7+97+rqqpKCxcu1N69e73PvJg1\na5by8vKUlpam22+/XYWFhfJ4PFq9erWWLFmiPXv26O6771Z+fn6A9wBmWa3BMRQ1AODUBUVoWLFi\nhbp376477rhDkjRjxgw5nU653W6lpaVJkoYOHaqVK1fKbrcrOztbkpSamiqPx6OSkhIlJiYGrH6Y\nZwuz6q6r+6t1YtTpN8ItGAAQEH4PDfn5+VqwYIHPtKSkJEVEROiPf/yjvvzySz388MN65pln5HA4\nvMvExMSoqKhIkZGRSkhI8E6Pjo6W0+kkNDQjWT1TAl0CAOA0+D005OTkKCcnx2fa/fffrwsvvFCS\ndPbZZ2v79u1yOBxyOp3eZVwul+Lj4xUeHi6Xy+UzPTY21j/FN0ONcTkgJaX6/XXERNQ6PzExusay\ndb0+0xokyREb2ejto2lxnJovjh1OFhSXJ7KysrRs2TJdcsklKiwsVLt27RQTEyO73a6ioiKlpaVp\nxYoVys3NVVhYmJ5++mndcsst2rNnjwzD8DnzAF8ez5mfyi8uLpUkOV3ltc4vKSmrsWxdr8+0Bkly\nlh5t9PbRdFJSYjlOzRTHrnlrisAXFKFh7NixmjVrlsaPHy9Jmj17tqTqjpDTpk2Tx+NRdna2MjIy\nJFWHjPHjx8swDM2cOTNgdbcEl56THugSAABBIihCg91u17x582pMz8zM1OLFi2tMz83NVW5urj9K\na/Gy+6cGugQAQJAIinEaELzat4oJdAkAgCBBaECdurWPl8XCuAoAgGqEBpwxcgUAtAyEBgAAYAqh\nAQAAmEJoQKO7dXTvJm2fQaQBIDAIDWh0FtHJAQBCEaEBAACYQmhoAa7I7hToEgAAIYDQ0AJcNaxL\noEsAAIQAQgMAADCF0AAAAEwhNAAAAFMIDQAAwBRCAwAAMIXQEOIac5glBm0CgJaN0ADTjCAZwNkI\njjIAoMUhNAAAAFMIDQAAwBRCAwAAMIXQgDPy+E1nB7oEAICfEBpwRjq2jfX7Ni3cxAEAAUFoAAAA\nphAaAACAKYQGNLpgGc8BANC4CA3wERVhU1yM/ZTWYaRIAGgZCA0h7lS/8/dIi1daSkyT1AIAaN4I\nDWh2GEYaAAKD0IBGx+UKAAhNhAYAAGAKoQEAAJhCaAAAAKYQGuDDYrHQ0RAAUCtCA+pGf0YAwEkI\nDagbZxwAACchNKAGniIJAKgNoQEAAJhCaIAPg16QAIA6EBoAAIAphAYAAGAKoQEAAJhCaAAAAKYQ\nGgAAgCmKKx7WAAAOdklEQVSEBjQ6g1GhACAkERoAAIAphAYAAGAKoQEAAJhCaIBpFh57CQAtGqGh\nhXNEheuK7E5Kb+1ocNnhme3Uu2OiHpwwwA+VAQCCDaGhhWvXKkZXDesiq4lHW0ZH2vTAhAHq1THR\nD5XVjedjAEBgEBoAAIAphAYAAGAKoQEAAJhCaAAAAKbYAl2AJDmdTt13330qKytTRESEnnrqKSUn\nJ6ugoEDz5s2TzWbTkCFDlJubK0nKy8vTsmXLZLPZ9PDDDysjIyPAewAAQOgLijMNb731lnr27KnX\nX39dI0eO1IsvvihJmjVrlubPn6+FCxdq3bp1Kiws1IYNG7R69WotWbJE8+fP1xNPPBHg6gEAaBmC\n4kxDjx49tHXrVknVZx3Cw8PldDrldruVlpYmSRo6dKhWrlwpu92u7OxsSVJqaqo8Ho9KSkqUmBjY\n2wABAAh1fg8N+fn5WrBggc+0mTNnauXKlRo1apQOHTqkhQsXyuVyyeE4MeBQTEyMioqKFBkZqYSE\nBO/06OhoOZ1OQkMdrFaLUlJi65wfHh6mlJRY2cKrTzrZI2yyuD0+8+qTkhKrw+VVPq9jYw/6vG4M\nJ7fjcEQ2evtoWhyn5otjh5P5PTTk5OQoJyfHZ9rdd9+t2267TePGjdPGjRuVm5urhQsXyul0epdx\nuVyKj49XeHi4XC6Xz/TYWH6p6+LxGCouLq1zvttdpeLiUlUeCwoV5ZVyV3l85tWnuLhUJSVlPq9L\nS4/4vG4MJ7fjdB5t9PbRdFJSYjlOzRTHrnlrisAXFH0a4uPjvWcVkpKSvGcZ7Ha7ioqKZBiGVqxY\noaysLA0YMEArVqyQYRjavXu3DMPwOfMAAACaRlD0afj1r3+tGTNmaOHChaqsrNScOXMkVXeEnDZt\nmjwej7Kzs713SWRlZWn8+PEyDEMzZ84MZOkIAKuVB2cBQCAERWho3bq1/vSnP9WYnpmZqcWLF9eY\nnpub6739EsGnX5dkWS0Wjb+oW5O0P6RfWxVs2afLB3dskvYBALULitCA0BIXbdeL0y9ssvYj7Tbd\nP+6sJmsfAFC7oOjTAAAAgh+hAQAAmEJoAAAAphAaAACAKYSGEMfNiQCAxkJoAAAAphAaWjrDCHQF\nAIBmgtAAAABMITQAAABTCA0h5Ly+bWpMa/Dig+XMu0o2QhMAgGaA0AAAAEwhNAAAAFMIDZAkGQ1f\nyAAAtHCEBviw0EEBAFAHQgMAADCF0AAfBoM9AQDqQGiAJMnCUyoAAA0gNAAAAFMIDQAAwBRCQwg5\nr1/bQJcAAAhhhIYQkRQXoX6dkwNdBgAghBEaQkR4GIcSANC0+KQBAACmEBoAAIAphAYAAGAKoQEA\nAJhCaAAAAKYQGgAAgCmEBgAAYAqhAQAAmEJoAAAAphAaAACAKYQGAABgCqEBAACYQmgAAACmEBoA\nAIAphAYAAGAKoQEAAJhCaAAAAKYQGgAAgCmEBgAAYAqhAQAAmEJoAAAAphAagpgl0AUAAHASQkMQ\nO7dvm0CXAACAF6EBAACYQmhoYdJSHIEuAQDQTBEaWhgLHSUAAKeJ0AAAAEwhNAAAAFMIDQAAwBRC\nAwAAMCVgoeH999/X1KlTva/Xrl2rcePGaeLEicrLy/NOz8vL09ixYzVhwgStW7dOklRSUqJJkybp\n+uuv1/3336/y8nK/1w8AQEsTkNAwd+5c/f73v/eZ9vjjj2v+/PlauHCh1q1bp8LCQm3YsEGrV6/W\nkiVLNH/+fD3xxBOSpOeff15jxozRa6+9pl69eulvf/tbIHYDAIAWJSChYeDAgZo1a5b3tdPplNvt\nVlpamiRp6NChWrlypdasWaPs7GxJUmpqqjwejw4cOKCvvvpKw4YNkyQNHz5cn332md/3AQCAlsbW\nlI3n5+drwYIFPtOefPJJjRw5Ul988YV3msvlksNxYtChmJgYFRUVKTIyUgkJCT7TnU6nXC6XYmNj\nvdNKS0ubcjeaBaPOGXXOqXW+4fNzA+sCAFqUJg0NOTk5ysnJaXC542HgOJfLpfj4eIWHh8vlcnmn\nO51OxcXFeZdPSkryCRD1efeZK09vJ5qZuvazof3Pe3DEaW8nJSW2Sd7fn28DzRfHr/ni2OFkQXH3\nhMPhkN1uV1FRkQzD0IoVK5SVlaUBAwZoxYoVMgxDu3fvlmEYSkhI0MCBA7V8+XJJ0vLlyzVo0KAA\n7wEAAKGvSc80nIrZs2dr2rRp8ng8ys7OVkZGhiQpKytL48ePl2EYmjlzpiRpypQpmj59ut544w0l\nJibqmWeeCWTpAAC0CBbDaOiiNwAAQJBcngAAAMGP0AAAAEwhNAAAAFMIDQAAwJSguXuiKRiGoVmz\nZmnjxo2y2+2aO3eu0tPTA10Wjrnmmmu8g3qlpaVp8uTJeuihh2S1WtW9e3c9/vjjkqQ33nhDixcv\nVnh4uCZPnqwLLrhA5eXleuCBB7R//345HA799re/VWJiYiB3p8VYu3atnn76af31r3/Vjh07zviY\nFRQUaN68ebLZbBoyZIhyc3MDvIeh6+Rj99133+mOO+5Qp06dJEkTJkzQyJEjOXZBqLKyUo888oh2\n7dolt9utyZMnq1u3boH52zNC2H/+8x/joYceMgzDMAoKCowpU6YEuCIcV15eblx99dU+0yZPnmx8\n+eWXhmEYxsyZM43333/fKC4uNkaPHm243W6jtLTUGD16tFFRUWG88sorxnPPPWcYhmH84x//MObM\nmeP3fWiJ/vznPxujR482xo8fbxhG4xyzK6+80igqKjIMwzBuu+0247vvvgvAnoW+nx+7N954w3jl\nlVd8luHYBac333zTmDdvnmEYhnHo0CHjggsuCNjfXkhfnlizZo33GRWZmZn65ptvAlwRjissLFRZ\nWZkmTZqkm266SWvXrtWGDRu8A3UNHz5cq1at0rp165SVlSWbzSaHw6FOnTqpsLBQa9as0fDhw73L\nfvrpp4HcnRajY8eOev75572vv/3229M+Zp999lmtz51ZtWqV/3esBajt2H388ce6/vrrNWPGDLlc\nLo5dkBo5cqTuueceSVJVVZXCwsLO6P/LMzl+IR0anE6nzxDTNptNHo8ngBXhuMjISE2aNEkvvfSS\nZs2apWnTpsk4aciQ2p4zIknR0dHe6ccvbfx8GHI0nUsuuURhYWHe12dyzEpLS2t97gzPkmkaPz92\nmZmZevDBB/Xaa68pPT1deXl5Nf7P5NgFh6ioKO+xuOeee3TfffcF7G8vpEODw+HweXaFx+OR1RrS\nu9xsdOrUSVdccYX354SEBO3fv9873+VyKS4uTg6Ho8ZzSY5PP35szT5/BI3v5L+n0zlmtT13Ji4u\nzn870IJdfPHF6tOnj/fnwsJCxcbGcuyC1J49e3TjjTfq6quv1qhRowL2txfSn6ADBw7UsmXLJEkF\nBQXq0aNHgCvCcW+++aZ++9vfSpL27t0rp9Op7Oxs79NPly9frqysLPXv319r1qxRRUWFSktLtXXr\nVnXv3l0DBgzwHttly5bx/JEA6dOnj7788ktJp3fM6nruDJrepEmTtH79eknSp59+qr59+3LsgtS+\nffs0adIkPfDAA7r66qslSb179w7I315IDyNtnHT3hFT9WO7OnTsHuCpIktvt1sMPP6zdu3fLarXq\ngQceUEJCgmbMmCG3262uXbtqzpw5slgsWrJkiRYvXizDMDRlyhRdfPHFOnr0qKZPn67i4mLZ7XY9\n88wzSk5ODvRutQi7du3S1KlTtWjRIm3fvl2PPfbYGR2zdevWae7cud7nztx7772B3sWQdfKx27Bh\ng37zm98oPDxcKSkpeuKJJxQTE8OxC0Jz587VP//5T3Xp0kWGYchisejRRx/VnDlz/P63F9KhAQAA\nNJ6QvjwBAAAaD6EBAACYQmgAAACmEBoAAIAphAYAAGAKoQEAAJhCaABQQ69evSRVD8V+1113NVq7\nN9xwg/fn44PUAGg+CA0AarBYLJKkgwcPqrCwsNHaPT7ipyQtXbq00doF4B+2QBcAIHjNnTtXP/30\nk+6++24999xzevvtt/WXv/xFhmGob9++mjlzpux2u84991z169dP+/fv15IlSzR79mxt3rxZ+/fv\nV+fOnfXcc8/pqaeekiSNHz9eixcvVq9evVRYWKijR49qxowZ2rhxo6xWq26++WZdddVVWrp0qT75\n5BMdOnRIRUVFys7O1uOPPx7gdwRo2TjTAKBOM2bMUOvWrfXcc89py5YtWrJkiRYtWqSlS5cqKSlJ\nL7/8sqTqMxKTJ0/W0qVLVVBQILvdrkWLFuk///mPjhw5ouXLl2vGjBmSpMWLF0s6cTbj2WefVWJi\not599129+uqrysvL06ZNmyRVPzMmLy9P77zzjj766CNt3rw5AO8CgOM40wDAlM8//1w//PCDxo8f\nL8MwVFlZqb59+3rnZ2RkSJIGDRqkhIQEvf7669q2bZt27Njh87TZ2tqdN2+eJCkxMVEXX3yxvvji\nC8XExGjAgAGKioqSJKWnp+vQoUNNuIcAGkJoAGBKVVWVRo4cqUcffVSSdOTIEVVVVUmqPmtgt9sl\nSR9++KGee+453XTTTbr22mtVUlJSb7s/f/yNx+NRZWWlJHnbrGtZAP7F5QkANRz/cLbZbN5gcM45\n5+iDDz7QgQMHZBiGHn/8cb366qs+y0vVj1m+/PLLddVVVykpKUlffvmltw2bzSaPx+OzzuDBg5Wf\nny9JOnDggD788EMNHjzYL/sJ4NQQGgDUcLy/QXJystq2basbb7xRvXr10p133qkbb7xRY8aMkWEY\nuv32232Wl6Rx48bp3Xff1TXXXKN77rlHZ511lnbu3ClJGjFihK688kpVVFR417nrrrt08OBBjRkz\nRjfccIOmTJmi3r1711kTgMDh0dgAAMAUzjQAAABTCA0AAMAUQgMAADCF0AAAAEwhNAAAAFMIDQAA\nwBRCAwAAMOX/A1qNdTZa3zEBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1206c9c88>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "advi_fit = pm.variational.advi(model=model, n=20000, accurate_elbo=False)\n",
    "trace_advi = sample_vp(500, model, advi_fit)\n",
    "plt.plot(advi_fit.elbo_vals)\n",
    "plt.ylabel('ELBO')\n",
    "plt.xlabel('Iteration')\n",
    "plt.ylim(-1000, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparison of posterior samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x11e94f358>"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2cAAAFFCAYAAABosJVxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xm0ZVdd9vvvand39jnVpipFEhBCJyYBQ+N7GUAoCJLQ\nKSQIGBhwI9gCIq2El1zEEdEXERRQkBevBKIXJa9EhQRChNgRRDEQQkiKtFR7+t2vdt4/1t6n6lR3\nqtnd2ev5jJGR1G7WmrOqduZ59m82ljHGICIiIiIiIiNlj7oBIiIiIiIionAmIiIiIiIyFhTORERE\nRERExoDCmYiIiIiIyBhQOBMRERERERkDCmciIiIiIiJjQOFMRERERIbmPe95D695zWtG3QyRsaRw\nJiIiIiIiMgYUzkRERERERMaAwpnIGh73uMfxN3/zN7zyla/k/PPP5wUveAG333471113HRdddBEX\nXnghb33rW4miCIDrr7+eJzzhCauucbTHREREJsX111/PpZdeynnnncfOnTv5kz/5k5XnPve5z/Gc\n5zyHJz3pSbz97W+n3W6PsKUi480ddQNE1oMPf/jDXHPNNTz84Q/nne98J294wxs4//zz+dSnPsV9\n993Hb/3Wb/GUpzyFV7ziFViWhWVZq95/tMdEREQmwQ9/+EOuvvpq/uiP/oif+qmf4o477uDtb387\n55xzDsYYPvCBD3D11Vfz5Cc/mS984Qv8+Z//OU996lNH3WyRsaRwJnICXv7yl/OsZz0LgBe/+MX8\n7u/+Lu973/vYsWMH5557Lo9//OO55557RtxKERGR4XvooYewbZsdO3awfft2tm/fzl/8xV+wfft2\n3vSmN/GSl7yEyy67DIC3vvWtfPOb3xxxi0XGl6Y1ipyAc845Z+W/y+XyyiDUUygUCMNwFE0TEREZ\nqWc84xmcf/75vPSlL+Vnf/Znef/7308URWzfvp177rnniGn9F1xwwYhaKjL+VDkTOQGuu/qjcrJT\nFOM47mdzRERExkahUODaa6/l+9//Prfeeiv//M//zHXXXcev//qvY1kWxphVr/c8b0QtFRl/qpyJ\n9JnneSRJQhAEK4/df//9o2uQiIjIAP3bv/0bH//4x3nCE57Ar/7qr3Ldddfxile8gi9/+cs8/vGP\n5zvf+c6q199xxx0jaqnI+FPlTKTPLrjgAizL4o//+I951atexe23387f/d3fjbpZIiIiA+F5Hh/7\n2MeoVqs8+9nPZnZ2lttuu40nPelJ7Ny5kze/+c2cd955PPOZz+RLX/oS//Vf/8WFF1446maLjCVV\nzkTWcCJTGA99zdlnn8373vc+brzxRi699FL+9m//lne84x2DbKKIiMjIPOUpT+Gaa67h85//PC98\n4Qt54xvfyNOe9jSuuuoqnvOc5/B7v/d7fO5zn+MlL3kJ3/3ud3n5y18+6iaLjC3LHD4RWERERERE\nRIbuhKY1fvKTn+SWW24hiiJe9apX8bKXvWzQ7RIRERl7Gh9FRKSf1gxn3/rWt/jOd77DX//1X9Nq\ntfj0pz89jHaJiIiMNY2PIiLSb2tOa/zQhz6EZVncc889NJtN3vGOdxxxXoWIiEjeaHwUEZF+W7Ny\ntri4yJ49e/jEJz7BQw89xK/+6q9y4403DqNtIiIiY0vjo4iI9Nua4WzDhg086lGPwnVdfuInfoJC\nocDCwgKbNm066uuNMSd9QK+IiMh6o/FRRET6bc1wduGFF3Lttdfy2te+lv3799PpdNi4ceMxX29Z\nFrOz9b42cr3YurWa275Dvvuvvuez75Dv/m/dWh11E0ZK4+PJyftnJa99h3z3X33PZ9/h1MfINcPZ\nRRddxLe//W0uu+wyjDFcffXV+uZPRERyT+OjiIj02wltpf+2t71t0O0QERFZdzQ+iohIP9mjboCI\niIiIiIgonImIiIiIiIwFhTMREREREZExoHAmIiIiIiIyBhTORERERERExoDCmYiIiIiIyBhQOBMR\nERERERkDCmciIiIiIiJjQOFMRERERERkDLijboCIyKkwxlCv10Zy72p1GsuyRnJvERGRtWiMXL8U\nzkRkXarXa3z1tl2UypWh3rfdanLx085lenrmpN532WUvwnVd/vIv/5pCobDquTe+8Zc566xzeOc7\nr+IZz3gK//N/vp/nPe/5R1zj0Oee8YynYFkWxpgjXmdZFq973et53eteTxiGXHvtX3DzzTexf/8+\nyuUy5513Aa997et57GMfd3KdFxGRdUFj5PodIxXORGTdKpUrlCvVUTfjhFiWxZ49u/nkJz/GG9/4\nW6d9vRtuuGnlv6+77lq+8Y1b+MQn/l8gG4hKpTIA11zzPn70o3v4rd96B2ef/QhqtSWuu+5afuM3\nXs///t/Xcs45jzjttoiIyPjRGJlZb2Ok1pyJiAzJjh0P42//9v/jjju+d9rX2rhx08o/pVIJ23bY\nuHHjymPFYpFWq8ktt3yVX/u1N/GUp/wM27dv5zGPeRzvfe/72bRpMzfc8Hd96JWIiMjp0xiZUTgT\nERmSSy55IeeddwEf+MD7iaJoCHe0sCyL2277d9I0XXnUtm0+8pE/44orXjuENoiIiKxNY2T3/iO5\nq4hIDlmWxbve9T/Zu3c3n/70Jwd+v3K5zM///OV84Quf5+d//hJ+93ev5h/+4e/Yv38f27dvZ8OG\nDQNvg4iIyInQGJnRmjMRkSE666yzufLKX+ETn/goz372c3jMYwa74Pg3f/NtPOEJP8U//MMN3HLL\nV/nKV74MwLOetZN3ves9VCpTA72/iIjIidIYqXAmIjJ0r3jFL/L1r3+Na675HT71qc+ses5xHIxJ\nj3hPb8cp1z35/21ffPHzufji5xMEHb773f/mlltu5ktf+nts2+Z977vm1DohIiIyAHkfIzWtUURk\nyGzb5rd/+708+OADfOYzn171XLU6TaPROOI9vfNqpqenT/g+3/nOf/Lxj39k5deFQpGnPOVneOc7\n38OrXvUa/v3f//UUeyAiIjIYeR8jFc5EREbgJ37ikbzmNa/j2mv/gj17dq88/tjHPo7vfe/2I15/\n++3fwbbtk5ri0Ww2+au/+iz33PPDI56bmppi06ZNp9Z4ERGRAcrzGKlpjSKybrVbzXV9zyuueC3f\n+MYt7Np1z8pjr3jFFbztbW/iU586m4svfj62bfODH3yfP/uzj/LSl778pL4VfPrTn8ETn/jTvOMd\nb+HKK3+ZJz3pwu60jdv57Gf/kt/8zbf1rS8iIjJeNEYe37iOkQpnIrIuVavTXPy0c0d275NnHfGI\n67q8+91X84Y3vBar+/STn/xU/tf/+jCf/exfcv31f0MYBpx55g5+4RdexS/8wi+e3B0tiw9+8CN8\n7nOf4fOfv46PfOSDgMWjH/0Y3v3u9/KMZ1x0Cv0QEZFxpzHyBO44pmOkZXor6Ppodrbe70uuC1u3\nVnPbd8h3/9X3fPYd8t3/rVuro27CupPXvyugz0pe+w757r/6ns++w6mPkVpzJiIiIiIiMgYUzkRE\nRERERMaAwpmIiIiIiMgYUDgTEREREREZAwpnIiIiIiIiY0DhTEREREREZAwonImIiIiIiIwBhTMR\nEREREZEx4I66ASIip8IYQ71eG8m9q9VpLMsayb1FRETWojFy/VI4E5F1qV6vccuuWymVy0O9b7vV\nYue5z2R6euaUr7GwMM/P//ylnHPOw7n22s+veu6yy17E/v37Vn7teR5btmzlWc/ayete93rK5TKd\nTocXv/hnednLXs4v//KvH3H9IOjwohf9LFde+QZ+4Rd+kcsvfzEvetHP8ZrX/N+n3GYREVk/NEau\n3zFS4UxE1q1SuUylWhl1M07aTTd9mR07HsYDD9zPd7/735x//hNXnrMsiyuueC0vf/krAWi329x1\n1w/46Ef/iDvu+C5/8iefoFgs8uxnP4ebb/7KUQeeb3zjn4jjiOc//wVD65OIiIwXjZHrc4zUmjMR\nsvJ/rba86h9jzKibJRPqxhv/gec+92d59KMfyw03/J8jni+VSmzcuImNGzexY8fD2LnzuXzgA3/I\n97//Pf7xH28A4NJLX8z+/Xu5447vHfH+m276Mv/jfzydmZkNA++LiIhIP+V9jFQ4EyEr/3/1tl38\ny/f28i/f28tXb9s1srnaMtnuuutO7rvvXp785KfxrGc9m69//Ws0Go013/eYxzyO889/Il/72lcA\nuOCCJ7Jjx1ncfPNNq163sDDPt799Gy984UsG0n4REZFB0RipcCayolSuUK5UKVeqlMrrbxqArA9f\n+tLfs3HjJi644Ins3HkxQRBw443/cELvfeQjH8W99+5a+fWll76Qf/qnm1dVeb/ylS+zceMmnva0\n/6vvbRcRERkkjZEKZyIiQxPHMV/72ld49rOfA8BZZ53NYx7zuKNO2ziaanWaZrO58uvnP/8FLC0t\n8p//+R8rj33lK1/mkkteiG3rf+8iIrJ+aIzMjG/LREQmzD//89ep1+tcdNFzVh7bufO53H//fXzv\ne7ev+f5ms8nUVHXl12ecsY0LL3zqyrSNe+/dxa5d93DJJS/sf+NFREQGSGNk5oR2a3zpS1/K1NQU\nAGeddRbXXHPNQBslIjKJvvzlfwTgN3/z147YcOaGG/4P5513wXHff/fdd/HoRz921WMveMGL+OAH\nP8Db3vbb3Hjjl3jCE87jnHMe3t+GyzFpfBQR6Q+NkZk1w1kYhgB85jOfGXhjREQm1cLCPN/61r/z\n0pdezkte8rJVz330ox/m61//Gm9601uP+f577rmbO+74Lldd9f+sevwZz7iIP/zD3+db3/om//RP\nN/Pa1/7SIJovR6HxUUSkPzRGHrRmOLvrrrtotVpceeWVJEnCW97yFi644PjJVUREVrvxxi9hjOGV\nr3w127ZtX/XcL/7ia/iP//gmN92UfWvYbrdZWJgHoNPpcOedd/Bnf/ZRnvSkC3ne8y5Z9V7f93nO\nc57Hpz/9SZaXl9m58+LhdEg0PoqI9InGyIPWDGfFYpErr7ySyy+/nPvvv5/Xv/713HTTTWO9kE5E\n8qHdaq2be9500z/y9Kc/84hBB+Cnf/rJnHvuo/n7v/87wOJzn/tLPve5vwSgXC6zbduZvOQlL+Pl\nL38llmUd8f4XvOBFfPGLX+CSS15IqVQ6yt2PfI+cPo2PIsdmjFk5ksb3U2q1OtXq9FH/HyaDoTEy\ns97GSMuscdJuGIYYYygUCgBcfvnlfPSjH2Xbtm1DaaDIMCwvL3Pztx6gMjUNQLNR47lPfTgzMzMj\nbpkcS3Zw+GjOopue1g8YovFR5HiWl5e54et3Uu4eTdNqNXnxRT+pcXVINEauX2tWzr7whS9w9913\nc/XVV7N//36azSZbt2497ntmZ+t9a+B6snVrNbd9h/Xd/1qtTqMZkNIBoNUMmJurE4Yn9g34eu77\n6Rpt30dToZibO3ggZt7/7PNM4+PJyftnJW99r9XqpMYlxac6VaRxkuPqpMjzGJnHv/eHOtUxcs1w\ndtlll/Hbv/3bvOpVr8K2ba655hpN2RARkdzT+CgiIv22ZjjzPI8PfvCDw2iLiIjIuqHxUURE+k1f\n8YmIiIiIiIwBhTMREREREZExoHAmIiIiIiIyBtZccyYiIiIiMgiHnocG6Cw0yT2FMxEREREZiXq9\nxi27bqVULtNutdh57jOZntZZaJJfmtYoIiIiIn213Iz49t1LLDWCNV9bKpepVCuUyuUhtExkvCmc\niYiIiEhf3XHfEvfvb3HTbQ8xu9QedXNE1g2FMxERERHpq+VWBEAQJXzxGz9ioR6OuEUi64PCmYiI\niIj0Va2ZhbPzHrWZMEq49Xvz1LuBTUSOTeFMRERERPqq1g1iP/mIDZz/qA3EieH7981Sqy1jjBlx\n60TGl8KZiIiIiPRVrRXh2BZx0KJpPQDAf+9+gFt23bpq63wRWU3hTERERET6qtYMKfo2lmUxvTE7\nuanecrAdi3q9puqZyDEonImIiIhI36Spod6OKfoOAIWCwbZT5hZS7t1T45++fa+qZyLHoHAmIiIi\nIn1Ta4UYAyU/+zHTsqBYSggCF9ctUixVRtxCkfGlcCYiIiIifbPcyLbN71XOAIrFGIxFEOhHT5Hj\n0SdERERERE6bMYZabZk9BxYAKPoHf8wslhIA2m3nqO8VkYw76gaIiIiIyPpXr9f46m272Luc/dol\nWXmuWIqBbjibyR5LTUo77gy7mSJjTeFMZA3GmCMWLler01iWNaIWiYiIjKdSuUKyHABQOLRyVsyC\nWqtj+IH9LX7w/a8zHywQpTFP2nAej6s+eiTtFRk3Cmcia+h9E1gqZwuY260mFz/tXKanZ0bcMhER\nkfHTDrIqWdE7GM48L8VxUgJ3iX32/fiBx9bSFvY093EgmONxKJyJgNaciZyQUrlCuVKlXKmuhDQR\nERE5UivIqmQF7+AMk96OjbFXB+AXH/Uy3v3Ut1B0CixH2lZfpEfhTERERET6ph3E2BZ4zurp/8Vi\ngl3JwtmO8nYsy+Jh5TNpxE3CJBpFU0XGjsKZiIiIiPRNO4gpeNYRa7OLpQS7XMNNSkx52SyUsypn\nArAYLA29nSLjSOFMRERERPrCGNMNZ0f+iOmXW1h+gBscXLP9sHIWzhY6i0Nro8g4UzgTERERkb4I\nohRjVq83W1HO1paZ1vTKQ73KmcKZSEbhTERERET6ohOmABSPEs4ipwlAUDsYzjb4M/i2x0JH0xpF\nQFvpi4iIiEifdMLuTo3ukd//B1YLgPZylfnFLIw1GnVmvGlmg3nCNBxeQ0XGlMKZyFGEccrsckdn\nmYmIiJyEdrdyVvCPrJx1aELsQ1Tgnx+4g23tIvOzc1TdMrPMsxRqS30RTWsUOYr/vHuJ3/ur7/Mf\ndx0YdVNERETWjV7lrHhY5SwmIrYi3KgMQBCXqFQrlEolpt0qgM47E0HhTOSollsRaQp/9sU7+M+7\nF0bdHBERkXVhZVrjYWvOOmRTGgtpCYDlQ3LYjJOFs6VoeQgtFBlvCmcihzHG0OokbKh4FH2Xz958\nHw/sb426WSIiImOvtyHI4eGst96sTBFYHc4Ktk/RKbCkypmIwpnI4VpBTGrgEdunePsrn0ix4PBf\nu5ZJUjPqpomIiIy1dphgWeC7h1XOeuHMKuI4yapwZlkWm4obaScdGlFzmM0VGTsKZyKHabQiADZP\n+zxi+zQXPHIDSWqoN7WLlIiIyPF0wpRSwcWyjpzW6BgXx7j4hYhGA+L44Jeem4obANjd2jvU9oqM\nG4UzkcM02lk42zRdAGD7pmx+/GIjGFmbRERExp0xhk6YUCqs3gw8NhGxFVKkgoVFwc/G2eXawXC2\nsZCFs32t2eE1WGQMKZyJHKYXzspuTK22zIZSNnjMLtQxRlMbRUREjqYVJKSGI8JZ23Q3AyH7srNQ\nyGaiLC2lK68pe9lztag+jKaKjC2dcyZymN60xh/Vv09nz25qaRsos3tpjp9qTY+2cSIiImOq1h0/\nywUHOPhlZmjaAPgm2wzEL2SvW1xOmc4yGSW3u4ujzjqTnFPlTOQwK9MaNxaoVCts3FzBdVOC0Btx\ny0RERMZXO8i20fc9Z9XjgekA4JEtFyh0w9nS0sEAV3QKWFiqnEnuKZyJHKbejih6FvYhn45yJSUM\nHOIkPfYbRUREciyIsnDmOat/vOyFM78bzhwnpeDD0vLBMdWyLIpOgeVQ4UzyTeFM5BBJmtLqxJQL\nqz8apXI24Cw341E0S0REZOwF3TPOXHf1GBqaDpaxcPEBsCyYmYZ6w5AkB19XtIs0ogap0Rehkl8K\nZyKHaLaz8FXyDw9n2UDRm08vIiIiqx2vcubhY3Fwe/3paTAGGq2Dry05RVIMNVXPJMdOKJzNz89z\n0UUXcd999w26PSIj1VtvVvItOp2ATqdNp9OhXMrCmSpnInI4jZEimU6UjZXeIZWzyAQkxCvrzXpm\nZrJ/15sHX1t0stcsBcsDbqnI+Fpzt8Y4jrn66qspFovDaI/ISPV2avSsiF0/XmRT4NBs1HG9IlBh\nuanKmYgcpDFS5KAgzCpnrmOvbNbYMlkVzDOrw9mG7ubH9YYNZKGu5GSfo6VAOzZKfq1ZOfv93/99\nXvnKV3LGGWcMoz0iI1XvVs6KHnieh+8X8TwfxwG/kFBrqXImIgdpjBQ5KOhWzlzHot1q0G41WOxk\nh0r7h1XOpnvh7GiVs44qZ5Jfxw1n119/PZs3b+bpT3+6Dt+VXFiZ1niUXfNLpZQgSumEyZFPikju\naIwUWa235iyNOzzQ+QF7uJeHOvcAR1bOioXsn1rdpvfxKdq9ypnCmeTXcac1Xn/99ViWxb/+679y\n11138c53vpM//dM/ZfPmzce96Nat1b42cj3Jc99h/fbf91OmKgu0g0Uc22LDtE/TeBSKHmHoYlkO\nU1XD8hKEqcWWLVVmZlb3db32vR/y3HdQ//PqVMbIvP9dyXP/c9H37hk001MFppIpSlNVFoN9AJT9\nCr7tYaUJdmKoVkucc3bC3btC4qRItVpkSzQN89CxWhP1+zVJfTlZee77qTpuOPvsZz+78t+vfvWr\n+Z3f+Z01gxnA7Gw+d9nZurWa277D+u5/rVan0QxYboRUii7tdkRoRwSdiDCIsay0e6imz/75FnNz\ndcLwYOF5Pff9dOW575Dv/ud90D2VMTKvf1dAn5U89L3WCAAIw4gwjHGCKDvjzAChQ2hFREFMnMY0\nGh3O3O5x9y649/6Ybds7JO1sN8d9y3MT8/uVlz/7o8lz3+HUx8gT3krfsqy1XySyjkVxShAlTJWP\nMqcRKJd01pmIHJ3GSJHDNgTpiq0QFw/7KD9ynrXDwbIM+w44ANiWzZRbYVkbgkiOrblbY89nPvOZ\nQbZDZOSaQTaoTJU8oHPE84ViimXBss46E5HDaIwUyTYEsSywu99VpKTEVkSZo1cQCgWLzRsS5hZd\nWq1sM5Fpf4rZzgLGGH3pIbmkQ6hFulqdrCI2dbTdQMim0ldLLrVmTKrF/yIiIqsEUYLrWCuhKiIA\n68idGg+1bWv2xeiDP+6GM69KlEa04/bgGywyhhTORLqanW7lrOwf8zUzFY8kNSzWw2E1S0REZF3o\nROmqKY0R2Rq044Wz7VuzL0bvfyAkCNqU0Flnkm8KZyJdK+HsGJUzgOlyNhN474K+0RMRETlUGCV4\nzsGpiOEJhLNyyVCtxOzdb7jnoTpzB7IvPxe1nb7klMKZSFdrZc3ZsZdizlS64Wxe4UxERORQnTDF\nPSScRdba4QzgjC0Rxlh02lNU3BkAlhXOJKcUzkS6gijFgu6W+Uc33d3Jcd/CkRuGiIiI5FWcpCSp\nWRXOwu7mWl53quKxnLEl22hracnDNyUA9tX2U6stU6st65B3yZUT3q1RZNKFUXaWmX2c3aEqRQfH\ntjStUURE5BCdo2yjHxFgGwcXl5j0mO+dria4Xkqt5kI7gircvW8/xfm9tFtNLn7auUxPzwy8DyLj\nQJUzka4wTil4x/9IWJbFdNll/2KHJD32QCMiIpInnTDb2KNXOTMYIgK89NibbPVYFszMJESRTRRk\nUyAjJ6RcqVIqVwbXaJExpHAmAhhjCKOUgn/sKY090xWXJDUcWFT1TEREBA6tnGXhLLFijGVwzdrh\nDGB6Qxbuag0Px7i008ZgGioy5hTORIB2mGCAwnHWm/XMdNed7Z5tDrhVIiIi60Nw2LTGxMrWkZ1I\n5QxgZiYLZwtNQ4GywpnklsKZCNDqbqN/opUzgB/PauAQERGBIytnsd2d5miOfTzNofyCoVhMWGqB\nb0qEpkNi4sE0VmSMKZyJAI1ONgCcVOVsTpUzEREROLjmzFuZ1phVzk40nAFMz0SkBqw4291R1TPJ\nI4UzEaDVC2cnUDkr+jalgqNpjSIiIl3Hqpw5JxPOprP3RO1sU5B2qnFW8kfhTARonkTlzLIsztxU\nYv9iiyhOBt00ERGRsXf4VvqnUjmrVmMsoN3M1qm103p/GymyDiiciXBy4Qxg+6YixsDe+dYgmyUi\nIrIuBNFRKmfGwjEnfqSu48BMCVqNLNB1Uo2xkj8KZyIcEs5OYFojwJmbSoB2bBQREYEjzzlLrAgP\nDwvrpK6zsWJBnFXOOkbhTPJH4UyEg+GseIKVs144+/GcFiuLiIh0goPTGlOTklgxLie2jf6hNpTB\nRNmaM1XOJI8UzkQ4ua30AbarciYiIrKi053W6DkWAW2wwDuFcFb2wUTZ+wJVziSHFM5EOFg580+w\ncjZVcpmu+ApnIiIirN6tsWOysfFUKmeeAw4upI4qZ5JLCmciZOHMdSwc+8Tnxj9sS4X5Wmdlnr2I\niEheHbrmbCWcmZMPZ5ZlUS44mMhXOJNcUjgTIQtnBe/kPg47NlcA7dgoIiIShAkW4NgW7W44O5Vp\njQDlgo2JfALTxhjTx1aKjD+FM8k9YwytTozvHv3jsGzN0+LIs1Z2bM3C2Z45TW0UEZF864QJvmdj\nWRYdTn1aI/TCWQFDSkTQz2aKjD2FM8m9MEqJEoN/lMrZfLyXffb97LPuP+K5HZvLAOyZVzgTEZF8\n64Qxxe6mWp01KmfGQKfTodNpE4QdgqADhxTIyr59yKYgncE2XGTMnPjJgCITqtGOACgcVjkL7Q57\nw/uy/7Y6xGb12rIdW7rTGuc0rVFERPItCBOKfjaOtk0TTLaxR0x4xGvjKOTuBxcoVyoszdbpBC02\nn+HjWFm4KxccTDvbTj+kPbxOiIwBVc4k93rh7NDKWUzEgeJDGFKKJgth7cOmNlbLPtWyp2mNIiKS\ne50woeD1KmctHHP8A6g9z8f3i/h+Ac/LqmTGGNqtJo7prFTOtCmI5I3CmeReo9MNZ4dUzu4t3UFs\nh2xzz2FLeibA0dedba4wu9Qm6J7vIiIikjdJmhLGKQUvO4C6Qws3PfnJWUG7ze5kF4ve/SsHUTfD\n5X43V2SsKZxJ7jVaqytnraTOgrcXPynxMO+RFJkCAy0aR7x3x5YKBvjx/iODm4iISB4EYQpAwXey\nA6gxOMY7pWv5xSLlqRIu2ftDSxuCSL4onEnuraw564azpWQWgHJcxbJsHBwKlGjTIDXpqvf21p09\npHAmIiK7nRNgAAAgAElEQVQ51TvjrODZh5xxdnrbGhSc7oYgaEMQyReFM8m9Znv1tMZeOPPT4spr\nSkxhrJSaWVj13l44e1DhTEREcqoTZlP7C56zEs6c9NQqZz1FN3t/J1U4k3xROJPcO3xDkKX4KOHM\nTAGwmO5f9V5VzkREJO96666LnrNyAPVpV84KLia1VDmT3FE4k9w7uJV+tqvUYjKLm/o4hwwsJbrh\nzBxY9d7pskel6CqciYhIbnWCbFqj79krB1CfbuWsVAITFYgshTPJF4Uzyb1DK2dhGtBKa1TS6VVb\nAHsUcIzHYro6nFmWxY4tFfbONYni1evRRERE8qDTq5z5ziFrzk5zWmPRQOwT2wHGmLXfIDIhFM4k\n9xrtCM+xcB2b5d5mIMn0qtdYWJSZokOTVlqnXq9Rqy1Tqy2zY3OF1MD+BZ3FIiIi+dNbc+Z7Nm3T\nwsLCNs5pXbNYMtl2+lZKkB55kLXIpDq9CcEiE6DRjigXs49CbzOQSjJNTLTqdWWq1Flkf+chvvFf\nIRs2babdarKpmgW53XNNzjpjariNFxERGbFeOCt6Dp12kyKV4x5AfSJ8n5WDqBtRkzNOu5Ui64Mq\nZ5J7jXZEpRvOFrubgVTS6SNeV6YKQM2ap1gqU65UKZUrbN+YbRyyZ645pBaLiIiMj6AbzjwPAloU\nrcppX9OyWDkrrRFpfJX8UDiTXIuTlE6YrISzpWQWB5dieuTAkn0TaLNsza16fPumbjib1+AhIiL5\n0zvnLHUCDIaiVe7LdXsHUc+3an25nsh6oHAmudY746xSdElNQi1ZYMbZfNTpGDY209YmGiyTkn1L\naIzBStqUCg67Z+vUastauCwiIrnSm9YYW22AvoUzz8rC2YGmwpnkh8KZ5Fq9G87KRYe6WcKQssHd\neszXz1ibMVZKk2UA2q0mt37nIQqezfxywFdv20W9rkFERETyoxfOom44K1Dqy3V9O5vVosqZ5InC\nmeRar3I2VXSpmQUANjjHCWf2ZgDqLK48ViyVma4UiRKDV+jPgCQiIrJe9KY1hnTDWb8qZ91wVtea\nM8mRNXdrTNOU97znPdx3333Yts373vc+zj333GG0TWTgemecOVbMQrgPbChGFZqtOmbqyOmJ09YW\nAGrWwqrHK6Vs6kU70FlnInmh8VHyzhhDvV6j0QoAqHWyLy4LVr8qZ9nY2owVziQ/1qyc3XLLLViW\nxV/91V/x5je/mQ996EPDaJfIUPSmNf64cR/z7AEDi+F+Hgp+SBxGR7y+am3AMtaqyhnAVDectYJk\n8I0WkbGg8VHyrl6v8dXbdnFgKauY/ejAPqB/0xoLroMxFp203ZfriawHa1bOnvvc57Jz504Adu/e\nzczMzMAbJTIsvWmN5YpLYHcoWmXKpSp+sXjU1zuWS4UZGiyRmoNVsl7lrDfvXkQmn8ZHESiVK6Sm\ng+tYpG73vLN+TWt0DUQ+oaNwJvlxQodQ27bNu971Lm6++Wb++I//eNBtEhma3rRGvA5pmlCy1z5E\numo20rCXqCcHpzZWVDkTySWNjyIQxSmuYxNYHSxj4VHoy3Vd12A6PrGrcCb5cULhDOADH/gA8/Pz\nXH755XzpS1+ieIzKAsDWrdW+NG49ynPfYf31P+4tKyu3oQHT/gYKRQ/Pd3E8l0LRIwxdLMvBsW2m\nKgU2J1vYG99H211iprIR2/aykzeBxMCWLVVmZtbX78PpWm9/7v2W9/7nncbHE5fn/k9i330/Zaqy\nQJIafM8htgIKlKhOFfHrLn7Bw0qzLy39QjaWWmlCnGbPFYoefsFd+fXhz6WJC40Cxq4zvbFAwfVH\n3ONTM4l/9icqz30/VWuGsy9+8Yvs37+fN7zhDRQKBWzbxraPv1RtdrbetwauJ1u3VnPbd1if/Z9d\naAGwFGZVsEJaIehERGFMaiyCTkQYxFhWimM5NAjwO1VwYV9zN06jTBgss3V7Fs4Wl9vcd99uduyw\nsKwjz0qbROvxz72f8tz/vA+6Gh9PTt4/K5PY91qtTqMZEEYJvmfTMe1s6n8zIAxjnCAiCmKKZWdl\nLI2CmCiMCYOIwMvG2N6vj3iuE2MnWSC7d89etpQ2jbjHJ29S/+xPRJ77Dqc+Rq65Icjznvc87rzz\nTq644gp+6Zd+iauuugrfX5/fXIgcrt4KcWyLRtoAoGyv/UGaYgMYi8XkAJ12k/vbd/JAeCeWZaiF\nbf7lwW/qrDORHND4KJLt2BgnBs9PSa2EAseuHJ8KJ82uVws0rko+rFk5K5VKfPjDHx5GW0SGrt6K\nqJQcakkdnyKOtfZMXweXCtMsxbMYDH6pSHmqQqEAcexQKuusM5E80PgoAkmarQ9w/BAAv087Nfa4\naYEEmGst88gNfb20yFjSIdSSa/VWRLkaEJuEkrX2ZiA9VTYSE9GxD5694hcMUWSR6qgzERHJiSjJ\nwpm1Es76VzkzxuAk2Zeme+YPUKstY8yRZ5CKTBKFM8mtMEoIogS3ms2HLp9MODMbAWjaB6dZFAoG\nsOgE+VhrJiIiEnd31rK87CDqgulf5SxotwnjDgA/XLifW3bdqmUDMvEUziS36q1sG31TWgI4ycpZ\ntii56SyvPOZ3dw5udxTOREQkH+KkO13Ey0JUPytnAJ6TbbgVWQmlcn/OTxMZZwpnklu1VjYFI/IX\nAShSOeH3VulWzpzDK2fQUTgTEZGc6E1rxM0qZ/0OZy4OAEEc9vW6IuNK4UxyK6ucGVr2AhWnjGM5\nJ/xeF48pewMNZ4mU7FtDvxvOVDkTEZG8iLvhLHWyylmhzxuCeHY2pgZp0NfriowrhTPJrXorxCo2\nSYiYcU58SmPPOf5jSK2EpptNiyx0pzVqzZmIiORFb1pjYg9mWqNjg0ltIqPKmeSDwpnkijGGWm2Z\nWm2ZucU6diWbllh1T/6gwEcWz8MyFjVvAWMMvq/KmYiI5EuvcpbYbRzj4ax9StNJcZ0UE/kkKJxJ\nPvT3EyQy5ur1GrfsupVSucwPZ8GuZBt6lLrnqJyMkj3Fpng7895eavEiBb8KGK05ExGR3OitOYus\ndt8PoAZw3RRin8Rrrv1ikQmgypnkTqlcplKtkKQuVrdyNmWf2g5Q28OHA7C38yCWDZ5vVDkTEZHc\niJMUrJSYTt+nNAJYloHYAzuh1WnonDOZeApnklvtTopdrlGxTm4zkENNJRvxkyIL0QHCtIPvZ5Uz\nDR4iIpIHcWLADcEC3/Q/nMVRCEm2nf4PHjpAo1Hv+z1ExonCmeRWO21gOQnTp7AZSI+FRTXaBBhm\n4934fkpqLJqduH8NFRERGVNxYlYOoB5E5QzATrNwZnuamSKTT+FMcqtjZ7ssTp/CZiCHqsQzuJbH\nXLwHt5CFsj0HllY2HlEVTUREJlUUpwfD2QAqZwC2ycJZnJ7s6nCR9UfhTHIr9rNwNuOcXjizsdlW\nOIuYiGTqxwDcdtcs//K9vXz1tl3U67U1riAiIrI+ZZWzbCdFv89nnPX0wlloFM5k8imcSS7FiYHy\nMhj7lDcDOdSZxYdjYRFWdwMGY/mUK1VK5crpN1ZERGRMDWNao9vdXDw2WjIgk0/hTHKp2Y6xSnW8\neBrbOv2PgW8X2ORsJ3E72BsP0Ar07Z6IiEy+OEmxC1k4KwxoWqPb3bQrOulDb0TWH4UzyaXZxhKW\nbSilG/t2zW3eOQB4Z95LO9S3eyIiMvmixOD4vWmNgwlnnpVVzhKFM8kBhTPJpfnOIgBTzoa+XbNk\nV6ikM9hTy9Stub5dV0REZFz1pjVaWHgUBnIP38kqZyn64lMmn8KZ5NJSmIWzDe6mvl53k9kGQGfD\nPX29roiIyDiKEwNeQMEqYzGYre49x8EYSGyFM5l8CmeSS410ERN7TBdPfzOQQ5WYwjSnMdP7aaeN\nvl5bRERknCSpIUlTjBtQtAe3AZbnGoh9jB0N7B4i40LhTHInSEJCq0XamKFU6u9HwMLCqm8GYDbc\n09dri4iIjJMgSsDO/vGMT7vVoN1q0Go26OcRn65rMJGPcRTOZPK5o26AyLAtRcsApM0ZSkWLsNnf\n6zvBFAlwINzDFufM/l5cRERkTARhunLGWSdsste6DwubeCkitfq3eYftGOh4UG6QmrRv1xUZR6qc\nSe4sht1w1thAsdD/+fFeXMEYi/l4X9+vLSIiMi6CKAG3u1OjU6RYLlMslymU+rtro2WBlfoAtJJO\nX68tMm4UziR3epUz05rB9/t/fdcB06pSN7MkRtv+iojIZOpEBytnLt5A72Wl2fXbSXug9xEZNYUz\nyRVjTFY5C0sUXR/L6n/lzHUT0sYGjJVSM/N9v76IiMg4CMIEq1s5cwa8UsY22bep9bA10PuIjJrC\nmeTKQrBEZCLS5sxApjQCuG5K2sjOT1tKZwdyDxERkVELomRolbPe9ZcChTOZbApnkit7Wtk6sKQ+\nQ7E4qHCWrISzRXNgIPcQEREZtU6Yrqw5G3g4s7LrLyucyYRTOJNcme1k0wzT9tRAw5kJSlhJgcVU\n4UxERCbToZWzQU9r9LrhrB4pnMlkUziTXOmFM9OpUDrJaY3GGFrNxiHnuNRpt5oYVh/mYtsJYGG3\nNtKhyZ6FPdRqy5h+HvoiIiIyYkGUrqw5G3TlzHOyNWcNhTOZcDrnTHJlLljAwsIEJRw3ptNJCMIO\nURyBqRz3vUG7zY/i72LZFhY2y51Z9gb3U/Krq14XdTo4TkzUmMKuwq0HvskZC1vYee4zmZ6eGWT3\nREREhiYIs8qZbRxsHFIGdwaZb3d3a0y1lb5MNoUzyQ1jDHOdeXxTpoXFYr3Brt0pS7N1OkGLzWf4\n+IXjn81SKBaxXQfLsilWyvjFo7/e8wxxbSP2mdC02pTK5UF0SUREZOiMMdTrNerNNhRDXFPADGal\nwIqim4WzTqqt9GWyaVqj5EYjatJJArwkq5AVSza+X8T3C3hefw8887yUuJ5VyRbCpb5eW0REZJTq\n9Rq37LqVPY15LDckjVOicLAVLc+zMbFLYIKB3kdk1BTOJDf2NbPNOayoBIDnDm4NmOsmkLoUTIWl\nYJnEDG6qh4iIyLCVymVSCywnXdlJcZA8z2BijwhNa5TJpnAmufHgwkMAtGuF7AErHNi9XC8LY348\nTYqhHtUHdi8REZFRCE1vM5DBr5JxPSD2SexAG2zJRFM4k9yY6+3UGGQbeHje4KpZnptd24mzKZS1\nuDGwe4mIiIxC1A1njnEGfi/HARP7YBk6iapnMrkUziQ3ZoMsnCWNKWw7xRngWOJ2wxlBFs4aCmci\nIjJhYjOcM856rCSbPtkItZ2+TC6FM8mNuc4CrvGIOj6eP9g1YK6XAGDaWZWuHjUHej8REZFhS6xe\n5WxY4Sy7TyPSF54yuRTOJBeSNGEhWKRkqsSxM9ApjXCwcpYEPgXbp67KmYiITJjEioDhhTPbZJWz\n5Y7WccvkUjiTXJjvLJKYlEIy+PVmh14/iiyq3hStpE2YDG4DEhERkWFL7eGtOQOw0ywEzjaXh3I/\nkVFQOJNcONCaBcDtbtAx6GmNtm2wbUMU2kz7UwDMdjckERERWe+SBHCHO63R7VbOFlq1odxPZBSO\n+2mK45h3v/vd7N69myiK+JVf+RV27tw5rLaJ9E0vnNlhFpS87pqwQbEs8PyDlTOA/Z05Hj/Qu4rI\nMGmMlDyLY7C84Yaz3sYjSx2FM5lcx/003XDDDWzcuJE/+IM/YHl5mZ/7uZ/TwCPr0v5uOKNTBgY/\nrRHA9w2NukXVzcLZgfbcwO8pIsOjMVLyLIrBckMwFtaQJmJ53R9ba6HWnMnkOm44u+SSS3j+858P\nQJqmuO5wvhkR6bcDrSwYJZ0KkAwlnHk+GGNRoBvOOrMDv6eIDI/GSMmzOAa8ECvxsbCGck/XdjBG\nx9PIZDvuSFIqlQBoNBq8+c1v5i1vectQGiXSb/tbs2zwpwnCbAAZ9JozyCpnAGno49se+1U5E5ko\nGiMlz+Ju5cxOSkO7p+cYiH1ats45k8m15td8e/fu5Td+4ze44ooruPTSS0/oolu3Vk+7YetVnvsO\n49V/Ywy1Wo1OHLAc1njMzCN4KEoAQ7li4zgOhaKHX3CJUxe/4FEoZouNPd/F8VwKRY8wdLEsBytN\nsB0H27OxrOy9vdcBK6/1fA/f9yiVs2keluWzoTDNbHuBmY0FfNcf1W/JwIzTn/so5L3/eXayY2Te\n/67kuf+T1HffT7Ox0ElwEh+/4K4aH600WRlXrTRb4+0X3COeO3wMPt5zhaJHsWRjIp/Qa62r38/1\n1NZ+y3PfT9Vxw9nc3BxXXnkl733ve/mZn/mZE77o7Gw+5wJv3VrNbd9h/Ppfqy1zy65bCdwYgHYn\npNYJcF2HOIpJ4pTAiwiDmCiMCYOIwMvObInCmNRYBJ3sectKiYIY202xUwfLyt7bex2w8toojAjD\nCNtOAJeFxZByuYxhnjsevI+zqztG9VsyEOP25z5see5/3gfdUxkj8/p3BfRZmaS+12p1FutNsMBK\nXcJo9fgYHTKuRkFMseysGksPHXPDw157rOcCL8KYCBMViGmwZ98CnuON+rdiTZP2Z38y8tx3OPUx\n8rgrOD/xiU9Qq9X4+Mc/zqtf/Wpe85rXEIY6q0nWj1K5TOxl4WxDeYYkcYey3gwOTmvsdKDqZVv4\n723uG8q9RWTwNEZKnnXi7k6NDC8gOU6KiQoA1EKtO5PJdNzK2VVXXcVVV101rLaIDEQjagJQoIQx\nNq43nB+e/O7sxU4HNnV3bNzXPDCUe4vI4GmMlDwL0mwsdYcYzmzbYMW9cFZnc2nj0O4tMiw6hFom\nXjPKFg7bcbZoeViVM69bOWt3YGolnO0fyr1FREQGKUyzZQDu2tsX9I1lgWeycFbXdvoyoRTOZOL1\nwpkJu+HMHewB1D29ylm7AwXbp+SU2KtwJiIiEyAy3WmN1nDXfflkY7nOOpNJpXAmE68Zt/Adn7CT\nfbs3rMqZbWf3ajbBsiy2lbYw254nSqKh3F9ERGRQIrJw5g85nJXsMgDLgcKZTCaFM5loxhhaUYsp\nt0y7nT3mesOpnAH4hZRWC4yBM0pbMBj2t3QYtYiIrG9xN5x59nAPXy87WeVsvlUb6n1FhkXhTCZa\nkIYkJqXslWl1w9mwKmcAhYIhNdnUxjOKWwGtOxMRkfUvsbqVM2d44cwYQ8FkawbmmwsYY4Z2b5Fh\nUTiTidZOskRW8Q5WzrwhVs4KhSwINluwrbQFQOvORERk3UutEJNauLYztHsG7Q61cAljYG9jP/W6\nqmcyeRTOZKK1kg6QhbNWGywrxXGG902bXzw0nGWVsz0KZyIiss6ldgixj+NYQ71vpeRCVCC24qHe\nV2RYFM5korXi1ZUz10mwhjiOFApZEGy1YMqtUPHKOohaRETWPeNEkPhDv6/nGUzkE5lg6PcWGQaF\nM5lovWmNJadMpwPOkLbR7zl0WqNlWZxZ2cZce4FQOzaKiMg6FacJODH2CMKZ76WYqEBqJYRJOPT7\niwyawplMtFY3nLlJEQO47nCnQfiHhDOAMyvbuzs2HhhqO0RERPqlETUBsM2IKmdxdt9G3Bz6/UUG\nTeFMJlo76eDZHlGQncPiDrly5jhQKBwazrYB2hRERETWr/l2dsaYyygqZwYTFQCoRwpnMnkUzmRi\nGWNoJW0qXplmK1v7NezKGUClDK02pMYonImIyLq30MrCmWeNonKWQtStnCmcyQRSOJOJ1UraJCbJ\ndmrshTNnuJUzgKkKpKlFvRUfEs60KYiIiKxPS50GAAXHG/q9XYeVylkjagz9/iKDpnAmE2spWAag\n4h5aORt+OKuUs38v1AOq/hRTXoW9DVXORERkfVoOsopVwRl+5cyywLe60xq15kwmkMKZTKzFsBvO\nvDL1RhbOPG/4uyRWKtm/F+rZrlI7KtuZ6ywQaJcpERFZh+phFopKbmEk9/esXuVM4Uwmj8KZTKyV\nyplXpl432PaIKme9cFbLzmQ5cyqb2rhP685ERGQdanYrViV/+NMaAYpOEYB6qGmNMnkUzmRiraqc\n1VOmKgz1AOqe3rTGxW7lTJuCiIjIetY7Q7TiD39aI0C54GOMxbLCmUwghTOZWEvdcOamJYIQpqZG\n047DpzWeWdkOKJyJiMj61DFtjLGoFEZTOSuXbIh8baUvE0nhTCbWYrCEYzkELRfIdk0cBc+18H3D\nQr07rVGVMxERWcci04HIp1gcwXQUoFS0MFFhZXqlyCRROJOJtRjWKDtF6t1ZD9URVc4gm9q4WA8x\nxlDxykz7VYUzERFZl2I7wMQeI1pyRqlkYSKfhJhOHIymESIDonAmE6kdt+kkHUpOiXo926lxVNMa\nAcoliBLD7v3z1GrLbC1sZqGzSCfqjK5RIiIiJylOY4wdYSU+1igWctMLZ9mOjbWwPpI2iAyKwplM\npIXOEgBlp0StngKjm9YIUOhu4X/rvbfz73u+TZpmbfrR7L2ja5SIiMhJ6m1fb6fD3wzEGAiCDo4d\nYKLs/rWgNvR2iAySwplMpPn2AsDqytkIw1m5mLUhNgUq1QpbpzYDsLetqY0iIrJ+1IJsrYDL8MNZ\nHAXct3uJfUt1TJhtp793ed/Q2yEySApnMpFWKmdukVrDUKlYOM7wp18YA51OB9fL5sT3DsPeVNwI\nwI+be4feJhERkVM138wqVaMIZwCu51MuF1bC2WL3TFORSaFwJhNprjMPQMEq0WoZpqdGMy8+CgPu\nfnCBxXo2mC0vZ9MbZwrTOJbDQ809I2mXiIjIqZhrZmHIs0YTzgBsG+w4C2e9Y3NEJoXCmUykuXYW\nzgiyE6Cr1dGEMwDP8ymXHcBQ606Nty2bGW+a/e1Z7TQlIiLrxmI7G8gKzoi2auzyTK9ytjTSdoj0\nm8KZTKTZ1jxFp0jQyr7Zm66O9q+64xgKhYj5RUjTbGrjRm8Gg+Gh+o9H2jYREZETtdzJ1pwVndFV\nzgBcx8bEHgua1igTRuFMJk5qUuY6C2wubKDRPZ9ylJWznlIpIElgfiHbqXGjPwPA/bWHRtksERGR\nE1YLu+HMHW048zyDCUoshcsYY0baFpF+UjiTibMc1IjTmE2FjSvhbNSVM4BSOZu+eGA2C2cbPIUz\nERFZX5pxNrBWvBGHMzfFBCViE69s7y8yCUb/E6tIn/XWm20+JJxVR7QhyKFKpdXhrOQUmXIrPKBw\nJiIi60Q7aWGMRbkw4jVnnlnZsXG+szDStoj0k8KZTJzZ7hlnm4tZOCsWwPdHH848L6ZYgP2zKcYY\nLMvi7MoOFoMllnWIpoiIrANB2obYo1gY7bjqeilpUALgxwu7qdU0vVEmg8KZTJzZ9hwAG7wNNFtQ\nHYMpjQCWBVu2QKtlaLayAeTsyg5AUxtFRGT8GWOI6GAiH0ybTqcDo8pDaRvTDWe3z32fW3bdSr2u\nLzpl/RuPn1pF+qg3rdGNpzDGYnoMNgPp2bo5+/f+A9nUxrNWwtmDo2qSiIjICVlaXiS1Q0zs8+CB\nGj98YI4oikbSFre75gwgsmNK5fJI2iHSbwpnMnHm2vO4tkunmc2HH4edGnu2bMn+3Vt31gtnWncm\nIiLjrhm3ALCTAoVCEW+Em4K4booJs3DW1IYgMkEUzmSiGGOYbc+zpbSZ+XoIQHVqfP6ab9oAjg0H\nZhMASm6RbeWtPFD7MalJR9w6ERGRY+uFMycd7U6NAI6TQuJC4tKM2qNujkjfjM9PrSJ90IxbtOMO\nW0ub2LeQ/c96w8z4VM4cx2LzZpuFRUMUZ489YvocOkmHA63Z0TZORETkOOphVqFyTXHELcnWcbtu\nAmGJZtTUZiAyMRTOZKL01pttKW3mwQMtLMuwceN4/TXfttXGGJhfMNTrNbb52VzHO/f/ULtNiYjI\nWDLGsHcx+xLRTryR7QNyKMdNSIMSsUkI09GsfRPpt/H6qVXkNM21emecbWLPXIuZKrjO+FTOAM44\nI/vY7dsX8S8PfpNauw7AN/d/W7tNiYjIWKrXa9wxew8AcRQShZ0RtwhcJyXtZOvO2ommNspkOKFw\ndvvtt/PqV7960G0ROW29M86sqEKUGDZuHHGDjmLbGQ62DQ/t9SgUS2zfuI2SW2Q2WKBYKo26eSJy\nkjRGSl7E3Z8aPXv0a84gm9bYO+uspXAmE8Jd6wWf+tSn+OIXv0ilUhlGe0ROS++Ms1YtGzg2bRhl\na46uWLB43KNd7vxhzP0PuWw5w2JHZTs/Wr6fxWh51M0TkZOgMVLyJEgicMDBGXVTgCyc9XZsbCWj\nr+SJ9MOalbOHP/zhfOxjHxtGW0RO21x7HguL2dlsKuPGMQxnABec5+E4hh/+yCOODWdWtgFwoDM3\n4paJyMnQGCl5EplsXZdnjUc4c9xk5awzTWuUSbFmOLv44otxnPH4EIqsZa49z6biBh7a18K2YcP0\nqFt0dKWSxSPPjggCizvvitlePgMLiwOBdmwUWU80RkqeRCbEGHDt8diywHUTTJDtHNmKFc5kMqw5\nrfFUbN1aHcRl14U89x1G2/9OFLAc1nnspkfx/f11Hra5xMxMwlQ1+x+36xp838X3PWzPxrIcCkUP\nv+ASpy5+waNQzA6u9nwXx3MpFD3C0MWyHKw0wXacVe/tvQ5Yea3ne/i+d+R7E8PUVJFSKWvP+T9p\n8cBuw/fujHnyT8+wfWorexsHCGjg+1Wmp6exrPHazORY9Pc+3/2XE5f3vyt57v9677vvpyR2DLFP\nuWSvGksPHx+tNFkZV600O9fTL7hHPHf4GHy85w4dc3vPFYtA4mEbl8B02LKlyszM+P0+r/c/+9OR\n576fqhMOZyezvffsbP2UGrPebd1azW3fYfT9v3t/tovU0qzN/9/evQdZVd0LHv/ufd7P7qZfPAVs\nSdREMeArNxeMudepIWHmxlEMiJBK8kcqlYcT8/SfxD803qmpVM0dR3NN/lBnbt0yUefO3MS5lWhU\nUjE+sBGkAaGhoWmapl/n/dhnv9b8caClBRpsaPZ5/D5VVjW996nzW/76rHV+e6+1tuUobMtgcrKI\n0mT1qYUAABq1SURBVKp/5oZhYJo2ut9Cd31omkslYGFWbCzTxqxYVALVKRuWaeMqjYpRPa5pLlbF\nRve701576jxg6lzLtDBN64zX2q5NoWBg29XzTbPCVcsVew8EefPtAvMWtDLCGL9+600WsIw7brmK\nZLLFg/+TH43XefdaM7dfBt0PXOgY2ax/KyCflXpvey6Xx8FEWSFcx8I01dR4+OHx0TptXLUqNuGo\nb9p4ePqYa37o3HMdO33MPTW+Vj93Gj4nQsEqMT6ewzRr467eKY2Q+9lq5rbD7MfIC/4Lrpcr+KJ5\nTRppAHxOdS5jeyLkZTgX5MplNqEQ7Ntv0xHsBCAbmCASlc0FhKgnMkaKRme5Nspnoaxg9eHPNeBU\nHJpVfdaZIZuCiAZwQcXZokWLePbZZ+c6FiEuSqpSLc4qxeq0wWS09teB+H1wzcf9VEwYH47jUwEm\n3GF5ELUQdUTGSNEMMmZ1N2HNjlAjS87QdYVPZ2o7/VQl43FEQly8Gvl4CXHxxo3qA6iLmSC6BvGw\nhmFUMIzyyf8M8LDmUao6tXIqnopBuWxw5VIbXYf390PUSWBiMFo6Sj6fI5fLkstlpVgTQgjhqZRR\nLXx8TtjjSD6gaRCNglWsFmeTJy/SClHP5mRDECG8MFwaQcdHNhWkNRHCqpQYOpZmXqV6B61YyGPb\nFuDNg55ty+TA0RTRk89DGhkcwef307XIT3tnmPHRIKHJOHSnOEwfvRMWoXyEcqnE565aWxfrz4QQ\nQjSmE8UUAH63tpYMxOMwWogRAiaMlNfhCHHR5M6ZaAiWazNaHifqtuE6Gu3J6pW9QCBAMBgmGAwT\nCAQ9jhICgeBUPH5/AP/Jfy+5onq8OLYElEbZXyAajxFLxIhEo94GLYQQoumNFqqFT0CrseIsBsqo\njpMTFSnORP2T4kw0hOOFERzlEjCrT51ub6mdaRcXIhqDZEuFSjlB0E5QpkjZloXNQgghasPEyWmN\nYV9tja/xGCcfRK3JnTPREKQ4Ew1hMHcMADNbnTKYCLmUS0WUl4vMPqLOruoDNFW2umvjifKYl+EI\nIYQQU7JWFuVqhAO1tSImHgfQCbhRKc5EQ5DiTDSEgdRhALKpCD6fS0rvZ6iyH9u0PI7swsUTFprm\nYIx3A1KcCSGEqB0llUeZEcKh2rroGT/55BmfFaXslClYRW8DEuIiSXEmGkJ1MxAdK5+kpRUi8SjB\ncG1NvTgfTYNwuIyZTxJUEcaMSWzX9josIYQQTa7imNiagapECNZocaaM6g9jpQkPoxHi4klxJuqe\n5ViMlscJunFAp6XV9TqkWYuEq1MbA5UOXOUyWhr3OCIhhBDNLmWc3KLejOD311ZxFghohEJgFaub\ngozJuCnqnBRnou4dL57AUS5aOQlAS2ttDRwfRSRcAsBNV9edDRdGvAxHCCGEYLJ8ahv9MJrmcTBn\nkYjrGHm5cyYagxRnou4dzVc3AzEyLeg+RSxWv8VZMFjB53MpjrcT1AMcL5yQB1ALIYTwTL6Q5419\nuwHQrADZXMHjiM6UiGs4pVPFmdw5E/VNijNR947mhgEop1pJJBy0Ov6r1jSIxUwqho/2QCdlxyBr\n5bwOSwghRJMqlw0mzOqUe78KYdm1t3QgkdDACuHXAoyV5c6ZqG91/DVWiKqh/DF8+FDlGIlk/W+g\nEYubAASNLgBOGHIVUAghhHeKbh6AIAGPIzm7RFwHNGJakrHSBK6qvQJSiAslxZmoa5ZjMVw8Qdhp\nA3QSScfrkC5a/GRxZkx0oKMxWpHiTAghhHfKKo9ydYL+2vraqBQYhkEoVB03A3YCy7XIVmTGiahf\ntfUpE+IjGi6M4CoXMx9H1xWxWP0XZ+Gwjc/vMj7qpzPaQdbKkTPzXoclhBCiSZl6EVWJEArW1hpo\n2zI5cDTFZL46RrrlCCCbgoj6JsWZqGv94wcBKKbihEJlHKd+Hjp9LpoGyaRDsQTt/uoDqQ/kDnkc\nlRBCiGZUcSq4uomqRAgEau8CaCAQJB4PAQqrdLI4K8uME1G/pDgTdW2oeBwAt9hCLGZ6HM2lk2ip\nrp3TCtV1Z/uzUpwJIYS4/DJmFqBmizMATYdg0KGYDQNy50zUNynORN1yXId92X50J4wqx4lEDa9D\numSSLdUBMDUaJuqL0J8dwHbrf7MTIYQQ9SVtVYsznx1Dr+FvjdOLM7lzJupXDX/MhJjZgcwhSnYZ\nMvPx6TqRSMXrkC6ZSMQlHIITJxTdoU4qrsmhzBGvwxJCCNFk0pVqcRZSMY8jmVkw6IATJOKLyp0z\nUdekOBN1a8foewCUx7roaAnW9BW9j0rToKsLSmVFQnUC0De5z+OohBBCNJvxUhqACHGPI5lZMFSd\ncZLwtTJhpHDc2pyCKcT5NNDXWdFMHNdh10QfYS2Km2+jqzXkdUiX3PzqcjPMTBsBPcCeyf3eBiSE\nEKLpTJ68cxbzJTyOZGbBYHXqf0glcZXLeHnS44iEmB0pzkRdOpA+RNEqEassBrSGLM66TxZnE+M+\nrkosY7Q0xoQMNkIIIS6jjJVFOT7igYjXocwoGKzeKQtU5gEwkD3iYTRCzJ4UZ6LuKKV4c3g7AJlj\n7QR8GkEMVG09fuWixWMQi2qMTcDHWnoA6Jt43+OohBBCNAvLtSk4OVQlQizs9zqcGYVCDroG+fHq\nHb4D6QGPIxJidqQ4E3Unk02za3IPQS1EYSJJNF7iYHYXttU4W+kDaJrGgvk6FVNjnrsYgF3jfR5H\nJYQQolnsT/XjajZOtoN4pLaLM59PsbgjzLGjGjF/jP7MIVSjXbUVTUGKM1F3DuWPYCmbpFoIaLS1\nQygS9jqsObFgvg+AkVHoaVlGf2aAlJH2OCohhBDN4N3x3dUfMvNJ1HhxBtCzIIqroCuwmEwlK+vO\nRF2S4kzUnV2pvQDYk/MBSLY27vO/FsyvfkT7h/PcsmA1CsX2E+96HJUQQohG57gO743vQZkhklon\nmqZ5HdJ59SyIVn8oVNed9acPeRiNELMjxZmoK8cLJ3h3cjcxX5TJoSSxKIRCjTttIR7TiccUB4/n\nWdlxHX7dz1snemWqhhBCiDnVnxmgZJdx0t3Mi9fHpluLO8KEgj7Gj1WfyXYgI8WZqD9SnIm6oZTi\nhf7folAs9V2NaWrM764+E6yRdXWAYbqMT9qs7PgEo6VxBvNDXoclhBCiwfQP9tN3qI++Q338sf81\nAJxUN63xgLeBXSCfrnH1klbGR33E/XH607LuTNQfKc5E3eib3Mf76X5WJK+kONYBwIL5Hgd1GXRX\nn0HNjgPj3LJgNQBvjezwMCIhhBCNaDg3QiZcIB3Kc6g8iE+FcPNttMWDXod2wa5dNg/QmKcvJGvm\nGStPeB2SEB+JFGeiLtiuzf/u/x26pvP5xX/DsWENvx8WNkFxtnA+REM+/rTrOD3JK0kGE/SO7sRy\nG3etnRBCCO9MlFNUXBN/YT4+3UcyVvt3zpSCfD7P0s5qrGY6Cci6M1F/pDgTdeG1Y68zVp5gzaJb\nGTiQo1DU6OywmUynvA5tTigFhmFgGGUcp8zqq5LkSxbb901wU/enKNol9kzKM8+EEEJcekOFYQCK\nI510JANUygVKxUJNP0/UNA3++OZ+du0fIhTQOHa4uovzASnORJ2R4kzUNKUUbx/t5f8e+jcivjBr\nOm7hjb3Vgqy1QzGRNTyOcG5YZoUDR1McHM4ycDyHU8mia/DyO8e4ef4qAF4ffsvjKIUQQjQa27U5\nmj+GDx9Wdh7BaIEh+wD9mZ1YVsXr8M7JLBuM+Y8x7h8inqxglWJE9Aj9mQFZdybqihRnoqbtGt7N\n/zr4HBqwqvU6do3u5XjeQddd2uY1dmcbCAQJBsMEg2Fak3Guv7KNobECxXSEFa1Xsje1n72T+70O\nUwghRAPZmzpA2TZoMReB0uno9BOORevieaKhcJhwLEp7pw5oBMqd5Mw8rx+Xi5mifkhxJmrW4exR\nnjn4GxSKNYs+zRUdSygaYRwnSEuLgd4kf71KKcqlAjdeVd0a+N/ePMy6hZ9DQ+P5/n/FlrVnQggh\nLoGCU2Rf6gARf5hAehkAbS3exjQb7Z0u4bDF2PtLCWphnt3/L7x19B1yuezUf3I3TdSq2n/cu2ga\nR44dYTh3HIBjxgh/yfZiK5sb225gQawbgIEjDgCtrQYQ9yrUy6pSLjNo7COpFWhribB7IE3nwuPc\n0rmKN8d72XbsL/zNFWu9DlMIIUSd21XYi6tcPtV5PTv3hgBFSwsU62x5t88Hy5anGT6ygHzfDUQ/\n8Q7/dPB5PtNxE23BVsqlEp+7ai3JZB1WnqLhNcm9B1EP8kYeN6mxx+7nT5m3UChWBq5mYaRamLmu\n4sigja47xBOmx9FeXsFImFgixvWfDKHQ2NEXYW3XXxPzR/l/h18iW8l7HaIQQog61jexj+PmKF2R\nDpbEF5LN68RjCn+dXsYPBl02374Qf2Ue1qGVuLi8MfkOO/N7mCBD2S57HaIQZyXFmagZjnJ4Y+Qd\ndk3sIeqP8LdXrGW+r2Pq+NFjDmUDEvF8wz94+sNO7d64cIFJz3JIZzT+54uH+UzbTRhOhWf3vkA2\nm5FpGkIIIT4y0zF5rv9f0dBY3b2SQgFsR6Ml6Xod2qwppUgGLe7722XY6S7MQ9eDE+Bo/hjvZnbz\n9+89xp8O/0WmOYqaI8WZqAklq8Tvxv7IYH6I9vA8/t3S25kXbps6btuKt9+x0DRobcl4GKk3bMvk\nwNEUh47naOvOEQ7lOTxmsu0Nl6AT5r30Xn6x+ylyuazXoQohhKgzvx34PRPlSVZEltMaamFo+OQS\ngqTjcWSzVykbvH1iB8XwIVZflyJU6ib/zm3oA7cSzc/Ddh1+ffj/8Iu9T/NS/6vk8zmvQxYCkOJM\n1IDJcpqf7/gFI5VRlsQX8rkla4j4p+8KtXO3RaGouO5aP6Fgc01pPOXU7o3hcJjuruOEQgaZVAz3\n8M3EfXGGzBF+O/QHufonhBDigh3ODvLq0J/pinTwydjHsSzFrt0Wfp9i8aL63nAqHIkQS8To7gzw\nVzfmmNduU5xoJXfwBj6pX0drqIWjpWHeKb5HzpTlAaI2SHEmPHUgfZD/2vsYJ4qjXJ+4hs8svAW/\n7pt2Ti4PfXtt4jGNldcFPIq0tui6YuHCEdrmueRTUVI7biboxnlj7B3++f3nKZhFr0MUQghR4yzH\n4p/2PYdCsfmaDfg0H337bIwK9Cw1CQW9jnD2lIJKxcAwylRMA8cx6FlRZtESB8sKsnNHF9eHb2ZZ\nfAk5u8A/7nuGoxNDMsVReO68yzyVUjz00EPs37+fYDDII488wpIlSy5HbKKBucrl90de4cXDL6Fp\nGhtW/B3tlSRZbXpRoRT07gLXhVtvChAINNlisxn4fC7XXmczctxmeDBCdudNhK/Zzl9GtrP9xE7+\nev7NrOm+ha62brRmW6QnxGUg46OoZ65y+ZdDL3KiNMZti/+Kq1qXc/DgUXbvsQmHoGepBdTv2GFb\nFQ4PGxSJkRnPY1RKtHcFWd7jUC5Mkkp38YeXFe1dVxFpL5GOT/Lf33uKlc5t/N0tN0zt5KiUOmPK\nYyKRlHFVzJnzFmcvv/wypmny7LPPsmvXLh599FGeeOKJyxGbaECuctmX6uelwVfpzwzQFmrlq5/c\nzJUtS9l9cPe0cydTLn/aGSGd11iy2McVS+p0y6g5pGmQTBYJrsgxdjxMZu9NaB0jqAUDvDryOq8e\n/wtXxBfxyc6rWdHWwxWJRYT9tf8gUSHqgYyPol5lKlme2ftrDqQP0hFp5z9euQ6A3QddbBtu/FSg\nbndpPJ3/5HKAYDCEc9ozQVuSKYIhl0y2m4nREIzeSGTxAOWF/bzJ79n32kFuaLuZaxcvpCup8crh\nV3ACLgW7RNbI0jNvOVd3fYzlySuIBqIetlA0ovN+9Hp7e1mzZg0AK1eupK+vb86DuliO61B2jPOe\nZ1QcHPfS3bq2/DapXLNOJ1OYPutk+xWO6+AoF1vZ5MwcmUqW8fI4707sImtWN624pu1q7u75T0T0\nCKmcQabgcCLnkM8rxsYd+g85KOVjySLFmk/X8dyKyyASCbBwYYmFiw2KpW6yQ0vIBY7im3eCQY5x\ntHgMjrwMQFugnQXRBXRG25kXaqMt3ErYHybsCxLUg+iajq7paJr+oWumZ79K2Nx/95eu/X6fRij4\nwZTeiC+M70NTfEVtqcfx8XIxHRPTtab9LlzRGG3STYsubz+pTvtJ4SoXpRQlu8yEMcFYeZw/j7xO\nyS5zdevVfGHxf2D7nkl6D4zz/qAiHtP4+Ao/6fHLFK5H4vEiPR+3mBzXGDzs4J7oQXfCqPn7ycb3\n81rlAK/ujaKFymj69F0rB0eGeWXkzwB0RTpZmriCpYmldITbCfvChP1h/JofTQPTb5HOlbxooucu\n9/eDaNh3zjuaQT1A0Fcf3yXPW5wVCgUSicQHL/D7cV0XXa/d5Wr/7d1/ZCA76HUY4iyU48OZXIw9\nvoQdxSQ72HHOcxNRxbWLi3Qv03AsKJ4c541CHtu2wLUoF/NUKjZGoYhRKaPpPhyremKxUMDv95NL\nV5+eWSoWKBYKgEL3+0DTcSyLQj4742tLxQJoOnalgq7r01576rzsZArTtKa9x6n3PPXayxVfKGqy\n9uYOTGsZBwcWkZ+IMGYdR4un0WM5UrEMaWsSmvM7Ut1YHF/Igzf/Z6/DEDOox/HxcjheOMF/2f4P\n2Kp+d/prdMrVsY5ey7tjS3iXD2attMYVn7rWwUgVMDIGWlRHuYpCeoJiNofjOphaANfvmzYefnh8\nPH3MsysVLDuI46izjoenj3EzHTt9zM2lUzMeAy7o3HwmRTAAHe0TdNlL6F5wDaXyIuLdBr3Z3eS1\nHKqSwCxGUJUoyoihjCj4bPR4Bj2eZtRJMVYeZ/tYr1fpFBfAr/t58Kb7mR/r9jqU8zpvcRaPxykW\nP6h6L2Tg6exMzHh8rv39v/+xp+8v5tZ9a77odQhCCFGX4+Pl0NmZ4J+X/w+vwxDiIn3J6wBEkzrv\n5b1Vq1axbds2AHbu3MnHPvaxOQ9KCCGEqHUyPgohhLjUNHWe/UJP340K4NFHH2X58uWXJTghhBCi\nVsn4KIQQ4lI7b3EmhBBCCCGEEGLuNfeqZSGEEEIIIYSoEVKcCSGEEEIIIUQNkOJMCCGEEEIIIWrA\nJSnOXnrpJb73ve+d9dgjjzzCXXfdxdatW9m6dSuFQuFSvGVNman9v/nNb7jrrrvYuHEjr7322uUN\nbA5VKhW+853vsHnzZr7+9a+TTqfPOKfRcq+U4qc//SkbN25k69atDA0NTTv+yiuvcPfdd7Nx40ae\ne+45j6KcO+dr/9NPP8369eun8n3kyBFvAp1Du3btYsuWLWf8vtFzD+duezPkfbakn5R+UvrJDzR6\n7qF5+0nbtvnhD3/I5s2bueeee3jllVemHW/k3J+v7bPKvbpIDz/8sFq3bp164IEHznp806ZNKp1O\nX+zb1KyZ2j8+Pq7Wr1+vLMtS+XxerV+/Xpmm6UGUl95TTz2lHnvsMaWUUi+++KJ6+OGHzzin0XL/\nhz/8Qf34xz9WSim1c+dO9Y1vfGPqmGVZ6o477lD5fF6ZpqnuuusuNTk56VWoc2Km9iul1Pe//321\nZ88eL0K7LH71q1+p9evXqy996UvTft8MuT9X25Vq/LxfDOknpZ+UfrKqGXLfzP3kCy+8oH72s58p\npZTKZDLqs5/97NSxRs/9TG1Xana5v+g7Z6tWreKhhx46V+HH4OAgP/nJT9i0aRMvvPDCxb5dzZmp\n/e+99x6rV6/G7/cTj8dZtmzZ1JbL9a63t5e1a9cCsHbtWt54441pxxsx9729vaxZswaAlStX0tfX\nN3Xs0KFDLF26lHg8TiAQYPXq1Wzfvt2rUOfETO0H2LNnD08++ST33nsvv/zlL70IcU4tXbqUxx9/\n/IzfN0Puz9V2aPy8XwzpJ6WflH6yqhly38z95Lp167j//vsBcF0Xv98/dazRcz9T22F2ufef/5Sq\n559/nmeeeWba7x599FHWrVvH22+/fdbXlEoltmzZwle+8hVs22br1q1cd911dfmgztm0v1AokEgk\npv4djUbJ5/NzGudcOFvbOzo6iMfjAMRisTOm4jRS7k/5cD79fj+u66Lr+hnHYrFYXeZ6JjO1H+AL\nX/gCmzdvJh6P881vfpNt27Zx2223eRXuJXfHHXcwPDx8xu+bIffnajs0ft4vlPSTVdJPSj8p/eSZ\nGj3vkUgEqOb5/vvv57vf/e7UsUbP/Uxth9nl/oKLs7vvvpu77777Iwe8ZcsWQqEQoVCIW2+9lfff\nf78uB57ZtD8ej08bjIvFIslk8lKHNufO1vZvf/vbFItFoNqu0z940Fi5PyUej0+1GZg24DZKrmcy\nU/sBvvzlL099Eb3tttvYu3dvQw0+59IMuZ9Js+b9w6SfrJJ+UvrJs2mG3M+kGfI+MjLCt771Le67\n7z4+//nPT/2+GXJ/rrbD7HI/p7s1Hj58mE2bNqGUwrIsent7+cQnPjGXb1lTrr/+enp7ezFNk3w+\nz8DAACtWrPA6rEti1apVbNu2DYBt27Zx4403TjveiLk/vc07d+6c9gWqp6eHwcFBcrkcpmmyfft2\nbrjhBq9CnRMztb9QKLB+/XrK5TJKKd588826z/e5KKWm/bsZcn/Kh9veTHmfDeknpZ+UfrKqGXJ/\nSjP2kxMTE3zta1/jBz/4AXfeeee0Y42e+5naPtvcX/Cds4/i6aefZunSpdx+++188YtfZMOGDQQC\nAe688056enrm4i1ryunt37JlC/feey9KKR544AGCwaDX4V0SmzZt4kc/+hH33nsvwWCQn//850Bj\n5/6OO+7g9ddfZ+PGjUB1Wuvvfvc7yuUyGzZs4MEHH+SrX/0qSik2bNhAV1eXxxFfWudr/wMPPDB1\nF+DTn/701FqbRqNpGkBT5f6Us7W9WfI+G9JPSj8p/WTz5P6UZuwnn3zySXK5HE888QSPP/44mqZx\nzz33NEXuz9f22eReUx8u8YUQQgghhBBCXHbyEGohhBBCCCGEqAFSnAkhhBBCCCFEDZDiTAghhBBC\nCCFqgBRnQgghhBBCCFEDpDgTQgghhBBCiBogxZkQQgghhBBC1AApzoQQQgghhBCiBkhxJoQQQggh\nhBA14P8DbWzSXg3uzLAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11df15470>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(15, 5))\n",
    "plt.subplot(1,2,1)\n",
    "sns.distplot(trace_nuts['mu'], label='NUTS')\n",
    "sns.distplot(trace_advi['mu'], label='ADVI')\n",
    "plt.legend(loc=0, fontsize=15)\n",
    "plt.title('mu', fontsize=15)\n",
    "plt.subplot(1,2,2)\n",
    "sns.distplot(trace_nuts['sd'], label='NUTS')\n",
    "sns.distplot(trace_advi['sd'], label='ADVI')\n",
    "plt.legend(loc=0, fontsize=15)\n",
    "plt.title('sd', fontsize=15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}