SI_7_parameter_recovery.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/felixmolter/anaconda3/lib/python3.6/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  from ._conv import register_converters as _register_converters\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import os\n",
    "from bambi import Model\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.patches as mpatches\n",
    "%matplotlib inline\n",
    "\n",
    "from plotting_functions import plot_correlation, cm2inch, show_img\n",
    "from pymc3 import plot_posterior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load parameter recovery results:\n",
    "recovery_dir = os.path.join('results', 'parameter_recovery')\n",
    "files = [file for file in os.listdir(recovery_dir)\n",
    "         if file.endswith('.csv')]\n",
    "\n",
    "results = []\n",
    "\n",
    "parameters = ['v', 'gamma', 's', 'tau']\n",
    "\n",
    "dataset = np.array(39 * ['krajbich2010'] + 30 * ['krajbich2011'] + 24 * ['folke2016'] + 25 * ['tavares2017'])\n",
    "\n",
    "for file in files:\n",
    "    _, _, subject, _, _ = file.split('_')\n",
    "    subject = int(subject)\n",
    "    result = pd.read_csv(os.path.join(recovery_dir, file), index_col=0)\n",
    "    results.append(pd.DataFrame(dict(subject=np.ones(len(parameters)).astype(int) * subject,\n",
    "                                     dataset=[dataset[subject]] * len(parameters),\n",
    "                                     parameter=parameters,\n",
    "                                     generating=result.loc[[parameter + '__0_0'\n",
    "                                                            for parameter in parameters],\n",
    "                                                           'generating'],\n",
    "                                     recovered=result.loc[[parameter + '__0_0'\n",
    "                                                           for parameter in parameters],\n",
    "                                                          'MAP'])))\n",
    "results = pd.concat(results).sort_values(['subject', 'parameter']).reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>subject</th>\n",
       "      <th>dataset</th>\n",
       "      <th>parameter</th>\n",
       "      <th>generating</th>\n",
       "      <th>recovered</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>krajbich2010</td>\n",
       "      <td>gamma</td>\n",
       "      <td>0.790000</td>\n",
       "      <td>0.980000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0</td>\n",
       "      <td>krajbich2010</td>\n",
       "      <td>s</td>\n",
       "      <td>0.010750</td>\n",
       "      <td>0.010854</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0</td>\n",
       "      <td>krajbich2010</td>\n",
       "      <td>tau</td>\n",
       "      <td>0.490000</td>\n",
       "      <td>0.810000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0</td>\n",
       "      <td>krajbich2010</td>\n",
       "      <td>v</td>\n",
       "      <td>0.000084</td>\n",
       "      <td>0.000074</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>krajbich2010</td>\n",
       "      <td>gamma</td>\n",
       "      <td>-0.130000</td>\n",
       "      <td>0.280000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   subject       dataset parameter  generating  recovered\n",
       "0        0  krajbich2010     gamma    0.790000   0.980000\n",
       "1        0  krajbich2010         s    0.010750   0.010854\n",
       "2        0  krajbich2010       tau    0.490000   0.810000\n",
       "3        0  krajbich2010         v    0.000084   0.000074\n",
       "4        1  krajbich2010     gamma   -0.130000   0.280000"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def figure_si_6(results, figsize=cm2inch(18, 5), colors=None, fontsize=7, plot_legend=False):\n",
    "\n",
    "    datasets = ['krajbich2010', 'krajbich2011', 'folke2016', 'tavares2017']\n",
    "\n",
    "    if colors is None:\n",
    "        colors = np.array(['C{}'.format(i)\n",
    "                           for i in range(results['dataset'].unique().size)])\n",
    "\n",
    "    parameters = ['v', 'gamma', 's', 'tau']\n",
    "    parameter_labels = [r'$v$', r'$\\gamma$', r'$\\sigma$', r'$\\tau$']\n",
    "    parameter_limits = [(0, 0.000225), (-1.5, 1.25), (0, 0.0225), (0, 5.5)]\n",
    "    parameter_ticks = [[0, 0.0001, 0.0002], [-1, 0, 1], [0, 0.01, 0.02], [0, 1, 2, 3, 4, 5]]\n",
    "    fig, axs = plt.subplots(1, 4, figsize=figsize)\n",
    "\n",
    "    for p, (parameter, label, limits, ticks) in enumerate(zip(parameters, \n",
    "                                                              parameter_labels,\n",
    "                                                              parameter_limits,\n",
    "                                                              parameter_ticks)):\n",
    "        results_parameter = results.loc[results['parameter'] == parameter]\n",
    "        color_idx = pd.Categorical(results_parameter['dataset'], categories=datasets, ordered=True).codes.astype(int)\n",
    "        axs[p].scatter(results_parameter['generating'],\n",
    "                       results_parameter['recovered'],\n",
    "                       marker='o',\n",
    "                       color='none',\n",
    "                       edgecolor=colors[color_idx],\n",
    "                       linewidth=0.5,\n",
    "                       s=30)\n",
    "        axs[p].scatter(results_parameter['generating'],\n",
    "                       results_parameter['recovered'],\n",
    "                       marker='o',\n",
    "                       color=colors[color_idx],\n",
    "                       alpha=0.5,\n",
    "                       linewidth=0,\n",
    "                       s=30)\n",
    "        axs[p].set_xlim(limits)\n",
    "        axs[p].set_ylim(limits)\n",
    "        axs[p].set_title(label, fontsize=fontsize)\n",
    "        axs[p].set_xlabel('Generating', fontsize=fontsize)\n",
    "        axs[p].set_ylabel('Recovered', fontsize=fontsize)\n",
    "        axs[p].set_xticks(ticks)\n",
    "        axs[p].set_yticks(ticks)\n",
    "        axs[p].set_xticklabels(ticks, fontsize=fontsize)\n",
    "        axs[p].set_yticklabels(ticks, fontsize=fontsize)\n",
    "\n",
    "    if plot_legend:\n",
    "        patches = [mpatches.Rectangle((0, 0), 0.5, 1, fc=\"C{}\".format(i))\n",
    "                   for i in range(len(datasets))]\n",
    "        axs[0].legend(patches,\n",
    "                      ['Krajbich 2010', 'Krajbich 2011', 'Folke 2016', 'Tavares 2017'],\n",
    "                      loc='lower right',\n",
    "                      handlelength=1, handleheight=1,\n",
    "                      fontsize=5,\n",
    "                      frameon=False)\n",
    "\n",
    "    # Labels\n",
    "    for label, ax in zip(list('abcd'), axs.ravel()):\n",
    "        ax.text(-0.3, 1.1, label, transform=ax.transAxes,\n",
    "                fontsize=fontsize, fontweight='bold', va='top')\n",
    "        ax.spines['top'].set_visible(False)\n",
    "        ax.spines['right'].set_visible(False)\n",
    "\n",
    "        # Font sizes\n",
    "        ax.tick_params(axis='both', labelsize=fontsize)\n",
    "\n",
    "        # Plot diagonal\n",
    "        ax.plot(ax.get_xlim(), ax.get_xlim(), color='black', alpha=1, linewidth=1, zorder=-1)\n",
    "        \n",
    "    fig.tight_layout()\n",
    "\n",
    "    return fig, axs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfcAAACGCAYAAADAfMP3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXecFdXZx79n5va6vbK79LIIUuwSEVtEEtSosbdoYtRYMK+x5DUxUWNeW9AYYzSJsUbUGCtGDEWaKIJI77BsY9l69/Yyc94/dkHpC+zubJnv53Ph3rln7vzOnufMc2bmnOcRUkpMTExMTExMeg6K0QJMTExMTExM2hfTuZuYmJiYmPQwTOduYmJiYmLSwzCdu4mJiYmJSQ/DdO4mJiYmJiY9DNO5m5iYmJiY9DBM525iYmJiYtLDMJ27iYmJiYlJD8N07h2AEOImIcSpre//LoTwGCzJpAMQQkwQQvy29f1vhBDDjNZk0jEIIVxCiOeEEI8LIbYLIS4zWpNJ+yOEOEEI8YoQYnPr/8VGazpcLEYL6KGsAEYKIXTgCyllyGhBJu2PlHK2EOIKIcRRQFxKucZoTSYdxk3AW1LKGUKIEinla0YLMml/pJSLhBA7gGullPcZredIMK/cO4YVQClwLfBXg7UYghDiGiHEdqN1dBJTgMeMFmHSoRwFLBRCOIGo0WKMRAixSAhxv9E6OpARwHKjRRwpnerchRDnCSFWCyEqhBB/6cxjdyZSyiZgPPCelDJltB4DcQohPhZCVAshjjdaTAdRDbwmpUwYLcQIhBClQoglrbeq/89oPR3I28AzwBPAgwZr6XSEEAOEEBuEEJ8DeUbr6WAagR8KIbr1xW9n35aXgAoI4CdCiF9KKes6WUOnIKU8ymgNXQAXcDHwPHAncKGxcjqEPGCO0SIM5FfANinlWCGE12gxHYWU8j3gPaN1GMiNQBNwMlBusJYORUo5F5hrtI4jpbNHJj8BPgWuaP1sTjTr2SSBFBCjZWDX45BSXi+l1IzWYSAC0Fvf6wcqaNLt2dmP40YLMTk4ne3c/wX8ELi1k49rYgxJ4N/AmbTczjTpeTwA9G+dX/Fro8WYdBh/AfKBxZhztboFwsznbmJiYmJi0rMwR2AmJiYmJiY9DNO5m5iYmJiY9DB6hHM/++yzJS0TPczXfl4NDQ2ytLRUTp06tVs/hzHb+uCv9957TxYUFJjt3MNfmqbJCy+8UF566aVmW/fw17p162ReXt4htXOPiFBXV9cjV9O1G7FYjPPOO4+JEydy++23Gy3niDDb+sAsWrSI6667junTpxst5Ygw2/nASCmZMmUKDQ0NfPTRR0bLOSLMtj4w27dv5+yzz+bhhx8+pP16hHM32T+6rnPllVdSUFDAo48+arQckw5k/fr1nHfeebz44osce+yxRssx6UAee+wxZs+ezbx587Db7UbLMekggsEg55xzDtdeey3XXnvtIe1rOvcejJSSO+64g9raWj7++GMUpUc8hTHZBzU1NUycOJGHHnqIc845x2g5Jh3Ia6+9xtNPP82CBQtIS0szWo5JB5FMJrnooos45phjuO++Qw9zbzr3HswTTzzBf//7X+bPn2+O7nswoVCISZMmcdVVV3HdddcZLcekA5k1axZTpkxh1qxZ9OnTx2g5Jh2ElJIf//jHWCwWnnnmGYQQh/wbpnPvofzzn//kySef7BajeyFEf+CXgF9K2RND1LYZTZcIQFFaOvNnm+p49fNtJDWJzSK45sS+jO2bsav8ztH96NGj+dWvfmWQapPOYPny5VxyySW8+eabDB8+3Gg5B6Q15fUDwCrgdSnlHEMFdTE2NGxgbeNaxuaOpcBTsNf39913H2vWrGHWrFlYLIfnpk3n3gOZPXs2t912GzNnzqSoqMhoOQdFSrkZuE4I8ZbRWowimtB48MPV7AjGQUJemoPrx/Xjr/O2cFShj5Yor5JHPl7L2JJ0RhSmcfqwHH56ww0oisKf//znwxrdm3QPtm3bxqRJk/jjH//I+PHjjZbTFiQQAhxAhcFaugy61Hlw0YOEEiG8Ni+zt81mSMYQbjj6hl1l/vKXvzBt2jQWLFiA2+0+7GO1m3MXQrhpyZqUAOZIKV9t3X4UcE9rsYellCuFEL+nJalIREp5955laAlbehcthvG1lLInZ5tqV1asWMHFF1/MtGnTGDFihNFy2gUhxE9oyUtAcXGxwWo6ht9NX4PTqlKa7wMgHE9x++tfcWy/TEAgpeSdZZXsCMZZtz2Ex1HF/a/OJrVmI59+Mv2wR/dHSlv7PWafPmwaGxuZOHEiU6ZM4eKLLzZaTluZJ6X8VAiRS0vo6cu//WVv6NP74pOtn5DSU/T19wUg05nJ+sb1bG7aTP+0/rz33nvcf//9zJs3j5ycnCM6VnvOsPoB8JaU8sfA5G9tvw24GbgJuEUIUQxYpZS3AqoQomjPMlLKdVLKH0kpLwN6aqrQdqe8vJxJkybx1FNPMWHCBKPltBtSyueklMdIKY/Jzs42Wk67I6Vke3MUt/0bB+22WwhEU7s66NJtjTSEE6gCUppOPBikJhjnzNsfx+MxNP9Sm/q92acPj53LWL/73e9yxx13GC2nzUgpdyYRagT2mvDT0/v0/vhyx5fkunJ325brymVB5YJdy1jfffddBg4ceMTHak/n3odvUgF+O0uWX0rZJKUMAF6g8FvltrXut2cZAIQQlwAz9nUwIcRPhBBfCiG+rK2tbcdqdE+ampqYOHEit912G5dcconRcg4JIUSmEOJZYLQQ4p6D7tADEex9S70ow8mm2hAJTWfp1gbCcY1IQieVSlETjJCWnUdtzACxu9PWfg+YffpQ0HWdq666iry8PB577DGj5RwSQogfCCH+ArwMPG20nq5CobuQYCK427bmRDO2sI3zzz+ff/zjHxx33HHtcqz2dO4VtHT0PX83IITwCyF8QBCo/Fa5otb99iyz8yRQIqV8dl8H660jv30Rj8c577zzOOOMM7rV6H4nUsp6KeVPpZQDpJSHFqmhByCEoG+Wm+Zocte2pkiCowr8XHJcMe8uqySutQSnsioQjcXx+vwkdYFdMfw5e1v7vdmnD5Gf//zn1NTU8OKLL3a7ZaxSyrellDdIKS82J9N9w0VDLqI8VE44GQagKdZEU6iJey+9lwceeIBJkya127Ha02LeBi4QQvwZeF8I8XLr9ieBp2gZvf1RSrkNSAohngA0KWX5nmWEEKOBx4B+reV6LeF4ir/O28x9767kg6+r0PXdIxDuHN1nZ2fzxBNPmJOquil3fncIdqvC+pog62uCeB0Wppw5mNHF6RzdJ418vxObCvF4DKfDTkoXpHTJuaMLjZbepn5v9ulD44knnuCTTz7hnXfeweFwGC3HpJ1wW908/J2H8dg81EZrybXn8tmvP+OKK67g+uuvb9djtdssHCllGPh2CJ1XW7evBK7eo+w9e3zeqwzfXA30WoKxJFOmLaMo3UWGy8q8DXXM31jH7y8YuavMnXfeSXV1NTNmzOh2o3uTb7CqCndPHLbXdq11MJfmUKjYVoPDn40UCgjBcX3T+c5gY69wD6XfY/bpNjFt2jT+8Ic/sHDhQtLT042WY9LO+Gw+bjz6RpLJJOeeey6jjhrF/fff3+7HMZfCdWFe+qyM4gwXXocVgDy/g631YVZWBjiq0M/UqVP56KOPmD9/vjm678ZUrlnM6wvXsMU6CGGx88NjinY57WyvnfToVuxbP6bAn09NZi5xixtNh3NG5uOymV24JzFnzhxuueUW/vvf/3aLZawmh4eUkp/+9KcAPPvssx1yx9U8M3RhqpqiZHl2n2jqc1hZXRVg9YKPefzxx1mwYAEZGRn7+QWTLk0ySuqdm1m1OcqJHg/Dgp+xKfdsXv4sxYDAQtx6gLeWNTB466cEXR6cjkpE4DWqfKO45rqbKcpwGV0Dk3ZkxYoV/PCHP2TatGmMHDny4DuYHJxYMzRsgowB4PAZrWYXv/nNb/j666+ZM2cOVqu1Q45hOvcuzJA8LxtqgmS4v3HwDeE4csdGbvnZz/jkk0961RrR7oSmS2auqWHjjhAThuYwLH8fJ5b5f2Bb3Ivw+hBWiUfRKamaTo5wUr6wlsaEhT6BbVRJH1lpbpKaZGB+HiWNy0zH3sOoqKhg0qRJPPnkkz1qGauhfPYnKPsMbC5IRqDoBDjpZ0ar4vnnn+fll19m4cKFHbqM1XTuXZhLji3m528sI5LU8Dus1DTH8RHl1h9dyj//+U+OPvpooyWa7INoQmPKtGX4nFb8TgtPz9rIgBw3d5w5ZPeCTdtI4CGt8WvUVAirVFA1yXZ7PzZZvBCpp0hGGOyEKqUATYdQIkVSKuiahqKqxlTQpF3ZuYz1lltu4dJLLzVaTs+gZjWUfwG53wrTW7EYtq+APOOCe3344Yfcd999zJ07l9zc3IPvcASYzr0LY7MoPHXpaOZtrGVDTYjTSuxc/v3zeeKJJzj99NONlmeyH15atJU8vwO/s+V228AcD+trglQ0RuiT/q0rbtXGwNBiViRjWBMS6mN4kzHs3jg2n0ZNVCHl9pOhJNBSSSIJqGjQcKkW07H3EOLxOOeffz6nnXYa//M//2O0nJ7Dyrcgvd/u2zL6w4p/GebcFy9ezDXXXMP777/P4MGDv/ki2gQr3255P+ICcPjb5Xjm9OoujhCCUwblcOGITG664gfcfPPNXH755Qff0cQwttZFdjn2naS7rHy1rWn3gsf9BNFcjSuRhI1xqNMRTTqZFc1YayKUZDhIV+NUJD00RTV0XaNWc/OB/Xu8sGBLJ9bIpCPQdZ2rr76azMxMcxlre+MrhERo922JEPj3TtLSGWzcuJHJkyfzt7/9jRNOOOGbL8q/gLd/DDtWt7z+dT1UfNkuxzSdezdg5+h+/Pjx3HnnnUbLMTkIuV47oXhqt23NsRSl+d7dC+YdhcwdTrRGQSgOUCyg2tF0ia1ZkJWsIVevYYk2mNflGXysjKd60GWMHTmCRZvrSaR0TLovv/jFL6isrOSVV15BNe/EtC9HXwJNZaAlWj5rCWjcCkdftlfRQDzA4u2LqY10TFTEHTt2MHHiRO6//34mT568+5efPQMFo8GV2fIqGN0yV6AdMG/Ld3F0Xefaa68lLS2NqVOnmqP7bsC14/pxx7Sv6JPuwu+0UtkUJdNtY0COl6qmKBIoTHOiNTfTWJaBVl9GvZqGhyCpZATd6iNfayTNrqDZ3JwS3UqRN5OmYd+cmCyKQmMkQa7PXALZHXnyySeZPn26uYy1o7C54ftPwfw/QLwabF743pNg330C2yurX2Hx9sW4bW5iyRhFviJuH3N7u51nw+Ew3/ve97j44ou54YYb9i6gJ2C30NPimwHJEWI69y7O3XffzbZt2/jkk0/M0X03we+08sfLxvDG4nK21Ef43sh8ju6Txr0vfERe0xKU5iDFa9ZSWrmFspz+hFN9SavfTrXbTaZTkG+J4cp0ouT2RQluJ+ErxhMqY/qqTeTn5tIvy01K18n27JWPw6Qb8NZbb/Hoo4/2+GWsrRkD5wK/llJ+0OkCfPlwziP7/boiWMGXNV8yJGPIbtsWVC1gXOG4Iz58KpXi4osvprS0lAceeGDfhZR9uOB9bTsMTOfehZDJJDumPkli8wZIRVndHGP618v5dMECnE6n0fJMDgGXzcI1J38zoeePf/8HpzYvwBVoJv2rVVgakmhRKKxZR11mFht1L86Qyvb0IrLUJXgG2CFUSyz/WJrLJLGkRlzVWFkZYNm2Ju77fimK8XHlTQ6RuXPnctNNNzFjxgxKSkqMltPR3AW8YbSI/TF722zy3fm7bSv0FDK3Yu4RO3cpJTfddBPJZJLnn39+/3cCRl4My16FrCGAhNr1MOaqIzr2Tkzn3oXY8cQfkJvn4hARmiJh0muref+Mk8nMzDRamskRkls3nzJnI2mbNmINJHEmJFIK9CjYQ3WkZeWT8JQQOTqT2lQ9w7KtUHQsX26sI9cHzaqkT2YWcV1BUcS+182bdGlWrVrFRRddxGuvvcaoUaOMltOhCCHOAFYD+3zm0BXyuRd4C9jQtAGv7Zu5MOFkmGznkYd0fvDBB/nyyy/59NNPDxykZuikFse+9MWWz2f+FrKOPN0rmM69S5FY/RkOr05TzMbqzRsYOXIkzoYNEAu02/IIE2PYIsppiiU4pT5BXEDYpuDTQCCJJQXSEabAuoWk14XTNhgaZkG0ESVSSKPwsznvuwzKarmFW9UUZUNNaPdldSZdmsrKSs455xwef/xxzjjjDKPldAYTADdQCkSFENO/leMdKeVzwHMAxxxzjNz3T3Qs4/uM592N7+Kz+XBZXSS0BFubt3LD0ft4Nn4IvPDCC7zwwgssXLgQr9d78B2yBsJZ+7ltfwSYzr0LoSbqCKfSWL1qBUOHDcPj8SCTDbBpNgw/z2h5JodBJJHioZfnk742xrBYEFdYosQEURdYbOBKQMQG1e4UwSwn+bpkTL8sEBdD5RK2Zl/ARnUAXtc3j2WaYylG9DEHe92FQCDAOeecw4033sgVV1xhtJxOQUr5SwAhxDVA3bcde1fBolh4aNxDvLDyBarD1XisHn594q/x2w/et+qidcytmEuuK5cTC05EES0Lzz766CPuuecePv30U/Ly8jq6CgfkgM5dCHEzsNuoSkr5TIcq6sXoPjvrVy1nQP8BZKSnIxMJrOkW8BhrJG2lt9mLlJLVVc3UhxIc1z8Dh3XvCY+Pz1hP+uqvsGpeUqKRiizoUwWuuKDJBTvSdRaOEWzJVRgk4hRYw6g7n89Z7Jw7aTJ3vLEMf1LgdahUNcUYXZy2V86B9qS3tWN7Eln6FY3TXodkCuHzkvaTn/CDK69k3Lhx3HXXXUbLOyIOxy6klP/oSE1Hitfm5dYxtx7SPu9sfIdZ22aR68plac1S3lj3Bg+Oe5B1y9dx1VVX8e677zJkyJCD/1AHc7Ar95XAVcDHQBxov0zyJmyoaaaiMcbYvumQiHLj8kqm9nfidYIM7sCWYcU3sgSKjzdaalvpNfYSjqe461/LUYTAZlF48bOtXHliCacOydmt3PZAjKNjITRLFlG9nKQHyvpIrFFYOUCwtY9CzGtFah5Uh4214XLOYADoKbC5cdstPHP5WGauraGsLswVx5fQL7vj4lG30mvasT2Jb91K/V//in3YUAQCPZFgxvnnk5afz1NPPdUTlrH2eruIJCPM3jaboRlDAUgnnUgywuPzHue5K5/jueee46STTjJYZQsHdO5Syk+FEOdJKd8AEELsnXDa5JBJajr3vr2CaFLDZVN5ddEWVkx/mbGjJ9D/VzfD53+BVBzhz4fv/By6yUmhN9nLH2dtIM/vwN2acrUwzcmri8oYNzALi/pNbChbLIS1cT3Dg1uo9ySpsAmCbogWaHw9XMUiFUoSEqdVRYZ2YBUSGre0ZLOa2LKMR1UEZ5V23t2b3tSOh8vGHUH+Nn8L8aRGltfBzRMGEnv9dewDByJa1y3PmvspwViM5266qUcsYzXtAjY1bcJtde++MQnTPpzGL3/5S84//3xjhO2DtjxzLxNC/AtIAAs6WE+v4MUFW/BuXEX/qi1IIShrChPNHsbvHrke4bbDuU8ZLfFI6BX2UhOIUZy5eyd32S1srgszONdL5OuvWfn4/ZxRsZ64S0CjJLdOQ08T1ORINhQrqBL6xVVSMpOGpBWH1csQmwbH/hgKx4JiaADJXtGOh0N1IMrvPlxDaYEfVRFEEynufHMZ98fiJDdtIllZSWNDI/6mRo476ywssbjRktuTXm0XBZ4CYqnYrs/JZJLXXn+NAf0HcPPNNxuobG/a4tz/CiyQUi4WQhR2tKAeR2MZfPEXiIdabq+PvJTYu29zVKgOq8vF5s2bGRVoZsD4gXxdEWD8Hrd1uyG9wl4sqoKUcrdbrYlUS2CZ2ObNrJ76ILXl6/FGdHIaIGmBgFvgQsOlgE21MDAmEJqHlMWGT9cYa3Mw3pMGRccaWLNd9Ip2PBxeXLiVwXle1NY4A06bBZ/TRnlUx//lEpKAFg7TLzub1JKluB96yFjB7UuvtotMZyZFviIqQ5XkOfN48503cRQ7+NO17RMytj1py6XBVOCa1vf3dpyUHkjdBph+Jyg28OZD2SLkf+6mqGIDisNJZWUldXV1DDlqOAVb1tKnZ+To7hX2cvnxxayubkbTWyYBVzVFqW2O8X8fr2Xez39FbNUacup1XFFQNbAlwBYTxIVCelJhYL0Fm3Qw3NKHCx1DON89nFx9EEI9wJrYzqVXtOPhEIylsFt2P3XaVYGyaAEpp5NQOIzP70dJJBB2O6mqaoOUdgi93i5uH3M7pxWdxnv/fY9UPMWbP36TEn/XC0jUliv3INDY+j52oIImezDztxCohEAlutVHqMJCorqWggo7/40oBOMJ+g0aRkoK8mJNDOj4iVKdQa+wl1HF6dz53SG8vGgbsaTGqqoAJ/TLxF25hUh9I564ilXTEHrLlImU2uLkg6qCQ0gCqpcqJYNTdS82xUosqVGiV8KIa42u2k56RTseDuMHZ/Pxqu2tcQYkcuF8SlYvw9VQS1M8hrtvX5zp6QirlWRlpdFy25tebxdCCD5/9XNq36hl7ty5+HxdM6BUW5x7LXCKEOIxQOtgPT2HyiVQuRTSimmOJWn+fDORgIOYPZ1Km4cBVWuo7D+SppjEocY5bXC3vx2/k15hL0k9yabIfLKLNpAIDcamFuJauZjcmdPxN9Xg0DRUnZaFQxJQQbNA2K0QVxXcLg238wfMUXPRKpbi1EMUDL+GwgFn7zukV+fTK9rxcDizNJevy5tYVxMk/4s5jF0+hxw1SSgRw+N242huRuTmgC6RmoZ98CCjJbcnvd4uXnrpJZ577jkWLlzYZR07HHyduwA+Ad4GFCnl6k5R1RNY8g9IK6Y2ECTQlICQC6lItoatLFP9jPIKimJNaNm5BOIpcn9+j9GKj5jeYi/RVJR75t1Duj2dNHsaC8rf5vgvG5nwdRnWhihqQtuV50kCKQVq/VBWYqE234lw9ePWs6/k7zsG89bSCgYUfZfSAh81CY3/fWclj110tJHV6zXteLgIIbjz7KFU/OMlAp+/hSUvj5XrN9Hf7sCmqshkkmT1doSq4r/ggp6wBA4w7QJgxowZ3HnnncyZM4eCAmNyw7eVgy2Fk0KICVLK/afWMdk3WhJZfAI7Fn2MU0uR1DQimsrqVB+kx8MGRx7JjCxsp5/JxqjgJ+npRis+YnqLvby57k1yXbm7IlmdpJVRWLMDS1hHSWq7JrLsjPaRssCWEjcrT8oiLSVpdl+OY/g5VKxfwXmjCned/L0OhepAjKdnbqAqECXba+eak/qR7rZ1av16SzseCQ2vvErk739DlTqB9esosVhIGz6c+KaNKE439v79cR13LJnXXWe01Hajp9hFY6yRUDJEH0+fQxp4LV26lMsvv5x///vfDBu29yrA2dtm80nZJ0gkQzOGcmXplbsi1xlBW27LTxZCnA4009K+P+xgTT2D7KEkG8spc4zAZw0RVyupj9uIKi6sFgs2BT7LHMjxbg9+2aOWyhyyvbSmhnyGluU1c6SUr3awxiOirLmMTGdLMh9NS+Ighr85jhoRiD2CbOoKRO0QyRCM1j0UJJ3EcyUWVUHT5V4nl5WVASQwJNdLOJHijjeW8cQPR3W6g8fs9/slVVdHeNEiVI+HmkAAuxD4/H60hgYsGRk4Bg+hz1NPGi2zo2iTXbSugb8NyAJmSin/3Ika90lSS/Lw5w/TFG9CVVQSWoKfjfkZg9IO/thk69atfP/73+fZZ59l3Li9M8b9Z/N/mF0xm2JvMUIIyoPlTF0ylTuOuaMjqtIm2uLcTwEKpZTlQggzgfRBCMWSRBIaGWN/QuOt3ye7KYZU7WxN+linFJLusFOfhDU5Q6h2prO1LsyD5x1ltOz25HDs5QfAW1LK94UQ04Au7dyLfcVUhapY07CGQKCctPgO+lo0VCzoAqQEQcuVe8ICzR5Jji1FvtqHEb5t5EycCMCgXC9VTVHSXS2Oe3sgSkLTGZzTMrHSbbMwIMfDPxZuYcqZnR7O0uz3+yC2cSNbr78e2RQglkoR0HUGFxYiQyFS9fWoaWnkP/Bbo2V2JG2yCynlGuCnQggFeL7T1B2A55c/j021MTC9JeuaLnWe/upppp469YBX8PX19Zx99tncfffdXHDBBfssM6tiFn19fXd9znBksK5xHdFUFKfFmHTdbXHuU2nJ7nMd8BRwZClzeii1wTi3vLaULXVhJJKTtn3OjyOCbJ+Nz2ttbNDzSKZn8XH/YzlhWD597Q7sgSh/vmKs0dLbm8Oxlz7Aitb3e03SMTo9pC51NjRuwGVxUeQr4odDfshlH1yGLjXS4hE0i52gN0LSIhFSoGigC0gJ2JAPFYU2MtPzUUNVOCff0pIxBrh+XD9+88Fq1tcEsVkUNtaEGNUnbbcTjdtmoTZoyKRks9+3oumSfyzcQtN7HzD+w7/h0hJIwColAxQF1WJBpqdDKkXRX55F9ffopD5ttgshxGTgbuDpfXzX6X16W3Abxb5vjqUIBZtioyJUQZG3aJ/7RKNRJk+ezOTJk7nlllv2+9tS7p3YThUqsVSsSzt3DShrfR/oQC3dmv99ZwXbm2P0zXIj4s2M1VdQnVRwRUPkRSoYmKfxZVTF1deH8LjZVhfi1gk9ahbtTg7HXipocfDL2EfsBSPTQ25o2sBTS57CZXWhSQ1VqNzpHsngigbqSeKTzRTWaWQ2COI2sCqSkCqIK1CZoTK/1MHG3Fwys89gdO53+Whlkv9rCUuNRVV44NyjqA/FaYwkcNpUHvxgzW7Hrw3GjcoAZ/b7Vu5/byWEQoyb+SbOVscuAKV1ECZVFdXvo8/UqVizjzwXeBenzXYhpXwPeE8I8SHw2h7fdXqfFkLsFXhKkxoe676XIGuaxuWXX05JSQm///3vD/jb2c5sIskILmtLrBJd6iAh3WHcXKq2OPc4MEwI8TOg+8/66gCaIgm2NURIc1qxRms5rulD+tqrIJWiMZJkVOkIbIkAfleKGS4ruG3ceOrRHZrZy0AOx17eBp4WQkwC3u8wZYf9caYmAAAgAElEQVTBn776E0MyhuyaGBOtW8+0Jc+RrXiwKBXUKzrKdpXydHBENXKCEFFUNpfAnOMEZb58ijOPoT5ci9Wi0BhJ7nWMTI+dzFZbGFWcxldlTWR5bATjKayqwkVj931V0cH02n4vpWTptibW1zQzOMfD9h07mPjeIxRG6tjr5q2U6JEwJX//W29w7NBGuxBCnErL4zY7ML1zpB2Ys0rO4qMtH1HiK0EIQXOiGb/Nv08HLKXktttuo6mpiY8++gjlIKGgfzb6Zzzw2QNUh6tbrti1GFPGTumoqrSJtjj3/wXOoOWKyvBJEV0Ri6qgIFCDNZwaeJ/+1kp86VFEbZQ8n4rNakWPSFxuhZsuPtlouR3NIduLlDIMdJnoLTtpiDWgChVFKCS0BCk9RbhsDTtSTo6NxCjLjBFB4I5JrCmBLgTSKnEKnSypM0hVKUiFqY42AXZ0qaMoB56de8MpA6gORFm4qY5BOV5G9knrnMruTa/r9ws21vGvJRUs3tqA32nh6KJ0/vLeUr63+m0GlpUj9nN9mX7llVhzekycioPRJruQUs4B5nSOpLYxoXgCADPKZqBLnWJvMT8f+/N9ln3kkUeYO3cu8+bNw24/+EWYy+ri4VMepjpUTTQVpZ+/n+FLINvi3F8F1gGvSil7ZdCCg6HNnsmdL/4GK0lcxEm6LYi+CZSBNrKbQmiJGqx+Df+DjxsttTPoMfbisrhI6knmVcwjnAyT1CRaMsZxuhu/Uk+aJlFUgeaQFGwHqxTEbAoqYFNVcus0NhbZcEQqUe1jaAinGFl48Fvs+X4nF4wx5Gr92/SYdmwLc9bu4I0lFcSTKfxOK4oimLt0K2etms1pG5ZhSe3Hs9tsZF19deeKNZZubRcTiifscvL749VXX+WZZ55h4cKF+A9x/kS+J/9I5LUrB3XuUsqLhRADgCuEEKdKKQ/8l+llRCurqLr7bnJjsV3ByKzxJEKDfHeEtKEe8Lng5FvA2SPCyx6QnmQvDouDxmgjMS1GpjOTsvowEWuM7Xo9W7YnSS9TER7YnC8YsKWl7RXA5YAci6Q5lMKhxgnGIvRxn8SQPC9Xn9jX4Fq1jZ7Ujm3hnWVVDMpxM2vtDuwWhbK6IAPqVjOmejm2VGqv8hqAaqV00WcIY7P3dSo93S5mzpzJHXfcwaxZsygs7N55cQ7q3IUQw4ErgP7AGx2uqBsQjqe4792VfLGlnpNXzObiRAoHtD6PkyiaQIR1otvdMDYHLp8Gdq+xojuJnmYvhbYcnIs2YalfT98kbPO7KalLousteda9AXBFBIuGSQq3+RiogdMeoVCx0z8tm+/5M3GfejUMOsXoqhwSPa0dD0bLbGeB12FlXXUT34nM45SGxfQJ7tjrObsEVub3J3Dvo4xw9YhkT22mJ9vF119/zaWXXsqbb77J8OHDjZZzxLTltvwFwHNSyi0HKrS/QCRCiKOAnbFVH5ZSrhRC/B5wAREp5d37KfMX4CQp5YhDr1bHkUjp3PjKErbUNNGXasJ2iWhdBtHyb6uLT4LUgImP9BrH3kqb7KWrs6RmCe9tfA/v+zMoCYKlSeIPKYzSgkgpqPMrJNMtWKTEmtCozvQQdA6kNuRmQKqaLBnk6CIP7v6nw6AzjK7O4dAR/b5L9mmA/DQH4XiK0nwfjZuWckbV56SFA6h8czt+57uA1cUX597EBf16xQS6Pema/Xvth7B2OiBh8EQo/f4h7V5WVsakSZN4+umnGT9+fMdo7GTa4tw/Bn4jhJDAn6WUi/ZTbn+BSG4DbqalbzwihHgIsEopbxVCPCqEKNqzDHCDlPIGIcRbh1+1juGFBZtZt62S8alF9BG1pCuNWIS+K3DJrlG+ruC75XHoM8ZAtYbQVnvpssytmMs7G9+hnz0fb0RSr2mUhCDeugTfmZSkhXUCboHXpmHVJS6HwFqqsSExmi/9F1KvW8iXDl475kSDa3PYtGu/pwv3aYDbzxjM7a9/xcYV6zm7fhFFO2qQFgFOEOGWwEQagoRqZUXWAMgv5OQBWUbLNoKu17+/fAHKFkJGv5bPaz+AUDUc95M27d7Q0MDE8yYy7n/HsTRvKSsXrOSK0isYnD64A0V3PG1x7jcC19Pit54B9teY+wtE4pdSNgEIIbxAIVDe+t221v32LHNQjAiCIKXk9S/KGRNdxciG9dgjCTJqm4m6rDgjSYQuUWhdA1tUhO+sszpFVxejrfbSZZm+eToD0waiBesJqS0z4eXO0ZsCkdbJszZNR1cUcKtEPF5O0DwUKquY5SwlW7FQVh8hqelY1W75TLa9+/1BMSxY0drpKKs/5LhZdfTb4ea44FrUkAYKSCvobpAJQTNukg4XQ08azfUXHW34bGiD6Hr9e9NsyPlWrPf0vrB5bpuceywW49zzziX/mnyOO+Y4HBYHutR5aulT3HPcPRR6u+9z97aedSSgw97LPL/FzkAke/5uQAjhF0L4aMkFXPmtckWt++1Z5uCCpHxOSnmMlPKY7E5aX1oViFFk1xm1dRX++hCe5iiOaALdYiHkddJo9xGyOok6PfR/6cVO0dRFaYu9dFk0qSGlpFwPsiUT0oOStBCkN4MtAWU5VirS3JTlevD4wFmaRmmkEItiJUOJIOLNCCHQJVgOsvSti9Oe/f7gBzOgT7PgKVj9Lm+uipIbrOf44GpcsThJoYAGQgqkWyHqtZM5bBBHXXURpzxwF2r3btcjpWv1b7n3hEfaMJFf0zSuvPJK/Ef5OWHMCTgsLcmWFaEwMG0gr699vb2VdiptuXL/M/BC6/u9wgh+i90CkQghXpZSXgk8SUuYQgE8IqXcJoRICiGeAOKtMYp3KwPQevt+tBDiWeA2KY3PrmKtqWLcx69Q0LwDnyVMxGlDFwJUhUa7lyZnGnZFsvK4s1hfnuSSrrMqojNpq710OaSU/HHmBuat9bDQugWro54xKStBZwKbT+KOgqZAVFXYMGwgwWyFCXEnZw/tz7QvthFNaCSxo1vd1AZjDMv3dueru3bt99AF+3T1crQlL/BpzEZjwEdJbQP+aByR0kFA0qqQElakw4rvrFMZ+vNHDZXbReh6/dvmAT0FSqs701Ngcx9wFykld9xxB7W1tdz+4O0sb1i+2/dWxUog0b0DM7bFuatSyisAhBAn7K/QPgKRvNq6fSVw9R5l79nj877K/BL4ZRv0dTwNW5GL/0rDX+fiT/jwyhjeZIy0ZJiQ204iYmOTo5CYMwNbfj5pp01g4aZ6Ljmu8+OgdwHaZC9dkd//Zy3vfVVJhns82+Nf0hyKY22yUe2GOodAKgJF6lRmWtmQY2W4byKW8vdRE0FOGpDFyo1bmJMsZXtEUpzh4pnLu3XegI7o912nT3/5D1j/MfPiTSy1WsiNNZJ0QSopsKYUdKmiImlM87DDlcXFU/7PaMVdha7Xv0+7Dz76RatDFxAPoZ/9O6pDlaTZ03Bb93b0jz/+ODNnzmT+/PkIp2BmxUyynN/MoagOV3Na0WmdWIn2py3O/XxgYev7yXSFZyydScUSmPsISSWXJkXhKH8FoaRKJGlDERKpKSzOGcLSrNEMHN6fjKGDoTWGcS+lW9qLlJL3l1VRnOkCBJHAWAKxUqLBBC7LV6h2GyRjZCKwKgqnxZJEIktZNeRGvutaTF+tioJzLyXXeSI3ex30yzrwlUM3oFu2Y5tIJWDjJ2z3FLJUtaKKFM6UIOoXqJogO6VDShC12GnwZvHFCddzclOMoozetextP3Qpu5BS8uyyBMuSP8cbrCLLrnPCKfm8vewP2FU7KT1FsbeY28fevusu2j//+U+eeuopFixYQFpaSwTIyQMn8/6m93FZXcS1OAP8Aw4a7Kar0xbnntsatEACBR2sp8uRmPEkTUt3UL2jEurjSEcKxSdJRd0omk7Y58BWAmPOOguLqgIQjCUpyey1J4JuaS/xlA67UoJInOEyjkpuJDNDkNXoRgoXPplAU3VEmhMrmfS1NnBx9nI4uSVnsw041sA6tDPdsh3bwpatG2naVEl1cisBv400koTSJc6gQjBbkubVcKZA5lloGHoMen4hgWgSw2MGdg3aZBdCiPOASUAO8Ccp5YyOEPPa59soqw8zNN8H+Iglk/xq9lNcPHrkLmdeG6nlXxv+xYWDL2T27NncdtttzJw5k6Kib1r0zJIzObXPqWxt3kquOxefzdcRcjuVtsaWv7n1/a87UEuXo6quAfmfJTRKB5rVS9ImSTXHsKRDzO8GVRDsn4UjLYtl1UEc1pZVsdluG7dOHHbQ3++hdEt7sakK6R47wViS3GQF4/RPSaWi2DJdlGZnsawiikxX6VPi5vxh43Ba7NgtKjRsNFp6R9Et23G/SAlfvYK2ZR6VqzeToQUoUBzkayliFgjkSFRNx18r8MsUtjyNUE4Oit1GMJZiWH73P9m3E22yCynlO8A7Qoh04DGg3Z27lJLpK6sJx1PUNMcozffREN+BRTgIxlP4HFYAsl3ZfL3ja6wNVqa8NIUfP/djSoaU7PV7VtXKoPSek6mzLc59CBCQUv5OCHEJ36T767FEYjF+9a+l9F/9d47RHTj1OFY9heKw06D4UK06deluIhluNLuDTenj+MPFo9Al2CwKHntb/qw9lm5pL4oiuPKEEl5asJXiwBJSqSSKzYnLZkFk96HU04BLOvjO6G89R481QfreJ4keQrdsx/3y2TPEq1eysEqQKYOoaFhTIU4KpfjKoxNQIZwnKfUn8Cck22UGCauHcELntsmDevvs+G9zqHbxv8Cf9tzYHsseH/5oLWV1EQrSHNSHEkxfWc2xA2zs+URUlzpratfwxvw3OOu0s/AWeLl77t388vhfduulbgejrRHqqlrfHwN07/UBB0JLwSe/Yt2adZxcH2dwpAw9IbGggUXDJaDCmsO/004mmuXCqiqEXIPIUdNJc9mMVt9V6Lb2csmxxVhCNVR/tA2Px4PP7SQc11i/PchIT4TR/XKgqQz8xRCth+aqlsk8PZNu2457oumSxQtn8WWzH7ui49AVHCSxAFHNz4hmHa9oRJUqVlS24SUh7IyxVXDcDT9FOA3LzNcVaZNdiJZ74r8HPpJSLt3z+yPN517eEKGyMcKgXA87gnE8dguqIti8HaRNYLWkACtSSlbWrGTdijWc0W8Aw3OtaMKKw+Lgt5/9lj+e/kdsas88d7dlnXsUQAihAJkdK8dg5jwE21dgCVYxuGoLnm0RLEmNVINCIqISkG42qEWsShvIFtsQ1qmDcHu8DM3vVeFlD0a3tZeKigruuOd+ajylDMlxU5DmYlCul6G5bhw2G/Yf/AmG/wDiAcgdARe+APYemwyo27bjt0lqOj99+UtmNBXisKikgLXOODPT4sz2x9nuiOAngk1aCEs362UfKmQuudY41nN+bzr2vWmrXdxCS2rYC4UQP21vEasqA3gdVkb28ZPmtFIXihOIJqkLJXj5gl+h6Rpbm7eyrXkb6//9JWf6JaVZgnDtWhaueJmKujWUh8q5ffbtbGjc0N7yugRtuXL/ELgVeJfWkVaPJBmDr15DhnYwICxpbvAgFYFw6TSpbkgIImkO3k0fR8rlYWiOl+ZYkoJ0J+MG9cowlPujW9pLU1MTEydO5JSLbsNy9PGsqZ5BQXAjqtCpE1lszb+QUULA4LNaXj2fbtmOe/LOV5U0hBM068UUJL4gzbmeckcYm5BIKVjrSuLVVZzJdOJYcROj2NpMn9NvOOT45L2ENtmFlPIpWuIcdAgji9J49+sqcn12/N5KnLGPKUoEOMMeJ33ZCRTkFOJQHcz800z+xxdmY0Y+wu5jY/MW7M4MZLgBJaMfpZmlPL/8eR4Z/0hHSTWMAzp3IcQvgFG03Ho5G+i5C7eXvYYWqiNSbSFQ5iQZtKBYJKpPw2bTiNgczPaNosmRjk0IYikNXUp8DgvfGdgrE0jsRZexl+qvYet8KBgLxcfDQQLJxONxzj//fCZMmMBdt1zGbz9YjX/A92nWNZAa5YEU547sURPGD0iXacd2YNamL1kVXEa+t4oq61oaHGF8UiczJZFCYJEaO2ySkgSsln3JcFk5YbAFTrzJaOldjq5kFwVpTgblevjPlunkaPNJtwYpt0b5rz3Ogs3/Qm8ayIINISr6VmAtKeJkVw6zozWE0UGmEHqCMRlDUYRCTIsZVY0O5WBX7sOklJcJIeYDZ0gpe+ZfAUhtXEjzJpXmMg9Sh1RMRbHqyGZBKkOwUSlgif8kVF0wINtFpsfOzacO4ph+GUZL70oYby//uRdCNeDJaXHyX70Ek5+G/eTc1nWda665hoyMDP7whz+gqirjB2czZ10tbruFWFKjf5abUwb3qgGc8e3YDvxpwe9YEpyGO1OlUY3zL8CKA7cuOTYao09Koggdn66jCp0hooqYuxRsqtHSuypdyi6uHpfF6sgKRNUatitgU1W+0CV+BH3LV7NFaJSUlvCnYCU/idv5Udoo/lQ9m6TVTrbFgupseapgUXrmBOiD1SpLCHEOUAecJlqCs0zvBF2dysrNNax7aQdHB33IGCgWiWLTkEkFXYWkbqGx7wDGH11CVUTw2EVH421dZmGyG8baS/liCO2AzIEtn+0+aK6ANe/D8HP3uctdd91FeXk5n3zyCWprnILLji9h8qhC1tcEKc5wkeWxd1YNugrdvt/P2fg+r615iZhqRQqNnWHQ41KiKYLVdhuFqQhWXaFvXBAUFvLUMNnaV3DCPwzV3oXpUnaxPbKdpvovabYJdEVF0TQUqRECKtUUo3LzsbqzCcgkX0Xq8Qe3cH76CD4MboD0fmgCtjZt4vSi042qQodysAl1bwHZwL9b/+9xly/JRJJZt96Mp7GGcNBGKq6QiikIKVCdOpasFKGRRYwsVRi65UUKPKrp2PePsfayZS749gjo7yuEsvn7LP7kk0/ywQcf8N577+F0Onf7zmO3MKY4vTc6djC6HduB/5t7LzFFRYo9JmILgQCsUjIwmeQHwSgRmUG9yCAsnGSOOBtyhhiiuRtgmF2E4im21oVJafo322rXktQ1NKGg6RJdQFwoqDoErHasQodIPUmLHW//M5mlBTjK6ueckdcSyBlCIB7g6tKrmTRgUmdVo1M54JW7lLJHpzZrjiWZOuV8RiaqSEvGQJdougK6JCnAateozsyi3F1EekxnqCfKhSVrgNFGS++SGG4veSNg9Tvg/1Z0wHAd5I3cq+hbb73Fo48+yoIFC8jIMB+tfBvD2/EICMQD3D3nN4SlRvIAa9NLkklGxZJ4pBW3Cp78UkS4DMZc0YlquxdG2IVc8RbL57xFRdTKdvy86/WSU6QzIq+YWGgbRYqd5VJHoGOR4NUlds1KRLWgq1ZC4RqKC47F4s0l0ucYOOVhRgAjOrsiBtAzHza0AZlI8MoNv6B/fS1ZgTgKEs0qIQVCKiStVlKjFGy5aYwqTqM03w+JMASrDv7jJsYw8AxYPq31mXsuROohuB2++7vdis2bN4+bbrqJGTNmUFLSY4PQ9DqklFz51oPUxgNE2MdVOyCkxKPrnBpKYNX8NGAhpTqoa6jjxLGnQOYAA5Sb7JPq5VR//ibbRR7eNFgYWUt6XGHx+oE0lC8klgxxLDEGqTZ2WFRsmoZFsTA6kk4imaLJlmSoLQt/1nAqghWc0ucUo2vUqfRa57720ksZXL4Dhx7GldCxaJC0Q8oKKJJUsRUt3Y5FizI4t3Ude9M2OP0aA1WbHBBFgfOfbXHwVV9B9hA46wGwOnYVWb16NRdeeCGvvvoqo0aNMlCsSXvzxZYGdoTDBJUV7GuOlEODvgnJ90IxToqqxC1eNjuG0yxc5E+6l5yh/TpftMn++eoV1sZz8DksLK7dRp1H4G+KcGpoDbhCWPwRcmNWInaVtGSCZtVKSdTJid7RpHSdhuQSFuUW0RjYyrjCcUzsP9HoGnUqvc+5S0ns5SmoW1eSrVtJoaHuyhciSDokqk+QyLURlS5OSk9hCWyDRAj6fQeyBhpdA5MDoVph9BUtrz2orKxk4sSJPPbYY5x55pkGiDPpSF746kMSSiWKunMSd0sSIACbDgNDfoZGMyiJNVAv4qRn9mFCQQac8whYnfv9XRODEAoWVdAUTRJFUlgeoqQyTEjVQGhU50kac1VGBi0sdqSTdCT42plDRWwtt2Rmctqpf+TCPt067fIR0auc+8KqhXwy+3ecPnMlmUk7VqmhKpCwgNAh5oZQiZ1wppeg189Rx52O47RroWEzZA4Ch5k8orsSCAQ455xzuPHGG7nyyiuNlmPSAdTzGSmZZKdDb0GgAj5d5/RkBoPTB2EVGp8HVCZMmgLFZq63LsvoKxhT+xD/rvSRE7Ugq5MklBbHriDpu10SSItT5QThzcCT0HGk5TNgwHd4LrKZqYVj6M0ZAdoSfrZH8P6m9/nbrMexri2nMegibBcIKVE0sOsgVQvN/bJ5v/j7POG9Avdlf2XgpDvAmQ6FY03H3kEIIe4QQiwVQhzVUcdIJBJccMEFjBs3jrvuuqujDmNiIDO2zmBLZDFYGkBXd/l3AShSZ0DYRproQzAW5+uwn6oBl9DPdOwdhhCivxDib0KItw73N/Sc4dT0PZdj3PWojRp6JB1NKkgkQoI/JbHFodquYbeqOFQYOXQwbrsVq7BSE6lpzyp1O3rFlfuKHat4cP6TeEIpBjdLpD1F3C2w6RJrQkEiCKa52dKnhM3OUnKLBnJGaZ7RsnsFUsonhBAdNnLSdZ0f/ehHKFn9cZxyPVOmLSPLY+fWMwbtSglp0n2JpqI8vvhxZpfPJqWn0MJRrG43kihWXeKSGt/RLNwz4Tf8MzCc5Q1Rxg3K4oZhuUZL79FIKTcD1x2uc69sjHDv2ytQ1QF4c35JZN0cjgkswqfVUu8TaALqbCpbfYJq1U6z1kDS6yRPi+NSXGhSw2PtsXkf2kSPd+4b6rdw0wdTCcZtaLYmKlwabruKlqljVSSeZIpGh4fyIbmk2VOcfdJobjh1MOIgIUtNOo8jSQ957733snZ7MyMu+xmF6S6EEESTGne9tZxnLh9jtnM359cLfk1FsIIsWxY19TWoHhW33YrARjKVoJ+nDw99/59gsfNjo8Wa7OJAfbo+FGfSk/PQ9RSj9FUc1biEozZtIqrpyIYEzohgeT/BkkEKDW4rocgQLM56SrNyWFC1gDE5Y8h35+Oxmc69RzNlxsM0pNah2gPEhc6KQiioFdR4FLItGnpcsLqgL19njeHBO27G4eu2CbC6FUKIq4AxwKMHK3u46SGffvpp3nnnHSb/5mWy07y7HLnTqmK3KCwrb2J0cfrhVcDEcNbWrwUBNsXG0nVLsWfaSXOlkdSTZDoycVldnNb3LLD0ykBEXZr99elYUuPKv35OJJnieH0FY9U1jK5YjyMURxWSuAaqBuVZbmr8KnFUdIfEIUuIpOoJJUIUeAq49v/bu/vgqMp7gePf3242IYQQIiFCEiAEFEgCVUoEpWL1euWlvIhacS4vUgXH/qHYTpWpxRlvaeWl1ul04KLXtreCvX2RoRWlYUrr9YrIdbggyItCfEFBIkoSQoAkm9393T920QghN5tk9+ye/X1mAju7z3Oe357fefY5Z8/Z55R/x7H3lihcfc5996f7+ejMHry+OvCEQKA2S9hyFdTkQ9NlPuqH9WZb/kTez7uJjGybzCReVHWdqj5E+LaQ04DHRKTb5pbYuHEjy5cvp7KyEk9aBmkXzC3v83qoO9fSXc2ZOFK/n/rNm2l47nl8DU1Uba/Cl+2jrKAMgJZgC3XNdUwomMD80vkOR5uaRKSviDwNXC0iP+xovb+89QnnWgJ4QgFGeT4gO3SO9NMteFuCZJ0Kkd0ARZ/DjW81c07TCIoPPGfJ75XDTYNuojSvlHml8/CIq4e2DnHtkfu+Y/XcvenHkNmItvrmVUU4k+Xl40HC1+rS2ekfztGsUh6bMsK+onVAZNarbp35avv27dx///1s2bKFIUOGMO7UMXa8X8OAnC9/7lR71s+1JfYtTbIJ1NVx/NEf4cvLY0BmDwr+4yWG9PVSPHkqB2sP4vP6KO5ZzGPjH6M0r9TpcFOWqtYAUd3H/e1jp3jujSN8Wt9ISEP0oYHLP6/BWx+Ec4IH8HoU9Qj59QFGHPPzziDweJoZXdiH+uZ6CrIKSPemx+ZNJRlXDu4hDfHAy0sg/T3a/C1ECOoC/XiH4Zwpu4s/zLyWnJ52cZUbvPvuu9x+++2sX7+eMWPGAHDrVYUcPnGGQycaSPcKjf4g868tJtPu/pV0ap79FRklQ/D0yGTHtm3UVAe4Jbs3B/2NFOcUk98zn0cqHrEjtyRz+MRpfvmP9xg9MIcjn9VzEzu57EgV/Q6dJdSkCJFfPoQAr9KSLlz5qfJBSQYeby1ng3XkZw3hvq/d5/A7SRyuG9xDoRAT13+ds+mBVicdzs9SE37cOygMHPwIt064jivOzz5nkl51dTVTpkxhxYoVTJo06YvnRYQlk0dwuqmFkw3NDO6bhbedecdN4grW1uIrLGTv3r3s2r2Le++5F19dLRW+b9CvYgJ5mXlOh2g64fkdHzOify+Ofl7HzYH/4dadm8n7zE9a8MupiITwP+oL/3z5dI6HnPQcxg8czxPXP9F+AynIdYP7sh3LaAg2E7rE/btzAgFWVTxMxdWT2nzdJKfTp08zdepUFi5cyIIFC9os07uHz37+luQ8vXrxweFDbN26lQUL7iY7O5um6mquGDWRtEy7ODJZNQeCeIMBrqpayxX795BT3RI5WtcvpiQKRQ7fVZSaAg8NRYWMyBvGpCH2Wd4W1313taFqA6G2zp0reBRuvmIWFVffG//ATMz4/X7uuOMOxo0bx6OPPup0OCaGqsddw4HNm5k9fRp9L+uL/6Mj9Bh+JWm5NrAns7FFWQx6cxkT979O7od+vCieyMD+xae5V/H1DlBzbYD9YzP5PLcPYweM5IaiGxyMPELIDYIAAAlYSURBVHG56sj9UO0hPBpqc3DPCLaw5PplfHv4tx2IzMSKqrJo0SJ69OjB6tWr7aJIFzty5AjTFixg9cqVDGpsJHDiM/rMnk1WRYXToZkumvjeemp3vI3/qI/WN/M7/5V8wKP0GXKWvJImcvsPoHjQ9Sy9bZVdPNcOVw3uH9Z/SLpC04Wf76pkiNjA7kJLly7l0KFDvPLKK6SluWpzNq3U1NQwefJklixZwqy5ds91NwkFAxz47WYGHPegevHOeUiUU6MDlI/y4M0pIevrC+g3/rtgO/LtctWn4dj+Y+nlTScY9NPS6oKp3GCQZ275jYORmVhYu3YtL7zwAtu3b6dnz55OhxM3f3/nU/68+xNaQkqPNC+Lri9hVFGO02HFTGNjIzNmzGD69Ok8+OCDTodjutnZI+8xIDQQDV08WAc9SktpC/UjSqkpHU3+zJ+Eb+1s/l+uWkueRg8nX2mmb0DpHQzRMxSiyO9n8eU3MnLQdU6HZ7rRiy++yLJly6isrKRfv35OhxM373/WwJ92HmNYfjYj+/emuG9Pntp6iDPNAadDi4lgMMicOXMYPHgwK1eudDocEwPelhAELp54UoGThcq7xaW8VLyUXcMesIE9Cq45cm9qamLmzJlMGjeVh+f8gDcP/J4Sfwtlo+YiOQVOh2e60Y4dO1i4cCGVlZUMHTrU6XDi6g87jzIs/8s5s0WEotye/PXt49xZEd28+4lOVVm8eDGnTp0KzzRoH+yuJKHwwC5AKPKcepSjJcreCflUNd9ISD2Mt0mnouKawX3u3LkUFRWxatUqPB4PM675ntMhmRhoampi1qxZrFu3jrFjxzodTtypXmJepg7PuJ88Vq1axWuvvca2bdvIyLD54d1KkcgMNeGvkgNe5diIEPuuymdf42Rqc8cwv+xycrPs4rloJOzgLiJZwL8BfuBVVf3dpcoePXqU7OxstmzZYnv3LldVVcXatWuZMmWK06E4YnZFEU/+7TAj+p+/S67ySd05Hpk83NG4OiKaPl1bW8uaNWt44403yMlx7/UEbhVNrlvES5P4yPAE8PhCBAf7qBk6mKE55RSOncu00YU2g2gniGpi7vKLyDzglKq+JCJ/VNXZlyqbmZmp1dXV9OnTJ44RJrWkvcy0oKBAjx8/7nQYjtqyr5pNe4/TEgyRkebhnm8MYczgNm96lFB5jqZP+3w+3bNnD2VlZXGMMKklba6H9srRJ8vH4Qv62VVwJSdLL2f8NdcwZ9b0OEacNDqc54Q9cgeKgH2Rx8ELX2x9P2CgOTc3d3+8AmtHHnAyCWLYoqqT4xFMd6uurj4pIh9FWS3eeYlne3lrLt1WouU5qj5dXl5ufbrjMSR1rm97c2sk1/8Nm2A10IEfPHYmN53NZzzrtVenw3lO5MH9GOENZA9tXNXf+n7AIvK/qur4CdhEiCMRYoglVY360vh4r5N4tpdk+bY+naQxdELMc92ZevFsq7P1uivfiTy4bwRWi8i3gJecDsYY02XWp1OH5dphCTu4q+pZ4DtOx2GM6R7Wp1OH5dp5brm0/N+dDiAiEeJIhBgSTbzXSTzbc2u+E+V9JUIciRBDLHX2/XWmXjzb6my9bsl3wl4tb4wxxpjOccuRuzHGGGMiknpwF5EsEXlORJ4VkTnxWr6IlIvI7yJ/5ZHnVojIL0VkRTtlnhGRfW23FvsYRGS4iPxGRP5TRJZ0fQ0lPhH5vojsPr+OYtxWTLfHNtorEZFfi8iGWLcVLw6sw2j6V9T9t7vjcFMf7myuO7vdi8itkbZeFJFbOlhnpIg8LSIbROS7UbaXJSK7RGRaFHW+KSLbIm1+M5r2LqKqSfsHzAOmRx7/MV7LB54F+gA5wDPAIODnkdd+Bgy8sEyruhucjiFSZqPT+YvjdvI4UJ7s22M77Ua1TSXyX7zXYUf7V6zXdbRxRF5L6j7c1Vx3NhdALvDrKOt4OlHnx8ASYFoUdW4AKoHfAsO6sn6T+sid8O8oj0YeXzRRQgyXn6Oqp1S1HsgGCluV+zhS78IyCRODiNwF/K0LMZm2xXp7TAXxXocd7V8JFYdL+rBT/WUpsKajhUVkBvA68I8o6twMHARORBnbNlWdQnin4F+jrPsVyT64n58oAWLzXi61/HoRyRGR3kAD8EmrcgMj9S4skxAxRD4UBqvq012IKeGJyHwR+YWIFMax2Vhvj6kg3uuwo/0rYeJwUR+Oa64lbCVQqaq7O1pPVTep6nVANKeJbgTGA/8CLBKRDr0/VT1/Y7w6oEt3S0rqq+UlfHOC1UAT8Lq2c3OC7lg+MFlV50XOwT1MeJ7fVaq6X0SWE05Gs6r+8BJlfgrcBWwFFqtqczxjAHyEJ5R4GTinqt/vlhWVwETkbuAB4H3gJ6oak3OmkbZiuj220V5f4KfAPwO/UtXlsWwvHhxYh9H0r6j7b3fHgYv6cGdz3dntXkQeBO4GdgJ7OrJzFDnvfRvhz9W3VbXDR/yR+guAk6r6cgfL3wZMInwqZq2qvhpNe19ZVjIP7sYYY4y5mH11aIwxxriMDe7GGGOMy9jgbowxxriMDe7GGGOMy9jgfgERSY/M8rZaRNaIyD1dXJ608dyGyP8FIvK9rizfdI7lOXVYrlOD5fmrEvaWrw5aRPh3kJUAIpImIj8C+hGeROIh4PfAq0A58BTwEeEJBwRoUNWlInIAWA/8VUTuJPzThr3ADqBURB4HNgCFIlIMrAM2AWXAPcCdwD8BZ4CQqv4g1m88xVieU4flOjVYnluxI/eLlQE7RcQjIr8ANgMTgVOAHxgZKfcU8HPgW4QnKsgkPPFAiYikA8dVdQVQRXgnqha4U1X3AwdV9XHCyT/vgKo+GSnXP1L2PsIbo+l+lufUYblODZbnVuzI/WIHgIrI3t9DIvJn4MNIQgEQkSZVDYlIC+HJDTzAZlXd1KpMfeThVMIbxDoR+a/Ic21NLnA28v/5ZdJOWdN1lufUYblODZbnVmxwv9izwJMSvpNPAHgTyBORpwjv4T3RRp3ngdUicj2QrqqLW732FrBCRAYA3shzJyR857bn24njTyLyDNBMeI/QdC/Lc+qwXKcGy3MrNkNdghKRMcDtQAGwXFUPOxySiQHLc+qwXKeGRMmzDe7GGGOMy9gFdcYYY4zL2OBujDHGuIwN7sYYY4zL2OBujDHGuIwN7sYYY4zL2OBujDHGuMz/AbXwWjnCSVWwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 510.236x141.732 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "figure_si_6(results);\n",
    "plt.savefig('results/figures/si_figure_6_parameter_recovery.png', dpi=330, bbox_inches='tight')\n",
    "plt.savefig('results/figures/si_figure_6_parameter_recovery.pdf', dpi=330, bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bias Analyses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>subject</th>\n",
       "      <th>dataset</th>\n",
       "      <th>parameter</th>\n",
       "      <th>type</th>\n",
       "      <th>value</th>\n",
       "      <th>is_recovered</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>krajbich2010</td>\n",
       "      <td>gamma</td>\n",
       "      <td>generating</td>\n",
       "      <td>0.790000</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0</td>\n",
       "      <td>krajbich2010</td>\n",
       "      <td>gamma</td>\n",
       "      <td>recovered</td>\n",
       "      <td>0.980000</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0</td>\n",
       "      <td>krajbich2010</td>\n",
       "      <td>s</td>\n",
       "      <td>generating</td>\n",
       "      <td>0.010750</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0</td>\n",
       "      <td>krajbich2010</td>\n",
       "      <td>s</td>\n",
       "      <td>recovered</td>\n",
       "      <td>0.010854</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0</td>\n",
       "      <td>krajbich2010</td>\n",
       "      <td>tau</td>\n",
       "      <td>generating</td>\n",
       "      <td>0.490000</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   subject       dataset parameter        type     value  is_recovered\n",
       "0        0  krajbich2010     gamma  generating  0.790000           0.0\n",
       "1        0  krajbich2010     gamma   recovered  0.980000           1.0\n",
       "2        0  krajbich2010         s  generating  0.010750           0.0\n",
       "3        0  krajbich2010         s   recovered  0.010854           1.0\n",
       "4        0  krajbich2010       tau  generating  0.490000           0.0"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = results.melt(id_vars=['subject', 'dataset', 'parameter'], value_vars=['generating', 'recovered'], var_name='type').sort_values(['subject', 'parameter']).reset_index(drop=True)\n",
    "df['is_recovered'] = (df['type'] == 'recovered').astype(float)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "v\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/felixmolter/anaconda3/lib/python3.6/site-packages/pandas/core/generic.py:4384: FutureWarning: Attribute 'is_copy' is deprecated and will be removed in a future version.\n",
      "  object.__getattribute__(self, name)\n",
      "/Users/felixmolter/anaconda3/lib/python3.6/site-packages/pandas/core/generic.py:4385: FutureWarning: Attribute 'is_copy' is deprecated and will be removed in a future version.\n",
      "  return object.__setattr__(self, name, value)\n",
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using advi...\n",
      "Average Loss = 8,825.4: 100%|██████████| 50000/50000 [00:40<00:00, 1227.10it/s]   \n",
      "Finished [100%]: Average Loss = 8,768.7\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [value_sd, is_recovered|dataset_offset, is_recovered|dataset_sd, 1|dataset_offset, 1|dataset_sd, is_recovered, Intercept]\n",
      "Sampling 2 chains: 100%|██████████| 61000/61000 [42:04<00:00, 12.55draws/s]\n",
      "The number of effective samples is smaller than 25% for some parameters.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                                            mean        sd  hpd0.95_lower  \\\n",
      "is_recovered|dataset_sd             6.951740e-06  0.000009   4.149151e-11   \n",
      "is_recovered|dataset[folke2016]    -4.179241e-07  0.000007  -1.428727e-05   \n",
      "is_recovered|dataset[krajbich2010] -2.107000e-07  0.000007  -1.333805e-05   \n",
      "is_recovered|dataset[krajbich2011]  1.790401e-07  0.000007  -1.327700e-05   \n",
      "is_recovered|dataset[tavares2017]   4.142537e-07  0.000007  -1.297363e-05   \n",
      "1|dataset[folke2016]               -3.744099e-05  0.000015  -6.945447e-05   \n",
      "1|dataset[krajbich2010]             1.270730e-05  0.000015  -1.855124e-05   \n",
      "1|dataset[krajbich2011]             2.604384e-06  0.000015  -2.838402e-05   \n",
      "1|dataset[tavares2017]              2.132381e-05  0.000015  -1.064517e-05   \n",
      "is_recovered                       -4.258305e-09  0.000007  -1.208974e-05   \n",
      "Intercept                           6.499155e-05  0.000015   3.434362e-05   \n",
      "value_sd                            3.022045e-05  0.000001   2.750301e-05   \n",
      "1|dataset_sd                        2.958834e-05  0.000011   1.228271e-05   \n",
      "\n",
      "                                    hpd0.95_upper  effective_n  gelman_rubin  \n",
      "is_recovered|dataset_sd                  0.000021        14559      0.999996  \n",
      "is_recovered|dataset[folke2016]          0.000013        35239      1.000005  \n",
      "is_recovered|dataset[krajbich2010]       0.000013        33201      0.999998  \n",
      "is_recovered|dataset[krajbich2011]       0.000013        33747      0.999986  \n",
      "is_recovered|dataset[tavares2017]        0.000014        35476      0.999993  \n",
      "1|dataset[folke2016]                    -0.000008        40785      1.000004  \n",
      "1|dataset[krajbich2010]                  0.000043        40210      1.000018  \n",
      "1|dataset[krajbich2011]                  0.000033        40828      1.000024  \n",
      "1|dataset[tavares2017]                   0.000051        40301      1.000032  \n",
      "is_recovered                             0.000012        34719      0.999990  \n",
      "Intercept                                0.000095        40214      1.000021  \n",
      "value_sd                                 0.000033        47475      0.999984  \n",
      "1|dataset_sd                             0.000052        41971      0.999984  \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xtc1uXh//HXxUGWlQ0zHIilSypREUNN3VKTblErzK+Vtn5fsJk6arN1mLOcU8qlK93KaQczN0xLS9Nceerk12UppaKVZmqZQSQpqKWkotfvDz7cA29AUJELez8fj/vBfV+f6/S5OLz5HLgx1lpERERcE1TbExARESmPAkpERJykgBIREScpoERExEkKKBERcZICSkREnKSAEhERJymgpE4zxrxX23M4k4wxg4wxUyrZHmaMmWuM2WaMWWOMaXbmZidyeimgpE6z1nY5lfbGmODTNRdH5jAYKLDWtgD+Dvz1NPYtckYpoKROM8Z8732MNMasNMZkGWM+NsZcXVkbY8xDxpg1QGdjTIIx5v+MMWuNMcuMMZFevRbGmDeNMRuMMeuMMZeaYo95Y3xkjBng1Z1rjOlTaox/GWP6G2OCvfofGGM2GmOGedu7G2PeMca8AHzklf0/Y0ymtw/PlASXMeZ2Y8xnxpj/A35xgiXpC2R4z+cBicYYcxJLK1LrFFBytvgVsMxaGw+0BbIqqXsu8LG19ipgDfAP4CZrbQIwA/iLV282MNVa2xboAuQC/wOUjHEt8JgXaHOAkrCqByQCiyk+otlnre0AdACGGGOae/13BEZZa2ONMS299r/w9uEocJvXdzrFweQDYk+wDk2ArwCstUXAPuDCE7QRcVJIbU9A5DT5AJhhjAkFFlprKwuoo8B87/nlQGvgDe9AIxjINcacDzSx1i4AsNb+AGCM+SXworX2KLDLO6rpACwBJhtjwoBewEprbaExpicQZ4y5yRvvAiAGOAxkWmu/8MoTgQTgA28e5wB5wFXACmvtt974c4HLKtm38o6W9IabUicpoOSsYK1daYzpClwHPG+MecxaO7OC6j94AQPFP9A/sdZ2Ll3BGNOggrblni6z1v5gjFkBJFF8JPRiqfq/s9YuO67/7sCB4/rNsNY+cFy9G6lewGQDTYFsY0wIxYGYX432Is7QKT45KxhjLgHyrLXPAs8BV1ax6RbgImNMZ6+fUGNMK2vtfop/yN/olYcZY+oDK4EB3rWli4CuQKbX1xzgduBqoCSQlgFp3pEdxpjLjDHnljOPt4CbjDERXr2G3j6tAbobYy70+rj5BPuzCEj1nt8EvG31LwukjtIRlJwtugN/MMYcAb4HUqrSyFp72Dv9NtkYcwHF3xOPA58A/ws8Y4x5CDhCcTgsADoDGyg+shlhrf3G6245MBNYZK097JVNB5oB67ybFb4FbixnHpuMMX8Clhtjgrzx7rLWrjbGjAXep/ga2DqKT0NW5DmKjyC3UXzkNLAq6yDiIqNfrkRExEU6xSciIk7SKT45a3l/5xR2XPH/Wms/qo35nE7GmFEEXo962Vr7l/Lqi9RFOsUnIiJOqs4RlJJMREROlxO+w4muQYmIiJMUUCIi4iQFlIiIOEkBJSIiTlJAiYiIkxRQIiLiJAWUiIg4SQElIiJOUkCJiIiTFFAiIuIkBZSIiDhJASU1puhY9d6+sbr1ReTsVp13M9dPD6m2Cet3V7nuyHaNanAmIuIYvVmsiIjUTQooERFxkgJKREScpIASEREnKaBERMRJCigREXGSAkpERJykgBIREScpoERExEkKKBERcZICSkREnKSAEhERJymgRETESQooERFxkgJKREScpIASEREnKaBERMRJCigREXGSAkpERJykgBIREScpoERExEkKKBERcZICSkREnKSAEhERJymgRETESQooERFxkgJKAjRr1ozdu3ef1j5zNm3g8Vu68lhyBxY9+gDW2oA61lqGDx9OixYtiIuLY926df5tGRkZxMTEEBMTQ0ZGhr987dq1tGnThhYtWjB8+HB/v/n5+fh8PmJiYvD5fBQUFJxwDBFxiwJKzoiF4/9Av1GTuP/VTPbs/JzP3nsroM6SJUvYunUrW7duZdq0aaSlpQHFYZOens6aNWvIzMwkPT3dHzhpaWlMmzbN327p0qUATJgwgcTERLZu3UpiYiITJkyodAwRcY8C6iyxY8cOrrjiCu644w5at27NbbfdxptvvskvfvELYmJiyMzMJD8/nxtvvJG4uDg6derExo0bAdizZw89e/akXbt2DBs2rMzRzaxZs+jYsSPx8fEMGzaMo0ePVntu+7/9hkMHvuOSth0wxtDu+gFsemdJQL1XX32VlJQUjDF06tSJvXv3kpuby7Jly/D5fDRs2JDw8HB8Ph9Lly4lNzeX/fv307lzZ4wxpKSksHDhQn9fqampAKSmppYpL28MEXGPAuossm3bNu6++242btzIp59+ygsvvMC7777LxIkTeeSRRxgzZgzt2rVj48aNPPLII6SkpACQnp7OL3/5S9avX09ycjI7d+4EYPPmzcydO5dVq1aRlZVFcHAws2fPBmDAgAHEx8cHPGbOnBkwr/3ffkODiCj/6wsiItmXFxgKOTk5NG3a1P86OjqanJycSsujo6MDygF27dpFZGQkAJGRkeTl5VU6hoi4J6S2JyCnT/PmzWnTpg0ArVq1IjExEWMMbdq0YceOHXz55ZfMnz8fgB49erBnzx727dvHypUreeWVVwC47rrrCA8PB+Ctt95i7dq1dOjQAYDCwkIiIiIAmDt3bpXnVd71JmNMletVt/x0zEVEap8C6iwSFhbmfx4UFOR/HRQURFFRESEhgZ/ukh/OFQVGamoq48ePD9g2YMAAtmzZElB+7733ctttt5GQkADAhR2v5aqbBrE/72t/nX15uTS46GcBbaOjo/nqq6/8r7Ozs4mKiiI6OpoVK1aUKe/evTvR0dFkZ2cH1Ado3Lgxubm5REZGkpub6w/WisYQEffoFN+PSNeuXf2n6FasWEGjRo1o0KBBmfIlS5b4b0BITExk3rx5/tNj+fn5fPnll0DxEVRWVlbAIyUlheDgYP9rX9pIGlz0M+rVP4+dGz/EWsv61+bSsnuvgPklJyczc+ZMrLWsXr2aCy64gMjISJKSkli+fDkFBQUUFBSwfPlykpKSiIyM5Pzzz2f16tVYa5k5cyZ9+/b191Vyt19GRkaZ8vLGEBH36AjqR2Ts2LHcfvvtxMXFUb9+ff8P8DFjxnDrrbdy5ZVX0q1bNy6++GIAYmNjGTduHD179uTYsWOEhoYydepULrnkkmqPfeODjzFvzO84cugHLuvSg8t/cS0Aa+b9C4CrbhpEnz59WLx4MS1atKB+/fr885//BKBhw4aMHj3af6rxz3/+Mw0bNgTgqaeeYtCgQRQWFtK7d2969+4NwMiRI7nlllt47rnnuPjii3n55ZcBKhxDRNxjyjsnX4EqVxQpMWF91f+eamS7RjU4ExFxzAkv/uoUn4iIOEkBJSIiTlJAiYiIkxRQIiLiJAWUiIg4SQElIiJOUkCJiIiTFFAiIuIkvZOE1Kg5a4eVeT0w4ZlamomI1DU6ghIREScpoERExEkKKBERcZICSkREnKSAEhERJymgRETESQooERFxkgJKREScpIASEREnKaBERMRJCigREXGSAkpERJykgBIREScpoERExEkKKBERcZICSkREnKSAEhERJymgRETESQooERFxkgJKREScpIASEREnKaBERMRJCigREXGSAkpERJykgBIREScpoERExEkKKBERcZICSkREnKSAEhERJymgRETESQqoOqhZs2a0adOG+Ph42rdv7y/Pz8/H5/MRExODz+ejoKCg3PYZGRnExMQQExNDRkYGAIcOHaJXr160bt2aJ5980l936NChrF+//qTnmv/el2wetYzNo5aR/96XJ9XHv/71Ly666CLi4+OJj49n+vTple5LwBwqWJf58+fTqlUrrr76avbs2QPA9u3bGThw4EnNU0ROLwVUHfXOO++QlZXFhx9+6C+bMGECiYmJbN26lcTERCZMmBDQLj8/n/T0dNasWUNmZibp6ekUFBSwbNkyEhIS2LhxI9OmTQNgw4YNHDt2jHbt2lVrbvn5+f6Pu/69mZgHryHmwWvY9e/NFO7fW2G7Y8eOsW/fvnK3DRgwgKysLLKysrjjjjsq3ZfjVbQukyZNYvXq1aSkpPDCCy8A8Kc//YmHH364WvsrIjVDAXUWefXVV0lNTQUgNTWVhQsXBtRZtmwZPp+Phg0bEh4ejs/nY+nSpYSGhlJYWEhRUZG/7ujRo3nooYeqNHZRURGLFi0iOTmZfv36+cc6LzaCkHPrEXJuPc6LjWDLe28HtN2bm82bTz/K5Zdfzrvvvlvl/a1oX45X0boEBQVx6NAhDh48SGhoKP/5z3+IjIwkJiamynMQkZoTUtsTkOozxtCzZ0+MMQwbNoyhQ4cCsGvXLiIjIwGIjIwkLy8voG1OTg5Nmzb1v46OjiYnJ4ebb76Z559/nquuuooRI0awaNEiEhISiIqKqnQu27Zt47nnnmPevHl06dKF++67j27duvnHqhde31+3Xvg57M/LBaDoyGE2r1jCBwtm8X3Bbq68fgDvv/8+jRo1Knec+fPns3LlSi677DL+/ve/07Rp0wr35XgVrcuYMWNISkoiKiqKWbNmccsttzBnzpxK91dEzhwFVB20atUqoqKiyMvLw+fzccUVV9C1a9cqtbXWBpQZYwgJCfGf5jpy5AhJSUksWrSIe++9l507d5KSkkJycnKZdvPnz2fAgAGMGjWKdevWcf7551dpLICpt/k4drSI/mMnc3GbBIAKw+mGG27g1ltvJSwsjKeffprU1FTefvvtSvuvCp/Ph8/nA4qvZfXp04ctW7YwceJEwsPDeeKJJ6hfv/4JehGRmqJTfHVQyVFNREQE/fr1IzMzE4DGjRuTm1t8hJKbm0tERERA2+joaL766iv/6+zs7ICjpCeffJLU1FTef/996tWrx9y5cxk3blxAXz6fjyeeeILXX3+d/v378+KLL/LDDz+UGetwwUH/68MFhZx/0c8A6Df6b1wc156X/pTGksfTyfv8swr398ILLyQsLAyAIUOGsHbt2irvS1XW5eDBg2RkZHDnnXfywAMPMGPGDBISEpg9e3aFcxKRmqeAqmMOHDjAd99953++fPlyWrduDUBycrL/TraMjAz69u0b0D4pKYnly5dTUFBAQUEBy5cvJykpyb+9oKCA1157jZSUFA4ePEhQUBDGmDLBU6JBgwbcddddfPjhh/z1r3/l3XffpWXLlowYMcI/1vef5FF04DBFBw7z/Sd5XNb5GgAubpNA/zFP8LsX36ZRs0uZ/9Dv6dSpE+vWrQsYpyRcABYtWkTLli2rtC8lTrQujz76KHfffbf/OpwxhqCgIA4ePBjQl4icQdbaqj7EAdu3b7dxcXE2Li7OxsbG2nHjxvm37d692/bo0cO2aNHC9ujRw+7Zs8daa+0HH3xgBw8e7K/33HPP2UsvvdReeumldsaMGWX6//3vf29XrFhhrbW2sLDQ+nw+GxsbaydPnlyl+RUWFtrXX3/d/zo69Upb76Jzbb2LzrVNB11px6/7tsLHpk2b7KZNmwL6HDlypI2NjbVxcXG2e/fudvPmzSfcl8GDB9sPPvig0nWx1tqcnBx73XXX+V+/9NJLNjY21nbp0sXm5eVVaZ9F5KScMHeMLec8fkVZVpNBKWen+On9y7wemPBMhXVHtiv/GpSInJVOeMFYp/hERMRJCigREXGSAkpERJykgBIREScpoERExEkKKBERcZICSkREnKSAEhERJymgRETESQooERFxkgJKREScpIASEREnKaBERMRJCigREXGSAkpERJykgBIREScpoERExEkKKBERcZICSkREnKSAEhERJymgRETESQoocUbRMXtG2ohI3RBS2xMQKRESZJiwfne12oxs16iGZiMitU1HUCIi4iQFlIiIOEkBJSIiTlJAiYiIkxRQIiLiJAWUiIg4SQElIiJOUkCJiIiTFFAiIuIkBZSIiDhJASUiIk5SQImIiJMUUCIi4iQFlIiIOEkBJSIiTlJAiYiIkxRQIiLiJAWUiIg4SQElIiJOUkCJiIiTFFAiIuIkBZSIiDhJASUiIk5SQImIiJMUUCIi4iQFlIiIOEkBJSIiTlJAiYiIkxRQIiLiJAWUiIg4SQElIiJOUkCJiIiTFFAiIuIkBZSIiDhJASUiIk5SQImIiJMUUCIi4iQFlIiIOEkBJSIiTlJAiYiIkxRQIiLiJAWUiIg4SQElIiJOUkCJiIiTFFAiIuIkBZSIiDhJASUiIk5SQImIiJMUUCIi4iQFlIiIOEkBJSIiTlJAiYiIkxRQIiLiJAWUiIg4SQElIiJOUkCJiIiTFFAiIuIkBZSIiDhJASUiIk5SQImIiJMUUCIi4iQFlIiIOEkBJWU88cQTtG7dmlatWvH444/7y8eOHUuTJk2Ij48nPj6exYsXA7Bq1Sri4uLo0KED27ZtA2Dv3r0kJSVhrS13jGlD+pK9Kcv/uuDrnTx+89UAfP7hKsZ2/TmTb72Gv/1PF9585rGA8kn9OvHM4BvYvHJ5jayBiLghpLYnIO74+OOPefbZZ8nMzKRevXr06tWL6667jpiYGADuuece7r///jJtJk2axPz589mxYwdPPfUUkyZN4uGHH+bBBx/EGHNS82gW34lBk1/gcOEBJg+8hpZX9yxTDvD1lo94/t5U3mrdmMTExFPYaxFxlY6gxG/z5s106tSJ+vXrExISQrdu3ViwYEGlbUJDQyksLOTgwYOEhoayfft2cnJy6Nat2ynPp94559KkZVv2ZH8RsC3q8jYkDrmPKVOmnPI4IuImHUGJX+vWrRk1ahR79uzhnHPOYfHixbRv396/fcqUKcycOZP27dszadIkwsPDeeCBBxg6dCjnnHMOzz//PPfffz8PP/zwCceaO+o3hIb9BICjR45gggJ/VzqwN5+dH62lx5D7OFCwJ2B7VMs43njp6VPYYxFxmQJK/Fq2bMkf//hHfD4f5513Hm3btiUkpPhLJC0tjdGjR2OMYfTo0dx3333MmDGD+Ph4Vq9eDcDKlSuJiorCWsuAAQMIDQ3lyJU/ENrgJwFjDfjL00THxgPF16Ay7r7Nv21H1mom33oNxgTR/fbhNL70Cj7/cFXghCu4xiUiZwcFlJQxePBgBg8eDMCDDz5IdHQ0AI0bN6bomCU4yDBkyBCuv/76Mu2stYwbN465c+fy29/+lvT0dHbs2MEb//g9kf1aVWsOpa81VebrTz+iZcuW1epbROoOBZSUkZeXR0REBDt37uSVV17h/fffByA3N5fIyEgmrN/Nu7NmEdwkhgnrd/vbrV30ImFtu/HMjqNk5e7l2U/3kVQ/iGOHj9bIPHM/+4S3p/+Nl2fOqJH+RaT2KaCkjP79+7Nnzx5CQ0OZOnUq4eHhAIwYMYKsrCx2/3CM8Kim3Dhqor/N4cKDrHttLr+e+jIAv7ztN8z+w+2saHAOjfo1P21zKzn1d+SHQs4Lb8QNf3hEd/CJnMVMRX+rUg6d8JcyR00nMrJdI+Kn9y9TNjDhmUrrV6f/kjYiUied8O9QdJu5iIg4SQElIiJOUkCJiIiTFFAiIuIkBZSIiDhJASUiIk5SQImIiJMUUCIi4iQFlIiIOEkBJSIiTlJAiYiIkxRQIiLiJAWUiIg4SQElIiJOUkCJiIiTFFAiIuIkBZSIiDhJASUiIk5SQImIiJMUUCIi4iQFlIiIOEkBJSIiTlJA1RGffvopnTt3JiwsjIkTJ55UH+PHj6dFixZcfvnlLFu2zF/erFkz2rRpQ3x8PO3btz9dU65R88YOZ1xiS1q3bl3ttmvXrqVNmza0aNGC4cOHY60FYOzYsTRp0oT4+Hji4+NZvHjx6Z62iFRDSG1PQKqmYcOGTJ48mYULF55U+02bNjFnzhw++eQTvv76a6699lo+++wzgoODAXjnnXdo1KjR6ZxyjUq4YSCdBwzmnUfurnbbtLQ0pk2bRqdOnejTpw9Lly6ld+/eANxzzz3cf//9p3u6InISdARVR0RERNChQwdCQ0MDts2aNYuOHTsSHx/PsGHDOHr0aECdV199lYEDBxIWFkbz5s1p0aIFmZmZZ2LqNaJ5QhfqXxBepmz79u306tWLhIQErr76aj799NOAdrm5uezfv5/OnTtjjCElJeWkQ19EapYCqo7bvHkzc+fOZdWqVWRlZREcHMzs2bMD6uXk5NC0aVP/6+joaHJycgAwxtCzZ08SEhKYNm3aGZv76TZ06FD+8Y9/sHbtWiZOnMidd94ZUCcnJ4fo6Gj/69LrADBlyhTi4uL49a9/TUFBwRmZt4iUT6f46ri33nqLtWvX0qFDBwAKCwuJiIgIqFdynaU0YwwAq1atIioqiry8PHw+H1dccQVdu3at2YmfZt9//z3vvfceN998s7/s0KFDAfUqW4e0tDRGjx6NMYbRo0dz3333MWPGjJqbtIhUSgHlsKlTp/Lss88CsHjxYqKiogLqWGtJTU1l/PjxZcoXLFhAeno6ANOnTyc6OpqvvvrKvz07O9vfX8nHiIgI+vXrR2ZmZp0LqGPHjvHTn/6UrKysMuVHjx4lISEBgOTkZNLS0sjOzvZvL70OjRs39pcPGTKE66+//gzMXEQqolN8DrvrrrvIysoiKyur3HACSExMZN68eeTl5QGQn5/Pl19+Sb9+/fxt27dvT3JyMnPmzOHQoUN88cUXbN26lY4dO3LgwAG+++47AA4cOMDy5ctP6s642tagQQOaN2/Oyy+/DBQH94YNGwgODvavw0MPPURkZCTnn38+q1evxlrLzJkz6du3L1B8farEggUL6uQ6iJxNdARVR3zzzTe0b9+e/fv3ExQUxOOPP86mTZuIjY1l3Lhx9OzZk2PHjhEaGsrUqVO55JJLyrRv1aoVt9xyC7GxsYSEhDB16lSCg4PZtWsX/fr1A6CoqIhf/epX9OrVqzZ2sVpefGAoX6xdReHefKKjo0lPT2f27NmkpaUxbtw4jhw5wsCBA2nbtm1A26eeeopBgwZRWFhI7969/XfwjRgxgqysLIwxNGvWjGeeeeZM75aIlGLKOydfgSpXlLPXhPW7q1x3ZLtGxE/vX6ZsYELFP/RHtmtUrf5L2ohInWROVEGn+ERExEkKKBERcZICSuq0omPVO/Nc3foiUnt0k4TUaSFBptrXxUSkbtARlIiIOEkBJSIiTlJAiYiIkxRQIiLiJN0kIQFK3jz1VD0AUPbvdNkw5JXK61fD+HXfVndKIlKH6AjqR0y3XIuIy3QE9SNW3Vu0RUTOJB1BiYiIk3QEJQEqurZT3TdzLbf+qGrWP82KjllCgqp+ja269UXk9FFAnSX0g7Rqqnta8/62F1arf30eRE6fKv+7jfT09KWA3iemWBTwdW1Poo7RmlWP1qv6tGbVU9vrtXvMmDGV//M5a60e1XyMHTvW1vYc6tpDa6b10pq59agL66WbJERExEkKqJOTXtsTqIO0ZtWj9ao+rVn1OL9e1fmX7yIiImeMjqBERMRJCigREXGSAkpERJz0owwoY0xDY8wbxpit3sfwCuqlenW2GmNSS5UnGGM+MsZsM8ZMNt7bf1fUryk22au/0RhzZam+jhpjsrzHopre9+owxvQyxmzx5j2ynO1hxpi53vY1xphmpbY94JVvMcYknahPY0xzr4+tXp/1TjSGaxxZr0HGmG9LfU3dUbN7fWrO8Jr91iuzxphGpcor/P50kSNr1t0Ys6/U19mfa2Rna/s+99p4AI8CI73nI4G/llOnIfC59zHcex7ubcsEOgMGWAL0rqxfoI9XzwCdgDWlxvm+ttejgjUKBrYDPwfqARuA2OPq3Ak87T0fCMz1nsd69cOA5l4/wZX1CbwEDPSePw2kVTaGaw+H1msQMKW218PRNWsHNAN2AI1KjVHh96drD4fWrDvwWk3v74/yCAroC2R4zzOAG8upkwS8Ya3Nt9YWAG8AvYwxkUADa+37tvgzNbNU+4r67QvMtMVWAz/1+nFZR2CbtfZza+1hYA7F+1Fa6f2dByQaY4xXPsdae8ha+wWwzeuv3D69Nj28PiBw7cobwzWurFddcsbWDMBau95au6OcedSl709X1uyM+LEGVGNrbS6A9zGinDpNgK9Kvc72ypp4z48vr6zfivoC+Ikx5kNjzGpjjEs/ZCqbc0Ada20RsA+4sJK2FZVfCOz1+jh+rIrGcI0r6wXQ3ztVNc8Y0/RUdqqGnck1O9V5uMKVNQPobIzZYIxZYoxpVZ2dqKqz9s1ijTFvAj8rZ1Ml76ddtotyymwl5SfTF8DF1tqvjTE/B942xnxkrd1exTnWpKrsZ3XXqLxfiE60piez3rXBlfX6N/CitfaQMeY3FP8m3aPcGde+M7lmpzoPV7iyZuuAS6y13xtj+gALgZgTtKm2s/YIylp7rbW2dTmPV4FdJYfw3se8crrIBkr/9hlN8RsrZnvPjy+nkn4r6gtrbcnHz4EVFJ/zdUGFcy6vjjEmBLgAyK+kbUXluyk+rRJyXHllY7jGifWy1u6x1h7yyp8FEk5pr2rWmVyzU52HK5xYM2vtfmvt997zxUBo6ZsoTpuavsjl4gN4jLI3MzxaTp2GwBcU3yAR7j1v6G37gOKLqSU3SfSprF/gOspehM30ysOBMO95I2Arx13wrMU1CqH4xpDm/PfCaavj6txF2YuxL3nPW1H2YuznFF+IrbBP4GXKXvS/s7IxXHs4tF6RpcbrB6yu7bVxZc1K9bmDshf8y/3+dPHh0Jr9jP++E1FHYGfJ69O6v7W94LX0Sb4QeIviQHiL/wZPe2B6qXq/pvhC4jbg9lLl7YGPKb7zZUqpT1RF/Rpgqlf/I6C9V97Fe73B+zi4ttfmuHXqA3zmzXuUV/YQkOw9/wnFPyi3UXxn489LtR3ltduCd5djRX165T/3+tjm9Rl2ojFceziyXuOBT7yvqXeAK2p7XRxas+EUHy0UUXyEMN0rL/f709WHI2v221JfZ6uBLjWxr3ovPhERcdJo9xI0AAAAMElEQVRZew1KRETqNgWUiIg4SQElIiJOUkCJiIiTFFAiIuIkBZSIiDhJASUiIk76/wC46I+iKXTkAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gamma\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/felixmolter/anaconda3/lib/python3.6/site-packages/pandas/core/generic.py:4384: FutureWarning: Attribute 'is_copy' is deprecated and will be removed in a future version.\n",
      "  object.__getattribute__(self, name)\n",
      "/Users/felixmolter/anaconda3/lib/python3.6/site-packages/pandas/core/generic.py:4385: FutureWarning: Attribute 'is_copy' is deprecated and will be removed in a future version.\n",
      "  return object.__setattr__(self, name, value)\n",
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using advi...\n",
      "Average Loss = 135.3:  37%|███▋      | 18653/50000 [00:15<00:25, 1238.69it/s] \n",
      "Convergence achieved at 18700\n",
      "Interrupted at 18,699 [37%]: Average Loss = 577.64\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [value_sd, is_recovered|dataset_offset, is_recovered|dataset_sd, 1|dataset_offset, 1|dataset_sd, is_recovered, Intercept]\n",
      "Sampling 2 chains: 100%|██████████| 61000/61000 [03:08<00:00, 323.82draws/s]\n",
      "There were 83 divergences after tuning. Increase `target_accept` or reparameterize.\n",
      "There were 153 divergences after tuning. Increase `target_accept` or reparameterize.\n",
      "The number of effective samples is smaller than 10% for some parameters.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                                        mean        sd  hpd0.95_lower  \\\n",
      "is_recovered|dataset_sd             0.087479  0.098996       0.000001   \n",
      "is_recovered|dataset[folke2016]    -0.007391  0.078848      -0.179695   \n",
      "is_recovered|dataset[krajbich2010] -0.009652  0.075945      -0.184323   \n",
      "is_recovered|dataset[krajbich2011]  0.007401  0.075575      -0.153483   \n",
      "is_recovered|dataset[tavares2017]   0.006528  0.077864      -0.153656   \n",
      "1|dataset[folke2016]               -0.099469  0.094485      -0.299020   \n",
      "1|dataset[krajbich2010]             0.085934  0.090272      -0.082528   \n",
      "1|dataset[krajbich2011]             0.034938  0.088481      -0.145946   \n",
      "1|dataset[tavares2017]             -0.025622  0.088999      -0.209345   \n",
      "is_recovered                       -0.033902  0.077722      -0.186047   \n",
      "Intercept                           0.136625  0.084162      -0.045565   \n",
      "value_sd                            0.389910  0.011266       0.368768   \n",
      "1|dataset_sd                        0.137915  0.087983       0.000003   \n",
      "\n",
      "                                    hpd0.95_upper  effective_n  gelman_rubin  \n",
      "is_recovered|dataset_sd                  0.269891         3305      1.000018  \n",
      "is_recovered|dataset[folke2016]          0.157398         7607      1.000174  \n",
      "is_recovered|dataset[krajbich2010]       0.145050         6556      1.000179  \n",
      "is_recovered|dataset[krajbich2011]       0.171411         6041      1.000210  \n",
      "is_recovered|dataset[tavares2017]        0.178103         6568      1.000037  \n",
      "1|dataset[folke2016]                     0.070924         6295      1.000455  \n",
      "1|dataset[krajbich2010]                  0.280246         6147      1.000389  \n",
      "1|dataset[krajbich2011]                  0.221465         6232      1.000630  \n",
      "1|dataset[tavares2017]                   0.156183         6387      1.000545  \n",
      "is_recovered                             0.121633         5934      1.000245  \n",
      "Intercept                                0.300313         5152      1.000713  \n",
      "value_sd                                 0.406005        36382      0.999996  \n",
      "1|dataset_sd                             0.309430         7887      1.000036  \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xlc1VX+x/HXYdHCLXNLRLLcRRF3HXNpXMdKrUwtE3OyzGxsc8wZ6+da0mjLWFraMkq525jUpJYaWZaaC5FWglOmCKbmQgkoV87vD66MBCoYwkHez8eDh9zvPef7/XyPl/u+3/P93nuNtRYRERHX+BR1ASIiIrlRQImIiJMUUCIi4iQFlIiIOEkBJSIiTlJAiYiIkxRQIiLiJAWUFGvGmM+LuobCZIy5xxjz8nnuL22MWWyM2W2M2WSMqVV41YkULAWUFGvW2j/8nv7GGN+CqsWRGu4Fjlpr6wAvAM8W4LpFCpUCSoo1Y8yv3n+rG2PWG2NijDE7jDEdztfHGDPJGLMJaGeMaWGM+cQYs9UYs9oYU93bro4xZo0x5itjzDZjTG2TaZp3G18bYwZ42y42xvQ6axtzjTG3G2N8ve2/NMbEGmOGe+/vbIz52BizAPjau+xuY8xm7z7MPhNcxpihxpg4Y8wnQPsLDEkfYJ7392VAF2OMuYihFSlyCii5XNwFrLbWhgFNgZjztC0D7LDWtgE2AS8B/ay1LYA3gae97eYDM621TYE/AEnAbcCZbXQFpnkDbRFwJqxKAV2AD8g8ojlurW0FtALuM8Zc511/a2CctbaRMaaht3977z6cBgZ51z2RzGDqBjS6wDjUAPYBWGs9wHGg0gX6iDjJr6gLECkgXwJvGmP8gXettecLqNPAO97f6wONgY+8Bxq+QJIxphxQw1q7HMBamwZgjLkBWGitPQ385D2qaQWsBGYYY0oDPYH11tpUY0x3INQY08+7vQpAXeAUsNla+4N3eRegBfClt44rgYNAGyDaWnvIu/3FQL3z7FtuR0v6wE0plhRQclmw1q43xnQEbgLeMsZMs9ZGnqN5mjdgIPMJfae1tt3ZDYwx5c/RN9fpMmttmjEmGuhB5pHQwrPa/8Vau/o36+8MnPjNeudZa//2m3Z9yV/AJAA1gQRjjB+ZgXgkH/1FnKEpPrksGGOuBQ5aa18D3gCa57HrLqCKMaaddz3+xpgQa20ymU/yfb3LSxtjAoD1wADvuaUqQEdgs3ddi4ChQAfgTCCtBkZ4j+wwxtQzxpTJpY61QD9jTFVvu6u9+7QJ6GyMqeRdxx0X2J8oYIj3937AOquvLJBiSkdQcrnoDPzVGJMO/AqE56WTtfaUd/pthjGmApl/Ey8CO4HBwGxjzCQgncxwWA60A74i88hmjLX2gHd1HwKRQJS19pR32etALWCb92KFQ0DfXOr4xhjzJPChMcbHu72R1tqNxpgJwBdkngPbRuY05Lm8QeYR5G4yj5wG5mUcRFxk9OJKRERcpCk+ERFxkqb45LLlfZ9T6d8sHmyt/boo6ilIxphx5DwftdRa+3Ru7UWKI03xiYiIk/JzBKUkExGRgnLBTzjROSgREXGSAkpERJykgBIREScpoERExEkKKBERcZICSkREnKSAEhERJymgRETESQooERFxkgJKREScpIASEREnKaBECpAnI/8fWXkxfURKgvx8mrn+ikTyIGL74Xy1H9us8iWqRMRp+rBYEREpnhRQIiLiJAWUiIg4SQElIiJOUkCJiIiTFFAi56FLwEWKjl9RFyDiMj8fk6/LxnXJuEjB0RGUiIg4SQElIiJOUkCJiIiTFFAiIuIkBZSIiDhJASUiIk5SQImIiJMUUCIi4iQFlIiIOEkBJSIiTlJAiYiIkxRQIiLiJAWUiIg4SQElIiJOUkCJiIiTFFAiIuIkBZSIiDhJASUiIk5SQImIiJMUUFKi1KpVi8OHDxfoOvd/8xUv9u/ItN6tGDVqFNbaHG2stUT9429M692Kf/bvxP5vvwLgaOI+WrRoQVhYGCEhIbz66qtZfcaNG0fNmjUpW7ZsrttdtmwZxhi2bNkCQHp6OkOGDKFJkyY0bNiQqVOnFuh+ihQ2BZTI7/Tu1L9y67jnGL1iM/Hx8cR9vjZHm10b1vDz3u8ZvWIztz75HO9OHQNAuSrV+Pzzz4mJiWHTpk1ERESQmJgIwC233MLmzZtz3eYvv/zCjBkzaNOmTdaypUuXcvLkSb7++mu2bt3K7Nmz2bNnT8HvsEghUUCJ8/bs2UODBg0YNmwYjRs3ZtCgQaxZs4b27dtTt25dNm/ezJEjR+jbty+hoaG0bduW2NhYAH7++We6d+9Os2bNGD58eLajm7fffpvWrVsTFhbG8OHDOX36dL5rSz50gJMnfuHapq0wxhAeHs43H6/M0e7b6FU0u3kAxhiCQ1uS9stxkg8dwM+/FKVLlwbg5MmTZGRkZPVp27Yt1atXz3W7Tz31FGPGjOGKK67IWmaM4cSJE3g8HlJTUylVqhTly5fP9z6JuEIBJcXC7t27efjhh4mNjeW7775jwYIFfPbZZ0yfPp1nnnmG8ePH06xZM2JjY3nmmWcIDw8HYOLEidxwww1s376d3r17s3fvXgC+/fZbFi9ezIYNG4iJicHX15f58+cDMGDAAMLCwrJ+ZgzszIyBndn2/uIcdSUfOkD5qoFZt4OCgjh+MClHu+MHk7iq2v/aVagaSPKhAwDs27eP0NBQatasyRNPPEFgYGCO/mfbvn07+/bt4+abb862vF+/fpQpU4bq1asTHBzM6NGjufrqq/MyvCJO8ivqAkTy4rrrrqNJkyYAhISE0KVLF4wxNGnShD179vDjjz/yzjvvAPDHP/6Rn3/+mePHj7N+/Xr+/e9/A3DTTTdRsWJFANauXcvWrVtp1aoVAKmpqVStWhWAxYuzB1HE9nOfs8rtfJMxJreWOduR2a5mzZrExsaSmJhI37596devH9WqVct1exkZGTz66KPMnTs3x32bN2/G19eXxMREjh49SocOHejatSvXX3/9OesXcZkCSoqFM9NgAD4+Plm3fXx88Hg8+PnlfCifCYrcAsNay5AhQ3K9kGDAgAHs2rUr6/bBVA8AN9w9grA/9ePlQV0AaNipJ2363UPywcSstgkJCZSvck2OdVaoGsixn/7X7vjBRMpVyR5CgYGBhISE8Omnn9KvX79cRiHz3NOOHTvo3LkzAAcOHKB3795ERUWxYMECevbsib+/P1WrVqV9+/Zs2bJFASXFlqb45LLQsWPHrCm66OhoKleuTPny5bMtX7lyJUePHgWgS5cuLFu2jIMHDwJw5MgRfvzxRyDzCComJibrZ9SiaEYtiqb5zQPw8fXNut1txFjKV7mGUgFl2Ru7BWstkZGRNOzcM0d9DTv1YPv7i7HWsjd2C1eULU/5Ktdw/KdEUlNTATh69CgbNmygfv3659zPChUqcPjwYfbs2cOePXto27YtUVFRtGzZkuDgYNatW4e1lhMnTrBx40YaNGhQcIMsUsh0BCWXhQkTJjB06FBCQ0MJCAhg3rx5AIwfP54777yT5s2b06lTJ4KDgwFo1KgRU6ZMoXv37mRkZODv78/MmTO59tpr873tvn+fxrLxfyH9ZBp39rmJGu27ArBp2VwA2vS7h/o3dGPXZ2uY3qc1/ldcSb8JMwA4+EMcbcbcjTEGay2jR4/OmsocM2YMCxYsICUlhaCgIIYNG8aECRPOWcfIkSMZOnQojRs3xlqbNR4ixZXJbQ79HPLcUORycr5zUL81tlnlfLU/00ekBMrtZG02muITEREnKaBERMRJCigREXGSAkpERJykgBIREScpoERExEkKKBERcZICSkREnKRPkhApJIu2Ds92e2CL2QB4Mix+Phd8z2KW/LYXKa4UUCJFzM/H5PvTKkRKAk3xiYiIkxRQIiLiJAWUiIg4SQElIiJOUkCJiIiTFFAiIuIkBZSIiDhJASUiIk5SQImIiJMUUCIi4iQFlIiIOEkBJSIiTlJAiYiIkxRQIiLiJAWUiIg4SQElJYonwxZ1CSKSR/rCQilR9OWAIsWHjqBERMRJCigREXGSAkpERJykgBIREScpoERExEkKKBERcZICSkREnKSAEhERJymgRETESQooERFxkgJKREScpIASEREnKaDECbt27SIsLCzrp3z58rz44osAHDlyhG7dulG3bl26devG0aNHz7me5ORkatSowUMPPQTAyZMn6dmzJ40bN2bWrFlZ7f49+TESv4u96Hq3vreI6X1aM71Pa7a+t+ii1vHp26/wwu3tCQ0N5fXht3E0cV/WfSv/OYkX7+jAi3d0IHb18nOuY8mSJTRq1IiQkBDuuusuIHMsW7RoQdOmTfniiy8A8Hg8dO3alZSUlIuqVaQoKKDECfXr1ycmJoaYmBi2bt1KQEAAt956KwARERF06dKF+Ph4unTpQkRExDnX89RTT9GpU6es26tXr6ZFixbExsYyZ84cAJLidmBtBoENQvNVY8rxo1n/rp0znQcjVzPyrQ9ZO2c6qcnHztkvIyODtF+ScywPrN+EkW9/RGxsLI273sLKf04E4LtPPyTxu1j+svBjHoxcxfrImaT9+kuO/vHx8UydOpUNGzawc+fOrECfPXs2ERERLFu2jOnTpwPwyiuvMHjwYAICAvK1zyJFSQElzlm7di21a9fm2muvBWDFihUMGTIEgCFDhvDuu+/m2m/r1q389NNPdO/ePWuZv78/qampeDyerGUfzYqg24ixearF4/HwzSeriHzkbt5+PLOGuC8+pm6bTgRUqMiV5a+ibptO7Pp8XY6+x5ISWPPqP3j+1rbsidmY4/7arW6g1JWZgRHcpAXJBxMBOPh9HNe1aIevnx+lrixD9XohxH2+Nkf/1157jZEjR1KxYkUAqlatmm2fU1JS8Pf359ixY7z33nuEh4fnaZ9FXKGAEucsWrSIO++8M+v2Tz/9RPXq1QGoXr06Bw8ezNEnIyODxx9/nGnTpmVb3q1bNw4cOECbNm0YM2YMUVFRBDZsSvkq15y3hsN7v2fVS5Np2LAhO9a+xw2DR3D/61EAJB9MosI1gVlty1cLJPlgEgCnTp3i649W8OaDdxD52GCuKFeeEXNX0qBD91y3c8aX786nXvsuAFxTL4S4DWs5lZrCiaM/898tGzj+U2KOPnFxccTFxdG+fXvatm3LqlWrABg5ciTPP/88DzzwAH//+9+ZNGkS48aNwxhz3hpEXKMvLBSnnDp1iqioKKZOnZqvfrNmzaJXr17UrFkz23I/Pz8WLFgAQHp6Oj169KDT5Dd5/7mnOH4ggWY3D6BRp57Z+uxY+x4Lx95H5z8/yrZt25i5+2T2jdmc38p75sm/ZcuW/PTrSW6fMIPgJi3yVPvbb7/N/m++4v7XVwBQr92N7N+5nVeH9qJMxcoEh7bEx9c3Rz+Px0N8fDzR0dEkJCTQoUMHduzYQXBwMNHR0QDs3r2bxMREGjRowODBgzl16hSTJ0+mXr16eapNpCjpCEqcsnLlSpo3b061atWyllWrVo2kpMwjlKSkpKyprLN98cUXvPzyy9SqVYvRo0cTGRnJ2LHZp/FmzZrFkCFD2Bu7BT9/f+6MeJ2PX38+x7rqtOnMzaOfZtdnH3H77bcTs+rfpJ9My7q/fLVAjh/43xFN8k+JlPMekb322msEh7ZkyZMjWPniRA5+H3fe/d296ROefvppwl98C79SpbOW3zjsMUYtiubeV5aBtVQKvj5H36CgIPr06YO/vz/XXXcd9evXJz4+PlubcePGMXnyZGbMmMGgQYOYOHEiEydOPG9NIq5QQIlTFi5cmG16D6B3797MmzcPgHnz5tGnT58c/ebPn8/evXvZs2cP06dPJzw8PNvFFEePHuX9998nPDyc9LRUjPEBY/CcPJljXVeULUe7Affy0Pw1PPvss/y4fSMv3N6elS9mPrHXa3cj8RujSU0+RmryMeI3RlOv3Y0AtGnThtvH/5O/LFxH5Vq1eWfSI8wK78n+b7/KsZ3E72JZ/vRooqKiKHt1lazlGadPc+LYEQCS4nZyIP4b6ra9MUf/vn378vHHHwNw+PBh4uLiuP76/wXZJ598Qo0aNahbty4pKSn4+Pjg6+urK/mk2NAUnzgjJSWFjz76iNmzZ2dbPnbsWPr3788bb7xBcHAwS5cuBWDLli28+uqrvP766xdc96RJk3jyyScxxlC33Y18seRN/tm/I236DTlvv2bNmtHnb/8g/WQa//3yUwACKlTkj8Me4+W7uwHwx/seJ6BCxWz9SgeUpVXfu2nV9+7/HUX9Jhc+eHEip1JOcMcdd3Aw1cNV1wQR/uLbnPakM+feWzLXU6Yc/afMwtcv80/1o1ciaHRzR3r37k2PHj348MMPadSoEb6+vkybNo1KlSoBYK1lypQpLFmyBID777+fQYMG4fF4eOWVVy44XiIuMDaX+fRzyHNDEZdFbD+c57Zjm1UusPaLtg7Pdntgi9kXvQ2Ry8AFr9rRFJ+IiDhJASUiIk5SQImIiJMUUCIi4iQFlIiIOEkBJSIiTlJAiYiIkxRQIiLiJAWUiIg4SQElIiJOUkCJiIiTFFAiIuIkBZSIiDhJASUiIk5SQImIiJMUUCIi4iQFlIiIOEkBJSIiTlJAiYiIkxRQIiLiJAWUiIg4SQElIiJOUkCJiIiTFFAiIuIkBZSIiDhJASUiIk5SQImIiJMUUCIi4iQFlEgx48mwl7S9iCv8iroAEckfPx9DxPbDeW4/tlnlS1iNyKWjIygREXGSAkpERJykgBIREScpoERExEkKKBERcZICSkREnKSAEhERJymgRETESQooERFxkgJKREScpIASEREnKaBERMRJCigREXGSAkpERJykgJJiTd91JHL50vdBSbGm70YSuXzpCEpERJykgBIREScpoERExEkKKBERcZICSkREnKSAEhERJymgRETESQooERFxkgJKREScpIASEREnKaBERMRJCigREXGSAkpERJykgBIREScpoERExEkKKBERcZICSkREnKSAEhERJymgRETESQooERFxkgJKREScpIASEREnKaBERMRJCigREXGSAkpERJykgBIREScpoERExEkKKBERcZICSkREnKSAEhERJymgRETESQooERFxkgJKREScpIASEREnKaBERMRJCigREXGSAkpERJykgBIREScpoEQuc54Me0nbi1wqfkVdgEhh2bBgNo0HL+RQqodWt97NDYMeAGDNq//gy+VvUaZiJQC6PzSOBjd0Y0/MJkKH/I0jp30Z+MxsKgdfT+ovx1n4xDCGzlyCMSbHNjp37kzj+54kqFEYAEcT9zLv4UE8svRTft11iB9mfkGpymWw6Rms6TONrsP/SnR0NBNu6c3VNa4lPS2VsldXoeOQv9CwY/cC2W8/H0PE9sN5bj+2WeUC2a7I76WAkhLhwO5v+XL523wfu5UXdibzr4cG0KBDNyoH1wag/aAH6Bg+Mlufz956heXvvEPEmq/YtGwuNz02iXWvPUfnex/JNZzyokydylw/6g+cPukh5tmlNOzQHWr4UiusLffMWABA4q6veeuxIfiXvoI6bTr+vh0XKcY0xSclwqEf4qjZpAUBAQH4+vlxXYs/sHPdB+ft4+PnR2pqKulpqfj6+fHzvh9IPpjE9S3a/+56fEv7UaNhU35O+CHHfYH1m9Dlvsf5Yskbv3s7IsWZjqCkRKhWuyGrZz7Dzz//zKnUFHZ9toagRk2z7v9i8Rtsf38JNRo15abHJnFl+avo/OeHuf/++zl02o/+k2fxwQvj6fbg2Atua/G4B/AvfQUAp9PTMT45Xwd6fj3J/q9j+eN9jwOeHPcHNgxlfeTMi99hkcuAAkqc4cmw+Plc3NTZhVS9vh6d7vkL3bp14yilqV4vBB/fzId/mzvuyQwKY/ho1lT+8/z/0W/CDALrN2Hjxo1EbD/MD1s/p1yVa7AWFjwxDF8/f3o9NpFylarm2NaAp1/NcQ7qjBO7D7Nr0lqMMfQYOo5qtRvA8R05C7a6UEFEASXOyO/JfMjfCf1Wfe/mnYmPELH9MKtfmkL5aoEA2UKm9W2DswUKgLWWdW+8wF0Rr7Hi2bF0feAJjibu5fOFr9HjoXH5qvfMOSiANi3uOWe7xO++pup19fK1bpHLjc5BSYnx65FDABxLSmDnx/8hrOdtACQfOpDVZue6DzKPas6y7b1FNLihK1eWv4r0tFSMjw/Gx4f0tNRLUmdS3E7Wvf48bfv/+ZKsX6S40BGUlBjzRw/l3yeTOerxofcTz3Jl+asAWPnPSSTF7cBgqBhYk77jpmf1SUlJYdv7i/nzzKUA3DDoAeb/dSi+fv4MnDqnwGrbE7ORGXfemHmZecXK3PLXZ3QFn5R4CigpMYa/+T5jm1XOMY04YMqsc/YJCAjgvjnvZt2+rnk7Hlmy/pzto6Ojs62/YmAwjyz9FICy9atQtn6VHH06d+7MhPXf53k/REoKTfGJiIiTFFAiIuIkBZSIiDhJASUiIk5SQImIiJMUUCIi4iQFlIiIOEkBJSIiTlJAiYiIkxRQIiLiJAWUiIg4SQElIiJOUkCJiIiTFFAiIuIkBZSIiDhJASUiIk5SQImIiJMUUCIi4iQFlIhk48mwhdJH5EL8iroAEXGLn48hYvvhfPUZ26zyJapGSjIdQYmIiJMUUCIi4iQFVDFirWXUqFHUqVOH0NBQtm3blmu7cePGUbNmTcqWLZtt+d69e7nxxhtp1qwZoaGhfPDBB1n3xcbG0q5dO0JCQmjSpAlpaWkAbN26lSZNmlCnTh1GjRqFtZnnGiZMmECNGjUICwsjLCws27qkZNi1YS3P3dqWab1bERERkeP+9evX07x5c/z8/Fi2bFnW8piYmKzHWmhoKIsXL8667+WXX6ZOnToYYzh8OPs0Y3R0NGFhYYSEhNCpUycA0tLSaN26NU2bNiUkJITx48dfor2VoqCAKkZWrlxJfHw88fHxzJkzhxEjRuTa7pZbbmHz5s05lk+ZMoX+/fuzfft2Fi1axIMPPgiAx+Ph7rvv5tVXX2Xnzp1ER0fj7+8PwIgRI5gzZ07WdletWpW1vkcffZSYmBhiYmLo1avXJdhjcVXG6dNEPTuWoS8t4tF3NrBw4UK++eabbG2Cg4OZO3cud911V7blAQEBREZGsnPnTlatWsUjjzzCsWPHAGjfvj1r1qzh2muvzdbn2LFjPPjgg0RFRbFz506WLl0KQOnSpVm3bh1fffUVMTExrFq1io0bN17CPZfCpIAqRlasWEF4eDjGGNq2bcuxY8dISkrK0a5t27ZUr149x3JjDMnJyQAcP36cwMBAAD788ENCQ0Np2rQpAJUqVcLX15ekpCSSk5Np164dxhjCw8N5991381yvruy6fO3bsY1KQbW4OqgWfv6lGDhwICtWrMjWplatWoSGhuLjk/1ppl69etStWxeAwMBAqlatyqFDhwBo1qwZtWrVyrG9BQsWcNtttxEcHAxA1apVgczH9JmZgvT0dNLT0zHGFOi+StHRVXzFyP79+6lZs2bW7aCgIPbv359rGOVmwoQJdO/enZdeeokTJ06wZs0aAOLi4jDG0KNHDw4dOsTAgQMZM2YM+/fvJygoKMf2znj55ZeJjIykZcuWPPfcc1SsWDHb9vJ7NZiuBCs+kg8lUeGaGlm3g4KC2LRpU77Xs3nzZk6dOkXt2rXP2y4uLo709HQ6d+7ML7/8wsMPP0x4eDgAp0+fpkWLFuzevZuRI0fSpk2bfNchbtIRVDFy5vzP2fLzanHhwoXcc889JCQk8MEHHzB48GAyMjLweDx89tlnzJ8/n88++4zly5ezdu3a825vxIgR/Pe//yUmJobq1avz+OOPX/yOSfHzOx+LAElJSQwePJh//etfOY6yfsvj8bB161b+85//sHr1aiZPnkxcXBwAvr6+xMTEkJCQwObNm9mxY0e+6hB3KaAcN3PmzKwLEQIDA9m3b1/WfQkJCVnTdHnxxhtv0L9/fwDatWtHWloahw8fJigoiE6dOlG5cmUCAgLo1asX27ZtIygoiISEhFy3V61aNXx9ffHx8eG+++7L9ZyXXL7KVw3k+IH/HU3n97GYnJzMTTfdxJQpU2jbtu0F2wcFBdGzZ0/KlClD5cqV6dixI1999VW2NldddRWdO3fOdp5UijcFlONGjhyZdSFC3759iYyMxFrLxo0bqVChQp6n9yDzpPXatWsB+Pbbb0lLS6NKlSr06NGD2NhYUlJS8Hg8fPLJJzRq1Ijq1atTrlw5Nm7ciLWWyMhI+vTpA5Dt3Nfy5ctp3Lhxwe64OC0opBmH9/3Akf0/4kk/xaJFi+jdu3ee+p46dYpbb72V8PBw7rjjjjz16dOnD59++ikej4eUlBQ2bdpEw4YNOXToUNYFFqmpqaxZs4YGDRpc9H6JWxRQxUivXr24/vrrqVOnDvfddx+zZs3Kui8sLCzr9zFjxhAUFERKSgpBQUFMmDABgOeee47XXnuNpk2bcueddzJ37lyMMVSsWJHHHnuMVq1aERYWRvPmzbnpppsAeOWVVxg2bBh16tShdu3a/OlPf8raRpMmTQgNDeXjjz/mhRdeKLyBkCLn6+dH7yem8ubI/rxwe3v69+9PSEgI//d//0dUVBQAX375JUFBQSxdupThw4cTEhICwJIlS1i/fj1z587Nmh2IiYkBYMaMGVlH7qGhoQwbNgyAhg0b0rNnT0JDQ2ndujXDhg2jcePGJCUlceONNxIaGkqrVq3o1q0bN998c9EMihQ4k9t5hnPQJVmSb/m9SOJiPmLnUm6jINsv2jo82+2BLWYXeU0F0f5MH5F8uuBJSx1BiYiIkxRQIiLiJAWU5JneeCvnkt/Hhh5Lkhd6o67kmd54K+eix4ZcCjqCEhERJymgRETESQooERFxkgJKREScpDfqlmCeDIufT/4+4PNcJ8L/1rxKQZR0ebv9N7ffKZIqCszUbYeyfs/vm3tHN62Ur8fexTxWxXkX/A/VVXwlmK68kqKix57khab4RETESQqoy4Te+CiXs/w+vtMv4u9Bf0Pu0RSfo/I7557fKRMo2GmTs89HnL1+fVjGTqtnAAAEJklEQVTseYxzsKaLaF8YLmZKML/7MLpppXy113mxS08BVUgudeBojl7k98nv31x+Ay09w+Kfj+eA/LaHyy8083wV38SJE1cBl+OzYCCQWNRFFFMau4ujcbt4GruL59rYHR4/fnzP87aw1pbonwkTJtiirqG4/mjsNG4au+LzUxzHThdJiIiIkxRQMLGoCyjGNHYXR+N28TR2F6/YjV1+PklCRESk0OgISkREnKSAEhERJymgRETESSUuoIwxVxtjPjLGxHv/rXietuWNMfuNMS8XZo2uysvYGWPCjDFfGGN2GmNijTEDiqJWFxhjehpjdhljdhtjxuZyf2ljzGLv/ZuMMbUKv0o35WHsHjPGfON9jK01xlxbFHW66EJjd1a7fsYYa4xpWZj15UeJCyhgLLDWWlsXWOu9fS6TgU8KpariIS9jlwKEW2tDgJ7Ai8aYqwqxRicYY3yBmcCfgEbAncaYRr9pdi9w1FpbB3gBeLZwq3RTHsduO9DSWhsKLAP+UbhVuimPY4cxphwwCthUuBXmT0kMqD7APO/v84C+uTUyxrQAqgEfFlJdxcEFx85aG2etjff+nggcBEril0W1BnZba7+31p4CFpE5fmc7ezyXAV2MMZfP59RcvAuOnbX2Y2ttivfmRiCokGt0VV4ed5D54vsfQFphFpdfJTGgqllrkwC8/1b9bQNjjA/wHPDXQq7NdRccu7MZY1oDpYD/FkJtrqkB7DvrdoJ3Wa5trLUe4DiQvw94uzzlZezOdi+w8pJWVHxccOyMMc2Amtba9wuzsItxWX5YrDFmDXBNLneNy+MqHgQ+sNbuK2kvaAtg7M6spzrwFjDEWptRELUVM7k9cH77psO8tCmJ8jwuxpi7gZZAp0taUfFx3rHzvvh+AbinsAr6PS7LgLLWdj3XfcaYn4wx1a21Sd4n0YO5NGsHdDDGPAiUBUoZY3611p7vfNVloQDGDmNMeeA/wJPW2o2XqFTXJQA1z7odRM4P6jzTJsEY4wdUAI4UTnlOy8vYYYzpSuYLp07W2pOFVJvrLjR25YDGQLT3xfc1QJQxpre1dkuhVZlHJXGKLwoY4v19CLDitw2stYOstcHW2lrAaCCyJIRTHlxw7IwxpYDlZI7Z0kKszTVfAnWNMdd5x2QgmeN3trPHsx+wzuqjXSAPY+edppoN9LbW5vpCqYQ679hZa49baytba2t5n982kjmGzoUTlMyAigC6GWPigW7e2xhjWhpjXi/SytyXl7HrD3QE7jHGxHh/woqm3KLjPaf0ELAa+BZYYq3daYyZZIzp7W32BlDJGLMbeIzzX1FaYuRx7KaRObux1PsY+234l0h5HLtiQ5/FJyIiTiqJR1AiIlIMKKBERMRJCigREXGSAkpERJykgBIREScpoERExEkKKBERcdL/A0hLZpnHvnNUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "s\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using advi...\n",
      "Average Loss = -1,090.7:  91%|█████████▏| 45648/50000 [00:37<00:03, 1214.60it/s]  \n",
      "Convergence achieved at 45700\n",
      "Interrupted at 45,699 [91%]: Average Loss = 3.5883e+06\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [value_sd, is_recovered|dataset_offset, is_recovered|dataset_sd, 1|dataset_offset, 1|dataset_sd, is_recovered, Intercept]\n",
      "Sampling 2 chains: 100%|██████████| 61000/61000 [14:00<00:00, 72.60draws/s]\n",
      "There was 1 divergence after tuning. Increase `target_accept` or reparameterize.\n",
      "The acceptance probability does not match the target. It is 0.890592701128019, but should be close to 0.8. Try to increase the number of tuning steps.\n",
      "The acceptance probability does not match the target. It is 0.887196398793421, but should be close to 0.8. Try to increase the number of tuning steps.\n",
      "The number of effective samples is smaller than 25% for some parameters.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                                        mean        sd  hpd0.95_lower  \\\n",
      "is_recovered|dataset_sd             0.000413  0.000507   1.175473e-08   \n",
      "is_recovered|dataset[folke2016]     0.000017  0.000392  -7.983618e-04   \n",
      "is_recovered|dataset[krajbich2010]  0.000021  0.000379  -7.571990e-04   \n",
      "is_recovered|dataset[krajbich2011] -0.000074  0.000390  -9.202820e-04   \n",
      "is_recovered|dataset[tavares2017]   0.000037  0.000392  -7.435034e-04   \n",
      "1|dataset[folke2016]               -0.002512  0.001158  -4.813972e-03   \n",
      "1|dataset[krajbich2010]             0.000235  0.001150  -2.105002e-03   \n",
      "1|dataset[krajbich2011]            -0.000627  0.001154  -2.960165e-03   \n",
      "1|dataset[tavares2017]              0.002898  0.001158   5.292771e-04   \n",
      "is_recovered                       -0.000285  0.000384  -9.891638e-04   \n",
      "Intercept                           0.008369  0.001141   6.124589e-03   \n",
      "value_sd                            0.001721  0.000080   1.567019e-03   \n",
      "1|dataset_sd                        0.002333  0.000750   1.096532e-03   \n",
      "\n",
      "                                    hpd0.95_upper  effective_n  gelman_rubin  \n",
      "is_recovered|dataset_sd                  0.001286        20949      0.999994  \n",
      "is_recovered|dataset[folke2016]          0.000801        30883      0.999996  \n",
      "is_recovered|dataset[krajbich2010]       0.000785        28394      0.999991  \n",
      "is_recovered|dataset[krajbich2011]       0.000643        30992      1.000021  \n",
      "is_recovered|dataset[tavares2017]        0.000857        30220      0.999984  \n",
      "1|dataset[folke2016]                    -0.000172        10257      1.000071  \n",
      "1|dataset[krajbich2010]                  0.002523        10179      1.000030  \n",
      "1|dataset[krajbich2011]                  0.001686        10416      1.000061  \n",
      "1|dataset[tavares2017]                   0.005192        10366      1.000057  \n",
      "is_recovered                             0.000446        29690      0.999991  \n",
      "Intercept                                0.010725         9904      1.000045  \n",
      "value_sd                                 0.001882        55253      0.999984  \n",
      "1|dataset_sd                             0.003825        17448      1.000036  \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XtYlVWix/Hvgo2UFea1QLJUNBVFCE3tpuWAljNkZWpjqd3TJudYOdrxlJec1KzJnC5m2ahjCZPOpJm31Dxd1USQvKSQmoGUqSinAU1gnT/Y7gFBBQNZ4O/zPPth7/dda+21Frh/+13vu7fGWouIiIhr/Kq6AyIiIqVRQImIiJMUUCIi4iQFlIiIOEkBJSIiTlJAiYiIkxRQIiLiJAWUVGvGmC+qug9nkzFmsDHmlVPsDzTGJBhj0owx64wxV5y93olULAWUVGvW2mt+TX1jjH9F9cWRPtwPZFlrw4CXgMkV2LbIWaWAkmrNGPOz92ewMeYTY0yyMWazMeb6U9Uxxow3xqwDuhhjoo0x/2uMSTTGLDfGBHvLhRljVhpjNhljNhpjmptCU7zP8bUxpp+3bIIx5pYizzHLGHOHMcbfW/4rY0yKMeZh7/5uxpiPjTHvAl97t91tjFnvHcMbx4PLGHOvMWaHMeZ/gWtPMyW3ArO99+cD3Y0x5gymVqTKKaCkpvg9sNxaGwm0B5JPUfYCYLO1thOwDvgr0MdaGw28DfzZW+4d4FVrbXvgGiATuB04/hy/AaZ4Ay0eOB5WtYDuwBIKj2gOW2s7Ah2BB40xTb3tXw2Mtta2Mca09ta/1juGfGCAt+1xFAZTDNDmNPPQGPgewFqbBxwG6p+mjoiTPFXdAZEK8hXwtjEmAHjfWnuqgMoHFnjvXwm0BT7yHmj4A5nGmIuAxtbafwFYa48AGGOuA+ZZa/OBH71HNR2BpcA0Y0wg0BP4xFqba4yJBSKMMX28z1cHaAH8Aqy31u7ybu8ORANfeftxPrAP6ASssdb+5H3+BKDlKcZW2tGSvnBTqiUFlNQI1tpPjDE3AL2Avxtjplhr55yk+BFvwEDhC/oWa22XogWMMUEnqVvqcpm19ogxZg3Qg8IjoXlFyj9mrV1+QvvdgH+f0O5sa+1TJ5TrTfkCJh24DEg3xngoDMSD5agv4gwt8UmNYIy5HNhnrX0TmAlcVcaq24GGxpgu3nYCjDHh1tpsCl/ke3u3BxpjagOfAP2855YaAjcA671txQP3AtcDxwNpOTDEe2SHMaalMeaCUvqxCuhjjGnkLVfPO6Z1QDdjTH1vG3eeZjyLgEHe+32A1Vb/ZYFUUzqCkpqiGzDCGHMM+BkYWJZK1tpfvMtv04wxdSj8NzEV2ALcA7xhjBkPHKMwHP4FdAE2UXhk8ydr7Q/e5lYAc4BF1tpfvNveAq4ANnovVvgJ6F1KP7YaY/4HWGGM8fM+36PW2rXGmLHAlxSeA9tI4TLkycyk8AgyjcIjp/5lmQcRFxm9uRIRERdpiU9ERJykJT6psbyfcwo8YfM91tqvq6I/FckYM5qS56Pes9b+ubTyItWRlvhERMRJ5TmCUpKJiEhFOe03nOgclIiIOEkBJSIiTlJAiYiIkxRQIiLiJAWUiIg4SQElIiJOUkCJiIiTFFAiIuIkBZSIiDhJASUiIk5SQImIiJMUUCKnkFdQvq+gLG95ETm58nybuf7lyTlpUtL+MpcdFdWgEnsiUqPoy2JFRKR6UkCJiIiTFFAiIuIkBZSIiDhJASUiIk5SQImIiJMUUCIi4iQFlIiIOEkBJSIiTlJAiYiIkxRQIiLiJAWUiIg4SQElIiJOUkCJiIiTFFAiIuIkBZSIiDhJASUiIk5SQImIiJMUUCIi4iQFlIiIOEkBJSIiTlJAiYiIkxRQIiLiJAWUiIg4SQElIiJOUkCJiIiTFFBS7V1xxRXs37+/QttMTEykXbt2hIWFsej5p7DWlihjrWXR808xJa4jL/ftSsa2Tb59s2fPpkWLFrRo0YLZs2eX2u6wYcN87Y4YMYJWrVoRERHBbbfdxqFDhwDYvXs3559/PpGRkURGRvLII49U6DhFXKaAEinFkCFDmDFjBqmpqRzYs5MdX6wqUWb75ys5sGcnTy5cz23/8yLvT/wTAAcPHmTcuHGsW7eO9evXM27cOLKyskq0m5qayrJlywCIiYlh8+bNpKSk0LJlSyZOnOh7nubNm5OcnExycjLTp08/C6MXcYMCSqrE7t27adWqFQ888ABt27ZlwIABrFy5kmuvvZYWLVqwfv16Dh48SO/evYmIiKBz586kpKQAcODAAWJjY4mKiuLhhx8udnQzd+5crr76aiIjI3n44YfJz88vd98yMzPJzs6mS5cuGGOI+m0/tn68tES5bWuWEfXbfhhjaBLRgSP/d5jMzEyWL19OTEwM9erVo27dusTExLBs2bIS7Q4cOJD3338fgNjYWDweDwCdO3cmPT39TKZVpEZRQEmVSUtL449//CMpKSl88803vPvuu3z22We88MILPPfcc4wZM4aoqChSUlJ47rnnGDhwIADjxo3juuuuIykpibi4OPbs2QPAtm3bSEhI4PPPPyc5ORl/f3/eeecdAPr16+dbJit6mzNnTol+ZWRkEBoa6ntcp1Ewh/dllih3eF8mF18SUqRcCBkZGWRkZHDZZZf5toeGhvq2F233+PYTvf3229x8882+x7t27SIqKoquXbvy6aeflnl+Rao7T1V3QM5dTZs2pV27dgCEh4fTvXt3jDG0a9eO3bt3891337FgwQIAbrrpJg4cOMDhw4f55JNP+Oc//wlAr169qFu3LgCrVq0iMTGRjh07ApCbm0ujRo0ASEhIKHO/SjvfZIwprWSp5U5Wvyzt/vnPf8bj8TBgwAAAgoOD2bNnD/Xr1ycxMZHevXuzZcsWgoKCyjgakepLASVVJjAw0Hffz8/P99jPz4+8vDzfkldRx1/QSwsMay2DBg0qdv7muH79+rF9+/YS2x9//HEGDBhAdHQ0AHFxcQwZMqTYEtvhfZkENby0RN06jUI49OPeIuX2EhISQmhoKGvWrPFtT09Pp1u3boSGhhZrNz09nZCQ/xyBzZ49m8WLF7Nq1Srf+AIDA33zEh0dTfPmzdmxYwcdOnQo0R+RmkZLfOKsG264wbdEt2bNGho0aEBQUFCx7UuXLvVdgNC9e3fmz5/Pvn37gMKLFb777jug8Ajq+IUGRW8DBw7E39/f93j8+PEEBwdz0UUXsXbtWqy1JC1OoHW3niX617prD5IWJ2CtZU/KBs67MIjg4GB69OjBihUryMrKIisrixUrVtCjR48S7c6ZM4dbb70VgGXLljF58mQWLVpE7dq1fc/x008/+c6j7dy5k9TUVJo1a1ZJMy7iFh1BibPGjh3LvffeS0REBLVr1/Zdrj1mzBjuuusurrrqKrp27UqTJk0AaNOmDRMmTCA2NpaCggICAgJ49dVXufzyy8v93K+//jqDBw8uXCbs0I0rr/0NAOvmzwKgU5/BXHldDNs/W8kLt15NwHnn02fsNADq1avH008/7VtqfOaZZ6hXr16Jdm+++WbfuaY//OEPHD16lJiYGKDwQonp06fzySef8Mwzz+DxePD392f69Om+tkRqOlPauvhJlLmgSE0yKansn7EaFdWgEnsiUqOUdmK3GC3xiYiIkxRQIiLiJAWUiIg4SQElIiJOUkCJiIiTFFAiIuIkBZScU/IK9GkJkepCH9SVc4rHz+hzTSLVhAJK5AzEJz5c7HH/6DeqqCciNZeW+ERExEkKKBERcZICSkREnKSAEhERJymgRETESQooERFxkgJKREScpIASEREnKaBERMRJCigREXGSAkpERJykgBIREScpoERExEkKKBERcZICSkREnKSAEhERJymgRETESQooERFxkgJKREScpIASEREnKaBERMRJCigREXGSAkpERJykgBIREScpoERExEkKKBERcZICSkREnKSAEhERJymgRETESQooERFxkgJKKs1LL71EeHg4bdu25a677uLIkSMA3H///bRv356IiAj69OnDzz//XKLuRx99RHR0NO3atSM6OprVq1cDcPToUXr27Enbtm157bXXfOUfeughkpKSzriviR/E88KtV/PCrVeT+EH8GbWxK/ELrrrqKjweD/Pnzy+2b8+ePcTGxtK6dWvatGnD7t27S9Tfs2cPN954I1FRUURERLBkyRIAPv/8cyIiIujYsSNpaWkAHDp0iB49emCtPaO+ilQHCiipFBkZGUybNo0NGzawefNm8vPziY8vfOF/6aWX2LRpEykpKTRp0oRXXnmlRP0GDRrwwQcf8PXXXzN79mzuueceAJYvX050dDQpKSnMmDEDgE2bNlFQUEBUVFS5+phzOMv3c9WMFxg6ZzmP/n0Fq2a8QG72oZPWKygoID/nWIntFweHMmvWLH7/+9+X2Ddw4EBGjBjBtm3bWL9+PY0aNSpRZsKECfTt25ekpCTi4+MZOnQoAC+++CILFizgueee4/XXXwfg2Wef5b//+78xxpRrzCLViaeqOyA1V15eHrm5uQQEBJCTk0NISAgAQUFBAFhryc3NLfVFtmjYhIeHc+TIEY4ePUpAQAC5ubnk5eX59j/99NNMnz69zH3a+r/L2PCvuRz5OZuH3lrEji8/pkWnrtSuUxeAFp26sv2L1UT2vL1Y3UOZ6WxY+C7JS+cT1PtSgiKCi+2vG9KEiIgG+PkVf9+3detW8vLyiImJAeDCCy8stW/GGLKzswE4fPiwb76OjzknJ4eAgAC+/fZbMjIy6Nq1a5nGLFJd6QhKKkXjxo158sknadKkCcHBwdSpU4fY2Fjf/nvvvZdLL72Ub775hscee+yUbS1YsICoqCgCAwOJiYnhhx9+oFOnTvzpT39i0aJFREdH+17MTyYtLY2nnnqK1q1bs3nVB1x3zxAeemsRANn7Mqlz6X/qB10SQva+TAB++eUXvv5oIW8PvZM5j9/DeRcFMWTW0hLhdCo7duzg4osv5vbbbycqKooRI0aQn59fotzYsWOZO3cuoaGh3HLLLfz1r38F4KmnnuKhhx5i6tSp/OEPf2D06NE8++yzZX5+kepKASWVIisri4ULF7Jr1y727t3Lv//9b+bOnevb/7e//Y29e/fSunVrEhISTtrOli1bGDlyJG+88QYAHo+Hd999l6SkJO68806mTp3KE088weOPP06fPn1YtGhRiTYWLFhAq1atqFWrFhs3bqTv+FdpFn3tfwqUch7n+FFdhw4dWDn9eX4zZBTD5n3MdQMe4YK69cs1F3l5eXz66ae88MILfPXVV+zcuZNZs2aVKDdv3jwGDx5Meno6S5Ys4Z577qGgoIDIyEjWrl3Lxx9/zM6dOwkJCcFaS79+/bj77rv58ccfy9UfkepCASWVYuXKlTRt2pSGDRsSEBDA7bffzhdffFGsjL+/P/369WPBggWltpGens5tt93GnDlzaN68eYn9r732GoMGDeLLL7+kVq1aJCQkMGHChBLlYmJiePnll/nwww+54447SF72T44dPeLbH3RJCId/2Ot7nP3jXi5qeCkAb775Jk0iOvCP/xnC0qnj2LdzR7nnIjQ0lKioKJo1a4bH46F3795s3LixRLmZM2fSt29fALp06cKRI0fYv3+/b7+1lgkTJvD0008zbtw4xo0bx9133820adPK3SeR6kABJZWiSZMmrF27lpycHKy1rFq1itatW2Ot9V2JZq3lgw8+oFWrViXqHzp0iF69ejFx4kSuvfbaEvuzsrJYvHgxAwcOJCcnBz8/P4wxvisFiwoKCuLRRx9lw4YNTJ48me+S1vLSHdeydOo4AFp2uZHUtWvIzT5EbvYhUteuoWWXGwHo1KkTd4x5mcfmrabBFc1ZMP6/eG1gT3K+yyrzXHTs2JGsrCx++uknAFavXk2bNm1KnbNVq1YBsG3bNo4cOULDhg19+2fPnk2vXr2oW7eub8x+fn7k5OSUuS8i1Ykpx2Wqup5VymXMmDEkJCTg8XiIiorirbfeIiAggOuvv57s7GystbRv357XX3+doKAgFi1axIYNGxg/fjwTJkxg4sSJtGjRwtfeihUrfFe/DR8+nN69e9O1a1eOHDlCXFwcGRkZPPLII6c9pzUpaT/Hjh7h268+pdV1hRcubHj/HT5+eyoAN94/nA63Fl6JNyqqAZOS9herv2/nDpZsHcN5wUG+bf2j3+D7LUksHHUvWVlZnHfeeVx66aVs2bIFKLxs/oknnsBaS3R0NDNmzKBWrVo888wzdOjQgbi4OLZu3cqDDz7Izz//jDGG559/3nfeLicnh169erFixQoCAgL49NNPGTp0KLVq1WLevHm0bNny1/yqRKrCaS9BVUDJOefEwDmV0gIKID7x4WKP+0e/4SsvImVy2oDSEp+IiDhJASUiIk5SQImIiJMUUCIi4iQFlIiIOEkBJSIiTlJAiYiIkxRQIiLiJAWUiIg4SQElIiJOUkCJiIiTFFAiIuIkBZSIiDhJASUiIk5SQImIiJMUUCIi4iQFlIiIOEkBJSIiTlJAiYiIkxRQIiLiJAWUiIg4SQElIiJOUkCJiIiTFFAiIuIkBZSIiDhJASUiIk5SQImIiJMUUCIi4iQFlIiIOEkBJSIiTlJAiYiIkxRQIiLiJAWUiIg4SQElIiJOUkCJiIiTFFAiIuIkBZSIiDhJASUiIk5SQImIiJMUUCIi4iQFlIiIOEkBJSIiTlJAiYiIkxRQIiLiJAWUiIg4SQElIiJOUkCJiIiTFFAiFSivwJ6VOiLnAk9Vd0CkJvH4GSYl7S9XnVFRDSqpNyLVm46gRETESQooERFxkgJKREScpIASEREnKaBERMRJCigREXGSAkpERJykgBIREScpoERExEkKKBERcZICSkREnKSAEhERJymgRETESQooERFxkgJKREScpIASEREnKaBERMRJCigREXGSAkpERJykgBIREScpoERExEkKKBERcZICSkREnKSAEhERJymgRETESQooERFxkgJKREScpIASEREnKaBERMRJCigREXGSAkpERJykgBIREScpoERExEkKKBERcZICSkREnKSAEhERJymgpNp6+eWXadu2LeHh4UydOtW3fezYsTRu3JjIyEgiIyNZsmQJAJ9//jkRERG8cncM+/fsBCD3/w7z9tA7sdaW+hzdunUjfWuy73HW3j1MvfN6AH7e/hNfD1vE9vGr+Mvt17DyjSkA7NzwOWNvaMa0u27kxds688b9v2PbJysqZQ5EajJPVXdA5Exs3ryZN998ky/WrqP2eYH07NmTXr160aJFCwCGDx/Ok08+WazOiy++yIIFC5i0chPr5s+i1+PjWf3mi3S7/78wxpxRPy4Ia0CzYddwe5u/MK3/jSQm9gXgisjODJ72LgB7t3/N3x8fREDgeYR1uuFXjFrk3KKAkmpp27ZtdO7cmaALL2BS0n5My44Me2UuXQc/xmeZOdSqbchL2l+sTtrPBeTm5nLsSC7+Hg8Hvt9F9r5MmkVf+6v7U+v8C2jcuj3ffvstEFhsX8iV7ej+4BN8+Y+ZCiiRclBASbXUtm1bRo8ezYEDB/glN4ftn60ktE173/4vE2aStPgfNG7Tnl6Pj+f8oIvpdt8feeihh/gp30PfZ19jyUtjiBk66rTPlTD6EQICzwMg/9gxjF/JlfF/HzrInq8TCQ+fwMov0krsD2kdwSdzXv0VIxY59yigpFpq3bo1I0eOJCYmhiwCCW4Zjp9/4Z9zpzsHc9ODT4AxfPTaRD78yzP0GTuNkCvbsXbtWiYl7WdX4hdc1PBSrIV3Rz6AvyeAWx4fx0X1G5V4rn5/nk5om0ig8BzU7D8O8O37d9p+to9fxU8X3Em3e4cRHh4OpQQUJznHJSInp4skpNq6//772bhxIw/P/IDaQRdTv0kzAC6q3wg/f3/8/Py4+vZ7SN+SVKyetZbVM1+i+4NPsGrGFH7zyEgib+nDF/PeLHcfLghrwJXPdOexd1fRqc/gk5bb+83XNGrastzti5zLFFBSbe3btw+AQ5npbPn4QyJ73g5A9k8/+MpsWb2ES5q3KlZv4wfxtLruN5wfdDHHjuRi/Pwwfn4cO5JbKf3M3LGF1W/9hc5976uU9kVqKi3xSbV1xx13cODAAbLy/IgbOZnzgy4GYOnL48ncsRmDoW7IZfQe/YKvTk5ODhsXJ3Dfq+8BcN2AR3hnxL34ewLoP3FGhfVtd/Japt11I8eO5HJh3Qb8bsRzukBCpJzMyT7/UQotoouTJp1wtd6pjIpqUCHl4xMfLva4f/QbZ9T+8Toi56DTfrZDS3wiIuIkBZSIiDhJASUiIk5SQImIiJMUUCIi4iQFlIiIOEkBJSIiTlJAiYiIkxRQIiLiJAWUiIg4SQElIiJOUkCJiIiTFFAiIuIkBZSIiDhJASUiIk5SQImIiJMUUCIi4iQFlIiIOEkBJSIiTlJAiYiIkxRQIiLiJAWUiIg4SQF1jrLWMmzYMMLCwoiIiGDjxo2llktMTKRdu3aEhYUxbNgwrLUAHDx4kJiYGFq0aEFMTAxZWVmnbXfkyJG0bduWtm3bkpCQ4Nu+evVqrrrqKtq2bcugQYPIy8urxJG7bdmyZVx55ZWEhYUxadKkEvuPHj1Kv379CAsLo1OnTuzevdu3b+LEiYSFhXHllVeyfPnyMrf52GOPceGFF/oez5o1i4YNGxIZGUlkZCRvvfVWxQ5SpIwUUOeopUuXkpqaSmpqKjNmzGDIkCGllhsyZAgzZszwlV22bBkAkyZNonv37qSmptK9e3ffC9/J2v3www/ZuHEjycnJrFu3jilTppCdnU1BQQGDBg0iPj6ezZs3c/nllzN79uyzMwmOyc/P59FHH2Xp0qVs3bqVefPmsXXr1mJlZs6cSd26dUlLS2P48OGMHDkSgK1btxIfH8+WLVtYtmwZQ4cOJT8//7RtbtiwgUOHDpXoS79+/UhOTiY5OZkHHnigcgcuchIKqHPUwoULGThwIMYYOnfuzKFDh8jMzCxWJjMzk+zsbLp06YIxhoEDB/L+++/76g8aNAiAQYMGFdteWrtbt26la9eueDweLrjgAtq3b8+yZcs4cOAAgYGBtGzZEoCYmBgWLFhwFmei6uUVFB6Vrl+/nrCwMJo1a0atWrXo378/CxcuLFZ24cKFDLhnIAB9+vRh1apVWGtZuHAh/fv3JzAwkKZNmxIWFsb69etP2WZ+fj4jRozg+eefP7sDFikjBdQ5KiMjg8suu8z3ODQ0lIyMjBJlQkNDSy3z448/EhwcDEBwcDD79u07Zbvt27dn6dKl5OTksH//fj7++GO+//57GjRowLFjx9iwYQMA8+fP5/vvv6+cQTvK42eYlLSfVz7bxsHaDZmUtJ9JSftJyqvDouQ03+NJSfvZ9O0eml7epLCex0OdOnU4cODASef9VL/nV155hbi4ON/vsagFCxYQERFBnz59zrnfh7jDU9UdkKpx/FxSUcaYcpcpa7uxsbF89dVXXHPNNTRs2JAuXbrg8XgwxhAfH8/w4cM5evQosbGxeDzn6J/lr/idnGx7QUFBqdv37t3Le++9x5o1a0rs/93vfsddd91FYGAg06dPZ9CgQaxevbocAxGpGDqCOoe8+uqrvhPfISEhxd4Zp6enExISUqx8aGgo6enppZa55JJLfEuCmZmZNGrUyFfnZO2OHj2a5ORkPvroI6y1tGjRAoAuXbrw6aefsn79em644Qbf9nNNUKMQDv/wn6PY7H17CWp4abEydRoF++Y3Ly+Pw4cPU69evZPO+8m2JyUlkZaWRlhYGFdccQU5OTmEhYUBUL9+fQIDAwF48MEHSUxMrLQxi5yKAuoc8uijj/pOfPfu3Zs5c+ZgrWXt2rXUqVOnxFJPcHAwF110EWvXrsVay5w5c7j11lsBiIuL813MMHv27GLbS2s3Pz+fAwcOAJCSkkJKSgqxsbEAvuXBo0ePMnnyZB555JGzMh+uCQ2PYv/3uziY8R15x35h0/L3ad21Z7Eyrbv29M37/PnzuemmmzDGEBcXR3x8PEePHmXXrl2kpqZy9dVX07FjR1JTU9m1axe//PIL8fHxxMXF0atXL3744Qd2797N7t27qV27NmlpaQDFzkUuWrSI1q1bn71JECniHF1LkVtuuYUlS5YQFhZG7dq1+dvf/ubbFxkZSXJyMgCvv/46gwcPJjc3l5tvvpmbb74ZgFGjRtG3b19mzpxJkyZNeO+9907Z7rFjx7j++usBCAoKYu7cub6lvClTprB48WIKCgoYMmQIN91001mbB5f4ezzEjZzI24/2xRYU0CHuLi5p3oqPXp9E4zaRtOnakw69B5Dy4nDCwsKoV68e8fHxAISHh9O3b1/atGmDx+Ph1Vdfxd/fHyg819SjRw/y8/O57777CA8PP2U/pk2bxqJFi/B4PNSrV49Zs2ZV9tBFSmVKW7s+iTIXFDmbJiXtL3PZUVENKqR8fOLDxR73j37jjNo/0z6J1ACnPqGNlvhERMRRCigREXGSAkpERJykgBJnHP9GBRER0FV84pDj36hQHrpgQKTm0hGUiIg4SQElIiJOUkCJiIiT9EFdccqkpP08dVXDqu7G6d1xwuMK/B9CJm786ZT7dd5Nagh9UFekpinv1Y66OlKqK13FJ5Umr8Di8TvtmyQpp/Je7agjLqmuFFBSafRCKiK/hgJKnHO6czBFVdSXv/7q8qPPrP0zrVMe5T2S1ZGvuEIBdQ4r7wvRsQJLgF64qp3yHsk+2b5+uZ9DoSaVocxX8Y0bN24ZUNVrMCHA3iruQ02geaw4msuKo7msONVhLvePGTOm5ylLWGurzW3s2LG2qvtQE26aR82lizfNpebyxJsuMxcRESdVt4AaV9UdqCE0jxVHc1lxNJcVp0bMZXm+SUJEROSsqW5HUCIico5QQImIiJMUUCIi4qQqDyhjTD1jzEfGmFTvz7onKTfIWybVGDOoyPZoY8zXxpg0Y8w0Y4w5od6TxhhrjKnqz3BVusqaS2PMFGPMN8aYFGPMv4wxF5+tMZ1txpiexpjt3jkYVcr+QGNMgnf/OmPMFUX2PeXdvt0prTGvAAADvElEQVQY06OsbdZEFT2PxpjLjDEfG2O2GWO2GGP+ePZGU7Uq42/Su8/fGJNkjFlc+aM4Q1V9nTvwPDDKe38UMLmUMvWAnd6fdb3363r3rQe6UPjV7UuBm4vUuwxYDnwHNKjqsVbXuQRiAY/3/uTS2q0JN8Af+BZoBtQCNgFtTigzFJjuvd8fSPDeb+MtHwg09bbjX5Y2a9qtkuYxGLjKW+YiYEdNn8fKmssi9R4H3gUWV/U4T3ar8iMo4FZgtvf+bKB3KWV6AB9Zaw9aa7OAj4CexphgIMha+6UtnPE5J9R/CfgT587/ZVUpc2mtXWGtzfPWXwuEVuYgqtDVQJq1dqe19hcgnsI5LaroHM8HunuPNG8F4q21R621u4A0b3tlabOmqfB5tNZmWms3Alhr/w/YBjQ+C2OpapXxN4kxJhToBbx1FsZwxlwIqEustZkA3p+NSinTGPi+yON077bG3vsnbscYEwdkWGs3VUanHVUpc3mC+yg8uqqJTjY3pZbxhvZhoP4p6palzZqmMubRx7uEFQWsq8A+u6qy5nIqhW/eCyq+yxXnrHxZrDFmJXBpKbtGl7WJUrbZk203xtT2th1bxvarjbM9lyc892ggD3injM9V3Zx2Dk5R5mTbS3sTWNOP6CtjHgsrGXMhhf9/8X9Za7PPuIfVR4XPpTHmt8A+a22iMabbr+xfpTorAWWt/c3J9hljfjTGBFtrM73LTPtKKZYOdCvyOBRY490eesL2vUBzCtdcN3nP84cCG40xV1trf/gVQ6lyVTCXx9seBPwW6O5dAqyJ0ik8b3lcsTk4oUy6McYD1AEOnqbu6dqsaSplHo0xARSG0zvW2n9WTtedUxlzGQfEGWNuAc4Dgowxc621d1fOEH6Fqj4JBkyh+In950spUw/YReFJ/bre+/W8+74COvOfE/u3lFJ/N+fGRRKVMpdAT2Ar0LCqx1jJ8+eh8KKRpvznhHT4CWUepfgJ6X9474dT/IT0TgpPcJ+2zZp2q6R5NBSeF51a1eOr7nN5Qt1uOHyRRNV3oHCtdBWQ6v15/MWyA/BWkXL3UXiSLw24t8j2DsBmCq9QeQXv1zed8BznSkBVylx6y30PJHtv06t6rJU4h7dQeIXYt8Bo77bxQJz3/nnAe945WQ80K1J3tLfedopfTVqizZp+q+h5BK6jcNkqpcjfYYk3ozXxVhl/k0X2Ox1Q+i4+ERFxkgtX8YmIiJSggBIREScpoERExEkKKBERcZICSkREnKSAEhERJymgRETESf8PBHOdrnNII60AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tau\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using advi...\n",
      "Average Loss = 286.87:  29%|██▉       | 14571/50000 [00:11<00:28, 1260.59it/s]\n",
      "Convergence achieved at 14600\n",
      "Interrupted at 14,599 [29%]: Average Loss = 530.41\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [value_sd, is_recovered|dataset_offset, is_recovered|dataset_sd, 1|dataset_offset, 1|dataset_sd, is_recovered, Intercept]\n",
      "Sampling 2 chains: 100%|██████████| 61000/61000 [04:20<00:00, 233.95draws/s]\n",
      "There were 337 divergences after tuning. Increase `target_accept` or reparameterize.\n",
      "There were 79 divergences after tuning. Increase `target_accept` or reparameterize.\n",
      "The number of effective samples is smaller than 10% for some parameters.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                                        mean        sd  hpd0.95_lower  \\\n",
      "is_recovered|dataset_sd             0.171427  0.191797       0.000004   \n",
      "is_recovered|dataset[folke2016]     0.014241  0.155317      -0.303411   \n",
      "is_recovered|dataset[krajbich2010]  0.000210  0.147585      -0.318546   \n",
      "is_recovered|dataset[krajbich2011]  0.015668  0.150932      -0.313711   \n",
      "is_recovered|dataset[tavares2017]  -0.030665  0.154272      -0.383646   \n",
      "1|dataset[folke2016]                0.753607  0.349399       0.070554   \n",
      "1|dataset[krajbich2010]            -0.329672  0.345228      -1.017057   \n",
      "1|dataset[krajbich2011]             0.026229  0.346175      -0.685019   \n",
      "1|dataset[tavares2017]             -0.413127  0.351035      -1.124929   \n",
      "is_recovered                        0.045551  0.149053      -0.248105   \n",
      "Intercept                           0.926238  0.340206       0.255271   \n",
      "value_sd                            0.743730  0.034783       0.676146   \n",
      "1|dataset_sd                        0.637339  0.249556       0.248613   \n",
      "\n",
      "                                    hpd0.95_upper  effective_n  gelman_rubin  \n",
      "is_recovered|dataset_sd                  0.534239         9787      1.000023  \n",
      "is_recovered|dataset[folke2016]          0.366882        17229      1.000024  \n",
      "is_recovered|dataset[krajbich2010]       0.316446        12566      1.000196  \n",
      "is_recovered|dataset[krajbich2011]       0.336622        12536      0.999998  \n",
      "is_recovered|dataset[tavares2017]        0.274286        10708      1.000165  \n",
      "1|dataset[folke2016]                     1.466746         1012      1.001111  \n",
      "1|dataset[krajbich2010]                  0.363144          963      1.001259  \n",
      "1|dataset[krajbich2011]                  0.709418         1001      1.001170  \n",
      "1|dataset[tavares2017]                   0.283241          940      1.001320  \n",
      "is_recovered                             0.345949        13224      1.000063  \n",
      "Intercept                                1.626750          918      1.001313  \n",
      "value_sd                                 0.812373        31846      0.999993  \n",
      "1|dataset_sd                             1.143974         1456      1.000528  \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd8VFXex/HPCQkBRKoGiMhSFUgIQ0dBqlJEEZSVACqwdlFXeBRQVATRwCKorFmxPouForACiiBSIsojIiXSlB6aIDWAtLTz/DHDbEJ6IOQEvu/Xa16Zufece3/3MMx3bpkZY61FRETENQEFXYCIiEhGFFAiIuIkBZSIiDhJASUiIk5SQImIiJMUUCIi4iQFlIiIOEkBJYWaMeb/CrqGi8kY088Y81YW84ONMdOMMVuMMT8ZY6pevOpELiwFlBRq1tobz6e/MabIharFkRruB45Ya2sCrwNjLuCyRS4qBZQUasaYP31/KxljlhhjYo0x64wxN2XVxxgz0hjzE3CDMaaRMeY7Y8xKY8w3xphKvnY1jTELjDG/GGNWGWNqGK+xvnWsNcb09LWdZoy5NdU6/m2MucsYU8TX/mdjzBpjzMO++W2MMYuNMZOBtb5p9xhjlvu24Z2zwWWM6W+M2WSM+Q5okc2Q3AFM8t2fDrQ3xpg8DK1IgVNAyaWiN/CNtdYD1Adis2h7BbDOWtsM+An4J9DDWtsI+BB4xdfuUyDaWlsfuBHYC9wJnF3HzcBYX6BNBc6GVVGgPfA13j2ao9baJkAT4EFjTDXf8psCw6y1dY0xdXz9W/i2IRno41v2CLzBdAtQN5txuAbYBWCtTQKOAuWz6SPipMCCLkDkAvkZ+NAYEwTMtNZmFVDJwAzf/euBcOBb345GEWCvMeZK4Bpr7RcA1trTAMaYlsAUa20y8Idvr6YJMBeYYIwJBjoBS6y1p4wxHYAIY0wP3/pKA7WABGC5tXa7b3p7oBHws6+O4sB+oBkQY6094Fv/NOC6LLYto70lfeGmFEoKKLkkWGuXGGNaAV2Aj40xY621H2XS/LQvYMD7gr7eWntD6gbGmFKZ9M3wcJm19rQxJgboiHdPaEqq9k9Ya785Z/ltgBPnLHeStfbZc9p1I3cBsxu4FthtjAnEG4iHc9FfxBk6xCeXBGPMX4D91tr3gA+AhjnsuhG42hhzg285QcaYMGvtMbwv8t1804ONMSWAJUBP37mlq4FWwHLfsqYC/YGbgLOB9A3wqG/PDmPMdcaYKzKoYyHQwxgT4mtXzrdNPwFtjDHlfcv4azbbMxvo67vfA1hk9ZMFUkhpD0ouFW2AZ4wxicCfwH056WStTfAdfptgjCmN9//EG8B64F7gHWPMSCARbzh8AdwA/IJ3z2awtXafb3HzgY+A2dbaBN+094GqwCrfxQoHgG4Z1LHBGPM8MN8YE+Bb3wBr7TJjzEvAj3jPga3CexgyMx/g3YPcgnfPKTIn4yDiIqM3VyIi4iId4hMRESfpEJ9csnyfcwo+Z/K91tq1BVHPhWSMGUb681GfW2tfyai9SGGkQ3wiIuKk3OxBKclERORCyfYbTnQOSkREnKSAEhERJymgRETESQooERFxkgJKREScpIASEREnKaBERMRJCigREXGSAkpERJykgBIREScpoERExEkKKJEsJKXk7isoc9teRDKXm28z1/88uSyNXn0wx22HNrgqHysRuaToy2JFRKRwUkCJiIiTFFAiIuIkBZSIiDhJASUiIk5SQImIiJMUUCIi4iQFlIiIOEkBJZcVfdODSOERWNAFiFxMgQFG3wwhUkhoD0pERJykgBIREScpoERExEkKKBERcZICSkREnKSAEhERJymgRETESQooERFxkgJKREScpIASEREnKaBERMRJCigREXGSAkpERJykgBIREScpoERExEkKKBERcZICSi4ZVatW5eDBnP8YYUaSEs4wecgDjO3ahOj7OhIXF5dhu41LFzKue3PGdm1CzP++6Z/er18/qlWrhsfjwePxEBsbC0BMTAylS5f2Tx85cmSa5SUnJ9OgQQNuu+02/7Tt27fTrFkzatWqRc+ePUlISDivbRMpbBRQIqn8PPNTipcqwzOzf6Zln0cYMmRIujYpycnMHjOU/v+cysAZS/ll3hf8sW2jf/7YsWOJjY0lNjYWj8fjn37TTTf5p7/44otplvnmm29Sp06dNNOGDBnCwIED2bx5M2XLluWDDz64wFsr4jYFlBSouLg4ateuzQMPPEB4eDh9+vRhwYIFtGjRglq1arF8+XIOHz5Mt27diIiIoHnz5qxZswaAQ4cO0aFDBxo0aMDDDz+Mtda/3E8++YSmTZvi8Xh4+OGHSU5OzlE9v8bMpeFtPQEIb387CxcuTLNcgF3rVlG+clXKVa5KYFBR6nfsxq8xc/M8Brt372bOnDk88MAD/mnWWhYtWkSPHj0A6Nu3LzNnzszzOkQKIwWUFLgtW7bw97//nTVr1vDbb78xefJkfvjhB1577TVeffVVhg8fToMGDVizZg2vvvoq9913HwAjRoygZcuWrF69mq5du7Jz504Afv31V6ZNm8bSpUuJjY2lSJEifPrppwD07NmTCZFt0t1WfTUNgGMH9lGm4jUAFAkMpHTp0pyMP5ym3mMH9lLa1wagVEgoR/fv9T8eNmwYERERDBw4kDNnzvin//jjj9SvX5/OnTuzfv16//SnnnqKf/zjHwQE/Pe/46FDhyhTpgyBgYEAVK5cmT179pz/YIsUIoEFXYBItWrVqFevHgBhYWG0b98eYwz16tUjLi6OHTt2MGPGDADatWvHoUOHOHr0KEuWLOE///kPAF26dKFs2bIALFy4kJUrV9KkSRMATp06RUhICADTpk1j9OrMz1Odu7cEgDHnNsqgibdNVFQUFStWJCEhgYceeogxY8bw4osv0rBhQ3bs2EHJkiX5+uuv6datG5s3b+arr74iJCSERo0aERMTk2Ud5tw6RC5xCigpcMHBwf77AQEB/scBAQEkJSX59yJSO/tindGLtrWWvn37EhUVlW5ez549+T52fbrpLe95lIa39aR0SCXi9+2hdIVQkpOSOHr0KCVKl03TtlRIKEf3/Xdv5tj+3yl1dUUAKlWq5N+m/v3789prr3n7lCrlb3/rrbfy2GOPcfDgQZYuXcrs2bP5+uuvOX36NMeOHeOee+7h448/Jj4+3r/9u3fvJjQ0NJMRFLk06RCfOK9Vq1b+Q3QxMTFcddVVlCpVKs30uXPncuTIEQDat2/P9OnT2b9/PwCHDx9mx44dgHcP6smpMeluZ8871WndyX+4b93CL2nXrl26EKwc1oCDu7ZzeM8OkhIT+OWbmdRp3QmAvXu9h/qstcycOZPw8HAA9u3b598rWr58OSkpKZQvX56oqCh2795NXFwcU6dOpV27dnzyyScYY2jbti3Tp08HYNKkSdxxxx35MLoi7tIelDjvpZdeon///kRERFCiRAkmTZoEwPDhw+nVqxcNGzakdevWVKlSBYC6desyatQoOnToQEpKCkFBQURHR/OXv/wl23U17taHz154jLFdm1CidFm+mz2dz456z03NGPkU/f85lSKBgXQdEsWHA+7GpqTQuGsvKtSoDUCfPn04cOAA1lo8Hg8TJ04EYPr06bz99tsEBgZSvHhxpk6dmu0huzFjxhAZGcnzzz9PgwYNuP/++89nGEUKHZPhMfeM5bihiMuyOgd1rqENrsp1exHJkWxPquoQn4iIOEkBJSIiTlJAiYiIkxRQIiLiJAWUiIg4SQElIiJOUkCJiIiTFFAiIuIkfZOEXBZ6z+jtv7/hyBkiG72TL+tJSrEEBuTuS13z0kfkcqCAErmAAgNMrr55AvTtEyKZ0SE+ERFxkgJKREScpIASEREnKaBERMRJCigREXGSAkpERJykgBIREScpoERExEkKKBERcZICSkREnKSAEhERJymgRETESQooERFxkgJKREScpIASEREnKaBERMRJCigREXGSAkpERJykgBIREScpoERExEkKKBERcZICSkREnKSAEhERJymgRETESQooERFxkgJKREScpIASEREnKaBERMRJCigREXGSAkryxenTp2natCn169cnLCyM4cOH++f169ePatWq4fF48Hg8xMbGpuu/ePFi/3yPx0OxYsWYOXMmAH369CEiIoLnnnvO3/7ll19m1qxZea5349KFjOvenLFdmxDzv2/meTkxMTFMiGzD6z1a8u4DXdPMS0lOZkKvtvz7yd4Z9p04cSL16tXD4/HQsmVLNmzYAMDSpUuJiIigSZMmbNmyBYD4+Hg6duyItTbPtYq4LrCgC5BLU3BwMIsWLaJkyZIkJibSsmVLOnfuTPPmzQEYO3YsPXr0yLR/27Zt/cF1+PBhatasSYcOHVizZg0Aa9as4aabbuLQkXgSTp9i+fLlvPDCC7mq8dSxeIqXKkNKcjKzxwzl/n99TqkKoUTf04E6rTtRofr12fZNM+34UR577DHue30yZSpV5s/DB9LMXzrlXUKqXcfpP49nuMzevXvzyCOPADB79mwGDRrEvHnzGDduHDNmzCAuLo63336bcePG8fLLL/Pcc89hjMnVNosUJgooyRfGGEqWLAlAYmIiiYmJeX4xnT59Op07d6ZEiRIEBQVx6tQpUlJSSEhIIDgokM73D6Z574GMXn0w02VsOHIGgDplinL81/1MnfEwu9at4pnZP7Nr3SrKV65KucpVAajfsRu/xsxNF1B/Hj7Aqq8+Y9WXU7mh5/0069EvzfzYuTO48847KVmpMgAly13tn3f0j9/Z+P23tL1/IN9/8naGNZYqVcp//8SJE/7xOrvNJ0+eJCgoiK1bt7Jnzx5at26dswEUKaQUUJJvkpOTadSoEVu2bGHAgAE0a9bMP2/YsGGMHDmS9u3bM3r0aIKDgzNdztSpUxk0aBAAderUoUqVKjRs2JB7772XLVu2YK0ltHZElrUkxp/i8NIdbP95N6ZCSW7q+yh3v/wvAI4d2Evpitf425YKCWXXupUApKSksHHpQlbM/IT92zZRv/Od9H9rGqUrhKZbx8EdWzlSNpDJD97BmRN/0qL3QzS8rScAX702jM5/H86Zk39mWWd0dDTjx48nISGBRYsWAfDss8/y0EMPUbx4cT7++GOefvppXn755SyXI3IpUEBJvilSpAixsbHEx8fTvXt31q1bR3h4OFFRUVSsWJGEhAQeeughxowZw4svvpjhMvbu3cvatWvp2LGjf9obb7zhv3/77bdzyxNRLH5/PHs3r6dmszY0vfPeNMvYtW4VG4bMo3zLqtz88s1sTYbajTr8t0EG53HO7r1069aN735awZ0vvM51N7bLci8wJTmJlSt/od+EaSSePs3b/Tpzbb3GHNyxlSvKXc01deuzbcXSLMdswIABDBgwgMmTJzNq1CgmTZqEx+Nh2bJlACxZsoTQ0FCstfTs2ZOgoCDGjRtHhQoVslyuSGGkiyQk35UpU4Y2bdowb948ACpVqoQxhuDgYPr378/y5csz7fvZZ5/RvXt3goKC0s2bNWsWjRs3JuHUCfZt/Y3eYz5g9ZzPSDh1Mk27irXqcm3fhpzee4wlY5ZwaMn2NOeBSoWEcnTfHv/jY/t/p9TVFQGIiooivP3tfPmPZ5k1egi71q/OtNbSFULp1KkTRYtfwRVly1Ot4Q3s27SOHb/8xK/fzWNMl4ZMefZBtq34gWnDHs1yzCIjI/0XhZxlrWXUqFG88MILjBgxghEjRnDPPfcwYcKELJclUlgpoCRfHDhwgPj4eABOnTrFggULqF27NuDdKwLvC+7MmTMJDw/PdDlTpkyhV69e6aYnJiby5ptv8swzz5B4+pR/z8baFJKTEtO0DQouRrkb/0LNwa1p/nhzzhw4wT97t2Pa848BUDmsAQd3befwnh0kJSbwyzczqdO6EwBhYWHc/swrPDX9B6o3vIH50a/y5t2t2fTj4nQ11W3dme+//57kpCQSTp1k17pVXF3tOjo98QLPzlvDkDmr6BX1HtUbt6TnK+nPQ23evNl/f86cOdSqVSvN/EmTJtGlSxfKli3LyZMnCQgIICAggJMnT567KJFLgg7xSb7Yu3cvffv2JTk5mZSUFO6++25uu+02wHuZ+IEDB7DW4vF4mDhxIgArVqxg4sSJvP/++wDExcWxa9euDC8GiI6Opm/fvpQoUYKKtcKw1vLG3a24vsXNFL+ydKZ1XVnpSkLvCuevI6PZ+MO3ABQJDKTrkCg+HHA3NiWFxl17UaFG7TT9AoOKEtGxOxEdu3Pk912cjD+cbtkh1a+jU6dOvNazNSYggMbd+lCxZp0sx+nbt0dT97ZWdO3albfeeosFCxYQFBRE2bJlmTRpkr/dyZMnmTRpEvPnzwdg0KBB3HXXXRQtWpQpU6ZkuQ6Rwsrk4nMU+sCFOCmrq/fOmrryYQDqlg1mw5EzRDZ6J0fLHtrgqhwtP6/tz/YRuQxle1mvDvGJiIiTFFAiIuIkBZSIiDhJASUiIk5SQImIiJMUUCIi4iQFlIiIOEkBJSIiTlJAiYiIkxRQIiLiJAWUiIg4SQElIiJOUkCJiIiTFFAiIuIkBZSIiDhJASUiIk5SQImIiJMUUCIi4iQFlIiIOEkBJSIiTlJAiYiIkxRQIiLiJAWUiIg4SQElIiJOUkCJiIiTFFAiIuIkBZSIiDhJASUiIk5SQImIiJMUUCIi4iQFlIiIOEkBJSIiTlJAiYiIkxRQIiLiJAWUiIg4SQElIiJOUkCJiIiTFFAiIuIkBZSIiDhJASUiIk5SQImIiJMUUCIFLCnF5mt7kcIqsKALELncBQYYRq8+mOP2QxtclY/ViLhDe1AiIuIkBZSIiDhJASUiIk5SQImIiJMUUCIi4iQFlIiIOEkBJSIiTlJAiYiIkxRQIiLiJAWUiIg4SQElIiJOUkCJiIiTFFAiIuIkBZSIiDhJASUiIk5SQImIiJMUUCIi4iQFlIiIOEkBJSIiTlJAiTOSUmxBlyAiDgks6AJEzgoMMIxefTBXfYY2uCqfqhGRgqY9KBERcZICSkREnKSAEhERJymgRETESQooERFxkgJKREScpIASEREnKaBERMRJCigREXGSAkpERJykgBIREScpoERExEkKKBERcZICSkREnKSAEhERJymgRETESQooERFxkgJKREScpIASEREnKaBERMRJCigREXGSAkoKpaWT3+GNv95EWFgYP3w60T99wcR/ENWxHhMi2zAhsg2//fAtACe2HGLjSwv4Zsg3nNn/JwCnjh/lw8f+irU2w3W8++AdrFixwv/4yO87eeOvNwGwbcVSXmpVnQm92jL+zhtZ8M5YAGJiYvzTx3Vvzjv3386vS+bnyxiIXOoCC7oAkdzat+VXfv7iEx776BueaxpK7RbtqX3TLVxVpQYALfo8Qqv7BqTpc2D+Zqo+2pyKZxJYH7MNOsOi98bR5v6nMMbkqY6qnub0mzCZhFMnmBDZljo3dYBrivinA/y+cS0fD+pLUHAxajZrdX4bLnKZ0R6UFDoHtm/i2nqNKFq8BIGBgVRrdCPrF32ddacihpTEZJLPJEORAA7t2s6x/Xup3qjFeddTtPgVXFOnPod2b083L/T6erR/8H/48bMPzns9Ipcb7UFJoVOhRh2+iX6VE/GHOXmyBBt/WEDluvX983+c9gGrv/qMa+rWp8ugkRQvVYYKna9n90eriC8RxFX3NWR+9Kvc8tjQbNfVp08fjhMEQHJiIiYg/Xu6E/GH2bl2Je0e/B8gKd380DoRLPkoOu8bLHKZUkBJoRNS/Tpa93uCDx/rwXdXl6HSdWEEFPE+lZv9tZ83KIzh239FMWf8i/R4aQLFq5Sh1nNtqVs2mOU/7aHc1VWxFiYPeYAigUHcOmgEV5YPSbeuTz/9lAVFqgLec1CT/t7HPy8udhkTerXFmADa9H+SCjVqw9F16QvO5ByXiGRNh/ikUGrS7R6emLyIJUuWUKJUGcpXqQ7AleVDCChShICAAJreeS+7169O089ayx9zfqP9g//DwnfHcvMjQ/Dc2oP/m/Jermuo6mnOk1MW88TkhTTr0S/Tdr//tpaQatflevkilzsFlBRKfx4+AMDOnTtZv3gOnk53AnDswD5/m/WLvvbu1aSyPWY7pepVpHipMiSePoUJCMAEBJB4+lS+1Ll303oWvT+e5nf/LV+WL3Ip0yE+KZQ+fbo/J48eYcqVxeg6ZAzFS5UBYO6bI9m7aR0GQ9nQa+k27DV/n5QzSWyP2U7Fx28EoGWfR/j0mf4UCQwiMurdC1bb2UN/iadPUbLsVdz+zKu6gk8kD0xmnwHJgA6kS74bvfpgrtoPbXBVjvpMXfkwAHXLBrPhyBkiG71zQZef1/Z5XYfIJSDbz3foEJ+IiDhJASUiIk5SQIkUMkkpuTvantv2Iq7QRRIihUxggNE5K7ksaA9KREScpIASEREnKaBERMRJCigREXGSAkpERJykgBIREScpoERExEkKKMk3+oCoiJwPfVBX8o0+UCoi50N7UCIi4iQFlIiIOEkBJSIiTlJAiYiIkxRQIiLiJAWUiIg4SQElIiJOUkCJiIiTFFAiIuIkBZRjrLU8+eST1KxZk4iICFatWpWuzcmTJ+nSpQu1a9cmLCyMoUOHpmszffp0jDGsWLECgMTERPr27Uu9evWoU6cOUVFR/rZ/+9vfCAkJITw8PM0yPv/8c8LCwggICPAvRy4N8+bN4/rrr6dmzZqMHj063fyJEydSr149PB4PLVu2ZMOGDWnm79y5k5IlS/Laa6/5p1WtWtXfp3Hjxv7pv/zyCzfccAP16tXj9ttv59ixYwDExcVRvHhxPB4PHo+HRx55JJ+2VgorBZRj5s6dy+bNm9m8eTPvvvsujz76aIbtnn76aX777TdWr17N0qVLmTt3rn/e8ePHmTBhAs2aNfNP+/zzzzlz5gxr165l5cqVvPPOO8TFxQHQr18/5s2bl24d4eHh/Oc//6FVq1YXdiOlQCUnJzNgwADmzp3Lhg0bmDJlSroA6t27N2vXriU2NpbBgwczaNCgNPMHDhxI586d0y178eLFxMbGpnlD88ADDzB69GjWrl1L9+7dGTt2rH9ejRo1iI2NJTY2lokTJ17gLZXCTgHlmFmzZnHfffdhjKF58+bEx8ezd+/eNG1KlChB27ZtAShatCgNGzZk9+7d/vkvvPACgwcPplixYv5pxhhOnDhBUlISp06domjRopQqVQqAVq1aUa5cuXS11KlTh+uvvz4/NlMK0PLly6lZsybVq1enaNGiREZGMmvWrDRtzj43AE6cOIExxv945syZVK9enbCwsBytb+PGjf43ObfccgszZsy4AFshlwMFlGP27NnDtdde639cuXJl9uzZk2n7+Ph4vvzyS9q3bw/A6tWr2bVrF7fddluadj169OCKK66gUqVKVKlShaeffjrDUJJLX06fY9HR0dSoUYPBgwczYcIEwBtWY8aMYfjw4enaG2Po0KEDjRo14t133/VPDw8PZ/bs2YB3T37Xrl3+edu3b6dBgwa0bt2a77///oJto1waFFCOsTb9T1SkfveaWlJSEr169eLJJ5+kevXqpKSkMHDgQMaNG5eu7fLlyylSpAi///4727dvZ9y4cWzbtu2C1y/uOfdnT7J7jp1tP2DAALZu3cqYMWMYNWoUAMOHD2fgwIGULFky3TKWLl3KqlWrmDt3LtHR0SxZsgSADz/8kOjoaBo1asTx48cpWrQoAJUqVWLnzp2sXr2a8ePH07t3b//5KRHQz204ITo6mvfeew+AJk2apHmHuXv3bkJDQzPs99BDD1GrVi2eeuopwHvuad26dbRp0waAffv20bVrV2bPns3kyZPp1KkTQUFBhISE0KJFC1asWEH16tXzd+OkwJ37syc7TpRkyfqt/mkxyzcCpf2Pz/3Zk8jISP+50J9++onp06czePBg4uPjCQgIoFixYjz++OP+52lISAjdu3dn+fLltGrVitq1azN//nwANm3axJw5cwAIDg4mODgYgEaNGlGjRg02bdqU5gILubxpD8oBAwYM8J8o7tatGx999BHWWpYtW0bp0qWpVKlSuj7PP/88R48e5Y033vBPK126NAcPHiQuLo64uDiaN2/O7Nmzady4MVWqVGHRokVYazlx4gTLli2jdu3aF3MzxRGVwxpwcNd2Du/ZQVJiAr98M5M6rTulabN582b//Tlz5lCrVi0Avv/+e//z66mnnuK5557j8ccf58SJExw/fhzwHgacP3++/6rQ/fv3A5CSksKoUaP8V+sdOHCA5ORkALZt28bmzZv1hknSUEA55tZbb6V69erUrFmTBx98kH/961/+eR6PB/DuVb3yyits2LCBhg0b4vF4eP/997Nc7oABA/jzzz8JDw+nSZMm9O/fn4iICAB69erFDTfcwMaNG6lcuTIffPABAF988QWVK1fmxx9/pEuXLnTs2DGftloupiKBgXQdEsWHA+7m9btaEHFLVyrUqM23b49mw3feqznfeustwsLC8Hg8jB8/nkmTJmW5zD/++IOWLVtSv359mjZtSpcuXejUyRt6U6ZM4brrrqN27dqEhobSv39/AJYsWUJERAT169enR48eTJw4UedFJQ2T0fHoTOj3uy9zSSmWwICMz4dlJre/qJub9rnpM3XlwwDULRvMhiNniGz0Tr7UlJ/bcDHbi1wE2b6Y6ByU5Jh+wl1ELiYd4hMREScpoERExEkKKBERcZICSkTSOPeDvfnVRyQ7ukhCRNLI7cUwoAtiJH9oD+oypne9IuIy7UFdxnTZuIi4TB/UvczlJqCebXh1PlaSz+465/Fl/osPUasOZDovrx82FsmlbD+oq0N8InLecnu4WIeXJSd0iO8SkZevIRK5UHS4WPKDAuoSoSuvRORSo4BylIt7RNZap77U9LzWMSyfl58Lrn1Z7MWQ2+d3YoolKBftXfz/I7mngMqj/P4PpkMmcinLy/M7N+2frl8+1zUp1NyT46v4RowYMQ9w7VUwFPi9oIvII9VeMApr7YW1blDtBcX12g8OHz68U5YtrLWF9vbSSy/Zgq5BtReuW2GtvbBuAW+RAAAF6ElEQVTWrdpV+/ncdJm5iIg4qbAH1IiCLuA8qPaCUVhrL6x1g2ovKIW5diB33yQhIiJy0RT2PSgREblEKaBERMRJCigREXGS8wFljClnjPnWGLPZ97dsBm3aGmNiU91OG2O6+eb92xizPdU8j0u1+9olp6pvdqrp1YwxP/n6TzPGFHWpdmOMxxjzozFmvTFmjTGmZ6p5F3XcjTGdjDEbjTFbjDFDM5gf7BvDLb4xrZpq3rO+6RuNMR3zs86M5KD2QcaYDb4xXmiM+UuqeRk+dy6WHNTezxhzIFWND6Sa19f3/NpsjOnrWN2vp6p5kzEmPtW8gh7zD40x+40x6zKZb4wxE3zbtsYY0zDVvAIb8zwp6Ovcs7sB/wCG+u4PBcZk074ccBgo4Xv8b6CHy7UDf2Yy/TMg0nd/IvCoS7UD1wG1fPdDgb1AmYs97kARYCtQHSgK/ALUPafNY8BE3/1IYJrvfl1f+2Cgmm85RS7iOOek9rapns+Pnq09q+eOQ7X3A97KoG85YJvvb1nf/bKu1H1O+yeAD10Yc9/6WwENgXWZzL8VmIv35yyaAz8V9Jjn9eb8HhRwBzDJd38S0C2b9j2Audbak/laVc7ktnY/Y4wB2gHT89L/Asi2dmvtJmvtZt/934H9QEH8aFRTYIu1dpu1NgGYirf+1FJvz3SgvW+M7wCmWmvPWGu3A1t8y7tYsq3dWrs41fN5GVD5ItaXlZyMe2Y6At9aaw9ba48A3wJZf6vAhZPbunsBUy5KZTlgrV2C9014Zu4APrJey4AyxphKFOyY50lhCKgK1tq9AL6/Idm0jyT9k+kV367u68aY4PwoMhM5rb2YMWaFMWbZ2UOTQHkg3lqb5Hu8G7gmf8tNI1fjboxpivfd6NZUky/WuF8D7Er1OKOx8rfxjelRvGOck775Kbfrvx/vu+OzMnruXCw5rf0u3/NgujHm2lz2zQ85XrfvcGo1YFGqyQU55jmR2fYV9HM915z4slhjzAKgYgazcvid0/7lVALqAd+kmvwssA/vi+e7wBBgZN4qzXCdF6L2Ktba340x1YFFxpi1wLEM2l3QD61d4HH/GOhrrU3xTc7XcT+3hAymnTtWmbXJSd/8lOP1G2PuARoDrVNNTvfcsdZuzah/PshJ7V8CU6y1Z4wxj+Ddi22Xw775JTfrjgSmW2uTU00ryDHPCVef67nmREBZa2/ObJ4x5g9jTCVr7V7fC+H+LBZ1N/CFtTYx1bL3+u6eMcb8L/D0BSn6v8s/79p9h8ew1m4zxsQADfD+KHkZY0yg7x1/ZS7wFz9eiNqNMaWAOcDzvsMJZ5edr+N+jt3AtakeZzRWZ9vsNsYEAqXxHibJSd/8lKP1G2NuxvvGobW19szZ6Zk8dy7Wi2W2tVtrD6V6+B4wJlXfNuf0jbngFWYsN//mkcCA1BMKeMxzIrPtK8gxz5PCcIhvNnD2apO+wKws2qY7Vux7cT17TqcbkOGVL/kk29qNMWXPHv4yxlwFtAA2WO9ZzcV4z6ll2j8f5aT2osAXeI93f37OvIs57j8DtYz3qseieF9Uzr26KvX29AAW+cZ4NhBpvFf5VQNqAcvzsdZzZVu7MaYB8A7Q1Vq7P9X0DJ87F63ynNVeKdXDrsCvvvvfAB1821AW6EDaIx/5KSfPF4wx1+O9mODHVNMKesxzYjZwn+9qvubAUd8bxoIc87wp6Ks0srvhPU+wENjs+1vON70x8H6qdlWBPUDAOf0XAWvxvkB+ApR0qXbgRl99v/j+3p+qf3W8L5ZbgM+BYMdqvwdIBGJT3TwFMe54r1zahPed7DDftJF4X9QBivnGcItvTKun6jvM128j0LkAnuPZ1b4A+CPVGM/O7rnjUO1RwHpfjYuB2qn6/s3377EF6O9S3b7HLwGjz+nnwphPwXvFbCLevaL7gUeAR3zzDRDt27a1QGMXxjwvN30Xn4iIOKkwHOITEZHLkAJKREScpIASEREnKaBERMRJCigREXGSAkpERJykgBIRESf9P7GTtA1of/TLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "bias_analysis_results = {}\n",
    "models = {}\n",
    "n_samples = 30000\n",
    "parameters = ['v', 'gamma', 's', 'tau']\n",
    "for parameter in parameters:\n",
    "    print(parameter)\n",
    "    result_path = 'results/mixed_effects_models/si6_parameter_recovery_{}_bias_analysis.csv'.format(parameter)\n",
    "    plot_path = 'results/mixed_effects_models/si6_parameter_recovery_{}_bias_analysis.png'.format(parameter)\n",
    "    if (not os.path.isfile(result_path)) or (not os.path.isfile(plot_path)):\n",
    "        model = Model(df.loc[df['parameter'] == parameter].reset_index(drop=True).copy())\n",
    "        models[parameter] = model\n",
    "        result_fit = model.fit('value ~ is_recovered',\n",
    "                           random=['is_recovered|dataset'],\n",
    "                           categorical=['dataset'],\n",
    "                           samples=n_samples,\n",
    "                           chains=2)\n",
    "        # save summary\n",
    "        result = result_fit.summary(ranefs=True)\n",
    "        print(result)\n",
    "        result.to_csv(result_path)\n",
    "        # save posteriorplot\n",
    "        plot_posterior(model.backend.trace,\n",
    "                       varnames=['is_recovered'],\n",
    "                       round_to=5,\n",
    "                       point_estimate='mode',\n",
    "                       ref_val=0)\n",
    "        plt.savefig(plot_path, dpi=330)\n",
    "        plt.show()\n",
    "    else:\n",
    "        print('Results loaded from: {}'.format(result_path))\n",
    "        result = pd.read_csv(result_path)\n",
    "        print(result)\n",
    "        show_img(plot_path)\n",
    "        \n",
    "    # bias_analysis_results[parameter] = result "
   ]
  }
 ],
 "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.5"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "30px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false,
   "widenNotebook": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}