062 Stan simulation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 24,
   "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>father</th>\n",
       "      <th>son</th>\n",
       "      <th>count</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>farm</td>\n",
       "      <td>farm</td>\n",
       "      <td>703</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>farm</td>\n",
       "      <td>unskilled</td>\n",
       "      <td>1478</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>farm</td>\n",
       "      <td>skilled</td>\n",
       "      <td>1430</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>farm</td>\n",
       "      <td>professional</td>\n",
       "      <td>1109</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>unskilled</td>\n",
       "      <td>farm</td>\n",
       "      <td>58</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>unskilled</td>\n",
       "      <td>unskilled</td>\n",
       "      <td>1756</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>unskilled</td>\n",
       "      <td>skilled</td>\n",
       "      <td>1630</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>unskilled</td>\n",
       "      <td>professional</td>\n",
       "      <td>1568</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>skilled</td>\n",
       "      <td>farm</td>\n",
       "      <td>63</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>skilled</td>\n",
       "      <td>unskilled</td>\n",
       "      <td>1453</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>skilled</td>\n",
       "      <td>skilled</td>\n",
       "      <td>2068</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>skilled</td>\n",
       "      <td>professional</td>\n",
       "      <td>2483</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>professional</td>\n",
       "      <td>farm</td>\n",
       "      <td>61</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>professional</td>\n",
       "      <td>unskilled</td>\n",
       "      <td>749</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>professional</td>\n",
       "      <td>skilled</td>\n",
       "      <td>1183</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>professional</td>\n",
       "      <td>professional</td>\n",
       "      <td>3315</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          father           son  count\n",
       "0           farm          farm    703\n",
       "1           farm     unskilled   1478\n",
       "2           farm       skilled   1430\n",
       "3           farm  professional   1109\n",
       "4      unskilled          farm     58\n",
       "5      unskilled     unskilled   1756\n",
       "6      unskilled       skilled   1630\n",
       "7      unskilled  professional   1568\n",
       "8        skilled          farm     63\n",
       "9        skilled     unskilled   1453\n",
       "10       skilled       skilled   2068\n",
       "11       skilled  professional   2483\n",
       "12  professional          farm     61\n",
       "13  professional     unskilled    749\n",
       "14  professional       skilled   1183\n",
       "15  professional  professional   3315"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy import stats as sts\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "data = pd.read_csv(\"socialmobility.csv\")\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "95% CI that the son is skilled while father is unskilled\n",
      " [0.31220254 0.33826378]\n",
      "\n",
      "95% CI that the son is professional while father is unskilled\n",
      " [0.30010122 0.32565323]\n"
     ]
    }
   ],
   "source": [
    "subset = data[data.father == \"unskilled\"]\n",
    "\"\"\"\n",
    "Setting a uniform distribution over \n",
    "the probability vector parameter of the multinomial.\n",
    "\"\"\"\n",
    "alpha_priors = np.repeat(1, len(subset))\n",
    "\"\"\"\n",
    "based on the conjugate prior table, we can compute the posterior \n",
    "hyperparameters of the Dirichlet posterior distribution by\n",
    "\"\"\"\n",
    "alpha_posteriors = [alpha_priors[i] + subset['count'].iloc[i] for i in range(len(subset))]\n",
    "\n",
    "posterior = sts.dirichlet(alpha_posteriors)\n",
    "\n",
    "samples = posterior.rvs(size=10**5)\n",
    "\n",
    "\"\"\"\n",
    "What is the posterior 95% probability interval over the probability that a son will\n",
    "become a skilled laborer if his father was an unskilled laborer?\n",
    "\"\"\"\n",
    "print(\"95% CI that the son is skilled while father is unskilled\\n\",\n",
    "      np.percentile(samples[:,2], [2.5, 97.5]))\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "What is the posterior 95% probability interval over the probability that a son will\n",
    "become a professional laborer if his father was an unskilled laborer?\n",
    "\"\"\"\n",
    "print(\"\\n95% CI that the son is professional while father is unskilled\\n\",\n",
    "      np.percentile(samples[:,3], [2.5, 97.5]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          father           son  count\n",
      "3           farm  professional   1109\n",
      "7      unskilled  professional   1568\n",
      "11       skilled  professional   2483\n",
      "15  professional  professional   3315\n",
      "\n",
      "95% CI that father works in a farm if son is a professional\n",
      " [0.12379266 0.13817723]\n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "What is the posterior 95% probability interval over the probability that a father \n",
    "works on a farm if his son works as a professional?\n",
    "\"\"\"\n",
    "\n",
    "subset2 = data[data.son == \"professional\"]\n",
    "print(subset2)\n",
    "alpha_priors = np.repeat(1, len(subset2))\n",
    "\n",
    "alpha_posteriors = [alpha_priors[i] + subset2['count'].iloc[i] for i in range(len(subset2))]\n",
    "\n",
    "posterior2 = sts.dirichlet(alpha_posteriors)\n",
    "samples2 = posterior2.rvs(size=10**5)\n",
    "\n",
    "print(\"\\n95% CI that father works in a farm if son is a professional\\n\",\n",
    "      np.percentile(samples2[:,0], [2.5, 97.5]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generating 2 samples:\n",
      "[[19 15 18 18 15 15]\n",
      " [16 16 15 18 16 19]]\n",
      "Value of probability mass function at samples:\n",
      "[1.28115926e-05 1.53138568e-05]\n"
     ]
    }
   ],
   "source": [
    "mult = sts.multinomial(100, [1/6, 1/6, 1/6, 1/6, 1/6, 1/6])\n",
    "num_samples = 2\n",
    "print('Generating', num_samples, 'samples:')\n",
    "samples = mult.rvs(size=num_samples)\n",
    "\n",
    "\n",
    "print(samples)\n",
    "print('Value of probability mass function at samples:')\n",
    "print(mult.pmf(samples))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generating 5 samples:\n",
      "[[0.16894931 0.40277231 0.20100879 0.22726959]\n",
      " [0.07946066 0.11558479 0.27263706 0.53231749]\n",
      " [0.14153622 0.15326053 0.59219034 0.11301291]\n",
      " [0.22383789 0.08118641 0.33564232 0.35933338]\n",
      " [0.00160571 0.41745923 0.39209485 0.18884022]]\n",
      "Sum of each sample (should be 1):\n",
      "[1. 1. 1. 1. 1.]\n",
      "Value of probability density function at each sample:\n",
      "[ 5.77690819 39.18902484  2.34595107 12.83249877 13.06957865]\n"
     ]
    }
   ],
   "source": [
    "dirichlet = sts.dirichlet([1, 2, 3, 4])\n",
    "num_samples = 5\n",
    "print('Generating', num_samples, 'samples:')\n",
    "samples = dirichlet.rvs(size=num_samples)\n",
    "print(samples)\n",
    "print('Sum of each sample (should be 1):')\n",
    "print(samples.sum(axis=1))\n",
    "print('Value of probability density function at each sample:')\n",
    "print(dirichlet.pdf(samples.transpose()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using the Stan simulation [example from class]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pystan\n",
    "# Code the type of work with integers\n",
    "farm = 0\n",
    "unskilled = 1\n",
    "skilled = 2\n",
    "professional = 3\n",
    "\n",
    "# Matrix of counts such that, for example, counts[farm][skilled] is the\n",
    "# number of times a father does farm work and a son skilled work.\n",
    "counts = [\n",
    "    [703, 1478, 1430, 1109],\n",
    "    [58, 1756, 1630, 1568],\n",
    "    [63, 1453, 2068, 2483],\n",
    "    [61, 749, 1183, 3315]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_4288838636dda9ea0168308ec88e8b29 NOW.\n"
     ]
    }
   ],
   "source": [
    "# This is a partial Stan model. Look up how to define the multinomial\n",
    "# and dirichlet distributions in your likelihood and prior.\n",
    "\n",
    "stan_code = \"\"\"\n",
    "\n",
    "data {\n",
    "    int<lower=1> C;    // number of categories in the multinomial\n",
    "    vector[C] alpha;   // prior hyperparameters for Dirichlet\n",
    "    int counts[C];     // the data set: counts for each category\n",
    "}\n",
    "\n",
    "parameters {\n",
    "    simplex[C] theta;  // the probability vector over categories\n",
    "}\n",
    "\n",
    "model {\n",
    "\n",
    "    theta ~ dirichlet(alpha); // prior\n",
    "    counts ~ multinomial(theta); // likelihood\n",
    "}\n",
    "\"\"\"\n",
    "\n",
    "stan_model = pystan.StanModel(model_code=stan_code)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate results\n",
    "\n",
    "# NOTE: The easiest option is to generate one set of posterior samples\n",
    "# for each row in the counts matrix.\n",
    "\n",
    "stan_data = {\n",
    "    'C': 4,\n",
    "    'alpha': [1, 1, 1, 1],\n",
    "    'counts': counts[unskilled],\n",
    "}\n",
    "\n",
    "results = stan_model.sampling(data=stan_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Inference for Stan model: anon_model_4288838636dda9ea0168308ec88e8b29.\n",
       "4 chains, each with iter=2000; warmup=1000; thin=1; \n",
       "post-warmup draws per chain=1000, total post-warmup draws=4000.\n",
       "\n",
       "           mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat\n",
       "theta[1]   0.01  2.3e-5 1.5e-3 9.0e-3   0.01   0.01   0.01   0.01   4198    1.0\n",
       "theta[2]   0.35  1.1e-4 6.8e-3   0.34   0.35   0.35   0.35   0.36   3656    1.0\n",
       "theta[3]   0.33  1.1e-4 6.6e-3   0.31   0.32   0.33   0.33   0.34   3786    1.0\n",
       "theta[4]   0.31  1.0e-4 6.6e-3    0.3   0.31   0.31   0.32   0.33   4006    1.0\n",
       "lp__      -5762    0.03   1.23  -5765  -5763  -5762  -5761  -5761   2086    1.0\n",
       "\n",
       "Samples were drawn using NUTS at Thu Oct 17 12:47:43 2019.\n",
       "For each parameter, n_eff is a crude measure of effective sample size,\n",
       "and Rhat is the potential scale reduction factor on split chains (at \n",
       "convergence, Rhat=1)."
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Inference for Stan model: anon_model_4288838636dda9ea0168308ec88e8b29.\n",
      "4 chains, each with iter=2000; warmup=1000; thin=1; \n",
      "post-warmup draws per chain=1000, total post-warmup draws=4000.\n",
      "\n",
      "           mean se_mean     sd   2.5%    50%  97.5%  n_eff   Rhat\n",
      "theta[1]   0.01  2.3e-5 1.5e-3 9.0e-3   0.01   0.01   4198    1.0\n",
      "theta[2]   0.35  1.1e-4 6.8e-3   0.34   0.35   0.36   3656    1.0\n",
      "theta[3]   0.33  1.1e-4 6.6e-3   0.31   0.33   0.34   3786    1.0\n",
      "theta[4]   0.31  1.0e-4 6.6e-3    0.3   0.31   0.33   4006    1.0\n",
      "\n",
      "Samples were drawn using NUTS at Thu Oct 17 12:47:43 2019.\n",
      "For each parameter, n_eff is a crude measure of effective sample size,\n",
      "and Rhat is the potential scale reduction factor on split chains (at \n",
      "convergence, Rhat=1).\n"
     ]
    }
   ],
   "source": [
    "print(results.stansummary(pars=['theta'], probs=[0.025, 0.5, 0.975]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAswAAAF1CAYAAAD8/Lw6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xt8XWWd9/3vL2maHtKDJZkeQktxWp+2WhtpUmQUFVGs6EB5hhkRppSOTIGK46PO7cDoDIrO4DjPrd4DFOkwVuiIyDADd/VGhBFQdIokxUDtASlMpU0PptQe0kOaw+/+41pr75XdnZWdNslOm8+7r/3K3mtd61rXXtlJv/nta61t7i4AAAAA+ZUUewAAAADAYEZgBgAAAFIQmAEAAIAUBGYAAAAgBYEZAAAASEFgBgAAAFIQmHFKMrNvmtnf9FFf08ysxcxKo8dPm9m1fdF31N8PzWxJX/V3ujKzb5vZl4s9jqHAzG4ws93R6/6MPOvnmdkvzeygmf1vM3tDSl/vMLOXo74W9e/IM/v87kDtK7HPrWY2PWX9l81sj5ntOon+33ei4yug/xYze2N/9d8bab9jk78vzewaM/tZYp2b2YwT3GdmWzP7mpldfyL9YOgiMGPQif7jOBL9Z73PzP7LzK43s8zr1d2vd/cvFdhX6n9C7v6au1e4e0cfjP0LZvavOf1/0N3vPdm+AUkys+nRf/7DTnD7Mklfk3RR9Lp/PWd9iaT/kHSXpAmSDkXtu3OrpDuivh45kTH1hpm9VdI8Sf+7v/dVKDObKukzkua4+6QC2g/4H4fR9+fVgdzniRig35f/KOlzZja8n/eD0wiBGYPVH7r7GElnSfqKpL+S9C99vZMTDR1DBcfntDRR0ghJG7pZf6GkMndf6e5tkv5a0kfMbGw37c9K6SvVCb6+rpP0HR9cn7p1lqTX3f23xR4IP7M9c/edkjZLuqTYY8Gpg8CMQc3d97v7GkkfkbTEzN4ida3QmFmlmf0gqkbvNbNnzKzEzFZLmibp+9HbkZ9NVOc+ZmavSXqym4rd75vZc2a2P3pLekK0r/eY2fbkGOMqtpktVDZctJjZC9H6zNuP0bg+b2a/MbPfmtl9ZjYuWhePY4mZvRa9vfu57o6NmY2Ltm+O+vt81H95dCzekmhbFVXtfy96/GEza0xU8N+a83z+ysxelHQo33/AZjbLzJ6IjvdLZvYniXUfsvB2/gEz22ZmX8jZ9p3RPvdF669JrH6Dmf2f6N2FX5jZ76c8/7z9dHdconXXmNnPzezr0XavmtkfRMu3Rd+TJYl9fNvC9J8nojH9xMzOSqz/AzOrj14n9Wb2B4l1T5vZl6L9HTSzx82sMrH+7Ynxv2Bm7ylw259GX/dFr7Pz8hybcjP7hpntiG7fiJa9SdJLie2fzHNoz5X08/iBu2+V1CJpbp79vCLpjcr+jJWb2RQzWxO9NraY2Z8n2n/BzB4ys381swOSrsnTZ+oxl/RBST9JtJ8Rtdkf/cx8L7HuhL8/hbLwDtYTkqZEx+Db0fJ/M7Nd0b5/amZvjpYvk3SVpM9G7b+f6K7GzF6MtvmemY1I7Odkf2aTUxIuNrON0fNuMrO/7Oa5dXnHzHJ+V6YdQzMbEX2fX4/GXG9mE/PsY3L0nP8y0WePU+Ki19r/b+F35e7oNTMysf5/mNnO6PX/Z3m6eFrSh3raD5Dh7ty4DaqbpK2S3pdn+WuSbojuf1vSl6P7t0n6pqSy6Ha+JMvXl6TpklzSfZJGSxqZWDYsavO0pCZJb4na/Lukf43WvUfS9u7GK+kLcdvE+qclXRvd/zNJWxRCRoXCW9+rc8b2z9G45klqlTS7m+N0n8Lb0mOibX8t6WPRum9J+rtE249Leiy6f46k3yoEo1JJS6LnUJ54Po2SpkoamWe/oyVtk7RU0rCovz2S3pw4RnMV/iB/q6TdkhZF66ZJOijpo9H36gxJNYnv6V5JC6J+vyPpgW6ee1o/acflGknt0dhLJX1Z4XV1p6RySRdF/VYkxnRQ0rui9f9L0s+idRMk/U7S4mi8H40en5H4vr8i6U3R9/NpSV+J1lVLel3SxdFxen/0uKqAbacr8Xrt5vjcKulZSb8nqUrSf0n6UiHbS7pD0j/kLGuQ9EeF/LwqhNkVClXsGknNki5M/Hy0SVoUPe98r6+0Yz46GntVov13JX0u6m+EpHee7Pcn5XlO72bde3T874U/U3gNlkv6hqTGnOf45Tz9PydpSjT2TZKu74uf2aiNS5oR3d8p6fzo/hskndPNNl9Q4vdZ7msn7RgqvBPwfUmjojHPlzQ2+TtR2Z/PZd38vrwm/t7neQ7fkLQmOlZjon3dFq1bqPB7J/4dfn9y26jN/yvp+e6+39y45d6oMONUskPhl2OuNkmTJZ3l7m3u/oy79/R27Rfc/ZC7H+lm/Wp3/5W7H5L0N5L+xKKTAk/SVZK+5u6vunuLpJslXZFTEfqiux9x9xckvaAQnLuIxvIRSTe7+0EPVcD/qRAOpPAfxEcTm1wZLZOkP5d0t7v/wt07PMwXbJX09kT7f3L3bd0cnw9L2uruq9y93d2fV/ij4nJJcven3X29u3e6+4sKgebdief/n+7+3eh79bq7Nyb6/g93f87d2xUCc03KcTyunwKOiyT9dzT2DknfUwgZt7p7q7s/LumYpOSJRf/H3X/q7q0Kwew8C3NWPyTpZXdfHR2H7yq8zfuHiW1Xufuvo+P4YOL5/KmkR9390eg4PaEQSi8uYNtCXBU9p9+6e7OkL+Ycg558MqoK7jOzfYXuOzou75T0V+5+NPre3pOz77Xu/kj0vLv7+evumI+P1h9MtG1TmBIxJdpnfJLYyXx/Tpq7fyt6DbYqBM95Fr2blOKf3H2Hu+9VCIDxeE72ZzZXm6Q5ZjbW3X8X/QyfqO6OYZvCH7IzojGvc/cDie3mKITjW9x9ZW92aGamcEw+5e573f2gpL+XdEXU5E+iccW/w7+Qp5uDyr6egB4RmHEqqVaoQOb6R4Wq7eMW3mK/qYC+tvVi/W8Uqpi9frs2jylRf8m+hynMK40lz7I/rFCJzlUpaXievqqj+09KGmlm50ZvZ9dIejhad5akz+QEoqnR2GJpx+csSefmbH+VpEmSFO3zKQtTIvZLul7ZYzdVoSLVnUKee1o/PR0XKVSeYkckyd1zlyX3mzkW0R85exWOVe73Mt++uns+Z0n645xj+E6FP/x62rYQ+V5nU7ppm2unpBXuPj6+SVofLS9kv3GASe47eUx6+tnr0ibnmO+LFo9JtP2sJJP0nJltSLz9fjLfn5NiZqVm9hUzeyWaerI1WtXT75C018vJ/Mzm+iOFP85+E01nOW5aTy90N+bVkn4k6YFoWsRXLZxwGrtK4Z28h05gn1UKlet1iePxWLRcCscl93d4rjHKvp6AHhGYcUowszqF/+h+lrsuquJ8xt3fqFA9+rSZXRiv7qbLnirQUxP3pylUS/YoXDFgVGJcpcr+ki6k3x0K//kl+25X1xBXiD3KVtaSfTVJkrt3KlR7PqpQXf5BIsRsU5iuMT5xGxVV4Ap5Htsk/SRn+wp3vyFaf7/CW6VT3X2cwnQZS2zb7bzkXuiun9TjcoIyrwUzq1B4l2OHjv9e9mZf2xTexUgew9Hu/pUCti3kZLd8r7MdBWwnSb9QosIehZyJkn5V4H4nmFky0OYek0LGn/eYR9XCeApA6Mx9l7v/ubtPUZgGsCKaq3sy35+TdaWkSyW9T9I4hakHUvbnoLcnLJ7sz2wX7l7v7pcqTNl5ROF3RT5dft8p+qO4wH20ufsX3X2OpD9QeGfq6kSTLyj8vN5/Au/e7VH4w/bNieMxzt3jsL5Tx/8OzzVb4R08oCAEZgxqZjbWzD4s6QGFuXTr87T5sIUTf0zSAUkd0U0KQfRErj36p2Y2x8xGKcwHfSh6C//XkkZYOLGtTNLnFeYoxnZLmm6JS+Dl+K6kT5nZ2VEQ+HtJ34umIBQsGsuDkv7OzMZEVeRPS0pe0u5+hekJVyk7HUMKc6SvjyrBZmajo+eTDDlpfiDpTWa22MzKoludmc2O1o9RqDIeNbMFCuEh9h1J7zOzPzGzYWZ2hpmdyNvgefsp8Lj01sUWTjAcLulLkn7h7tskPapwHK6MxvARhbeZf1BAn/8q6Q/N7ANRNXKEhRNKzyxg22ZJnUp/XX9X0uctnOxZKelvVfgxeFLSTDP7oJmVK0xJ+lHO2+l5RcflvyTdFj2nt0r6mML3qze6O+ZSOO7xFB+Z2R8njtvvFIJjh07u+3OyxihMmXhdIXD+fc763v5eOtmf2QwzG25mV5nZOA9XQYl/Z+bTKOldFq5VP05hClmh+7nAzOZGYfiAwh+yyf20SfpjhTnGq1N+Zx4nKgj8s6SvW/ZE5moz+0DU5EFJ1yR+h9+Sp5t3S/phofsECMwYrL5vZgcVKiufU7gO7NJu2s6U9J8KZ/KvVXg7+elo3W0KwWGfdXMmeDdWK5yYs0vhRKK/kMJVOyQtV5iX2aRQgUleNePfoq+vm1m+eYHfivr+qaT/lnRU0id6Ma6kT0T7f1Wh8n5/1L+isf4iWj9Fif8Y3L1BYf7fHQoBY4vyXK2gO1Gl+iKF+YI7FI7RPyj7h8NySbdG37+/VaJ65e6vKbwV/BmFt9kblWeOdgFjSOsn9bicgPsV/sPdq3Di0lXRGF5XqJp9RiEYfVbSh919TwHj36ZQgfxrhQC8TdL/UAG/k939sKS/k/Tz6HX99jzNvqwwJ/pFhekUz0fLehSFkY8q/OzsVZjO8+lCto18VKGiukNhGtAt0Rzt3sh7zCMrJV0V/YEsSXWSfmFmLQrvbHzS3f/7ZL4/feA+hWkATZI2KpyAmfQvCnOI95lZj9euPtmf2TwWS9oaTRe5XmFOfb79PqEwz/9FSevUuz82JilMtzigcALjT5TzR5u7H1M4+e73JH2rN6FZ4VKjWyQ9Gz2P/5T0/0T9/lDhpMAnozZdrgZjZpMV/njq9+uG4/QRX0kAAJDDwiXCtrv754s9lqGikGNuZvdLetAH4INSEvvcKuk9Hk4kxSnMzP6npFfcfUWxx4JTBxc4BwCcUtz9yp5bAfm5+2eKPQacepiSAQBAz74hrqoADFlMyQAAAABSUGEGAAAAUhCYAQAAgBSD7qS/yspKnz59erGHAQAAgNPcunXr9rh7VU/tBl1gnj59uhoaGoo9DAAAAJzmzCzfR6cfhykZAAAAQAoCMwAAAJCCwAwAAACkIDADAAAAKQjMAAAAQAoCMwAAAJCCwAwAAACkIDADAAAAKQjMAAAAQAoCMwAAAJCCwAwAAACkIDADGFyeuq3YIwAAoAsCM4CswRBWL7i52CMAAKALAjOArBMJq30ZsgdDYAcAIAeBGUB+qy7O3s8XZJ+6LdxOtiKc7JvqMgBgEDJ3L/YYuqitrfWGhoZiDwPAyeiLIA0AQD8zs3XuXttTOyrMALoXV5F7WpfbJjcsM9UCAHAKIzAD6Dr9Inn/gpvDLV/gTVvXk+6meAAAMAgRmAFISx/NBuWlj4bwmgzOSasuDuu//pbwOA7OyXU9ueDmbP8EZQDAIEdgBhDC6/Tzs48vuDn7ODfQTj8/rP/Ur45fH6/L3WbrM13bPnVbaBs/TobntLAOAEAREJiBoe6p20JVOQ61sWQVOGnrM9lQmwy8qy7uuk1cgU62iSvRyX1tfSbsP7nfZJgGAKDICMzAUHfBzSHc5laUv/4Wadf6cD8OyU/dJu17TfrZ18Lyxu9kw+/SR6XbpmbDd81VoX0yEMd971ovPbsiO4ZVF4e+4v3G4wIAYBAYVuwBABgEju6XfvrVcL/xO9kw+/bl4evSR0OQHT8ttM2177UQlm/elq1YJ6d5xFXnOJzH/V5ws7QqEahvmyqNGBdC9qqcyjMAAEVCYAYgTZobwu3WZ0IolrKhd3t9WH50fwjGI8ZJ4+eGNjVXhXUHmqR3fTaE5bjqvGt9aH90f9jmU78K65OBOw7ScbgeMS70+eyKEL4BABgEmJIBIPjpV0No3fdauH3xDeHrsBHha6xld/gaT73Y95o0tjpUi+OpGlIIv3E4btkdwvJPviK1H81O49heH25SCNotu8O6SXOzUzMAACgyAjMw1N02NYTW4RXZ0Dp+WgjBR/dLx1rC1/aj2QD82toQgH/z89D2QFMI2MNGhO2nn5+tLB9rkTqOhb6tRKqYGEL2T78qdbZLZ9aFavT4aaHdb34e1sfBHACAIiMwA0PdiHEhxLYeCMFXylZ9Ww9I084LXzvbQ+W341gI08kqccmwbF/DRoRKshRCsHdKpcOl/YkpFnHwLhmWPbFwe30I1OOmhsdn1vXfcwYAoBcIzABCUC4fG6rMne0hFB/dL531jmx4LhmWvX90f6gElw4P7aUQog80hXBtJeFrbpvhFSE4H2sJ94eNCO1aD4R9vuuzYf3R/aGKDQDAIEBgBoa6o/uzFeLWAyH4xnatD6HWol8Vne0hWLcekN59Uwi53hna7N8WQnD52NDfWe8I28Vt4hCdNGJctu/S4dnKdO44AAAoooICs5ktNLOXzGyLmd2UZ/31ZrbezBrN7GdmNidaPt3MjkTLG83sm339BACcpHiaRVwZjqdOHGsJ85ZbD4QA3HEsLG8/Gr7G4TZWOjy7TgrVZe/Mrovvl4/N3j/QFO6PmxrCeHJKxn6ukgEAGBx6vKycmZVKulPS+yVtl1RvZmvcfWOi2f3u/s2o/SWSviZpYbTuFXev6dthA+gz+14LgTauBJePDcuT1eA4LHtn9n68Tenw0LZ8bKgY5wbduMoci/tN9p/cpmV3CM3xfGoAAIqskArzAklb3P1Vdz8m6QFJlyYbuHvyfdbRkrzvhgig38XB10qyleV4eax8bLZNMgTHX1sPhOBrOb9WvDNMz8gV95O7rONY6KeEy8QDAAaHQgJztaRkyWh7tKwLM/u4mb0i6auS/iKx6mwz+6WZ/cTMzs+3AzNbZmYNZtbQ3Nzci+EDOGlxdbfjWJiDPLa6a1CWQpBtPxqmTSSDbOnw4/tLzj22khC0f/PzbOVayp5gmKxoS12nbQAAMEgUEpgtz7LjKsjufqe7/76kv5L0+WjxTknT3P1tkj4t6X4zO+5/Qndf6e617l5bVVVV+OgBnLw4nJaPDaE4DtDvvim73DtDWI6nZMRBNzdYS12vn+ydx0/BiEN2vqkZUpiOkbsMAIAiKiQwb5c0NfH4TEk7Uto/IGmRJLl7q7u/Ht1fJ+kVSW86saEC6BfHWrIhNZ6aUT42fGpfXFmWstdTLh/bNdTGAXjc1ONDdOnwrlVoKwlX1Gg/mp2OkZyWEZ90OG4q12FGUa1oXNHv25zIPgAURyGBuV7STDM728yGS7pC0ppkAzObmXj4IUkvR8uropMGZWZvlDRT0qt9MXAAfWRsdbYqnKz+JivKUvbqFcdaQqgtH9s1IO/fFraLQ3VcuU6eJDi8ImwfO+sdXadlxFfM4FP+MIDyBdflNct7vW2h25xoewDF02Ngdvd2STdK+pGkTZIedPcNZnZrdEUMSbrRzDaYWaPC1Isl0fJ3SXrRzF6Q9JCk6919b58/CwAnbvy08El/Z70jVH87jmXDcfnY7Dzk/dtCOyl7ZYz2o/kryMkTB+M+/qY5bBN/YMnwiux2ycvRjZ8W1gMDJDe4Flr5XdG4IrNtvm3iZWn9rWhcoaWPLS10qACKpKDrMLv7o+7+Jnf/fXf/u2jZ37r7muj+J939ze5e4+4XuPuGaPm/R8vnufs57v79/nsqAE7I9Ohc3NfWhqB71jtCaJ12XtdqcOnwsDy+JrOVhABdMTFbUY5P+ItP3us4Fi41VzJMum1qmBc9YlxYN2Jc+GCUcVO79rG9PlwbOv7IbKAfrWhccVygTQboOMzmtonDclpQjvtJqyQvr1muVQtXHbctgMGFT/oDhrqtz3S9DvKu9eHazLvWh6CbrATve01656dDyB1bHT4lsGV3qBCPGCd96leh3btvyobguGI8YlzY14Em6e3Lw7ZvzxMk3vnpsJ9864ABFofZ+l31mWXJynLu19yqc76gnW9ZjGkawOBEYAaGuqWPSl/YHwLwLb8Lwfbo/hCC/ya6zOOkueEkvPHTwsmANVdlA+87Px1ukrTq4vC18TthfVwpvnlb2FYK+9n6TNh26zPZcdwchfZnV4T+L7h5YJ4/hqzuwmzuNIkVjStUN6n7k1CTIbh+V72WPra0S5W5u6px3C4OyUsfW8r0DGCQIjADkJ66LYTU+GsckG+L5jJPPz87dSOeXzxiXAjGW58JXz/1q1AZHjEuhOO3L88+/lJV2H7po6H/6eeH7eLHsUlzs5Xop24b2GOAIWd5zXI9suWRvNXi5DSJZHU52S536kUyWOdWnHOryEsfW9plH0sfW6q6SXWpwRw4FTXffkexh9AnzH1wfShfbW2tNzQ0FHsYwNCx6uIQYBu/Ex7H0yok6etvCeE3DrhxtTgOs1ufCesuuDm0rbkqVIhv3pbtVwp9x8H42RVdp2/EVel9r2Wr1nGYBvpJMsgmp0Fc9NBFWjRjUerJfLHVG1dr8ZzFmbbxtlLXynL9rnrVTarrNkQnA3kyRAPof2a2zt1re2rHZ88CQ93SR0NoHT8thNY4wE4/Pxtqn7ottEsG5X2vdV0f37/g5vB432tdp1XEoToO08n9JyWDNtBP4sD6yJZHMsviAHz3C3drec1ynXf/eVo8Z3GX7R7Z8ogev/xxLX1sqdZeuVZLH1uquffOVe3E2kxYfmTLI5l2Fz10kR6//HGdd/95mf1K2SAe73/RjEW6Z/09/f68gd5qvv0OVX3ixmIPo+ioMAMIATcZbpOPV12cDdVxaI7XxfeT0yeSVeek26Zmq9O56/L1DfSzix66SHuO7FHlyEotmrFIqzeulqRM1fic1eeocmSldh/arevmXZepBNdNqtM96+9ReWl5JlDH7SXp2rnXannN8kz/8WMpTL1oammSpEwl+6KHLlJ1RbU2792sw22H9cKSFwb6UGAIKkYQHozhu9AKM4EZQNdQnAy73YXj5LJYPI2ip9Db3dzkfH0TntGPzrv/PI0ZPkaStOfIHklS5chKSdLOQztVVlKm8tJyzZowS1KYLjH33rkqUYkmjg7XJK+uqFZTS1MmVK/euDrT/oXmEHzLS8sz+zl47GAmZNfvqtfmvZs1ZvgYHTx2UJI0ZvgYPX754wPx9AGo8MDMSX8AstMilj56/DSKfOLQe8HN2dv087uG4eS0i3ibZAiOt5OyV8tIVqsJy+hH8ZSInYd2atGMRbp27rVq62zTniN7tOfIHk0ePVmVIyvV2tGqzXs3q6mlSXPvnStJKi0plZQN2XuO7NGoslF6ZMsjWnvlWjW1NGnz3s2ZPls7WjP7Pdx2OFPJbmppUmtHqxbNWJQJ2dUV1QN2DDB0neiJeGnb9abPU/FEQCrMAAqXr+rbXVW4kApx2lQQoB/F85fjSvK8qnl6fvfzKi0pVVtnm0pUok51qnZirRp2N6ispExtnW2aPHqy9hzZo/LScklSS1v4cJ+4vSRVlFWopa1FFWUVau1ozWxXXVGtF5pfUEdnR5e2Y4aP0c5DOzV59GSqyyiqrYuv1vTV9x23vDdTKQptO1imZ1BhBtA7hVzGLV+Y7a4inW+ecr5t48pzd/0D/WB5zXLtPLRTtRPD/5MNu0Ohpq2zTWUlZZJCCE6GZUnafWi35lXNO66/OABLyoTllrYWtXW2qaKsQnuO7FHD7gZVjqzMVKhvmHdDJixLododnxwIFEMclnMrwLnBNnf91sVXZ5antW2+/Y7M40PPPXfcssGMCjOAgZdWSabKjAEy9965mXnKcSU4Drrx8pa2FpVEtaVOdapEJZkq9OTRkzOV4WTV+XDb4eOqzXF/8bpk/0llJWV6fvHzA3sgMCTF4TYZcrcuvlqjFyzoVeU3NyTH4Tm3n2S7fPsuFk76AzC4EIQxiKxoXKG7Xrgr87hEJRpVNioTnGOTR0/W7kO7JYW5y/EJfAePHewSdpMBOFmRzhWH647ODo0qG3XclI4SlXCVDAyYl997oWY++eNMhTdfdThf6E1WhNMCb/Ptd2jfww+r48ABzWqo77ZdMRGYAQBIEZ/EJ4XpEXe/cHemMhyH3oqyCh1uOyzp+ApzXD2WlLlyxu5Du7tMz4j7ibeLQ3VcnY63javOkqgwY0DEYTkprg5PX31fJhxvrq3TrIb6zPzmePnL771Q4y+7rMv28TQLKVSY42Xxdrn3c9sUA3OYAZx6+DhsDJCljy3NBFRJuuuFu7oE3bhC3NLWonMmnnNcWJak1o5WlZWUZarTOw/tzATfirKKLpXmTnV2uWJGHJbjdfE+u6tMA30pt7IshWpwvtA6YvZsSSEAx2F56+KrNfPJH2vfww/r0HPPZW6jFyzQ0U2bNH31fTr03HPa9/DDmr76Pr383gsldQ3icVje9/DDamtqGvTzmPmkPwCDB1M2MECe3/18lxBcObLyuGkWsfh6yp3qVGdn11BdO7FWz+9+PrNdXDnODb7JkwAlHTcvOjk9A+hPuZXh3Mrv6AULMnOZ4xCbDLNbF1+dCbgzn/yxNs19qyqXLcusLx07VlL2BMKX33uhyqrD5RLjanVyWVl1ddGqy73BlAwAwJA09965XS4Bl5x/nHvSX7L6LHW9dJyUPdEvt21ZSVnmMnLxZeZyLysXbxuHb+Ywoz/lTrGIA/Tr994rb2nR7M2btHXx1Tqybp1mb9zQpap8pLFRs9e/qM21dZnKc1tTk9qbmyVJVh5OfI1DcbKKHU/XSE7DSFaZx192WVFOAGQOMwAA3Vj62FK90PxCQVMgcsNxspq858ieHvsoKylT5cjKLtMwck8MjIN2RVmF1l659kSeElCQrYuv1tFNmzRi9my1NTVp5pM/zgTkYZMmZSq/Rzdt0hlLlmjfww+rfdcu2ahRKh07Vu27dmnk/PmZsCuF8HuksVFSCM1++HAmbMcn/cWVZylUlduawkfEtzc3a2RNTa+vztFXmMMMAEA3Nu/drPLSclWUVWjy6Mkqif7FoTjpcNthlahErR2tmXnLUrgmc0dnR2aRo14WAAAf3klEQVR6Rbx88ujJmW3jdQePHVRZSVmXedPx+mRYTn4qINAfjqxbJz98WEcaG9W+Y4c2zXmz2pqaNHL+/LC+sTETlvesXBk2Kg3XDm/ftUuKpiXtufPOzNzlI42NUkeHKpctkx8+rJHz52tzbZ32rFypsurqELR37ND4yy5T+44dmXDdceCApFClLvbl5XpCYAYADDmL5yzW4bbDmStgxB8mEps8enIm3I4qG6VOdaqjs0OVIytVXlreZV28bfJDSqQQoEeVjcqc7Deval7m+s1x23MmnqNzJp6jyaMnq6WtJXMtZ6DfdHZKnZ0aWVOjkXV1kkKV90hjY6ge19TIW1v1+r33Sh0dat+xQ+rokCTZqFEaWVenI42NGjZlio6sW6cj69ZpZE2NKm+4IUy76OzMhHIpVKrLqqtlFRUhgJeUSG1tam9u1ojZs1W5bFkI4oMcgRkAMCSVlpTqunnX6eCxg6ocWalzJp6jWRNmHXfJtzHDx+iGeTdo4uiJevzyxzVrwqzMPOS1V65VR2dHplI9ZviYTOjt6OzIVKc7OjvU1NKkirKKzHWdD7cdVlNLk+om1engsYNFOw4YWqwivIsSV5KHTZoktbWFW2dnZqpE6dixslGjpLIyqbNTfviwvKUlVIfb2kKQlkJArq/XnpUrdaQ+utZyaWnoV9IZS5bo6KZN8pYWDauqkkpLwxja2jLbxRXswYzADAAYcpbXLNfzi5/XPevv0awJs1RdUa3NezdLkvYc2aMb5t2ga+deqxeWvKCdh3aqfle9Dh47qPPuP091k+p0zsRzVFZSpnn3ztPE0RMzfRw8dlCtHa2qKKvQdfOu08TREzVx9ESNKhul6oowN3Ti6ImZAF5dUa3VG1frcNvhzEdlA/3JW6KrsbS1yVtassG3pEQqKclUlNt37JC3toYgLYXKdFmZrLw8BN6yskwojperLJpylKhM77nzztBPWVlmCoa3tkolJRo2ZUqoXrcN/sspEpgBAEPOisYVWtG4osuy1o5Wbd67WZUjK3XP+nu0vGa5LnroIt0w7wZt3rtZY4aP0awJs/TIlkckSfOq5umcied0qQ63tLXo2rnXataEWbpn/T1aNGOR9hzZo1kTZmnz3s1aPGexHr/8cdXvqteiGYskSbMmzNJ1865T/a56PX754wN3EICyxJz60lKNnD8/hNjS0vBVyoTjuCrsLS2Z6Rbx1TEkZSvIJSXZIB1fhjGqYHtLSzYgd3aG7aMgHl+rebAiMAMAhpzlNcu1vGZ5lznD86rmafGcxVo0Y5HKS8u19LGlWjRjkZbXLM9Ufusm1am6olpNLeFt6/gazZK0auEq3TDvhky7ypGVWr1xtSpHVqpuUp3WXrlWqzeulqRMNbuppUmrFq7KhPDcEA/0p/gycMkpEu27dkltbaEa3NGRqTLHIVmS1NkZto0rw9Ec5/YdO0IQ3rGjaxhPSk6/iLdvazvuUwcHGwIzAGBIWvrYUq29cq1WLVyluknh5Kf6XfV6ZMsjmjVhliRp9cbVWtG4QotmLNKiGYtUvyvM0Vw0Y5FWLVylypGVmbZSCOL1u+q1vGa5Fs1YlFm3vGa5lj62NPN47ZVrVb+rXtUV1ZlgvmrhqoF8+hjihk2ZIm9p0ci6uuw0jQRvaQnV33zTJUpKwvo4FMeV5GRI7maaRRzSR9bVZdpbRQWf9AcAwGATV3JXNK7Q8prlkpQJrMllSUsfW5pps/SxpVrRuELVFdVatXCVlj629Li2dZPqtGrhKq1oXJHZNreCnNtfHMiB/hZfmeJIfb2soqJraC4r07Cqquz8ZimE4pISqbNTwyZN6nL1jHj5cSE5Xp4Q7ydzgmBZWahiD3J8cAkAYMhKhuM4zMaP49C7vGZ5pl0causm1WXCbVydzheyc9unVZGTgRzoL5tmzZZKSmSjRh0XktXRkQ24cZvW1q7LIyPr6nRk3brMB5rEH07SJWTHcoJz5cc/rj133pkN6iUlmr1xQ3883R4V+sElVJgBAMgjt+KcG2iT1eJkqE7qLkQnt4nvE5YxIMrKsifglZSEOcXRiXjDpkzJfjhJZ6e8tTVUmnftyoRbq6jIfPBJHLrbo77aowA+sq5OR+rru/RnFRWZ+dB77rorTOs4fFgqKdHI+fMzH6U9WDGHGQAwZCWnQOQLt8nluYE2PnGwEPG0jNztcrfnpD8MhOT8YSl8IImk4y4xF19vedikSZmPth4xe/bx0y9yKshHGhulkhJ1HDiQ6Tsz7SK6rrONGiUbNUrDJk3SkcbGQR2WJQIzAGAIyxeCpeODayFBNi08L69Znpm6kewvt99CAzhwwjo6dHTTpsw0CysvD5eEmzJFw6ZMyQTc+ANGrKJC7c3NmTB9pL4+hN74xD+pS1i2iorMiX3JS9ANq6qSOjq6rIuncFh5+aA/6Y/ADABAjmRwzTfPOXkd57SQm+/kwtz9dBfSgf6QrBZLykzNKKuuVntzc5iGEX/6X7S+y+XnSkJ0jIN1HKqHTZmSmWYxYvbs7Cf6dXZqZF2dOg4c0Mj587uMpX3XrsxVOqo+ceNAPP0TRmAGACBFMugmp1IUEnR7muZRSFugL3UcONAlvFpFhUbOn68jjY2hCqzoChplZRoxe3aYviFp9uZNKh07Nny4SVVVCMUKUzT88GGVVYdPshw5f77ampo0sqYms/2R+nqNmD1bR9atkxQuaRfvd/rq+06Jy8oVFJjNbKGZvWRmW8zspjzrrzez9WbWaGY/M7M5iXU3R9u9ZGYf6MvBAwBQbIUEXcIwBosRs2frjCVL1NbUJKuo0BlLlmj0ggUaWVOj9uZmVS5bptkbN8jKy8PUDUl++LC2Lr5akjR6wQKVVVdr9IIFqvz4x9XW1KTZGzfoyLp1qrzhhsw2bU1NmVscukfOn6/SsWNVVl2d2e/L771QsxrqT/0Ks5mVSrpT0gclzZH00WQgjtzv7nPdvUbSVyV9Ldp2jqQrJL1Z0kJJK6L+AAA47fRmWgVTMFAM01ffp0PPPafxl12mEbNna8/KlTr03HOavvo+DauqUtUnbtTWxVdrxOzZGjF7dgjT8+dngrIUQrMkvX7vveo4cEDNt9+hkfPn69Bzz+mMJUs088kfq6y6WmXV1eo4cCDTPmnfww9Lksqqq9V8+x2DvsJcyGXlFkja4u6vSpKZPSDpUkkb4wbufiDRfrSk+OLOl0p6wN1bJf23mW2J+lvbB2MHAGBQ6U0lubtL0QH96eX3Xqiy6upQ0b39Dk1ffZ8kqfn2OzT+ssskherw+Msuy4TaZODNrQRnHkeBd8+dd2aWTV99n7Yuvlp77rorzIuWNP6yy1T1iRvVfPsdOvTcc2pratKelSs1e/2L/fek+0Ahgbla0rbE4+2Szs1tZGYfl/RpScMlvTex7bM521af0EgBABgE+jLkEpYx0GY++WNtXXx1pqL78nsvzATlWPw4/hoH3Omr7+tSCY6Dcbzs9Xvv1ci6Om1dfHUmiI9esCATwONt4ukdybA+2BUSmC3PsuM+HtDd75R0p5ldKenzkpYUuq2ZLZO0TJKmTZtWwJAAACgOQi5OZckwK2XDcG74TVaS4/XdheX4fjIMJ5fHFe2ti6+Wbr8j7xSNZPvBqJDAvF3S1MTjMyXl+dzDjAck3dWbbd19paSVUvho7ALGBAAAgF6Kw3IcnJNBdXNtnWY11GeCcTI85wbt3GAd95Nsk9kmmqOcXJc0mINyzNzT86mZDZP0a0kXSmqSVC/pSnffkGgz091fju7/oaRb3L3WzN4s6X6FectTJP1Y0kx37+huf7W1td7Q0HByzwoAAAAFS1aRuwuw+SrMJ7O/vujnZJnZOnev7aldj1fJcPd2STdK+pGkTZIedPcNZnarmV0SNbvRzDaYWaPCPOYl0bYbJD2ocILgY5I+nhaWAQAAMPCS0yqk/POKqz5xY+bWndzqdHfr0vYzGPVYYR5oVJgBAAAGn8E+z/hE9FmFGQAAAKen3lR484XltIry6YQKMwAAAE7YqVx5psIMAACAfneiYflUqkoTmAEAANDnegrEp1JVmsAMAACAPncqBeKeEJgBAACAFARmAAAAIAWBGQAAAEhBYAYAAABSEJgBAACAFARmAAAAIAWBGQAAAEhBYAYAAABSEJgBAACAFARmAAAAIAWBGQAAAEhBYAYAAABSEJgBAACAFARmAAAAIAWBGQAAAEhBYAYAAABSEJgBAACAFARmAAAAIAWBGQAAAEhBYAYAAABSEJgBAACAFARmAAAAIAWBGQAAAEhBYAYAAABSEJgBAACAFARmAAAAIEVBgdnMFprZS2a2xcxuyrP+02a20cxeNLMfm9lZiXUdZtYY3db05eABAACA/jaspwZmVirpTknvl7RdUr2ZrXH3jYlmv5RU6+6HzewGSV+V9JFo3RF3r+njcQMAAAADopAK8wJJW9z9VXc/JukBSZcmG7j7U+5+OHr4rKQz+3aYAAAAQHEUEpirJW1LPN4eLevOxyT9MPF4hJk1mNmzZrboBMYIAAAAFE2PUzIkWZ5lnreh2Z9KqpX07sTiae6+w8zeKOlJM1vv7q/kbLdM0jJJmjZtWkEDBwAAAAZCIRXm7ZKmJh6fKWlHbiMze5+kz0m6xN1b4+XuviP6+qqkpyW9LXdbd1/p7rXuXltVVdWrJwAAAAD0p0ICc72kmWZ2tpkNl3SFpC5XuzCzt0m6WyEs/zax/A1mVh7dr5T0DknJkwUBAACAQa3HKRnu3m5mN0r6kaRSSd9y9w1mdqukBndfI+kfJVVI+jczk6TX3P0SSbMl3W1mnQrh/Cs5V9cAAAAABjVzzzsduWhqa2u9oaGh2MMAAADAac7M1rl7bU/t+KQ/AAAAIAWBGQAAAEhBYAYAAABSEJgBAACAFARmAAAAIAWBGQAAAEhBYAYAAABSEJgBAACAFARmAAAAIAWBGQAAAEhBYAYAAABSEJgBAACAFARmAAAAIAWBGQAAAEhBYAYAAABSEJgBAACAFARmAAAAIAWBGQAAAEhBYAYAAABSEJgBAACAFARmAAAAIAWBGQAAAEhBYAYAAABSEJgBAACAFARmAAAAIAWBGQAAAEhBYAYAAABSEJgBAACAFARmAAAAIAWBGQAAAEhBYAYAAABSEJgBAACAFAUFZjNbaGYvmdkWM7spz/pPm9lGM3vRzH5sZmcl1i0xs5ej25K+HDwAAADQ33oMzGZWKulOSR+UNEfSR81sTk6zX0qqdfe3SnpI0lejbSdIukXSuZIWSLrFzN7Qd8MHAAAA+lchFeYFkra4+6vufkzSA5IuTTZw96fc/XD08FlJZ0b3PyDpCXff6+6/k/SEpIV9M3QAAACg/xUSmKslbUs83h4t687HJP2wN9ua2TIzazCzhubm5gKGBAAAAAyMQgKz5VnmeRua/amkWkn/2Jtt3X2lu9e6e21VVVUBQwIAAAAGRiGBebukqYnHZ0rakdvIzN4n6XOSLnH31t5sCwAAAAxWhQTmekkzzexsMxsu6QpJa5INzOxtku5WCMu/Taz6kaSLzOwN0cl+F0XLAAAAgFPCsJ4auHu7md2oEHRLJX3L3TeY2a2SGtx9jcIUjApJ/2ZmkvSau1/i7nvN7EsKoVuSbnX3vf3yTAAAAIB+YO55pyMXTW1trTc0NBR7GAAAADjNmdk6d6/tqR2f9AcAAACkIDADAAAAKQjMAAAAQAoCMwAAAJCCwAwAAACkIDADAAAAKQjMAAAAQAoCMwAAAJCCwAwAAACkIDADAAAAKQjMAAAAQAoCMwAAAJCCwAwAAACkIDADAAAAKQjMAAAAQAoCMwAAAJCCwAwAAACkIDADAAAAKQjMAAAAQAoCMwAAAJCCwAwAAACkIDADAAAAKQjMAAAAQAoCMwAAAJCCwAwAAACkIDADAAAAKQjMAAAAQAoCMwAAAJCCwAwAAACkIDADAAAAKQjMAAAAQIqCArOZLTSzl8xsi5ndlGf9u8zseTNrN7PLc9Z1mFljdFvTVwMHAAAABsKwnhqYWamkOyW9X9J2SfVmtsbdNyaavSbpGkl/maeLI+5e0wdjBQAAAAZcj4FZ0gJJW9z9VUkyswckXSopE5jdfWu0rrMfxggAAAAUTSFTMqolbUs83h4tK9QIM2sws2fNbFGvRgcAAAAUWSEVZsuzzHuxj2nuvsPM3ijpSTNb7+6vdNmB2TJJyyRp2rRpvegaAAAA6F+FVJi3S5qaeHympB2F7sDdd0RfX5X0tKS35Wmz0t1r3b22qqqq0K4BAACAfldIYK6XNNPMzjaz4ZKukFTQ1S7M7A1mVh7dr5T0DiXmPgMAAACDXY+B2d3bJd0o6UeSNkl60N03mNmtZnaJJJlZnZltl/THku42sw3R5rMlNZjZC5KekvSVnKtrAAAAAIOaufdmOnL/q62t9YaGhmIPAwAAAKc5M1vn7rU9teOT/gAAAIAUBGYAAAAgBYEZAAAASEFgBgAAAFIQmAEAAIAUBGYAAAAgBYEZAAAASEFgBgAAAFIQmAEAAIAUBGYAAAAgBYEZAAAASEFgBgAAAFIQmAEAAIAUBGYAAAAgBYEZAAAASEFgBgAAAFIQmAEAAIAUBGYAAAAgBYEZAAAASEFgBgAAAFIQmAEAAIAUBGYAAAAgBYEZAAAASEFgBgAAAFIQmAEAAIAUBGYAAAAgBYEZAAAASEFgBgAAAFIQmAEAAIAUBGYAAAAgBYEZAAAASEFgBgAAAFIUFJjNbKGZvWRmW8zspjzr32Vmz5tZu5ldnrNuiZm9HN2W9NXAAQAAgIHQY2A2s1JJd0r6oKQ5kj5qZnNymr0m6RpJ9+dsO0HSLZLOlbRA0i1m9oaTHzYAAAAwMAqpMC+QtMXdX3X3Y5IekHRpsoG7b3X3FyV15mz7AUlPuPted/+dpCckLeyDcQMAAAADopDAXC1pW+Lx9mhZIQra1syWmVmDmTU0NzcX2DUAAADQ/woJzJZnmRfYf0HbuvtKd69199qqqqoCuwYAAAD6XyGBebukqYnHZ0raUWD/J7MtAAAAUHSFBOZ6STPN7GwzGy7pCklrCuz/R5IuMrM3RCf7XRQtAwAAAE4JPQZmd2+XdKNC0N0k6UF332Bmt5rZJZJkZnVmtl3SH0u628w2RNvulfQlhdBdL+nWaBkAAABwSjD3QqcjD4za2lpvaGgo9jAAAABwmjOzde5e21M7PukPAAAASEFgBgAAAFIQmAEAAIAUBGYAAAAgBYEZAAAASEFgBgAAAFIQmAEAAIAUBGYAAAAgBYEZAAAASEFgBgAAAFIQmAEAAIAUBGYAAAAgBYEZAAAASEFgBgAAAFIQmAEAAIAUBGYAAAAgBYEZAAAASEFgBgAAAFIQmAEAAIAUBGYAAAAgBYEZAAAASEFgBgAAAFIQmAEAAIAUBGYAAAAgBYEZAAAASEFgBgAAAFIQmAEAAIAUBGYAAAAgBYEZAAAASEFgBgAAAFIQmAEAAIAUBQVmM1toZi+Z2RYzuynP+nIz+160/hdmNj1aPt3MjphZY3T7Zt8OHwAAAOhfw3pqYGalku6U9H5J2yXVm9kad9+YaPYxSb9z9xlmdoWkf5D0kWjdK+5e08fjBgAAAAZEIRXmBZK2uPur7n5M0gOSLs1pc6mke6P7D0m60Mys74YJAAAAFEchgbla0rbE4+3Rsrxt3L1d0n5JZ0TrzjazX5rZT8zs/JMcLwAAADCgepySISlfpdgLbLNT0jR3f93M5kt6xMze7O4HumxstkzSMkmaNm1aAUMCAAAABkYhFebtkqYmHp8paUd3bcxsmKRxkva6e6u7vy5J7r5O0iuS3pS7A3df6e617l5bVVXV+2cBAAAA9JNCAnO9pJlmdraZDZd0haQ1OW3WSFoS3b9c0pPu7mZWFZ00KDN7o6SZkl7tm6EDAAAA/a/HKRnu3m5mN0r6kaRSSd9y9w1mdqukBndfI+lfJK02sy2S9iqEakl6l6RbzaxdUoek6919b388EQAAAKA/mHvudOTiqq2t9YaGhmIPAwAAAKc5M1vn7rU9teOT/gAAAIAUBGYAAAAgBYEZAAAASEFgBgAAAFIQmAEAAIAUBGYAAAAgBYEZAAAASEFgBgAAAFIQmAEAAIAUBGYAAAAgBYEZAAAASEFgBgAAAFIQmAEAAIAUBGYAAAAgBYEZAAAASEFgBgAAAFIQmAEAAIAUBGYAAAAgBYEZAAAASEFg7gdff+LXxR4CAAAA+giBuRe+/sSvCwrDn3r/m3rs50TWAQAAYOANK/YABpOvP/HrTNhNBtdPvf9NeYNs7rJ828aeffV1fe+6845b9/Unfp1Z95G71+rtbzxDH7l7rSTp7W88o8fwDQAAgP5l7l7sMXRRW1vrDQ0NRdn33FsekySNHVkmSWo+2Kob3jNDtz/5siRp9PBSHTrWoU6XhpeajnV0PXbV40dIkpr2Hc08ju9L+bcZU16qg60dqh4/Qs0HWzPL2ztdw0pC+zHlpVr/xYV9/GwBAACGNjNb5+61PbWjwpwQh+GDrR2ZZf/rxy9n7ieX5wZfSV3Ccb7H+baJ+8xtm2wfB3gAAAAMPOYwJ3QOrmJ7Rr4wDQAAgIFBYAYAAABSEJhPAeeePaHYQwAAABiyCMwJ5549QSUWTsT75IUzde7ZEzS81DKPP3nhTFWPH6ESCyfwDS81nXv2hEy7c8+eoOrxI467P6a8VNXjR2S23/qVD6l6/AiNKS/VuWdP0CcvnKnhpaYSU6ZN3Ickfe+684p8ZAAAAIYurpLRC/El4Qq91Ftu++Rl6wAAAFBchV4lg8AMAACAIanQwMyUDAAAACBFQYHZzBaa2UtmtsXMbsqzvtzMvhet/4WZTU+suzla/pKZfaDvhg4AAAD0vx4Ds5mVSrpT0gclzZH0UTObk9PsY5J+5+4zJH1d0j9E286RdIWkN0taKGlF1B8AAABwSiikwrxA0hZ3f9Xdj0l6QNKlOW0ulXRvdP8hSReamUXLH3D3Vnf/b0lbov4AAACAU0Ihgbla0rbE4+3Rsrxt3L1d0n5JZxS4LQAAADBoFRKYLc+y3EtrdNemkG1lZsvMrMHMGpqbmwsYEgAAADAwCgnM2yVNTTw+U9KO7tqY2TBJ4yTtLXBbuftKd69199qqqqrCRw8AAAD0s0ICc72kmWZ2tpkNVziJb01OmzWSlkT3L5f0pIcLPK+RdEV0FY2zJc2U9FzfDB0AAADof8N6auDu7WZ2o6QfSSqV9C1332Bmt0pqcPc1kv5F0moz26JQWb4i2naDmT0oaaOkdkkfd/eOfnouAAAAQJ/jk/4AAAAwJPFJfwAAAEAfIDADAAAAKQjMAAAAQAoCMwAAAJBi0J30Z2bNkn5T7HH0gUpJe4o9iCGKY188HPvi4dgXD8e+eDj2xXU6HP+z3L3HDwEZdIH5dGFmDYWcdYm+x7EvHo598XDsi4djXzwc++IaSsefKRkAAABACgIzAAAAkILA3H9WFnsAQxjHvng49sXDsS8ejn3xcOyLa8gcf+YwAwAAACmoMAMAAAApCMwnwcy+ZWa/NbNfdbPezOyfzGyLmb1oZucM9BhPVwUc+/eY2X4za4xufzvQYzxdmdlUM3vKzDaZ2QYz+2SeNrz2+0GBx57Xfj8wsxFm9pyZvRAd+y/maVNuZt+LXve/MLPpAz/S00+Bx/4aM2tOvO6vLcZYT1dmVmpmvzSzH+RZNyRe98OKPYBT3Lcl3SHpvm7Wf1DSzOh2rqS7oq84ed9W+rGXpGfc/cMDM5whpV3SZ9z9eTMbI2mdmT3h7hsTbXjt949Cjr3Ea78/tEp6r7u3mFmZpJ+Z2Q/d/dlEm49J+p27zzCzKyT9g6SPFGOwp5lCjr0kfc/dbyzC+IaCT0raJGlsnnVD4nVPhfkkuPtPJe1NaXKppPs8eFbSeDObPDCjO70VcOzRT9x9p7s/H90/qPBLtDqnGa/9flDgsUc/iF7LLdHDsuiWexLQpZLuje4/JOlCM7MBGuJpq8Bjj35iZmdK+pCke7ppMiRe9wTm/lUtaVvi8Xbxn9tAOi96C++HZvbmYg/mdBS99fY2Sb/IWcVrv5+lHHuJ136/iN6WbpT0W0lPuHu3r3t3b5e0X9IZAzvK01MBx16S/iiaAvaQmU0d4CGezr4h6bOSOrtZPyRe9wTm/pXvLyz+Kh4Yzyt83OU8SbdLeqTI4zntmFmFpH+X9P+5+4Hc1Xk24bXfR3o49rz2+4m7d7h7jaQzJS0ws7fkNOF1308KOPbflzTd3d8q6T+VrXjiJJjZhyX91t3XpTXLs+y0e90TmPvXdknJv3LPlLSjSGMZUtz9QPwWnrs/KqnMzCqLPKzTRjSP8N8lfcfd/yNPE177/aSnY89rv/+5+z5JT0tamLMq87o3s2GSxompY32qu2Pv7q+7e2v08J8lzR/goZ2u3iHpEjPbKukBSe81s3/NaTMkXvcE5v61RtLV0RUD3i5pv7vvLPaghgIzmxTPoTKzBQqv9deLO6rTQ3Rc/0XSJnf/WjfNeO33g0KOPa/9/mFmVWY2Pro/UtL7JG3OabZG0pLo/uWSnnQ+7OCkFXLsc86RuERhfj9Okrvf7O5nuvt0SVcovKb/NKfZkHjdc5WMk2Bm35X0HkmVZrZd0i0KJyPI3b8p6VFJF0vaIumwpKXFGenpp4Bjf7mkG8ysXdIRSVecjj/ARfIOSYslrY/mFErSX0uaJvHa72eFHHte+/1jsqR7zaxU4Y+QB939B2Z2q6QGd1+j8MfMajPbolBhu6J4wz2tFHLs/8LMLlG4ksxeSdcUbbRDwFB83fNJfwAAAEAKpmQAAAAAKQjMAAAAQAoCMwAAAJCCwAwAAACkIDADAAAAKQjMAAAAQAoCMwAAAJCCwAwAAACk+L/zWb5QL1ntIwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "samples = results.extract()\n",
    "plt.figure(figsize=(12, 6))\n",
    "for i in range(4):\n",
    "    plt.plot(sts.uniform.rvs(loc=i+1-0.1, scale=0.2, size=samples['theta'].shape[0]), \n",
    "             samples['theta'][:,i], ',', alpha=0.5)\n",
    "plt.title('Distribution over each component of θ for p(son | father is unskilled)')\n",
    "plt.show()"
   ]
  },
  {
   "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}