docs/source/notebooks/sampler-stats.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "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": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "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": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "As a minimal example we sample from a standard normal distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "model = pm.Model()\n",
    "with model:\n",
    "    mu1 = pm.Normal(\"mu1\", mu=0, sd=1, shape=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 2000/2000 [00:01<00:00, 1361.25it/s]\n"
     ]
    }
   ],
   "source": [
    "with model:\n",
    "    step = pm.NUTS()\n",
    "    trace = pm.sample(2000, tune=1000, init=None, step=step, njobs=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "NUTS provides the following statistics:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'depth',\n",
       " 'diverging',\n",
       " 'energy',\n",
       " 'energy_error',\n",
       " 'max_energy_error',\n",
       " 'mean_tree_accept',\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": {
    "deletable": true,
    "editable": true
   },
   "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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "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": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f234bfd5588>]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYVNW97vFv9UTP0NDFPMqwFBEVREFBEdEomuPNUWOi\nGUw0ecwlOTE3k8abPJmOiRoOiZico4keb3KcohHFOEaEgOCACIiKS+apGRq66YEe6O6q+0cNVHVX\nd1V1V3XtDe/neXyeqtq79n7dVP161Vpr7+3x+/2IiIh7ZWU6gIiI9IwKuYiIy6mQi4i4nAq5iIjL\nqZCLiLhcTm/vsLKyrtvTZMrKCqmubkhlnJRwai5wbjblSo5yJedEzOX1lng6W+aqFnlOTnamI8Tk\n1Fzg3GzKlRzlSs7JliuhFrkx5h5gVnD9X1lrn4lYNhe4C2gDXrTW/iIdQUVEJLa4LXJjzMXAJGvt\nDOBy4LftVrkPuAa4ALjMGDMx5SlFRKRTiXStrACuCz4+AhQZY7IBjDGnAFXW2t3WWh/wInBJWpKK\niEhMcbtWrLVtwNHg05sJdJ+0BZ8PBiojVj8IjO1qe2VlhT3qJ/J6S7r93nRyai5wbjblSo5yJedk\nypXwrBVjzNUECvllXazW6ahqSE9Gkr3eEior67r9/nRxai5wbjblSo5yJedEzNXVH4BEBzs/BdwJ\nXG6trYlYVEGgVR4yLPiaiIj0kkQGO/sC9wJXWWurIpdZa3cApcaY0caYHOAq4NV0BBURkdgSaZFf\nD5QDfzXGhF57HdhorV0MfAN4PPj6k9baT1KeEjhc08SL7+zmkrOG0ifPmXNERUQyIZHBzgeBB7tY\nvgKYkcpQsbxrD/L061sYUV7IGacMSPfuRERcwzVndvqCN8Dw+XQjDBGRSK4p5CIiEpsKuQhQXdfM\nd3+/itXva9KVuI8KuQjwxsZ9VNc186v/tybTUUSSpkIuIuJyKuQiIi6nQi4i4nKuK+SafCgiEs01\nhdwT/3pcIt2mT5e4mWsKuYiIxKZCLoK67MTdVMhFRFxOhVxExOVUyEXQYKe4mwq5iIjLqZCLiLic\n+wq5pheIiERJ9ObLk4DngIXW2vvbLZsPfAFoA9611t6W8pQiItKpRG6+XAQsApbGWFYKfB+YZa2d\nCUw0xkxPeUoREelUIl0rzcA8INYV948F/ys2xuQAhUBV6uKJiEg8idx8uRVoNcbEWtZkjPkZsA1o\nBJ6w1n7S1fbKygrJyclOOmhxcR8A+vYtwOstSfr96ebETCFOzeakXEVFfcKPnZQrknIl52TKlVAf\neWeCXSs/AiYAtcDrxpgzrbUbOntPdXVDt/ZVX98MQE1NI5WVdd3aRrp4vSWOyxTi1GxOy3X0aHP4\nsZNyhTjteIUoV3J6kqurPwA9nbVyGrDNWnvIWnsMWAlM7eE2RXqdR2cEiYv1tJDvAE4zxhQEn58D\nbO7hNkVEJAlxu1aMMVOBBcBooMUYcy2wBNhurV1sjLkXWGaMaQVWW2tXpjOwXxPJRUSiJDLYuRaY\n3cXyB4AHUpgpJv30lXTyq30gLua+MztFRCSKCrkI+sUn7qZCLiLicirkIiIup0IuIuJyKuQiIi6n\nQi4i4nLuK+Sa7ysiEsU1hVyzw0REYnNNIRcRkdhUyEVEXE6FXETE5VTIRURcToVcRMTlVMhFRFxO\nhVxExOVcV8h1PpCISDT3FHJdMFpEJKa4t3oDMMZMAp4DFlpr72+3bATwOJAHvGetvTXlKUVEpFNx\nW+TGmCJgEbC0k1UWAAustecCbcaYkSnMJyIicSTSIm8G5gE/bL/AGJMFzAI+D2CtnZ/SdCK9xKOu\nO3GxuIXcWtsKtBpjYi32AnXAQmPMFGCltfaOrrZXVlZITk520kGLi/sAUFpagNdbkvT7082JmUKc\nms1JuYqK+oQfOylXJOVKzsmUK6E+8i54gGHA74AdwAvGmCuttS909obq6oZu7ai+vhmA2tpGKivr\nurWNdPF6SxyXKcSp2ZyWq76+KfzYSblCnHa8QpQrOT3J1dUfgJ7OWjkE7LTWbrXWthHoRz+9h9sU\nEZEk9KiQB7tdthljxgdfmgrYHqfqgl8TyUVEosTtWjHGTCUwM2U00GKMuRZYAmy31i4GbgMeCQ58\nbgSeT0dQDUVJOmmwU9wskcHOtcDsLpZvAWamMJOIiCTBPWd2iohITCrkIiIup0IuIuJyKuQiIi6n\nQi4i4nIq5CIiLufCQq4zgkREIrmnkOt8DUkjfbzEzdxTyEVEJCYVchERl1MhF0EjL+JuKuQiIi6n\nQi6CBjvF3VTIRURcznWFXDeWEBGJ5ppCrp++IiKxuaaQi4hIbCrkIiIul1AhN8ZMMsZsNcZ8s4t1\nfmWMWZ6yZCIikpC4hdwYUwQsApZ2sc5E4MIU5hIRkQQl0iJvBuYBFV2sswC4MyWJREQkKTnxVrDW\ntgKtxpiYy40xNwH/BHYkssOyskJycrITTxhUXJwPQGnfArzekqTfn25OzBTi1GxOylVU1Cf82Em5\nIilXck6mXHELeVeMMf2BrwBzgWGJvKe6uqFb+6qvbwKgtqaRysq6bm0jXbzeEsdlCnFqNqflOtrQ\nHH7spFwhTjteIcqVnJ7k6uoPQE9nrcwBvMBKYDEwxRizsIfbFBGRJPSoRW6tfRp4GsAYMxp4xFr7\nnRTk6sDj0SlBIiKxxC3kxpipBAYzRwMtxphrgSXAdmvt4vTGExGReBIZ7FwLzE5gvR2JrCciIqml\nMztFRFxOhVxExOVUyEVEXE6FXETE5VxXyHVfCUkHj654Ly7mukIuIiLRVMhFAL9+64mLqZCLiLic\nCrmIiMupkIugwU5xNxVyERGXUyEXEXE5FXIREZdzXSH3+zVNTEQkkmsKue4rISISm2sKuYiIxKZC\nLiLicgnds9MYMwl4Dlhorb2/3bKLgV8BbYAFbrHW+lIdVEREYovbIjfGFAGLgKWdrPIgcK219gKg\nBLg8dfFERCSeRLpWmoF5QEUny6daa/cEH1cCA1IRTEREEhO3kFtrW621jV0srwUwxgwBLgNeTF08\nERGJJ6E+8niMMQOB54H/ba093NW6ZWWF5ORkJ72P4uJ8AEpLCvB6S7oTM62cmCnEqdmclKuoqE/4\nsZNyRVKu5JxMuXpcyI0xpcBLwJ3W2lfjrV9d3dCt/RytbwKgtq6Rysq6bm0jXbzeEsdlCnFqNqfl\nOnq0OfzYSblCnHa8QpQrOT3J1dUfgFRMP1xAYDbLyynYloiIJClui9wYM5VAsR4NtBhjrgWWANuB\nV4AvAeONMbcE3/KYtfbB9MQVEZH24hZya+1aYHYXq/TpYpmIiKSZzuwUEXE5FXIREZdTIRcRcTkV\nchERl3NdIdd9JUREormnkOvOEiIiMbmnkIuISEwq5CLoB5+4mwq5iIjLqZCLiLicCrmIiMupkIuI\nuJwKuYiIy7mukPvRGUEiIpFcU8g1O0xEJDbXFHIREYlNhVwE/eITd1MhFwGNvIirxb3VG4AxZhLw\nHIGbLN/fbtlc4C6gDXjRWvuLlKcUEZFOxW2RG2OKgEXA0k5WuQ+4BrgAuMwYMzF18UREJJ5Eulaa\ngXlARfsFxphTgCpr7W5rrQ94EbgktRFFRKQrcbtWrLWtQKsxJtbiwUBlxPODwNiutldWVkhOTnYy\nGQEoLskHoKSkAK+3JOn3p5sTM4U4NZuTchUX54cfOylXJOVKzsmUK6E+8iTEHfyvrm7o1obr65sA\nqKttpLKyrlvbSBevt8RxmUKcms1puY4GP1+Ao3KFOO14hShXcnqSq6s/AD2dtVJBoFUeMowYXTAi\nIpI+PSrk1todQKkxZrQxJge4Cng1FcFERCQxcbtWjDFTgQXAaKDFGHMtsATYbq1dDHwDeDy4+pPW\n2k/SlFVERGJIZLBzLTC7i+UrgBkpzCQiIknQmZ0iIi6nQi7Si1rbfPj9uiCApJYKuUiKrdq4jxff\n2tnh9cbmVr5+73JuvnsZ+w4fzUAyOVG5rpCrLSNp4UnN9Q9b23w89MImnl6+laZjrVHLdh04Pn/4\nL6/YlOxPBFJ/QlDa6DKj4gYbtx4OPz5c28yw8uNfsc17asKPfT4/H++s5uV3duHtW8CO/bX84Iaz\nye3GWc8irinkImmVon7rg0caw4837znCsPKi8PNnVmwLP/5kTw33PL4u6r1/+vsm1nx8EIAF8y+g\nrKRPSjLJiU+FXCSFnnx9S/jx8nV7mTV5CF+7Z3nUOhNHl/HRjuoO7w0VcQC7u5oV6yu4fs54Rg12\n5jVDxDlc10eeaa1tPm79zXIWPLEu/Ly+sQWAhqZW1tpKzUo4yfj8fr5+73Ju/c3y6Nd9UHmkqcP6\nWRH98SMGFnPzlad1WOfBJR/x8a4j/OyRNXyy+0jKMwfy+fnHu7tZsmo7La1tadmH9A61yGOormtm\nW0UNU83ADstWbdzHsVYfH+6o5pd/fpfdB+tpafVx/22z+NPfN7F+yyGuu3gsnzp3JEvX7uHx1zbz\n3c+dxemj+2fg/0QS1oPBzv2HG2ht84Wf5+dlU1qUR13jMY42tUSte9t1Z1KYn8OmndVMnziIL1xm\n6JOXzdGmVpa9t4fPzhnHor9tjHrPb55Yx4yzhnc7X2fWba7k8dc2A/Dsyu08+P3Z5GQn1rb7aEcV\ndtcRJowZQF1tY/w39LKS0hrH5fJ4PMyampeWbauQx/D7xRvZVlHLyIHF/PSr50Yt27L3+IDVtora\n8GO7+wjrtxwC4KllW3lq2dbwsv94cj0/v/k8fvKntzn/jMHcfKXuvXEi2b6vNup5ed98cnOyOVjd\nyJsf7A+//p/fvYg+uYHBzD/+4OKo91w2bQSXTRuB3+9n1uQhfLyrmuYWH7VHj9Ha5mf/4aPkpjCz\n3VXNE0s3R732owff4hv/axJjhpR2+j6/388DSz7knU3BbqDVO1KY6sS3/0gTV58/KuXbPWkLeUur\njwef/5Cq2mZ+/OVzopaFCvSug/UsX7eX2WcPCy87UNVIlsfDpdOG88o7u8Ov3/9MdCsqkt8PP/7T\n2wCs2rifaacOZPLY8vDyxuZWcnOyyM7ysOtAPTVHjzFhRF/y87r+51m1cR//WLOb+f96Bt5+BYn/\nz0tKfbwrur+7pDCPQzWB1uDr7+0F4LbrJoeLeFc8Hg9fmXe8q+Xlt3fx12Vb+HhnFWeMKks6m9/v\np7G5lZfe3sWWPTXcdMWptLT5uPuxQNfg5LEDuH7OOO55bB2Hapq49/F1fO9zZzOsvIgVGyqYOXkI\nBX0Cn8Ote2t41x4MF/ELzxzKxLHlNDQ0J50r3UqK86mr79itlUlZHg9zzh1Fa3NL/JWT5L5CnoLu\n5/VbDnHf0++Hnx880sjAYCFsafWRneWhzRfY0Z9fsUyZ4KW0KI8DVQ1s2VvDwLICrp8znrFD+/KH\nZz8IxArm+vevncedfwwU7Umn9OcbV09i/sIVUfv/7VPvU943n1lnDuXc0wZyxwNvxcx5w9zxzD1n\nRMxlofnKAD/8rze5+OxhfPFTMW/+IWm2auP+qOdDy4v4/CXj+cnD74RfmzRmQLe2fcYp/fnrMnh+\n5TbOGDU1qff6/X7+8OwHrLXH7/1yx4PRn7WvzjuN0qI8Fn5rJn9fvYNnVmzjl39+N7x83eZK+pfm\ns7WilgNVgXsJ5OVk8eObpjGsvOiEvO53OpWV5lNZeRIXck+MPsyGplaeWLqZK6aPZMiAohjvCmjz\n+fjj8x9x1rhypp8+OKqIQ6C1/MD3ZgNwoKqBNp+fCSP6hQeZXn5nFys3VHC0KXCCx+D+hQCcc+pA\nHr59Doue2ci6Tyo5c+wAhgwo4p5bZ1Df1MLowYGfqDddcSr/WLObr316Is+u3M76LYc4VNPE4hXb\nWBwxJa29x17bzPTTB1Nc0PFH9eoPoovHsnV7OW/iICaM6Nfp9iQ1/H4/x1p99MnNprml4yDh0PIi\nhg8s5oa543nstc0U5eeQldW9Pvhh3mLGDCll294amo+10Scvfqt+1cZ9vPLOLiaO7h9VxCOde9pA\nbrlqYlSf+JUzRuHz+Xn9vT3UNgSKzce7ogdaP3PhKcyYOIhy/QJ0FNcU8lheXbOLNzbuw+6u5u5b\nz+90vV0H6nlnU+An4fTTB1OUnxMuyhBohYfsqawH4Bzj5TOzxnD3Y+t4+e1dUdtr34f4zevO4p9r\nd3HBpMA9Nsr7FVDO8Q/6hWcO5cIzhwJww6Xjw33pkeZMGUZVbTOnjSpj9JASFv1tI/WNLdz96Hv8\n/OZzw3/IqmqbuP+ZjezYX0deTha/vOU87O4jPPTCJn796HtcOWMUV88cw57Ket768AC5eTk0Nh5L\n6Hj2poKCPEflal+wOuPz+fnP5wKt3KtnjmHaqYEB8SkTvLz3SaBoDh0Q+EM/95wRDO5fyODg8+4a\nN6wv2/fVsvNAHVW1TXy0s5ovX27Izspi044qBvUvpH9p4FZ11XXN4V9qeyoDlwH43JxxnDXBywfb\nDrNifQVf+/REhnmLO+zH4/HwLzPHMG/GKNZsOkjjsVb+59XAValvuuJUzj1tYNzuPskMV/+rHAsW\n4CP1xwtCa1ugaySyBV8TsfzD7VVRRTykzedj1cb9PPLSx0CgJWRGRvdJeoDzJg7ismnR3R0D+xcy\n+6xhJKK8bwEP3z6H+sYW7vrLWvx+P9///NnhL2LIwm9dwG33vcHeQ0fZuK2KA1UNvLZ2d9R0tvn/\negbl/Qro3zef5ev2srWilhfe3MkLb3a8zockJl6Ld9UH+8Kt3Ofe2E7/4Ek7IwcVHy/kEScBTTql\ne10qkcYMCcwj//Wj74Vfe3/LIeZNH8UTwXnrv/3WTEqL8vhg++EO77902gg8Hg9zpgxnzpT4s19y\nsrOYEWyUnDWunH4lfaKmTIrzuLqQZwd/roZa1B9ur2LBk+u5bvZYrph+fGS4ruF4IV/w5Hog8IWd\necYQ9h0+ykc7qqmpPxYu4gDDvIEv4xXTR/LSW7u4/LyRfPbicSnLXlyQy11fn97F/1sW37pmMr9+\n9D1++9SGqGVzpw5n3oxR9CsOFJEsj4c7v3QOL721k6eWH58tM2F4X+Z/9mxqarp3n9R06l9WRFW1\ncy4cdbimid89/T4XTB7a5XpbI2YtAXyyJ9CSH1Cazx1fmMK+ww2UFKZ2ilmsPwa1DS3hIg5w26I3\nmDSmP7k5ga6Sa2ePpbG5latnjonZLZmo9g0McSZXF/L2c15DRfqp5VsxI8s4ZWigC6SusePgwvzP\nTGLSmAE8tXwLH+2o5nt/WB21vDT4ZbzivFGUlwYGJnvb+OF9KS7Ipb6xhUH9C/nCZRPIz83mlKGl\nMb+cV0wfxRXTR1Hf2MKxljb6l+YHBn3ynHfel9dbQmGOc1p5+Qn0PQfWi/7KhAY6B5TmM354P8YP\nT/0YRXFBLj+5+Tx+/tDbZHk8XHT2UJYFZ8Pk5mSFGzIfbK8CoKQwl8vPG6lW9EnE1YU8Oyuy+yR6\nCtQv//wuD98+B4D6ho6FfGhwcHTGxMG89NbxPvBBZQXM/8wZ4efFBblcnMDP0XTweDz87t9mBua0\nDypO+IJKxQW5EGOAVOKLd1ZuW1vs5aFGQ7pMmziYh354MbUNLfQtyuOLl5nwH+zahmP8/JHjM01m\nnD5YRfwkk1AhN8YsBKYTmPz3bWvtmohl84EvAG3Au9ba29IRNBaf7/iXavHK7R2WP79qe9Trk8cO\n4P3g1elKiwIt7uEDi7nx0gk8+o/AoM63rzszPCvFCTweD2OH9c10jBOeJ8Hra9YHz9T84Q1nh+di\nA+QlMEe8pzweD32LjnfbhP5g9y/N5+Hb51DbcIwNWw4xfeKgtGcRZ4n7m9sYcxEw3lo7A7gZuC9i\nWSnwfWCWtXYmMNEY03nHb4pFfnlWbKjosLx9cf9ScJ51dpYnqlvmkqnDefj2OTx8+xxHFXFxliP1\nzbz90QEgMHMp9Fn5wefPzmSssNLCPGZNHqpL4Z6EEmmRXwI8C2Ct3WSMKTPGlFpra4Fjwf+KjTH1\nQCFQlba0BE5QCI2ot/l8cdaOFmq5iHRH5BmcebnZ3PX16fj8fnVjSMYlUsgHA2sjnlcGX6u11jYZ\nY34GbAMagSestZ90tbGyskJyutFiKCkJjJ6/ayvxegPTsXJyE+/iv2Dy0PD70iGd2+4pp2ZzUi5/\n8DPpJ3Yun8/P4brjs58ykd1JxyuSciUnHbm6M9gZbn4Eu1Z+BEwAaoHXjTFnWms3dPbm6uruTYWr\nqzs+fzp06u2RJK5uduPccWk7ZdeppwODc7M5LVd1bcfPV6QVGyr427LAdL+B/Qp6PbvTjleIciWn\nJ7m6+gOQSCGvINACDxkK7As+Pg3YZq09BGCMWQlMBTot5Kng8/nJyvJwrCW6a+Vzl4zn7PHlvPnB\nfp59I9A/fuvVp3PuaRr8kZ6xEWd+XnX+6MwFEYkhkQnGrwLXAhhjpgAV1trQn5QdwGnGmND56OcA\nmztsIcVqjgZ+4h5rd52Lgf0K8PYr4MrzR3H+pMH82zWTVcQlJbKzj/eDJzrnXKS3xG2RW2tXG2PW\nGmNWAz5gvjHmJqDGWrvYGHMvsMwY0wqsttauTG/k4/cAaG6NbpGfOS5wBlx2Vha3XKVrfkvq9IkY\n12k6prvpiLMk1Edurb293UsbIpY9ADyQylDxHKhqoF9xnw4t8p6ciiwCdHqZ5KXv7Qk/nmq8vRRG\nJDHOO3c7AaETMSKvWnj7jVMyFUdOIJVHGlm+bm/4DM9VG/fxz/V7o9YJ3WhBxClc/YnctDMwr1dz\nwyVVPtx2mA+3HWZA33zOOGVA+JKwIfd9e1aGkol0zjUtcnWaSG+qPdrxWunZWR6K8l3d9pETlGsK\neWeyu3nnFZGuvPnh/g6vtfn8GocRR3JtIb87eJH9RG5oK5Ksj3ZUdxhMF3Eq1xZyG7yfZkNzx7v9\niCQrVkv7WLvprYncL1MkE1xbyENGD3bm9RTE/dp3r/RN8Z1/RFLF9YU83Rf0l5PX469Fn6T8k5um\nZSiJSNdcX8gP1TTFX0mkh+aeM5xCzVgRh3J9IQ/d8UcknSaP7XgDZBGncH0hv/HSCZmOICe4Gy+d\nwKQxKuTiXO4p5J1M3x05qLh3c8hJZ5zumSoO55pC3tnNcfsV9+nlJHKyyclxzddETlKu+YR2dkKd\nt19B7AUiKdIn1zVfEzlJueYT6u/k8qIiqdDVmfflfdVYEGdzTSH3qZKLiMSkQi4i4nKuKeR+X8dC\nfta48gwkERFxloROVTPGLASmE7gR1rettWsilo0AHgfygPestbemI2hTjCvR3TTv1HTsSiSsUHcD\nEheI2yI3xlwEjLfWzgBuBu5rt8oCYIG19lygzRgzMvUxob6hpcNrpbqIkaRIrLHO00eXcf93Luz1\nLCLJSqRr5RLgWQBr7SagzBhTCmCMyQJmAUuCy+dba3elI2h2dvRX7YuX6YxOSZ/8vGy++7mzMx1D\nJCGJ/G4cDKyNeF4ZfK0W8AJ1wEJjzBRgpbX2jq42VlZWSE5O8td1LiyMPvFnwpgBeL3OuYStk7K0\n59RsTsqV0yc36nmfvGxH5QNnHa9IypWcdOTqTgegp93jYcDvgB3AC8aYK621L3T25urqhm7sEhoa\nmqOe19U2UVlZ161tpZrXW+KYLO05NZvTctXUR3++crI8jsrntOMVolzJ6Umurv4AJNK1UkGgBR4y\nFNgXfHwI2Gmt3WqtbQOWAqd3K2Uc7WcfjtINJSSNuvOrUSRTEinkrwLXAgS7TyqstXUA1tpWYJsx\nZnxw3amATUfQSFONlwLNJpBUandq5/VzxmUoiEjy4lZDa+1qY8xaY8xqwAfMN8bcBNRYaxcDtwGP\nBAc+NwLPpzMwwA1zNdApKRbxk2/m5CE6R0FcJaFmrbX29nYvbYhYtgWYmcpQsYS+ZznZHspKdMVD\nSa3InrucrC4uvCLiQO45szP8VdOXTFLvWKsv/DhLhVxcxjWFXCSd+hUdP7lMhVzcxj2FPNgg7+py\noyLdlZd7fJZKtgq5uIx7CrlIL1GLXNzGNYVcPeTSW/rkag65uItrCnlLcDAqJ9s1kcWlpk7wZjqC\nSFJcUxVDX66rZ47JcBI50fUvzc90BJGkuOb0yFNHlfHUXVdSV9uY6ShyglMfubiNa1rkAPk6LV96\nQa6678Rl9IkVaUctcnEbNXFFgr50uaG2sTXTMUSSpkIuEjT7rGGOvY61SFfUtSIi4nIq5CIiLqdC\nLiLicirkIiIup0IuIuJyKuQiIi6nQi4i4nIq5CIiLufx+/3x1xIREcdSi1xExOVUyEVEXE6FXETE\n5VTIRURcToVcRMTlVMhFRFxOhVxExOVcc2MJY8xCYDrgB75trV3Ti/ueDTwFfBh8aSNwD/AXIBvY\nB3zRWttsjLkRuA3wAQ9aax9KU6ZJwHPAQmvt/caYEYnmMcbkAo8Ao4A24CvW2m1pyvUIMBU4HFzl\nXmvtCxnIdQ8wi8Bn/lfAGpxxvNrn+hcyfLyMMYXB7Q4C8oFfABvI8PHqJNe1OODzFcxXAHwQzLWU\nXjxermiRG2MuAsZba2cANwP3ZSDGP621s4P/fQv4OfB7a+0sYAvwVWNMEfATYC4wG/iOMaZ/qoME\n97OIwIclJJk8NwBHrLUzgX8nUEDSlQvgjohj90IGcl0MTAp+fi4HfoszjlesXJDh4wV8GnjXWnsR\n8FngP3DA8eokF2T+eIX8X6Aq+LhXj5crCjlwCfAsgLV2E1BmjCnNbCRmA0uCj58n8I9zHrDGWltj\nrW0EVgEXpGHfzcA8oKKbeS4BFgfXfS2FGWPliqW3c60Args+PgIU4YzjFStXdoz1ejWXtfZJa+09\nwacjgD044Hh1kiuW3v53xBhzKjAReCH40mx68Xi5pZAPBiojnlcGX+tNE40xS4wxbxhjLgWKrLXN\nwWUHgSExcoZeTylrbWvwgxApmTzh1621PsBvjMlLUy6AbxpjXjfGPGGMKc9ArjZr7dHg05uBF3HG\n8YqVq40MH68QY8xq4DECXQEZP16d5AJnHK8FwP+JeN6rx8sthbw9Ty/vbzPwM+Bq4MvAQ0SPL3SW\np7dzxtvOaJojAAACIUlEQVRvJnL+BbjdWjsHWA/8NIn9pzSXMeZqAgXzmz3cfzpzOeZ4WWvPJ9Bn\n/z/ttp3R49UuV8aPlzHmS8Cb1trtSe4nZbncUsgriG6BDyUwgNArrLV7gz/r/NbarcB+At07BcFV\nhgUzts8Zer031CeRJ/x6cKDFY609lo5Q1tql1tr1wadLgDMykcsY8yngTuAKa20NDjle7XM54XgZ\nY6YGB88JZskB6jJ9vDrJtTHTxwu4ErjaGPMWcAvwY3r58+WWQv4qgdFpjDFTgAprbV1v7dwYc6Mx\n5nvBx4MJjJr/N3BNcJVrgJeBt4Fpxph+xphiAn1dK3sp5mtJ5HmV432znwaWpSuUMeZvxphTgk9n\nExjV79Vcxpi+wL3AVdba0GBUxo9XrFxOOF7AhcB3g3kGAcU44Hh1kuuBTB8va+311tpp1trpwJ8I\nzFrp1ePlmsvYGmN+TeAf0gfMt9Zu6MV9lxDok+sH5BHoZlkH/JnANKidBKYMtRhjrgW+T2Ca5CJr\n7aNpyDOVQJ/caKAF2AvcSGAKU9w8xphsAh+48QQGKG+y1u5OU65FwO1AA1AfzHWwl3N9ncBP7k8i\nXv5ycF+ZPF6xcv03gS6WTB6vAgLdhyOAAgKf93dJ8PPey7nqCUwFztjxapfxp8AO4BV68Xi5ppCL\niEhsbulaERGRTqiQi4i4nAq5iIjLqZCLiLicCrmIiMupkIuIuJwKuYiIy/1/kKXM8+1RSQgAAAAA\nSUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f239c54c7b8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(trace['step_size_bar'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "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": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f23417d0dd8>]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWkAAAD4CAYAAAAuNhccAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmUG9Wd6PGvutWt3hd3ywteMTbXLGEJe0gAOySQQMLk\nTQiPyUbCvGQmy0wymUkyk7wEZvIIgSHLmcxLIGyHbV4gmQQIJBDA7LuNzWK4xsa722713q3ullrL\n+0NLbyV1SaqSruTf5xwfq0ulql9dVf30062qK088HkcIIYSZqkodgBBCiMwkSQshhMEkSQshhMEk\nSQshhMEkSQshhMG8Ti8wEBjO+3KR9vYG+vtHnQzHERJXbiSu3EhcuTE1LigsNr+/2WM13ahK2uut\nLnUIliSu3EhcuZG4cmNqXOBObHNW0kqpJuA2oB3wAVdqrR9yPBIhhBCz2KmkLwO01not8HHgZ65G\nVIHkhqG5SRsJYc1Oku4BOpKP25N/C5sGg2Eu/9F6HnhuZ6lDMVZ4IsrlP1rPXY9sLXUoQhjHY6eC\nUUr9CVhFIklfoLV+PtO8kUg0bnKfUbE9sXEv/37nBgDuv+6iEkdjpp1dQ3z139cD0kbikGZ54tBO\nn/SngN1a6/OVUscDNwEnZ5q/kLOufn8zgcBw3q93SyFxDQ2PpR87vW2V0l79fcH0Yze3p1Laq1gk\nrtwVEpvf32w53U53x5nAQwBa683AYUopKZWFEKII7CTpbcBpAEqp5cCI1jrqalRCCCEAezezXA/c\nrJR6Ijn/37gbkhBCiJQ5k7TWegT4RBFiEUIIMYNRdxwKIYSYTpK0EEIYTJK0EEIYTJK0KD3LS/iF\nECBJWgghjCZJWgghDCZJWgghDCZJWgghDCZJWgghDCZJWgghDCZJWgghDCZJWpScXCYtRGaSpIUQ\nwmCSpIUQwmCSpIUQwmCSpIUQwmB2fpkFpdQngW8CEeB7WusHXI1KCCEEYKOSVkp1AN8H3gtcCFzk\ndlBCCCES7FTS5wKPaK2HgWHgC+6GJIQQIsVOkl4BNCil7gPagSu01o+6GpUQQgjAXpL2AB3Ax4Dl\nwHql1HKtddxq5vb2Brze6rwD8vub836tm/KNq6V5sOBlZFMJ7TUandyV3N6eSmivYpK4cud0bHaS\n9EHgWa11BNiulBoG/EC31cz9/aN5B+P3NxMIDOf9ercUEtfQ8Fj6sdPbVintNXWfcXN7KqW9ikXi\nyl0hsWVK7nYuwXsYWKeUqkqeRGwCevKKQgghRE7mTNJa633Ab4DngT8CX9Vax9wOTAghhM3rpLXW\n1wPXuxyLEEKIGeSOQyGEMJgkaSGEMJgkaSGEMJgkaSGEMJgkaSGEMJgkaSGEMJgxSXrrngE2vmV5\nEyNPv9pFdwF3Mm7dM8DrO3rzfr0d0ViM9Rv3MhQM5/X64PgEj23cS3giavn8tr2DvLrdehvi8TiP\nb9pH39D4rOeefb2Lg335t52Vg32jPPLynmnxPr/lAPt7grPmHQtFeGzjXsbDkYzLs3qdHfsCI7yw\n5WDG5+PxOE9u3k/PwFjGeYQwna3rpIvhrj9vJRiKcO3fvmfa9F0Hhrn5wTep8ni48Vtr81r21Xdu\nBODmb68rOM5Mntzcxe0Pb+XZNw7wnU+fnPPrb/3jW2zQAYaCYf7ifStnPX/VHRsA623YtK2H2/6k\n6Wip49ovTbZfV2+QG//wZsbX5eufb3g+/XgoGObsExZzw31bLNfzm8e3s/6VfewNBPnMecpyeb/4\n/et5xfG/b3oRgKOWt9PSWDvr+a17Brj1j2/R4PPy86+fldc6hCg1YyrpSCxOeGL2jYzDY4nKNBa3\nHM/JGD2DiWpt14H87tvf0z0CQHd/7lXf4EiijXpnVNLB8czVq1O6+8eyVsldvYkqeX9gxLUYwhHr\nbx+DyW81oyH320EItxiTpD0AWCRis3PzpHScnlJGYRyPJ9Eebo4j4JE2FxXMmCQtx1llSuZo4oZ/\nExLCVOYkacDqOC6XQzsVp0c+bKapSjaImzk6U5vL54KoBMYkaQ/lk5CzkRw9nSedpCvh3RWi+IxJ\n0uCxLH3K5tgulziLLFXlxkrQPnF5U0QFMCZJV0w3QaVsh0OqilBJeypm5xFiNmOSNGQqRsujGpKq\nzdrkicMSrFzeElEBjEnSHsqoayMLuRxsOumTFqIwtpK0UqpeKbVdKXWZa5F4wKr0KZdju1ziLLZS\nVtLylohKYLeS/i7Q52YgFVOBVshmOCV9M4t8igmRlzmTtFJqDXA08ICbgQyNhhkLRXnujQPc9chW\nfnL3Zq66fQOv75j8bBgcCQEQCke58+GtHMgwcNCfX9rDyxkGa0rZ3xPkzj9vZSLDLcU9A2Pc8bBm\ndHwizy2a27Z9g9yzflvBXQEvvpl5kKGZJiJR7vzz1lmDGq1/ZR/PvX6goDisVM1RSTuRvDO2n3wu\niApgZ4Cl64CvAJ+1s8D29ga83uqcA+kfTiTgX92/Zdr0bfsG04/vfuId/uWyU/nNY2/z6Ma9vLaz\nj5u+84FZy/qvR98G4P7rLpo23e9vTj/++s+fZnAkzKpl7Vz43tkDGl1910a27h6gpbmOy5fOm/Za\nK/X1iQF+qjyeafO2NE/GP3MZn7/6MQDWnbac6urE56XPV5N1XTOf8/ubeWv3gOXzvcGJWdPvf+od\nHt2wlw06wO1Xnp9+/vaHNAAfXbs622bO4vPVMG9e46z1pP6v89UAUFXlsdyujXr6h+lc7Wxl3rwm\n/O31s6Y3t8xul3yWXwwSV25MjQucjy1rklZKfQZ4Tmu9QynrEcxm6i9gSNG59PSPEggM09OXqAJ7\nB8YIBDIPaDTzual/pwYl6ukLWi6jNzm8ZSC5rmzrARgdTSwvHo9Pm3doeHLApEzL6OkZIRpNjG4R\nCk3Y3ia/vznrNvYPjM6a3t2bGOhoYCRkuZ65tnOmUGiCvr7JqjwQGJ4WVyg5+FI0GrNcdveMgZdy\nXT9Ab+8IRGYPojQ0NL3trdrLBBJXbkyNCwqLLVNyn6uSvgBYqZS6EFgChJRSe7XWj+QVRYFS32on\nb5Bw8/tsfrczm37Jbqm6hjOt1omrPuTyR1HJsiZprfUlqcdKqSuAnaVK0DB5MKZPMhZhPAi7OaRc\nEkWxL4UrymdWpi7p8nhLhMjKmOuk7ZhZSbtaR6fXketa3E1LhSbZkiUuN08cFrwEIcxl+5dZtNZX\nuBiHLZNJ2v36zJPrJ0GRMkWcwj4Gij5aXxHWI0laVLIyq6RT3R3uS60j10rP9dgKzEiTbVjczvNM\n30gcOa8g/RqigpVXkk7+X4wqMNdqvVhpwqm+72JV0nOuxokcnWm65G5RAcorSc/slM42T4HyHWLT\n7eRX6OalKlcn47QTU8b7TZx4uyQZiwpWVkk6dTBWZUkwTlVPntwv7ygPLvTrZ9/07Otx4ptBxkq6\nbN4UITIrqyRtp6p1rDsgtTzDjvNC40l/GSk8lNkLLe5LpyzDsDdJCAeVVZJO1UxVWbs7srw6h4M5\n94s7itcrXdirU90dxamk51qNq/lVcreoAB6nq5BAYDivBabGsTBZe7MvPcaI23w11cxvr2dP98is\n56qrPEQzfK248D0r+MOzOx2NpaWxlqFg2Na8Kw9rocZbjd7d72gMAJ/64JHc8fBWx5ebi9VLWnl7\n7+DcMxagxlvFRCTm6jqmynW/LtZx4AEWdjTQ1evOUBMdLT5qvNUZB2rL1ZplbVz792cXclu4ZUkj\nSVoIccjKVvDk4/7rLnI8SZdZd4cQQjhnxSJzR9NLkSQthDhkZTu/lY+Yg1V5iiRpIcQhy/Ek7cKZ\ncEnSQohDVlW2my7yIJW0EEI4yOEcLUlaCCGc5HG6kpbuDiGEcI6cOBRCCIM5naSdvOY6RZK0EOKQ\n5fiJQxe6O2z9MotS6hrgfcn5f6i1/m/HIxFCiCKriBOHSqm1wLFa6zOA84GfOh6FEEKUgPOX4Dm6\nOMBed8eTwMXJxwNAo1Kq2vlQhBCiuJzuk+4ZGHN0eWAjSWuto1rrYPLPy4EHtdbRTPO3tzfg9zfn\n/M8ULY21ltNrvdOb6rhVnbPmsZrmtvnzGqb9Pa+lLuv8all7+nFTfU3G+Rrra/jwe1bMuf4jlrSm\nHy9b2MwpRy9g6YKmWfP5au19rmdrQ2/13DVFpm1qb/bltJypqgustj70nhV846/ebWveE4708+Ov\nnZXzOo4+fF7Or0lZsagFgCOXtfHuNfPnnL+2Jv8abVFHI8evtn+c/OXaVenHZ7xrEeefsSL9d423\nsFNqC+Y18OH3ruRDU5Y5VUPdZG9wvc+6Z3ixf/q+vnRh7rlvrhxoexQ8pdRFwL8AH9RaZxyrsdBR\n8L5/2SlceetLAFyybhXnnbps2nw//vUmXt/Rx8rDWvjuZ062XNbe7hG+d/OL6b9v/va6vOO5+oun\nc8yRCwgEhtPTvnbxcRx3RPET8szYfvK1s2mtM+9Ljd/fnHUksFT8f/G+w/nomYfPev7lt7r5v79/\nHcjvvcs3rlKRuHJjalxQWGyZRsGze+LwPOA7wPnZErQTpn77sPz8SP+qVebPAicvULf6OuTCuYG8\nyM9DCVH55kzSSqlW4FrgXK11n/shTbJKxKmkme0LgJPnAixPLBiSGyv1V6Oc/NUYIcqdnUr6EqAT\nuFsplZr2Ga31bjcCmnqAZrvmMGuSdjBLWyUMU35Tz5Q4nCY5WohJcyZprfUNwA1FiAWYu7sjXUln\nKWedPGNrlfBNSY2mxOE0SdJCTDLujsOpx2e2SjFbEenkQW7Z22FIBRsv3s/gFZXH2d8yF6KsGZek\np2ZYq1yY/hXvrH3SLlfSZuToij1xKJW0EJOMS9JTj0+rPmmPje4OJ1klfFNSoykfFk6TJC3EJPOS\ntGcyUWerpLNlSievDrBM0oZkR1PicJpc3SHEJAOTtCdrtZw6fLNd+eFon7RFC5mSG02Jw2mSooWY\nZF6SJnu/s50qy8lKzPISPEM6PEyJw2lSSQsxybgkXVtTzWJ/IwDNDbPH0ehoTYxN4W+rz7iMmmp3\nThzW1iSaq8GXecyLYqqrtXXDqHEWJscbabF4f2H6mAlCHOqqr7jiCkcXODoazmuBJ6zq5JgjOjli\nUQvHHdFBc0MN5560dNbgNkcsbsFXU81fnLUSX4aBXmq81bQ3+6irreaLHz0m46BJ2Ry1vJ0jl7ax\nfGEzjY0+RkfDnLJmPvNafJx57MKSVnsnKz8L5jVwzklLGR0NlyyOTFLtlckJqzppaaxh3buXWF49\n09ZUS0NdDZesW5XXe5dvXKUiceXG1LigsNgaG31XWk23PcCSXfkOsATmDpwiceVG4sqNxJUbU+MC\ndwZYMq67QwghxCRJ0kIIYTDHuzuEEEI4RyppIYQwmCRpIYQwmCRpIYQwmCRpIYQwmCRpIYQwmCRp\nIYQwmCRpIYQwmCRpIYQwmCRpIYQwmCRpIYQwmCRpIYQwmCRpIYQwmCRpIYQwmCRpIYQwmCRpIYQw\nmCRpIYQwmCRpIYQwmCRpIYQwmCRpIYQwmNfpBQYCw3n/aGJ7ewP9/aNOhuMIiSs3ElduJK7cmBoX\nFBab39/ssZo+Z5JWSjUBtwHtgA+4Umv9UF5RzMHrrXZjsQWTuHIjceVG4sqNqXGBO7HZ6e64DNBa\n67XAx4GfOR6FEEIIS3aSdA/QkXzcnvxbGOpg3yhv7x0odRhCCId44vG5u5CVUn8CVpFI0hdorZ/P\nNG8kEo2b/HWk0n3kG/cCcP91F5U4EiFEjvLuk/4UsFtrfb5S6njgJuDkTPMX0qHv9zcTCAzn/Xq3\nlGNcpYy3HNurlCSu3JgaFxQWm9/fbDndTnfHmcBDAFrrzcBhSikplYUQogjsJOltwGkASqnlwIjW\nOupqVEIIIQB710lfD9yslHoiOf/fuBuSEEKIlDmTtNZ6BPhEEWIRQggxg9wWLoQQBpMkLYQQBpMk\nLYQQBpMkLYQQBpMkXaHs3EkqhDCfJGkhhDCYJOkKJXW0EJVBkrQQQhhMknSlklJaiIogSVoIIQwm\nSbpCxaWUFqIiSJIWQgiDSZKuUHKZtBCVQZK0EEIYTJK0EEIYTJK0EEIYTJJ0hZI+aSEqg52fz0Ip\n9Ungm0AE+J7W+gFXoxJCCAHYqKSVUh3A94H3AhcCF7kdlHCClNJCVAI7lfS5wCNa62FgGPiCuyGJ\n3z6xncWdjZx+zMJSh1I2Xt/Ry0Yd4NPnKTweT6nDEcIxdpL0CqBBKXUf0A5cobV+NNPM7e0NeL3V\neQfk9zfn/Vo3FTOuB57bBcBHzlk957yZ4urobMZXk//7UKhiv48/vvoxAC5au5pVS9oyzif7V24k\nrtw5HZudJO0BOoCPAcuB9Uqp5Vpry+/T/f2jeQfj9zcTCAzn/Xq3lCquudaZLa6ewDC1JUrSpXwf\n+/qCBHzW2y37V24krtwVElum5G7n6o6DwLNa64jWejuJLg9/XlGIOTn1iyrSIy1EZbCTpB8G1iml\nqpInEZuAHnfDEkIIATaStNZ6H/Ab4Hngj8BXtdYxtwM7VDlWAUspLURFsHWdtNb6euB6l2MRIMlV\nCDGN3HFoGKfGgZbxpIWoDJKkhRDCYJKkDePUmBuH6tgdch+LqDSSpIUQwmCSpA1zqFbAQghrkqSN\n49CJw0M02XuQ/g5RWWxdglcM96zfRiQOl65bxVgowlV3bOAj71nBqUctsL2M8XCEq27fyAVnLOe0\no7O/7tXtPdy9fjvfvPREWhprc4r1/md3snXPAP/wiePzHsxnKBjmmv96hU+sPYLjjujMON9YKMIP\n79jAR848nFPWzM9rXSnPvNbFn1/ew7986qT0LeO9g+P8+O5NfOqDiqOWtxe0/ExC4ShX3bGB809b\nxhnHLCQWj3PNXa9w4upOzjt1WcHL/9MLu23PG56I8jfXPTFr+sJ5DRzomxzS4LgjOvjaxcfzzv4h\nfnDby5bLetfKDt7eO8B4OEprUy1nHL0Qvaefhroa3tjRB8Cnz1OsPXExAHf9eSs7uobYvn8IgLOO\nX8STm7s4+agFfOmiY9LLffrVLm5+8E3OPWkJu7tHOHF1J4s6GvnpPZvTsQ2PhtnRlbj9uKm+hpGx\nCQBOWNXJ3338OJ7fcoAb7tsCwIdPX47e08/2fUNzts8VnzuFK255adb0j521kodf3E1wPMI3Lz2R\nNcvbeWtXP7c/rPn6JxLt9Mt73+DE1Z1EY3Fe3d4LwDcuOYF6n5cf3PYyS/yN/Ovlp/HCloP84dmd\nfO3i49nRNcR9z+xg4bwGPB4PL73VDcC3/upEhkcn+PVj26iqgvccu4iNWwPsCwRZvayN/qFxuvvH\nAPjomSu475md0+Jd1NHAh09fzk0PvJmeduXnT+XGP2xJ74efT473ssTfyN5AEIDGOi+joQhtTT76\nh0MAfO7Da+juH0uPqeOtriISjfG5D6/hvqd30DsU4vILjuLMdy2as33zYUyS3ry9l5GxCS5dt4pN\n23rYFwjyy3vfyClJv7q9l72BEa6/7405k/RP73kVgKde3c8FZ6zIKdbfPfkOkKhW8z1R9dSr+9nf\nE+Sn97zKzd9el54+swLeuDXA3kCQX/z+dU6ZMt/cZpfSqR12654Bjl3ZAcBDL+6mq3eU//zv1/j5\n18/KeTvseH1HH3u6R/jV/Vs445iFDAXDbN0zwNY9A44k6bvXb7M979Y9A5bTpyZoIJ1kbn7wTavZ\nAXjtnd7048GRMH96cfaHxe0P6XSSfmTD3mnPPbm5C4CX3zwIU5J0ap2p+WfGnIotJZWgATZtS9wM\nnErQAA8+vyvjNsx0zV2vWE5P7fOJuPezZnk7//m71wiOR3johT08ujER6ytvT78Z+bpfb+KwzkaA\ndCK88Q9biMbibNwa4L8efXvacym/uPcNhoLh9N/3Pr0j/Vjv6p8278wEDdDVOzotQQNcfecGxkLR\n9H6YMnXdwfEIQDpBA9zy4FvTlhOJxmZNv+mBN11L0kZ1d6QSVLl8YS3kWuRM3REzJ7vdbVGcXpHp\naylpV0yOO5dTY6mUCztFR6zANonGEq/P2rYutHu5vpXGJOnEvlFYli72m+DK+py6BC/bkxbte6hc\nuiZ91tmZMha3wYdW0RmTpPGU7yddPuxuqut3DhahzY16X3PIQYdaFQ1QVcQcna11TS6Ais2YJO2Z\n9tiMT/O5uLEjOXZbeJbFlEv75mOuQjCXLS/TY7ogplTSbijXoRKMSdIweVCUz35SwJueIYvOmux6\nIV38HbeUFWpOu1bcsG8BRVDMYy97l7QrFVBZMihJl19/R5mFO8myT7psPhkLk8N2lmvlVYhK3g/K\n9d00Jkl7POXXiIXEa79P2l3FaHOT3tecujviZsVeDMXsk87GlUI6PvVx+byz5iTpqY/L5dPc5R2p\nsOVkXlCZtK4rymXXKkQhCciUY8+dFDq51DLK0eYkaTi0rpPO/NJiX0dYhFXMOCLK5QCJxymfYKco\nJOKi9klnv77D+fVNraTL6DuSrSStlKpXSm1XSl3mWiQeKLcvlyaf28i2HKvj0JACynW5VYrltT+m\nFRC2MZW02yfMy+ittVtJfxfoczMQDx4mIrH0LZcpE5EYE5EY8XicaCzxXCweJxKdPq1cpOK13Sft\n8t7k9r5q9f5YVTEmvo8TEWfqLTvb5uS+PBHNfzn57G9z3YE485hOvy6W+XWuXCYdt35sujnH7lBK\nrQGOBh5wM5BdBxODxXzh2sc5YdXkgENf/PfHp81387fX8dc/Wg9AlcdDLB7nxm+tpaoEFUCub/TI\n2AR/97OnWHvi4oyDOjm282QtpWe3lRut99Sr+7nlwbd4/7uXTH9iRmz3Pb2D3z+9g6u/eDrz2xsK\nWudclWAuu8lXfvqkI98w/tc1j3Pdl8/MOs8v7n2Dl9/q5pffOLvg9f2txQBSdgUGxnN+zfpX9mV9\nPjUQ0ky/feIdy+mQGJTLaVM/clP5xkmbtvXwAX+z48u1M8DSdcBXgM/aWWB7ewNeb3VBQaUGibHi\nn9IIqU/wefOaqPFW0dIyaDlfNo2NPlvzWs3T0dlEU32NrfUA7H87ACR26ks/qCyX7ZsyqIzf30xT\nU13WGDLF3tHRRFuzz/K5trb69Ovq6hLxV1VV2W4zO/z+5vTAN8+83jVterSqatrfv08OnrMzMMox\nR9ofUMtKe3tD1u1ob2vMaXlOfWgeHAxlff7l5Ohv1T77+1Op1Plq8PubqcrjUhAn97FCvLazf+6Z\ncvTE5v184IzDHd/GrElaKfUZ4Dmt9Q6lVLZZ0/r7R+eeqQAHu2cPtxgIDFHjrWZoaHzKNHuflMFg\naM55/f5my3l6eoYZq7N/UPUPTLZNMDh50E5d9vBoeNr04ZHM25QpLoCe3hEmxsOWzw0OjKVfNz6W\nmCcWj9tus7mk4kp9zZ2a5wKBYXqn7CPTtn14vOAY+vqCNNVY9+L5/c0MDLq7f2YyNGxdTc7U0zvi\nciSFGw9NEAgMZ+2uyKTb4vgtlqkfuKNB62Mjk8Y6b3qEvExCocTz+e7DmZL7XJX0BcBKpdSFwBIg\npJTaq7V+JK8oHGC1Y5SqOzPnKiufqqych+7IcLVOpnUW49pVM06LZVFGfaX5MGXzyunqjqxJWmt9\nSeqxUuoKYGcpEzRYH8iFDp1YLHaSkGNbku06aas7Dp1a77QQ4pbryzhMqwMbP+cyDM/S5bIv582Q\nzXOjmfP4YmGLUddJ22FVNZdqx8618ps6d8aXxrP+6Tg3m25y0dMzY6Z2K8a7WMmDS5UDUypYV6Jw\n6WCy/cssWusrXIkgR1YJOZ++MSfk3NthEbvdroBc5bwcF3JXppuTMlfShW/9XEnAlCSRSYl25aIx\n5YuCG11rbvW6ll8lbdndUYJAIOdMOP2Opzxe5Ar3lp+xuyPj/K6FMvfKTWFKFqtw5TRedfklacsT\nh+VSSVtMtFtmOrGu1CqLdJ10rhWzE1XuXM1negqs9BxtysBGbsTh1raVX5K2aIeSvfE590lP65XO\nME/2v51WjD7pWR8KNvvjXWFGjsjI8PAKZsqNpWVUSJdhkq6QSnqyv3bmSbXMr3FDxkTqxLIzBJ/p\nRK8TJ4DnrqTNToOmVJrZFBKiKbf/x13IGVJJJ1kl5GjJru7IdX7zD0AnpT+IinkJ3lwnDg1/C0yP\nDwqrGKOGnBl1Iwq3Ns321R2m2KC7Z0176c1ujljcypu7JseAen7LAXw11dTWVNPg8zI6HmEsFKHO\nV01N9eRn08BImI1bAwyNhqn1VrFqcSvtzT6GRycAeHvvILV1fRCJgodpt4HvOjBMz+AYDT4vI2MT\nNNbV4PFAY30NkUiMvuEQE5EY/vZ6RkYn2LJr8lbUzclb32PxOC++eZBYPM685jr6hsenzbOne/IO\ntHf2D+GrrWYoGGZwJERtXR+eaDSxjXVe9vcEp722oc5Lc0MtQ8HwtArmZd3NwEiIUDjKc68fAKB/\nOMQzr3VRV1tNvc/LmuXt7O0eIR6H3qFxQuEoSxckbr8PDIwRi8WJx6GztY7RUIRIJEYkOe2M5jp2\ndA0RmkiMv5BqS0h8yAYGJu++23Vg8u6sxzftY9mCZmq9VRzsH6WlsZbR8QgtjbX42+rZfXCY6qoq\nlsxvpKWhls3be+hsqZ+2LzyxaT9joQjtzXW0NtbSOzjO8FiYaDTO6Y0+V8ZssGPb3sG5ZyLxPphu\nZ9cQz7zWNecdeFZ2dJWm/Wey+36k2NnWngF7d5XmyuN0dRcIDOe1wM9f/ZijcRSis7WOnsHcB5qp\nJCcrPy/rgOPLnd9WT7dLO7MQpXb/dRcVclu4ZZ9j2XV3FMOhnqABVxI0IAlaiBxJkhZCCIe4cd5J\nkrQQQjjEjZOHkqSFEMIhblwOLElaCCEc4sZgb5KkhRDCIVJJCyGEwSRJCyGEwaS7QwghDOZGJW3r\ntnCl1DXA+5Lz/1Br/d+ORyKEEGWuJN0dSqm1wLFa6zOA84GfOh6FEEJUgFJ1dzwJXJx8PAA0KqWq\nHY9ECCHKXCRagu4OrXUUSA2vdjnwYHKapfb2BrxeyeFCCPOtWtrGtj0Dji2vd3CMow/vcGx5kMNQ\npUqpi0h6o7n3AAAMW0lEQVQk6Q9mm6+/f7SggM49aQmPbNib/nvZgiZ2H0wM13nRew/nrV396GSj\nXnruag5f1MILWw7yaPI1px41n9VL2nhrdz8dLXU0N9SwqKOR197ppcHnpa3JR0Odl3e6hgiOTbB0\nfhP7AkHwwML2Btqafdz6x7cA+NhZK9mxf4iWJh+rF7fwTtcQ3qoqdnQNsW3f5FCHHztrJc31NYyG\nEsMZtjTUEovHqamuoq62mq6+RJuMh6MctayNnQeH6RkYZ16Lj1gcvFUe9gZGaPDVEInFiMcTYzC3\nNtYSjcU5rKMRgP29QVoba+kbDlHv8zI2EePIw1qIxeMMBsME+sfY0z3MGzsTQ6JefsFRVFd70LsH\niMbinHSkn76hcd7c1U9LYy0ej4fBkRALOxrxt9axNxBk2YImBoPhRGyxOL2D44yHowTHI8xr8VFd\n5SE0EQPitDX5CE/E6B0aJxaL09bsYyISY/Xyeew/OET3wBgTEzFWLGpmKBimq3eU41d1EInGGRgJ\nMRQMc9TydkZDEd7eO8jASIhl85vxt9Xx9KtdLOpspLHOS1N9DR6Ph7raaprqa6jyeBgemyAwMEZz\nfQ1N9TVUV3vYumeAwMA4DXVe5rfV09Fax1AwzL6eIB0tdRyzys+2XX2EI1F8NdUExyOEJ6Is8Tcx\nNBpmdDxCNBYnNBElGoun94+RsQlGxibSw7j2D4eIx+HIpa2Mh6P0D4cYGAlRV+vlsM4GotE43QNj\n1HiraKqr4UD/KAvaG4hGY6xY1MIRi1vZ8FY30Vic0fEIJx29kO17+ghNRKn3eQmFo4QnYrQ1+3jt\nnV7mt9ezdH4T1R4P+3uDVHk86f3VV1tNeCLGO11D+Nvq2N8TZIm/iR1dQ0Siceprq+loraOu1ovH\nA3u6RzjYP8oRh7Vy+KIWNmwNsHpxKzXeKiLRGHsDQRp8XqKxOEsWtjA0PMZ4OMqB3lEWzGvghFWd\nbNrWg97dz6olrby6vZcTV/upq61mLBRhwbwGdh8cpq7WS2drHdVVHroHxnh1ey+LOxuZ317PeDjK\n4Yta2LZvkFg8zpFL2tigu2lt8rGzayi9bY31NRyxuJWR0TCjoQgL2hsYDUUYnYjhq4JdB0dYOr+J\nXQeH2X1wmKXzmwCYiMQYGA6xdH4zKxY1s/vgCJ7kMMO13ipGxiMMjIQ489iFrDyslfuf3cmOriGq\nqzysWNjM23sHmdfiY4m/iZ7BcSYiMfb3Bnn3aj8AewMjdLbWsad7hKFgGH9bPYs6G2lrrGXN8nmF\njIJnOd3WUKVKqfOAfwPO11r3ZZu30KFKr/rC6dz92DY2beth2fwmzjttGb+6fwsAN397HV29Qb7z\nqxfSf898/dRp+Uot66ZvrcXj8eD3N89q+O/f/CJ7ukc4YVUnf/fx4wpeZz6s4npq835uSX7IONEW\n+bCKywQSV24krtwVElumoUrnrKSVUq3AtcC5cyVop8WZ/aseVVVu/GSqtWw/KeXCr005wo2fwRJC\nlI6d7o5LgE7gbqVUatpntNa73QjIw2QCjMdJfLWbYubfpZL6bUI3zuYWwpDmEUI4xM6JwxuAG4oQ\ni6WZlaExSdqMMGYxNS4hRH4Mv+MwzsycU8zujmxSydC4SnpWiwkhypmRSTpVPcexqqRLEJCFdFxm\n5WippIWoMOYl6cTVRQnx2UnZlBNjqSjc+LmcgpjRPEIIh5iXpGFaopmZlA3J0dOqfZOY0mcvhHCG\nmUk6yeoSPFMqaaZcgSKEEG4xMklP7UowtU861XCmdXdIJS1EZTEySU8tn6uqZj5lRhJKd3eYlaON\n6Q4SQjjDmCR9+KIWIDHuxalr5gNw1vGHsbgzcT/+Cas6AaitSYS8qKNh2uudzE2+mrkHiHrPsQsB\nOCP5vymWLUjc/3+y8pc4EiGEE2yN3ZGLfMfuiMZiNDbXMx4MATA0GqaloRaAkbEJGuq86a/ywfEJ\nfDXVeKsnP2PCE1Fi8Th1tbbHjMpoIhJNDE7jSywr0/34U2MshUxxDY+GaUwORFQKpo6tIHHlRuLK\nXUnG7iiW6qoqmhtq00l6avJrqq+ZNm9j3fS/AWptVL921XirqbHRMqVM0Nk0GxqXECJ3xnR3CCGE\nmE2StBBCGMzxPmkhhBDOkUpaCCEMJklaCCEMJklaCCEMJklaCCEMJklaCCEMJklaCCEMJklaCCEM\nZsxt4UqpnwCnkxhG+u+11i8Vef3XAO8j0SY/BD4KnAT0Jme5Vmv9gFLqk8DXgBhwg9b6JhdjOge4\nB3gjOek14BrgdqAa6AI+rbUOFTmuy4FPT5l0MvAy0AgEk9O+obXeoJT6J+BiEu/rlVrrB12K6Vjg\nXuAnWuufK6WWYrOdlFI1wK3AciAKfE5r/Y6Lcd0C1AATwKe01geUUhPAM1Ne+n4SRVSx4roVm/t7\nkdvrHiA1Wtg84HngKhLHwobk9IDW+mKlVCtwF9AKjAB/pbXucyiumfnhJYq0fxmRpJVSZwOrtdZn\nKKWOAm4Gziji+tcCxybX3wG8AjwG/LPW+g9T5msEvgecCoSBl5RSv3NqR8jgCa31x6fEcAvwn1rr\ne5RSVwGfV0rdVsy4kh8ANyXjORv4BHAMiZ3v9SmxHg78TxLvZSvwlFLqIa111Ml4ku/LfwCPTpn8\nr9hsJ+AjwIDW+pNKqQ+SOAgvcSmuH5A4eO9WSn0Z+Afgm8Cg1vqcGa//VBHjApv7O0VsL631xVOe\nvxm4cfKp6e1FIjk+rrW+Vin1BeBbyX+FxmWVHx6lSPuXKd0d7wd+D6C1fhNoV0q1FHH9T5Ko9gAG\nSFSEViM2nQa8pLUe1FqPkah8zixOiGnnAPclH98PnFviuL4H/FuG59YCf9Rah7XWAWAXcLQLMYSA\nDwP7p0w7B/vt9H7gd8l5H8G5trOK60vAb5OPA0BHltcXMy4rJrQXAEopBbRprV/M8vqpcaXecydY\n5YdzKNL+ZUqSXkhih00JJKcVhdY6qrVOfU2/HHiQxNeSryilHlNK/T+lVKdFnN3AIpfDO1opdZ9S\n6mml1AeARq11aMb6SxEXSqlTgD1a6wPJSf+qlHpSKXW9Uqq+WHFprSPJg2KqXNopPV1rHQPiSqmC\nhxK0iktrHdRaR5VS1cCXSXw9B6hTSt2llHpGKfUPyWlFiyvJ7v5e7LgA/p5ElZ2yUCn1G6XUs8ku\nBmbE69i+liE/FG3/MiVJz1SSgZCVUheReBO+QqK/6dta63XAJuAKi5e4HefbwJXARcBnSXQxTO2i\nyrT+YrXfX5PoawP4GfBPWuuzSPTHfbmEcdldb0naL5mgbwce01qnvtr/I/AF4IPAJ5VSJxc5rkL2\nd7fbqxZ4r9Z6fXJSL/C/gUtJnDv6N6XUzITseEwz8oOddTnSXqYk6f1Mr5wPI9EZXzRKqfOA7wAf\nSn5deVRrvSn59H3AuyziXMzcXxnzprXep7X+tdY6rrXeDhwg0RVUP2P9RY1rinOAZ5Ox/i4ZIyS+\n/hW9vWYYyaGd0tOTJ3k8Wuuwi7HdAryttb4yNUFr/Uut9UiyYnuUGe3ndlw57u/Fbq+zgXQ3h9Z6\nWGt9i9Z6QmvdQ+Kk9ZoZ8Tq6r83MDxRx/zIlST8MfBxAKfVuYL/Wumg/vZA8K3wtcGHqZJtS6rdK\nqZXJWc4BXgdeAE5RSrUppZpI9C095WJcn1RK/WPy8UJgAYkD/C+Ts/wl8Kdix5WM5zBgRGsdVkp5\nlFKPKKXakk+fQ6K9HgMuUErVJudfDGxxM64pHsF+Oz3MZJ/jR4D1uCT51Tystf7+lGkq2dXhUUp5\nk3G9UeS4ctnfixZX0inA5imxrlVK/Tj5uBE4Adg6I67Ue14wq/xAEfcvY4YqVUpdDaS/KmutN8/x\nEifX/QUSX++2Tpl8C4mvNaMkLuf5nNa6Wyn1ceCfSFxS9h9a6ztdjKuZRJ9lG1BLouvjFeA2oI7E\nibjPaa0nihlXMraTgB9orT+U/PsTJM6kB4F9wOVa61Gl1FeBTybj+u6Ur/dOx3IdsILEZW37kuu8\nFRvtlOx+uBFYTeLk1WVa6z0uxTUfGAeGkrNt0Vp/SSn1I2Adif3/Pq31/ylyXP8BfBsb+3uR4/of\nJPb7p7XWv07O502uX5E4wf8LrfUtycR4B4mTsQMkLm8cdCAuq/zw2WQMru9fxiRpIYQQs5nS3SGE\nEMKCJGkhhDCYJGkhhDCYJGkhhDCYJGkhhDCYJGkhhDCYJGkhhDDY/wf7Sm3O5SK5/QAAAABJRU5E\nrkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f234edab3c8>"
      ]
     },
     "metadata": {},
     "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": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7f23480810f0>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAD4CAYAAAATpHZ6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEf5JREFUeJzt3X2QXXddx/F3yA6Th93STb2wYakWRvxCZUag1FJC6KIJ\nj2UYDAXHWKuFsdTKkI7AdIQppMyIoLVI6SiRjK1lcMA4SDKUhAnyEFuaCSgPinylQCs01Sx0jRub\nLk0T/zhnO9tw7+7N3fvQ/eX9mtnpvb/z9J1vdz/39HfOPV124sQJJElleNygC5AkdY+hLkkFMdQl\nqSCGuiQVxFCXpIIMDfLgk5PTxd96Mzq6iqmpBwZdxkDZg4p9qNiHymL60GiMLGu1zDP1HhsaWj7o\nEgbOHlTsQ8U+VHrVB0NdkgpiqEtSQRacU4+INwCXzhl6HrAO+AvgBPCNzLyyXvdtwCX1+NbMvK3r\nFUuSWlow1DNzO7AdICIuAl4HfAB4S2YeiIiPRcTLgW8Dvw5cCDwB2BcRezLz4Z5VL0l6lFOdfrkW\neB/w1Mw8UI/tAjYALwY+k5k/ycxJ4B7g3K5VKklaUNu3NEbE+cAPgGPA1JxFh4C1wI+BySbj32y1\nz9HRVafFlfBGY2TQJQycPajYh4p9qPSiD6dyn/obgZubjLe6X7LlfZSzTod7VRuNESYnpwddxkDZ\ng4p9qNiHymL6MN+HwalMv0wAd1CdjZ81Z3wcOFj/jDUZlyT1SVuhHhFPBo7U8+UPAd+OiBfWi38N\n2A38I/DKiHh8vf448K1eFC1Jaq7d6Ze1VHPks7YAH46IxwH7M3MvQET8FfAlqlsar8zM490sVpIe\nK77wtXsXtf0lG5/RpUoera1Qz8yvAi+f8/5bwPom690I3Ni16iRJp8RvlEpSQQx1SSqIoS5JBTHU\nJakghrokFcRQl6SCGOqSVBBDXZIKYqhLUkEMdUkqiKEuSQUx1CWpIIa6JBXEUJekghjqklQQQ12S\nCmKoS1JBDHVJKoihLkkFMdQlqSCGuiQVZKidlSJiM/B24BhwLfAN4FZgOXAfcGlmztTrbQGOA9sy\nc3tPqpYkNbXgmXpEnAW8C3ghcDHwauA64KbMXA/cBVweEaupAn8DMAFcHRFrelS3JKmJds7UNwB7\nM3MamAZ+NyK+D7ypXr4LeCuQwIHMPAwQEbcD6+rlkqQ+aCfUzwFWRcROYBR4N7A6M2fq5YeAtcAY\nMDlnu9nxlkZHVzE0tPwUS156Go2RQZcwcPagYh8qJfRhZHjFovfRiz60E+rLgLOA1wA/B3y+Hpu7\nvNV285qaeqCNwy9tjcYIk5PTgy5joOxBxT5USunD9JEHF72PTvsw34dBO3e//DdwR2Yey8zvUk3B\nTEfEynr5OHCw/hmbs93suCSpT9oJ9c8CvxIRj6svmg4De4FN9fJNwG5gP3B+RJwZEcNU8+n7elCz\nJKmFBUM9M+8FdgB3Ap8B3kx1N8xlEbEPWAPckplHgWuAPVShv3X2oqkkqT/auk89Mz8MfPik4Y1N\n1ttB9QEgSRoAv1EqSQUx1CWpIIa6JBXEUJekghjqklQQQ12SCmKoS1JBDHVJKoihLkkFMdQlqSCG\nuiQVxFCXpIIY6pJUEENdkgpiqEtSQQx1SSqIoS5JBTHUJakghrokFcRQl6SCGOqSVJChhVaIiAng\n74B/q4e+CbwfuBVYDtwHXJqZMxGxGdgCHAe2Zeb2XhQtSWqu3TP1L2bmRP3zZuA64KbMXA/cBVwe\nEauBa4ENwARwdUSs6UXRkqTmOp1+mQB21q93UQX5BcCBzDycmUeB24F1i65QktS2BadfaudGxE5g\nDbAVWJ2ZM/WyQ8BaYAyYnLPN7HhLo6OrGBpafmoVL0GNxsigSxg4e1CxD5US+jAyvGLR++hFH9oJ\n9e9QBfkngKcBnz9pu2Uttms1/oipqQfaOPzS1miMMDk5PegyBsoeVOxDpZQ+TB95cNH76LQP830Y\nLBjqmXkv8PH67Xcj4r+A8yNiZT3NMg4crH/G5mw6DtzZUcWSpI4sOKceEZsj4q316zHgScBfA5vq\nVTYBu4H9VGF/ZkQMU82n7+tJ1ZKkptqZftkJfCwiXg08HrgS+BfgbyLiCuAe4JbMfCgirgH2ACeA\nrZl5uEd1S5KaaGf6ZRp4VZNFG5usuwPY0YW6JEkd8BulklQQQ12SCmKoS1JBDHVJKoihLkkFMdQl\nqSCGuiQVxFCXpIIY6pJUEENdkgpiqEtSQQx1SSqIoS5JBTHUJakghrokFcRQl6SCGOqSVBBDXZIK\nYqhLUkEMdUkqiKEuSQUZameliFgJ/CvwHuBzwK3AcuA+4NLMnImIzcAW4DiwLTO396ZkSVIr7Z6p\nvxO4v359HXBTZq4H7gIuj4jVwLXABmACuDoi1nS5VknSAhYM9Yh4BnAu8Ol6aALYWb/eRRXkFwAH\nMvNwZh4FbgfWdb1aSdK82pl+uR74feCy+v3qzJypXx8C1gJjwOScbWbH5zU6uoqhoeXtV7tENRoj\ngy5h4OxBxT5USujDyPCKRe+jF32YN9Qj4reAL2fm9yOi2SrLWmzaavxRpqYeaGe1Ja3RGGFycnrQ\nZQyUPajYh0opfZg+8uCi99FpH+b7MFjoTP2VwNMi4mLgKcAMcCQiVtbTLOPAwfpnbM5248CdHVUr\nSerYvKGema+ffR0R7wbuBl4AbAI+Wv9zN7Af+EhEnAkco5pP39KTiiVJLXVyn/q7gMsiYh+wBril\nPmu/BtgD7AW2Zubh7pUpSWpHW/epA2Tmu+e83dhk+Q5gRxdqkiR1yG+USlJBDHVJKoihLkkFMdQl\nqSCGuiQVxFCXpIIY6pJUEENdkgpiqEtSQQx1SSqIoS5JBTHUJakghrokFcRQl6SCGOqSVBBDXZIK\nYqhLUkEMdUkqiKEuSQUx1CWpIIa6JBVkaKEVImIVcDPwJGAF8B7g68CtwHLgPuDSzJyJiM3AFuA4\nsC0zt/eobklSE+2cqb8K+EpmXgS8Dvgz4DrgpsxcD9wFXB4Rq4FrgQ3ABHB1RKzpSdWSpKYWPFPP\nzI/PeXs28EOq0H5TPbYLeCuQwIHMPAwQEbcD6+rlkqQ+WDDUZ0XEHcBTgIuBvZk5Uy86BKwFxoDJ\nOZvMjrc0OrqKoaHlp1TwUtRojAy6hIGzBxX7UCmhDyPDKxa9j170oe1Qz8wXRMSzgY8Cy+YsWtZi\nk1bjj5iaeqDdwy9ZjcYIk5PTgy5joOxBxT5USunD9JEHF72PTvsw34fBgnPqEXFeRJwNkJlfo/og\nmI6IlfUq48DB+mdszqaz45KkPmnnQumLgD8AiIgnAcPAXmBTvXwTsBvYD5wfEWdGxDDVfPq+rlcs\nSWqpnVD/S+CJEbEP+DRwFfAu4LJ6bA1wS2YeBa4B9lCF/tbZi6aSpP5o5+6Xo8BvNFm0scm6O4Ad\nXahLktQBv1EqSQUx1CWpIIa6JBXEUJekghjqklQQQ12SCmKoS1JBDHVJKoihLkkFMdQlqSCGuiQV\nxFCXpIIY6pJUEENdkgpiqEtSQQx1SSqIoS5JBTHUJakghrokFcRQl6SCGOqSVJChdlaKiPcD6+v1\n3wscAG4FlgP3AZdm5kxEbAa2AMeBbZm5vSdVS5KaWvBMPSJeDDwrMy8EXgZ8ALgOuCkz1wN3AZdH\nxGrgWmADMAFcHRFrelW4JOmntTP98iXgkvr1/wCrqUJ7Zz22iyrILwAOZObhzDwK3A6s62q1kqR5\nLTj9kpkPA/9Xv30DcBvw0sycqccOAWuBMWByzqaz45KkPmlrTh0gIl5NFeovAb4zZ9GyFpu0Gn/E\n6OgqhoaWt1vCktVojAy6hIGzBxX7UCmhDyPDKxa9j170od0LpS8F3gG8LDMPR8SRiFhZT7OMAwfr\nn7E5m40Dd86336mpBzqreglpNEaYnJwedBkDZQ8q9qFSSh+mjzy46H102of5PgzauVD6BOBPgIsz\n8/56eC+wqX69CdgN7AfOj4gzI2KYaj59X0cVS5I60s6Z+uuBnwE+ERGzY5cBH4mIK4B7gFsy86GI\nuAbYA5wAtmbm4R7ULElqoZ0LpduAbU0WbWyy7g5gRxfqkiR1wG+USlJBDHVJKoihLkkFMdQlqSBt\nf/lIkkryha/dO+gSesIzdUkqiKEuSQUx1CWpIIa6JBXEUJekghjqklQQb2mU1JFObwkcGV7B9JEH\nmXj2+ECOXzpDXTpNGYplcvpFkgrimbqkgfC/FHrDM3VJKoihLkkFMdQlqSCGuiQVxAul0hLlhUY1\n45m6JBWkrTP1iHgW8Cnghsz8UEScDdwKLAfuAy7NzJmI2AxsAY4D2zJze4/qliQ1seCZekSsBm4E\nPjdn+DrgpsxcD9wFXF6vdy2wAZgAro6INV2vWJLUUjvTLzPAK4CDc8YmgJ31611UQX4BcCAzD2fm\nUeB2YF33SpUkLWTB6ZfMPAYci4i5w6szc6Z+fQhYC4wBk3PWmR2XJPVJN+5+WXaK448YHV3F0NDy\nLpTw2NZojAy6hIGzB5Vu9mFkeEXX9tVvS7n2burF30WnoX4kIlbW0yzjVFMzB6nO1meNA3fOt5Op\nqQc6PPzS0WiMMDk5PegyBsoeVE7uw+l6S+Lso3dFx38X830YdHpL415gU/16E7Ab2A+cHxFnRsQw\n1Xz6vg73L0nqwIJn6hFxHnA9cA7wUES8FtgM3BwRVwD3ALdk5kMRcQ2wBzgBbM3Mwz2rXBqgTs6y\nPUNVP7RzofSrVHe7nGxjk3V3ADsWX5YkqRN+o1SSCmKoS1JBDHVJKoihLkkF8dG7Oi2drveIq3ye\nqUtSQQx1SSqIoS5JBTHUJakgXijVkuSFTqk5Q10DYShLveH0iyQVxFCXpIIY6pJUEENdkgrihVJ1\n5FQudPo/h5D6x1BfohZ798jEs8e7VImkxxJD/TTlLYVSmZxTl6SCGOqSVBBDXZIK4pz6gDinLakX\nuh7qEXED8HzgBPCWzDzQ7WNIkprr6vRLRFwEPD0zLwTeAHywm/uXJM2v22fqvwr8A0Bm/ntEjEbE\nGZn5v10+zpKZvvCLN5L6qduhPgZ8dc77yXqsaag3GiPLOj3QJRuf0emmkvSY0GiMdH2fvb77pePQ\nliSdum6H+kGqM/NZTwbu6/IxJEktdDvUPwu8FiAingsczMzpLh9DktTCshMnTnR1hxHxx8CLgOPA\nVZn59a4eQJLUUtdDXZI0OD4mQJIKYqhLUkF89kuXzPd4hIh4MfBe4GEggTdm5vGBFNpj7TwmIiLe\nC1yYmRN9Lq9vFvh9OBv4W+DxwD9n5psGU2VvLdCDq4DfpPqb+EpmbhlMlb0XEc8CPgXckJkfOmnZ\nBuCPqPpwW2a+Z7HH80y9C9p4PMI24LWZuQ4YAV7W5xL7op3HRETEuVQX0ovVRh+uB67PzF8GHo6I\nn+13jb02Xw8i4gzgbcD6zHwhcG5EPH8wlfZWRKwGbgQ+12KVDwKbgHXAS+q/j0Ux1LvjUY9HAEbr\nX9xZ52XmD+vXk8BZfa6vXxbqA1SB9o5+F9ZnLfsQEY8D1gM76+VXZeZ/DqrQHprvd+En9c9wRAwB\nq4D7B1Jl780Ar6D6Ds+jRMTTgPsz8wf1f7nfRtW3RTHUu2OMKqxnzT4eAYDZZ99ExFrgJVT/8ko0\nbx8i4reBLwJ397Wq/puvDw1gGrghIv6pnooqUcseZOaDwFbge8A9wP7M/I++V9gHmXksM4+2WHxy\njw4Baxd7TEO9N37q8QgR8URgF/B7mfnj/pc0EI/0ISLWAL9DdaZ+ull20utx4M+Bi4DnRMQrB1JV\nf839XTgD+EPgF4CnAhdExC8NqrDHkK48VsVQ7455H49Q/xJ/BnhnZn62z7X103x9+BWqs9R9wCeB\n59YX0ko0Xx9+BNyTmd/NzIep5lp/sc/19cN8PXgm8L3M/FFm/oTqd+K8Ptf3WHByj8ZpMk1zqgz1\n7ljo8QjXU1353j2I4vqoZR8yc0dmnpuZzwdeQ3XXx9WDK7Wn5uvDMeB7EfH0et3zqO6IKs18fxN3\nA8+MiJX1++cB3+l7hQOWmXcDZ0TEOfW1hYup+rYofqO0S05+PALwHOAwsAeYAr48Z/WPZea2vhfZ\nB636kJmfnLPOOcDNhd/S2LIPEfHzwM1UJ1XfBK4s8RbXBXpwBdV03DHgjsx8++Aq7Z2IOI/qpO4c\n4CHgXqqL5N+v+/Ai4H316n+fmX+62GMa6pJUEKdfJKkghrokFcRQl6SCGOqSVBBDXZIKYqhLUkEM\ndUkqyP8DELoGG86E/1sAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f2353ca90b8>"
      ]
     },
     "metadata": {},
     "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": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.80329863296935067"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "accept.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Find the index of all diverging transitions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([], dtype=int64),)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trace['diverging'].nonzero()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "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": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f234168ada0>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8nHd96PvP8zyzaaSRNJJGu+Xdv8TO4iQksUkggaRA\nS+6hbOWc9kBT0styyTlQbhcoPbdwUuiltE0J5dLQjbYcKHApKUsasoCzkDiLE+/2T7YlWfuu0Tqa\n7XnOHzNj7ZZsS5pF3/frpVjzbPOVMvOdn77PbzEcx0EIIURhMbMdgBBCiNUnyV0IIQqQJHchhChA\nktyFEKIASXIXQogC5Mp2ABkDA+Nr1m0nGPQzMjK1VpdfExLz+snHuCXm9ZEPMYdCAWOx7Rui5e5y\nWdkO4ZJJzOsnH+OWmNdHPsacsSGSuxBCbDSS3IUQogBJchdCiAIkyV0IIQqQJHchhChAktyFEKIA\nSXIXQogCJMldCCEKUM6MUBVCbAwHDnet6vXu3NuwqtcrFJLcxYb0XNfBC98HxnyMj08DcHvDvmyF\nJMSqkuQuhChoyWSSP/uzz9Pd3UUikeC3f/sj/OM//i0333wrr776CuFwmC9+8UFqa2t5+OGvcvTo\nYWw7ybve9Wv8+q+/l89//rO4XG7GxsJ8+tN/zB/90e8TjUbZv/82fvSjR/jkJ/+AJ598jP/xPx4A\n4Itf/BNuu+0N3H77HVn9uaXmLoQoaE888RiVlVV85SsP86d/+hc89NBfAFBcXMyXv/w19u17Pc88\n8zOOHHmNvr5evvrVv+XLX/4b/umf/oHp6dRfdKWlpXz+81/iscd+zJYt2/ja1/6ekpIAjuNwyy37\nOHHiBNFoFNu2OXbsCLfe+vps/siAtNyFEAXu+PGjHDnyGkePHgYgGo1SVJTg+utvAKC6uprR0VGO\nHTvCiRPHuP/+DwHgODYDAwMA7N69B4C2tjZuuOEmAG6//Y1861v/jGVZ3Hbb7Rw8+AsqK6u47rq9\nuN3u9f4xF5DkLoQoaC6Xmw984IP80i+97cK2++//EJY1M+Oj4zi43W7uuecdvP/9v3VheygUuHCN\n9JGYZmqGXcOYmWn3bW97O9/85j9RV1c/53myScoyQoiCtnv3NTz33NMAjIwM8/DDX13yuF/84lls\n2yYajfLgg3+24Jj6+kZOnz4FwMGDz1/YvnOnYnBwgFOnTrB3741r8FNcOmm5CyHW1Xp3XXzzm+/m\n1Vdf5iMf+SDJZJIPfvBDHDt2ZMFx1157PTfccBMf/vBvAQ7vfOd7FxzzK7/yf/DpT3+S++//EDff\nfCumOdM+vvnmW5mamprTos8mw3GWXwBJKfUgsA9wgI9rrV+etc8HPAzs0Vq/bt55RcBx4AGt9Tcu\n9hxruRJTKBRgYGB8rS6/JiTmtTWnK2Qg/7pC5tPvOqMQYu7t7eH8+TZuvXU/x48f5e///mEefPCr\nOI7DJz7xMX7v9z5NY+Om9Y5x0U+TZVvuSqk7gJ1a6/1KqauBfwD2zzrkS8BhYM8ip/8RMHzp4Qoh\nRO4pLi7hO9/5X3zjG3+L48AnPvG79PR085nP/D5vfvPd657YL2YlZZm7gEcAtNanlFJBpVSp1nos\nvf8PgUrgN2afpJS6CtgN/GQV4xVCiKwJBAL85V/+9YLt//AP38xCNBe3kuReCxya9XggvW0MQGs9\nrpSqXOS8vwDuB35zJYEEg/41Xa8wc9c7n0jMaycw5pv7OJB6nC/xQ37FmiExr5/LuaG67N0CpdQH\ngBe01q1KqRVddC1XGC+EWl8+yKeYMzV2mFtzz5f48+l3nSExr42lPnxWkty7SbXUM+qBnmXOeTuw\nTSl1D9AIRJVSnVrrJ1fwfEIIIa7QSpL748DngIeVUjcC3Vrri36Uaa3fl/leKfVZoE0SuxBCrJ9l\nk7vW+nml1CGl1POADXxMKXUvMKq1/oFS6nvAJkAppQ4AX9daf2stgxZC5K/Z3VBXQ750X53tvvve\nz5/8yRd58snHueGGG7nqqt189KP3sXnzFu6//3e4//4Pcfvtb+QjH7n/sp9jRTV3rfWn5m06Mmvf\nwp7+c8/97KWHJcSVWWzOcJn3W+Sa97//XgB6e3uJx+P80R99jsOHX2XTpk1XlNhBRqgKIQrcYlP+\n3nTTzRdGmV5s2t977/1Nbr31jlWZ9vev/upLHD9+jKamzSQScQA+//nPcuedd/Hooz+kq6uTL3zh\nc5w9e4b+/l7+5m/++ooSvMwtIzasocgwT3U8y4HWF5iIT2Y7HLFGlpryF5af9vdrX/sa0eiVT/vb\n2trCsWNH+frXv8GHP/wx2tvPz4nx/vt/h6amzfzhH/4x99//CfbuvVFa7kLApS3dFkvGeKXvMGfC\nLQD0Tw1wdqiNXcEdvL7+FkxD2jyFZLEpf+PxVMt5uWl/bdtmcHAQuLJpf9vaWti9+xpM06Smppb6\n+rUvEUpyFxvO95p/yJlwCwFPCa+r3ovhtnmx8zCnR87wUu+r7Kt73fIXEXljsSl/M5ab9nd2P/cr\nmfbXcbhwDqQ+NNaaNFFEwUraC+eiOz18hud7XqLcW8Yvb76L2uJqdlRu4a5Nb8QyLB459yiRRCQL\n0Yq1stIpfzPHzp7294EHHlhwzOVM+9vUtBmtT+M4Dr29PfT0dK/Gj3ZR0nIXBWVsMsYhPcDgaIRI\nNEljqJhrtlVQHfQTTcb49unvY2Bwa+2NWOZMq63Y7Wd3peLY4En+o/Up3rXzniz+FIVtvbsuLjbl\n71LmT/v7gQ+8f8ExlzPt744dO9m2bTsf/vBvsWlTEzt37lq1n28pK5rydz3IlL9zScyX5sDhLoZG\np3nqUCfTsSR+nwuPyyQ8EQPgJhWi7trzPNX+DL/UdCdVRRUXzs1MP5C0kzzV8Swj0TCfueWT1BZX\nZ+VnWQl5fayPxWLOtWl/L3vKXyHyQXgiyuMvdRBP2uzbXcOupnIcx6F/JMLPj7Rz6Fwn/tLnKPb4\nKfeWLXoNy7R4146387fH/4Un25/mv1590SEcYoPKl2l/JbmLvOc4Di+e6COetLn9ujq21ZcCqZtd\nNRV+rr0+zvHe8ziGTQ3bcZlLzz56XWgPQW85r/Yf4b273oHX8qzXjyHyRL5M+ys3VEXee+FEL30j\nETZVl1xI7LN5ihK4atpx4h6aX67hxLkwzR2pr/lMw2Rf3U1EkzEO9x9bj/CFWBOS3EVem5pO8N2f\nncUyDW6+evEa+WCiB9tIUBJvJB5zcbbZxcVuNWW6Qr7Q8/LSBwmR46QsI/LO7AFLJ1qHGZuKs3dH\nJSVF7gXH2k6S/kQ7BibbK+qJldkMDZoMDxlUVi2e4auKKtlZvo0z4RYGI0NUFS22Fo0QuU1a7iJv\n2Y6Dbg9jmQaqKbjoMd3xFmJOlCpXHW7Tw45dCcChrcWFc5FxJJnW+4s9h5Y+SIgcJsld5K3O/gkm\nInG2N5Ti9Sx+k7Q1egKAkCs13NtfDDV1NpEpg77epV/+e0PX4rE8vNj7KrnSXViISyHJXeStU20j\nAFy1efFW+1RyjN74eYrNUorMkgvbm7YkMU2H820W8fjiidvn8nJt5dUMTQ/TNbHcwmNC5B6puYu8\nFJ6I0jcSoa7ST3mJd9Fj2mKpIeJVrvo5271eaGi06Wi3aD6bYN8tM/tmLySR6QZ5bPAkjYG51xAi\n10nLXeSl1u4xAHY2Lj4gyXFsWqMncOEmaC3sRVPfmGq9Hz+ZIJlcvPVeV1yLicHRwROrF7gQ60SS\nu8g7juPQ2jOOyzJorC5Z9Ji+RAdT9jibvLuwjIV/oLo9UFtvMznlcLo5tug1PJaban+I9vEuRqYX\n9okXIpdJchd5Z3B0molInKaaAC5r8ZdwW/QkAFu9e5a8TkNjEtOEV16LYC8ygyRAQ0kdkCrNCJFP\nVlRzV0o9COwDHODjWuuXZ+3zAQ8De7TWr5u1/c+AN6Sf40+11v+2moGLjStTktlat3A0asv0MRJO\nnM7YWXyGn3B8YMEMfRleH2zfanHmXJLObpumxoU9bhpK6jjUf4Sjgyd5Y+PrF7mKELlp2Za7UuoO\nYKfWej9wH/DQvEO+BByed86bgGvS57wN+KvVCVdsdEnbpq13HK/boq7Sv+gx4eQADjYVrtolE3vG\nnqtTA59Onk6tzJOZliDz1dUbY1NJPc0j54gkplf3hxFiDa2kLHMX8AiA1voUEFRKzW4y/SHwg3nn\nPANkptQLA8VKqaVnaxJihc52jjIdS7K5tmTOyjazDSV6Aai0ape9XkXQpKHeRXePTXh08VFN14b2\nkHSSnBpuvvzAhVhnKynL1AKzh+kNpLeNAWitx5VSc8Zna62TQGbF4fuAR9PblhQM+nG51i7/h0KB\nNbv2WpGYF/rhC6mFhXdtriBQ4luw30nGmZgKU+qqIOBfPpZAwMf115p0dU9wtsUh1LhwCoPbt+/l\n0dYnaJls4a2h2678h1gl8vpYH/kYM1xeP/eL/507i1LqHaSS+1uWO3ZkZOoyQlmZQlkkINetR8wv\nHO3GZRmU+V2MTywsk/ROdQIQNKuJTseXvd74+DTbt3op9hucPB3ldeVxXPPeFYFkkBJ3Ma92Hae/\nf2zZUs96kNfH+siHmJf68FlJWaabVEs9ox5YdsieUuqtwGeAX9Zaj67geYS4qJ6hSfpGItRXFWOZ\nC1+6juMwlOjFwFy0b/tSTNNA7XSRSMBg/8LrPt/9EpW+CkZjY/y45XGe6zo4Z7CTELloJcn9ceA9\nAEqpG4FurfVFP8qUUmWkbrTeo7UevuIohQAOnx0EYNMSfdtHkn1EnQjlVtWifdsvZud2C8OA3p7F\n3xL1JTUA9Ez2XtJ1hciWZd8BWuvnlVKHlFLPAzbwMaXUvcCo1voHSqnvAZsApZQ6AHwdKAGqgO8q\npTKX+oDWun0NfgaxQRw+M4gBNISKF93fHtUAVLiWv5E6X3GxSWO9SUcXTEwYlJTM7fde688k9z52\nV6rFLiFETllR80Zr/al5m47M2rfUQpNfv9yghJhvIhLnbNco2xvK8HkWvmwdx6Yj1oyFi1KzYpEr\nLG/XThcdXTH6ekxKds69/+9zeanwBRmIDBFPxnFbC2+8CpFLZISqyAunz4/gOHDNtsUT90Cii2ln\niqCrGtO4vJf1pgYLt8ehv8/EXqRXZH1xDQ4OfVMDl3V9IdaTJHeRF060pW7d7NmyeHJvj6VLMlbN\nZT+HaRpU19gkEwbDQwt7xNQVp67dLXV3kQdkyl+R0zJL6h3SA7hdJuf7xhcMXko6CTpjZykySigx\ny6/o+ULVNl0dFgN9FlWhxJx9Fb4KPKabnsk+WcBD5DxpuYucNz4VYyISp67Sv+io1L54O3Enyibv\nzivug15c4lDkdxgeMkjMze2YhkFtcQ1TiQhjsdzu+yyEJHeR87oHUwPclppLJlOSafJceS8Ww4Dq\nmiSOYzA0sPDtMVOa6bvi5xJiLUlZRuS8nqHUTBZ1lQu7QCacGN2xFkrMcsqtakYS/Zd07eaOMF6f\ne85o1lC1zflW6O8zqambe2c1k9x7JLmLHCctd5HTbMehd2iKkiI3Af/C7oddsRaSJGjyqFWbFsBX\nBIFSm9GwQTQ6d1+Ry0fQW8ZAZJDpRHTxCwiRAyS5i5w2Mh4llrCprfAvmrw7MiUZ7+oOLKqusQFj\n0ekI6oprsR2bM+Fzq/qcQqwmSe4ip/UNp+rtNRVFC/ZF7Qi98XbKrWoCVnBVn7cqZGMYqT7v82VK\nMyeH9Ko+pxCrSZK7yGn9IxEAaoILb6Z2xs7iYK/KjdT53B4oDzpMTpgL5nmvKqrAbbo5MaSlS6TI\nWZLcRc5yHIe+4QjFPhcli9TbMyWZTd6da/L8qdIMnGud2yfSNExq/SGGpofpjwyuyXMLcaUkuYuc\n1T04STSepKZiYau9OfIqA4kuis1SemNttEwfo2X62Ko+f0WVjWk6tLQmF7TQ64pTk5NJaUbkKknu\nImc1d4QBqAkurLePJocAKLOq1uz5LQsqq2zGJxwGhxbvEinJXeQqSe4iZ+lMcl+k5T6aTJVDytcw\nuQOE0qWZlta5s0T63UXUF9dyJnyOWHL5FZ+EWG+S3EVOchwH3RHG57EW9G+3nSRjyWE8hg+fsfjc\n7qulPOjg8UDL+SS2Pbc0s6fyKuJ2gjPhljWNQYjLISNURc7ITBIGqfnbRydiNNWULOjfPpjoJkmC\nCqt2zdczNU3Y0mTRfDZJX79NXe3MIu67K3fxRPsBTg6dZo8s4CFyjLTcRU4aCKe6QIbKF9bbu2Ot\nAJRblesSy/atqTbQ/F4z28q24LU8nByWurvIPZLcRU4aDE8DUFXum7PdcRx64i2YWJSYqztwaSk1\n1SZFRQZt7UmSyZnSjMt0oYI76Z8aZDAytC6xCLFSktxFThoIRzAMqCydm9wn7DAT9iilVsVlr7h0\nqUzTYNtmi1gMunrm9prJrKd6QnrNiByzopq7UupBYB/gAB/XWr88a58PeBjYo7V+3UrOEeJikrbN\n8FiUioAXlzU3gffHOwAotS5vndTLtW2rxYnTCVpaEzQ1puruz3UdZCo+deF7yzC5vWHfusYlxFKW\nbfoope4Admqt9wP3AQ/NO+RLwOFLPEeIJQ2PRbEdh6pF6u0DidRN18A6lWQyqipNAiUG7Z1JEomZ\n0kyx20+pJ0Df1ABJO3mRKwixvlbyd+1dwCMAWutTQFApVTpr/x8CP7jEc4RYUqbeHlqk3t4f78Rn\nFOM1Fib+tWQYBlu3WCQS0N45N4nXFdeQdJIMSN1d5JCVlGVqgUOzHg+kt40BaK3HlVLzuy1c9JzF\nBIN+XC5rqd1XLBQKrNm118pGizlQkkrm4ckYAJvrywmUeDk9kfrDMGJPEnWmqPLU4ivyXHmws3h9\nC+euuRBXIBXXtXtcHD0+Rkenw/XXznzwbLc3oUfOMhgfXNf/Zxvt9ZEt+RgzXF4/98vpWLzsOSMj\nU5dx2ZUJhQIMDOTXmpcbMebxiVSLvWdwEq/bwnRsxiemL6ySNBQfAMDvlM1ZOelKzV+Jab5DJ3sv\nfO8vdtHSFmNwKILXk3pZl1CKZVicH+lat/9nG/H1kQ35EPNSHz4rKct0k2p1Z9QDPWtwjhBEogkm\nInGqyn0LBiiN2yMABKzybIQGpJbgcxyD8+0zpRnLtKjxhxiLjTM8PZK12ISYbSXJ/XHgPQBKqRuB\nbq31ch9ll3OOEAyOZurtc2vqjuMwYYdx4cFrLL5Q9noIVafnmmmbO6ApM5GYdIkUuWLZ5K61fh44\npJR6nlSvl48ppe5VSr0TQCn1PeBfU9+qA0qpX1/snLX7EUQhyYxMrSqbezM16kSIOzECVvmaTzlw\nMb4iCARsenptIpGZXjP16SmAT0lyFzliRTV3rfWn5m06Mmvfe1d4jhDLujAydV5yv1CSWecukIup\nqrEZHzd58fA49Y0zg5oC7hJOjZwhnozjtpa+QSvEepARqiJn2I7D4GiEshIPHvfcnlMTyVEASrJY\nb88IhWzAYWDe4tkNJXXEkjH0yNnsBCbELJLcRc4YnYiSSDqEyhb2YZ+yxzCx8GWx3p7h8UJZucP4\nmMl0ZGZ7Y0kdAEcHT2QpMiFmSHIXOWNgicnCkk6CaWcKvxnIar19tswiHrNb75VFlZS4izk2eArb\nsZc6VYh1Icld5IylpvmdslMdrYrN3BnkXFVlYxhzSzOmYXBt1W7GYuOcH+vIYnRCSHIXOWQwPI3L\nMigrmTv6NJPc/WbujBR0uSFY4TA1aTI5MfPXxPWhPQAcGZDSjMguSe4iJ0xNxxmdjFFVVoQ5r/Qy\neaHlnjvJHRYvzajgTjymm6ODJ7MVlhCAJHeRI1p7Ugl8fr0dUjdTLVx41nmysOVUVNqYpsPggImT\n7vLusdxcXanom+qnb7I/uwGKDU2Su8gJ57pTXR3n19tjdpSoE8mpm6kZlpVK8NMRg6nJmdiuq9oN\nIK13kVWyQLbICS3dqQlD5w9eCidTrd9cupk6W2XIYXAABgdMuDq1aEc0GcUg9X2RK/XzyCIeYr1J\ny11kneM4tHSPUVLkpsg7t70xnOgDcutm6mzBilRpZmhgpuXutbyEiqoYnB4mkpjOYnRiI5PkLrKu\nPxxJzQRZtrDePpJuuedqcne5oDzoMDVlEh6d6dvekB7Q1DUhk6GK7JDkLrKupStVkplfbwcYSfTh\nwo3HWJj4c0VVKJXUXzk6QXNHmOaOMPHxEkCSu8geSe4i6y7U2+f1lInaESbtsZy8mTpbRWVqrpnh\noZm3k9csosxTSu9UP3E7sfTJQqwRSe4i6851j+KyDCpKvXO2z5RkcvNmaobLDaVlDhPjBrHYzPbG\nQD22Y9M72Ze94MSGJcldZFUsnqSjf4KmmgCWOfflOJK+mZprg5cWk2q9G4wMz/wMmYnEOie6sxSV\n2MgkuYusau+bIGk7bKtb2DofSeRHyx0yyR1GhmbKR0FvOX5XEd0TvSTt5FKnCrEmJLmLrMoMXtrW\nsDCBDyf78Bl+3IZnwb5cU+QHr89hZNjETneaMQwjNce7HefcaGt2AxQbjiR3kVXn0jdTt9eXzdk+\nbU8SsScIumpy+mZqhmGkWu/JpMHY6Ey8mS6RMpGYWG+S3EVWtXaPEvC7F/Rxz5Rkgq7qbIR1WSoq\n0qWZWXX3an8Id3oiMcdxljpViFW3oukHlFIPAvsAB/i41vrlWfvuBr4AJIFHtdYPKKVKgH8GgoAX\n+JzW+qerHbzIbyPjUYbGouzdUbWgdZ7pKVNh1RCxJ7IR3iUrLXcwDIfwsAHbU9ssw6S+uIbz4510\nTvSwKVCf3SDFhrFsy10pdQewU2u9H7gPeGjeIQ8B7wZuA96ilNoN3AtorfWbgPcAX17NoEVhyPRv\n31q/SL093VMmn1rulpXqEjk5ac7tElmSSuiy/J5YTyspy9wFPAKgtT4FBJVSpQBKqW3AsNa6Q2tt\nA4+mjx8EKtPnB9OPhZijpSd1M3X7Isl9JNFPkVmCzyxe77CuSDBdmhkdmXlr1RXXYBkWR6XuLtbR\nSsoytcChWY8H0tvG0v8OzNrXD2zXWn9FKXWvUuosqeT+9uWeJBj043JZyx122UKh3O8rPV+hx9ze\nP4lhwM3X1uP3uQmUpOruU4lxpp1Jmnw7CZT48OJeq3Av8PpW5zlCNdDWAmNjFoFA5j6Cj2tqFEd6\nT2IUx6nyV6zOcxX46yNX5GPMcHlT/l6s64IBoJT6r0C71vptSqnrgb8HXnexi46MTF1GKCsTCgUY\nGBhfs+uvhUKPOZ6w0edH2BQqYXJ8msnxacYnUjModsdS64+WUsX4xDTR6fiaxQypxL5az+HxgMvt\nZmgQxsYiF+4l7Azs4EjvSX5x5jVeX3/LFT9Pob8+ckU+xLzUh89KyjLdpFroGfVAzxL7GtLbbgN+\nCqC1PgLUK6XWrlku8k5rzxiJpM3OTeUXtrVMH7vwBTBtT134Pl8YBpSX28SiBqNjM71jIvHU4t/P\ndL3Ac10Hea7rYLZCFBvESpL746RuiqKUuhHo1lqPA2it24BSpdQWpZQLuCd9/Fng1vQ5m4EJrbUM\n0RMXnOkMA6BmJfeMXF0zdaXKK1JJvbtn5iUf8JRQ7PLTN9mPLV0ixTpYNrlrrZ8HDimlnifVM+Zj\n6Xr6O9OHfBT4NvAs8B2tdTPwMLBFKfU08C3gI2sSvchbzR2pm6k75yV3x3GYssfwGD5ceTAydTHl\nwdRN1a6emfndDcOgtriamB1neHokW6GJDWRFNXet9afmbToya98zwP55x08Av3bF0YmCZNsOZ7vC\n1ASLKCuem8DjTpQEccrNhS36fOHzQVGRQ29fEtt2MM1U3b2uuIZzo230TPZSVbQ6N1WFWIqMUBXr\nrqN/gkg0ya5FSzKpvu/5WpLJKA/axOPQPzjTeq/xV2Ng0DPZn8XIxEYhC2SLdXHgcNeF70+1pcoS\niyX3qXS9PR9mgryY8gqbnm6L7p4ktdWpvgQey01lUQVDkSFiydgyVxDiykjLXay7vnS314u13HN1\nzdSVKit3MAzonlV3B6jzV+MAvVPSehdrS5K7WFeO49A/EsHvdS2YLCx1M3Ucr1GEy1j7gUtryeWC\nUJXJwKBNLDbTO6a2uAaAvsmBpU4VYlVIchframwyznQsSXWwaMFkYTFnmiSJvG+1Z9TXmTgO9PTO\ndIms8JXjMl30RSS5i7UlyV2sq/50SaamomjBvpmSTH7X2zPq61K19tldIk3DpLqokvHYBOHoaLZC\nExuAJHexrvpGUiM1q4P+Bfum8nzw0nzVVSZu99zBTJCa4x3gzEhLNsISG4Qkd7Gu+kcieNwm5SUL\nByhNFcjN1AzTNKitsRgbd5iYmGm9Z5J788i5bIUmNgBJ7mLdTEbiTETiVAf9C+rtSTvJpD2OzyjG\nMgqnh25DXeotNrs0E/SW4zbdnAlLchdrR5K7WDeZkkxNcGG9vWuiB5skJWbZgn35LFN3n12aMQ2D\n6qJKBiJDjEyHsxWaKHCF00QSOS9zM7U6ndxnD2w6M30YgBKrsJJ7WalBsd+guzeJ4zgX/mKp9ofo\nmuyleeQct9bdlOUoRSGSlrtYN/0jEVyWQWWpb8G+wURqFulCa7kbhkF9nUk0CkPDM/3dazI3VcNy\nU1WsDUnuYl1E40nCEzGqyoouTKSV4TgOg/FuXHjwGAtLNvlusdJMubcMv6tIbqqKNSPJXayLgQtd\nIBcm7yl7jGlnkhKrbMGN1kJQX5vp7z6T3A3DYHv5Voamh6XuLtaEJHexLvouktwLtSSTUVRkUBE0\n6O+3SSRmSjM7yrcCcC7cmq3QRAGT5C7WRf9IBAMIlS+W3LsBKMnjOdyX01BnkbShr3+mS+TO8m0A\nnBmV5C5WnyR3sebiiSRDo9MES724XQtfckOJbixc+M2SLES3PmamIpgpzTSW1OOxPJyVlrtYA5Lc\nxZpr7RnHdpxFSzIxO8pocogKVw2GUbgvx5pqE8uce1PVMi22l22hd7KP8dhEFqMThahw300iZ2QW\nw65ZZD6ZgUQnAFWuhnWNab25XAY11SbDIw6RSKru/lzXQdxmaqjJj1t+ynNdB7MZoigwKxrEpJR6\nENgHOMBKt83mAAAfaElEQVTHtdYvz9p3N/AFIAk8qrV+IL39N4DfBxLA/6O1/skqxy7yxJnO1OyH\ni9Xbe+PnAah1NzGWHF7XuNZbQ71Fd69NZ3eSndtTb71QURUAA5EhNgUK+wNOrK9lW+5KqTuAnVrr\n/cB9wEPzDnkIeDdwG/AWpdRupVQl8MfA7cA9wDtWNWqRN2zH4WznKAG/G79vblvCcRx64+dxGx4q\nXHVZinD9NG1K1d3b2mdKM5W+IKZh0j81mK2wRIFaSVnmLuARAK31KSColCoFUEptA4a11h1aaxt4\nNH383cCTWutxrXWP1vpDaxO+yHXdA5NMRRNUL9Jqn7DDTNljVLuaMAu43p5RVmoSDBp0dScvrM5k\nmRaVviDhaJhYMp7lCEUhWUlZphY4NOvxQHrbWPrf2UvK9APbAT/gV0r9EAgCn9VaP3WxJwkG/bhc\n1iWEfmlCofybRrYQYn75TKpF2lRXSqBk7rQDHeOpLpBbAjsIlPjwkr2l9by+1X/uQGDhNAtqp8PB\nlyIMDJlctcsLQGN5LQORISaN8Uv6f14Ir498kI8xw+VNHHaxIYTGrH8rgXcCm4GfK6U2a62dpU4c\nSU8qtRZCoQADA+Nrdv21UCgxv3aqD4DSIhfjE9Nz9rVNnAWgPFnP+MQ00enstFy9PveaPPf4+PSC\nbfW1qbfAKT1NQ13q+zIr1b+/fah7xf/PC+X1kevyIealPnxWkty7SbXQM+qBniX2NaS3TQLPa60T\nwDml1DgQItWyFxvImc4wJUVuSovnLs6RdBIMxDsptSrwW/nZMlpOc8fi0wr4/S46Oh1OtkZwuWBr\nQwUGBv2RoXWOUBSylRQ6HwfeA6CUuhHo1lqPA2it24BSpdQWpZSL1M3Tx9Nfb1ZKmembqyWA3DHa\nIA4c7uLA4S5+/EIbQ2NRggHvgjljBhPdJElQ696cnSCzqDLk4DgGw4Opt5/bdBP0lTMcGZa6u1g1\nyyZ3rfXzwCGl1POkesZ8TCl1r1LqnelDPgp8G3gW+I7Wullr3QX8/8BB4D+A/5a+4So2kItNFtYT\nawPYkMm9uibVW6avd+btV11UhY1D21h7tsISBWZFNXet9afmbToya98zwP5FznkYePiKohN5rT+c\nTu7zeso4jkNX/Cwuw1Pwg5cWU+SH0jKb0bDJdOpXRMhfyemRM5wNt7AruD27AYqCUPj9z0TWDI1O\nYxhQUeqds30k2ceUPU69e1tBrZd6KWpqU3/I9vWmeoiFiioBZJ4ZsWokuYs1kbQdhsdT9XbLmvsy\n64ylesls8uzMRmg5oSpkY5oO/b0mjuPgtbyUeUppGT1Pwk5kOzxRACS5izURHo9i286CJfUcx6Ez\ndgYXHmrcTVmKLvssF1RV20SjBl09qVZ8tb+KuB2nY7xrmbOFWJ4kd7EmhkZTfbyryuYm93Cyn0l7\njHrP1g1bksmob0gl9ZOnUj1kMvPMSGlGrAZJ7mJNDI6lknvlvOSeKck0buCSTEZJwKG0zKaz2yY8\nalPtT9Xdm8Oyrqq4chu76STWzNDoNJZpMOI6Qzg9UNNxHFqjJzCxiCTHaZk+lt0gc0B9Y5KxUZOT\npxO8/tYiavwhzoVbSdpJLHPtpuMQhU9a7mLVJZI24YkoFaVeZs8HNmWPE3UilFlVmIYkLoDKSoeS\nYoOz5xJEow67gjuIJmOcH+/Mdmgiz0lyF6tueCyK40BV2dz+7cPJXgAqXTXZCCsnGSZcfZWLRBL0\n2cSFPu7NI2ezHJnId5LcxarL3EydXW93HJvhRB8u3JSaFdkKLSft2uHC5YJTpxNsK90KgB6Ruru4\nMpLcxaobHE0Nu5zdU2bMHiFBnGCBr5V6Obweg53bXUxOOTS3TtFQUkfLaBtxmWdGXAF5l4lVNzQ6\njdtlEvDPzJE+nEiXZCwpySxm91Wpvg1PvNKBCu4gYSdoHTuf5ahEPpPkLlbV1HScsak4lWW+CzNB\nJp0E4eQAXqMIv1ma5QhzT3NHmL7RMYKVNue6xujvTC0kLqUZcSUkuYtV1dabWtigatbI1HByEBub\nCqtmwdS/YkZDY2q2yN72IkzDlJuq4opIcherqrVnDJh7M3UokVrbpcJVu+g5IqWs3CEY8NLZO029\nv562sQ4iiYWrOQmxEpLcxapq60m13DPJPWpHGLdHKDHL8Jn+bIaW8wwDdm8J4jhgToawHVta7+Ky\nSXIXq6q1dwyfx6LYl7pBONO3vS6bYeWNLXWlFHkt2s+kPhxPDuksRyTylSR3sWpGxqcZHoteuJnq\nOA6DiR5MTIJWdbbDywuWaXBVU5DISAA3Xk4ON+M4S64rL8SSJLmLVXMmvSB0pn/7YKKbmDNNuVW9\n4WeAvBQ7N5VjmRaJcAXD0yP86NXj2Q5J5CFJ7mLVnGlPJfdMvb0tehKAKinJrFjL9DG67ZOEahNE\nh1OzRJ6YOpjlqEQ+WlFzSin1ILAPcICPa61fnrXvbuALQBJ4VGv9wKx9RcBx4AGt9TdWMW6Rg850\njACplnvCidERO4PH8FFilmc5svxT35ik97XU/O5jyaEsRyPy0bItd6XUHcBOrfV+4D7goXmHPAS8\nG7gNeItSavesfX8EDK9SrCKHOY7DmY4wlaU+fB4XnbGzJIlT6aqTvu2Xwe+HYKkHeyrAeDJMLBnL\ndkgiz6ykLHMX8AiA1voUEFRKlQIopbYBw1rrDq21DTyaPh6l1FXAbuAnaxG4yC1DY9OMTcbYWhcA\nZkoylZb0bb9cDQ1JkuEqHMOmWUariku0kuReCwzMejyQ3rbYvn4gU2D9C+CTVxqgyA+Z/u1b60qZ\nSI4ykOgi5GrEaxYtc6ZYSlnQwR1NlWYO9Z7IcjQi31xOF4aL/Y1tACilPgC8oLVuVUqt6KLBoB+X\na+0WcAiFAmt27bWSTzH3vdgOwN6ra/jxmccAuKr0epIkshnWinl97uUPyoJNoXLaE26ODpygsuo3\nMWfNqJlPr48MiXn9rCS5dzPTUgeoB3qW2NeQ3vZ2YJtS6h6gEYgqpTq11k8u9SQjI1OXEvclCYUC\nDAyMr9n110K+xXzi3CCGAaVek+aJY7hwU+Vspj16OtuhLcvrcxOdzs3pdSsroa07xHRlNwebT7Cz\nYguQf68PkJjXylIfPispyzwOvAdAKXUj0K21HgfQWrcBpUqpLUopF3AP8LjW+n1a65u11vuAvyPV\nW2bJxC7ym+04nO8bpyFUQvtkG1P2OI2enbiM3GwN5xPLgjIjBMDjM53UhFjWsslda/08cEgp9Typ\nnjEfU0rdq5R6Z/qQjwLfBp4FvqO1bl6zaEVO6hueIhJNsnNTOU93/QKA7b7rshxV4dhUUY6TNNFj\np2W0qlixFdXctdafmrfpyKx9zwD7L3LuZy8rMpE3MjNB1tYZ/Hv/aSqsWipkndRV4/dblIYbGPd0\n8FLrOW7dtiPbIYk8ICNUxRU7151K7oPWaRwcdviuz3JEhWdfY+p3+thpKc2IlZHkLq5YS9cYLpfD\na4OHCLhLaPRIy3K13b3rJnAMepPn6A9Hsh2OyAOS3MUVicaTdPRPULV1iMn4FLc13CqThK2BEk8x\n9d4mzOIxHnnpWLbDEXlAkru4Iud7x7Edm1j5GSzT4vb6W7MdUsG6a9s+AF4bOEx4PJrlaESukyaW\nuGTPdc3MUnjsVBwz2EfEGOVNm19P0FcOTGYvuAJ14HAXCacCw3FBsIsH//UQ9//qtdkOS+QwabmL\nK3KufRp3fQsAFdN7OHC4K8sRFS6X4aHRsx3TF+FYz1nGp2QyMbE0abmLKzLhDGMVj1FuVdMTbc/Z\nkZ75rmU6VWf3mumFx8s7+ckL5/nPd+3MYlQil0nLXVy2yUkbQqlWe517c5aj2RhKzSAuPFiVvfzs\ntXaGRqezHZLIUZLcxWU72zeIVTqCJ1aJ38zPyZXyjWGYVLhqMFxx7EAfP3i2JdshiRwlyV1cttap\n1EwT1aa02tdTyFUPQMmmTp4/3supNlkPRywkyV1clpHpUSLufuzxckIlsozeevKZxZS7Kon7BjH9\nY/zTY5pYPJntsESOkeQuLsvxAQ2AZ2wrpryK1l2trwmApj1D9IcjPPJsa5YjErlG3pbiko3HJuic\n6sSeDBD0VGQ7nA0p6A5RVVTJkHWOUKXJYy+1c/jsYLbDEjlEkru4ZCeHU632RM82yqUikxWGYXBH\nw34SdoKb9k/jdpn83Y9OMiDzzog0Se7ikgxFhmkdbceIFsNoLSUBmV88W/bV3YzP8vHqyEv857u3\nMhVN8JXvHyUSzY+lDcXakuQuLslPz/8cB4do53aqQ5bU27PI7y7iTZtuYzw+QbKijTff2EDnwCSf\n/5dX+NmrnRw43CUjhjcweWuKFRueHuFgzyv4jGKSQ3XU1azdguZiZd686Q0UuXw8cf4A77qziYZQ\nMd2DU7x4sk9WbdrgJLmLFXv8/AGSThLf2A7AoLFBkns2Pdd1kFf7j7K9bAsT8Um+efq7NKkRggEv\nZzpHOdEq/d83MknuYkVGpsO80P0SVb5KRlpqKCoyqKwwsh2WAFRwB27TzanhMzhWnLtuasDvdfFq\n8yBtvePZDk9kyYomDlNKPQjsAxzg41rPLMOulLob+AKQBB7VWj+Q3v5nwBvSz/GnWut/W+XYxTp6\nov0ACSfJjeX7+fdogl07TAxDknsu8Fgedlfs4sjgCXrirewqvZE339TAYy+289zRHoq9LkLBogvH\n37m3IYvRivWybMtdKXUHsFNrvR+4D3ho3iEPAe8GbgPeopTarZR6E3BN+py3AX+1umGL9RSOjvKL\n7peo9FUQH6gDYJOUZLKuuSN84cucDOE1iuhPdDGaGKKi1McdextwHIefvdol0wNvQCspy9wFPAKg\ntT4FBJVSpQBKqW3AsNa6Q2ttA4+mj38GeG/6/DBQrJSSbJCnnjz/NAk7wVu3vInjLSOYJtTXyf/O\nXGIaJo2enYDD4amncRyHhlAxt15dQzSe5MBr3SSTdrbDFOtoJcm9FhiY9XggvW2xff1AndY6qbXO\nLMdzH6lyjUx+kYfC0VGe6z5IhS/IruJraO+boLbGxO2WkkyuKTMrKTUr6U900BU/C8CupnJ2NpYx\nMh7l5dP9WY5QrKfLWazjYu/qOfuUUu8gldzfstxFg0E/LtfatQZDofybkjYXYv7Owe8TtxP82rX3\n0NwxAcCunT4CgdSiEV6fe87x8x/ni3yMe7GYt3uu5ujoC7w2dYCtZdvxWX7efHMTw+NnaO4Ypamu\nLKuvq1x4TV+qfIwZVpbcu5lpqQPUAz1L7GtIb0Mp9VbgM8DbtNajyz3JyMjUSuK9LKFQgIGB/Oo1\nkAsxt46288z5Fwl6yxgZneBHzw1iGFBf4zA+nlokYvbKS16fOy9XYsrHuJeK2cTDnqL9HI08x7MD\nP2VfyS8D8Ibr6vjJ8+c5cKiDO6+ro76qeL1DzonX9KXKh5iX+vBZSVnmceA9AEqpG4FurfU4gNa6\nDShVSm1RSrmAe4DHlVJlwJeAe7TW0tk2DzmOw/fP/BCAG6uvZ2zMYWjYobHexOeTkkwu2+W7gQqr\nlo5YM52xVHmmtNjD/mtrSSQdvvbIcaIxqZIWumWTu9b6eeCQUup5Uj1jPqaUulcp9c70IR8Fvg08\nC3xHa90MvA+oAr6rlDqQ/mpamx9BrIWXel+ldaydTYEGqv1VnGtJzVeyfassu5vrDMPk5pJfwsTi\n0ORTTCbHANhSG+CqpnK6Bif55uM6y1GKtbaid6rW+lPzNh2Zte8ZYP+8478OfP2KoxNZMRgZ5rvN\n/47H8rA3dA2O43CuNYnLBU2bpJdMPii1Ktjrv4NXp37GCxM/4U2l78UyXNx0VYhoPMkvjveyq6mc\nN1xXn+1QxRqREapijqSd5BsnvsV0cpr37fpVStzFdPfaTEw6bGmycLmkJJMvtnmvYYvnakaS/bw2\n9TQAlmnykXdcg9/r4n893kznwESWoxRrRZK7mOPHrY/TOtbOzTU3cGvtTQCcPJW6cXe1cs0ZONPc\nEc5mqGIJLdPHaJk+Rmv0OEGrmiKjhNbocU5HUgPLQ+VFfPDtVxNL2Px/PzjOdEymCC5EktzFBQc6\nf8Hj539Ola+C96l38vSRbo7oMB1dSQKlNiOR3O41IBYyDYsd3uvwmwGORZ7n7HSqonrjrhBvuXkT\nvcNT/PNjWmaQLECS3AWQmmHwe83/TsBTwv91/QcpcqX6sXd3WYBBfaOMbsxXHtPHVs8eXHh4beoA\n3zjxLZ7rOkj9rkG215dy8GQfT7zSme0wxSqT5L7BOY7DE+cP8K/6B5S4i/nvez9ETXE1ANFYkv5e\nE4/HobJKkns+85l+dvn24sLNy32HOTJwAsOAj/7qNZQVe/jOz85wRNZgLSiS3DewaDLGP574Fo+c\ne5RST4Db62+lZbSN57oO8lzXQX5xVpNMGjRsSsqKSwWgyCxB+W6ixF3MyWHNCz0vU1Js8t/fcx0u\ny+RvfniCtt6xbIcpVom8ZTeozvFuvvTKVzjUf4RtZZv5g5s/TtA3s9r1VMShuzPVaq+tl1Z7ofCZ\nft7SdCdVvgrOj3fy54e+SnFZlP/znt3EYkn+8jtH6JIeNAVBkvsGYzs2T7Y/zZde+Qo9k33c0fh6\nPn7Dhynzzh3CfPR4HNs22LQ5iSVd2wvK+Z4Im4zrCLka6Jro4fMvfpmTo8fZd00tE5E4X/jmIX7w\nbEu2wxRXSIYbbiC9k/1889T3aB07T8BTwvuv/jX2VF614LjRMZvTzQm8PoeaOmm1FyLTMGnyKLZ6\nr+HQ5FO8OPkYTeWKm66+gUOnwvz0xXb2bKlg16by5S8mcpIk9w0gaSd5qv0ZftL2BAk7QVOgkZuq\nr2dkOsxzXQfnHOs4Ds+/GMO2Yev2hNTaC1zSiXOV73W0xU7SHtO4S1vZvXc/p474+PN/Pcx/uXsn\nd+6tl1W38pAk9wL3o3OP8VLvqwxHw/gsL/vqb2JTYOll1lrbkvT02jTWm1RWSd/njcBn+lHeG+lN\nnKcnfp5Wz9PUvq6BYb2Vf/mp5pkj3dxydTW/fOvmbIcqLoEk9wI1Mh3mJ61PcLDnZRxga2kTN1Rf\nh9fyLHlOZNrhxVdiWBbsu8VDb3h6/QIWWWUYJnXurQStGgYT3fQl2jGv6qJkoo6O9ia6n53AZZnc\nubcBt0v+nMsHktwLzFQ8whPtB/h5x7PE7QRlnlJuqL6WuuKai57nOA7PPR8lMg033+imNGDSK7ML\nbDg+088bAr9Kf6KdE5EXGSrpwbu7B2e6mO8db+HRlzfz1ht38fprayn1L91QENknyb1APN3xC86E\nWzgxpInZMfyuIm6svo4tpZsxV1AvPdWcoKPLpq7W5Jrd8rLYyFqjxwHY7L6KSquWgUQ3Yd8A7iZN\n1GnmkY7D/NvhTTQEatgeqqO+shjTnHmN3bl36bKfWD/yLs5ztmPzSt9hftL6BJOJKdymm+ur9rAr\nuAOXubI+jD29SV56JY7XC3fc5pGbZwIAwzAIWEECVpCEE8cVGONsuI3Rin6sin56oz66ehpxnW5i\nc2U1W+oCVJcXZTtskSbJPY+dGmrmkXOP0jnRjWmYXBXcye7KXXgt74qvER61eerpKI7jsPOqBJ1D\nMRhaw6BFXnIZbnYFt7MruJ3h6RHOhltpG+vEbDyL03CW1nCIM6c2443V0D8S4Zara9hSG5CGQhZJ\ncs9DPzz3HxweOEHfVGo1+y2lTVxXtZtit/+SrhMetfnpk1FiMdh5VZLyoPSOEcur8AW5pTbIDdXX\n8fL5ZgbiXUwFB7CCA9iREp5s3cxPX64nVFbMLVfXcMvVNTSG1n/N1o3OyJWpPgcGxtcskHxY5Ha+\nxWIejAzzo5bHeKXvMAC1/mr2hq6ZM23ASg0OJXn8qSjT0dQNVG/p5BXHnI8LTUN+xp1rMU8mx+hP\ndDCc7AccLMdLoq+R6Z5NEPdRV+nn9r0NbK8NsL2hFCtPBlDkQ+4IhQKL/nkkLfc8EI6O8uT5p3mm\n6wWSTpKgt5y9oWuoTc/eeCkcx+Hk6QQvvxrHtuH1t3q4apeL5o41CFxsGMVWKVutPTTYO6iv8fJc\n90Ema8/hr20lEGtisK2a7z01CRj4vS72bK1ANZWzo6GMxlDJnBuyYnVIcs9RjuPQNtbOs50Hebnv\nNZJOkkpfBf9p+9uIxCOXVcvs6Uty6LU4/QM2Pi+84TYvmxpk4hixejyml9Lxa3lrydWcj57mTPQw\nY542XLvaKKYEz0QDE72VvHw6zsunU2VFr8diW10pjaES6qr81FcWU1vhJ+B3S83+CqyoLKOUehDY\nBzjAx7XWL8/adzfwBSAJPKq1fmC5cxaz0csySTvJ8HSY9vEOWsfaOT58ioHJ1J3NGn+Iuza9kVvr\nbsJluhZMGbAUx3EYn3Do6ExyrjXJ4FBqnpjKKpvtuxJ4Vrmbcq6VClYqH+POl5gdx2HCDjOc6GXE\nHiDppJb085geip0qArEmhvt8DPS6wJ7b1jQNg4pSL+UlXgJ+Nz6PC6/Hwus28botvB4Ln9uitXcM\nt8vCbZm4XCYel8mdNzRQ4nNf8V8E+ZA7lirLLJvclVJ3AL+ntb5HKXU18A9a6/2z9p8E3gp0AU8D\nHwZCFztnMZeb3G3HZjI+RdJJ4jgOtmNjOw62k8Qm9bi8vIjB4XGSTpJ4MkHCSRC3EyTsmX8zX0kn\nmfnVzPovGBgXHhjM/C7ntywMDGzHJm7HiduJ1L/J+IXnidnx9PPGiSdT/07GJxmLTeAw8ysocvkI\n+ULUFzVS7U0NQHKcWV+AbUMi4ZBIQDzuEI3B5KTDxKTN+IRDeNQhEnHScUJjvcne69xrtlxeviSc\n+fIx7nyM2e01GZzqJ5wYYNQeIu7E5uz3UoIr6Ye4j2TcJJEwScQN4jEDbDP9/nPASL9PDGfmcXqf\ncWFbar/bZeC2LDyWC48r9eV1u/C53fjcLtymC5dl4TItXKYLt2XhtlyYhomFRSDgJzIRT1/fSF07\nkwGc9LZZecE0DCzTwDItLNPCbZrpf11YloXLNLFME5dpUVLkwcDAY3kuOnJ8OVdSc78LeARAa31K\nKRVUSpVqrceUUtuAYa11B4BS6tH08aGlzrnsn2AJ/3jiW7zaf3S1L7suLMPEMizclpuqogr8Lj9B\nXzmVvnLqyuv4538ZQ8cALm8aAL/fYMtmC9MTpSpk4/HASESmFBDZYRoWZVYVZVYVjuMQdSJM2GGm\n7Aki9gRRJ8Kk1Q8W4Js5z30Fz+kAsfTXxTdmj8t08Zlbfodqf2h1r7uCY2qBQ7MeD6S3jaX/HZi1\nrx/YDlRd5JxFLfXps5xPvemjl3NaXvhPf5rtCIQQ+epy+iNdLAkvtU/uigghxDpaScu9m1SrO6Me\n6FliX0N6W+wi5wghhFhjK2m5Pw68B0ApdSPQrbUeB9BatwGlSqktSikXcE/6+CXPEUIIsfZW2hXy\n/wXeCNjAx4AbgFGt9Q+UUm8Evpg+9Pta6z9f7Byt9ZE1iF8IIcQicmb6ASGEEKsnPyZ4EEIIcUkk\nuQshRAEq6Lll0qNrvwd8UGv94/S2A0AxkJkG8f/WWh9a/Arrb4mYrwe+Rmr4xVGtdU527ldK3Qs8\nAJxLb3pCa/357EV0cZc6RUa2KaXuJPXaOJHedExr/d+yF9HFKaWuAf4deFBr/ddKqU3Av5AaptQD\nvF9rHc1mjPMtEvM3gJuYWeXgS1rrn2QrvktRsMldKbUd+CTwi0V2/5bW+vg6h7Ssi8T8V6STj1Lq\nW0qpX9Za/8f6R7gi39Fa/262g1hO+kN0p9Z6f2aKDOCiU2TkiKe11u/JdhDLUUoVA18Bnpq1+X8C\nX9Vaf08p9QXgg6QaLTlhiZgBPp1paOWTQi7L9ADvAkazHcglWBCzUsoDbJ3VqvwRcHcWYis0c6bV\nAIJKqdLshlRQosCvkBr3knEn8MP097n4Ol4s5rxVsC13rfUUgFJqsd3/UylVBZwCPqG1jqxnbEtZ\nIuYqYGTW436gbh3DulR3KKUeIzUlyO9qrV/LdkBLuNi0Grlst1Lqh0AF8Dmt9RPZDmgxWusEkJj3\nWi6eVYbJudfxEjED3K+U+iSpmO/XWg+ue3CXoSCSu1Lqt4Hfnrf5j7XWP13k8C+TqlufU0p9jVS/\n/T9f6xjnu8SYZ8uJqRyWiP/bwGe11j9RSu0H/hm4dt2Duzw58Xtdxhngc8B3gW3Az5VSO7TWOTIF\n1iXJh983pO4RDGmtDyulPgV8Frg/uyGtTEEkd6313wF/t8JjfzDr4Y+A961JUMvHsdKYB4DKWY8z\nUzxk1XLxa61fUEqFlFKW1jq51HFZdLFpNXKS1roL+E764TmlVC+p10Nr9qK6JBNKqaL0X8o58Tpe\njtZ6dv39h+TQPYLlFHLNfQGllKGUelIplVl09E4g526szqa1jgOnlVK3pze9C3gsiyEtSSn1+0qp\n/5L+/hpgIEcTO+ThFBlKqd9QSv1u+vtaoIbUOgr54kng3env302Ovo5nU0p9Pz21OeRBvpitYEeo\nKqXeDvwecBWp1m+P1votSqlfA/6AVFfILuC+TK072y4S827gYVIfxi9qrT+ZxTCXpJRqJPVnrEnq\nr8Lf0Vq/lN2olpZvU2QopQLAt4BywEOq5v5odqNanFLqJuAvgC1AnNR77TeAb5Carf08qV5rObPi\nyBIxfwX4FDAFTJCKuT9bMV6Kgk3uQgixkW2osowQQmwUktyFEKIASXIXQogCJMldCCEKkCR3IYQo\nQJLchRCiAElyF0KIAvS/AR6EmHEcWB+2AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f23416fef28>"
      ]
     },
     "metadata": {},
     "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": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "If the overall distribution of energy levels has longer tails, the efficiency of the sampler will deteriorate quickly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "model = pm.Model()\n",
    "with model:\n",
    "    mu1 = pm.Bernoulli(\"mu1\", p=0.8)\n",
    "    mu2 = pm.Normal(\"mu2\", mu=0, sd=1, shape=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 10000/10000 [00:02<00:00, 3695.88it/s]\n"
     ]
    }
   ],
   "source": [
    "with model:\n",
    "    step1 = pm.BinaryMetropolis([mu1])\n",
    "    step2 = pm.Metropolis([mu2])\n",
    "    trace = pm.sample(10000, init=None, step=[step1, step2], njobs=2, tune=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'accept', 'p_jump', 'tune'}"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trace.stat_names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Both samplers export `accept`, so we get one acceptance probability for each sampler:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  1.00000000e+00,   2.99250475e-05],\n",
       "       [  2.50000000e-01,   2.58944314e-04],\n",
       "       [  1.00000000e+00,   5.89596665e-03],\n",
       "       ..., \n",
       "       [  2.50000000e-01,   7.71917301e-02],\n",
       "       [  1.00000000e+00,   2.43450517e-01],\n",
       "       [  1.00000000e+00,   2.40661294e-02]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trace.get_sampler_stats('accept')"
   ]
  }
 ],
 "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": 2
}