{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Model comparison\n",
    "\n",
    "To demonstrate the use of model comparison criteria in PyMC3, we implement the **8 schools** example from Section 5.5 of Gelman et al (2003), which attempts to infer the effects of coaching on SAT scores of students from 8 schools. Below, we fit a **pooled model**, which assumes a single fixed effect across all schools, and a **hierarchical model** that allows for a random effect that partially pools the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import pymc3 as pm\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set_context('notebook')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data include the observed treatment effects and associated standard deviations in the 8 schools."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "J = 8\n",
    "y = np.array([28,  8, -3,  7, -1,  1, 18, 12])\n",
    "sigma = np.array([15, 10, 16, 11,  9, 11, 10, 18])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pooled model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Average ELBO = -42.09: 100%|██████████| 200000/200000 [00:17<00:00, 11675.35it/s] , 9446.54it/s]\n",
      "100%|██████████| 2000/2000 [00:01<00:00, 1747.24it/s]\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as pooled:\n",
    "    mu = pm.Normal('mu', 0, sd=1e6)\n",
    "    \n",
    "    obs = pm.Normal('obs', mu, sd=sigma, observed=y)\n",
    "    \n",
    "    trace_p = pm.sample(2000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA04AAACCCAYAAABxYpjTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXl8VPW5/z9zZt+zTPaNJCQQdpHFilKU4oJSFyggt1iX\nWuut9upP+7K2Fq1V0V+19mp/tde2ehVcsCAudWmL4lIVBDRAICyBJGTfZpLZ13N+f0xmMjOZ5cyW\nScjzfr14MZk5y3O+Z3ue77MJOI7jQBAEQRAEQRAEQUSEybQABEEQBEEQBEEQ4x0ynAiCIAiCIAiC\nIGJAhhNBEARBEARBEEQMyHAiCIIgCIIgCIKIARlOBEEQBEEQBEEQMSDDiSAIgiAIgiAIIgZkOBEE\nQRAEQRAEQcSADCeCIAiCIAiCIIgYkOFEEOOEvXv3Yt26dbjjjjtw2WWX4ZprrsFHH32EG2+8EcuW\nLcOjjz6KvXv34sorrwxaJ/BvgiAIgkgUeg8RRHREmRaAIIgRDh8+jO3bt2PGjBn44Q9/iOeeew4v\nvfQSzGYzli5ditmzZ2daRIIgCOIsht5DBBEZMpwIYhxRWlqKGTNmAADKy8uhVqshkUiQk5MDpVKJ\noaGhDEtIEARBnM3Qe4ggIkOhegQxjpBIJEF/i0TBcxu1tbXgOM7/t8vlGhO5CIIgiMkBvYcIIjJk\nOBHEBEKj0aCzsxMDAwPgOA67du3KtEgEQRDEJILeQ8RkhkL1CGICwTAM1q9fj9WrVyMvLw/Lli3L\ntEgEQRDEJILeQ8RkRsAF+lsJgiAIgiAIgiCIUVCoHkEQBEEQBEEQRAzIcCIIgiAIgiAIgogBGU4E\nQRAEQRAEQRAxSJvhxLIsNm3ahHXr1mHjxo1obW0dtYxer8ell14Kh8MBALDb7bjjjjuwYcMG3HLL\nLdDr9ekSjyAIgiAIgiAIgjdpq6q3a9cuOJ1ObNu2DfX19Xjsscfw7LPP+n//7LPP8OSTT6Kvr8//\n3auvvora2lrccccdePfdd/HHP/4R999/f9T99PWZ0nUIBEEQxBiTl6fOtAhxk4r3UHa2AgaDNQXS\npJ+JIivJmXomiqwTRU5g4sg6UeQEkpc12nsobR6nAwcO4MILLwQAzJs3Dw0NDcE7Zhi88MILyMrK\nCrvO0qVL8eWXX6ZLPIIgCIIYN4hEwkyLwJuJIivJmXomiqwTRU5g4sg6UeQE0itr2jxOZrMZKpXK\n/7dQKITb7fZ3oF6yZEnYddRqr5WnVCphMpE3iSAIgiAIgiCIzJM2w0mlUsFisfj/ZlnWbzTxWcdi\nsUCj0aRLPILICCzLYcBoh8Xugs3hgcPlCfpdAEDICMAwAggZAYQME/DZ+71AAAgEI/8zAkCjkEAi\nnjizQQRBEET8cByHvkEbslRSeuYTRAZIm+E0f/587N69GytXrkR9fT1qa2t5rfPJJ59gzpw5+PTT\nT3HuueemSzyCGBM8LIvDp/WoP9mPli4juvRWuNxsWvalkouRlyVDdYkW50zVYVp5NhhGkJZ9EQRB\nEGOP3ujA6S4j5BIR5k7VZVocgph0pM1wWrFiBT7//HOsX78eHMfh0UcfxQsvvIDy8nIsX7487DrX\nXXcd7r33Xlx33XUQi8V48skn0yUeQaQVt4fF7q878P7eVgyanQAAiYhBca4SRToFNAoJZBIhpGIh\nBIIR44bjOLAcB4+Hg4cN+cxy8LAsuOHlOA7gOK9xNmRxwmByoK3XjOYuE3btb0d+lhxXnj8F588u\nBCMgA4ogCGKi43R7oxRsTneGJSGIyUnaDCeGYfDQQw8FfVddXT1quY8++sj/WS6X4+mnn06XSAQx\nJjR3GfHnd46iW2+FTCLExfNLcN7MQlQVadLuAXK5PWhqH8Keoz348kg3nn+vEZ8c7MAPr5yBgmxF\nWvdNEARBEARxNpM2w4kgJiOfHezES/84Dg/L4aL5Jbj6gkqoFZIx279YJETdlBzUTcnBVRdUYttH\nTdh3rBe/+d/9+PHVMzGrMnfMZCEIgiAIgjibSFs5coKYbLy3pxUvvH8McqkId6+fh42XTBtToymU\nHI0Mt109CzdfUQenm8V//+0Q9h/rzZg8BEEQBEEQExkynAgiBXz0dTu2f3wKuRopfrHxXMyckpNp\nkfwsmV2Eu9fNhVjE4Nm3GlB/sj/TIhEEQRAEQUw4yHAiiCRpaB7Ay/86AY1CjHvWn4PCnPGXSzSt\nPBt3r5sHsYjBn95qwOlOY6ZFIghiAnC604gDx/vAcVymRSEIgsg4ZDgRRBL0D9rwP28dgZAR4I41\nc1AwDo0mH9UlWtx21Sy43CyefbMBZpsr0yIRBDHO6R20wuXxIBm7yeZwo1tvTZ1QBEEQGYIMJ4JI\nEA/L4rm/H4XF7sb3L5mG6mJtpkWKydypOnz3gkoMGO14/t1GmkUmJiTt7e34+OOP4fF40NbWlvB2\nXC4Xfvazn2HDhg1Ys2YNPvzwQ7S2tuK6667Dhg0b8MADD4Bl09N3bTJx+PQAWrqNGLI4My0KQUwK\n3B4WNgeVrE8HZDgRRIK8t+cMmtqHsHB6Pi6cU5RpcXiz6vwpqKvIRn1TP/65L3GlkyAywXvvvYfb\nbrsNDz/8MAYHB7F+/Xq89dZbCW3r7bffRlZWFl555RX85S9/wW9+8xts3rwZd955J1555RVwHIcP\nP/wwxUcwMeGQ+CQLOzxB4/GQEUoQY8HBpn4cPNUPN91zKYcMJ4JIgM5+C975vBlalQTXXzYtqInt\neIdhBPjRqhnQKCXY/vEpnOkxZVokguDNn//8Z7z66qtQqVTIzc3Fzp078dxzzyW0rcsuuwz/9V//\nBcDbVFooFOLIkSNYtGgRAGDp0qX44osvUib7RIac0wQxcXANG0wUVZJ6qI8TQcQJx3F46YNjcHs4\nbLxkGpQycaZFihutSoofXlGH371+EM+/24j7f7AAIiHNoxDjH4ZhoFKp/H/n5+eDYRK7dpVKJQDA\nbDbjpz/9Ke688048/vjj/okQpVIJkyn2xEJ2tgIikTAhGQLJy1MnvY1Uo1EPAQB0OjXEopFxjkdW\n3zZyc1XQZclTK2AMxuOYhoOvnA4OMFjdca2Tavjut7PfDLPVhdry7DRLFJ6Jcu6B1Mvqu+fy8jRB\n922yTOYx9UGGE0HEyf7jfTjRPoRzanSYX5uXaXESZlZVLi6YXYR/H+7C+3tasWpJZaZFIoiY1NTU\nYOvWrXC73WhsbMQrr7yC6dOnJ7y9rq4u/OQnP8GGDRuwatUq/Pa3v/X/ZrFYoNFoYm7DYEi+8EFe\nnhp9fePP+2s02QAAfX0mvwIWr6y+bQwMmMG5xi7vYryOaSjxyKnXW4LOyVgTj6wHjnYDAORCQCYZ\nW3Vzopx7ID2yhrtvk2UyjWk0o4ummAkiDlxuFn/b3QQhI8Dai6dmWpykWb98KrQqCd7+vAUdfeZM\ni0MQMdm0aRN6enoglUrxi1/8AiqVCg888EBC2+rv78dNN92En/3sZ1izZg0AYMaMGdi7dy8A4NNP\nP8WCBQtSJns82BzucZbcnXzIT7SQ5l6DFQaTI+l9pBuzzYWOfkumxZhQnM3RYizHoaXbCKudqtSm\nE7eHxZEW/bgoMEOGE0HEwYcH2tE/ZMfyc0tRkD1+S4/zRSET4weXToeH5fD8e8fAsmfxG444K1Ao\nFLj77ruxY8cO7Ny5E/fee29Q6F48/OlPf4LRaMQf//hHbNy4ERs3bsSdd96JZ555BuvWrYPL5cKl\nl16a4iPgx8FT3uTu8UKqlN8hi3PUJI3LzeJ0lxHH2wyp2UkaaWgeQFuvCVb7eDJqE4flOJhtrrC5\nMAaTA3Zn+OO0O92882fOZsOpf8iObr0VDc36TItyVtNjsMFkdaKxNfPjTKF6BMETk9WJd75ogVIm\nwpXnT8m0OCljXo0Oi2cUYO/RHnz4dTtWLCjLtEgEEZHp06eP8lzk5eXh008/jXtb999/P+6///5R\n32/dujVh+VINx3nr2TEZLkCTKuXXp/jkZ8shHs4Lm4gTNp4YZeo5joPR6oJaLgbD8Dt3VrsbDpcH\n2WppKkTkRXuvGZ0DFlQXa5EXkH/mdHn8hux5MwqD1jFZnTjSokdelnxCtOEIh83hRku3CZVFat5h\nhBzHweHyBC3vu3bZCDeI3emGRCTkfQ2knol3b4UlxgPI5nBDJhGOSaEuXh6nW265Be+//z5cLnJF\nEpOXtz9vgc3hxqollVDJJ15BiGhct7wGSpkIb3xyGgND9kyLQxAROXbsGBobG9HY2IhDhw7hd7/7\nHS6//PJMi5U26k/2Y19jb6bFSKoceTgC7Y5ISud4pjVGNdK+ITsaW/Vo7jLy3uah0/043mYYNR4O\nlydtHq5egzcXxhIQatZrsEaV2xdC2jdo47mX8Xd+m7uMGLI40NLFPw+mtceE+qZ+DJr5hZS63Czq\nm/px+PRAomKGheO4uMuMuz3shJyg8BFO8hNtgzh0yns+Dp7qR0v32ORf8TKcfvSjH+Gzzz7DpZde\nil//+tc4dOhQuuUiiHFF14AFH3/TgfxsOS6eX5JpcVKORinBuotr4HB5sOWfx6mEKTEhEIvFuPzy\ny7Fnz55Mi5I2HG5Pyo2WREj1IyHwmCbi48Zsiz6RbBn+na+SHUTIeHxzsg+HTqcnbNM9bMGKAipT\nnu4ywhBF7nhn9cfj6fXJFM+7rlvvLQLDt4WHy+0dW1uEcMdEOdKsx/7jvTG9nsDIvbX/eC/2H0/9\nBIzT5cGeo93oSUGBnGj4jkOAkWtPb7LD6nDDOJz3lG4ZfPDyTy5cuBALFy6E3W7HBx98gJ/+9KdQ\nqVRYs2YNNmzYAIlEkm45CSKj/G33KXhYDt9bNvWsLdu9ZHYhvjzSjUOnBrDvWC8W1RVkWiSCGMWb\nb77p/8xxHE6ePAmx+OzyAFvGYaJ5Y6shpVVEA/XVYCOKm1B98dKB3mSHTju6bDubxrERxBFKFiiC\nh2UhTLAdQCAcx8HDchPi/WrlWbQl3lPl9rCjjt/tYdHSZUKxTgmFzKuym4efD24Ph3iGKx2eXV9B\nl+YuY9ryvj0sC6fbAyD8mI51GDPvId+7dy8eeughPPXUU7jwwgvxy1/+EgMDA7jtttvSKR9BZJxj\nrQbUN/WjtlSL+bW6TIuTNgQCAa6/bBrEIgav/OtEzBlVgsgEe/fu9f/76quvAABPPfVUhqVKLQ6n\nJ9MijMKnuCRDpNn9ICNqPLonxpgeQ/gQuHSGWvnsJj77CJz1N1pivyf4nNNDpwaw/3hv2Gukd9CW\ndIVJjuPGdUho36AN+4/3on8o+Nx3DVjRb7Th+JnECqfEe8RjEW3i9rBxTQ59c6I/eljoGM+z8PI4\nXXTRRSgtLcXq1auxadMmyGQyAMDixYuxevXqtApIEJmE5Ti89tFJAMC65TVn/UxoQbYCV11Qie0f\nn8Lru5tw08q6TItEEEFs3rw50yKknVQ2rAzE5fbAw3Jj3lMnHIEKGhvymRlrTSiNuDxsyrxoHk/6\nldpoCm2vwQq5Mj2FK3zhbByC9WCr3YXTnUMQQIDFMxKPgjh8Wg+rw+UvdOHbx1iZUkMWJ8w2F0p0\nyrC/+3LN+gZHvI02hxs9w+GB7jAGbarvEreHxf7jvSjMUaCiQM3vmk1AiIbTethdbpxTkwepOHbj\ncHdASKIwjGd0rJ8WvJ6eL774IpRKJXJzc2G329Ha2oqKigowDIOdO3emW0aCyBhfNnTjTI8Z580s\nQGVR7EaYZwOXLCzD3qM9+PehLnxrZiHqKjLT9Z0gArn44oujvsg//PDDMZRmfDAwZAcHLmxYVzgO\nnOgDMLpKWipo6zVj0OzArMqciOcpcDLb7vRAIROP+v5szK/sHbTFF8YUYQjGwmMS6dzZHG6c7jJi\nyO5BjiK+0NhkzqnD5VWak83zszqSi6Dw5Tcliq+aZEG2PHw44vCwD1kc/pC9wHYECRsHcQybfdjT\n3a23oltvxfyaPEgCDBtPijye9uEm2C43y8twCsRXiTOQsa55wctw+vjjj7Fz507s3LkTAwMD+PGP\nf4wbbrgB69atS7d8BJExHC4P3vj0NMQiBquXVmdanDFDJGRww+XT8fBL+/HiB8fw65sWxf1wI4hU\ns2XLlkyLMGbw0QM4jsPJjkEA4G048cXl9uBE+xAK41D2O/q9vZlYjoMwkuEUcGSuwKpgQR6nOIVN\nAR39FhiMdsyMYvSFopBGNx4CN2OyulAQx/wTBw6t3SbosmRQykb2k0mb0mf8uNxsyouVBBpWHMfB\nYHZCq5KAEQhiFkBgWQ4NzQMozFEgP429FVu6+VdHjEakcxjkZXO4oVGkpnZAR78FVcUjk74cx6F3\n0IYctXSUERLqzDFanf5ni9Plwb/rOyDkWGiUEhRkK8AwgqS8PYk4YcOtM9bVAnnFA7z++ut4+eWX\nAQAlJSV44403xlWfC4JIB//YewYGkwOXLCxDrlaWaXHGlMoiDS5ZWIZegw1/292UaXEIAiUlJSgp\nKUFeXh6OHj2Kffv2Yd++fdizZw+2b9+eafHSjtMVnGN0rHUk5yHe0sSBmG2uUfkDBrMTJqsTTR1D\nQd+73B50DVjQ3DmE9j5zUvsN1L0D9Z5klSCr3YX2XnNcXo62XhPMdpe/Clo01HKvQhtPOGW8Hhez\nzYUuvSVMGevkFcSm9qGoVeF4KbNxihFr8Y5+i//zwJAdx9sMON3hNVRieTlMViesw96wUNgweU1D\nw9UC0x2qd6bHFFcDa0GEz1HX4bFg76A1KGey12BDc5cRJ9uD722HywOD2Rn0XeDQWYbL4RvMDrT2\nmNA1YEGypCrEbqx7ZPHyOLlcrqDKeWdbBSOCCKV/yIZ397RCq5Rg5XkVmRYnI1y7tAoNp/X46OsO\nzJ2qw+yq3EyLRBC4/fbbYbPZcObMGSxYsAD79u3DvHnzMi1W0nAcB5eb9SrkYbS5br0V5QVqsCwH\nhhFgyOocvVCYbTpdLKSSyB7jhmavcp6tlvpDiEaUymBBTrYPwWh1QmN1w2iywWJzYVp5sCslUNky\nWZ0RjZHALYd6HPjCchyau4zIy5L7Z+gPDRsbaqUE+by35CWenI4hC/8y46nyFKViYr3f6DWSywvU\nyW8shBNtg7DYXJg7VReszMaQW28cGUvTcFEivckOQBvbkI5yzhpbDbA73DgnoBpk4xlDUKiqy8Xi\nRNsgyvJVkEtTl/vXmaRhkcriTIH3ceewkRpagObQqYFR3r1o96UvhDIht1EShLuXRMJxWFXvO9/5\nDn7wgx9g69at2Lp1K2666SZcfPHF6ZaNIDLGto+a4HKzWLOsOqUP04mEWCTELatmQMgI8Py7jTDx\nUNQIIt00NzfjpZdewooVK/DDH/4Qf/vb39Dbm/kGscnS2KLHgRO9cLo8YfVMjvMm5391rGdU5a1I\nnOow4pumPl5KWKBiFEkNCa1sFqsp65EWPU60DwbsJHh/LMuhR28NMq7iMQ4GTQ70DdpwtEUf9reR\nbXIwJvD8cqWgkiCQurytvkFbxG3ZHG6c7jRG9AIeadGjNaRBqN44utl5uHPfO2iD0TpyDYWTgGU5\n6E12ONwe9EaogMayXFj5A72pvp9ZjoOHZWNeY9FUZpPVCZeHjVpUw+5yQ2+yh/VYjSkhBohvQiMa\nfC+rwAkBh7+sd/D+woVERtt8IuGavQZryH0Y3eBhOQ4NYZoHsxyHUwHe8FjNqFMNL8PpZz/7GTZu\n3Ijm5ma0tbXh+uuvx1133ZVu2QgiIxxt0ePA8T5MLdHiW7NSn0Q9kSgvUOPapVUYsjjxwnvHzsrE\nbWJikZubC4FAgMrKShw/fhwFBQVwOie+Ud83XFXL6WbDakQcOHQNeBPUQ0PoIt2WPu+ChYfhFGiw\nRPK8xNOvh0/56K4BC5q7jWgOyB+JVACB5Ti09ZrhGFayWZYLNspCt60fmfFvah/C0Ra9P0yLD539\nFhw40RdkXPTorQlNIKUqBaPHYEVn/+giBS3dRhw81Y/eQSs6+oI9HV0DFgyaHTBZnUFjojfag8Yv\n2qP9dOdQUI5PqLfQYncFnTd3iJeRG17uYFM/jraMLqsddM4DPpptbv/5jnjtBVyqQ2YHDhzvQ3uv\nOWiRxtbYpbxdLtZvdMZThCP0TmnvNeN0Z/xGWOB2Uv2WTdQfE20Yov3Gshw6+i1BBjHLcTjdZQyZ\n5Ih+pGaby9+zKnCdQZMDfREmjwaG7EmXro8F76n06upq6HQ6v+K0b98+LFy4MG2CEUQmcHtYvLLr\nJAQA/mNF7Zg3VhuPXLqoHA3NetQ39eODr87g8sWTM3SRGB/U1NTgN7/5Da677jrcc8896O3thct1\n9vQci6S0cVwSsfwhq4VrNprQpEgUcY6F6TsTuAeWA+yu0R6dSGJ09lnQ0W+G0eLEzMoc9IfxlkTC\nG/YFvxIemZGd+8pD600OsCyH1h5TcEELeA2pgpzwBQlC0+aPNOshEAAzpuTwljscLs/oYwis+BaY\nE+RysxFn40PHIra/MXDhkaW7Bqww2Zwozh0psz3qGuY42BweONwev8cjEoHr+irRAYBULPSHqUYQ\nBXqTAy6PB+39ZpwTsEy4inqh59LucvtLcXfrrTi3Ni9sBbdYtPebYy8UhlB5ksFgCp4gCOvB5mGe\nRX0mcN5zZQvjEezSW9HWa8Kg2YGZw9e72Tr6HESTwO1hw3qSvXJFXs9XMKe8NH3VgHkZTr/+9a+x\ne/dulJWV+b8TCAR46aWX0iYYQWSC975sRWe/BcvmFaOiMPUx4BMRhhHg1u/OxIMvfIXtH59CVZFm\nVF4DQYwVDz74IL755htMnToVd9xxB7788ks8+eSTmRYrZXg8XNi0gR6D1V+YYDTxGT0n2gYxaHbi\n3GkjuR+BkTp8w3DC1dTyKTUxCy1E8ixFcM9Yh2eR/aFoCRh68RiePgWeEQjQ1DkUdpnmbiM8LIfi\nCL15AjHZvJ6qjj4zSvJUvOWIN1ywd9CKKYVqMIzAbzCmmsCR9x1X/9DIvkJPjcnmQntfcDGTsCW5\nEVlZtzpcqG/qx/yAfCUg2LiKh0geLJ8RarK6kKPh02MoNf6hoP5ZYTbpZkf3A4u05+NtIZMWERaM\ndh6A6J5SDhzO9JjClmn3eZrsjmHvMMfhaJjzFO0WPnwqfKgix415WtUoeBlOn3/+OT744AN/41uC\nOBtp6TbinS9akK2WYs2yyVN+nA8apQS3XT0L//eVb/DsW0fwwA0Lka1OTyNEgojGHXfcge9+97tw\nOp1Yvnw5li9fnmmRkiZwhjjcrLqPRBWG0NUMwyFrgcZNkLGUpohcLsLnoGUiaFM+T4pP0eMjYmhu\nV6ChZ3O4YbW7IQxILPftun/IBqc/FyT6Ps70mqCUiaBVSWFzuCGVCKNGKrQNG04cx+FMjxlqhRh5\neZEn6RqaIxsG35zsQ01p1qjvewxWFOUqRxUAiEo8hmjAooxAAHa4sImP/iFb0MRjqLeBZTnYPW4I\nGQZiEQOJWAi7rwFuFDGcbg96DVbotPKw9wjfEDuO46BWiP1GXzj4bqutN/n8mqa2yCGngXTrrWEb\nwPrgOC5smFq4iRCHy+P3sE0pjNCjMsoYcFzk3lYjkw4jcsVLNM/kWPQziwavgOWysjLKbSDOalxu\nD/7690Z4WA43XVHnb8xIjFBTmoXvXTQVRosTz+w4FN9LmSBSxNq1a7Fr1y585zvfwS9/+Uvs3bs3\n0yIlTWDejCfKuzaSysRx3m3YHG70GqyjQnUiEbirQE8P77e9wJs/FJxDwld5Df99pJClUGWMD23d\n4ZXa3kEbDp7qx8mOwbAhhYE5ZHzCtR0uD8w2Fw6e6kdTe3jvVCAHjvehs9+CLr0lap6Wb9vRfmsJ\nU9TA6WKhN9rjquzW1pdYiJlvfAKV89BzGzqELMehvqkfB054i7roAtp9xNI1T3cZg8qXB8JXTeWz\nWLRqfpGuiUQV+o6QsbeMyuvx0tpjCi5iMbw7D8tiz9Fu7G3s8VeVDCSaWNEa+wY3pg7+LZo309fe\nwOcdi7T/ROwKm9ONQVNmc1p5eZy0Wi2uuOIKnHPOOUFlyTdv3pw2wQhiLNn5aTM6+i24aH6JPyaX\nGM2KBaVo7zXj34e78Oe/H8V/XjOL8sCIMWXZsmVYtmwZ7HY7Pv74Yzz++OMwGAzYvXt3pkVLmMB7\niGU5IEKEUDQ140hIPsCcKh2PPY9skeW8JYp7B22QiPkVgRBgpADFqC1HKg3o/xj+aJq7jNAqJZBJ\nvOpJr8GKvkG73+Pk9zSErN7Rb0FJSMicSBT8bBIIhivsmcMrXuEk4vN46xqwoijXm+vEJzzO5fEE\nGSoGoz1mo9dIjE6eB4RCQUyDLDKxlVlLQF4LwwiAENHdLBtkfIYWG4lmlPAZhUjJ/7x7gHGxjazA\n30MrMrIcF/Z6i6cHmdnmgkgoCBvSyrdKnM9YtTnSM4nJchx6DFYoEqwsHOveSdQdE+mZM1bwGo0L\nL7wQF154YbplIYiM8M3JPnzw1RnkZ8uxdtnUTIszrhEIBLj+smnoH7Lh6xN92P7xKay9iMaMGFua\nmprw7rvv4oMPPkBRURGuv/76pLZ38OBBPPHEE9iyZQtaW1vx85//HAKBADU1NXjggQfAxFFNLhEE\nAW6Ulm4jlHF6vMMpIM0BM9Onu4xBym649ViWw/G2wbCJ9D7srtjVqlIRnGK2ufyGU2iZ6EgV/9p6\nTaMUWe84jgikNzpwsn0o5Q0zbU53UokXh5r6IUphfGSikvBtaNw7OOKliHTYgSXzQ5dxRykPzqdq\nod/LFXKx8c3p4lUYIeDzyTChdG29pqCS9+HkibRlm8MdVG5co5bzWC8yoQ2sR+0xwUvLw3L+50he\nVvwyxuyLFiJXW68Zdqc7bPjpeIKX4XTNNdegvb0dTU1NuOCCC9DV1RVUKIIgJirtfWb85e9HIREx\n+M+rZ0VtFEl4EQkZ/OTa2XjkpQP4YO8ZaJUSXLqoPNNiEZOEVatWQSgU4qqrrsKLL76I/Px425wG\n8+c//xlvv/025HKvYrB582bceeedWLx4MTZt2oQPP/wQK1asSIXoEQnV4yOF6kTSfzzhFN6QbfYY\nRofkBIeCaYCUAAAgAElEQVTicEHlg/kQTi+KpqMFFhDoNdiQrQmfJxnNi80wAvQP2vzFIqLhCwnz\n7394ptoTh2eAb9hjIImEIKW0kAMPIy6cV2P/8V7Mqoyv0XkiEQejqinGOVwCHvMY/VGMiXhPT6Tw\n0dAcKT6XldXhzdNKJeHu7WASs5wCQw/7Bm0xDbzGFn3QpESsS+PYGQMWTM/3X0MdwxUJK1NYYTAd\n8JpGe++993DbbbfhkUcewdDQENavX4+33nor3bIRRFrpG7ThyW31sDk8uHFlXVo6qZ+tKGVi3LV2\nLrJUEmz7qAmfHerMtEjEJOGJJ57Am2++iRtvvDFpowkAysvL8cwzz/j/PnLkCBYtWgQAWLp0Kb74\n4ouk9xELsYifRyuSQh6uNw6fmfvA7Z1oH4Q77nCxKInqYZS1wJAnN8vC6YpfQWJZDk2dQzyURaB7\nwMrLS+YnzPDanPzWDxxLp5uNaPyOdzoj5A9FIqZXAaMNlcDwtI5+S9xqfd+gDS43Gz10NUyuj4+D\nTf0xvU4t3ca4+n4B/EL1Glv1GIijlH40+BqA0Tx80YjlyQplyOr0F54JJJKcLMdBPzyZkuo6CjGr\neiYBL4/Tn//8Z7z66qv4/ve/j9zcXOzcuRM33ngjrrrqqrQJRhDpZMjixJOv1WPI7MR1y2uweEZB\npkWacORlyXH3unl47OWv8b/vH4NCKg4qb0wQ6WDatGkp3d6ll16K9vZ2/9+BJX+VSiVMptj5BtnZ\nCogS6Pni36dIiF6jnlfIjkadWL5BOHJzVdAMJJYv4JNVIwmWR5erAgd+4UcalQRcmDBImULqrzSn\nUQcXW8jKksMdwWDLy1MHLW+xueIKg8rVqSCXikbtkw9ZWQpozF5jSa2RA91m/76zNNKwxxlIsuFa\ngcjkkoS35+YhS+DvMokQIkn0azL0nGVlKTBk8xqkQzY3ivOU0MSZp2OwuVE3JQcadeRms9GOQ6aQ\nQuOKrqx3GOyonpLLeyxzclXQ9MZnePpI5HxJ5BK091tirtuuj+4tCr1vohGPnKrhipEutweargiF\nWkxOODkBqku1/m3rdGpoOpOrVmi1u6JWq0wGXk9ghmGgUo30HcjPz097zDdBpAuDyYHfvV6P3kEb\nrjx/ClYspLDTRCnJU+HOtXPxxKv1+J+3G3Dn9+Ym3eCRIDJJ4LvNYrFAo4lQqjcAQ5KhNwPD+SBG\n09gmPff1mxPap0Ytj7jev/Y0894O5/GE9Yz1SxgoRAL06K2j9hNN3r4+U9Dv0eQMR3+/CVKxMKEx\n6ell/Os1nOiF2+n2hxMKWBZGS2TvRbxyxiKd11GorJHNlhGEHAtjgCeiX8wEbyMBeVmXGzqlOOy6\nUrEQUpkk6nb57vOT/Wd4eW8BoLfXmPL7KRr7GlJznjs6B3ntP145PS43DjZ2g0P08TaabHA6nP5l\nQu/jROA473YSJZrRxcv6qampwdatW+F2u9HY2Ihf/epXmD59esICEUSm6Oi34JEt+9HRZ8F3FpTi\nmgsrMy3ShKe6WIvbV88GADyz4zBORWgWSRATgRkzZvhLnH/66adYsGBB+neaoW4fTQlXXksNoQrp\n1BItAKBzwAKXm0VzNx+1fIRkw304LvFTEVjpTW+yw+UO3lK4ZsGThdDTMhTFiOQLy0Uu8hCthHu8\n8DWaAEzYFh2BeYCpxGJ34UyviVevq2hl0ROBbxPvROBlOG3atAk9PT2QSqX4xS9+AZVKhQceeCDq\nOizLYtOmTVi3bh02btyI1tbWoN9ff/11XHvttVi7dq2/jOzg4CAWL16MjRs3YuPGjXjxxRcTPCyC\nGM2xVgMe23oAeqMDq79dheuW1/CKzyZiM3NKDm797kw43R78/vWDo/pSEESq6OjowI033ohLLrkE\nvb29uP7664NC7ZLl3nvvxTPPPIN169bB5XLh0ksvTdm2I5GpLomRkt4zhSwg5MvX4ycexlO7SZdn\nYirR6SD0tPAp7BELm8M9qrFupmmaoJOGmW4oO5oUyJPGQ+IVqqdQKHD33Xfj7rvv5r3hXbt2wel0\nYtu2baivr8djjz2GZ599FgDQ19eHLVu2YMeOHXA4HNiwYQOWLFmCo0eP4sorr8SvfvWrxI6GIMLA\ncRz+8VUbtn98CgIBcPMVdVgyuyjTYp11nDstHzdcPh0vvHcMT2yrx33fPxf5CZQwJYhobNq0CTff\nfDOefPJJ5OXl4corr8S9996Ll19+OeFtlpaW4vXXXwcAVFZWYuvWrakSlxfjTW3JBOfW5ied0J1o\nL6QgEjwZ0XRPRpDeGfDxAiMQhFfC06CYc+DQGKZ5MUEA6TUGeXmcpk+fjrq6uqB/S5cujbrOgQMH\n/L2f5s2bh4aGBv9vhw4d8jfTVavVKC8vx7Fjx9DQ0IAjR47g+9//Pn7605+itzf+GSeCCMRqd+PZ\nNxvw+u4mqBVi/Oy6c8hoSiMXzinG+ounYsjsxJOvfZNQKV+CiIbBYMAFF1zgL+Kwdu1amM0T28PZ\nPZDaMJWJiFjEJNMKCUDi1cMCOd0ZX3igj8Eo4WeTIbJBIhJi7tTwTZdtEzSEjcgM35xMPnQwnU40\nXh6nY8eO+T+7XC7s2rUL9fX1Udcxm81BBSWEQiHcbjdEIhHMZjPU6pHEK6VSCbPZjKqqKsyaNQvn\nn38+3n77bTz88MN4+umn4z0mggAAnGwfxHNvH8WA0Y7aUi1+fPUsZKnC9w0hUscli8phdbjx9uct\n+N22etz7H/OhksfX0JMgIiGTydDd3e1XRvfv3w+JRJJhqZLD6nCNqk43magq8hbgSKQnUCDx9GeK\nhK/XUyqZBHYTakq1kIrDV5Z0uslwSjdqhSSufKzxTCq8RSzHpS2rMO4ntVgsxuWXX44//elPUZdT\nqVSwWEbKMrIsC5FIFPY3i8UCtVqNOXPm+JsQrlixgowmIiE8LIt3Pm/BO1+0AABWnT8Fq5ZMgUhI\nlSDHiqsuqITV7sauA+146vWDuGf9PMilk1cxJFLHz3/+c9x66604c+YMrrrqKgwNDeH3v/99psVK\nislcOKCiQI38bAWA5A2MeBv4hpKKXjIihhnVDytVfXvGM5PBq5YpJCJhTOMzE4UpakuzcCLDBWYi\nwXHRuswlBy9N5s033wwQhsPJkychFkefQZ4/fz52796NlStXor6+HrW1tf7f5syZg9///vdwOBxw\nOp04deoUamtrce+99+KSSy7BypUr8eWXX2LmzJkJHhYxWekbtOG5d47gVIcRuRopblk1E7VlWZkW\na9IhEAiw/js1sDnc+LyhG8/sOIS71s6FOIleNwQBeN8f27dvR0tLCzweD6qqqia8x2ky65wa5ci5\nS3YcklXiUuCwSqCJMDEREECQsTy1wNsinGEOAEJm9M0jZJjU5P1NQNQKMSym9BiTvAwnX2lWH9nZ\n2XjqqaeirrNixQp8/vnnWL9+PTiOw6OPPooXXngB5eXlWL58OTZu3IgNGzaA4zjcddddkEqluPvu\nu/GLX/wCr776KuRyOR5++OHEj4yYdHzZ0I0t/zwOu9ODRXX5uP7SaVDIKEQsUzACAW5YOR02pwdf\nn+jDs28ewX9eM4s8f0RC3HfffVF/37x58xhJknp8IWoF2QrkZcnR1mPCkNUJqViY0tLK0RALhRmp\nBCcK6gmZWQvSYhs/VdpK81Ron4DVSacUatASZxn5VJOjliH5un3BFOTIU14ymzcBt4UuSxZWDqFw\nbO8diUg4rr2MCpkYFlN6PL28DKdEXkgMw+Chhx4K+q66utr/ee3atVi7dm3Q72VlZdiyZUvc+yIm\nN1a7G1v/eRx7jvZAKhHi5ivqcP6swnF9U08WhAyDW787E/+9/SDqm/rxwnuNuPnKGUnnMhCTj0WL\nFmVahLTBBMwWq+RiTC3Ngt3phlohgdXuxtEWPW9Phkwigt0Zn9pYVaSBLkuOrxp74lovGlqlFIU5\nChxvi175TCoZ8UKLRYlNqpTlqdCWAiMj3r5R6UQuEUGnkacl5yoUhVQMqyM5o9EzXNqeCeP5CESr\nkGAojbk4IoaBQiqC0RH/JECuRhY2rFKnkaM8X43SPBVOdxqhT5NCHonAUN5IYb3h3qmpCD2NRG1Z\nVtJVMCcqvAyniy++OKwS6qtq9OGHH6ZcMILgw/EzBvz13Ub0D9lRVazBj1bN8MfLE+MDsYjB7dfO\nxpOv1ePLIz2QSUT4/iW1ZNgScXHNNdf4Pzc2NmLPnj0QCoVYsmRJ0KTcRKS2LAtGuwc6lddDLhYx\nEIu8IWwKmQjn1Oqw7xi/KrNSEQN7DL2UEQigVUpgMHsrweVoZCmdzCjIVqCySANzAh6cbJXULxdf\nstTSlBhO4wmBAGPmgJtVmYOvjiVnNCt5FgCaVpGdMgM9XO6Pm2WDQj5z1DLeho48pEBLqU6F9n4z\n8rLlYBgBGAgiFsBIKwHHI4owuRDtfaqUiWGxjx9vairIpEeWl+G0atUqiMVirF27FiKRCO+88w4O\nHz6Mu+66K93yEURYHC4PdnxyCrv2t0MgAK48fwq+SwUgxi0yiQh3rp2Lx1/+Bru/6YBcKsKaZRNb\n2SUyw/PPP4/XXnsNy5cvh8fjwW233YZbb70Vq1evzrRoCaOSi1FZnoO+PlPY34UMA7VCAovNFbPi\nlFIuhtHqipqPIRIyaZ24cA97H6LNeNeVZ0MSRglNZJI8GaMvnQpYUY4SbpZF32ByXqNk82sKshXo\nMUQOM0vFpeB798baVLK7ysuS+8ezRKcc5SVUKyRBB5TosZ1TkwepWIiCHEXCntBwyMQi2F3xeYQD\nD4ERRAirjXLjaJUS5GfJU+5R5Tu2fDzCQoaBAPxyBItzlSjNU6FrwJqRHC5ehtNnn32GN954w//3\nD37wA1x77bUoKSlJm2AEEYmmjiH89e9H0WOwoTBHgZuvrEN1sTbTYhExUMrEuHv9PDy29QDe29MK\nuVSIK741JdNiEROMbdu24Y033vC3u/jJT36C6667bkIbTnyYOSUHbg+L/ccje56qirXIUUshFjE4\n02OGkBHwUkRSbUP5dLhoxpk2QmuIRMyDWIZTtKpkMrEwbsOkokCN1p7wRm4gLMeBS4FeJxAk15dm\nSqE6huE0evwiNrONQUzDKYmLraJAjVyNbMQQDbMpsZBJiaPO51lKpdEEADqtDAqZaFQhk2g5hqFD\nFm5CIlcrixoCWZCjyFgoarjwzdB7kuM4eHheb6k+J/HCe+9ffPGF//Pu3buhVCrTIhBBRMLl9uBv\nu5uweesB9BpsuGRhGR68cSEZTRMIrVKCe9afgxyNFDs+OY2Pvm7PtEjEBEOr1fpbWwCAQqGYNO+j\naB51EcMgP0sOkZBBUa4SC6fnRwyfClVjUl0OXTact6SUJdCCIBELIYb4sypzwn4vYhjkaGRxG46y\ngJCu2tLIVVsdTg8k4sSUPEGQ1yS585PI+upEe++lMbyQYQTB4xJu94JgQ6MgJ77QfaVMDEWaW2dk\nqUdPGkS/ToKPNNSgrSrWjnmKApdkn6RkJgIyHebP6+p46KGHcO+996K/39vNt6qqCo8//nhaBSOI\nQFq6jfjL3xvR2W9BfpYcN11RR2XGJyi5WhnuWX8OHtt6AFv/eQJyiQjfmlWYabGICUJZWRnWrVuH\nK664AiKRCP/617+gUqnwhz/8AQBw++23Z1jCzBCqh0RL0heE/pZiPaQo12vICgQCKKQiWB3BoUk5\nalnEdSPlcEQjUPzp5dnQGx3oHYxdAa2mLCtmMQNg9Ox44CrR+tMJBN5QwM4BS8Rl+JDM6ZHwbAGx\nqK4AeqMdDheLtl4TsjXRPRiRSNYIDwzFC0XICIKN3DAKdOjxahTxtSqIZGSnDEH481mer0LjmZFC\nKnOqcnGyfQg2p3vUYYYaTvlZ8rC74uP5zRQ6rQxd+sTuC60ys+0neBlOs2bNwrvvvgu9Xg+pVDpp\nZveIzOP2sPj7Fy34+xetYDkOy+eXYs2y6qBKTMTEozBHgbvXn4PHX/4af323EVKJEPNr8zItFjEB\nqKysRGVlJZxOJ5xOJ5YsWZJpkcYFYfOJIkzrhvZ84ZsjVJyrjGgE+JLwdVp5UChNqAiL6gqi7m9K\nodpf2UwhFaG8QI1THUMozVehuSt2qJFAIEBVsSbIcIq0O3a4cVOs2e/Q8YqmiAYaihw32oBVysQo\n0Slj9pwSCLweu/4hb6+rRCu5TS/nN8HICATQab0KeI5aCrlUlPay4tPKskdVXYxmeDECgf/akYiE\nCC5kL0BRrgLFOmXCOWUSCf8S25H6KfEhdB9Zw+MdiEIm9l+34hBPc6JhlIkUXokEB4y6sSoK1Og1\n2GALqOpZV549auIE8PZZ6tLH3k9+lgLFOgXqm7yOm4XT8yFkMhuqx8tw6ujowP3334+Ojg68/PLL\nuO222/Doo4+itLQ03fIRk5iOfgv+8vejaO02IUcjxU0r6zBjSppng4gxoyxfhbvWzsUTr9XjT281\n4L/WzMXMdM/2EROeyepRSiWBhotCGj0ka2qJFk0dQ8PLRlYZKovUkIgZlOiCJ1ZHecJiKKZikdCv\n4AmFDLJUUpw7LR+2MMqXj8BNht/86C9FDAO1QjwsY3QlNPTXSIcwuyoXZpvLb+CFU25nV+XCYOKn\nvBbmKCAVC72G03H+htPcah0OnvIqmokUTIrmRYtF4NhUFWlwOoqxG2+FOpVcDIYR4JypeRCJBEHj\nKBEzKC9QAwDcnuBxL9Gp0NEfuThBRYEaHg+HPG1kT6iPRIp0qBUSmIa9d6GXTlGOEtPKszEwMFo+\nn10WeA5ZDpgxJQddAxYMGO1BHrZzavLQNWCJ2G+qpiwrpS0HQsnReMcvMP9PrZCENZz4kK2SoqpY\nE/Qdn0me2VW5Ce2PL7zuqE2bNuHmm2+GQqGATqfDlVdeiXvvvTetghGTF5bl8MHeM/j1C/vQ2m3C\nktmFeOimxWQ0nYVUl2hxx+rZAAR4eschHGnhMQVFTGpefPFFLFq0CHV1dairq8P06dNRV1eXabEm\nFDyi0wKWHVlYHEXRFYuEmFKogTg0NCyBXIaKQjW0CgmqikaUpugT7CMyxlKsRAyDsjwVFkzP529U\nhOxbIBBgWlk2Kos0fu+aiGGglImD5IzkFdCqJJBLRMjVyEIaAAfsY3g/ORrZKI9XLPgaPjOn5GDh\n9Py4th0Po64FHkQ7fT7vnVQiHK7CFrBwwFCrhvOzyvK8BWSyw+QUAd5mvQqpCHlZcpTmq9IW0hbt\nmqwoVEMmFYU9bt/1E+S15Dio5GLUlGZhwbR8zKvR+X+SioXIzxrJdfIZeb61+XqWZ1XyMDy40Uag\nOFy1ziSGNFwYbaxzNK0sG0pZgvl5POH11DAYDLjgggsAeIVeu3YtzOazq2cCMT4YNDvwf1/9Bq/v\nboJCKsQdq2fj5itmQJFIkjExIZgxJQe3XzsbHAc8vf0QGpoHMi0SMY558cUX8eabb6KxsRGNjY04\nduwYGhsbMy3WmDG1hH8xnGi2Bp/mmAIIoFFKMG+qDnUVOXEr8Ikik4hQNyUnxAAIllcmHvktUJcK\nJ2Lg7wum56NkWKHmS6iXQSDwKuMF2QqIhAxmVeZi7lSvAhuo7EWK5GIEAsydqkNNlMISwSW1Rz7H\nG6YUTc8UCZmUhz0FFs6IpacLBF7jbU6VjpdSPzpkMvxy2Wopzptd5D/PkZYrzFFgTrUu4TYmedrw\nuUWhBIkdQRixSAixUBjynVcuSWDoa8DvIiEzatxSYfup+BYGCdjXubX5YBjBqKa4vgmAREhkPaEw\n/c8oXleLTCZDd3e3/yD2798PiSSzyVnE2UdzlxEP/e8+nGgbxPzaPDz0w8U4p4byXiYDc6pzccdq\nn/F0GA2nyXgiwlNdXQ2dThd7wbMUnVaO82aMLqYSzg6aUqgOG17HwRvyA4xWRgFgTpUO580oxOIZ\nBRAJGcgkImiVkoSUsvxsfsplLPgXvwj/vVwiiqmcl0UyqEJ2Hqpoq+Riv5Kr04yEe/EyTgPDDANk\nj+/oxg+Bk5yRZJ1eno3CHAXkUhHUCknwOlEOMFSRFkYxeALDANM1ZtUlWv5GBg/K8oPDXGtKtSjK\nUaJIF1AxL8YlFXb8xuCi8T1HfHmDI/Lw23m4WyWR500y1fr4wmsa/7777sOtt96KM2fO4KqrrsLQ\n0BD++7//O92yEZOILxu68cL7x+BhWay9aCouXVQ2LivBEOljdlUufrpmNp7ZcRhP7ziM26+dhTnV\nk1dBJsKzceNGrFq1CnPnzoUwYIZ28+bNKdsHy7J48MEHcfz4cUgkEjz88MOoqKhI2fbHCoVMjDnV\nOuw52h30PQeAG1Zwwtsf8Wkf0RL6i3VKyKUiHG8zoDg38cJSspCCQGqFGFaHy7v/QI/TsC5dnq/G\noG0kt2JOdeTwo6IcJbr0FmhVUohEDHr01qC8jMDRmFqijZqbwzACVBdrcapzyB8iNqcqF00dRsyY\nkh31GPOz5SP9lni8/qqKtRgyO/zFNMIT4pEI6FmVTONgPgRWSCzOVfr3l6WSIiukj5e/AlyAvCq5\nGBqFJGJBEv7V1dJ3nHKJCGabK+oygcp8VElCzodMIkJFoTdvSyUTw2x3xfSopFpvqizUwGRzoX/I\nFnQ+ZRJh0HFH9WwOP2QCC2rw8V7zvT4Dr2k+kxXJwstwGhgYwPbt29HS0gKPx4OqqiryOBEpgWU5\nbP/4FD746gzkUhHuuGp22hP7iPHLrMpc/HT1HDy94xCe2XEYN1w+HUtmF2VaLGIc8cgjj2DVqlVp\nbcC+a9cuOJ1ObNu2DfX19Xjsscfw7LPPpm1/qSCucGZuuG+MDZDFUQggkhoTqzVEtlqKBdPiyCkK\ng5BhsLiuACzHQW90IEsl8RsZgcq2L1TMd3yAVwGLplBWFKpRrFNCLGKgkotRkK2A3mj3V74TCRk4\n3R7kqGX+ynPRyMuSQykT+UMNvQZspPeaV678LAVUcjF6fEXmQvQ/rVIKqVgIfYCRJJMIwaikQYZT\nrHLg507L8zdRTkTHDq1oq5KLUZavDvquskgDl4sNMnZ9hRsi4VN8A5XzLJU0Zjn3aOXLA7fuQ6uU\nYsiSeGW5UL28olANhUwEk9UFvckOqUiIOVNzIWQYnGwfxJDZGaTM+8a8okA9KqwtGlXFWvQN2pAX\nofS4D4mIQbZKiiy11D8BIImjxP+84ZBTX5VMtUIMk9VrIIlF3nvQ7eEgFjHIEjLI1chQlKv031+B\nYa0Vw+c8VyuDzemBTivzFy3h0waAL4U5ioRLmycCryfmb3/7Wyxbtgw1NTXploeYRFjtLvzp7SNo\nOK1HYY4Cd6ye7e//QUxeZlbm4O518/D09kP467uNGDQ7sPK8CvJAEgAAiUSS9sp6Bw4cwIUXXggA\nmDdvHhoaGtK6v1QQWn0qGhzHYUqhGnKpCIVxNggNJVzYYDiSMZp8CAQCCAUC5GXJ4fZEVzoDDQg+\nSpo4RLnMCQi5K8pVgGW5mEprIAqeCeqBj7VAnTz0+OoqvN4qgynQSAJ0WXIo5WL0DdrQOWCJaQwF\nnodEnqmh5zFcIYGC4Was8cz+15Zmoa3XDJ1WhrZeU+wV4sATED5WV5E9ygObDL6G01bHkP87X96Y\nL4eta8AyqidWvLqOQjbifYqGQCDAtHLvtZKlkkJvtEMX47pVysSoKdXCw3L+iYeppVo4XSrIJKKg\nIhMCgQBikfe6YRjBqDw93ykXMUxQP7eyfG8YrFQkhMPtGeW1DVetMB6P6IyKHPQN2qAZgx5PvAyn\nsrIy3HfffZg7dy5kspGHydVXX502wYizm64BC57ecRg9eitmV+Xi1u/O4P2iIc5+asuycN/Gc/HU\n6/XY8clpDJqcuO47NSmdpSImJueffz4ee+wxLF26FGLxyDNj4cKFKduH2WyGSjWS7yIUCuF2uyES\nhX9lZmcrIEqgglgoeXmxFSMfGrVXUZtRmYscbfTKa75lfagUYhQXZaG4KCvscrk6ddjcDYvNBU2f\n18tToFPD5nDHJXMqcXtYaDq9CnZenhqaDqP/MwBAJEKP0QGNWp6wjCUFNpisTpSVZEEdZyNVvswX\niXC0eQAza/JgsrrQb/Iq2FnZSuSFMWo1XSa/p0KnU/sVRRcEMDtZMIzAOx7D5zIvTw2xiAn62/c5\nP18d1aANvW4Ab08pi3pEEY81tpp2I6/l8vLUmFadB47jcKrbW3wsJ0c56lhC0VtdcHi8nrDQ331/\n6zgOdhYoyFEgVyuPur1YDNrdsLm9HpfA9QeG5ZCFkUOnU8FgdQcdU7jj54RC//lP1X1VWhz8d7hz\nqlKIUVYSOYzU6PDAxQmgkIljyhVpfHxclK2Ew+mGSiFBU3dwkbk5NTocOtnv/zs3d2Sswp0zTacR\nHg+H7GwFqkuzUB2yr3Q9m6IaTj09PSgoKEB2tndADx48GPQ7GU5EIhw61Y//efsIbA4PLl9cjtXf\nriaFmBhFiU6JX3z/XDz1t4P48Ot2DBjtuGXVjKR6jBATn6NHjwIAjhw54v9OIBDgpZdeStk+VCoV\nLJaR0A+WZSMaTQBgMITvmxIPeXlq9PXxn2k3mobDk9xu6MP0gAmkqkAJm8ODvkEb9CY73E5X2H35\ntjnQb4ItzEQWy3IQeDyYOiUXIo4FxyEumVOJh2X98vb3m/yfffL4QtqMJlvCMhZlSaGVC2G3OGBP\nIrwrFt+eX4q+PhMGDFb/cXCu8OfIZLLD6fYAACwmGxxWr1zGIRuMJhukYiH6+kxBYyMSMkHjM3Ke\nzVHfvf5rbJja0ixUl2ahtWMQdqcbpXmqmGNrNNkgE4sSurYNeiEUQgFqirylwsNtwzDoHTOpSBj0\ne+j9lKeSgHW6g44/kevCYrLDaLJBLZcErW8YPncOsTDqvaXXiyANGXKfrAajPSnZ+BB6TgHA43JH\n3Z9CJIBWLkJteVZMuXzjIBaGHwcfNotjlCwOa/B3Tq3Uv41w42I02uFhWRjEDPqkwRNX8T5PQ4lm\ndFGze6YAABkKSURBVEXVQH784x9j586d2Lx5M55//nncdNNNCQtBEBzn7c+0/eNTEIkY3LJqBr41\nk1+YBzE5ydHIcN9/zMf/29mA+qZ+/ObF/RTSOcnZsmVL2vcxf/587N69GytXrkR9fT1qa2vTvs94\nUSsk/Mriwpv3I5OI0D80rJQnmD/NMALUTclB3rDCPF6iZ9MVxisSMtCkydMUloDzEqlMeKChE+gt\nKsxVwOH2xFeAI85h8/WUmludC47jFwK5qK4g6dIM0fpBaZUS9A3akMujeW0qKMpVggOQHyH8LVaO\nWTQi9f1KJb7cpXgQCb2Nrfn05YrnEAKLRQDB9/GsytyUVixMJVENp8D41HfeeYcMJyJhnC4P/vf9\nY9hztAfZailuv3Y2Kov4x+QTkxeFTIz/s24utn98Cv/4qg2/eXE/fnDZdCyeUZBp0YgMsH//fvz1\nr3+F1WoFx3FgWRadnZ346KOPUraPFStW4PPPP8f69evBcRweffTRlG07VcxMoCG4P4E7gnJTmqdC\nj94a1ItnvJKMgjpeSUZtFgkZVBfz6/Hlq0KWaFU9gUDA22hOd+U+nVYOpUw8qupiumAYAUrDla3n\nefKiGflZKinkEhFKdOmbGKwp1cLlVuPrk30jMqVw+9lqKXoHrSjg0YZgfm0e2vvM6BywQMR4e1LV\nlmZBJhGO69SNqE/HwBM8FiX+iLMTvdGOP7xxGC3dJlSXaPCTa2aPKkVKENEQMgzWXVyDikI1Xnz/\nOP7n7SNoaB7Af6yonRBKHpE67r//ftxyyy3YuXMnNm7ciE8//RQzZsxI6T4YhsFDDz2U0m2OB3xv\n9HCJ2IDXcAqrFE5AJqLK4kuYV0VRGn3nMFeTuIdldlVuxGsgEr6CD+OReMO387MU6WvmnMRmRULG\n30g5XQgEAkiilNNPlmy1FPNr8njtg2EEKMlTwulmUZTrvb5y4riu5VIhzDbWW0FzDOF9tVFFKyIR\njp8x4Nm3jsBoceKC2UXYeOm0UdWLCIIv580oxJRCDf7n7SP4/HA3TrYNYeOl0zCzMv7Zd2JiIpPJ\nsHr1anR0dECj0eDhhx/Gtddem2mxJgSxPE4TihgqyUSc7M1WS1Fbmr5CFD7iKl2fxDqJ4utZFK1f\nVjLEU4GSL7EM0aklWrR2m5GrOfsnjeMxzIQMg6kl/DylopDw1drSLPQN2cfcqI96J5w8eRLLly8H\n4C0U4fvMcRwEAgE+/PDD9EtITEg4jsP7e89gxyenwAgEuG55Db6zoJQMcCJpCnMU+OXGc7Hzs9P4\nYO8ZPLmtHufNLMD6i2vGpBQpkVmkUikGBwdRWVmJgwcP4lvf+has1uSLM0wGinUKmKxOVKZBcRxr\nGIEAZXkqf0hPca4yqCfVWOSLpIN4ZtzjYcG0fLBs4mMylm/uaeVZMJidY5a3lEoijZNOK+fVAywT\nxFNmP1Msqhsdmi8RC9Ma1hiJqIbTP/7xj7GSgziLGDI78ML7x3Do1ACyVBLcdvWsUbX+CSIZREIG\n31s2FYumF+ClfxzDniM9ONQ0gNXLqrF0blHExGpi4nPDDTfgrrvuwjPPPIM1a9bgnXfewaxZszIt\n1oRAJhGlPRRoLCkJCCsMbbA6Qe2mmPjmHmMdX115Nsw2l7+AhEjIAAk6cHSasVX6xSJhxOILRGqY\nU5ULm8MDjVLMq+hDpkl3rlw8RDWc0tmZnTg7+aqxB1v+cRwWuxszpmTjR6tmkheASBsVhWr8cuMC\n7P6mAzs+OYUt/ziOf+1rwzVLq3DutLxx9bAlUsPll1+Oyy67DAKBAG+88QZaWlowffr0TItFjDNy\nNFKYnCym8iyYcLahVUmhTVEu8dTSyTmG8eAzZCfKK0chE4/rAgzjGcqqJlJCr8GK1z5sQn1TPyQi\nBv+xohYXzS8hxZVIOwwjwPJzSzG/Ng/vfN6MTw924dk3G1BRoMY1S6swuyqHQkTPEnbv3o2pU6ei\nrKwMu3btwvbt21FXV4fa2low5GUkAhCLhPjW7KKM9ZkiCOLshN40RFJY7W7s+OQU7v/LV6hv6sf0\n8iz8+qZFWH5uKRlNxJiSrZbi+sum45FbFmNRXT5ae0z4/d8OYtPzX+GzQ51wudnYGyHGLX/961/x\nhz/8AQ6HA8eOHcM999yD5cuXw2q14vHHH8+0eARx1pGrkSE/a/xW0xufkN5ztkMeJyIhbA43dh1o\nxz/2noHV4Ua2Wop1F0/Fwun5NLtPZJSCHAV+fNUsrDzPhA/2nsFXjb144b1j2PHJaSybV4zzZxdR\n/PwE5K233sK2bdsgl8vxxBNP4OKLL8b3vvc9cByHlStXZlo8ghgjxu79SrnJBDEaMpyIuNAb7dh1\noB2f1HfC5nBDKRNhzbJqLJ9fCukYNaAjCD6UF6jxo+/OxJpl1f5r9u3PW/D25y2YWqrF+TMLsWB6\n/rjtTk4EIxAIIJd7Dd69e/diw4YN/u8JgiAyia/8PT2Ozn7IcCJiwrIcGpr1+OxgJ+qb+uFhOWiU\nEly+uArLzy2Nu/kcQYwlORoZ1l40Fd9dMgUHjvfhi4ZuHGs1oKl9CFv/eQJTSzSYM1WH6eXZKC9Q\n+atQEeMLoVAIo9EIq9WKxsZGLFmyBADQ0dEBkYieQcTkYKSq3llaNnCC4utdlK7eU8T4gd42RER6\nDVZ80dCNfx/ugt7oAACU5avwnQWlOG9GITWyJSYUMokIS2YXYcnsIuiNduw52oNvTvThZPsQTrQP\nAQDEIgYVhWpMLdZiSpEa5QVq5GfLKV9vHPCjH/0IV199NdxuN9asWYP8/Hy89957eOqpp/CTn/wk\n0+IRBDGJKctXQcgIUJRLOWFnO2Q4EUEMWZzY19iDPUd7cLrTCACQSYRYNq8YF84txpRCNYXGEBOe\nHI0MK8+rwMrzKmCyOnGkWY+T7UM41eH91zRsSAGAVCJEeb4K5QVqlBeoUFGgRrFOSZ6pMeayyy7D\nOeecA4PB4C8/rlQq8fDDD2Px4sUZlo4giMmMSMiM6iVGnJ2Q4UTAZHWi/mQ/vjrWi6MtenCcNxxg\nVmUOFs8owIJp+ZS/RJy1qBUSnDezEOfNLAQA2J1utHSZ0Nrj/Xemx4ymjiGcDDCmREIBinVKlBeo\nMaVQjWnl2SjOVdCkQpopKChAQcFIB/lvf/vbGZSGIMYeesIQRGYhw2mS0j9ow9cn+/HNiT6caB/0\nN2+rKtZg8YwCLKorgJYa1xKTEJlEhOkV2Zheke3/zuHyoL3PjDM9ZpzpMeFMjwltvRac6THj34e6\nAABalQR1Fdmoq8jGrMpcZKtT03ySIAiCIIjxARlOkwSW49DabcLhUwP4+mQfzvSYAXhnr6pKNJhf\nm4dza/OQn03xuQQRilQsRHWxFtXFWv93HpZF14AVpzuNaGw1oLHVgD1HerDnSA8AoESnxMzKHMyq\nysG0siyIReS1JQgiNVBtCILIDGQ4ncUMWZw40jyAhtN6HGnRw2R1AQCEjACzqnIwvyYP82p0yFLR\nzDhBxIuQYVCap0JpngpL5xaD4zh09FtwtMWAhuYBnDgziH/ua8M/97VBLGIwrSwLsypzMLMql8L6\nCIJICCHjfW4wDD0/CCITkOF0lsBxHPqG7GhqH0RT+xBOdgyho8/i/12rkuCC2UWYVZWDWZW5UMjo\n1BNEKhEIBH5D6pKFZXC5PTjRNoSG5gE0NOv9//BRE7LVUsyqzMH0imyUF6hRlKMgRYggiJhUFWvR\n1mtGeYEq06IQxKQkbdozy7J48MEHcfz4cUgkEjz88MOoqKjw//7666/jtddeg0gkwm233YaLLroI\ner0e99xzD+x2O/Lz87F582Z/w0PCi9vDwmByQG+0o8dgQ1uvGR19ZrT3WWC2ufzLScQM6iqyMbsq\nF7Mqc1CSp6QZboIYQ8QiIWZW5mBmZQ7WATCYHGhoHsCRZj2ONOvx2aEufDacHyURMSjJU6FEp4Qu\nSwadVgadVo4slQQKmRhyqRBChqr4pYt//etf+OCDD/Dkk08CAOrr6/HII49AKBTiggsuwO23355h\nCQnCi1QixNRSbewFCYJIC2kznHbt2gWn04lt27ahvr4ejz32GJ599lkAQF9fH7Zs2YIdO3bA4XBg\nw4YNWLJkCf74xz/iyiuvxLXXXovnnnsO27Ztww033JAuEcFxHE53GWG1u4f/BgDOHzs8/Cc4jHzB\nITC2mAuKM+YQ3JSO860bYXsc5/3nZlm43CzcbhYuz/BnDwuH0wOz3Q2LzQXL/2/v/mOqLP8/jj/P\nuQ+/5EdqyIYTSkw+wwoNqelC+gNdZQbLYoI/mEMTWWVSGkkpNo6os6gNbdPhnEPDEGvt09L6J2WE\nOYepAdIPQy0tC1LhHETg3NfnD/IohPGrc+5zvt/3Y2Ny7gPH13W9ua+L69w3993eSYu9g2u2Dnqf\n2mwCxowM4D+RI5k4biQTx91FRJjcyFMITzIq2I8ZsWOZETsWXVec+62Vs5eu/XWxie6LTjT+2nLH\n7/f31RjhbyHA14LFYsaimfDRzGiaGR+t+7FFM3d/WMxYzCbn11k0M5rZhGY2o2kmLGYTmnNbz8/N\nZhO6rtCVQte7x7Tuz29t0//apnSFQ1foqvtG2ZPuHeV1l+S1Wq1UVVURExPj3Jafn09xcTEREREs\nW7aM+vp6Jk2aZGBKIYQQnsBlC6eamhpmzJgBwJQpU6itrXU+d/r0aR566CF8fX3x9fUlMjKShoYG\nampqyMrKAiAxMZGioiKXLpzOXmqhsLTGZa//b/Lz0Qge4UN0xEhGh/hz911+hN4VwLi/3qWWy4UL\n4T3MZhNRY0OIGhvi3NbZpdN07TrN19ppammn6Wo7LfYO2m500dbe+de/XVy13aDLoehy6Dh0z/oL\n8QfGj+aVeVOMjjEocXFxzJw5kw8//BAAm81GR0cHkZGRACQkJFBdXS0LJyGEEK5bONlsNoKCbp2D\nq2kaXV1dWCwWbDYbwcG33pUMDAzEZrP12B4YGEhra2u//8+YMUN/d3PMmGD+O2XckL9fCCH+TWPD\n5RQcV9m/fz+7d+/usa2wsJDZs2dz7Ngx57bec1dgYCA///zzP772cOYhV7yOO3hLVsn57/OWrN6S\nE7wnq7fkBNdlddnCKSgoCLv91sUJdF3HYrH0+Zzdbic4ONi53d/fH7vdTkhIyN9eVwghhBis1NRU\nUlNT+/26vuYnmYuEEEIAuOyPYOLi4qisrAS6/9A2Ojra+VxsbCw1NTXcuHGD1tZWzp49S3R0NHFx\ncRw5cgSAyspKpk6d6qp4QgghxN8EBQXh4+PDhQsXUEpRVVVFfHy80bGEEEJ4AJcdcZo1axZfffUV\naWlpKKUoLCxk165dREZGkpSUxKJFi5g/fz5KKXJycvDz8yM7O5vc3FzKy8sZNWqU8wpHQgghhLu8\n9dZbrFq1CofDQUJCApMnTzY6khBCCA9gUkruPy2EEEIIIYQQ/0SuVy2EEEIIIYQQ/ZCFkxBCCCGE\nEEL0QxZOQgghhBBCCNEPl10cQgyOUorExETuvfdeoPumwa+++qqxoYZA13XWr1/Pd999h6+vL1ar\nlXvuucfoWMPyzDPPOO/rMm7cODZu3GhwoqE5deoUb7/9NqWlpZw/f57XX38dk8nExIkTyc/Px2z2\nvvdRbm9TfX09WVlZzn0oPT2d2bNnGxtwgDo7O8nLy+PixYt0dHSQnZ3Nfffd57U16qs94eHhXlsf\nV/LEMXMw9du6dSuHDx/GYrGQl5dHbGysW7P2Hp/nzZvHhg0b0DSNhIQEXnzxRY/o448++oiPP/4Y\ngBs3bnDmzBmKiorYvHkz4eHhALz00kvEx8cblnUgc0Rf9Xb3fHJ7zjNnzlBQUICmafj6+rJ582ZC\nQ0OxWq2cOHGCwMBAAN5//306OztZtWoV7e3thIWFsXHjRgICAlyWs3fWO81RntanOTk5NDU1AXDx\n4kUmT57Mu+++S3Z2NleuXMHHxwc/Pz9KSkrcmnMw86RL+1QJj3Du3DmVlZVldIxh+/zzz1Vubq5S\nSqlvvvlGLV++3OBEw9Pe3q5SUlKMjjFsO3bsUHPmzFGpqalKKaWysrLU119/rZRSau3ateqLL74w\nMt6Q9G5TeXm52rlzp8GphqaiokJZrVallFJXrlxRjz32mFfXqK/2eHN9XMkTx8yB1q+2tlYtWrRI\n6bquLl68qObOnevWnH2Nz8nJyer8+fNK13W1dOlSVVdX53F9vH79erVv3z5VVFSkDh061OM5o7IO\nZI64U73dOVb1zrlgwQJVX1+vlFKqrKxMFRYWKqWUSktLU83NzT2+t6CgQB04cEAppdT27dvVrl27\nXJazr6yD2YeM7NObrl69qpKTk9Xly5eVUko9+eSTStf1Hl/jzpwDnSdd3afe8fbl/wN1dXVcvnyZ\nRYsW8fzzz/PTTz8ZHWlIampqmDFjBtB91Ky2ttbgRMPT0NDA9evXyczMJCMjg5MnTxodaUgiIyMp\nLi52Pq6rq+ORRx4BIDExkerqaqOiDVnvNtXW1nL48GEWLFhAXl4eNpvNwHSD88QTT/Dyyy8D3Uef\nNU3z6hr11R5vro8reeKYOdD61dTUkJCQgMlkYuzYsTgcDv7880+35ew9Ph8/fpyOjg4iIyMxmUwk\nJCRQXV3tUX387bff8uOPPzJv3jzq6uo4cOAA8+fPZ9OmTXR1dRmWdSBzxJ3q7c6xqnfOoqIiYmJi\nAHA4HPj5+aHrOufPn2fdunWkpaVRUVEB9NzX3DGmDmSO8sQ+vam4uJiFCxcSFhZGU1MTLS0tLF++\nnPT0dL788kvAvb9LDHSedHWfysLJAPv372fOnDk9PkJDQ1m2bBmlpaVkZWWxevVqo2MOic1mc542\nAaBpGl1dXQYmGh5/f3+WLFnCzp07nfd28cb2PP7441gst87MVUphMpkACAwMpLW11ahoQ9a7TbGx\nsbz22mvs3buXiIgItm3bZmC6wQkMDCQoKAibzcaKFStYuXKlV9eor/Z4c31cyRPHzIHWr3d2d/+c\n9h6f16xZ0+PUq5t5PKmPt2/fzgsvvADAo48+ytq1a9m7dy9tbW3s27fPsKwDmSPuVG93jlW9c4aF\nhQFw4sQJ9uzZw+LFi2lra2PhwoVs2bKFkpISPvjgAxoaGrDZbAQHB7slZ19ZB7MPGdmnAM3NzRw9\nepS5c+cC3afJZWZmsm3bNrZu3crGjRtpbm52a86BzpOu7lNZOBkgNTWVTz/9tMfHgw8+SFJSEgDx\n8fH8/vvvKC+8xVZQUBB2u935WNf1v+2Q3mT8+PEkJydjMpkYP348I0eO5I8//jA61rDdfl6v3W4n\nJCTEwDT/jlmzZvHAAw84P6+vrzc40eD8+uuvZGRkkJKSwtNPP+31NerdHm+vj6t46pg5kPr1zm63\n252/mLpD7/E5ODiYq1ev9sgTEhLiMX3c0tJCY2Mj06ZNA+DZZ58lIiICk8lEUlJSn31qVNa+xp87\n1dvoseqzzz4jPz+fHTt2MHr0aAICAsjIyCAgIICgoCCmTZtGQ0NDj/xG5BzMPmR0nx46dIg5c+ag\naRoAoaGhpKWlYbFYuPvuu4mJiaGxsdHtOQcyT7q6T2Xh5CG2bt3K7t27ge7TD8LDw50rY28SFxdH\nZWUlACdPniQ6OtrgRMNTUVHBpk2bALh8+TI2m40xY8YYnGr4Jk2axLFjxwCorKwkPj7e4ETDt2TJ\nEk6fPg3A0aNHuf/++w1ONHBNTU1kZmayevVqnnvuOcC7a9RXe7y5Pq7kiWPmQOsXFxdHVVUVuq5z\n6dIldF1n9OjRbsvZe3y+fv06I0aM4MKFCyilqKqqIj4+3mP6+Pjx40yfPh3oPqKTnJzMb7/9BvTs\nU0/I2tf4c6d6GzlWffLJJ+zZs4fS0lIiIiIAOHfuHOnp6TgcDjo7Ozlx4oSzb48cOeLMOXXqVLfl\nhMHtQ0aP/0ePHiUxMdH5uLq62nmanN1u54cffiAqKsqtOQc6T7q6T03KGw9r/B907do1Vq9eTVtb\nG5qmsW7dOiZMmGB0rEG7efWi77//HqUUhYWFXtmOmzo6OlizZg2XLl3CZDKxatUq4uLijI41JL/8\n8guvvPIK5eXlNDY2snbtWjo7O4mKisJqtTrfWfImt7eprq6OgoICfHx8CA0NpaCgoMfhek9mtVo5\nePAgUVFRzm1vvPEGVqvVK2vUV3tWrlzJli1bvLI+ruSJY+Zg6ldcXExlZSW6rrNmzRq3/oLX1/hs\nNpspLCzE4XCQkJBATk6Ox/RxSUkJFouFxYsXA1BVVcV7772Hv78/EyZM4M0330TTNMOyDmSO6Kve\n7p5PbuYsKytj+vTphIeHO48ePPzww6xYsYKSkhIOHjyIj48PKSkppKen09TURG5uLna7nVGjRvHO\nO+8wYsQIl+W8Pes/zVGe1Kfl5eUAPPXUU5SVlfU4KrNhwwZOnTqF2Wxm6dKlzJw50605BzNPurJP\nZeEkhBBCCCGEEP2QU/WEEEIIIYQQoh+ycBJCCCGEEEKIfsjCSQghhBBCCCH6IQsnIYQQQgghhOiH\nLJyEEEIIIYQQoh+ycBJCCCGEEEKIfsjCSQghhBBCCCH68T8CEXhQgkniZgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f2e3733af98>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.traceplot(trace_p);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hierarchical model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Average ELBO = -43.175: 100%|██████████| 200000/200000 [00:23<00:00, 8435.73it/s]3, 8478.86it/s]\n",
      "100%|██████████| 2000/2000 [00:03<00:00, 477.25it/s]\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as hierarchical:\n",
    "    \n",
    "    eta = pm.Normal('eta', 0, 1, shape=J)\n",
    "    mu = pm.Normal('mu', 0, sd=1e6)\n",
    "    tau = pm.HalfCauchy('tau', 5)\n",
    "    \n",
    "    theta = pm.Deterministic('theta', mu + tau*eta)\n",
    "    \n",
    "    obs = pm.Normal('obs', theta, sd=sigma, observed=y)\n",
    "    \n",
    "    trace_h = pm.sample(2000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1kAAACCCAYAAAC0EAJCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl4VOXZP/DvmX2fTJJJgGwkYUkCRAEBFUUQEVEoIBQp\nr+BSreWtWv3ZVqsWl+JWtdaqtcXWBVAWUVxel7YgFquCgoQlC1tIyL5NltnX8/tjMpOZyczkTGYm\nk+X+XJeXJJk55z5nzsw893me534YlmVZEEIIIYQQQgiJCV6iAyCEEEIIIYSQ4YSSLEIIIYQQQgiJ\nIUqyCCGEEEIIISSGKMkihBBCCCGEkBiiJIsQQgghhBBCYoiSLEIIIYQQQgiJIUqyCCGEEEIIISSG\nKMkihBBCCCGEkBiiJIuQQezgwYO44YYbcNddd+Gaa67B8uXL8cUXX+CWW27B3Llz8eSTT+LgwYNY\nvHix33N8fyaEEEJijb6fCAlPkOgACCHhHT9+HLt27UJRURFuu+02bNq0CZs3b4bBYMCcOXMwZcqU\nRIdICCFkBKLvJ0JCoySLkEEuMzMTRUVFAIDs7GwolUqIRCIkJydDLpejs7MzwRESQggZiej7iZDQ\naLggIYOcSCTy+1kg8L83MmHCBLAs6/3ZbrcPSFyEEEJGNvp+IiQ0SrIIGeJUKhXq6+vR1tYGlmWx\nZ8+eRIdECCGE0PcTGdFouCAhQxyPx8Pq1auxYsUKaLVazJ07N9EhEUIIIfT9REY0hvXtxyWEEEII\nIYQQEhUaLkgIIYQQQgghMURJFiGEEEIIIYTEECVZhBBCCCGEEBJDlGQRQgghhBBCSAwNy+qCLS36\nRIdACCEkRrRaZaJDiJlYfD9pNDK0t5tiEE38DZVYh0qcwNCJleKMvaES61CJE4hNrKG+o6gnixBC\nCBlCBAJ+okPgbKjEOlTiBIZOrBRn7A2VWIdKnEB8Y6UkixBCCCGEEEJiaFgOFyRkJGFZFkaLAx0G\nKwwmO4wWO/g8HsRCHmQSIUalyCAWDp27SoQQQkg82R0u6PQWaJOkiQ6FDGOUZBEyBDicLrR0mNHS\nYUFrp7nn3x1mtHSaYbY6Qz6XATB2tArTJqTisimjoVaIBy5wQgghZJA5W9eJDqMVLheL9DRVosMh\nwxQlWYQMMgazHadrOlDVqEd9mxH1rUY06cxwsWyvx4qFfKQmSaBVS6FRiqGQCiGXCuFysbDanegy\n2VDTbEBlXRfONXThw/+ew5XTMvGj2WMhkwgTcHSEEEJIYhktdgCAxRb6BiUh0aIki5AEY1kWdS1G\nfFfRjGNnWlHTbIBvOiUVC5A7RolRyTKkJUmh9flPKROCYZg+92G02HGgtAmfHzyPf31fg+8rmnHb\n4iIU5mjid2CEEELIoNT39yaJnU6DFXanC6nqkTU8k5IsQhLEZHHg6xMN+PJIHRra3OVDBXweJmYn\nYWK2BvkZKmSkKpCkEHFKpMKRS4SYPz0Tcy4YjU8PnMf/fVOF57eX4H+unoB5UzNicTiEEEIIIX4s\nNgfKz7cDACVZhJD46jTa8Mm3VfjqaAOsdicEfB4umqjFRQVpKM5PgUQUv7elUMDH0styMWlsMl56\n/xi2/PMkTBY7rrtkbNz2SQghhITicrGwO1wQiwauQJPnvmWQUfgkxupajIkOIWEoySJkgJitDnx6\noBr/PlQDm90FjVKMxZfm4PILxkAlEw1oLOMy1Xho7XT8YdsRvPefSoiEfCy4KGtAYyCEEEJKq3Qw\nWuyYPkEL4RBaX2mwsdmdsNidA96eIKFRkkXIADh8sgVv//skOgw2qBUi3DBvLC6/YAwE/MQtVZem\nkeE3a6bhqa2HsX3PaWjVUlw4PjVh8RBCCBl5PEUorHbXgCVZngH4LIZPV1bJmVa4WJZTsupwutBh\nsCJFJYl6OkKfRvD0N1qMmJA46jLa8NJ7x/DK7uMwmO340eyxePqOSzBvWmZCEyyPtCQp7l5RDKGA\nh799VIrzTfpEh0QIIWQY0ZtscDhdiQ5jwBnMdjhd8T3u+lYjyqt0AOCtQOx09Z04nq3rxJm6TjR3\nmOMa30iX+FYeIcNU6TkdNrz+HY6cbsXErCQ8dutMLLs8b9AtDJw7WoXblxTBanfixV3H0GW0JTok\nQmKmtrYWX375JZxOJ2pqahIdDiEjisXmQGmVDodONic6FH89XVlxoTfZcOJcG07VdMZnB93ON+vR\nGZDEMhy6jvSm7hL2YdbYJNGL23BBl8uFRx99FCdPnoRIJMLGjRuRk5Pj/fvOnTuxfft2CAQCrF+/\nHvPmzcMTTzyBiooKAEBLSwtUKhV27tyJjRs34ocffoBcLgcA/OUvf4FSqYxX6IRExeVi8f7+Snx6\noBp8HoNV88bh6plZ4MW7Sz4K0yemYcUVeXjvP5X4xyfl+OWPiwd1vIRw8emnn+LVV1+F2WzGjh07\nsHr1avzmN7/B0qVLEx0aISMCl16VRIpXdEaLAwDQabTGaQ/+XL7nmb66B424JVl79uyBzWbDjh07\nUFJSgqeffhqvvvoqAHcCtWXLFrz33nuwWq1Ys2YNZs+ejYceeggAYLfbsWbNGvz+978HAJSWluLv\nf/87kpOT4xUuITFhstjxt4/KcLyyDWlJUtyxdBJyRw+N1eQXXZyDivMdOF7Zhj2HanH1DCqEQYa2\n1157Ddu2bcONN96IlJQU7N69G7fccgslWYQMENcgTbK49PZEY6CHR7p8yiRSjjV4xG244OHDh3H5\n5ZcDAC688EKcOHHC+7djx45h6tSpEIlEUCqVyM7O9vZgAcDWrVsxe/ZsTJw4ES6XC9XV1diwYQNW\nr16NXbt2xStkQqLSpDNh4+bDOF7Zhsl5ydhw80VDJsECAB7D4LbrCqGSCfHuvjOobqT5WWRo4/F4\nUCgU3p/T0tLA4/X/a89ut+PXv/411qxZg5UrV2Lv3r2orq7GT37yE6xZswaPPPIIXHGegxGO1eZE\nZX3XiJz/QganQZpjxZ0nuRyoESG+HzuD7ZTHO6EdzOLWk2UwGPy+3Ph8PhwOBwQCAQwGg99wP7lc\nDoPBAACw2WzYvn27N5kymUy48cYbccstt8DpdGLdunWYPHkyCgoK4hU6IRGrauzCCzuPQm+y45qZ\n2Vg5Nx883tD7YFErxPjp4iK8sPMoNn1cikdungHRIJtDRghX48ePx9atW+FwOFBeXo533nknqu+O\njz76CElJSXj22WfR0dGBZcuWoaCgAPfccw9mzZqFDRs2YO/evViwYEEMj4K7s/Wd6DK551TmjRk6\nN3jI8DVYe7LiPSfL07EU98p93VxBFvzSdVlgsTkxJlU+IDFwwbLsgJ2TwSBuPVkKhQJGY88CZC6X\nCwKBIOjfjEajN+n69ttvMWPGDO/PUqkU69atg1QqhUKhwMUXX+zX60VIopVW6fDMO0dgMNmxduFE\nrLpy3JBMsDym5KXgqumZaGgz4cP/nkt0OIT024YNG9DU1ASxWIwHH3wQCoUCjzzySL+3d8011+CX\nv/wlAHdjgc/no7S0FDNnzgQAzJkzB998801MYo+G2eZIdAhxw7IsrHaarD/QmnQmfFfeBLsjsl5S\ndoSu9uspDZ+QlkD3KT9V24HzzfqEvwbNHaaE7j+R4taTNW3aNOzbtw/XXnstSkpKMGHCBO/fiouL\n8ac//QlWqxU2mw1nz571/v2bb77BnDlzvI+tqqrCPffcgw8++AAulws//PADli9fHq+wCYnI0TOt\neGX3cQDA+mWTcVFBWoIjio0VV+Tj6NlWfP7deUybqEX+GHWiQyIkYjKZDPfddx/uu+++mGzPU3zJ\nYDDg7rvvxj333INnnnnGe2dWLpdDrw8/zFajkUEQg7WAtFr3jUiDyQapWAA+n4ekDgvA50MlF3n/\nPhjEMpbjZ1qh67JgRlE6ZBJhzLYLxDbOeBvoWMtqOqFQSMAXC6FNlnF+niZZDlWnu/hDqJhVSncF\nvpQUBVTy3gvpduitcLEsklWSfkQeXFKLESaLA+okadjY+qvD4oDJzkIo4MV8277b85y71FQFVK3u\nZCYlVQGpWOD9W3KKoteSMap6PRxOF5I0srhdS57teuJwx6mM+U1oh9OF0so2ZI9SQqPs3zUSr3MQ\ntyRrwYIF+Prrr7F69WqwLIsnn3wSb7zxBrKzszF//nysXbsWa9asAcuyuPfeeyEWiwEA586dw7Jl\ny7zbyc/Px9KlS7Fq1SoIhUIsXboU48ePj1fYhHB27GwbXtl9HDyGwd0ri1E0dvgUZhGL+Lj12kI8\n884RvPFpBR65+aIBWySSkFgpKCjoNTRFq9Vi//79/d5mQ0MDfvGLX2DNmjVYsmQJnn32We/fjEYj\nVKrww/Ta26O/q6vVKtHSoofBbMeJc21QyUQoGpuMzk4zuoxWsA4nWloGx5xKT6yxUlXXDgCoqeuI\naaM71nHGUyJi7dK711PS6QzgObn1JGq1SrS1GbzPDRWz5+9tbQZYTb0T5wNljQCA/DFqqGQiiEXR\nfxd1dpphtjkgZFgAKTE/n+3tJnTpzRDy+THdduBr7zl3ra0Gn3/rIREJvD83NnX1Wjqmq8sMh8sF\nuZCHFnHsv9t94/TEAbivgb6SLBfLwu5wcV7uplFnwvnGLpyv70CqWorc0UrwI5h7G4v3U6gkjVOS\ndfvtt+P666/HVVddBaGQ250jHo+Hxx9/3O93+fn53n+vWrUKq1at6vW8TZs29frdbbfdhttuu43T\nfgkZCCcq2/Dy+8fBDMMEy2NitgZXTsvAFz/U4aOvq7Diivy+n0TIIOI7tNxut2PPnj0oKSnp9/Za\nW1tx6623YsOGDbjkkksAAEVFRTh48CBmzZqF/fv34+KLL446bq7MVvewQM88rFhyudghPex5ODNZ\n7GhoM2J0ysDPtYm4iIHPSLVo5+Ocre8Ej2EwszC939vw8IQRr5F0njlSA/UW8p2TFXhMiZ4XJxUJ\nIhrCfPJ8BzqNVlw4LhUSUWR9Qa2dZsglgoS8N4LhlOr97Gc/w1dffYWFCxfisccew7Fjx+IdFyGD\n1olzbfjze8fBMBi2CZbHyrn5SFVL8OmBapxr6Ep0OIT0m1AoxKJFi3DgwIF+b+Ovf/0rurq68Je/\n/AVr167F2rVrcc899+Cll17CDTfcALvdjoULF8Yw6vDi1XSqbTHg+4pmGC32OO2BRONweTOqm/To\njGDh+OZ2Ew6UNcJqi24+W6Q5UqzrXAYr8DDYuFxsT2IzUPcpBv9pAdAzVy0Yh9MFu8PlXVvMzHGh\n5MA5Z4PpEuGUIs6YMQMzZsyAxWLB559/jrvvvhsKhQIrV67EmjVrIBL1HkNLyHB0qqYDL73nnoN1\n14opmDSMEywAkIgEuGVRAZ7dXoLXPynHhptnQCiIW70cQmLqgw8+8P6bZVmcPn2a82iMYB5++GE8\n/PDDvX6/devWfm9zoFhsDpgsjl5D7LpMNlRUt6MwRwOlzP1dXtvirvbbZbRBHuN5T6R/HE4XbHYX\nZBKBN9FwRFCEorL7JplOb4nJXX4Xy+JMbSfSNFIkKcTe35+q6UCSUoy07rlOvi1eFiNjDafvKpq8\n/462kp7d4UJtiwGZWnnYIfvh8opEJx2+uw8Xy6GTzXGPZaBxbi0dPHgQjz/+OF544QVcfvnleOih\nh9DW1ob169fHMz5CBo2aZgNe3HUMLheLO6+fjMm5KYkOaUAUjk3G3KkZqGs14uNvqhIdDiGcHTx4\n0Pvfd999BwB44YUXEhxV/AW7W1xyphWnajt6VearaTbAxbKobTb0eo7TOYhuCY9wped0OFbZGnVl\nxWgb/Z7ndxps0OktqDjf7v2b1e6ETm9BZX1PoQO/K4h1D3EdHD2k7uPoiqA3sP976b/6ViOa2k04\nWdPB+Tm9E5nh/T4ezJ9TnHqy5s2bh8zMTKxYsQIbNmyAROK+EzZr1iysWLEirgESMhi0dJjxx50l\nMFsduH1JEYrzUxMd0oD68dx8HD/bik+/rcbU8alDapFlMnI99dRTiQ4hvgJaU8EadEaLHY1tPcU2\nnE4X4DOhPNhyQTyGgYtl4RysaxzFmNFiB6szDepeFs+cFt8kq7XTjKpGPYrzUziPMIh2iaKeuUwc\nr42Ahx2vbIOLZTGrMH1A15ByOlm/cyQS8GCyAg6XK+icpS6TDR16K7LTo6s6J4mySIfnhomlr6Fz\nie6uCifOodW29r5BNFhwSrLeeustyOVypKSkwGKxoLq6Gjk5OeDxeNi9e3e8YyQkobqMNjy/owSd\nBht+Mn88Lpk0KtEhDTipWIBbri3Ec9tL8Pf/K8Ojt8ygaoNk0LryyivDNuD27t07gNEkjotlUVHd\nDruzZ1gZl/YOn8eDy+kMOf+lorodNocz6ptNR063QC0Xh1w42elyobndDG2StFcJ6v4KVnzheGUb\nVEopJmZEVpUsIXxeknaDe+6KTm9BuoZbWfVoE5tIe9L8CjKADTunKpbrOdnsToi6byYcO9sGi82B\nmYXp4HUfv8An4eroPo++yqp0AIBUtSSipQICz4/vfqoau6BRSqAOUqY+WiznnyLXabCipcOC/AxV\n1NdPf15im92J880GZKUpOFUc7CtEo8WO45VtmJilgUYpDv/gKHFKsr788kvs3r0bu3fvRltbG37+\n85/j5ptvxg033BDX4AhJNIvNgRd2HkVzuxnXXZKDBTOyEh1SwhSNTcb8aZnY+0Mt3t9fiRuupKUU\nyOC0ZcuWRIcwIMK1V1wsi+/Km8I8IvS2+DwGdmfoYTgdxt6N0kh5FhVu7jCFTLLqWoyobzNCb7Jj\nQlaS39+adCYkKcXeRjMXx862wWZ3hlzPMFGdAXaHC7ouC7Qaqd/xsCyL+jYTkn0agkFDjCDuaKvd\nnWvo4pzQBfI9v1znZ3kqaEaiw2BFxfl2jEmRIztdCUt3L6DLxYLHd+/Vd9/h4oi0M/dEpS7o7w1m\nOxp1JjTqTLi4KPY3asNdu9Fe1+XdQ0K1SRKoFWJYbU5Y7E7OyWK4YhdcVDXqodNbYHe4UJijgcPp\nnqc2KlkWcfVBAN6e/aqGLmiU2qhi6wunWzY7d+7E22+/DQDIyMjA+++/PyQm+hISDZeLxV8/LEV1\nkx6XF4/G9XPyEh1Swq2cl490jRT/+q4GJ33G4hMymGRkZCAjIwNarRZlZWX4/vvv8f333+PAgQPY\ntWtXosMbEKGKIgQ2uII1MPndDVGnKza14dr1Vhwoa4TB3DMXh0vDz2Z379/U3dC2+xxTp8mG78qb\nIur9MFntcHA4prN1nThe2cZpmy6WRafRBpZlYTDbYXf09GRwje1MbQfONXahSee/hlqn0YaaZr1/\nLFG2mCNJSiMVbMvRNvD7M3/L0zNV32ZETZC5hoC7dyRQsF42i82BI6daUNts8HttQ7EHriHWvUmn\nk/t7qbndhEMVzTBFcOydPjc+uJxyz7E26IxwOF0oq9J556exLIu6FgNMFv8E17PdI2daUF6t454A\nB07M48jz/vEM5/ScQ0+yWnKmlfO2AHePnO/nkJXD6xktTkmW3W73qyAYTXUmQoaKbXtP49jZNkzO\nTca6ayYO2PjxwUws5OOni4sABvjHJ+X9ustIyEC58847sXnzZrzwwgv46quv8OKLL+Ls2bOJDiv+\n2OiGhXka4pHcxe8y2qAPsV5XdZN7oU/fuWFcSnEz3S0UtjuQYI1c33ljLMviTF0n2vX962nzhNTS\naXbP0/JU8XO6QpY+r27Uo7xah/o2E06ca8OR0z0Nv+8rmkPejGJZFueb9DBZ7DBaes+38o3H91xx\nfU3O1HbiQFkjzjV0BU0oosV5SpbPA41m/6Shy2SDI0jy0Wns+X3E63IFqPOZr+MJxWJzoNP3WmXc\nBWC+K2/y9np5NLaZYHU4UdtqQMkZbol3MH29bueb9DhY1gSny4XzTQY4XC60dXG/jht9EvTA1yZY\nsu97TTXpTOgy2VBWrUNplQ66LitqWgw4Vtnqd2MjUOm54L12sRb4Uebbwx7qZoh7Tql/7PWtRgCI\naM2uaHFKsq666ircdNNN2Lp1K7Zu3Ypbb70VV155ZbxjIyRh/n2oBnsP1yJDK8f6ZZMH/zj9ATQu\nQ41Fs3LQ2mnBzn1nEh0OISGdO3cOmzdvxoIFC3Dbbbfh3XffRXPzMCoTHKbhFqzxGlR3C0ZvsgWd\nm8KVp4EWdBfd//cMG7LZnWhuN3v/brU7oeuyeH8ObGD3NAh7N7h9G5B6sx2tnWacrGkPuwBrdaPe\nu92K6tA98o7uxlzJ6VYcOdMStLHa0Z3QdXWfO0+sLOuef9RusOJUTQcqqtu9jTwA6DDYUN9mxLEw\nPWbBej/OdyesvoIdaWuX+/w2tZtQcd6/Ml273ormDnOQZ8Web2zlPgmn3mhDWZWuV9U8XZcF5dU6\nHDrZjLoWQ78KdfSVAAb20FisTm8yFrj2mMGnN8npcvW7dzcwefMwWx0wWx2obzOCBesXm+81HKpA\nRyz4blVvsvklbIdPhf685NIrHLj9aHo2vZ8CPtdEqJ7OktOt+L4i8Z/1nFqOv/71r7F27VqcO3cO\nNTU1WLduHe699954x0ZIQpScbsX2vaehlotwz8oLIBVHPuZ3uFt6WS4ytQr8p6QeP5xqSXQ4hASV\nkpIChmGQm5uLkydPIj09HTZbfEs2D6TWTovfz76Nj6NnIxtKA8CvHDfgHoLkm6w5nC7uyZtfXO7A\nPA2sivMdON/ckywcOd2CU7Ud3b06dhw62YwDZY1o7jD5PS9Yg9tvno/Pv33XKgrUoHMnOy3t5oD5\nZf4tQO96VN2NybDtw4DgfB+r01vQYbTifHPPAsJcGutnfEqhe/TnLrzJ6t8QPVnTjsr6zqh6uDjP\nswnxMM8Q0MDeT98hpTUthqA9ssGSDb3JhpPn2/u4Pv2Hn3mcimDo+/cVzSFHcARLwj2/0YXolTp6\nttXvvVpapfNed4HVPU9FUMY9WAxchTqHXJO8LpMtZGn8/uRYgZ8fXHrpbcGGAgZ5XnmVDoYQPfCx\nwPn2fH5+PhYtWoSrrroKarUa33//fdyCIiRRqhv1+NtHpRDyebh7ZTFS1JK+nzQCCQU8/OxHRRAK\neHj9k3K0dg7MXVFCIjF+/Hj8/ve/x6xZs/Dmm29i06ZNsNsHwxo90XG5WHQarNCbI28cBDYEA5sd\ngQ2sBp8hfodONgddMLSvxldg2e/ARr+H3eFCpyH0MQW7Cx5NtbrAoVAuV0A1vF7jrkJvy/c8nqrp\nCDkHp7za3dsXbG4Uy7rLskeayFpsDpRW6VDdqEeZTyO9L6dqQzfaDWa7d1HqoDjswmx1oCuCBqzZ\n5kR9m9Hvd8EKdXxX0YTTPrE7XS6UVunQbrCiqb3v76KwSxNwOC5PYljV2IXm9p73R22LMdRToJBy\nn2bjef0C31dcC85Ect0GE+r8nKrt4JSYl1XpcMTnc8L/Rkjvbbd2mv2S69Dczw01LNnD9zOrL50m\nm3e4bjxwukX/2GOPYd++fcjK6qmsxjAMNm/eHLfACBloui4LXtx1FDa7E/+7fAqtBdWHTK0C/7Ng\nAt78rAJ/+6gU96+ZFrMyy4TEwqOPPoojR45g3LhxuOuuu/Dtt9/i+eefT3RYUatu7PIbeuXhab9E\nM6joh1MtkPuUrPY0+MMlLH1NQPe0k+0OV6+hWr5Y9BV777+arc6eCmMBMTa3m5GezL0S3pEz/r3y\ngYfcoDMhI1Xu/8sgSYBOb4FcErp5FSoJauu0oKm70T59AveqZ57hXZ7GZ7hz7CvUPDMAOHHOPYxR\nwOdB0I+ShFWNXXA0GiJKsiqD9NyFGi/Y1mXBeLiP2XeYak1z7+GUHjGrHskwYFnWe97Tuqsttust\n4Z4CILLCI+7rxP/x1Q1dOHG6GdMnasNOYwg1PJGLcOeJWzLUo8Ng9SsGEtiB63K551FGoq9rKmgv\nFvyLg/gSCnhwBhYsiRFOSdbXX3+Nzz//3LsIMSHDjdXuxEvvHUeHwYZV88Zh+sT4lvUcLi4vHo3y\n6nYcLGvC7q8q8eO54xIdEiFed911F370ox/BZrNh/vz5mD9/fqJDiolghR1Kq3TeRna4imye9pOn\nkRh4d9zFsn7tOk+jMNwE+FCNGg/P8B6DxY5jlaETMrPVEbaRHOwG+8ma9pAlsTuNtoiSrGB85y7V\nNOt7J1khGqQ1YXqBnE7/tbo8wxF957j4FtCI1Mnz3IaV2UP0mPkOh6tq7Ar6mMAev8AhXI06E1RK\nKac4wukrJeF6rIC7h6a10wxbmGuZkxALdQdNTqJI7ILto6qhCy6WhdHigEoWvIS6i/W/8eG7FYfT\n5VeAJrjIKgA6fBZ69p3jWFal65UQea4bg9kOHgNON2Y9l5bDyXrnU3J18nw7UtXSsKOS4lnSjNNt\n56ysrJguEkfIYMKyLN74tNxbqn3hzJG7FlakGIbBuoUTkaaR4rMD5zmXPSZkIKxatQp79uzBVVdd\nhYceeggHDx5MdEgxEWpOSiRaOszeqn9cRNMw5Tp/J3COmS+7w9lnO8QcpmcmeFzhma2O4L0rUTp8\nqhln+7h7z3XIXzC9yohHgGVZTvP5Kht6ki9PlblIIuZ8eH20gLkWXwDc85/O1HX6VRzsFRfHbUXa\nJg47RDGEcMNw27usIQvVNAYMufQ92XUtRtSGOf6Ah/f5t/Lqdhw+1ey9CeN70yZYj5PnPJw414Zj\nlW19vmdZlvUWwLE5nN75lFy1G6w4XdcR8edjrHDqyVKr1bjuuuswdepUv1LuTz31VNwCI2Sg/N+3\n1fiuvBnjMtW48Woq1R4pqViA9Usn44kth/Dax2V45OYZNJeNDApz587F3LlzYbFY8OWXX+KZZ55B\ne3s79u3bl+jQohJYCrs/PGtQ9YXfPVQsVIMvcP6Q3eEEn8cDz2eIWbheMF+hGq4OlwuHT7UgLSl0\nr1Sn0Ray16W/QsXdYbCiqkGPorEabxIQcXGBGK1BFq12vRUna9pRNDYZDocrbKIbiqEfa1kNVtWN\neozqq/eTYXoNewunw2D1DgONRLhErkFnDJlwmK3+iYvnrdvWaQk5ZM5vv9xD9CZSNofT25sVTuDn\niGeOYjAFAnavAAAgAElEQVR1rUacqu2ALAbFx8IOoY1jk49T5Jdffjkuv/zy+EVBSIL8cKoFu/dX\nIkUlxp3Lp3D6kCC95YxSYvX88dj6r1N4+f3jeODGaRAL+YkOixCcOXMGn3zyCT7//HOMHj0a69at\nS3RIiRXhvC2JyP0+DnYn/kBZIwQB80IOn2qBRCjAheNTvb8LXP8pFFMf6+6FKrBTcqY16BwUtTz4\ncCrA3VPUV2/EuRBJ28nzHWDBciqyMNh5FuttbDNBF2ZOEScRtM45PzSG84O47Y7ltAhwsJ7GYKEa\nLXa0nh/Y6yQwtpM17ZhRkIbTddyGVoZ7X4TqlXa6WM7rZnItg+8Z9tzX5wIXod7LQPRrsYXDKcla\nvnw5amtrcebMGVx22WVoaGjwK4IRjMvlwqOPPoqTJ09CJBJh48aNyMnJ8f59586d2L59OwQCAdav\nX4958+aho6MDCxcuxIQJEwD0rM8V7LGERKum2YDXPi6DSMjDXSuKoQrzhUz6Nm9qBqob9fjqWAPe\n+qwCty8pol5BklBLliwBn8/H0qVL8dZbbyEtLS3qbR49ehTPPfcctmzZgurqajzwwANgGAbjx4/H\nI488At4gX1Mv0mFoZ+o7kZokDd2TFaTBZLHHp1pXqNhDTfLnhynYcL5JH1ERAo+QvXJDdkqFO+5I\nEqzBMn3EU5wj1uzOvo+P6ykINczWZnfih9Phlz/p71kOPpyYe0Ia7thCjWCsqG7n9NnCsiyqG8MP\nVxxOOCVZn376KV599VVYLBZs374dq1evxm9+8xssXbo05HP27NkDm82GHTt2oKSkBE8//TReffVV\nAEBLSwu2bNmC9957D1arFWvWrMHs2bNRVlaGxYsX43e/+513O6Ee6ztskZBIdZls+POuY7Danfjf\nZZORna5MdEhDHsMwuPHqiahvM+JAWROy0hVYNCun7ycSEifPPfccJk6cGLPtvfbaa/joo48glbon\n9D/11FO45557MGvWLGzYsAF79+7FggULYra/eDhb14mCHE3EDeX+zBEqOdOKrAgad7HGAmhoM0Kn\nt4IfpDDDmBR58CeGcfRMq/dufl2rwXsXPNL5YIkSrmw7VwfLe69BZrU5ua+dhfBVAH31Z5hdtDqC\nFJYJ5Ftc5lxDF7LSFEGzolBpfFtXlL2GYQS78dEcpte1d8X30K9jqDmKXD8fdHrroFvyJZ73gjnd\ncnvttdewbds2yOVypKSkYPfu3di0aVPY5xw+fNg7xPDCCy/EiRMnvH87duyYd36XUqlEdnY2Kioq\ncOLECZSWluLGG2/E3Xffjebm5pCPJaS/HE4X/rL7BNq6LFh6WS4uKoj+7jZxEwp4+MXyKdAoxdi1\n7ywVwiAJFcsECwCys7Px0ksveX8uLS3FzJkzAQBz5szBN998E9P9xYPd6cLxyrZe6xGFw7Js2IqF\nwTS3m2CxOdAUwZo1scYCqG7SQ2+yBV1jyFOCOxKBDVhPg7SvCouDVSyGYgHu8veRVn7jIp7JSCh9\nFVew2Zw465NsNLWb0Nhm4tyDa7Y6OBWdiWWHYbieynCFQGJtsCVY8capJ4vH40GhUHh/TktL63NI\nhMFg8HsOn8+Hw+GAQCCAwWCAUtnTcyCXy2EwGJCXl4fJkyfj0ksvxUcffYSNGzdi/vz5QR9LSH+w\nLIu3/30Kp2o6cNFELZbMHpvokIadJIUYd14/BU9t/QF//bAUD66d3rvsMSFD0MKFC1FbW+v92bd0\ntVwuh17fd8NJo5FBIIhyvmJNZ0zKY4fju321Rg5TbVdE+2w12L2Pj3esIWPQ2yLad6Li7I/BGKsT\nveMajHEGE0mccqUEqoDiEsF+F8q5ZiPn/Zns7Ig4p4mm1cZnNBOnJGv8+PHYunUrHA4HysvL8c47\n76CgoCDscxQKBYzGnrsBLpcLAoEg6N+MRiOUSiWKi4u9wzAWLFiAP//5z1i6dGnQxxLSH1/8UIf/\nlNQjO02Bn15X1K8x+aRvuaNVuOXaArz2cRn+tPMoHlo3HUkKcaLDIiSmfG82Go1GqFR9L2DeHuXw\nJ8+wnC79wN0R/tc3lf1+rkopHdBY+2uoxAkMnViHa5yM09mrPPlAHedwPaeJ1tISXS9sqCSN03DB\nDRs2oKmpCWKxGA8++CAUCgUeeeSRsM+ZNm0a9u/fDwAoKSnxFrMAgOLiYhw+fBhWqxV6vR5nz57F\nhAkT8PDDD+Of//wnAODbb7/FpEmTQj6WkEiVVemwbc9pqGRC3LWiGGIRVb+Lp0smjcLyOXlo67Lg\nxXePRbUCPSH9UVdXh1tuuQVXX301mpubsW7dOr+eqGgVFRV5197av38/LrroophtO6TBUXOAkBGr\nM0FrLpGhh1NPlkwmw3333Yf77ruP84YXLFiAr7/+GqtXrwbLsnjyySfxxhtvIDs7G/Pnz8fatWux\nZs0asCyLe++9F2KxGPfddx8efPBBbNu2DVKpFBs3boRWqw36WEIi0agz4S+7T4DHA35x/RRax2mA\nLL4kB22dFuw/Wo9XPyjF3SungD/Iq6+R4WPDhg346U9/iueffx5arRaLFy/G/fffj7fffjsm27//\n/vvxu9/9Dn/84x+Rl5eHhQsXxmS74USzSC0hhJCBw7AcSgwVFBT0KsWs1Wq9PVWDTbTdfmR4MVrs\n2Lj5MJp0Jvz0ukLMnjI60SGNKE6XC3/edRzHK9sw54LRuOma3p8nhITT3/Hy119/Pd5//30sW7YM\nH3zwAQBg6dKl+PDDD2MZXkSi/X6y2ByobDIOmaE4Q2XY0FCJExg6sVKcsTdUYh0qcQLAnOnZsJmj\n650M9R3FqSfLt5qf3W7Hnj17UFJSElVAhAwETyXBJp0Ji2ZlU4KVAHweD+uXTcIzbx/B/qMNSFFL\nseTSsYkOi4wAEokEjY2N3qT+0KFDQ375j1Dr7hBCCBlcIh63IxQKsWjRIhw4cCAe8RASU9v2nkZ5\ndTsuHJeKFVfkJzqcEUsiEuCeHxcjRSXB7v2V+O+xhkSHREaABx54AHfccQeqqqqwdOlS/OpXv8JD\nDz2U6LCiwoZaDZQQQsigwqknyzPMAnCXrD19+jSEQmHcgiIkFvYersW+H+qQqZXj9iVF4PFoiFoi\nqRVi3LvqAjy19TDe+KwcEhGf1igjcVVcXIxdu3ahqqoKTqcTeXl5Q74nazDmWKOT5X2uLUQIISMN\npyTLUz3JQ6PR4IUXXohLQITEQum5nkqCd68shlTM6VIncTYmVY7/d8OFeHbbEfzto1KIRXxMyUtJ\ndFhkmPntb38b9u9PPfXUAEUSexymUQ+4ZJWYkqwRhscwVISFkD5wankO5S8kMvJUN+rxyu7j4PGA\nO68vRqp66CyINxLkjlbhlyuL8cedR/HK+8fx/264EBOykhIdFhlGZs6cmegQ4oZLw3bsKBWqGrsG\nIBo3uUSImQXpaO4wD+h+SeIMlfwqWSmBTm9JdBhkhOKUZF155ZVBq4F5Vrvfu3dvzAMjpD+a2014\nYWcJrDYn7lg6CeMy1YkOiQQxMVuDXyyfjJfeO44/vXsUv1o9FXlj+l7IlRAuli9f7v13eXk5Dhw4\nAD6fj9mzZyM/f2jPzXRyGC8o4A/w0GgGCVvYXSYWYEyKHGfqOxOy/5GKHaAF24R8PuxOZ7+fzx/o\n90IflFIR9FFWsutL4E0WhUQIg8Ue132S4DgVvliyZAmWL1+Obdu24d1338W6deswdepUbNmyBZs3\nb453jIRw0mGw4vkdJegy2fE/V0/AzML0RIdEwijOT8XtS4pgtTvx3PYjOFNLjSQSW6+//jp++ctf\norm5GbW1tVi/fj3ee++9RIcVFYN58DWWmID/DyShgI/UpNCjFdI1spju74L8VO+/k5W03mIsXTQx\nDQqJ/3z/eOfuaUmxvT76kpWuwMyC+LZNRiWHPqac9N6lxpmEvHMHlkKamDoSnJKsr776CnfeeSfS\n0tKQnJyMm266CZWVlcjIyEBGRka8YySkT51GG57ddgQtHRb8aPZYXDktM9EhEQ5mFqbjjh9Ngs3u\nwvM7SnDyfHuiQyLDyI4dO/D+++/j/vvvx4MPPoh3330X//jHPxIdVlSytArIJOEHoQx0o2kwrHs3\nbkzwUQtJCnFM9+M7vzc7XRHTbffX9AnDo4CQgM9DXsDrGO2V1dd7ITVJglmTR0W1D669uKkqKZRS\nYUKLcAn4vZv9vvHI+pi/7ns+ZWJB0BsNeWPUUCuGdoGhWOFcwv2bb77x/nvfvn2Qy+VxCYiQSHV1\nJ1gNbSYsnJmFpZflJjokEoGZhen43+WT4XC68MLOoyit0iU6JDJMqNVqCAQ9jQaZTDbkv7tEQj4K\nxyb3+v34jOjmNY4dNbSH66YmSSHgRbwqTVTCzUsK1ZvAj0OMA3zYcSUS+h9MvBN4IZ8X9T7GZ/b9\n3ssfo8a4THVMj0ebJMXo5PCfZ4HvawGfB7GQ7/c7BkBBtgZquRiTc1M4J43aJCmEgt4Xn2cqUazl\npCsxKchnHxdqeWxvtnDF6a35+OOP44knnsCsWbMwa9YsbNq0CU888US8YyOkT+16K57ddgT1rUYs\nuCgLq+aNGxR3VUlkpk3Q4q4VU+BigRffPYYjp1sSHRIZBrKysnDDDTdg06ZNeP3117Fu3TooFAq8\n/PLLePnllxMdXr8pZL3vEitl/R8Oo5aJwg4x0oYZjudLJY/93WuJMLrKsNF8HeRGkXiGSlqnT9T2\ne5uRkIm5XQ/Bho8NBAGPB2WQ69i3pyV3lArjMtVRXwOhpGtkkIoFUfWWycQCJClESFFJoPUpsuWb\n/KSqpdAoY9fIz0pTQiISYOwoJTLT5EhWSkIOiw0cJsfnM0Ff8ySFGIU5GvB4DPhhetoC5+IFmyPa\n1/ksyNZgVj+mczAME/Sa4SJDm5iba5ySrMmTJ+OTTz7BZ599hi+++ALbtm1DdnZ2vGMjJKzaFgM2\nbj6Euu4Ea/V8SrCGsuL8VPxyZTF4PODl949j/9H6RIdEhrjc3FwsWLAANpsNJpMJs2fPxvTp0xMd\nVlyIfO5OR1qUIFjPmK9RyTIUZGtQkK3pdZfbd+6rVCzARRPTcHHRKEzJS4l42KLvfCePUMOOlNLe\nv89Ki+3wvTRN7CvT+p4/sYAf5pHc9bfSX7gkcoJP70xhtqZ/OwhjdIqszwYoj8dAIRViSn7/ei+A\n8El2uMQn3E0HX2NS5GAYBuMzk/y2l+UzlHRchrrXMD2uyULgHDUAyEiV48JxqeDzeODzeJiQlRRy\nzlFgNVIBj9drCG1ggbBwPVmqgLg9S0r49iLzGKbPGy6xaKsJ+dzeP3wez++Y4nE9h8Lp9kBdXR0e\nfvhh1NXV4e2338b69evx5JNPIjOT5r2QxCiv0uHl3Sdgtjqwcm4+Fs3KpgRrGJiUm4xf/2QqXnz3\nGN78rAKdRhsWX5JDry3plzvvvDPRISRELMtrF2ZrIPdp6F04LhXHK3Xeim+BDTJPY1IuEWJWUTrK\naqIvaDN1nBZHzvTdu52mkeJcd1U1XnfDN5JPDo1CjHaD1ftzrD93cke7k5rCbA3AMKhvNcLq6H/l\nPE95ct+eh56S5X1fBAwT+mGBw8DUMhE6TbGrihfJuY00WRcL+JzOq7dgS5DNp2tkaNSZQj5XKRO5\ne376eY0U5WjAssC5xi44nSwytQocq2wF4E7wPPtOUoo5VQYMjEImFsBkdUAU8Dry+UyvOWGBSRcT\npidrfGYSDp9q9v7s6u7JEgh4cNhcANzJcc4oFU5XtcFqd78OxXmp3uMb6G9zTR/zMiViAWxxqvjI\nqSdrw4YN+OlPfwqZTIbU1FQsXrwY999/f1wCIiQclmXx+cHzeH7HUdjsTty+pAjXXkyN8OEkf4wa\nv71xGlJUEuzeX4m3/33K+0FOSCTeeustzJw5E4WFhSgsLERBQQEKCwsTHVbchRvuE6jYZzHwCUHm\nlqgDGigiIR/5GdyH0Xl62GRiAWRiISbn9uwvcIhTqI/xYPM+gvH9Hpg2QevuWeD43TAxRwOJKLJh\naaGSWU+vwoyCNFw0sacohaT7XKgVYqgjHFrpic337v34TDVmFaaDYRgU5SRjYlbPHXqWBaaN1/Z5\ntz/UR2uo71SFRIhkpQRFOeF7l7K0PT05KpkIConQLyFhGPT52vT3W310KsehYd37789NCT7DhE2w\neAyD3NEqFIXoJWYYd7KTP0aNCVlJfsVsfIeajknhdiyenqPMVPd5LxqbjCl5Kb2uaWGQwhfBYg8l\n8L3ouX58l41IUojB4zHe5E3A40EmEWB8RhLUMhGU/RxWHMuF2PO6b3hMzNL0mqMWS5w+udrb23HZ\nZZcBcF8Yq1atgsFgiFtQhARjtjrw6oel2LnvDJQyIX79k6m4ZFJ0VYHI4DQ6RY4H105HplaOL36o\nw8vvH4fF5kh0WGSIeeutt/DBBx+gvLwc5eXlqKioQHl5eaLDijsBn4dUVe+hbklysV+jbXSyHDKf\nXqp43KzybFIuFaI4PwUKqRCZWgUmZCYhd7TKLzGQiATISFX0fr5PWGNS5FD3c16Gh1jADzoMK5hp\n47WYPiENMrHA26OXmaqAkM+DRBy8ceaZnM/n8SDg85A7SgWxkA9FH/Pmws094vMYXFw0ym8YGsMw\n3tdMJRdBoxT75S0iIR+yEDF6hJrL59fQ9k2OeAwmZCX1Gg4WGLvvUEulTITJeSl+Q1qBCJKoEA/0\nJLO+xRCm5KX4DfUL1y73HJZIyEdWmjKi4jG5HNZ1TNfIeg2vixjHkyQS8jGrMB2Z3UNmBXye93qV\nd58nbZKUU2XDSG7SeIYj8n2uEc8+Aj9OUtQSFI5NjnpNvYifH+ThaRoZZhWmx3SuXDCcbttIJBI0\nNjZ638yHDh2CSETlGcnAOV7Zhjc/q0C73orxmWqsXzY55qV5yeCiUYrxwP9Mwyu7T6DkTCue3voD\n7l5ZjGQVrU1DuMnPz0dqau95PrHkcrnw6KOP4uTJkxCJRNi4cSNycnLiuk8upBIB0OX/u4IcDQxm\nO+rbjACAnFH+E+C5tq1UMhHkEiHneSuBMn16OQLvjGelKSAU8LyLqTIBg8Wyuyftl54LX4XU0w4L\ndkhajRQM4DcMK1TVP09iUOwzXywzTeFtzAbft/9e05NlSA9yrgJjy89QodNoQ7JSjGOVbUG37XS6\nQu7Xvc3u3pnun/u696+SiVA0NhllAVVdA9uxnu2EukR8T5+7Z6D3I31/w2OYkB1ZarkYnUart3cn\n1D4LsjUwWRx+85siSRB8G+sZ3b1fp+uCP9Z3UeSLiwbu5i4Dd3LU0mHu+7EhTmh2ugJyqRCpHL87\nI6kwz3Z3ZQXbdyyHLfvi8xi4nCz4fAZ2DqNtAyPzXssDMAKKU5L129/+FnfccQfOnz+PpUuXorOz\nEy+++GLY5/T1xbNz505s374dAoEA69evx7x581BfX48HH3wQTqcTLMvi8ccfR15eHt588028++67\nSE5236147LHHkJeXF8Vhk6HCaLFj+97T+Pp4I/g8Bj+aPRaLLx0bdK0HMvzIJELcu+oCvPPvU/iy\npB6/33wId68o9s5tICSctWvXYsmSJbjgggvA9xk29dRTT8VsH3v27IHNZsOOHTtQUlKCp59+Gq++\n+mrMth9KX3NkQjUf5BIB0pJkSFb1vknFtdHB4zGY4jPMMBoKqRA56Uq/Est+jXEeEzQukZAHmEMP\nJQw3j4dl/QsCZKQqkJokQXVtz2MinRwv4PHgcIVPgPoikwiglIlgC9Jy9IQbrJpbUN0P43LXXxJk\nuFQ0azkV5mj6HqUZ5u8Ts5NgtTn91iQLRsDneXvUJmZp0GGweofHpWtkMPYxlylcjL5/GzdGDR6P\nwanajrDbi1ZOuhIWm/u1nzpeC7vDBYZxDynkkmSFwufxkMaxQigA8CNoX3neRzweg+K8VL9k23PN\nRnMtZaYqUNsafOScgOt2EzibhFOS1dbWhl27dqGqqgpOpxN5eXl99mSF++JpaWnBli1b8N5778Fq\ntWLNmjWYPXs2XnzxRdx444246qqr8NVXX+GPf/wjXn75ZZw4cQLPPPMMJk+eHP0RkyGj5HQr3vpn\nBToNNmSnK3DrtYXeu5hk5BDweVi7cCJGpcixY+9pPPP2D7htcREuKhgeC3CS+HniiSewZMkSZGRk\nxG0fhw8fxuWXXw4AuPDCC3HixIm47cvXhOwkGMwOlFcH79HxvaOfqpYiuXtYDMMwyAs11CkOjREu\nd7NHh5l3EqohPHaUEmIhP/RzmYD/wz08skFnhFouQk1LT8MtK03hl8gJeLxec9FCkYmFMFnt0KjE\nUTWE3aH2HmallImg90mmk5USdFmMoW80BZyvsaOVON/EIE0jRUWIxd6DNYJ9kzOlTNhnj5jJ6h7O\nnaKSQCISwOmTcEpFfO+xmLuHfTNwD9lqN1ihlIn85qjxGMYvwQpMsieNTe7V86hRiv2GfnnOT2W9\nf3duT2GQCApqMD2JR7i12KLtGfG9lsVCflznCoWTk66E0+mCSMQPek1naRWoaTEgSSFGk879dx7D\n9Fok3ZuARXFexmjl0OmtMFntvc6HVCKAw8UiWSn29s4Hk6r2TzDj1cMWDKck69lnn8XcuXMxfvx4\nzhsO98Vz7NgxTJ06FSKRCCKRCNnZ2aioqMD9998PpdLdiHY6nRCL3W+Y0tJSbNq0CS0tLZg7dy7u\nuOMOznGQocdgtuOdf5/CgbImCPgMls/Jw6JZ2dR7NYIxDIOrZ2QhTSPF3z4qxV8+OIHr5+ThOqo8\nSMIQiURxrzBoMBigUPQMHePz+XA4HH6LIPvSaGQQxKB096h0NViWRV13I0erVUKldFfyS06RQyUX\nQyBxz38Ktq5WMGKjDaq2nkbVBeO1SIpyzsKZRgNUSimSNXJotdxuktnAoM3o7oVITVH4HZvvNsaM\n7j2HxvO4NK0SDMNAZLB6j+miKWNgszshEvJhBwNGZ0JWutK7zaQkGYx2FwQCHudYr0iWw2x1QMDn\n4bvSRozPSoJWy62UfGOnFS6fRrtWqwSPx8Bmd0LV4E4CU5OkYPhmyKVCb0zZmUkhP/faTHbYXIBE\nxPc+PivD3StX327xe2xKigLa7ob9ZTIxxCI+vi9rAgCkpymhatB3/1uFhg4LGD4fSUqxd7uec+0r\nNcX9OjtdLE7W6aFSSlEwTguGYZCcLMd/u5fmSElRYFSKHOPGpnDq6fDdV14O917UdrMDFqe7VZ2k\nFCM7XYljZ9xV7rRahXdOYuAxaVOVUDUZvbGmqqUQSYRITZL6zWP0xQgFUHVY/LYXK564Irk2uWwP\nCB5rZkYSXC4Wze0mnKxu93usVqvEBd0LDte1WyCyO5GcLOu1HbVaCpvL3UMbbB/Brh9fPB6D9DQV\nNBo5dF0WpCfLwDAMUltNMJrtSNcqcfEYNewOJwzHGgAAhbnJaOuwoLndBLGIjxlFo7w3nDz702jk\n0AYszxDr18uDU5KVlZWF3/72t7jgggsgkfSM6Vy2bFnI54T74jEYDN5kCgDkcjkMBoN3OGBlZSWe\neeYZvPLKKwCA6667DmvWrIFCocCdd96Jffv2Yd68eZEdKRkSDlU0Y+u/TqLLZEfuaCVuvbYQGRy/\nsMjwd+G4VDx443S8uOso3t9fidoWA265tjBhd/zI4HbppZfi6aefxpw5cyAU9jSMZsyYEbN9KBQK\nGI09d1FdLlfIBAsA2ttDl4XmSqtVoqXF3QAeo5FALOSjpUWPLr07mWhrM8JhsSNFJoTZaIXZaA23\nOa8uk827DQCwW2xosURf2rhLb4aEz6BFyq16n05n8sbhssnQ0qKHxWyDQir0Hne4fQFAa/cQI73P\nMfk+Vy3hwy4XQi5k0NKih1arhK7diC69GUI+r8/9BFOUpQbAcn5uZ6cJXT4l41tb9WAYBnaHyxuz\nSsJHl94MMb8nft/XP1BHh/vcWQX8Xo/xfW0BoE0nAt+nx8lhtfudP9/z1tFpht5kA+N0erfru72s\nNCVqmvXIH6VAS4ve24vRpTd7XwsAYJ1O6E02mI0ytEQwvNJ3X5G8Nu3tPddSUZYabTqj3zFKxQK/\n8+n9W1vP8bfrjGAcTkj5DIx6C4x6S5A9Ae16a9BrLRY8271sWlZMts31fPL7eGxnpwl2pwtSAeP3\nN61Wia4uM7r0ZjBOUdB9BF6PgWYVpnufx0fPezpdJUKD3QEJzx2P7/uFcTgh4bu3PT4jCbq2nmuv\n53XXA46eQlrh3k9chUrSwn7iNTU1IT09HRqN+y7I0aNH/f4eLskK98UT+Dej0ehNug4cOIDHHnsM\nf/jDH5CXlweWZXHTTTd5/37FFVegrKyMkqxhpstow9Z/n8KhimYI+Dz8eF4+rp6RFXIyMhm5stIU\n+N1NM/DK7uP4rrwZTToz7loxhQpikF7KysoAuEdDeDAMg82bN8dsH9OmTcO+fftw7bXXoqSkBBMm\nTIjZtrmIZQEgNg5LJeSMVqJVZ0Cqun/vT0+vwbQJ2liGBQGfF6ZwR2J6x4P1To1JlUMi5EMTZA5d\ntMJ2IIWuFhD04Rmpcm/xiGBP9yjIToLZ6gy5eG5fol3AWehTajyauULBRFJ0o780SglaOKybFUsT\nMpNCzkdTykTQ6S2QS3qnE9lpCvAYBmO4ltQPEKq3ViIShJ2XrZAKwxYniW7mZGTCJlk///nPsXv3\nbjz11FN4/fXXceutt3LecLgvnuLiYvzpT3+C1WqFzWbD2bNnMWHCBBw4cABPPPEE/v73v3vH0BsM\nBixevBiffvopZDIZDh48iBUrVvTzcMlgdOR0C978rAJ6kx3jMtS45dqCsGP0CVHLRfjNT6Ziyz9P\n4qtjDXj8ze/xi+unYHyQdX7IyLVly5a472PBggX4+uuvsXr1arAsiyeffDLu+4yXeAy9HZOqAL8w\nPeqyzf0VyTF55mokegSy7/55DIPUCIoWBFYXDCfcaxLtKWAYBoW5yb16ffg8HhTSyG+eatVStHSa\nI7n86fAAABmmSURBVA6M7T4Tou7kTCYRIn+MGgwQdgREf45f3D33LNLFkwe7ZJUE0ydog76X8sao\nkGqUBC2FLhTwB6RIVaTv11iut9WXsEmWbyAff/xxRElWsC+eN954A9nZ2Zg/fz7Wrl2LNWvWgGVZ\n3HvvvRCLxXjyySdht9vxwAMPAAByc3Px+OOP495778W6desgEolwySWX4Iorrujn4ZLBxGx1YMcX\np7H/aAMEfB5uuHIcFlyUFfO7S2R4EvB5uHlRAbLSFNi+9wz+8M4RrF04EXMuGJPo0MggcejQIfzj\nH/+AyWQCy7JwuVyor6/HF198EbN98Hg8PP744zHbXjSUUhH0ZlvQanGcni8TYkyKPOwk8v5IVIIV\nKU+DfGhEG5znVHNpRwZrNCtlIm9VO18sp7TNX5pGhhYHhxrbHGiUYrR0mpGmiWzZAM+157ugvTaC\npBXgnqiLhXxMzNJA2sfaZEORMEQPooDP6/coknSNDB0GK6xc6rDHgKdSYdRrl0UgbJLle2FFmvkF\n++LJz8/3/nvVqlVYtWqV398/+uijoNtatmxZ2KGJZOg5U9uJ1/6vFC0dFmSlKXD7kiK/tVMI4YJh\nGFx1URbGpMrx6gcn8OZnFahpNmD1/HE01JTg4Ycfxu23347du3dj7dq12L9/P4qKihIdVtwUjtXA\nZnd676hHimEYZKcrMSpZxr1U+CAX0T27QXLI8brTLhMLYLO7vKXmg50b34V9falkIhjMdqj6WFA5\nXjy9KaEa+6EIuwtmRVpev7/3BeK1uO2UvBS/RHE48PRyHShr9Pu9Wi7GmJT+rcEXTmaaAmO08gG9\n6cNtFioGZtEuMvw5nC58+N9z+PRANcAC116cg6WX5YZc64QQLorGJuN3N12El947jr2Ha1HbbMDP\nl07iXIaZDE8SiQQrVqxAXV0dVCoVNm7ciOuvvz7RYcUNj2G86wRFQ5TgQjLRtDYyUxUQCnu+TyIZ\nGTGQi5SGIxTwkaKSRPn51btB7lnb7IdTrbA7nRG9zplpCqjlIu+6VAAwKlmGRp1pwBqtkSZYgDs5\nq201ICut/9XjBsPgGnmIiobDSbpGhqZ2E7LSFP2es9eXge5VD/tpfPr0acyfPx+AuwiG599sd+nG\nvXv3xj9CMmzUtxrx2sdlqG7SI1UtwW2LizAhi+bQkNhI08jw4NrpeP2Tchw+1YJH3/we/7tsMs3T\nGsHEYjE6OjqQm5uLo0eP4pJLLoHJFH11PzJ4Zab5j4iIqEfbm2XFLh4uZGIhJuf59yD193PLU3wh\n2IKynuSxOD8ZJoujz8V+ffEYplfSNxSGgcokAswsSI9qGkKik+7hzpNc5aQrkTNKOSSuK67CvsP+\n+c9/DlQcZBhzsSz2Hq7Fri/Pwu5w4bLi0fjJ/PERfcATwoVULMD/Lp+Mz787j11fnsUf3jmCVVeO\nw1XTM+mLcgS6+eabce+99+Kll17CypUr8fHHH9Oi9iMMv7uaHJdlHgZ6TpbvZ1KsGpZjUuWwO10Y\nE6Z4lFDAh1rR9/mYnJsStmJeT6GQwf3Z2r8Eq+c5g/zwhrzc0aoBKZCRCGFbuZ4Kf4T0V2unGW9+\nVoGyqnYopEL8bMkkTJ8Y21K8hPhiGAaLZuUgd5QKf/3wBLbtOY2K6nbctKhgQCe8ksRbtGgRrrnm\nGjAMg/fffx9VVVUoKChIdFhkAPEYBheOS+VUXtuT6Azlhe8FfB7yx6hjsi2uQ7aGZw7SM9xysCeR\n0YjF8GISGp1dEhculsWXR+rw7pdnYbU5UZyfglsWFdAcGTJgCnI0eOSWmXjt41IcOd2Ks/Xf4dZr\nC1Ccn5ro0MgA2LdvH8aNG4esrCzs2bMHu3btQmFhISZMmAAeFUUZUbg2JLPTFWBZIDONlhDhoj8V\nB4cK3xoTg2FOVrzkDYMepAGsyB4x+qYhMdfUbsIf3jmCrf86BQGPwW2LC/HLlcWUYJEBp1GK8auf\nTMWqeeNgstjxp3ePYfM/T8Jic/T9ZDJk/eMf/8DLL78Mq9WKiooK/OpXv8L8+fNhMpnwzDPPJDo8\n0odEtZmEAj7GZaoTcHd/ELcSR5hxY9RIUUn8hpcO556s4bBkzmB+eagni8SM2erAJ99W41/f18Dh\ndGHq+FSsXTgRSZRckQTiMQyumZWNSbnJ2PRxKb48UoejZ1qx5qrxmBZigUUytH344YfYsWMHpFIp\nnnvuOVx55ZX48Y9/DJZlce211yY6PEKGhcGyeHMspSZJey3+PJyObzgS8HkozEnu9/qA8URJFoma\n3eHE/qMN+Pjrc+gy2aFRinHDleMwoyCNGrBk0MhKU2DDTRfh42+q8dmBaryy+wSm5KVg5dx8ZKXR\nGm3DCcMwkErdDaWDBw9izZo13t8TMlh45oAN/etyqMcf3KhkGTr0Noj6UTqeDCy1fHDOt6Yki/Sb\n2erAV0fr8dl359FpsEEk5GHZ5blYODObUyUnQgaaUMDH9XPycMmkdGz91ykcr2zDico2zCxKx49m\nj8XoMBW5yNDB5/PR1dUFk8mE8vJyzJ49GwBQV1cHgYC+9ga74dlk7y0rXQGH04Xs9P6v4TQYDIMR\nZ0GNHaUCRiU6iviK18LXxI2+bUhEWJZFVaMe/ympx8GyJljtTohFfCy6OBsLZ2T7LVRIyGA1OkWO\nX62+ECfO6fDef87iYFkTDpY1oTg/BQsuykLhWM2wWqtjpPnZz36GZcuWweFwYOXKlUhLS8Onn36K\nF154Ab/4xS8SHR4hANxl5QtyNIkOo98ytXJYbU5kpdNIgKEmI1WBhjYjZBJKA+KJzi7hpK7FgIPl\nzfi+vAlN7WYAQIpKgusuycHcqRlxW52bkHhhGAZT8lIwKTcZR0614J/f1+DY2TYcO9uGVLUEl04e\nhUsmjUJ6sizRoZIIXXPNNZg6dSra29u9Jdvlcjk2btyIWbNmJTg6QoYHoWBoJ4kjWVaagobJDwBK\nskhQDqcLZ2o7cayyDUfPtKKhzQQAEAl5mFGQhksmj0JxXsqwqExDRjYew2D6xDRMn5iGyvou7Puh\nFodOtuCjr6vw0ddVGJ0iwwXjUjFpbDLyxqhoEe0hIj09Henp6d6fr7jiigRGQwghZKSh1gLx6jRY\ncbxSh2NnW1FapYPZ6gQAiAQ8TB2fillF6bggPxViEc23IsNT3hgV8sYU4X+uduDwyRYcPtmCsiod\nPj94Hp8fPA+GATK1CozLUCMrXYExKXKMSpHRIseEEEII8UNJ1ghmMNtx8nwHTp5vR8X5DtS2GLx/\nS1VLcOmk0Sgel4KJWUkQUSELMoJIRALMnjIas6eMhs3uRMX5dpyq6cSZuk6ca+hCTbPB7/EKqRAp\nKgnUChFUchHU3f+p5CKoZCIoZUIoZSIopELq/R3i/v3vf+Pzzz/H888/DwAoKSnBE088AT6fj8su\nuwx33nlngiMkhBAyGFCSNUK4XCwa2oyoatTjXEMXTtV0+iVVQgEPhTkaFOenoDg/BaOSZcOgrCwh\n0RMJ+SjOT0VxfioA91DammYD6lqMaNAZ0dhmQkObCQ06I6qb9GG3xQCQS4XepEstFyFJIUaSUgSN\nQowkhRgapfv/1GM8+GzcuBH//e9/UVhY6P3dI488gpdeeglZWVn42c9+hrKyMhQVFSUwyqGPvnoI\nIcMBJVnDiMvFwmC2o9NoQ0uHGY06ExrbTGjUmVDTYoDV5vQ+1pNUTcxOQkG2BrmjVRAKeAmMnpCh\nQcDnIXe0CrmjVb3+ZrY60GW0ocNgRafRhk6jDXqTHQaT+/96kw16sx1dRpt3nmMoUrEASQqRN+nq\nScBESFKKoVGIoZKLIODT+3agTJs2DVdddRV27NgBADAYDLDZbMjOzgYAXHbZZfjmm28oyYpSqlqK\n/9/e3cc0dfZ9AP+W05aXttyCwD0WAXWTxC0Dh2y6iJiFEDbGIGPy6kv24CYSp3MbDGXRYUDAx8n+\nEHESzWLQDRC2LE+y4LZkQpjMGJhjoLg3xIVtTBQfaHlpOed6/qicSS1QnrvtOcXfJzbtOT2YL7+e\nXlev04tzhkZMeIhOOkMIcWEOG2QJgoDCwkJcu3YNarUaxcXFCAkJEZ+vq6tDTU0NlEolcnJy8Oyz\nz+L27dvIzc3F2NgYAgICUFpaCk9PT6vbzgeMMZgmBIyZeIwbzbfJx2NGHsYJHkYTj3GTAKPJvDxu\nFKasHzdOYHjEJH6YE6xc88BNoUDgQi8sDtRh8UPmD4dBAVoaVBFiZ57uSni6K206IyEvCBgymHBH\nP26+DY9jUG/EnWHz8uDddTMNxhQAtF4quKs4qFUc1Eo3821yWXV3WclBpXKDinODanJZ6TblZsu6\nB+W09mfPnsWpU6emrCspKUF8fDwuXrwortPr9dBq/zlDl0ajwe+//z7j/+3j4wWlHS5u6u/vOtdW\n+v9k/fe/7z+I4WjzvaZSoJz25ypZXSUn4LisDhtkff311zAajaitrcXly5dRVlaGY8eOAQBu3ryJ\n6upqNDQ0YHx8HJmZmVizZg0qKyuRkJCA5ORkVFVVoba2Fi+88ILVbdVqx/2h+f/qx9Hz1zAYY2AM\n4r1gZXlyHc8L5kGPiRdv9w6Q7h1AjZt4jBknMGbkYY/rwLmrOfxLo4a/j6f4dyB+3h54yNcLDy30\ngv8CTzraTYjMcG5u8NGZv52aidHE447hnsHXlHvzt2VGE4+hu/fGCcGBmRVQi4O1+wdgCsXdC8kq\nFOIFSn107viv+OUu1QalpKQgJSVl1u20Wi0MBoO4bDAY4O098+BgcHDmbzBt4e+vw82bM09NlQtX\nyeoqOQHXyUo57c9VsrpKTsA+WacbpDlskNXW1oa1a9cCAFasWIHOzk7xuY6ODjz55JNQq9VQq9UI\nDg5Gd3c32trakJ2dDQCIjo5GeXk5goKCrG4bFhbmqOio+p8ruNo7aPf/113FwV3NwUPFQfsvT/Hx\nvffiYxUHtZqDu9J8RNpd5Xb33nx0evLItbvK/GGHEDI/qVUcAhZ4ImCBp03bC4xhYkKAcUIQB12T\n9ybxxouPjfeu4y3XzfwzI2MmmHgBRpNgPgAFhrv/RDovFcaMPLSerjPIspVWq4VKpcKNGzcQFBSE\nlpYWOvEFIYQQAA4cZFlOo+A4DhMTE1AqldDr9dDp/hn1aTQa6PX6Kes1Gg2Gh4en3XYm/+nXfv+9\nM/o/+nlCCCEPhv379yM3Nxc8zyMqKgrh4eEzbm+vaSk0Fcf+XCUn4DpZKaf9uUpWV8kJuOB0Qctp\nFIIgQKlUWn3OYDBAp9OJ6z08PMRpF9NtSwghhDjbqlWrsGrVKnF5xYoVqKurkzARIYQQOXLY/I2I\niAg0NzcDMF9HJDQ0VHwuLCwMbW1tGB8fx/DwMH799VeEhoYiIiICTU1NAIDm5masXLly2m0JIYQQ\nQgghRI4UjNnj1Av3mzy74E8//QTGGEpKStDc3Izg4GDExMSgrq4OtbW1YIwhOzsbcXFxGBgYQH5+\nPgwGA3x8fHD48GF4eXlZ3ZYQQgghhBBC5MhhgyxCCCGEEEIIeRDNv9M9EUIIIYQQQoiEaJBFCCGE\nEEIIIXZEgyxCCCGEEEIIsSOHncLd1Xz11VdobGzE4cOHAZjPiHjgwAFwHIeoqCjZX2CSMYbo6Ggs\nXrwYgPm0wm+//ba0oWYxeXKUa9euQa1Wo7i4GCEhIVLHmpOXXnpJvB7cokWLUFpaKnGi2f3www94\n//33UV1djd7eXuzevRsKhQLLli3De++9Bzc3+R57uTf7lStXkJ2dLe7zGRkZiI+PlzbgNEwmEwoK\nCtDX1wej0YicnBw8+uijLlN7a/kDAwNdpv7zhRzbzLnsGxUVFTh//jyUSiUKCgoQFhbm1KyW7XVa\nWtp9/bwcavzpp5/is88+AwCMj4/j6tWrKC8vx8GDBxEYGAgA2LFjByIjIyXLaks/Yu31dnafc2/O\nq1evoqioCBzHQa1W4+DBg/Dz80NxcTHa29uh0WgAAJWVlTCZTMjNzcXY2BgCAgJQWloKT0/bLghv\nj6zT9W9yq+mbb76JgYEBAEBfXx/Cw8PxwQcfICcnB4ODg1CpVHB3d8eJEyecmnMufa5Da8oIKyoq\nYnFxcWzXrl3iusTERNbb28sEQWCvvvoq6+rqkjDh7K5fv86ys7OljjEn586dY/n5+Ywxxr7//nu2\nbds2iRPNzdjYGEtKSpI6xpxUVVWxhIQElpKSwhhjLDs7m3333XeMMcb27t3LvvzySynjzcgye11d\nHTt58qTEqWxTX1/PiouLGWOMDQ4OsnXr1rlU7a3ld6X6zxdybDNt3Tc6OzvZpk2bmCAIrK+vjyUn\nJzs1p7X22lo/L7caFxYWspqaGlZeXs4aGxunPCdVVlv6keleb2e2e5Y5N2zYwK5cucIYY+yTTz5h\nJSUljDHG0tPT2a1bt6b8bFFREWtoaGCMMXb8+HH20UcfOSyntaxzeQ9JWdNJd+7cYYmJiay/v58x\nxtjzzz/PBEGYso0zc9ra5zq6pvI8bOpkERERKCwsFJf1ej2MRiOCg4OhUCgQFRWFCxcuSBfQBl1d\nXejv78emTZvw2muv4bfffpM60qza2tqwdu1aAOZv3jo7OyVONDfd3d0YHR1FVlYWNm/ejMuXL0sd\naVbBwcE4cuSIuNzV1YWnn34aABAdHS3r/dwye2dnJ86fP48NGzagoKAAer1ewnQze+655/DGG28A\nMH/rzHGcS9XeWn5Xqv98Icc209Z9o62tDVFRUVAoFHj44YfB8zxu377ttJyW7fWlS5es9vNyqvGP\nP/6IX375BWlpaejq6kJDQwMyMzNRVlaGiYkJybLa0o9M93o7s92zzFleXo7ly5cDAHieh7u7OwRB\nQG9vL/bt24f09HTU19cDmPpec0b7bEv/JseaTjpy5Ag2btyIgIAADAwMYGhoCNu2bUNGRga++eYb\nAM79vGFrn+vomj5Qg6yzZ88iISFhyq2jowPx8fFQKBTidnq9XpxSAAAajQbDw8NSRLbK2u/h5+eH\nrVu3orq6GtnZ2cjLy5M65qws68xxHCYmJiRMNDceHh7YsmULTp48if379yM3N1f2+ePi4qBU/jNL\nmDEm7vty288tWWYPCwvDO++8gzNnziAoKAhHjx6VMN3MNBoNtFot9Ho9du7ciV27drlU7a3ld6X6\nzxdybDNt3Tek7lct2+s9e/ZMmf41mUdONT5+/Di2b98OAFizZg327t2LM2fOYGRkBDU1NZJltaUf\nme71dma7Z5kzICAAANDe3o7Tp0/jlVdewcjICDZu3IhDhw7hxIkT+Pjjj9Hd3Q29Xg+dTueUnNay\nzuU9JGVNAeDWrVtobW1FcnIyAPNUvaysLBw9ehQVFRUoLS3FrVu3nJrT1j7X0TV9oP4mKyUlBSkp\nKbNup9VqYTAYxGWDwQBvb29HRpsTa7/H6OgoOI4DAERGRuLvv/+espPIkWWdBUG4780rZ0uWLEFI\nSAgUCgWWLFmCBQsW4ObNm+KceVdw7xxjue3ns4mNjRXzxsbGoqioSOJEM/vzzz+xfft2ZGZm4sUX\nX8ShQ4fE51yh9pb5h4aGXKr+84Fc20xb9o2YmJj7+tXJD7HOYNle63Q63LlzZ0oeb29vjI2NyaLG\nQ0ND6OnpwerVqwEAL7/8sljTmJgYnDt3DjqdThZZrfUj1j5H6XQ6yfucL774AseOHUNVVRV8fX3B\n8zw2b94sDrhXr16N7u5uMb+Hh4ckOa31b9O9h6SuaWNjIxISEsTPoH5+fkhPT4dSqcTChQuxfPly\n9PT0OD2nLX2uo/fTB+qbLFtptVqoVCrcuHEDjDG0tLQgMjJS6lgzqqiowKlTpwCYp0UEBgbKeoAF\nmKdpNjc3AzCfaCQ0NFTiRHNTX1+PsrIyAEB/fz/0ej38/f0lTjU3jz32GC5evAgAaG5ulv1+fq8t\nW7ago6MDANDa2orHH39c4kTTGxgYQFZWFvLy8rB+/XoArlV7a/ldqf7zhRzbTFv3jYiICLS0tEAQ\nBPzxxx8QBAG+vr5Oy2nZXo+OjsLLy+u+fl4uNb506RKeeeYZAOZvihITE/HXX38BmFpTOWS11pZN\n93pL2e59/vnnOH36NKqrqxEUFAQAuH79OjIyMsDzPEwmE9rb28XaNjU1iTlXrlzptJzA3N5DUvcl\nra2tiI6OFpcvXLggTtUzGAz4+eefsXTpUqfmtLXPdXRNpT8EJlOT0794nkdUVBTCw8OljjSjrVu3\nIi8vD01NTeA4ziXOchcbG4tvv/0W6enpYIyhpKRE6khzsn79euzZswcZGRlQKBQoKSmRxVHlucjP\nz8fevXtRXl6OpUuXIi4uTupINissLERRURFUKhX8/Pxk/U3Khx9+iKGhIVRWVqKyshIA8O6776K4\nuNglam8t/+7du1FSUuIS9Z8v5Nhm2rpvaLVaREZGIi0tDYIgYN++fU7Naa29dnNzu6+ff+KJJ2RR\n456eHixatAgAoFAoUFxcjNdffx0eHh545JFHkJqaCo7jZJHVWj/CcZzV11uqPofneRw4cACBgYHY\nsWMHAOCpp57Czp07kZSUhNTUVKhUKiQlJWHZsmXIyclBfn4+6urq4OPjI5552lms9W/TvYek7sd7\nenrEQSsArFu3Di0tLUhNTYWbmxveeust+Pr6OjWnrX2uo/dTBWOM2e23IoQQQgghhJAHHE0XJIQQ\nQgghhBA7okEWIYQQQgghhNgRDbIIIYQQQgghxI5okEUIIYQQQgghdkSDLEIIIYQQQgixIxpkEUII\nIYQQQogd0SCLEEIIIYQQQuzo/wDXvtrR6JVm9gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f2e1432bf28>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.traceplot(trace_h, varnames=['mu']);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAg8AAAF8CAYAAABIe1hQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt4VPWdx/HPTIZAYsAAJuLiZeUi2uTpQ2HqJS5WbNAo\nIShKXQKsqMvFar0giIEIEYaLgFtNQHmgGLtBAiWb0JTVqFG6XgIuk2epJLvIrhFYkOL0MVxCAkmY\n2T9cU2IwzoE5OSeH9+t5fB6Y2/nMzzMnn5yZ+eIKhUIhAQAAhMltdQAAANC5UB4AAIAhlAcAAGAI\n5QEAABhCeQAAAIZQHgAAgCGUB8AGCgoKdMcdd2j06NGaPn26jhw50nLdjTfeqNGjR7f8V1paKkl6\n//33NXLkSN1zzz3atWtXy+2zs7NVUVHR7vY+//xz/epXv9KoUaOUkZGhCRMmyO/3n/fzWLt2rZ59\n9llJ0pw5c1RRUaEDBw7oJz/5yVlvn5eXp/nz5xvaRnFxsaZOnfqDt9u0aZPeeOMNQ499Ltp7foBT\neawOAFzotm/frjVr1uh3v/ud+vTpo82bN2vu3LnKzc1VTU2NLr74Yv3+979vc7+8vDy9/vrrOnjw\noNasWaPc3Fx9+umnOn78uFJSUr53ezU1NXrggQe0ePFiDRs2TJK0bds2TZs2TYWFhRo4cGBEntfC\nhQslffPD1QqVlZURey4AWqM8ABarrq5WSkqK+vTpI0m6/fbblZ2drcbGRv3Hf/yH3G63Jk6cqCNH\njuiOO+7QI488oqioKEVHR+vkyZNqaGhQly5dFAqFtGzZMi1ZsqTd7a1Zs0b33ntvS3GQpJtuukkv\nvviiunXrpgMHDmj8+PHq37+/Dh48qIKCAh04cEDLly9XQ0ODXC6XfvWrX2n48OFqamqSz+dTRUWF\nevfurd69e6t79+6SpIkTJ2r8+PFKTk5WMBjUnDlzVF1dLY/Ho+zsbA0ePLhVrsOHD2v+/Pk6dOiQ\nmpqaNHLkSE2bNq3d55KXl6eDBw8qEAjo4MGD6tWrl37961/r008/1fvvv6+PP/5Y3bp10/jx4/Xq\nq6/qnXfeUTAYVN++fTVv3jxdeumlmjhxoi6++GLV1NTo/vvv1yuvvKIPP/xQ0dHROn36tIYPH67X\nXntNdXV1WrZsmRobGxUIBJSSkqJFixa1yvP5559rzpw5amxsVCgU0n333afx48eHvS8AnQXlAbDY\nj3/8YxUUFOjgwYPq27eviouL1dTUpCNHjuj06dO6+eab9cwzz+jkyZOaMmWK4uLiNGnSJM2cOVNP\nPfWUunXrpueff15FRUW68cYb1bdv33a3V1VVpRkzZrS5/Gc/+5mkb84U/PnPf9aLL74or9ero0eP\nKisrS2vXrtXll1+uw4cP6xe/+IUGDRqkd999V3v37tW//uu/qrm5WRMmTGgpD2c6efKkbr75Zi1c\nuFAffvihnnzySb3zzjutbjNz5kxNmjRJt912m06dOqXJkyfryiuv1F133dXu8/H7/dq8ebPi4uI0\nbdo0bdy4UY8//rjee+89DRw4UOPHj9fmzZu1Z88ebdq0SR6PRxs3blR2drbWrFkjSerRo4fefPNN\nSdK7776r999/X2lpafroo4/Ut29fDRgwQNOnT9fjjz+uG264QSdOnNDPf/5zVVVVKT4+viXL2rVr\nddttt2nKlCkKBAJatGiRxo0bJ7ebd4jhLJQHwGI//elP9eijj+qxxx6Ty+XSvffeq/j4eHXp0kW/\n+MUvWm4XHR2tBx98UAUFBZo0aZK8Xq+KiookSUePHtWmTZtUUFCgVatWaefOnRowYMBZS4LL5VIw\nGGw3k8fjaTkzsHPnTgUCAT366KOtHuOzzz7Ttm3blJ6erujoaEVHR2vUqFH67LPP2jxejx49WkrA\nsGHDFAqFVFNT03J9fX29duzYoaNHj+rll19uuWz37t0/WB6uv/56xcXFSZJ+9KMf6ejRo21us3Xr\nVu3atUv33nuvJCkYDKqhoaHleq/X2/LnsWPHqqSkRGlpaSouLtbYsWMlSUuWLNEHH3ygVatWqaam\nRidPnlR9fX2r8jBixAjNmjVLn376qW666SZlZ2dTHOBIlAfAYnV1dbr++utbfkj95S9/UW5uruLj\n47V582Zde+21uvbaayVJoVBIHk/bl+3LL7+sqVOn6ssvv9S2bdv029/+VtnZ2dq2bZtuuummVrcd\nPHiwdu7cqeHDh7e6fMWKFbryyis1ZMgQRUdHt2zn9OnT6t+/vzZt2tRy28OHD6tXr17auHFjq8eI\nioo663P87g/QUCikLl26tPw9GAwqFAppw4YNiomJkSR9/fXX6tq16/cv3P/r1q1by59dLpfO9s/1\nBINB/eM//qMyMzMlSY2Nja1KRmxsbMuf09LStHjxYn3++efasWNHy9tA48eP17XXXqthw4bpzjvv\n1J/+9Kc22xo+fLjefvttVVRUaNu2bVq5cqU2bNigK6+88gefB9CZUIkBi3311VeaOHGi6urqJEmv\nvPKKRo4cKZfLpf/+7/9Wbm6uTp8+rZMnT+qNN95o85v47t279eWXX+rnP/+5GhsbW37ou93uVr9d\nf+vhhx/Wpk2b9NFHH7Vc9sEHH6igoKClpJxp8ODB2rdvn3bs2CFJ+q//+i/dcccd+uqrrzRs2DBt\n3rxZp06d0qlTp1pO/X/XkSNHtHXrVknffEuka9euuuqqq1quj4uL0+DBg5Wfny9JOnbsmMaNG6f3\n3nsv7HX8rqioKDU3N0uS/u7v/k5FRUUta/zyyy/rmWeeOev9unbtqpEjR+rZZ5/V7bffrpiYGB09\nerTl7Z7bb79dhw8f1v79+9ucwXn66af15ptvauTIkZo3b57i4uJ06NChc34OgF1x5gGwWL9+/TRl\nyhSNHTtWwWBQQ4cO1dy5cyVJjz32mObPn69Ro0apublZaWlpLWcovrVkyRLNmzdPkjRo0CD17t1b\no0aN0hVXXNHqQ5Hfuuqqq7Rq1Sq99NJLeuGFFxQMBtWrVy+9+uqruuaaa9p8O6JXr17Kzc3V0qVL\nderUKYVCIS1dulR9+/bV3//932v//v1KT09XfHx8q0Jwpt69e+udd97RSy+9pJiYGOXl5bU5g7J8\n+XItWLBAo0aNUmNjo9LT05WRkXHO63rLLbdowYIFkqTJkye3fFbD5XLpsssua/eDpWPHjtW6deuU\nk5MjSbr44os1ZcoU3XPPPYqPj1fPnj01ZMgQ7du3T1dccUXL/X75y19qzpw52rhxo6KiopSamqrr\nr7/+nJ8DYFcu/kluAABgBG9bAAAAQygPAADAEMoDAAAwhPIAAAAMccy3LQKB4x26vZ49Y1VbW9+h\n2+yMhg5Nltvt0o4du374xhc49qnwsVbhY63Cwzq1lZDQdlrstzjzcI48nrMPw0FrlZVV2rt3r9Ux\nOgX2qfCxVuFjrcLDOhlDeQAAAIZQHmAqny9HWVlZVscAAEQQ5QGmKikpUmFhodUxAAARRHkAAACG\nUB4AAIAhlAcAAGAI5QEAABhCeYCpmPMAAM5DeQAAAIZQHmAqn485DwDgNJQHmIo5DwDgPJQHAABg\nCOUBAAAYQnkAAACGUB4AAIAhlAeYijkPAOA8lAcAAGBIu+Xh1KlT2rRpkyQpLy8v7K/cHTlyRH/4\nwx/Cuu0nn3yim266SUuXLm257N1339XTTz8tSfL7/Ro9erSeeuqpsB4P9uLzMecBbfn9buXmRsvv\n5/cXoDPytHdlIBDQpk2bNHbsWEMP+tlnn+n999/XqFGjwrr9jTfeqGeeeUaS5PP59NFHH+m6666T\nJHm9Xs2ePVsbNmwwlAH2UFJSJLfbpenTZ1sdpVPKzIxReXm7L9NOrut53Ld7xFJ0pNTUZq1f32B1\nDOC8tHtUWrVqlf7nf/5HK1askCS99957Kisr05EjR/TEE0/otttu01tvvaXXX39dbrdbQ4cO1YwZ\nM7Rq1Srt3r1bGzdu1E9+8hMtWbJEp0+fVm1trXJycjRkyJDv3eaQIUOUmpqqjRs3GnoiPXvGyuOJ\nMnSf85WQ0DkOXsnJUnW1VVvfJ0lKTLRq+51N59incO7Kyz1KTOzo/8/sV+HpvOuUlCRVVXXc9tot\nD9OmTdOePXv02GOPKS8vT5deeqkWLlyoTz75RL/5zW80ZMgQ5eXl6V/+5V8UExOjmTNn6uOPP9a0\nadO0YcMG3X///XrzzTc1a9YsDRo0SH/4wx9UXFzcbnm466679Mknnxh+IrW19Ybvcz4SErorEDje\nods8V1u3WrftoUOT5Xa7tGPHLutCdBKdaZ86H36/WxkZsWpudsnjCam0tF5eb9DQY1woaxUJrFV4\nnLBOgUBkH6+9X5ANnQ9NSkqSJF1yySU6efKk9u/fr6+//lpTpkyRJJ04cUL79+9Xv379Wu6TmJio\nV155Rd26ddOJEycUFxd3Ls8BgEN4vUGVltarosKjlJRmw8UBgPXaLQ9ut1vB4F9f2C6Xq9X1l19+\nuS677DK99tpr6tKli4qLi3Xdddeprq6u5X4LFy7U8uXL1b9/f+Xm5urgwYMmPA0AnYnXG5TX22h1\nDADnqN3y0Lt3bzU1NWnZsmXq1q1bm+t79eqlSZMmaeLEiTp9+rT69u2rO++8U8eOHdOePXv0+uuv\nKyMjQ0888YR69OihPn36qLa21rQnA/uprKxyxOlAAMBftVseunbtqt///vdtLu/fv78KCgokSaNH\nj9bo0aNbXR8TE6O33nqr5e8PPvigoVA33HCDbrjhBkP3AQAAHcMWX7Levn17qzkPZ/L7/Vq0aFEH\nJ0Kk+HzMeQAAp3GFQqGQ1SEioaNPi3MqPjx82yJ87FPhY63Cx1qFh3Vqq71vW9jizAMAAOg8KA8A\nAMAQygMAADCE8gAAAAyhPMBUlZVV2rt3r9UxAAARRHkAAACGUB5gKp+POQ8A4DSUB5iqpKRIhYWF\nVscAAEQQ5QEAABhCeQAAAIZQHgAAgCGUBwAAYAjlAaZizgMAOA/lAQAAGEJ5gKl8PuY8AIDTUB5g\nKuY8AIDzUB4AAIAhlAcAAGAI5QEAABhCeQAAAIZQHmAq5jwAgPNQHgAAgCGUB5jK52POAwA4DeUB\npmLOAwA4D+UBAAAYQnkAAACGUB4AAIAhlAcAAGAI5QGmYs4DADgP5QEAABhCeYCpfD7mPACA01Ae\nYCrmPACA83isDiBJxcXFys3N1X333Se/36/Tp08rFApp/vz52rNnj1566SWlpqZqxowZVkcFLlh+\nv1sVFR6lpDTL6w1aHQeAhWxRHiQpPT1d//u//6sJEyYoNTVVH374of7pn/5JK1asUH19vWpqaqyO\nCIQtMzNG5eXn8vLqHvEskdfV6gD/7/zWKjW1WevXN0QoC3BhsU15kKRZs2ape/dvDginT59W167h\nH6R69oyVxxNlVrSzSkjomAN9crJUXd0hmzLBPklSYqLFMYDvKC/3KDGxM5S1SHDW80xKkqqqIv+4\nHXVMdwJblYdevXpJkmpqavTCCy9o5cqVYd+3trberFhnlZDQXYHA8Q7Z1tatHbIZUwwdmiy326Ud\nO3ZZHcX2OnKfMsrvdysjI1bNzS55PCGVltZb+taFndfKbpy6VoFAZB/Pqet0PtorU7YqD5K0fft2\nPf/881q6dKn69etndRycp8rKKl6UDuD1BlVaWs9nHgBIsll52L59uxYuXKjf/OY36tu3r9VxAJzB\n6w3K6220OgYAG7BVeVi0aJGampr07LPPSpKuvvpqzZ8/3+JUOB8+X45iY6M1ffpsq6MAACLEVuWh\ntLTU6giIsJKSIrndLsoDADiIbYZEbdmyRfn5+W0uLysr0+rVqy1IBAAAzsYWZx7GjBmjMWPGnPW6\ntLQ0paWldXAiAADwfWxz5gEAAHQOlAcAAGAI5QGmqqys0t69e62OAQCIIMoDAAAwhPIAU/l8OcrK\nyrI6BgAggigPMFVJSZEKCwutjgEAiCDKAwAAMITyAAAADKE8AAAAQygPAADAEMoDTMWcBwBwHsoD\nAAAwhPIAU/l8zHkAAKehPMBUzHkAAOehPAAAAEMoDwAAwBDKAwAAMITyAAAADKE8wFTMeQAA56E8\nAAAAQygPMJXPx5wHAHAaygNMxZwHAHAeygMAADCE8gAAAAyhPAAAAEMoDwAAwBDKA0zFnAcAcB7K\nAwAAMITyAFP5fMx5AACnoTzAVMx5AADnoTwAAABDKA8AAMAQW5SH4uJi3XrrrXr11Vf1wAMPKDMz\nU4888ojq6upUVlamtLQ0LV++3OqYgOP5/W7l5kbL77fFoQGATXmsDvCt9PR0ff3117rnnnt09913\nKy8vT0VFRZo0aZLq6+tVU1NjdUSgRWZmjMrLzXj5dDfhMc9FV6sDhCH8tUpNbdb69Q0mZgEuLLYp\nD5I0e/ZshUIhBYNBHTp0SH/zN38T9n179oyVxxMV8UzJyVJ19fdda5cDvZ3tkyQlJlocAxe08nKP\nEhMv5Ners557UpJUVRX5x01IcNY6mclW5cHlcqm5uVmjR4/WqVOn9Oijj4Z939raelMybd169ssT\nErorEDhuyjadhrUKj9Xr5Pe7lZERq+ZmlzyekEpL6+X1Bi3L0x6r16ozcepaBQKRfTynrtP5aK9M\n2ao8SFKXLl305ptvqqKiQrNmzdK6deusjoTz4PPlKDY2WtOnz7Y6Cn6A1xtUaWm9Kio8Sklptm1x\nAGA9W30qKicnR9u3b5ckXXTRRXK5XBYnwvlizkPn4vUG9fjjjRQHAO2y1ZmHiRMnKicnRytXrpTb\n7VZOTo7VkQAAwHfYqjz0799fBQUFVscAAADtsM3bFlu2bFF+fn6by8vKyrR69WoLEgEAgLOxxZmH\nMWPGaMyYMWe9Li0tTWlpaR2cCAAAfB/bnHmAM1VWVmnv3r1WxwAARBDlAQAAGEJ5gKl8vhxlZWVZ\nHQMAEEGUB5iKOQ8A4DyUBwAAYAjlAQAAGEJ5AAAAhlAeAACAIZQHmIo5DwDgPJQHAABgCOUBpvL5\nmPMAAE5DeYCpmPMAAM5DeQAAAIZQHgAAgCGUBwAAYAjlAQAAGEJ5gKmY8wAAzkN5AAAAhlAeYCqf\njzkPAOA0lAeYijkPAOA8lAcAAGAI5QEAABhCeQAAAIZQHgAAgCGUB5iKOQ8A4DyUBwAAYAjlAaby\n+ZjzAABOQ3mAqZjzAADOQ3kAAACGUB4AAIAhlAcAAGAI5QEAABhii/JQXFysW2+9Vfn5+ZKkf//3\nf9fPfvYzSVJZWZnS0tK0fPlyKyPiHDHnwVn8frdyc6Pl99vi0AHAIh6rA3wrPT1dDz74oA4dOqT8\n/Hw1NzdLktLS0lRfX6+amhqLEwLhy8yMUXn5uby8ukc8izm6Wh1AkVir1NRmrV/fEIEswIXFNuVB\nkk6dOqV58+ZpwYIFGjNmjKH79uwZK48nyqRkZ5eQYO2BPjlZqq62NIIBneWHIi4k5eUeJSZeCPvm\nhfAczy4pSaqqCu+2Vh/TOxNblYf58+froYce0qWXXmr4vrW19SYk+n4JCd0VCBzv0G1+19atlm4+\nLEOHJsvtdmnHjl1WR7E9O+xT7fH73crIiFVzs0seT0ilpfXyeoOWZLH7WtkJayUFAj98G9aprfbK\nlG3Kw9GjR+X3+7V//36tXLlSR48e1VNPPaVf//rXVkcDIMnrDaq0tF4VFR6lpDRbVhwAWM825eHi\niy/W22+/3fL3m2++meIA2IzXG5TX22h1DAAW4yPTAADAENuWh48//tjqCAAA4CxsUx62bNnSMufh\nTGVlZVq9erUFiRAJzHkAAOdxhUKhkNUhIqGjPyXLJ3PDx1qFh3UKH2sVPtYqPKxTW+1928I2Zx7g\nTD5fjrKysqyOAQCIIMoDTFVSUqTCwkKrYwAAIojyAAAADKE8AAAAQygPAADAEMoDAAAwhPIAUzHn\nAQCch/IAAAAMoTzAVD4fcx4AwGkoDzAVcx4AwHkoDwAAwBDKAwAAMITyAAAADKE8AAAAQygPMBVz\nHgDAeSgPAADAEMoDTOXzMecBAJyG8gBTMecBAJyH8gAAAAyhPAAAAEMoDwAAwBDKAwAAMITyAFMx\n5wEAnIfyAAAADKE8wFQ+H3MeAMBpKA8wFXMeAMB5KA8AAMAQygMAADCE8gAAAAyhPAAAAEMoDzAV\ncx4AwHkoDwAAwBCP1QEkqbi4WLm5uXrggQe0atUqXXPNNZKk1NRURUVF6bXXXtPkyZM1btw4i5PC\nKJ8vR7Gx0Zo+fbbVUQAAEWKL8iBJ6enpGjRokNLT0/Xcc8+1uq62ttaiVDhfJSVFcrtdlAeL+P1u\nVVR4lJLSLK83aHUcAA5hm/IgSVVVVaqurtaECRPUq1cvZWdnKzEx0epYuEBkZsaovNzql0R3kx63\nq0mPa53U1BitX99gdQzggmT1kbKVfv36KTk5WSkpKSotLZXP51Nubm5Y9+3ZM1YeT5TJCVtLSPjm\nQJ+cLFVXd+imO5F9kiQ6ICKtvNyjxESzypbTdMw6JSVJVVUdsilTfHtMxw+zVXm48cYbFRMTI0ka\nMWJE2MVBkmpr682KdVYJCd0VCByXJG3d2qGb7lSGDk2W2+3Sjh27rI5ie2fuU5Hg97uVkRGr5maX\nPJ6QSkvrHfPWRaTXysk6eq0CgQ7bVESxT7XVXpmy1bctsrOz9fbbb0uStm3bpqSkJIsTAZ2X1xtU\naWm9srNPOao4ALCerc48PP3005o9e7YKCwsVExMjn89ndSScp8rKKhq9hbzeoLzeRqtjAHAYW5WH\nK664QgUFBVbHAAAA7bDN2xZbtmxRfn5+m8vXrVunkpISCxIhEny+HGVlZVkdAwAQQa5QKBSyOkQk\ndPRpcU7Fh4cPTIaPfSp8rFX4WKvwsE5tdZoPTAIAAPujPAAAAEMoDwAAwBDKAwAAMITyAFNVVlZp\n7969VscAAEQQ5QEAABhCeYCpfD7mPACA01AeYKqSkiIVFhZaHQMAEEGUBwAAYAjlAQAAGEJ5AAAA\nhlAeAACAIZQHmIo5DwDgPJQHAABgCOUBpvL5mPMAAE5DeYCpmPMAAM5DeQAAAIZQHgAAgCGUBwAA\nYAjlAQAAGEJ5gKmY8wAAzkN5AAAAhlAeYCqfjzkPAOA0lAeYijkPAOA8lAcAAGAI5QEAABhCeQAA\nAIZQHgAAgCGUB5iKOQ8A4DyUBwAAYAjlAaby+ZjzAABOQ3mAqZjzAADO47E6gCQVFxcrNzdX999/\nv7744gsdOHBATU1Neu655/Tll1/qpZdeUmpqqmbMmGF1VAAALni2KA+SlJ6erubmZg0cOFBLly7V\n7t27tXv3bt19992qr69XTU2N1REBx/H73aqo8CglpVleb9DqOAA6CduUB0n66KOPdOedd+rhhx/W\nRRddpHnz5lkdCWglMzNG5eVmvmy6m/jY7elq0XbPx/evVWpqs9avb+jALMCFxVbloba2VseOHdPa\ntWu1efNmvfDCC1q6dGlY9+3ZM1YeT5TJCVtLSDB2oE9OlqqrTQpjW/skSYmJFsfABaW83KPERKuK\nmB11rrVISpKqqjp+u0aP6RcyW5WH+Ph43XbbbZKk4cOHa/Xq1WHft7a23qxYZ5WQ0F2BwHFD99m6\n1aQwNncua3Uh6uh18vvdysiIVXOzSx5PSKWl9Z3mrQv2qfB11rUKBDp2e511nczUXpmyVXkYOnSo\n/u3f/k3JycnasWOHBgwYYHUkwLG83qBKS+v5zAMAw2xVHqZOnars7Gzdf//98ng8euGFF6yOhPPk\n8+UoNjZa06fPtjoKzsLrDcrrbbQ6BoBOxlblIT4+XitWrLA6BiKopKRIbreL8gAADmKbIVFbtmxR\nfn5+m8vLysoMffYBAACYyxZnHsaMGaMxY8ac9bq0tDSlpaV1cCIAAPB9bHPmAQAAdA6UBwAAYAjl\nAaaqrKzS3r17rY4BAIggygMAADCE8gBT+Xw5ysrKsjoGACCCKA8wVUlJkQoLC62OAQCIIMoDAAAw\nhPIAAAAMoTwAAABDKA8AAMAQygNMxZwHAHAeygMAADCE8gBT+XzMeQAAp6E8wFTMeQAA56E8AAAA\nQygPAADAEMoDAAAwhPIAAAAMoTzAVMx5AADnoTwAAABDKA8wlc/HnAcAcBrKA0zFnAcAcB7KAwAA\nMITyAAAADKE8AAAAQygPAADAEMoDTMWcBwBwHsoDAAAwhPIAU/l8zHkAAKehPMBUzHkAAOehPAAA\nAEMoDwAAwBCP1QEkqbi4WLm5uTp27JiSkpIkSYFAQD169FBGRoZee+01TZ48WePGjbM4KQAAsEV5\nkKT09HTNmDFDktTU1KTMzEwtWLBAgwYNUm1trcXpAESC3+9WRYVHKSnN8nqDVscBcI5sUx7OtG7d\nOt18880aNGiQ1VFwniorq5SQ0F2BwHGro3RqmZkxKi+35cv1HHWNwGN0j8BjREZqarPWr2+wOgbQ\nYWx3NGpsbNSGDRtUVFRk6H49e8bK44kyKdXZJSTY5+B1PpKTpepqs7fijLUyH+vUGZWXe5SYaOf/\nd3bOZid/XaekJKmqysIoNme78rBt2zb99Kc/Vffuxnb22tp6kxKdnZN+m9661bzH9vlyFBsbrenT\nZ5u3EYdw0j51Nn6/WxkZsWpudsnjCam0tP6c37pw+lpFEmsVnrOtUyBgURibaO8XZNuVh4qKCt1y\nyy1Wx0CElJQUye12UR4grzeo0tJ6PvMAOIDtysMXX3yhu+++2+oYAEzg9Qbl9TZaHQPAebJdeVi9\nerXVEQAAQDtsMyRqy5Ytys/Pb3P5unXrVFJSYkEiAABwNq5QKBSyOkQkdPQHgvgQUniGDk2W2+3S\njh27rI5ie+xT4WOtwsdahYd1aqu9D0za5swDnKmyskp79+61OgYAIIIoDwAAwBDKA0zl8+UoKyvL\n6hgAgAiiPMBUJSVFKiwstDoGACCCKA8AAMAQygMAADCE8gAAAAyhPAAAAEMoDzAVcx4AwHkoDwAA\nwBDKA0zl8zHnAQCchvIAUzHnAQCch/IAAAAMoTwAAABDKA8AAMAQygMAADCE8gBTMecBAJyH8gAA\nAAyhPMBJv3CUAAAHN0lEQVRUPh9zHgDAaSgPMBVzHgDAeSgPAADAEMoDAAAwhPIAAAAMoTwAAABD\nKA8wFXMeAMB5KA8AAMAQygNM5fMx5wEAnIbyAFMx5wEAnIfyAAAADKE8AAAAQygPAADAEMoDAAAw\nxBblobi4WLfeeqtWrFihCRMmaPz48frlL3+phoYGlZWVKS0tTcuXL7c6Js4Bcx4AwHlsUR4kKT09\nXceOHdOdd96pN954QwMHDlRRUZHS0tI0ZcoUq+MBOIPf71ZubrT8ftscQgB0II/VAc503XXX6c9/\n/rMkqa6uTn369LE4Ec6Xz5ej2NhoTZ8+2+oolsnMjFF5ebgvte6mZom8rhZuO/JrlZrarPXrGyL+\nuIDT2Ko89OnTRy+++KK2bNmixsZGPfbYY2Hft2fPWHk8USamayshwfwDfXKyVF1t+mZM9KIkackS\ni2MAYSgv9ygxsbMVuHA48TmZ4fvXKSlJqqrqwCg2Z6vysHTpUi1evFjDhg3TH//4R82aNUurV68O\n6761tfUmp2stIaG7AoHjpm9n61bTN2GqoUOT5Xa7tGPHLquj2F5H7VPny+93KyMjVs3NLnk8IZWW\n1svrDXZohs6yVnbAWoUnnHUKBDoojE209wuyrcpDjx491L37N2ETExN17NgxixMB+C6vN6jS0npV\nVHiUktLc4cUBgPVsVR6ee+45zZ8/X8FgUKFQSHPnzrU6EoCz8HqD8nobrY4BwCK2Kg8DBgzQP//z\nP1sdAwAAtMM237PasmWL8vPz21xeVlYW9uceYD/MeQAA53GFQqGQ1SEioaM/EMSHkMLHWoWHdQof\naxU+1io8rFNb7X1g0jZnHuBMPl+OsrKyrI4BAIggygNMVVJSpMLCQqtjAAAiiPIAAAAMoTwAAABD\nKA8AAMAQygMAADCE8gBTMecBAJyH8gAAAAyhPMBUPh9zHgDAaSgPMBVzHgDAeSgPAADAEMoDAAAw\nhPIAAAAMoTwAAABDHPNPcgMAgI7BmQcAAGAI5QEAABhCeQAAAIZQHgAAgCGUBwAAYAjlAQAAGEJ5\nAAAAhnisDtBZhUIh3XLLLfrbv/1bSdLgwYP19NNPWxvKRoLBoHJycvTZZ58pOjpaPp9PV111ldWx\nbOuee+5RXFycJOnyyy/X4sWLLU5kP3/605+0fPlyFRQUaN++fXr22Wflcrk0cOBAzZs3T243vwtJ\nrdfpP//zPzV16tSW49S4ceN01113WRvQBpqamjR79mwdPHhQjY2NeuSRRzRgwAD2KQMoD+do//79\nSkpK0qpVq6yOYkvl5eVqbGzUxo0btXPnTi1ZskSvvvqq1bFs6dSpUwqFQiooKLA6im2tWbNGpaWl\niomJkSQtXrxYTz75pG644QbNnTtX7733nkaMGGFxSut9d52qq6v14IMP6qGHHrI4mb2UlpYqPj5e\ny5Yt05EjR3T33Xfr2muvZZ8ygFp1jqqrq3X48GFNnDhRkydPVk1NjdWRbKWyslLDhg2T9M1Zmaqq\nKosT2dfu3bvV0NCghx56SP/wD/+gnTt3Wh3Jdq688krl5eW1/L26ulrXX3+9JOmWW25RRUWFVdFs\n5bvrVFVVpT/+8Y8aP368Zs+erbq6OgvT2UdaWpqeeOIJSd+cRY6KimKfMojyEIZNmzYpPT291X+X\nXHKJpkyZooKCAk2dOlUzZ860Oqat1NXVtZyGl6SoqCg1NzdbmMi+unXrpocfflhr167V888/rxkz\nZrBW33HHHXfI4/nridJQKCSXyyVJuuiii3T8+HGrotnKd9fpxz/+sZ555hm98cYbuuKKK7Ry5UoL\n09nHRRddpLi4ONXV1enxxx/Xk08+yT5lEG9bhGHs2LEaO3Zsq8saGhoUFRUlSfJ6vfrqq69a7XwX\nuri4OJ04caLl78FgsNVBDX919dVX66qrrpLL5dLVV1+t+Ph4BQIBXXbZZVZHs60z34s+ceKEevTo\nYWEa+xoxYkTL2owYMUILFiywOJF9HDp0SI8++qgyMzM1atQoLVu2rOU69qkfxpmHc7RixQr99re/\nlfTNaefLLruM4nCGIUOG6IMPPpAk7dy5U9dcc43FieyrqKhIS5YskSQdPnxYdXV1SkhIsDiVvf3o\nRz/SJ598Ikn64IMP5PV6LU5kTw8//LA+/fRTSdK2bduUlJRkcSJ7+Mtf/qKHHnpIM2fO1H333SeJ\nfcoo/lXNc3T06FHNnDlT9fX1ioqK0ty5c9W/f3+rY9nGt9+22LNnj0KhkBYtWsT6fI/GxkZlZWXp\nyy+/lMvl0owZMzRkyBCrY9nOgQMHNH36dP3ud7/TF198oeeee05NTU3q16+ffD5fy5nAC92Z61Rd\nXa0FCxaoS5cuuuSSS7RgwYJWbydeqHw+n9566y3169ev5bI5c+bI5/OxT4WJ8gAAAAzhbQsAAGAI\n5QEAABhCeQAAAIZQHgAAgCGUBwAAYAjlAQAAGEJ5AAAAhvwf71cEF2i41fwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f2e26c14dd8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.forestplot(trace_h, varnames=['theta']);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Deviance Information Criterion (DIC)\n",
    "\n",
    "DIC (Spiegelhalter et al. 2002) is an information theoretic criterion for estimating predictive accuracy that is analogous to Akaike's Information Criterion (AIC). It is a more Bayesian approach that allows for the modeling of random effects, replacing the maximum likelihood estimate with the posterior mean and using the effective number of parameters to correct for bias."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "90.832077760613302"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pooled_dic = pm.dic(trace_p, pooled)\n",
    "\n",
    "pooled_dic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "124.82406598767808"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hierarchical_dic = pm.dic(trace_h, hierarchical)\n",
    "    \n",
    "hierarchical_dic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Widely-applicable Information Criterion (WAIC)\n",
    "\n",
    "WAIC (Watanabe 2010) is a fully Bayesian criterion for estimating out-of-sample expectation, using the computed log pointwise posterior predictive density (LPPD) and correcting for the effective number of parameters to adjust for overfitting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "61.121151094964134"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pooled_waic = pm.waic(trace_p, pooled)\n",
    "    \n",
    "pooled_waic.WAIC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "61.570270537072844"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hierarchical_waic = pm.waic(trace_h, hierarchical)\n",
    "    \n",
    "hierarchical_waic.WAIC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PyMC3 includes two convenience functions to help compare WAIC for different models. The first of this functions is `compare`, this one computes WAIC (or LOO) from a set of traces and models and returns a DataFrame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>WAIC</th>\n",
       "      <th>pWAIC</th>\n",
       "      <th>dWAIC</th>\n",
       "      <th>weight</th>\n",
       "      <th>SE</th>\n",
       "      <th>dSE</th>\n",
       "      <th>warning</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>61.1212</td>\n",
       "      <td>0.674471</td>\n",
       "      <td>0</td>\n",
       "      <td>0.555905</td>\n",
       "      <td>2.20084</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>61.5703</td>\n",
       "      <td>1.08163</td>\n",
       "      <td>0.449119</td>\n",
       "      <td>0.444095</td>\n",
       "      <td>1.96714</td>\n",
       "      <td>0.0259316</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      WAIC     pWAIC     dWAIC    weight       SE        dSE warning\n",
       "1  61.1212  0.674471         0  0.555905  2.20084          0       0\n",
       "0  61.5703   1.08163  0.449119  0.444095  1.96714  0.0259316       0"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_comp_WAIC = pm.compare((trace_h, trace_p), (hierarchical, pooled))\n",
    "df_comp_WAIC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have many columns so let check one by one the meaning of them:\n",
    "\n",
    "\n",
    "1. The first column clearly contains the values of WAIC. The DataFrame is always sorted from lowest to highest WAIC. The index reflects the order in which the models are passed to this function.\n",
    "\n",
    "2. The second column is the estimated effective number of parameters. In general, models with more parameters will be more flexible to fit data and at the same time could also lead to overfitting. Thus we can interpret pWAIC as a penalization term, intuitively we can also interpret it as measure of how flexible each model is in fitting the data. \n",
    "\n",
    "3. The third column is the relative difference between the value of WAIC for the top-ranked model and the value of WAIC for each model. For this reason we will always get a value of 0 for the first model.\n",
    "\n",
    "4. Sometimes when comparing models, we do not want to select the \"best\" model, instead we want to perform predictions by averaging along all the models (or at least several models). Ideally we would like to perform a weighted average, giving more weight to the model that seems to explain/predict the data better. There are many approaches to perform this task, one of them is to use Akaike weights based on the values of WAIC for each model. These weights can be loosely interpreted as the probability of each model (among the compared models) given the data. One caveat of this approach is that the weights are based on point estimates of WAIC (i.e. the uncertainty is ignored).\n",
    "\n",
    "5. The fifth column records the standard error for the WAIC computations. The standard error can be useful to assess the uncertainty of the WAIC estimates. Nevertheless, caution need to be taken because the estimation of the standard error assumes normality and hence could be problematic when the sample size is low.\n",
    "\n",
    "6. In the same way that we can compute the standard error for each value of WAIC, we can compute the standard error of the differences between two values of WAIC. Notice that both quantities are not necessarily the same, the reason is that the uncertainty about WAIC is correlated between models. This quantity is always 0 for the top-ranked model.\n",
    "\n",
    "7. Finally we have the last column named \"warning\". A value of 1 indicates that the computation of WAIC may not be reliable, this warning is based on an empirical determined cutoff value and need to be interpreted with caution. For more details you can read this [paper](https://arxiv.org/abs/1507.04544)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second convenience function takes the output of `compare` and produces a summary plot in the style of the one used in the book [Statistical Rethinking](http://xcelab.net/rm/statistical-rethinking/) by Richard McElreath (check also [this port](https://github.com/aloctavodia/Statistical-Rethinking-with-Python-and-PyMC3) of the examples in the book to PyMC3)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdkAAAFXCAYAAADu/TSqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFqlJREFUeJzt3X90lnX9x/HXvY2bDaZfVo7CgAHlQA9HFGJGuVYZkS3g\nDBScxY7K18NRQwwUtQ7mKQ6OlRrxB2R0MhAXc56ZBAsOcTwjBISxgfwYRsE8y8TFOMVuRtvuXd8/\n+DoD9gvYe/fnM5+Pv7jvm933mw83vPa5dr3uKxQEQSAAANDt4mI9AAAAvRUhCwCAEUIWAAAjhCwA\nAEYIWQAAjBCyAAAYSejuJ6ytPd3dT9mulJR+OnXqTI+93scN62urq+u7Zs0qSdKsWf9rPVKvwvvX\nFuv7kdTUq9p9zOudbEJCfKxH6NVYX1usry3W1xbr2zVehywAAC7r9sPFALoXh4kBf7GTBQDACCEL\nOG7Hjm3asWNbrMcAcBkIWcBxR48e0dGjR2I9BoDLQMgCAGCEkAUAwAhnFwMOa2ho0N695Tp58qT6\n9bta2dlTlJiYGOuxAHRRl3ay+/bt06xZs6xnwWUqKSlWVtYEDRqUoqysCSopKY71SOgGFRXlysgY\noz179qi5uVlFRYUaP/5GVVSUx3o0AF3U6U72V7/6lV5//XUlJSX1xDy4RCUlxZoz577W24cPH2y9\nnZNzR6zGwhVqaGhQXl6uCgqe1+23Z7feX1q6QXl5udq9ez87WsADoSAIgo5+w6ZNmzRy5EgtXLhQ\nRUVFnT5hd3928bhxo9t9LC4upJaWDsfv9d5//x9qamq66P4+ffro058edEXPzfra6mh9I5GIzpyJ\nKDV14EWP1dZ+oH79+qt///7WI3qN96+tnl7f8vIDPfZal6qjzy7udCc7adIk1dTUdPnFUlL6detn\nWsbFha7o8d6urYD98P7uWJuP+/paa299o9FmhcNhxcWFNG7cOElSefm5w8ThcFjRaDN/N13AGtnq\nyfXtKMhc1u0nPnX3VRl273673cdSU6/q0av+uCgra4IOHz540f033DBab7zx5hU9N+trq6P1ffXV\nIhUVFWrdupLWq/CsXPmiJGnmzBzNmJGr6dNn9NSoXuL9a6un19flv8teexUeSI88sqDN++fNm9/D\nk6A7ZWdP0aFDB1VauuG8+0tLN+jQoYPKzp4So8kAXAoqPJ778OSmZcue0zvvVCk9fZTmzZvPSU+e\nS0xM1OrVhcrLy9WAAf+jwYOH6A9/OBewq1cXctIT4IlOT3y6VD25pedwkC3W11ZX1rehoUE/+MF8\nnTxZp6lTp9OTvQS8f22xvh+5ohOfAMROUlKSxo49d+ITP4MF/EPIAo7jerKAvzjxCQAAI4Qs4Diu\nJwv4i5AFHMf1ZAF/EbIAABghZAEAMELIAgBghJAFAMAIPVnAcfRkAX+xkwUAwAghCziOnizgL0IW\ncBw9WcBfhCwAAEYIWQAAjBCyAAAYIWQBADBCTxZwHD1ZwF/sZAEAMELIAo6jJwv4i5AFHEdPFvAX\nIQsAgBFCFgAAI4QsAABGCFkAAIzQkwUcR08W8Bc7WQAAjBCygOPoyQL+ImQBx9GTBfxFyAIAYISQ\nBQDACCELAIARQhYAACP0ZAHH0ZMF/MVOFgAAI4Qs4Dh6soC/CFnAcfRkAX8RsgAAGCFkAQAwQsgC\nAGCEkAUAwAg9WcBx9GQBf7GTBQDACCELOI6eLOAvQhZwHD1ZwF+ELAAARghZAACMELIAABghZAEA\nMEJPFnAcPVnAX+xkAQAwQsgCjqMnC/iLkAUcR08W8BchCwCAEUIWAAAjhCwAAEYIWQAAjNCTBRxH\nTxbwFztZAACMELKA4+jJAv4iZAHH0ZMF/EXIAgBghJAFAMAIIQsAgBFCFgAAI/RkAcfRkwX8xU4W\nAAAjhCzgOHqygL8IWcBx9GQBfxGyAAAYIWQBADBCyAIAYISQBQDACD1ZwHH0ZAF/sZMFAMAIIQs4\njp4s4C9CFnAcPVnAX4QsAABGCFkAAIwQsgAAGCFkAQAwQk8WcBw9WcBf7GQBADBCyAKOoycL+IuQ\nBRxHTxbwFyELAIARQhYAACOELAAARghZAACM0JMFHEdPFvAXO1kAAIwQsoDj6MkC/iJkAcfRkwX8\nRcgCAGCEkAUAwAghCwCAEUIWAAAj9GQBx9GTBfzFThYAACOELOA4erKAvwhZwHH0ZAF/EbIAABgh\nZAEAMELIAgBghJAFAMAIPVnAcfRkAX+xkwUAwAghCziOnizgL0IWcBw9WcBfhCwAAEYIWQAAjBCy\nAAAYIWQBADBCTxZwHD1ZwF/sZAEAMELIAo6jJwv4i5AFHEdPFvAXIQsAgBFCFgAAI4QsAABGCFkA\nAIzQkwUcR08W8Bc7WQAAjBCygOPoyQL+ImQBx9GTBfxFyAIAYISQBQDACCELAICRTis8LS0tevrp\np3XkyBGFw2EtXrxYaWlpPTEbAMedOPEPnT17Vmlpw2M9CuCkTkN2y5Ytamxs1Lp161RZWan8/Hyt\nWLGiJ2YDILd7svv3V+if//xAQ4akKS6OA2PAhTr9V1FeXq7MzExJ0k033aQDBw6YDwXAfb/97Qv6\ny1+qdOpUnX7zm5WxHgdwUqc72fr6eiUnJ7fejo+PV3NzsxIS2v7SlJR+SkiI774JO5GaelWPvdbH\nEetrqyvru2XLFknS17/+detxLklzc1Prrxsb/6Pk5AQlJSXFcKKL8f61xfp2rtOQTU5OViQSab3d\n0tLSbsBK0qlTZ7pnsi5ITb1KtbWne+z1Pm5YX1tdXd/9+9+WJI0Zc4v1SF0WjUYVH58g6T+SpCAI\nVFq6WZmZX4vtYP+F968t1vcjHX2z0enh4rFjx6qsrEySVFlZqfT09O6bDICXjh49ojNnIufdd/Dg\nfp06VRejiQA3dRqyEydOVDgc1l133aVnnnlGTz75ZE/MBcBhb79dedF9LS0t2r79jZ4fBnBYp4eL\n4+Li9OMf/7gnZgHgiTvuuDvWIwBe4Jx7AACMcD1ZwHEu92QBdIydLAAARghZwHFcTxbwFyELOI7r\nyQL+ImQBADBCyAIAYISQBQDACCELAIARerKA4+jJAv5iJwsAgBFCFnAcPVnAX4Qs4Dh6soC/CFkA\nAIwQsgAAGCFkAQAwQsgCAGCEnizgOHqygL/YyQIAYISQBRxHTxbwFyELOI6eLOAvQhYAACOELAAA\nRghZAACMELIAABihJws4jp4s4C92sgAAGCFkAcfRkwX8RcgCjqMnC/iLkAUAwAghCwCAEUIWAAAj\nhCwAAEboyQKOoycL+IudLAAARghZwHH0ZAF/EbKA4+jJAv4iZAEAMELIAgBghJAFAMAIIQsAgBF6\nsoDj6MkC/mInCwCAEUIWcBw9WcBfhCzgOHqygL8IWQAAjBCyAAAYIWQBADBCyAIAYISeLOA4erKA\nv9jJAgBghJAFHEdPFvAXIQs4jp4s4C9CFgAAI4QsAABGCFkAAIwQsgAAGKEnCziOnizgL3ayAAAY\nIWQBx9GTBfxFyAKOoycL+IuQBQDACCELAIARQhYAACOELAAARujJAo6jJwv4i50sAABGCFnAcfRk\nAX8RsoDj6MkC/iJkAQAwQsgCAGCEkAUAwAghCwCAEXqygOPoyQL+YicLAIARQhZwHD1ZwF+ELOA4\nerKAvwhZAACMELIAABghZAEAMELIAgBghJ4s4Dh6soC/2MkCAGCEkAUcR08W8BchCziOnizgL0IW\nAAAjhCwAAEYIWQAAjBCyAAAYoScLOI6eLOAvdrIAABghZAHH0ZMF/EXIAo6jJwv4i5AFAMAIIQsA\ngBFCFgAAI4QsAABG6MkCjqMnC/iLnSwAAEYIWcBx9GQBfxGygOPoyQL+ImQBADBCyAIAYISQBQDA\nCCELAIARerKA4+jJAv5iJwsAgBFCFnAcPVnAX4Qs4Dh6soC/CFkAAIwQsgAAGCFkAQAwQsgCAGCE\nnizgOHqygL/YyQIAYISQBRxHTxbwFyELOI6eLOAvQhYAACOELAAARghZAACMELIAABihJws4jp4s\n4C92sgAAGCFkAcfRkwX8RcgCjqMnC/iLkAUAwAghCwCAEc4uBnBZGhoatHHjeh0/fkzDhg1XdvYU\nJSYmxnoswCmd7mRbWlr01FNPaebMmZo1a5aqq6t7Yi6vlJQUKytrggYNSlFW1gSVlBTHeiTAVEVF\nuTIyxqioqFCNjf9RUVGhxo+/URUV5bEeDXBKpzvZLVu2qLGxUevWrVNlZaXy8/O1YsWKnpjNCyUl\nxZoz577W24cPH2y9nZNzR6zGQi/iWk+2oaFBeXm5Kih4Xrffnt16f2npBuXl5Wr37v3saIH/FwqC\nIOjoNzzzzDO68cYblZ197h9TZmamtm1rv05QW3u6eyfsQGrqVT36em3Jypqgw4cPXnT/DTeM1htv\nvHnJzzdu3OjuGKtbxMWF1NLS4dsDV8DX9Y1EIjpzJqLU1IEXPVZb+4H69euv/v37x2Cy8/m6vr7w\neX3Lyw906/Olpl7V7mOd7mTr6+uVnJzcejs+Pl7Nzc1KSGj7S1NS+ikhIf4yxrw8Hf3hesI771S1\ne//lzBYXF7rSkbqVa/P0Nl1Z33HjxkmSysvdOBQbjTYrHA63OXs4HFY02uzM+8aVOXorX9e3J3Oj\n05BNTk5WJBJpvd3S0tJuwErSqVNnumeyLnBhJ5uePqrNnWx6+qjLmm337re7Y6xu4cL69mZdXd81\na1ZJklaufNF4oq559dUiFRUVat26kosemzkzRzNm5Gr69BkxmOx8vH9t+by+3T13R6Hd6YlPY8eO\nVVlZmSSpsrJS6enp3TdZL/DIIwvavH/evPk9PAnQM7Kzp+jQoYMqLd1w3v2lpRt06NBBZWdPidFk\ngHs63clOnDhR27dv11133aUgCLRkyZKemMsbH57ctGzZc3rnnSqlp4/SvHnzOekJvVZiYqJWry5U\nXl6uXnxxlcaMuVn79lXo0KGDWr26kJOegP/S6YlPl+rjduJTb8b62rrUw8UunmW8ceN6VVcfV1ra\nMOd6srx/bbG+H7miE58AoC1JSUlO/OwVcBkhCzjOtR0sgK7js4sBADBCyAKO43qygL8IWcBxXE8W\n8BchCwCAEUIWAAAjhCwAAEYIWQAAjHT7Jz4BAIBz2MkCAGCEkAUAwAghCwCAEUIWAAAjhCwAAEYI\nWQAAjHh1qbucnBwlJydLkgYPHqzs7Gz97Gc/U1JSkjIzM/Xggw/GeEK//fKXv9TWrVvV1NSk3Nxc\nZWRk6IknnlAoFNJ1112nH/3oR4qL4/uyy3Xh+t55552SpCVLlmj48OHKzc2N8YR+u3B9R48erZ/8\n5CeKj49XOBzW0qVLdc0118R6TC9duLY333yzFi1apCAINGzYMC1evFgJCV7FSc8JPHH27Nlg6tSp\nrbej0WiQlZUVvPvuu0EQBMGCBQuC3bt3x2o87+3cuTOYM2dOEI1Gg/r6+uAXv/hFMGfOnGDnzp1B\nEATBokWLgs2bN8d4Sn+1tb4nT54MZs+eHdx2223Byy+/HOsRvdbW+n7nO98JDh06FARBEBQWFgZL\nliyJ8ZR+amttH3jggeCtt94KgiAIHn/8cf5v6IA333pUVVWpoaFB9913n5qbmzV37lxdffXVGjJk\niCRp7Nix2rt3rz7/+c/HeFI//fnPf1Z6eroeeugh1dfXa+HChSoqKlJGRoYk6ctf/rK2b9+uiRMn\nxnhSP7W1vpFIRHPnzlVZWVmsx/NeW+s7c+ZMDRw4UJIUjUbVt2/fGE/pp7bW9sEHH1R8fLwaGxtV\nW1vbeoQRF/MmZBMTEzV79mzdeeedOn78uO6//361tLTor3/9q4YNG6aysjKNGjUq1mN669SpU3rv\nvfe0cuVK1dTU6IEHHlAQBAqFQpKk/v376/Tp0zGe0l9tre8f//hHDRkyhJDtBu2tryTt3btXL730\nktauXRvjKf3U3tr+/e9/17333qvk5GT+7+2ANyE7fPhwpaWlKRQKafjw4RowYICeeOIJPf300wqH\nw0pPT1dKSkqsx/TWgAEDNGLECIXDYY0YMUJ9+/bV+++/3/p4JBLR1VdfHcMJ/dbW+tbV1emTn/xk\nrEfrFdpb3127dmnFihV64YUX9IlPfCLWY3qpvbX9zGc+o82bN+uVV15Rfn6+li5dGutRneTNWSzF\nxcXKz8+XJJ04cUL19fXauXOnfv3rX2vVqlV699139cUvfjHGU/pr3Lhx2rZtm4Ig0IkTJ9TQ0KAJ\nEyZo165dkqSysjIOxV+BttZ3wIABsR6r12hrfcvKyvTSSy9pzZo1rT9WwqVra21/+MMf6vjx45LO\nHeXihMj2eXOBgMbGRj355JN67733FAqF9Oijj+ro0aNau3atEhMTNXnyZH33u9+N9ZheKygo0K5d\nuxQEgb7//e9r8ODBWrRokZqamjRixAgtXrxY8fHxsR7TWxeub2ZmpiRp+fLluuaaazi7+ApduL4L\nFizQoEGDWo/AjB8/Xg8//HCMp/TThWvbv39/FRQUqE+fPkpKStLixYtbf/6N83kTsgAA+IY9PgAA\nRghZAACMELIAABghZAEAMELIAgBgxJsPowB6i5qaGn3zm9/UZz/7WUnS2bNnNXLkSD311FOX/AH2\ny5Yt0+jRo3XbbbdZjArgClHhAXpYTU2N8vLytHXrVklSEAR67rnnVF5erpdffjnG0wHoTuxkgRgL\nhUKaO3euvvSlL6mqqkplZWUqLS1VNBrVrbfeqscee0z5+fkaOHCgZs+eLUl6+OGH9e1vf1tbt25V\nRkaGpk2bpueff147duzQv/71L6WkpGj58uVKTU3VrbfeqkmTJqm8vFzx8fH6+c9/riFDhujNN99U\nfn6+giDQtddeq2effVZJSUkqKCjQW2+9pWg0qmnTpumee+6J7QIBHuNnsoADwuGw0tLSVFVVpQMH\nDqi4uFivvfaaTpw4oddff11Tp07Vhg0bJEn19fXau3evvvKVr7R+fXV1tf72t7/pd7/7nTZt2qSh\nQ4dq/fr1kqTa2lpNmDBBr732msaPH6+1a9eqsbFRjz76qJYuXar169dr5MiRKikpUVFRkSSppKRE\nxcXF+tOf/qQ9e/b0+HoAvQU7WcARoVBIq1evVl1dnaZNmybp3M9rr732Wk2dOlWNjY2qrq5WRUWF\nvvrVryocDrd+bVpamh5//HG98sorOnbsmCorKzV06NDWxz/8CMfrrrtOe/bs0ZEjR/SpT31K119/\nvSRp/vz5ks7tkA8fPqydO3dKks6cOaMjR47wudXAZSJkAQc0Njbq2LFjuuWWWzR58mTde++9kqR/\n//vfrZ8XPWXKFG3cuFEVFRW6//77z/v6AwcOaMGCBbrnnns0adIkxcXF6b9Pt/jwWqqhUEhBEKhP\nnz7nff3p06cViUQUjUb12GOP6Rvf+IYkqa6uTv369TP7cwO9HYeLgRhraWnR8uXLNWbMGE2fPl2/\n//3vFYlE1NzcrIceekibNm2SJE2ePFkbN25UdXX1RTvL3bt3KyMjQ7m5ufrc5z6n7du3KxqNtvua\nw4cPV11dnY4ePSpJWrVqlQoLC/WFL3xBRUVFampqUiQS0d133619+/bZ/eGBXo6dLBADH3zwgaZO\nnSrpXMhef/31evbZZzVgwABVVVVpxowZikajyszMVE5OjiRp0KBBSklJ0U033aRQKHTe833rW9/S\n9773PU2ePFl9+vTRyJEjVVNT0+7r9+3bVz/96U+1cOFCNTU1aejQoSooKFA4HFZ1dbVycnLU3Nys\nadOm6ZZbbrFbCKCXo8IDAIARDhcDAGCEkAUAwAghCwCAEUIWAAAjhCwAAEYIWQAAjBCyAAAYIWQB\nADDyf2IvK/9TovU3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f2e70199978>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.compare_plot(df_comp_WAIC);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The empty circle represents the values of WAIC and the black error bars associated with them are the values of the standard deviation of WAIC. \n",
    "\n",
    "The value of the lowest WAIC is also indicated with a vertical dashed grey line to ease comparison with other WAIC values.\n",
    "\n",
    "The filled black dots are the in-sample deviance of each model, which for WAIC is  2 pWAIC from the corresponding WAIC value.\n",
    "\n",
    "For all models except the top-ranked one we also get a triangle indicating the value of the difference of WAIC between that model and the top model and a grey errobar indicating the standard error of the differences between the top-ranked WAIC and WAIC for each model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Leave-one-out Cross-validation (LOO)\n",
    "\n",
    "LOO cross-validation is an estimate of the out-of-sample predictive fit. In cross-validation, the data are repeatedly partitioned into training and holdout sets, iteratively fitting the model with the former and evaluating the fit with the holdout data. Vehtari et al. (2016) introduced an efficient computation of LOO from MCMC samples, which are corrected using Pareto-smoothed importance sampling (PSIS) to provide an estimate of point-wise out-of-sample prediction accuracy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/osvaldo/anaconda3/lib/python3.5/site-packages/pymc3/stats.py:208: UserWarning: Estimated shape parameter of Pareto distribution is\n",
      "        greater than 0.7 for one or more samples.\n",
      "        You should consider using a more robust model, this is\n",
      "        because importance sampling is less likely to work well if the marginal\n",
      "        posterior and LOO posterior are very different. This is more likely to\n",
      "        happen with a non-robust model and highly influential observations.\n",
      "  happen with a non-robust model and highly influential observations.\"\"\")\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "61.591020562529295"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pooled_loo = pm.loo(trace_p, pooled)\n",
    "    \n",
    "pooled_loo.LOO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/osvaldo/anaconda3/lib/python3.5/site-packages/pymc3/stats.py:208: UserWarning: Estimated shape parameter of Pareto distribution is\n",
      "        greater than 0.7 for one or more samples.\n",
      "        You should consider using a more robust model, this is\n",
      "        because importance sampling is less likely to work well if the marginal\n",
      "        posterior and LOO posterior are very different. This is more likely to\n",
      "        happen with a non-robust model and highly influential observations.\n",
      "  happen with a non-robust model and highly influential observations.\"\"\")\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "61.655191059793324"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hierarchical_loo  = pm.loo(trace_h, hierarchical)\n",
    "    \n",
    "hierarchical_loo.LOO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use `compare` with LOO."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>LOO</th>\n",
       "      <th>pLOO</th>\n",
       "      <th>dLOO</th>\n",
       "      <th>weight</th>\n",
       "      <th>SE</th>\n",
       "      <th>dSE</th>\n",
       "      <th>warning</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>61.591</td>\n",
       "      <td>0.909406</td>\n",
       "      <td>0</td>\n",
       "      <td>0.508021</td>\n",
       "      <td>2.13514</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>61.6552</td>\n",
       "      <td>1.12409</td>\n",
       "      <td>0.0641705</td>\n",
       "      <td>0.491979</td>\n",
       "      <td>1.9971</td>\n",
       "      <td>0.0258501</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       LOO      pLOO       dLOO    weight       SE        dSE warning\n",
       "1   61.591  0.909406          0  0.508021  2.13514          0       1\n",
       "0  61.6552   1.12409  0.0641705  0.491979   1.9971  0.0258501       1"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_comp_LOO = pm.compare((trace_h, trace_p), (hierarchical, pooled), ic='LOO')\n",
    "df_comp_LOO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The columns return the equivalent values for LOO, notice that in this example we get two warnings. Also notice that the order of the models is not the same as the one for WAIC.\n",
    "\n",
    "We can also plot the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdkAAAFXCAYAAADu/TSqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFUxJREFUeJzt3X2QlmW9B/Dv7sLKmx1oxA4WopagjaMmiTFJ1BRZkTJg\nqTQjo6eapsyXMK1OozZFijRajn9o6czpgIpuNJgeIRu1BkMzXEEDAcdJmaFOygmmZMGWfTl/eKRQ\nkGPutc/tvZ/PXz73ss/zey5Gvnvde3+fu6m3t7c3AECfa270AABQV0IWAAoRsgBQiJAFgEKELAAU\nImQBoJBBff2EW7a80NdPWWujRg3Ltm07Gj3GgFPXdV+06OYkydlnf67Bk+xdXde96qx7WaNHH7jP\nr9nJNtigQS2NHmFAsu6NYd0bw7o3jpAFgEL6/HQx0DhVPU0MA5WdLAAUImShRh5++ME8/PCDjR4D\n+D9CFmrk6ac35umnNzZ6DOD/CFkAKETIAkAhri6Gmti5c2cee6w9f/7znzNs2FsyffppGTJkSKPH\nggHt/7WTffzxx3P22WeXnqWSli5dkqlTJ2fMmFGZOnVyli5d0uiR4FVWr27PpEnH5dFHH01XV1fa\n2hbnxBOPzerV7Y0eDQa0/e5kb7rpptx1110ZOnRof8xTKUuXLskXvvBvux+vX79u9+OZMz/VqLFg\nDzt37sycObOzYMH38/GPT999fPnyezJnzuysWvWEHS00SFNvb2/va/2Be++9NxMmTMill16atra2\n/T5hf3x28cSJxxR/jST505/+O7t27XrV8cGDB+df/3VMn7xGc3NTenpe86+AAuq07h0dHdmxoyOj\nRx/8qq9t2fJ8hg0bnuHDhzdgsler07q/mVj3PbW3r+3T53utzy7e7072lFNOyebNm//fLzZq1LDi\nn5PZ3NxU9PlftreAffl4X87QX++HPdVl3bu7u9La2prm5qZMnDgxSdLe/tJp4tbW1nR3d1XqvVZp\nloHEuv/da4ViX+vzC5/6404Pq1b9rvhrJMnUqZOzfv26Vx1/97uPya9+9VCfvMbo0Qe6c1ED1Gnd\nf/rTtrS1Lc4ddyzdfReeG2/8cZLkzDNn5owzZuf0089o4IR/V6d1fzOx7nvq67VwF55/0kUXXbzX\n4xdeOLefJ4F9mz79tDz55LosX37PHseXL78nTz65LtOnn9agyQAVntfw8sVN1113bZ56akPGjz8q\nF14410VPVMqQIUOycOHizJkzOyNH/kve8Y6x+a//eilgFy5c7KInaKD9Xvj0ejkl8fo4jdMYdVz3\nnTt35t//fW7+/OetmTHj9Er2ZOu47m8G1r2sN3ThE/DmMHTo0JxwwksXPlXld7Aw0AlZqBH3k4Vq\nceETABQiZKFG3E8WqkXIQo24nyxUi5AFgEKELAAUImQBoBAhCwCF6MlCjejJQrXYyQJAIUIWakRP\nFqpFyEKN6MlCtQhZAChEyAJAIUIWAAoRsgBQiJ4s1IieLFSLnSwAFCJkoUb0ZKFahCzUiJ4sVIuQ\nBYBChCwAFCJkAaAQIQsAhejJQo3oyUK12MkCQCFCFmpETxaqRchCjejJQrUIWQAoRMgCQCFCFgAK\nEbIAUIieLNSInixUi50sABQiZKFG9GShWoQs1IieLFSLkAWAQoQsABQiZAGgECELAIXoyUKN6MlC\ntdjJAkAhQhZqRE8WqkXIQo3oyUK1CFkAKETIAkAhQhYAChGyAFCInizUiJ4sVIudLAAUImShRvRk\noVqELNSInixUi5AFgEKELAAUImQBoBAhCwCF6MlCjejJQrXYyQJAIUIWakRPFqpFyEKN6MlCtQhZ\nAChEyAJAIUIWAAoRsgBQiJ4s1IieLFSLnSwAFCJkoUb0ZKFahCzUiJ4sVIuQBYBChCwAFCJkAaAQ\nIQsAhejJQo3oyUK12MkCQCFCFmpETxaqRchCjejJQrUIWQAoRMgCQCFCFgAKEbIAUIieLNSInixU\ni50sABQiZKFG9GShWoQs1IieLFSLkAWAQoQsABQiZAGgECELAIXoyUKN6MlCtdjJAkAhQhZqRE8W\nqkXIQo3oyUK1CFkAKETIAkAhQhYAChGyAFCInizUiJ4sVIudLAAUImShRvRkoVqELNSInixUi5AF\ngEKELAAUImQBoBAhCwCF6MlCjejJQrXYyQJAIUIWakRPFqpFyEKN6MlCtQhZAChEyAJAIUIWAAoR\nsgBQiJ4s1IieLFSLnSwAFCJkoUb0ZKFahCzUiJ4sVIuQBYBChCwAFCJkAaAQIQsAhejJQo3oyUK1\n2MkCQCFCFmpETxaqRchCjejJQrUIWQAoRMgCQCFCFgAK2W+Fp6enJ9/61reycePGtLa2Zt68eRk3\nblx/zAb0keee+++8+OKLGTfu8EaPAgPKfkP2vvvuS2dnZ+64446sWbMm8+fPzw033NAfswGv0756\nsk88sTr/8z/PZ+zYcWludgIL+st+Q7a9vT1TpkxJkhx//PFZu3Zt8aGAvvOf//mjdHRsT5L8x3/c\nmNbWVh9aAf1kvyG7ffv2jBgxYvfjlpaWdHV1ZdCgvX/rqFHDMmhQS99NOACMHn1go0cYkOq47vfd\nd1+S5CMf+cjuY11du3b/d2fn39LaOrih772O6/5mYN0bY78hO2LEiHR0dOx+3NPTs8+ATZJt23b0\nzWQDxOjRB2bLlhcaPcaAU9d1f+KJ3yVJjjvupCRJd3d3WloGJflbkqS3tzeHHfbOhr33uq571Vn3\nsl7rB5j9/nLmhBNOyIoVK5Ika9asyfjx4/tuMqCop5/emB07OvY4tm7dE9m2bWuDJoKBZb8hO23a\ntLS2tuass87KVVddlW984xv9MRfQB373uzWvOtbT05OVK3/V/8PAALTf08XNzc359re/3R+zAH3s\nU5/6TKNHgAHNtfwAUIj7yUKNqOZAtdjJAkAhQhZqxP1koVqELNSI+8lCtQhZAChEyAJAIUIWAAoR\nsgBQiJ4s1IieLFSLnSwAFCJkoUb0ZKFahCzUiJ4sVIuQBYBChCwAFCJkAaAQIQsAhejJQo3oyUK1\n2MkCQCFCFmpETxaqRchCjejJQrUIWQAoRMgCQCFCFgAKEbIAUIieLNSInixUi50sABQiZKFG9GSh\nWoQs1IieLFSLkAWAQoQsABQiZAGgECELAIXoyUKN6MlCtdjJAkAhQhZqRE8WqkXIQo3oyUK1CFkA\nKETIAkAhQhYAChGyAFCInizUiJ4sVIudLAAUImShRvRkoVqELNSInixUi5AFgEKELAAUImQBoBAh\nCwCF6MlCjejJQrXYyQJAIUIWakRPFqpFyEKN6MlCtQhZAChEyAJAIUIWAAoRsgBQiJ4s1IieLFSL\nnSwAFCJkoUb0ZKFahCzUiJ4sVIuQBYBChCwAFCJkAaAQIQsAhejJQo3oyUK12MkCQCFCFmpETxaq\nRchCjejJQrUIWQAoRMgCQCFCFgAKEbIAUIieLNSInixUi50sABQiZKFG9GShWoQs1IieLFSLkAWA\nQoQsABQiZAGgECELAIXoyUKN6MlCtdjJAkAhQhZqRE8WqkXIQo3oyUK1CFkAKETIAkAhQhYAChGy\nAFCInizUiJ4sVIudLAAUImShRvRkoVqELNSInixUi5AFgEKELAAUImQBoBAhCwCF6MlCjejJQrXY\nyQJAIUIWakRPFqpFyEKN6MlCtQhZAChEyAJAIUIWAAoRsgBQiJ4s1IieLFSLnSwAFCJkoUb0ZKFa\nhCzUiJ4sVIuQBYBChCwAFCJkAaAQIQsAhejJQo3oyUK12MkCQCFCFmpETxaqRchCjejJQrUIWQAo\nRMgCQCFCFgAKEbIAUIieLNSInixUi50sABQiZKFG9GShWoQs1IieLFSLkAWAQoQsABTi6mKouZ07\nd2bZsrvz7LPP5LDDDs/06adlyJAhjR4LBoT97mR7enpy+eWX58wzz8zZZ5+dTZs29cdcWbp0SaZO\nnZwxY0Zl6tTJWbp0Sb+8LtTJ6tXtmTTpuLS1LU5n59/S1rY4J554bFavbm/0aDAg7Hcne99996Wz\nszN33HFH1qxZk/nz5+eGG24oOtTSpUvyhS/82+7H69ev2/145sxPFX1teDP7x57szp07M2fO7CxY\n8P18/OPTdx9fvvyezJkzO6tWPWFHC4U19fb29r7WH7jqqqty7LHHZvr0l/4nnTJlSh58cN8VgS1b\nXnjDQ02dOjnr16971fF3v/uY/OpXD73h5/9HEyce06fP93o1Nzelp+c1/wooYCCse0dHR3bs6Mjo\n0Qe/6mtbtjyfYcOGZ/jw4f0600BY9yoaiOve3r62315r9OgD9/m1/e5kt2/fnhEjRux+3NLSkq6u\nrgwatPdvHTVqWAYNavknxvy7p57asM/jr/Vm/hnNzU19+nxv1hkGojqu+8SJE5Mk7e3t6e7uSmtr\n617fZ2tra7q7uxqyBnVc9zeDgbbufZ0V/6z9huyIESPS0dGx+3FPT88+AzZJtm3b8YaHGj/+qL3u\nZMePP6pPdsr/aNWq3/Xp871eo0cf2Ofvif2r67ovWnRzkuTGG3+cn/60LW1ti3PHHUtf9efOPHNm\nzjhjdk4//Yx+na+u6151A3Hd+/P9vlag7/fCpxNOOCErVqxIkqxZsybjx4/vu8n24aKLLt7r8Qsv\nnFv8taEupk8/LU8+uS7Ll9+zx/Hly+/Jk0+uy/TppzVoMhg49ruTnTZtWlauXJmzzjorvb29ufLK\nK4sP9fLFTdddd22eempDxo8/KhdeONdFT/A6DBkyJAsXLs6cObPz4x/fnOOOe08ef3x1nnxyXRYu\nXOyiJ+gH+73w6fUaaKck3qiBeBqnCuq67i+fLn7lVcbLlt2dTZuezbhxhzW0J1vXda86617WG7rw\nCXhzGzp0aL//7hV4iZCFGnE/WagWn10MAIUIWagR95OFahGyUCPuJwvVImQBoBAhCwCFCFkAKETI\nAkAhff6JTwDAS+xkAaAQIQsAhQhZAChEyAJAIUIWAAoRsgBQiFvd9bMf/vCHeeCBB7Jr167Mnj07\nkyZNyte//vU0NTXlyCOPzBVXXJHmZj/79LVXrvunP/3pJMmVV16Zww8/PLNnz27whPXzyjU/5phj\n8p3vfCctLS1pbW3N1VdfnYMOOqjRY9bOK9f9Pe95Ty677LL09vbmsMMOy7x58zJokH/6+4t/zfvR\nI488ktWrV2fx4sVZtGhR/vSnP+Wqq67KRRddlNtuuy29vb25//77Gz1m7ext3bdu3ZrPfe5zeeCB\nBxo9Xi3tbc2/+93v5rLLLsuiRYsybdq03HTTTY0es3b2tu7XXntt5s6dm9tvvz1J8stf/rLBUw4s\nfpzpR7/+9a8zfvz4nHfeedm+fXsuvfTStLW1ZdKkSUmSD3zgA1m5cmWmTZvW4EnrZW/r3tHRkfPP\nPz8rVqxo9Hi1tLc1P/PMM3PwwQcnSbq7u3PAAQc0eMr62du6f+lLX0pLS0s6OzuzZcuWjBgxotFj\nDihCth9t27Ytf/zjH3PjjTdm8+bN+eIXv5je3t40NTUlSYYPH54XXnihwVPWz97W/ec//3nGjh0r\nZAvZ15onyWOPPZZbbrklt956a4OnrJ99rfsf/vCHnHvuuRkxYkSOOuqoRo85oDhd3I9GjhyZk08+\nOa2trTniiCNywAEH7BGqHR0dectb3tLACetpb+u+devWRo9Va/ta82XLluWKK67Ij370o7z1rW9t\n9Ji1s691f/vb355f/OIXmT17dubPn9/oMQcUIduPJk6cmAcffDC9vb157rnnsnPnzkyePDmPPPJI\nkmTFihV573vf2+Ap62dv6z5y5MhGj1Vre1vzFStW5JZbbsmiRYsyduzYRo9YS3tb929+85t59tln\nk7x0tsyFlf3LDQL62YIFC/LII4+kt7c3X/nKV/KOd7wjl112WXbt2pUjjjgi8+bNS0tLS6PHrJ1X\nrvuUKVOSJNdff30OOuggVxcX8Mo1v/jiizNmzJjdZ2tOPPHEXHDBBQ2esn5eue7Dhw/PggULMnjw\n4AwdOjTz5s3b/btxyhOyAFCI8wYAUIiQBYBChCwAFCJkAaAQIQsAhfjEJ+hnmzdvzsc+9rG8853v\nTJK8+OKLmTBhQi6//PLX/YH51113XY455ph8+MMfLjEq8Aap8EA/27x5c+bMmbP75gS9vb259tpr\n097enttuu63B0wF9yU4WGqypqSnnn39+3v/+92fDhg1ZsWJFli9fnu7u7px88sm55JJLMn/+/Bx8\n8MH57Gc/myS54IIL8slPfjIPPPBAJk2alFmzZuX73/9+Hn744fzlL3/JqFGjcv3112f06NE5+eST\nc8opp6S9vT0tLS35wQ9+kLFjx+ahhx7K/Pnz09vbm0MOOSTXXHNNhg4dmgULFuS3v/1turu7M2vW\nrJxzzjmNXSB4E/M7WaiA1tbWjBs3Lhs2bMjatWuzZMmS3HnnnXnuuedy1113ZcaMGbnnnnuSJNu3\nb89jjz2WD37wg7u/f9OmTfn973+f22+/Pffee28OPfTQ3H333UmSLVu2ZPLkybnzzjtz4okn5tZb\nb01nZ2e++tWv5uqrr87dd9+dCRMmZOnSpWlra0uSLF26NEuWLMn999+fRx99tN/XA+rCThYqoqmp\nKQsXLszWrVsza9asJC/9vvaQQw7JjBkz0tnZmU2bNmX16tX50Ic+lNbW1t3fO27cuHzta1/LT37y\nkzzzzDNZs2ZNDj300N1ff/ljJI888sg8+uij2bhxY972trfl6KOPTpLMnTs3yUs75PXr1+c3v/lN\nkmTHjh3ZuHGjz9SGf5KQhQro7OzMM888k5NOOimnnnpqzj333CTJX//6192fZX3aaadl2bJlWb16\ndT7/+c/v8f1r167NxRdfnHPOOSennHJKmpub84+XW7x879ampqb09vZm8ODBe3z/Cy+8kI6OjnR3\nd+eSSy7JRz/60STJ1q1bM2zYsGLvG+rO6WJosJ6enlx//fU57rjjcvrpp+dnP/tZOjo60tXVlfPO\nOy/33ntvkuTUU0/NsmXLsmnTplftLFetWpVJkyZl9uzZede73pWVK1emu7t7n695+OGHZ+vWrXn6\n6aeTJDfffHMWL16c973vfWlra8uuXbvS0dGRz3zmM3n88cfLvXmoOTtZaIDnn38+M2bMSPJSyB59\n9NG55pprMnLkyGzYsCFnnHFGuru7M2XKlMycOTNJMmbMmIwaNSrHH398mpqa9ni+T3ziE/nyl7+c\nU089NYMHD86ECROyefPmfb7+AQcckO9973u59NJLs2vXrhx66KFZsGBBWltbs2nTpsycOTNdXV2Z\nNWtWTjrppHILATWnwgMAhThdDACFCFkAKETIAkAhQhYAChGyAFCIkAWAQoQsABQiZAGgkP8FgTvf\n9H5Up+UAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f2e31cd5b70>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.compare_plot(df_comp_LOO);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpretation\n",
    "\n",
    "Though we might expect the hierarchical model to outperform a complete pooling model, there is little to choose between the models in this case, giving that both models gives very similar values of the information criteria. This is more clearly appreciated when we take into account the uncertainty (in terms of standard errors) of WAIC and LOO.\n",
    "\n",
    "## Reference\n",
    "\n",
    "[Gelman, A., Hwang, J., & Vehtari, A. (2014). Understanding predictive information criteria for Bayesian models. Statistics and Computing, 24(6), 997–1016.](http://doi.org/10.1007/s11222-013-9416-2)\n",
    "\n",
    "[Vehtari, A, Gelman, A, Gabry, J. (2016). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing](http://link.springer.com/article/10.1007/s11222-016-9696-4)"
   ]
  }
 ],
 "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}