03 Chapter 3/01 H1,2,3,4,5.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pymc3 as pm\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.style.use('seaborn-whitegrid')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "births = pd.DataFrame(dict(\n",
    "    first = [\n",
    "         1,0,0,0,1,1,0,1,0,1,0,0,1,1,0,1,1,0,0,0,1,0,0,0,1,0,\n",
    "         0,0,0,1,1,1,0,1,0,1,1,1,0,1,0,1,1,0,1,0,0,1,1,0,1,0,\n",
    "         0,0,0,0,0,0,1,1,0,1,0,0,1,0,0,0,1,0,0,1,1,1,1,0,1,0,\n",
    "         1,1,1,1,1,0,0,1,0,1,1,0,1,0,1,1,1,0,1,1,1,1],\n",
    "    second = [\n",
    "         0,1,0,1,0,1,1,1,0,0,1,1,1,1,1,0,0,1,1,1,0,0,1,1,1,0,\n",
    "         1,1,1,0,1,1,1,0,1,0,0,1,1,1,1,0,0,1,0,1,1,1,1,1,1,1,\n",
    "         1,1,1,1,1,1,1,1,1,0,1,1,0,1,1,0,1,1,1,0,0,0,0,0,0,1,\n",
    "         0,0,0,1,1,0,0,1,0,0,1,1,0,0,0,1,1,1,0,0,0,0],\n",
    "))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_births = births['first'].count() + births['second'].count()\n",
    "n_boys =  births['first'].sum() + births['second'].sum()\n",
    "\n",
    "first_births = births['first'].count()\n",
    "second_births = n_births - first_births\n",
    "first_boys = births['first'].sum()\n",
    "second_boys = n_boys - first_boys\n",
    "\n",
    "births_after_girl = births[births['first'] == 0]['second'].count()\n",
    "boys_after_girl = births[births['first'] == 0]['second'].sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**H1** Using grid approximation, compute the posterior distribution for the probability of a birth\n",
    "being a boy. Assume a uniform prior probability. Which parameter value maximizes the posterior\n",
    "probability?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/petrushev/etc/jup_env/lib/python3.7/site-packages/pymc3/tuning/starting.py:61: UserWarning: find_MAP should not be used to initialize the NUTS sampler, simply call pymc3.sample() and it will automatically initialize NUTS in a better way.\n",
      "  warnings.warn('find_MAP should not be used to initialize the NUTS sampler, simply call pymc3.sample() and it will automatically initialize NUTS in a better way.')\n",
      "logp = -7.0748, ||grad|| = 11: 100%|██████████| 6/6 [00:00<00:00, 700.18it/s]\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as model_h1:\n",
    "    # uninformed prior for the gender\n",
    "    p = pm.Uniform('p', lower=0, upper=1)\n",
    "\n",
    "    # both first- and second born are modeled with the same p\n",
    "    boys = pm.Binomial(\n",
    "        'boys',\n",
    "        n=[first_births, second_births],\n",
    "        p=[p, p],\n",
    "        observed=[first_boys, second_boys],\n",
    "    )\n",
    "\n",
    "    start = pm.find_MAP()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5549999989296014"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.asscalar(start['p'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**H2** Using the sample function, draw 10,000 random parameter values from the posterior distribution you calculated above. Use these samples to estimate the 50%, 89%, and 97% highest posterior\n",
    "density intervals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "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: [p]\n",
      "Sampling 2 chains: 100%|██████████| 12000/12000 [00:05<00:00, 2040.28draws/s]\n"
     ]
    }
   ],
   "source": [
    "with model_h1:\n",
    "    model_h1.trace = pm.sample(5000, chains=2, tune=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "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>interval</th>\n",
       "      <th>lower</th>\n",
       "      <th>upper</th>\n",
       "      <th>mean_</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.50</td>\n",
       "      <td>0.526177</td>\n",
       "      <td>0.573165</td>\n",
       "      <td>0.554412</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.89</td>\n",
       "      <td>0.498953</td>\n",
       "      <td>0.610093</td>\n",
       "      <td>0.554412</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.97</td>\n",
       "      <td>0.479315</td>\n",
       "      <td>0.629175</td>\n",
       "      <td>0.554412</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   interval     lower     upper     mean_\n",
       "0      0.50  0.526177  0.573165  0.554412\n",
       "1      0.89  0.498953  0.610093  0.554412\n",
       "2      0.97  0.479315  0.629175  0.554412"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_sample = pm.trace_to_dataframe(model_h1.trace, varnames=['p']).p\n",
    "\n",
    "pd.DataFrame(\n",
    "    [[interval] + list(pm.hpd(p_sample, alpha=1-interval)) + [p_sample.mean()]\n",
    "     for interval in (.50, .89, .97)],\n",
    "    columns=['interval', 'lower', 'upper', 'mean_']\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**H3** Use rbinom to simulate 10,000 replicates of 200 births. You should end up with 10,000 numbers, each one a count of boys out of 200 births. Compare the distribution of predicted numbers of boys to the actual count in the data (111 boys out of 200 births). There are many good ways to visualize the simulations, but the dens command (part of the rethinking package) is probably the easiest way in this case. Does it look like the model fits the data well? That is, does the distribution of predictions include the actual observation as a central, likely outcome?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 10000/10000 [00:07<00:00, 1282.00it/s]\n"
     ]
    }
   ],
   "source": [
    "with model_h1:\n",
    "    ppc = pm.sample_posterior_predictive(\n",
    "        model_h1.trace,\n",
    "        vars=[boys], samples=10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "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>first</th>\n",
       "      <th>second</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>8142</th>\n",
       "      <td>53</td>\n",
       "      <td>49</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4504</th>\n",
       "      <td>58</td>\n",
       "      <td>49</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6408</th>\n",
       "      <td>51</td>\n",
       "      <td>54</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>385</th>\n",
       "      <td>64</td>\n",
       "      <td>68</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5733</th>\n",
       "      <td>55</td>\n",
       "      <td>49</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      first  second\n",
       "8142     53      49\n",
       "4504     58      49\n",
       "6408     51      54\n",
       "385      64      68\n",
       "5733     55      49"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "boys_ppc = pd.DataFrame(ppc['boys'], columns=['first', 'second'])\n",
    "boys_ppc.sample(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiAAAAFLCAYAAADvUjcYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAMTQAADE0B0s6tTgAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X9wVPW9//FXkl1EAiZZh2yAZFqLjlRL1VFDYouEYIJLQggRi7UGpf5oq9JeBhEw1TRgkGKKPwa1l0v7xaZyr2hAEpOFXipojJNQbusVpRWl/thg2M2yBJNozO7mfP+g7iUSSNTkbJJ9Pmacyfmxn/N5f/yweeXsOWejDMMwBAAAYKLocHcAAABEHgIIAAAwHQEEAACYjgACAABMRwABAACms4S7A2fS1dWlQCCg6OhoRUVFhbs7AACgDwzDUFdXlywWi6Kjez7XMagDSCAQ0P79+8PdDQAA8BVMnjxZI0aM6HHboA4gn6emyZMnKyYmJsy9MU9jY6OSk5PD3Y2wivQxoP5T6z98WLroohM/HzggTZgQho6ZiDkQ2fVLQ3sMgsGg9u/ff9qzH9IgDyCff+wSExMTUQEkOjo6ourtSaSPAfWfWn90tNTe/n8/D/fhYQ5Edv3S8BiDM10+wUWoAADAdAQQAABgOgIIAAAwHQEEAACYjgACAABMRwABAACm61MAcTqdysnJUXZ2ttavX3/KdrfbrcLCQjkcDi1YsEBHjx7ttr2trU3XXHONGhoaQuvKy8vlcDiUnZ2tZ5999muWAQAAhpJeA0hzc7PWrl2r8vJyVVdXa9++faqtre22T0lJiQoKCuR0OpWXl6fS0tJu21etWqWPP/44tHzgwAFt2bJFFRUV2rp1q/74xz/q0KFD/VQSAAAY7HoNIHV1dUpLS5PNZpPValV+fr5qampC2/1+vxoaGpSbmytJys/P1549e+T3+yVJNTU1io2N1YUXXhh6ze7duzVz5kyNGjVKo0eP1syZM+V0Ovu7NgAAhpzGxkZlZmaGuxvdDESfen0SqsfjUWJiYmg5MTFRbrc7tNzS0qLY2FhZrdYTDVosGj16tHw+n4LBoJ5++mk9/fTTuv3220Ovcbvdmjx5cmjZbrfrjTfeOG0fGhsbz/g41+GmtbVVLpcr3N0Iq0gfA+o/tf6mphhJ4yVJH330kQwjGIaemYc5MPD1BwJSc/PAPWl07NigLF/heeNHjhxRIBAYVHPg8z71tT9dXV297tPr0PTUyMmPVj3TQYqKinT//fdr5MiR3dYbhnHGNr8oOTl5yD+O9stwuVxKSUkJdzfCKtLHgPpPrf/kt4jx48driH5FRp8xBwa2/kBA+s53pLffHrBD6MILpTffVK8h5Le//a0qKysVExOj733ve7rxxhsVDAb16KOP6siRI0pJSdHq1asVFxen9evXa+fOnbJYLJo8ebJWrlypYDCosrIy1dfXKxAI6Nprr9Vdd92lhoYGrV27VoZh6Jvf/Kb27t2riooK2e12ffLJJ7rmmmv03//933rvvfe0evVqffrppxozZoyKi4s1ceJEHThwQEVFRZKkSZMmyWKx9Pn/STAYlNfrPeM+vQaQpKQk7d27N7Ts8XiUlJQUWrbZbGpra1MgEJDFYlEgEFB7e7uOHTumf/7zn6HOf/jhh/rlL3+pkpIS2e12eTye07YJACfbsEE6duz/lp95RkpI+Gpt3XFH//QJ6A8vv/yydu3apYqKClmtVi1atEi1tbU6evSo8vPzNWvWLK1du1ZPPPGEFi1apM2bN6u2tlZRUVEqKSlRU1OTXnnlFfn9fm3dulWBQEB33nmndu/erVGjRum9997T7t27FRcXp9WrV8vpdOqWW27Rrl279L3vfU8jRozQfffdpyeffFLJycn6n//5Hy1dulRbt27VsmXLdO+992rq1Kl64oknut1I0h96DSDp6el6/PHH5fV6FRcXp8rKSv3whz8MbbdarUpNTVVVVZXmzp2rqqoqpaamatKkSXr55ZdD+xUWFuruu+/WlClTNGbMGBUVFemWW26RJO3YsUNr1qzp18IAADgdi+XE2YkjRwbuGElJvZ/9qK+vV25urs4++2xJ0nXXXacXXnhB3/jGN0KXKuTl5Wn58uUaM2aMJk6cqHnz5mn69Om66aabNG7cONXV1enAgQP6y1/+Ikn69NNPdfDgQV166aX61re+pbi4OEknrtH81a9+pVtuuUUvvviiFixYoPfee08ffvih7rrrrlCffD6fjh49KrfbralTp0qSCgoKVFFR0a/j02sAsdvtWrp0qRYuXKjOzk5lZmYqKytLRUVFyszM1IwZM1RcXKwVK1Zo48aNiouLU1lZ2RnbnDx5sq6//nr94Ac/UCAQ0A033KCLL76434oCAKA3FovC/lFeT5cxBAKBU657tPwryfzhD3/QX//6V9XW1urWW29VWVmZgsGg7r33XmVnZ0uSjh07ppEjR+qNN94IBRtJuuiii/Tpp5/qzTff1Lvvvqv09HS98847SklJ0fbt2yWduETC7XYrOjq62+USlq9yMUsv+tSiw+GQw+Hotu7kW23HjRunTZs2nbGN8vLybsuFhYUqLCzsYzcBABh+0tLS9MQTT2j+/PmyWq2qqKjQlVdeqZdfflkHDx5USkqKnn/+eV111VVqbGzUz372Mz333HO6/PLL5Xa79fbbbystLU1btmzR9OnTFQgEdMstt2jRokUaM2bMKcebM2eOHnjgAc2cOVMxMTH61re+pePHj+svf/mLrrzySlVVVen3v/+9XnjhBU2YMEG7du3SNddco+rq6n6vvf8jDQAA6JPp06fr73//u+bNm6dAIKCrrrpKM2bM0JYtW7R582Y9+uijuuCCC7R48WLFxsZq1qxZmjt3rs4++2yNHz9ec+fO1VlnnaUPPvhAc+fOld/vV05OzikP//zc7NmztW7dOq1evVqSNGLECD322GNavXq1Ojo6NGrUqNCnGA8//LBWrFih9evX69JLL+332qOMnm5JGSSCwaBef/11XXrppdwFE2EifQyov3v9n1+Eunz5ieU1a4b/RajMgciuXxraY9CX39+R83ANAAAwaBBAAACA6QggAADAdAQQAABgOgIIAAAwHQEEAACYjueAAAAi1oYNA9f2ULnlO1w4AwIAwDDX2NiozMzML/WazMxMNTY2DlCPOAMCYAB83b8qfb5Y2Wz905cv6u+/ePkrF/hqCCAAAISJz+fTkiVLdPz4cUnSnXfeqXPOOUfr1q1Ta2urPv30Uy1ZskQ5OTlavny5Ro4cqddff10tLS1avHixdu3apb///e+aPn26ioqKtHXrVu3cuVNtbW3yeDyaNm2aioqKTjnmAw88oI8++kiSdPfddyszM1MtLS1aunSpmpqadP755+uzzz4b0NoJIAAAhElVVZUuvPBCLV++XIcOHdKWLVvU1NSkVatWaeTIkTp8+LBKS0uVk5MjSXK73XrhhRe0bds2rVq1Sjt37tRZZ52lq6++WosWLZIkvfHGG6qsrFRcXJwWLFigHTt2aPLkyaFjlpaWKi8vT9nZ2fL5fJo/f74uueQSPfHEE5o0aZL+4z/+Q/v27ZPT6RzQ2gkgAACEyZQpU3TrrbeqsbExFCKsVqt2796tv/3tb3rvvffU3t4e2j8jI0OSNH78eF1wwQU699xzJUnx8fH6+OOPJZ34gruxY8dKkmbNmqW//OUv3QLIq6++qnfeeUdPPPGEJCkQCOjQoUPau3dv6IvorrjiigH/HhoCCAAAYTJp0iTt3LlTr7zyinbv3q3f//73io2NVWpqqi644ALNnDlT99xzT2h/q9Ua+tli6flX+Mnru7q6FB3d/X6Trq4u/eEPf1B8fLwkyePxyGazKSoq6rTtDATuggEAIEw2bNigjRs3atasWfrVr36l999/X++//77+7d/+TVOmTFFdXZ2CweCXarO2tlYff/yxPvvsM1VXV+v73/9+t+1paWnavHmzJOn9999Xbm6ujh8/rvT0dG3btk2S9Oabb+qDDz7onyJPgzMgAICIFe67mObNm6d77rlHs2fPlsViUVFRkQ4fPqycnByNGDFCV155pTo6Orp9DNObsWPH6ic/+Yl8Pp9yc3OVkZHR7XbaX/7ylyouLtbs2bNlGIZKS0t17rnn6uc//7mWL1+unJwcfetb3xrwj2CiDMMwBvQIX0MwGNTrr7+uSy+9VDExMeHujmlcLteA/48f7CJ9DIZ6/V//NlyfbF+4D/fYMWn58hM/r1kjJSR8vWP0l4H6BTbU58DXFen1S19tDLZu3aq9e/dqzZo1A9SrvunL728+ggEAAKbjIxgAAIaJgoICFRQUhLsbfcIZEAAAYDoCCAAAMB0BBAAAmI4AAgAATEcAAQAApiOAAAAA0xFAAACA6QggAADAdH1+EJnT6dT69evl9/uVl5enu+++u9t2t9ute+65R16vV2PHjtUjjzyic889Vx999JGWL1+ulpYWnXXWWVq5cqW+/e1vS5KmTp3a7XHLv/3tbzVu3Lh+Kg0AAAxWfToD0tzcrLVr16q8vFzV1dXat2+famtru+1TUlKigoICOZ1O5eXlqbS0VJK0Zs0azZ49W5WVlVq0aJFKSkokSU1NTUpISND27dtD/xE+AACIDH06A1JXV6e0tLTQ2Yr8/HzV1NRo6tSpkiS/36+GhgY99thjoe2rV6+W3+/Xo48+GmqnsbFR55xzjiRp//79CgQCuvHGG9XR0aGf/vSnys7O7tfiAPTd1/0COQD4MvoUQDwejxITE0PLiYmJcrvdoeWWlhbFxsbKarWeaNRi0ejRo+Xz+WS32yVJ2dnZ+uijj/TUU09JOhFapk2bpiVLlqipqUk33XSTJk6cqIkTJ55y/MbGRkVHR87lKq2trXK5XOHuRlhF+hiEo36fL9bU451JR0eHfD5ft3XHj0dJOvEVuC0txzRYvsjb5er716R/GfwbiOz6paE9Bl1dXb3u06cA0lNDUVFRvR7o5NDwpz/9SW+99ZZuvfVW7dixQzk5OcrJyZEkpaSkKCsrS3V1dT0GkOTk5NN+ne9wxNdQMwbhqP+ky7HCzufzdbs+TJJOestRfHyCEhJM7tRppKQMzMDxbyCy65eG9hgEg0F5vd4z7tOn0wpJSUlqbm4OLXs8HiUlJYWWbTab2traFAgEJEmBQEDt7e2Kj4/Xrl271NHRIUm6+OKLNWHCBLlcLu3YsUP/+Mc/Qm0YhhFRIQMAgEjWpwCSnp6u+vp6eb1e+f1+VVZWKiMjI7TdarUqNTVVVVVVkqSqqiqlpqbKarXqueeeU0VFhSTp4MGDOnr0qCZOnKgPPvhATz31lAzDkMfj0Z///OdubQIAgOGrTx/B2O12LV26VAsXLlRnZ6cyMzOVlZWloqIiZWZmasaMGSouLtaKFSu0ceNGxcXFqaysTJJUXFys++67T1u2bNFZZ52ldevWadSoUVq4cKHuv/9+5eTkyDAM3XvvvZowYcKAFgsAAAaHPj8HxOFwyOFwdFv3+a22kjRu3Dht2rTplNeNHz++x/UjRozQr3/96773FAAADBuRc2sJAAAYNAggAADAdAQQAABgOgIIAAAwHQEEAACYjgACAABMRwABAACmI4AAAADTEUAAAIDpCCAAAMB0BBAAAGA6AggAADAdAQQAAJiOAAIAAExHAAEAAKYjgAAAANMRQAAAgOkIIAAAwHQEEAAAYDoCCAAAMB0BBAAAmI4AAgAATEcAAQAApiOAAAAA0xFAAACA6Szh7gAADGUbNvRve3fc0b/tAYMVZ0AAAIDp+hxAnE6ncnJylJ2drfXr15+y3e12q7CwUA6HQwsWLNDRo0clSR999JEWLFigvLw8XX/99fr73/8ees26det07bXXaubMmXrppZf6oRwAADAU9CmANDc3a+3atSovL1d1dbX27dun2trabvuUlJSooKBATqdTeXl5Ki0tlSStWbNGs2fPVmVlpRYtWqSSkhJJ0q5du7R//369+OKLKi8vV2lpqY4fP97P5QEAgMGoTwGkrq5OaWlpstlsslqtys/PV01NTWi73+9XQ0ODcnNzJUn5+fnas2eP/H6/Hn30UV133XWSpMbGRp1zzjmSpJdeekl5eXmyWCxKTExUamqqdu/e3d/1AQCAQahPF6F6PB4lJiaGlhMTE+V2u0PLLS0tio2NldVqPdGoxaLRo0fL5/PJbrdLkrKzs/XRRx/pqaeeknTiI5vPt33e5pEjR3o8fmNjo6KjI+dyldbWVrlcrnB3I6wifQzCUb/PF2vq8c6ko6NDPp+v27rjx6MkJUiSWlqOyTCMMPRs4Llc7ZL4NxDp9UtDewy6urp63adPAaSnhqKiono90Mmh4U9/+pPeeust3XrrrdqxY0ePbx6nCxnJycmKiYnpS1eHBZfLpZSUlHB3I6wifQzCUb/NZurhzsjn88n2hQ6d9Jaj+PgEJSSY3CmTpKScqJt/A5FdvzS0xyAYDMrr9Z5xnz6dVkhKSlJzc3No2ePxKCkpKbRss9nU1tamQCAgSQoEAmpvb1d8fLx27dqljo4OSdLFF1+sCRMmyOVyyW63y+PxhNpobm7udkYEAAAMX30KIOnp6aqvr5fX65Xf71dlZaUyMjJC261Wq1JTU1VVVSVJqqqqUmpqqqxWq5577jlVVFRIkg4ePKijR49q4sSJysjIUGVlpQKBgLxer+rr63XVVVf1f4UAAGDQ6dNHMHa7XUuXLtXChQvV2dmpzMxMZWVlqaioSJmZmZoxY4aKi4u1YsUKbdy4UXFxcSorK5MkFRcX67777tOWLVt01llnad26dRo1apSys7P15ptvas6cOQoGg1q8eLHGjh07oMUCAIDBoc9PQnU4HHI4HN3WfX6rrSSNGzdOmzZtOuV148eP73F9VFSUlixZoiVLlvS9twAAYFiInFtLAADAoEEAAQAApiOAAAAA0xFAAACA6QggAADAdAQQAABgOgIIAAAwHQEEAACYjgACAABM1+cnoQIYXDZsCHcPAOCr4wwIAAAwHQEEAACYjgACAABMRwABAACmI4AAAADTEUAAAIDpCCAAAMB0BBAAAGA6AggAADAdAQQAAJiOAAIAAExHAAEAAKYjgAAAANMRQAAAgOkIIAAAwHQEEAAAYDoCCAAAMJ2lLzs5nU6tX79efr9feXl5uvvuu7ttd7vduueee+T1ejV27Fg98sgjOvfcc3Xs2DEVFRWpsbFRhmHopz/9qXJyciRJ8+fPV3t7u2JiYiRJK1eu1CWXXNLP5QEAgMGo1wDS3NystWvXqqKiQmPGjNHtt9+u2tpaTZ06NbRPSUmJCgoKNHfuXD3//PMqLS3VunXr9Pjjj+uiiy7Sk08+qebmZs2dO1dTpkxRXFycDh8+rFdeeUXR0ZyEAQAg0vT627+urk5paWmy2WyyWq3Kz89XTU1NaLvf71dDQ4Nyc3MlSfn5+dqzZ4/8fr+uvvpqzZ8/X5I0duxYxcfHy+v16uDBg7JYLPrxj3+svLw8bd68eYDKAwAAg1GvZ0A8Ho8SExNDy4mJiXK73aHllpYWxcbGymq1nmjQYtHo0aPl8/k0ffr00H7V1dXq7OzU+eefr3379ik9PV3FxcX65JNPtGDBAp133nlKT0/vsQ+NjY0RdaaktbVVLpcr3N0Iq0gfg77U7/PFmtQb83V0dMjn83Vbd/x4lKQESVJLyzEZhhGGng08l6tdEv8GIr1+aWiPQVdXV6/79BpAemokKiqq14OcHBi2b9+uhx9+WBs3bpTFYlFaWprS0tIkSSNHjtS8efO0Z8+e0waQ5OTk0LUikcDlciklJSXc3QirSB+DvtRvs5nUmTDw+XyyfaHAk952FB+foIQEkztlEqfzRN09jcGXdccd/dGj8Ij09wBpaI9BMBiU1+s94z69nlZISkpSc3NzaNnj8SgpKSm0bLPZ1NbWpkAgIEkKBAJqb29XfHy8JGnDhg167LHH9PTTT2vSpEmSpNdee0379u0LtdHV1RVRAQMAgEjXawBJT09XfX29vF6v/H6/KisrlZGREdputVqVmpqqqqoqSVJVVZVSU1NltVq1detWbdu2Tc8++6wmTpwYeo3P51NZWZk6OzvV1tambdu2KSsrq/+rAwAAg1KvH8HY7XYtXbpUCxcuVGdnpzIzM5WVlaWioiJlZmZqxowZKi4u1ooVK7Rx40bFxcWprKxMkvTII48oKipKt912W6i9lStXKicnR2+99ZbmzJmjrq4u/ehHP9Jll102cFUCAIBBpU/PAXE4HHI4HN3WlZaWhn4eN26cNm3adMrramtrT9vmsmXLtGzZsj52EwAADCeRc2sJAAAYNAggAADAdAQQAABgOgIIAAAwHQEEAACYjgACAABMRwABAACmI4AAAADTEUAAAIDpCCAAAMB0BBAAAGA6AggAADAdAQQAAJiOAAIAAExHAAEAAKYjgAAAANMRQAAAgOkIIAAAwHSWcHcAiBQbNvR9X58vVjbbwPUFAMKNMyAAAMB0BBAAAGA6AggAADAdAQQAAJiOAAIAAExHAAEAAKYjgAAAANMRQAAAgOn6FECcTqdycnKUnZ2t9evXn7Ld7XarsLBQDodDCxYs0NGjRyVJx44d05133qm8vDzNnj1b1dXVodeUl5fL4XAoOztbzz77bD+VAwAAhoJeA0hzc7PWrl2r8vJyVVdXa9++faqtre22T0lJiQoKCuR0OpWXl6fS0lJJ0uOPP66LLrpIlZWV+v3vf6+HHnpIXq9XBw4c0JYtW1RRUaGtW7fqj3/8ow4dOjQwFQIAgEGn1wBSV1entLQ02Ww2Wa1W5efnq6amJrTd7/eroaFBubm5kqT8/Hzt2bNHfr9fV199tebPny9JGjt2rOLj4+X1erV7927NnDlTo0aN0ujRozVz5kw5nc4BKhEAAAw2vX4XjMfjUWJiYmg5MTFRbrc7tNzS0qLY2FhZrdYTDVosGj16tHw+n6ZPnx7ar7q6Wp2dnTr//PO1efNmTZ48ObTNbrfrjTfeOG0fGhsbFR0dOZertLa2yuVyhbsbYTUcx8Dni+3zvh0dHfL5fAPYm8Gtp/qPH4+SlCBJamk5JsMwwtAz8/THHHC52vupN+Ybju8BX9ZQHoOurq5e9+k1gPTUSFRUVK8HOTkwbN++XQ8//LA2btwoi8XS4xvHyW1+UXJysmJiYnrr6rDhcrmUkpIS7m6E1XAcgy/z5XI+n0+2CP42up7qP/ktIj4+QQkJJnfKZP0xB1JShu4cGo7vAV/WUB6DYDAor9d7xn16Pa2QlJSk5ubm0LLH41FSUlJo2Wazqa2tTYFAQJIUCATU3t6u+Ph4SdKGDRv02GOP6emnn9akSZMknTjj4fF4TtsmAAAY3noNIOnp6aqvr5fX65Xf71dlZaUyMjJC261Wq1JTU1VVVSVJqqqqUmpqqqxWq7Zu3apt27bp2Wef1cSJE0OvmTZtmnbu3Kn29na1t7drx44dmjZtWv9XBwAABqVeP4Kx2+1aunSpFi5cqM7OTmVmZiorK0tFRUXKzMzUjBkzVFxcrBUrVmjjxo2Ki4tTWVmZJOmRRx5RVFSUbrvttlB7K1eu1CWXXKLrr79eP/jBDxQIBHTDDTfo4osvHrgqAQDAoNJrAJEkh8Mhh8PRbd3nt9pK0rhx47Rp06ZTXvfF23VPVlhYqMLCwj52EwAADCeRc2sJAAAYNAggAADAdAQQAABgOgIIAAAwHQEEAACYjgACAABMRwABAACmI4AAAADTEUAAAIDpCCAAAMB0BBAAAGA6AggAADAdAQQAAJiOAAIAAExHAAEAAKYjgAAAANNZwt0BAMDA2LChf9u7447+bQ+RjTMgAADAdAQQAABgOgIIAAAwHQEEAACYjgACAABMRwABAACmI4AAAADTEUAAAIDpCCAAAMB0BBAAAGC6PgcQp9OpnJwcZWdna/369adsd7vdKiwslMPh0IIFC3T06NFu21977TXdfPPN3dZNnTpVc+bMCf3X1NT0FcsAAABDSZ8CSHNzs9auXavy8nJVV1dr3759qq2t7bZPSUmJCgoK5HQ6lZeXp9LSUklSMBjU7373Oy1evFhdXV2h/ZuampSQkKDt27eH/hs3blw/lgYAAAarPgWQuro6paWlyWazyWq1Kj8/XzU1NaHtfr9fDQ0Nys3NlSTl5+drz5498vv9euedd/Tee+9p1apV3drcv3+/AoGAbrzxRhUUFOhPf/pTP5YFAAAGsz59G67H41FiYmJoOTExUW63O7Tc0tKi2NhYWa3WE41aLBo9erR8Pp8mTZqkBx98UA0NDd3a9Pv9mjZtmpYsWaKmpibddNNNmjhxoiZOnNgfdQH9or+/TRQAcEKfAsjJH518Lioq6ozbJSk6+vQnWHJycpSTkyNJSklJUVZWlurq6noMII2NjWdsa7hpbW2Vy+UKdzfCarCMgc8XG5bjdnR0yOfzheXYg0FP9R8/HiUpQZLU0nJMhmGEoWfmGYxzwOVqN+1Yg+U9IJyG8hicLhecrE8BJCkpSXv37g0tezweJSUlhZZtNpva2toUCARksVgUCATU3t6u+Pj407a5Y8cOffOb39SkSZMkSYZhKCYmpsd9k5OTT7ttOHK5XEpJSQl3N8JqsIyBzRae4/p8PtnCdfBBoKf6T/qbR/HxCUpIMLlTJhuMcyAlxbz+DJb3gHAaymMQDAbl9XrPuE+fTiukp6ervr5eXq9Xfr9flZWVysjICG23Wq1KTU1VVVWVJKmqqkqpqamhj2R68sEHH+ipp56SYRjyeDz685//3K1NAAAwfPUpgNjtdi1dulQLFy5Ubm6uLrzwQmVlZamoqEh//vOfJUnFxcXavn27cnJy9Nxzz+n+++8/Y5sLFy7UyJEjlZOTo5tvvln33nuvJkyY8PUrAgAAg16fPoKRJIfDIYfD0W3d57faStK4ceO0adOm075+ypQpmjJlSmh5xIgR+vWvf/0lugoAAIaLyLmyEwAADBoEEAAAYDoCCAAAMB0BBAAAmI4AAgAATEcAAQAApiOAAAAA0xFAAACA6QggAADAdAQQAABgOgIIAAAwHQEEAACYjgACAAC6FHGNAAARaElEQVRMRwABAACmI4AAAADTEUAAAIDpCCAAAMB0BBAAAGA6AggAADAdAQQAAJiOAAIAAExHAAEAAKYjgAAAANMRQAAAgOkIIAAAwHQEEAAAYDoCCAAAMB0BBAAAmK7PAcTpdConJ0fZ2dlav379KdvdbrcKCwvlcDi0YMECHT16tNv21157TTfffHO3devWrdO1116rmTNn6qWXXvqKJQAAgKGmTwGkublZa9euVXl5uaqrq7Vv3z7V1tZ226ekpEQFBQVyOp3Ky8tTaWmpJCkYDOp3v/udFi9erK6urtD+u3bt0v79+/Xiiy+qvLxcpaWlOn78eD+WBgAABqs+BZC6ujqlpaXJZrPJarUqPz9fNTU1oe1+v18NDQ3Kzc2VJOXn52vPnj3y+/1655139N5772nVqlXd2nzppZeUl5cni8WixMREpaamavfu3f1YGgAAGKwsfdnJ4/EoMTExtJyYmCi32x1abmlpUWxsrKxW64lGLRaNHj1aPp9PkyZN0oMPPqiGhoZubbrdbtnt9m5tHjlypMfjNzY2Kjo6ci5XaW1tlcvlCnc3wmqwjIHPFxuW43Z0dMjn84Xl2INBT/UfPx4lKUGS1NJyTIZhhKFn5hmMc8DlajftWIPlPSCchvIYnPyJx+n0KYD01FBUVFSvBzpTaOjpzeN0+ycnJysmJqa3bg4bLpdLKSkp4e5GWH3VMdiwoX/7YbP1b3t95fP5ZAvXwQeBnuo/6S1H8fEJSkgwuVMmG4xzwOnsv/7ccceZt/M+OLTHIBgMyuv1nnGfPp1WSEpKUnNzc2jZ4/EoKSkptGyz2dTW1qZAICBJCgQCam9vV3x8/GnbtNvt8ng8oeXm5uZuZ0QAAMDw1acAkp6ervr6enm9Xvn9flVWViojIyO03Wq1KjU1VVVVVZKkqqoqpaamhj6S6UlGRoYqKysVCATk9XpVX1+vq6666utVAwAAhoQ+fQRjt9u1dOlSLVy4UJ2dncrMzFRWVpaKioqUmZmpGTNmqLi4WCtWrNDGjRsVFxensrKyM7aZnZ2tN998U3PmzFEwGNTixYs1duzYfikKAAAMbn0KIJLkcDjkcDi6rfv8VltJGjdunDZt2nTa10+ZMkVTpkwJLUdFRWnJkiVasmTJl+guAAAYDiLn1hIAADBoEEAAAIDpCCAAAMB0BBAAAGA6AggAADAdAQQAAJiOAAIAAExHAAEAAKYjgAAAANMRQAAAgOkIIAAAwHQEEAAAYDoCCAAAMB0BBAAAmI4AAgAATEcAAQAApiOAAAAA0xFAAACA6QggAADAdAQQAABgOgIIAAAwHQEEAACYzhLuDgAAIs+GDWfe7vPFymbre3t33PH1+gPzcQYEAACYjgACAABMRwABAACmI4AAAADT9ekiVKfTqfXr18vv9ysvL0933313t+1ut1v33HOPvF6vxo4dq0ceeUTnnnuu/H6/iouL9frrrysqKkqrV6/WJZdcIkmaP3++2tvbFRMTI0lauXJlaBsAABjeeg0gzc3NWrt2rSoqKjRmzBjdfvvtqq2t1dSpU0P7lJSUqKCgQHPnztXzzz+v0tJSrVu3Ts8884wMw1BNTY3effdd3XXXXaqurpZhGDp8+LBeeeUVRUdzEgYAgEjT62//uro6paWlyWazyWq1Kj8/XzU1NaHtfr9fDQ0Nys3NlSTl5+drz5498vv92r17t+bOnStJOv/882W32/W3v/1NBw8elMVi0Y9//GPl5eVp8+bNA1QeAAAYjHo9A+LxeJSYmBhaTkxMlNvtDi23tLQoNjZWVqv1RIMWi0aPHi2fzye32y273d7ttUeOHNHYsWOVnp6u4uJiffLJJ1qwYIHOO+88paen99iHxsbGiDpT0traKpfLFe5umOaZZ2JPWdfREaORI31h6M3g0NHRIZ+P+k92/HiUpARJUkvLMRmGEYaemYc58OXqd7naB7A34TGUfxd0dXX1uk+vAaSnRqKiono9SHR0dI9vENHR0UpLS1NaWpokaeTIkZo3b5727Nlz2gCSnJwculYkErhcLqWkpIS7G6bp6WFDPp9Pti/zFKJhhvpPrf+ktx3FxycoIcHkTpmMOfDl6k9JGX5jNZR/FwSDQXm93jPu0+tphaSkJDU3N4eWPR6PkpKSQss2m01tbW0KBAKSpEAgoPb2dsXHx8tut3d7bXNzs+x2u1577TXt27cvtL6rqyuiAgYAAJGu1wCSnp6u+vp6eb1e+f1+VVZWKiMjI7TdarUqNTVVVVVVkqSqqiqlpqbKarUqIyNDFRUVkqRDhw7pww8/1He/+135fD6VlZWps7NTbW1t2rZtm7KysgamQgAAMOj0+hGM3W7X0qVLtXDhQnV2diozM1NZWVkqKipSZmamZsyYoeLiYq1YsUIbN25UXFycysrKJEk33XSTSkpKlJOTE7oNd8SIEcrJydFbb72lOXPmqKurSz/60Y902WWXDXixAABgcOjTc0AcDoccDke3daWlpaGfx40bp02bNp3yuhEjRnTb73NRUVFatmyZli1b9iW7CwAAhoPIubUEAAAMGgQQAABgOgIIAAAwHQEEAACYjgACAABMRwABAACmI4AAAADT9ek5IAAADGYbNvRve3fc0b/t4VScAQEAAKYjgAAAANMRQAAAgOkIIAAAwHRchIovrb8v9gIARB7OgAAAANMRQAAAgOkIIAAAwHQEEAAAYDoCCAAAMB0BBAAAmI4AAgAATEcAAQAApuNBZBGAB4cBwJfTn++bfLNuzzgDAgAATEcAAQAApiOAAAAA03ENyCDFdRsAgOGMAAIAwAD6qn9Q+nyxstlOXT9cLmrt00cwTqdTOTk5ys7O1vr160/Z7na7VVhYKIfDoQULFujo0aOSJL/fr/vuu0+zZs1STk6O/vd//zf0mvLycjkcDmVnZ+vZZ5/tp3IAAMBQ0OsZkObmZq1du1YVFRUaM2aMbr/9dtXW1mrq1KmhfUpKSlRQUKC5c+fq+eefV2lpqdatW6dnnnlGhmGopqZG7777ru666y5VV1fr4MGD2rJliyoqKtTV1aUf/vCHuuKKKzRx4sQBLRYAgKFuuNwi3GsAqaurU1pammz/Og+Un5+vmpqaUADx+/1qaGjQY489Ftq+evVq+f1+7d69W3fddZck6fzzz5fdbtff/vY37d27VzNnztSoUaMkSTNnzpTT6dTdd9/d7diGYUiSgsFgP5Xb3f/7fwPS7Nd27NjZSkgYmJqHipiYLkVHR+4YUP+p9cfESLGx//dz9DC/hJ45ENn1S+aMwQD9eg393v7893hPeg0gHo9HiYmJoeXExES53e7QcktLi2JjY2W1Wk80aLFo9OjR8vl8crvdstvt3V575MgRud1uTZ48ObTebrfrjTfeOOXYXV1dkqT9+/f31s2v5LLLBqTZfvJBuDswCET6GFD/F02fHoZuhBVzAAM7Bq+/PqDNh36P96TXANLTi6OionptPDo6usfkc7r1J7cZ6pzFosmTJys6OrrH7QAAYPAxDENdXV2yWE4fM3oNIElJSdq7d29o2ePxKCkpKbRss9nU1tamQCAgi8WiQCCg9vZ2xcfHy263q7m5Wd/4xjcknbiexG63y263y+PxnLbNz0VHR2vEiBF9qxYAAAwZvX6Kmp6ervr6enm9Xvn9flVWViojIyO03Wq1KjU1VVVVVZKkqqoqpaamymq1KiMjQxUVFZKkQ4cO6cMPP9R3v/tdTZs2TTt37lR7e7va29u1Y8cOTZs2bWAqBAAAg06UcaYrRP7F6XTqySefVGdnpzIzM7Vs2TIVFRUpMzNTM2bMUFNTk1asWKHm5mbFxcWprKxM48ePV2dnp0pKSvT6668rKipKRUVFSk9Pl3TiNtz/+q//UiAQ0A033KCFCxcOeLEAAGCQMBBWL7zwgjFr1ixj1qxZxpo1awzDMIy3337buP76642ZM2caixYtMtrb28Pcy4H17//+70Z2draRm5trPPnkk4ZhRMYYtLa2Grm5uYbL5TIMwzAaGhqMvLw8Izs723jggQcMv99vGIZhHDlyxLjpppuMa6+91igsLDS8Xm84u91vvli/YRhGZ2encfPNNxv19fWhdcN5LnxxDCorK43c3FwjNzfXuPPOO42WlhbDMCJnDmzatMlwOByGw+Ew1qxZY3R1dRmGEVlz4HN//OMfjZtuuim0PBznAAEkjD755BPjyiuvNLxer+H3+4158+YZdXV1Rl5entHQ0GAYhmE8+uijRllZWZh7OnDq6uqMnJwc4+OPPzYCgYDxk5/8xNi5c+ewH4O//vWvRm5urnHxxRcbLpfL+Oyzz4yMjAzjn//8p9HV1WXce++9xn/+538ahmEYP/vZz4ytW7cahmEYzz33nLF48eJwdr1ffLF+wzCMd99915g/f74xefLkbgFkuM6FL47BkSNHjKlTpxpHjx41DMMw1q1bZ6xatcowjMiYAwcPHjSys7ONTz75xAgEAsb8+fON2tpawzAiZw587p133jGmTp3aLYAMxzkwzO+kH9yCwaC6urrU0dGhYDCoYDAoi8Wi1tZWpaamSpKuv/561dTUhLmnA+fAgQOaOnWqxowZo5iYGE2dOlXl5eXDfgy2bNmi4uLi0C3ub7zxhpKTk3XeeecpKipK8+bNU01NTeg5O7m5uZJOPGdnz5498vv94ez+1/bF+iXp+eef12233aZLLrkktK6pqWnYzoUvjkFUVJRWrlwZeubSt7/9bTU1NUXMHLjgggv04osv6uyzz9bHH3+strY2nXPOORE1BySps7NTDzzwgH7xi1+E1g3XOUAACaPRo0frF7/4hRwOh66++mqNHz9eVqv1lGennPzcleHm4osv1quvvqqWlhZ99tlneumllyJiDB566CFdccUVoeWenpnjdrvP+JydoeyL9UvSsmXLdM0113Rbd7pxGQ6+OAaJiYmhC/w//fRTbdiwQddcc01EzQGr1arNmzdrxowZGjt2rCZNmhRRc0CSfvOb3+i6665TcnJyaN1wnQMEkDD6xz/+oYqKCu3evVu1tbWKjo5WXV3dKfsN52egpKenq6CgQIWFhbrtttt0+eWXKxAInLLfcB4D6fTP2znTc3YiQW/PIRqOjh07pltvvVUXXXSR5s6dG3Fz4MYbb9TevXtls9m0fv36iJoDdXV1ampq0nXXXddt/XCdA0O790Pcq6++qvT0dJ177rkaMWKECgoK1NDQoObm5tA+zc3NPT4jZbhoa2tTVlaWqqqqVF5errPPPlsTJkyIqDGQTjxvp6eaT37OjqRuz9mJBKcbl+Hq8OHD+uEPf6jLLrtMq1atkqSImQONjY16/V+P5bRYLJo9e7befvvtiJoDL774ot555x3NmTNHv/zlL/Xmm2/q5z//+bCdAwSQMJo0aZJeffVVtbe3yzAMvfTSS7riiit09tlnhx7+VlFRMayfkXL48GH97Gc/k9/vV2trq5577jnNmzcvosZAki655BK9//77OnTokKQTNWdkZJzxOTuRYPz48REzFz777DPdeuutuuGGG7R06dLQX/mRMgeOHTumpUuXqq2tTV1dXXI6nbriiisiag489NBDcjqd2r59ux588EF95zvf0eOPPz5s50CvT0LFwPn+97+vAwcOqKCgQCNGjNB3vvMd3XHHHXI4HLr//vvV2tqq5ORk/eY3vwl3VwfMhRdeqNzcXM2ZM0fBYFC33HKLLr/8cv3mN7+JmDGQpBEjRujhhx/WkiVL1NHRoe9+97u68cYbJUnFxcVasWKFNm7cGHrOTiSJlLmwfft2ffjhh9q2bZu2bdsmSbrooov00EMPRcQcmDx5shYsWKD58+crOjpaqampuuWWWyRFzhw4k+E4B/r0IDIAAID+xEcwAADAdAQQAABgOgIIAAAwHQEEAACYjgACAABMRwABAACmI4AAAADTEUAAAIDp/j9P2w9EwbG3KgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(8, 5), dpi=80)\n",
    "ax.hist(boys_ppc['first'] + boys_ppc['second'],\n",
    "        bins=23, density=True, alpha=.4, color='b', label='sampled')\n",
    "ax.axvline(x=n_boys, color='b', label='observed')\n",
    "ax.legend()\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**H4** Now compare 10,000 counts of boys from 100 simulated first borns only to the number of boys\n",
    "in the first births, birth1. How does the model look in this light?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhkAAAFLCAYAAAB7tJ7qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAMTQAADE0B0s6tTgAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X9YlHW+//HXAGMmuMKYgAbZZq5urfYbwWoFTAtBRNNq2/TEZVZXWacfUhkaakEdssyiTsfjKcvTj12TFALWsnQjPGJutWm2/bBMWJBfAySkMcPM9w+/zkqoYM0HBnk+rqvrmns+93zu982ne+bl577nHovb7XYLAADAy/y6uwAAAHByImQAAAAjCBkAAMAIQgYAADCCkAEAAIwI6O4CJMnlcsnpdMrPz08Wi6W7ywEAAJ3gdrvlcrkUEBAgP7/28xY+ETKcTqd27NjR3WUAAICfYdSoUerTp0+7530iZBxOP6NGjVJlZaUiIiK6uSL8VHl5OePiozoam3/+UzrnnEOPd+2STj+9iwrr5ThmfBdj4z2tra3asWPHUWcxJB8JGYdPkfj7+8vPz0/+/v7dXBF+inHxXR2NjZ+f1Nz8r8cMY9fgmPFdjI33HetSBy78BAAARhAyAACAEYQMAABgBCEDAAAYQcgAAABGEDIAAIARhAwAAGAEIQMAABhByAAAnDTKy8sVHx/f3WW0cayampqaNH36dE2ZMkUvvfSSli9f3uk+33vvPb344ouSpGeeeUbPPPOM1+r1Jp+44ycAoGdxOqV9+8z0HR4uBfSCT6fPP/9c/v7+euONN074tZ999pmBiryvFwwjAMCbnE7pd7+TvvjCTP8jRkg7d3YcNJ5//nnl5eXJ399fl156qdLS0iRJP/74o+666y598803ioyMVFZWlgYMGKCcnBxt2LBBLpdLF110kZYsWaLW1lYtXbpUW7duldPp1FVXXaXbb79dpaWlys7Oltvt1plnnqlt27Zp7dq1CgsL0w8//KArrrhC77zzjr799ltlZWXpwIED6t+/vzIyMjRs2DDt2rVL6enpkqSRI0e2q72urk4PPvigamtrNWfOHCUkJGjbtm167LHHFB8fr1GjRukf//iHVq1apUcffVR79+6VxWLRtddeqwsuuECvv/66JCk8PFyS9Omnn+qaa65RU1OTrrnmGt14443H/BtVVlZq9uzZGjhwoCwWi66++mq9//77ampqUllZmX77298qOzv7qD94dqIIGQCMW7HCu/3dfLN3+0PP89e//lUbN27U2rVrZbVadccdd+j111/XuHHjVFdXpxtuuEEXX3yxsrOz9eyzz+qOO+7Qq6++quLiYpWXl+uFF15QZWWl3n//fTkcDuXm5srpdOq2227Tpk2b1K9fP3377bfatGmTBgwYoKysLBUVFenGG2/Uxo0bdemll6pPnz568MEH9dxzzykiIkJ/+9vflJaWptzcXN1///267777dPnll+vZZ59VaWlpm/oHDhyoRx55RDk5Ofrv//5v5ebmtmm/7LLLtHz5cm3cuFFut1vr1q1TfX29srKydN111+m6666TJM2YMUPPPPOMqqur9dprr8nlcmnatGkaM2aMqqurj/k32rNnj1asWKGhQ4cqNzdXH3/8sd566y3169dPM2bM0AcffOCV006EDADACQkIODTT0J2nS7Zu3aqkpCSdeuqpkqSrr75a69at07hx4zR06FBdfPHFkqTk5GQ98MAD6t+/v4YNG6bp06froosu0g033KDBgwerpKREu3bt0ocffihJOnDggL788kudf/75OuusszRgwABJUkpKihYtWqQbb7xRb731lmbNmqVvv/1We/fu1e233+6py263q66uTlVVVbr88sslSdOmTdPatWtP6G9wwQUXSDr06+SZmZmaPXu2xo0bp/vvv/+o6yckJKhfv36SpLi4OG3btk379u075t/IZrNp6NChntdfeOGF6t+/vyRp+PDhamxsPKF6j4WQAQA4YQEBUnf+WrrL5Wr3nNPplKR2Pzse8P8Ty8svv6yPPvpIBQUFmj17tpYuXarW1lbdd999mjhxoiSpvr5effv21aeffur5cJakc845RwcOHNDOnTv19ddfKyYmRl999ZUiIyO1fv16SZLb7VZVVZX8/Pzkdrvbbf9E9O3bV5IUFhamoqIilZSUqLi4WFOnTlVBQUG79Y/cZ7fbrYCAgOP+jY7cN0k65ZRTPI8tFkub+n8Jvl0CAOhxoqOj9dZbb+nAgQNyOp1au3atLrnkEknSnj17tHPnTknSG2+8obFjx6q8vFzJyck699xzlZqaqksvvVRffPGFoqOj9ec//1kOh0MHDhzQjTfeqJKSkqNuc8qUKXrooYd05ZVXyt/fX2eddZYaGxs9syD5+fm69dZbFRISotNPP10bN26UpKOGgs7Kz8/XokWLNH78eC1YsED9+vVTZWWl/P39PYFBkt5++221tLSosbFRmzZtUnR09HH/Rl2FmQwAQI8TFxenzz//XNOnT5fT6dTYsWM1a9YsVVVV6YwzztB//dd/ac+ePRo+fLjuvvtuBQYGatKkSZo6dar8/f115plnaurUqTrllFP03XffaerUqXI4HEpMTNQVV1zR7hoKSZo8ebKefPJJZWVlSZL69Omj5cuXKysrSwcPHlS/fv20dOlSSdLjjz+u+fPnKycnR+eff/7P3s+JEyfq3XffVWJioqxWq6688kqNGDFCTU1NSktLU0hIiCTp9NNP1x/+8AcdOHBAN998s4YNG6Zhw4Yd82/UVSxub82J/AKtra365JNPdP7556uiokKRkZHdXRJ+oqysjHHxUR2NTXm5dLi5rKx7prh744WfHDO+i7HxniM/v/39/du1c7oEAAAYQcgAAABGcE0GgB7H26dfpJ5xCgboaZjJAAAARhAyAACAEYQMAABgBNdkAAB+FhPXxkhcH3MyYSYDAIDjKC8vP+EfC4uPj1d5eXm75x988EFNnDhR69at05w5czrdX1lZmR588EFJUmlpqWbOnHlC9XQXZjIAAOgiubm5+vTTT9WnTx+lpKR0+nUVFRUqKyszWJkZhAwAQI9jt9t17733en4t9LbbbtMVV1yhbdu26cknn9SPP/6oxsZG3XvvvUpMTNQDDzygvn376pNPPlFtba3S0tK0ceNGff7554qLi1N6erpyc3O1YcMGNTU1qbq6WuPGjVN6enq77T700EOqqKiQJM2dO1fx8fFqaGhQWlqaKisrdfbZZ+vHH39sV/Ott94qt9utGTNmKCMjQ/PmzdN7772nBx54QPX19dq7d6/+/d//XV9//bU2bNiggIAAjRo1SkuWLNGSJUtUUVGhhx56SImJiaqvr9dNN92kffv2afTo0Vq0aJH69OmjTZs26amnnpLL5VJkZKSWLFmi0047TfHx8Ro1apT+8Y9/KDMzU0uWLNG5556rzz77TH379tUTTzxh5C6onC4BAPQ4+fn5GjFihHJzc/X44497fqTsf//3f/Xwww/rzTffVFZWlp5//nnPa6qqqrRu3Tqlpqbq4Ycf1qJFi7Ru3TqtXbtW33//vSTp008/1VNPPaWCggLt3LlTf/nLX9psNzMzU8nJycrNzdXKlSv16KOPqq6uTk8//bRGjhypt956SzfccINqa2vb1Xy4lvXr1ys0NLRNW//+/VVUVKRLL71Ur776qqcui8WiyspKPfTQQ/rd736nJUuWSDp0+mTBggXKz8/XDz/8oNdee011dXVauHChcnJylJ+frwsvvNCzviRddtll2rBhg8LDw/Xll1/q+uuvV15enkaPHq1XXnnFC6PSHjMZAIAeZ8yYMZo9e7bKy8v1+9//XnfccYekQz9MtmnTJr399tv6+9//rubmZs9rYmNjJUmhoaEaPny4Bg4cKEkKDg72hIy4uDgNGjRIkjRp0iR9+OGHGjVqlKePDz74QF999ZWeffZZSYd+On337t3atm2b58fRLr744hOeFbjgggskHQobw4YN0/Tp0xUXF6cbbrhBgwcP1t69e9usf9FFF+nMM8+UdOiH23Jzc3XGGWdo9OjRnm1fe+21WnHE1bmHtyFJAwcO9OzXb3/7W23fvv2E6u0sQgYAoMcZOXKkNmzYoPfff1+bNm3SCy+8oKKiIl1//fWKiopSdHS0YmJiNG/ePM9rrFar53FAwNE//o583uVyyc+v7YS/y+XSyy+/rODgYElSdXW1bDabLBbLMfvpjFNPPdXz+OWXX9ZHH32k4uJizZ492xNejnRkXW63WwEBAXK5XG3WcbvdcjgcnuW+fft6Hp9yyimexxaLRaZ+K5XTJQCAHmfFihVauXKlJk2apEWLFslut6uxsVF79uzRXXfdpXHjxqmkpEStra0n1G9xcbG+//57/fjjjyooKNBll13Wpj06OlqvvvqqJGnPnj1KSkpSY2OjYmJi9Oabb0qSdu7cqe++++5n7Vd5ebmSk5N17rnn6q677tKll16qL774Qv7+/nI6nZ71Pv74Y/3zn/+Uy+XSunXrNHbsWJ133nn69NNPPReI/ulPf1JUVNTPqsNbOh21ioqKlJOTI4fDoeTkZM2dO7dNe1VVlebNm6fa2loNGjRIy5Yt08CBA9XS0qLs7Gx9+OGHcjqdmj9/frtBAwD0PN15P4vp06dr3rx5mjx5sgICAnTHHXcoODhYM2bMUGJiooKCgnTeeefp4MGDbU6ZdGTQoEG65ZZbZLfblZSUpNjY2DZfRV2wYIEyMjI0efJkud1uZWZmauDAgbrzzjv1wAMPKDExUWedddbPvogyIiJCkyZN0tSpU3XqqadqyJAhmjp1qpxOp5qamnTPPffo2muv1fDhw7Vw4UJVVVVpzJgxmj59uvz9/bVkyRLNnTtXTqdT4eHhysrK+ll1eIvF3Yk5kpqaGl1zzTVau3at+vfvrzlz5mj27Nm6/PLLPevcdtttmjBhgqZOnao33nhDW7Zs0ZNPPqnnnntOu3fv1tKlS/X1118rNTVV77//fpupniN/j76iosLIFa74ZcrKyhgXH9XR2JSXS4eby8qkiIguKuwIpm7a5E3e/sDkmPFdxxqb3Nxcbdu2TY899lg3VNUzHfn57e/v3669U6dLSkpKFB0dLZvNJqvVqpSUFBUWFnraHQ6HSktLlZSUJElKSUnR5s2b5XA4VFRUpDlz5shisWj48OFatWqVsXM/AADAd3TqdEl1dXWbr9uEhoaqqqrKs9zQ0KDAwEDPRTUBAQEKCgqS3W7Xd999p7/97W9KT0+Xv7+/7rnnHp199tlH3U55ebmam5t75A1HTnb79+9nXHxUR2NTWekvaYikQzf0cbtP7By1N9jtgV2+zRNVVtb5KfXO4JjxXccam0suuUSXXHIJ43YCfnqx6U91KmQcrZMjr6Q91kb8/PzU2tqq8vJyvfHGG/riiy900003qaioSP3792+3fkREBKdLfBRTv76ro7E58qL3IUOGdMvpEput67d5oiIjvVskx4zvYmy8p7W19aj3BDmsUyEjPDxc27Zt8yxXV1crPDzcs2yz2dTU1CSn06mAgAA5nU41NzcrODhYp512mhISEmSxWDRy5EiFh4fr22+/1ejRo3/BbgGAd3n7upGEBO/2B/REnbomIyYmRlu3blVtba0cDofy8vI8NzWRDn33OCoqSvn5+ZIO3YktKipKVqtVcXFxKioqknTodEhlZaV+/etfe39PAACAT+lUyAgLC1NaWppSU1OVlJSkESNGaMKECUpPT9e7774rScrIyND69euVmJioNWvWaOHChZKkefPmqa6uTomJibrlllu0ZMmSo54qAQAAJ5dO3ycjISFBCT+Z/8vMzPQ8Hjx4sFatWtXudUFBQcrOzv75FQIAgB6JO34CAAAjCBkAAMAIQgYAADCCX2EF0EZPuAU4gJ6BmQwAAGAEIQMAABhByAAAAEYQMgAAgBGEDAAAYAQhAwAAGEHIAAAARhAyAACAEYQMAABgBCEDAAAYQcgAAABGEDIAAIARhAwAAGAEIQMAABhByAAAAEYQMgAAgBGEDAAAYAQhAwAAGEHIAAAARhAyAACAEYQMAABgBCEDAAAYQcgAAABGEDIAAIARhAwAAGAEIQMAABhByAAAAEYQMgAAgBGEDAAAYAQhAwAAGEHIAAAARgR0ZqWioiLl5OTI4XAoOTlZc+fObdNeVVWlefPmqba2VoMGDdKyZcs0cOBAHThwQGPHjtUZZ5zhWTc3N1f+/v7e3QsAAOBzOpzJqKmpUXZ2tlavXq2CggJt375dxcXFbdZZvHixpk2bpqKiIiUnJyszM1OS9Nlnnyk6Olrr16/3/EfAAACgd+gwZJSUlCg6Olo2m01Wq1UpKSkqLCz0tDscDpWWliopKUmSlJKSos2bN8vhcGjHjh2qqqrSjBkzdN1112n79u3m9gQAAPiUDk+XVFdXKzQ01LMcGhqqqqoqz3JDQ4MCAwNltVoPdRgQoKCgINntdlksFl111VWaM2eOdu3apVtuuUX5+fkKCQk56rbKy8vV3NyssrKyX7pf8LL9+/czLj6qo7GprPSXNESSVFFRIbe79bj92e2B3iyv1+KY8V2Mjfe4XK7jtncYMo7WgcVi6XADfn5+uvHGGz3L5557rkaNGqWPPvpI48ePP+prIiIiVFFRocjIyI7KQhcrKytjXHxUR2NzxOGqIUOGKCLi+P3ZbF4qrJfr37+ZY8ZH8X7mPa2traqtrT1me4enS8LDw1VTU+NZrq6uVnh4uGfZZrOpqalJTqdTkuR0OtXc3Kzg4GCtWbNGlZWVnnXdbrcCAjp1rSkAAOjhOgwZMTEx2rp1q2pra+VwOJSXl6fY2FhPu9VqVVRUlPLz8yVJ+fn5ioqKktVq1Y4dO/Tyyy9Lkr7++mvt2rVLF110kZk9AQAAPqXDaYWwsDClpaUpNTVVLS0tio+P14QJE5Senq74+HiNHz9eGRkZmj9/vlauXKkBAwZo6dKlkqS7775b8+fPV2Jiovz8/JSdna2goCDjOwUAALpfp85dJCQkKCEhoc1zh7+mKkmDBw/WqlWr2r0uJCREzz///C+rEAAA9Ejc8RMAABhByAAAAEYQMgAAgBGEDAAAYAQhAwAAGEHIAAAARnD7TQAw4JVXAr1+i/abb/Zuf4BpzGQAAAAjmMkAeriO/sVcX3/kutIxfgQZALyOmQwAAGAEIQMAABhByAAAAEYQMgAAgBGEDAAAYAQhAwAAGEHIAAAARhAyAACAEYQMAABgBCEDAAAYQcgAAABGEDIAAIARhAwAAGAEIQMAABhByAAAAEYQMgAAgBGEDAAAYAQhAwAAGEHIAAAARhAyAACAEYQMAABgBCEDAAAYQcgAAABGEDIAAIARhAwAAGAEIQMAABjRqZBRVFSkxMRETZw4UTk5Oe3aq6qqNHPmTCUkJGjWrFmqq6tr097U1KQrrrhCpaWl3qkaAAD4vA5DRk1NjbKzs7V69WoVFBRo+/btKi4ubrPO4sWLNW3aNBUVFSk5OVmZmZlt2h9++GF9//333q0cAAD4tA5DRklJiaKjo2Wz2WS1WpWSkqLCwkJPu8PhUGlpqZKSkiRJKSkp2rx5sxwOhySpsLBQgYGBGjFihKFdAAAAviigoxWqq6sVGhrqWQ4NDVVVVZVnuaGhQYGBgbJarYc6DAhQUFCQ7Ha7Wltb9dJLL+mll17SnDlzOiymvLxczc3NKisr+zn7AoP279/PuPiogwf9Zbfbj9ne2GiRFCJJamiol9vt7qLKereDBw8ed1x+jrKyZq/211vxfuY9LpfruO0dhoyjdWCxWDq1gfT0dC1cuFB9+/btaDOSpIiICFVUVCgyMrJT66PrlJWVMS4+qm9fu2w22zHbjzhcFRwcopCQLigKstuPPy4/R2Skd/vrrXg/857W1lbV1tYes73D0yXh4eGqqanxLFdXVys8PNyzbLPZ1NTUJKfTKUlyOp1qbm5WfX29vvnmG6Wnp2vKlCnauXOnFixYoC1btvyS/QEAAD1EhyEjJiZGW7duVW1trRwOh/Ly8hQbG+tpt1qtioqKUn5+viQpPz9fUVFRGjlypP76179q/fr1Wr9+vX73u9/pkUce0dixY43tDAAA8B0dhoywsDClpaUpNTVVSUlJGjFihCZMmKD09HS9++67kqSMjAytX79eiYmJWrNmjRYuXGi8cAAA4Ns6vCZDkhISEpSQkNDmuSO/pjp48GCtWrXquH2sXr36xKsDAAA9Fnf8BAAARhAyAACAEYQMAABgBCEDAAAYQcgAAABGEDIAAIARhAwAAGAEIQMAABhByAAAAEYQMgAAgBGEDAAAYAQhAwAAGEHIAAAARhAyAACAEYQMAABgBCEDAAAYQcgAAABGEDIAAIARhAwAAGBEQHcXAPQmK1Z0dwUA0HWYyQAAAEYwkwEAPYS3Z8Juvtm7/QE/xUwGAAAwgpABAACMIGQAAAAjCBkAAMAIQgYAADCCkAEAAIwgZAAAACMIGQAAwAhCBgAAMIKQAQAAjCBkAAAAIwgZAADACEIGAAAwotMho6ioSImJiZo4caJycnLatVdVVWnmzJlKSEjQrFmzVFdXJ0mqqKjQrFmzlJycrBkzZujzzz/3XvUAAMBndSpk1NTUKDs7W6tXr1ZBQYG2b9+u4uLiNussXrxY06ZNU1FRkZKTk5WZmSlJeuyxxzR58mTl5eXpjjvu0OLFi72/FwAAwOd0KmSUlJQoOjpaNptNVqtVKSkpKiws9LQ7HA6VlpYqKSlJkpSSkqLNmzfL4XDoqaee0tVXXy1JKi8v169+9SsDuwEAAHxNQGdWqq6uVmhoqGc5NDRUVVVVnuWGhgYFBgbKarUe6jQgQEFBQbLb7QoLC5MkTZw4URUVFfrP//zPY26nvLxczc3NKisr+1k7A3P279/PuHiB3R7o9T4PHjwou91+zPbGRoukEElSQ0O93G6312tAex2Niy8oK2vu7hK6Be9n3uNyuY7b3qmQcbROLBZLhxvx8/vXRMnbb7+tzz77TLNnz9Zf/vIXBQcHt1s/IiJCFRUVioyM7ExZ6EJlZWWMixfYbN7v0263y3acjo84VBUcHKKQEO/XgPY6GhdfEBnp2/WZwvuZ97S2tqq2tvaY7Z06XRIeHq6amhrPcnV1tcLDwz3LNptNTU1NcjqdkiSn06nm5mYFBwdr48aNOnjwoCTp3HPP1emnn06CBACgF+hUyIiJidHWrVtVW1srh8OhvLw8xcbGetqtVquioqKUn58vScrPz1dUVJSsVqvWrFmjtWvXSpK+/PJL1dXVadiwYd7fEwAA4FM6dbokLCxMaWlpSk1NVUtLi+Lj4zVhwgSlp6crPj5e48ePV0ZGhubPn6+VK1dqwIABWrp0qSQpIyNDDz74oP785z/rlFNO0ZNPPql+/foZ3SkAAND9OhUyJCkhIUEJCQltnjv8NVVJGjx4sFatWtXudUOGDDnq8wAA4OTGHT8BAIARhAwAAGAEIQMAABhByAAAAEYQMgAAgBGEDAAAYAQhAwAAGEHIAAAARhAyAACAEYQMAABgBCEDAAAYQcgAAABGEDIAAIARhAwAAGAEIQMAABhByAAAAEYQMgAAgBGEDAAAYAQhAwAAGEHIAAAARhAyAACAEYQMAABgBCEDAAAYQcgAAABGEDIAAIARhAwAAGAEIQMAABgR0N0FAL5sxYrurgAAei5mMgAAgBHMZABAL2Vipu7mm73fJ3ouZjIAAIARhAwAAGAEIQMAABhByAAAAEYQMgAAgBGEDAAAYESnQkZRUZESExM1ceJE5eTktGuvqqrSzJkzlZCQoFmzZqmurk6SVF9fr9tuu03JycmaPHmyCgoKvFs9AADwWR2GjJqaGmVnZ2v16tUqKCjQ9u3bVVxc3GadxYsXa9q0aSoqKlJycrIyMzMlSU8//bTOOecc5eXl6YUXXtCjjz6q2tpaM3sCAAB8Socho6SkRNHR0bLZbLJarUpJSVFhYaGn3eFwqLS0VElJSZKklJQUbd68WQ6HQ7///e917bXXSpIGDRqk4OBgQgYAAL1Eh3f8rK6uVmhoqGc5NDRUVVVVnuWGhgYFBgbKarUe6jAgQEFBQbLb7YqLi/OsV1BQoJaWFp199tnH3FZ5ebmam5tVVlb2s3YG5uzfv79XjovdHtjdJXTo4MGDstvtx2xvbLRICpEkNTTUy+12d1FlvVtH43KyKitr7u4SOtRb389McLlcx23vMGQcrQOLxdLhBvz8/jVJsn79ej3++ONauXKlAgKOvcmIiAhVVFQoMjKyo7LQxcrKynrluNhs3V1Bx+x2u2zHKfSIw1XBwSEKCemCotDhuJysIiN9f5976/uZCa2trcc9Q9Hh6ZLw8HDV1NR4lqurqxUeHu5ZttlsampqktPplCQ5nU41NzcrODhYkrRixQotX75cL730kkaOHPmzdwQAAPQsHYaMmJgYbd26VbW1tXI4HMrLy1NsbKyn3Wq1KioqSvn5+ZKk/Px8RUVFyWq1Kjc3V2+++ab+9Kc/adiwYcZ2AgAA+J4OT5eEhYUpLS1NqampamlpUXx8vCZMmKD09HTFx8dr/PjxysjI0Pz587Vy5UoNGDBAS5culSQtW7ZMFotFN910k6e/JUuW6LzzzjO3RwAAwCd06qfeExISlJCQ0Oa5w19TlaTBgwdr1apV7V7306+6AgCA3oM7fgIAACMIGQAAwAhCBgAAMIKQAQAAjCBkAAAAIwgZAADACEIGAAAwgpABAACMIGQAAAAjCBkAAMAIQgYAADCCkAEAAIwgZAAAACMIGQAAwAhCBgAAMIKQAQAAjCBkAAAAIwgZAADACEIGAAAwgpABAACMIGQAAAAjCBkAAMAIQgYAADCCkAEAAIwI6O4CAAAnjxUrvNvfzTd7tz90LWYyAACAEYQMAABgBCEDAAAYwTUZOGl4+1wwAOCXYSYDAAAYQcgAAABGEDIAAIARhAwAAGAEIQMAABhByAAAAEZ0OmQUFRUpMTFREydOVE5OTrv2qqoqzZw5UwkJCZo1a5bq6uratG/ZskX/9m//9ssrBgAAPUKnQkZNTY2ys7O1evVqFRQUaPs4kenRAAAKOUlEQVT27SouLm6zzuLFizVt2jQVFRUpOTlZmZmZkqTW1lb9z//8j+6++265XC7v7wEAAPBJnQoZJSUlio6Ols1mk9VqVUpKigoLCz3tDodDpaWlSkpKkiSlpKRo8+bNcjgc+uqrr/Ttt9/q4YcfNrMHAADAJ3Xqjp/V1dUKDQ31LIeGhqqqqsqz3NDQoMDAQFmt1kOdBgQoKChIdrtdI0eO1COPPKLS0tIOt1NeXq7m5maVlZWd6H7AsP379/v8uNjtgd1dQrc4ePCg7Hb7MdsbGy2SQiRJDQ31crvdXVRZ79bRuKBzysqavd5nT3g/6yk6OkPRqZBxtE4sFkuHG/HzO7HrSiMiIlRRUaHIyMgTeh3MKysr8/lxsdm6u4LuYbfbZTvOzh9xqCo4OEQhIV1QFDocF3ROZKT3/4Y94f2sp2htbVVtbe0x2zuVAsLDw1VTU+NZrq6uVnh4uGfZZrOpqalJTqdTkuR0OtXc3Kzg4OCfWzcAAOjhOhUyYmJitHXrVtXW1srhcCgvL0+xsbGedqvVqqioKOXn50uS8vPzFRUV5Tl9AgAAep9OhYywsDClpaUpNTVVSUlJGjFihCZMmKD09HS9++67kqSMjAytX79eiYmJWrNmjRYuXGi0cAAA4Ns6/VPvCQkJSkhIaPPc4a+pStLgwYO1atWqY75+zJgxGjNmzIlXCAAAeiTu+AkAAIwgZAAAACMIGQAAwAhCBgAAMIKQAQAAjCBkAAAAIwgZAADACEIGAAAwgpABAACM6PQdPwFvW7GiuysAAJjETAYAADCCkAEAAIwgZAAAACMIGQAAwAhCBgAAMIKQAQAAjOArrAAAn2Xiq+4JCd7vE0fHTAYAADCCkAEAAIwgZAAAACMIGQAAwAhCBgAAMIKQAQAAjCBkAAAAI7hPBjrllVcCZbN1dxUAgJ6EkAEA6FW8/Y+mm2/2Xl8nG06XAAAAIwgZAADACEIGAAAwgpABAACMIGQAAAAj+HbJScrEzyMDAHAimMkAAABGEDIAAIARnC4BAOAXMHF6+mS5wVenQkZRUZFycnLkcDiUnJysuXPntmmvqqrSvHnzVFtbq0GDBmnZsmUaOHCgHA6HMjIy9Mknn8hisSgrK0vnnXeekR0BAOBk4e3g0l2hpcPTJTU1NcrOztbq1atVUFCg7du3q7i4uM06ixcv1rRp01RUVKTk5GRlZmZKkl555RW53W4VFhZq+fLluu++++R0Os3sCQAA8CkdzmSUlJQoOjpatv9/o/eUlBQVFhbq8ssvlyQ5HA6VlpZq+fLlnvasrCw5HA5t2rRJt99+uyTp7LPPVlhYmD7++GNdcsklbbbhdrslSa2trXK5XGptbfXeHvYAL77o/T79vHy1jb+/S35+vWtceoqOxsbfXwoM/Ndjb/+/gaPjmPFdvXFsTH2sHv68Pvw5/lMdhozq6mqFhoZ6lkNDQ1VVVeVZbmhoUGBgoKxW66EOAwIUFBQku92uqqoqhYWFtXntvn372m3D5XJJknbs2CFJqq2t7XDHTiYXXNDdFXTWd91dAI7p+GMTF9dFZeAnOGZ8V+8am08+Mdv/4c/xn+owZBzthRaLpcOO/fz8jpps/I7yz6iAgACNGjVKfn5+bfoGAAC+y+12y+VyKSDg6HGiw5ARHh6ubdu2eZarq6sVHh7uWbbZbGpqapLT6VRAQICcTqeam5sVHByssLAw1dTUaOjQoZIOXd9x5MzGYX5+furTp88J7xwAAPBdHZ6djYmJ0datW1VbWyuHw6G8vDzFxsZ62q1Wq6KiopSfny9Jys/PV1RUlKxWq2JjY7V27VpJ0u7du7V3716NHj3azJ4AAACfYnEf62qNIxQVFem5555TS0uL4uPjdf/99ys9PV3x8fEaP368KisrNX/+fNXU1GjAgAFaunSphgwZopaWFi1evNjzFdb09HTFxMR0xX4BAIBu1qmQYcrSpUv13nvvyWKxaPr06UpNTdW2bduUmZmpgwcPKjo6WgsXLjzmuR6Y9R//8R+qr6/XY489pi+//FILFizQ999/r9/85jd67LHH1K9fv+4usde555579Nlnn6lv376SpLlz52ro0KGMTTd77733lJOTox9++EGXXXaZFixYwHuZD3jttdf0+uuve5YrKioUGxurOXPmcMx0FXc32bRpk/uPf/yj2+l0ug8cOOCOi4tz79692x0bG+v+5ptv3C6Xy33fffe5X3vtte4qsVfbsmWLe8yYMe7777/f7Xa73cnJye7S0lK32+12P/XUU+6lS5d2Z3m91oQJE9z19fVtnmNsutfevXvdl112mbuystLd0tLi/sMf/uB+5513eC/zMd988407Li7OvW/fPo6ZLtRt35iPjY3Viy++KH9/f9XV1am1tVV2u10RERH69a9/7ZndKCws7K4Se62GhgYtW7ZMt956qySpsrJS+/fvV1RUlCRpxowZjEs3qK+vl91u13333afJkycrJyeHsfEB77zzjiZNmqTw8HBZrVYtW7ZMv/rVr3gv8zFLlizRnXfeKZfLxTHThbp17u7wAblq1SolJCRo37597e6rceQ9OdA1HnroId19992qrKyUpKPe74Rx6Xp1dXUaO3asFi9erFNOOUW33nqrAgICGJtu9t1336lPnz6aPXu2qqurFR8fr+HDhzMuPmT79u2qra3VlClT9Pe//52x6ULdfu+/u+++W//3f/+niooK7dmzp107983oWmvWrNHgwYPbXKDb0b1S0DXOPvtsPf300woJCVG/fv00c+ZMbdmypd16jE3Xam1tVXFxsbKzs7VmzRrt2LFD5eXl7dZjXLrPa6+9ptTUVFksFt7Puli3zWR89dVXcrlcGjFihPr166crr7xSGzZsaDPYNTU1be7JAfMKCwtVU1OjKVOmqLGxUT/88IMsFotqamo86zAu3WPHjh2qqalRfHy8pH+FP8ame5122mmKiYnRwIEDJUnjx4/nvcyHOBwObdmyRYsXL5Z06N5PHDNdp9tmMnbv3q1FixbJ4XCopaVF77zzjq677jrt2bNHu3fvliStXbu2zT05YN6LL76ot956S+vXr9edd96p+Ph4Pfroozr11FM9N2Vbu3atxo0b182V9j4Oh0NZWVlqampSS0uLXn/9dV1zzTWMTTeLi4vTli1b1NjYqNbWVn3wwQeaNGkS72U+4ssvv9TQoUMVFBQkSRoyZAjHTBfqtpmMq666Srt27dKUKVPk7++vq666SpMmTdJpp52me++9VwcPHtTo0aN1/fXXd1eJOMITTzyhhQsXav/+/YqIiNATTzzR3SX1OhdeeKGuv/56zZgxQ62trbryyiuVlJSk3/zmN4xNNzrvvPM0Z84c/fGPf5TD4VBMTIxmzJihs846i/cyH7B3714NHjy4zXO8n3Wdbr1PBgAAOHl1+4WfAADg5ETIAAAARhAyAACAEYQMAABgBCEDAAAYQcgAAABGEDIAAIARhAwAAGDE/wNBCvhj0FTe7QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(8, 5), dpi=80)\n",
    "ax.hist(boys_ppc['first'],\n",
    "        bins=22, density=True, alpha=.4, color='b', label='sampled firstborn')\n",
    "ax.axvline(x=first_boys, color='b', label='observed firstborn')\n",
    "ax.legend()\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**H5** The model assumes that sex of first and second births are independent. To check this assumption, focus now on second births that followed female first borns.\n",
    "\n",
    "Compare 10,000 simulated counts of boys to only those second births that followed girls.\n",
    "To do this correctly, you need to count the number of first borns who were girls and simulate that many births, 10,000 times. Compare the counts of boys in your simulations to the actual observed count of boys following girls. How does the model look in this light? Any guesses what is going on in these data?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf4AAAEfCAYAAACgd6QCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAANEgAADRIBtYA3dAAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XtclHXe//EXRyVRXFFEwtjWUtk8bZmVYRm4eUpFhGLXLNtN+WmUZneZuYRK6XqnecokrXTDQ4mwiqV5gLb1tJapiYbbXasjiHJMcyQ5DPP7g5wNAR2cUdTr/Xw8enRdw/e6vp/ryzhvrsNcl4vVarUiIiIihuDa0AWIiIjI1aPgFxERMRAFv4iIiIEo+EVERAxEwS8iImIgCn4REREDcW/oAn6psrKSc+fO4e7ujouLS0OXIyIics2yWq1UVFTQuHFjXF3t34+/poL/3LlzZGVlNXQZIiIi143g4GBuuukmu9tfU8Hv7l5VTnBwMB4eHg1czfUlJyeHwMBAh9bxzTcQGlo1nZEBv/2tEwozAGeMvVwejX3D0vjb5+zZb/j666oP165dM2jSxPEP15ycHFq3bk1WVpYtO+11TQX/+cP7Hh4eCv56cnV1dXjMXFyguPi/0/oV2McZYy+XR2PfsDT+9nF3d8FqLbZNO2PMfjn29T01rov7REREDETBLyIiYiAKfhEREQNR8IuIiBiIgl9ERMRAFPwiIiIGouAXERExEAW/iIiIgdgd/Onp6QwaNIi+ffsyceJESktLa21XVlbGyJEj2bp1q+01i8XCm2++SXh4OH379uXdd991vHIREYMbMWJEtc/aa0ldtYWGhl4Tt2a3WCzExMTQt29fUlNTmTx5MocPH74ifS1a9CG7dzet8XpOTg7du3e/In1ejF137issLCQuLo7Vq1cTGBjIlClTSExMZNy4cdXaZWZmEh8fz/fff8/jjz9ue33p0qUcPHiQ1atXU1ZWxpAhQ+jRowddunRx7taIyBWx+KvFDi0/+q7RTqrk6rFYzvLTT99f0T68vG7Dzc3+e6yL8+Tl5bF792727t2Lq6srb731VrXccqYxY6LZs2fGFVn35bAr+Ldv3063bt1s92SOjo7mmWeeqRH8K1euZPz48SxZsqTa62lpaUyePBlPT088PT1ZtmwZzZs3d9ImiIg4l8Vyln/96zeUl+df0X48PPy4997/4ObW5KLtli9fzqpVq3B1dSUoKIj4+HhatWoFwMaNG3n77bcpLS0lKiqKkSNHcvbsWV566SWOHz+Oi4sLvXv3Zty4cVitVhYuXEh6ejqVlZXcfvvtxMfH07RpU0JDQ+nSpQuHDx9m5MiR7Nq1i3nz5gGwa9cuZs2aRUpKClu2bOGdd97BYrHg4eHB9OnTue2228jLy+Pll1+moKCAwMBATp06Vef2fPDBBxw+fJiKigrGjBnDgAEDAJgzZw5btmzBzc2NTp062fbCJ02axObNm3FxcSEvL49hw4aRkZHBwoULycjIwMPDg8DAQGbOnImXl1e1vjIyMli8eDFlZWUUFxczfPhwnnzySZ5++mlKS0sZOnQoXbp0IT8/n/HjxzN//nyaNWvG1KlTyc/Px2KxMHz4cB599FF2797NjBkz8PT0pLKyktWrV1d7Kl5ycjLvvvsuN910E3fddRcZGRlkZGQQH/8WXl7N6dPnFD17DufBB3vz7bffkpiYeFnvG0fZFfz5+fn4+/vb5v38/MjLy6vRbsaMqr9oLgx+k8lEVlYWCxcu5NSpU0RERDBy5Mg6+8vJyanXIwYFzGYzJpPJoXXk5noAAT9P5+LjU+6Eym58zhj7a11RUZFDy1+p8blSY19Z+RMWS6XT13shi6WSY8eycXX1qrPN3r17Wbp0KfPmzaNp06YsX76cCRMm8Nprr3Hu3DmKi4uZNWsWJSUlPPPMM7Ru3Zrc3FwqKiqYN28e586d48033yQrK4vt27eTlZXF7NmzcXNzY/ny5cTHxzNu3DgqKiq4+eabef755zl79iyzZs3i0KFDeHt7s3z5ch566CF27NjBnDlzmDVrFk2aNGH37t3ExMTw/vvvM2XKFG6//XamTJnCsWPHGDNmDPn5+TV+PxUVFVgsFubOnUt+fj7PPPMMLVu2ZO/evezYsYO5c+fi6enJ3LlzmTJlCs8++yyurq6sX7+erl27snLlSnr37s3XX39NcnIyK1aswNXVlffee4/PP/+c4OBgW19Wq5XExEQmTJhAmzZtyM3N5c9//jO9e/cmPj6eMWPGMH/+fAA+//xzJk6cSOPGjXnuuef44x//yJ133klJSQkTJkygadOmlJeX83//938kJSXh6+tLdna2rS+TycSbb77J22+/za9+9Svmzp1LRUUFJpOJn376ifN/j5SWltGjRw/+53/+h5MnT1JZWXlZ72Gz2UxOTk69lwM7g7+ysuY/gPoEc0VFBf/+97957733OHPmDCNGjKBt27aEhYXV2j4wMFAPfqgnk8lEUFCQQ+s4ffq/0wEBATi4OsNwxthf63wLfR1a/kqNz5Uc+1tuOXpNHOr/6KOPCA8Pp1OnTgA8++yz9OzZk5tvvpnGjRvz1FNP8Zvf/AaAvn378v333zN48GD+9re/MW3aNHr27ElcXByBgYEsWrSIrKwsnn/+eaDqs7lFixYEBQXh7u5OWFiYbTzDwsI4ePAggwYNYu/evcyYMYP169dTVFTEpEmTgKprukpLS/Hx8WH//v28/vrrP392BNGjRw/8/Pxq/H7c3d2JiYkhKCiIoKAg7r77bnJzc8nKyuKPf/wj7du3ByAmJoZnn32WoKAgRowYwT//+U8GDRpERkYG7777LoGBgdxyyy08//zzPPjggwwbNoyuXbvWGL+lS5fy2Wef8eWXX/Ldd99hsVjw9/fHy8vLdgTlfF1t2rShVatWHDx4kGXLlrFs2TLbdv7www+0a9eOwMBA7rzzzhr9fPbZZ/Tu3Zvf/e53AIwePZrY2FiCgoJqHIV4+OGHad26NW5ubtVqqA+TyURAQADF55+sVg92Bb+/vz+HDh2yzV94BOBS/Pz8GDhwIB4eHrRo0YIHHniAAwcO1Bn8IiINzc2tCd7eDX8d0oU7XpWVlVgsFtu8m5ubbdpqteLh4cEtt9zC5s2b2bVrFzt37uTRRx9l0aJFWK1WxowZQ1RUFAAlJSWcO3fOtvwvA2rYsGHMnTuXxo0b07NnT5o1a0ZlZSW9evVi1qxZABw9epRGjRrh4+ODi4tLtVp/WdeFLqzZ3d291u2sqKgAYMiQISxYsIB//OMftj8sAFatWsW+ffvYuXMnL7zwAk888QRPPPGEbR0lJSVERETw8MMPc9dddxEREcGnn36K1Wqts7bzdSQnJ+Pp6QlUHfHy9vZm//79NG7cuNblXF1dq633YjvHF/4hcLXZtdseEhLC3r17bYc1kpOT6xXaffr0IS0tDavVSklJCbt27aJz586XV7GIiIHcf//9fPLJJ5z++ZDcihUr6NGjh+0Z7KmpqVitVoqLi9m6dSshISGsWLGCqVOn8tBDD/HKK69w2223ceTIEXr27ElycjIlJSVYrVYSEhKYPXt2rf3efffdFBUV8cEHHxAZGQnAPffcw7Zt2zh27BgAW7ZssZ227dWrF2vWrAEgOzubPXv21LlNKSkpQNVe6/79+7n77ru5//77SU1N5dy5c1itVlatWkXPnj0BaNq0KX369CEhIYHHHnsMgKysLMLDwwkODiY2Npbw8PAaV+UfPXqUs2fP8txzz9G7d2+2bdsGUO0Pp/Pc3NyoqKjA29ubzp07k5SUBEBBQQFDhgwhMzPzYr8m7r//fnbs2EFhYaHt93KtsmuPv2XLliQkJDB27FjKy8tp3749M2bMIDMzk/nz59c4p3+hCRMm8Ne//pWBAwdisVh45JFH6NOnj1M2QETkRtarVy+io6MZPnw4FouFtm3bMnPmTNvPmzRpQkREBKWlpYwdO5Y77riDoKAgdu3axYABA2jUqBEdO3ZkwIABuLu7c/LkSaKiomwX973++uu19uvi4kJERAQfffQR9957LwDt27fn1VdfJTY21ranPn/+fFxcXIiPj2fSpEkMHDiQNm3a0KFDhzq36aeffiI8PByLxcLrr79O69atefTRRzl+/DjDhg2joqKCTp06ER8fb1tm6NChZGRk2LIjODiY3r17Ex4eTpMmTWjWrBmvvfZatX46duxISEgI/fv3p2nTprRr145bbrmFY8eO4etb/fRVnz59mDBhAjNmzGD27NlMnTqVtWvXUl5ezpgxY+jevTu7d++uc5vatWvHCy+8wJNPPkmjRo1o165dg+/Z18XFerFjHldZeXk5Bw4coEuXLjrHX0/OONd54ACcP0X29degb1vaxwjn+K/Vr/MZYeyvZVdr/K1WK3PnzsVqtTJhwoQr3t/lOHbsGBs2bCAmJgYXFxeWLVvG/v37mTt3LmbzAfbsqfpw7d79a6ecQjp/jv9yMtOuPX4REZGG0rt3bwICAhrs62/2aNOmDdnZ2QwaNAgXFxf8/PzqPJrS0BT8IiJyTfv8888buoRL8vDwuGaD/kL6sryIiIiBKPhFREQMRMEvIiJiIAp+ERERA1Hwi4iIGIiCX0RExEAU/CIi16kRI0awdevWhi6jVnXVFhoaSlZWVgNUVJ3FYiEmJoa+ffuSmppqewTw5crJyaF79+5OrPDK0ff4RURqcfYsfH9lH87HbbfBTRd/OJ9cIXl5eezevZu9e/fi6urKW2+9xeOPP97QZV0VCn4RkQucPQu/+Q3k51/Zfvz84D//gSZNLt5u+fLlrFq1yvYI1/j4eFq1agXAxo0befvttyktLSUqKoqRI0dy9uxZXnrpJY4fP46Liwu9e/dm3LhxWK1WFi5cSHp6uu1e/fHx8TRt2pTQ0FC6dOnC4cOHGTlyJLt27WLevHkA7Nq1i1mzZpGSksKWLVt45513sFgseHh4MH36dG677Tby8vJ4+eWXKSgoIDAwkFOnTtW5PR988AGHDx+moqKCMWPGMGDAAADmzJnDli1bcHNzo1OnTra98EmTJrF582ZcXFzIy8tj2LBhZGRksHDhQjIyMvDw8CAwMJCZM2fWuD9+RkYGixcvpqysjOLiYoYPH86TTz7J008/TWlpKUOHDqVLly7k5+czfvx45s+fT7NmzZg6dSr5+flYLBaGDx/Oo48+yu7du5kxYwaenp5UVlayevXqak/hs1gsTJo0iW+++QYvLy+mTJlCx44dOX36DIsXt+HECU+8vCYwdOhjjBo1irfffpucnBymT58OwIYNG/j44495++236/1eqg8Fv8h14Fq9V/6NysXl2ulr586dLF++nI8++ggfHx8WLlzI5MmTWby46j1x7tw5kpOTMZvNRERE0LVrV0wmE40aNWLt2rX89NNPTJ48GbPZzObNm8nOzmbNmjW4ubnx1ltv8cYbbzBt2jQAOnXq9PO95c3MmTOHH3/8kWbNmrF27VqGDRvGkSNHWLhwIUlJSTRt2pT169cTGxvLp59+yrRp07j77rsZO3Ys33//PeHh4XVuU5MmTfj73/9Obm6ureadO3fy1VdfkZqaSqNGjYiLi2P27NnEx8fTuHFjvvjiC+655x5SU1MZPHgwRUVFpKam8vnnn+Pq6srs2bP59ttv6Xr+gSNU3eN/6dKlvPHGG7Rt25Zjx47Rv39/nnzySRYvXkx4eDjr1q0DYMeOHcydO5cOHTowYsQIxowZQ8+ePTGbzQwfPpzbb78dgO+++4709HRat25dY7tKSkro3bs3M2bMYOvWrYwfP56NGzfyxhtL8fUtZ/ToE3Ts+CH/7//9lcDAQCIjI3nkkUeYPHkyTZo0ITk52fa0wytJwS8icoGbbqraE78WDvVv376dgQMH4uPjA1SdO1+0aJHtWfWPPfYYbm5u+Pj4EBYWxo4dOxg8eDBvvvkmf/7zn+nZsycTJkzA29ubf/zjHxw8eJCIiAgAKioqaNGiha2vO++8EwBvb2969+7Nxo0bGTRoEP/85z+ZPHky69ev58SJE7ZD4mVlZZjNZk6dOsXOnTuZPHkyUPWkurvvvrvObfrDH/4AQEBAAHfddRdffvkl27dvJyIiwva8++HDh/Pss8/atjE5OZkePXqQmprKu+++i5+fHzfffDMRERE8+OCD9OnTp1roQ9UTBhMTE/nss89Yu3Yt3333HRUVFZSWltZZW0lJCXv27Kn2BESz2UxWVhbt2rXj5ptvrjX0oepJtn379gWqnvY3adIkTpw4wb/+9TX/8z+nfh7bmwgPD2fbtm3079+fe++9lw0bNnDfffdhMpno1atXnbU5i4JfRKQWTZpcG0+orKysrDH/y+fJu7m52aatViseHh7ccsstbN68mV27drFz504effRRFi1ahNVqZcyYMURFRQFVIXfu3Dnb8r88TD5s2DDmzp1L48aN6dmzJ82aNaOyspJevXoxa9YsoOp5940aNcLHxwcXF5dqtf6yrgtdWLO7u3ut23n+j5shQ4awYMEC/vGPfxAQEGB7IuCqVavYt28fO3fu5IUXXuCJJ57giSeesK2jpKSEiIgIHn74Ye666y4iIiL49NNPudhDac/XkZycjKenJwBFRUV4e3uzf/9+2x8mtfnlYf+LbZvVarVtW3R0NAsWLLCdwrhwHVeCruoXEbmG3X///XzyySecPn0agBUrVtCjRw/c3av221JTU7FarRQXF7N161ZCQkJYsWIFU6dO5aGHHuKVV17htttu48iRI/Ts2ZPk5GRKSkqwWq0kJCQwe/bsWvu9++67KSoq4oMPPiAyMhKAe+65h23btnHs2DEAtmzZYjs03atXL9asWQNAdnY2e/bsqXObUlJSgKpHy+7fv5+7776b+++/n9TUVM6dO4fVamXVqlX07NkTgKZNm9KnTx8SEhJ47LHHAMjKyiI8PJzg4GBiY2MJDw+vcVX+0aNHOXv2LM899xy9e/dm27ZtANX+cDrPzc2NiooKvL296dy5M0lJSQAUFBQwZMgQMjMzL/ZrAiA/P9/Wx8cff8zNN9+Mn58f997blc8/bw7A2bM/sW7dOtu23XfffRQXF/P3v//dNs5Xmvb4RUSuYb169SI6Oprhw4djsVho27ZttcPQTZo0ISIigtLSUsaOHcsdd9xBUFAQu3btYsCAATRq1IiOHTsyYMAA3N3dOXnyJFFRUbaL++p6opyLiwsRERF89NFH3HvvvQC0b9+eV199ldjYWNve7Pz583FxcSE+Pp5JkyYxcOBA2rRpQ4cOHercpp9++onw8HAsFguvv/46rVu35tFHH+X48eMMGzaMiooKOnXqRHx8vG2ZoUOHkpGRQZ8+fQAIDg6md+/ehIeH06RJE5o1a8Zrr71WrZ+OHTsSEhJC//79adq0Ke3ateOWW27h2LFj+Pr6Vmvbp08fJkyYwIwZM5g9ezZTp05l7dq1lJeXM2bMGLp3787u3bsv+rvy9/cnLS2NWbNm0bRpU+bMmQPAiy8+xcSJnzJ1ahAeHhMZNGiY7RoIFxcXBg8ezKFDh+o8heBsLtaLHfO4ysrLyzlw4ABdunTBw8Ojocu5rphMJtvhr8t14ACcP0X29dfXxmHO64Ezxv5SGvrivobuvy5XY+ylbldr/K1WK3PnzsVqtTJhwoQr3p+zmc0H2LOn6sO1e/ev8fb+74drWVkZTz/9NKNHjyYkJMTudZpMJgICAi4rM3WoX0RErmm9e/fmiy++4M9//nNDl+JUhw8fJiQkhKCgoHqFvqN0qF9ERK5pn3/+eUOXcEV07NiRL7744qr3qz1+ERERA1Hwi4iIGIjdwZ+ens6gQYPo27cvEydOrPMGCGVlZYwcObLWhzOcOXOGPn36XLMPlRAREbnR2RX8hYWFxMXFsWjRIjZt2oSXlxeJiYk12mVmZhIdHc2+fftqXc/kyZM5c+aMYxWLiIjIZbPr4r7t27fTrVs3AgMDgao7DT3zzDOMGzeuWruVK1cyfvx4lixZUmMdSUlJ3HzzzbRv394JZYuIXHmOfo3xUvQMBWkIdu3x5+fn4+/vb5v38/MjLy+vRrsZM2bwwAMP1Hg9MzOTjRs3XpffvxQRaUi7d+9myJAhl73MiRMnGDhwYK13tqvL5SzTEDp06MCPP/5IamoqY8eOvSp9jhgx4ro/XW3XHv+F9xmGmvckrsuZM2eIi4tjwYIFdt9gICcn56rcr/hGYjabMZlMDq0jN9cDCPh5Ohcfn3InVHbjc8bYX0pRUZFDyztaX0P3X5crPfaObvel2FN7Xl4eZWVl9drOXy6zZcsWmjdvbnv0qz3rsXeZq/Hev5Ts7GyKioooKSm5KrWcO3eO/Pz8evVVVpZrmz5xIhdPTx+H6zCbzeTk5FzWsnYFv7+/P4cOHbLNX3gE4GK2bdvGjz/+SGxsLADHjh1j+vTpmM3mOh/bGBgYqDv31ZMz7qD1863AAX5+EIaDRRnE1bh7mW+h76UbXYSj9TV0/3W50mPv6HZfij21nzx5krKyMmbPns3Ro0eprKzkf//3f/ntb3/Lvn37mDVrFqWlpRQWFvL73/+eyZMnc/LkSTw9PSkqKmL58uWcOXOG1157rcZp2A0bNpCYmIjVaqVly5a8+uqr/PDDD3UuU1FRQXx8PJmZmbi5uREUFMScOXNwcXHhww8/5KOPPsJqteLv78+UKVPw9/cnPz+fuLg4jh07hpubG6NGjWLIkCEcPnyYhIQEfvzxRzw8PBg3bhwPPvggqamppKenY7FYMJlMtGrVirlz59KiRQv27NlDQkICLi4u/O53vwOgbdu2+Pr6cubMGaZMmUJhYSG//vWvSUhIoHnz5hftZ82aNZSUlODv788dd9xBbm4ueXl55OTk0KFDB954440aD+Vp3LgxX331FWvWrKG0tJSoqCjb8wqWL1/OqlWrcHV1JSgoiPj4eMrLyxk27BVef92FRo2stGjRksGD/0haWhobN25k9erVuLu707x5c2bOnEmrVq3seu+cv3NfcXGxXe1/ya7d6pCQEPbu3Ut2djZQ9dSisLAwuzoYMGAAGRkZrFu3jnXr1tGpUydeeeWViz6rWURE/isvL4+YmBjS0tLo27cvixdXXXuQlJTEhAkTWLNmDWlpaaxZs8b2AB2oeszuc889xz333FMj9L/77jumT5/OO++8w/r16xkyZAjPPvssv/vd7+pcZv/+/Xz33Xe2vqDqCO2//vUvPv30Uz766CPWrl1Lv379bI/onTZtGsHBwWzcuJH333+fhQsXYjabGTt2LKNHj2b9+vW88cYbTJw4kePHjwPYAn7jxo00adKE1atXU1ZWxrhx44iLi2Pt2rV06tSpWm0mk4lp06axfv162rRpw5tvvklFRcVF+zl69CgrV660Xay+d+9eFixYwMaNGzl+/Dhbtmyp9fdx7tw5kpOTWblyJUlJSbYnBC5fvpyVK1eyfv16goODmTx5MgEBAXTt2oGvvmoKQHr6brp3707z5s2ZPXs2K1euJDU1lQceeMCuBwE5g13B37JlSxISEhg7diz9+vWjoKCA2NhYMjMzGTVq1JWuUUTE0IKCgrjjjjuAqru9nT8F8de//pXi4mIWLVrEtGnTqKiooKSkxK517t69m5CQENq0aQNAeHg4J0+e5OTJk3Uuc/vtt3P69Gn+8Ic/kJiYyNChQ2nbti2ff/453333HVFRUQwZMoT33nvPdih8586dtscA+/n5sXnzZlsfDz74IADt2rXjzjvv5MsvvwSgU6dOtj3f4OBgiouL+fbbb2ncuDHdu3cHICIiwvbYXKh6mFHbtm1tP9uxYwdHjx69aD8dO3bkpptusq3jnnvuoUmTJri5uXH77bfXuTf92GOP4ebmho+PD2FhYezYsYPt27czcOBAfHyqDuOPGDGCnTt3UlFRwbBhD7N9e9Xra9emEx0djZubGw8++CDDhg1jxowZdOjQgdDQ0Ev81pzD7lv2hoaG1iiqc+fOdV7BX5eL/UxERGo6/wheqHqa23nDhw8nODiYBx54gEceeYTdu3df9Fnzv1TbtVuVlZW258TXxsfHh48//pgvv/ySnTt38vLLLzNt2jSsVivDhg3j+eefB6ru53Lq1KkatQMcOXKk1sfinu/b1dWVRo0aVdteq9Vq+/8vX3dzc7PN//K6sPNPDrxUPxcexq+t39r8sl+r1YqHh0eN8aysrLT1HxLyO157zZ1vvrmJgoJi29MO586dS1ZWFtu3b2fGjBncd999TJo0qdY+nUlX0ImIXIdOnz7NoUOHePHFF+nTpw9Hjx4lPz+/1rCrzb333su2bds4ceIEAGvXrsXX19f2te3aZGRkMHr0aHr06MELL7xA9+7d+fbbb+nZsycff/yxbQ/5nXfe4aWXXrL1s3btWqDqYskRI0bQtGlTrFar7R7833//PXv27OGee+6ps+/27dtTXl7Orl27gKrrE3766Sfbz3fu3Gk7krBq1SoeeOABbr311nr3Y4/U1FSsVivFxcVs3bqVkJAQ7r//fj755BNO/3yx1IoVK+jRowfu7u64urpy//2nWb68NUOGhOLi4kJxcTGhoaG0bt2aUaNG8dRTT121b1DoIT0iInW4lr9n7+Pjw9NPP82QIUNo3rw5fn5+dOnSpdZnzdfm9ttvZ9KkScTExGCxWGjevDlvv/12tSMKF3rggQfIyMhg4MCB3HTTTTRr1ozHH38cHx8fnnjiCZ544gmg6pD+X//6VwDi4uKIj49n8ODBALbz3m+99Ravv/46s2bNAmD69Om0bdvWdhj+Qh4eHixcuJApU6ZgsVjo2rUrTZo0sf28Y8eOvPjii/zwww907NiRl156CU9Pz3r3Y48mTZoQERFBaWkpY8eOtZ2GiY6OZvjw4VgsFtq2bcvMmTNty9x334988okvgwc/BECLFi3405/+xPDhw/Hy8qJRo0bExcVddk314WK197jQVVBeXn5ZzxYW51zdfOAAdK16ZDRffw1duly8vVS5Glf1O3ojGUcDrKH7r8vVeh681E7jb58zZ75m7twH+fe/vUhM3IS3t+Mfruev6r+czNQev4iIyBX05z/HUVhbZoDhAAAbl0lEQVTYnLFjjzd0KYCCX0RE5Ip6//3X2LOna0OXYaPgF5Erzhn3vL+Wz7eLXE90Vb+IiIiBKPhFREQMRMEvIiJiIAp+ERERA1Hwi4iIGIiCX0RExEAU/CIiIgai4BcRETEQBb+IiIiBKPhFREQMRMEvIiJiIAp+ERERA1Hwi4iIGIiCX0RExEAU/CIiIgai4BcRETEQu4M/PT2dQYMG0bdvXyZOnEhpaWmt7crKyhg5ciRbt261vZaXl0dMTAyDBw9m4MCBLFu2zOHCRUREpP7sCv7CwkLi4uJYtGgRmzZtwsvLi8TExBrtMjMziY6OZt++fdVenzp1Kj169CAtLY0PP/yQVatW8eWXXzpnC0RERMRudgX/9u3b6datG4GBgQBER0eTlpZWo93KlSsZP348Xbp0qfb6gAEDiIyMBKBp06YEBQVx/PhxR2sXERGRenK3p1F+fj7+/v62eT8/P/Ly8mq0mzFjBgBLliyp9vojjzxim96xYwf79+9n+vTpdfaXk5ODq6suP6gPs9mMyWRyaB25uR5AwM/Tufj4lDuhshufM8b+UoqKihxa3tH6HO3fGWrbhqsx9lI3jb99yspybdMnTuTi6enj8DrNZjM5OTmXtaxdwV9ZWVnjtcsJ5g0bNjBt2jTmzZtHy5Yt62wXGBiIh4dHvddvZCaTiaCgIIfWcfr0f6cDAgJwcHWG4YyxvxTfQl+Hlne0Pkf7d4batuFqjL3UTeNvH7P5NOcPcrdpE4C3t+NjZjKZCAgIoLi4uN7L2hX8/v7+HDp0yDZ/4REAe8yfP5+UlBTef/99fvvb39avShEREXEKu4I/JCSEN954g+zsbNq2bUtycjJhYWF2d/LWW2+xefNmkpOT8fPzu+xiRa5Xi79a3NAliIgAdl7c17JlSxISEhg7diz9+vWjoKCA2NhYMjMzGTVq1EWXPXv2LImJiZjNZkaNGsWQIUMYMmQIa9eudcoGiIiIiP3s2uMHCA0NJTQ0tNprnTt3rnEhH0BSUpJtukmTJhw8eNCBEkVERMRZdOm8iIiIgSj4RUREDETBLyIiYiAKfhEREQNR8IuIiBiI3Vf1i4g0pNruhVBUVGT3XQVH3zXa2SWJXJe0xy8iImIgCn4REREDUfCLiIgYiIJfRETEQBT8IiIiBqLgFxERMRAFv4iIiIEo+EVERAxEwS8iImIgCn4REREDUfCLiIgYiIJfRETEQBT8IiIiBqLgFxERMRAFv4iIiIHYHfzp6ekMGjSIvn37MnHiREpLS2ttV1ZWxsiRI9m6davttZKSEsaPH0///v3p378/u3btcrxyERERqTe7gr+wsJC4uDgWLVrEpk2b8PLyIjExsUa7zMxMoqOj2bdvX7XX582bh6+vLxs3biQxMZGJEydy5swZ52yBiIiI2M2u4N++fTvdunUjMDAQgOjoaNLS0mq0W7lyJePHj6dLly7VXk9PTycyMhKAoKAgOnfuTHp6uqO1i4iISD2529MoPz8ff39/27yfnx95eXk12s2YMQOAJUuWVHs9Ly+P1q1b2+Zbt27NyZMn6+wvJycHV1ddflAfZrMZk8nk0Dpycz2AgJ+nc/HxKXdCZTc+e8a+qKjoKlVTO0ffGw1df11KS0vtrs3RMZCanPG5YwRlZbm26RMncvH09HF4nWazmZycnMta1q7gr6ysrPFafYLZarXWa/nAwEA8PDzsXr9UfagFBQU5tI7Tp/87HRAQgIOrMwx7xt630PcqVVM7R98bDV1/XYqKivD1ta82R8dAanLG544RmM2nOX68arpNmwC8vR0fM5PJREBAAMXFxfVe1q709vf3p6CgwDZ/4RGAK728iIiIOIddwR8SEsLevXvJzs4GIDk5mbCwMLs7CQsLIzk5GYDs7Gz2799PSEjIZZQrIiIijrAr+Fu2bElCQgJjx46lX79+FBQUEBsbS2ZmJqNGjbrk8s8++yw//PADAwcOJCYmhilTptCiRQuHixcREZH6sescP0BoaCihoaHVXuvcuXONC/kAkpKSqs17e3sze/bsyyxRREREnEWXzouIiBiIgl9ERMRAFPwiIiIGouAXERExEAW/iIiIgSj4RUREDETBLyIiYiAKfhEREQNR8IuIiBiIgl9ERMRAFPwiIiIGouAXERExEAW/iIiIgSj4RUREDETBLyIiYiAKfhEREQNR8IuIiBiIgl9ERMRAFPwiIiIGouAXERExEPeGLkDkalj81WKHlh9912gnVSIi0rDs3uNPT09n0KBB9O3bl4kTJ1JaWlqjzerVq+nfvz8PP/wwM2fOpLKyEoDi4mJiYmIYNGgQgwYNYt26dc7bAhEREbGbXcFfWFhIXFwcixYtYtOmTXh5eZGYmFitzeHDh0lMTGTVqlVs2LCBI0eOkJqaCsDbb7/Nr3/9a9avX897773H1KlTKS4udv7WiIiIyEXZFfzbt2+nW7duBAYGAhAdHU1aWlq1Nunp6YSFhdG8eXPc3d2JjIy0tbFarZSUlNj+7+7ujouLi5M3RURERC7FrnP8+fn5+Pv72+b9/PzIy8ur1iYvL4+goCDbfOvWrW1tnnnmGf7whz8QEhLC6dOneeGFF/jVr35VZ385OTm4uuq6w/owm82YTCaH1pGb6wEE/Dydi49PuRMquzYUFRU5tPzFxtaesXe0f0c5+t5o6PrrUlpaandtjo6B1OSMzx0jKCvLtU2fOJGLp6ePw+s0m83k5ORc1rJ2Bf/5c/W/dGEwW63WOtskJCTQr18/xo8fz8mTJ3nyySfp2LEj9913X639BQYG4uHhYU9p8jOTyVTtD6/Lcfr0f6cDAgJwcHXXFN9CX4eWv9jY2jP2jvbvKEffGw1df12Kiorw9bWvNkfHQGpyxueOEZjNpzl+vGq6TZsAvL0dHzOTyURAQMBlnTa3K/j9/f05dOiQbf7CIwDn2xQUFNTa5rPPPmPjxo24uLjQpk0bwsLC+Ne//lVn8IuIczn6rQYRuXHYdTw9JCSEvXv3kp2dDUBycjJhYWHV2oSGhrJ161aKi4uxWCykpKTY2txxxx18+umnQNXhiV27dtGlSxdnboeIiIjYwa7gb9myJQkJCYwdO5Z+/fpRUFBAbGwsmZmZjBo1CoDg4GBiYmIYMWIEAwYMoFWrVkRHRwMwc+ZMdu7cyYABA4iOjuaRRx6p8YeDiIiIXHl238AnNDSU0NDQaq917tyZJUuW2OajoqKIioqqsWxgYGC1diIiItIwdOm8iIiIgSj4RUREDETBLyIiYiAKfhEREQNR8IuIiBiIgl9ERMRAFPwiIiIGouAXERExEAW/iIiIgSj4RUREDETBLyIiYiAKfhEREQNR8IuIiBiIgl9ERMRAFPwiIiIGouAXERExEAW/iIiIgSj4RUREDETBLyIiYiAKfhEREQNxb+gCRESuB4u/WuzQ8qPvGu2kSkQcY/cef3p6OoMGDaJv375MnDiR0tLSGm1Wr15N//79efjhh5k5cyaVlZUAWCwW3nzzTcLDw+nbty/vvvuu87ZARERE7GZX8BcWFhIXF8eiRYvYtGkTXl5eJCYmVmtz+PBhEhMTWbVqFRs2bODIkSOkpqYCsHTpUg4ePMjq1atJSUlh1apVHDhwwPlbIyIiIhdlV/Bv376dbt26ERgYCEB0dDRpaWnV2qSnpxMWFkbz5s1xd3cnMjLS1iYtLY2YmBg8PT3x9vZm2bJl3HrrrU7eFBEREbkUu4I/Pz8ff39/27yfnx95eXnV2uTl5VVr07p1a1sbk8lEVlYWTzzxBIMHDyY9PZ2mTZs6o34RERGpB7su7jt/rv6XXF2r/81gtVrrbFNRUcG///1v3nvvPc6cOcOIESNo27YtYWFhtfaXk5NTY/1ycWazGZPJ5NA6cnM9gICfp3Px8Sl3QmXXhqKiIoeWv9jY2jP2jvYvtSstLbV7bB3993El30PXK2d87hhBWVmubfrEiVw8PX0cXqfZbCYnJ+eylrUr+P39/Tl06JBt/sIjAOfbFBQU1NrGz8+PgQMH4uHhQYsWLXjggQc4cOBAncEfGBiIh4dHvTfGyEwmE0FBQQ6t4/Tp/04HBATg4OquKb6Fvg4tf7GxtWfsHe1faldUVISvr31j6+i/jyv5HrpeOeNzxwjM5tMcP1413aZNAN7ejo+ZyWQiICCA4uLiei9r1251SEgIe/fuJTs7G4Dk5OQaoR0aGsrWrVspLi7GYrGQkpJia9OnTx/S0tKwWq2UlJSwa9cuOnfuXO9iRURExDF2BX/Lli1JSEhg7Nix9OvXj4KCAmJjY8nMzGTUqFEABAcHExMTw4gRIxgwYACtWrUiOjoagAkTJuDl5cXAgQMZOnQoYWFh9OnT58ptlYiIiNTK7hv4hIaGEhoaWu21zp07s2TJEtt8VFQUUVFRNZb18vJi6tSpDpQpIiIizqAr6ERERAxEwS8iImIgule/XBccvU+6iIhU0R6/iIiIgSj4RUREDETBLyIiYiAKfhEREQNR8IuIiBiIgl9ERMRAFPwiIiIGouAXERExEAW/iIiIgSj4RUREDETBLyIiYiC6V7+IHS72rICioiJ8C32vYjUiIpdPe/wiIiIGouAXERExEAW/iIiIgSj4RUREDETBLyIiYiAKfhEREQNR8IuIiBiI3d/jT09PZ+7cuZSVldGtWzemTZtGo0aNqrVZvXo1S5cuxWKxEBYWxosvvoir63//tjhz5gxDhw7l5Zdfpk+fPs7bChGRS7jYvRhEjMSu4C8sLCQuLo7Vq1cTGBjIlClTSExMZNy4cbY2hw8fJjExkdTUVLy9vYmNjSU1NZXIyEhbm8mTJ3PmzBnnb4WIyDXO0T88Rt812kmViNHZdah/+/btdOvWjcDAQACio6NJS0ur1iY9PZ2wsDCaN2+Ou7s7kZGR1dokJSVx88030759eyeWLyIiIvVh1x5/fn4+/v7+tnk/Pz/y8vKqtcnLyyMoKMg237p1a1ubzMxMNm7cyN/+9jf+9Kc/XbK/nJycaqcI5NLMZjMmk8mhdeTmegABP0/n4uNT7oTKnKOoqKihS6hTaWnpNV3fjcxIY+/ov+8rwRmfO0ZQVpZrmz5xIhdPTx+H12k2m8nJybmsZe0K/srKyhqvXRjMVqu11jZnzpwhLi6OBQsW4OHhYVdRgYGBdreVKiaTqdofXpfj9On/TgcEBODg6pzqWr4XflFREb6+1259NzIjjb2j/76vBGd87hiB2Xya48erptu0CcDb2/ExM5lMBAQEUFxcXO9l7Qp+f39/Dh06ZJu/8AjA+TYFBQU12mzbto0ff/yR2NhYAI4dO8b06dMxm82Eh4fXu2ARERG5fHYdTw8JCWHv3r1kZ2cDkJycTFhYWLU2oaGhbN26leLiYiwWCykpKYSFhTFgwAAyMjJYt24d69ato1OnTrzyyisKfRERkQZgV/C3bNmShIQExo4dS79+/SgoKCA2NpbMzExGjRoFQHBwMDExMYwYMYIBAwbQqlUroqOjr2jxIiIiUj92f48/NDSU0NDQaq917tyZJUuW2OajoqKIioq66HqSkpLqWaKIiIg4iy6dFxERMRAFv4iIiIEo+EVERAxEwS8iImIgCn4REREDUfCLiIgYiIJfRETEQBT8IiIiBqLgFxERMRAFv4iIiIEo+EVERAxEwS8iImIgCn4REREDUfCLiIgYiIJfRETEQBT8IiIiBqLgFxERMRAFv4iIiIEo+EVERAzEvaELkGvf4q8WO7T86LtGO6kSERFxlPb4RUREDMTu4E9PT2fQoEH07duXiRMnUlpaWqPN6tWr6d+/Pw8//DAzZ86ksrISgLy8PGJiYhg8eDADBw5k2bJlTtsAERERsZ9dwV9YWEhcXByLFi1i06ZNeHl5kZiYWK3N4cOHSUxMZNWqVWzYsIEjR46QmpoKwNSpU+nRowdpaWl8+OGHrFq1ii+//NL5WyMiIiIXZVfwb9++nW7duhEYGAhAdHQ0aWlp1dqkp6cTFhZG8+bNcXd3JzIy0tZmwIABREZGAtC0aVOCgoI4fvy4M7dDRERE7GDXxX35+fn4+/vb5v38/MjLy6vWJi8vj6CgINt869atbW0eeeQR2+s7duxg//79TJ8+vc7+cnJycHXV5Qf1YTabMZlMDq0jN9cDCPh5Ohcfn3IAioqKHFqvo3U5o4YrqbS09Jqu70ZmpLF3xr8jZ3PG544RlJXl2qZPnMjF09PH4XWazWZycnIua1m7gv/8ufpfujCYrVbrJdts2LCBadOmMW/ePFq2bFlnf4GBgXh4eNhTmvzMZDJV+8Prcpw+/d/pgIAAzq/Ot9DXofU6WpczariSioqK8PW9duu7kRlp7J3x78jZnPG5YwRm82nOH+Ru0yYAb2/Hx8xkMhEQEEBxcXG9l7Ur+P39/Tl06JBt/sIjAOfbFBQU1Nlm/vz5pKSk8P777/Pb3/623oXK9cvRrwOKiIjz2HU8PSQkhL1795KdnQ1AcnIyYWFh1dqEhoaydetWiouLsVgspKSk2Nq89dZbbN68meTkZIW+iIhIA7Jrj79ly5YkJCQwduxYysvLad++PTNmzCAzM5P58+ezZMkSgoODiYmJYcSIEVRUVHDvvfcSHR3N2bNnSUxMpGXLlowaNcq2zqeeeorw8PArtmEiIiJSk9137gsNDSU0NLTaa507d2bJkiW2+aioKKKioqp34O7OwYMHHSxTREREnEGXzouIiBiI7tUvImIAeuaGnKc9fhEREQNR8IuIiBiIgl9ERMRAFPwiIiIGouAXERExEF3VLyJyHdCtr8VZtMcvIiJiINrjNwB79xRy/q8FEAnAmm/W8K/y+j/1SURErm3a4xcRETEQBb+IiIiBKPhFREQMROf4RUTkkmq7VqioqAjfQl+7lte9/q8d2uMXERExEAW/iIiIgSj4RUREDETBLyIiYiAKfhEREQPRVf1XgaP32NbVsCIi4iza4xcRETEQu/f409PTmTt3LmVlZXTr1o1p06bRqFGjam1Wr17N0qVLsVgshIWF8eKLL+Lq6kpJSQmvvPIK//73vwF49dVXue+++5y7JTcwe44Y1Of7tCIiYlx2BX9hYSFxcXGsXr2awMBApkyZQmJiIuPGjbO1OXz4MImJiaSmpuLt7U1sbCypqalERkYyb948fH192bhxIyaTiREjRvDJJ5/QtGnTK7ZhzqTHYYqIOMYZn6M67ekcdgX/9u3b6datG4GBgQBER0fzzDPPVAv+9PR0wsLCaN68OQCRkZF88MEHREZGkp6ezoIFCwAICgqic+fOpKenEx4eXq0fq9UKQHl5ueNb5kTu18GlEJ4ung7X2cjNlRYtym3T18N2XwucMfZyeTT2Detqj//7X71/1fqqzYguIy5ruYoKKy4uLWzTzsi4yspK23rOZ6e97PqN5efn4+/vb5v38/MjLy+vWpu8vDyCgoJs861bt7a1ycvLo3Xr1tV+dvLkyRr9VFRUAJCVlVWPTbjyutK1oUu4tBaOr6Lrb6Df5gM/zwX9/J9ckhPGXi6Txr5hGWz8Dxw4cOlGdfD23gzA999XAJe/nl8qLq56dHpFRQWenp52L2dX8FdWVtZ4zdW1+nWBtf3Fcb7NxX72S40bNyY4OBh3d3dcXFzsKU1ERMSQrFYrFRUVNG7cuF7L2RX8/v7+HDp0yDZ/4RGA820KCgpqbXP+Zy1atLD9rFu3bjX6cXV15aabbqrXBoiIiBhVffb0z7Pr63whISHs3buX7OxsAJKTkwkLC6vWJjQ0lK1bt1JcXIzFYiElJcXWJiwsjOTkZACys7PZv38/ISEh9S5WREREHONitfOqgIyMDObMmUN5eTnt27dnxowZ/Oc//2H+/PksWbIEqPqDYNmyZVRUVHDvvfcSFxeHu7s7ZrOZ+Ph4Dh8+jNVqZcKECfTp0+eKbpiIiIjUZHfwy7Vj9erVJCUl4eLiQosWLUhISKBNmzZMnz6dnTt3YrVaeeaZZxg8eHBDl3rDqW3sAwIC6NGjh+1bLwBxcXF07969ASu98SxZsoS///3vAHTp0oUpU6ZQWVmpe4RcJbWN/5kzZwgLC+PWW2+1tVuwYAG33HJLQ5V5Q1u+fDnJycmsW7fOofvj6Hsw15lvvvnGdr+E5s2bs2LFCiZPnkyfPn04efIkGzZs4NSpU0RFRdG1a9dq37QQx9Q19i+//DK33nora9asaegSb1h79uxh7dq1pKam0qhRI8aNG8fy5cspKCi4ru8Rcr2oa/yDgoJ46KGHmDdvXkOXeMM7ePAg77zzju1aOUfuj6Nb9l5nmjRpQkJCgu1+CZ07dyY3N5f09HQiIiJwdXWlRYsWhIWFsX79+gau9sZS19jv27eP0tJS/vjHPzJkyBBWrFjRwJXeeLp3787atWtp3LgxZ8+e5YcffsDHx4f09HQiIyOB6vcIEeeqa/z37dvHiRMnePTRR4mIiGDTpk0NXeoN6cyZM8THx/PCCy/YXnPkva/gv84EBQVx//33A1BWVsacOXPo169frfdKuPBeC+KYusbearXSu3dvPvjgA5YsWUJSUhJbt25t4GpvPB4eHiQnJ/PQQw/xww8/8Pvf/97ue4SI42obfw8PD/r3789HH33E7NmzmTZtGt98801Dl3rDmTx5MmPGjKFNmza21xx57yv4r1OnTp1i9OjRNG7cmHHjxtl9rwRx3IVj//jjj/PCCy/g7u6On58fjz32mPY6r5CoqCi++OILevXqxcSJE/W+v8ouHP/nn3+ep556ChcXF2699Vb69+9PRkZGQ5d5Q/nggw9o1apVjQviHXnv61/Idejo0aM89thj3Hbbbbz11lt4eHhc9D4K4jy1jX1KSgrff/+9rY3VasXDw6MBq7zx/Oc//7HdNc3FxYWhQ4eSlZWl9/1VUtf4L126lMLCQls7vfedLy0tjd27dzNkyBD+8pe/cOTIEaKjox167yv4rzMFBQU8/vjjDB8+nL/85S+4ubkBVfdKSElJwWKxUFxczNatWwkNDW3gam8sdY19VlYW77zzDlarldOnT5OSkkK/fv0auNoby/Hjx3nppZc4e/YsAB9//DE9evTQPUKukrrG/4svviApKQmAEydOsGnTJn1V28nWrFnDxx9/zLp163jttde49dZb+fDDDx167+uq/utMUlISp06dIiUlhZSUFAC8vLxISkrCZDIxePBgKioqGDVqFB06dGjgam8sdY39u+++S3x8PAMHDsRisTB8+HB69uzZwNXeWHr16kVkZCRRUVG4ubnRoUMH4uLicHNzs4291WplypQptquexXnqGv9z584RFxfHI488YvtqZbt27Rq6XEN49tlnL/u9r+/xi4iIGIgO9YuIiBiIgl9ERMRAFPwiIiIGouAXERExEAW/iIiIgSj4RUREDETBLyIiYiAKfhEREQNR8IuIiBjI/wepP2qI+yrzWwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 595x340 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(7, 4), dpi=85)\n",
    "ax.hist(boys_ppc['second']/2, bins=26, density=True, alpha=.4, color='g', label='half of secondborn boys')\n",
    "ax.axvline(x=boys_after_girl, color='y', label='observed boys after girl')\n",
    "ax.axvline(x=(n_boys - boys_after_girl - first_boys), color='b', label='observed boys after boy')\n",
    "ax.legend()\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Model the second birth **separately, conditioned** on the first gender"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "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: [p2_g, p2_b, p1]\n",
      "Sampling 2 chains: 100%|██████████| 12000/12000 [00:08<00:00, 1450.50draws/s]\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as model_h5:\n",
    "    p1 = pm.Uniform('p1', lower=0, upper=1)\n",
    "\n",
    "    boys1 = pm.Binomial('boys1', p=p1, n=first_births,\n",
    "                        observed=first_boys)\n",
    "\n",
    "    # two uninformed priors for second born\n",
    "    p2_b = pm.Uniform('p2_b', lower=0, upper=1)\n",
    "    p2_g = pm.Uniform('p2_g', lower=0, upper=1)\n",
    "\n",
    "    # use second _b, _g for inference\n",
    "    boys2_b = pm.Binomial('boys2_b', p=p2_b, n=(n_births - first_births - births_after_girl),\n",
    "                          observed=(n_boys - first_boys - boys_after_girl))\n",
    "    boys2_g = pm.Binomial('boys2_g', p=p2_g, n=births_after_girl,\n",
    "                          observed=boys_after_girl)\n",
    "\n",
    "    model_h5.trace = pm.sample(5000, tune=1000, chains=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "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>mean</th>\n",
       "      <th>sd</th>\n",
       "      <th>mc_error</th>\n",
       "      <th>hpd_2.5</th>\n",
       "      <th>hpd_97.5</th>\n",
       "      <th>n_eff</th>\n",
       "      <th>Rhat</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>p1</th>\n",
       "      <td>0.509666</td>\n",
       "      <td>0.048672</td>\n",
       "      <td>0.000434</td>\n",
       "      <td>0.412515</td>\n",
       "      <td>0.602438</td>\n",
       "      <td>13192.347632</td>\n",
       "      <td>0.999901</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>p2_b</th>\n",
       "      <td>0.415151</td>\n",
       "      <td>0.066587</td>\n",
       "      <td>0.000487</td>\n",
       "      <td>0.286184</td>\n",
       "      <td>0.541124</td>\n",
       "      <td>15413.593704</td>\n",
       "      <td>0.999970</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>p2_g</th>\n",
       "      <td>0.784031</td>\n",
       "      <td>0.057677</td>\n",
       "      <td>0.000500</td>\n",
       "      <td>0.666565</td>\n",
       "      <td>0.887501</td>\n",
       "      <td>16633.953240</td>\n",
       "      <td>0.999917</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          mean        sd  mc_error   hpd_2.5  hpd_97.5         n_eff      Rhat\n",
       "p1    0.509666  0.048672  0.000434  0.412515  0.602438  13192.347632  0.999901\n",
       "p2_b  0.415151  0.066587  0.000487  0.286184  0.541124  15413.593704  0.999970\n",
       "p2_g  0.784031  0.057677  0.000500  0.666565  0.887501  16633.953240  0.999917"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pm.summary(model_h5.trace)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 10000/10000 [00:13<00:00, 754.44it/s]\n"
     ]
    }
   ],
   "source": [
    "with model_h5:\n",
    "    ppc = pm.sample_posterior_predictive(model_h5.trace, samples=10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "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>boys1</th>\n",
       "      <th>boys2_b</th>\n",
       "      <th>boys2_g</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>5154</th>\n",
       "      <td>50</td>\n",
       "      <td>18</td>\n",
       "      <td>32</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>869</th>\n",
       "      <td>43</td>\n",
       "      <td>22</td>\n",
       "      <td>40</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7868</th>\n",
       "      <td>52</td>\n",
       "      <td>12</td>\n",
       "      <td>38</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1931</th>\n",
       "      <td>64</td>\n",
       "      <td>25</td>\n",
       "      <td>37</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4417</th>\n",
       "      <td>47</td>\n",
       "      <td>14</td>\n",
       "      <td>34</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      boys1  boys2_b  boys2_g\n",
       "5154     50       18       32\n",
       "869      43       22       40\n",
       "7868     52       12       38\n",
       "1931     64       25       37\n",
       "4417     47       14       34"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "boys_ppc = pd.DataFrame(ppc)\n",
    "boys_ppc.sample(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsIAAAEsCAYAAADJtNW2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X1cVHXe//H3wIShg81SaaFhglDt9iAFN+sy282bEqyulU1AAjNbLXdrC63UvEdBMTG3WlBTUcc7KN3dau3aLjcvNbsqpZBsM8MtXR02KR1jIECY+f3h1fxE0BEYhDyv51+cu+/5nMN3jm++njnH5Ha73QIAAAAMxq+tCwAAAADaAkEYAAAAhkQQBgAAgCERhAEAAGBIBGEAAAAYEkEYAAAAhuQ1CLtcLs2YMUOJiYlKTU3VoUOHGqxz/Phx3XPPPaqurpYk1dXVae7cuUpKSlJ8fLy2bdvm+8oBAACAFvAahLdu3aqamhrl5+dr4sSJmj9/fr3lO3fu1JgxY1RWVuaZ95e//EW1tbXauHGjcnNzGw3PAAAAQFsye1uhsLBQAwYMkCT17t1b+/btq7fcz89PeXl5+vWvf+2Z9+677yoiIkLjxo2T2+3W9OnTG20XAAAAuBhiYmIazPMahJ1OpywWi2fa399ftbW1MptPb9q/f/8G25w4cUKHDx/W0qVLtXv3bk2ZMkXr1q1rsN61117bpAPAhSkvL1dQUFBbl4Efmeb0m26vdJMkHR17tDVKQjvHtca4unULkSQdPWpv8rb0GzRHS/tNaWlpo/O9BmGLxaKKigrPtMvl8oTgc7FarfrlL38pk8mkW2+9VV999VWj64WEhHjbPZrBbrdzbtFkLek39Ddj4lqD5vz+6Tdojpb2m3MFYa/3CEdHR2vHjh2SpKKiIkVGRnrdWUxMjLZv3y5J2r9/PyO/AAAAaHe8jggPGTJEu3btUlJSktxutzIzM5WXl6fQ0FANGjSo0W0SEhI0c+ZMJSQkyO12a/bs2T4vHAAAAGgJr0HYz89P6enp9eaFh4c3WO+dd97x/BwQEKB58+b5oDwAAACgdfBCDQAAABgSQRgAAACGRBAGAACAIRGEAQAALlBqaqoOHjzYpjVs2LBBL730Ur15mzdv1sKFCy96LXv37tWQIUOUnZ2t3bt3a//+/T5rOyMjQ3Z7/WdVHzlyRAkJCT7bB0EYAAAAzbJz506NGjVKEydO1KZNm3Ts2DGftT116tRWf+a016dGAAAAtKVhw6QtW3zbZlyc9Ne/nnv5qVOnNGXKFB05ckR1dXV6+OGHFRcXJ0l68cUXdeLECQUEBGjBggWSpKeeekput1vV1dWaPXu2brrpJtlsNr355psymUyKi4vTqFGjNHnyZDkcDjkcDvXs2VO33nqrhg8frrKyMj366KPavHmzsrOztWfPHrlcLo0ePVqxsbHas2ePMjMz1blzZ/n7+6t3794Nai4qKtJDDz0kp9OpJ554Qr/85S+1a9cuLV68WB06dJDValVmZqaWL1+url276sEHH9TJkyf18MMPa/ny5Y0eww+cTqemTp2q8vJyHTt2TMnJybr55pu1efNmXXbZZTp16pR27typTz/9VL169dLevXu1atUq+fn5KSYmRk8//bReeuklffzxx6qsrFRGRobnKWRVVVV69tlndezYMV177bXavXu33n33XaWmpmrWrFnasmWL3n//fdXV1SkjI8OHvYAgDAAAmmnZsqZv43B0lNUqjRvn+3p8KT8/X8HBwVq4cKGcTqfi4+N12223SZLuvvtuDRs2TOvWrdPSpUt1++23y2q1asGCBSopKVFlZaVKSkq0ZcsWrV+/XpL08MMP64477pAk3XbbbRo9erRKSkqUnp6u4cOH6y9/+Yvi4+O1fft2HTlyRBs2bFB1dbUSEhLUv39/zZ49Wy+++KJ69uypmTNnNlpzYGCgli1bpuPHj2vEiBEaMGCApk+frg0bNqhr165avXq1cnNzlZycrAkTJujBBx/Um2++qfvuu0/FxcUNjuFMhw4d0rBhw3T33Xfr66+/Vmpqqt5++20NHz5cV111lUaOHKkDBw4oLi5OHTt21EsvvaRNmzYpMDBQzzzzjHbt2iVJCgsL07Rp0xqc6+7du+vFF1/UwYMHde+99zY4ttDQUM2bN09Hjhxp2S/2LARhAADQrp1v5La1HDx4UP/xH/8hSbJYLAoPD9e//vUvSVLfvn0lnX777vbt2zVp0iR99dVX+u1vfyuz2azx48frwIEDstvtGj16tCTp5MmTOnTokCSpZ8+ekqRevXqprq5OR48e1ZYtW7Rq1Srl5+fr008/VWpqqiSptrZWR48e1TfffOPZLjo6WocPH25Qc0xMjEwmk6688koFBQXp5MmTslgs6tq1qyTp5z//uRYtWqTrrrtOnTp1UklJid544w3l5OTIarU2OIYzXXXVVVq9erXefvttWSwW1dbWnvPcHT58WMePH9e4//trp6KiwlPvD8dw9rm+8847JZ1+V0VwcHCDda677rpz7q8luEcYAADgLOHh4dqzZ4+k07cFHDhwQN27d5ckffLJJ5KkPXv2KCIiQh988IG6dOmilStXavz48Vq0aJHCwsLUq1cvrVmzRjabTfHx8brhhhskSSaTybOfBx54QM8//7x69eqlzp07KywsTP369ZPNZtPq1asVGxur6667Tl27dvV8Se+H/Z/th/llZWWqrKzUT37yEzmdTs99ux9++KGuv/56SaffApyTk6OuXbsqODi40WM408qVK9W7d28tXLhQQ4cOldvtbrB/k8kkt9ut7t2769prr9XKlStls9mUkpLiuZXDz69h9IyMjNTHH38s6XSIPnHiRIN1GtvOFxgRBgAAOEtCQoKmT5+ukSNHqrq6Wo8//riuvPJKSdLWrVu1evVqderUSVlZWXK5XJowYYI2bNig2tpa/e53v9ONN96o22+/XSNHjlRNTY2ioqI8I7NnGjp0qDIyMpSbmytJGjhwoD788EMlJyersrJSgwcPlsViUXp6up599llZLBZ16tRJV1xxRYO2qqqqNGrUKFVWVio9PV0mk0lz587VE088IZPJpCuuuMLz5t/BgwcrPT1dzz//vCTpxhtvbHAMZ7rrrrs0d+5cbdmyRUFBQfL391dNTU29dW655RYtXLhQixcv1ujRo5Wamqq6ujp169ZNsbGx5zzXDzzwgCZPnqwHH3xQISEh6tChQxN+Uy1jcjcW6S+CwsJCxcTEtMWuL3l2u73Vv2WJS09z+o1p9ulRDffMNrmMoI1xrTGuHwY0ly5t+rYOh0NWq7Xd3yN8qfv++++VkpKiV199tdVGWy/URx99pMrKSt1xxx366quv9Jvf/EZbt26tt05Lrzfnyp2MCAMAABjIRx99pJkzZ+p3v/tdm4dg6fT9vxMmTNDLL7+s2tpazZgx46LtmyAMAABgINHR0XrjjTfaugyPq6++WjabrU323fZ/BgAAAABtgCAMAAAAQyIIAwAAwJAIwgAAADAkgjAAAMAFSk1N9bzYoq1s2LBBL730Ur15mzdv1sKFCy96LXv37tWQIUOUnZ2t3bt3a//+/U1uo61qlwjCAAAAaKadO3dq1KhRmjhxojZt2uR5i92PBY9PAwAA7dqw9cO05YstPm0zLiJOf03+6zmXnzp1SlOmTNGRI0dUV1enhx9+WHFxcZKkF198USdOnFBAQIAWLFggSXrqqafkdrtVXV2t2bNn66abbpLNZtObb74pk8mkuLg4jRo1SpMnT5bD4ZDD4VDPnj116623avjw4SorK9Ojjz6qzZs3Kzs7W3v27JHL5dLo0aMVGxurPXv2KDMzU507d5a/v7/nlcVnKioq0kMPPSSn06knnnhCv/zlL7Vr1y4tXrxYHTp0kNVqVWZmppYvX66uXbvqwQcf1MmTJ/Xwww9r+fLljR7DD5xOp6ZOnary8nIdO3ZMycnJuvnmm7V582ZddtllOnXqlHbu3KlPP/1UvXr10t69e7Vq1Sr5+fkpJiZGTz/9tF566SV9/PHHqqysVEZGhsLDwy+49ssvv1yLFi1qtPbNmzc3ux8QhAEAAM6Sn5+v4OBgLVy4UE6nU/Hx8brtttskSXfffbeGDRumdevWaenSpbr99ttltVq1YMEClZSUqLKyUiUlJdqyZYvWr18vSXr44Yd1xx13SJJuu+02jR49WiUlJUpPT9fw4cP1l7/8RfHx8dq+fbuOHDmiDRs2qLq6WgkJCerfv79mz56tF198UT179tTMmTMbrTkwMFDLli3T8ePHNWLECA0YMEDTp0/Xhg0b1LVrV61evVq5ublKTk7WhAkT9OCDD+rNN9/Ufffdp+Li4gbHcKZDhw5p2LBhuvvuu/X1118rNTVVb7/9toYPH66rrrpKI0eO1IEDBxQXF6eOHTvqpZde0qZNmxQYGKhnnnlGu3btkiSFhYVp2rRpTa79xRdfPGftLUEQBgAA7dr5Rm5by8GDB/Uf//EfkiSLxaLw8HD961//kiT17dtX0ukXU2zfvl2TJk3SV199pd/+9rcym80aP368Dhw4ILvdrtGjR0uSTp48qUOHDkmSevbsKUnq1auX6urqdPToUW3ZskWrVq1Sfn6+Pv30U6WmpkqSamtrdfToUX3zzTee7aKjo3X48OEGNcfExMhkMunKK69UUFCQTp48KYvFoq5du0qSfv7zn2vRokW67rrr1KlTJ5WUlOiNN95QTk6OrFZrg2M401VXXaXVq1fr7bfflsViUW1t7TnP3eHDh3X8+HGN+7/3aFdUVHjq/eEYmlr7LbfcIpvN1mjtLcE9wgAAAGcJDw/Xnj17JJ2+LeDAgQPq3r27JOmTTz6RJO3Zs0cRERH64IMP1KVLF61cuVLjx4/XokWLFBYWpl69emnNmjWy2WyKj4/XDTfcIEkymUye/TzwwAN6/vnn1atXL3Xu3FlhYWHq16+fbDabVq9erdjYWF133XXq2rWr50t6P+z/bD/MLysrU2VlpX7yk5/I6XR67tv98MMPdf3110uSEhISlJOTo65duyo4OLjRYzjTypUr1bt3by1cuFBDhw6V2+1usH+TySS3263u3bvr2muv1cqVK2Wz2ZSSkuK5leNcr3T2VvvevXvPWXtLMCIMAABwloSEBE2fPl0jR45UdXW1Hn/8cV155ZWSpK1bt2r16tXq1KmTsrKy5HK5NGHCBG3YsEG1tbX63e9+pxtvvFG33367Ro4cqZqaGkVFRXlGN880dOhQZWRkKDc3V5I0cOBAffjhh0pOTlZlZaUGDx4si8Wi9PR0Pfvss7JYLOrUqZOuuOKKBm1VVVVp1KhRqqysVHp6ukwmk+bOnasnnnhCJpNJV1xxhebNmydJGjx4sNLT0/X8889Lkm688cYGx3Cmu+66S3PnztWWLVsUFBQkf39/1dTU1Fvnlltu0cKFC7V48WKNHj1aqampqqurU7du3RQbG3ve8+2t9g4dOuiFF15otPaWMLkbi/QXQWFhoWJiYtpi15c8u92ukJCQti4DPzLN6Tem2adHNdwz2+QygjbGtca4fhjQXLq06ds6HA5ZrVb93/+ao418//33SklJ0auvvnrOUdr25MzrTXNqP1fubP9HDgAAAJ/56KOPlJCQoLFjx/4oQvCZfF271xZcLpdmzJihxMREpaamem70PtPx48d1zz33qLq6ut78gwcPKiYmpsF8AAAAtI3o6Gi98cYbGjp0aFuX0mS+rt1rEN66datqamqUn5+viRMnav78+fWW79y5U2PGjFFZWVm9+U6nU1lZWQoICPBJoQAAAIAveQ3ChYWFGjBggCSpd+/e2rdvX/0G/PyUl5cnq9Xqmed2uzV9+nRNmDBBgYGBPi4ZAAAAaDmvT41wOp2yWCyeaX9/f9XW1spsPr1p//79G2zz8ssv6xe/+IVuvPHG87Ztt9ubWi8uQHl5OecWTdaSfkN/MyauNUZ2+ktLDoejyVtWVVXJ4XDIbq/0vjLwf1rreuM1CFssFlVUVHimXS6XJwSfy+uvv65rrrlGmzZtUllZmcaMGaN169Y1WI9vG7cOvsmN5mhJv6G/GRPXGpz5v8EX6oenRoSENH1bGFdLrzelpaWNzvcahKOjo7Vt2zbFxcWpqKhIkZGRXnf23//9356fBw4cqJUrVzahVAAAgPqWLfNte+3h8W2pqamaNWuWwsPDva6blpampKQk9evXzzNv8uTJiouL05133tmaZTawdu1arVu3Tk888YTKy8sVHx+vyy67zCdtP/7443r55Zfrzdu8ebOKi4s1a9Ysn+zjTF6D8JAhQ7Rr1y4lJSXJ7XYrMzNTeXl5Cg0N1aBBg3xeEAAAANqvt99+W4sXL9YNN9yggQMH6le/+pXP2j47BLc2r0HYz89P6enp9eY19pfLO++80+j255oPAADQXn355ZeaMmWKzGazXC6XsrOz1aVLF82YMUP//ve/dezYMQ0cOFBpaWmaPHmyzGaz7Ha7ampqFBcXp23btqm0tFQ5OTkqLS3VkiVL5Ofnp7KyMiUmJurBBx/07Ku8vFxTp07ViRMnJEnTpk3TDTfcoHXr1unVV1/V1VdfrW+//bbROtevX68VK1aorq5OGRkZ6tGjh1auXKm//vWvMpvN6tu3r5555hklJSVpzpw5ioiI0Pbt27Vt2zbdd999ysrKktlsVmBgoP7whz/U+17YgQMHNH/+fNXV1enEiROaNWuWvvjiC/3jH//Q1KlTlZiYqLKyMqWlpSknJ0fZ2dnas2ePXC6XRo8erdjYWKWmpio4OFgnT57UihUr5O/vL0k6dOiQ57x169ZNR48elc1mU//+/bVr16562w0bNqzVfs8/rqcoAwAAXATvvfeeoqKilJeX57kFoLS0VL1799aKFSv02muvaePGjZ71u3XrppUrVyosLExHjhzRK6+8orvvvtszIPj1118rNzdXBQUFWrVqVb1gu2TJEt12222y2WyaM2eOZs2apW+++UZr1qxRQUGBcnJydOrUqUbrjI6O1urVqzV27Fg9//zz+vzzz/XWW29p48aN2rhxow4dOqRt27ZpxIgR+tOf/iRJ2rRpk0aMGKGtW7cqNjZWa9eu1ciRI/Xdd9/Va7ukpESTJk3ytL9582YlJibqpptuUlZWlkaMGKGrr75aL7zwgrZv364jR45ow4YNWrNmjZYsWeJp795779WqVas8IViSFixYoMcee0w2m03R0dGNHltj2/ma1xFhAAAAo3nggQf0yiuv6De/+Y2CgoKUlpYmq9WqTz75RO+//74sFotqamo86//0pz+VJHXu3FlhYWGen39Yp0+fPp53K0REROjw4cOebQ8cOKD3339fb731liTp5MmTOnz4sHr16uXZJioqqtE6+/bt62l/wYIF+uc//6lbbrnFc89u37599cUXXyglJUXx8fF65JFH9PXXX+tnP/uZunfvriVLluihhx5S165dG+yjS5cuysnJ0eWXX66Kiop6o8VnO3DggD799FOlpqZKkmpra3X06FFJUs+ePRusf/DgQfXp00eSFBMTozfeeKPBOo1t52uMCAMAAJzl73//u2JiYrR69WoNHTpUy5cv1+bNmxUUFKTs7GyNGTNGVVVVcrvdkiSTyXTe9j777DPV1dXp+++/V0lJiXr06OFZFhYWptGjR8tms2nx4sW6//77df3116ukpERVVVWqq6vTZ5991mi7xcXFkqQ9e/YoIiJCYWFhKi4uVm1trdxut3bv3q2ePXuqY8eO6tevnzIyMnT//fdLOv2Ur+HDh8tmsykiIkIFBQX12s7IyNDvf/97ZWVlKTIy0nOsZzKZTHK5XAoLC1O/fv1ks9m0evVqxcbG6rrrrjvnuYmMjNTHH38sSdq7d2+jx+btnPoCI8IAAABnufnmmzVp0iTl5ubK5XJpypQpCggI0MSJE1VUVKSAgAD16NFDx44du6D2amtrNXbsWDkcDo0fP17BwcGeZY899pimTp2qgoICOZ1OPf744woODtbYsWOVlJSk4ODgc76gbO/evRo1apRMJpMyMzPVrVs3xcbGauTIkXK5XIqJidHgwYMlSQkJCUpOTvY8fSEqKkrTpk1TYGBgo98Ju//++/Xkk0+qc+fOuuaaazz3MJ+pb9++GjdunNasWaMPP/xQycnJqqys1ODBg887gvz000/rueee08qVKxUUFOT10bytxeRuLN5fBIWFhYqJiWmLXV/yeLYnmqM5/cY0+/Rf6+6ZbXIZQRvjWmNcPwzULV3a9G1/eI5we3h82cXywQcfaOPGjXrhhRfatI7i4mKtXbtWCxYsaNM6pNOj0bfccot69OihV199VR999JHmzZt3zvVber05V+5kRBgAAOASt3btWr322mtavHhxW5ciSbr22muVlpbmGY3OzMxskzoIwgAAAK2oX79+9V6E0RZSUlKUkpLSpjWc6ec//7k2b97c1mXwZTkAAAAYE0EYAAAAhkQQBgAAgCERhAEAAGBIfFkOAAC0e8sKl/m0vXExbf/8ttTUVM2aNUvh4eFe101LS1NSUlK9L91NnjxZcXFxuvPOO1uzzAbWrl2rdevWeV49HR8f73mT3YVqq9rPxogwAAAALtjbb7+txYsXKy4uTkuXLpXL5WrrkpqNEWEAAICzfPnll5oyZYrMZrNcLpeys7PVpUsXzZgxQ//+97917NgxDRw4UGlpaZo8ebLMZrPsdrtqamoUFxenbdu2qbS0VDk5OSotLdWSJUvk5+ensrIyJSYm6sEHH/Tsq7y8XFOnTvW8uW3atGm64YYbtG7dOr366qu6+uqr9e233zZa5/r167VixQrV1dUpIyNDPXr00MqVK/XXv/5VZrNZffv21TPPPKOkpCTNmTNHERER2r59u7Zt26b77rtPWVlZMpvNCgwM1B/+8Id6b4M7cOCA5s+fr7q6Op04cUKzZs3SF198oX/84x+aOnWqEhMTVVZWprS0NOXk5Cg7O1t79uyRy+XS6NGjFRsbq9TUVAUHB+vkyZNasWKF/P39m1X7E088oZCQEE/tP7wdr6UYEQYAADjLe++9p6ioKOXl5XluASgtLVXv3r21YsUKvfbaa9q4caNn/W7dumnlypUKCwvTkSNH9Morr+juu+/WO++8I0n6+uuvlZubq4KCAq1atapesF2yZIluu+022Ww2zZkzR7NmzdI333yjNWvWqKCgQDk5OTp16lSjdUZHR2v16tUaO3asnn/+eX3++ed66623tHHjRm3cuFGHDh3Stm3bNGLECP3pT3+SJG3atEkjRozQ1q1bFRsbq7Vr12rkyJH67rvv6rVdUlKiSZMmedrfvHmzEhMTddNNNykrK0sjRozQ1VdfrRdeeEHbt2/XkSNHtGHDBq1Zs0ZLlizxtHfvvfdq1apV9UJwU2v/29/+Vq92X2FEGAAA4CwPPPCAXnnlFf3mN79RUFCQ0tLSZLVa9cknn+j999+XxWJRTU2NZ/2f/vSnkqTOnTsrLCzM8/MP6/Tp00cBAQGSpIiICB0+fNiz7YEDB/T+++/rrbfekiSdPHlShw8fVq9evTzbREVFNVpn3759Pe0vWLBA//znP3XLLbd47tnt27evvvjiC6WkpCg+Pl6PPPKIvv76a/3sZz9T9+7dtWTJEj300EPq2rVrg3106dJFOTk5uvzyy1VRUVFvtPhsBw4c0KeffqrU1FRJUm1trY4ePSpJ6tmzZ4trz83N1bfffuup3VcYEQYAADjL3//+d8XExGj16tUaOnSoli9frs2bNysoKEjZ2dkaM2aMqqqq5Ha7JUkmk+m87X322Weqq6vT999/r5KSEvXo0cOzLCwsTKNHj5bNZtPixYt1//336/rrr1dJSYmqqqpUV1enzz77rNF2i4uLJUl79uxRRESEwsLCVFxcrNraWrndbu3evVs9e/ZUx44d1a9fP2VkZOj++++XJL3++usaPny4bDabIiIiVFBQUK/tjIwM/f73v1dWVpYiIyM9x3omk8kkl8ulsLAw9evXTzabTatXr1ZsbKyuu+66856bptTep0+ferX7CiPCAAAAZ7n55ps1adIk5ebmyuVyacqUKQoICNDEiRNVVFSkgIAA9ejRQ8eOHbug9mprazV27Fg5HA6NHz9ewcHBnmWPPfaYpk6dqoKCAjmdTj3++OMKDg7W2LFjlZSUpODgYAUGBjba7t69ezVq1CiZTCZlZmaqW7duio2N1ciRI+VyuRQTE6PBgwdLkhISEpScnOy5vzYqKkrTpk1TYGCg/Pz8lJ6eXq/t+++/X08++aQ6d+6sa665xnMP85n69u2rcePGac2aNfrwww+VnJysyspKDR48+LwjyE2tfdiwYXryySd9dm/wD0zuxuL9RVBYWKiYmJi22PUlz263KyQkpK3LwI9Mc/qNafbpv/LdM9vkMoI2xrXGuH4Y4Fu6tOnbOhwOWa1WjWv7p5ddNB988IE2btyoF154oU3rKC4u1tq1a7VgwYI2raM53nnnHf3Xf/1Xs2s/V+5kRBgAAOASt3btWr322mtavHhxW5fSZGvXrtWGDRv0xz/+0edtE4QBAABaUb9+/eq9CKMtpKSkKCUlpU1raK6UlBQNHDiwVf4Hii/LAQAAwJAIwgAAADAkgjAAAAAMiSAMAAAAQyIIAwAAwJC8BmGXy6UZM2YoMTFRqampOnToUIN1jh8/rnvuuUfV1dWSpPLycj322GNKSUlRYmKiPv74Y99XDgAAALSA1yC8detW1dTUKD8/XxMnTtT8+fPrLd+5c6fGjBmjsrIyz7y8vDzddtttWrt2rebNm9fgTSUAAABAW/P6HOHCwkINGDBAktS7d2/t27ev3nI/Pz/l5eXp17/+tWfe6NGjFRAQIEmqq6tThw4dfFkzAAAA0GJeg7DT6az3rmh/f3/V1tbKbD69af/+/Rts07lzZ0lSWVmZnnnmGT333HONtm2325tVNM6vvLycc4sma0m/ob8ZE9caIzv9YgOHw9HkLauqquRwOGS3V/q6KFzCWut64zUIWywWVVRUeKZdLpcnBJ/P559/rgkTJujZZ5/Vrbfe2ug6vKO+ddjtds4tmqwl/Yb+Zkwn7BkVAAAXbUlEQVRca2C1Wpu8jcPhkNVqVUhI07eFcbX0elNaWtrofK/3CEdHR2vHjh2SpKKiIkVGRnrdWUlJiZ588kllZ2frF7/4RRNLBQAAAFqf16HdIUOGaNeuXUpKSpLb7VZmZqby8vIUGhqqQYMGNbpNdna2ampqlJGRIen0qHJubq5vKwcAAABawGsQ9vPza/DUh/Dw8AbrvfPOO56fCb0AAABo73ihBgAAAAyJIAwAAABDIggDAADAkAjCAAAAMCSCMAAAAAyJIAwAAABDIggDAADAkAjCAAAAMCSCMAAAAAyJIAwAAABDIggDAADAkAjCAAAAMCRzWxcA4MdnWeGy8043x7iYcS1uAwCApmBEGAAAAIZEEAYAAIAhEYQBAABgSARhAAAAGBJBGAAAAIZEEAYAAIAhEYQBAABgSARhAAAAGBJBGAAAAIZEEAYAAIAhEYQBAABgSARhAAAAGJLXIOxyuTRjxgwlJiYqNTVVhw4darDO8ePHdc8996i6ulqSVFVVpSeeeELJyckaO3asjh8/7vvKAQAAgBbwGoS3bt2qmpoa5efna+LEiZo/f3695Tt37tSYMWNUVlbmmbdhwwZFRkZq/fr1+tWvfqWcnBzfVw4AAAC0gNcgXFhYqAEDBkiSevfurX379tVvwM9PeXl5slqtjW5z55136n//9399WTMAAADQYmZvKzidTlksFs+0v7+/amtrZTaf3rR///6NbhMUFCRJ6tSpk8rLy31VLwAAAOATXoOwxWJRRUWFZ9rlcnlC8IVsU1FRoc6dOze6nt1ub0qtuEDl5eWcWzRZU/qNw+E473Rz0Gd/fLjWGFmIpOZ99quqquRwOGS3V/q6KFzCWut64zUIR0dHa9u2bYqLi1NRUZEiIyO9NhodHa3t27crKipKO3bsUExMTKPrhYSENL1ieGW32zm3aLKm9BtrqbX+tNV6jjUvHH32x4drDZrz2Xc4HLJarQoJafl1A8bR0utNaWlpo/O9BuEhQ4Zo165dSkpKktvtVmZmpvLy8hQaGqpBgwY1us3IkSM1adIkjRw5Updddpmys7ObXTgAAADQGrwGYT8/P6Wnp9ebFx4e3mC9d955x/NzYGCgXnzxRR+UBwAAALQOXqgBAAAAQ/I6Igzg4llWuMznbY6LGefzNgEAuBQwIgwAAABDIggDAADAkAjCAAAAMCSCMAAAAAyJIAwAAABDIggDAADAkAjCAAAAMCSCMAAAAAyJIAwAAABDIggDAADAkAjCAAAAMCSCMAAAAAyJIAwAAABDIggDAADAkAjCAAAAMCSCMAAAAAzJ3NYFAADQ3ixb5pt2xo3zTTsAWgcjwgAAADAkgjAAAAAMiSAMAAAAQ+IeYQCXrGWFPrrR8wzjYrjpEwAuFYwIAwAAwJAIwgAAADAkgjAAAAAMyes9wi6XS7NmzdLnn3+ugIAAzZ07Vz169PAsLygo0MaNG2U2mzV+/HjdddddstvtevbZZ+V2u3XFFVcoOztbgYGBrXogAAAAQFN4HRHeunWrampqlJ+fr4kTJ2r+/PmeZWVlZbLZbNq4caNWrFihRYsWqaamRqtWrVJsbKzWrVuniIgIvfbaa616EAAAAEBTeQ3ChYWFGjBggCSpd+/e2rdvn2dZcXGx+vTpo4CAAAUFBSk0NFT79+/XTTfdpO+++06S5HQ6ZTbzcAoAAAC0L14TqtPplMVi8Uz7+/urtrZWZrNZTqdTQUFBnmWdOnWS0+nUNddco+zsbL355puqqanR448/3mjbdrvdB4eAs5WXl3Nuf6QcDofP27zQvtCUfnN2nb6ouzX6bFueTyO4lK81DkdHn7Rjt1f6pJ32J0RS8z5jVVVVcjgcl/C5QWtoreuN1yBssVhUUVHhmXa5XJ4R3rOXVVRUKCgoSDNmzNC8efM0YMAA/c///I8mTZqkZY28uD0kJMQXx4Cz2O12zu2PlLXU6vM2L7QvNKXfnF2n1dryulujz7bl+TSCS/la44MuLUkKCfF9H2xPmvPZdzgcslqtl/y5gW+19HpTWlra6Hyvt0ZER0drx44dkqSioiJFRkZ6lkVFRamwsFDV1dUqLy/XwYMHFRkZqc6dO3tGirt06eK5TQIAAABoL7yOCA8ZMkS7du1SUlKS3G63MjMzlZeXp9DQUA0aNEipqalKTk6W2+1WWlqaOnTooOnTpys9PV0ul0tut1szZsy4GMcCAAAAXDCvQdjPz0/p6en15oWHh3t+TkhIUEJCQr3lvXr10po1a3xUIgAAAOB7vFADAAAAhkQQBgAAgCERhAEAAGBIBGEAAAAYEq98AwDAIBp5pD9gaIwIAwAAwJAIwgAAADAkbo0AAAAXna9u0xg3zjftwJgYEQYAAIAhEYQBAABgSARhAAAAGBJBGAAAAIZEEAYAAIAhEYQBAABgSARhAAAAGBJBGAAAAIZEEAYAAIAh8WY5AADaOV+9hQ1AfYwIAwAAwJAIwgAAADAkgjAAAAAMiSAMAAAAQyIIAwAAwJAIwgAAADAkHp8GAAB+tHz1aLlx43zTDn5cGBEGAACAIXkNwi6XSzNmzFBiYqJSU1N16NChessLCgoUHx+vhIQEbdu2TZJUWVmpZ599VsnJyRoxYoSKi4tbp3oAAACgmbzeGrF161bV1NQoPz9fRUVFmj9/vnJzcyVJZWVlstls2rRpk6qrq5WcnKz+/ftrxYoVioiI0IIFC7R//37t379fUVFRrX4wAAAAwIXyOiJcWFioAQMGSJJ69+6tffv2eZYVFxerT58+CggIUFBQkEJDQ7V//369++67uuyyy/TII48oJyfHsz0AAADQXngdEXY6nbJYLJ5pf39/1dbWymw2y+l0KigoyLOsU6dOcjqdOnHihL777jutWLFCf/7zn5WVlaUFCxY0aNtut/voMHCm8vJyzu1Z1n621udtptyU4vM2HQ6Hz9u80L7QlH5zdp2+qLs1+mxbnk8juJSvNQ5HR5+0Y7dX+qQdX9XjO1ZJzfuMVVVVtcpns6V89btC62it643XIGyxWFRRUeGZdrlcMpvNjS6rqKhQUFCQrFarBg4cKEm66667tOwcX+kMCQlpUfFonN1u59yexVpq9XmbrXGO27LOpvSbs+u0Wlte96V2Po3gUr7W+KBLS5JCQnzTkK/q8bXmfPYdDodPrhm+5qvfFVpHS683paWljc73emtEdHS0duzYIUkqKipSZGSkZ1lUVJQKCwtVXV2t8vJyHTx4UJGRkYqJidH27dslSbt371avXr2aXTgAAADQGryOCA8ZMkS7du1SUlKS3G63MjMzlZeXp9DQUA0aNEipqalKTk6W2+1WWlqaOnTooEcffVTTpk1TYmKizGazsrKyLsaxAAAAtCmea/zj4jUI+/n5KT09vd688PBwz88JCQlKSEiot9xqterll1/2UYkAAACA7/FmOQAAWomvRgcBtA7eLAcAAABDIggDAADAkAjCAAAAMCSCMAAAAAyJIAwAAABDIggDAADAkAjCAAAAMCSCMAAAAAyJIAwAAABDIggDAADAkHjFMnCJW1Z4Ye94dTgcspZaW7kaAADaD0aEAQAAYEiMCANoFy505BoAAF8hCAMAAMNbxt/ihsStEQAAADAkgjAAAAAMiSAMAAAAQyIIAwAAwJAIwgAAADAkgjAAAAAMicenAUAba41nKI+LGefzNgHgUsOIMAAAAAyJIAwAAABDIggDAADAkAjCAAAAMCSvQdjlcmnGjBlKTExUamqqDh06VG95QUGB4uPjlZCQoG3bttVb9uGHH+oXv/iFbysGAAAAfMDrUyO2bt2qmpoa5efnq6ioSPPnz1dubq4kqaysTDabTZs2bVJ1dbWSk5PVv39/BQQEqLS0VHl5eaqtrW31gwDaQmt80x8AAFw8XkeECwsLNWDAAElS7969tW/fPs+y4uJi9enTRwEBAQoKClJoaKj279+v6upqzZw5U7NmzWq1wgEAAICW8Doi7HQ6ZbFYPNP+/v6qra2V2WyW0+lUUFCQZ1mnTp3kdDqVnp6uMWPGqGvXrudt2263t6B0nEt5eTnn9iwOh6OtS2j3qqqqmn2ejHR+W+Oz1Rrn72JcAy7la43D0bGtS2jnrJKa13dbcq0xEru9sq1LaFda63rjNQhbLBZVVFR4pl0ul8xmc6PLKioqdNlll2nPnj06fPiw/vjHP+rkyZNKS0vTCy+80KDtkJAQXxwDzmK32zm3Z7GWWtu6hHbP4XDIam3eeWrudj9GrfHZao3+eTGuAZfytcZAXbpFmvPZb8m1xkhCQjhHZ2rp9aa0tLTR+V6DcHR0tLZt26a4uDgVFRUpMjLSsywqKkqLFy9WdXW1ampqdPDgQUVFRelvf/ubZ53+/fs3GoIBAACAtuQ1CA8ZMkS7du1SUlKS3G63MjMzlZeXp9DQUA0aNEipqalKTk6W2+1WWlqaOnTocDHqBgAAAFrEaxD28/NTenp6vXnh4eGenxMSEpSQkHDO7Xft2tWC8gAAAIDWwQs1AAAAYEgEYQAAABgSQRgAAACGRBAGAACAIRGEAQAAYEgEYQAAABgSQRgAAACGRBAGAACAIXl9oQYA4P9bVrisrUsAAPgII8IAAAAwJIIwAAAADIkgDAAAAEMiCAMAAMCQCMIAAAAwJJ4aAQC4ZCzjoR4AmoARYQAAABgSQRgAAACGRBAGAACAIRGEAQAAYEgEYQAAABgST41Aiywr9P1XtMfFjPN5mwAAAGdjRBgAAACGRBAGAACAIRGEAQAAYEgEYQAAABiS1y/LuVwuzZo1S59//rkCAgI0d+5c9ejRw7O8oKBAGzdulNls1vjx43XXXXfJbrfrueeeU11dndxut9LT0xUWFtaqBwIAAAA0hdcgvHXrVtXU1Cg/P19FRUWaP3++cnNzJUllZWWy2WzatGmTqqurlZycrP79++sPf/iDUlJSNHjwYO3cuVOLFi3Syy+/3OoHAwD4cVrm+wfQAIBXXoNwYWGhBgwYIEnq3bu39u3b51lWXFysPn36KCAgQAEBAQoNDdX+/fs1adIkBQUFSZLq6urUoUOHViofl6LWeCQbAADA2bwGYafTKYvF4pn29/dXbW2tzGaznE6nJ/BKUqdOneR0OhUcHCxJ+uc//6msrCz98Y9/bIXSAQAAgObzGoQtFosqKio80y6XS2azudFlFRUVnmD8/vvva/bs2VqwYME57w+22+0tKh6NKy8vv2jn1uFwXJT9oPVVVVU1+/dJP2h/fH0NWPvZ2gbzqqqqdPlnlze7zZSbUjw/Oxwdm90O2oJVUvM++y251hiJ3V7Z1iW0K62VbbwG4ejoaG3btk1xcXEqKipSZGSkZ1lUVJQWL16s6upq1dTU6ODBg4qMjNT777+vjIwMLV++XN26dTtn2yEhIb45CtRjt9sv2rm1llovyn7Q+hwOh6zW5v0+m7sdWo+vrwGNfdZb0mek+jXShX6cmvP7b2m/MYqQEM7RmVqabUpLSxud7zUIDxkyRLt27VJSUpLcbrcyMzOVl5en0NBQDRo0SKmpqUpOTpbb7VZaWpo6dOigzMxMnTp1SpMnT5Yk9ezZU+np6c0uHgAAAPA1r0HYz8+vQYgNDw/3/JyQkKCEhIR6y19//XUflQcAAAC0Dl6oAQAAAEMiCAMAAMCQvN4aAQC4NDXlJRY7GvkC+/eVAQrsKN05wHc1AcDFRBAGgEvQhbyYprFwCwBGQhAGALTIjp3N3LDQp2UAQJNxjzAAAAAMiSAMAAAAQyIIAwAAwJAIwgAAADAkgjAAAAAMiSAMAAAAQ+LxaQAAAO1MU154cz7jxvmmnUsVI8IAAAAwJIIwAAAADIkgDAAAAEMiCAMAAMCQCMIAAAAwJIIwAAAADInHpxnIskIfPYsFAADgEkAQBgC0iR2Vvv/j/M6OPDQVwIXj1ggAAAAYEiPC7VRLbmNwOByyllp9WA0AAPgx4g1150cQBoAfmR0727qC9ovbLQA0BUEYAC4SAiwAtC/cIwwAAABDIggDAADAkLzeGuFyuTRr1ix9/vnnCggI0Ny5c9WjRw/P8oKCAm3cuFFms1njx4/XXXfdpePHj+vpp59WVVWVunTponnz5ikwMLBVDwQAWgu3NADApclrEN66datqamqUn5+voqIizZ8/X7m5uZKksrIy2Ww2bdq0SdXV1UpOTlb//v2Vk5Oje++9V/Hx8Vq2bJny8/M1evTo1j4WAKiHAAtf4At4gG+ePtEenzzhNQgXFhZqwIABkqTevXtr3759nmXFxcXq06ePAgICFBAQoNDQUO3fv1+FhYV69NFHJUl33nmnFi1a1K6CMG9YA1qPr8LnnQN80w7QHhGugfbBaxB2Op2yWCyeaX9/f9XW1spsNsvpdCooKMizrFOnTnI6nfXmd+rUSeXl5Y22XVhY2NL6myVGMW2y34uGRwijOZrRb/bcu6fBvJh2FGDbUy1A67t4/6buafjRB7xqaewrLS31TSFn8BqELRaLKioqPNMul0tms7nRZRUVFQoKCvLMv/zyy1VRUaHOnTs3aDcm5hIPowAAAGjXvD41Ijo6Wjt27JAkFRUVKTIy0rMsKipKhYWFqq6uVnl5uQ4ePKjIyEhFR0dr+/btkqQdO3YQegEAANDumNxut/t8K/zw1IgDBw7I7XYrMzNTO3bsUGhoqAYNGqSCggLl5+fL7Xbr0Ucf1T333KNvvvlGkyZNUkVFhX7yk58oOztbHTt2vFjHBAAAAHjlNQjjx2Hv3r1auHChbDabDh06pMmTJ8tkMikiIkIzZ86Unx+PjMb/d+rUKT333HM6evSoampqNH78ePXq1Yt+g3Oqq6vTtGnT9OWXX8pkMmn27Nnq0KEDfQYX5Ntvv1V8fLxWrlwps9lMv8F5DR8+3PP9tO7duysxMVEZGRny9/fXHXfcoccff9xn+6LnXQJeeeUVTZs2TdXV1ZKkefPm6amnntL69evldrv197//vY0rRHvz+uuvy2q1av369Vq+fLnmzJlDv8F5bdu2TZK0ceNGPfXUU3rhhRfoM7ggp06d0owZM3T55ZdL4t8onF91dbXcbrdsNptsNpvmzZunmTNnKjs7Wxs2bNDevXv1j3/8w2f7IwhfAkJDQ/XSSy95pj/99FPdeuutkk4/vu69995rq9LQTg0dOlRPPvmkJMntdsvf359+g/MaPHiw5syZI0my2+3q3LkzfQYXJCsrS0lJSerSpYsk/o3C+e3fv1/ff/+9xowZo1GjRmn37t2qqalRaGioTCaT7rjjDp/2GYLwJeCee+7xPMlDOh1sTCaTpPM/vg7G1alTJ1ksFjmdTv3+97/XU089Rb+BV2azWZMmTdKcOXN033330Wfg1ebNmxUcHOx5H4HEv1E4v8svv1yPPPKIVqxYodmzZ2vKlCn13k7s6z5DEL4EnXmv1bkeXweUlpZq1KhR+s///E/dd9999BtckKysLP3tb3/T9OnTPbdjSfQZNG7Tpk167733lJqaqs8++0yTJk3S8ePHPcvpNzhbz549df/998tkMqlnz54KCgqSw+HwLPd1nyEIX4J++tOf6oMPPpB0+vF1ffv2beOK0N588803GjNmjJ555hk98MADkug3OL8///nPWrp0qSQpMDBQJpNJN998M30G57Vu3TqtXbtWNptNN910k7KysnTnnXfSb3BOr732mubPny9J+vrrr/X999+rY8eOOnz4sNxut959912f9hmeGnGJOHLkiCZMmKCCggJ9+eWXmj59uk6dOqWwsDDNnTtX/v7+bV0i2pG5c+fqrbfeUlhYmGfe1KlTNXfuXPoNGlVZWakpU6bom2++UW1trcaOHavw8HCuNbhgqampmjVrlvz8/Og3OKeamhpNmTJFdrtdJpNJTz/9tPz8/JSZmam6ujrdcccdSktL89n+CMIAAAAwJG6NAAAAgCERhAEAAGBIBGEAAAAYEkEYAAAAhkQQBgAAgCERhAEAAGBIBGEAAAAYEkEYAAAAhvT/AHNq9LUOsJg2AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(12, 5))\n",
    "ax.axvline(x=boys_after_girl, color='b', label='observed boys after girl')\n",
    "ax.hist(boys_ppc['boys2_g'],\n",
    "        bins=23, density=True, alpha=.4, color='b', label='sampled boys after girl')\n",
    "ax.axvline(x=(n_boys - first_boys - boys_after_girl), color='g', label='observed boys after boy')\n",
    "ax.hist(boys_ppc['boys2_b'],\n",
    "        bins=23, density=True, alpha=.4, color='g', label='sampled boys after boy')\n",
    "ax.legend()\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsIAAAEtCAYAAAAC6AYTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X1UlGX+x/EPMA6rDEq25YkSE5XKWlK0tNR+Wz6krLVHN8lI1B7U3NXOkrU+pRIhmQ9ktWlLmbljCpTuqTY727qH1ZNlrnQItdSFNjWHTSzHZTAZYe7fH53mpKKjMreI1/v1T8x9Ddd8ry/3gY9X99wTYVmWJQAAAMAwkU1dAAAAANAUCMIAAAAwEkEYAAAARiIIAwAAwEgEYQAAABiJIAwAAAAjOUI9IRAIKCsrS7t27ZLT6VROTo46dOgQHC8qKlJBQYEcDocmTpyo22+/XV6vV3feeaeSkpIkSQMGDNCYMWPsWwUAAABwlkIG4fXr18vv96uwsFClpaWaN2+eli5dKkmqqqqS2+3WmjVrVFtbq/T0dPXp00eff/65hg4dqlmzZtm+AAAAAOBchLw0oqSkRP369ZMkdevWTdu3bw+OlZWVqXv37nI6nYqNjVVCQoJ27typ7du3a8eOHRo1apQeffRRHThwwL4VAAAAAOcg5I6wz+eTy+UKPo6KilJdXZ0cDod8Pp9iY2ODYzExMfL5fEpMTNQNN9ygW2+9Ve+8845ycnL0wgsvHDdvSUlJGJcBAAAAnFqPHj1OOhYyCLtcLtXU1AQfBwIBORyOBsdqamoUGxur5ORktWzZUpI0cODAk0Lwj6644oqzWwEaVF1dfdw/SBB+9Nhe9Ndeofp75ZXxkqT9+z3nq6SLCuev/eixvUzob2VlZYPHQwbhlJQUFRcXKzU1VaWlpcE3wElScnKyFi9erNraWvn9flVUVCgpKUlTp07VoEGDlJqaqo8//ljXX399g3PHx8ef43LwUx6Ph17ajB7bi/7a60z7y8/g3HD+2o8e28uE/p5zEB44cKA2bdqkkSNHyrIs5ebmavny5UpISFD//v2VkZGh9PR0WZalzMxMRUdHa8qUKZoxY4ZWr16tli1bKicnJ+wLAgAAABojZBCOjIxUdnb2ccc6deoU/DotLU1paWnHjbdv315utztMJQIAAADhxwdqAAAAwEgEYQAAABiJIAwAAAAjEYQBAABgJIIwAADARSQjI0MVFRVNWsPq1av14osvHnds7dq1WrhwYRNV1DCCMAAAAIwU8vZpAAAAOL1f/Upaty68c6amSu+9d+rxY8eOafr06fr6669VX1+vBx54QKmpqZKkF154QYcOHZLT6dT8+fMlSb///e9lWZZqa2v11FNP6brrrpPb7dbatWsVHR2t1NRUjR49WtOmTZPX65XX61XHjh118803a9iwYaqqqtKECRO0du1aLVq0SFu3blUgENDYsWM1ZMgQbd26Vbm5uWrdurWioqLUrVu3k2ouLS3VmDFj5PP5NHnyZP3yl7/Upk2btHjxYkVHRysuLk65ubl69dVX1a5dO91///06fPiwHnjgAb366qsNrqExCMIA0Izk55/993i9rRQXF/5aADStwsJCtW3bVgsXLpTP59Pw4cPVu3dvSdKgQYP0q1/9Sm+88Yb+9Kc/6ZZbblFcXJzmz5+v8vJyHTlyROXl5Vq3bp1eeOEFxcfH64EHHlDfvn0lSb1799bYsWNVXl6u7OxsDRs2TG+//baGDx+uDRs26Ouvv9bq1atVW1urtLQ09enTR0899ZReeOEFdezYUXPmzGmw5pYtWyo/P1/fffedRowYoX79+mnWrFlavXq12rVrpxUrVmjp0qVKT0/XY489pvvvv19//etfddddd6msrOykNTQWQRgAAKCRTrdza5eKigrdeuutkiSXy6VOnTpp3759kqSePXtKklJSUrRhwwZNnTpVX331lX7729/K4XBo4sSJ2r17tzwej6ZMmaLo6GgdPnxYe/bskSR17NhRktS5c2fV19dr//79WrdunV5//XUVFhZqx44dysjIkCTV1dVp//79OnjwYPD7UlJStHfv3pNq7tGjhyIiInTppZcqNjZWhw8flsvlUrt27SRJN910k/Ly8tS+fXvFxMSovLxc7777rpYsWaK4uLiT1tBYXCMMAADQDHXq1Elbt26VJPl8Pu3evVtXXXWVJGnbtm2SpK1bt6pLly765JNPdPnll+u1117TxIkTlZeXp8TERHXu3FnPPfec3G63hg8frmuuuUaSFBEREXyde+65RwsWLFDnzp3VunVrJSYmqlevXnK73VqxYoWGDBmi9u3bq127dsE36f34+if68XhVVZWOHDmiSy65RD6fTwcOHJAkbdmyRVdffbWkHz69eMmSJWrXrp3atm3b4Boaix1hAACAZigtLU2zZs3Sfffdp9raWk2aNEmXXnqpJGn9+vVasWKFYmJi9OyzzyoQCOixxx7T6tWrVVdXp9/97ne69tprdcstt2jy5MmSpOTk5ODO7E8NHjxYc+fO1dKlSyVJd9xxh7Zs2aL09HQdOXJEAwYMkMvlUnZ2tv7whz/I5XIpJiZGbdq0OWmuo0ePavTo0Tpy5Iiys7MVERGhnJwcTZ48WREREWrTpo2eeeYZSdKAAQOUnZ2tBQsWSJKuvfbak9bQWBGWZVmNnuUclJSUqEePHk3x0hcdj8ej+Pj4pi7jokaP7UV/z9y5XSPsVdxpLhKeMOGH/zbNX4Pmj/PXfvTYXhdqf7///nuNGjVKb775piIjG3cRw6lyJ5dGAAAA4ILy6aefKi0tTePGjWt0CD4dLo0AAADABSUlJUXvvvuu7a/DjjAAAACMRBAGAACAkQjCAAAAMBJBGAAAAEbizXIAAABhcC63Nzyd8ePDO9+5yMjIUFZWljp16hTyuZmZmRo5cqR69eoVPDZt2jSlpqbqtttus7PMc8aOMAAAAIzEjjAAAEAz9J///EfTp0+Xw+FQIBDQokWLdPnll2v27Nn673//qwMHDuiOO+5QZmampk2bJofDIY/HI7/fr9TUVBUXF6uyslJz5szRvn379PLLLysyMlJVVVW69957df/99wdfq7q6WjNnztShQ4ckSU8++aSuueYavfHGG3rzzTd12WWX6dtvv22wzlWrVmnZsmWqr6/X3Llz1aFDB7322mt677335HA41LNnTz3xxBMaOXKknn76aXXp0kUbNmxQcXGx7rrrLj377LNyOBxq2bKlnn/+eblcrrD1kB1hAACAZuijjz5ScnKyli9frsmTJ6u6ulqVlZXq1q2bli1bprfeeksFBQXB51955ZV67bXXlJiYqK+//lqvvPKKBg0apI8//liS9M0332jp0qUqKirS66+/flywffnll9W7d2+53W49/fTTysrK0sGDB/XnP/9ZRUVFWrJkiY4dO9ZgnSkpKVqxYoXGjRunBQsWaNeuXXr//fdVUFCggoIC7dmzR8XFxRoxYoT+8pe/SJLWrFmjESNGaP369RoyZIhWrlyp++67T//73//C2kN2hAEAAJqhe+65R6+88ooefvhhxcbGKjMzU3Fxcdq2bZs2b94sl8slv98ffH7Xrl0lSa1bt1ZiYmLw64MHD0qSunfvLqfTKUnq0qWL9u7dG/ze3bt3a/PmzXr//fclSYcPH9bevXvVuXPn4PckJyc3WGfPnj2D88+fP19ffvmlbrzxRrVo0SI4/u9//1ujRo3S8OHD9dBDD+mbb77R9ddfr6uuukovv/yyxowZo3bt2p3yNc4VO8IAAADN0D/+8Q/16NFDK1as0ODBg/Xqq69q7dq1io2N1aJFi/Tggw/q6NGjsixLkhQREXHa+b744gvV19fr+++/V3l5uTp06BAcS0xM1NixY+V2u7V48WLdfffduvrqq1VeXq6jR4+qvr5eX3zxRYPzlpWVSZK2bt2qLl26KDExUWVlZaqrq5NlWfrXv/6ljh07qlWrVurVq5fmzp2ru+++W5L0zjvvaNiwYXK73erSpYuKiorC0bogdoQBAACaoRtuuEFTp07V0qVLFQgENH36dDmdTk2ZMkWlpaVyOp3q0KGDDhw4cEbz1dXVady4cfJ6vZo4caLatm0bHHvkkUc0c+ZMFRUVyefzadKkSWrbtq3GjRunkSNHqm3btmrZsmWD83722WcaPXq0IiIilJubqyuvvFJDhgzRfffdp0AgoB49emjAgAGSpLS0NKWnpysrK0vSD7vMTz75pFq2bKnIyEhlZ2c3rmkniLB+/GfCeVZSUqIePXo0xUtfdDwej+Lj45u6jIsaPbYX/T1z53J7Jq/Xq7i4uFOOT5jww3+b5q9B88f5az96bC+Px6N9+/apoKBAzz33XJPWUlZWppUrV2r+/PlhnfdUuZMdYQAAADS5lStX6q233tLixYvP22sShAEAAAzXq1ev4z4IoymMGjVKo0aNOq+vyZvlAAAAYCSCMAAAAIxEEAYAAICRCMIAAAAwEkEYAAAARuKuEQAQRudyn18AQNNgRxgAAABGYkcYAHBe2LlbPn68fXMDuHixIwwAAAAjEYQBAABgJIIwAAAAjEQQBgAAgJEIwgAAADASQRgAAABGIggDAADASARhAAAAGIkgDAAAACMRhAEAAGCkkEE4EAho9uzZuvfee5WRkaE9e/YcN15UVKThw4crLS1NxcXFx41t2bJF//d//xfeigEAAIAwcIR6wvr16+X3+1VYWKjS0lLNmzdPS5culSRVVVXJ7XZrzZo1qq2tVXp6uvr06SOn06nKykotX75cdXV1ti8CAAAAOFshd4RLSkrUr18/SVK3bt20ffv24FhZWZm6d+8up9Op2NhYJSQkaOfOnaqtrdWcOXOUlZVlW+EAAABAY4QMwj6fTy6XK/g4KioquMvr8/kUGxsbHIuJiZHP51N2drYefPBBtWvXzoaSAQAAgMYLeWmEy+VSTU1N8HEgEJDD4WhwrKamRi1atNDWrVu1d+9evfTSSzp8+LAyMzP13HPPnTS3x+MJxxqMV11dTS9tRo/tdTH11+tt1dQlnOTo0aPyer2neUacJPt/J9vZG4/niG1zh3Ixnb8XKnpsL5P7GzIIp6SkqLi4WKmpqSotLVVSUlJwLDk5WYsXL1Ztba38fr8qKiqUnJysv/3tb8Hn9OnTp8EQLEnx8fFhWAI8Hg+9tBk9ttfF1N+4uKau4GRer1dxZ1CY3T8DO3sTH990jb+Yzt8LFT22lwn9raysbPB4yCA8cOBAbdq0SSNHjpRlWcrNzdXy5cuVkJCg/v37KyMjQ+np6bIsS5mZmYqOjg578QAAAEC4hQzCkZGRys7OPu5Yp06dgl+npaUpLS3tlN+/adOmRpQHAAAA2IMP1AAAAICRCMIAAAAwEkEYAAAARiIIAwAAwEgEYQAAABiJIAwAAAAjEYQBAABgJIIwAAAAjEQQBgAAgJEIwgAAADASQRgAAABGIggDAADASARhAAAAGIkgDAAAACMRhAEAAGAkgjAAAACMRBAGAACAkQjCAAAAMBJBGAAAAEYiCAMAAMBIBGEAAAAYiSAMAAAAIxGEAQAAYCRHUxcAAEBj5efbO//48fbOD6BpsCMMAAAAIxGEAQAAYCSCMAAAAIxEEAYAAICRCMIAAAAwEkEYAAAARiIIAwAAwEgEYQAAABiJIAwAAAAjEYQBAABgJIIwAAAAjEQQBgAAgJEIwgAAADASQRgAAABGIggDAADASI6mLgAAGpKfb9/c48fbNzcAoPlgRxgAAABGIggDAADASARhAAAAGIkgDAAAACMRhAEAAGAk7hoBAJBk7506AOBCFDIIBwIBZWVladeuXXI6ncrJyVGHDh2C40VFRSooKJDD4dDEiRN1++2368CBA3riiSd07NgxtWnTRgsWLJDL5bJ1IQAAAMDZCHlpxPr16+X3+1VYWKgpU6Zo3rx5wbGqqiq53W4VFBRo2bJlysvLk9/v1yuvvKJhw4Zp1apV6tq1q9566y1bFwEAAACcrZA7wiUlJerXr58kqVu3btq+fXtwrKysTN27d5fT6ZTT6VRCQoJ27typGTNmyLIsBQIBVVZWKj4+3r4VAAAAAOcg5I6wz+c77rKGqKgo1dXVBcdiY2ODYzExMfL5fIqIiFB9fb2GDh2qTz75RL1797ahdAAAAODchdwRdrlcqqmpCT4OBAJyOBwNjtXU1ASDcYsWLbRu3Tp99NFHmjp1qlauXHnS3B6Pp9ELgFRdXU0vbUaP7dVQf73eVra9nsdzxLa57az7XB09elRer/c0z4iTpBDPMdvpzhl+P9iPHtvL5P6GDMIpKSkqLi5WamqqSktLlZSUFBxLTk7W4sWLVVtbK7/fr4qKCiUlJSkrK0uDBw9W7969FRMTo4iIiAbn5pKJ8PB4PPTSZvTYXg31Ny7OvteLj7dvcjvrPlder1dxZ1DYmTzHVKc7Z/j9YD96bC8T+ltZWdng8ZBBeODAgdq0aZNGjhwpy7KUm5ur5cuXKyEhQf3791dGRobS09NlWZYyMzMVHR2tjIwMZWVl6aWXXlJkZKSysrLCvR4AAACgUUIG4cjISGVnZx93rFOnTsGv09LSlJaWdtK42+0OU4kAAABA+PHJcgAAADASQRgAAABGIggDAADASARhAAAAGIkgDAAAACMRhAEAAGAkgjAAAACMRBAGAACAkQjCAAAAMBJBGAAAAEYiCAMAAMBIBGEAAAAYiSAMAAAAIxGEAQAAYCSCMAAAAIxEEAYAAICRCMIAAAAwEkEYAAAARnI0dQEAcL7l5zd1BQCACwE7wgAAADASQRgAAABGIggDAADASARhAAAAGIkgDAAAACMRhAEAAGAkgjAAAACMRBAGAACAkQjCAAAAMBJBGAAAAEYiCAMAAMBIBGEAAAAYiSAMAAAAIxGEAQAAYCSCMAAAAIzkaOoCADRP+fnhm8vrbaW4uPDNBwDAmWBHGAAAAEYiCAMAAMBIBGEAAAAYiSAMAAAAIxGEAQAAYCSCMAAAAIxEEAYAAICRCMIAAAAwEkEYAAAARiIIAwAAwEgEYQAAABiJIAwAAAAjOUI9IRAIKCsrS7t27ZLT6VROTo46dOgQHC8qKlJBQYEcDocmTpyo22+/XR6PRzNmzFB9fb0sy1J2drYSExNtXQgAAABwNkLuCK9fv15+v1+FhYWaMmWK5s2bFxyrqqqS2+1WQUGBli1bpry8PPn9fj3//PMaNWqU3G63JkyYoLy8PFsXAQAAAJytkDvCJSUl6tevnySpW7du2r59e3CsrKxM3bt3l9PplNPpVEJCgnbu3KmpU6cqNjZWklRfX6/o6GibygcAAADOTcgg7PP55HK5go+joqJUV1cnh8Mhn88XDLySFBMTI5/Pp7Zt20qSvvzySz377LN66aWXGpzb4/E0tn5Iqq6uppc2o8cn83pbhW2uo0ePyuv1hm0+HC90f+MkiZ/BaXg8R045xu8H+9Fje5nc35BB2OVyqaamJvg4EAjI4XA0OFZTUxMMxps3b9ZTTz2l+fPnn/L64Pj4+EYVjx94PB56aTN6fLK4uPDN5fV6FRfOCXGcM+0vP4NTi48/dW/4/WA/emwvE/pbWVnZ4PGQ1winpKRo48aNkqTS0lIlJSUFx5KTk1VSUqLa2lpVV1eroqJCSUlJ2rx5s+bOnatXX31Vv/jFL8K0BAAAACB8Qu4IDxw4UJs2bdLIkSNlWZZyc3O1fPlyJSQkqH///srIyFB6erosy1JmZqaio6OVm5urY8eOadq0aZKkjh07Kjs72/bFAAAAAGcqZBCOjIw8KcR26tQp+HVaWprS0tKOG3/nnXfCVB4AABe3/Hz75h4/3r65gYsBH6gBAAAAIxGEAQAAYKSQl0YAAGC6012+4PW2CutdVACcP+wIAwAAwEgEYQAAABiJIAwAAAAjEYQBAABgJIIwAAAAjEQQBgAAgJEIwgAAADASQRgAAABGIggDAADASARhAAAAGIkgDAAAACMRhAEAAGAkgjAAAACMRBAGAACAkQjCAAAAMBJBGAAAAEYiCAMAAMBIBGEAAAAYiSAMAAAAIxGEAQAAYCSCMAAAAIxEEAYAAICRCMIAAAAwEkEYAAAARiIIAwAAwEgEYQAAABiJIAwAAAAjEYQBAABgJIIwAAAAjEQQBgAAgJEIwgAAADCSo6kLAAAA9sjPt3f+8ePtnR+wGzvCAAAAMBJBGAAAAEYiCAMAAMBIBGEAAAAYiSAMAAAAIxGEAQAAYCSCMAAAAIxEEAYAAICRCMIAAAAwEkEYAAAARuIjloEmZPfHnwIAgFMLuSMcCAQ0e/Zs3XvvvcrIyNCePXuOGy8qKtLw4cOVlpam4uLi48Zef/11LVy4MLwVAwAAAGEQckd4/fr18vv9KiwsVGlpqebNm6elS5dKkqqqquR2u7VmzRrV1tYqPT1dffr0USAQ0MyZM7Vt2zYNGjTI9kUAAAAAZyvkjnBJSYn69esnSerWrZu2b98eHCsrK1P37t3ldDoVGxurhIQE7dy5U7W1tRo2bJgeeeQR+yoHAAAAGiHkjrDP55PL5Qo+joqKUl1dnRwOh3w+n2JjY4NjMTEx8vl8atOmjfr27au1a9eedm6Px9OI0vGj6upqemkzu3rs9bYK+5zN0dGjR+X1epu6jItW6P7GSRI/g3Nk8vnr8Rw5L6/D3zl7mdzfkEHY5XKppqYm+DgQCMjhcDQ4VlNTc1wwDiU+Pv5sasUpeDweemkzu3ocFxf2KZslr9erOJphmzPtLz+Dc2Py+Rsff37Wzd85e5nQ38rKygaPh7w0IiUlRRs3bpQklZaWKikpKTiWnJyskpIS1dbWqrq6WhUVFceNAwAAABeqkDvCAwcO1KZNmzRy5EhZlqXc3FwtX75cCQkJ6t+/vzIyMpSeni7LspSZmano6OjzUTcAAADQKCGDcGRkpLKzs4871qlTp+DXaWlpSktLa/B7hw8f3sjyAAAAAHvwyXIAAAAwEkEYAAAARiIIAwAAwEgEYQAAABiJIAwAAAAjEYQBAABgJIIwAAAAjEQQBgAAgJEIwgAAADASQRgAAABGIggDAADASARhAAAAGIkgDAAAACMRhAEAAGAkgjAAAACMRBAGAACAkQjCAAAAMBJBGAAAAEZyNHUBwIUuP1/yelspLq6pKwEAAOFEEAYAAOckP9++ucePt29u4EdcGgEAAAAjEYQBAABgJIIwAAAAjEQQBgAAgJF4sxwAALjg/PSNeHbcuYc340FiRxgAAACGIggDAADASARhAAAAGIkgDAAAACMRhAEAAGAkgjAAAACMRBAGAACAkQjCAAAAMBJBGAAAAEYiCAMAAMBIBGEAAAAYiSAMAAAAIxGEAQAAYCSCMAAAAIzkaOoCgHDIz2/qCgAAzYmdfzfGj7dvboQXQRjnBUEVAABcaLg0AgAAAEYiCAMAAMBIBGEAAAAYiSAMAAAAI/FmOQTxhjYAABqPO1I0HyGDcCAQUFZWlnbt2iWn06mcnBx16NAhOF5UVKSCggI5HA5NnDhRt99+u7777js9/vjjOnr0qC6//HI988wzatmypa0LAQAAAM5GyEsj1q9fL7/fr8LCQk2ZMkXz5s0LjlVVVcntdqugoEDLli1TXl6e/H6/lixZoqFDh2rVqlXq2rWrCgsLbV0EAAAAcLZC7giXlJSoX79+kqRu3bpp+/btwbGysjJ1795dTqdTTqdTCQkJ2rlzp0pKSjRhwgRJ0m233aa8vDyNHTvWnhVcYJri8gKvt5Xi4s7/6wIAgPPLjpxxvnLEhXhZR8gg7PP55HK5go+joqJUV1cnh8Mhn8+n2NjY4FhMTIx8Pt9xx2NiYlRdXd3g3CUlJY2t/4LTo0dTVwAAZ2fr1qauAIAJLsTYFzIIu1wu1dTUBB8HAgE5HI4Gx2pqahQbGxs8/rOf/Uw1NTVq3br1SfP2IDECAACgCYW8RjglJUUbN26UJJWWliopKSk4lpycrJKSEtXW1qq6uloVFRVKSkpSSkqKNmzYIEnauHEjoRcAAAAXnAjLsqzTPeHHu0bs3r1blmUpNzdXGzduVEJCgvr376+ioiIVFhbKsixNmDBBd955pw4ePKipU6eqpqZGl1xyiRYtWqRWrVqdrzUBAAAAIYUMwriwHDt2TNOmTdP+/fsVGRmpp59+Wg6HQ9OmTVNERIS6dOmiOXPmKDKSz0o5V36/X9OnT9e+ffvkcrk0e/Zseb1ezZ07V1FRUerbt68mTZrU1GU2S5999pkWLlwot9utPXv2NHje/vGPf9Q///lPORwOzZgxQ8nJyU1ddrPx0/7+KDc3Vx07dtR9990nqeFbXuLM/LS/X3zxhZ5++mlFRUXJ6XTq2Wef1c9//nP62wg/7W95eblmzZoly7J09dVXKycnRw6Hg/42UkO/I959912tXLkyeIcv43psoVn5+9//bj366KOWZVnWhx9+aE2aNMmaMGGCtXnzZsuyLGvWrFnWBx980JQlNntut9t68sknLcuyrIqKCuvBBx+07r77bmvPnj1WIBCwHn74YWvHjh1NXGXzk5+fbw0dOtQaMWKEZVlWg+ft9u3brYyMDCsQCFj79++3hg8f3pQlNysn9vfbb7+1HnroIat///7WqlWrLMuyrAMHDlhDhw61amtrrf/973/BrxHaif29//77rc8//9yyLMtavXq1lZubS38b4cT+Tpw40dqyZYtlWZY1depU64MPPqC/jXRijy3Lsnbs2GGNHj06eMzEHrNt2Mx07NhR9fX1CgQC8vl8cjgc2rFjh26++WZJP9yu7qOPPmriKpu38vJy3XbbbZKkxMREbdu2TX6/XwkJCYqIiFDfvn3p8TlISEjQiy++GHzc0HlbUlKivn37KiIiQvHx8aqvr9d3333XVCU3Kyf2t6amRpMnT9avf/3r4LGf3vIyNjY2eMtLhHZif/Py8nTddddJkurr6xUdHU1/G+HE/r744ou66aab5Pf7VVVVJZfLRX8b6cQeHzp0SHl5eZoxY0bwmIk9Jgg3M61atdL+/fs1ZMgQzZo1SxkZGbIsSxEREZJOf7s6nJnrrrtOxcXFsixLpaWlqq6uPu4ad3pJQ3wEAAACo0lEQVR8bu68887gHWckNXjenni7Rnp95k7sb/v27XXjjTce95xT3fISoZ3Y38svv1yS9Omnn2rlypUaO3Ys/W2EE/sbFRWl/fv3a+jQoTp06JCuvfZa+ttIP+1xfX29Zs6cqenTpysmJib4HBN7HPL2abiwvP766+rbt6+mTJmiyspKjRkzRseOHQuOn+p2dThzv/nNb1RRUaH09HSlpKTo2muv1ffffx8cp8fh8dPr2H/s6aluyYjwoL/htW7dOi1dulT5+flq27Yt/Q2zK6+8Uh988IHefPNNzZs3T4MGDaK/YbJjxw7t2bNHWVlZqq2tVXl5uebOnavevXsb12N2hJuZ1q1bB0/KNm3aqK6uTl27dtUnn3wi6Yfb1fXs2bMpS2z2tm3bpltuuUWrV6/W4MGDdfXVV6tFixbau3evLMvShx9+SI/DoKHzNiUlRR9++KECgYA8Ho8CgYDatm3bxJVePE51y0ucvbffflsrV66U2+1W+/btJdHfcHrkkUf01VdfSfphVzIyMpL+hlFycrLee+89ud1u5eXlqXPnzpo5c6aRPWZHuJkZO3asZsyYofT0dB07dkyZmZm64YYbNGvWLOXl5SkxMVF33nlnU5fZrHXo0EHPP/+8Xn75ZcXGxmru3LmqrKzU448/rvr6evXt2/ek/+WMszd16tSTztuoqCj17NlT9957rwKBgGbPnt3UZV5ULrvsMmVkZCg9PV2WZSkzM1PR0dFNXVazU19fr7lz5+qKK67Q5MmTJUk33XSTHn30UfobJuPHj9e0adPUokULtWzZUjk5OZy/54GJPeb2aQAAADASl0YAAADASARhAAAAGIkgDAAAACMRhAEAAGAkgjAAAACMRBAGAACAkQjCAAAAMBJBGAAAAEb6fz4du0aP4JcwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(12, 5))\n",
    "ax.axvline(x=n_boys, color='b', label='observed boys')\n",
    "ax.hist(boys_ppc.apply(np.sum, axis=1),\n",
    "        bins=26, density=True, alpha=.4, color='b',label='sampled boys')\n",
    "ax.legend()\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}