{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sampler statistics\n",
    "\n",
    "When checking for convergence or when debugging a badly behaving\n",
    "sampler, it is often helpful to take a closer look at what the\n",
    "sampler is doing. For this purpose some samplers export\n",
    "statistics for each generated sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sb\n",
    "import pandas as pd\n",
    "import pymc3 as pm \n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a minimal example we sample from a standard normal distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = pm.Model()\n",
    "with model:\n",
    "    mu1 = pm.Normal(\"mu1\", mu=0, sigma=1, shape=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [mu1]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "        <style>\n",
       "            /* Turns off some styling */\n",
       "            progress {\n",
       "                /* gets rid of default border in Firefox and Opera. */\n",
       "                border: none;\n",
       "                /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "                background-size: auto;\n",
       "            }\n",
       "            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "                background: #F44336;\n",
       "            }\n",
       "        </style>\n",
       "      <progress value='6000' class='' max='6000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [6000/6000 00:07<00:00 Sampling 2 chains, 0 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 9 seconds.\n"
     ]
    }
   ],
   "source": [
    "with model:\n",
    "    step = pm.NUTS()\n",
    "    trace = pm.sample(2000, tune=1000, init=None, step=step, cores=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NUTS provides the following statistics:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'depth',\n",
       " 'diverging',\n",
       " 'energy',\n",
       " 'energy_error',\n",
       " 'max_energy_error',\n",
       " 'mean_tree_accept',\n",
       " 'model_logp',\n",
       " 'step_size',\n",
       " 'step_size_bar',\n",
       " 'tree_size',\n",
       " 'tune'}"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trace.stat_names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- `mean_tree_accept`: The mean acceptance probability for the tree that generated this sample. The mean of these values across all samples but the burn-in should be approximately `target_accept` (the default for this is 0.8).\n",
    "- `diverging`: Whether the trajectory for this sample diverged. If there are many diverging samples, this usually indicates that a region of the posterior has high curvature. Reparametrization can often help, but you can also try to increase `target_accept` to something like 0.9 or 0.95.\n",
    "- `energy`: The energy at the point in phase-space where the sample was accepted. This can be used to identify posteriors with problematically long tails. See below for an example.\n",
    "- `energy_error`: The difference in energy between the start and the end of the trajectory. For a perfect integrator this would always be zero.\n",
    "- `max_energy_error`: The maximum difference in energy along the whole trajectory.\n",
    "- `depth`: The depth of the tree that was used to generate this sample\n",
    "- `tree_size`: The number of leafs of the sampling tree, when the sample was accepted. This is usually a bit less than $2 ^ \\text{depth}$. If the tree size is large, the sampler is using a lot of leapfrog steps to find the next sample. This can for example happen if there are strong correlations in the posterior, if the posterior has long tails, if there are regions of high curvature (\"funnels\"), or if the variance estimates in the mass matrix are inaccurate. Reparametrisation of the model or estimating the posterior variances from past samples might help.\n",
    "- `tune`: This is `True`, if step size adaptation was turned on when this sample was generated.\n",
    "- `step_size`: The step size used for this sample.\n",
    "- `step_size_bar`: The current best known step-size. After the tuning samples, the step size is set to this value. This should converge during tuning.\n",
    "- `model_logp`: The model log-likelihood for this sample."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the name of the statistic does not clash with the name of one of the variables, we can use indexing to get the values. The values for the chains will be concatenated.\n",
    "\n",
    "We can see that the step sizes converged after the 1000 tuning samples for both chains to about the same value. The first 2000 values are from chain 1, the second 2000 from chain 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fbc86827550>]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD4CAYAAADlwTGnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAYrklEQVR4nO3df5Bd9X3e8feDZBEnAoOljQZYDciGjL3xaABvcRiDpdrjiWQnEET/ALspxEkYR1Z/xNXU0jDDOMKMTCoa7MLYo7oqVn6YgGJnCFEKFITp1Alm+SGBjAVrGaqVaFibqCnFRQY9/eN+Vxxdr7R3d6/26ug8r5kdnfM933PO55xd3WfP+d6zV7aJiIjmOanXBURERG8kACIiGioBEBHRUAmAiIiGSgBERDTU7F4XMBnz58/3Oeec0+syIiJq5fHHH/+R7b729o4CQNIy4EvALOBrtr/YtvxsYBPQB7wC/HPbI5XlpwLfA/7S9qrSNge4DVgKHASut/0XR6vjnHPOYWhoqJOSIyKikPTieO0T3gKSNAu4HVgODABXSxpo67YB2Gx7MbAOWN+2/Ebgkba264GXbf9S2e63J6olIiK6p5MxgIuAYdu7bR8A7gQub+szADxUprdVl0t6P7AAuL9tnU9RgsL2Qds/mnz5ERExVZ0EwFnAnsr8SGmr2g6sKNNXAKdImifpJOAWYHW1s6TTyuSNkp6QdLekBePtXNJ1koYkDY2OjnZQbkREdKJb7wJaDSyR9CSwBNgLvAmsBLZWxwOK2UA/8B3bFwJ/S+s20s+wvdH2oO3Bvr6fGcOIiIgp6mQQeC+wsDLfX9oOsb2PcgUgaS5wpe39ki4GLpW0EpgLzJH0KrAWeA34ZtnE3cBvT+dAIiJicjoJgMeA8yQtovXCfxXwiWoHSfOBV2wfpPXivgnA9icrfa4FBm2vKfN/ResdQA8BH6H1LqGIiJghE94Csv0GsAq4D3gWuMv2TknrJF1Wui0Fdkl6jtaA700d7PtzwOcl7QB+E/i3U6g/IiKmSHX6c9CDg4OeynMAd/yPH/LK/z1wDCqKaPnAu+bxwXPn97qMiHFJetz2YHt7rZ4Enqo/++7/5PmXX+11GXGCsuF9u17m3n95aa9LiZiURgTA/b+/pNclxAnsdzcPMfIPP+l1GRGTlj8GFxHRUAmAiIiGSgBERDRUAiCiC+r0brqIMQmAiIiGSgBETJN6XUDEFCUAIiIaKgEQEdFQCYCIiIZKAERENFQCIGKalFHgqKkEQEREQyUAIiIaKgEQ0QV5EDjqKAEQEdFQCYCIaVKeBY6aSgBERDRUAiAioqE6CgBJyyTtkjQsac04y8+W9KCkHZIeltTftvxUSSOSbhtn3XskPTP1Q4joPZNR4KifCQNA0izgdmA5MABcLWmgrdsGYLPtxcA6YH3b8huBR8bZ9gogn9YeEdEDnVwBXAQM295t+wBwJ3B5W58B4KEyva26XNL7gQXA/dUVJM0FPgt8YWqlRxwf8iRw1FUnAXAWsKcyP1LaqrYDK8r0FcApkuZJOgm4BVg9znZvLMteO9rOJV0naUjS0OjoaAflRkREJ7o1CLwaWCLpSWAJsBd4E1gJbLU9Uu0s6Xzg3ba/NdGGbW+0PWh7sK+vr0vlRkTE7A767AUWVub7S9shtvdRrgDKrZ0rbe+XdDFwqaSVwFxgjqRXgReBQUkvlBp+UdLDtpdO83gieiJPAkcddRIAjwHnSVpE64X/KuAT1Q6S5gOv2D4IrAU2Adj+ZKXPtcCg7bF3EX2ltJ8D3JsX/4iImTXhLSDbbwCrgPuAZ4G7bO+UtE7SZaXbUmCXpOdoDfjedIzqjTjuZBA46qqTKwBsbwW2trXdUJneAmyZYBt3AHeM0/4C8L5O6oiIiO7Jk8AREQ2VAIjogowBRx0lACIiGioBEDFN+XPQUVcJgIiIhkoAREQ0VAIgogucR4GjhhIAERENlQCImK6MAUdNJQAiIhoqARAR0VAJgIguyBBw1FECICKioRIAEdOUMeCoqwRARERDJQAiIhoqARDRDRkFjhpKAERENFQCIGKalA8FjprqKAAkLZO0S9KwpDXjLD9b0oOSdkh6WFJ/2/JTJY1Iuq3M/7ykv5b0fUk7JX2xO4cTERGdmjAAJM0CbgeWAwPA1ZIG2rptADbbXgysA9a3Lb8ReKR9HdvvAS4APihp+RTqj4iIKerkCuAiYNj2btsHgDuBy9v6DAAPlelt1eWS3g8sAO4fa7P9mu1tZfoA8ARw2FVDRJ1kDDjqqJMAOAvYU5kfKW1V24EVZfoK4BRJ8ySdBNwCrD7SxiWdBvw68GCnRUdExPR1axB4NbBE0pPAEmAv8CawEthqe2S8lSTNBr4BfNn27iP0uU7SkKSh0dHRLpUb0T0ZAo66mt1Bn73Awsp8f2k7xPY+yhWApLnAlbb3S7oYuFTSSmAuMEfSq7bHBpI3As/bvvVIO7e9sfRjcHAwV9oREV3SSQA8BpwnaRGtF/6rgE9UO0iaD7xi+yCwFtgEYPuTlT7XAoNjL/6SvgC8A/id6R9GRERM1oS3gGy/AawC7gOeBe6yvVPSOkmXlW5LgV2SnqM14HvT0bZZ3iZ6Pa3B4yckPSUpQRC1lc8Ejjrq5AoA21uBrW1tN1SmtwBbJtjGHcAdZXqE3DqNiOipPAkcMU15EDjqKgEQEdFQCYCIiIZKAER0QYaAo44SABERDZUAiJimjAFHXSUAIiIaKgEQEdFQCYCILsiDwFFHCYCIiIZKAERMUz4TOOoqARAR0VAJgIiIhkoARHSB8yxw1FACICKioRIAEdOUIeCoqwRARERDJQAiIhoqARDRBXkSOOooARAR0VAdBYCkZZJ2SRqWtGac5WdLelDSDkkPS+pvW36qpBFJt1Xa3i/p6bLNLyuPU0Zd5Sc3amrCAJA0C7gdWA4MAFdLGmjrtgHYbHsxsA5Y37b8RuCRtravAL8LnFe+lk26+oiImLJOrgAuAoZt77Z9ALgTuLytzwDwUJneVl0u6f3AAuD+StsZwKm2/862gc3Ab0z5KCIiYtI6CYCzgD2V+ZHSVrUdWFGmrwBOkTRP0knALcDqcbY5MsE2AZB0naQhSUOjo6MdlBsREZ3o1iDwamCJpCeBJcBe4E1gJbDV9sjRVj4a2xttD9oe7Ovr6061EV2WdwFFHc3uoM9eYGFlvr+0HWJ7H+UKQNJc4Erb+yVdDFwqaSUwF5gj6VXgS2U7R9xmRF0oo8BRU50EwGPAeZIW0XqRvgr4RLWDpPnAK7YPAmuBTQC2P1npcy0waHtNmf9HSb8CPAr8C+A/TvtoIiKiYxPeArL9BrAKuA94FrjL9k5J6yRdVrotBXZJeo7WgO9NHex7JfA1YBj4AfA3ky8/IiKmqpMrAGxvBba2td1Qmd4CbJlgG3cAd1Tmh4D3dV5qRER0U54EjohoqARAxDTlGfaoqwRARERDJQAiIhoqARAR0VAJgIgucB4FjhpKAERMU8aAo64SABERDZUAiIhoqARARERDJQAiuiBDwFFHCYCIacqTwFFXCYCIiIZKAERENFQCICKioRIAEV2QB4GjjhIAEdOUzwSOukoAREQ0VAIgIqKhEgAREQ3VUQBIWiZpl6RhSWvGWX62pAcl7ZD0sKT+SvsTkp6StFPSpyvrXC3p6bLOf5U0v3uHFTGznGeBo4YmDABJs4DbgeXAAHC1pIG2bhuAzbYXA+uA9aX9JeBi2+cDHwDWSDpT0mzgS8A/LevsAFZ144AiIqIznVwBXAQM295t+wBwJ3B5W58B4KEyvW1sue0Dtl8v7SdX9qfy9QuSBJwK7JvyUUT0UP4URNRVJwFwFrCnMj9S2qq2AyvK9BXAKZLmAUhaKGlH2cbNtvfZ/inwe8DTtF74B4D/PN7OJV0naUjS0OjoaIeHFRERE+nWIPBqYImkJ4ElwF7gTQDbe8ptnnOBayQtkPQ2WgFwAXAmrVtAa8fbsO2NtgdtD/b19XWp3IiImN1Bn73Awsp8f2k7xPY+yhWApLnAlbb3t/eR9AxwKfBiaftBWecu4GcGlyPqIk8CRx11cgXwGHCepEWS5gBXAfdUO0iaL2lsW2uBTaW9X9Lby/TpwCXALloBMiBp7Ff6jwLPTvdgIiKicxNeAdh+Q9Iq4D5gFrDJ9k5J64Ah2/cAS4H1kgw8AnymrP5e4JbSLmCD7acBJP0B8Iikn9K6Iri2q0cWMUMyCBx11cktIGxvBba2td1Qmd4CbBlnvQeAxUfY5leBr06m2IiI6J48CRwR0VAJgIguyBhw1FECICKioRIAEdOWUeCopwRARERDJQAiIhoqARDRBXkSOOooARAR0VAJgIhpypPAUVcJgIiIhkoAREQ0VAIgoisyChz1kwCIiGioBEDENGUMOOoqARAR0VAJgIiIhkoARHRBngSOOkoAREQ0VAIgYpryJHDUVUcBIGmZpF2ShiWtGWf52ZIelLRD0sOS+ivtT0h6StJOSZ+urDNH0kZJz0n6vqQru3dYERExkQk/FF7SLOB24KPACPCYpHtsf6/SbQOw2fbXJX0YWA/8JvAScLHt1yXNBZ4p6+4Drgdetv1Lkk4C3tndQ4uIiKOZMACAi4Bh27sBJN0JXA5UA2AA+GyZ3gb8JYDtA5U+J3P4FcengPeUfgeBH02h/ojjQsaAo446uQV0FrCnMj9S2qq2AyvK9BXAKZLmAUhaKGlH2cbNtvdJOq30vbHcIrpb0oLxdi7pOklDkoZGR0c7PKyIiJhItwaBVwNLJD0JLAH2Am8C2N5jezFwLnBNeaGfDfQD37F9IfC3tG4j/QzbG20P2h7s6+vrUrkR3aM8Cxw11UkA7AUWVub7S9shtvfZXmH7Alr39rG9v70P8AxwKfBj4DXgm2Xx3cCFUzmAiIiYmk4C4DHgPEmLJM0BrgLuqXaQNL8M5AKsBTaV9n5Jby/TpwOXALtsG/grYGlZ5yMcPqYQERHH2ISDwLbfkLQKuA+YBWyyvVPSOmDI9j20XsjXSzLwCPCZsvp7gVtKu4ANtp8uyz4H/LGkW4FR4Le6eFwRM8p5FDhqqJN3AWF7K7C1re2GyvQWYMs46z0ALD7CNl8EPjSZYiMionvyJHDENOVJ4KirBEBEREMlACIiGioBENEFGQKOOkoAREQ0VAIgYpoyBhx1lQCIiGioBEBEREMlACK6IA8CRx0lACIiGioBEDFNyqPAUVMJgIiIhkoAREQ0VAIgogvy56CjjhIAERENlQCIiGioBEBEREMlACIiGioBENEFGQKOOkoAREQ0VEcBIGmZpF2ShiWtGWf52ZIelLRD0sOS+ivtT0h6StJOSZ8eZ917JD0z/UOJ6I08CBx1NWEASJoF3A4sBwaAqyUNtHXbAGy2vRhYB6wv7S8BF9s+H/gAsEbSmZVtrwBenfZRRETEpHVyBXARMGx7t+0DwJ3A5W19BoCHyvS2seW2D9h+vbSfXN2fpLnAZ4EvTL38iIiYqk4C4CxgT2V+pLRVbQdWlOkrgFMkzQOQtFDSjrKNm23vK/1uBG4BXjvaziVdJ2lI0tDo6GgH5Ub0QEaBo4a6NQi8Glgi6UlgCbAXeBPA9p5ya+hc4BpJCySdD7zb9rcm2rDtjbYHbQ/29fV1qdyIiJjdQZ+9wMLKfH9pO6T8Vr8CDt3audL2/vY+ZbD3UqAPGJT0QqnhFyU9bHvpFI8jomeUTwWOmurkCuAx4DxJiyTNAa4C7ql2kDRf0ti21gKbSnu/pLeX6dOBS4Bdtr9i+0zb55S25/LiHxExsyYMANtvAKuA+4Bngbts75S0TtJlpdtSYJek54AFwE2l/b3Ao5K2A98GNth+usvHEBERU9DJLSBsbwW2trXdUJneAmwZZ70HgMUTbPsF4H2d1BFxvMoYcNRRngSOiGioBEDENOVJ4KirBEBEREMlACIiGioBENEF+UzgqKMEQEREQyUAIqYpY8BRVwmAiIiGSgBERDRUAiAioqESABFdkPcARR0lACKmKU8CR10lACIiGioBEBHRUAmAiIiGSgBEdEH+EkTUUQIgYpqUUeCoqQRARERDJQAiIhqqowCQtEzSLknDktaMs/xsSQ9K2iHpYUn9lfYnJD0laaekT5f2n5f015K+X9q/2N3DioiIiUwYAJJmAbcDy4EB4GpJA23dNgCbbS8G1gHrS/tLwMW2zwc+AKyRdObYOrbfA1wAfFDS8mkfTUSPOM8CRw11cgVwETBse7ftA8CdwOVtfQaAh8r0trHltg/Yfr20nzy2P9uv2d421gd4AuifzoFE9EqGgKOuOgmAs4A9lfmR0la1HVhRpq8ATpE0D0DSQkk7yjZutr2vuqKk04BfBx6cfPkRETFV3RoEXg0skfQksATYC7wJYHtPuTV0LnCNpAVjK0maDXwD+LLt3eNtWNJ1koYkDY2Ojnap3IiI6CQA9gILK/P9pe0Q2/tsr7B9AXB9advf3gd4Bri00rwReN72rUfaue2NtgdtD/b19XVQbkREdGJ2B30eA86TtIjWC/9VwCeqHSTNB16xfRBYC2wq7f3Aj23/RNLpwCXAH5VlXwDeAfxOl44lomf+308P8tH/8O1elxEnsHv/1SWcPHtWV7c5YQDYfkPSKuA+YBawyfZOSeuAIdv3AEuB9ZIMPAJ8pqz+XuCW0i5a7/x5ugTD9cD3gSfKk5S32f5aV48uYgZ8fPEZjOz/Cc7fg4hjSMfg7Qaq0w/t4OCgh4aGel1GREStSHrc9mB7e54EjohoqARARERDJQAiIhoqARAR0VAJgIiIhkoAREQ0VAIgIqKhEgAREQ1VqwfBJI0CL05x9fnAj7pYTrekrslJXZOTuibnRK3rbNs/88fUahUA0yFpaLwn4XotdU1O6pqc1DU5Tasrt4AiIhoqARAR0VBNCoCNvS7gCFLX5KSuyUldk9OouhozBhAREYdr0hVARERUJAAiIhrqhA8AScsk7ZI0LGlND/b/gqSnJT0laai0vVPSA5KeL/+eXtol6cul1h2SLuxyLZskvSzpmUrbpGuRdE3p/7yka45RXZ+XtLect6ckfayybG2pa5ekX620d+17LWmhpG2Svidpp6R/Xdp7er6OUlevz9fPSfqupO2lrj8o7YskPVr28eeS5pT2k8v8cFl+zkT1drmuOyT9sHK+zi/tM/ZzX7Y5S9KTku4t8zN7vmyfsF+0PsLyB8C7gDnAdmBghmt4AZjf1vaHwJoyvQa4uUx/DPgbWh+f+SvAo12u5UPAhcAzU60FeCewu/x7epk+/RjU9Xlg9Th9B8r38WRgUfn+zur29xo4A7iwTJ8CPFf23dPzdZS6en2+BMwt028DHi3n4S7gqtL+VeD3yvRK4Ktl+irgz49W7zGo6w7gn43Tf8Z+7st2Pwv8GXBvmZ/R83WiXwFcBAzb3m37AHAncHmPa4JWDV8v018HfqPSvtktfwecJumMbu3U9iPAK9Os5VeBB2y/YvsfgAeAZcegriO5HLjT9uu2fwgM0/o+d/V7bfsl20+U6f8DPAucRY/P11HqOpKZOl+2/WqZfVv5MvBhYEtpbz9fY+dxC/ARSTpKvd2u60hm7Oderc9G/zjwtTIvZvh8negBcBawpzI/wtH/sxwLBu6X9Lik60rbAtsvlen/BSwo072od7K1zGSNq8pl+KaxWy29qKtcbl9A67fH4+Z8tdUFPT5f5XbGU8DLtF4gfwDst/3GOPs4tP+y/H8D82aiLttj5+umcr7+SNLJ7XW17f9YfB9vBf4dcLDMz2OGz9eJHgDHg0tsXwgsBz4j6UPVhW5dxx0X78U9nmoBvgK8GzgfeAm4pRdFSJoL/AXwb2z/Y3VZL8/XOHX1/HzZftP2+UA/rd9C3zPTNYynvS5J7wPW0qrvn9C6rfO5maxJ0q8BL9t+fCb32+5ED4C9wMLKfH9pmzG295Z/Xwa+Res/xt+P3dop/75cuvei3snWMiM12v778h/3IPCfeOuydsbqkvQ2Wi+yf2r7m6W55+drvLqOh/M1xvZ+YBtwMa1bKLPH2ceh/Zfl7wB+PEN1LSu30mz7deC/MPPn64PAZZJeoHX77cPAl5jp8zWdAYzj/QuYTWuwZhFvDXT98gzu/xeAUyrT36F13/Dfc/hA4h+W6Y9z+ADUd49BTedw+GDrpGqh9dvSD2kNhJ1ept95DOo6ozL9+7TucwL8MocPeu2mNaDZ1e91Oe7NwK1t7T09X0epq9fnqw84rUy/HfjvwK8Bd3P4oObKMv0ZDh/UvOto9R6Dus6onM9bgS/24ue+bHspbw0Cz+j56uqLy/H4RWtU/zla9yOvn+F9v6t8c7YDO8f2T+ve3YPA88B/G/tBKj90t5danwYGu1zPN2jdHvgprXuFvz2VWoBP0RpsGgZ+6xjV9cdlvzuAezj8Be76UtcuYPmx+F4Dl9C6vbMDeKp8fazX5+sodfX6fC0Gniz7fwa4ofJ/4Lvl2O8GTi7tP1fmh8vyd01Ub5freqicr2eAP+GtdwrN2M99ZbtLeSsAZvR85U9BREQ01Ik+BhAREUeQAIiIaKgEQEREQyUAIiIaKgEQEdFQCYCIiIZKAERENNT/B0sq/EaGebaJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(trace['step_size_bar'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `get_sampler_stats` method provides more control over which values should be returned, and it also works if the name of the statistic is the same as the name of one of the variables. We can use the `chains` option, to control values from which chain should be returned, or we can set `combine=False` to get the values for the individual chains:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fbc8630a690>]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD4CAYAAADFAawfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO2deZwdVZ3ov6c73Ul39qRDAkmaJmyKUVkCAoIKigKOg844DvNGHGccM/rm+fSNo5OZYZ76EFRUxo2RQUEBHUCEESVsARKSQEjohix09qVDlk53p/d9Pe+Pu+Qutd5bVffc3N/38+lP160653d+Z6lf/ep3TlUprTWCIAiCuZQVWgFBEATBGTHUgiAIhiOGWhAEwXDEUAuCIBiOGGpBEATDmRSG0JqaGl1XVxeGaEEQhJOShoaG41rreVbHQjHUdXV11NfXhyFaEAThpEQpddDumGvoQyk1RSm1SSm1RSnVqJT6RrDqCYIgCE548aiHgau11n1KqQpgvVLqKa31KyHrJgiCIODBo9Yx+uI/K+J/oT3O+Kl7N3H32n38j5+9wr3rD/DktmY+cMeLTEykF/ntp3bypYde50fP7+Gvf7GJq7+/hqffOEb34Cjvuu05vvXUDt7//TWMZ+QbHhunbsVK6las5B9+szmr/K6BES6+NZY/s9y71+5L5l23p40Lb1nFoY4BPvLj9dyz/kDyWNPx/mSehoOd1K1YydcefwOAv72vnh8/vyd5/Ofr9vPJn28E4L82vsmf/MdLHDjez0W3rOJI1yDffWYnX3jwdYZGx3n3t1/gz+56meX319M7NErdipVc98N1SVkPbGjiE3dt8NzW33l6J1986PWs/Tf/bhs3/24bf/6fGzj35qeoW7GSz94fC2XddM9Gfr5uv6W8va19vP3rz3DuzU+xelcrAP/vD9v5yiNb0tINjY5zxXdeYN2etuS+lp4hLrplFXtb+5L7jvcNs+ybq9jR3JPc90+/3co3/tDouY6Z/NldL/PAK7E7zK89/gb//Ng2ugdGueTW5zj35qd4YENTMu1Hfrye3zYctpTzkR+v5+8eqOeq761hdHzCscyvPLKFLz70OnUrVvKOrz9jmaZ/eIzLv/U8r+xvB+CZxmNc/b01jLnIXn5/PXUrVnKoYyBt/yf+cwP3p9TFidaeoeTYbekZSu7/m1++yn+s2csLO1uoW7GSi25Zxbu//QKPbz7iSS7AFx6M1fue9QfS9n/9943JMm97cgc3/GQ9AH/605f59cZY/yTG2l0v7uPTv9jkWtbAyBiX3PocdStW8sCGpqT81t6hZJlf/W1sLP5s7X5uumdjMu8Nd76UTP9qUwef/1UD/75qd3JfKsNj47zn9tV866kd1K1Yyc2/28Yd8bTX/mCt57bxg6cYtVKqHGgAzgLu1FpvtEizHFgOUFtbm7NCa3e3sXZ37AR+eV87VRXlDI6OMzg6ztTJJ9S968V9WXlXPLaVOz7xTlp6hvnPF2PGpG94jJlVFck0hzsHk9uPvXaEOz5xfpqMV/Z30NZ7Iv/I+ARTysoBuO3Jncl0N90TGziPbz7CtiPdbDvSnTz2cP0h/unatwDw76t2A3DfhoN844alPLejhed2tPCF958NwDdX7kjm+5f/3gbAg5vepL1/hD9sOcqdq2P1/Nx7l3Cka5AjXTH9Nx/qAkgzYv/2uD8D9tM1Mdk/vPGCtP2/euXNrLSrtrcAsG7PcdbtOc7fXrkkK839G5roHRoD4GuPN3LVV0/h3pdiJ+h3/+ydyXRN7f0c7hzklie28+z/eS8AT25rpr1/hAc2NPGNG5YC8MLOVo73jXDP+gN8L57/4fpDMfkfeZuvuiZ4tamTV5s6uenS07lvQ8wgXHXuPFp7h4FYG950WR0A245084+PbOHjFy3KkpPa5539I5wyY4ptmY+kGPueePtksqO5h6PdQ9z+9E4e+5/v5p8f20ZH/wjdg6PMnTbZVvaz8X753etHkmMKYNOBDjYd6OBT8bo48UxcBsQuEIk8L+xs5YWdrZw2M1a39v4RAL708GZuOH+hq1yAP2w5CsAtT2znM1eckdz/y5ebktt3rz1x4W842EnDwU7+8l2nJ8eaV3Y096T1Y4LVO1v584trk2Xe/vF3cuuTO9LybomfTwDfe2YXGw908NQbxyzLOdw5yJsdA0kbkXq+HO8b8ayvHzwtz9Naj2utzwcWAZcopZZapLlba71Ma71s3jzLiUtBEISTGqXCketrHbXWugtYDVwbjjr25NoAYTVcFKi0bfMr4lfD1DqZX7toUDkO2DDHeaZO5vZVMJrl05ZhtY2XVR/zlFKz4ttVwDXATudcgiAIpUdYF0wvMepTgfviceoy4Dda6yfCUUcQBKF4CevO19VQa623Ahe4pTtZkdd1C4LgFSNi1IWkUAZTu6xEFEN+ciIf1BByoWAxalMo1Gnjdr7K6Xxy4tTvYsSFqCkeQ53jyVHM55RO2za/In41TK2T+bWLBhPHeaZO5vZVMJrl05a5rtpxo2gMdaFw67NivhAI9jh1q/S5EDVFY6hzPjfyPKncPJwwPV1ZR104nPo9TDst66jzwYB11DKZmGO+PE8r8ahLE2ePWjo9DE6Gdi15Q52rvc2372UysTRxnEyMTo2SGl8ngZ0O7c63aAx1rp5x3n0vLnVJ4jTeouzyUhpeJ0NVS96jzjn0kedId11HnZd0wVScPeroVkEUw2qfoChU6CPIcmUddcT5BEGIFlOflTCB4jHUBVpfWshOlHXUZpI5JoIcI1nj3KPsk2Eddf51CCY8Kuuo8yB3jzrcVR8TIZwhmV+zEaLHT7eGeREtpZEQxrnkBS9OoFdHseRDH4XCdR11CGOrlE7OKPFzV+ZnMjFcbzY82aZhcl096yaTiblmzLPccMVby4wLlQdegsWXl1ygycTMW2ev3vrJ8MBL/ncmuWmWWapVWxbYTheRoS7Q8jzXddSheNR5hmtMdk0KiJ9W8fMIeZjtXUpdWaiIn5c29hz6KPUYdcEeeHFdnhf86PI2cPLLX4r4Cn34eIQ81NBHLnmKdAAULEbtoZXFo/ZIoSYTi5HSq3FxEbYhLVI7jZ4oULkBtpc88BKQR+1bjqHL80x8F4Xpy/OCC33kG5ryLttrWanJgm7LyJbn5R+oDEYPCzFp7VuA06t4DHVgMWp/clwnE8OIUXsQ6nSbWKQOVegEN5mYu1yr/LnqYZ+nOEeA0THqDLfJDnnXR84edaaHEmy5YZwUniRKjNo3vi72PtrXrxMReugjVOnhYXSM2qNHXfKhj6DI5/bc8ngoHnXwMk2iYN+/9GWnvbvUso46GAplqP1SCC2LxlDnPJmY50mV72tOcwnZJPI4XZwdH8gokE/l1Zk4sU7c5A8HOBzLcx21U/og1lEHbe8i+3BA3nrnuI46o1zLddSePeqSD30ENFEQcYw6txij+75iXp5XqAuJP4/auxy/YzOoWLltniINfhQsRu0pjbbczkSW5wW06sNUinUCSPBP2Ia0WIdSaugjyvMh0NeclmqM2k/FFdmzrlm3qS59klleLp2YKsJrdr+eeSE8aj9t4dRvXm8dVcb/fLEzkFa3q74eePGrh0WGTBUSP6MwV8pmO2rS1lUEWHG3lRiePGqZTHQm0SjeltC4G+YgDSIW5SX0cDpuLSd720mO4615SKe3n3i93yVPyTpri32etPNQno0gK6PsHPrIbyWRtcx02Vbt4SV/LvpYjb102fnF5L2S+tbI3MqwzuX8gi3t6ZwP6yLiFaMMtbMXE0zXBW3EgopRpxkoK8ORsctxHXVYHnVQcgyORSbTOk4m5iM5uFi5fZ7ijH2kX2yiqYOOeXce0nmNUZfAZGIYt/OFXkfttTi3K3bWBceXIQmGQk3oBoUprzn1U/+gJqOLgbQYdQ75c5t49dYfXj3qkgh9ON/OBy/TW/48DXEOMWrrVR/eI6NheSOl5FE7L6POLxTg5UKcu/Ro4tphkD6ZmH08KKcoU6anMJ1H4SWx6iOK253gPWqX43msozaNoLqnULULLOSQp0dtKdI2fp6/7GLB6A8beZ0DKIV11H4mcDzLNLnzU7AyzmmrR3zJKgymfzjAFJzP89xaIv2Bl2BHQGQPvLisVQ7lmQWy+yMfW1siHnX2vkSjeekDq+V5mTndvNWs5XkeyrXSI5k/h9CHl+MFWZ7nozWcl+c5TMaEuDzPTn3r5XnexQTxro/Mffksz/Obx2p5XiHW9buG/1zzW6dwmuDT2ltdPT/wUgoxaisKvzzPLS5mcdLZbHslWWdbqf5uzYMiqOV5Ovk/e/Io1OV5dsu3LJfnhTmZaL8ve3meR5leb81d9LEqN7Llea7nZm4xarfXLWRdeC3j487HE5SGRx3CdGK295Nf/kzyHVwn0jnL9LU8z9AYd4KCTSb6iVE7XmicnYFcZNvGqHPpyyBi5vmL8I3rZKJL/lxXyHh1Ar3qEQZmGeoQbuf9esh+y/XnZXozrm5eulu5oYU+ApNbGEvtp1RHlyFzTPmtj6URsvP2/YnOSR/LcqPvo/TleVbhIef8uT9r4SH04fKcQwJ5KVPO+fK9bXO9jjsf9XjLlHYsb486HILy1AvnUXsvOKhvJnoOq9h51Dl6iflSiC5yj1H78Ir8lBmgR12w0IdSarFSarVSartSqlEp9cWQdInIo84vv+/jNtt+0/l5cCcsbyio5VOFCsz48qgDal+vYY5c4qu2ZfrOYSGjAJ3k9sCLn3PNK1Yxast0XmPUIVnqSR7SjAFf1lq/ppSaDjQopVZprbcHrYzpsVWTKERbyRv+hDA5GYZXWI+QK78nn1LqceAnWutVdmmWLVum6+vrfStz7s1PMTxm/SniSWWKsbhLV1lexsi4dbryMsV4ius3qUxRlnKZm9A6KSchK5WxiYk0zzE1v12ZViTkpuapKFeMjmvL46nHrEitP6TX00qWlwGTSJ/ZBnb1tNLfKV9qP6WmT+0Dq3byss+qfDc0Ok3/hKzMMVNZXpaV1k4OZI8xp7Rg3T+ZbZLQzUk2WLeRk+5uMhJ5rPTOTOMFK/2syoT0MW51PriN68xzO1Pf1LGTeq4AjnXNLNupnEuXzOGh5Zc5yrJDKdWgtV5mdcyLR50qqA64ANhocWw5sBygtrbWt5IAf/eeJWxq6uC0mVUc7hpk8exqqivL2XyoiyvOruGna/YB8MlLT6elZ4i+4TEWzJhCx8AIzd2DXFg7m6mTJ7FqewuXnDGHbYe7ueLsmqxyHtr0JgDnL57FW06dkXX82cZjXHLGXBqPdvPus07kb+keYldLL41He/jce8/kmcZjXLt0ARv2tfPWU6fzalMnI2MTXP/2U5O3QGPjE/zipSY+esFC5k2fzNbDXdRMm8xps6piMnuGONI5yIWnz6ajb4S9bX1cXDcnKfvNjgF6h8Z422kzWL2zlXMXTGdwZJwzT5nGw68e4t1n1bBodkzW8d5hmtr7WVY3x1N7v9k+QN/wGOedlt4Ge1r6UAq6B0bRaHoGx6iZXsk7Fs3itYOdLJxdxfwZU7LkaR1ru8kV5bzrjDlUVZZzoK2fkfEJzl0wPS3t6p2tXHbmXKZUlCf3PfPGMT74tgVpt4/PNB7jmvPmJ43V7mO9TCpXLJk3zVMdM3n1QAdL5k1l7rTJ7G3tQ2vN2fOn82zjMWZUVXDWvGnUTJ8MkOzXWdWVWXJe2d9O7Zxq9rb28Z5z5jmW+WbHABVlihd2tnLO/OlcfIZ1/7ywo5Urzq6hclIZ/cNjNBzsdJV9oK2fpxuP8fn3nZm2v76pg7q5U5N1cePpN44BcO3SBcl9Ww51ccr0ycyqruTlfcfRGqory1m6cCYzqio8yW082sOell4uP7OGU2ac0GVfax/Pbm8B4MaLF7OrpZdLl8xl04EOzoz3T2KsaQ3Huoc4v3aWa3kv7GhFo7lsyVxe3N1GmVJ8KF6n5q5BxiY0i+dU09I9xNHuQS6onQ3ApgMdNHcNopTiuqUL2NXSy6zqSg6299M7NJbWLgBrdrVxQe0s/mvjm1y3dAGTysvoGhjh7686y1O7+MWzR62Umga8CNyqtX7MKW2uHrUbdStWArDlax9kpseBIgiCmSTO56Zvf7jAmpiBk0ft6f5FKVUBPAr82s1IR0FYAXtBEAQT8bLqQwH3ADu01neEr5I7YqcFQSglvHjU7wZuAq5WSm2O/10fsl6CIAhCHNfJRK31egxzYsN6+kcQBMFEiubJxFTETAuCUEoUp6EWSy0IQglRnIZafGpBEEqI4jTUYqcFQSghitJQC4IglBJFaajFoxYEoZQoSkMtCIJQShSloZbJREEQSoniNNRipwVBKCGK01AXWgFBEIQIKU5DLS61IAglRHEa6kIrIAiCECHFaajFUguCUEIUpaEWBEEoJYrSUEuMWhCEUqIoDbUgCEIpIYZaEATBcMRQC4IgGI4YakEQBMMRQy0IgmA4YqgFQRAMRwy1IAiC4YihFgRBMBwx1IIgCIYjhloQBMFwxFALgiAYjhhqQRAEwxFDLQiCYDhiqAVBEAxHDLUgCILhiKEWBEEwHDHUgiAIhiOGWhAEwXDEUAuCIBiOGGpBEATDEUMtCIJgOK6GWil1r1KqVSn1RhQKCYIgCOl48ah/CVwbsh6CIAiCDa6GWmu9FuiIQBdXKidJpEYQhNJjUlCClFLLgeUAtbW1QYlNY80/vo+jXYOhyBYEIVqe+MIVVFeWF1qNoiAwQ621vhu4G2DZsmU6KLmpnDaritNmVYUhWhCEiFm6cGahVSgaJJYgCIJgOGKoBUEQDEdp7RylUEo9CLwPqAFagK9pre9xydMGHMxRpxrgeI55w0T08ofo5Q/Ryx8no16na63nWR1wNdRRo5Sq11ovK7QemYhe/hC9/CF6+aPU9JLQhyAIguGIoRYEQTAcEw313YVWwAbRyx+ilz9EL3+UlF7GxagFQRCEdEz0qAVBEIQUxFALgiAYjhhqQRAEwxFDLQiCYDhiqAVBEAxHDLUgCILhiKEWBEEwHDHUgiAIhiOGWhAEwXDEUAuCIBiOGGpBEATDEUMtCIJgOGKoBUEQDEcMtSAIguFMCkNoTU2NrqurC0O0IAjCSUlDQ8Nxu28muhpqpdQUYC0wOZ7+t1rrrznlqauro76+PhddBUEQShKllO0Hwb141MPA1VrrPqVUBbBeKfWU1vqVwDQUBEEQbHE11Dr2CZi++M+K+F8on4VZvauVxbOrOOuU6fQMjdI3NMapM6ewbs9x3rFoJm29w5w9f3oy/b62PvqGxpg2ZRJLaqayu6WPcxfEju9u6eWMmqm82THAmfOmAXC4c4D9bf28feFMBkfHOW1WVVLW0a5BBkbGGRmbYHB0nMVzqphUVkZH/wh1c6sZGZ+gvW+ExXOqs/Sub+rgtFlV9A2PcU6KfnYc6hig8Wg3l59Vw4wpFcn9e1p6OXPeNMrKFIc6Bpg7rZLBkXHa+0cAqJ1TTUvPEKfPnQpAe98wu471cvlZNQA0Hu1mSkU5S2qmsqe1j3PmT6ezf4TR8QmOdA1SM20yi+dUMzAyxnM7WvngefM52jXI4jnVVJSfmK5oONjJgplTOHi8n7a+Yd57zjxae4cZHp1g2pRJzJ8xOa0t2vuG2dPaF5dfxfDYBH1DY+w81sPs6koWzqqiZ2iMmVUVTK4oS/brntY+tCbZZ6ltMDA6zrHuIXqHRlk0u5p50ydntePmQ128ZcF0lIK1u49zxVk11B/sYN70ycysqqBcKSrKy1AKRsYmOGXGFA53DtDeN8KhzgGufdsCJsXrfbC9n1lVlXQOjLBodhVPbG1m8ZxqDncOsGh2FQtmVrEwZby09gxROamMWdWVyfzbj/Ywd9pkaqZVsiQ+5gD6h8fo6B+hurKcCU1WXbTW7GntY2RsgoPtA1x5Tvq42N3Sy9mnTGN0XHO4c4DJFeW09gzxzkWzKCtTdPaP8MS2Zi6sncVbFsygZ3CU0fHYOJ4/YwpjE5qugREWzU4fu797/QhL5k2lvExRM20y82dMyWrjva29zKqOjcMdzbH6jY1PcEbNVIZGJ+gcGOH0udVsO9JN7Zxq5s+YwpSKcrYe7mJmVUVyrNrROzRK49Ee3rloFlWV5ZZphuJjYcHMKbT0DHHK9CnsPNbD4jnV1Ew70ZZNx/vpHhylurKcwdFxZldXUl1Zztx4mqHR8bTzJ0Fz9yDHuoeY0Jq3njqDo12DVJaX0zkQG+P9w2NMaE2ZUsydVkl15aRknyXO99aeISrKy5g9tdKxvrniKUatlCoHGoCzgDu11hst0iwHlgPU1tbmpMznf9XApy6r41+ufyvX/WAdR7oGufnDb+WbK3ck0zzxhStYunAmAO///ovJ/V/50Ll895ld3PvpZUxMwN/eX8+cqZV09I/w0oqrWTiriiu+szqtvKZvfzi5ffm3X7DVa/l7lrDlUBcbD3Sk5QHY29rHx+/akPz94Gcv5bIz5zrW88rbT+iRkLftcDcf+cl6/unat/D5953JlbevZtnpsxkdn2DL4W4AZlVX0DUwyu5vXkflpDIu+uZzADz3D++ldk41H/7RegBu//g7+Opvt/LLv76YT//i1aw6v+u25+kdGuPGixfz0KuH+Mt31XLrx94OwPo9x/nkPVndm8YldXPY1HSiLRJ6AHz0/NPYdKCDo91DlnkXzqrK6te7PnkR1y5dwPajPVz/o3V8+ZpzeLrxGI1He7LaKcHe1l4+eudLfPbKM9iwv503jvRgh1KgdUxG6hj46Pmn8YMbL0BrzXu/uya5/5Yb3sa/Pd6YJmPu1Eoa/u2aE21w2/OUKdj/rQ8zNj6Rlj9T37/42StsjfehVV1+vu4Atz55YoxPKlPsve16AF5/s5OP/cfL3Pzht7K7pZff1B9OpvvyNefwhfefzQW3rEru+/urzuTO1fuSv685bz4Hjvezt7UvrdwntzXzpYc32+qc4AN3rOXtC2eyo7mHsQl3/+x9587jTy5cxP9+8HUADnzrepRStumv/9E6DnUMcuXZNTzwmXdZpvnyI1tYubWZK86qYf3e41x+5lxe3teepvPw2Djv+94ay/yJNJ/7VQNrdrVl6XTZt+zP/UwurpvNI5+7nPtebuLrf9jOI5+7jIvr5nDJbc+nlRU0nlZ9aK3HtdbnA4uAS5RSSy3S3K21Xqa1XjZvnmU83BXFicY70jUIkDRSCY7G92fy8r7jAOxr7WdPa+wGoCPuiXbG/+fKq00dbDzQYXmspSfdIDW19+dUxuHOAQC2HOpK7qs/2JlW/66BUQAmMr5z2do7xNjERPL39riB299mrUvv0BgAL+5uA+CV/e3JY/uP91nmSWVTk3VbALy0r93WSIN1v+5t7QVO9O3mQ11pRtqKnngdXm3qdDTSEDPSVjy3o9Xy+O6W7DZotxhDCbs17vLd0a0ZYziTzYe70n6nGsQ3O+Lj4nB30jgl8x1Kzwexu6FU1u5uY29rdn12Hut11CmVbUe6PRlpgDW72th17ER/uH2S9VBHrM/X7Tlum2ZdfJyu3xtLk9kOAKPj7vqt2dXmSScnXm2Kte/WI7E+Pdg+kLswH/hanqe17gJWA9eGoUzM85GP7QrRIiMuPExsWxN1csPVUCul5imlZsW3q4BrgJ1hKSR2WogacQ7Cw8S2NVEnN7zEqE8F7ovHqcuA32itnwhDGUXuVztT2r5QeuRSriltlkmh1dI+Nci7HUOscKHbMioiN74RF+dl1cdW4IIIdEEpZazxEE5eZMiFh4lta6JObhj1CHnMo86tGR0mliMlCj2sLmam1D8Ioq5KZnuqqDUIsbhCD4uoHC+nlSWZBKJTxA1rlKFGmXs7bhJWF7OcQh9F6VsET2Y7RB76OIkxcYyZqJMbRhnqQl/9vWDCRISVCrloZUBVjCCzHfy2SzGe+FER1Rjzc14W47g3y1AXwf27CZ1spUIuF5C858BMaIwQ8Furk7QZioqTvQsMM9Tmn/xRaOfWBlbHC+FRF7Krgiw7f49asCM6jzqctKZglKEG8we9qRcSQ9UqCUwdEyXFSd4FRhlqhfkGx029IPR3k2F5OKdy81M2rK7yJje40rNjzD4nE/NXIDTyER3EBSiq+L2fcgLRKWI7ZZahVirnRjTdwPshiotBEHIK6UmaFPo4WYnC6SgEJurkhlmGmtwbMap5SDf9gtDD1QBarvoowGRinvnt8NKEQZad5U/7jVHnq4yh66iDaOOobKKvGHUQBZbyOmqlzA81RXEr5+pRB7WOOk8Lc/JMJua3jtr4QZsjgYQ+IhokfkopxjkFoww1mP8IeRThAleHOqh11DnkSc9fyNBHkDHqTNl+8xs+aK3wMg6DKCYAGZ7K8bOOOkQ9wsIoQ10Ey6jzxst4cjvxA1tHXczL84KUlRmjzjN/MeBF5WKKUfvzqENTIzTMMtSA6de7vI1bMGpEJldwR9q+8BSj8fWDUYYazG9wV2/XNWwRROgj9xh1al5TY3VheXu29c33gRdTr94OooO4s8tLgYDxpWsxxXTiGGWoVR4vZTLxCSjL/EHoYLnP4+pjbb2dkx4FDX0EF+rJeimTz4qZeblzxkv7BRL6iM5S+0hafD1mlqHG+zrqQnmDeU/AefFkcphMJM0A2wvQNtu5UNABn4tHbbe/FGPUReLYeC7HT9oi7C+zDLUPj7oYG9sr7pOJZlS+2CYT7S5gZrRmtEQ2mZi/iMAxUSc3zDLUeG/EzHTRPfDirKGbHoHccrosz3PKHmSMOqwB7+mBl0A96nxDH3m2RAEeeIkqRh3ZOmo/oY8gdCrtB178fKXhJA595HDca+w50NBHIR8hDzRG7fzbgzIWZZnnt6VdpCOLUUeDv3d9FB9GGWowP6Rhqn65TCZOmFqZIsPtwimEz8ne3kYZ6tgj5OktnuljJ7zurNBHPKVS2eGHfMMiadm1wzEP2A2ohI5e3smdeVhlPNGZaEOreqe2byJP6p2Mn/pYaek1f2q6RPl++umE7j7y2KzuyPt91C6hqCBQeBvXTt97TNMzgDs7O1J1CCTK4KGTCzWZGFUExChDDWS1eNZtqd3JxYn9Qb8NLX18a9tjXsqy83wT+bT2EvrI1JekwSMAAA7eSURBVCF9T6osu3LS06XeEnsnHwOVrm96n3qa6HLQwTaPzdhy61N3XbJz+Ap9eDSaXsa1U138hr1yDd+kOQMRxbn9PUIeXEwnKkfeKEPt56VMmY09MRG8Ppbl5mv0vZyUbsbeJSY64VHHvEMfBbzdDOaFQYmN9P1+2yUKjzoI/E4kB1KHiBrCV5eZ2DkumGWoUZ5PQDuP2j59ML2T9wRcAIncY6L2AoKMURf0pUy55LG7C8tTeLHEqK3uuhzTB+d4GoWJOrlhlqEO6TWnJp40TkR10fHqedvrEYgaORaeS5bMGHX6/5Od9HkMLxmCLdMUTNTJDbMMNd4b0U8cWqOdj/uJb+W7jtpTvC2840F6wWGNd28fDvBfuvcYtd/Qh0WM2o+MiNZRW00kOxFIfDki/9XflEAAOsk66nBw6pqA531CzW8rN8cLnJA/bpO2QviY8rRuWJhlqHE3ZCeW51nPblstz4utBHGI23rQK1WW3TEveFqe51Nm5jtSnJauWcnOeXmeRWVyudRGtzzPekfQq4TCIOjleWE+8FKQ5Xl+nK0A+7c0l+dZrCG2+53V2CknndUtrrNH7RITTtu2vkCckOUoylaRtOV5PtdRZ4Z2UpcqZue1uE3PdXmex32WeS1WIPhanpdDfNnuk1uZIgJZ9ZFvQ1okyXd5nlveHNSyyZfbeLKVF/AKlUDstI+xGgRmGWq8Gz4/J5fW1jHqpIHwrl7ePePJk8lBRuoep0nCIAdWPp5JIVYG2o2tTEPgd5LVuj/Mc8v9tnlRfTMxwHkmEzHKUCvIOpvsHxDx5jkkjjl9EDbSGLUXT8YljZun7Dxx6l6+V/IxRvkasiA+Pabt9kftUUdEengsWC/VVkZE7eDLozawb9wwy1Cr7PdRe/WoHWPQFuGQVBm+XugSgSeY0zcTveYP1KXOI2sBPOosJ8BGCf8etbd9hcbv8rxiMmjFpGsumGWo8RaHs9rv37Rlx0a9EIkn6OpRZ1dep1tqB9HBjeh8JOVtqHPIbxe/zZTlP0ZtdbdmnuVIGyKehmEQoY+8RXgtyXtK87rGFbMMtctLhDIOeCbMjglDtnuM2l/6sMgrRl2AhY52TkDW6wiK8ET2gt/XnAYxsIyM1RuokxtmGeqMt8BBcIbQKfSRr5xU3D8cEBbeJBeDNxHxswRZBPHNRF8SCvDhgJONXMd1znc+pf3Ai/3t6Yk0Nuuo9QkZWeuoLdKn5nE1vg76+H+KzaaMtNec+pORGYNPbIa+jtpCWlGto86QlStRTCYGso46dTtEh7qYXnOar37GrKNWSi1WSq1WSm1XSjUqpb4YpkJuM/D2rzk9kd/NkJ3Io9P+2+rkqE92Wc6y7OLlJ/67G//Mi1nma07tY++BrqPO4y4l/WRJ19fTRJeDDvZl2rR9xm//L6vK83bNo9HMdx11VJOJBVlH7Wcc2Gz7wsdYDYJJHtKMAV/WWr+mlJoONCilVmmtt4ehUPZJ4y2d86oPa9OX0/I8mwuEdwH+y3DVIeNCFNk66jzyFuKFUPYedUaM2ucrc70YTCPw6UkW1TcTfU0masttk3E11FrrZqA5vt2rlNoBLAQCN9RKKY52DfJow+Hkvu1He9LSvLyvnd6hMfpHxtL2H+4cBKD+YAfjGSfai7vb2NvWl1Xe714/QkV5GcNjzmfm7mO9ye2n3zjGnKmVyd+7WnrT0m460M7kSfY3Kj1Do2m/E3XdcrgLgMbmbp7c1uyoz7PbW9h6uDv5e/3e48ysqkj+3nks1mb1Bzuy8j6x5WjWvn1t/Uk9Nh3IzmPHU9uamV1dmbavrXfYU96EjgANBzt5tOEwjfG+3tHck5U+dUzAifZqau/3rO/Krdnt+mjDYVozdD7UOWCZP1OHxL5jPUNZ+/+w5SjVldanV6acxqPdtmkSfdh4pDs5xhNsb+7JkrXrWPp4HEkZ27/fciSp0+tvdrrq1TU4mpXGC6+/2ZXcfqaxhVOmdzmkti8/Qc/QmOX+1DzN3YOuaax0Gs/BUD/acJg34n22cX97Wvjjia1H+aN3nOZbphvKzxVFKVUHrAWWaq17Mo4tB5YD1NbWXnTw4EHfytx0z0bW7TnuO58gCIIJ1EybTP3NH8gpr1KqQWu9zPKYV0OtlJoGvAjcqrV+zCntsmXLdH19vW9Fh0bHae2JeTfJ7/5lvHAobaIiY79GJ48ntjP3TSovYyzuctvJSk2Tmi5VVip2+tnhVB+rssrLFWPjE7b1GR9Pj7Wnpsv+BmVs35SKcoZHJyzrZNcW5WUqGRbI1CNTfoLyMsW4RZwjNZ1TG1i1U1J2uaJncJTrfrgOgLVfuSo5wWYVZ04tM7UuqeUm/k+pKGNgZJxpUyYxNDpuMWGbrrtGp9U1Vd8fPL+bx147wicvreWzVy5xHEOJ9vbTJon+HBodtz1frPq5clIZI2MTlv2QoLxcMT6uUQomlStGxiayzo9UMvV0Ox8S76kpc5kwzOwfu7GTwOk8txvviXPJKYzi1i/lZYpFs6sd62Ir28FQe4lRo5SqAB4Ffu1mpPNhSkU5tXNzq6RQmqSGfEwdO7VzYnrNmTqZ0+dOLbA2QjHiZdWHAu4Bdmit7whfJUHwTqmsExZKGy/rqN8N3ARcrZTaHP+7PmS9BMETIX5rQhCMwcuqj/WI4yIYipc5AUEodox6MlEQ/FJMHnURqSoYhhhqQYiI4ni0QjARMdRCUVNMHrUg5IoYakEQBMMRQy0UNTKZKJQCYqiFokZCH0IpIIZaKGrETgulgBhqoajx8lJ5UygeTQXTEEMtFDXFZPxkeZ6QK2KohaKmiBxqQcgZMdRCUVNMoQ9ByBUx1IIgCIYjhloQBMFwxFALgiAYjhhqQRAEwxFDLQiCYDhiqAVBEAxHDLUgCILhiKEWBEEwHDHUgiAIhiOGWhAEwXDEUAuCIBiOGGpBEATDEUMtCIJgOGKoBUEQDEcMtSAIguGIoRYEQTAcMdSCIAiGI4ZaEATBcMRQC4IgGI4YakEQBMMRQy0IgmA4YqgFQRAMRwy1IAiC4YihFgRBMBwx1IIgCIbjaqiVUvcqpVqVUm9EoZAgCIKQjheP+pfAtSHrIQiCINjgaqi11muBjgh0EQRBECwILEatlFqulKpXStW3tbUFJVYQPPGxCxYWWgVbrjr3FACufsspBdZEKFaU1to9kVJ1wBNa66VehC5btkzX19fnp5kgCEIJoZRq0Fovszomqz4EQRAMRwy1IAiC4XhZnvcgsAE4Vyl1WCn1mfDVEgRBEBJ4ilH7FqpUG3Awx+w1wPEA1QkK0csfopc/RC9/nIx6na61nmd1IBRDnQ9KqXq7gHohEb38IXr5Q/TyR6npJTFqQRAEwxFDLQiCYDgmGuq7C62ADaKXP0Qvf4he/igpvYyLUQuCIAjpmOhRC4IgCCmIoRYEQTAcYwy1UupapdQupdRepdSKiMterJRarZTarpRqVEp9Mb7/60qpI0qpzfG/61Py/HNc111KqQ+FqFuTUmpbvPz6+L45SqlVSqk98f+z4/uVUupHcb22KqUuDEmnc1PaZLNSqkcp9aVCtZfVO9NzaSOl1F/F0+9RSv1VSHp9Vym1M172fyulZsX31ymlBlPa7q6UPBfFx8DeuO4qBL18913Q56yNXg+n6NSklNoc3x9JeznYhmjHl9a64H9AObAPWAJUAluA8yIs/1Tgwvj2dGA3cB7wdeAfLdKfF9dxMnBGXPfykHRrAmoy9t0OrIhvrwC+E9++HngKUMClwMaI+u4YcHqh2gt4D3Ah8EaubQTMAfbH/8+Ob88OQa8PApPi299J0asuNV2GnE1xXVVc9+tC0MtX34VxzlrplXH8+8D/jbK9HGxDpOPLFI/6EmCv1nq/1noEeAi4IarCtdbNWuvX4tu9wA7A6b2ZNwAPaa2HtdYHgL3E6hAVNwD3xbfvAz6asv9+HeMVYJZS6tSQdXk/sE9r7fQkaqjtpa3fme63jT4ErNJad2itO4FV5PnBDCu9tNbPaq3H4j9fARY5yYjrNkNr/YqOnfH3p9QlML0csOu7wM9ZJ73iXvEngAedZATdXg62IdLxZYqhXggcSvl9GGdDGRoq9krXC4CN8V3/K34Lc2/i9oZo9dXAs0qpBqXU8vi++Vrr5vj2MWB+AfRKcCPpJ0+h2yuB3zYqhI5/Q8z7SnCGUup1pdSLSqkr4/sWxnWJQi8/fRd1e10JtGit96Tsi7S9MmxDpOPLFENtBEqpacCjwJe01j3AT4EzgfOBZmK3XlFzhdb6QuA64O+VUu9JPRj3GgqyxlIpVQn8MfBIfJcJ7ZVFIdvIDqXUvwJjwK/ju5qBWq31BcA/AP+llJoRoUpG9l0Kf0G6QxBpe1nYhiRRjC9TDPURYHHK70XxfZGhlKog1hG/1lo/BqC1btFaj2utJ4CfceJ2PTJ9tdZH4v9bgf+O69CSCGnE/7dGrVec64DXtNYtcR0L3l4p+G2jyHRUSn0a+CPgL+MnOfHQQnt8u4FY/PecuA6p4ZFQ9Mqh76Jsr0nAnwAPp+gbWXtZ2QYiHl+mGOpXgbOVUmfEvbQbgd9HVXg8/nUPsENrfUfK/tT47seAxGz074EblVKTlVJnAGcTm8AIWq+pSqnpiW1iE1FvxMtPzBr/FfB4il6fis88Xwp0p9yehUGal1Po9srAbxs9A3xQKTU7ftv/wfi+QFFKXQt8FfhjrfVAyv55Sqny+PYSYm20P65bj1Lq0vg4/VRKXYLUy2/fRXnOfgDYqbVOhjSiai8720DU4yvX2dCg/4jNlu4mdmX814jLvoLYrctWYHP873rgAWBbfP/vgVNT8vxrXNdd5DkL76DXEmKz6VuAxkS7AHOB54E9wHPAnPh+BdwZ12sbsCzENpsKtAMzU/YVpL2IXSyagVFisb/P5NJGxGLGe+N/fx2SXnuJxSoT4+yueNo/jffxZuA14CMpcpYRM5z7gJ8Qf6I4YL18913Q56yVXvH9vwQ+l5E2kvbC3jZEOr7kEXJBEATDMSX0IQiCINgghloQBMFwxFALgiAYjhhqQRAEwxFDLQiCYDhiqAVBEAxHDLUgCILh/H9WU56w6T86pAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sizes1, sizes2 = trace.get_sampler_stats('depth', combine=False)\n",
    "fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True, sharey=True)\n",
    "ax1.plot(sizes1)\n",
    "ax2.plot(sizes2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7fbc862fdd50>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAOi0lEQVR4nO3df6zdd13H8eeLlYHKjw56aZa2szOU6IIBlpttBKPAhIxq1iXCwi8pS2MTMgwKUaf+gb/+gBiZLCFoZYSOOGCiuAanuIyRRWMrdw7GfohcJttaBy1jq5IFdPD2j/OpudR7e87tPT/aT5+P5OR8v5/v55zv+35y+7rffr7f8z2pKiRJfXnKrAuQJI2f4S5JHTLcJalDhrskdchwl6QOrZt1AQAbNmyorVu3zroMSTqt3Hnnnd+sqrnltp0S4b5161YWFhZmXYYknVaSPLjSNqdlJKlDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ6fEJ1QlqQc3Hnho1a9548XnTaCSEY/ck3wtyZeSfCHJQmt7TpJbk3ylPZ/T2pPkuiSLSe5OcuFEKpckrWg10zKvqKoXV9V8W78GuK2qtgG3tXWA1wDb2mM38MFxFStJGs1a5tx3AHvb8l7giiXtN9TAfmB9knPXsB9J0iqNGu4F/H2SO5Psbm0bq+qRtvx1YGNb3gQ8vOS1B1ubJGlKRj2h+lNVdSjJ84Bbk/zr0o1VVUlqNTtufyR2A5x33mROKEjSmWqkI/eqOtSeDwOfAi4CvnFsuqU9H27dDwFblrx8c2s7/j33VNV8Vc3PzS17r3lJ0kkaGu5JfiTJM48tA68G7gH2ATtbt53AzW15H/CWdtXMJcDRJdM3kqQpGGVaZiPwqSTH+t9YVX+X5PPATUl2AQ8CV7b+twDbgUXgCeCqsVctSTqhoeFeVQ8AL1qm/VHg0mXaC7h6LNVJkk6Ktx+QpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdWjkcE9yVpK7kny6rZ+f5ECSxSSfSHJ2a39aW19s27dOpnRJ0kpWc+T+DuD+JevvBa6tqucDjwG7Wvsu4LHWfm3rJ0maopHCPclm4OeAD7X1AK8EPtm67AWuaMs72jpt+6WtvyRpSkY9cv9j4NeB77f15wKPV9WTbf0gsKktbwIeBmjbj7b+PyDJ7iQLSRaOHDlykuVLkpYzNNyT/DxwuKruHOeOq2pPVc1X1fzc3Nw431qSznjrRujzMuDyJNuBpwPPAt4PrE+yrh2dbwYOtf6HgC3AwSTrgGcDj469cknSioYeuVfVb1bV5qraCrwe+GxVvQm4HXht67YTuLkt72vrtO2fraoaa9WSpBNay3XuvwG8M8kigzn161v79cBzW/s7gWvWVqIkabVGmZb5P1X1OeBzbfkB4KJl+nwHeN0YapMknSQ/oSpJHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUoeGhnuSpyf55yRfTHJvkt9t7ecnOZBkMcknkpzd2p/W1hfb9q2T/REkSccb5cj9u8Arq+pFwIuBy5JcArwXuLaqng88Buxq/XcBj7X2a1s/SdIUDQ33Gvh2W31qexTwSuCTrX0vcEVb3tHWadsvTZKxVSxJGmqkOfckZyX5AnAYuBX4KvB4VT3ZuhwENrXlTcDDAG37UeC5y7zn7iQLSRaOHDmytp9CkvQDRgr3qvpeVb0Y2AxcBPz4WndcVXuqar6q5ufm5tb6dpKkJVZ1tUxVPQ7cDrwUWJ9kXdu0GTjUlg8BWwDa9mcDj46lWknSSEa5WmYuyfq2/EPAq4D7GYT8a1u3ncDNbXlfW6dt/2xV1TiLliSd2LrhXTgX2JvkLAZ/DG6qqk8nuQ/4eJI/AO4Crm/9rwc+mmQR+Bbw+gnULUk6gaHhXlV3Ay9Zpv0BBvPvx7d/B3jdWKqTJJ0UP6EqSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDo3yHqiSdcW488NCsS1gTj9wlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6NDTck2xJcnuS+5Lcm+Qdrf05SW5N8pX2fE5rT5LrkiwmuTvJhZP+ISRJP2iUI/cngXdV1QXAJcDVSS4ArgFuq6ptwG1tHeA1wLb22A18cOxVS5JOaGi4V9UjVfUvbfm/gPuBTcAOYG/rthe4oi3vAG6ogf3A+iTnjr1ySdKKVjXnnmQr8BLgALCxqh5pm74ObGzLm4CHl7zsYGs7/r12J1lIsnDkyJFVli1JOpGRwz3JM4C/BH6lqv5z6baqKqBWs+Oq2lNV81U1Pzc3t5qXSpKGGCnckzyVQbD/eVX9VWv+xrHplvZ8uLUfArYsefnm1iZJmpKh36GaJMD1wP1V9b4lm/YBO4H3tOebl7S/PcnHgYuBo0umbyRpJk7370RdrVG+IPtlwC8CX0ryhdb2WwxC/aYku4AHgSvbtluA7cAi8ARw1VgrliQNNTTcq+ofgKyw+dJl+hdw9RrrkiStgZ9QlaQOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalD62ZdgCSt1o0HHpp1Cac8w13Sik4mRN948XkT34eGM9wljZVhfWpwzl2SOmS4S1KHDHdJ6pDhLkkdMtwlqUNDwz3Jh5McTnLPkrbnJLk1yVfa8zmtPUmuS7KY5O4kF06yeEnS8kY5cv8IcNlxbdcAt1XVNuC2tg7wGmBbe+wGPjieMiVJqzE03KvqDuBbxzXvAPa25b3AFUvab6iB/cD6JOeOq1hJ0mhOds59Y1U90pa/Dmxsy5uAh5f0O9ja/p8ku5MsJFk4cuTISZYhSVrOmk+oVlUBdRKv21NV81U1Pzc3t9YyJElLnGy4f+PYdEt7PtzaDwFblvTb3NokSVN0suG+D9jZlncCNy9pf0u7auYS4OiS6RtJ0pQMvXFYko8BLwc2JDkIvBt4D3BTkl3Ag8CVrfstwHZgEXgCuGoCNUuShhga7lX1hhU2XbpM3wKuXmtRksbPuzWeWfyEqiR1yHCXpA75ZR3SacppFp2IR+6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ14KKZ0CvKxR4+aRuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQl0JKJ+AlijpdeeQuSR0y3CWpQ4a7JHXIcJekDnlCVWcUT5DqTGG465RwMqH7xovPm0AlUh8Md522PAqXVma4ayIMXmm2PKEqSR3yyP00t9ojZOeppTOD4a6hnGKRTj+G+xnGoJbODM65S1KHPHKfEI+QJc2SR+6S1KGJHLknuQx4P3AW8KGqes8k9gPT+2SjR+KSTidjD/ckZwEfAF4FHAQ+n2RfVd037n2dLINaUu8mMS1zEbBYVQ9U1X8DHwd2TGA/kqQVTGJaZhPw8JL1g8DFx3dKshvY3Va/neTLE6jlVLIB+OasizgFOA6OwTGOA/CmtY3Dj660YWZXy1TVHmDPrPY/bUkWqmp+1nXMmuPgGBzjOAxMahwmMS1zCNiyZH1za5MkTckkwv3zwLYk5yc5G3g9sG8C+5EkrWDs0zJV9WSStwOfYXAp5Ier6t5x7+c0dMZMQQ3hODgGxzgOAxMZh1TVJN5XkjRDfkJVkjpkuEtShwz3MUtyWZIvJ1lMcs0y29+Z5L4kdye5LcmK16meroaNwZJ+v5CkknR5Odwo45Dkyvb7cG+SG6dd4zSM8G/ivCS3J7mr/bvYPos6JynJh5McTnLPCtuT5Lo2RncnuXDNO60qH2N6MDiB/FXgx4CzgS8CFxzX5xXAD7fltwGfmHXd0x6D1u+ZwB3AfmB+1nXP6HdhG3AXcE5bf96s657ROOwB3taWLwC+Nuu6JzAOPw1cCNyzwvbtwN8CAS4BDqx1nx65j9fQWy9U1e1V9URb3c/gcwA9GfX2E78PvBf4zjSLm6JRxuGXgA9U1WMAVXV4yjVOwyjjUMCz2vKzgf+YYn1TUVV3AN86QZcdwA01sB9Yn+TctezTcB+v5W69sOkE/Xcx+Gvdk6Fj0P7LuaWq/maahU3ZKL8LLwBekOQfk+xvd1PtzSjj8DvAm5McBG4Bfnk6pZ1SVpsdQ/llHTOS5M3APPAzs65lmpI8BXgf8NYZl3IqWMdgaublDP4Hd0eSn6yqx2da1fS9AfhIVf1RkpcCH03ywqr6/qwLO5155D5eI916IcnPAr8NXF5V351SbdMybAyeCbwQ+FySrzGYX9zX4UnVUX4XDgL7qup/qurfgX9jEPY9GWUcdgE3AVTVPwFPZ3AzrTPJ2G/bYriP19BbLyR5CfCnDIK9xznWE45BVR2tqg1VtbWqtjI473B5VS3MptyJGeU2HH/N4KidJBsYTNM8MM0ip2CUcXgIuBQgyU8wCPcjU61y9vYBb2lXzVwCHK2qR9byhk7LjFGtcOuFJL8HLFTVPuAPgWcAf5EE4KGqunxmRY/ZiGPQvRHH4TPAq5PcB3wP+LWqenR2VY/fiOPwLuDPkvwqg5Orb612CUkvknyMwR/yDe3cwruBpwJU1Z8wONewHVgEngCuWvM+OxtDSRJOy0hSlwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1KH/BSAIv/HIuBBEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "accept = trace.get_sampler_stats('mean_tree_accept', burn=1000)\n",
    "sb.distplot(accept, kde=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8062428977787616"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "accept.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find the index of all diverging transitions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([], dtype=int64),)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trace['diverging'].nonzero()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is often useful to compare the overall distribution of the\n",
    "energy levels with the change of energy between successive samples.\n",
    "Ideally, they should be very similar:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fbc86b0f450>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXxb1Zn4/88jyZL3fUvsJM4KZCcrEAiUHYZCGWgJtAxd+AGdMl+mnf469MvvRSnTTqfTzgylpQttKQxrgBZIS2jYyhaS4CQYQvbESRzb8b7Jm2RJ5/fHlYPjLUpi+0ry8369/JJ079HVI0V+fHLuuc8RYwxKKaXil8PuAJRSSo0uTfRKKRXnNNErpVSc00SvlFJxThO9UkrFOZfdAfSXm5trSkpK7A5DKaViypYtWxqMMXmD7Yu6RF9SUsLmzZvtDkMppWKKiBwaal9EQzcicrmI7BaRfSJy9yD7V4rIVhEJiMj1g+xPF5FKEfnFiYWulFLqVB030YuIE3gIuAKYDdwoIrP7NasAvgw8NcRh/g145+TDVEopdbIi6dEvA/YZY8qNMX7gGeCavg2MMQeNMR8Dof5PFpHFQAHw6gjEq5RS6gRFMkZfBBzu87gSWB7JwUXEAfwX8CXg4mHa3QbcBjB58uRIDq2UikI9PT1UVlbS3d1tdyhxKzExkeLiYhISEiJ+zmifjP1HYK0xplJEhmxkjHkYeBhgyZIlWnxHqRhVWVlJWloaJSUlDPc7r06OMYbGxkYqKyuZOnVqxM+LJNFXAZP6PC4Ob4vE2cB5IvKPQCrgFpF2Y8yAE7pKqdjX3d2tSX4UiQg5OTnU19ef0PMiSfSlwEwRmYqV4FcBN0VycGPMF/sE+GVgiSZ5peKbJvnRdTKf73FPxhpjAsCdwDpgJ/CsMWa7iNwvIleHX3ipiFQCnwd+IyLbTzgSpZRSoyKiMXpjzFpgbb9t9/a5X4o1pDPcMR4FHj3hCJVSMeupTRUjeryblutkjZMRdVfGKmWH/glJE4rqFQgEcLliO1VqUTOlVFx54oknWLZsGQsXLuT2228nGAySmprKPffcw4IFCzjrrLOora0FoL6+nuuuu46lS5eydOlS1q9fD8B9993HzTffzIoVK7j55pupr6/nkksuYc6cOdx6661MmTKFhoYG7r33Xh544IGjr33PPffws5/9zJb3PRxN9EqpuLFz505Wr17N+vXrKSsrw+l08uSTT9LR0cFZZ53FRx99xMqVK/ntb38LwF133cU3v/lNSktL+eMf/8itt9569Fg7duzg9ddf5+mnn+b73/8+F154Idu3b+f666+nosL6H+BXv/pV/vd//xeAUCjEM888w5e+9KWxf+PHEdv/H1FKqT7eeOMNtmzZwtKlSwHo6uoiPz8ft9vNVVddBcDixYt57bXXAHj99dfZsWPH0ee3tbXR3t4OwNVXX01SUhIA7733Hi+88AIAl19+OVlZWYBVhDEnJ4cPP/yQ2tpazjzzTHJycsbmzZ4ATfRKqbhhjOGWW27hRz/60THbf/rTnx6dluh0OgkEAoDVC9+4cSOJiYkDjpWSkhLRa9566608+uij1NTU8NWvfvUU38Ho0KEbpVTcuOiii3j++eepq6sDoKmpiUOHhqzey6WXXsrPf/7zo4/LysoGbbdixQqeffZZAF599VWam5uP7rv22mv561//SmlpKZdddtlIvI0Rpz16pSI02FRBnZ0zvLH+fGbPns0PfvADLr30UkKhEAkJCTz00ENDtn/wwQf5xje+wfz58wkEAqxcuZJf//rXA9p973vf48Ybb+Txxx/n7LPPprCwkLS0NADcbjef+cxnyMzMxOl0jtp7OxViTHSVllmyZInRhUfUWItkeqUm+uPbuXMnZ5xxht1hjDifz4fT6cTlcrFhwwa+/vWvH+39h0IhFi1axHPPPcfMmTPHJJ7BPmcR2WKMWTJYe+3RK6XUcVRUVPCFL3yBUCiE2+0+Omtnx44dXHXVVVx77bVjluRPhiZ6pZQ6jpkzZ/Lhhx8O2D579mzKy8ttiOjE6MlYpZSKc5rolVIqzmmiV0qpOKeJXiml4pyejFVKjZ7NfxjZ4y35ysgezyb33XcfqampfPvb3+bee+9l5cqVXHzxxbz77rvccccdJCQksGHDBu69917Wrl3LlVdeyU9+8pOTfj1N9EopdRyjWar4/vvvP3r/ySef5Lvf/e7RwmgPP/wwTU1Np3whliZ6pfpp6+7h52/s5d29DTR1+pmel8KV8yZgjEFEmF7x3KeNndnWbZz0NOPBE088wYMPPojf72f58uX88pe/xOl0kpqayl133cVf/vIXkpKSeOmllygoKKC+vp477rjjaEXKBx54gBUrVnDfffexf/9+ysvLmTx5Mg8++CA33XQT1dXVnH322bz22mts2bKFBx98kOzsbP75n/8ZsEoV5+fnc9dddx0T1w9/+EMee+wx8vPzmTRpEosXLwbgy1/+MldddRUtLS08++yzrFu3jldeeQWv10t7ezuLFy/mu9/9LjfccMNJfyY6Rq9UHx9XtvCz1/fyX6/toasnyNTcFLZXt3HXM2U8/UEFXf6g3SGqYQxVphiwtVTxli1beOaZZygrK2Pt2rWUlpYOiP3WW2/l6quv5ic/+QlPPvkka9asISkpibKyslNK8qA9eqWOemt3Ha/uqGVSVhKPfnUZ0/NSAQiGDA+/U85P1+2mq+cQP5gEDl3/OioNVaYYsLVU8bvvvsu1115LcnLy0eOOJU30SgGlB5p4dUctCydlct2i4qNJHsDpEL5+wXT21nr504dVrHFn87nCJhujVUMZqkwxQEJCwrgtVaxDN2rc++BAEy+WVbEwvZ3/N28jsyqft2aL9JsxsnhKFvOKMlhdlUd5p8emaNVwTrRMMYxNqeKVK1fy4osv0tXVhdfr5c9//vMJv7dTEVGPXkQuB34GOIHfGWP+o9/+lcADwHxglTHm+fD2hcCvgHQgCPzQGLN65MJX6tR0+AJ8+7mPyExO4JvTqnANMyQjInxuYRHlNU386Ugu355eNXaBxqoxPkk9VJniKVOmDPmcsShVvGjRIm644QYWLFhAfn7+0aGlsXLcMsUi4gT2AJcAlUApcKMxZkefNiVYyfzbwJo+iX4WYIwxe0VkIrAFOMMY0zLU62mZYjWW/r8Xt/HkpgpuPXcaF3euPbp9+dR+s2lCQf781vsEHQmUbt3KU3WT+O85B7n2jNRj241z8VqmGKKrVPFolCleBuwzxpSHD/YMcA1wNNEbYw6G94X6PtEYs6fP/WoRqQPygCETvVJjpexwC09srOBr505lam4KDCw3Dw174YPfwo4X+Wx7LQCfA/7Rk0P5oVmQuwDyTh/TuJU9YrlUcSSJvgg43OdxJbD8RF9IRJYBbmD/iT5XqZFmjOHf/rKD3FQP37xkFmvKqgc2aquC310EPd0w61I2OpcAwoT6d2loqGd+z0ewaQPknQEl50JudP6Sq5ERy6WKx2TWjYhMAB4HbjHGhAbZfxtwG8Dkybpijxp9f/n4CFsONfOf180n1TPIr4G3Bjb+ChIz4PZ3IKuE8vAKU2ICVKa4WbVjMr/IeZ6LGtchD53N20t+wYV/t2qM30n06b2wTI2Ok1kVMJJEXwVM6vO4OLwtIiKSDrwM3GOM2ThYG2PMw8DDYI3RR3pspU6GLxDkx3/dxewJ6Vy3uHjAfkfQB5t+AyLwD2sgq2RAm+IkPxOSAvx399XkzZjBaYee5Pwtd4J757Htx9nYfWJiIo2NjeTk5GiyHwXGGBobGwedCjqcSBJ9KTBTRKZiJfhVwE2RHFxE3MALwP/2nqBVym7ffu5jKpu7+PI5BawuPTxgf2HTJuhugRX/DLkzhjzO2VleVlfnccRkwpQvMvvgH3B+8DCsuAtSC0bzLUSt4uJiKisrqa+vtzuUuJWYmEhx8cAOynCOm+iNMQERuRNYhzW98hFjzHYRuR/YbIxZIyJLsRJ6FvBZEfm+MWYO8AVgJZAjIl8OH/LLxpjBJ6oqNcq6/EHe2lVHSU4yM/NTB+x3BruZ0LARCuYO2pPv65ysNlZX57GhOZ3PFgTYNeVLLDz0CHy8Gs7+J+t/BONMQkICU6dOtTsM1U9EY/TGmLXA2n7b7u1zvxRrSKf/854AnjjFGJUaMY9vPIjXF2DVssmDDi0UNm7AFepmW+o5dB5oYn9wsKk44baJPUxL7mJDUxqfLWjC586C066Cbauh+kMoWjSab0WpiOmVsWrc6PIHefidcmbkp1rTKftxBTqZ0LiJxvQz6EwsjOiYZ2V52d+ZRKM/3GeavBzSi2HnSxDwjWT4Sp00rXWjxo1nSitoaPdz7ZmDj2/mtXyIM+SnKu/8o9uOKUk8iEUZ7TxVlU9ZWwoX5baCOGDu38P7D8L+N+CsO0b0PSh1MrRHr8YFXyDIb94uZ9nU7EF78wC5LdvwJhXTlZgf8XGLE/3kJPRQ1tpnvD97GkxYCAfeAV/7qYau1CnTRK/GhT9uqaKmrZt/unDwWTSZbbtJ9tXRmDHvhI4rAgszOtjWlkyg78TgqSsh0A2f6GQzZT8dulFxLxAM8au397GgOINzZ+Ty9AcDp1SWVL9MCAeNGbNP+PgL09t5oyGTPe1JrOjdmDUV0ibC2/8JodCxM3DG2dx6ZT/t0au4t+ajag43dXHnhTMHv4jHhJhy5BVaU6cTcEVWf7yvuemdODGUtfUZvhGBkhVWGYWW4cvkKjXaNNGruBYMGR762z5OL0zjotMHH3vPb9pCSncNjZknNmzTK9kZ4rTULj5q7fdHomgxuDxwaP1JHVepkaKJXsW113bUsL++g298ZgaOIdb/KznyMj3OZJrTTjvp15mf3sHBrkQafX1ew5UIRUutOfX+jpM+tlKnShO9imuPrD9IcVYSV86bMHgDE6Ko9i2q81cSciSc9OvMTbcS+ft17mN3TFoGoQDU7RjkWUqNDT0Zq+LW9upWPjjQxBVzCwetaQOQ07qdJH8jVXnn4wx2nfRrTUvuJskR5P06N5+d1OdCqYxiqwJmzTYoHttVhZTqpYlexZWnNn1asuCPWypJcApLpmQP2X5i3duEcFCddy6Tal476dd1CsxJ62R9XfKxO8Rh1c2pLIWgH5zuwQ+g1CjSoRsVlzp8AT6qbGHR5CyS3APX8OxVVP8ODVkL8bszT/k156Z3UtHh5HBHv1+rwnlWkm/Ye8qvodTJ0B69iksfVbYQCBmWT80ZsK+3rEFCTxvZbTupyL/ouKUOIjEvzRqnX1/nZtXU7k935MywTszWbIOCOaf8OkqdKO3Rq7i0taKZoswkCjOGXqAh02v1sFvSZo3IaxYl+slPDLK+/wlZhwvyz4DaT2DgAmtKjTpN9Cru1LR2U93SzZmThx+OyWrfS3dCJl2e3BF5XRFYkd/D+3VuQv3XSSucB/52aD44Iq+l1InQRK/iztaKZpwiLCgeOtFLqIf09nJa0maO6AIh5+T7afQ52N3a77xA3mzrxKxOs1Q20ESv4krIGMoOt3BaYRopgy36HZbeWYHTBGhJHXqpwJOxIt8PMHD4JiERMidD474RfT2lIqGJXsWVQ42dtPsCLJg0/LBNRvt+QuLEm1Iyoq8/MTnEtNTAwEQP1knZlgotXazGnCZ6FVd2VLfidAizBlkPtq+M9v14kyef0tWwQ1lR4GdTfQI9/c+75sywTsYe3jTir6nUcDTRq7hhjGHHkTZm5KXiSRh67nxCTxvJvnpaU6ePShwr8nvoDDooa+r3RyRrqjVOf/C9UXldpYai8+hVTOt7JeyR1i6aO3u4YNbwK0RltJcD0DJKif6sPD+CYX1dAktzez7d4fJY4/Sa6NUY0x69ihs7jrQhwOkT0oZtl9m+H78rlS5P5EsGnohMt2FeVmBggTOwhm+qt+o4vRpTESV6EblcRHaLyD4RuXuQ/StFZKuIBETk+n77bhGRveGfW0YqcKX623XEy6TsZNIShx53FxMkvaPcGrYZwWmV/Z2T72drYwIdgX6vkTPDqmap4/RqDB030YuIE3gIuAKYDdwoIv3XW6sAvgw81e+52cD3gOXAMuB7IpJ16mErdaxOX4Dqli5mFQx/Eja7dQcJwS5aU6aNajwr8v0EjPBB/SDj9A6XDt+oMRVJj34ZsM8YU26M8QPPANf0bWCMOWiM+RjoP8/gMuA1Y0yTMaYZeA24fATiVuoY5Q0dGGB63vCJfkLDegzQmjq6iX5pbg9uhxk4fOPywMRFmujVmIok0RcBfYt5V4a3RSKi54rIbSKyWUQ219fXR3hopT61v74dt8tBcVbysO0KG96nI3HiSa0NeyISnbA4p4f36gYZRpq8HI58BAHfwH1KjYKoOBlrjHnYGLPEGLMkLy/P7nBUDNpf387UnBScQywXCJDQ00pu80ejNq2yvxX5fna2Jhy7vCBYC5AEfVDzyZjEoVQkib4KmNTncXF4WyRO5blKRaS1q4eGdj/T84bvpRc2bMJBaNSmVfbXWw5hQ//hm96VpipLxyQOpSJJ9KXATBGZKiJuYBWwJsLjrwMuFZGs8EnYS8PblBox++utqYrTjjs+/z5+VxrtycVjERbzsgKkuUIDyyGkT4T0Ik30aswcN9EbYwLAnVgJeifwrDFmu4jcLyJXA4jIUhGpBD4P/EZEtoef2wT8G9Yfi1Lg/vA2pUZMeX0HyW7nsLXnMYYJDe9Rk3OWdXXqGHA5YHlez8BEv/kPkJIL5X+z7m/+w5jEo8aviK6MNcasBdb223Zvn/ulWMMygz33EeCRU4hRqWEdauxgSk4KjmHmxWe07yelu5ZPZpwDpn+x+JGz6cCx/Zhil+H1jgL+vMtLvse6Snb51GzInGKdkPV5wTP8BV5KnaqoOBmr1Mnq8AVo7PAzOXv42TYTGtYDcCR3xViEdVTv8oLbvP3iyyqxbpsPjWk8anzSRK9i2uGmToDjJ/r69bSkTqczacJYhHVUUaKfrIQePmnrd6I4o9gaQmo5OKbxqPFJE72KaRVNnTgEijKThmzjDHaR37yFmtxzxjAyiwjMTevkE2/yscsLOt3WCVnt0asxoIlexbSKpk4mZCThdg39Vc5v2owz5Kc699wxjOxTc9M6aQu4qOjyHLsjc4q1EIkuGK5GmSZ6FbMCwRCVzV1MOs6wzcT69QQcidRlLx6jyI41P90ap/+4//BN1hTrwilvjQ1RqfFEE72KWbtrvfiDISZnDz1sA1DYsJ667MWEnJ5h242WbHeA4kQfH/VP9JmTrduWioFPUmoEaaJXMevDihYAJmcPfUVsSmclGR0HOZJnz7BNrwXpHexqT6I72GcKaEoeuBI10atRp4lexayPK1tIdjvJSh66/vyEhvcBqB7jaZX9LchoJ2Ac7GzvM8wkDqtXr4lejTJN9CpmfVLVRlFmEjLMhVITGtbTnjQRb0rJ2AU2iDNSu0iQ0ODDN95q6OmyJzA1LmiiVzHJFwiyp9bLxGGmVTpCPRQ2bLIukhrF1aQi4XYYZqd1Dp7oTQhqttkTmBoXNNGrmLS7xksgZIZN9HnNW0kIdlBt8/h8r/npHVR3e6jq7PNrlxE+IVu11Z6g1LgQUa0bpaLNJ1VtwPAXShXVvU1InCR3HmF6xXNjFdqQFqR38Djwbq2bVVO7rY1JmeDJgKottsam4pv26FVM+qS6lfRE17AnYifWvUNbSgkhp3vINmOpONFPdkIP79T0iydzkiZ6Nao00auY9ElVK3OLMoY8EZvWcZD0zkO0pM4a48iGJmL16t+rcxPoezFs5hRo2g9dzbbFpuKbJnoVc3qCIXYd8TK3KGPINhPr3gGgOW3mWIUVkQXpHbT1OPiouc+oae+FU9Uf2hOUinua6FXM2Vvbjj8YYs7E9CHbFNW/Q0vqDPzuzDGM7PjmpncgGN6t7TN8kxFebVNPyKpRoolexZwdR6wTsXMmHtujn17xHNMrnmPWgcfJbyylMzHfjvCGleYKMT8rwDs1fcoxuJMhZ4YmejVqNNGrmLOn1ovb5aAkZ/BiZhnt5TgI0RxF4/N9nV/op6zJRau/z/mFiYv0hKwaNZroVcx4alMFT22q4K3ddeSkuHl2c+Wg7TK9e+hxJo3ZIuAnamWBnxBy7FqyRYuhvQbaqu0LTMUtTfQq5tS2+ShIH2IhcBMis30frakzxmwR8BO1MLuHNFeId/qO0xctsm61V69GQXT+Jig1hO6eIK1dPRSkDV5yOLWrmoRgJy1RNtumL5cDzsn3826t+9N1ygvngcOl4/RqVESU6EXkchHZLSL7ROTuQfZ7RGR1eP8mESkJb08QkcdEZJuI7BSR745s+Gq8qW2zrigdqkef6d2DQWhJmT6WYZ2wlYV+qjqd7Pc6rQ0JSZA/W3v0alQcN9GLiBN4CLgCmA3cKCKz+zX7GtBsjJkB/A/w4/D2zwMeY8w8YDFwe+8fAaVORl2bD4D8oRJ9+168yZMIuoZfjMRuKwv8AP2GbxZDdRmEdGlBNbIi6dEvA/YZY8qNMX7gGeCafm2uAR4L338euEisSxYNkCIiLiAJ8ANtIxK5GpdqvN24nQ4yByl94O5pJaW7lpa06Jxt09eklBDTUgMDE72v1bpKVqkRFEmiLwIO93lcGd42aBtjTABoBXKwkn4HcASoAH5qjGnq/wIicpuIbBaRzfX19Sf8JtT4UdvWTX66B8cgpQ8yvXuB6LsadigrC/1srHfjC4Y36AlZNUpG+2TsMiAITASmAv8iItP6NzLGPGyMWWKMWZKXlzfKIalYVtfmoyBtqPH5vXQnZNLtzh3jqE7Oefl+uoPC5obw/07yToeEFD0hq0ZcJIm+CpjU53FxeNugbcLDNBlAI3AT8FdjTI8xpg5YDyw51aDV+NTuC9DuC1CQPnDGjSPoI73jgDXbxuZFRiJ1Vn4PCWI+Hb5xOGHCAqjabG9gKu5EkuhLgZkiMlVE3MAqYE2/NmuAW8L3rwfeNMYYrOGaCwFEJAU4C9g1EoGr8adumBk3BU2bcZoALakzxjqsk5biMizJ7eHtvuP0k5bCkY91aUE1oo6b6MNj7ncC64CdwLPGmO0icr+IXB1u9nsgR0T2Ad8CeqdgPgSkish2rD8YfzDGfDzSb0KND7XeoWfcTKh/j5C4bF8bNhKbDjQd/ZnqbmFXa8LRP2JMWg6hHmv2jVIjJKIVpowxa4G1/bbd2+d+N9ZUyv7Pax9su1Ino7atm8QEB+mJA7+2E+vfpS1lCiHH0AuRRKMF6R08VQXv7m3gusXFVqIHOLwJppxtb3AqbuiVsSpm1LV1U5CWOGCxkdSOw+FFRmJjtk1fk5N8ZLgCvLM3PNssJReyp8PhD+wNTMUVXTNWxQRjDLVtPuYNstjIxIZ3AWJqfL6XQ6xFw9/dUUmotAyHAMnZcOBtKH0Eln7V7hBVHNAevYoJdV4fXT3BQWfcTKx/l7bkKfg82TZEduoWpHfQ5HewvSXc78qaBv526GiwNzAVNzTRq5iwp9YLDDwR6wx2k99YypG8c+0Ia0TMT+8A+pRDyC6xbpsP2BOQijua6FVM2F1jJfr+Uyvzm0pxhXxU551nR1gjIiMhyJzMHt6pCSf61AKryJkmejVCNNGrmLCn1kuKx0Wq59jTShPr3yPgSKQ2O7avw1tZ4GdLYwLtPWLV0c8sgeaDdoel4oQmehUT9tS2H61B37s27PSK55hy5BW8yZOYWtX/Gr7Ycl6Bn4ARNtSHp4dmTwVvDXS12BuYigua6FXUM8awt9Y7YNjG42sk0d9srSYV44y3lkRHkGf3GjYdaILsaYCBio12h6bigCZ6FfWqWrro8AcHJPrM9n0AtKTFfqJ3OWBuWidlrSnWqlOZU6wVpw6+a3doKg5ooldRr3fGTf+plZnt++hy5+BzZ9kR1ohbmNFBvd9Ntc8NzgTIKtFEr0aEJnoV9XbXtAOQ36c8sSPUQ3rHwbjozfc6M8N6n2WtKdaGnBlWgTMdp1enSBO9inp7a70UpieS5HYe3ZbecQCHCcbk1bBDyXUHKE70Udaaam3ImYE1Tr/B1rhU7NNEr6Le7lovswrTjtmW0V5OUFx4k6fYFNXoWJjRzo72JDoDWOP0Tg8c0OEbdWo00auoFgwZ9tW1c1pB6jHbMzrK8aZMwTjiq1zTwvQOAsbBhrrwOP2kZTpOr06ZJnoV1SqaOvEFQsws+LRH7+5pI8nXQGvKVBsjGx2np3bhcYR4q7ccQsl5ULMNuprtDUzFNE30Kqr1lj44rU+iT28vB6AtdcDywzEvwWGYm9bBWzUea5plybmAgUPv2x2aimGa6FVU651aObPP0E1GRzk9zhQ6PQV2hTWqFmZ0cLjDSXm7E4qXgCsJyt+yOywVwzTRq6i2p9bLpOwkkt3hsXhjyGg/QGvq1JhZBPxELQxXs3yrxg0uj9Wr3/eGzVGpWKaJXkW1PbXeY4ZtMr17SAh20JoSf8M2vfI9PUxPC1iJfvMfwJ0KTfvhnZ9aj5U6QZroVdTyB0KU13cwq0+iL2y0ar/E4/h8XxcU+tlU76YjIJB/urWxfpe9QamYpYleRa2DjR0EQubYRN+wgS53Lv6EdBsjG30XT/DhD4lVoz4lH5KyoG6n3WGpGKWJXkWt3hk3vYneEfSR37SF1jjvzQMsze0hyx1iXbXHOheRdzo07oVQ0O7QVAyKKNGLyOUisltE9onI3YPs94jI6vD+TSJS0mfffBHZICLbRWSbiCT2f75Sg9lb68XpEKblWbVf8lrKcIW6rROxcc7lgIsm+HjjiJueEFaiD/h0MRJ1Uo6b6EXECTwEXAHMBm4Ukdn9mn0NaDbGzAD+B/hx+Lku4AngDmPMHOACoGfEoldx7fWddWQlu/nT1iqe2lRBYcNGQuKiLbnE7tDGxGVFPrw9DjbWJ0DuLGvlKR2nVychkh79MmCfMabcGOMHngGu6dfmGuCx8P3ngYtERIBLgY+NMR8BGGMajTH6f08Vkdq27mNKExc2bqAhcz4hp2eYZ8WP8wr8JDkN66o81hqyWSU6Tq9OSiSJvgg43OdxZXjboG2MMQGgFcgBZgFGRNaJyFYR+c5gLyAit4nIZhHZXF9ff6LvQYYjInoAABxDSURBVMWh7p4gTR3+o4uNuP0tZLfuoCbnLJsjGzuJTji/0Mer1R5CBsifDW2V0Fpld2gqxoz2yVgXcC7wxfDttSJyUf9GxpiHjTFLjDFL8vLyRjkkFQv21bVj4GiiL2jchGCoyT3b3sDG2GUTfdR1O/moyQUFc62Nu9faG5SKOZEk+ipgUp/HxeFtg7YJj8tnAI1Yvf93jDENxphOYC2w6FSDVvFvb114VanwguATGjfgd6XSmDHXzrDG3IUT/LjEWLNvUgsgJQ92v2J3WCrGRJLoS4GZIjJVRNzAKmBNvzZrgFvC968H3jTGGGAdME9EksN/AM4HdoxM6Cqe7a5pxylCTqoHjKGwYSO12UvjrizxUDYdaGLTgSZ2VTUyO7WTlw66MIjVqz/wDnS32R2iiiHHTfThMfc7sZL2TuBZY8x2EblfRK4ON/s9kCMi+4BvAXeHn9sM/DfWH4syYKsx5uWRfxsq3uyp9ZKX5mFW5fPM3v8wqV1VBB2JTK94zu7QxtySTC9HfB72e51Wog/1wH6tfaMiF1H3yBizFmvYpe+2e/vc7wY+P8Rzn8CaYqlUxPbUeskPz7jJCJclHg/z5wezNLOdRw7DuioPM04rgaRsa/hmzrV2h6ZihF4Zq6JOuy9AZXPX0ROxGe3l+BLS6Xbn2ByZPbLdAWamdLG2ygMOJ8y6DPasg2DA7tBUjNBEr6LO3treE7GJYEKkdxy0qlXGaVniSJyd1cb2lgRr+Oa0K6G7BQ69Z3dYKkZooldRZ29tOwAF6R5SuqrDZQ/iv77NcM7O8iIY/nzYAzMuhoRk2P6i3WGpGKGJXkWd3bVeEhMcZKW4yegILxsYh+vDnohsd4DleT2sOZyISUiyhm92/lmHb1RENNGrqLOn1svM/DQcImS0H6AjsZCAK8XusGx39aRuyr0utle3wezPQWcDHFpvd1gqBmiiV1FnT62XmQWpuAKdpHYdpnWc9+Z75QdrcGL48Su7WN16BgFnEuzQ4Rt1fJroVVRp7eyhts3HaQVp5DdtxmFC4358vleaK8SZGe2UHW7BL4lU5Z1nDd9ojXp1HJroVVTZEy59MKswjcLGDYTEhTd5ss1RRY+VOa14fQH21bVTUXgpdNTDofftDktFOU30Kqr0XVWqsGED3uTJGEeCzVFFj0UZHaQ6g+zZ9TGJvgZwJMBb/2F3WCrKaaJXUWV3jZdUj4uJ0khm+35adNjmGAkOw4rsNja3pOINJULBbKj5SIdv1LA00auosrvGy+mFaUj52wC0aaIfYGVOKz3GwfvN6TBhIfi8ULHB7rBUFNNEr6KGMYZdNW2cVpgG+9+ky51Dp6fA7rCizvTkbiYndfNGQ6a1GIkjQS+eUsPSRK+ixpHWbtq6A5xekALlb1mLjIzjsgdDEYGLc1s40JnINm8K5J8BO9fo8I0akiZ6FTV6T8Qu8lRBZ8O4W03qRJyb3YZbQjx9IMkavmmvhYqNdoelopQmehU1doUT/bS2DwA4kqOJfigprhBnZ7fxUoWH9uy54ErUi6fUkMbHcj0qJuyqaWNiRiJJFW9D/hy6E3X94OFcnNvC242ZvFCdwc25s+CjZyDvdJBw/23JV+wNUEUN7dGrqPDUpgo2lTeR6wkSPPg+O5MX2x1S1JuZ0s2CrB7+sC+J0ISF4GuDpnK7w1JRSBO9igqBUIh6r4/zPHtxmh5qcs+xO6SoJwJfmdlJudfFehZas2+qP7Q7LBWFNNGrqNDg9RM0huWhjwg63NRlL7I7pJhwZbGPvMQgvz+QBQVz4IhePKUG0kSvokJNWzcAs7s2U591JkFnks0RxQa3A740rYu3ajwcyVwM/nZo2m93WCrKaKJXUaG2rZtCaSG3cz9HdNjmhNw8vYskp+GBxiXgdOvwjRogokQvIpeLyG4R2Scidw+y3yMiq8P7N4lISb/9k0WkXUS+PTJhq3hT09rNZz1bAHD1dDC94jmmVzxnc1SxIdtjWDW1iz8eTqczd54O36gBjpvoRcQJPARcAcwGbhSR2f2afQ1oNsbMAP4H+HG//f8NvHLq4ap4VdPWzQXObfQ4k+lM1LIHJ+rWWZ0AvBQ4C3o6oWGPzRGpaBJJj34ZsM8YU26M8QPPANf0a3MN8Fj4/vPARSLWtesi8jngALB9ZEJW8aa1s4fWLj/zQzutRUa07MEJK0oO8bnJ3fywejE9Dg91ezbx1KYKntpUYXdoKgpEkuiLgMN9HleGtw3axhgTAFqBHBFJBf4V+P5wLyAit4nIZhHZXF9fH2nsKk7srvUyVw6QZry0pk63O5yY9U9ndNJlPGx2LCC7bReOUI/dIakoMdonY+8D/scY0z5cI2PMw8aYJcaYJXl5ejXkeLO7po2LnVsxCC2pM+wOJ2ZNSQ1yQW4Lv+taiSvUTWGDli5WlkhKIFQBk/o8Lg5vG6xNpYi4gAygEVgOXC8i/wlkAiER6TbG/OKUI1dxY2eNly86P6Q9qYiAK8XucGLKpgNNxzy+boKLbzfOpZ1kphz5K9X5K22KTEWTSHr0pcBMEZkqIm5gFbCmX5s1wC3h+9cDbxrLecaYEmNMCfAA8O+a5FV/9VUHmCMHaE6bZXcoMS/HHeCCvHbWBpZSVPsmjqDP7pBUFDhuog+Pud8JrAN2As8aY7aLyP0icnW42e+xxuT3Ad8CBkzBVGowxhiK698BoEUT/Yj4XGEj68xy3MEOJja8Z3c4KgpEVL3SGLMWWNtv27197ncDnz/OMe47ifhUnKts7uK8UCkNiRPo8uj5mZGQkRAkM3cCTS2p5Fe8AmgVy/FOr4xVttp9uIYVju0czD5Pp1WOoCsLW3nNLGdKw9vg77Q7HGUzTfTKVh0738AjPTQWXWR3KHEl1RWiquhykuhm//t/sjscZTNN9MpW+YdfoU3SaM5bancocSdn9oU0kkHdhqcxxtgdjrKRJnplH38n8zvWsz3jfEKOBLujiTsJbjc7Mj/Dwu4PeHvbAbvDUTbSRK9s07ZtLSl00zT1s3aHErc6Z15Nkvh5/5WnCIa0Vz9eaaJXtvGVPUu9ySB79oV2hxK3mnIW0Z2Yx+L2v/FSWf/rHNV4oYle2WPjr8iqfJOXg8uZXb9WSxKPEiNOPPP/ns84P+JX68rwBbR88XikiV7Zo+YTXKaHD9zLyXDrkMJokrnX4aaHOd71PLlRq1mOR5rolT2qtlBLNmROsTuS+Fe8FNKL+If0Lfzib/vwdmtVy/FGE70ae82HMPW7eCZwPrOzQnZHE/8cDphzLQv9Wwl0NPG7d3UGznijiV6NvS2PAvBM4EIWZGvvckzM+XscoR6+PXkvv323nKYOv90RqTEUUa0bpUZMwA8fPs7B5Pkc6c5hfpYuNDNaek9wb6oAjGFBQiaXel/ke/4z+ebqMi6bU3i07U3LJ9sUpRoL2qNXY2vXn6GjnpccFzEtNaAnYseKCE3ps8nv2s85E2FDeSOdvoDdUakxooleja3SR/AmFfFYy3yK3B1sOtA0YPEMNToaMufhIMT/k7GZnkCI9/Y12B2SGiOa6NXYqS6DQ+/xyYTraA64mZHSZXdE40pXYgHtiRNZ0PBn5k5M533t1Y8bmujV2Hn7PyExg1eTrgRgRkq3zQGNP/VZC8ny7mVVcSP+QIj39muvfjzQRK/GxpGPYffLcNY32N/mwCmGKUm6zN1Ya8yYS8DhYWnLWuZOTGfD/kY6/dqrj3ea6NXYePvH4MmA5bdzuLmLkqRuEhx6InasBZ2JHC68hJLqtVw6MwNfIMT6fY12h6VGmSZ6NfpqtsGuv8BZX6fHnU5lcyczdXzeNvuLr8Ud8LKk613mTExnQ3kDbXq1bFzTRK9GlzHwyr9CUhacdQefVLXSEzScnqaJ3i512UtoSynhtENPccHMPLp7Qjy+4ZDdYalRpIleja5tz8Gh9XDR9yApi9KD1lTKM1J1HVPbiIPdU75ITusnLJA9zCpI5ZH3DtDl18qW8SqiRC8il4vIbhHZJyJ3D7LfIyKrw/s3iUhJePslIrJFRLaFb7Xw+HjS3Qrr7oGixbDoFgA+ONBEToqbzARNKnY6UHQ1flcapx18nAtm5dPY4efpD7SyZbw6bgkEEXECDwGXAJVAqYisMcbs6NPsa0CzMWaGiKwCfgzcADQAnzXGVIvIXGAdUDTSb0JFqb/9CDrq4cwvwdbHCBko3ZfLovQ2uyMb13pLIzRkzGNyzWtckXY6ZVOv5OF3yvniWZPxuJw2R6hGWiQ9+mXAPmNMuTHGDzwDXNOvzTXAY+H7zwMXiYgYYz40xlSHt28HkkTEMxKBqyhXsw0++A1MORsyrToqe9uctPY4OCNNh22iQW2OtSB7YdMH3PmZGdS0dfOnrboKVTyKpKhZEXC4z+NKYPlQbYwxARFpBXKwevS9rgO2GmMGTJ4WkduA2wAmT9biSjFr8x+sWxOC938OCUlw2t8d3f1BgxuAM1L1RGw08Cdk0Jgxh/zmLRRONMwvzuBXb+3n84uLcTn19F08GZN/TRGZgzWcc/tg+40xDxtjlhhjluTl5Y1FSGo0VW6G5gNw+mfBnXJ086b6BAqTguS5dSpftKjKW4kjFGD38/czvyiDiqZO/u8L2+wOS42wSBJ9FTCpz+Pi8LZB24iIC8gAGsOPi4EXgH8wxuw/1YBVlPN3ws41kFUCk5Yd3Rw08F6dmxX5fkTsC08dq9uTS33mAmZWrGZRZif5aR7e2l1PKKQXs8WTSBJ9KTBTRKaKiBtYBazp12YNcEv4/vXAm8YYIyKZwMvA3caY9SMVtIpiu9eCvwPmXg/y6dfr4yYXLX4H5xfqghfRpipvJRjD/P2/4YLT8qnz+nh1R63dYakRdNxEb4wJAHdizZjZCTxrjNkuIveLyNXhZr8HckRkH/AtoHcK5p3ADOBeESkL/+SP+LtQ0aH1sDVnvuRcyCg+ZtfbtW4Ew3n5muijjd+dyb7JX2Ba1Yucl15Ldoqbh/62D2O0Vx8vIlphyhizFljbb9u9fe53A58f5Hk/AH5wijGqWBAKwbbnrTH5064YsPvtGg8LsgNkeTR5RKNPZtzOlOq1nL39fi6Y+QB/KqvhzV11XHRGgd2hqRGgp9bVyPjwcWg5BLOvhoTkY3Y1+4SPmlycX6C9+Wjlc2exZfa/ktv6MTeaV5iam8KP/7qLoI7VxwVN9OrUdTbB6/dB9jQoWjpg93t1bkII5xdqWeJodmjClVTlnceZ+37BveemsKe2nT9uqbQ7LDUCNNGrU/fG961yB3OvZ7ApNa9Ve8hyh1iQrXXPo5oIpXPuxeDggo+/w9LiZP7rtd106CpUMU8TvTo1lVtgy2Ow/HZInzhgd1cAXq92c3mRD6dOq4x6nUmFbFjwQ6R6K7/KepI6bzc/fXW33WGpU6SJXp28UBBe/hak5sMF3x20yVs1HjqDDq4q1mUDY0VlwUWw8jvk7n2On0/bzKPvH6TscIvdYalTENGsG6UG2PwHayrlkTI482b45I+DNvtLpYdcT4jleXo1bDTrLXR2VEkhzLqCv9vzAKXJ/4fvPJ/KS984lyS3FjyLRdqjVyfH1w67XoacGTBx0aBNOgPwxhEPVxR349JvWmwRB1z/CFJyLveFfs7Mhte558VtOrc+RmmPXp2cXX+GQDfMvW7QE7AAr1Z76A4K05z1bDqghcxizser4fSrkNZKHmz+Bd8p8/P4pEz+4ewSuyNTJ0j7WerElb8NhzfB1PMhbcKQzR7fn0RJaoDTtVplzNl0oMn6OdzB5oIv4E2Zwn+5f039yz9gTZmWMo41mujVielug5e+ASl5g14B2+uTZhdbGt18aVoXDp1tE9OCTg97Jt9EYOIS/sX1HME/3s5fPyy3Oyx1AjTRqxOz7v9CWxUsvAmc7iGbPb4/iSSn4fMlOtsmHhiHE9eZX8R33t1c43yPKS9cw7OvvmN3WCpCmuhV5HassUodrLgLsqYO2azZJ7x0OJHPTe4mw60n7+LFpoPN/DH1i7x55i8odjZxxfovsPq3/0G3Xy+oinaa6FVkqrbCn26zFvoeYs58r1/vTsYXhK/M0CUD41Fd4UpeX/k8VYkzuKHqR7z/o7/jly9v4qlNurh4tNJZN+r4Wirg6VWQmgc3PgOuoZf9PdLp4NF9yVw7pZtZGcExDFKNhb7z7VunX8PGqmLOa3mVuR98gRdTV+HjdDy9U+2XfMWeINUAmujV8Gq3w9M3Qk83/MMa6yrYYTy4MwUDfHN2x9jEp+wjDqR4KVsyp1Bc8RK3dfyal/96ARMX/x1nFiYcbda/p3/Tcl0Xeqzp0I0a2o6X4HeXQMAHN78A+acP23xro4vVBxL54rQuJqWExihIZTdJzaf69K9Qlno+l5l3KCr9Ib9/aycH6tvtDk2FaY9eDXT4A3jzB3DgbShaAjc8AelDz5cH6AgI3ypNJ8cdYGVqJZsOaKIfT4zDhW/K+fRkLMNsXs3X2n/Dlp//jddO/xcoXkaqR1ONnfTTH682/+HT+yYELYfB4bTWfK3bAcm5cNm/w9Jbhx2TBzAG7i9L5VC7k3tnVZDs1CQ/XiVlF5F0yT/TVl7KrL3rWLzndt7bNY+NE28mZ+4leBI05dhBP/XxKtgDDXuhdps1Du9rA3HClHPg8h/DmV8CT2pEh3pgRwqrDybxjdM7mJ2iV8GOe+Igffpy+NxPaPzbQ8zZ9CvOrfkOB45MYHvuFfin/RPuvGl2RzmuSLQVKVqyZInZvHmz3WHEp2AAyt+yapjseBGCfnB6rLH3gnmQf4a15muEjIFf7krmJ9tT+XxJFz9e7KX0YNPoxa9iyv7J1jLSjqCPtH0vMvHQi5wZ/ASAxuRpZE1dhGPiAvCk6QydESAiW4wxSwbbpz36eGcMVH8I256zFu/uqIPEDGs+fOE8yJkFzhP/GrT1CP+6OY1XqhK5elI3P1rk1VIHalAhp4fW026gZdYX2FH2FL6Gcs5p30DO9ucJbf8TJmcmTocLzrgKkrLsDjcuaY8+HoWCULMN9r5m9d4b91rlCmZeCvNvgFmXQdlTJ3Xo7iCsPpDEz3ak0Noj3FhUz1X5TUMVsFRqgKUl2bxa7eYvO5o4rWMz1zjfZ7LUEZIEgtM+Q8L86+H0K62evorYcD36iBK9iFwO/AxwAr8zxvxHv/0e4H+BxUAjcIMx5mB433eBrwFB4P8YY9YN91qa6CNgjLVGa0cDdNRbP+210LgfGnZD1RZrP8CUcyG9CCYsAHfySb1cVwA21rt544ibNYcTaetxcHaen3vmt9PRUjeCb0yNB8unZgPW13hrk4un9ydy+MgRLjQbucq5kSJppAcXdRnz8RYshwnzSSqaQ1rhDNweD26ngwSnINq7OMYpJXoRcQJ7gEuASqAUuNEYs6NPm38E5htj7hCRVcC1xpgbRGQ28DSwDJgIvA7MMsYMecnkiCZ6Y8I/QauX2/c2GIBQj3VSMhT49NaEet8UIBHe9mmPsS4u6umCQJd129MZfuyz9h/9zM2xsfbdZoz1vI76TxN6e92n90ODrNiUkAK5M60hmakroeQ8a1pk3xk2vYc2EAxBjxECIfAFhSa/g8ZuB01+oaLdye42F7taXRxsdxI0QoKEWJ7l5cLcVmandmovXo2YgIFd7cns87pxdRxhlu8T5pndzJMDOOXT35N2k0grKbSaVNpIwUcCIRwEcRIS69Y4nIg4MQ4X4nCSmpyMJHiQhCScCYk43cm4PEm4PMm4E5NxepIxrkRwJWFciYjDBQ4n4nCAOBGnExEH4nAiDic4HIg4EYcLI0IIByFxWLcIIZy4nE6cTlf4VnA5BKej763j6GPHCI15nuoY/TJgnzGmPHywZ4BrgB192lwD3Be+/zzwC7H+3F4DPGOM8QEHRGRf+HgbTuaNDKujAR6YZyXq3mRu4mCan9Nj9cTdadYsmMxi66SpJzW8LbzdnWrdl/A1cP4O2PPXo4f5+oZ03jjiIWggaI7/xRIMk1OCzMoIckWRj6W5PTg6anE7omuoT8UHl8DctE7mpnUCScBSWoPLWOsLYjqbcXc3UJDQiSvQiSvYiSfQSZ7fi4t2xJijKdZhrFsJhZCQwUEQlz+I2/TgwX/MH42xEDJCCCGI9Ydgv5nIVf5/H/j+w38AzpycyTO3nT3icUSS6IuAw30eVwLLh2pjjAmISCuQE96+sd9zi/q/gIjcBtwWftguIo1AQyRvwGa5xHGcB4ExLkQb15+nDTTOkTUCcTYAVw25dw+w+vaTPviUoXZExawbY8zDwMO9j0Vk81D/BYkmGufI0jhHlsY5smIlzsFEUuumCpjU53FxeNugbUTEBWRgnZSN5LlKKaVGUSSJvhSYKSJTRcQNrALW9GuzBrglfP964E1jneVdA6wSEY+ITAVmAh+MTOhKKaUicdyhm/CY+53AOqzplY8YY7aLyP3AZmPMGuD3wOPhk61NWH8MCLd7FuvEbQD4xnAzbvp4+PhNooLGObI0zpGlcY6sWIlzgKi7YEoppdTI0nr0SikV5zTRK6VUnIuaRC8iPxGRXSLysYi8ICKZffZ9V0T2ichuEbnM5jg/LyLbRSQkIkv6bC8RkS4RKQv//Doa4wzvi5rPsz8RuU9Eqvp8jlfaHVMvEbk8/JntE5G77Y5nOCJyUES2hT/DqKkpIiKPiEidiHzSZ1u2iLwmInvDt7ZXNhsizqj9bh5P1CR64DVgrjFmPtZ1A98FCJdRWAXMAS4Hfhkuy2CXT4C/Z/BrifYbYxaGf+4Y47j6GzTOKPw8B/M/fT7HtXYHA0dLgTwEXAHMBm4Mf5bR7DPhzzCa5n4/ivW96+tu4A1jzEzgjfBjuz3KwDghCr+bkYiaRG+MedUYEwg/3Ig15x76lFEwxhwAesso2MIYs9MYs9uu14/UMHFG1ecZQ46WAjHG+IHeUiDqBBhj3sGamdfXNcBj4fuPAZ8b06AGMUScMStqEn0/XwVeCd8frATDgDIKUWKqiHwoIm+LyHl2BzOEWPg87wwP4T0SDf+ND4uFz60vA7wqIlvCJUaiWYEx5kj4fg1QYGcwxxGN383jGtMSCCLyOlA4yK57jDEvhdvcgzXn/smxjK2vSOIcxBFgsjGmUUQWAy+KyBxjTFuUxWm74eIGfgX8G1ai+jfgv7D+8KsTc64xpkpE8oHXRGRXuJca1YwxRmSMK49FLma/m2Oa6I0xFw+3X0S+jFXx5yLz6QT/MS+jcLw4h3iOD/CF728Rkf3ALGDUToSdTJxEQVmKSOMWkd8CfxnlcCJl++d2IowxVeHbOhF5AWvoKVoTfa2ITDDGHBGRCUBULnJgjKntvR9l383jipqhm/DiJt8BrjbGdPbZFRNlFEQkr/ekpohMw4qz3N6oBhXVn2f4F73XtVgnlaNBJKVAooKIpIhIWu994FKi53McTN8SKrcAUfm/0Sj+bh6fMSYqfrBOCh4GysI/v+6z7x5gP7AbuMLmOK/FGp/1AbXAuvD264Dt4di3Ap+Nxjij7fMcJO7HgW3Ax1gJYILdMfWJ7UqsGWH7sYbHbI9piDinAR+Ff7ZHU6xYCxEdAXrC38+vYZU0fwPYi7U4UXaUxhm1383j/WgJBKWUinNRM3SjlFJqdGiiV0qpOKeJXiml4pwmeqWUinOa6JVSKs5poldKqTiniV4ppeLc/w9zBLNNbfyKkAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "energy = trace['energy']\n",
    "energy_diff = np.diff(energy)\n",
    "sb.distplot(energy - energy.mean(), label='energy')\n",
    "sb.distplot(energy_diff, label='energy diff')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the overall distribution of energy levels has longer tails, the efficiency of the sampler will deteriorate quickly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multiple samplers\n",
    "\n",
    "If multiple samplers are used for the same model (e.g. for continuous and discrete variables), the exported values are merged or stacked along a new axis.\n",
    "\n",
    "Note that for the `model_logp` sampler statistic, only the last column (i.e. `trace.get_sampler_stat('model_logp')[-1]`) will be the overall model logp."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = pm.Model()\n",
    "with model:\n",
    "    mu1 = pm.Bernoulli(\"mu1\", p=0.8)\n",
    "    mu2 = pm.Normal(\"mu2\", mu=0, sigma=1, shape=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "CompoundStep\n",
      ">BinaryMetropolis: [mu1]\n",
      ">Metropolis: [mu2]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "        <style>\n",
       "            /* Turns off some styling */\n",
       "            progress {\n",
       "                /* gets rid of default border in Firefox and Opera. */\n",
       "                border: none;\n",
       "                /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "                background-size: auto;\n",
       "            }\n",
       "            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "                background: #F44336;\n",
       "            }\n",
       "        </style>\n",
       "      <progress value='22000' class='' max='22000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [22000/22000 00:13<00:00 Sampling 2 chains, 0 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 2 chains for 1_000 tune and 10_000 draw iterations (2_000 + 20_000 draws total) took 13 seconds.\n",
      "The number of effective samples is smaller than 10% for some parameters.\n"
     ]
    }
   ],
   "source": [
    "with model:\n",
    "    step1 = pm.BinaryMetropolis([mu1])\n",
    "    step2 = pm.Metropolis([mu2])\n",
    "    trace = pm.sample(10000, init=None, step=[step1, step2], cores=2, tune=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'accept', 'accepted', 'p_jump', 'scaling', 'tune'}"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trace.stat_names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both samplers export `accept`, so we get one acceptance probability for each sampler:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[2.50000000e-01, 3.96972963e-05],\n",
       "       [2.50000000e-01, 1.32994734e-01],\n",
       "       [1.00000000e+00, 1.71805501e-02],\n",
       "       ...,\n",
       "       [1.00000000e+00, 4.23641715e-02],\n",
       "       [2.50000000e-01, 2.22210183e-03],\n",
       "       [2.50000000e-01, 1.80279632e-01]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trace.get_sampler_stats('accept')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "numpy   1.18.5\n",
      "seaborn 0.10.1\n",
      "pymc3   3.9.0\n",
      "pandas  1.0.4\n",
      "last updated: Mon Jun 15 2020 \n",
      "\n",
      "CPython 3.7.7\n",
      "IPython 7.15.0\n",
      "watermark 2.0.2\n"
     ]
    }
   ],
   "source": [
    "%load_ext watermark\n",
    "%watermark -n -u -v -iv -w"
   ]
  }
 ],
 "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}