lectures/Week-4-Thu.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### ASTR-598, Winter 2018, Connolly & Ivezic, University of Washington\n",
    "https://github.com/dirac-institute/uw-astr598-w18/tree/master/lectures/Week-4-Thu.ipynb\n",
    "\n",
    " \n",
    "### The objectives of this lecture\n",
    "- A gentle (re)introduction to MCMC (from a mathematical and intuative perspective)\n",
    "- An overview of probablistic programming"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### For further reading\n",
    "\n",
    "These lectures have been adapted from a range of different materials including\n",
    "- [A Practical Python Guide for the Analysis of Survey Data](http://press.princeton.edu/titles/10159.html) Chapter 5.  \n",
    "- [Thomas Wiecki: ``MCMC sampling for dummies by Thomas Wiecki\"](http://twiecki.github.io/blog/2015/11/10/mcmc-sampling/)\n",
    "- [Probabilistic Programming and Bayesian Methods for Hackers by Cam Davidson Pilon](http://nbviewer.jupyter.org/github/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers/blob/master/Chapter1_Introduction/Ch1_Introduction_PyMC2.ipynb)\n",
    "- [Abraham Flaxman's (UW) discussion of MCMC](https://healthyalgorithms.com/)\n",
    "- [Astro 324 special topics course at the University of Washington](https://github.com/uw-astr-324-s17/astr-324-s17)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### For those who want to dive deep:\n",
    "- [Andrieu et al. ``An Introduction to MCMC for Machine Learning\" (includes a few pages of history)\"](http://www.cs.princeton.edu/courses/archive/spr06/cos598C/papers/AndrieuFreitasDoucetJordan2003.pdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Software needed\n",
    "\n",
    "For the examples we will be using PYMC3 \n",
    "\n",
    "- pip install git+https://github.com/pymc-devs/pymc3\n",
    "\n",
    "To demonstrate how MCMC techniques sample we will use Jake Vanderplas's animation code\n",
    "- pip install JSAnimation\n",
    "\n",
    "Aside: PYMC3 should install Theano (a deep learning package) by default as it is used in manny of the definitions of distributions. We will call Theano a couple of times in the later exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Testing it all works \n",
    "\n",
    "If you want to ensure that you have everything in place try to run the following piece of code (we will walk through what these lines mean as we go through the class). You should see something like the following figure\n",
    "\n",
    "<center><img src=\"figures/pymc3-gaussian.png\" alt=\"pymc3 test\" style=\"width:600px;\"/></center>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from scipy.stats import cauchy\n",
    "import scipy.stats as stats \n",
    "\n",
    "from JSAnimation import IPython_display\n",
    "import pymc3 as pm\n",
    "import theano as thea\n",
    "import theano.tensor as T \n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 445.119047\n",
      "         Iterations: 11\n",
      "         Function evaluations: 16\n",
      "         Gradient evaluations: 16\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 10500/10500 [00:05<00:00, 1871.75it/s]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean</th>\n",
       "      <th>sd</th>\n",
       "      <th>mc_error</th>\n",
       "      <th>hpd_2.5</th>\n",
       "      <th>hpd_97.5</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>mu</th>\n",
       "      <td>10.065646</td>\n",
       "      <td>0.153347</td>\n",
       "      <td>0.001589</td>\n",
       "      <td>9.771033</td>\n",
       "      <td>10.366624</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sigma</th>\n",
       "      <td>2.194541</td>\n",
       "      <td>0.109839</td>\n",
       "      <td>0.001060</td>\n",
       "      <td>1.982834</td>\n",
       "      <td>2.413016</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            mean        sd  mc_error   hpd_2.5   hpd_97.5\n",
       "mu     10.065646  0.153347  0.001589  9.771033  10.366624\n",
       "sigma   2.194541  0.109839  0.001060  1.982834   2.413016"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAEYCAYAAABBWFftAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXe8FOXV+L+HC1IUEKSjoCKgiA0VUVQudixgiy15LW8S\nTdSoSUxi/CURY968mqrG5I0aY9BELNiwRUW5sWAQBUREKaKASDOK9H5+fzwz3Lmzs7uzZbbce76f\nz352dvaZmfPMzs485zlNVBXDMAzDMAzDMAyjcJqVWwDDMAzDMAzDMIzGgilYhmEYhmEYhmEYRcIU\nLMMwDMMwDMMwjCJhCpZhGIZhGIZhGEaRMAXLMAzDMAzDMAyjSJiCZRiGYRiGYRiGUSRMwTIMwzAM\nwzAMwygSpmAZhmEYhmEYhmEUCVOwDMMwDMMwDMMwioQpWIaRICLykYhcKyLviMhqEblbRLqIyLMi\nskpEXhCR9iIyTEQWRWx7TLlkNwzDMJoG9qwyjOJiCpZhJM+ZwLFAP2Ak8CxwHdAJqAGu8tppWaQz\nDMMwDHtWGUbRMAXLMJLnD6r6maouAV4FJqvqDFXdBDwOHFRe8QzDMAzDnlWGUSxMwTKM5FkWWF4f\n8Xmn0opjGIZhGCnYs8owioQpWIZRGawFdvQ/iEgN0Ll84hiGYRhGCvasMowYmIJlGJXBHKCliIwQ\nkebAT4AdyiyTYRiGYQSxZ5VhxMAULMNIlnAwcGRwsKquAq4A7gE+AVZ774ZhGIaRNPasMowiIqrJ\nJYMRkZbAK7jZjebAOFW9MaLd7cAInOn5YlWdnphQhmEYRpNERO4BTgWWqer+3roOwENAb+Bj4BxV\n/TJi2+8CXwe2Ae8Cl3jB/4ZhGIbRgEQtWKq6ERiuqgcBBwIjRGRwsI2IjAD6qGpf4DLgz0nKZBiG\nYTRZ7gVODK27Dpigqv2Bl4EfhzcSkR7Ad4BBnmLWHDgvYVkNwzCMKiVxF0FVXecttsQ9lMIms1HA\nfV7byUB7EematFyGYRhG00JVXwO+CK0eBYzxlscAp6fZvAbY0Ys7aQN8moiQhmEYRtWTuIIlIs1E\nZBqwFHhRVaeEmvQEglXBF3vrDMMwDCNpuqjqMgBVXQp0CTdQ1U+B3wILcc+olao6oaRSGoZhGFVD\n86QPoKrbgINEpB3whIgMUNVZue5HRKxyuGEYRpWgqlJuGfIk5VkjIjvjLF29gS+BcSJygao+ELUD\ne14ZhmFUD0k8r0qWRdDLPDMROCn01WJgt8DnXb11UftoNK8bbrih7DJYfxLsz9VXu1cF9MV+n8p/\nNbb+VBnLfLd0EekGLI9ocxwwX1U/V9WtwGPAEZl2Wu7foJJfje16t3Nk56hSX3aOsr+SIlEFS0Q6\niUh7b7k1cDzwQajZeOBCr80QnOvFMgzDMAyj+Ij38hkPXOwtXwQ8GbHNQmCIiLQSEQGOBd5PUkjD\nMAyjeknagtUdmCgi04HJwPOq+qyIXCYilwKo6rPARyIyD7gTuDxhmQzDMIwmiIg8AEwC+onIQhG5\nBLgZOF5EZuMUp5u9tt1F5GkAVX0TGAdMA97BKWh3laELhmEYRhWQaAyWqr4LDIpYf2fo85VJylGJ\n1NbWlluEomL9qWysP5VNY+tPpaKqF6T56riItktwNbP8zzcCKXUcjdyx6z07do6yY+coO3aOykei\nhYaLiYhotchqGFxzjXu/9dbyymEYZUBE0OpNclEw9rwyDMOoDpJ6XpUsyYVhGIZhGIZhGEYmNm+G\nLVvKLUVhmIJlGIZhGIZhGEZF8NRTMHFiuaUojMTrYBmGYRiGYRiGYcRh48ZyS1A4ZsEyDMMwDMMw\nDMMoEqZgGYZhGE0CEblHRJaJyIzAug4i8oKIzBaR5/3ajRHbtheRR0TkfRF5T0QOK53khmEYRjVh\nCpZhGIbRVLgXODG07jpggqr2B14Gfpxm29uAZ1V1H+AArNCwYRhFYuVKGDu23FI0HjZtKn+SDFOw\nDMMwjCaBqr4GfBFaPQoY4y2PAU4Pbyci7YCjVPVebz9bVHVVkrIahtF0WLeu3BI0Lp58Ev71r/LK\nkKiCJSK7isjLnjvFuyJyVUSbYSKyUkSmeq+fJCmTYRiGYQTooqrLAFR1KdAlos0ewGcicq/3nLpL\nRFqXVErDMIw82bSpcSSOiMuWLbBmTXllSNqCtQX4nqruCxwOXCEie0e0e0VVB3mvXyQsk2GUhG0K\nVmvUMKqOqH9tc2AQ8EdVHQSsw7kWGoaRMGPHNv5nadL9e+EFeOaZZI9hNCTRNO3ebOBSb3mNiLwP\n9AQ+CDUtegVlwyg1X34Jjz8Ozz4LI56HVavhmtth551h//1h+HA44ww44IByS2oYRoBlItJVVZeJ\nSDdgeUSbT4BFqvqW93kc8KNMOx09evT25draWmpra7MKMnYsDB0KvXrFlNwwDCMGa9fCtm3llqK0\npFNa6+rqqKurS/z4JauDJSK7AwcCkyO+PlxEpgOLgR+o6qxSyWUYhaAKb74Jd97plKvhw2HkSBjZ\nGtq3h6tugxUrYPp0N4M0ciR07w5XXAHnnAMtW5a7B4bR5BAaTuqNBy4GbgEuAp4Mb+ApX4tEpJ+q\nzgGOBTI+p4IKVi6sXGkKlmHEYds2WLXKTWKmY+NG1661OfQaHuEJrxtvvDGR45REwRKRnXAzfler\natgr8m2gl6quE5ERwBNAv6j95DMjaBhJsGIF/OMf8Ne/uuDUSy+F2bOhix+9Md17F7fuhBPc65Zb\nnJn+j3+Ea6+Fq66CK690yphhVCulmhH08eKfeqnq7By3ewCoBXYRkYXADcDNwCMi8t/AAuAcr213\n4G5VPdXb/CrgHyLSApgPXFKMvhiGkR/z58OUKXD++enbvPyym7TI1MYwkiBxBUtEmuOUq/tVNWpm\ncE1g+TkR+ZOIdFTVz8Nt850RNIxi8OGHzgr13HPwyivOGnX77XD00dAsZjRjTY1n4RoJs2bBzTdD\nnz5w+eVwzTXQsWOyfTCMJCjVjCCAiJwG/AbYAdhDRA4Efq6qI7Ntq6oXpPnquIi2S4BTA5/fAQ7N\nS+gMfP65/e+NxsV//uMyuJ15ZrLH2brVvX/5JXzyCey7b2qbUiR2WL0aZs6Eww9vuL6S/tvVGMO2\ncaPz/jnwwPy2L3efS5Gm/a/ALFW9LepLEekaWB4MSJRyZRilZutWeO01Z2nq2xeOPBL+/W847zxY\nsADuuw9qa+MrV2EGDHD7mDwZPv3UHeO66+CLcBJpwzCCjAYGAysBVHU6LstfVfL8884ini9Ll8LC\nhcWTp1hs2NC0spZFMXNm+WvxFJONG501KBuffVaa394fQM+dCzNmZG6bJJ9+Ch9/nLr++eerJ/36\nG2+4cIdK4/0qrjaYdJr2ocBXgWNEZJqX3vYkEblMRC71mp0tIjNFZBpwK3BukjIZRjZWrICf/Qx6\n9HCxUjvtBA8/7G6iY8bABRcU16WvTx/4y19g2jQ349W/P9x2m0urahhGCptV9cvQuiqcn61HFd57\nL79tJ02C118vrjzF4Omn3QCzKfPuu07ZSJK1a2Hz5tT1r7ziBs3F5I03nAeHEZ+XXqpPLlFui0om\nPv7YuVzG4fMKMIG89VblT14kqmCp6uuqWqOqB6rqQV4a9n+q6p2qepfX5o+qOtD7/ghVjUqCYRiJ\ns3Ur3HEH7LMPLFvmHlDvvAOjR8NBB4EknOuyVy+46y7nM/7cczBwIJQwrMUwqoX3ROQCoEZE+orI\nH4BJ5RaqUD78MPP3S5cWZunKh5Ur878Hbd7srFhGbmzblttAfPz4aMvD4sWwaFHx5IJ6l7xKI+ln\nM+Q/4blmTb0C/OCDbqK2mlm5sjImTubOdQlOgqxcWVlKbClcBA2j4lmxAk48ER56CF591WUF7N+/\nPLIMHAj//Cf87nfw1a/CD35grjaGEeA7wL7ARmAssAq4pqwSRbBkSfy2IvUDg3SDxYkT3aRPUixY\n4NLE+2zZ4gaD6fqxfLlLMJAvq1alDpAaI8uX53b/fuIJZ5XMhUrzdvCv5YcfLt4+H3wwv8FzsQbc\njz5anAmDtWsL30chFDpJU0mp3sO/7XPPRbtqlgtTsIwmzwcfwKGHulddnbNgVQKnnuosaPPmweDB\nLimGYTR1VHWdqv4/VT1UVQ/xlivOVjJ1avY2QYXGp1wzsF+GnC4fecQpB1HfgXMnmjcv/+M980xx\n3c02b3YTUwATJmS2CG7YUDo3p/fec/fxbHzxhYs3iRvnFKSSZu2DFNvilc0lbN684rtFpjv+2LHx\nFedS/z4rVqS3Mk2YUFpZSk3wmiv3/8IULKNJM2MGHHOMcwP83/91Wf4qiU6d4LHH4DvfcQk1Xnqp\n3BIZRnkRkYki8nL4FXPbe0RkmYjMCKzrICIviMhsEXleRNJGWIpIMy+WeHwx+uKjmv9gIO4gVjX3\ngHt/xv7ZZ52Fq5JZu7Y+QdCKFc5FLh1TphTHzWnNmnhJieK4sM2Z4zKmNWZWrEjWAjJvXqoFI0n3\nwajYt2Lz4ovxJlcXLHD/U3DuxOEJhDj3l3IrJI0NU7CMJsvcuc4t8Pe/h4svLrc06RGBb3zDzShf\ncAH87W/llsgwysq1wA+8109xVefeirntvcCJoXXXARNUtT/wMvDjDNtfTZYCw/nw2muwfn1+28YN\n9J4/H55MKZSSStANKjjgiuu2tnlzZSRCWLy4Xv5Vq+Dtt+u/K9Yg/+WX661m+eL/fqWe3Fu0yKVT\nz5eFC51lM67b3Pz5znoSN5FCFEFlKQnrUVTbFSvqLc3h7+fMST7BzGefOQuor8iPHRttfV2yJNrS\nXC0Uy801/Bs9+2z6WMSXXkr2nMVSsERkv+REMIzSs2SJU65uugnOrZK8lcOGudoiP/853HijzTYZ\nTRNVfTvwel1Vv4crHhxn29eAsM1hFDDGWx4DnB61rYjsCpwM/CUvwTNQihjLbAMY/37y+OOFHWft\n2swubhs3utIUkHw8h69ALFzoBsOff55bbNy2bU75TUfce3Cmdo88UnhcTKb9r1kTrQS99lr975AP\nr7/uBq9xvSpWr3bvxfrNH3us3qqai5Vq0qRoBeWLL1ycF7jz6StVmSyU8+fHK5FQjFjD4HFWr3aK\n+UMPFb7ffCnmf/fLL12MW7HHNKpu30uXNlw/dqxLZLZ8uctGmBRxLVh/EpE3ReTyTO4ThlENbNgA\np58OF17oLEPVxN57Ox/zcePghhvKLY1hlB4R6Rh4dRKRE4FCnktdVHUZgKouBbqkafd7nNUs1jBg\n1ar0s82ZKKZLUyGuh0nx2WfprRjbtmUfsM6bl9ldMRjXtmZNQ1fASZPiZ0VctMhNaEXNfm/bll6G\nKVNyt0YWmjxh+fL0CvRTT6UqQdkG/Js2OQ+POGzc2ND1dPnyhrWLJk3KnAxly5b83TXzOW8LFrhX\nOJtf8DcL/mfi/B+ztSmGa3/4f7xxo7sON2xw8Yx+8oyNG4t7D5kzJ9rd9qGHXHHnYuCf+6TSrkfd\nA33LlR9nmgSxFCxVPQpXz2o34G0ReUBEjk9OLMNIBlW47DLo3bt6FZSuXd0N+9FHnSXLMJoYb+Nc\nAt8G3gC+D3y9iPtPeRyLyCnAMq+osXivWKxZk7ou3wH1pk3prQ6rV6cqA+PHu+Lo0HCQsXmzU0Si\nZAtSjuB83+Xqww+jk4BMmRJ/1nnFingKbpQF6cMP62e+w1avd991ikNUTNu8eW52PEhwwDtjRup2\nIvEVmnRExeL5xw0PXJ95xr2ruvPzr381/H7RosznOHjdbNzY0PV05syG1/eCBfDRR6n7WLXKZexd\nu9bJsGZN/m6yUXz+eXr3vU8/Te1zOoK/3fr10f+JoPtpqfnyS3cufUUhTkxm1P8qiiefdH2bNi36\n+3nz6ieS4hIVt5Zvspywcu8Ttq6Vq7xA87gNVXWuiPwE92C7HThIRAS4XlUfi9rGc6m4D+gKbAPu\nVtXbI9rdDowA1gIXew8xwyg6f/mLy+7173+XpnZGUnTp4vz/hw93/fjZz8otkWGUBlXdo8i7XCYi\nXVV1mYh0A6LmNIcCI0XkZKA10FZE7lPVC9PtdNy40YCrUTRqVC21tbUsX14/m33OOfkJO38+HHZY\n6vpp09xM8/nnu89r17rBh+8mFyx46w+QZ81yGUrzIer+OX8+tGhR/7mQgU0m16zwIHfsWBgwILVd\n3IHbhAnOZbxjx+jv6+rc77V1q1NWsiUd2LCh4SDPl3f9epdVsEWLwrPVPvMMHH105jZxFORPPsm9\nNtNTT6X/Lu5z9dNP3bF9ZfSpp2DHHWHkyOz7Dg+go6yJvlWsVat4MqaT+4MP6pdfesn99zp0SC9j\nUoR/y1wzTeZDnKQ4uSrF48a5cIcePerXZavXNnUq7Lwz7Llnw/W+kuzf83zCrpMffwyHH17/edas\nOiZNqku8JlksBUtE9gcuAU4BXgROU9WpItIDN4MYqWABW4Dvqep0EdkJZ/16QVW3X7IiMgLoo6p9\nReQw4M/AkPy7ZBjRzJgB11/vasnsuGO5pSmcrl2dknX00dC5M3z72+WWyDCSQ0TOzPR9uom+qF3R\n0AI1HrgYuAW4CEhJBaGq1wPXe3IMA76fSbkCOPvs0QAMHeqKiENDl5pCawSpNnRvCbvxjPfyHPqx\nL0HiDoKjAsDr6lyfgmze7BS4yZNTB5+bN7tBFcBJJ6Xu74svXKKI888vbNIrTqY1f5AapXgEB+2P\nPprqcvfvfzv3xUwKgM+0aQ23//BDp1A9/bT7nEs/t26NToCxalVqkooNG1wMXXjAmY5Vq8o/0Ri0\naAQV8rFj4eyzoXnEKDU8gI5jEY76zbdudRMR7dpFb7NtW+r/Z+PG+sQmmc5dutT8y5dD376p6+PU\nb/KvqdWrM7vShi2omQgrqy+9BMceG3/bl0P5W1VdfNxZZ6XfLuqeBO732LAB2rZtuH72bPfuK1hx\nY8/STTAMGFDLoEG120tpPPZYMq5AcWOw/gBMBQ5Q1StUdSqAqn4K/CTdRqq61LdGqeoa4H2gZ6jZ\nKJyVC1WdDLQXka459cIwsrBmjUtm8dvfVk6dq2LQrZvL2HXTTfUDKsNopJyW4XVqnB2IyAPAJKCf\niCwUkUuAm4HjRWQ2cKz3GRHpLiJPF70XWYjrlvfYY6mDm7iEB4aPPtrQ7SuT9WjJktSB/dSp6WOb\ngoPmqCKrUTPx775bOrceX0kNuiZGxTP5robplLPwQDU8s/904EoK9y38ewQ/P/xwfPdB//y++mq8\n9kHefDN1kP/00/GSOEDuJQAy4V+LmzfXD67B/Ub5XhcPPljvGun/NrNmuXXr10craYVYOIKWryDp\nzmecRCe+jO++mzn7XXDiZevWzNauiRMbZsJMF5MUdU+IcvfbssX9f1Trz8HHHzdM5JOuRuBLLzX8\nn6SjkoodZyKui+ApwHpV3QquFgjQyiv4eH+cHYjI7sCBQNiDvCcQNBAu9tbloIMbRmauuAKGDHGJ\nLRobffrAE0/AKae4h0W+Lj+GUcmo6iVF2McFab46LqLtEiIUN1X9FxAzgiN34rrc5JrWODhoDw/o\nN21qOGDLlnp84sSG7jrBGJ9cXZeC/fCTI8ycGe2G5St2+dYfirI4+ErJunVuAPvmm9Hb+oP+qMB+\nf0Z9993r12UaBM6YAfvum1XcBmnmowgOWidPhp49G8oYjMFatAh22y39sT780L1WrKh3lVy92ik1\nr78OhxySOTmGH4cVdBH1yVUp8l0QVRvGvy1enD1uMBO+5cS/5vzr9oknottHKapx/p+bN6cq4nEs\nrFC/nUjhMWmbNrkkFe++m77N559nTy6xerW7Jxx8MLRp03DbML5SOnmyi73zE3Ptsw/slyUfuX99\nFaJAJZkVMFfiKlgTcA8g/9JuA7wAHBFnY889cBxwtWfJyovRo0dvX66tdT7thpGNMWNcYHSmTEbV\nzuDB8Ne/wqhRLv1unz7llshoStTV1VEXNz1bEfCSTuwLbI+wUNWfl0yAHHj99XoXwTj47j0rVrjB\n0dChhc2kQ+pgL2p2vZB6MMFZ+VwTYwRns6P6uWSJizmtqYEXXqhfP3Omi7vKJe16NuLInkm5C1qA\nsg0SJ02CI9KMoHxLhR9DFpbrjTfce/DcLVmS/lxs2uSeC0OHZldQ5s2DnXZKXb9yZTxXvEKL7wb7\nGvV75JqVMxPprEyZCFrU0uG7xAYJugzOn++U8WbNnPISTCM+frxz/x8ypGFykHwygm7bll552rw5\nWhkG5545fHjq+jjJPHwFJ+yy/P77xfntcjkHUf9Bf/t0VrRiElfBahVUjFR1jYi0ybSBj4g0xylX\n96tqVJnDxbjshD67eutSCCpYhhGHWbPg2mvdjGtjiLvKxGmnuWQXI0a4h3enTuWWyGgqhCe8bkww\nvaWI/Bk3yTccV5PqbCCN3aEyePNNF6MUNXBNx4IFTnHZd9/UjGeZXPiWLnWuw0Fmzaq3dqxe7SwV\nUJr4m6DytWkTtGwZ3S6ddaSuzllPwnEr777rBqm5WEeiBme51iAr1jlbsKBewQrvc+NGp8xEzcZv\n2hQvXieK119PfQ5G9SdKCcs305tPMDX/ypXRmek2bKjPepmOTN/7/4ugkldoCvxM5Fscd/Jkl7Sh\nY0en5AXPbTAxTXD/69YVVylYuDDzRGw+ytDWrann5L336pfjxoZlUuTSKcVRLqpR9wa/1lkpiBuD\ntVZEBvkfRORgIK7x8q/ALFW9Lc3344ELvf0OAVb6NUkMoxDWrXNxVzffDAMHllua0vDtb8OZZzpL\nVpIPFsMoI0d4CSa+UNUbgcOBfmWWKSPLlzsLUVQ9mTAiDQeIzz2X2iaTC9/Eianr1q2rj/EIxjjM\nm+fc1ZIkOFh66630FhQ/PiYdhVhG/G3DMWBbt6ZPQZ2OqJTjUeRSgDZKyQkOGH35VV28XCGEz0FU\nXFwS+M+jDz+MvqZ9wspjPqUC8rFM+YwdW/5n56pVzuUyWE9s0aLUgrmQXeFPN4Hw5pvOYpzu/KZL\n0pGJqMQ9+dxf0inzr70G09PkGM/mQFGKYu5h4lqwrgEeEZFPcdmXugHnZttIRIbi6me9KyLTcPVF\nrgd6A6qqd6nqsyJysojMw6VpL9jP3jAArr4aDjgA/vu/yy1JafnlL+GCC+CSS+CBB8qfJcowiow/\nubfOy2T7H6B7GeUpKmvWRLsYFUIm60M4YUUmgspSPkVBt26NPwsfdFlcsCDamvPUU6kZDaPwB9z+\nYNJ3HQoPCNPF4gQpJAYoTLqCy2H8JAHBGK9Sy1AsconR27w5t4x4PoU+8/L9jXORN5uMxUoakklZ\njFsLLB+SSESRLp27avbz9dhjcHyJq/fGUrBUdYqI7A3091bNVtWs80mq+joQkWA0pd2VceQwjLg8\n8IBLx/7WW01PwWjWDO69F445xhVT/nlFRqYYRt48LSI7A7/GZbdV4O7yipSZUlkJkmbOnPpl30KT\nS5HRXAgO0DJlWEtXTDYT+Qzak8C3UMS5PjLFWKUjHyW4kgi6l+VCoRko0yU6yUYuEyNvveUG/IUm\nsghbxYMWwPfeK934J1yLrJTX3osvxrNwv/hi8rIEiV1oGDgU2N3bZpCIoKr3JSKVYRTAnDnOevXi\ni6n1FJoKrVs7F4MhQ2CvvRpn9kSjaaKqN3mLj3pp1FupagEpGpKnktMKR7kdGaWlMSdgKoRsBWjT\nWS3iZuxLR7o6TcXks8+KMzkRTgwTdO2bMwd23bXwY8Rh0qTc2n/2mXMfPOqowo+dixW+lMSKwRKR\n+4HfAEfiFK1DgUMSlMsw8mLDBjjnHGe1OfDAcktTXrp0cfEW116brCuAYZQSEZkhIteLSB9V3ZiL\nciUi94jIMhGZEVjXQUReEJHZIvK8iLSP2G5XEXlZRN4TkXdF5Kpi9ScJ4tTUMQyj+smUOh+KW5+s\nmHzyibMkF9sdupKIm+TiEGCoql6uqt/xXhX9gDGaJt//PvTrB9/6VrklqQwGDHDukuec09C9xzCq\nmNOALcDDIjJFRK4VkbiJ0O8FTgytuw6YoKr9gZeBH0dstwX4nqrui0uqcYXnNl+RTJhQbgmMQnn+\n+XJLYDQGipnWvpgEE3g0VuIqWDNxiS0Mo2J55BGXXevuu5te3FUmjjsOfvELV4i4Uk3phhEXVV2g\nqr9S1YOBC4D9gVi53VT1NSCc5HwUMMZbHgOcHrHdUlWd7i2vAd4HeubXA8MwDKOxEzcGqxMwS0Te\nBLYnO1TVkYlIZRg5MmsWXH65U7Dapzj4GN/8JsydC2ec4WLT0tWiMYxqQER64zLZngtsBX5YwO66\n+KVBVHWpiHTJcuzdgQOByQUc0zAMw2jExFWwRicphGEUwpdfOsXh17+Ggw8utzSVy803w1e+At/4\nBtx3n1n5jOpERCYDLYCHga+oarETTaetvCMiOwHjgKs9S1Zaxo0bvX15wIBaBgyoLZJ4hmEYRr7M\nmlXHrFl1iR9HNGYVN2/GsK+qThCRNkCNqpYg18r242tcWY2mw7Ztrqju7rvDH/5QbmkCXHONe7/1\n1vLKEWLdOhg2DE47DX72s3JLYzRWvCyziajwItJfVWcXsH1v4ClV3d/7/D5Qq6rLRKQbMFFV94nY\nrjnwNPCcqt6W5Rj6wAP2vDIMw6h0LrggmedV3CyC38TN2t3preoJxCjJZxjJ8vOfOwvW735Xbkmq\ngzZtYPx4uOcel/zCMKqNQpQrD/FePuOBi73li4An02z3V2BWNuXKMAzDMOImubgCGAqsAlDVuUBG\nP3WITokb+n6YiKwUkane6ydxBTcMX1F45BFo0aLc0lQP3bu79O3XXJNfkU7DqFZE5AFgEtBPRBaK\nyCXAzcDxIjIbONb7jIh09+psISJDga8Cx4jINO95dVJ5emEYhmFUOnFjsDaq6ibxgjY8V4k4/g/3\nAn8AMhUkfsWSZRi5Mnu2iyV66ino2rXc0lQf++3n4rDOOgteecWltjeMxo6qXpDmq+Mi2i4BTvWW\nXwdqEhTNMAzDaETEtWD9S0SuB1qLyPHAI8BT2TZKkxI3jIXaGzmxerVLavHLX8Jhh5VbmurlpJPg\nppvg5JMHjpTaAAAgAElEQVStMKlRPYhIGxH5qYjc7X3uKyKnllsuwzAMw/CJq2BdB6wA3gUuA54F\niuXOd7iITBeRZ0RkQJH2aTRSVOHrX4ehQ50FyyiMb34Tzj3XJb2o1IrvhhHiXly5kMO9z4uBX5RP\nHMMwDMNoSCwXQVXdBtztvYrJ20AvVV0nIiNwiTPSOiuNHj16+3JtbS21tbVFFseodH7/e/joI3j1\n1XJL0nj4xS/g44/ha19z8Ww15ghl5EhdXR11dXWlOlwfVT1XRM4H8J4f5glhGIZhVAyx0rSLyEdE\nxFyp6p4xtm2QEjfGcQ5W1c8jvrM07U2cf/3LWVsmT4bevcstTRYqNE17OjZudC6DBx7olFjDKISE\n07RPwiWjeF1VB4lIH2Csqg5O4nj5ICL60vUTyi2GYRiGkYVjf3lcIs+ruEkuDgkstwK+AnSMuW04\nJW79FyJdVXWZtzwYp/ClKFeGsXgxnH8+3H9/FShXVUjLlvDYY8718tZb6/VDw6hAbgD+CewmIv/A\nZbi9uKwSGYZhGEaAuC6C/wmtulVE3gYylir1UuLWAruIyELcg3EHt0u9CzhbRL4NbAbWA+fmJr7R\nFNi0Cb7yFbjySjj++HJL03jp0AGeew6OOAJ69YIzzyy3RIaRiqq+KCJTgSG4yburVfWzONuKyD24\nzIDLAoWGOwAPAb2Bj4FzVPXLiG1PAm7FxS7fo6q3ZDrWsoHHxu6TYRiG0biIpWCJyKDAx2Y4i1bW\nbTOkxPW//yPwxzgyGE2X738fOneG664rtySNn969XX2xk06CHj1gyJByS2QYjtBzCGCJ995LRHqp\n6tQYu4kqHXIdMEFVfyUiPwJ+7K0LHrsZcAfONfFTYIqIPKmqH+TRFcMwDKORE9dF8LeB5S14s3xF\nl8YwQvz97/D88zBlCjSLm/PSKIiDD4a//c2lwn/1Vdhrr3JLZBhAw+dQGAWOybYDVX3NiwsOMgoY\n5i2PAeoIKVjAYGCuqi4AEJEHve1MwTIMwzBSiOsiODxpQQwjzIwZ8N3vwssvQ/v25ZamaXHKKXDD\nDa5G1qRJ0KlTuSUymjoJPoe6+LHAqrpURLpEtOkJLAp8/gSndBmGYRhGCnFdBL+X6XtV/V1xxDEM\nx8qVLgbo9tthv/3KLU3T5FvfcunbR42CCROgdetyS2QYICKtgMuBI3GWq1eBP6vqhiIdoijpaseN\nG719ecCAWgYMqC3Gbg3DMIwCmDWrjlmz6hI/Ti5ZBA8FxnufTwPeBOYmIZTRtNm2Df7rv5wV5fzz\nyy1N0+aXv4SvfhUuvBAeesjcNI2K4D5gNS6WCuAC4H5cdtt8WOZntBWRbsDyiDaLgV6Bz7t669Jy\n9tmj8xTHMAzDSIrwhNdjj92YyHHiKli7AoNUdTWAiIwGnlHVryUildGk+eUv4Ysv4Ne/LrckRrNm\nLh7r+OPhhz+E3/ym3BIZBgNVdUDg80QRmZXD9uHSIeNxad5vAS4CnozYZgqwlxe/tQQ4D7DpH8Mw\nDCOSuPPRXYFNgc+bvHWGUVSefx7+7//gkUdghx3KLY0BrkbWE0/A00/DHXeUWxrDYKqIbM9vKSKH\nAW/F2dArHTIJ6CciC0XkEuBm4HgRmY3LEniz17a7iDwNoKpbgSuBF4D3gAdV9f0i9skwDMNoRMS1\nYN0HvCkij3ufT8dlWzKMovHxx3DRRU656t693NIYQTp2dDWyhg6F3XZzcVmGUSYOBiZ5tRXBue7N\nFpF3cTUW90+3YYbSIcdFtF2Cq5nlf/4n0D9vqQ2jStllF/hPuBqqYRgZiZtF8H9E5DngKG/VJao6\nLdt2UUUdI9rcDowA1gIXq+r0WJIbjYp16+Css1ytq6OOyt7eKD177AFPPukyC3bvDoMth5pRHk4q\ntwCG0ZTo3dsULMPIlVxC1tsAq1T1NuATEdkjxjb3Aiem+1JERgB9VLUvcBnw5xzkMRoJqvCNb8A+\n+8DVV5dbGiMThx4K99wDp58O8+eXWxqjKeLVoloFtAd28V+qusCvU2UYlUSvXtnbFJv+RbS1alHy\najZtrJ5k0yNumvYbcJkE++OUphbA34GhmbZLU9QxyCic+yGqOllE2vvZnOLIZTQOfvUrmDsXXnkF\nRLK3N8rLyJGwaJGzZL3+unMfMYxSISI34ZJSfEh9SvVYhYYNoxy0bVv6Y+63H8yeXdpj9u4NC2yK\no6R06wZLl5ZbCiOKuBasM4CRODc+VPVToBi3jHDxxsXeOqOJ8MwzrtbV449bnaVq4oor4LTT4Iwz\nYNOm7O0No4icg/N8qFXV4d6r4pWrIUOytwnSqlXm7zNZRVq0yO1YjYXDDy+3BNF07lz6YxZrsvKU\nU+K3tdjp9CQ1eWzWxcolroK1SVUVb7ZQRHZMTqT0jB49evurrq6uHCIYRWTmTLjkEpfUYtddyy2N\nkSu33AI772xunQbU1dU1uD8nzExg56QPUkzyGVy1b5/5+6Fp/EeGD8++bRxatoxef/zxhe87X2pr\nM39fqZln8/n9Bw5s+PmMM4ojSyb23Td1Xbt2yR+3KRDnGjgpj+jSYvzXAQ45pDj7yYd8/h9HHll8\nOYpN3CyCD4vIncDOIvJN4L+Bu4tw/MXAboHPGYs3luDBbZSIhQthxAi47TY44ohyS2PkQ7Nm8Pe/\nw2GHwV13waWXllsio1zU1tZSGxj93nhjMoUbPf4XmCYiM4GN/kpVHVnITkXkauAb3se7VfX20Pft\ncK7xvYAa4Leq+rc4+z7vvNLFLBZrwNWhQ6rr0c475+dpcP75MHZs4TJVgwv58OEwcWLDdfnIHXa9\nzmbRDJPPMXdOM21RLVaStm1h9Wro0gWWR5ULT5h27WDVqujvdtgBdt/dZUv2KUZ2xp49Yc6cwvYB\nhf23Dj4Y3n47+rs998x+7+vRAxZnKNs+dKgLRwhSDddkLAuWqv4GGAc8iovD+pmq/iHmMcJFHYOM\nBy4E8OqarLT4q8bP55+7mZrvftc9eI3qpV07l1nwJz9JvQEaRkKMwRUFvhn4beCVNyKyL/B1XKzx\ngcCpIrJnqNkVwHuqeiAwHPitiMSdpKRHj9zlOuig3LfJh9NPT10XNdgeMcJNrJSLmprUddmsVjvt\nFH//ffvmJk8UUTGpIrnJUernYpcu7l0Emkdc0X36xN+X77q6997p2+yYxgcq18mB/fZr+HnYMPfe\npk1u+wnTO1PmgAxkOm6LFqnXatRkxYEHQtccq8xGWR5LSSYr52GHZd8+m3KXzfU0l/9WKcl6qxSR\nGhGZqKovquoPVPVaVX0xzs6jijqKyGUicimAqj4LfCQi84A7gcsL6ItRBaxf72J3Tj4Zvve9cktj\nFIN+/eBvf4OvfAU++aTc0hhNgHWqeruqTlTVf/mvAve5DzBZVTd6RYVfAc4MtVHqY4/bAv9R1S1x\nD9CqlUsQU4mDgaiBXroZ4kJiZfv1y39biB78B4my8jRrFj1g3nPP1PbBgV6+M/o1Nc6qAPVKaosW\nucUyFYNc5D/2WKec9EwTAb/DDnBi2nzQ0aRzMQX3P8j3OjrnHKeEQKobJThXynQlRFq3ziyXT9R+\no8h2PQYJ/h777OOsPmHFSMR9F5yM6dAh+77juMZGxWzm62bn38P891z/K+HJm+B5zBY/evLJ8X7D\nSiCrguU9bLaJSM6OB6p6gar2UNWWqtpLVe9V1TtV9a5AmytVdS9VPUBVp+Z6DKN6WL/eFajdc0+X\nOdBoPJx8Mlx1lXu4rV9fbmmMRs6rIvK/InK4iAzyXwXucyZwlIh0EJE2wMk0dF8HuAMYICKfAu8A\nOUcf7rhjw0HVbuEjhAgrOb6lwefUU0mhWK4z2SxVu+/e0BIRlg2cRSg4u53OBQ2irVNBunfPPtjs\n2DHeOnCDunB2v3bt6mPbfKUsnTXp+OMbWiyGDXNtmzWDo492g/SBA90zL84gOUxw0Jou3i4uffq4\n3ytMUMHs0SPzb57uPIb3l26wHR44R7Xbbbfs10FNTeasjK1apd/HqadG/2d8cq3BGUxe0qtXfAuo\niJts6NgRzj67fr2vsPTt65TeIH361FsSw1ahTp2ijxNsF6V0plOog0T9r8O/XbrjpyNsqQzGf512\nWvZtzzwznuzpGDCgfrnQSZ9MxNW/1wDvisiLeJkEAVT1qkSkMhodvnLVqRPce2953UyMZPjRj2D6\ndLjsMhgzpjriJYyqxHecC+blKyhNu6p+ICK3AC/innfTgK2hZicC01T1GBHpA7woIvur6pqofY4b\nN3r7cvfuDWPUfPbe25U8iEIku7KUZPrvdG5cQUaOdLFVzZq5AWE4zqpt2/jB8+eckzlOK52SEj5H\nHTrAxo2ueD24Cb0PP4wng0hmZej00+GJJ9zyLru4gdpbb7nP4UF92H0t12decPAetnoOG+YUokLj\n2rLdo7MpO337OmvCokXu92/ZEj76KL9jgTtn3bvDiyEfqXBcW76TCEFLSVS8VLdu8WUNy+ErwXHi\nDYM1sVq0cCETrVvXy1dTk6rYiLjJmW7dnDI3eXL9d0EFp3Xr+knOPfds2M7njDNg7dp412TPnqnx\nbDvu6GLdfKU512d9uH1Q+Y5rFaypcfJv25bbsQFmz67jgw/qmDkzeuKhWMRVsB7zXoaRM0Hl6r77\ncjOrG9WDCPz1r+5Bc+utLsbOMIqNqg5PaL/34uo8IiL/Q8MSIgCX4BJsoKofishHwN7AW1H7O/vs\n0YCbdU7nshRF795u4Lf33rB5c259KBfZBuI+/sCqe3dYsgROOAFeeCHeAM13gwoOIMFZAZYFIreP\nP96dt8cfd5/btMm+/zZtnEKWTakMu7XttVe9glVsmjd3Fr+VK1PljxPP16xZw+QA+Ux4ZbNatWjh\nlKKwMhlFIZbV8PWVbl+5ZJHce29YswZmzYrXPnzd5Urw/Ievs3RK/VFHubb//Kf73KZNvauf/x8K\n/6677JLdVb9Vq/hJU3r2hGnT6j+3bevk2rDB/e/efLNehkMOccku0v0+gwe79jU1LqTgkUcaft+h\nQ/p7SdSEU7ZrukMH+OKLhuuOOAJ6964Fahk71mUm/d3vkknKlFF/FZFeAKo6JuqViERGo8JXrjp3\nNuWqKdCmjZvh/dWv4Nlnyy2N0VgRkVNE5Ici8jP/VYR9dvbee+FqPz4QarIAOM5r0xXoB2TNDejH\ni8SleXM3o929uxtMnXNObtu3aJF5MOsPYLK5J2YbvASPEbfYuD+o8wfCcQd5I0dC//71y3371lsb\nwq5nNTUNB2lRVr6BA51y5PehXz/ndpQumD7KrUyk/hx16hTPhS5f0il+2ZJhBONukvAoiNpnXGW7\nEILX3p6BVDS5xOZ07w4HHNBwXaZzFBy7+BMmSWey23XXeuUrLJvvzhhcX0j9u3BfBqVxuh4xwp2L\nnXZKTW6x117xfoNWraLHgkGX2/D/KepcB/sexxJ15JGlLQmUzUD4hL8gIo8mLIvRyFi/3j0MO3d2\nLmOmXDUNeveGxx6Diy6C114rtzRGY0NE/gycC3wHl6H2K0Ceeb8a8KiX+v1J4HJVXRVMygT8AjhC\nRGbgXAl/qKqf53qQXAdluQxYzz8/epAVLHKca60oX5HJRBw3wLPOyi+TIjgFw3dnatbMHW94AXbM\n/fZzsRxBBS/TwDDbwPWgg5It7pzLb9a5c7I1s9q2zXyuslkBi11Xy5clrCwVEoaQSdHq0qVeEcj0\nm3fu7K6LpOqIRd0X8nUZ3jOcL5XUa+7ww11/g8ft3Nndc/zz5bsxpmP33Z2bbVQSkcGDG8YaxkkG\ndPDBDeU77zw3URIkqLRFxfglGcqQ7RIMHjriJzCMaHzlqksXU66aIocfDv/4h7vZvfNOuaUxGhlH\nqOqFwBeqeiNwOM6aVBCqerSqDlTVg1S1zlu3PSmTqi5R1RNVdX/vlVcETCliE8NK3B571C/vumvD\ntOzpEih07uxm+bMpMm3bxssI5w/YDj00ntJWCg4/PHVd0A2rlPipubOlpA4TZSHNtWZWOqImAzp1\nSh3EBtl334aJG/xkDf51X4xsilFyBRMXAJx7brx9nXlmvXU0E/41EUxZf9hhsP/+0e379nVuiKXO\nHpmOTMpg27bZJ352373h7xrGV2j79XOWosMOa2iF6tzZKTetWzdUfv17U58+mSeTTjoptRBzWDEU\nyc2KedxxuafEz4VsCpamWTaMtJhyZYCLr7jjDudSMHduuaUxGhF+JMQ6EekBbAZyHJYmjz94DytU\nwQFArtasKOUsUzatQw+tH+AFXWNat67fVzplYscdXXxCNk49Nbd7/F57xW9frHTM/ftHxwlFyVFT\n45TOpNy/dtst2tLgu2wGz3kcGfbZJ17dp3wU+3zOv0jDwbyfrCHcl27d6t3fjjwyt5ThwX3F6Vcm\nF9aWLeNdj/vv76w1u+9ef/wddoj+LUeOTC0NUKiFM04/g+clrmJ3+OENlcbwuYhrCQwqtLvt5pQf\nP7X/iSemv5dETbYcd1xD6xS4ayWfbJyZ/kOdOyc74ZXtsjpARFbhLFmtvWW8z6qqCRk/jWolqFzd\nd19p/LGNyuWcc1zA/vDhLlA3bn0Rw8jA0yKyM/BrYCpu8u/u8oqUyu67wxtvFGdffmayqMHA0Uen\nZi0bNMj973IpEBumFJa28DFOOcVZXx591FlCdtzRBcVn44ADYOnS9N937OhePXrkl3WsWAwe7AZ1\nr77acP3QoS7zYSnp0SM1AUCQUaOyKwWFXCN+CvFt23J358tF+a2pyV1R9Ps1fLhLErJpU/q2Xbs2\nzAoIqW6Sp5ziFLHPc3Yojk/Xru7lHyOua2I4dql5c9gSqO5XU5N7HGiYdPGJe+8drWAFM2hC9Wad\nzii2qtaoajtVbauqzb1l/7MpV0YD1q1zNQy6djXlyqjnG99wSS+OPTY6Zaxh5IKq3qSqK1X1UVzs\n1d6qWnCSi1LhD/bCtWB80g1ajzsue2IKn86dU5Ur30Uv06DYH1Dvv3/2+3cxLTy+TO3a5R4jBs49\n7JhAkv50snXsGF2zJ53L2UEHpa6P63oWRZ8+ro/h4/XqFS1Dr17Fc2EKuw0OG5a5fZs2xY0rSxeb\nlc/guXPn+LWXRJwyd8IJDdfHuX7juLK2bOksxZlo167wCYts2x9zTOb4p2z4NerCBaVFkhvLHXRQ\nPHfWXI8fbO/vP1t9rSRIXC8UkZNE5AMRmSMiP4r4fpiIrBSRqd7rJ0nLZBSftWudq0iPHs4t0JQr\nI8gFF7gU7qedllrjxDDiICKHiki3wOcLgYeBm0QkwRxuxcWPYTn55Ibrs6VyL9SdJexyE2bUqHol\nJd+BWjC2KxPVMCO9885uhh0auk01a5Y9e182wu5j6Rg4sKHiWMgAOug9kEk5KZblslOnhsrygAEu\n0Uk24hy/bVuXjh/iKUqtWmV2Ewy7xWVzGcx3cqEUVuHwfytuHbouXdx13aaNUxh9l+KkMyUWypAh\nmRXc4cPdfSlO0oxik2h0jIg0A+4AjgU+BaaIyJOq+kGo6SuqOjJJWYzkWLPGKVd77AF/+YspV0Y0\np5ziXH/OPhtGj4Zvf7vcEhlVxp3Up0k/GrgZl0nwQOAuIEMIdmUR5bLkzyDnSzaXoGyDxjZt0rtC\nZZuh92nd2tWZWbs2c3Kbnj2dRcGPHQpbrVq1iq9QRiXYKLTwaZjmzQtXqoLst5+rIzRvXvxzC86y\nGJVU4Ygj3P4yZW1t1sxZJ9q0KU4SjAMOyKwodulSr1B16ZKqcMUlmKAlH9L9tkHFYe+9nbXQb7vD\nDi4jHbhzG3SZKxdR/RgyJFVZPuSQhgk/dtvN/R9zIezyWMlEXR/Bc1WsOM58SDr9wGBgrqouABCR\nB4FRQFjBKoFebyTB6tVu4Ny3L9x9d3XMTBrl46ijnE/7qFFuAHb77fk9dI0mSU0gLfq5wF2em+Cj\nIjK90J2LyNXAN7yPd6vq7RFtaoHfAy2AFUkVPc6Vbt3iu0z5DBoUfyCVy4Crd2/nLp5JwRKptyhE\nWb2CacYzueSdeWa0G1uSKdOLTTEGs+GYlXQUs1ZXOGtfJvxMgrmw226waFF2BStfC0vwudOsmbNw\n+PF5wQF6rlkdS0nUuWnZsqFS0apVw0yTuT5vq3lMFye7aZIkfep6AosCnz/x1oU5XESmi8gzIpLD\n39YoJ6tWuQxxe+9typURn732csH/S5a4uJIlS8otkVEl1IiIPyl4LPBy4LuCJgtFZF/g68AhOIvY\nqSKyZ6hNe+CPwKmqOhBXfyvGvtN/Vyz3m+HDo7PkZaJ169T4nmINSOJktQseM9NxMz1XWrbM/H24\nNlIl0aNHMsWJo85HMZ/NSbu5hS2xSRzv/PMLm9jr1i27S2+xKYZn0E475ZaJb9So6pwA9a+ZdCUo\nSkUlJNB+G+ilqutEZASuuHFkTZPRo0dvX66traU2Tg5ZIxE++cRZro48Ev7wB1OujNxo1w4efxx+\n8QsXG3L//fnNchqVQV1dHXV1dUkfZizwLxH5DJeq/VUAEdkL+LLAfe8DTFbVjd4+XwHOBH4TaHMB\n8KiqLgZQ1c8KPGbO+AOHQrIDZqJly/iucJUemwGpGdLSUY6+9OyZOcV+Ppx8cqpL1AknpC9Am0+/\nc7WU5kom62NUHNWAAcWrqxb3fNTUJPcfTEcx3DpPPjk3hTWXiZJKokUL59ZZ7hJBSR9+MRCssrGr\nt247qromsPyciPxJRDoGXEG2E1SwjPIxfbpLVvCd78APflCawE2j8dGsGfzsZ05J/9rX4NJL4ac/\ntRi+aiQ84XXjjTcW/Riq+j8i8hKu5tULqtuHQ81wsViFMBP4hYh0ADYCJwNTQm36AS1EZCKwE3C7\nqt6faadduuR+PWe6n/qZ2Eo9e54vcVNFJ8HZZ5fOVbAYcR7FUPKiMlNmSu6QD6V63sc9TsuWmd34\nKml8km+c5fHHF8faWW3P1k6d4LM8prFOOKEyJoCSVrCmAHuJSG9gCXAe0GB+TES6quoyb3kwIFHK\nlVEZPPMMXHKJKyBbaG0EwwCXJWvqVPjqV11NnzFjqivI1igdqvrviHVzirDfD0TkFuBFYA0wDdga\natYcGAQcA+wIvCEib6jqvKh9+hOCr75aPI+Lrl3rg+8rnWImhciHUilX55xTnIFrOYPxK42DD3Zx\nZQsXlluS4rLnnu6VK0lbDSuVfD2jsrk6l8jjIlkFS1W3isiVwAu4WcZ7VPV9EbnMfa13AWeLyLeB\nzTi3jwKqTBhJsXEj/OQn8OCD8MQTLrOOYRSLbt1c+vY77nCV5f0sg+Z6apQKVb0XuBdARP6HhvHD\n4GKIP1PVDcAGz43wACCjglVsCp2RL1bcT7q6Rk2NYlkFevd2Fk8D+oWCRKopaYlRPJL63UvhcQEl\niMFS1X8C/UPr7gws/xEXOGxUKG+84YrF9u0L06Y13dkUI1maNYOrrnKphC+6yCnzt98eXezTMIqN\niHRW1RUi0gs4AxgSavIk8AcRqQFaAocBvyvW8f14hyStt8WyKp11VvnjGxobIqWPecnVjWrYsOIV\nPs6FQkoY5DIhUVPjsmsalcGQIW5yv1qxW6SRlg8/dFarV16B3/3OuUJUkj+z0Tjp39+lcr/nHpel\nctQolwwjbipiw8iTR72CxZuBy1V1VdDbwnMjfB6YgXMfvEtVZxVywOAAt3Xr8rvVxaUaM4sZqXTu\n7EqtxKVHj+RkiUOu449+/XJPDtG/f/Y2RmnYYYfqvteYA47RAFWYNMnFwxx2GOyzD8yZ42qRmHJl\nlIqaGpf04v333cBz771dQpVPPy23ZEZjRVWPVtWBqnqQqtZ56+70XNn9Nr9R1X1VdX9V/UPZhDWM\nInDkkW4Sqxro3Dn3rIsHHwz77puMPIaRDVOwDACWLYNbb4WBA+Hii51b1vz5Lsub+dob5aJDB3dd\nTp8Omze76/Ob34QpUyojS5BhFINyF8Q0miYilT9x6st33HG5FTc2jHJjClYTZuFCV8Nq2DBnFn/7\nbfjTn2D2bLj22vKm2DWMILvt5hStOXNcbZvzz3eFVX/zG1iwoNzSGUZhnH56uSUwDMMwionFYDUR\nVOHjj12SiokT4YUX4PPPXeG5a691dRaKUcjOMJKkUyf4f/8Prr8eXnvNpXQ/5BDnOjJypHsNGmTZ\nBw3DiE+1FlRNggMOgE2byi2FYVQ/pmA1Qtavh5kz4Z133Gv6dJgxw1VzP+AAV2vowQfdsg1EjWpE\nBI46yr3uvBP+/W948klXsHjVKjj2WBg+HGprYY89Kt8NxjCM8lDKgsTVQN++5ZagnmbNYKedyi2F\nYeSHaJUEMoiIVouspUIVli5tqEi98w589JHLnnPggU6J8l+WXr2EXHONe7/11vLK0QSZO9dZaSdO\nhLo6l4Vo2DAYPBgOPdT9F8xamywigqo2WbU27vNqxQqYMKF6sgcahmE0NpJ6XiVuvxCRk0TkAxGZ\nIyI/StPmdhGZKyLTReTApGWqBHKtIr15s7NK/f3vLpvaCSe4ehQDB8Kvf+2SVJx0EowdC19+6RSt\nMWPge99zs/lJK1elqIpdSqw/lU2m/vTt6zIQjh3rsg6+8ILLljVjBnzrW67Q6qBBLlnGrbfC88/D\nokXlTZrR2H6fakRErhaRd73XVRnaHSoim0XkzEKP2blz9WRxKyZ2vWfHzlF27Bxlx85R+UjURVBE\nmgF3AMcCnwJTRORJVf0g0GYE0EdV+4rIYcCfSS3w2Oioq6trUEnaZ+1aNwM/e7Z7zZkDs2bBBx+4\nQH/fGnX11e69Z8/KcH9K159qxfpT2cTtj4hL4NK/v1O6wLnQzpgBU6e6/9bTT7v31audcta7d+qr\nVy/YZZfkXGob2+9TbYjIvsDXgUOALcBzIvK0qs4PtWsG3Aw8X6xjF1JEtVqx6z07do6yY+coO3aO\nykfSMViDgbmqugBARB4ERgEfBNqMAu4DUNXJItJeRLqq6rKEZUsMVRckum6de61e7axKwdekSfDT\nnz3/DqoAACAASURBVLoZ9sWL69/XrYO99nIufv37u9SkV17pMqZZunTDKJzWrV2Nt8MOa7j+iy/c\n5MbChS4z4UcfORdD//Pq1U7J6tLFWR46d06/3KkTtG9vsR1VxD7AZFXdCCAirwBnAr8JtfsOMA44\ntLTiGYZhGNVE0gpWT2BR4PMnOKUrU5vF3rqiKVjz5rn046r1bkD+ci4vcK56vuIUfK1d2/Bz8+Yu\nM1GbNi5Is337hq/Vq93ga8gQZ4Xq0cO9d+pUGRYpw2hqdOjg4rQGh+9QHps2wWefubiZFStg+fL6\n5alTG67//HNYudLFeu28s3t16OD++y1b1leo919vvw1r1rj/frNmybyLuPvY0KHO+m00YCbwCxHp\nAGwETgamBBuISA/gdFUdLiJprhLDMAzDSDjJhYicBZyoqpd6n78GDFbVqwJtngL+V1UneZ8nAD9U\n1amhfVmGC8MwjCqh2pJciMglwBXAGuA9YKOqfi/w/cPAb1T1TRG5F3haVR9Nsy97XhmGYVQJSTyv\nkrZgLQZ6BT7v6q0Lt9ktS5uqe1gbhmEY1YOq3gvcCyAi/0NDzwpw8VkPiogAnYARIrJZVcdH7Mue\nV4ZhGE2YpLMITgH2EpHeIrIDcB4QfhiNBy4EEJEhwMpqjr8yDMMwqg8R6ey99wLOAB4Ifq+qe3qv\nPXBxWJdHKVeGYRiGkagFS1W3isiVwAs4Ze4eVX1fRC5zX+tdqvqsiJwsIvOAtcAlScpkGIZhGBE8\nKiIdgc045WlV8FkVamsugIZhGEZaqqbQsGEYhmEYhmEYRqWTeKHhXClHscckidMfEakVkWkiMlNE\nJpZaxlzI1h8RaSci472i0e+KyMVlEDMtInKPiCwTkRmBdR1E5AURmS0iz4tI+zTbZi2aXWry7Y+I\n7CoiL4vIe9n+a6WkkN/Ha9tMRKaKSEW4bhV4vbUXkUdE5H3vdzosql0pKbA/3/XucTNE5B+e23ij\nohLvEaUg3f0k07UhIj8Wkbne9X1CYP0g7xqZIyK3lqM/SRK+R9k5akjUfc/OUUOi7qVN/Rzl+mzK\n9Zx45/hBb5s3xLmSZ0ZVK+YF7AvMAFoCNTjXwj0j2jUDXgKeBs4st9yF9Adoj8tY1dP73KncchfY\nnx/jskKCCwT/D9C83LIH5DsSOBCYEVh3Cy5zJcCPgJvTXHPzgN5AC2A6sHcV96cbcKC3vBMwu5r7\nE2j7XeDvwPhy96XQ/gB/Ay7xlpsD7aq1P0APYD6wg/f5IeDCcvenyOemIu8RJep75P0k3bUBDACm\nedf17t558z1qJgOHesvP4jIRl72PRTxXDe5Rdo5Szk/4vtfezlGD8xN1L72oqZ+jXJ5N+ZwT4NvA\nn7zlc4EHs8lUaRas7cUeVXUr4Bd7DOMXe1xeSuHyIE5/LgAeVdXFAKr6WYllzIU4/VGgrbfcFviP\nqm4poYwZUdXXgC9Cq0cBY7zlMcDpEZtuL5qtqpsBv2h2Wcm3P6q6VFWne8trgPdx9efKSgG/DyKy\nK65+0V8SEzBH8u2PiLQDjlKX2Q5V3aKqq5KUNQ6F/D64SZkdRaQ50Ab4NBEhy0dF3iNKQZr7ya6k\nvzZG4gYoW1T1Y2AuMFhEugFtVdWvQXYf6a+nqiPNPcrOkUea+96X2DkKE7yXtsZl3m7S5yjHZ1M+\n5yS4r3HAsdlkqjQFayZwlGfWa4O7EQVTuAeLPf4fUOmpcLP2B+gHdBSRiSIyRUT+q+RSxidOf+4A\nBojIp8A7wNUlljEfuqiXuVJVlwJdItpEFc0uu0KShjj92Y6I7I6b+ZmcuGT5Ebc/vwd+QOUnIIjT\nnz2Az0TkXs+d6C4RaV1SKeOTtT+q+inwW2AhbjCwUlUnlFTK5Kmme0RiBO4n/wa6prk2wudqsbeu\nJ+68+TS2cxh1j7JzVE/Ufa8Ndo62E3Ev/dK7l9o5SiXdsymfc7J9G8/AsFJcUqS0VJSCpaof4Ex6\nL+JMc9OAraFmt+JMfT4Vq2TF7E9zYBAwAjgJ+KmI7FVKOeMSsz8nAtNUtQdwEPBHEdmppIIWTqUP\n0HMlbX+832YccLU381wNpPRHRE4Blnmz6EIF3xciiPp9/PvCH1V1ELAOuK6kUuVP1O+zM24GsDfO\nxWUnEbmg1IIZyRJxPwlfC43t3hqbiHtUOprsOSL1vrcWd9+z68gj4l66o4h8FTtHcSjmOck6xqgo\nBQtcsUdVPURVa4GVwJxQE7/Y40fA2bgB/MgSixmbGP35BHheVTeo6n9wbncHlFjM2MTozyXAY17b\nD4GPcL74lcwyEekK4JmIo1xP4xTNrhTi9AfPvWAccL+qPllC+XIlTn+GAiNFZD4wFhguIveVUMZc\niNOfT4BFqvqW93kcbuBRicTpz3HAfFX93Jv9eww4ooQyloJqukcUnTT3k3TXxmIaej/45yrd+sZA\n+B51jIjcDyy1c7Sd8H3vUdx9z66jesL30sdx91I7R6kU85xs/05EanAx0Z9nOnjFKVjSyIo9ZusP\n8CRwpIjUeKbww3D+6xVJjP4swN0A8C7sfriAzEoibOEYD1zsLV+E+03CxCmaXS7y6Q/AX4FZqnpb\ncqLlRc79UdXrVbWXqu6J+21eVtULkxY0Jvn0ZxmwSET6eauOBWYlKGMu5HO9LQSGiEgrERFcfyr2\nPpcnlXyPKAVR95N018Z44DwvM9cewF7Am54bz5ciMti7Ti4k/f2rqkhzj/ov4CnsHAFp73vvYddR\nkKh76SzsHEH8Z1M+52S8tw+ArwAvZ5VGKyD7R/CFs+DMxLmf1XrrLgMujWj7Vyo4i2Dc/gDX4m4i\nM4DvlFvmQvoDdAee9/oyAzi/3DKH5H8AF1y/EXejugToAEzAZb56Adg50JenA9ue5LWZC1xX7r4U\n0h/cbOpWXKazacBU4KRq7U9oH8OonCyChVxvB+AG7dNxFp/2Vd6fG3BK1QxcsHCLcvcngfNTcfeI\nEvU78n4CdIy6NrxtfozL3vU+cEJg/cHAu945vK3cfUvofG2/R9k5Sjk3Kfc9O0cp5yjlXtrUz1Eu\nz6Z8zgkue/bD3vp/A7tnk8kKDRuGYRiGYRiGYRSJinMRNAzDMAzDMAzDqFZMwTIMwzAMwzAMwygS\npmAZhmEYhmEYhmEUCVOwDMMwDMMwDMMwioQpWIZhGIZhGIZhGEXCFCzDMAzDMAzDMIwiYQqWYRiG\nYRiGYRhGkTAFyzAMwzAMwzAMo0iYgmUYhmEYhmEYhlEkTMEyDMMwDMMwDMMoEqZgGYZhGIZhGIZh\nFAlTsAzDMAzDMAzDMIqEKViGUWRE5Mcicle55TAMwzCMdNizyjCSQ1S13DIYhmEYhmEYhmE0CsyC\nZRiGYRiGYRiGUSRMwTKMAhCRH4nIJyKySkTeF5HhInKDiNwfaHOhiHwsIitE5Cci8pGIHON9d4OI\nPCwi93v7eEdE+orIdSKyTEQWiMhxgX1dLCKzvLbzROTScvTbMAzDqB7sWWUYpcUULMPIExHpB1wB\nHKyq7YATgY+9r9VrMwD4I3A+0B1oD/QI7epUYAywMzAdeB4Qr91NQNBHfhlwsne8S4Dfi8iBxe6b\nYRiG0TiwZ5VhlB5TsAwjf7YCOwADRaS5qi5U1Y9Cbc4CxqvqG6q6BfhZxH5eVdUJqroNeAToBNys\nqluBB4HeItIOQFWfU9WPveVXgReAo5LonGEYhtEosGeVYZQYU7AMI09U9UPgGmA0sFxEHhCR7qFm\nPYBFgW3WA/8JtVkWWF4PfKb12WfW42YIdwIQkREi8oaI/EdEvgBG4B5yhmEYhpGCPasMo/SYgmUY\nBaCqD6rqUUAvb9UtoSZLgF39DyLSGtgln2OJyA7AOOBXQGdV7QA8h3uoGYZhGEYk9qwyjNJiCpZh\n5ImI9PMChXcANuFm8LaGmo0DThORISLSAjeDmC87eK/PVHWbiIwATihgf4ZhGEYjx55VhlF6TMEy\njPxpCdwMrAA+BToDPw42UNVZwHeAh7w2q4DlwMYcjqPevtYAVwGPiMjnwHnAk4V1wTAMw2jk2LPK\nMEpMSQoNi0gz4C3gE1UdGfH97Tj/3LXAxao6PXGhDKMMiMiOwEpgL1VdUG55DMMwDCOMPasMozBK\nZcG6GpgV9YVnOu6jqn2By4A/l0gmwygJInKqiLT2Hli/BWbYA8swDMOoJOxZZRjFI3EFS0R2BU4G\n/pKmySjgPgBVnQy0F5GuSctlGCVkFM7l4hOgD/+/vTMPk6K8HvV72AZEEBcEFAFFiaIiIIuCy7iL\nGy64Jnr1lxizmGhMjInGiEZvjDe5cYvXJWqiRlxQARVFRUYUFZFdkE1BWQIoCgoIDMy5f3xdTE1P\nVXf1Wt0z532eerqqazv9dS3f+c7m3CUMwzAMo5Swd5Vh5IlmRTjH34FrcUXrgtgTX2pQYHniu1XB\nmxtGeaGqlwOXxy2HYRiGYYRh7yrDyB8FVbBE5FRglarOEJFKckjRKSKFDxYzDMMw8oKqNtqUzPa+\nMgzDKB8K8b4qtIvgYOAMEfkUGAEcIyKPJW2zHNjLt9w58V09VLWkpptuuil2GUp5svYJma66yk3W\nRpEma6PyayOj9N5XpTSV2vVaipO1kbWRtVFxpkJRUAVLVa9X1S6qug/Ol/dNVb0kabMxwCUAInIY\nsFZVzT3QMAzDMAzDKBobNsC25AphhpEFsdTBEpErROTHAKo6FlgsIouAB4CfxSGTYRiGYRiG0XgZ\nMwZmzoxbCqMhUIwkFwCo6lvAW4n5B5LWXVksOfJJZWVl3CKUNNY+6bE2So+1UXqsjYxyopSv14kT\nYeBAqKiIV45SbqNSoVBtVF1dkMPGgl1H8VGUQsP5QES0XGQ1jJRcfbX7vPPOeOUwjAIhImgjT3Jh\n76vyZMQIOPZY6JCHYjHr18P778Pxx+d+LKM4jBgB++zjlGyjcVCo91UsLoKGYRiGUSqISGcReVNE\n5ojIbBH5Zch2lSIyXUQ+EpEJxZbTKC++/BK++CJuKQzDiIOiuQgaRmNk0ybn0/3aa86ve/Nm+ONX\nsGdn6PCpGykzDCN2tgLXqCspsiMwVUReU9V53gYishPwD+BEVV0uIrvFJaxRl88/d25d3bun3q66\nGpo1A2m0tlXDMIpFQS1YIlIhIpMTI36zReSmgG2OFpG1IjItMf2hkDIZRjFYvx5uuw06d4aHHoLe\nveHuu+GJJ6BfP6jeAv37wxVXwLffxi2tYTRuVHWlqs5IzK8HPsYVvPdzEfCcqi5PbPdlcaU0wpg8\nGT74wM3PmgXjxwdvN3IkzJ9fPLlMkSs9Nm2C6dPjlsJoDBQ6Tftm4BhV7QP0BoaIyICATSeqat/E\ndGshZTKMQrJ5s1Ok9tsP5syBd9+F11+HK6+Eww+HXr2gWzc46ihYtAhqaqBvX9cpMAwjfkSkG+59\nNTlpVQ9gFxGZICJTROTiYstWqqxbB6UScrZsGaxeHb5+w4biyWKUHitWwLx56bfLhq1bYcuWwhzb\nKD8K7iKoqhsTsxWJ8wU9hm2cxyhrtm6Fxx+H4cPh4IPh1VfhkENS77Pzzs669cQTcMIJMGqUU8IM\nw4iHhHvgSOCqhCXLTzOgL3As0Bp4T0TeU9VFQccaPnz49vnKysoGnc1r7Fj37OrWLW5JSkfRKwW+\n/RbatIlbivyxdKkbjDz11LglCeatt1zc3fnnxy1JabB+Pey4Y9xS1KeqqoqqqqqCn6fgCpaINAGm\nAt2Bf6jqlIDNDheRGcBy4FpVnVtouQwjH2zd6hSkW2+FvfaCJ5+EwYMzO8YPfgC77gpDh8JLL8GA\nIBuvYRgFRUSa4ZSrx1V1dMAmy4AvVXUTsElEJgKHAGkVrMbA1q1xS5A/GoJrX02Ne59ceGHckuSP\nlSvhm2/iliKcDRtcu6u6a2jtWmjXLm6p4uPFF+G886Bp07glqUvygNfNN99ckPMUw4JVA/QRkbbA\nKBHpmaRATQW6qOpGERkCjMK5YtSjMY0IGqXNN9/Av/7l3AG7dIGHH4ajj87+eEOGuGMMHQpvvw37\n7ps3UQ2j4BRrRLDAPALMVdW7QtaPBu4RkaY4j4yBwP8NO1h1NTRvnn8hjeKxbRtMnWqDXkZmvPCC\ns+bOnw+nnAI77RS3RPHRmC3KxSw0/E0ire3JwFzf9+t986+IyH0isouqfpV8jMY2ImiUHvPnwz33\nOEvViSfCY4/BoEH5Ofbpp8PNN8PJJ8OkSfmpw2IYxaBYI4KFQkQGA98HZovIdJwr+/VAV0BV9UFV\nnSci44BZwDbgwVTeFhMnwnHHFUH4IrJ6NaxZAwccELck2ZNJh2/9evjkE1OwGjNffw2zZ7u46XR4\nls/Nm919As6ilQnz58O0ac617rTTGoY1tbFSUAUrkca2WlXXiUgr4ATg9qRtOqjqqsT8AFzx43rK\nlWHEhapLs37XXW4088c/dg/cPZNzjOWBH/8Yli93PuZVVaXpv2wYDQ1VnQSkdWRR1b8Cf41yzE2b\ncpUqXmbPdh3Fnj1hhx1qv1u9OljBaswj1XFTXe3av0WLuCVpeCxf7qZi4Slm69c7C2ozK6aUMZMm\nQatWLoFYnBT6r+sE/DsRh9UEeFpVx4rIFSRGBYFhIvJToBr4DrDwQKNkGD8errvOZQb61a/g+eeh\nZcvCnnP4cJcJ6/zzYfRoe8AaRjlS7grHggXuubdwYenH8ZR7W+fKG284hf6ss4pzvpoaaFLQHNSF\no1jXSmO/Jj3iaIfPP4eKigauYKnqbFzWpeTvH/DN/wNXvNEwSoY1a+Cqq1ya9dtvh2HDivdCEYH7\n73dWrF/8Au67z9wEDMMwGhrbtrnkIBUVuR1n/frwJCNewoV88vTT0KcP7L9/fo9r1EXVDbZ27hy3\nJEY2lOkYhGEUjilT3MujfXvnEnPeecUfrWve3BXFnDQJ/hrJIckwDAARaSUi34tbjm+/hXHj3PzC\nhfBVFo7vNTWuE26UBqr5G5H/7jt45hnnFVGOZJLNb+xY93tzxa8obtxYuJpmURXSLVsKe39u2uSS\nXnlUVxfuXA2NUrAgmoJlGD7GjHFZf+68E/7+d2jdOj5Z2raFl192mQr/85/45DCMckFETgdmAK8m\nlnuLyJi45PGUqg8/dIM1qaiuro2/8Jg40aU6NkqDd991ykI+WLgwP8cpB9atc1MmpOsgjx3r3o+F\nJpWy9dxz8MEHwevy0cH3H+Prr92ga0Nh8+a4JSg8pmAZRoLnnoPLL4dXXoGzz45bGsdeezl5fvMb\neOqpuKUxjJJnODAAWAugqjOAveMUKCqzZrlkOn7Wrs3PyH+ufPZZ+pH6UhgxjkIuyUfWrClMHabF\ni2HevPwfF8rnf/GzalX69111dWlYdzduLM55SlEh2bw5PEtiqgGlbdvisdyuWJGdJ0G2mIJlGLiM\nfT/7mXPp6dcvbmnqctBBruP1q18533fDMEKpVtXksfKS62KuXOnSf/vJNJ1zMXn3Xfjvf/N3PE9J\nGT8e3nknf8eF9Arp0qXh69L9B5nEMq1f7xSFKEyf7qaobNqUX8Wpqiq64rlli/sspOK2fn36bRoa\nNTXOZTTqtqWgcD3/fPB1qwoffRS+X1zPurfecs+yYlFQBUtEKkRksohMF5HZInJTyHZ3i8hCEZkh\nIr0LKZNhJDNvnsvYN2IE9C7Rq+/gg53yd/XVTk7DMAKZIyIXAU1FZD8RuQco4is1Gh9+GO5aFIV0\nsRhr16buAM+fD59+mv35/Xz5Zeb7vPyyk3H1aqdsekyf7jqZuVgFcolTyWen9YMP6ivR+eKFF5zF\nM5W83v//xRfh1oTp012c73//m3pkf+7cWmvRjBnZyeyxdGnm7oJxkquVLKpSvm5d/XOF3cMzZpRO\n7F7QvfrGG3WXt2yB99+vXS5G0q4oAwBvv11Yi1ZBFSxV3Qwco6p9gN7AkEStq+2IyBCgu6ruB1wB\n3F9ImQzDz4YNcM45cNttcOyxcUuTml69nCXruuvg1lvL0/XDMArML4ADgc3ACOAb4OpYJcqBoI7I\n6tXpYzFeeQWWLAlfP22amzI9dxATJ9Zd3rzZjV6PHp16P68zqQqTJ7vO/vz57vv586OdO9941hmP\nTJ+x/u0L3YmcO9f9z+mYN6/WmrB+fV3l5rPPXErrdMycWdsRzdX68M47boAhmS1bnDIYhSAFLVk5\nWbvWZVZcssR1+N99N1p7JZOJS2ku/3lUV+Dp0wvjppoJc+fWDtAE3SPJgy5r1jg32GQ2bgzPflkM\nli0rbI2zgrsIqqqn31bg0sIn/x1DgccS204GdhKRDoWWyzDAuQX27w8//GHckkTj4INdZ2T0aLjk\nktJwEzCMUkFVN6rqDaraX1X7JeZLouTvihW180EdMX9HZfLk8I5v1A5friPvyc+W5AQcyXjyz5zp\nLCYbN6a2JnnxZlu3us5aLtaefMXipEtgkaoDvXKlixv68kuXkCCdsuB3ucx2sCxKp9xvIXzjDfcb\nPVe0Ypf/SFZg/cyeXd/ykYwn79df11/3zDN1rbKvvOKsfMuWuf9i9WqndIWxenX09+m2bS5LaBzM\nm5feXTeqhfDbb13W5Nmz67plzpuXWombOdO1LeSmoIweXdeSr1r3Gpk5M7eYSVX33/utVIXKPBlE\nwRUsEWkiItOBlcDrqjolaZM9Ab9X9PLEd4ZRUEaPhvfeg3/8o7zqTHXq5HyJN22C445zLwbDMEBE\nJojIm8lT3HKlIqhznaxwjBpVmID+TDpHc+dmrgh4nbBMyYfyAa5DnYtLddh7YeZM9zliRK1Fx+ug\nvv66UxSC/q9XX3UFnCGaa9LEie4dlQuefCNG1CpkQbJFaecRI3L7P156KfxcYccNUmbC2iQ5diuT\ne2b8+Fr3x2++qbXCeHKNGwcTJrj37qRJtb8lnVJWXR3NSuid57PPosscxtix0axcn38OixY5C+eL\nL9Zak6ZPd268X3xRe72mI5VlM1X/yt9+c+a4ZGMec+c6d9hsef99dy/6FfeamrqDDoWkGBasmoSL\nYGdgoIj0zPZYw4cP3z5VVVXlTUaj8bFuHVx5JTz0ULyp2LNlhx1cwotjj4UBA3L3izeMXKiqqqrz\nfI6R3wDXJqYbcSnbAxyS6iIinRPK2JxEvPAvU2zbX0SqRSQvuUafesq98L0Olpf0QbU2vuG77zKP\nLUo3aFRdXd+9L9/4R6O3bUttwfDL643Az5+fPptcKqtTKotFFMI6/XPnpt4mzEL59dd1LZnJeJ12\nj+XLo3XOk9m0yVlBM2HixLrXXhipkoSk4ptvsvO4mDOnfrtEJeh/SBX76FnAvvjCKcp+vvrK3adv\nvlk37shvLQqyan3+uVPI0rnMBpEuFtFT9IOIolwmt0/yM2bWLJg6NXjfZOtpWAKuKPeHZ6XyrEtb\nttR1HVy4sL5b6YgR6S11ntU9uS38xy5kqEWzKBuJyMGqmqaKR2pU9RsRmQCcDPgeTywH9vItd058\nV4+YX9xGA+J3v4MhQ+Doo+OWJHuaNIFbboEDD4QTToD773fxZIZRbCorK6msrNy+fPPNN8cih6om\ndwcmiUiUdBJbgWtUdYaI7AhMFZHXVLVO8mwRaQLcDozLj8QO/8i7Z71Jl4Hu669h550zO8/YsXDI\nIZntk47kTrO/87RkiRsAatrUxcAsWwYXXhj92P44se++c8uDB9fdxt/JKnRcapCSG/Wc/t+SrPh5\no/9BHfSamsz/65Urs7OELF3qOtSffw5nnRW8jT92zpuaRBiqD+uopyNfxXWzLXeQ/P+mc30NY+NG\nd2936hR9n3R2hLlzU9/Pq1Y5D5eDDw5en4nnzqZN0LJl8LpU98DixbDTTuHrV650Vqphw2q/GzUK\ndt21dtlTrpIzPG/cWPfYEya4Z41H2O/zf59q0CdXolqw7hORD0TkZyKSoqnqIiK7eduLSCvgBCC5\n2sMY4JLENocBa1U1YnJTw8ict9925vA77ohbkvxw/vnO9eSqq+C+++KWxjDiQ0R28U27ichJQNp3\nlqquTNTMQlXXAx8T7Kr+C2AkUDDH3HQd9kmT3Oerr2Z+7HXrUrvHfPKJGxlO1+nwlBr/yHCY3N73\nUVJv+4/h/U6P1avTW3OSrV1+S1M+8MeDeEqSJ7O/0xbUAfe2/+9/6ydbSBfon+l/nYvL+4IF7nd6\nbnhhVpQvv3TWoGxcuIJ+r1/mbBXlYiR+ysWSHKTkjRvnFNt0smfz2z7+ODxd+tq19a+TZMXfv/6F\nF9K7HQbJGDUjqKesg3uuBMV8hrkheorkypXpY12DzlsoIilYqnok8H2cpWmqiDwpIidE2LUTMEFE\nZgCTgXGqOlZErhCRHyeOPRZYLCKLgAeAn2XzQwwjCps3u2LC99wD7drFLU3+OPRQF5f117/Cn/8c\ntzSGERtTcS6BU4H3gF8DGaWwEZFuuKy3k5O+3wM4U1X/H5BT1GZyR2XKlFr3pELEWi1dGi0o33Of\nCutEegqG12lbtKju+hUrwpNVpOvsJK8PU6aiZlAbObKw6cA9JamQHTR/myxdGty2b7xRX2FJ5ToW\nFS8LpdepTU5A4bk8esr4V19Fb++vv65fJsAfkzYlOVKf7Ns5l6yH+ajHlS4j5ldfpXYbzRTv3n31\n1drrIvmeWb3aXb/Jz4Rki1nyPZluIGDr1mCXXa+IdiorZvL/G/QcTN5fxO335ptuSsYvf1gGyWXL\nwmXKlUguggCqulBE/oB7ed0N9BERAa5X1cCM/Am3wr4B3z+QtHxlRlIbRpbccw/06BHu/lDOdO/u\nfOhPOMGNRN1+e3kl7zCMXFHVvXPZP+EeOBK4KmHJ8nMncJ1/82zOkUuHPLmDo5r+Hn/vvdQp24MI\nk9HrkIV1PD/8MHzfoMxvfqIqli+/DAMHuhijVM/xqJaGxYthjz2goiLa9nHhxUd17173+6BshfnM\nlBZ2fa1fX9dtcVzCafbCC9380UeHu5QFyRhWTy2T+yVqPFxUXnwRunWLtm2yq2yQsjl5MuyzmdpV\nzgAAIABJREFUT/37Np2MI0aAzwM7JX4rlHdPvfwydOzoLDx9+9a6qyYPkKQ6fxS++y74d3tuzwsW\nuMHgbFmxwg28dOnillXrp4v3Pqur67aFf97f3rlkKUxH1BisXsBlwKnA68DpqjotMaL3HlAiJc8M\nI5wvv4S//KX2RdUQ2WMPF6x84olw000uRsswGjrpEk6EDQImHaMZTrl6XFWDQtL7AU8lBhZ3w9V1\nrFbVMUHHGzlyOODVdqqkZ89Kpk3LvsbT4sX1k9k89ZSLJW3XznUsamrqxiBAeuVqxAg4/XTYccfo\nsgRlwPv662Al6dlnoXPn9MdMpWAtWFBXAfLcHHPJMOad8/33XY3BHj2yO0YUBcBzYcoUb4Q+Kl9+\nmVldoUJZ3776ylkGunWDZpGH8euzdGl6xdxPdXXqWKFCsXlz/cK/YW27eLG75i68sNa6tnZtfopS\n19TUdUP0Kw/ePZOu/l2UcyST6jrK5zW2caNzHfYULMhOQfrnP6uoqakKLcCdL6Je+vcA/8RZq7b/\nfaq6ImHVMoyS5+ab4YIL4Hvfi1uSwrLrrs4cfsQR0KED/PzncUtkGAXn9BTrlGiDgI8Ac1X1rsCD\nqO7jzYvIo8CLYcoVwLBhwwEXI+ll2EpWrjKJ5whLqf7dd07BmjXLxRx5SSREgl2PwuIk/J2NbFzM\nliyB5s2D1+XqhjN1KgwaVLucaaetpiY4EYNndaipSV+8OYwosgS5L2V77FT1tRYsiJ7cIh8ZJD0r\nVJDFbMoUp2iFKQ5R2i3qYKhnnVi0yE1t2rjlhQvTK3jbttXWBcuWZOUqFX63PM/lLV2NueT9/Mya\n5QYIwGVc9MdcRY1/CiPZsrZsWa2lMiqZlIIIKkacimytk/vtV8mFF1Zut8w9/3xhkjJFVbBOBb5T\n1W2wPZNSy0RRx8cLIplh5BEv3e/HH8ctSXHYfXf3IDzySGjfHs47L26JDKNwqOpluewvIoNxccaz\nE3UbFbge6OoOrw8mnzLqscPSF0NmnfqwjrXXUfVcc7wR7K1bXVxmNsePkhp84kT3bPGTS8a3TIL8\n811gPRc3oVxifLIhVXbJTDIHLl/uyn2kYsQIOOOM8PVe/FXYezWVVSbsek6VmCQo/u7bb+tbufzK\nSJhFb+tWp3ylSwFfKCvfmjWZKUBhbTxnjlM0Djww//dFsgKTTUxa0D5vv13/u2zaefp02LOEq+ZG\nVbDeAI4HvKbaAXgNGBS6B662CPAY0AGoAR5S1buTtjkaGA14IY/Pq+qtEeUyjEj89rdu2m23uCUp\nHnvv7XyvTzjBWbLKOSW9YURFRE4FDgS2OwqpakpnWVWdBDRNtU3S9v+TtYAFwuvojxrlPjNJ8JBN\nCuvlyzMbnfbIROnz47f+ZdoZW7s2OM2514EMSzMeFhjvpxBJSVKRiQtgOqK0Y1hsVK5s3OgUoYUL\n6/63Ycp9WIHYl15yMXmZ8uyzcPbZ6ZWcbAv/plO8X3vNufRHJZWcH31U120uXyRbJsNkSDWIlEx1\ndd1C4rmwdi20bZufYxWCqApWS3/Ar6quF5E0Yx9AxNoiwERVTTFOYhjZU1XlzOiZPAQaCoccAk8+\n6Vwj338funaNWyLDKBwicj9uAPAYnFv7MCBKHayyJayTHDWAHZw7V7ES4mSbMc0f95VpPNPy5Zm5\nNqlGVzqrq53VJRtlM4jNm93A2JAhwf9tPhWsKL8x25jBdKgGW4/CYq5SFTjOtKiyRyaufYUg10LY\nfmbMyO/xoP4gTSr31Kjk24som0LcUJyBkagK1gYR6auq0wBE5FAg7a2pqiuBlYn59SLi1RZJVrAs\n15lREGpq4Ne/dhn1ih34Wiocf7yz3p15pvNpb906bokMo2AMUtVeIjJLVW8Wkb8BEewQ5U8urkz5\nKuZaqiTXAlqypLbOUxCzZ9ftTC5YEJ5NLtOYlHR4nf58H7eQ+K+9qG5kmV6vqf6vcmLOnNr5XGOk\n/OQr1XumcVD5YkxoNGt98pE+3/8/FIqohYavBp4VkbdF5B3gaSCj1OphtUUSHC4iM0TkZRHpmclx\nDSMVTzwBLVpYDNLVV7tq7j/8YXGKMRpGTHgDfxsTWW6rcfUYjQTpigg3BpI76wsW1F1OHqmfOrWw\n9XKCCLMuLVxYGnKEETWeLZ/KhZE/kuuTZUKU/z5MscnEovTii9G3zVSOfBLJgqWqU0Rkf8DLvzZf\nVSOPeaWpLTIV6KKqG0VkCDAKCEyYOnz48O3zlZWVVEYtDGA0SjZuhBtucK6Bjb0elAg88ICLw/rL\nX+B3v4tbIqMhUVVVRVVylcp4eElE2gH/B5iGS0bxULwiFZZt29yoc9SBk1SuVkY4jXVgKlNF6PXX\nCyOHURxySdoyfnz+5GgIiEZ8aojIIKAbPqVMVR+LsF8z4CXglbD0t0nbLwYOVdWvkr7XqLIaBsCf\n/uRcQ0ou9urqq93nnXcW/dTLlrmA4AcfhFNPLfrpjUaCiKCqsQ5riEgFLn44g3QPeTu3Pvlkcd5X\n3bplXkjYyJwBA+CDBh3Nlz177x2fa5lh5MpFFxXmfRW10PDjQHdgBuAZ8hSXITAdKWuLiEgHVV2V\nmB+AU/oCyhgaRnRWrIC77oIPP4xbktKic2eXPenMM12q1IZeE8xoXIjILOAp4GlV/QTIc+Li0sOU\nq+KQjwD/hkq2iQYMoyETNclFP6BnpiakiLVFhonIT3G+8t8B52dyDsMI4g9/gB/9KDwwuTEzaBDc\ndptTsiZPLu00p4aRIafj3iHPiEgNLl74GVW1LqCRE2ahCafYqeoNoxyIqmB9BHQEMspeH6W2iKr+\nA/hHJsc1jFRMn+7qlxQqvWxD4PLLYdo0uPhieOGF8DowhlFOqOpnwB3AHSKyH3Aj8BcyqHFlGIZh\nGLkStVu1GzBXRMaJyBhvKqRghpENqnDNNTB8uFlm0nHXXa6a/C0pS7AaRnkhIl1F5Lc4V8H9gd/G\nLJJhGIbRyIhqwRpeSCEMI1889ZQrSPnDH8YtSenTogWMHAn9+0OfPjB0aNwSGUZuiMhkoDnwDHCu\nquaQdNgwDMMwsiNqmva3RKQrsJ+qviEiO2AuF0aJsWaNs16NHg3Nog4dNHI6dnRK1mmnQY8ecMAB\ncUtkGDlxiaqac7BhGIYRK5FcBEXkclwdqwcSX+2Jq1eVbr/OIvKmiMwRkdki8suQ7e4WkYWJYsO9\nowpvGH6uvdYVFB4wIG5JyouBA+GOO1zSi3VFT2htGPkjW+UqyrtKRC4SkZmJ6R0ROTh3iQ3DMIyG\nSNRx/p8DA4DJAKq6UER2j7DfVuAaVZ2RKDY8VUReU9V53gaJ4sLdVXU/ERkI3A8cltGvMBo9b77p\nitx99FHckpQnl10GU6fC978PY8ZY0guj0ZH2XQV8ChylqutE5GRcAWN7VxmGYRj1iNqN2qyqW7yF\nRPHgtCnbVXWlqs5IzK8HPsZZv/wMJVFPS1UnAzuJSIeIchkG330HV1wB994LbdrELU358ve/w7ff\nwk03xS2JYRSXKO8qVX3fV7T4/eT1hmEYhuERVcF6S0SuB1qJyAnAs8CLmZxIRLoBvUlYwXzsCSz1\nLS/HXlxGBlx7LfTrB6efHrck5U3z5q4I8b//Dc8/H7c0hpE5IrKDiNwoIg8llvcTkdMyPEY3gt9V\nfn4EvJKtnIZhGEbDJqqL4O+AHwKzgSuAscA/o54k4XIxErgqMTqYFcOHD98+X1lZSWVlZbaHMhoI\no0bByy+72ldG7uy+u1OuhgyB/faDgy3KxIhAVVUVVVVVcYsB8CgwFTg8sbwcNyD4UpSdo7yrROQY\n4DLgiFTHevW+S7bPH9L1EHp3tfBiwzCMuJnx2Qxmfjaz4OcR1bSefrmdwLkTvgS8oqp3Bay/H5ig\nqk8nlucBR6vqqqTttNCyGuXFokUweLBTsg4/PP32JcPVV7vPO++MV44UPPkk3HADTJ7slC7DyAQR\nQVUlhvN+qKr9RGS6qvZJfDdTVQ+JsG/Kd1Vim17Ac8DJqvpJimPp+OvfyO5HGIZhGEXjuP99fEHe\nV5EsWCKymICYK1XdJ8LujwBzw15YwBhcEo2nReQwYG2ycmUYyaxf77Le3XRTmSlXZcJFF8HcuXD2\n2S55SEVF3BIZRiS2iEgrEu8rEekObI64b8p3lYh0wSlXF6dSrjxWHXRcxNMahmEYDY2oLoL9fPMt\ngXOBXdLtJCKDge8Ds0VkOu6ldz3QFVBVfVBVx4rIKSKyCNiAc70wjFCqq+H88+Gww+CnP41bmobL\nLbfAsGHwk5/AI4+AFN0eYRgZcxPwKrCXiPwHGAxcmm6nKO8q4Ebce+8+ERGgWlWtKIRhGIZRj6xd\nBEVkqqoemmd5Up3PXAQNVF1K8dWrXUHh5s3jligLysBF0GPDBjjySDjrLLjxxrilMcqFuFwEE+fe\nFZc+XYD3VfXLGGTQJ5+095VhGEapc9FFhXlfRXUR7OtbbIKzaEW1fhlG3rj+evj4Y1f3qiyVqzKj\ndWuXRGTwYOjYES6/PG6JDKM+Se8ogP8mPruISBdVnVZsmQzDMIzGS1Ql6W+++a3AEuC8vEtjGCm4\n+2544QV45x3X8TeKQ6dOMG4cHHWUS3gxdGjcEhlGPf6WYp0CxxZLEMMwDMOIpGCp6jGFFsQwUvH0\n03DHHTBpEuy2W9zSND722w9efNGlb2/XDo4+Om6JDKMWe0cZhmEYpURUF8FrUq1X1f+bH3EMoz4T\nJsAvfgFvvAFdu8YtTeOlXz946imX+GLkSFOyjNJDRFoCP8PVqFLgbeB+Vd0Uq2CGYRhGo6JJxO36\nAT8F9kxMPwH6Am0SUyAi8rCIrBKRWSHrjxaRtSIyLTH9ITPxjYbO7NkuY+Azz0CvXnFLYxx3nLMm\nDhsGpVFX1jDq8BhwIHAPcG9i/vFYJTIMwzAaHVFjsDoDfVX1WwARGQ68rKo/SLPfo7gX3WMptpmo\nqmdElMNoRCxbBqeeCnfdBZWVcUtjeBx7rFN4zz3XfR5jzllG6XCQqvb0LU8QkbmxSWMYhmE0SqJa\nsDoAW3zLWxLfpURV3wG+TrOZVdcx6rFunVOufv5zuPDCuKUxkjnmGHj2WWddHDcubmkMYzvTEgXr\nARCRgcCHMcpjGIZhNEKiKliPAR+IyPCE9Woy8O88yXC4iMwQkZdFpGf6zY2GTk0N/OAHMGgQ/Pa3\ncUtjhFFZCaNGwcUXu5pkhlECHAq8KyJLRGQJ8B7QX0Rmh7mqG6lp2TL4+2MtL6NhGEYoUbMI3iYi\nrwBHJr66TFWn5+H8U4EuqrpRRIYAo4AeYRsPHz58+3xlZSWV5jfWILntNvjqK3juORCzb5Y0gwbB\n2LHO2rhli3MbNBofVVVVVJVGUN7JcQvQEOjWDZYscfM77ACbGnmKkOOOg/HjC3f8Vq3gu+8Kd/zG\nTMeOsHJl3FIYHh06wKpVcUtRHDIpFrwD8I2qPioi7UVkb1VdnMvJVXW9b/4VEblPRHZR1a+Ctvcr\nWEbD5JVX4P77YcoUaNEibmmMKPTrB6+9Bief7DpiF18ct0RGsUke8Lr55ptjkUNVPxORnYG98L3f\nrNBwelq0cIMkAN271ypYxaJvX5hWpv/SwQe7hExGao48Et5+u3jnO+ssmDevuArW3nvD4px6xqVB\nly7w+efRtu3ZE+ZGjHRVzV6mciOSi6CI3ARcB/w+8VVz4ImI5xBC4qxEpINvfgAgYcqV0fD59FO4\n9FKXCnyPPeKWxsiEQw6BN9+E3/8eHnoobmmMxoqI/AmYBdyNKz78N+CvccjSs0AO74ceWpjjhhVv\nD/s+3x2lZpkM9+aZdG2aypOifXto3jy/8pQSTaIGkkSg2P9xy5auhmMx6dYtf8caMCB/x4pK27bu\nM5//eyGI83kRlahNeBZwBrABQFVXkCI9u4eIPAm8C/QQkc9F5DIRuUJEfpzYZJiIfCQi04E7gfMz\n/gVGg2DjRjj7bLjhBjfKZZQfBxzgUrffeivcc0/c0hiNlPOA7qpaqarHJKa00UIi0llE3hSROYl4\nrV+GbHe3iCxMxA33Tn3MLH9BCLvs4j733Td4fa6d/G7dgjsthx8evH1c7tudOmW2fbrC9Geemf4Y\nYUpmt24uFq1z58xkSmbvvXPb34//P9xxx7rrvGsomaA4u44d3WfYb4/CEUc498ds2Gmn7M8LcPrp\n7rN16+yv1WxlL2dOOCHzfUSix2Tma2CmHDycoipYW1RVcYUbEZFIt5yqXqSqe6hqhap2UdVHVfUB\nVX0wsf4fqnqQqvZR1UGqOjm7n2GUM6rwk5/AgQe6gsJG+bLvvvDWW3DnnXDHHXFLYzRCPgLaZbHf\nVuAaVT0QOBz4uYjs798gESfcXVX3A64A7s9V2CC6d3effoWpVy8YPBiOPz58ZLlN2iHP1Oy/Pxx1\nVN3vmjZ1UxC77163M9a3L5x4Yvrz9OkTvq5HaAR2LZl2lg8+OPX6Vq3qHvOAA+rHku6wQ/C+PXq4\n/8NTQtKdK4wwpTkqBx3kPjt2TB0Hu9de0Y/pld/IRZHea6/Uin8q69Ipp2R/XqirXBbCLS1MMS/U\nwEOxlD1PcYnSZoMHu08RF1sVhaDjHnAA7LNPtP09vOfSrrsGr0/1PAzbJ99EVbCeEZEHgHYicjnw\nBmCOQEZeuO8+mDkTHnzQklo0BLp1g4kT4ZFH4E9/ilsao5HxZ2C6iIwTkTHelG4nVV2pqjMS8+uB\nj4E9kzYbSqKmY2IwcCe/m3sy3gi81/kNIsi60r+/+6yoqPv9jjs6d7RMyaZAu9eZO/DA8G1E6srf\npg20y0a19ZEP98dkBTRKZ8r/3mnSJNz9aOedw/eDaO5offvW/65Vq9xqPXqKXZjc3v+Uqdvq/vtn\nd/3kSjYZKk87zXnBHHdc+t85bFj9/7IQpLJsZmKBOfXU+s+DYhH2/zdvXtuGqe77vfd2sXCpfq+q\nyx6dDSee6MrFJBOmYDVvnp2VLhsiKViq+ldgJPAc8D3gj6pqTkBGzrz7LtxyCzz/fG6uCEZpseee\nzl3wySedy6BhFIl/A38Bbqc2ButvmRxARLoBvXHlSPzsCSz1LS+nvhK2Ha9zFaSknH22GwE/4ICg\n87tPb3Q4KmGd6913r7scZC3p16/ucps2zoIRxaLkp2nTwtct7NWrVgkNwnMN8xTcKK6Tu+wSLabj\nyCOdFdEjkzgQzz2voqKuYtq0aX0F73vfq7tvFOvAqafCwIHR5fEzcGCwK2ifPtGsXn6lMfn68lss\nPCUhnUulX5mIahlp08btt/vuLibYz6BBdZebN48WY+TvpGcz+JtKdq9dzj47/XHatq0f2+VZGAvB\noYfWKkTt2zsLbrIFbdiw2nnvGglSAps3dy6oRx/tloMsWJncR8nPKgj/L8OeryJ15S8UaS8xEWkq\nIhNU9XVVvVZVf6OqrxdeNKOhs3IlnHees3R4bjFGw6FjR5gwAf7zH1OyjKKxUVXvVtUJqvqWN0Xd\nWUR2xA0mXuXPcpsNt9wynJEjh3PLLcOZO7eqzrqKCtdh6dw5XAHwfx/FXadrV/eZrHwkx914671g\ndqjtoPg7SB07ZhbXFUXGVB3rsA5scudr553rduK9Tqy3vxdP1L599EG7XXap71oX5AKmWvs7e/Wq\n24aerEFJDioq4KSTgs993nlOdu+4Ik5h8b8TU7lV7plQ8du2re0Ue5a0o45yFh0/vX2Rg17s1x57\npE7O4Fea27Sp2zk977y6CmH//s69zzuv/7rwLB7ef+Vdb3361A5ENG+eWfzVnqFDHLV07VrfypLu\neh082CkF3jWbSlnylIdkdt45fD/Pxa2iItqgRHLsod9tNVkhDyP5fk4e/PEsOz16OLlbt3YK69Ch\n0WIVTzst/TZB7d6zZ/04xCDX5FTK/lln1b8fvTjCZKqqqrjtNvd8njBheEp5cyGt3qiq20SkRkR2\nUtV1mRxcRB4GTgNWqWqgoVFE7gaG4BJoXOq5aRgNm+pq92D+0Y/cyJvRMPGUrGOOcS/VG26IWyKj\ngfO2iPwZGANs9r6MkqZdRJrhlKvHVTWodPZyXPp3j86J7wK56abhPP206zyNGBFV/GByse57nVl/\nJx6cW9mkSXW3bdfOPZcLxZFHurTZmdC8OWzdGr5+773r1tXxfmMqK1dFhetApkpD3aqV67B9803t\nd02bOsXthBPqu3d6neTDD6+f4t6zqmVC//7wySe1yzvs4JJBJROWhASCFRX/aH/PnuGJL8Jo0qRu\nRz2oI+w/b8uWsD5kqKJDB+eG6B1vzhyn6AUp2127wrffuhqZfg4/PHqKcKhVtMIKaHt06eI+k5WS\ndu2c3N7+bdvWZj3eaadal1QRt+2xxwbf/wMHZiZ3Kg480CnXYaUCdt0V1qxxz5G1a2u/79XLtXkQ\ne+2VWcweOBn23z/4HvcGSjyXQ3/ZxKZN6ytDrVrVv2523jncVbJly/B40WS8siIjRjiL2EMPFaas\nSNQYrPXAbBF5OJFF6e6EYpSOR4GQcZviBQ0bpcc117iRsD/+MW5JjELTsaNL4f7EE66ItGEUkD7A\nYcD/JvM07Y8Ac1X1rpD1Y4BLAETkMGCtqhalZGby6O5ZZwXHHaSiX7/aGKchQ8I7T2GdlEw7W+kI\nszx41hHP3a1799TJIzyrjf/3JXPOOXWXmzSJlvnPH38ydGhtQox0mQmT8XfSwyx16eK7hgxJf+wo\n+JXstm0zdwP1khEkt0HY7wr63vsuWVkL48ILnZtfkJWtefP6LoFBeNebp2AOGuTuI4/DD3fun1Fc\n9vr2rW8xBSeHp0icfnr93+aXc4896rqbevgtUx07RrtOKypcrGfYfXL00c4S1blzcAxgMWjXzrV3\np07RMoEGubsecIBTfKNY04Lo1Cn3jJ+ZENXz8fnElBGq+o6IdE2xSZ2gYRHZSUQ6FOulZcTDP/8J\nr78OkyeXfq0FIz906uSULM+Sdf31cUtkNERUNavIBBEZDHwfN5A4HZcx93qgqzusPqiqY0XkFBFZ\nhPO4uCzVMZs0ydw6n5yFL9nq5JE8+p4qGYWH362uXbu6x03nktWsmetQf/11uDUiU446yo3s77OP\nq4Ho4XVKvXdD69Zum7DR+Y4dXYfLixHZtq3+Ni1auM7shg3w3nvZyRuWSTBT/P+dX8n05PcrJEOH\nwuiELbVU0lJ7cUndu9dPAx+Fk05y/2nPnsVJNOGRbK1p3rz2WuvRo1Z527y53q6Au0eSCxafeWZm\n8UPNm7tkJl98UX/diSfCa6+5/9yzeLVuXVu2JooLbs+e7n54/XUn26hRtes8l03vXguKT0rHuefC\ns89mto//ek5nNfSTHDsKtc+EsIyKbdq4Z5R33uQ2C0okU8jCxykvDRHpoqqfq+q/C3T+sKBhU7Aa\nKO++6zrXb7+de50Lo7zo1Kmuu+Dvf59+H8PIFBE5FTgQ2P46V9VbUu2jqpOAtA4mqnplJrL4YwJa\ntIAtW4KO6T7bt6+b8U41XMHy2G03+PLL1J3vdMH5UeI/vPikffd1nVS/a1OYS1CvXjBrVrRz+xUs\ncJ3B3XarqwxVVNTt/Prb09/hato0ON6pfXs3eccMatf+/aPF9KTjnHNgyhSn7C33OZEOG1bXqpGc\nFn/33ev+l/lS6vbbz10n6Tj7bJdwKsi90nN19dptn30yT60NdZN9ZMtJJ7lyIJs2ZX8MP6k6/p5F\nt3dvmD+/7rpUqdOT77vjj3f3d5MmwRac5Kx32SifTZq4+2bIkFrZ9t8/uK17B1TxS6csNmvmLG8r\nVkSXKdP4+tat3UAIOCv9vHkuy3QUDjvMFWcePdpdZ/m6f7Ilne49CugLICLPqeo5abYvKMOHD98+\n7/lQGuXDsmXuRf2vf0UPyjQaFn4lC0zJaihUVVVR5XeqjwkRuR/YATgG+CcwDPggVqESVFQ4hcCL\n1/A47DDXEU92f4rSAT32WGdRatPGxSkFWW9EguOqRDLP+NesWf2YnT596itYzZrVT4+eSdauoNHr\n5s3rKljnnBMe25YqrqhNm/DBveQMeNmWDWnRwiVJ+OKLunFT6dzhkhNS5Itu3VInsfDwknGEtV9F\nRXhK7p13Dr7+CkGmcWO54Ck+Ubxt/NskXzu77JL6GC1a1L0f01kHVZ3y4o/T8/D/R8nuxd261VUM\n993XKZHz52fuUZTO+tOkSfi9NmSIG2hPtognZ9Ts2dPF7CfHq/Xs6WLy/Hh1+5KfNc2bO8tgsUmn\nYPkvkSzGKtKSUdCwX8EyyotvvnHuMlddlXsBQaO88StZW7a4ODyrf1beJA943XxzYYKGIzBIVXuJ\nyCxVvVlE/ga8EpcwHkOGuJf8mICKXEGxTVEVn6ZNazswLVq4gatFi4K3yxd77RUsn7+zFVTsNp1y\nke4ZcPzxwbVyMn12eG6bn32Wftvdd3cWgGxp3x5OPjn7/QtBOlfDVMpLqvikoHiiTDnooOixfoV0\n7fLzve9Fs9Qdf3xtXNqJJ9ZVYgpVuiDdtX/++fWVpmQLmmet/N73oiXT6ds3emHsVPK1axe8PsiK\ntv/+9WP+olpPTz89XFndZ5/6A175JJ2+qiHzmSDUVdT8xBY0bBSPLVvcg3nwYLj22rilMUqBTp2c\ni8eYMfCTn6TOEmYYGfBd4nOjiOwBVAMRQqoLS7t2rvPSr1/qhA3JHHhg5sVhS4l81ppp1So/tRK9\nOjhROugtWqROkV4sPGvmEUe4z7CU71HYe2/nRpVvPOtBLhx8cO7FqtORqUIuEi3+rX372mNHKW6d\nKdkolJlYpKLeW23a5MeNNojzz69NEe+noiL7c6ayBA4cmF0cYVTSWbAOEZFvcApSq8Q8iWVV1bbh\nu4KIPAlUAruKyOfATUALsgwaNsqPmhr4n/9xN+U995ilwqilQweXqvWcc2rdfeL2mTbKnpdEpB3w\nf4BpuIHBh+IVqRavPlFUDjoo83N4z9iw2KhikmmGu0wZOrThJ0ryWz/CrIeZIFIYBaBq44wNAAAQ\nzUlEQVQYFMtqFUb37sVzT0zuK3kWseRYxHLmqKPq/pYmTVwbN5S6qCkVLFXNaTxCVS+KsE1GQcNG\n+aAKv/2tC2AePz6/bipGw6BNG3jpJaeEH3MMPPdccdOoGg0LVf1TYvY5EXkJaJlp/cZyx3MZTA6a\nLxX22CM4i1o25DIgs+uujTfRUrt2cMEFxTlXroOqYf2GQw/NXNHo1Cm4+GxUGQth+YtKq1ZOud6y\nxb0nS4FMMigGkVwYuKHRwMd+jLhQheuuc4rViy+mzrZjNG5atIDHH3c1MgYMgIkT45bIKDdEpL+I\ndPQtXwI8A/xJRIoYEh8/u+1WuJiPfNC2bW3qaT9hylKhBubato0vHjiTdNVh5JrivFjeJMlJFjLh\nlFNcEpggunbNvIZXZWX6mLpS97Lx3BUrKlyq9TiVvh12qF+TyrPwDRwYXMuqMZGj/mkY9VGF3/zG\nuX+NH1/cjD9GeSICv/udSx177rkulf8vf1n6LzujZHgAOB5ARI4Cbgd+AfQGHsRlEzQKSCr3raOP\nTr3veecFK1Inntgw3x+HHZZ73Gk+lLRisM8+LrV/NvXTki2MuaR2b0hkkm2w0CQPnnft6ibDFCwj\nz2zbBtdcA5MmwRtvFLeQoFH+nHyyq1NzwQXw8svw8MPRs0oZjZqmqvpVYv584EFVfQ7nKjgjRrka\nPVGsaWFWqnKNFUqHv8htthx+eHBdtYZMt27BLn4NEQupKH8K7iIoIieLyDwRWSAi1wWsP1pE1orI\ntMT0h0LLZBSGdetc0PHMmaZcGdmzzz6uIPVRR7mUsP/6V/zBzUbJ01REvAHD44A3fetsILEImLW5\nuFRUhMfZNdT/QqQw4QalFgt0yikuxswobwqqYIlIE+Be4CTgQOBCEQnygJ2oqn0T062FlMkoDAsW\nOLeHbt3g9dcLn2rVaNg0awZ/+IO7lu6+26Un/qAkysUaJcoI4C0RGY1L1f42gIjsCzSqJBdxcOyx\nhUvdbGRO586Nx9KTKxdeWHqJlXbaqfAZOI3CU2gL1gBgoap+pqrVwFNAUD3lBjre0vBRdRaGI45w\nroH33msPBiN/9O4NU6bA5Ze7JBjf/z4sXBi3VEapoaq3Ab8G/gUcobrd5tkEF4tlFJAOHRqu1aQc\n6dDBZWUtFSx2ymiMFFrB2hNY6ltelvgumcNFZIaIvCwiZVxWsXGxZImLmbn7bhg3znWCDSPfNG0K\nl14K8+e7avODBrkYrRkWWWP4UNX3VfUFVd3g+26Bqk6Lsr+IPCwiq0RkVsj6tiIyJvGumi0il+ZJ\ndMNo0Bx9tAsfKAeaNTOF0MgPpeCbPhXooqobRWQIMAoITL45fPjw7fOVlZVUVlYWQz4jiQ0b4M47\n3fTrX7vJrFZGodlxR/jjH+FXv4IHHnB+6vvvD1dc4axbXvpao/hUVVVRVVUVtxi58ihwD/BYyPqf\nA3NU9QwR2Q2YLyJPqGqO+eAKg1mUjFKhnBSWJk3g7LPjlsJoCIgWMHpcRA4DhqvqyYnl3wGqqn9J\nsc9i4FBfRijvey2krEZ6tm6FRx+F4cOdS+Btt8G++8YtVRly9dXu884745WjzNm8GUaNcsrWnDlw\n0UVw/vmu9oZ1LuNFRFDVsvsXRKQr8KKq9gpY9zugs6peKSJ7A+NUNXAwsBTeV5984uIWS7kmlmEY\nRtwU6n1VaBfBKcC+ItJVRFoAFwBj/BuISAff/ACc0vcVRslQUwNPPQUHHgj/+Q+88AI8/bQpV0a8\nVFQ4herNN+Htt10mqEsvdYlWrr3WxW7ZmIyRR+4FeorICmAmcFXM8hiGYRglSkFdBFV1m4hcCbyG\nU+YeVtWPReQKt1ofBIaJyE+Balz2p/MLKZMRHVUYPRpuvNFV7L73Xjj+eLMOGKVHjx5w883OuvrR\nR/DMMy4hxpYtcM45bjrsMOf+YRhZchIwXVWPFZHuwOsi0ktVA0uoxu3S3qVLUU9nGIZRFhTLpb2g\nLoL5pBRcLhoLqi5pxR/+4NwCb70VTj3VFKu8YS6CRUHVKVsjR8Jzz8HXXzvf+mHDnIurFXIsHA3U\nRfAl4M+qOimxPB64TlU/DNjW3leGYRhlQLm6CBplRlUVHHmkS7l+3XUwbRqcdpopV0b5IQIHH+ws\nWx99BOPHu9owV1/tavb87GcwYQJs2xa3pEYJIYSXDfkMOB62u7b3AD4tklyGYRhGGWEWLAOA995z\nroBLljg3qwsvtBH+gmEWrNhZtMhZtp55BlascJatc8+Fo46y6z4flKMFS0SeBCqBXYFVwE1ACxLu\n7CLSCVdnq1Nilz+r6oiQY9n7yjAMowwo1PvKFKxGTE0NvPQS/PWv8NlnziXw0kst5XrBMQWrpPCU\nrWefheXLnbJ11llO2Sqn9MKlRDkqWPnE3leGYRjlgSlY9sLKG0uWuKyAjz7qagtde62LS2lWClXR\nGgOmYJUsn3zilK0xY2D2bKdknXyymyxrZnRMwbL3lWEYRjlQtjFYInKyiMwTkQUicl3INneLyEIR\nmSEivQstU74ol8Ka334Lb7wBf/oTDB4M/fs7Jevhh+HDD+GCCwqjXJVL+8SJtVF6itlG3bu72MNJ\nk9w9csklLg7xqKOgc2c47zz4+9/h/fddHa5Swa4jo5yw6zU91kbpsTZKj7VRfBRUwRKRJrjaIScB\nBwIXisj+SdsMAbqr6n7AFcD9hZQpn8R54dbUwJo1MHeuS0zx9NNw991www3wox/BGWe4gqvdurnA\n/ptvhvXrnRvgihVw//0uk1ohk1fYjZ0ea6P0xNVGu+ziFKpHHnGug2+9BaefDgsXwk9/6tb36uVi\nt2680dWImzLF3V9btxZXVruOjHLCrtf0WBulx9ooPdZG8VFop7ABwEJV/QxARJ4ChgLzfNsMBR4D\nUNXJIrKTiHRQ1VUFli1WamrcCPjmzbBpU+20YQN89ZWb1qyp/VyzBlavhlWr3OcXX0Dr1tChQ91p\n991hwIC6y507Q4sWcf9iwyhfRJx1q3t3uPhi992GDbBgAcyb56YXX3TLK1a4+3XXXaFTJ2jf3rni\n+qeKCmc1btrUfUaZT7V+6VKYPDn7/Zs2tUyhhmEYhpEvCq1g7Qks9S0vwyldqbZZnviuZBWs556D\nO+5wnZoxY5yyFDRt21Z3ecuWWkWqutopPS1b1p1atXIj47vu6qZddoG99oJDDnHKkl9xMqXJMOKj\ndWvo08dNyWzd6gZC/vtfNxiyYYOzIHvT5s3u+bBli9t227a6n+m+S16/eLGLGctkH/9327bVKlp+\npUvVTc2bw5dfFr+NDcMwDKMcKWiSCxE5BzhJVX+cWP4BMEBVf+nb5kVcutt3E8tvAL9V1WlJx7KI\nYcMwjDKhsSe5iFsGwzAMIxqFeF8V2oK1HOjiW+6c+C55m73SbNOoX9aGYRhG+WDvK8MwjMZNobMI\nTgH2FZGuItICuAAYk7TNGOASABE5DFjb0OOvDMMwDMMwDMNomBTUgqWq20TkSuA1nDL3sKp+LCJX\nuNX6oKqOFZFTRGQRsAG4rJAyGYZhGIZhGIZhFIqyKTRsGIZhGIZhGIZR6hS80HC5IyIPi8gqEZkV\nsr6tiIxJFEmeLSKXFlnEWBGRziLypojMSfz+X4ZsV5bFpPNBlDYSkYtEZGZiekdEDo5D1riIeh0l\ntu0vItUicnYxZYybDO61ShGZLiIficiEYsvZmBGRk0VknogsEJHr4panWIRdmyKys4i8JiLzRWSc\niOzk2+f3iXfCxyJyou/7viIyK9GGd8bxewqJiDQRkWkiMiaxbG3kI1Gq59nEb54jIgOtjeoiIr9K\nPN9nich/RKRFY2+joL56Ptsk0cZPJfZ5T0T8+SWCUVWbUkzAEUBvYFbI+t/jsiAC7AasAZrFLXcR\n26cj0DsxvyMwH9g/aZshwMuJ+YHA+3HLXYJtdBiwU2L+ZGuj+m2UWNcEGA+8BJwdt9yl1kbATsAc\nYM/E8m5xy91YpsS1uQjoCjQHZgRdww1xCrs2gb/gsgIDXAfcnpjvCUzHhSl0S7Sb51EzGeifmB+L\ny0Qc+2/MY1v9CngCGJNYtjaq2z7/Ai5LzDdLPNOsjWrbZw/gU6BFYvlp4H819jYioK+ezzYBfgrc\nl5g/H3gqnUxmwUqDqr4DfJ1qE6BNYr4NsEZVtxZcsBJBVVeq6ozE/HrgY1wdMz91ikkDO4lIh6IK\nGiNR2khV31fVdYnF95PXN3QiXkcAvwBGAquLKF5JELGNLgKeU9Xlie2selXxGAAsVNXPVLUaeAr3\n7GvwhFybnXG//9+Jzf4NnJmYPwPXQdmqqkuAhcAAEekItFHVKYntHvPtU/aISGfgFOCfvq+tjRKI\nSFvgSFV9FCDx29dhbZRMU6C1iDQDWuEybzfqNgrpq+ezTfzHGgkcl04mU7By516gp4isAGYCV8Us\nT2yISDfcCMLkpFVhxaQbHSnayM+PgFeKIU8pEtZGIrIHcKaq/j+gUafBTnEd9QB2EZEJIjJFRC4u\ntmyNmOTn3DIa4XPOd22+D3TQRFZgVV0J7J7YLOydsCeu3TwaWhv+HbgWNzDrYW1Uy97AlyLyaMKN\n8kER2QFro+2o6grgb8DnuN+7TlXfwNooiN3z2Cbb91HVbcBaEdkl1clNwcqdk4DpqroH0Af4h4js\nGLNMRSfxm0cCVyVGMI0korSRiByDy6TZaOI3/KRpozup2y6NUslK00bNgL44t9yTgRtFZN8ii2g0\nUgKuzeQsWo02q5aInAqsSlj6Uj27Gm0bUfv8+oeq9sVllv4ddh1tR0Ta4awpXXHugq1F5PtYG0Uh\nn22Stv9hClbuXAY8D6CqnwCLcb7njYaEmXok8Liqjg7YJFIx6YZMhDZCRHoBDwJnqGoqt9QGSYQ2\n6gc8JSKLgWG4wYwziilj3ERoo2XAOFXdpKprgInAIcWUsRGzHPAHPjeq51zItbnKcwdPuN94rr1h\n74SG/K4YDJwhIp8CI4BjReRxYKW10XaWAUtV9cPE8nM4hcuuo1qOBz5V1a8SlpQXgEFYGwWRzzbZ\nvk5EmgJtVfWrVCc3BSsaQri2+hnugifxR/bABSA2Jh4B5qrqXSHrrZh0mjZKZKR5Drg4oag3RlK2\nkaruk5j2xnXkfqaqyYXLGzrp7rXRwBEi0jThWjMQFw9jFJ4pwL4i0lVEWgAX4J59jYWga3MMcGli\n/n/hrk/v+wsSmbn2BvYFPki48awTkQEiIrj3RuCAVLmhqterahdV3Qd3bbypqhcDL2JtBECiX7BU\nRHokvjoOl7THrqNaPgcOE5GWid92HDAXayOo31fPZ5uMSRwD4FzgzXTCFLTQcENARJ4EKoFdReRz\n4CagBYlCycCtwL98qSF/m06rbUiIyGDg+8BsEZmOM8FejzNfq1ox6UhtBNwI7ALcl7ixq1V1QFwy\nF5uIbeSn0bk/RLzX5onIOGAWsA14UFXnxiZ0I0JVt4nIlcBruMHLh1W1USi3Ka7NvwDPiMj/4AYj\nzwNQ1bki8gyuY1iNGyzx7umf4zLJtQTGquqrxfwtMXA71kZ+fgn8R0Sa4warL8MldbA2AlT1AxEZ\nicuCV534fBCXZK3RtlFIX/124Nk8tcnDwOMishCXLfyCtDLVHtMwDMMwDMMwDMPIBXMRNAzDMAzD\nMAzDyBOmYBmGYRiGYRiGYeQJU7AMwzAMwzAMwzDyhClYhmEYhmEYhmEYecIULMMwDMMwDMMwjDxh\nCpZhGIZhGIZhGEaeMAXLMAzDMAzDMAwjT/x/535FLhKeclQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x12bc70240>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#generate a set of data\n",
    "N = 200\n",
    "mu_0 = 10.\n",
    "sigma_0 = 2.\n",
    "y = np.random.normal(loc=mu_0, scale=sigma_0, size=N)\n",
    "\n",
    "#run MCMC\n",
    "nsamples = 10000\n",
    "with pm.Model() as model:\n",
    "    mu = pm.Uniform('mu', lower=-20, upper=20)     # a simple uniform prior\n",
    "    sigma = pm.Uniform('sigma', lower=0, upper=10)\n",
    "    y_obs = pm.Normal('Y_obs', mu, sigma, observed=y) # we use the canned distributions in PYMC3\n",
    "   \n",
    "    start = pm.find_MAP()\n",
    "    step = pm.NUTS(scaling=start)\n",
    "    \n",
    "    trace = pm.sample(nsamples, step, start, random_seed=123, progressbar=True)\n",
    "lines = {var:trace[var].mean() for var in trace.varnames}\n",
    "pm.traceplot(trace, lines= lines)\n",
    "pm.df_summary(trace)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**traceplot** provides\n",
    "- marginal posterior distribution for each parameter (left)\n",
    "- the timeseries of the parameter values from the chain (right)\n",
    "\n",
    "Right hand side shows broad oscillation: both inlying and extreme values occur frequently but at irregular intervals indicating the markov sampler chain was 'well-mixed' (sampling evenly around the optimal position)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bayesian Statistical Inference\n",
    "\n",
    "As you have seen in the earlier lectures (and will likely see many times again) in Bayesian inference, we evaluate the **posterior probability** by using\n",
    "** data likelihood** and **prior** information: \n",
    " \n",
    "$$p(M,\\theta \\,|\\,D,I) = \\frac{p(D\\,|\\,M,\\theta,I)\\,p(M,\\theta\\,|\\,I)}{p(D\\,|\\,I)},$$\n",
    "\n",
    "The prior can be expanded as \n",
    "$$p(M,\\theta\\,|\\,I) = p(\\theta\\,|\\,M,I)\\,p(M\\,|\\,I).$$\n",
    " \n",
    "It is often the case that the ** evidence** $p(D\\,|\\,I)$ is not evaluated explictly since the posterior probability \n",
    "can be (re)normalized. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**The Bayesian Statistical Inference process** is then\n",
    "* formulate the likelihood, $p(D\\,|\\,M,\\theta,I)$\n",
    "* chose a prior $p(M,\\theta\\,|\\,I)$, which incorporates *other information beyond the data in $D$*\n",
    "* determine the posterior pdf, $p(M,\\theta \\,|\\,D,I)$\n",
    "* search for the model parameters that maximize $p(M,\\theta \\,|\\,D,I)$ \n",
    "* quantify the uncertainty of the model parameter estimates (credible region)\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    " \n",
    " ### Estimating $p(\\theta)$\n",
    " \n",
    "We want to estimate $p(\\theta)$ give data, a likelihood, and priors. For the simplest (low dimensional) case we could simply undertake a grid search to evaluate $p(\\theta)$ for all $\\theta$. As you might expect this becomes quite slow for large dimensional data sets or for large samples. \n",
    "\n",
    "Once we have the posterior we wish to evaluate the multidimensional integral\n",
    "\n",
    "$$I(\\theta) = \\int g(\\theta) p(\\theta) \\, d\\theta$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Three examples of this are:\n",
    "\n",
    "- Marginalization:  if $P$ elements of $\\theta$ are the model parameters, and the next $k-P$ parameters are nuisance parameters (often used to capture systematics within the data) we integrate over the space spanned by nuisance parameters $\\theta_j$, $j=(P+1), \\ldots, k$. For this this case, $g(\\theta)=1$. \n",
    "- Estimating the posterior mean: $g(\\theta)=\\theta_m$ and the integral is performed over all other parameters\n",
    "- Estimating credible regions (interval enclosing $1-\\alpha$ of the posterior probability): we evaluate the integral up to an unknown normalization constant because the posterior pdf can be renormalized to integrate to unity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Monte Carlo Methods \n",
    "\n",
    "Monte carlo methods have been used prior to computers including Buffon's needle (Comte de Buffon) that was proposed as a method for estimating $\\pi$\n",
    "\n",
    "<center><img src=\"figures/220px-Buffon_needle.png\" alt=\"buffon\" style=\"width:300px;\"/></center>\n",
    "\n",
    "Probability a needle will cross a line\n",
    ">$P = \\frac{2l}{t\\pi}$\n",
    "\n",
    "The modern version was invented in the late 1940s by Stanislaw Ulam, while he was working on nuclear weapons projects (studying the distance neutrons traveled through different materials) with von Neumann developing techniques for pseudorandom numbers. The name Monte Carlo was given to the method by Nick Metropolis (inventor of the Metropolis sampler, which evolved into one of the most famous MCMC algorithms, the Metropolis-Hastings algorithm). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Assume that you can generate a distribution of M random numbers $\\theta_j$ uniformly sampled \n",
    "within the integration volume V. Then our interval can be evaluated as \n",
    "$$ I = \\int g(\\theta) \\, p(\\theta) \\, d\\theta = \\frac{V}{M} \\sum_{j=1}^M g(\\theta_j) \\, p(\\theta_j).$$\n",
    "    \n",
    "Note that in 1-D we can write a similar expression \n",
    "$$ I = \\int f(\\theta) \\, d\\theta = \\Delta \\, \\sum_{j=1}^M f(\\theta_j).$$\n",
    "\n",
    "where $ f(\\theta) = g(\\theta) \\, p(\\theta) $, and it is assumed that the values\n",
    "$\\theta_j$ are sampled on a regular grid with the step $\\Delta = V/M$ ($V$ here is the\n",
    "length of the sampling domain). This expression is the simplest example of\n",
    "numerical integration (\"rectangle rule\", which amounts to approximating $f(\\theta)$\n",
    "by a piecewise constant function)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The reason why we expressed $f(\\theta)$\n",
    "as a product of $g(\\theta)$ and $p(\\theta)$ is that, as we will see shortly,\n",
    "we can generate a sample drawn from $p(\\theta)$ (instead of sampling on a \n",
    "regular grid), and this greatly improves the performance of numerical integration.\n",
    "\n",
    "One of the simplest numerical integration methods is generic Monte Carlo. We generate a random set of  M values  $\\theta$,  uniformly sampled within the integration volume  $V(\\theta)$ as ${ I \\approx { V\\theta \\over M} \\sum{j=1}^M g(\\theta_j) \\, p(\\theta_j). }$ This method is very inefficient when the integrated function greatly varies within the integration volume, as is the case for the posterior pdf. This problem is especially acute with high-dimensional integrals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fractional integration error for pi:\n",
      "M:  10 0.0185916357881 0.145915590262\n",
      "M:  100 0.0949860084722 0.0122254380645\n",
      "M:  1000 0.00585924034078 0.0255944532842\n",
      "M:  6000 0.00373717443289 0.00910132430986\n",
      "M:  500000 0.000601177109208 0.00168215748269\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARgAAAEPCAYAAACZXduvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFARJREFUeJzt3XusnHWdx/H3Z3ootFDAXuzaQusG0LMQkWWxtgimcrEt\n0UBwA4WNLmRZ+cMSknUja7IGNG6MZjWusmqqLMZEtiCoW1RuKxwBbUvdBYrS0xZ0S2+2lJv0Bqdn\nvvvHPJXpcC5z5sxvnmdmPq/kSWfm/M4zn5yc8+3v+3suo4jAzCyFUt4BzKxzucCYWTIuMGaWjAuM\nmSXjAmNmybjAmFkySQuMpFsk7ZS0boQxX5O0SdITks5ImcfMWiv1DOZWYNFwX5S0BDgpIk4BrgW+\nlTiPmbVQ0gITEY8CL40w5GLge9nYNcBxkmamzGRmrZP3GsxsYEvV823Za2bWAfIuMGbWwXpyfv9t\nwIlVz0/IXnsTSb5oyiwnEaFGvq8VMxhl21BWAh8DkDQfeDkidg63o8hhG8y2g8AeehgEHmE+03i+\nju++MafUjW7tlteZW7M1LukMRtJtwEJgmqTnqPxkJwIREcsj4meSLpL0DLAXuDplnkZUV+CjOQjA\nOaxmOzPooVJ8AJ5jLku4j028s9URzQoraYGJiCvrGLMsZYZUJmb/HipAJ7GZfnoPG3MjMI81bGMO\n/8y/8ALTWxnRLHdql/vBSGqTpG/oozJ9AyhTmWwKWMOZfJj7Clhw+ngjcbvow5lTU8NrMC4wORkE\nyohHeB/TeZm/5k63V1ZQLjBtrwwcoIcjKHMFt3EXl+cdySzjAtNRylWPi9tOWfdwgelog8AWZnOA\nyaxmAf/Il11wrIVcYLpKGdjHJM7nQR5jft5xrOM1XmB8qUAbKgHHsJ9VLGAQ8X/M4hQ25B3L7E1c\nYNpYKdvmsoN+ehlEHETcy/lMY3fe8czcInWqMtDHOVzGj7xeY+PkFslqlIDzeJRdzGAzb+NW/taz\nGms5z2C6SBnYyyQu8OKwjYlnMFaHEjAlWxyurNVc4FmNJeUZTJfzWo2NzjMYa1D1Ws0gYhlfyTuS\ndRDPYOwwh676/m/O42+43bMaw2fyWhJlYA+TuZCfe1G4q7lFsgRKwLHsyxaFS3yE2/OOZG3GMxir\nWxnYxTTezy9975qu4hbJWqgMvMjxnM1qF5qu4BbJWqgETOdl+unldeTWyYblGYyNWxk4SIkP8VMe\nYHHecazpPIOxHJWAiZS5lyW+ktsO4wJjTVMCFvEgu5jBDXwu7zhWAG6RLIlD9xXu41wu44c+Ya+t\nuUWygjl0M6zzeIQdvJXbuNytUxfyDMZaZjszOJ2nPZtpO57BWBuYxfPsYgbPMdv3EO4SLjDWUiXg\nRLbTTy8Xcm/ecSwxt0iWmzLwPMdxLmt8RnChuUWyNlQCZvIK6+nlTi7xInAH8gzGCmMfPZzBbzyb\nKRzPYKwDTOYg/fT6bOAO4gJjhVJ9NrBv39n+3CJZYZWBXvrdMuXOLZJ1oBLQTy+vMpl5rM47jjXA\nBcYKrQQck32Wk+87036SFxhJiyX1S9oo6YYhvn6spJWSnpD0lKSrUmey9lMC7mApq/grLwC3kaRr\nMJJKwEbgfGA7sBZYGhH9VWM+DRwbEZ+WNB3YAMyMiIM1+/IajAHwAm9hAau8NtMyxV2DmQdsiojN\nETEArAAurhkTwJTs8RTghdriYlZtGi/RTy8Ps8CzmYJLXWBmA1uqnm/NXqt2M3CqpO3Ak8D1iTNZ\nBygB57Kanczw2kyB9eQdAFgEPB4R50k6CXhA0ukRsad24E1Vjxdmm3W3CVTWZhZznO8H3DR92TZ+\nqQvMNmBO1fMTsteqXQ18ASAinpX0e6AX+HXtzm5Kk9HaXAm4lyX8ivdyCT/x/WbGbSGH//f92Yb3\nlLpFWgucLGmupInAUmBlzZjNwAUAkmYC7wB+lziXdZgScA5r+D1zfK+ZAklaYCJiEFgG3A/8FlgR\nEeslXSvp49mwzwNnS1oHPAB8KiJeTJnLOtcU9tNPL1/iuryjGL5UwDpUGdjNdD7M3TzG/LzjtDl/\ndKzZYcpUpud7OIop7M87Tpsr7nkwZrk49Is9mQO8ylG+liknLjDW0SrXMr3GKhZwFd/OO07XcYtk\nXaMMLGCV12TGzC2S2ahKwCoW8HWuyTtK1/AMxrpOGbiTD3E5d+cdpU34KJLZmJSBvUziAh50yzQq\nFxizhrzC0RzPmy57s8N4DcasIVPYy6sc4cPYiXgGYwYcBE71DcaH4RbJbNy2MpMT+UPeMQrILZLZ\nuM1iJ68y0e1SE3kGY1bD7VItz2DMmqYHWE8vF3Jv3lHanguM2RAmULlLnj++dnzcIpmNoAz8Hcv5\nLn+fd5Qc+SiSWTK+SNJrMGbJHLpI8rO86YNJbRSewZjVqXvbJbdIZi1RBi5jBXdxed5RWsgFxqxl\num9NxgXGrKVeZRLHsi/vGC3iRV6zljqa/RxEvs/vKDyDMRuH7miX3CKZ5WYAOK2jr11ygTHL1XO8\njblszztGIl6DMcvVCezwdUtD8AzGrEnKwJOcyoX8gheYnnecJvIMxix3JeAveZrvd9VJeCNzgTFr\nsgt5kI9we94xCsEtklkCnXX42keRzApnD0cwhdfzjtEEXoMxK5zJDHADn8s7Rq48gzFLqAxcz5e5\nmX/IO8o4uEUyK6wy0NvWZ/oWuEWStFhSv6SNkoa8JZikhZIel/QbSQ+lzmTWSiXgSd7FKWzIO0rL\nJZ3BSCoBG4Hzge3AWmBpRPRXjTkO+BXwwYjYJml6ROweYl+ewVhb282xzOCVvGM0oLgzmHnApojY\nHBEDwArg4poxVwJ3RcQ2gKGKi1knmMof+RVnMY3u+RVPXWBmA1uqnm/NXqv2DmCqpIckrZX00cSZ\nzHJRAhbwPzzEOXlHaZkiHKbuAc4ElgCLgc9IOjnfSGbpnMaGrvnUyJ7E+98GzKl6fkL2WrWtwO6I\nOAAckPQw8G7gmdqd3VT1eGG2mbWbEnAPFzGTXQW9KLIv28Yv9SLvBGADlUXeHcBjwBURsb5qTC/w\ndSqzlyOBNcDlEfF0zb68yGsd5decznt4Mu8YdSjoIm9EDALLgPuB3wIrImK9pGslfTwb0w/cB6wD\nVgPLa4uLWSc6k3Udf1GkT7Qzy1EZWMw9PMDivKOMwGfymrWtQWAmzxd0PQYK2yKZ2egmAN9nad4x\nknCBMSuAC/l5Rx66dotkVhAHgT8rZKvkFsms7fUAP+JDecdoKhcYswJ5H2s66iZVbpHMCqYMTKBI\nv+1ukcw6SqecgOcZjFkBFetTCXyinVnHeYWjOZ49ecfALZJZB5rC3rb/vGsXGLOCKgH/xifbusi4\nRTIruPyPKrlFMmu6ct4Bqsxjdd4RGuIZjFkbeI0Ss9mZ02UEnsGYdbQjKbflFdeewZi1ifzOjUk4\ng5F0naS3NLJzM2ueEvAwZ7fV5yrV0yLNBNZKuiP7GNiGKpmZjd+RBF/nE3nHqFtdLVJWVD4IXA2c\nBdwB3BIRz6aNd1gGt0hmwKtM4lj2tfAdEy/yRqUK/SHbDgJvAe6U9KVG3tTMGnc0+7mKb+cdoy6j\nzmAkXQ98DNgNfAf4cUQMZB9svykiTkof0zMYs2pl4K0tu/td2hnMVODSiFgUET/IPsSeiChDh91+\ny6xNlIDb2uCwtQ9Tm7WpMtBLP5t4Z+J38u0azLrSDqYyixcSv4vP5DXrSjN5sdDXKXkGY9bmXuQY\npvFqwnfwDMasax3PnsLOYjyDMesAe5nIMbyWaO+ewZh1tUm8XsiT7zyDMesQg0BPkjvfeQZj1vUE\nhbvS2gXGrEOUgO9zed4xDuMWyayDpDm712fymllmC29lDjubuEevwZhZZja7CrMWk7zAZHfB65e0\nUdINI4x7j6QBSZemzmTWyUrADyjGn1HSFim7Z8xG4HxgO7AWWBoR/UOMewDYD/xHRPxwiH25RTKr\nU3PvF1PcFmkelZtSbc7uI7MCuHiIcdcBdwK7Eucx6wol4G4W5R0jeYGZDWyper41e+1PJM0CLomI\nb1I5lG9mTfBe/jfvCIVY5P0qUL024yJj1iR5XwTZk3j/24A5Vc9PyF6rdhawIvvkgunAEkkDEbGy\ndmc3VT1emG1mNrQS0Mf7OZHtY1yL6cu28Uu9yDsB2EBlkXcH8BhwRUSsH2b8rcDdXuQ1a54HOZfz\neXgce2h8kTfpDCYiBiUtA+6nUlBviYj1kq6tfDmW135Lyjxm3Wghj+T23j6T16zDDQATx/V/d3EP\nU5tZziZAbveK8QzGrAuM714xnsGY2QgEnMKG1r+vZzBm3WEHM5jV0MnynsGY2Shm8nzLr7J2gTHr\nEiXgjhZfZe0WyayLlIEJY17sdYtkZnVq5WKvC4xZFykBK/lwy97PLZJZl3kdmDWmm1G5RTKzOk0E\n/pVPtuS9PIMx60KbeRtvZ3udo/2xJWY2BmO7dMAtkpmNgYBlfCX9+3gGY9ad6j8nxjMYM2tA6ksH\nXGDMulQJWM41Sd/DLZJZFzvABCZxcJRRbpHMrAE9DCbdvwuMWRcrkfbaJBcYsy5WAn7KRcn27zUY\nsy73OiWOHLFV8hqMmTWoh3KyfbvAmBkf4fYk+3WLZGYMUGLisG2SWyQzG4fXErVJLjBmxmTStElu\nkcwMgIOII4acybhFMrNxUsMfLTs8FxgzAyr3iGn21dUuMGYGVIpBs+/V6wJjZn9yLr9o6v68yGtm\nfzL0ZQNe5DWzJmj2ZQMuMGZ2mHmsbtq+3CKZ2WFeYgpT+WPVKwVukSQtltQvaaOkG4b4+pWSnsy2\nRyW9K3UmMxveZPY1bV9JC4ykEnAzsAg4DbhCUm/NsN8B74+IdwOfB76dMpOZjWywiSfcpZ7BzAM2\nRcTmiBgAVgAXVw+IiNUR8Ur2dDUwO3EmMxvBUU1c6E1dYGYDW6qeb2XkAnINcE/SRGY2qmbdp7en\nKXtpAkkfAK4GzhluzE1Vjxdmm5k118PAIuazievHva+kR5EkzQduiojF2fN/AiIivlgz7nTgLmBx\nRDw7zL58FMmsRQ4/4a64R5HWAidLmitpIrAUWFk9QNIcKsXlo8MVF7OhpLuTrDXrhLukLVJEDEpa\nBtxPpZjdEhHrJV1b+XIsBz4DTAW+IUnAQETMS5nLOoPPEi0+n2hnZm9SBi5jBXdxOeNpkVxgzGxI\nA4iJlCnyGoyZtan9TBr3PlxgzGxI+zly3Ptwi2RmQyoDEwjcIplZEuM9o9cFxsyGVAJ+zCXj3oeZ\n2ZAmcHBc3+8CY2bD2sfEcX2/F3nNbFiVhV68yGtmxeMCY2bJuMCYWTIuMGaWjAuMmQ1rvAXCBcbM\nknGBMbNkXGDMLBkXGDNLxgXGzJJxgTGzZFxgzCwZFxgzS8YFxsyScYExs2RcYMwsGRcYM0vGBcbM\nknGBMbNkXGDMLBkXGDNLxgXGzJJxgTGzZFxgzCwZFxgzS8YFxsySSV5gJC2W1C9po6QbhhnzNUmb\nJD0h6YzUmcysNZIWGEkl4GZgEXAacIWk3poxS4CTIuIU4FrgWykztVJf3gHGqC/vAA3oyztAA/ry\nDtBCqWcw84BNEbE5IgaAFcDFNWMuBr4HEBFrgOMkzUycqyX68g4wRn15B2hAX94BGtCXd4AWSl1g\nZgNbqp5vzV4bacy2IcaYWRvyIq+ZJdOTeP/bgDlVz0/IXqsdc+IoYwBQU6O1xmfzDjBG7ZYXnLnI\nUheYtcDJkuYCO4ClwBU1Y1YCnwBulzQfeDkidtbuKCLasb6YdbWkBSYiBiUtA+6n0o7dEhHrJV1b\n+XIsj4ifSbpI0jPAXuDqlJnMrHUUEXlnMLMOVbhF3nY7MW+0vJKulPRktj0q6V155KzJNOrPOBv3\nHkkDki5tZb5hstTze7FQ0uOSfiPpoVZnrMky2u/FsZJWZr/DT0m6KoeYtZlukbRT0roRxoztby8i\nCrNRKXjPAHOBI4AngN6aMUuAn2aP3wusLnje+cBx2ePFeeatN3PVuJ8DPwEuLXpm4Djgt8Ds7Pn0\nguf9NPCFQ1mBF4CenH/O5wBnAOuG+fqY//aKNoNptxPzRs0bEasj4pXs6WryP8ennp8xwHXAncCu\nVoYbRj2ZrwTuiohtABGxu8UZq9WTN4Ap2eMpwAsRcbCFGd8kIh4FXhphyJj/9opWYNrtxLx68la7\nBrgnaaLRjZpZ0izgkoj4JsU4O6Cen/M7gKmSHpK0VtJHW5buzerJezNwqqTtwJPA9S3KNh5j/ttL\nfZjaMpI+QOUI2Tl5Z6nDV4HqdYMiFJnR9ABnAucBRwOrJK2KiGfyjTWsRcDjEXGepJOABySdHhF7\n8g7WTEUrME09Ma8F6smLpNOB5cDiiBhpCtoK9WQ+C1ghSVTWB5ZIGoiIlS3KWKuezFuB3RFxADgg\n6WHg3VTWQlqtnrxXA18AiIhnJf0e6AV+3ZKEjRn7316ei0pDLCJN4I3FsYlUFsf+ombMRbyx0DSf\nfBd568k7B9gEzM/751tv5prxt5L/Im89P+de4IFs7GTgKeDUAuf9d+DG7PFMKq3H1AL8frwdeGqY\nr435b69QM5hosxPz6skLfAaYCnwjmxEMRMS8gmc+7FtaHrI2QH2/F/2S7gPWAYPA8oh4uqh5gc8D\n3606JPypiHgxj7yHSLoNWAhMk/QccCOVAtnw355PtDOzZIp2FMnMOogLjJkl4wJjZsm4wJhZMi4w\nZpaMC4yZJeMCY2bJuMCYWTIuMJaMpLOyG21NlHR0diOoU/POZa3jM3ktKUmfAyZl25aI+GLOkayF\nXGAsKUlHUPl0if3A2eFfuK7iFslSmw4cQ+WubUflnMVazDMYS0rSfwH/Cfw5MCsirss5krVQoW7X\nYJ0lu23l6xGxQlIJ+KWkhRHRl3M0axHPYMwsGa/BmFkyLjBmlowLjJkl4wJjZsm4wJhZMi4wZpaM\nC4yZJeMCY2bJ/D+wSRwtDlIH0QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x12bf3f518>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "print (\"fractional integration error for pi:\")\n",
    "for M in [10, 100, 1000, 6000, 500000]:\n",
    "    x = stats.uniform(0, 1).rvs(M)\n",
    "    y = stats.uniform(0, 1).rvs(M)\n",
    "    z = stats.uniform(0, 1).rvs(M)\n",
    "    \n",
    "    r2 = x*x+y*y\n",
    "    r3 = x*x+y*y+z*z\n",
    "\n",
    "    print ('M: ', M, np.abs((np.pi-4.0*np.size(x[r2<1])/M)/np.pi), np.abs((np.pi-6.*np.size(x[r3<1])/M))/np.pi)\n",
    "\n",
    "fig, ax = plt.subplots(subplot_kw=dict(aspect='equal'))\n",
    "plt.plot(x, y, '.k', markersize=3, c='blue')\n",
    "plt.plot(x[r2<1], y[r2<1], '.k', markersize=3, c='red')\n",
    "ax.axis([0, 1, 0, 1], aspect='equal');\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Exercise: Extend this technique to >2 dimensions\n",
    "\n",
    "Estimate how many samples would be required to reach a 1% estimate of $\\pi$ (for the adventurous see Chapter 7 in \"A Practical Python Guide for the Analysis of Survey Data\" for formula for an n-dimensional hypersphere)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "##  Markov Chains\n",
    "\n",
    "A number of methods exist that are much more efficient than generic Monte Carlo integration. \n",
    "The most popular group of techniques is known as Markov Chain Monte Carlo (MCMC) . \n",
    "\n",
    "MCMC returns a sample of points, or **chain**, from the k-dimensional parameter space, with \n",
    "a distribution that is **asymptotically proportional** to $p(\\theta)$. The constant of \n",
    "proportionality is not important in the first class of problems listed above. In model \n",
    "comparison problems, the proportionality constant must be known and we will return to this \n",
    "point later.\n",
    "\n",
    "Given such a chain of length M, the integral can be estimated as\n",
    "$$ I = \\int g(\\theta) \\, p(\\theta) \\, d\\theta = \\frac{1}{M} \\sum_{j=1}^M g(\\theta_j).$$\n",
    "\n",
    "Again, here **the values of $\\theta$ are not sampled uniformly from the volume**;\n",
    "they are sampled **proportionally** to $p(\\theta)$! Note that there is no $p(\\theta_j)$ \n",
    "term next to $g(\\theta_j)$ because the proper weighting in the sum is taken care of \n",
    "by the sample itself! "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### What is a Markov process or chain?\n",
    "\n",
    "A Markov process is defined as a sequence of random variables where a parameter depends \n",
    "*only* on the preceding value.  Such processes are \"memoryless\".  \n",
    " \n",
    "Mathematically, we have\n",
    "$$p(\\theta_{i+1}|\\{\\theta_i\\}) = p(\\theta_{i+1}|\\,\\theta_i).$$\n",
    "\n",
    "\n",
    "Think of $\\theta$ as a vector in multidimensional space, and a realization of the chain represents a path through this space. \n",
    "\n",
    "For  equilibrium, or a stationary distribution of positions, it is necessary that the transition probability is symmetric:\n",
    "$$    p(\\theta_{i+1}|\\,\\theta_i) = p(\\theta_i |\\, \\theta_{i+1}). $$\n",
    "\n",
    "\n",
    "This is called the detailed balance or reversibility condition (i.e. the probability of a jump between two points does not depend on the direction of the jump)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The use of resulting chains to perform Monte Carlo integration is called *Markov Chain Monte Carlo* (MCMC).\n",
    "\n",
    "Given such a chain of length $M$, the integral can be estimated as\n",
    "$\n",
    "               I = {1 \\over M} \\sum_{j=1}^M g(\\theta_j).\n",
    "$\n",
    "As a simple example, to estimate the expectation value for $\\theta_1$ (i.e., $g(\\theta)=\\theta_1$),\n",
    "we simply take the mean value of all $\\theta_1$ in the chain.\n",
    "\n",
    "Given a Markov chain, quantitative description of the posterior pdf\n",
    "becomes a density estimation problem. To visualize the posterior pdf\n",
    "for parameter $\\theta_1$, marginalized over all other parameters, $\\theta_2, \\ldots, \\theta_k$,\n",
    "we can construct a histogram of all $\\theta_1$ values in the chain, and normalize its integral\n",
    "to 1. To get a MAP (maximum a posterori) estimate for $\\theta_1$, we  find the maximum of this marginalized pdf."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "In order for a Markov chain to reach a stationary distribution proportional to $p(\\theta)$,\n",
    "the probability of arriving at a point $\\theta_{i+1}$ must be proportional to $p(\\theta_{i+1})$,\n",
    "$$ p(\\theta_{i+1}) = \\int  T(\\theta_{i+1}|\\theta_i)  \\,   p(\\theta_i) \\,    d \\theta_i, $$\n",
    "where the transition probability $T(\\theta_{i+1}|\\theta_i)$ is called the jump kernel or\n",
    "transition kernel (and it is assumed that we know how to compute $p(\\theta_i)$).\n",
    "\n",
    "This requirement will be satisfied when the transition probability satisfies the detailed\n",
    "balance condition\n",
    "$$ T(\\theta_{i+1}|\\theta_i)  \\,  p(\\theta_i) = T(\\theta_i|\\theta_{i+1})  \\,  p(\\theta_{i+1}). $$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "##  Markov Chain Monte Carlo\n",
    "\n",
    "How do we build the chain? Need to choose how far we will be allowed to move at each step. Need to decide whether we will accept the move (e.g. of $p(\\theta_{i+1})$ is smaller)\n",
    "\n",
    "*Various MCMC algorithms differ in their choice of transition kernel*\n",
    "\n",
    "The most classic approach is the **The Metropolis-Hastings algorithm** \n",
    "\n",
    "**The Metropolis-Hastings algorithm** adopts acceptance probability\n",
    "$$ p_{\\rm acc}(\\theta_i,\\theta_{i+1}) = { p(\\theta_{i+1}) \\over p(\\theta_i) }, $$\n",
    "where the proposed point $\\theta_{i+1}$ is drawn from an *arbitrary* symmetric density distribution $K(\\theta_{i+1}\\,|\\,\\theta_i)$. A Gaussian distribution centered on \n",
    "$\\theta_i$ is often used for $K(\\theta_{i+1}|\\theta_i)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "** When $\\theta_{i+1}$ is rejected, $\\theta_i$ is added to the chain instead. **\n",
    "\n",
    "The original Metropolis algorithm is based on a symmetric proposal distribution,\n",
    "$K(\\theta_{i+1}|\\theta_i) =  K(\\theta_i|\\theta_{i+1})$, which then cancels out from\n",
    "the acceptance probability. **In this case, $\\theta_{i+1}$ is always accepted if\n",
    "$p(\\theta_{i+1}) > p(\\theta_i)$, and if not, then it is accepted with a probability\n",
    "$p(\\theta_{i+1})/p(\\theta_i)$.**\n",
    "\n",
    "**The key aspect** is that\n",
    "\n",
    "$\\frac{p(M,\\theta_{i+1} \\,|\\,D,I)}{p(M,\\theta_{i} \\,|\\,D,I)} = \\frac{\\frac{p(D\\,|\\,M,\\theta_{i+1},I)\\,p(M,\\theta_{i+1},|\\,I)}{p(D\\,|\\,I)}}{\\frac{p(D\\,|\\,M,\\theta_i,I)\\,p(M,\\theta_i\\,|\\,I)}{p(D\\,|\\,I)}} = \\frac{p(D\\,|\\,M,\\theta_{i+1},I)\\,p(M,\\theta_{i+1},|\\,I)}{p(D\\,|\\,M,\\theta_i,I)\\,p(M,\\theta_i\\,|\\,I)}\n",
    "$\n",
    "\n",
    "so there is no dependence on the evidence.\n",
    "\n",
    "This algorithm guarantees that the chain will reach an equilibrium, or stationary, distribution, and it will approximate a sample drawn from $p(\\theta)$! "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "##  Markov Chain Monte Carlo\n",
    "\n",
    "**In summary, the Metropolis-Hastings algorithm consists of these steps:**\n",
    "\n",
    "1) given $\\theta_i$ and $K(\\theta_{i+1}|\\theta_i)$, draw a proposed value for $\\theta_{i+1}.$ \n",
    "\n",
    "2) compute acceptance probability $p_{\\rm acc}(\\theta_i,\\theta_{i+1})$.\n",
    "\n",
    "3) draw a random number between 0 and 1 from a uniform distribution; if it is smaller than\n",
    "   $p_{\\rm acc}(\\theta_i,\\theta_{i+1})$, then accept $\\theta_{i+1}$.\n",
    "   \n",
    "4) if $\\theta_{i+1}$ is accepted added it to the chain, if not, add $\\theta_{i}$ to the chain.\n",
    "\n",
    "5) use the chain (of $\\theta$ values) for inference; e.g. a histogram of $\\theta$ is\n",
    "  an estimator of the posterior pdf for $\\theta$, $p(\\theta)$, and the expectation value for \n",
    "  $\\theta$ can be computed from \n",
    "  $$ I = \\int g(\\theta) \\, p(\\theta) \\, d\\theta = \\frac{1}{M} \\sum_{j=1}^M \\theta_j.$$\n",
    "\n",
    "where M is the number of elements in the chain (e.g. \n",
    "the expectation value for $\\theta$ is simply the mean value of chain elements). \n",
    "  \n",
    "  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Caveats\n",
    "\n",
    "Although $K(\\theta_{i+1}|\\theta_i)$ satisfies a Markov chain requirement that it\n",
    "must be a function of only the current position $\\theta_i$, it takes a number\n",
    "of steps to reach a stationary distribution from an initial arbitrary position $\\theta_0$.\n",
    "**These early steps are called the \"burn-in\" and need to be discarded in analysis.**\n",
    "There is no general theory for finding transition from the burn-in phase to\n",
    "the stationary phase; several methods are used in practice. Gelman and Rubin\n",
    "proposed to generate a number of chains and then compare the ratio of\n",
    "the variance between the chains to the mean variance within the chains (this\n",
    "ratio is known as the $R$ statistic). For stationary chains, this ratio will\n",
    "be close to 1.  \n",
    "\n",
    "When the posterior pdf is multimodal, the simple Metropolis--Hastings algorithm can\n",
    "become stuck in a local mode and not find the globally best mode within a reasonable\n",
    "running time. There are a number of better algorithms, such as Gibbs sampling, parallel\n",
    "tempering, various genetic algorithms, and nested sampling."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### How do we choose the next step\n",
    "\n",
    "How far should we step (small steps in parameter space or large). This impacts the efficiency of the process but not if we will reach equilibrium. We want our samples to be independent of one another which we can determine with the autocorrelation function - we can then adjust the step size to account for this (e.g. ensure that we reject 50% of the proposed steps)\n",
    "\n",
    "Other techniques: **NUTS (No-U-Turn Sampler)**, Metropolis, Slice, HamiltonianMC, and BinaryMetropolis\n",
    "\n",
    "[NUTS](http://arxiv.org/abs/1507.08050) is generally the sampler of choice for continuous parameters (it uses  the gradient of the log posterior-density). See the end of this notebook for animations of different strategies \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<iframe src=\"https://player.vimeo.com/video/19274900\" width=\"640\" height=\"480\" frameborder=\"0\" webkitallowfullscreen mozallowfullscreen allowfullscreen></iframe><p><a href=\"https://vimeo.com/19274900\">Metropolis in the Square</a> from <a href=\"https://vimeo.com/user3812935\">Abraham Flaxman</a> on <a href=\"https://vimeo.com\">Vimeo</a>.</p>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import HTML\n",
    "HTML('<iframe src=\"https://player.vimeo.com/video/19274900\" width=\"640\" height=\"480\" frameborder=\"0\" webkitallowfullscreen mozallowfullscreen allowfullscreen></iframe><p><a href=\"https://vimeo.com/19274900\">Metropolis in the Square</a> from <a href=\"https://vimeo.com/user3812935\">Abraham Flaxman</a> on <a href=\"https://vimeo.com\">Vimeo</a>.</p>')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<iframe src=\"https://player.vimeo.com/video/19274173\" width=\"640\" height=\"480\" frameborder=\"0\" webkitallowfullscreen mozallowfullscreen allowfullscreen></iframe><p><a href=\"https://vimeo.com/19274173\">Metropolis in Diagonal Region</a> from <a href=\"https://vimeo.com/user3812935\">Abraham Flaxman</a> on <a href=\"https://vimeo.com\">Vimeo</a>.</p>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "HTML('<iframe src=\"https://player.vimeo.com/video/19274173\" width=\"640\" height=\"480\" frameborder=\"0\" webkitallowfullscreen mozallowfullscreen allowfullscreen></iframe><p><a href=\"https://vimeo.com/19274173\">Metropolis in Diagonal Region</a> from <a href=\"https://vimeo.com/user3812935\">Abraham Flaxman</a> on <a href=\"https://vimeo.com\">Vimeo</a>.</p>')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}