10 Chapter 10/02 Graduate school admissions.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import pandas as pd\n",
    "import pymc3 as pm\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats\n",
    "import numpy as np\n",
    "\n",
    "plt.style.use('seaborn-whitegrid')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def summarize(alpha=0.05):\n",
    "\n",
    "    def _summarize(ser):\n",
    "        hpd_ = pm.hpd(ser, alpha=alpha)\n",
    "        return pd.Series({\n",
    "            'mean_': ser.mean(),\n",
    "            'lower': hpd_[0],\n",
    "            'upper': hpd_[1],\n",
    "        })\n",
    "\n",
    "    return _summarize\n",
    "\n",
    "\n",
    "def agg_prediction(trace_var, source_data, group_columns, aggregate_fc):\n",
    "    df = pd.DataFrame(trace_var).stack().reset_index()\n",
    "    merged = pd.merge(\n",
    "        df, source_data,\n",
    "        left_on='level_1', right_index=True\n",
    "    )\n",
    "\n",
    "    return (\n",
    "        merged.groupby(group_columns)[0].apply(aggregate_fc)\n",
    "        .unstack().sort_index().reset_index()\n",
    "    )\n"
   ]
  },
  {
   "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>dept</th>\n",
       "      <th>applicant.gender</th>\n",
       "      <th>admit</th>\n",
       "      <th>reject</th>\n",
       "      <th>applications</th>\n",
       "      <th>dept_i</th>\n",
       "      <th>is_male</th>\n",
       "      <th>ratio</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>A</td>\n",
       "      <td>male</td>\n",
       "      <td>512</td>\n",
       "      <td>313</td>\n",
       "      <td>825</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0.620606</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>A</td>\n",
       "      <td>female</td>\n",
       "      <td>89</td>\n",
       "      <td>19</td>\n",
       "      <td>108</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.824074</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>B</td>\n",
       "      <td>male</td>\n",
       "      <td>353</td>\n",
       "      <td>207</td>\n",
       "      <td>560</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0.630357</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>B</td>\n",
       "      <td>female</td>\n",
       "      <td>17</td>\n",
       "      <td>8</td>\n",
       "      <td>25</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0.680000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>C</td>\n",
       "      <td>male</td>\n",
       "      <td>120</td>\n",
       "      <td>205</td>\n",
       "      <td>325</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>0.369231</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>C</td>\n",
       "      <td>female</td>\n",
       "      <td>202</td>\n",
       "      <td>391</td>\n",
       "      <td>593</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>0.340641</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>D</td>\n",
       "      <td>male</td>\n",
       "      <td>138</td>\n",
       "      <td>279</td>\n",
       "      <td>417</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>0.330935</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>D</td>\n",
       "      <td>female</td>\n",
       "      <td>131</td>\n",
       "      <td>244</td>\n",
       "      <td>375</td>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>0.349333</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>E</td>\n",
       "      <td>male</td>\n",
       "      <td>53</td>\n",
       "      <td>138</td>\n",
       "      <td>191</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>0.277487</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>E</td>\n",
       "      <td>female</td>\n",
       "      <td>94</td>\n",
       "      <td>299</td>\n",
       "      <td>393</td>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "      <td>0.239186</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>F</td>\n",
       "      <td>male</td>\n",
       "      <td>22</td>\n",
       "      <td>351</td>\n",
       "      <td>373</td>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "      <td>0.058981</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>F</td>\n",
       "      <td>female</td>\n",
       "      <td>24</td>\n",
       "      <td>317</td>\n",
       "      <td>341</td>\n",
       "      <td>5</td>\n",
       "      <td>0</td>\n",
       "      <td>0.070381</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   dept applicant.gender  admit  reject  applications  dept_i  is_male  \\\n",
       "0     A             male    512     313           825       0        1   \n",
       "1     A           female     89      19           108       0        0   \n",
       "2     B             male    353     207           560       1        1   \n",
       "3     B           female     17       8            25       1        0   \n",
       "4     C             male    120     205           325       2        1   \n",
       "5     C           female    202     391           593       2        0   \n",
       "6     D             male    138     279           417       3        1   \n",
       "7     D           female    131     244           375       3        0   \n",
       "8     E             male     53     138           191       4        1   \n",
       "9     E           female     94     299           393       4        0   \n",
       "10    F             male     22     351           373       5        1   \n",
       "11    F           female     24     317           341       5        0   \n",
       "\n",
       "       ratio  \n",
       "0   0.620606  \n",
       "1   0.824074  \n",
       "2   0.630357  \n",
       "3   0.680000  \n",
       "4   0.369231  \n",
       "5   0.340641  \n",
       "6   0.330935  \n",
       "7   0.349333  \n",
       "8   0.277487  \n",
       "9   0.239186  \n",
       "10  0.058981  \n",
       "11  0.070381  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = pd.read_csv('../data/UCBadmit.csv', sep=';').reset_index(drop=True)\n",
    "\n",
    "dept_idx = (\n",
    "    data.dept.drop_duplicates().sort_values().reset_index(drop=True).reset_index()\n",
    "    .set_index('dept')['index'].rename('dept_i').to_frame()\n",
    ")\n",
    "\n",
    "data = pd.merge(\n",
    "    data,\n",
    "    dept_idx,\n",
    "    left_on='dept',\n",
    "    right_index=True)\n",
    "\n",
    "data = data.assign(\n",
    "    is_male=(data['applicant.gender'] == 'male').astype('int'),\n",
    "    ratio=data.admit / data.applications)\n",
    "\n",
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**single intercept**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [alpha]\n",
      "Sampling 2 chains: 100%|██████████| 8000/8000 [00:04<00:00, 1827.20draws/s]\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as model_single_intercept:\n",
    "    alpha = pm.Normal('alpha', mu=0, sd=1.5)\n",
    "    \n",
    "    p = pm.Deterministic('p', pm.math.invlogit(alpha))\n",
    "    \n",
    "    admit = pm.Binomial('admit', n=data.applications.values, p=p,\n",
    "                        observed=data.admit.values)\n",
    "    \n",
    "    model_single_intercept.name = 'single intercept'\n",
    "    model_single_intercept.trace = pm.sample(2000, tune=2000, chains=2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAEcCAYAAACBPmBXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAGYZJREFUeJzt3XuQldWZ7/FvazgIwTGUmjCTgzI4utJTzhlmOEwZLwFHJMiMMkFLPBwEvCDeYoi3oPESjaNmvASJShCFVjgG7yIJegQvhAgpOJYkmuAS5WLLRZEA0m5ElD5/7G3sYaB9UejVe/H9VFHdvffu9/31qqYfnvW8/VLT2NiIJEktbY/UASRJuycLkCQpCQuQJCkJC5AkKQkLkCQpCQuQJCmJL6UOIO2oEEJbYCLQFXgPOC/GuCiE8B3gZqC+8tKrgReBJ4B2wIgY4+9DCEcCR8QYf7Kd458FDAa2AG2AH8YYn/8cOc8GOgE/B66KMZ4bQlgKfCPG+EGT1w2rPDaqwDHrgCkxxqe28/wBwN/HGKftaN7POG8v4OwY4yk787javdkBqRoNBxpijIcB3wVurzzeHbg0xtir8mcW0IdyAToXOCOEUAN8D7htWwcOIZwCHAscE2PsRbkQTQoh7Pd5w8YYV8UYz/28n7+D/hk4ooXOJX0hdkCqRn8LPAkQY4whhNrK492BfwghjATmAT8AGih3P+2A94FBwGNNO5CtjAAujDFurhx/SQihW4xxTaX72Lfy51+AS4GjgD2BW2OMD1W6q9uAtcBHwG9DCF0ody2HVc4xrvLY28DQpicPIXy3krGx8jljthWy0pH8APiQcic4BbgRGAW0DyHMAZYAY4AaYA1wOvAPwE8qnzcT6BVjPLpyzF8CVwIHAedR7v4age9sde6JwN9U1vS2GOOk7ayl1Cw7IFWjBcC/hhBqQgiHAV8PIewJzKDcEX0L6ACcTfmH7NcoF5a7KP8w/V0IYVwI4dJtHPuvgMVNH4gxrmny4bMxxsOBw4C/jjEeCRwN/DCE8BVgLPC/Yoy9KReAbRkbY+wJLKXczQEQQvhbYCBwJOXC9m8hhNDMOhwInFjJcmmM8WPKRej+GOMTwHjK25O9gOmUCybAXjHGo2KM1wB7hRAODCH8JbBfjPEl4BDgXypf2x+BbzfJuDfl9R0A9AU+biaf1Cw7IFWjCUAtMBt4AXgxxvhxCGFCjHEdQAhhKnBijHELMLLy2GWUu5MrKBeqH4UQDokxvtbk2MuAzsD6Tx4IIXwb+H3lw1h5+3dA9xDC85WP2wBdgK81Od4LlDuFpj6MMf628v4cytt98ysfH0q5qDxT+bgjcHCTc27t5RjjR8BHIYSN23i+FrizUsPaAIu2+hoA7gGGAJsoz9UA3gHuDSE0AN8A5n7y4hjjhkqHeRfwF8Dk7WSTPpMdkKpRD+CZyr/QHwIWV2Y7vw8h/PfKa46hfAECACGErwIhxjgbaE/5X+6NwJe3OvYE4MoQwpcqn3cIcDef/kt/S+Xtq8Bzle7in4EHgTeA5U22BHtsI/t/CyF0q7x/FPBKk+ci8Afg6Mpx6/i08G3Ltm7kuIVP/15HYEjlWJcCv9zqa4Dy1t2/Uu4M7w8h7ANcA5wCnAlspLyFB0ClU+oeY/wO5W3I//hkraQdZQFSNVoEjAwhzAV+THlm00j5B+ajIYRZlIvM+CafcwVwXeX9O4H/C/wl8LumB44xTgF+C/wmhPBryl3B4BjjO1tlmAY0hBBmUy50jTHGDZS3+u4LITxDuZvZ2ibgu5Vjf5Vycfvk3L+j3P38JoTw/yh3P8uLLwsALwP9KxdTnFPJ8hvKW3P/pZjFGBsqa7Cwkv89yp3bXMod5kbK25KfWAV0qsyYZgA3V7owaYfVeDdsSVIKdkCSpCQsQJKkJCxAkqQkLECSpCQsQJKkJFrt9fsvvviil+dJUga6d+9es63HW/Nl2K022M62cOFCamtrP/uFuznXqXnjx49n5cqVXHXVVamjtHp+LxWzk9Zpm8UHWnEHJGnHDB8+nIULF6aOIRXmDEiSlIQdkJSJXr16USqVmDdvXuooUiEWICkTw4YNY8WKFaljSIVZgKRMDBs2zBmQqoozICkTmzdvZvPmzaljSIXZAUmZOPbYY50BqapYgKRMnHnmmSxfvqP/fZCUjgVIysTgwYOdAamqOAOSMlEqldi4cWPqGFJhdkBSJvr16+cMSFXFAiRl4pxzznEGpKpiAZIyMXDgQGdAqirOgKRMrF+/ng0bNqSOIRVmByRlon///s6AVFUsQFImLrjgAt56663UMaTCLEBSJgYMGOAMSFXFGZCUiXfffZe1a9emjiEVZgckZeKkk05yBqSqYgGSMnHRRRdRX1+fOoZUmAVIysTxxx/vDEhVxRmQlIlVq1axevXq1DGkwuyApEyccsopzoBUVSxAUiZGjRrFm2++mTqGVJgFSMpE3759nQGpqjgDkjJRX1/PypUrU8eQCrMDkjJx6qmnOgNSVbEASZm44oorWLZsWeoYUmEWICkTvXv3dgakquIMSMrE4sWLvROCqoodkJSJ008/nVKpRJ8+fVJHkQqxAEmZuOaaa5wBqapYgKRM9OzZ0xmQqoozICkTMUaWLFmSOoZUmB2QlIkRI0ZQKpXo169f6ihSIRYgKRPXX389S5cuTR1DKswCJGXi8MMPp2PHjqljSIU5A5Iy8corr7Bo0aLUMaTC7ICkTJx//vmUSiVOOOGE1FGkQixAUiZuuukmr4JTVbEASZno0aMHHTp0SB1DKswZkJSJBQsW+Iuoqip2QFImRo4cSalUYsCAAamjSIVYgKRMjB49msWLF6eOIRVmASrozTUlzrh3PotXv0/X/b/MPUN7cMC+7VPHkv6sW7dutG3bNnUMqTALUBOnTZzHc3H1Z75u0TsNfOum57b53NFhfyae9k87O5r0mebPn8+SJUuora1NHUUqJOsC1Oens3jt7YYWPedzcTVdRv3qc3zmYg75Wgee/n7PnZ5Ju4dLLrmEUqnEySefnDqKVEjWBWhn/jA/9tZZvLG6gS2NsEcNHLR/B2ZcuHOOv3DhQv/Vqi/s9ttvdwakquJl2AXdM7QHB+3fgT1rajho/w7cM7RH6kjSf3LooYdy8MEHp44hFZZ1B7QzHbBv+53W8Ui7wpw5c1i6dKndtKqGBUjKxOWXX06pVGLQoEGpo0iFWICkTIwbN4433ngjdQypMAuQlIkQAlu2bEkdQyrMixCkTMyaNYv58+enjiEVZgckZeLqq6+mVCoxZMiQ1FGkQixAUiYmTJjA66+/njqGVNgX3oILIaxq5rleIYQpX/Qckj5b165d6dy5c+oYUmF2QDvAG5KqNZs5cybLli3z94BUNQoXoBDCocCtwJ7AfsA5McY5TZ5/HngV+AZQAwysPHVwCOFJ4KvAtBjjj0IIPYGrKXdgHYBBMcbXvviX88V90RuSejNSpXLddddRKpU444wzUkeRCqlpbGws9MIQwkDgjzHGl0MIg4CjY4zDQwirYoydKgVoQozxvhDCuZQL0aPAHcA/Ui5cb8YY96s8/3iMcUUI4XKgJsb4703Pt3Tp0saNGzd+oS/u7Kn1LFu3+QsdoyUd+JU2/Ly/Wyjb88EHH7DXXnuljtFqrVy5kk2bNtGlS5fUUVo9v5eK2RnrVFtbW7O953ZkC245cGUIYSOwN/DeNl7zbOXtHKB/5f1XYoybAEIIHzU51pgQQgPwdeCFrQ+0M/4SzdrJWxG76oak3oy0GNepebW1ta5RQa5TMbt6nXbkIoQxwNUxxqHAy5S32bbWvfL2COAPlfe31WKNB06LMQ4DVmznWK2ONyRVa/bUU08xe/bs1DGkwnakA5oMPBRCWAu8RXkOtLVhIYQLgfeBU4G/a+ZYs0MI7wNvA3+1AzmS8Yakas1uvPFGSqUSZ511VuooUiGFC1CM8VbKFyFs/XinJh9eFmN8tcnHz1f+/KfXxhgv3NGgkpo3ZcoUXnutVVzLIxXiZdhSJjp16sTatWtTx5AK22kFKMbYa2cdS9KOmzZtGvX19Q7XVTW8GamUiVtuuYW6urrUMaTC3IKTMvHwww87A1JVsQOSMrHffvvRsWPH1DGkwixAUiYeffRRZsyYkTqGVJgFSMrEmDFjmDx5cuoYUmHOgKRMTJ06lRhj6hhSYXZAUib22Wcf9t5779QxpMLsgKRMPPDAAyxfvtzfA1LVsAOSMjF27FimTPE/IFb1sAOSMjF9+nReffXVz36h1ErYAUmZaN++Pe3atUsdQyrMDkjKxOTJk50BqarYAUmZuPvuu3nkkUdSx5AKswOSMjFjxgwWLlyYOoZUmB2QlIk2bdrQpk2b1DGkwuyApEzU1dWxYsUKZ0CqGnZAUibq6up4/PHHU8eQCrMDkjLx/PPPOwNSVbEDkiQlYQGSMjF+/Hgeeuih1DGkwixAUiYeeOABnnzyydQxpMKcAUmZmDlzpjMgVRU7IElSEhYgKRN33nknv/jFL1LHkApzC07KxLRp02hoaEgdQyrMAiRl4sknn3QGpKriFpwkKQkLkJSJ2267jUmTJqWOIRXmFpyUiWeeecYZkKqKBUjKxBNPPOEMSFXFLThJUhIWICkTN998MxMmTEgdQyrMLTgpE3PnzmXDhg2pY0iFWYCkTDzyyCPOgFRV3IKTJCVhAZIyceONNzJ+/PjUMaTC3IKTMrFgwQLee++91DGkwixAUiamTJniDEhVxS04SVISFiApEz/+8Y8ZO3Zs6hhSYW7BSZmIMbJ+/frUMaTCLEBSJiZPnuwMSFXFLThJUhIWICkTV111FT/72c9Sx5AKcwtOykR9fT3r1q1LHUMqzAIkZWLixInOgFRV3IKTJCVhAZIycdlll3HrrbemjiEVZgGSMrFmzRp/D0hVxRmQlIm77rrLGZCqih2QJCkJOyApExdffDFr1qxh4sSJqaNIhdgBSZnYuHEjmzZtSh1DKswOSMrEHXfc4QxIVcUOSJKUhB2QlImRI0fypz/9ifvuuy91FKkQOyBJUhJ2QFImRo8e7QxIVcUOSJKUhB2QlInzzjuPtWvXcv/996eOIhViByRlol27drRt2zZ1DKkwOyApEzfffLMzIFUVOyBJUhJ2QFImzjrrLNatW8eDDz6YOopUiAVIysS+++5LY2Nj6hhSYRYgKRM33HCDMyBVFWdAkqQk7ICkTJx22mmsW7eOxx57LHUUqRALkJSJzp07065du9QxpMIsQFImrr32WmdAqirOgCRJSdgBSZkYPHgw69evZ9q0aamjSIVYgKRMhBBYvXp16hhSYRYgKRNXXnmlMyBVFWdAkqQk7ICkTJxyyim89957TJ8+PXUUqRALkJSJbt268c4776SOIRVmAZIyMWrUKGdAqirOgCRJSdgBSZk48cQT2bBhA08//XTqKFIhFiApE9/85jd5++23U8eQCrMASZm4+OKLnQGpqjgDkiQlYQckZeKEE06goaGBZ599NnUUqRALkJSJY445xhmQqooFSMrE9773PWdAqirOgCRJSdgBSZk47rjjaGhoYPbs2amjSIVYgKRMHH/88axatSp1DKkwC5CUiXPPPdcZkKqKMyBJUhJ2QFImevfuzfvvv8/cuXNTR5EKsQBJmRg4cCArV65MHUMqzAIkZWL48OHOgFRVnAFJkpKwA5Iy0atXL0qlEvPmzUsdRSrEAiRlYtiwYaxYsSJ1DKkwC5CUiWHDhjkDUlVxBiRlYvPmzWzevDl1DKkwOyApE8cee6wzIFUVC5CUiTPPPJPly5enjiEVZgGSMjF48GBnQNphb64pcca981m8+n267v9l7hnagwP2bd8i53YGJGWiVCqxcePG1DFURU6bOI9v3fQci95p4OPGRha900Df237dYue3A5Iy0a9fP2dAKqTPT2fx2tsN23yu9OHHdBn1qz9/fMjX3uHp7/fcJTksQFImzjnnHGdAKqRpQTn21lm8sbqBLY2wRw0ctH8HZlxYfn7hwoXU1tbushxuwUmZGDhwIMcdd1zqGKoy9wztwUH7d2DPmhoO2r8D9wzt0WLntgOSMrF+/Xo2bNiQOoaqzAH7tv9zx9PSLEBSJvr37+8MSFXFAiRl4oILLuCtt95KHUMqzAIkZWLAgAH+HpCqihchSJl49913Wbt2beoYUmF2QFImTjrpJGdAqioWICkTF110EfX19aljSIVZgKRMHH/88c6AtMNS3gsuSQEKIdQBNUBnoAMwJMb4aoosUi5WrVrF6tWrd+lvris/fW/7NaUPPwb4873g/nht3xY5d01jY2OLnKipSgFaHGO8NoTQDzg7xnhC09csXbq0cXe5seIHH3zAXnvtlTpGq+c6NW/o0KFs2bKFSZMmpY7S6u3u30tnT61n2bpi/3nhgV9pw8/7d/7c56qtra3Z3nMpt+CerbydA/x06ye7dOnSomFS2tX3W8qF69S8a6+9ljfffNM1KmB3/16a1eRr313vBde98vYI4A8Jc0hZ6Nu3L0cddVTqGKoyu+u94I4LIfQH9gSGJcwhZaG+vp6VK1fu1v+y147bXe8FNzrG+FTC80tZOfXUU/09IFUVL8OWMnHFFVewbNmy1DGkwpIUoBjjsBTnlXLWu3dvfw9IVcV7wUmZWLx4sXdCUFVxC07KxOmnn06pVKJPnz6po0iFWICkTFxzzTXOgFRVLEBSJnr27OkMSFXFGZCUiRgjS5YsSR1DKswOSMrEiBEjKJVK9OvXL3UUqRALkJSJ66+/nqVLl6aOIRVmAZIycfjhh9OxY8fUMaTCnAFJmXjllVdYtGhR6hhSYXZAUibOP/98SqUSJ5xwwme/WGoFLEBSJm666SavglNVsQBJmejRowcdOnRIHUMqzBmQlIkFCxb4i6iqKnZAUiZGjhxJqVRiwIABqaNIhViApEyMHj2axYsXp44hFWYBkjLRrVs32rZtmzqGVJgzICkT8+fP5+WXX04dQyrMDkjKxCWXXEKpVOLkk09OHUUqxAIkZeL22293BqSqYgGSMnHooYey5557po4hFeYMSMrEnDlzeOmll1LHkAqzA5Iycfnll1MqlRg0aFDqKFIhFiApE+PGjeONN95IHUMqzAIkZSKEwJYtW1LHkApzBiRlYtasWcyfPz91DKkwOyApE1dffTWlUokhQ4akjiIVYgGSMjFhwgRef/311DGkwixAUia6du3Kpk2bUseQCnMGJGVi5syZzJkzJ3UMqTALkJSJ6667jnHjxqWOIRXmFpyUiUmTJrFo0aLUMaTCLEBSJjp37kxDQ0PqGFJhbsFJmXjqqaeYPXt26hhSYRYgKRM33ngjd999d+oYUmFuwUmZmDJlCq+99lrqGFJhFiApE506dWLt2rWpY0iFuQUnZWLatGk899xzqWNIhVmApEzccsst1NXVpY4hFeYWnJSJhx9+2BmQqoodkJSJ/fbbj44dO6aOIRVmAZIy8eijjzJjxozUMaTCLEBSJsaMGcPkyZNTx5AKcwYkZWLq1KnEGFPHkAqzA5Iysc8++7D33nunjiEVZgckZeKBBx5g+fLl1NbWpo4iFWIHJGVi7NixTJkyJXUMqTA7ICkT06dP59VXX00dQyrMDkjKRPv27WnXrl3qGFJhdkBSJiZPnuwMSFXFDkjKxN13380jjzySOoZUmB2QlIkZM2awcOHC1DGkwuyApEy0adOGNm3apI4hFWYHJGWirq6OFStWOANS1bADkjJRV1fH448/njqGVFhNY2Nj6gzb9OKLL7bOYJKkHdK9e/eabT3eaguQJClvbsFJkpKwAEmSkrAASZKS8DLsBEII7YDJwFeBDcDQGOPqbbyuPTAHGBVjfKplU6ZVZI1CCDcBR1L+Pr4rxji+xYMmEkLYA7gT+HtgE3BmjPH1Js8PB0YAHwHXxRh/mSRoQgXW6PvAKZUPp8cYr2n5lOl91jo1ec2vgKkxxp/vrHPbAaVxDvByjPEo4D7giu287g5gd71KpNk1CiEcDfxNjPGblIvQD0IIHVs+ZjL/BuxV+fpHAbd88kQIoRNwAXAE8G3ghhBC2yQp02pujboC/xs4HDgM6BNC+B9JUqa33XVq4jpgp//9sgClcSTwSUfzJNB76xeEEC6m3P38rgVztSaftUZzgdMr7zcCewKbWyZaq/Dn9Ykx/hb4n02e+yfghRjjphjjeuB1YHf84drcGtUDfWOMH8cYG4E2wActH7FVaG6dCCGcBGzh07+PO41bcLtYCOEM4PtbPfw2sL7y/gZgn60+5xjg4BjjiBDCEbs+ZVqfZ41ijB8AH4QQ2gD3Ut6Ca9jVWVuRv+DT9QH4OITwpRjjR9t47r+s325iu2sUY9wMvBtCqAFuAl6KMb6WJGV6212nEMKhwCDgJOCqnX1iC9AuFmO8B7in6WMhhEeBvSsf7g2s2+rTzgAODCE8D3wD+McQwqoY44JdHDeJz7lGVLbcHgaejzHesKtztjLv8en6AOxRKT7bem6b67cbaG6NCCHsBUygXKDPbeFsrUlz6zQE+DrwLNAF+DCEsHRnzaQtQGm8APQD5gHHAbObPhljHPTJ+yGEOmBKrsWnGc2uUeUihWeAW2KM/6fl4yX3AnA88GAI4TDg5SbPzQP+vfIDti1QC7zS8hGT2+4aVTqfqcCzMcafJMrXWmx3nWKMl37yfgjhR8CqnXlBlAUojbHAvSGE3wAfUm5xCSH8B/BwjHFeynCtRLNrRHnA3hUYXrniC+C0GOOSFGETeAw4NoQwB6gBTgshXAi8HmN8IoQwhnLR3gP4YWXLcnez3TWiPDPsCbQNIRxXef1lMca5aaIm1ez30q48sbfikSQl4VVwkqQkLECSpCQsQJKkJCxAkqQkLECSpCQsQJKkJCxAkqQk/j/w2Sfx/3GhpQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.forestplot(model_single_intercept.trace, rhat=False)\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**gender model**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [beta_m, alpha]\n",
      "Sampling 2 chains: 100%|██████████| 8000/8000 [00:06<00:00, 1147.55draws/s]\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as model_gender:\n",
    "    alpha = pm.Normal('alpha', mu=0, sd=1.5)\n",
    "\n",
    "    beta_m = pm.Normal('beta_m', mu=0, sd=1)\n",
    "\n",
    "    logit_p = alpha + beta_m * data.is_male\n",
    "\n",
    "    p = pm.Deterministic('p', pm.math.invlogit(logit_p))\n",
    "\n",
    "    admit = pm.Binomial('admit', n=data.applications.values, p=p,\n",
    "                        observed=data.admit.values)\n",
    "    \n",
    "    model_gender.name = 'gender'\n",
    "    model_gender.trace = pm.sample(2000, tune=2000, chains=2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaMAAAEcCAYAAABqCdtUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHTBJREFUeJzt3XuYVXW9x/H3ZBwu4jESi+pgHk2/zXk4xYlDj5kmJppSSqEJh1DA+y0jL4WKkuZRCiUkFRGEUXgMUlSkgxQoEgo+EElpDV9QLo5cFInbuEdCmfPHWuRuHPYacJjfrB+f1/PwzN5r7ctnzWLv73zXd82estraWkREREL6SOgAIiIiKkYiIhKcipGIiASnYiQiIsGpGImISHAqRiIiEtxHQwcQaSxm1hKYABwBbAUud/flZvYd4A6gKr3pUGAx8CTQGrjY3f9sZscBX3X3n+3m8S8C+gE7gRbADe7+7F7kvAToANwH3OTul5nZKuDz7v5O0e0GpMsGN+AxK4DJ7j5zN+sPA77o7tP3NG/G83YDLnH3Po35uLL/UWckMbkQqHb3Y4DvA3eny7sAP3L3bum/ucApJMXoMuB8MysDfgDcVd8Dm1kf4GTgJHfvRlKUJppZ+70N6+7r3f2yvb3/Hvo68NUmei6RPabOSGLyH8BTAO7uZlaeLu8C/JeZDQIWAj8Gqkm6otbA20Bf4PHizqSOi4Gr3H1H+vgrzayzu29Mu5JD0n/fBH4EHA8cAIxw90fSrusuYBPwLvCCmR1O0s0ckz7HmHTZG0D/4ic3s++nGWvT+4yqL2TaqfwY+DtJhzgZGAYMBtqY2XxgJTAKKAM2AucB/wX8LL3fbKCbu5+YPuZvgBuBI4HLSbrCWuA7dZ57AvC59Ht6l7tP3M33UuQD1BlJTJYA3zKzMjM7BviMmR0AzCLplL4GtAUuIXnD/SRJkbmf5I31T2Y2xsx+VM9jfxpYUbzA3TcWXX3G3Y8FjgH+3d2PA04EbjCzjwGjgf9x9+4kxaA+o939BGAVSZcHgJn9B9AbOI6kyH3bzKzE9+GzwJlplh+5+3skBelhd38SGEtyCLMbMIOkeAK0cvfj3f1moJWZfdbMPgW0d/cXgaOBb6bb9lfgG0UZDyL5/vYCTgXeK5FP5APUGUlMxgPlwDzgeWCxu79nZuPdfTOAmU0DznT3ncCgdNl1JF3LEJKi9RMzO9rdlxU99mqgI7Bl1wIz+wbw5/Sqp1//E+hiZs+m11sAhwOfLHq850k6iGJ/d/cX0svzSQ4JLkqvdyIpME+n19sBRxU9Z10vufu7wLtmVlPP+nLg3rSetQCW19kGgAeAc4HtJHM4gDeBB82sGvg8sGDXjd19W9p53g/8KzBpN9lE6qXOSGLSFXg6/cn9EWBFOgv6s5n9W3qbk0hOXgDAzD4BmLvPA9qQ/ERfCxxY57HHAzea2UfT+x0NjOP9DmBn+nUpMCftOr4O/Bp4FVhTdNiwaz3Z/8XMOqeXjwdeLlrnwF+AE9PHreD9Ilif+j5wcifvv94dODd9rB8Bv6mzDZAc3vsWScf4sJkdDNwM9AEuAGpIDvMBkHZQXdz9OySHKn++63sl0hAqRhKT5cAgM1sA/JRkxlNL8ub5mJnNJSk4Y4vuMwS4Nb18L/Bb4FPAn4of2N0nAy8Az5nZ70m6hX7u/madDNOBajObR1L0at19G8nhwIfM7GmSLqeu7cD308f+BEmh2/XcfyLpip4zsz+QdEVrGv5tAeAloGd6IsalaZbnSA7ffaCwuXt1+j2oTPNvJenoFpB0njUkhy53WQ90SGdSs4A70u5MpEHK9KndIiISmjojEREJTsVIRESCUzESEZHgVIxERCQ4FSMREQku178HsHjxYp0KKCKSI126dCmrb3neT+2uBaisrKS8vDzrts2etqN5iWE7xo4dy7p167jppptCR2kUMewT2K+3o95CBDnvjESktAsvvJDKysrQMUQyaWYkIiLBqTMSiVi3bt0oFAosXLgwdBSRklSMRCI2YMAA1q5dGzqGSCYVI5GIDRgwQDMjyQXNjEQitmPHDnbs2BE6hkgmdUYiETv55JM1M5JcUDESidgFF1zAmjV7+qePRJqeipFIxPr166eZkeSCZkYiESsUCtTU1ISOIZJJnZFIxHr06KGZkeSCipFIxC699FLNjCQXVIxEIta7d2/NjCQXNDMSidiWLVvYtm1b6BgimdQZiUSsZ8+emhlJLqgYiUTsyiuv5PXXXw8dQySTipFIxHr16qWZkeSCZkYiEXvrrbfYtGlT6BgimdQZiUTsrLPO0sxIckHFSCRiV199NVVVVaFjiGRSMRKJ2Omnn66ZkeSCZkYiEVu/fj0bNmwIHUMkkzojkYj16dNHMyPJBRUjkYgNHjyY1157LXQMkUwqRiIRO/XUUzUzklzQzEgkYlVVVaxbty50DJFM6oxEInbOOedoZiS5oGIkErEhQ4awevXq0DFEMqkYiUSse/fumhlJLmhmJBKxFStW6BMYJBfUGYlE7LzzzqNQKHDKKaeEjiJSkoqRSMRuvvlmzYwkF1SMRCJ2wgknaGYkuaCZkUjE3J2VK1eGjiGSSZ2RSMQuvvhiCoUCPXr0CB1FpCQVI5GI3XbbbaxatSp0DJFMKkYiETv22GNp165d6BgimTQzEonYyy+/zPLly0PHEMmkzkgkYldccQWFQoEzzjgjdBSRklSMRCI2fPhwnU0nuaBiJBKxrl270rZt29AxRDJpZiQSsSVLluiXXiUX1BmJRGzQoEEUCgV69eoVOopISSpGIhEbOXIkK1asCB1DJJOKUSN4bWOB8x9cxIoNb3PEoQfyQP+uHHZIm9CxROjcuTMtW7YMHUMkk4rRHhg4YSFzfEPJ2yx/s5qvDZ/zgeUn2qFMGPjlfRVNpF6LFi1i5cqVlJeXh44iUtJ+X4xO+cVclr1Rvc+fZ45v4PDB/9eAWyaHVI7+ZFt+98MT9m0oid61115LoVDg7LPPDh1FpKT9vhg1xhv+ySPm8uqGanbWwkfK4MhD2zLrqj1/3MrKSv0EK43q7rvv1sxIckGndjeCB/p35chD23JAWRlHHtqWB/p3DR1JBIBOnTpx1FFHhY4hkmm/74waw2GHtNmrTkhkX5s/fz6rVq1Sxy3NnoqRSMSuv/56CoUCffv2DR1FpCQVI5GIjRkzhldffTV0DJFMKkYiETMzdu7cGTqGSCadwCASsblz57Jo0aLQMUQyqTMSidjQoUMpFAqce+65oaOIlKRiJBKx8ePH88orr4SOIZKpSQ7Tmdn6Euu6mdnkpsghsr854ogj6NixY+gYIpnUGTUCfVCqNFezZ89m9erV+j0jafYatRiZWSdgBHAA0B641N3nF61/FlgKfB4oA3qnq44ys6eATwDT3f0nZnYCMJSke2sL9HX3ZY2Zd0805ENSQR+UKs3LrbfeSqFQ4Pzzzw8dRaSkstra2kZ7MDPrDfzV3V8ys77Aie5+oZmtd/cOaTEa7+4PmdllJEXpMeAe4EskRew1d2+frn/C3dea2fVAmbv/b/HzrVq1qrampoZ33nmHVq1a7XXuS6ZVsXrzjr2+/7722Y+14L6e+TnU8mH3R3MRw3asW7eO7du3c/jhh4eO0ihi2Cew/25HeXl52e7WNfZhujXAjWZWAxwEbK3nNs+kX+cDPdPLL7v7dgAze7fosUaZWTXwGeD5ug+06wX2YT9gdO6HPIShD0r9Z9qO5qO8vDyK7dgllm3RdnxQY5/AMAoY6u79gZdIDsXV1SX9+lXgL+nl+tqzscBAdx8ArN3NYzUL+qBUaa5mzpzJvHnzQscQydTYndEk4BEz2wS8TjI3qmuAmV0FvA2cA/xniceaZ2ZvA28An27krI1GH5QqzdWwYcMoFApcdNFFoaOIlNSoxcjdR5CcwFB3eYeiq9e5+9Ki68+m//7ptu5+VWNmE9kfTZ48mWXLgp33I9JgOrVbJGIdOnRg06ZNoWOIZGrSYuTu3Zry+UT2d9OnT6eqqiqKYbnETR+UKhKxO++8k4qKitAxRDLpMJ1IxB599FHNjCQX1BmJRKx9+/a0a9cudAyRTCpGIhF77LHHmDVrVugYIplUjEQiNmrUKCZNmhQ6hkgmzYxEIjZt2jTcPXQMkUzqjEQidvDBB3PQQQeFjiGSSZ2RSMSmTJnCmjVr9HtG0uypMxKJ2OjRo5k8WX9IWZo/dUYiEZsxYwZLly7NvqFIYOqMRCLWpk0bWrduHTqGSCZ1RiIRmzRpkmZGkgvqjEQiNm7cOKZOnRo6hkgmdUYiEZs1axaVlZWhY4hkUmckErEWLVrQokWL0DFEMqkzEolYRUUFa9eu1cxImj11RiIRq6io4IknnggdQySTOiORiD377LOaGUkuqDMSEZHgVIxEIjZ27FgeeeSR0DFEMqkYiURsypQpPPXUU6FjiGTSzEgkYrNnz9bMSHJBnZGIiASnYiQSsXvvvZdf/epXoWOIZNJhOpGITZ8+nerq6tAxRDKpGIlE7KmnntLMSHJBh+lERCQ4FSORiN11111MnDgxdAyRTDpMJxKxp59+WjMjyQUVI5GIPfnkk5oZSS7oMJ2IiASnYiQSsTvuuIPx48eHjiGSSYfpRCK2YMECtm3bFjqGSCYVI5GITZ06VTMjyQUdphMRkeBUjEQiNmzYMMaOHRs6hkgmHaYTidiSJUvYunVr6BgimVSMRCI2efJkzYwkF3SYTkREglMxEonYT3/6U0aPHh06hkgmHaYTiZi7s2XLltAxRDKpGIlEbNKkSZoZSS7oMJ2IiASnYiQSsZtuuolf/vKXoWOIZNJhOpGIVVVVsXnz5tAxRDKpGIlEbMKECZoZSS7oMJ2IiASnYiQSseuuu44RI0aEjiGSScVIJGIbN27U7xlJLmhmJBKx+++/XzMjyQV1RiIiEpw6I5GIXXPNNWzcuJEJEyaEjiJSkjojkYjV1NSwffv20DFEMqkzEonYPffco5mR5II6IxERCU6dkUjEBg0axN/+9jceeuih0FFESlJnJCIiwakzEonYyJEjNTOSXFBnJCIiwakzEonY5ZdfzqZNm3j44YdDRxEpSZ2RSMRat25Ny5YtQ8cQyaTOSCRid9xxh2ZGkgvqjEREJDh1RiIRu+iii9i8eTO//vWvQ0cRKUnFSCRihxxyCLW1taFjiGRSMRKJ2O23366ZkeSCZkYiIhKcOiORiA0cOJDNmzfz+OOPh44iUpKKkUjEOnbsSOvWrUPHEMmkYiQSsVtuuUUzI8kFzYxERCQ4dUYiEevXrx9btmxh+vTpoaOIlKRiJBIxM2PDhg2hY4hkUjESidiNN96omZHkgmZGIiISnDojkYj16dOHrVu3MmPGjNBRREpSMRKJWOfOnXnzzTdDxxDJpGIkErHBgwdrZiS5oJmRiIgEp85IJGJnnnkm27Zt43e/+13oKCIlqRiJROwrX/kKb7zxRugYIplUjEQids0112hmJLmgmZGIiASnzkgkYmeccQbV1dU888wzoaOIlKRiJBKxk046STMjyQUVI5GI/eAHP9DMSHJBMyMREQlOnZFIxE477TSqq6uZN29e6CgiJakYiUTs9NNPZ/369aFjiGRSMRKJ2GWXXaaZkeSCZkYiIhKcOiORiHXv3p23336bBQsWhI4iUpKKkUjEevfuzbp160LHEMmkYiQSsQsvvFAzI8kFzYxERCQ4dUYiEevWrRuFQoGFCxeGjiJSkoqRSMQGDBjA2rVrQ8cQyaRiJBKxAQMGaGYkuaCZkUjEduzYwY4dO0LHEMmkzkgkYieffLJmRpILKkYiEbvgggtYs2ZN6BgimVSMRCLWr18/zYwi99rGAuc/uIgVG97miEMP5IH+XTnskDahY+0xFSORiBUKBWpqakLHkA9p4ISFzPENmbdb/mY1Xxs+Z7frT7RDmTDwy40ZrdGoGIlErEePHpoZBXbKL+ay7I3qetasaPIsc3wDhw/+vz2+39GfbMvvfnjCPkj0PhUjkYhdeumlmhkFVt+beGVlJeXl5Y3y+CePmMurG6rZWQsfKYMjD23LrKv2beHYF3Rqt0jEevfuzWmnnRY6huxDD/TvypGHtuWAsjKOPLQtD/TvGjrSXlFnJBKxLVu2sG3bttAxZB867JA2ueyE6lIxEolYz549NTOSXFAxEonYlVdeyeuvvx46hkgmFSORiPXq1Uu/ZyS5oBMYRCL21ltvsWnTptAxRDKpMxKJ2FlnnaWZkeSCipFIxK6++mqqqqpCxxDJ1OBiZGYVwGR3n9mA2x4GfNHdp3+IbCLyIZ1++umaGUVOn01X2teBzwMqRiIBrV+/ng0bNjTab/tL02vo59JBvj+brqy2trZBN0w7o48DB5IUsfPd/RUz+z7QF6gFJgP3AH8B2gBXAFuAoSQnS7QF+rr7st08x0+AzwHtgUPSxzoTOBro7+4vFN9+1apVtTU1Nbzzzju0atWq4VvdTGk7mpcYtqN///7s3LmTiRMnho7SKPK6Ty6ZVsXqzfn+I4ef/VgL7uvZ8Z+W7en+KC8vL9vduj3tjOa7+zAz6wH83MyGAL2B49L1s4DfAsOAz7v7k2Z2GdDP3dea2fXAd4H/LfEcNe5+qpkNBnq4++lmNhDoA/xTMTr88MOBxv2cp5C0Hc1LDNtxyy238Nprr+V+O3bJ6z6ZWydzLJ9N15jbsaendv8+/TofMKAT8Fng6fTfIcBRde6zBhiVdlYnAi0ynuOP6dfNwF/Ty5uA/P04JBLYqaeeyvHHHx86huxD++tn032ZpBAdD7wMOMkhudPcvdbMfgj8GejG+4VuLHCku28zsweB3bZpqYYdNxSRTFVVVaxbty6X3YQ0zP762XTHmNkzJAXjPHdfbWZPA8+ZWUtgIUkn9BJwg5n9EZgEzDOzt4E3gE83XnwRKeWcc87R7xlJLjS4GLn7gN0sHw4Mr7P4RZLDeJCc1NDQ5/hJ0eX7ii4/ATzR0McRkcSQIUNYvXp16BgimYL80quZPUZyZl6xLe7eM0QekVh1795dv2ckuRCkGLl7rxDPK7K/WbFiBVVVVZoZSbOnjwMSidh5551HoVDglFNOCR1FpCQVI5GI3XzzzZoZSS6oGIlE7IQTTtDMSHJBf89IJGLuzsqVK0PHEMmkzkgkYhdffDGFQoEePXqEjiJSkoqRSMRuu+02Vq1aFTqGSCYVI5GIHXvssbRr1y50DJFMmhmJROzll19m+fLloWOIZFJnJBKxK664gkKhwBlnnBE6ikhJKkYiERs+fLjOppNcUDESiVjXrl1p27Zt6BgimTQzEonYkiVL9EuvkgvqjEQiNmjQIAqFAr166bOJpXlTMRKJ2MiRI1mxYkXoGCKZVIxEIta5c2datmwZOoZIJs2MRCK2aNEiXnrppdAxRDKpMxKJ2LXXXkuhUODss88OHUWkJBUjkYjdfffdmhlJLqgYiUSsU6dOHHDAAaFjiGTSzEgkYvPnz+fFF18MHUMkkzojkYhdf/31FAoF+vbtGzqKSEkqRiIRGzNmDK+++mroGCKZVIxEImZm7Ny5M3QMkUyaGYlEbO7cuSxatCh0DJFM6oxEIjZ06FAKhQLnnntu6CgiJakYiURs/PjxvPLKK6FjiGRSMRKJ2BFHHMH27dtDxxDJpJmRSMRmz57N/PnzQ8cQyaRiJBKxW2+9lTFjxoSOIZJJh+lEIjZx4kSWL18eOoZIJhUjkYh17NiR6urq0DFEMukwnUjEZs6cybx580LHEMmkYiQSsWHDhjFu3LjQMUQy6TCdSMQmT57MsmXLQscQyaRiJBKxDh06sGnTptAxRDLpMJ1IxKZPn86cOXNCxxDJpGIkErE777yTioqK0DFEMukwnUjEHn30Uc2MJBfUGYlErH379rRr1y50DJFMKkYiEXvssceYNWtW6BgimVSMRCI2atQoJk2aFDqGSCbNjEQiNm3aNNw9dAyRTOqMRCJ28MEHc9BBB4WOIZJJnZFIxKZMmcKaNWsoLy8PHUWkJHVGIhEbPXo0kydPDh1DJJM6I5GIzZgxg6VLl4aOIZJJnZFIxNq0aUPr1q1DxxDJpM5IJGKTJk3SzEhyQZ2RSMTGjRvH1KlTQ8cQyaTOSCRis2bNorKyMnQMkUzqjEQi1qJFC1q0aBE6hkgmdUYiEauoqGDt2rWaGUmzp85IJGIVFRU88cQToWOIZCqrra0NnWGvLV68OL/hRUT2Q126dCmrb3mui5GIiMRBh+lERCQ4FSMREQlOxUhERILL3andZvYd4Lvu3reedRcCFwPvAre6+2/MrD3wMNAaWAsMdPdCU2auk7E1MAn4BLAN6O/uG4rWnwoMTq+WAccBnYBWwG+A5em60e4+paly15W1HeltpgHtgR1AjbufZmafAyqAWuBl4HJ339mU2etkbMh2DCfZDx8F7nf3sWb2cWAZyTYAPO7udzVd8n9k+whwL/BFYDtwgbu/UrS+2b8moEHb8UOgT3p1hrvfbGZlwOu8/5pY4O7XNWHsejVgW+4i+f+0LV3UE2hBjvaJmXUGRhbd/Bjg28BC9vJ1kavOKN2Jt1NPbjPrAFwJfBX4BnC7mbUEbgIedvfjgRdJXpghXQq8lOZ5CBhSvNLdZ7p7N3fvRlJ8fubulUAXYMSudSELUarkdqSOAo5L856WLhsBDEnvV0byQgyp5HaY2YnA59z9KyRvID82s3bAl4BfFe2PJi9EqW8DrdJ8g4E7d63I0WsCSm/HEcD3gGNJ3vROMbMvAEcCfyzaB8ELUWq325LqAnyjKPcWcrZP3H1J0fvUPcBUd5/Jh3hd5KoYAfNJ3jzq82XgeXffnu7cV4AvkLyBzExv8xTQfZ+nLK1Beczs34BzgJvTRV2Ab5rZ783sATML/ec7S26HmX0S+Bgw3cyeM7Nvpau6AHN3d78AsvbHAuC89HItcABJp9cF6GJmc83sETP7VFOErcc/8rv7C8B/F63Ly2sCSm9HFXCqu7/n7rUkXcQ7JPvgM2Y2x8xmmJk1dejd2O22pN3GUcD9Zva8mZ1X9z7kY58AYGYHkrxH/SBdtNevi2Z5mM7Mzgd+WGfxQHefYmbddnO3fwW2FF3fBhxcZ/muZU1iN9vxRgPzXAX8wt23p9cXAuPcfbGZ3QAMBa5p5Mj12svt+BeSn6TuAj4OPG9mC4Gy9A1ld/fbZ/ZmO9z9HeAdM2sBPEhymK7azJYCi919tpl9D/glcNY+3YD61f1//56ZfdTd361nXfDXRAm73Q533wG8lR6WGw686O7L0s7vdnd/xMyOIznc2rXpo39AqX1yIMn/lREkP9jMMbM/kLN9UrTsfOARd38rvb7Xr4tmWYzc/QHggT2821aguFs4CNhctLymaFmTqG87zOwx3s9Zb570p6dvATcULX7c3Xfd9nGSndwk9nI71gP3pf9x3zSzFwEDiudDedkf7YBHgWfd/fZ08TPArmP6jwO37IvMDVD3//1Hit4smt1rooRS24GZtQLGk7xRX5Yu/gPJLAx3f87MPm1mxT/shFJqWwrAXbvmQWb2DMlMJnf7JPU9/rnY7PXrIm+H6UpZCBxvZq3M7GCgnGSI9jzQI73NacC8QPl2aUieTsBSd68pWvZbM/tyevkkYPG+i9ggWdvRHXgEwMzakmxTJfBiUXfb7PdHeoLD08B4d/9p0apxwJnp5ZD74x/5zewY4KWidXl5TUCJ7Ug7omnAn9z9Ynd/L101FBiU3uaLQFUzKERQep8cTXKU4IC02z4O+CM52yfpsoOBlu5eVbR4r18XufsEhvSN7BJ375Nevwp4xd2fTM8cuoikyN7m7lPT2cWDJBX+LaCvu78dJj2YWZs0z6eAv6d51pvZz4FH3X2hmX0X+Kq7Dyq635dIuqEdJF3HRe6+tem34B95GrIdI0kGzjuBn7v7E2Z2NDCW5DBeJXBh0ZtLk8vaDpLh/1BgSdHdBqZfx5OchPE2yZlG65oseKrojKcvpFkGkryB5OY1AaW3g+Rw1q+AF4ruch3JIaFJQFuSDulydw/+N9YbsE+uBc4meS0/5O735W2fpNvRFbjB3b9ddJ9/Zy9fF7krRiIiEp+YDtOJiEhOqRiJiEhwKkYiIhKcipGIiASnYiQiIsGpGImISHAqRiIiEtz/A/bDhU3xGULzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.forestplot(model_gender.trace, varnames=['alpha', 'beta_m'], rhat=False)\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/petrushev/etc/jup_env/lib/python3.7/site-packages/pymc3/stats.py:218: UserWarning: For one or more samples the posterior variance of the\n",
      "        log predictive densities exceeds 0.4. This could be indication of\n",
      "        WAIC starting to fail see http://arxiv.org/abs/1507.04544 for details\n",
      "        \n",
      "  \"\"\")\n",
      "/home/petrushev/etc/jup_env/lib/python3.7/site-packages/pymc3/stats.py:218: UserWarning: For one or more samples the posterior variance of the\n",
      "        log predictive densities exceeds 0.4. This could be indication of\n",
      "        WAIC starting to fail see http://arxiv.org/abs/1507.04544 for details\n",
      "        \n",
      "  \"\"\")\n"
     ]
    },
    {
     "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>WAIC</th>\n",
       "      <th>pWAIC</th>\n",
       "      <th>dWAIC</th>\n",
       "      <th>weight</th>\n",
       "      <th>SE</th>\n",
       "      <th>dSE</th>\n",
       "      <th>var_warn</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>gender</th>\n",
       "      <td>984.38</td>\n",
       "      <td>113.01</td>\n",
       "      <td>0</td>\n",
       "      <td>0.65</td>\n",
       "      <td>288.14</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>single intercept</th>\n",
       "      <td>1050.27</td>\n",
       "      <td>86.37</td>\n",
       "      <td>65.89</td>\n",
       "      <td>0.35</td>\n",
       "      <td>293.19</td>\n",
       "      <td>158.95</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                     WAIC   pWAIC  dWAIC weight      SE     dSE var_warn\n",
       "gender             984.38  113.01      0   0.65  288.14       0        1\n",
       "single intercept  1050.27   86.37  65.89   0.35  293.19  158.95        1"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pm.compare({\n",
    "        model_gender: model_gender.trace,\n",
    "        model_single_intercept: model_single_intercept.trace,\n",
    "    },\n",
    "    method='BB-pseudo-BMA',\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjQAAAGACAYAAAC6OPj9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAANEgAADRIBtYA3dAAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl4FEX+x/H3kIRwoyiXXKKYUjBgRMQFUUBlEURRRMXVVRcRFgFPXC9ARVgR1huPxf3pKqiIKyCKCLq6XvEgBgkepSLIoSCnyBEISf/+qAmZHMAQZtL08Hk9T56kuyvT3xQh+aSqujvkeR4iIiIiQVbJ7wJERERE9pcCjYiIiASeAo2IiIgEngKNiIiIBJ4CjYiIiASeAo2IiIgEngKNiIiIBJ4CjYjsF2PMfcaY98r5uc8aY14Kfxwyxgw0xqTEsLY2xpjTwx93NsZ4xpgqsXp9ETlwKNCIiJ+uAwaGPz4NeBJIiuHrzwCOC3/8MdDQWpsbw9cXkQNEst8FiMjBy1r7W8RmKA6n2PWa1todwKo4nENEDgAhPfpARCIZY9oD44B2uFHcT4FB1tpvw8f/ADwKtALeB5YCxlrb2RhzJXA1MBO4NfyS9wJZuNGXpsBbwJ+stbnGmGeBKuG2SyLK6GKtfa9EXe8BXwNnATWAtkBN4AHg1PDrLASus9Z+HG5/evjT/w08C7wLVA2fuyEwAfgjblRoBnB9iZAlIgGhKScR2cUYUxN4E/gQOB43DVQDuD98/PDw8Y+AE8IfX13iZdqFj50CPAiMBx4KtzsP6AZcUeJzlgN9wh8fiZseKstVQH/gXOAXYBZu1OWk8NtaXHACuABYAfwNN7UV+XWmAO8A9YAzgXNwAemZ3ZxXRA5wCjQiEqk6MBYYYa1dYq39HDey0Sp8/GJgE3CjdR4C3ijxGinAtdba74GJuNGPx6y1H1tr38GN6rSM/ARrbT6wPry5Ojw9VJa3rLXvh+uqBkzCjap8b63NCZ+vVfg11wP5wKYyRl26A0cBl1prF1hrP8KFrPONMWlR9JOIHGAUaERkF2vtKtwoxXXGmGeMMZnAPyhaqNsSyAkHkEKfl3iZ9dbajeGPt4XfL404vg1ILWeJu17HWrsFeBy41BjzT2PM+8BLRPdz7ThgibV2TcTrZYdra7nbzxKRA5YWBYvILsaYRsB84EtgLvAckAEMi2hWcvFuydGUnWW8dEGMStx1hZIxpgYuTG3ErX95FaiNCzVRv07E64VwX1ssr7ISkQqiERoRiXQ+sNVa291a+4C19l3cQt7CELMIOKHEvWIyYnTufb1CoTNuvc3p1tpx1to5QEPYFU729JrfAs2NMfUi9mXgFhbbfaxDRA4AGqERkUjrgIbGmD/ifrH3BP4KbAgffwm4C3jcGDMBdxXR+UBmDM69Ofz+RGPMF1HcL2YdLoBcaIz5COgIjAofS8WNwmwGjDGmTonPfRt3xdQUY8xwoCrwBPCetXbR/n8pIlLRNEIjIpFexi0CfhFYAPTFBZp6xphG1toNuKuUWoaP/xl4LEbnzsFd0v0u0GNvja21mbgA8wDwFXB9uNZ8ikaNHsNdXfV0ic8twF1xtQV3xdbruOmr3jH4OkTEB7oPjYiIiASeRmhEREQk8BRoREREJPAUaERERCTwFGhEREQk8BRoREREJPAq/D40WVlZKUAa7rktsbp7qIiIiCSuSkAd4Lu2bdvmldXAjxvrpeHuNioiIiKyL47H3XeqFD8CzXqA4447jpSUlL213W85OTmkp6fH/TwHG/VrfKhf40P9Gh/q1/hQv5aWl5fHN998A+EMURY/Ak0BQEpKSoUEmlAoVCHnOdioX+ND/Rof6tf4UL/Gh/p1j3a7VEWLgkVERCTwFGhEREQk8BRoREREJPAUaERERCTwFGhEREQk8BRoREREJPAUaERERCTwFGhEREQk8BRoRERE4qWgAGbPhmuvdW9vvAH5+RVy6ssvv5y33367Qs61r+JRmwKNiIhIPGzZAmecAZdcAr//Dps3w6WXQteubltiyo9HH0gJmzfDu+/Czp3QqRMcfrjfFYmIyH677TZYvx5+/LHoB/v69S7k/O1v8PjjMTvV5MmTefHFF6lUqRLNmjVj1KhRALz55ps8/vjjbN++nb59+3LllVeyZcsWbrnlFlauXEkoFKJz585cd911eJ7HxIkTeeeddygoKOCYY45h1KhR1KxZk65du9K6dWu+/fZbrrzySjIzM3n44YcByMzMZMKECfznP/9h3rx5PPXUU+Tn51OzZk1GjhxJixYtWL16Nbfeeitr1qyhcePGbNy4MWZfeyGN0PjskUfgiCOgf38YMgQaN3bf5wW7fVqFiIgc8LZtg2eegXHjiv+VWqcO3H8/PPccbN0ak1N9/PHHTJ48mRdeeIFZs2Zx3HHHcccddwCQm5vLtGnTeOGFF3j++efJzs5m3rx5pKamMmPGDF544QV++uknNm/ezPTp01m+fDmvvPIKM2fO5Mgjj2T8+PG7znP88cczZ84czjnnHD755BM2bdoEwIwZM+jTpw9Llixh4sSJPPPMM0yfPp2hQ4cyZMgQAO655x7atWvH66+/zvDhw1m6dGlMvvZIGqHx0ZQpLsA/+yxccAFUquRGavr1g5o14c47/a5QRETKZe1aN/zerl3pY+3auemoX3+FI4/c71N9+OGH9OzZk9q1awNufcoTTzxBeno6F198MUlJSdSuXZszzjiDjz76iHPPPZcHHniA/v3706FDB2688UZq1KjBe++9x6JFi7jgggsA2LlzJ3Xq1Nl1nhNPPBGAGjVq0LlzZ95880169erF+++/zx133MGsWbP45ZdfuOyyy3Z9zubNm9m4cSMff/zxrpB19NFH066sftlPCjQ+8Ty47z4YNQr69i3a37UrTJwIAwfCzTdDlSr+1SgiIuV0+OFQvTpkZUG3bsWPZWVB1apQt25MTlVQYki/oKCA/PDC46SkpF37Pc8jJSWFpk2bMnfuXDIzM/n444+56KKLeOKJJ/A8j7/+9a/0Df9S2rp1K7m5ubs+v2rVqrs+7tOnDw899BBVqlShQ4cO1KpVi4KCAjp16sSECRN2nW/VqlXUrl2bUChUrM7IumJFU04+2boVFi2C3r1LHzvvPDfN+t13FV+XiIjEQNWq8Oc/u2H4DRuK9m/cCLfeCpdf7gJPDHTs2JE33niD3377DYApU6Zw8sknk5yczKuvvorneaxfv563336bU089lSlTpnD33XfTpUsXbr/9dlq0aMGSJUvo0KED06ZNY+vWrXiex+jRo/nHP/5R5jnbtWvHunXreO6557jwwgsBaN++PR988AHLli0DYPr06Vx55ZUAdOrUiVdeeQWA5cuXM3/+/Jh87ZE0QuOTypUhJQVWr4a0tOLHfv3Vva9Zs+LrEhGRGBk3Drp3h6OPhj59IBSCV14BY9w6mhjp1KkTl1xyCX/605/Iz8+nSZMmjBs3juHDh1O9enUuuOACtm/fzuDBg2nVqhXNmjUjMzOTHj16kJqayrHHHkuPHj1ITk5m1apV9O3bd9ei4DFjxpR5zlAoxAUXXMDUqVM55ZRTAEhLS2PkyJEMGTIEz/OoXr06jzzyCKFQiFGjRnHbbbfRs2dPGjZsiDEmZl//rpo8z4v5i+5JVlZWfWBV69atSUlJifv5srOzycjIiPt5yqNfPzeNOmOGWz9TaNgw+Owz+OQT/2rbmwO5X4NM/Rof6tf4UL9GYedOeP119wbQsyf06gXJux9PUL+WlpeXx8KFCwEatG3bdnVZbTRC46OxY6FDB3cF38CBkJoKzz8Pc+fCf//rd3UiIrLfkpPd2oKy1hdITGkNjY+aN4fPP4c2bdyU6tChcMghbt/JJ/tdnYiISHBohMZnjRvDQw+5NxERESkfjdCIiIhI4EU9QmOMORcYA6QCmcBAa21uiTa3AFcCBUBWWW1EREREYi2qERpjTH1gEtDLWpsGbAFuL9EmA7gWOBlIB2oCg2NarYiIiEgZop1y6gZkWmuXhrefAi4r0aYSUBmoihv5qQpsi0GNIiIiInsU7ZRTI2BFxPZKoHFkA2ttljHmJWAZbgRnMfDM7l4wJyeHUCi0b9WWQ0FBAdnZ2XE/z8FG/Rof6tf4UL/Gh/p1z5I3bCBpwwa2H3VUsf2pS5aQX7s2OyOekxRJ/VpaNPfMizbQlDWSkx+5YYzpCZwONAU2AI8AT+LW1JSSnp5+0N9YL8jUr/Ghfo0P9Wt8qF/3YPt2d/+NX35xTx1u1crt//prGDzYPcfp88/LfGCf3/366KOPsmnTJu644w4GDBjAbbfdxlElQlmh33//nWHDhvHMM7sdvyjTrbfeyrHHHrvr0Qh7E3Fjvd2KdsppOdAwYvsIio/YAPQApltr11hrdwKPA2dG+foiIiKJIzUV+veHNWugSxf46isXZrp2dc+8+ctfAvH04UmTJu02zAD89ttv5OTkVGBFuxftCM1cYLwxprm1dgkwAJhZok0WcI0xZoK1dhtwPvBZ7EoVEREJkGHD3PvrroPjjy/a/8ADcMMNMTvNp59+yoQJEzjssMNYuXIlNWrU4P777+f222+nVq1aLFmyhAEDBtCqVSvGjBnDpk2bCIVCDB48mDPPPJO8vDzuvfdeMjMzqVOnDrVr16Zp06YAdO3alYkTJ3LccccxdepU/v3vf1OpUiUaNWrEuHHjuPPOO9myZQt9+/Zl2rRpzJ8/nwkTJrB9+3ZSU1O5+eabOemkk9i8eTN33HEH1lrq1atHfn4+xx57bMz6AKIMNNba1caYAcBMY0xlIAe4yhhzEnCPtbYHbr1Mc2CBMWY78APw15hWKyIiEiTDhsE338CTT7rt/v1jGmYKffXVV0yZMoWMjAz+9a9/MWLECADq1avHxIkT2blzJ7179+bBBx/kmGOOYe3atVx00UUcd9xxvP3226xcuZLZs2eTm5tLv379dgWaQt988w2PPPII06dPp169ejzwwANMmTKFe++9l969ezNt2jQ2btzIiBEjePbZZ6lfvz5LlizhiiuuYM6cOTzyyCPUrFmTOXPmsGbNGnrH4VEQUd+Hxlo7C5hVYvd83FQT1loPGBF+ExERka+/hunTi7Zfe81NPxWuqYmRli1b7lp3c9FFFzFhwgTatGnDiSeeCMCSJUv46aefuPnmm3d9jud5fPfdd2RmZnLOOeeQnJxMjRo16NmzJ+vWrSv2+p988gmnnXYa9erVA+DGG28EYMWKotUnCxYsYPXq1VxzzTXFPnfZsmVkZmYycuRIAOrWrUvXrl1j+vWDHn0gIiISH5FrZh54AJKS3PRTly7FFwrHQHLE07s9zyMUClGpUiWqVq0KuCunDj30UGbOLFotsnr1aurUqcPUqVMpKCjYtT8pKanU65fct3nzZjZs2FDsauWCggKMMbz44ou79q1atYq6devuOr6nc+wvPfpAREQk1rZvh7PPLgozN9zgpp8eftgtFO7eHXJjdyP9RYsW8cMPPwAwdepUTjvttGKhoXnz5iQlJTFnzhwAFi9ezNlnn83q1avp1KkTr732Gjt27GDbtm272kQ6+eST+eijj3aN3Dz99NM8/fTTJCcnk5/vLnpu06YNixcvZsGCBYAb1enVqxe5ubl06tSJV199Fc/z2LBhA++9917MvvZCGqERERGJtdRUeOIJ+P57NypTaNgwCIWgefOYXuVUr149xo8fz8qVK2nQoAFjx47lpptu2nW8cuXKPPbYY4wZM4aJEyfieR5jx46lcePGXHLJJSxbtoxzzjmHQw45hCZNmpR6/WOPPZbrr79+12XWzZo14+9//zvVqlWjZcuWnHPOOUybNo2HH36Y0aNHk5eXR1JSEo899hjVq1dn6NChjBo1iu7du1O3bl2OOeaYmH3thULR3KwmlrKysuoDq1q3bq370ASY+jU+1K/xoX6ND/VrfOxrv3766aeMHTu22HRSoom4D02Dtm3bri6rjaacREREJPAUaERERAKsffv2CT06Ey0FGhEREQk8BRoREREJPAUaERERCTwFGhEREQk8BRoREREJPAUaERERCTwFGhEREQk8BRoREREJPAUaERERCTwFGhEREQk8BRoREREJPAUaERERCTwFGhEREQk8BRoREREJPAUaERERCTwFGhEREQk8BRoREREJPAUaERERCTwFGhEREQk8BRoREREJPAUaERERCTwFGhEREQm85GgbGmPOBcYAqUAmMNBamxtxfDjwp4hPaQj8aq1Nj1GtIiIiImWKaoTGGFMfmAT0stamAVuA2yPbWGvHW2tPsNaeAHQLt7kytuWKiIiIlBbtlFM3INNauzS8/RRw2R7aPwQ8ba3N2o/aRERERKIS7ZRTI2BFxPZKoHFZDY0xJwCdgKv29II5OTmEQqEoT19+BQUFZGdnx/08Bxv1a3yoX+ND/Rof6tf4UL+W5nneXttEG2jKGsnJ303bocDD1trte3rB9PR0UlJSojx9+WVnZ5ORkRH38xxs1K/xoX6ND/VrfKhf40P9WlpeXh4LFy7cY5top5yW4xb5FjqC4iM2ABhjkoDzgMlRvq6IiIjIfos20MwFOhpjmoe3BwAzy2iXjruyaVUsihMprx07YMoU6NMHevWChx6C337zuyoREYmXqAKNtXY14RBjjPkWaADcZYw5yRgzO6LpMcCy2JcpEr3Nm6FrV7jpJmjeHDIy4MknoU0bWLrU7+pERCQeor4PjbV2FjCrxO75QI+INtOAabEpTaR8Ro+GTZvgm2/g0EPdvpEj4eKLYeBAeOstf+sTEZHYizrQiASB58G//gVPPFEUZgCSk+G++yAtDZYtg6ZN/atRRERiT48+kISSlwfr1kGrVqWPtWgBlSvDKq3wEhFJOAo0klAqV3brZj78sPSxzz93IzjNm5c+JiIiwaZAIwln6FAYNQq+/rpo35o1MGQIXHIJ1K3rX20iIhIfWkMjCee661yYSU+HLl2genWYOxfat4fHHvO7OhERiQeN0EjCqVQJJk2CBQvgtNOgdWuYPRvefRdq1fK7OhERiQeN0EjCSk93byIikvg0QiMiIiKBp0AjIiIigadAIyIiIoGnQCMiIiKBp0AjIiIigadAIyIiIoGnQCMiIiKBp0AjIiIigadAI/smPx+ysqi2cCHk5vpdjYiICKBAI/ti2jQ4+mg45RSOueYaaNwYHnjAPcJaRETERwo0Ep0334TLLoMbboCNG/nygw/ckx7vuQceftjv6kRE5CCnZzlJdEaPhptuco+yBkhJgUsugZ073f7Bg6FyZX9rFBGRg1ZijtAUFMC8eTB0KI3HjoUZM9wvXimf7dshMxP69St9rG9fWLMGvv664usSEREJS7xAk5sLPXvC+efDhg1U2r4d/vIX6NgRNm70u7pgqlQJkpLg999LH9uyxa2hqVKl4usSEREJS7xAc/fdsGQJ/PADTJ7MssLtpKSi6RLZNykpcN558OijpRcAP/YYtGwJxvhTm4iICIkWaPLzYdIkuPdeaNCgaH/t2vCPf8BLL8GGDf7VF2SjR8OcOXDRRfDuu+6y7aFDYcwY17ehkN8ViojIQSyxAs2mTbBuHbRrV/pYu3awYwesXFnxdSWCli3dOprkZDj7bNL694fvvoN33oHu3f2uTkREDnKJdZVTzZpw6KHwxRfQrFnxY1lZburkiCP8qS0RHHssvPgieB4LsrPJOPFEvysSEREBEm2EJjkZ+veHESPclTeFNm+G4cPdFTl16vhXX6IIhTTFJCIiB5TEGqEBtyg4OxtatIA+fWiyaRO8/767q+0jj/hdnYiIiMRBYo3QAFSrBm+9BS+8AMnJhHbudFfifPIJHHaY39WJiIhIHEQ9QmOMORcYA6QCmcBAa21uiTZdgfuAKsAK4HJr7brYlRulpCR3L5qePVmWnc1hGRkVXoKIiIhUnKhGaIwx9YFJQC9rbRqwBbi9RJsmwMvAVdba1sBXwD2xLVdERESktGinnLoBmdbapeHtp4DLSrS5EHjdWvtVePteYMJ+VygiIiKyF9FOOTXCTSEVWgk0LtHmGCDXGPMK0AI3QnP97l4wJyeHUAVcKVNQUEB2dnbcz3OwUb/Gh/o1PtSv8aF+jQ/1a2leybvUlyHaQFPWSE5+ie0U4GzgVOAn3Fqap4ALynrB9PR0UlJSojx9+WVnZ5OhNTQxp36ND/VrfKhf40P9Gh/q19Ly8vJYuHDhHttEO+W0HGgYsX0ExUdsAH4G3rfWLrXWesBzQPsoX19ERESk3KINNHOBjsaY5uHtAcDMEm1mAKcbYwqnos4DPt//EkVERET2LKpAY61dTTjEGGO+BRoAdxljTjLGzA63yQZuBGYbY74GTgf+Gp+yRURERIpEfR8aa+0sYFaJ3fOBHhFtpgJTY1OaiIiISHQS707BIiIictBRoBEREZHAU6ARERGRwFOgERERkcBToBEREZHAU6ARERGRwFOgERERkcBToBEREZHAU6ARERGRwFOgERERkcBToBEREZHAU6ARERGRwFOgERERkcBToBEREZHAU6ARERGRwFOgERERkcBToBEREZHAU6ARERGRwFOgERERkcBToBEREZHAU6ARERGRwFOgERERkcBToBEREZHAU6ARERGRwFOgERERkcBToBEREZHAU6ARERGRwFOgERERkcBLjrahMeZcYAyQCmQCA621uSXazAGaA9vCu/5trX0wRrWKiIiIlCmqQGOMqQ9MAtpba5caYx4HbgdGRrQJAW2Bo621m+JRrIiIiEhZop1y6gZkWmuXhrefAi4r0ea48PupxpgcY8zDxpiqMahRREREZI+inXJqBKyI2F4JNC7R5jDgHeCvwFbgeWAscENZL5iTk0MoFNqnYsujoKCA7OzsuJ/nYKN+jQ/1a3yoX+ND/Rof6tfSPM/ba5toA01ZIzn5kRvW2g+ADwq3jTH3Aa+wm0CTnp5OSkpKlKcvv+zsbDIyMuJ+noON+jU+1K/xoX6ND/VrfKhfS8vLy2PhwoV7bBPtlNNyoGHE9hEUH7HBGNPFGNM9YlcIyIvy9UVERETKLdpAMxfoaIxpHt4eAMws0eYQ4EFjTA1jTCXgemBabMoUERER2b2oAo21djXhEGOM+RZoANxljDnJGDM73GY68CLwOfAtkAuMjkvVIiIiIhGivg+NtXYWMKvE7vlAj4g29wD3xKY0ERERkejoTsEiIiISeAo0IiIiEngKNCIiIhJ4CjQiIiISeAo0IiIiEngKNCIiIhJ4CjQiIiISeAo0IiIiEngKNCIiIhJ4CjQiIiISeAo0Ptq5E379tfT+TZtg8+aKrydR5OfD6tWl92/e7PpWyqegAFatKr1/61bYuLHi60kUnge//FJ6f24urF9f8fUkit31644dsHZtxdeTKDwPfv659P68PFizpuLriaRA45P8fLj8cujQAVasKNr/22/wxz9Cz56wZYt/9QVVQQH07w+nnAJLlxbt37wZevSA7t3h9999Ky+wPA8GDYJ27eCHH4r2b90KvXrBWWcp1JSH58GwYXDiiWBt0f7cXOjdG7p2Vagpr1tugTZt4Kuvivbt2AEXXgidO/v/yzeoRoyA9HRYsKBoX14e9OsHp55a9h89FUWBxieVKkGTJrB4sfvPtWKFCzPdu8Mnn0CDBpCa6neVwRMKQdOmLsx06eLeF4aZDz5Qv5ZXYb+uWOH69YcfisLMf/8L9epB1ap+Vxk8hf26apX7OWBtUZh56y316/5o2tSFli5dXKgpDDOzZkHdulCtmt8VBlOTJi5kn3GGCzWFYeY//4HDD4fq1f2rLeqnbUtshUIwbpz7ePx49xdaaqr7hXHRRTBlCiTrX2efhUJw993u49Gj3YhCtWqwbBmcfz689BJUruxvjUF1553u/YgR0L491KwJP/3kwuKrryooltfw4e79LbfAH/4AtWu7IH7WWTBzpgJNeQ0d6t4PGwYdO7pftoV/QL7+ur+/eINs4ED3ftAgOP10Fw4XL3azDXPmuJ8LftGvTB8VhppNm+Cpp9y+jh0VZvZXYajZvBkefNDtO/FEhZlYuPNO2LYNxo51f6UZozATC8OHu5GZkSNhwwY3uqAws/+GDnUjMzff7EbA69ZVmImFgQPdGtAhQ9zvr+rV/Q8zoCkn323aVHwucsUKf+cgE8WWLTB/ftH2qlVlL2STfbN1q5sSLbR+PSxf7l89iSI3Fz76qGh70yY3qij7Z8cO+N//ira3bi2+tk7KJy8P3n23aHvHDjdK4zcFGh8Vrpn59FM3zXTDDW4Iv3BNjZRP5JqZ88+Hv/3NhZnCNTVSPpFrZnr0cFNPhWsUIhcKy76JXDNz1llwzz1ugXXhmhopn8g1M507w9//7v7QKVxTI+UTuWamQwe3ZCIvr2hNjZ80seGT/Hz3S+GTT4rWzCQluamm8ePdf8CsLDefLtErKHC/dAvDzEsvQUqKm2oaPdr9MMvKgjp1/K40WDwPLrigKMwUTjNVruyCTWG/1qvnd6XB4nnu/39hmCmcZqpSxa2pKfw5cMQRflcaPJddVhRmCqeZqld3a2q6dHEjuE2b+l1l8PzlL0VhpnCaqWZNt6bmjDPg88/hqKP8qU2BxidJSe7y4qZN4fnni9bMFC4UrloVatXyr76gqlQJrr7aLQCcMqVozUzhQuEdO+DQQ/2rL6hCIff9WqUKTJ1atGamcKHwunVufYLsm8J+LSiAadOK1swULhRessRdmSf77qqr3C0aXnmlaM1M4ULhL7+Exo39qy3IrrjCTeG/+mrRmpnChcIffQTNmvlXW8jzvAo9YVZWVn1gVevWrUlJSYn7+bKzs8nIyIj7ecrL89wPtZL7oPT+A4n6NT6C2K972n+gUL/Gh/o1PtSvpeXl5bFw4UKABm3bti3j1qkaofFdWf/4B/J/tKBQv8bH7vpQfbt/1K/xoX6NjwO1X7UoWERERAJPgUZEREQCT4FGREREAk+BRkRERAJPgUZEREQCT4FGREREAk+BRkRERAIv6vvQGGPOBcYAqUAmMNBam7ubthOAFtba3jGpUkRERGQPohqhMcbUByYBvay1acAW4PbdtO0FXBGzCkWalJrJAAAa1ElEQVRERET2Itopp25AprV2aXj7KeCyko2MMc2AkcDomFQnIiIiEoVoA00jYEXE9kqg2KO9jDEpwPPAYGBTTKoTERERiUK0a2jKCj75JbbHATOttZ8bY1rt7QVzcnIIVcCDHwoKCsjOzo77eQ426tf4UL/Gh/o1PtSv8aF+LS2aB2lHG2iWA20jto+g+IgNQD/gV2PM5UAdoLYx5hVr7YVlvWB6erqeth1g6tf4UL/Gh/o1PtSv8aF+LS3iadu7FW2gmQuMN8Y0t9YuAQYAMyMbWGsbFn5sjLkS6L27MCMiIiISS1GtobHWriYcYowx3wINgLuMMScZY2bHs0ARERGRvYn6PjTW2lnArBK75wM9ymj7LPDs/hQmIiIiEi3dKVhEREQCT4FG5ECxfTuhbdv8rkJEJJAUaET8lpMD550H1apxQseOcMop8MYbflclIhIoCjQifvrqKzj1VKhVCzIz+eaVV+Css+CCC+Dll/2uTkQkMKJeFCwicXDXXdCjBzz3HIRC5KakQJ8+cPjhcMst7uOkJL+rFBE54GmERsQvngevvQaDBkHJu2ZffTWsWAFffulPbSIiAaNAI+IXz4OdO6FKldLHUlLcyExeXsXXJSISQAo0In6pVAnOPBMmTy597OWXoXZtOOGEiq9LRCSAtIZGxE8jRsAZZ8Bhh8GQIYRyc+HZZ2HYMLj3XkhN9btCEZFAUKAR8dOpp8KsWXDjjXD33ZwA0LAh/P3vMHiw39VJBfE8ePttmDYNtm2D006DSy+F6tX9rkwkODTlJOK3bt3cvWis5ZupU+Gnn+Daa0svFJaEtHMn9O0LvXvD9u1usG7MGEhPh6VL/a5OJDg0QiNyIAiFIC2N3C1b3IJgOWg8/DB89hl88w00ber2jR8P/frBlVfCe+/5WZ1IcCjQiMg+WbvW3ch4+3Y4/XQwxu+Kgu2f/4Rbby0KM+Ay7QMPQPPmYK36WCQamnISkaiNGweNG7v1yk8+CccdB5dcArm5flcWXMuWQUZG6f1Nm7rpp59+qviaRIJIgUZEovLSSzB6NEyfDt99B1984aZJvvwSrr/e7+qC6+ij3ZRTSYsXw7p1cNRRFV+TSBAp0IhIVCZMgNtug7PPLlqvbAxMmgTPPON++cq+GzTIXdT23XdF+7Ztg6FD3W2KWrTwrzaRINEaGhHZK89zIzJPPVX6WMeOULkyfP01dOpU8bXt1pYt8H//B6+8Qtr69dC9OwwZAs2a+V1ZMX/9K2RlwfHHQ8+ecOih7kr+evVg3jy/qxMJDo3QiMhehUJQp07Z6znWrIGtW93xA8amTe5mLg89BN27s/6cc+Dzz92dl7Oy/K6umKQkl7s+/NBNP1WtCo89BtnZcMQRflcnEhwaoRGRqFx+Odx3n5tyqlrV7fM8t0C4TRto2dLf+ooZOxYKCtwCnxo1WJudTZOxY908ztVXu+GmA+g+P6EQnHyyexOR8lGgEZGojBjhBj3atoUBA6BWLbdQ+LPP3F1uD6B8AP/+Nzz4INSoUbQvFIKRI92wx6JF7s51IpIwFGhEJCp16sDHH7t1NC+95C7V7tLF3UeleXO/qyth7dqyLw+qVw9q1nTHRSShKNCISNRq1YLhw93bAe2EE2Du3NJzOJ9/Dps3H2DzYyISC1oULCKJ58Yb3V0A581zC33APRjp6qvdUx/r1/e1PBGJPY3QiEji6dcPfvzRXQd9zDEck5zs1s2ccw48/rjf1YlIHGiERkQS0x13wPLlcP31bPzjH92VTTNnQvXqflcmInGgERoRSVz168OAAazJzqZxmzZ+VyMicZR4IzTvvQf33196f/iOoVJOH30EY8YUrUcoNHkyvPCCPzUlgs8+g7vuKt2vL7/snicg5ZOdDXfeWbpfZ850l2VJ+eTkuEeDFxQU3z97trsboJTPN9+4lfYl+3XePPfYdYlKYo3Q7NwJAwe6h6Lk5bkhZ3Bh5uqr4fDD4Y9/dJdtSvQKCtwDZxYtgh073C9gcGHmiiugdm3o0QMOOcTXMgPH89x977/4wj2857773P6XX3YLV6tVg1693PetRM/z4NprITPTXdH04INu/8yZcOGF7jkNvXpBw4b+1hlEw4a5Pxp/+w0mTnT7Zs+G88939/k57zxo0sTXEgPpxhthzhz3QLSnn3b75s2Dc8+F/Hzo3VtPKY1CYgWa5GR44w3o3Nn9dQbU2bHDPSK4dm14802FmfKoVMk9XKZzZ7jnHgAOTUmBUaPceoTZsxVmyiMUghkz3M1cwqOKh9Sp44J4lSrue1lhZt+FQvCf/7jv14cfBqB206bwt79BSgq89prCTHlNneq+X598EoBaxx3nRhYKv5cVZspn8mQ444xdo7I127aFm292YeY//1GYiVLIKzkkuxvGmHOBMUAqkAkMtNbmRhyvBIwHzg7vehMYbq0tNoaWlZVVH1jVunVrUlJS9v8rKMsPP7gfZitXuu1DDnG3Mm3bNj7nO1gsXer6tfCBPjVrunt9nHKKn1UF3/Ll7pfE4sVuu3p1F74PqCc9BtAvv7jv18LHWFepAq+/7n5xSPn9+qv7fv36a7ddubIb/ere3d+6gm7dOve9+eWXbjspCaZPd6OJQl5eHgsXLgRo0LZt29VltYlqDY0xpj4wCehlrU0DtgC3l2h2KdAaSA+/7wD0KV/p+6lFC/jzn4u2u3ZVmImFI490U3eF/vAHaN/et3ISRpMmMHhw0Xbr1nDqqf7VkygaNoQbbijaPuYY97NA9k+9esXvrHjEEdCtm3/1JIrDDnPrkwrVqeNuOyBRi3ZRcDcg01q7NLz9FHBZZANr7WSgh7U2H6gD1AbWx6jOffN//1e0HgHg1VfdglbZP5Mnu2mmQnPnlr2gVfbNyy/DLbcUbWdmuh9s6tf9M3OmexhloZwcF3DUr/tn9my3VrHQ0qVuzVLJBa2yb+bNg6uuKtpes8b9Aal+jVq0a2gaASsitlcCjUs2stbmGWNGAcOBz4CPdveCOTk5hOLwNLs6M2bQdPRo8mvU4IcnnmBn9eqkDRxI5Tvv5Oeff2Z15AiDRO3QN96g2ahRFFStyg8TJ7KjTh3SBg4k9Z57+OWXX1g1aNAB9nTCYDhk7lyOvOMOCipXZvGjj5LboAFm0CBS77+f1atX8/OwYerXcqj93ns0v+UWvORkFj/6KFubN8cMHEiVhx/m119/ZeXNN6tfy6HWhx/S/KabIBTix0cfZYsxpA0aRNUnn2TN2rWsuPVWt+ZO9knNTz7hqBtugPx8ljz4IL+np5M2eDDVnnmGdevWsWzkyIO+X6NaHuN53l7f0tLSbk9LS3ssYvvwtLS0bXtoXzktLe3FyM8pfJs/f379+fPnezt27PBiLi/P89q187xDDvG8+fM9z/O8L774wvO+/97zGjXyvBYtPG/TptifN9Hl53tex46eV7Om52Vmep4X7tclSzyvWTPPO/JIz9uwwdcSA6mgwPO6dPG86tU97/33Pc8L9+uyZZ539NGe17ix561Z43ORAVRQ4Hndu3te1aqe9/bbnueF+/Xnnz0vLc3zGjTwvF9+8bnIgOrd2/NSUz3vzTc9zwv36+rVnteypefVret5y5f7XGBAXXyx56WkeN5rr3meF+7XtWs9r00bz6tTx/2sPcjt2LHDmz9/vjd//vz63m6yR7QjNMuByEUoR1B8xAZjTFsg11r7lbV2hzFmCnArFSk5Gd56C5Ytg8ibaLVo4S41rFJFVzmVR6VK7oqbH3+EjIyi/Uce6foVdJVTeYRC7oqbb7+Fk04q2t+kCbz7rrtEXlc57bvCq5wWLSr+cMqGDd3366ZN0KCBb+UF2osvwoIFxS8EqFfPfb+uXQuNSw3cSzSeew6uu86tSyx02GHwzjtucfuRR/pWWpBEG2jmAuONMc2ttUuAAcDMEm1OAi4xxnQDCoCLgfdiVWjUDj3UvZXUokWFl5JQatcuHmYK6T/a/qlRo3iYKaTLX/dPtWqln7QNLtToku3yq1Kl7Ksa69Vzb1I+lSsXDzOFDjvMvUlUopqUs9auJhxijDHfAg2Au4wxJxljZoebTQK+AL4EFgC/AaNjX7KIiPhm3Tq4/XZo0YLju3WDiy+G+fP9rkok+hvrWWtnAbNK7J4P9AgfLwBuCr+JiEiiWbsWOnaEWrVgxAiWbdjA0QsWuH3Tp7s7hov4JLHuFCwiIvFz331uSv/996FyZTZlZ8P118Oxx7p7KS1e7G4IJ+KDg/s6MBERid6LL7p7+VSuXHz/kCHuDsKffOJPXSIo0IiISLQ2b4b69Uvvr17dLXD//feKr0kkTIFGRESic8op7g7MJWVmuidw6xEz4iOtoRERkejccgucfTa0bAlXXun2zZ/vnp3Xvz/UretreXJw0wiNiIhE54wz4Jln3GXb9evTsmdPN2rTtSs8+KDf1clBTiM0IiISvT/9Cfr0gffeY+WiRRzVrx80auR3VSIKNCIiso+qVIHu3fmtfn2FGTlgaMpJREREAk+BRkRERAJPgUZEREQCT4FGREREAk+LgkVEJCF5Hvz3vzB1KmzZ4p6hefnlULOm35VJPGiERkREEk5+Plx6KZx7LmzfDg0awAMPQKtW8P33flcn8aARGhERSThPPAH/+x8sWgTNm7t948bBFVfAZZfBp5/6W5/EnkZoREQk4fzznzB8eFGYAUhOhgkT4IsvYOFC/2qT+FCgERGRhLNsGWRklN7fsKGbfvrpp4qvSeJLgUZERBLO0UeXPa20bBn8/LM7LolFgUZERBLO4MEwfjx89VXRvtxcGDoUOnVyDwyXxKJFwSIiknD+8heYPx/atIGzz4bDD4fZs+HQQ2HePL+rk3jQCI2IiCScUMhd6fTZZ240pnp1eOghtxi4SRO/q5N40AiNiIgkrBNPdG+S+DRCIyIiIoGnQCMiIiKBp0AjIiIigadAIyIiIoGnQCMiIiKBp0AjIiIigRf1ZdvGmHOBMUAqkAkMtNbmlmgzArgEKAC+B66x1q6NXbkiIiIipUU1QmOMqQ9MAnpZa9OALcDtJdr0BnoB7ay16bhAc39syxUREREpLdopp25AprV2aXj7KeCyEm1+Aq6z1m4Nb88Hmu13hSIiIiJ7Ee2UUyNgRcT2SqBxZANrbXbhx8aY2sAI4PHdvWBOTg6hUCj6SsupoKCA7OzsvTeUfaJ+jQ/1a3yoX+ND/Rof6tfSPM/ba5toA01ZIzn5ZTU0xjQFZgIfWmuf3N0Lpqenk5KSEuXpyy87O5uMjIy4n+dgo36ND/VrfKhf40P9Gh/q19Ly8vJYuHDhHttEO+W0HGgYsX0ExUdsADDGtAc+BV6y1g6O8rVFRERE9ku0gWYu0NEY0zy8PQA3CrOLMaYVMAe41lo7LnYlioiIiOxZVIHGWruacIgxxnwLNADuMsacZIyZHW52B5ACjDTGLAi/vRKXqkVEREQiRH0fGmvtLGBWid3zgR7h45fGsC4RERGRqOlOwSIiIhJ4CjQiIiISeAo0IiIiEngKNCIiIhJ4CjQiIiISeAo0IiIiEngKNCIiIhJ4CjQiIiISeAo0IiIiEngKNCIiIhJ4CjQiIiISeAo0IiIiEngKNCIiIhJ4CjQiIiISeAo0IiIiEngKNCIiIhJ4CjQiIiISeAo0IiIiEngKNCIiIhJ4CjQiIiISeAo0IiIiEngKNCIiIhJ4CjQiIiISeAo0IiIiEngKNCIiIhJ4CjQiIiISeAo0IiIiEngKNCIiIhJ4ydE2NMacC4wBUoFMYKC1NreMdpWB2cBj1toZsSpUREREZHeiGqExxtQHJgG9rLVpwBbg9jLanQR8DHSIZZEiIiIiexLtlFM3INNauzS8/RRwWRntrgXuBD7b/9JEREREohPtlFMjYEXE9kqgcclG1tqrAIwxt+7tBXNycgiFQlGevvwKCgrIzs6O+3kONurX+FC/xof6NT7Ur/Ghfi3N87y9tok20JQ1kpO/T9WUkJ6eTkpKyv68RFSys7PJyMiI+3kONurX+FC/xof6NT7Ur/Ghfi0tLy+PhQsX7rFNtFNOy4GGEdtHUHzERkRERMQ30QaauUBHY0zz8PYAYGZ8ShIRERHZN1EFGmvtasIhxhjzLdAAuMsYc5IxZnY8CxQRERHZm6jvQ2OtnQXMKrF7PtCjjLad968sERERkejpTsEiIiISeAo0IiIiErVt2+DRR6FDB2jdGq65Br7+2u+q9mHKSURERA5uW7ZAt27wyy8waBDUrw8zZkDbtvDaa3DWWf7VpkAjIiIiUXnkEVi7Fr74Ag45xO274gq4917o3x+WLIGkJH9q05STiIiIRGXyZLjhhqIwU+jGG2HjRnj/fX/qAo3QiIiIHBiys2HqVBovWQLnnQd9+kBqqt9VFbN+PTRtWnp/tWpQty6sW1fxNRXSCI2IiIifPA+GDoV27WDRIirl5rohj/R0+Oknv6srJiMD5swpvd9aN93k5xMbFGhERET89Nxz8MILboTm9ddZdtddLsi0agV/+pPf1RVz443w1FPw0ktQUOD2LVsGl18OvXvD0Uf7V5sCjYiIiJ8ef7xoRKZQaio89hh88gnk5PhXWwlnngkTJ8LVV8NRR7lBpebN4fDD4Zln/K1Na2hERET8tHgxtG9fen+jRtCkCfzwQ/Gw47Orr4a+feGNN+D33+GUU6BNG7+rUqARERHxV9OmsGCBG/6I9OuvsHIlNGvmT117ULs2XHqp31UUpyknERERP11zDfzjH/Djj0X78vPh5pvdrXj9XGkbIBqhERER8dOAAfDhh3D88dCnD40AMjMhNxfeeQdCIb8rDASN0IiIiPgpKQmef94tSqlWjeR169zd677+Gozxu7rA0AiNiIiI30Ih6NIFunThp+xs6miaaZ9phEZEREQCT4FGREREAk+BRkRERAJPgUZEREQCT4FGREREAk+BRkRERAJPgUZEREQCT4FGREREAk+BRkRERAJPgUZEREQCT4FGREREAk+BRkRERALPj4dTVgLIy8urkJN5nldh5zqYqF/jQ/0aH+rX+FC/xof6tbSI/tjtQEzI87yKqSYsKyurFbCoQk8qIiIiieD4tm3bflXWAT9GaL4DjgfWAwU+nF9ERESCpRJQB5chylThIzQiIiIisaZFwSIiIhJ4CjQiIiISeH6soakwxpjqwEpgurX2Kr/rCTpjzJHAYiAnvCsZWAcMstZ+41ddQWeMSQJuAi7B/ZGRDLwA/N1aqznhcirj+zUJ2Abca619za+6gq6Mfi201lp7ZsVXlDiMMR7uopn8iN2vWWtH+lRSoCR0oMH9gngbON8Yc4u1do3fBSWA3621JxRuGGNuBh4EuvtXUuA9CdQATrfW/m6MOQR4HagC6AfZ/in5/ZoOzDPGrLPWfuRjXUFXrF8lpjpZazf6XUQQJfqU0yDgJVyoGehzLQnHGBMCDgV+9ruWoDLGNAcuBQZaa38HCP8wuwJ438/aEpG1Ngd4BLje71pEJLYSdoTGGHMikIb7SzcPeNwYM85aq7sV7Z+axpgF4Y/r4kYWOvtXTuCdCHxlrd0UudNauxg3rC+x9yVwmd9FBFzkz4FC06y1Y3ypJrF8YIyJnHJqa63N321r2SVhAw1udGaatTbXGDMbqAxcCLzob1mBV3II/wJgrjHmqMIRBtknHok/Unqg8YAtfhcRcJpyih9NOZVTQv4gNcbUBPoBPY0xS4HvcYFmmI9lJSRr7au4XxAt/a4loD4DWoa/Z3cxxmQYYyb7VFOiO4nSC1pFJOASMtDghpOXAUdYa4+01h6JG9pvZ4xp52tlCcYY0wG3eHW3d2+U3bPWrgCmAv80xtQAMMYcjlvnscTP2hKRMaYtbvT2Yb9rEZHYStQpp4HAQ5GXvFprFxtjpgLXofnz/RE5dx4CdgJ9rLUbfKwp6K4B7gY+M8bk4S4vfh6439eqEkPk92sBbqrpCmvtlz7WlAjKWkMDcIa1dl2FVyOCHn0gIiIiCSBRp5xERETkIKJAIyIiIoGnQCMiIiKBp0AjIiIigadAIyIiIoGnQCMiu2WMudIYMyPO56htjJkXz3NEnOtpY0ybijiXiFQsBRoR8duhQEXd8PJM3P2TRCTB6D40IlKMMeYe3BPA1wFLgVRgKDARaIS78d9Ea+0kY0xnYBzwK3Ak8BvwZ2vtj8aYP4SPVQEaAK9aa68Pf86DwHbcH1WbgdOB+dba9saYLcDjwNlAVeBm3NPHWwELgYustfnGmPOB23A3CN0IDLHWfm2MeRbYBLQGmgHzrLXXhL+uW3F3YL4QODb8+flALnCNtfabmHamiFQYjdCIyC7GmPOAc4E2QCegTvjQ88Aj1tq24f3XhgMLQFtgrLU2HZgB/DO8fxhwm7X2ZFy46G+MOTp8rBVwQfjYX3APO2wfPlYN+N5aezwwGXgS6I97XtgJQEdjTBowEjjLWnsicBfwasSX0ho4C0jHPdOtg7V2JPAzcLG1Nge4D+hnrW0XPkfH/eg6EfFZoj76QETK50zcSMoWAGPMM8BVuBAzwRhT2K4WLlx8A3xhrc0M758EjDPGpOBGVXoaY+4AjgNSgBrhdkuttT/voY7CdTtLcCM368L1LMeFrHSgCfC/iJpqG2MKA9hb1to8IM8Y8wNQt4xzvAK8bYyZBcwBXtpjz4jIAU0jNCISyaP4z4WduCkmgPbW2hOstScApwDPRLQpFAq/RgHwAfBH4CtgBG4Kq3D9yta91LEj4uO8Mo5XAuZE1JMBnAwUPlNsW4mvqdS6GWvt34AzgG+B4cC0vdQkIgcwBRoRiTQH6GuMqWWMScY9yHUT8DluHQ3GmAbAl8BJ4c9pa4xpGf74GuBN3AhOW+Bv1toZQBpwBEXhKNLO3ezfk3eBsyOmsK4A3oni83YCycaYZGPMEgBr7SPAnbhpNhEJKAUaEdnFWjsbt14mC8jEjaqAWyR8hjFmIfAecK+19sPwsV+A+40xi4AuuMW1G3BPC19gjMkChgCfAS3KOO0vQLYxZpExplqUdS4CrgVeDdd0De6p73u7ymEGbmrpFOCW8Od/Ea51UDTnFpEDk65yEpFyC1+x9FB42kdExDcaoREREZHA0wiNiIiIBJ5GaERERCTwFGhEREQk8BRoREREJPAUaERERCTwFGhEREQk8BRoREREJPAUaERERCTw/h9KE++gR6jA5QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 680x425 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "tmp = agg_prediction(model_gender.trace['p'], data, ['dept', 'is_male'], summarize())\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 5), dpi=85)\n",
    "color = tmp.is_male.map(lambda is_male: 'b' if is_male else 'r')\n",
    "\n",
    "ax.scatter(data.index, data.ratio, marker='o', facecolors='none', edgecolors=color, label='observed',)\n",
    "ax.scatter(tmp.index, tmp.mean_, marker='x', c=color, label='predicted')\n",
    "\n",
    "ax.set(title='admit ratio', xlabel='departments',\n",
    "       xticks = range(0, 12, 2), xticklabels=tmp.loc[range(0, 12, 2)].dept)\n",
    "ax.legend()\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**department intercepts**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [beta_m, dept]\n",
      "Sampling 2 chains: 100%|██████████| 8000/8000 [00:10<00:00, 766.12draws/s]\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as model_dept:\n",
    "    alpha = pm.Normal('dept', mu=0, sd=1.5, shape=data.dept_i.nunique())\n",
    "\n",
    "    beta_m = pm.Normal('beta_m', mu=0, sd=1)\n",
    "\n",
    "    _alpha = alpha[data.dept_i.values]\n",
    "    _logit_p = _alpha + beta_m * data.is_male.values\n",
    "    \n",
    "    p = pm.Deterministic('p', pm.math.invlogit(_logit_p))\n",
    "\n",
    "    admit = pm.Binomial('admit', n=data.applications.values, p=p,\n",
    "                        observed=data.admit.values)\n",
    "    \n",
    "    model_dept.name = 'department'\n",
    "    model_dept.trace = pm.sample(2000, tune=2000, chains=2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjQAAAGACAYAAAC6OPj9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAANEgAADRIBtYA3dAAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8VcX9//FXSAKiIkoVUXFXRsWoiKg/961u4FKtW6vV1qqtdd/qjnWrVupW3KqtddcvVqS0alG/+tVqXMAouI1WwQUFFUEEBQO5vz/mhoQkwBVyuTnh9Xw88uCecyf3fDL6IG9m5swpy+VySJIkZVmHUhcgSZK0qAw0kiQp8ww0kiQp8ww0kiQp8ww0kiQp8ww0kiQp8ww0kiQp8ww0khZJCOGKEMLTC/m9fwsh3J9/XRZCOC6EUNmKtW0aQtgx/3qnEEIuhLBUa32+pLbDQCOplE4Gjsu/3gG4GShvxc9/GNgw//p5YJUY44xW/HxJbURFqQuQtOSKMX7V6LCsCJeY85kxxu+ACUW4hqQ2oMxHH0hqLISwFXAl0I80ivsi8KsY49v59/8f8CegN/AMMA4IMcadQghHAb8EhgFn5z/yUmAUafRlDeDfwE9jjDNCCH8Dlsq3HduojJ1jjE83qetp4E3gh8CyQF+gC3A1sF3+c0YDJ8cYn8+33zH/7XcAfwOeAjrnr70KMAjYgzQq9DBwSpOQJSkjnHKSNEcIoQvwKPAfYGPSNNCywB/y76+Yf/85YLP86182+Zh++fe2Bq4BrgKuzbfbD9gdOLLJ93wEHJh/vRZpeqglPweOBvYFPgWGk0Zdtsh/fUEKTgAHAB8DvyVNbTX+OSuBJ4HuwG7AAFJAun0e15XUxhloJDW2DHA5cEGMcWyM8WXSyEbv/PuHAFOB02JyLfCvJp9RCfwmxvgucANp9GNwjPH5GOOTpFGdjRp/Q4xxNvBl/nBifnqoJf+OMT6Tr2tp4FbSqMq7McYx+ev1zn/ml8BsYGoLoy57AusAP4kxvhpjfI4Usn4UQuhVQD9JamMMNJLmiDFOII1SnBxCuD2EUA38kYaFuhsBY/IBpN7LTT7myxjjlPzrb/N/jmv0/rdAp4Uscc7nxBinAzcCPwkh/DmE8AxwP4X9vbYhMDbG+Hmjz6vJ17bRPL9LUpvlomBJc4QQVgNGAq8BI4A7gT7ASY2aNV2823Q0ZVYLH13XSiXOuUMphLAsKUxNIa1/eQjoSgo1BX9Oo88rI/1srXmXlaTFxBEaSY39CPgmxrhnjPHqGONTpIW89SHmdWCzJnvF9Gmla3/fOxR2Iq232THGeGWM8TFgFZgTTub3mW8Da4cQujc614e0sDh+zzoktQGO0EhqbBKwSghhD9Iv9v7Ar4HJ+ffvBy4CbgwhDCLdRfQjoLoVrj0t/+fmIYRXCtgvZhIpgPw4hPAcsC0wMP9eJ9IozDQghBC6NfneJ0h3TN0TQjgT6AzcBDwdY3x90X8USYubIzSSGvsf0iLg+4BXgYNIgaZ7CGG1GONk0l1KG+Xf/xkwuJWuPYZ0S/dTwN4LahxjrCYFmKuBN4BT8rXOpmHUaDDp7qrbmnxvHemOq+mkO7b+SZq+2r8Vfg5JJeA+NJIkKfMcoZEkSZlnoJEkSZlnoJEkSZlnoJEkSZlnoJEkSZm32PehGTVqVCXQi/TcltbaPVSSJLVfHYBuwDt9+/atbalBKTbW60XabVSSJOn72Ji071QzpQg0XwJsuOGGVFZWLqjtIhszZgxVVVVFv86Sxn4tDvu1OOzX4rBfi8N+ba62tpa33noL8hmiJaUINHUAlZWViyXQlJWVLZbrLGns1+KwX4vDfi0O+7U47Nf5mudSFRcFS5KkzDPQSJKkzDPQSJKkzDPQSJKkzDPQSJKkzDPQSJKkzDPQSJKkzDPQSJKkzDPQSJJULHV18Mgj8JvfpK9//Qtmz14slz7iiCN44oknFsu1vq9i1GagkSSpGKZPh113hUMPha+/hmnT4Cc/gV12ScdqVaV49IGamDYNnnoKZs2C7beHFVcsdUWSpEV2zjnw5Zfw/vsNf7F/+WUKOb/9Ldx4Y6td6u677+a+++6jQ4cOrLnmmgwcOBCARx99lBtvvJGZM2dy0EEHcdRRRzF9+nTOOussxo8fT1lZGTvttBMnn3wyuVyOG264gSeffJK6ujrWX399Bg4cSJcuXdhll13YZJNNePvttznqqKOorq7muuuuA6C6uppBgwbx97//nccff5xbbrmF2bNn06VLFy688ELWW289Jk6cyNlnn83nn39Oz549mTJlSqv97PUcoSmx66+HVVeFo4+GE06Anj3T/+d183xahSSpzfv2W7j9drjyyrn/ldqtG/zhD3DnnfDNN61yqeeff567776be++9l+HDh7Phhhty3nnnATBjxgyGDBnCvffey1133UVNTQ2PP/44nTp14uGHH+bee+/lgw8+YNq0aQwdOpSPPvqIBx98kGHDhrHWWmtx1VVXzbnOxhtvzGOPPcaAAQN44YUXmDp1KgAPP/wwBx54IGPHjuWGG27g9ttvZ+jQoZx44omccMIJAFx88cX069ePf/7zn5x55pmMGzeuVX72xhyhKaF77kkB/m9/gwMOgA4d0kjNYYdBly5w/vmlrlCStFC++CINv/fr1/y9fv3SdNRnn8Faay3ypf7zn//Qv39/unbtCqT1KTfddBNVVVUccsghlJeX07VrV3bddVeee+459t13X66++mqOPvpottlmG0477TSWXXZZnn76aV5//XUOOOAAAGbNmkW3bt3mXGfzzTcHYNlll2WnnXbi0UcfZZ999uGZZ57hvPPOY/jw4Xz66accfvjhc75n2rRpTJkyheeff35OyFp33XXp11K/LCIDTYnkcnDFFTBwIBx0UMP5XXaBG26A446DM86ApZYqXY2SpIW04oqwzDIwahTsvvvc740aBZ07w0ortcql6poM6dfV1TE7v/C4vLx8zvlcLkdlZSVrrLEGI0aMoLq6mueff56DDz6Ym266iVwux69//WsOyv9S+uabb5gxY8ac7+/cufOc1wceeCDXXnstSy21FNtssw3LLbccdXV1bL/99gwaNGjO9SZMmEDXrl0pKyubq87GdbUWp5xK5Jtv4PXXYf/9m7+3335pmvWddxZ/XZKkVtC5M/zsZ2kYfvLkhvNTpsDZZ8MRR6TA0wq23XZb/vWvf/HVV18BcM8997DllltSUVHBQw89RC6X48svv+SJJ55gu+2245577uF3v/sdO++8M+eeey7rrbceY8eOZZtttmHIkCF888035HI5LrnkEv74xz+2eM1+/foxadIk7rzzTn784x8DsNVWW/Hss8/y4YcfAjB06FCOOuooALbffnsefPBBAD766CNGjhzZKj97Y47QlEjHjlBZCRMnQq9ec7/32Wfpzy5dFn9dkqRWcuWVsOeesO66cOCBUFYGDz4IIaR1NK1k++2359BDD+WnP/0ps2fPZvXVV+fKK6/kzDPPZJllluGAAw5g5syZHH/88fTu3Zs111yT6upq9t57bzp16sQGG2zA3nvvTUVFBRMmTOCggw6asyj4sssua/GaZWVlHHDAATzwwANsvfXWAPTq1YsLL7yQE044gVwuxzLLLMP1119PWVkZAwcO5JxzzqF///6sssoqhBBa7eefU1Mul2v1D52fUaNGrQxM2GSTTaisrCz69WpqaujTp0/Rr7MwDjssTaM+/HBaP1PvpJPgpZfghRdKV9uCtOV+zTL7tTjs1+KwXwswaxb885/pC6B/f9hnH6iY93iC/dpcbW0to0ePBujRt2/fiS21cYSmhC6/HLbZJt3Bd9xx0KkT3HUXjBgB//u/pa5OkrTIKirS2oKW1heoVbmGpoTWXhtefhk23TRNqZ54Iiy/fDq35Zalrk6SpOxwhKbEevaEa69NX5IkaeE4QiNJkjKv4BGaEMK+wGVAJ6AaOC7GOKNJm7OAo4A6YFRLbSRJklpbQSM0IYSVgVuBfWKMvYDpwLlN2vQBfgNsCVQBXYDjW7VaSZKkFhQ65bQ7UB1jHJc/vgU4vEmbDkBHoDNp5Kcz8G0r1ChJkjRfhU45rQZ83Oh4PNCzcYMY46gQwv3Ah6QRnPeA2+f1gWPGjKGsrOz7VbsQ6urqqKmpKfp1ljT2a3HYr8VhvxaH/Tp/FZMnUz55MjPXWWeu853GjmV2167MavScpMbs1+YK2TOv0EDT0kjO7MYHIYT+wI7AGsBk4HrgZtKammaqqqqW+I31ssx+LQ77tTjs1+KwX+dj5sy0/8ann6anDvfunc6/+SYcf3x6jtPLL7f4wL5S9+uf/vQnpk6dynnnnccxxxzDOeecwzpNQlm9r7/+mpNOOonbb5/n+EWLzj77bDbYYIM5j0ZYkEYb681ToVNOHwGrNDpelblHbAD2BobGGD+PMc4CbgR2K/DzJUlqPzp1gqOPhs8/h513hjfeSGFml13SM29+8YtMPH341ltvnWeYAfjqq68YM2bMYqxo3godoRkBXBVCWDvGOBY4BhjWpM0o4NgQwqAY47fAj4CXWq9USZIy5KST0p8nnwwbb9xw/uqr4dRTW+0yL774IoMGDeIHP/gB48ePZ9lll+UPf/gD5557Lssttxxjx47lmGOOoXfv3lx22WVMnTqVsrIyjj/+eHbbbTdqa2u59NJLqa6uplu3bnTt2pU11lgDgF122YUbbriBDTfckAceeIA77riDDh06sNpqq3HllVdy/vnnM336dA466CCGDBnCyJEjGTRoEDNnzqRTp06cccYZbLHFFkybNo3zzjuPGCPdu3dn9uzZbLDBBq3WB1BgoIkxTgwhHAMMCyF0BMYAPw8hbAFcHGPcm7ReZm3g1RDCTOC/wK9btVpJkrLkpJPgrbfg5pvT8dFHt2qYqffGG29wzz330KdPH/7yl79wwQUXANC9e3duuOEGZs2axf77788111zD+uuvzxdffMHBBx/MhhtuyBNPPMH48eN55JFHmDFjBocddticQFPvrbfe4vrrr2fo0KF0796dq6++mnvuuYdLL72U/fffnyFDhjBlyhQuuOAC/va3v7HyyiszduxYjjzySB577DGuv/56unTpwmOPPcbnn3/O/kV4FETB+9DEGIcDw5ucHkmaaiLGmAMuyH9JkqQ334ShQxuO//GPNP1Uv6amlWy00UZz1t0cfPDBDBo0iE033ZTNN98cgLFjx/LBBx9wxhlnzPmeXC7HO++8Q3V1NQMGDKCiooJll12W/v37M2nSpLk+/4UXXmCHHXage/fuAJx22mkAfPxxw+qTV199lYkTJ3LsscfO9b0ffvgh1dXVXHjhhQCstNJK7LLLLq3684OPPpAkqTgar5m5+mooL0/TTzvvPPdC4VZQ0ejp3blcjrKyMjp06EDnzp2BdOfUCiuswLBhDatFJk6cSLdu3XjggQeoq6ubc768vLzZ5zc9N23aNCZPnjzX3cp1dXWEELjvvvvmnJswYQIrrbTSnPfnd41F5aMPJElqbTNnwl57NYSZU09N00/XXZcWCu+5J8xovY30X3/9df773/8C8MADD7DDDjvMFRrWXnttysvLeeyxxwB477332GuvvZg4cSLbb789//jHP/juu+/49ttv57RpbMstt+S5556bM3Jz2223cdttt1FRUcHs2emm50033ZT33nuPV199FUijOvvssw8zZsxg++2356GHHiKXyzF58mSefvrpVvvZ6zlCI0lSa+vUCW66Cd59N43K1DvpJCgrg7XXbtW7nLp3785VV13F+PHj6dGjB5dffjmnn376nPc7duzI4MGDueyyy7jhhhvI5XJcfvnl9OzZk0MPPZQPP/yQAQMGsPzyy7P66qs3+/wNNtiAU045Zc5t1muuuSa///3vWXrppdloo40YMGAAQ4YM4brrruOSSy6htraW8vJyBg8ezDLLLMOJJ57IwIED2XPPPVlppZVYf/31W+1nr1dWyGY1rWnUqFErAxM22WQT96HJMPu1OOzX4rBfi8N+LY7v268vvvgil19++VzTSe1No31oevTt23diS22ccpIkSZlnoJEkKcO22mqrdj06UygDjSRJyjwDjSRJyjwDjSRJyjwDjSRJyjwDjSRJyjwDjSRJyjwDjSRJyjwDjSRJyjwDjSRJyjwDjSRJyjwDjSRJyjwDjSRJyjwDjSRJyjwDjSRJyjwDjSRJyjwDjSRJyjwDjSRJyjwDjSRJyjwDjSRJyjwDjSRJyjwDjSRJyjwDjSRJyryKQhuGEPYFLgM6AdXAcTHGGY3ePxP4aaNvWQX4LMZY1Uq1SpIktaigEZoQwsrArcA+McZewHTg3MZtYoxXxRg3izFuBuyeb3NU65YrSZLUXKFTTrsD1THGcfnjW4DD59P+WuC2GOOoRahNkiSpIIVOOa0GfNzoeDzQs6WGIYTNgO2Bn8/vA8eMGUNZWVmBl194dXV11NTUFP06Sxr7tTjs1+KwX4vDfi0O+7W5XC63wDaFBpqWRnJmz6PticB1McaZ8/vAqqoqKisrC7z8wqupqaFPnz5Fv86Sxn4tDvu1OOzX4rBfi8N+ba62tpbRo0fPt02hU04fkRb51luVuUdsAAghlAP7AXcX+LmSJEmLrNBAMwLYNoSwdv74GGBYC+2qSHc2TWiN4qSF9d13cM89cOCBsM8+cO218NVXpa5KklQsBQWaGONE8iEmhPA20AO4KISwRQjhkUZN1wc+bP0ypcJNmwa77AKnnw5rrw19+sDNN8Omm8K4caWuTpJUDAXvQxNjHA4Mb3J6JLB3ozZDgCGtU5q0cC65BKZOhbfeghVWSOcuvBAOOQSOOw7+/e/S1idJan0FBxopC3I5+Mtf4KabGsIMQEUFXHEF9OoFH34Ia6xRuholSa3PRx+oXamthUmToHfv5u+ttx507AgTXOElSe2OgUbtSseOad3Mf/7T/L2XX04jOGuv3fw9SVK2tb9Ac++98OtfQ11dw7lcDgYOhGuuKV1dWmxOPDH9537zzYZzn38OJ5wAhx4KK61UutokScXRvtbQzJoFV10Fr74Ks2enW1tyObjoIrj4YlhzTfjlL6FLl1JXqiI6+eQUZqqqYOedYZllYMQI2GorGDy41NVJkoqhfQWaior0m2vXXeHWWwFYJZeD225LYebppw0zS4AOHdJ//pNOgqFD07qaU06BnXaCxfC0DUlSCbSvQANpPuHJJ+eEmh4Aq6+ewsxaa5W2tqz77ru0SKWxWbNSSigvL01N81FVlb4kSe1f+1tDA7DiirDddg3HG2zgfbqL6uKL0/zN1183nJs1C448Mn3NntejvSRJKr72F2jq18zcdBMsvTSzl1kGHn8cfvWruRcKq3B1dTB6NDz/POy5Zwo19WHm3nvhnXdg+vRSVylJWoK1rymnpguAn36aN99/n6pTTpmzpoabb06LLFS4Dh3gvvvgsMPg73+HH/6QdTp1gmeegX790rql5ZYrdZWSpCVY+/rNPns2jBzZsAB4rbWYtcIKaU1NVRWMGQPfflvqKrOpsjKFmv33hxdfpOszz6SpvBEjYPnlS12dJGkJ175GaCoq4KGH0qYjPXs2nK9fKNyxY7qHt4347jv47LO5S4W0021FBXTtWpq65qmsLAWbeuXlbXIxsCRpydO+RmgAOnVqnhAghZo2lBBmz4aDD4ZttoH33284P2kS7LYb7L57emp0m1G/ZmbIENhsM6ZutRW88UbDmhpJkkqo/QWajCgvh802g48+SvujvP9+Q5h59VXYaCPo3LnUVebV1TUsAO7XD556iveuvx4OPLBhoXCbSl+SpCVN+5pyypiLLkp//u53sMkmKeRMnQpHHZX2AmwzszllZWkPn/oFwMsv37Cm5rDD0qLhTp1KXaUkaQnmCE2JXXRRevbQ9OkpzOywQxsLM5ACzaWXwv/9X9ppedQolh49Os2b3Xcf3HPP3GtrSuzbb+Gtt5qfnzABPvlk8dcjSSo+A02JTZoEzz7bcPzee/DBB6WrZ57KyuCf/4R114Wtt2b9Y49Na5X+9Ke0grmNqKuDH/0Itt0Wamoazk+YALvskqb3Jk8uWXmSpCIx0JRQ4zUzRx0FF1wA48c3rKlpUx59FA4/HE49FaZM4bVnn01Perz4YrjuulJXN0eHDmlpz+TJ6ekXNTUNYeatt2CPPbzLXJLao7bzT+slzOzZaS1tfZipn2bq0CGtqdlpp7RtTpu5MeuSS+D009OjrCFNMR16aLr76fTT4fjjmz/nqUSOOSb9eeyxsPnmDedPOAGuv94HVEpSe9Q+A01dXdp35h//oOeECfDTn8KAAW1qaqS8HM44Iz2V4ZZbGtbMXHRR+oXbuXMbCjMzZ0J1dSq0qYMOgp/9DN58M9221UYccwxMnJhGvSBNQRlmJKn9an9TTjNmQP/+aSHF5Ml0mDkTfvGL9BttypRSVzeXQw5peQHwwIFw1lmlqalFHTqkIlvab2b69PTIiaWWWvx1zceECeku83pvvplGwyRJ7VP7CzS/+x2MHQv//S/cfTcf1h+XlzdMl+j7qayE/fZLC4BzubnfGzw4bZoTQmlqa0HjNTMnnJAGlhqvqZEktT9tZw6mNcyenR5CefPN0KNHw/muXeGPf0wLU669FlZYoWQlZtYll6RRroMPhuOPZ+kPP4S//hX+/GcYNqzNzOXU1aXZxfowUz/NVFaW1tTsthu8+y5061bqSiVJral9BZqpU9OtQ/36NX+vX7/08KTx4w00C2OjjdI6mt/9Dvbai16zZqUhjyefhO22K3V1c3ToAFdemW7KuuqqhpxVv1C4ttYwI0ntUfsKNF26pLDyyivpiduNjRqVpk5WXbU0tbUHG2yQNtLL5Xi1poY+jW8hakN23TV9NVUfaiRJ7U/7WkNTUQFHH51ubfn884bz06bBmWemO3L85/miq5/DkSSpjWhfIzSQpkRqamC99eDAA1l96lR45pm0q+3115e6OkmSVATta4QGYOml4d//TvfsVlRQNmtWuhPnhRfgBz8odXWSJKkICh6hCSHsC1wGdAKqgeNijDOatNkFuAJYCvgYOCLGOKn1yi1QeXnai6Z/fz6sqeEHffos9hIkSdLiU9AITQhhZeBWYJ8YYy9gOnBukzarA/8D/DzGuAnwBnBx65YrSZLUXKFTTrsD1THGcfnjW4DDm7T5MfDPGOMb+eNLgUGLXKEkSdICFDrltBppCqneeKBnkzbrAzNCCA8C65FGaE6Z1weOGTOGssVwp0xdXR01bg/b6uzX4rBfi8N+LQ77tTjs1+ZyTXepb0GhgaalkZzZTY4rgb2A7YAPSGtpbgEOaOkDq6qqqKysLPDyC6+mpoY+rqFpdfZrcdivxWG/Fof9Whz2a3O1tbWMHj16vm0KnXL6CFil0fGqzD1iA/AJ8EyMcVyMMQfcCWxV4OdLkiQttEIDzQhg2xDC2vnjY4BhTdo8DOwYQqifitoPeHnRS5QkSZq/ggJNjHEi+RATQngb6AFcFELYIoTwSL5NDXAa8EgI4U1gR+DXxSlbkiSpQcH70MQYhwPDm5weCezdqM0DwAOtU5okSVJh2t9OwZIkaYljoJEkSZlnoJEkSZlnoJEkSZlnoJEkSZlnoJEkSZlnoJEkSZlnoJEkSZlnoJEkSZlnoJEkSZlnoJEkSZlnoJEkSZlnoJEkSZlnoJEkSZlnoJEkSZlnoJEkSZlnoJEkSZlnoJEkSZlnoJEkSZlnoJEkSZlnoJEkSZlnoJEkSZlnoJEkSZlnoJEkSZlnoJEkSZlnoJEkSZlnoJEkSZlnoJEkSZlXUWjDEMK+wGVAJ6AaOC7GOKNJm8eAtYFv86fuiDFe00q1SpIktaigQBNCWBm4FdgqxjguhHAjcC5wYaM2ZUBfYN0Y49RiFCtJktSSQqecdgeqY4zj8se3AIc3abNh/s8HQghjQgjXhRA6t0KNkiRJ81XolNNqwMeNjscDPZu0+QHwJPBr4BvgLuBy4NSWPnDMmDGUlZV9r2IXRl1dHTU1NUW/zpLGfi0O+7U47NfisF+Lw35tLpfLLbBNoYGmpZGc2Y0PYozPAs/WH4cQrgAeZB6BpqqqisrKygIvv/Bqamro06dP0a+zpLFfi8N+LQ77tTjs1+KwX5urra1l9OjR821T6JTTR8AqjY5XZe4RG0IIO4cQ9mx0qgyoLfDzJUmSFlqhgWYEsG0IYe388THAsCZtlgeuCSEsG0LoAJwCDGmdMiVJkuatoEATY5xIPsSEEN4GegAXhRC2CCE8km8zFLgPeBl4G5gBXFKUqiVJkhopeB+aGONwYHiT0yOBvRu1uRi4uHVKkyRJKow7BUuSpMwz0EiSpMwz0EiSpMwz0EiSpMwz0EiSpMwz0EiSpMwz0EiSpMwz0EiSpMwz0EiSpMwz0EiSpMwz0EiSpMwz0EiSpMwz0EiSpMwz0EiSpMwz0EiSpMwz0EiSpMwz0EiSpMwz0EiSpMwz0EiSpMwz0EiSpMwz0EiSpMwz0EiSpMwz0EiSpMwz0EiSpMwz0EiSpMwz0EiSpMwz0EiSpMwz0EiSpMyrKLRhCGFf4DKgE1ANHBdjnDGPtoOA9WKM+7dKlZIkSfNR0AhNCGFl4FZgnxhjL2A6cO482u4DHNlqFUqSJC1AoVNOuwPVMcZx+eNbgMObNgohrAlcCFzSKtVJkiQVoNBAsxrwcaPj8UDPxg1CCJXAXcDxwNRWqU6SJKkAha6haSn4zG5yfCUwLMb4cgih94I+cMyYMZSVlRV4+YVXV1dHTU1N0a+zpLFfi8N+LQ77tTjs1+KwX5vL5XILbFNooPkI6NvoeFXmHrEBOAz4LIRwBNAN6BpCeDDG+OOWPrCqqorKysoCL7/wampq6NOnT9Gvs6SxX4vDfi0O+7U47NfisF+bq62tZfTo0fNtU2igGQFcFUJYO8Y4FjgGGNa4QYxxlfrXIYSjgP3nFWYkSZJaU0FraGKME8mHmBDC20AP4KIQwhYhhEeKWaAkSdKCFLwPTYxxODC8yemRwN4ttP0b8LdFKUySJKlQ7hQsSZIyz0AjtRUzZ1L27belrkKSMslAI5XamDGw336w9NJNKl+6AAAXdklEQVRstu22sPXW8K9/lboqScoUA41USm+8AdttB8stB9XVvPXgg/DDH8IBB8D//E+pq5OkzCh4UbCkIrjoIth7b7jzTigrY0ZlJRx4IKy4Ipx1VnpdXl7qKiWpzXOERiqVXA6GDUvTTU13zT7wQPjoI3jttdLUJkkZY6CRSqWuDmpr4fjj09RTvcmTYf/9U+CZNKl09UlShhhopFIpL4cQUoDZZRd44w3Kp06F3XeHUaOgY0fYdttSVylJmWCgkUrptttSsPnsM9h8czbeYw8YORIqKuCKK2DppUtdYTZ98AH89a/Nz7/4Ijzi5uZSe+SiYKmUttsu/YI9/HD4/PP0L4yKCrj6ajjhhFJXl025HBx6KLzwAnz9NZx8cjr/4otp9Ku2Ft57D1ZZZf6fsxjlcvDEEzBkCHz7LeywA/zkJ7DMMqWuTMoOR2ikUuvXD9Zcs+G4a9c0BdV0obAKU1aWRr66d4dTToHrrmPpMWNSmPn6a7jlljYVZmbNgoMOSsumZs6EH/wALrsMqqpg3LhSVydlh4FGKqXJk9Mv2pEj4ac/5ePTTksLgfNrarSQeveG//3fOaEmHHlkCjN33AFHHFHq6uZy3XXw0kvw1lupvGuvhXffhc03h6OOKnV1UnYYaKRSyeVgn33mhBnuuIPPDz88/Ub77LMUar74otRVNvPFF+kX75//DDGWupr56N07rUOqt/nmbS7MQOrHs8+GNdZoOFdZmWYdn322jfex1IYYaKRSKSuDc86BX/wiJYT6DfROPjmFmhNPTBvstSFXXgk9e8Kll8LNN8OGG6blKjNmlLqyFrz4YppyqjdqVBoOaWM+/BD69Gl+fo010vTTBx8s/pqkLDLQSKXUvz/85S/NdwM++WQ4//zS1DQP110HF14IQ4fCO+/AK6+kaZLq6vSkhjalfgHw11/DnXfy1pAhc62paUvWXTdNOTX13ntp9nGddRZ/TVIWGWgkLVAuB+eeC999l6ZA6tcrf/55+nr0UXj77dLWOEcuB6efPteamRnrrtuwpub88+HTT0td5Ry/+hX8/vcpJNb79ts0QLfbbrDeeqWrTcoSb9uWVJBvvoFu3eDUU9PxFlvAnnumO3M6dUrBZoMNSlsjkNLWQw/Bk0+mBT877kivL79Mxd59N3Tu3Kbucvr1r+Hll2HjjdOA3QorwPDhsNJK8O9/l7o6KTscoZG0QGVlaT3HBRekQY5TT4Xtt08jCTfdlLZ26dat1FU2stRSMGhQWou05558OWBASg0HH5wCTRty551pdObRR9P0U+fO8Kc/wV57wRlnpNu6JS2YgUZSQY44Au69d+6lPVtuCa+/DptuChttVLramrn88vSsrNdeg3PO4YuDD4annkp3k/3yl2laqg2o31DvhRfSlN7558PgwWl90tVXp3XMPs5LKoxTTpIKcsEF0Ldvw8a7kH4R19Sk24vb1D6Ad9wB11wDyy7bcK6sLK1qXnXVlMKqqkpXX6OS7rgjZa/770/rmPv0Sbdyr7tuymArr1zqKqVsMNBIKsibb6Z1MpB+2ZaXw0cfpWmn555LGx63GV980fLtQd27Q5cubWp/n4oKuOuu9Pr++9PM2IorpjCz+uqlrU3KEqecJC1QLgdnnZXCy333wX//m+52evnllBEuuijtBdhmbLYZjBjR/PzLL8O0aW1sfiyFw8brlDt2TLlLUuEMNJIWqKwMhg1LX4cc0nC+/gkD//53CjZtxmmnpV0AH3+8Yb3MuHFp/cxPftKm5nFyOfjtb9MM2ZprpueVfvIJ7LEHTJlS6uqk7HDKSVJBVloJBgxofr5378VfywIddhi8/366D3r99Vm/oiKtmxkwAG68sdTVzVEfZq66qmHNzCqrpAXY99+fQs2IEel5pZLmz0AjqX0677w0IvOPfzDl3XdZ9s470+1YbUhZGSyzTEOYqV8zU7+m5pNPmm8iLallBhpJ7dfKK8Mxx/B5TQ0921iYqTdwYNrXZ7nlGs7VLxSeOTMFHkkLZqCRpBJrHGbqVVSkL0mFcVGwJEnKPAONJEnKvIIHNEMI+wKXAZ2AauC4GOOMRu93AK4C9sqfehQ4M8ZY13rlSpIkNVfQCE0IYWXgVmCfGGMvYDpwbpNmPwE2Aaryf24DHNh6pUqSJLWs0Cmn3YHqGOO4/PEtwOGNG8QY7wb2jjHOBroBXYEvW6lOSZKkeSp0ymk14ONGx+OBnk0bxRhrQwgDgTOBl4Dn5vWBY8aMoWwxPM2urq6Ompqaol9nSWO/Fof9Whz2a3HYr8VhvzaXq9/xez4KDTQtjeTMbqlhjPF3IYTfA3cAg4ATWmpXVVVFZWVlgZdfeDU1NfTp06fo11nS2K/FYb8Wh/1aHPZrcdivzdXW1jJ69Oj5til0yukjoNGj01iVuUdsCCH0DSH0BogxfgfcA2xWcLWSJEkLqdBAMwLYNoSwdv74GGBYkzZbAINDCJUhhHLgEODpVqlSkiRpPgoKNDHGieRDTAjhbaAHcFEIYYsQwiP5ZrcCrwCvAa8CXwGXtH7JkqSSmTQJzj0X1luPjXffPT1+feTIUlclFb4PTYxxODC8yemRwN759+uA0/NfkqT25osvYNtt07MaLriADydPZt1XX03nhg6FvfcudYVagvmkEElSYa64AlZYAZ55Bjp2ZGpNDZxyCmywARx/PLz3no8HV8n46ANJUmHuuy89Grxjx7nPn3ACfPYZvPBCaeqSMNBIkgo1bRqsvHLz88ssA8suC19/vfhrkvIMNJKkwmy9NQxreoMrUF0NX30Fffsu/pqkPNfQSJIKc9ZZsNdesNFGcNRR6dzIkfCzn8HRR8NKK5W0PC3ZHKGRJBVmlVXSnUznnAMrr8xG/funUZtevVLIkUrIERpJUmF+8xt4+mn41a9g330Z/8YbrLPqqvCLX8D//i/stx+svnqpq9QSykAjSSrMfffBzjvDzTdDRQV1G2yQwkwuB3//u2FGJWWgkSQVpkcPeOqpFGoGD2Y9gMpKeOghGDCg1NVpCecaGklS4Xr0SGto6vXsCf37l64eKc9AI0kq3IgRcOyxDcdjx8JJJ6VpJ6mEDDSSpMKMGAH77pvCy/DhjBkxIj32YPBgQ41KzjU0kqTCDB7csAB4wABm1dQ0rKm5//60T00bWRj8ySfw8MMQAjzwAEyfnp6huckm6SkNBxxQ6grV2gw0kqTCPPAAvPJKSgb16hcKf/FFmwkzAEccke4kr6yEww5LZf7+9ynolJfDO+/AWmuVukq1JqecJEmF6dx57jBTr0cP2HjjxV/PfGy3HXToALW1ad+/ww+Hb76BujpYbTXDTHtkoJEktTtDh8Jvf5vCy/nnw+abw5dfwh/+AB9/DKNHl7pCtTYDjSSp3fnwQ9h9d7j22oZza64JZ56ZBpQ++KB0tak4XEMjSWp31l037fd3zz0N5z74IAWaTz5J76t9cYRGktTuDBiQbsr68ku45RZ4911YdVUYNAjWWMNnabZHjtBIktqdF15Id5iXlcHw4VBdDTNmQEUFfPopjBvnwuD2xkAjSWp37r8fnnwS1lmnYR+awYNh001h4kTDTHtkoJEktTsrrAA//nF6vfnmc7/ndFP75BoaSZKUeQYaSZKUeQYaSZKUeQYaSZKUeQYaSZKUeQYaSZKUeQXfth1C2Be4DOgEVAPHxRhnNGlzAXAoUAe8CxwbY/yi9cqVJElqrqARmhDCysCtwD4xxl7AdODcJm32B/YB+sUYq0iB5g+tW64kSVJzhU457Q5UxxjH5Y9vAQ5v0uYD4OQY4zf545HAmotcoSRJ0gIUOuW0GvBxo+PxQM/GDWKMNfWvQwhdgQuAG+f1gWPGjKGsrKzwShdSXV0dNTU1C26o78V+LQ77tTjs1+KwX4vDfm0ul8stsE2hgaalkZzZLTUMIawBDAP+E2O8eV4fWFVVRWVlZYGXX3g1NTX06dOn6NdZ0tivxWG/Fof9Whz2a3HYr83V1tYyevTo+bYpdMrpI2CVRserMveIDQAhhK2AF4H7Y4zHF/jZkiRJi6TQQDMC2DaEsHb++BjSKMwcIYTewGPAb2KMV7ZeiZIkSfNXUKCJMU4kH2JCCG8DPYCLQghbhBAeyTc7D6gELgwhvJr/erAoVUuSJDVS8D40McbhwPAmp0cCe+ff/0kr1iVJklQwdwqWJEmZZ6CRJEmZZ6CRJEmZZ6CRJEmZZ6CRJEmZZ6CRJEmZZ6CRJEmZZ6CRJEmZZ6CRJEmZZ6CRJEmZZ6CRJEmZZ6CRJEmZZ6CRJEmZZ6CRJEmZZ6CRJEmZZ6CRJEmZZ6CRJEmZZ6CRJEmZZ6CRJEmZZ6CRJEmZZ6CRJEmZZ6CRJEmZZ6CRJEmZZ6CRJEmZZ6CRJEmZZ6CRJEmZZ6CRJEmZZ6CRJEmZV1FowxDCvsBlQCegGjguxjijhXYdgUeAwTHGh1urUEmSpHkpaIQmhLAycCuwT4yxFzAdOLeFdlsAzwPbtGaRkiRJ81PolNPuQHWMcVz++Bbg8Bba/QY4H3hp0UuTJEkqTKFTTqsBHzc6Hg/0bNooxvhzgBDC2Qv6wDFjxlBWVlbg5RdeXV0dNTU1Rb/OksZ+LQ77tTjs1+KwX4vDfm0ul8stsE2hgaalkZzZ36uaJqqqqqisrFyUjyhITU0Nffr0Kfp1ljT2a3HYr8VhvxaH/Voc9mtztbW1jB49er5tCp1y+ghYpdHxqsw9YiNJklQyhQaaEcC2IYS188fHAMOKU5IkSdL3U1CgiTFOJB9iQghvAz2Ai0IIW4QQHilmgZIkSQtS8D40McbhwPAmp0cCe7fQdqdFK0uSJKlw7hQsSZIyz0AjSVIpzZoFLd2mPWUKvPPO4q9nAb79Fv70J9hmG9hkEzj2WHjzzVJXZaCRJKm0fv7zlA5GjGg4N2UK7LEH7LADfPBB6WprYvp02G03uOYa2H9/OP10+Pxz6NsXHn+8tLUVvIZGkiQVwSGHwAMPwH77wbBhlC+1VAozL70Ehx0Gq61W6grnuP76FGBeeQWWXz6dO/JIuPRS+MUvYNw4KC8vTW2O0EiSVEoDBsBDD8Hs2bDHHmyy444NYebOO6Gi7Yw93HgjfPddmiVrrK4OJkyAJ58sTV1goJEkqfQGDICBAxuOu3SBW29tU2EGYNKkNAO2227wxRfp3BVXpNJzOXjvvdLV1rZ6SpKkJU0uB8cdlwJMva+/hl694PnnYc01S1dbEzvumALNa6+lUPPDH8KgQbDSSmkqavfdS1ebIzSSJJXSTTfBX/+aXh92GO//8Y9QWZnmcPZuttVbSZ1+Ovz3vynYvPZaCjMdO0KPHvCjH8G665auNgONJEmldO65af1Mfs3MVzvvnNbUlJen+6EffbTUFc6x225pHc3zzzec++476N4dbr+9dHWBgUaSpNLbZ5+5FwDXLxRefnmYMaO0tTXxxRdQWwtdu8J66819rpQMNJIkldI666T9ZpouAN5yy7TxSxtaQ3PFFXDOObDKKulGrBjTNjr1a2rqFwqXgoFGkqRSOvZY+OMf4f33G87Nng1nnJG24u3Tp3S1NVFRAauuCk8/ndYsd+gAt92WQk2HDumrZLWV7tKSJIljjoH//Ac23hgOPJDVAKqr01TTk09CWVmpK5zjjDPgl79s2FQPGkLNtGmw3HKlq80RGkmSSqm8HO66C/71L1h6aSomTYJTT00LgkModXXNNA4z9Tp0KG2YAUdoJEkqvbIy2Hln2HlnPqipoVsbmmbKCkdoJElS5hloJElS5hloJElS5hloJElS5hloJElS5hloJElS5hloJElS5hloJElS5hloJElS5hloJElS5hloJElS5hloJElS5pXi4ZQdAGpraxfLxXK53GK71pLEfi0O+7U47NfisF+Lw35trlF/zHMgpiyXyy2eavJGjRrVG3h9sV5UkiS1Bxv37dv3jZbeKMUIzTvAxsCXQF0Jri9JkrKlA9CNlCFatNhHaCRJklqbi4IlSVLmGWgkSVLmlWINzWITQlgGGA8MjTH+vNT1ZF0IYS3gPWBM/lQFMAn4VYzxrVLVlXUhhHLgdOBQ0j8yKoB7gd/HGJ0TXkgt/P9aDnwLXBpj/Eep6sq6Fvq13hcxxt0Wf0XtRwghR7ppZnaj0/+IMV5YopIypV0HGtIviCeAH4UQzooxfl7qgtqBr2OMm9UfhBDOAK4B9ixdSZl3M7AssGOM8esQwvLAP4GlAP8iWzRN/3+tAh4PIUyKMT5Xwrqybq5+VavaPsY4pdRFZFF7n3L6FXA/KdQcV+Ja2p0QQhmwAvBJqWvJqhDC2sBPgONijF8D5P8yOxJ4ppS1tUcxxjHA9cAppa5FUutqtyM0IYTNgV6kf+nWAjeGEK6MMbpb0aLpEkJ4Nf96JdLIwk6lKyfzNgfeiDFObXwyxvgeaVhfre814PBSF5Fxjf8eqDckxnhZSappX54NITSecuobY5w9z9aao90GGtLozJAY44wQwiNAR+DHwH2lLSvzmg7hHwCMCCGsUz/CoO8lR/sfKW1rcsD0UheRcU45FY9TTgupXf5FGkLoAhwG9A8hjAPeJQWak0pYVrsUY3yI9Atio1LXklEvARvl/5+dI4TQJ4Rwd4lqau+2oPmCVkkZ1y4DDWk4+UNg1RjjWjHGtUhD+/1CCP1KWlk7E0LYhrR4dZ67N2reYowfAw8Afw4hLAsQQliRtM5jbClra49CCH1Jo7fXlboWSa2rvU45HQdc2/iW1xjjeyGEB4CTcf58UTSeOy8DZgEHxhgnl7CmrDsW+B3wUgihlnR78V3AH0paVfvQ+P/XOtJU05ExxtdKWFN70NIaGoBdY4yTFns1Ej76QJIktQPtdcpJkiQtQQw0kiQp8ww0kiQp8ww0kiQp8ww0kiQp8ww0kuYphHBUCOHhIl+jawjh8WJeo9G1bgshbLo4riVp8TLQSCq1FYDFteHlbqT9kyS1M+5DI2kuIYSLSU8AnwSMAzoBJwI3AKuRNv67IcZ4awhhJ+BK4DNgLeAr4GcxxvdDCP8v/95SQA/goRjjKfnvuQaYSfpH1TRgR2BkjHGrEMJ04EZgL6AzcAbp6eO9gdHAwTHG2SGEHwHnkDYInQKcEGN8M4TwN2AqsAmwJvB4jPHY/M91NmkH5h8DG+S/fzYwAzg2xvhWq3ampMXGERpJc4QQ9gP2BTYFtge65d+6C7g+xtg3f/43+cAC0Be4PMZYBTwM/Dl//iTgnBjjlqRwcXQIYd38e72BA/Lv/YL0sMOt8u8tDbwbY9wYuBu4GTia9LywzYBtQwi9gAuBH8YYNwcuAh5q9KNsAvwQqCI9022bGOOFwCfAITHGMcAVwGExxn75a2y7CF0nqcTa66MPJC2c3UgjKdMBQgi3Az8nhZhBIYT6dsuRwsVbwCsxxur8+VuBK0MIlaRRlf4hhPOADYFKYNl8u3Exxk/mU0f9up2xpJGbSfl6PiKFrCpgdeD/GtXUNYRQH8D+HWOsBWpDCP8FVmrhGg8CT4QQhgOPAffPt2cktWmO0EhqLMfcfy/MIk0xAWwVY9wsxrgZsDVwe6M29cryn1EHPAvsAbwBXECawqpfv/LNAur4rtHr2hbe7wA81qiePsCWQP0zxb5t8jM1WzcTY/wtsCvwNnAmMGQBNUlqwww0khp7DDgohLBcCKGC9CDXqcDLpHU0hBB6AK8BW+S/p28IYaP862OBR0kjOH2B38YYHwZ6AavSEI4amzWP8/PzFLBXoymsI4EnC/i+WUBFCKEihDAWIMZ4PXA+aZpNUkYZaCTNEWN8hLReZhRQTRpVgbRIeNcQwmjgaeDSGON/8u99CvwhhPA6sDNpce1k0tPCXw0hjAJOAF4C1mvhsp8CNSGE10MISxdY5+vAb4CH8jUdS3rq+4LucniYNLW0NXBW/vtfydf6q0KuLalt8i4nSQstf8fStflpH0kqGUdoJElS5jlCI0mSMs8RGkmSlHkGGkmSlHkGGkmSlHkGGkmSlHkGGkmSlHkGGkmSlHkGGkmSlHn/Hyr446A68taRAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 680x425 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "tmp = agg_prediction(model_dept.trace['p'], data, ['dept', 'is_male'], summarize())\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 5), dpi=85)\n",
    "color = tmp.is_male.map(lambda is_male: 'b' if is_male else 'r')\n",
    "\n",
    "ax.scatter(data.index, data.ratio, marker='o', facecolors='none', edgecolors=color, label='observed',)\n",
    "ax.scatter(tmp.index, tmp.mean_, marker='x', c=color, label='predicted')\n",
    "\n",
    "ax.set(title='admit ratio', xlabel='departments',\n",
    "       xticks = range(0, 12, 2), xticklabels=tmp.loc[range(0, 12, 2)].dept)\n",
    "ax.legend()\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEcCAYAAAALEfkWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xt4VYWZ7/FvoFTAWC+YyswcMIXKa3roqT5MPNYqlxaoMsW06FQOBwW826pFvExqrbZClRksorUCcgvCeLDiBRkvI2qaqsjAONpqJ74gIRBBIFJAwg6YQs4fe0W3MXd2stZe+X2eh4edtdZe68fKZr9517v3TlZtbS0iIiLp1CXsACIiEj8qLiIiknYqLiIiknYqLiIiknYqLiIiknYqLiIiknZfCDuASGuZ2VHAIqAf8BHwY3ffYGY/AO4BKoJN7wDeAJ4GegBXufufzOxs4Fvu/s+N7P9KYDxwGOgG/Mzdf9+GnFcDvYE5wO3u/iMzKwdOdfcDKdtNDJYVtmCfRcAyd3++kfV9gW+4+8rW5m3muEOBq919bDr3K/GlzkUy0RVAlbufCVwHPBAsHwTc4u5Dgz8lwEiSxeVHwGVmlgX8BLivoR2b2VhgBPAddx9KssgsMbMT2xrW3be7+4/aev9W+jbwrQ46lkij1LlIJvoa8ByAu7uZ5QXLBwGnm9lkYC3wT0AVya6lB7AfGAc8mdo51HMVMMXda4L9bzKz09x9V9A19Ar+/ANwC3AO0BWY6e6PBV3RfcBu4K/AGjPLJdltnBkcY26wbAcwIfXgZnZdkLE2uM/9DYUMOol/Aj4m2cEtA6YDhUBPM1sNbALuB7KAXcClwOnAPwf3exEY6u7Dgn3+G/BzoD/wY5JdWy3wg3rHXgR8NTin97n7kkbOpXRi6lwkE70FfM/MsszsTODvzKwrsIpkJzMYyAauJvkEehLJovEQySfKP5rZXDO7pYF9/y1QlrrA3XelfPmyu58FnAl8xd3PBoYBPzOz44DZwP9x9+Ekn9wbMtvdhwDlJLswAMzsa8BFwNkki9b3zcyaOA8nAxcEWW5x90MkC8wj7v40MI/kJcOhwLMkiyFAd3c/x91/CXQ3s5PN7G+AE939TWAA8A/Bv+2/ge+mZDyG5PkdA5wLHGoin3Ri6lwkEy0E8oBXgNeAN9z9kJktdPc9AGa2ArjA3Q8Dk4NlPyXZVdxGsgj9wswGuPv6lH1vBvoAe+sWmNl3gT8FX3rw99eBQWb2++DrbkAucFLK/l4j+RN+qo/dfU1wezXJS3Drgq8HkiwYLwVfHw+cknLM+t52978CfzWz6gbW5wEPBvWpG7Ch3r8BYAFwCXCQ5BwLYCew2MyqgFOB1+s2dvd9QWf4EPAlYGkj2aSTU+cimSgfeCn4yfoxoCyYpfzJzP5HsM13SA7zATCzLwPm7q8APUn+xF0LHF1v3wuBn5vZF4L7DQDm8+lP6IeDv98FioOu4NvA74CNwNaUy3T5DWT/opmdFtw+B3gnZZ0DfwaGBfst4tOi1pCGPhjwMJ/+v3bgkmBftwD/Vu/fAMnLad8j2dE9YmbHAr8ExgKXA9UkL6sBEHQ4g9z9ByQvDf5L3bkSSaXiIploAzDZzF4HppKckdSSfDJ8wsxKSBaQeSn3uQ2YFtx+EPh34G+AP6bu2N2XAWuAV83sDyR/mh/v7jvrZVgJVJnZKySLWK277yN5+e1hM3uJZBdS30HgumDfXyZZuOqO/UeSXcurZvafJLuWrS0/LQC8DRQEL0y4JsjyKsnLZZ8rVO5eFZyD0iD/RyQ7rtdJdobVJC8V1tkO9A5mOquAe4LuSeQzsvSpyCIikm7qXEREJO1UXEREJO1UXEREJO1UXEREJO1UXEREJO0y+vXpb7zxhl7qJiISkkGDBmU1ti7TX4rcovClpaXk5eU1v2HIlDO9lDO9MiHnvHnz+OCDD7j99tvDjtKsTDif0GTORgsLZHjnIiKS6oorrqC0tDTsGIJmLiIi0g7UuYhIbAwdOpREIsHatWvDjtLpqbiISGxMnDiRbdu2hR1DUHERkRiZOHGiZi4RoZmLiMRGTU0NNTU1YccQ1LmISIyMGDFCM5eIUHERkdi4/PLL2bq1tb8CR9qDiouIxMb48eM1c4kIFRcRiY1EIkF1dXXYMdrNll0JLlu8jrLK/fTLOZoFE/Lp26tn2LEapOIiIrExatSoWM1cJi1aS7FXNrhuw84qBs8o/syyYZbDoklndES0Zqm4iEhsXHPNNRkxcxl5bwnrd1QBZWndb7FXklv4TIu2HXBSNi/cMCStx0/VpuJiZucCY919Yivv9wPgP9x9W8qyicCdwCzgYeARoAewDZgEXArcBEx39zltySsincNFF12UETOXF24Y0qYPrhwxs4SNlVUcroUuWdA/J5tVU9qvQByJjn6fy0+ALzWw/BF3nwncHtw+B3gTuMrdHwCKOi6iiGSqvXv3sm/fvrBjtJsFE/Lpe0JyxnK4FmoOHWbLrkTIqRrW4s7FzPKAhcD+4M/uYPk/AlOAQ8Cr7l5oZr8ATgW+DBwPXAccC5wGPGxmZ7v7xw0c5mzgruD2c8Hte1v/zxKRzqigoKDTzFwAynclIjt3ac1lsRnA7e6+ysz+CcgzsxOAXwJ/7+4JM1tiZiOC7RPu/m0z+58ku5FvmNlbwNWNFBZIdjV7g9v7SBakRpWXl7folSEHDhzIiFZZOdNLOdMrE3KOGTOGjz/+OPI5r15RweY9NaR75gItm7ucfFw35hT0adH+Gvu+N3dJrzXFZQBQ9+PAa0Ae8FUgB3jWzACOAfoH27wM4O5/NrPeLTzGR8E+qoO/9zS1cW5ubot2GoNfyhMpypleypk+eXl5GZGz5AhydvTcpa05WzNz+W/gm8Ht/ODvTUAFMMLdhwK/AdYE6wYBmNlAoO7lG4ebOeZrwKjg9nnAK63IJyKd3Icffsju3bvDjtGuFkzIp39ONl2zsuifk82CCfnN3ykErelcbgQWm9nNQCVwwN0rzWwmUGJmXYFy4HfB9qeb2UvA0cAVwbLVJGcuI939Lw0cY1pwjCuAD4Fxrf4XiUindeGFF8Zq5tKQvr16RvYVYqlaXFzcfSPJgXv95UuBpanLgktky+q/dNjdbwNua+IYO4BzW5pJRCTVjTfeSEVFRdgxhOh85P44M5vS0AozuxaY2LFxRCQTjR49mmHDhoUdQ2ind+i7+y9asW0RTbyPJXifywNHHEpEYm/79u1UVlZGfqDfGejjX0QkNsaOHRv7mUumUHERkdgoLCxky5YtYccQVFxEJEbOPffcyL+BsrOIykBfROSIVVRU8MEHH4QdQ1DnIiIxcvHFF2vmEhEqLiISG7fddhubN28OO4ag4iIiMTJ8+HDNXCJCMxcRiY2ysjK9Qz8i1LmISGxceumlJBIJRo4cGXaUTk/FRURi45e//KVmLhGh4iIisTFkyBDNXCJCMxcRiQ13Z9OmTWHHENS5iEiMXHXVVSQSCUaNGtX8xtKuVFxEJDbuuusuysvLw44hqLiISIycddZZHH/88WHHEDRzEZEYeeedd9iwYUPYMQR1LiISI9deey2JRILzzz8/7CidnoqLiMTGjBkz9GqxiFBxEZHYyM/PJzs7O+wYgmYuIhIjb731lt5EGRHqXEQkNiZPnkwikWDMmDFhR+n0VFxEJDZmzZpFWVlZ2DEEFRcRiZHTTjuNo446KuwYabNlV4LLFq+jrHI//XKOZsGEfPr26hl2rBZRcRGR2Fi3bh2bNm0iLy8v7CitNmnRWoq9stH1G3ZWMXhG8eeWD7McFk06oz2jtYmKi4jExs0330wikeCHP/xh2FEAGHlvCet3VDWxxZFfwiv2SnILn2n1/QaclM0LNww54uM3JhLFxcwmAncCs9x9pplNBnq7e6GZXQhMA55y98Iwc4pItD3wwAORmrk09eRdWlrabIc1YmYJGyurOFwLXbKgf042q6a0X0FIpyi9FPkRYLaZ/Svw47qF7r4cmB5aKhHJGAMHDuSUU04JO0baLJiQT98TkjOWw7VQc+gwW3YlQk7VMpHoXFJ0BxYDq4BTQ84iIhlm9erVlJeXx3LmAlC+K/G5uYtmLi3g7ruBF4LLZM0qLy+nurq62e0OHDiQEW+sUs70Us70yoSckydP5vDhw5x++ulhR/nE1Ssq2LynppG14cxcTj6uG3MK+rRo28a+780V8EgVl9bKzc1t0XYtubYZBcqZXsqZXpmQc8mSJWzcuDFSOUsaydLS8xn23KWt3/cozVxERI6ImfGVr3wl7BhptWBCPv1zsumalUX/nGwWTMgPO1KLZHTnIiKSqqSkhM2bN0eqczlSfXv1zJhXiKWKZHFx96KwM4hI5rnjjjtIJBJccsklYUfp9KJUXMaZ2U53n5m6MHifSyHwVDixRCRTLFy4kPfeey/sGEJEikvQqRQ1sm45sLwj84hIZurXrx8HDx4MO4aggb6IxMiLL77I6tWrw44hqLiISIxMmzaNuXPnhh1DiMhlMRGRdFiyZAkbNmwIO4ag4iIiMdKnTx+qqpr6FGLpKLosJiKx8fzzz/PKK6+EHUNQcRGRGJk+fTrz588PO4agy2IiEiPLli1j/fr1YccQVFxEJEZ69+7N7t27w44h6LKYiMTIypUrKS7+/O+Zl46n4iIisfHrX/+aoqKisGMIuiwmIjGyfPlyzVwiQp2LiMTGiSeeyPHHHx92DEHFRURi5IknnmDVqlVhxxBUXEQkRu6//36WLl0adgxBMxcRiZEVK1bg7mHHENS5iEiMHHvssRxzzDFhxxDUuYhIjDz66KNs3bqVvLy8sKN0eupcRCQ2Zs+ezbJly8KOIahzEZEYefbZZ3n33XfDjiGocxGRGOnZsyc9evQIO4agzkVEYmTp0qWauUSEOhcRiY358+fz+OOPhx1DUOciIsCWXQkuW7yOssr99Ms5mgUT8unbq2fYsVpt1apVlJaWhh1DUHER6bQmLVpLsVd+bvmGnVUMnpH82PphlsOiSWd0dLQ269atG926dQs7hqDiIhJrI+8tYf2Oqjbfv9gryS18JmVJ2WfWDzgpmxduGNLm/adbUVER27Zt08wlAiJRXMxsInAnMB8YTDJXFnAl8HVgGvCUuxeGlVEkE7X0iX/EzBI2VlZxuBa6ZEH/nGxWTfnsfUtLSyP/pF1UVEQikeDWW28NO0qnF6WB/iNAf+ABdx8K3AXc7e7LgelhBhOJu6kFA+naJQuArl2ymFowMOREbfP73/+exYsXhx1DiEjnkuJGYG9w+wvAgRCziMRaYzOXmkO1jJ235pOvM23uItGQVVtbG3aGustip9Zd9jIzA54Cvu/uXn99nfLy8trq6upm93/gwAG6d++e9tzpppzppZxw9YoKNu+paZd91zn5uG7MKejTrsdoqccee4yamhrGjRsXdpRmZfrjMy8vL6up+0Wtc8HMhgEPAhd7M5+dnZub26J9ZsK1YlDOdFNOKGnhfuMyc/nDH/7A/v37mTp1athRmpUJ5xPanjNKM5e6wnIfcK67/2fYeUQ6iwUT8umfk03XrCz652SzYEJ+2JHa5MUXX2ThwoVhxxCi17nMAr4ILE5eGcPd/apwI4nEX99ePT/XqYgciUgVF3f/RtgZRCRzPfjgg2zfvp0777wz7CidXpQui40zsyn1F5rZhYDe3yIizVq5ciXFxcVhxxAi0rm4exFQ1Mi65cDyjswjIpnpueee02eLRUSUOhcREYkJFRcRiY377ruPJUuWhB1DiMhlMRGRdHjppZeoqmr7B3VK+qi4iEhsPP3005q5RIQui4mISNqpuIhIbNxzzz16h35E6LKYiMTG66+/zr59+8KOIai4iEiMPP7445q5RIQui4mISNqpuIhIbEyfPp158+aFHUPQZTERiZG33nqLjz76KOwYgoqLiMTIsmXLNHOJCF0WExGRtFNxEZHYmDp1KrNnzw47hqDLYiISI+7O3r17w44hqLiISIwsXbpUM5eI0GUxERFJOxUXEYmN22+/nd/85jdhxxB0WUxEYqSiooI9e/aEHUNQcRGRGFm0aJFmLhGhy2IiIpJ2Ki4iEhs//elPmTlzZtgxBF0WE+nUtuxKcNnidZRV7qdfztEsmJBP3149w47VZrt27dL7XCJCxUWkE5m0aC3FXtngug07qxg8o/iTr4dZDosmndFR0dLioYce0swlIlRcRGJi5L0lrN9Rlbb9FXsluYXP1Fta9smtASdl88INQ9J2PImXSBQXM5sI3AnMBb4NfBH4CzAe+C4wDXjK3QvDyigSdW15oh8xs4SNlVUcroUuWdA/J5tVUxreT2lpKXl5eUcas13ddNNN7Nq1i0WLFoUdpdOL0kD/ESAHWOzu5wBvApe7+3JgeqjJRGJqasFAunbJAqBrlyymFgwMOdGRqa6u5uDBg2HHECLSuaS4Acgysy5AH2BzyHlEYqWpmUvNoVrGzlvzydeZOHP57W9/q5lLRGTV1taGnaHustip7l5oZt2APwLdgWHuvjl1fer9ysvLa6urq5vd/4EDB+jevXv6g6eZcqZXZ8t59YoKNu+pSUOiljv5uG7MKejTocdsTmf7vre3xnLm5eVlNXW/qHUuuHsN8DUzGw48DDR6ITk3N7dF+8yEa8WgnOnW2XKWtGEfcZu5TJ48mb/85S88/PDDYUdpViacT2h7zijNXDCzB81sWPDlPuBwmHlE4m7BhHz652TTNSuL/jnZLJiQH3YkiYmodS73A3PM7HaSheVHIecRibW+vXo22qlkolmzZmnmEhGRKi7u/i4wNOwcIiJyZKJ0WWycmU2pv9DMLgT0/hYRadaPf/xjpk6dGnYMISKdi7sXAUWNrFsOLO/IPCKSmXr06EEikQg7hhCR4iIikg733HOPZi4REaXLYiIiEhPqXEQkNq688kr27NnD7373u7CjdHoqLiISG7169SIKnzoiKi4iEiN33323Zi4RoZmLiIiknToXEYmNSZMmsWfPHp588smwo3R6Ki4iEht9+vShR48eYccQVFxEJEbuvPNOzVwiQjMXERFJO3UuIhIb48ePZ+/evaxcuTLsKJ2eiouIxIaZUVnZ8K9xlo6l4iIisfHzn/9cM5eI0MxFRETSTp2LiMTG2LFj+eijj3j22WfDjtLpqbiISGycdtpp7Ny5M+wYgoqLiMRIYWGhZi4RoZmLiIiknToXEYmNCy64gH379vHCCy+EHaXTU3ERkdj45je/yY4dO8KOIai4iEiM3HTTTZq5RIRmLiIiknbqXEQ6wJZdCS5bvI6yyv30yzmaBRPy6durZ9ixYuf888+nqqqKl19+OewonZ6Ki0iaTVq0lmJv/POtNuysYvCM4s8sG2b7WDTpjPaOFnvf+c53NHOJCBUXkRYaeW8J63dUtcu+i72S3MJnmt1uwEnZvHDDkHbJEAc/+clPNHOJiEgUFzObCNwJzHL3mWY2BFjq7n3M7EJgGvCUuxeGmVM6tyN5Uh8xs4SNlVUcroUuWdA/J5tVU5L7Ky0tJS8vL10xRSIhSgP9R4LC0geYAnQDcPflwPRQk4kcoakFA+naJQuArl2ymFowMORE8XTeeedx5ZVXhh1DiEjnUsfMugNzgCuBN0KOI9Imzc1cag7VMnbems8s08wlPUaPHs327dvDjiFAVm1tbdgZ6i6LnQqcCPyruxeb2XZ37526vv5lsfLy8trq6upm93/gwAG6d++e9tzpppzp1V45r15RweY9NWnfb0udfFw35hT06fDjdvbve7ples68vLyspu4Xpc7lBOAc4Ktmdgdwgpktc/exjd0hNze3RTvOlGvayple7ZWzpA37jMPMRTnTK+45o1Rc/uLuVvdF0Lk0WlhEMsmCCfmfe5+LpN/w4cPZv38/r7/+ethROr0oFReR2Orbq+cnnYq0n4suuogPPvgg7BhChItL3bxFRKSlrrjiCr3PJSKi9FLkcWY2pf7C4H0uen+LiEgGiUTn4u5FQFEj65YDyzsyj4hkpqFDh5JIJFi7dm3YUTq9SBQXEZF0mDhxItu2bQs7hqDiIiIxMnHiRM1cIiJKMxcRkSNSU1NDTU14b3CVT6lzEZHYGDFihGYuEaHiIiKxcfnll7N169awYwgqLiISI+PHj9fMJSI0cxGR2EgkErTkw2yl/alzEZHYGDVqlGYuEaHiIiKxcc0112jmEhEqLiISGxdddJFmLhGhmYuIxMbevXvZt29f2DEEdS4iEiMFBQWauUSEiouIxMb111/P+++/H3YMQcVFRGJkzJgxmrlEhGYuIhIbH374Ibt37w47hqDORURi5MILL9TMJSJUXEQkNm688UYqKirCjiGouIhIjIwePVozl4jQzEVEYmP79u1UVlaGHUNQ5yIiMTJ27FjNXCJCxUVEYqOwsJAtW7aEHUNQcRGRGDn33HM1c4kIzVza2ZZdCUbMLKH/T59lxMwStuxKhB1JJLYqKir44IMPwo4hqHM5YpMWraXYWzZA3LCzisEzihtcN8xyuOXMY9IZTaTTufjiizVziYhOXVxG3lvC+h1VYccAoNgrgyJV1ur7DjgpmxduGJL+UCIZ5rbbbmPz5s1hxxAiUlzMbCJwJzALuBV4J1j1JHAIuAmY7u5z0nncjnhCHjGzhI2VVRyuhS5Z0D8nm1VTGj5uaWkpeXl57Z5JJK6GDx+umUtERKK4BB4B/gT8P3e/LnWFmZ0YTqQjN7VgIBcv/A8OH6qla5csphYMDDuSSGyVlZVRUVGhH9IiIErFBWAQMMjMSoCdwPXuHunpXGtmLjWHahk7b02j6/P/rgeP6T+FSJtdeumlJBIJRo4cGXaUTi+rtrY27Ax1l8VOBV4H9rv7i2b2f4EfuPuFZvYLYHv9y2Ll5eW11dXVze7/wIEDdO/evcF1V6+oYPOemiP8F0TDycd1Y05Bn3Y/TlPnM0qUM70yIee6dev4+OOP+da3vhV2lGZlwvmExnPm5eVlNXW/qHUuLwN1r9V9kuQcplG5ubkt2mlTs4ySdu4U4jhzUc70Us70ycvLy4ickBnnE9qeM2rvc5kPXBDc/g7wRohZ0mLBhHz652TTNSuL/jnZLJiQH3YkkdhydzZt2hR2DCF6nUshsNDMfgTsBy4POc8R69urZ6Odioik11VXXUUikWDUqFFhR+n0IlVc3H0TMCzsHCKSme666y7Ky8vDjiFE67LYODObUn+hmV0LTOz4OCKSac466yxOP/30sGMIEelc3L0IKGpk3QPAAx2ZR0Qy0zvvvENZWVlGDMrjLhLFRUQkHa699loSiQTnn39+2FE6PRUXEYmNGTNm6NViEaHiIiKxkZ+fT3Z2dtgxhGgN9EVEjshbb72lD66MCHUuIhIbkydPJpFIMGbMmLCjdHoqLiISG7NmzaKsrPW/E0nST8VFRGLjtNNO46ijjgo7hqCZi4jEyLp163j77bfDjiGocxGRGLn55ptJJBL88Ic/DDtKp6fiIiKx8cADD2jmEhEqLiISGwMHDqRr165hxxA0cxGRGFm9ejVvvvlm2DEEdS4iEiO33noriUSCcePGhR2l01NxEZHYmDt3Lhs3bgw7hqDiIiIxYmYcPnw47BiCZi4iEiMlJSWsW7cu7BiCOhcRiZE77riDRCLBJZdcEnaUTk/FRURiY+HChbz33nthxxBUXEQkRvr168fBgwfDjiFo5iIiMfLiiy+yevXqsGMIKi4iEiPTpk1j7ty5YccQdFlMRGJkyZIlbNiwIewYgoqLiGSALbsSXLZ4HWWV++mXczQLJuTTt1fPz23Xp08fqqqqQkgo9am4iEgkTVq0lmKv/NzyDTurGDyj+JOvh1kOiyadAcDzzz/Pli1byMvL67Cc0jAVFxEJxch7S1i/48i7jGKvJLfwmZQlf8ddn/k6acBJ2bxww5AjPp60jIqLiISiNU/0I2aWsLGyisO10CUL+udks2rK5++/fft21q9fz+DBg9MZVdqgxa8WM7MiMzu3hdv2NbPRbY8lIvKpqQUD6dolC4CuXbKYWjCwwe169+5NTk5OR0aTRrRX5/Jt4FRgZTvtX0RirrGZS82hWsbOW/PJ16kzl5UrV1JRUaGZSwRk1dbWtmhDMysCTgCOJlmULnP398zsOmAcUAssA34L/BnoCVwL7AXuINklZQPj3H19I8f4BfBV4ESgV7CvC4ABwAR3X5O6fXl5eW11dXWz2Q8cOED37t1b9O8Mk3Kml3KmV7pzXr2igs17atK2v8acfFw35hT0affjtFamf9/z8vKymrpfazuX1e4+3cxGAf9iZrcBFwFnB+tXAf8OTAdOdfenzexHwHh332ZmtwL/CPyqiWNUu/u5ZlYIjHL30WY2CRgLfKa45Obmtih0aWlpRvwko5zppZzple6cJa3YV0tnLh9++CHr16/nrLPOSlvO9hL373tr36H/h+Dv1YABA4GTgZeCP72AU+rdZytwf9D5DAO6NXOM/wr+3gP8d3B7NxD9Ei8i7WLBhHz652TTNSuL/jnZLJiQ3+B2J554Iscff3wHp5OGtLZzOYNkYTkHeAdwkpfAznP3WjO7AfgTMJRPC9c8oL+77zOzxUCTrRTJy2siIp/o26tng51KfU888QTvv/9+RnQEcdfa4nKmmb1MsgBc6u6bzewl4FUzOwpYS7JTeRv4mZn9F7AUeMXM9gM7gL9NX3wRkU/df//9JBIJrr/++rCjdHotLi7uPrGR5TOAGfUWv0nyshkkh/wtPcYvUm7PSbn9FPBUS/cjIp3TihUrcPewYwghvYnSzJ4g+cqzVHvdvSCMPCISD8ceeyzHHHNM2DGEkIqLu48J47giEm+PPvooW7du1cwlAvT7XEQkNmbPns2yZS2+Ei/tSJ8tJiKx8eyzz/Luu++GHUNQ5yIiMdKzZ0969OgRdgxBnYuIxMjSpUs1c4kIdS4iEhvz58/n8ccfDzuGoM5FRGJk1apVlJaWhh1DUOciIjHSrVs3unVr7uMLpSOocxGR2CgqKmLbtm2auUSAOhcRiY2ioiKeekqfFBUFLf5lYVH0xhtvZG54EZEMN2jQoEY/5T6ji4uIiESTLouJiEjaqbiIiEjaqbiIiEjaxe7EBpVkAAAEFUlEQVSlyGZ2NPAIcDzwMTDB3bfW2+YK4Crgr8A0d/+3EHIeS/K3dH4J+CIwxd1fr7fNfcDZwL5gUYG7741gztDPZ708PwD+0d3HNbAu9HOakqWpnKGfUzPrQfJ7/2WS52uCu1fW22YFcCJQA1S7+3kdlK0L8CDwDeAgcLm7v5eyPvTzl5KluaxRekz+b+Cf3X1oveWjgdtJns+F7j6vuX3FsXO5AnjD3QeT/I9xS+pKM+sNXA98C/gucHfwK5o72hTgJXcfAkwEftvANoOA77r70OBPGA+4JnNG6HzW5bkPuJvGH9tROKdN5ozQOb0GeNvdzwEeBm5rYJtTgLODc9khhSXwfaC7u38TKAR+XbciQuevTqNZA1F5TN4CzAe611veDbgXGAkMAa40s5Oa21/siou7zwJ+FXzZF9hTb5MzgNfc/WDwTXwP+F8dGLHOvcDc4PYXgAOpK4Ofdk4BHjKz18zs0g7OV6fJnETnfNZZTfJJ8XMidE6hiZxE55yeDTwf3H4OGJ66MniCOQ5YaWavmtn3wsjm7muAv09ZF5XzV6fRrBF7TG4EGvpFjnnAe+6+290/Bl4FBje3s4y+LGZmlwE31Fs8yd3XmdnLwNeBEfXWfwlI/clgH3Bs+6VsNmdvkh3W5HrrjwZ+A8wEugLFZvaf7v6niOXs8PMJTWZ91MyGNnK3KJ3TpnJG5TG6IyVHQxm+SPKn8PtI/try18xsrbvvbM+sgfrn6JCZfcHd/9rAug55TDahqawd/phsjLs/bma5Daxq0/nM6OLi7guABY2s+7aZnQo8A/RPWfURkPpLto/h891NWjWW08y+DiwDbnL3knqrE8B97p4Itn2Z5DXbdnvQtTFnh59PaPp734TInNNmROIxamZPpORoKMN2YE7wJLnTzN4EDOiI4lL/HHUJcjS0rkMek01oKmuHPybboE3nM6OLS0PM7KfA++6+BKgCDtXbZC3wKzPrDhxFsuV7p2NTgpl9DXgMuMjd/9jAJgOAR83sdJKXL88GFndgRKBFOSNxPlsoEue0BaJyTl8DRgV5zgNeqbd+OHAdMMrMsoGBQEd9JPFrwGjgd2Z2JvB2yrqonL86TWXNhMdkKXCKmZ1A8jl1MHBPc3eKXXEBFgKLgza/KzAJwMymkLxu+LSZ3U/yP0oX4GfuXn+O0BHuJjk4u8/MAPa6e0G9nEuANSRfifOwu/85ojmjcD4bFcFz2qAIntPZJP8vvUrylZfjgpz/Aix39+fM7LtmtgY4DNzq7h92ULYngRFmthrIAiZF8Py1NGskH5NmNg7IdveHgrz/TvJ8Lqz/CtyG6ONfREQk7WL3ajEREQmfiouIiKSdiouIiKSdiouIiKSdiouIiKSdiouIiKSdiouIiKTd/wcLx9e9VpKXIAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.forestplot(model_dept.trace, varnames=['dept', 'beta_m'], rhat=False)\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/petrushev/etc/jup_env/lib/python3.7/site-packages/pymc3/stats.py:218: UserWarning: For one or more samples the posterior variance of the\n",
      "        log predictive densities exceeds 0.4. This could be indication of\n",
      "        WAIC starting to fail see http://arxiv.org/abs/1507.04544 for details\n",
      "        \n",
      "  \"\"\")\n",
      "/home/petrushev/etc/jup_env/lib/python3.7/site-packages/pymc3/stats.py:218: UserWarning: For one or more samples the posterior variance of the\n",
      "        log predictive densities exceeds 0.4. This could be indication of\n",
      "        WAIC starting to fail see http://arxiv.org/abs/1507.04544 for details\n",
      "        \n",
      "  \"\"\")\n"
     ]
    },
    {
     "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>WAIC</th>\n",
       "      <th>pWAIC</th>\n",
       "      <th>dWAIC</th>\n",
       "      <th>weight</th>\n",
       "      <th>SE</th>\n",
       "      <th>dSE</th>\n",
       "      <th>var_warn</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>department</th>\n",
       "      <td>108.47</td>\n",
       "      <td>9.41</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>14.38</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gender</th>\n",
       "      <td>984.38</td>\n",
       "      <td>113.01</td>\n",
       "      <td>875.92</td>\n",
       "      <td>0</td>\n",
       "      <td>292.81</td>\n",
       "      <td>308.1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              WAIC   pWAIC   dWAIC weight      SE    dSE var_warn\n",
       "department  108.47    9.41       0      1   14.38      0        1\n",
       "gender      984.38  113.01  875.92      0  292.81  308.1        1"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pm.compare({\n",
    "        model_gender: model_gender.trace,\n",
    "        model_dept: model_dept.trace,\n",
    "    },\n",
    "    method='BB-pseudo-BMA'\n",
    ")   "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}