docs/source/notebooks/howto_debugging.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How to debug a model\n",
    "\n",
    "There are various levels on which to debug a model. One of the simplest is to just print out the values that different variables are taking on.\n",
    "\n",
    "Because `PyMC3` uses `Theano` expressions to build the model, and not functions, there is no way to place a `print` statement into a likelihood function. Instead, you can use the `Theano` `Print` operatator. For more information, see:  theano Print operator for this before: http://deeplearning.net/software/theano/tutorial/debug_faq.html#how-do-i-print-an-intermediate-value-in-a-function.\n",
    "\n",
    "Let's build a simple model with just two parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "CompoundStep\n",
      ">Metropolis: [sd]\n",
      ">Metropolis: [mu]\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='24000' class='' max='24000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [24000/24000 00:14<00:00 Sampling 4 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 4 chains for 1_000 tune and 5_000 draw iterations (4_000 + 20_000 draws total) took 15 seconds.\n",
      "/env/miniconda3/lib/python3.7/site-packages/xarray/core/nputils.py:215: RuntimeWarning: All-NaN slice encountered\n",
      "  result = getattr(npmodule, name)(values, axis=axis, **kwargs)\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import pymc3 as pm\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import pandas as pd\n",
    "import theano.tensor as tt\n",
    "\n",
    "x = np.random.randn(100)\n",
    "\n",
    "with pm.Model() as model:\n",
    "    mu = pm.Normal('mu', mu=0, sigma=1)\n",
    "    sd = pm.Normal('sd', mu=0, sigma=1)\n",
    "\n",
    "    obs = pm.Normal('obs', mu=mu, sigma=sd, observed=x)\n",
    "    step = pm.Metropolis()\n",
    "    trace = pm.sample(5000, step)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 0., 0., ..., 0., 0., 0.])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trace['mu']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hm, looks like something has gone wrong, but what? Let's look at the values getting proposed using the `Print` operator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mu __str__ = 0.0\n",
      "sd __str__ = 0.0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Only 3 samples in chain.\n",
      "Sequential sampling (1 chains in 1 job)\n",
      "CompoundStep\n",
      ">Metropolis: [sd]\n",
      ">Metropolis: [mu]\n",
      "Sampling 1 chain for 0 tune and 3 draw iterations (0 + 3 draws total) took 0 seconds.\n",
      "/dependencies/arviz/arviz/data/base.py:146: UserWarning: More chains (3) than draws (2). Passed array should have shape (chains, draws, *shape)\n",
      "  UserWarning,\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sd __str__ = 1.753782660690967\n",
      "mu __str__ = 0.0\n",
      "sd __str__ = 0.0\n",
      "sd __str__ = 0.0\n",
      "mu __str__ = -1.179953715933347\n",
      "mu __str__ = 0.0\n",
      "sd __str__ = -0.7921743241812005\n",
      "mu __str__ = 0.0\n",
      "sd __str__ = 0.0\n",
      "sd __str__ = 0.0\n",
      "mu __str__ = -1.6164899771368246\n",
      "mu __str__ = 0.0\n",
      "sd __str__ = 1.7649365130019956\n",
      "mu __str__ = 0.0\n",
      "sd __str__ = 0.0\n",
      "sd __str__ = 0.0\n",
      "mu __str__ = 0.5633904409221188\n",
      "mu __str__ = 0.0\n",
      "sd __str__ = 0.0\n",
      "mu __str__ = 0.0\n",
      "sd __str__ = 0.0\n",
      "mu __str__ = 0.0\n",
      "sd __str__ = 0.0\n",
      "mu __str__ = 0.0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/dependencies/pymc3/pymc3/sampling.py:603: UserWarning: The number of samples is too small to check convergence reliably.\n",
      "  warnings.warn(\"The number of samples is too small to check convergence reliably.\")\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as model:\n",
    "    mu = pm.Normal('mu', mu=0, sigma=1)\n",
    "    sd = pm.Normal('sd', mu=0, sigma=1)\n",
    "    \n",
    "    mu_print = tt.printing.Print('mu')(mu)\n",
    "    sd_print = tt.printing.Print('sd')(sd)\n",
    "    \n",
    "    obs = pm.Normal('obs', mu=mu_print, sigma=sd_print, observed=x)\n",
    "    step = pm.Metropolis()\n",
    "    trace = pm.sample(3, step, tune=0, chains=1, progressbar=False) # Make sure not to draw too many samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the code above, we set the `tune=0, chains=1, progressbar=False` in the `pm.sample`, this is done so that the output is cleaner.\n",
    "\n",
    "Looks like `sd` is always `0` which will cause the logp to go to `-inf`. Of course, we should not have used a prior that has negative mass for `sd` but instead something like a `HalfNormal`.\n",
    "\n",
    "We can also redirect the output to a string buffer and access the proposed values later on (thanks to [Lindley Lentati](https://github.com/LindleyLentati) for providing this example):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Only 5 samples in chain.\n",
      "Sequential sampling (1 chains in 1 job)\n",
      "CompoundStep\n",
      ">Metropolis: [sd]\n",
      ">Metropolis: [mu]\n",
      "Sampling 1 chain for 0 tune and 5 draw iterations (0 + 5 draws total) took 0 seconds.\n",
      "/dependencies/arviz/arviz/data/base.py:146: UserWarning: More chains (5) than draws (2). Passed array should have shape (chains, draws, *shape)\n",
      "  UserWarning,\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeQAAAEGCAYAAAC0IuZwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOy9e3xdZZXw/12535ukTShtAuVSLi0ilwriyEVQB5lRQJSLzoi+Og4qo+NlHPw5L4Oo7wjeZ3SG11FeL8MICOKAIOooiqOgtFBKQ6G0XEyaQtLm5J6cXM76/bGfJ9k9PUlOkrNP9mnX9/PJJ/s8+9nPXnv3ss66PGuJqmIYhmEYxtJStNQCGIZhGIZhCtkwDMMwYoEpZMMwDMOIAaaQDcMwDCMGmEI2DMMwjBhQstQCGIXLihUrdM2aNUsthmEYRsGwadOmParalOmcKWRjwaxZs4aNGzcutRiGYRgFg4i8MNM5c1kbhmEYRgwwhWwYhmEYMcAUsmEYhmHEAFPIhmEYhhEDTCEbhmEYRgzIm0IWkfNF5GkR2SEi12Q4/04R6RaRze7nPW78JBF5SETaRGSLiFwWukZE5LMisl1EtonIB93434XW2SoikyLS6M7dLCJdIrI17f6fdutvFpGficgqN35haHyjiLw6dM2NTq5tIvLPIiJu/DJ3TZuI3BCa/+WQXNtFpNeNHy4ij7rxNhG5KnTNFSLyhFvvfhFZ4cbf6uamRGRDaP7rRGSTu2aTiJwbOlcmIt9w935KRC5x44eJyAMi8pi7zwUL+1M2DMMwFoyqRv4DFAM7gSOBMuBxYF3anHcCX8tw7THAWne8CtgN1LvP7wK+CxS5z80Zrn8j8MvQ57OAU4CtafPqQscfBG5yxzWAuOMTgafc8auA37pnKwYeAs4BlgN/BJrcvO8A52WQ62+Am91xGVAeut/z7llLgC5ghTt3I3CdOz4eOBb4FbAhtO7JwCp3fAKwK3TuU8Bn3HFRaN1vAO9zx+uA57P5cz311FPVMAzDyB5go87wf2q+LOTTgB2q+qyqjgG3Ahdmc6GqblfVZ9xxJ4GC8puq3wdcr6opd74rwxJXAN8Prfcg0JPhPv2hj9WAuvFB9xL3GXe/K3DKFCgFXiL40vGMqna7ef8NXDKbXKo6pqpJN17OtOdC3E+1s77rgE53zTZVfTrDczzm3hNAG1ApIuXu8/8C/snNS6nqntCz1LnjZf4ehmEUBr95ppvn9gwttRjGIsmXQl4NtIc+d7ixdC5xLtM7RKQ1/aSInEagAHe6oaOAy5wr+ScisjZtfhVwPnBnNkI693c78Hbg2tD4xSLyFHAvgVJDVR8CHiCw2HcDP1XVbcAO4FgRWSMiJcBFQGvafQ4HjgB+GRprFZEtBO/pBlXtVNVxgi8dTxAoyXXAt7J5FsclwKOqmhSRejf2aece/4GIHOLGrgP+QkQ6gPsIrPeZ3tF73fve2N3dPdM0wzDyyIdu3cxNv9o590Qj1sQpqeseYI2qngj8nMDVO4WIHAp8D3iXt4gJrMlRVd0A/Dtwc9qabwR+q6r7WcSZUNVPqmorcAtwdWj8LlU9jkC5ftrJczSB27iF4MvFuSJypqomCJTobcBvCNzPk2m3uhy4Q1WnxlW13T370cCVInKIiJS6tU4mcGFvAT6RzbOIyHrgBuCv3VCJk/V3qnoKgYv9C+7cFcC3VbUFuAD4nohk/Luhqt9Q1Q2quqGpKWP1N8Mw8shkSkkMj7F3KDn3ZCPW5Esh72JfK7HFjU2hqntDbttvAqf6cyJSR2CdflJVHw5d1gH80B3fRRDjDXM5IXf1PLiFDG5m5+4+0iVWXQw87Fzag8BPgDPcvHtU9XRVPQN4GtierVzO3bwVOBM4yY3tdG7z2wli17MiIi0E7+Mdquq/Nu8Fhpl+Xz8giKUDvNut7S3/CmDFXPcxDGPp6RsZRxUSw+NLLYqxSPKlkB8B1orIESJSRqCQ7g5PcBaw503ANjdeRqBcvquqd6St+yPgNe74bEKKT0SWubH/ykbANHf3hcBTbvzoUPb0KQRW+V6CxK2zRaTEWbJnh2Rudr8bgPcTfMHw9zkOaCCwUP1Yi4hUhq55NYEi3wWsExFvir7O32OW56gn+PJyjar+1o87hX4PQeIZwHnAk+74j+4zInI8gUI2f7RhFAA9Q2MAJNxvo3DJS3MJVZ0QkauBnxJkJN+sqm0icj1BxtndwAdF5E3ABEHS1Tvd5ZcSZEYvFxE/9k5V3Qx8DrhFRD4MDALvCd32YuBnqrpPpoOIfJ9AKa1wMdN/VNVvAZ8TkWOBFPAC4LceXQK8Q0TGgRHgMlVVEbkDOJcgvqvA/ap6j7vmqyLycnd8vaqGLeTLgVtDiWIQuL6/KCJKkMT1BVV9wsn7KeBBd/8X/HsRkYuBfyFIcLtXRDar6p8SuNqPBq4VER8Hf71LePt7Anf0VwgU7rvc+Y8C/+7eo7r3G5bPMIyYkhgOFHHPsCnkQkfs/11joWzYsEGt25NhLC0/a3uR935vEyKw47MXUFwkSy2SMQsissnlPe1HnJK6DMMwjHniLWTVIJ5sFC6mkA3DMAqYnqHx0LG5rQsZU8iGYRgFTCIUO05YHLmgMYVsGIZRwISzqy3TurAxhWwYhlHAJIbHqK8qnTo2ChdTyIZhGAVMz9AYRzXVuGNL6ipkTCEbhmEUML3D4xy6rILykiJ6zUIuaEwhG4ZhFDA9w2M0VpfRWF1mWdYFjilkwzCMAmUypfSNjNNQVUZDVZnFkAscU8iGYRgFim8s0VBVSkN1qTWYKHBMIRuGYRQo3kXdUO0sZHNZFzSmkA3DMAoU76KeiiGby7qgMYVsGIZRoExZyC6G3DcyzmTKGgYVKqaQDcMwChS/zSlwWZdag4kCxxSyYRhGgeILgTRWldFQXebGzG1dqJhCNgzDKFASw2NUlBZRWVZMo1PItvWpcDGFbBiGUaAkhsZoqAoUsf9tmdaFiylkwzCMAiUxHFLIZiEXPKaQDcMwCpSeobEpV3VjlY8hW1JXoWIK2TAMo0BJDI9PtV6sLCumvKTILOQCxhSyYRhGgZIYnraQISgQYjHkwsUUsmEYRgEyMZmaaizhsQYThY0pZMMwjALEN5ZIt5BtH3LhYgrZMAyjAPGdnXwM2R/3WsengsUUsmEYRgESbizhsQYThY0pZMMwjAIk3FjC4xtMTEymlkosYxGYQjYMwyhAEqFeyB5rMFHYmEI2DMMoQHwMuTFsIU9V6zKFXIiYQjYMwyhAwo0lPNZgorAxhWwYhlGA9AyN7WMdw3Q82bY+FSamkA3DMAqQ3uEx6tMVsrOQe81CLkjyopBF5DgReUhEkiLysVnmHSEivxeRHSJym4iUufEvi8hm97NdRHrTrqsTkQ4R+Zr7XBuav1lE9ojIV9KuuUREVEQ2uM+lIvIdEXlCRLaJyCdCc28WkS4R2Zq2xudF5CkR2SIid4lIvRtfIyIjofvfFLrmMje/TURuCI2Xu2fe4d7BmtC5E937a3PyVbjxz4pIu4gMpsl1uIj8wt3nVyLSMtv7cmNXuLW3iMj9IrJipj8nwzCWnnBjCY81mChs8mUh9wAfBL4wx7wbgC+r6tFAAng3gKp+WFVPUtWTgH8Bfph23aeBB/0HVR3w8901L4SvEZFa4EPA70NrvBUoV9WXAacCfx1Sit8Gzs8g78+BE1T1RGA78InQuZ0hGa5y910OfB44T1XXAytF5Dw3/91Awj37l927QERKgP8ArnLXnAP4f233AKdlkOsLwHedXNcD/zTb+3L3+CrwGnfNFuDqDOsahhETEsPj+2RYQ9BgoqLUGkwUKnlRyKrapaqPMK1I9kNEBDgXuMMNfQe4KMPUK4Dvh647FTgE+NkM6x4DNAO/CQ1/mkDhjYbFBKqdcqoExoB+J/+DBF8q0p/rZ6o64T4+DLSkz0njSOAZVe12n/8buMQdX0jwzBC8g/PcO3k9sEVVH3f33Kuqk+74YVXdneE+64BfuuMH3NrAjO9L3E+1u2cd0DnHsxjGAcftG9t5x81/WGoxsqJnaIyGUJUuT0OVlc8sVOIUQ14O9IYUXAewOjxBRA4HjsApGxEpAr4IzOgGBy4HblNVddecArSq6r1p8+4AhoDdwB+BL6jqfkp4Fv4X8JPQ5yNE5DER+bWInOnGdgDHOpd2CcEXjlZ3bjXQDuDeQR/BOzkGUBH5qYg8KiIfz0KWx4E3u+OLgVoRWT7T+1LVceB9wBMEingd8K1MC4vIe0Vko4hs7O7uzjTFMAqWh3fu5cHt3YxNxLuwxsRkiv7RfRtLeBqqyiyGXKDESSFnw+XAHd5CBN4P3KeqHXNc832YUuBfAj6aYd5pwCSwikDpf1REjsxGKBH5JDAB3OKGdgOHqerJwEeA/xSROlVNECi+2wgs9ufdPWejBHg18Hb3++KQm3smPgacLSKPAWcDu9x9Mr4vESl1cp1M8Pxb2Nf9PoWqfkNVN6jqhqampjnEMIzComsgCcCeweQSSzI7mRpLeKzBROFSEtXCIvIB4K/cxwtUdS4X6F6gXkRKnIXYQqBIwlwOfCD0+QzgTBF5P1ADlInIoKpe42R4OVCiqpvc/FrgBOBXgWeWlcDdIvIm4G3A/c5a7BKR3wIbgGfneM53An9OEBdWAFVNAkl3vElEdhJYuhtV9R6C2C8i8l6mFfIuAmu5w1nPy9w76QAeVNU97pr7gFOAX8wkk3vXb3bza4BLVLVXRDK+L+BOd91Od83twDWzPbdhHIh0DYy630lW1VcusTQz42PE6TFkP7ardyTfIhk5IDILWVW/HkpqmjMe6ZTZA8Bb3NCVwH/58yJyHNAAPBS65u2qepiqriGwCr/rlbFjn3izqvap6gpVXeOueRh4k6puJHBTn+vuVQ28EnhqNplF5Hzg426N4dB4k4gUu+MjgbU4xS4ize53A4HF+k132d3umXHv4JfunfwUeJmIVDlFfTbw5BxyrXDeAAgs3ZvneF+7gHUi4k3e1wHbZruHYRyIdDsL2f+OK74SV+YYcqkldRUo+dr2tFJEOgjct//gttzUuXP3icgqN/XvgY+IyA6C+Gk4jnk5cKu3QrPkUkIKeQ6+DtSISBvwCPD/VHWLk/H7BF8EjnWyv9td8zUCq/vnadubzgK2iMhmgtj0VaF49FdF5Engt8DnVHW7G/8WsNw9+0dwFqpzc3/JybQZeNTHv0XkRvdeq5xc17m1zgGeFpHtBAlcn53twd0Xpk8BD4rIFuAk4P9k+d4M44BgbCI1pei8pRxXMjWW8FiDicIlMpd1GFV9kRkykFX1gtDxs2TexoOqXjfHPb5NsD0pPDZrDFhVzwkdDxJsfco074oZxo+eYfxOnBt4HmuNznL//yDY+pQ+/nECCz19/A6ms9Uzkv6+VPUm4KaZ5hvGgU53KG7c1R9zC3lo/9aLnsbqsqkGE8tryvMtmrEICi2pyzAMIxLCburumCd1+Z7HmSzkeufGNrd14WEK2TAMA+jqD9zURRJ/C7l3eHy/xhKeRuv4VLCYQjYMw2DaKj66uSb+FnKGxhIeazBRuJhCNgzDILCKReC4lXV098c7qSsxtH9jCc9UT2RTyAWHKWTDMAyCvceNVWUcWl9B92CS+W3oyC+J4f0bS3i85Wwu68LDFLJhGAZBUldTbTnNtRWMTyq9MVZomRpLeKzBROFiCtkwDAPoHhilua6C5tpgq1BXjIuDBDHk/YuCeBqtwURBYgrZMAyDQAE31ZTTNKWQ4xlHnphM0TcyPmMMGaC+qsxiyAWIKWTDMA56Uillz2CS5rryKQs5ruUz+0YCV/pMMWR/zlzWhYcpZMMwDnp6R8YZn1Saa8tprqsA4uuynq2xhKehusySugoQU8iGYRz0ePd0U2051WXFVJYWx7Y4SM/QzI0lPA1VpRZDLkBMIRuGcdDj3dPNtRWICM115bEtDpKYpWymp6GqjP5RazBRaJhCNgzjoMdbwz6hq7m2fKqUZtyYrbGEJ9xgwigcTCEbhnHQ461hn9DVVBtfC3m2xhKeqWpdlthVUJhCNgzjoKerP0l1WTHV5UFH2ubaCrpjGkNODI3N2FjC4+PLPt5sFAamkA3DOOjpGhidcldDYCEPJCcYGZtcQqkykxgen7GxhKehyizkQsQUsmEYBz3dA0maayumPjfFeC9yYmhs1i1PEGrBaJnWBYUpZMMwDnq6B5I01U1byM0xrtbVMzw2a/wYQi0YzUIuKEwhG4Zx0OPLZnqaYlzPuneWxhIe32Aizg0yjP0xhWwYxkHN8NgEg8kJmvexkAP3dRxd1nM1lvBYg4nCwxSyYRgHNV7phi3kxuoyioskdi7rbBpLeKzBROFhCtkwjIOaqSpdddNJXcVFwvLqsthZyL1ZNJbwNFaXWQy5wDCFbBjGQU3XwL5FQTzNdeWxiyH3ZtFYwtNQXWYx5ALDFLJhGAc1vkRmU7pCrq2IXYMJX+hjrn3IwRxrMFFomEI2DOOgpnswSUmR7KfkmmriVz7TK9j6LJK66qvK6BuxBhOFhClkwzAOarr6k6yoKaeoSPYZb64rZ+9gksmULpFk++Nd1tnGkMEaTBQSppANwzio6RpI7ueuhsCFnVLYGyMrOZvGEh5rMFF4mEI2DOOgJiibub9Cbo5hcZBsGkt4rMFE4WEK2TCMg5qugeQ+RUE8TTEsDtIzNHdjCc9U+UxL7CoYTCEbhnHQMjGZYu/QvmUzPXGsZ907PHdjCY+PIfeay7pgyJtCFpHzReRpEdkhItdkOP9lEdnsfraLSG/o3GTo3N2h8d+ExjtF5Edu/BwR6Qudu3YuOSTgs+7e20Tkg6Hxf3bzt4jIKW78cBF51K3fJiJXhdb6rIi0i8jgDO/iEhFREdngPr89JOtmEUmJyEnu3BUi8oS79/0issKN3xaa/7yIbHbjrxORTe6aTSJybui+M631VvcMKS+TYRwM9AyNoQpNoaIgnjh2fOoZHssqoQuswUQhUpKPm4hIMfB14HVAB/CIiNytqk/6Oar64dD8vwFODi0xoqonpa+rqmeGrrkT+K/Q6d+o6p/PQ453Aq3AcaqaEpFmd9kbgLXu53Tg39zv3cAZqpoUkRpgq1urE7gH+BrwTIZ3UQt8CPh96DluAW5x518G/EhVN4tICfBVYJ2q7hGRG4GrgetU9bLQml8E+tzHPcAbVbVTRE4Afgqsnm0tYCvwZuD/pstrGAcyXRnKZnoqSoupqyiJXQy5paEqq7m+wYSVzywc5rSQReQYEfmFiGx1n08UkX+Y531OA3ao6rOqOgbcClw4y/wrgO9nu7iI1AHnAj9ahBzvA65X1RSAqna58QuB72rAw0C9iByqqmOq6v+llhN6l6r6sKrunkGGTwM3ADP5wa5wcgGI+6kWEQHqgM60ZxfgUtz7UtXH3JcCgDagUkTKZ1tLVbep6tMzyGMY80ZVuX/ri7HfAztdNnN/hQyBlRwnCzkxPJ5VYwlPY1UZCavWVTBk47L+d+ATwDiAqm4BLp/nfVYD7aHPHW5sP0TkcOAI4Jeh4QoR2SgiD4vIRRkuuwj4har2h8bOEJHHReQnIrI+CzmOAi5z9/mJiKyd6xoRaRWRLe78DSFFmBHn7m5V1XtnmXYZ08p1nOCLwhMEynMd8K20+WcCL6nqftY4cAnwqKoms1xrTkTkve4dbezu7p7v5cZBwqYXElz1H5t44Ol4/x3x8eFMWdbBeEVsLGTfWCLbGDIEW5/MQi4cslHIVar6h7SxiSiEcVwO3KGqk6Gxw1V1A/A24CsiclTaNekW9aPumpcD/8LcljMEVu6ou8+/AzfPdYGqtqvqicDRwJUicshMc0WkCPgS8NFZ5pwODKuq90aUEijRk4FVwBaCL0dhMnoT3JeQG4C/nsdac6Kq31DVDaq6oampab6XGwcJz+0ZAuCFvUNLLMns+NKYKzK4rMHXs45HUpdvLJHNHmRPQ5U1mCgkslHIe5wCVAAReQtB/HQ+7CKIz3pa3FgmLidNwajqLvf7WeBXhOLLLjHpNODe0Px+VR10x/cBpW7ebHJ0AD90x3cBJ2Yru7OMtxJYqzNRC5wA/EpEngdeCdydlkSV/uwnufV3qqoCtwOvCj17CUHs97bwjUSkxT3DO1R1ZzZrGUYuaU+MBL97hpdYktnpHkyyrLKUitLM+3qbagKXdfBPZmnxlq5ZyAcu2SjkDxAk+xwnIruAvyWwtObDI8BaETlCRMoIFM/d6ZNE5DigAXgoNNbgYqBe+f4J8GTosrcAP1bV0dA1K12cFBE5jeA5984hx4+A17jjs4Ht7vhu4B0u2/qVQJ+q7haRFhGp9DICrwZmjMOqap+qrlDVNaq6BngYeJOqbnRrFBHEgm8NXbYLWCci3hR9HbAtdP61wFOq2hF69nqCLyfXqOpv57GWYeSMDqeIvWKOK139mYuCeJrryhkdTzGQjNIpmB0+FpztPuRgbqnFkAuIObOsnVX6WhGpBopUdWC+N1HVCRG5miDjtxi4WVXbROR6YKOqeqV4OXCr7vt19Hjg/4pIikCxfi6cne2u+VzaLd8CvE9EJoAR4HK3ZkY53DWfA24RkQ8Dg8B73Ph9wAXADmAYeFdIri+KiBIkS31BVZ8AcBnMbwOqRKQD+KaqXjfHazoLaHfv27+3ThH5FPCgiIwDLxBkg4efPd1dfTWBC/1amd7u9frZ1hKRiwlc+03AvSKyWVX/dA55DWNG2hNOIcfcQu4aGM1YNtPjz3X1J6mryD6ZKgrm01jCE24wUVJsZSfizpwKOfSfuv8MgKpeP58bOdfxfWlj16Z9vi7Ddb8DXjbLuudkGPsawbajrORw473An2UYVwIvQfr4z5l2a6ef+zjw8ZlkziS3qv6KwI2dPu8m4KYZ1nhnhrHPAJ+ZYX7GtVT1LgIXt2HkhA5nGXckRlDVqf834kb3YJJTD2uY8XxzqFrX0c01+RIrI4l5NJbwTBUHGRmfMU5uxIdsvjINhX4mCfblrolQJsMwCpjkxCQv9o/SWF3GyPgke2Maw1RVuvozN5bwxKlaV2IejSU8DVatq6DIxmX9xfBnEfkCgcvXMAxjPzp7R1GFVx7ZyH1PvEh7z3AsrbP+0QmSE6kpKzgTcarWlRgao7K0OKvGEp7GqXrWFkcuBBYSVKgiyDQ2DMPYDx83PuPI5cHnmCZ2zVUUBGBZZSllJUWxUMg9Q+NTHZyypX6q45NZyIVANjHkJ3BbnggSoZqAecWPDcM4ePAJXWcc5RRyTBO7vBs6U9lMj4jQVFMei+IgiXk0lvA0Wk/kgiKbWtbhetATBFWhln4PgGEYsaS9Z4TSYuGIFTUsry6jIxFPhZyNhQzxKZ+ZmEdjCY+PN5tCLgxmVMgi0ugO07c51YkIqtoTnViGYRQqHYlhVtVXUlwktDRW0d4Tb5d10ywxZAgSu56PQcWx+TSW8FiDicJiNgt5E4GrOtN+BQWOjEQiwzAKmvbECK1OcbQ2VLJ1V98cVywNXQNJykqKqKuY3VHYXFfOH55fevujZ2hsXo0lPI1VZZbUVSDM+DdRVY/IpyCGYRwYdPQM8/r1QVn31sYqftr2IpMppbgoXnuRuweCKl1z7ZFuqqmgd3ic5MQk5SXZZzjnkonJFP2jE/OOIUOw9cm2PRUGWfVDdqUh1wJTvh1VfTAqoQzDKEyGkhPsDblWWxoqGZ9UXuofZVV95RJLty9dA6Ozls30+BjznsExVi/RM/jGEvONIftrrMFEYZBNP+T3AA8S7D3+lPt9XbRiGYZRiPgKXS0NgeLyrus4ZlrPVRTE47Owu/qXrjhIYqps5vwVcn2VNZgoFLLZh/wh4BXAC6r6GoJOS72RSmUYRkHiFW9rY9U+v+O4F7l7MDlrURCPt5CXMtPa7yOeT2MJT2NVqe1DLhCyUcijvpOSiJSr6lPAsdGKZRhGIeL3IHvLeFV9BSLxs5CTE5P0Do9nZSF7pb2Ue5F9x6aG6vkndTVUl9E/OsHEZCrXYhk5JpsYcodr6fcj4OcikiDoFGQYhrEPHYkRKkuLWVETWHLlJcWsrKuYUtRxYWoPchYKeXlNGSJLrZDnX8fa46+xBhPxJ5ta1he7w+tE5AFgGXB/pFIZhlGQtPcM09JQuU/mcmtD1VRsOS5kWxQEoLS4iMaqsli4rBekkH21rqExU8gxJ5ukrn8WkVcBqOqvVfVuVbWAhGEY+9GeGJmKG3taGivpiJnL2lu7TTVzx5DBV+tauqSu3uH5N5bwNE5V67K9yHEnmxjyJuAfRGSniHxBRDZELZRhGIWHqtLhLOQwLQ1V7O4fZWwiPjHM+VjIsPTlMxfSWMJjDSYKhzkVsqp+R1UvIMi0fhq4QUSeiVwywzAKir6RcQaSE1MJXZ7WhkpUobM3Pm7rroEkIrA8y329zbUVSx5DXkhRELAGE4XEfNovHg0cBxwOPBWNOIZhFCq+ZnVr474W8vTWp/i4rbsHRlleXUZJcXb/BTbXBRZyKqVzT46AhTSW8FiDicIhmxjyjc4ivh54Atigqm+MXDLDMAoK39UpvQHClEKOUZOJ7oHknE0lwjTVlDOR0qmKWfkmMTS2oIQuCBpMVJYWW3GQAiCbbU87gTNUdU/UwhiGUbhM7UFOS+paWVdBabHEqg1j10B2Vbo8PtbcNTC6YEt1MfQMjS04hgzQUFVqDSYKgGxiyP/XlLFhGHPR3jNCXUUJyyr3VRzFRcKq+spYVevq6k9mtQfZM10+M/9x5MU0lvA0VJeZy7oAmE8M2TAMY0baE8P7WceelobK2FTrSqWUPYPzU8jNdYF7eykyrRfTWMLTaAq5IDCFbBhGTmjPsOXJExQHiYdCTgyPMZHS+bmsa73LOv8KeTGNJTzWYKIwmDOGLCKHZRpX1T/mXhzDMAoRVaUjMcJrjm3OeL61sYo9g2MMj01QVZZV19fI6Joqm5l9Uld1eQlVZcV0LUFxkMU0lvBYg4nCIJt/GfcCCghBP+QjCPYjr49QLsMwCojuwSTJidSsLmsIal0fc0htPkXbj/kWBfE0L1FxkMU0lvCEG0xku9XLyD/ZJHW9TFVPdL/XAqcBD0UvmmEYhcJMe5A901uflusx0/gAACAASURBVN5tPV02c74KeWmKg/jY72JjyMCSbdsysmPeX5VU9VHg9AhkMQyjQOlIa7uYjh+PQ5OJhVrITbXl7FkChbyYxhIeH3+2OHK8ySaG/JHQxyLgFKAzMokMwyg4vOWbXhTEs6KmjIrSophYyKPUlJfMO5bdVFvOr7cvTVJXZWkxFaXzbyzh8fFniyPHm2z+RoYDPhMEMeU7oxHHMIxCpL1nhBU15TN2IxIRWhqqYlE+c75FQTzNdeUMJifynpiWGB5fdDESH3+2jk/xJpt+yJ/KhyCGYRQu7YmZtzx5WhsqY1E+s3uBCtnHnLsHkhy+PJ8KeWyqY9NCsXrWhUE2taw3iMhdIvKoiGzxP/kQzjCMwmC2oiCe1sZ4WMgLVci+OEi+E7t6hhbeWMLTYC7rgiCbpK5bgP8HXAK8MfSzaETkOBF5SESSIvKxWeZ9W0SeE5HN7uckN/529wXhCRH5nYi8PHTN+SLytIjsEJFrQuNHiMjv3fhtIlLmxq9y62wWkf8RkXVu/HUissmd2yQi54bWOtWN7xCRfxYRCZ37GxF5SkTaRORGN7ZGREZCz3FTaP5l7lnaROSGtOe/VESedOf+MzR+oxvbFr6/iFzh5NoiIveLyAo3/lY3P5Wpr7WIHCYig7P9WRhGOhOTKXb3jtI6p4VcxcDoBH1L7Dbt6h+dV5Uuz1KVz+wdXnhjCY9vMNFrFnKsyUYhd6vq3ar6nKq+4H9ydP8e4IPAF7KY+3eqepL72ezGngPOVtWXAZ8GvgEgIsXA14E3AOuAK7yCBW4AvqyqRwMJ4N1u/D/d1q6TgBuBL7nxPcAb3T2uBL4XkunfgL8C1rqf8939XwNcCLxcVdenPd/O0HNc5eYvBz4PnOfmrxSR89y5tcAngD9x5/7Wjb8K+BPgROAEgn7VZ4tICfBV4DWqeiKwBbja3Xsr8GbgwRne8ZeAn8xwzjAy8mL/KBMpzcJCDhT2UlrJQ8kJhsYm51UUxOOzsrvzXBwkFxYyBFufrMFEvMlGIf+jiHzTWV1v9j+5uLmqdqnqI8CC/pao6u9UNeE+Pgy0uOPTgB2q+qyqjgG3Ahc6C/Jc4A437zvARW6t/tDS1QTFUFDVx1TVZ5W3AZUiUi4ihwJ1qvqwqirwXb8W8D7gc6qa9M85x6McCTyjqt3u838TeCQgUPhf988ZWksJCrWUAeVAKfASQQEXAard89bhsuJVdZuqPp1JABG5iOALTtscshp55po7t/C9h55fajFmZGoP8gwZ1p6Wqa1PS6eQ/ZanhbisG6vKKC6SvLqsx11jicXGkAHqq0othhxzslHI7wJOIrD+vLv6z6MUagY+61ywXxaRTP+a3s20dbcaaA+d63Bjy4FeVZ1IGwdARD4gIjsJLOQPZrjHJcCjTtGudten3wPgGOBM5xr/tYi8IjTvCBF5zI2f6cZ2AMc6l3YJgWJvDa11jIj8VkQeFpHzAVT1IeABYLf7+alTuOMEXwieIFDE64BvZXiWKUSkBvh7YM4EPhF5r4hsFJGN3d3dc003FklyYpI7NnVw12O7llqUGZluuzi3yxqWti9y96Avmzl/hVxUJKyoKctrta7e4cU3lvAEFrIp5DiTTargK1T12MglmZ1PAC8SWIPfIFAe1/uTzkX8buDVi7mJqn4d+LqIvA34BwIXtb/HegJ39+uzWKoEaAReSeBKvl1EjiRQnIep6l4RORX4kYisV9WEiLwPuA1IAb8DjgqttRY4h8AD8KCIvAxYARzPtFfg507BP0ygkE8GngX+heD9fWYWea8jcOMPhsLgGVHVb+BCAxs2bNAs3oWxCLa/OMhEStm2e4DJlFJcNPufz1LQ0TNMkcChy2ZXyMuqSqmtKFlSl7WP/863KIgn39W6fMx3sTFkv0Yc9oEbM5ONhfy7UPx10Tgr1Cc1rcrmGlXdrQFJggSz00LrnQh8E7hQVfe64V1MW5gQKK1dwF6g3lmh4fF0bmXa/YyItAB3Ae9Q1Z2he7SErgmv1QH80Mn8BwIlu0JVk15GVd0E7CSwgFHVe1T1dFU9g6BW+PbQWner6riqPufG1wIXAw+r6qCqDhJ4B84g8GagqjudK/124FUZX+w0pwM3isjzBDHq/09Erp79EiMftHX2ATAyPslze4aWWJrMtCdGWFlXQVnJ3P+dtDZULalS8M0h5ls209NcW55XhZyLKl2eBmswEXuyUcivBDa7jGWf0bzgbU+q+vVQUlNWFb9cvBYXE72IIDnJd6L6IfCXqro9dMkjwFoJMqrLgMsJlJoSuHnf4uZdCfyXW2tt6Po/A55x4/UExVCuUdXfhp5jN9AvIq90cr3DrwX8CHiNu/4YAst+j4g0uYQznMW8lsCKRUSa3e8G4P0EXzL8Wue4cysIFPizwB9xSVwiUgqcDWwj+FKwTkSa3PWvc+MzoqpnquoaVV0DfAX4P6r6tdmuMfLDVqeQYVo5x42OxDAtcyR0eVobK2lfwvKZ3QNJSopkwQquKc8NJnzMdzGNJTzhBhNGPMnGZX1+VDcXkZXARoLEo5SI/C2wTlX7ReQ+4D1Oad/iFIwAm4Gr3BLXEsSF/9W5WidUdYOqTjgL76dAMXCzqvpkpb8HbhWRzwCPMR1fvVpEXkuQYJZg2l19NXA0cK2IXOvGXu+Sq94PfBuoJLBQfQz7ZuBmEdkKjAFXqqqKyFnA9SIyTmA1X6WqPe6ar8r0tq3rQ18wfgq8XkSeBCYJss33isgdBAlqTxAkeN2vqve49/opAtf2OPAC8E43fjGBC7sJuFdENqvqn87xx2QsIW2d/Ww4vIEtu/po6+znwpNWz31RnmnvGeFPjl6R1dzWhip+vb0bVWWu8EgUdA0kWVFTTtECXf/NteXsHUrmrWtSIscxZAgaTKxYoIfAiJZsKnXlaotTprVfZF+3b/jcBaHjc2eY8x7gPTOcuw+4L8P4s4Rc3qHxD82wzmeYIf6qqhsJthylj48Bf5Fh/E5mKDuqqlfMMK7AR9xPeHwS+OsZrrkJuCnD+F0ErvcZUdXrZjtv5I/JlLJtdz9vO+1wkhOpWFrIyYlJXhoYnTOhy9PaWMXoeIruweSCth4tlq6B5ILjxwBNdRWowt6hMQ6pi17+XLqsww0mTCHHE2uMaRgx5dnuQUbHU6xfVcf6VXVs3dVP8P0sPuxKjKA695YnT7gv8lLQPZBcUIa1J1w+Mx/korGExxpMxB9TyIYRU9o6g63xJ6xexvrVy+gbGWdX79LXgg7j48FzFQXxLHVf5O6B0QXtQfZ467orT8VBeoZzUxQEwg0mTCHHlfxVSDcMY160dfZRXlLEUU3VDI1NuLH+GVscLgVesWbrsl5KC3liMsXeoTGaFuEq99Z1vspn9g6P5yShC6ZjyNbxKb6YhWwYMaWts5/jVtZSUlzE8SvrKJJpqzkutCeGKS2WrOPBVWUlrKhZmv2we4fGUF1YURDPijy7rHuGFl/H2mMNJuKPKWTDiCGqytZdfaxbtQwImgMc1VRD2654JXZ1JEZYXV85r4IlS9UX2Vu1i3FZV5QWs6yyNG97kRM5aCzhqSgNGkwkTCHHFlPIhhFDOhIj9I9OcMLquqmx9avqYmchd/TM3XYxndbGqiUpn9k9GMR9F2MhQ373Iidy1FjC01hdZi7rGGMK2TBiiFe8652F7I9f7B9l72B+2//NRntiZN4x7daGSjp7R5hM5TdjfLps5uK2KwXVuqJP6splYwmPNZiIN6aQDSOGPNnZR3GRcNzK2qmx9asCazkuVvJQcoKeobGsE7o8LQ1VTKSUF/vz28bQu5lX1CzO4sxX+cxcNpbwWIOJeGMK2TBiyNbOfo5qqt5n/6m3lrfGpEDIVJen+VrIvi9ynhO7ugeS1FeVUl6yuD293mUd9Z7wRA4bS3gaqsrMQo4xppANI4a0dfZxQshdDUG3pJaGythYyFN9kOcbQ25Ymr3IXQOjC24qEaa5toLkROBOjhKffJXzGLJZyLHFFLJhxIzugSQv9SdZt6puv3PrV9XxZEwUcoezkP3e4mxZVV+JCHlvMrHYspkev0Z3xHFkb8nmOobcPzrBuDWYiCWmkA0jZvia1evTLGQ/9tyeIQZGlz5Ttr1nhMrSYpbP04IrKyni0LoKOpbAZZ2L+tneyo46jtwzFE0MGabj00a8MIVsGDHDu6QzWch+G9S23QN5lSkT7YlhWhsrF9S1qaUxv3uRVZWugeSi9iB7pi3kaBVyVDFkgF6LI8cSU8iGETOe7OznsMYqllXu76r0VnMcOj+19wzPO6HL09JQmde9yP0jE4xNpBa9BxmYKr0ZdfnMXDaW8Fi1rnhjCtkwYkZbZ9/UFqd0mmvLWVFTtuSJXapKR2Jk3gldntaGKl4aGCU5MZljyTLji4LkwkKuqyihrKSI7oj3g+eysYTHGkzEG1PIhhEj+kfHeX7v8IwKWURYv2oZW5e4hGbv8DiDyYl5J3R5WhurUIXO3vzsRc5F2UyPiAR7kSPeR50YGstZYwmPV/A+Pm3EC1PIhhEjtvkKXav3T+jyrF9Vx46uwbxZl5lon8qwXqiFnN+9yN6azUVSF7i9yBFbyInh8ZzGj2HaZW0WcjwxhWwYMWK6ZGZmCzk4t4yJlLL9xcF8ibUfHVN9kBduIQN5S+zKpYUMrlpX1DHkHDaW8FiDiXhjCtkwYkRbZz9NteWzWnLTJTSXzm093Qd5YRbyIXUVlBZL3hK7ugZGKS8poq4iNy3gm2sr8rDtKfcxZHDlM81CjiWmkA0jRsyW0OU5rLGK2vKSJS2h2Z4YZlllKXUVC4txFhcJq+sr82Yhd7uiIAvZopWJptpy+kbGIwsbjE+mGBidyLmFDEFil+1DjiemkA0jJoyOT/JM1+B+JTPTKSoSjl/iVoztPSMLdld7WhqqplzfUdM1kMxJ2UyP3z4V1V5krzBzndQFQRzZtj3FE1PIhhETtr80wGRK57SQIXBbb9vdn/cWhp72xML3IHtaGyvzVq2rK0dVujy+OEhUbusoioJ4rMFEfDGFbBgxYeuu/Xsgz8QJq5YxOp7i2e78J3alUovbg+xpaahi79AYQ8lomzTAtMs6VzTVVEytGwU9ETSW8FgLxvhiCtkwYkJbZx+1FSVZuYLXr1663sh7BpOMTaQWvAfZ4xV61G7r0fFJ+kbGc+uyjthC7o2gsYSnvqqUAWswEUtMIRtGTGjr7Gf9qrqsEo+OaqqhrKRoSTKtF9oHOZ187UX2VmwuLeTl1WWIQHdExUGiaCzhsQYT8cUUsmHEgInJFE+92J+VuxqgtLiI41bWLomFPN0HOTcWctSZ1r6AR672IAOUFBexvLossuIgUceQw/cw4oMpZMOIAc/uGWJ0PJVVQpfHl9BUzW9il7doF1qly7O8uozK0uLI9yL7Ah65TOqCoMlEVMVBEkNjVJXltrGEx1vIVhwkfphCNowY4F3PJ8xSMjOd9avq6B+dyNvWIU97Ypim2vJFKwsRoaWhko48Wci56PQUJsrymT0RVOny+Li0WcjxwxSyYcSAtl39lJcUceSK6qyvma7YlV+3dXvPyFT8d7G0NlbRHvEXiu7+UURyH4+NsnxmFI0lPNZgIr6YQjaMGNDW2c9xh9ZRUpz9P8njD62juEh4Ms+JXe2J4UVvefK0NgR7kaN0u3cNJFleXT6vd5sNzbXl7BlMkopgL3hPBI0lPBZDji95UcgicrOIdInI1hnOv11EtojIEyLyOxF5eejc8258s4hsDI1/2l2zWUR+JiKrQufOceNtIvLr0PiH3dhWEfm+iFS48XNF5FE3/h0RKclCrg+5+W0i8reh8etEZJe7/2YRucCNLxeRB0RkUES+lvb8l7n7tInIDaHxw9w1j7nzfq01IjISusdNoWs+KyLtIjKYdo/DReQXbp1fiUiLGz9JRB5y994iIpfN/Sdq5BJVzapkZjoVpcUc1VTN1jxayBOTKXb3jS56y5OntbGKgeQEfSPRWWvdA8mcu6shcFlPpDQSxdYbocvaGkzEl3xZyN8Gzp/l/HPA2ar6MuDTwDfSzr9GVU9S1Q2hsc+r6omqehLwY+BaABGpB/4VeJOqrgfe6sZXAx8ENqjqCUAxcLmIFAHfAS534y8AV84ml4icAPwVcBrwcuDPReTokGxfdvKepKr3ubFR4H8DHws/mIgsBz4PnOfkXSki57nT/wDcrqonA5e75/LsDN3jqtD4PU6udL4AfFdVTwSuB/7JjQ8D73D3Ph/4inuHRp7oSIzQPzoxZ8nMTKxftSyvW592940ymdJFb3ny+MSwKBO7ugaSOc2w9vgksSj2IkfVWMJjDSbiSV4Usqo+CPTMcv53qppwHx8GWrJYM2wWVAPeb/Q24Ieq+kc3rys0rwSodBZwFdAJLAfGVHW7m/Nz4JI55Doe+L2qDqvqBPBr4M1zyDukqv9DoJjDHAk8o6rd7vN/+/u7Z/Jm0zIn76yo6sOqujvDqXXAL93xA8CFbv52VX3GHXcCXUDTXPcpRPqGx9nc3rvUYuyHV6jztZD9NS/1JyOrGJXO1B7kHLmsvaUd5danroHRSCzkqIqDRNlYwtNQXWoWcgyJYwz53cBPQp8V+JmIbBKR94Ynevcs8HachQwcAzQ4t+wmEXkHgKruIrAS/wjsBvpU9WfAHqBERLz1/RagdQ65tgJnOjd0FXBB2jVXO/fvzSLSMMfz7gCOdW7oEuCi0FrXAX8hIh3AfcDfhK47wrmyfy0iZ85xD4DHmf7ScDFQ66zzKUTkNKAM2DnTIiLyXhHZKCIbu7u7Z5oWSz5z75O85d9+R1dExRwWytZd/RQXCceurJ33tX7fcr6s5A6/BzlHFvLUXuSIioOkUsqewbGcFgXx+Mpfuf4yFGVjCU9Qz9qSuuJGrBSyiLyGQPH9fWj41ap6CvAG4AMicpY/oaqfVNVW4BbgajdcApwK/Bnwp8D/FpFjnGK8EDgCWAVUi8hfaJBNcjnwZRH5AzAA7NNTLV0uVd0G3AD8DLgf2By65t+Ao4CTCBT/F2d7ZmeBvw+4DfgN8HxorSuAb6tqC4HS/55zse8GDnOu7I8A/ykic5lXHwPOFpHHgLOBXeHnFJFDge8B71LVGWvqqeo3VHWDqm5oaiocQ3owOcGPt+xmIqX88LFdSy3OPrR19nF0U82CthGty3OmdXtimCKBQ+tzs6c3aOFYEtnWrZ7hMSZTmtOymZ5pCzm3X/CiLArisQYT8SQ2CllETgS+CVyoqnv9uLNsvev5LjLHR29h2s3bAfzUuYj3AA8SxHlfCzynqt2qOg78EHiVW/shVT1TVU9z8737eja5vqWqp6rqWUDCX6OqL6nqpFNq/z6DvPugqveo6umqegbwdOj+7wZu9zICFcAKVU16WVR1E4FFe8wc9+hU1Tc7Jf5JN9brnrEOuBf4pKo+PJe8hciPH+9kZHySQ+rKuf2R9rwX05iNts7+qdrU82VZZSmtjZU8mS+F3DPMocsqKc1hxnKw9SkaC3mqKEhdbouCAFSVlVBTXpLzrU9RNpbwWIOJeBILhSwihxEoyL8MxXIRkWoRqfXHwOsJ3MWIyNrQEhcCT7nj/wJeLSIlzp18OrCNwFX9ShGpEhEBznPjiEiz+11OYAXfNJtcadccRuAK/k/3+dDQtIu9vHM8v1+rAXg/wRcAnMznuXPHEyjkbhFpEpFiN34ksBZ4do57rHDWNcAngJvdeBnBF53vquodc8laqNy2sZ21zTV89PXH8uyeITa9kJj7ojzQNTBK10Ay65KZmTghj4ldQZen3GRYe1obqiJzWUdRNjNMFMVBfGw3agvZGkzEj3xte/o+8BBBrLRDRN4tIleJiM8OvpYguepfZd/tTYcA/yMijwN/AO5V1fvduc+5bUdbCBT1h2DKnXw/sMVd801V3aqqvwfuAB4FniB4dp/N/Xciss1dc4+q+uSnmeQCuFNEniTIav6AtzaBGyXYJrUFeA3w4dB7eB74EvBO9x7WuVNfdWv9FvhcSPl/FPgr9/zfB97pXOxnAVtEZLN7pqtUtcfd40YXc65y97jOrXUO8LSIbHfv9bNu/FK33jtD26hO4gDimZcGeOyPvVy6oZU/e9mhVJcVc9sj7UstFjDtal5IQpdn/ao6nt87TP9o9DHB9sTwoktmptPaWElHYiQSr4XPF4giqQucQs6xhZzIRwzZrW0NJuJFST5uoqpXzHH+PcB7Mow/S+BuznTNJZnG3bnPE2wlSh//R+AfM4z/HfB32crlzmVMpFLVv5xFrjUzjGd8P6r6JPAnGcbvBO6c4ZqPAx/PMH4HgfJOH/8P4D9mkvlA4PaN7ZQUCRefsprq8hLe+PJV3P14J//4pvXUlOfln8CMeFfzukUp5MC63tbZz+lHLp9j9sIZHZ/kpf5kzhK6PK2NVSQnUq5ncW5dy/mwkHMdLshXDNnfK6p3Y8yfWLisDSMqxidT/PDRXZx3fDMrXGLPWze0Mjw2yb1b5txFFjltnX0cvryKuoqFW0P5KqG5qzc3XZ7SiXLrU1d/kpryEqrKovniFZTPzG1SV0+EjSU80+UzLY4cJ0whGwc0v9jWxd6hMS57xfSutFMOq+fo5ppYuK19D+TF0FxXQVNteeQK2cd5c7UH2dMaYXGQqKp0eZprKxgam2QoOZGzNRMRVuny+AYTvZZpHStMIRsHND/Y2M4hdeWctXZ6i5aIcOmGFh79Yy87ugaWTLb+0XFe2Du8qIQuz/pVdZEndvkmELl2WfuYdBRdn7ojqtLl8Wvnci9ylI0lPNZgIp6YQjYOWF7qH+WBp7u45JSW/RoLXHxyCyVFwu0bO5ZIuun48WItZL/GM12DjI5Pzj15gXQkhikrKcq5xVlZVsyKmvJILOSugdFIFbJ/F7ms1hVlYwmPNZiIJ6aQjQOWOx/tIKVw6Yb9C6811ZZz7nHN/PDRjiXb+jGdYZ0LC3kZkyll+0vRWfwdPSO01FdSVCQ5X7u1sTKaGPJAcqrmdBREURwkEXEdawgaTFSVFVsMOWaYQjYOSFSVH2zs4LQjGlkzQ4/hy17Ryp7BMX75VFfG81HTtquP5trynFhwvjHF1l3RxZHbE8OszlGXp3RaG3JfHGQoOcHw2GS0LusIymfmI4YMVq0rjphCNg5I/vBcD8/tGeKyDNax5+xjmmiuLecHG5cmuSsXCV2e1sZKaitKIo0jt/fkrg9yOi0NlXT2jjKRQ2+FdyNHmdTVUFVGSZHkzGWdj8YSHmswET9MIRsHJLdv7KCmvIQ3vGzljHNKiou45NQWHni6O+8NJ0bHJ9nRPcgJqxfvroYgUW3doXWRZVoPJidIDI/nPKHL09pYxWRK2d2Xuz8Hb7VG0VjCU1QkQXGQHClkb7E2RpzUBcGXiR4rDBIrTCEbBxwDo+Pc98Ru3vjyVXPuP710QyuTKeXOR/PbcOLpFweYTGnOLGSAE1Yv46kX+3NqZXqmtzxF57KG3O5F9nHdqAtfNNWW58xC9pWz6vPksrZtT/HCFLJxwPHjLbsZGZ/k0g1zttXmiBXVnLamkR9szG/Dia1TPZBzYyEHa9UxOp7i2T1DOVvTM6WQI7OQA0Wfy65PU40lIkzqCtbPXXGQfDSW8FiDifhhCtk44LjtkXaOOaSGk1rrs5p/6StaeXbPEBvz2HCirbOfuoqSqSpVuSDK3sheUUYVQ15VX0mRQEcOm0x0DyYpLRbqK6N1/zbVlrMnRw0m8tFYwmMNJuKHKWTjgGL7SwNsbg8aSQRNvebmgpetpKa8JK+Vu4KErmVZy5gNRzVVU15SRFsEmdbtiWGqyoppqIpGuZUWF3Hossqp4iO5oKs/yYqa8ki2aYVpqq1g79BYTkIFPb6OdT5iyO4elmkdH0whGwcUtz/STmmxcPHJq7O+pqqshDe+/FDu3bKbwRyWQJyJickUT+3OXYa1p6S4iOMiSuxq7xmhtaEqp18g0mlpqMxpG8augdFIM6w9zbXlqMKewcUrNh9DzpeFHL6nsfSYQjYOGMYmUtz12C5ee/whLK+Z33/Eb93Qysj4JD9+PPqGEzu7h0hOpFi/OrcKGaZLaOY6Ht6RGI4socvTkuO9yEHZzGjjx5Db8pn5aCzhsQYT8cMUsnHA8MunXmLv0FjGylxzcXJrPWuba7g9D3uSfYz3hBwmdHnWr6qjf3Qip8lRqkp7T+77IKfT2ljJS/3JnJX/jLqOtWe6fObiE7sSQ/kpCgKh8pmmkGODKWTjgOG2R9pZWVfBWcc0zT05jaDhRGteGk60dfZTUVrEkU01OV/7hAgSuxLD4wyNTUaW0OXxGdydvYv/MjE+mWLv0Fh+XNauh3Mutj4lhqNvLOGZjiGbyzoumEI2Dghe7Bvl19u7ueTU1RQvMInn4lNW56XhxNZdfRy3sm7Bcs7GsStrKS6SnJbQ9F2YWiMqm+nxCj8XiV17XTw3HxbyiprA0syJyzoPjSU81mAifphCNg4IfCOJt546f3e1Z0VNOecdH23DCVXlyQgSujwVpcUc3VSTUwvZd2GK3EJ2MepcJHZ593E+LOTykmLqq0pz5rLOxx5ksAYTccQUslHwqCq3b2zn9FkaSWSLbzjxi23RNJxo7xlhYHQiZyUzMxEkduXOQvaJVrncM52JQ2orKC2WnCR2TZfNjD6pCwLFnwsLOV+NJTzWYCJemEI2Cp7fP9fDC3uHuewVC7eOPWetjbbhRNtUha5oLGSA9auX0TWQzFlLwPaeYeqrSqmtiDa2WVQkrK6vpCMHfZF9PDcfLmt/n8XGkPPZWMJjDSbihSlko+C5fWM7teUlvOGEQxe9VklxEW85tYUHnu7ipQgaTmzt7KO4SDjmkNqcr+3xyj5XVnJ7YiSykpnptDbmZuuTL5vZNM/tbwulubZi6p4LJZ+NJTzWYCJemEI2Cpp+30jipFVUluVm7+ZbN7SS0iAunWvaqQ4RzwAAD2FJREFUOvtZ21wT6T7TdU4hP5kjhdzRE/0eZE9LQ1VOYsjdg6M0VJVSVpKf/+KaasvpHkwuav93YsgVBclTDBmCvchmIccHU8hGQfPjx3czOp5a0N7jmThiRTWnHdHIDzZ25LzAhi+ZGSV1FaUcvrwqJ4ldqZTSkVcLuZLE8PiiK6Z19ednD7KnubacsYkU/SMLl9tbyBZDPngxhWwUNLdtbOfYQ2p5eUtuldylG1p5bs8Qjzyfu4YTXf2jdA8kI40fe3KV2NU9mGRsMkVLxBnWHq/4Oxbptu4aSEbe5SlMUw6Kg+SzsYTHGkzEC1PIRsHy9IsDPN7ey6WvyL6RRLZE0XDCK8j8KORlvLB3mP7RxcUHp9su5sdlPbUXeZGJXfmq0uXJRfnMnuH8tV70NFqDiVhhCtkoWG7fOP9GEtniG07c98RuBhap1DzehbwuDwo5V3Hk6S1P+bKQF78XWVXpHkjmZQ+yx1vji8m09k0e6iPqqJWJemswEStMIRsFiW8k8bp1h0RmUVzqG05s2Z2T9do6+1mzvCry7UMQLqG5SIXsLNWo9yB7GqvLqCwtXlSmdd/IOGOTqfzGkOsW77LOZ2MJjzWYiBemkI2C5BfbXqJnaIy35jCZK52TctxwYmtnX+QJXZ6m2nKaa8tp27W4xK72nmGaa8vzpiREhNbGykW5rLvzvAcZoLa8hPKSokW5rPPZWMJjDSbihSlkoyC5baNrJLF2/o0kskVEuOwVrTz2x16eeWlxDSf6RsZp7xnJi7vak4vErvbEcOQlM9NpbahaVFKXdxvnM6lLRGiuW1xxkJ7h/JXN9ExZyBZDjgWmkI2CY3ffCA9u7+Ytp7ZE0qAhzEUn+4YTi7OSfSw3ypKZ6axftYwd3YOLamcYbHnKj7va09pYRUdiZMFbzqbLZubPQobgC8CiLOTh8bzGj2E6Xm0x5HhgCtkoOO7c5BpJbGiJ/F4rasp57fGH8MNHdzE2sfCtIfkomZnOCavrmEwpT7+4MOt+YjLF7r7RvFvILQ2VDCYnFqwkfBw3ny5rCKqCLcZCzmdjCY81mIgXppAPAkTkfBF5WkR2iMg1Gc6Xi8ht7vzvRWRN/qXMjlRKuX1jB688spHDly+ukUS2XPqKFvYOjfHLpxbecKKts59D6spZkadSjsBUvHrrAguE7O4bZTKleSsK4pluw7gwt3VXf5KK0iJqy0tyKdacNNeV07WIcqtLEUMGVxzEFHIsMIV8gCMixcDXgTcA64ArRGRd2rR3AwlVPRr4MnBDfqXMnt8/18Mfe3LTSCJbzlrbxCF15YtyW7flMaHL09JQSV1FyYLjyH7rUb4yrD0tU1ufFpbY1T0YFAXJ9d70uWiuLad/dGJBIYLxyRQDyfw2lvA0VJfaPuSYkN+vkMZScBqwQ1WfBRCRW4ELgSdDcy4ErnPHdwBfExHRXNeNdLzxX/5nwXHNnqExastLOH/94htJZItvOPFvv9rJ67706wWtsaN7kPPXr8yxZLMjIqxftYwfPbaLR57rmff1A6NBGci8J3W5+113Txtf+e/t875+V+8Ixx+av9CAx7vIL/jqb+ad2zDp/qk15LGxhKehqoyHnt274L/bByMNVWXcftUZOV/XFPKBz2ogbNp1AKfPNEdVJ0SkD1gO7ElfTETeC7wX4LDDDluQQEc1VTO2iFJ95x53SM4aSWTLlWesobN3lOTEwr5IHH9oHRdFUMBkLt579pGLaiW5sq6S1fX5tZDrKkr5m3OPZmf34IKuX3tIDX9+4qocSzU3Zx/TzMUnr17w35ETVy/j3OOacyzV3Fx5xhpqK0wVzIe6iGoJSERGkBETROQtwPmq+h73+S+B01X16tCcrW5Oh/u8083ZTyGH2bBhg27cuDE64Q3DMA4wRGSTqm7IdM5iyAc+u4BwwLXFjWWcIyIlwDJgb16kMwzDMABTyAcDjwBrReQIESkDLgfuTptzN3ClO34L8Muo4seGYRhGZixwcIDjYsJXAz8FioGbVbVNRK4HNqrq3cC3gO+JyA6gh0BpG4ZhGHnEFPJBgKreB9yXNnZt6HgUeGu+5TIMwzCmMZe1YRiGYcQAU8iGYRiGEQNMIRuGYRhGDDCFbBiGYRgxwAqDGAtGRLqBFxZ4+QoyVAKLASbX/DC55ofJNT8ORLkOV9WMjdxNIRtLgohsnKlazVJics0Pk2t+mFzz42CTy1zWhmEYhhEDTCEbhmEYRgwwhWwsFd9YagFmwOSaHybX/DC55sdBJZfFkA3DMAwjBpiFbBiGYRgxwBSyYRiGYcQAU8hGpIjI+SLytIjsEJFrMpwvF5Hb3Pnfi8iaPMjUKiIPiMiTItImIh/KMOccEekTkc3u59pMa0Ug2/Mi8oS758YM50VE/tm9ry0ickoeZDo29B42i0i/iPxt2py8vC8RuVlEukRka2isUUR+LiLPuN8NM1x7pZvzjIhcmWlOjuX6vIg85f6c7hKR+hmunfXPPAK5rhORXaE/qwtmuHbWf7sRyHVbSKbnRWTzDNdG+b4y/t+Qt79jqmo/9hPJD0G7x53AkUAZ8DiwLm3O+4Gb3PHlwG15kOtQ4BR3XAtszyDXOcCPl+CdPQ+smOX8BcBPAAFeCfx+Cf5MXyQobpD39wWcBZwCbA2N3Qhc446vAW7IcF0j8Kz73eCOGyKW6/VAiTu+IZNc2fyZRyDXdcDHsvhznvXfbq7lSjv/ReDaJXhfGf9vyNffMbOQjSg5Ddihqs+q6hhwK3Bh2pwLge+44zuA80REohRKVXer6qPueADYBqyO8p455ELguxrwMFAvIofm8f7nATtVdaEV2haFqj5I0LM7TPjv0HeAizJc+qfAz1W1R1UT8P+3d/cxdlR1GMe/T1+MsYW1AkHlJYBBFBRabRteutiYQsAXkvpKS7CIiWkiNvYfaWwiWBMFNDQgIAkvEdsVjIK1VoxVaGsFKy+tXaqt0hgSIHWrRFsL2sjy849zhp3e3rt72+7MvdDnk9zcmTNnZs49d+6cPWdm58evgIuqLFdErI6Il/PsBuD40drfoZSrTe38dispV/79fwq4d7T2165hzg21HGNukK1KxwHPluafY/+G79U8+eS1CziqltIBeYh8CvD7JovPkbRZ0i8knVFTkQJYLelJSZ9vsrydOq3SpbQ+UXaivgCOjYgdefpvwLFN8nS63q4kjWw0M9J3XoWr8lD63S2GXztZX73AQEQ83WJ5LfXVcG6o5Rhzg2yHLUkTgfuBL0XE7obFG0nDsmcB3wFW1FSsGRHxPuBi4AuSzq9pvyOS9AbgEuBHTRZ3qr72EWnssKv+l1PSYuBloK9Flrq/8+8C7wAmAztIw8PdZA7D944rr6/hzg1VHmNukK1KzwMnlOaPz2lN80gaB/QAL1RdMEnjST+4voh4oHF5ROyOiD15+kFgvKSjqy5XRDyf33cCPyENHZa1U6dVuRjYGBEDjQs6VV/ZQDFsn993NsnTkXqTdAXwEeCyfCLfTxvf+aiKiIGIGIyIV4A7WuyvU/U1DvgY8MNWeaqurxbnhlqOMTfIVqXHgVMlnZx7V5cCKxvyrASKuxE/ATzc6sQ1WvI1qruArRFxY4s8by2uZUuaTvqtVPqHgqQJko4opkk3BW1pyLYS+IySs4FdpaG0qrXsuXSivkrKx9A84KdN8vwSuFDSpDxEe2FOq4yki4AvA5dExEst8rTznY92ucr3HMxusb92frtVmAVsi4jnmi2sur6GOTfUc4xVcaeaX34VL9JdwX8h3bG5OKctIZ2kAN5IGgLdDjwGnFJDmWaQhpz6gT/k14eA+cD8nOcq4I+ku0s3AOfWUK5T8v42530X9VUul4Bbc30+BUyt6XucQGpge0pptdcX6Q+CHcD/SNfoPke65+Ah4Gng18Bbct6pwJ2lda/Mx9l24LM1lGs76ZpicYwV/03wduDB4b7zisu1LB87/aSG5m2N5crz+/12qyxXTv9ecUyV8tZZX63ODbUcY350ppmZWRfwkLWZmVkXcINsZmbWBdwgm5mZdQE3yGZmZl3ADbKZmVkXcINsZq95StGmVjVJnyrp5lKec0dxnydJmttsX2YHY1ynC2Bmrx+SxkbEYKfLUYiIJ4AiRN9MYA/waLvrSxoXQwEiGp0EzAV+0GRfZgfMPWQzG1HuDW6T1Cdpq6QfS3pTXvaMpOslbQQ+KWlOjle7RdL1pW3skbQ0x5l9SNIxOX2ypA0aihs8KacvyHFp+yXdl9OmS/qdpE2SHpV02gjlnilpVQ4UMB9YqBRHt1fSMZLul/R4fp2X17lW0jJJjwDL8mdfL2ljfhW97OuA3ry9heVeulL83BW57BsknVna9t2S1kr6q6QFo/Ud2WufG2Qza9dpwG0R8W5gNymWdeGFSA/8/w0p9u8HScELpkkqQtVNAJ6IiDOAdcA1Of37wNURcSbpCVJF+iJgSk6fn9O2Ab0RMQX4KvCNdgoeEc8AtwNLI2JyRKwHbsrz04CPA3eWVjkdmBURc0jPLb4gf75PA8Ww9CJgfd7e0oZdfg3YlMv+lfwZC+8iheqbDlyTn51s5iFrM2vbsxHxSJ5eDiwAvp3ni2AA04C1EfF3AEl9pGD0K4BXSvmWAw9I6gHeHBHrcvo9DEWT6gf6JK1gKHpUD3CPpFNJjzg8lMZsFnC6hsJvH6kU5QdgZUT8J0+PB26RNBkYBN7ZxrZnkBp5IuJhSUdJOjIv+3lE7AX2StpJCuXX9NnNdnhxg2xm7Wp8zm55/sVR2F6jD5Ma848CiyW9F/g6sCYiZudh6LUHsd/CGODsiPhvOTE30OXPsxAYAM7K6+yT/yDsLU0P4vOwZR6yNrN2nSjpnDw9F/htkzyPAR+QdLSksaQIUUXvdwwpoter60fELuCfknpz+uXAOkljgBMiYg1wNalnPDG/FyHtrjjA8v8bOKI0vxr4YjGTe8DN9AA7IoUrvBwY22J7ZeuBy/J2ZwL/iP1jbpvtww2ymbXrz6SA8FuBSaRA9/uIFApyEbCGFJHnyYgoQtW9CEyXtIV0jXlJTp8HfEtSP+m68xJSo7dc0lPAJuDmiPgXcAPwTUmbOPCe5c+A2cVNXaQh96n5xqs/MXSdutFtwDxJm0nXf4vecz8wKGmzpIUN61wLvD9/pusYCt1n1pKjPZnZiPLw8KqIeM8hbGNPREwcOafZ4ck9ZDMzsy7gHrKZmVkXcA/ZzMysC7hBNjMz6wJukM3MzLqAG2QzM7Mu4AbZzMysC/wfKOxWDwWsb6AAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from io import StringIO\n",
    "import sys\n",
    "\n",
    "x = np.random.randn(100)\n",
    "\n",
    "old_stdout = sys.stdout\n",
    "mystdout = sys.stdout = StringIO()\n",
    "\n",
    "with pm.Model() as model:\n",
    "    mu = pm.Normal('mu', mu=0, sigma=1)\n",
    "    sd = pm.Normal('sd', mu=0, sigma=1)\n",
    "\n",
    "    mu_print = tt.printing.Print('mu')(mu)\n",
    "    sd_print = tt.printing.Print('sd')(sd)\n",
    "\n",
    "    obs = pm.Normal('obs', mu=mu_print, sigma=sd_print, observed=x)\n",
    "    step = pm.Metropolis()\n",
    "    trace = pm.sample(5, step, tune=0, chains=1, progressbar=False) # Make sure not to draw too many samples\n",
    "\n",
    "sys.stdout = old_stdout\n",
    "\n",
    "output = mystdout.getvalue().split('\\n')\n",
    "mulines = [s for s in output if 'mu' in s]\n",
    "\n",
    "muvals = [line.split()[-1] for line in mulines]\n",
    "plt.plot(np.arange(0, len(muvals)), muvals)\n",
    "plt.xlabel('proposal iteration')\n",
    "plt.ylabel('mu value');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 0., 0., 0., 0.])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trace['mu']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that for each iteration, 3 values were printed and recorded. The printed values are the original value (last sample), the proposed value and the accepted value. Plus the starting value in the very beginning, we recorded in total `1+3*5=16` value above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "seaborn 0.10.1\n",
      "pandas  1.0.4\n",
      "pymc3   3.9.0\n",
      "numpy   1.18.5\n",
      "last updated: Fri Jun 12 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.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}