MBML_Chapter2_PeoplesSkills.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Model-based machine learning - Chapter 2 - Assessing People’s Skills"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook is based on the [Model-based machine learning - Chapter 2 - Assessing People’s Skills](http://mbmlbook.com/LearningSkills.html) and this [forum thread](https://forum.pyro.ai/t/model-based-machine-learning-book-chapter-2-skills-example-in-pyro-tensor-dimension-issue/464) using [Pyro](http://pyro.ai) and [PyTorch](https://pytorch.org).<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import pyro\n",
    "import pyro.infer\n",
    "import pyro.optim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pdb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'0.3.0'"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pyro.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert pyro.__version__.startswith('0.3.0')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "pyro.enable_validation(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "pyro.set_rng_seed(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Raw responses data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"Raw data collected from volunteers completing a real assessment test. This data consists of the self-assessed skills (S1-S7) and the answers to each question (Q1-Q48). The first row of data gives the correct answers to each question. Each subsequent row gives the data for one of the participants.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw_responses = pd.read_csv('http://mbmlbook.com/Downloads/LearningSkills_Real_Data_Experiments-Original-Inputs-RawResponsesAsDictionary.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "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>#</th>\n",
       "      <th>S1</th>\n",
       "      <th>S2</th>\n",
       "      <th>S3</th>\n",
       "      <th>S4</th>\n",
       "      <th>S5</th>\n",
       "      <th>S6</th>\n",
       "      <th>S7</th>\n",
       "      <th>Q1</th>\n",
       "      <th>Q2</th>\n",
       "      <th>...</th>\n",
       "      <th>Q39</th>\n",
       "      <th>Q40</th>\n",
       "      <th>Q41</th>\n",
       "      <th>Q42</th>\n",
       "      <th>Q43</th>\n",
       "      <th>Q44</th>\n",
       "      <th>Q45</th>\n",
       "      <th>Q46</th>\n",
       "      <th>Q47</th>\n",
       "      <th>Q48</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>ANS</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>2</td>\n",
       "      <td>4</td>\n",
       "      <td>...</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>4</td>\n",
       "      <td>2</td>\n",
       "      <td>5</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>P1</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>2</td>\n",
       "      <td>4</td>\n",
       "      <td>...</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>4</td>\n",
       "      <td>2</td>\n",
       "      <td>5</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>P2</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "      <td>...</td>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>4</td>\n",
       "      <td>4</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>P3</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>3</td>\n",
       "      <td>4</td>\n",
       "      <td>...</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>4</td>\n",
       "      <td>2</td>\n",
       "      <td>5</td>\n",
       "      <td>5</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>P4</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>2</td>\n",
       "      <td>4</td>\n",
       "      <td>...</td>\n",
       "      <td>2</td>\n",
       "      <td>5</td>\n",
       "      <td>5</td>\n",
       "      <td>4</td>\n",
       "      <td>2</td>\n",
       "      <td>5</td>\n",
       "      <td>4</td>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 56 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     #     S1     S2     S3    S4     S5     S6     S7  Q1  Q2 ...   Q39  Q40  \\\n",
       "0  ANS    NaN    NaN    NaN   NaN    NaN    NaN    NaN   2   4 ...     2    1   \n",
       "1   P1   True   True   True  True  False   True   True   2   4 ...     2    1   \n",
       "2   P2   True   True   True  True  False  False   True   1   4 ...     5    1   \n",
       "3   P3  False  False  False  True  False   True  False   3   4 ...     2    1   \n",
       "4   P4   True   True  False  True  False  False  False   2   4 ...     2    5   \n",
       "\n",
       "   Q41  Q42  Q43  Q44  Q45  Q46  Q47  Q48  \n",
       "0    5    4    2    5    2    2    2    1  \n",
       "1    5    4    2    5    2    2    2    1  \n",
       "2    1    4    3    3    4    4    3    1  \n",
       "3    5    4    2    5    5    3    3    1  \n",
       "4    5    4    2    5    4    2    3    1  \n",
       "\n",
       "[5 rows x 56 columns]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "raw_responses.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Skills for questions data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "skills_questions_mask = pd.read_csv('http://mbmlbook.com/Downloads/LearningSkills_Real_Data_Experiments-Original-Inputs-Quiz-SkillsQuestionsMask.csv', header=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "#skills_questions_mask.T.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs4AAACfCAYAAADpoV9XAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHwlJREFUeJzt3XmYJFWd7vHvCw1Iy1ouyCY9yCaU2mqLOiIy4uUqg+KCl/GiyGYrXtxRcfCRRVEUHJFhVEoUmNFpRkURcWvEC4KsDUKTiIJoMyCISCkNIkjTv/njRElSnXHyRHZFLd3v53nyqVzizXMy45eZJ6NORCoiMDMzMzOzvDWmugNmZmZmZjOBB85mZmZmZgU8cDYzMzMzK+CBs5mZmZlZAQ+czczMzMwKeOBsZmZmZlbAA2czs4qkAyRdUnPbfpIWdl0OSdtU58+Q9LEB2nuRpJsl3S/p1YP3fHqRdIOk3aa6H2ZmE80DZzNbrUjaRdKlku6VNCrpp5Ke1y8XEV+NiD0muDvHAqdExHoRcc4E3/ek6PWlISJ2iogLp6hLZmatmTXVHTAzmyySNgDOAw4FvgasDbwYeGiKurQVcMMgQUmzImLZBPfHzMwyvMXZzFYn2wFExIKIeCQi/hIRCyNica+FJZ0g6RJJG+amcYzLPFHSeZL+VG3RvljSCu+1km4Btga+U03VWEfSZpLOrXK/kvSWruWPlvQNSV+RtBQ4oMd9PqHKL5V0paSPjvVZ0pxqesmsruUvlHRI1+WDJN0o6Y+Sfihpq+p6SfqMpN9XW+oXSxqWNB/YD/hA9Ri+Uy2/RNLLqvPrSDpJ0h3V6SRJ61S37Sbpdknvq+77TkkHdvVnT0k/l3SfpN9KOrzf829m1iYPnM1sdXIT8IikMyW9QtLGvRaStIakLwLPBPaIiHsbtPE+4HbgScAmwD8DMX6hiHga8N/AK6upGg8BC6rsZsA+wMcl7d4V2xv4BrAR8NUebf8b8CCwKXBQdSpSzbH+Z+C1Vd8vrvoDsAewK+mLx0bAvsA9ETFS9eNT1WN4ZY+7PhJ4ATAXeBawM/DhrtufAmwIbA4cDPxb13r5EvDWiFgfGAZ+XPp4zMza4IGzma02ImIpsAtpIPtF4O5qC+0mXYutRRowDpEGtQ80bOZh0sB1q4h4OCIujogVBs7jSdqy6tsHI+LBiLgWOA14U9dil0XEORGxPCL+Mi6/JvA64CMR8eeI6ABnNuj3W4FPRMSN1RSQjwNzq63ODwPrAzsAqpa5s/B+9wOOjYjfR8TdwDHjHtPD1e0PR8T3gPuB7btu21HSBhHxx4i4psHjMTObcB44m9lqpRr0HRARW5C2Ym4GnNS1yDakLbvHRMRfB2jiBOBXwEJJv5Z0RGFuM2A0Iu7ruu5W0pbYMbdl8k8i7bfSvcythW1Dmm/92WqKyZ+AUUDA5hHxY+AU0hbtuySNVPPFS2w2rh+3VteNuWfcXO0HgPWq868D9gRulXSRpBc2eDxmZhPOA2czW21FxC+AM0gD6DE3AgcC35e0fa9cn/u8LyLeFxFbA68E3jtuukWdO4AhSet3XfdU4Lfdd5/J3w0sA7Yclx/z5+rv7K7rntJ1/jbStIiNuk7rRsSl1eM6OSKeC+xEmrLx/oI+jT2urcb16Y4+Gao2r4qIvYEnA+eQdug0M5syHjib2WpD0g7VjmhbVJe3BN4AXN69XEQsIM33/ZGkpzVsYy9J20gSsBR4pDplRcRtwKXAJyQ9TtIzSXN+e81l7pV/BPgmcLSk2ZJ2BN7cdfvdpEH4GyWtKekgoPuxfQH4kKSdqsexoaTXV+efJ+n5ktYiDcAf7HpMd5F2cqyzAPiwpCdJeiLwEeAr/R6PpLWVjp29YUQ8zKPPpZnZlPHA2cxWJ/cBzweukPRn0oC5Q9qh7zEi4kzScZZ/LGlOgza2BX5Emqt7GfC5Bsc0fgMwh7RF9lvAURFxfoO2DyNNc/gdaUv66eNufwtpS/E9pC3Hl47dEBHfAj4JnFUdtaMDvKK6eQPSnPA/kqZa3AOcWN32JdI85D9J6nUs6o8Bi4DFwPXANdV1Jd4ELKn68zbgjYU5M7NWqGCfFTMzm4EkHQAcEhG7THVfzMxWBd7ibGZmZmZWwANnMzMzM7MCnqphZmZmZlbAW5zNzMzMzAp44GxmZmZmVmDWVHcgw3NIzMzMzGwyqGSh6TxwZnR0tHFmaGiocW5oaGig9laH3PDwcP8Fx+l0Os4551yVAQbOzYT3CBjs8U33deecc9M5t7LvLatybtB1UMpTNczMzMzMCnjgbGZmZmZWoLWpGpJ2APYGNifNV74DODcibmyrTTMzMzOztrSyxVnSB4GzSBOtrwSuqs4vkHREG22amZmZmbWprS3OBwM7RcTD3VdK+hfgBuD4lto1MzMzM2tFW3OclwOb9bh+0+q2niTNl7RI0qKRkZGWumZmZmZm1lxbW5zfDVwg6Wbgtuq6pwLbAIfVhSJiBBgbMccgh6MzMzMzM2tDKwPniPiBpO2AnUk7Bwq4HbgqIh5po00zMzMzsza1dlSNiFgOXN7W/ZuZmZmZTSYfx9nMzMzMrIAHzmZmZmZmBRQRU92HOtO2Y2ZmZma2SlHJQt7ibGZmZmZWoLWdAyfC8PBw40yn02mc63Q6A7W3OuQGOSTg0NCQc841yk3Wa32yczPptT7Zuem+7pzL5/xan9rcTHqtT3Zu0HVQyluczczMzMwKeOBsZmZmZlbAA2czMzMzswKTPnCWdOBkt2lmZmZmtrKmYovzMVPQppmZmZnZSmnlqBqSFtfdBGySyc0H5gOceuqpLfTMzMzMzGwwbR2ObhPgfwN/HHe9gEvrQhExAoyMXTz55JPb6Z2ZmZmZWUNtDZzPA9aLiGvH3yDpwpbaNDMzMzNrTSsD54g4OHPb/22jTTMzMzOzNvlwdGZmZmZmBTxwNjMzMzMroIiY6j7UmbYdMzMzM7NVikoW8hZnMzMzM7MCbR1VY0IMDw83znQ6nca5TqczUHurQ26y1oFzq3dudHS0cW5oaGja54aGhgAGzs2E94hBczOlNp1zbjrmZtJrfbJzg66DUt7ibGZmZmZWwANnMzMzM7MCrQ2cJe0gaXdJ6427/uVttWlmZmZm1pZWBs6S3gl8G3gH0JG0d9fNH2+jTTMzMzOzNrW1c+BbgOdGxP2S5gDfkDQnIj5L4eE+zMzMzMymk7YGzmtGxP0AEbFE0m6kwfNWZAbOkuYD8wFOPfXUlrpmZmZmZtZcW3Ocfydp7tiFahC9F/BE4Bl1oYgYiYh5ETFv/vz5LXXNzMzMzKy5tgbO+wO/674iIpZFxP7Ari21aWZmZmbWmlamakTE7ZnbftpGm2ZmZmZmbfJxnM3MzMzMCnjgbGZmZmZWQBFRf6P0HaB2gYh4VRudGrv7Fu/bzMzMzGxM0eGS+81xPnECOmJmZmZmNuNltzhPsRgeHm4c6nQ6NM11Oh0A53rkJmsdOOfcqpibSa/1lcmNjo42yg0NDTXOONdObrq/hpyrz8DMeY+YCeMWJmKLs6TryU/VeGajnpmZmZmZzVD9pmrsNSm9MDMzMzOb5rID54i4ddA7lrRzuou4StKOwMuBX0TE9wa9TzMzMzOzqdJvqsYlEbGLpPt47JQNkQbFG9TkjgJeAcySdD7wfOBC4AhJz46I4yak92ZmZmZmk6TfFuddqr/rN7zffYC5wDqkn97eIiKWSjoBuALwwNnMzMzMZpSiH0CRdHCP647PRJZFxCMR8QBwS0QsBYiIvwDLM+3Ml7RI0qKRkZGSrpmZmZmZTYrSXw7cR9J+YxckfQ54cmb5v0qaXZ1/blduQzID54gYiYh5ETFv/vz5hV0zMzMzM2tfv6NqjHktcK6k5aS5y6MR8fbM8rtGxEMAEdE9UF4LePNAPTUzMzMzm0L9dg4c6rp4CPBt4BLgWElDEdHziOtjg+Ye1/8B+MOAfTUzMzMzmzL9tjhfTTqahrr+7lmdALZur2tmZmZmZtNHv4HzvsBtEXEngKQ3A68DlgBHt9ozMzMzM7NpRBG1v6iNpGuAl0XEqKRdgbOAd5AONff0iNinxb7Vd8zMzMzMbOKoZKF+W5zX7JrHvC8wEhFnA2dLunZlemdmZmZmNpP0HThLmhURy4Ddge5jxJUekWNgw8PDjTOdTqdxrtPpDNTe6pCbrHXgnHOrYm4mvdYnOzfd193qkhsd7bmPf9bQ0JBzNTm/t0x9btB1UKrf4HcBcJGkPwB/AS4GkLQNcG/jnpmZmZmZzVD9fnL7OEkXAJsCC+PRCdFrkOY6m5mZmZmtFvpOt4iIy3tcd1PThiT9e0Ts3zRnZmZmZjYdtDJPWdK5468C/kHSRgAR8ao22jUzMzMza0tbO/htAfwcOI1HfzhlHvDpltozMzMzM2vVGi3d7zzSrw4eCdwbERcCf4mIiyLiorqQpPmSFklaNDIy0lLXzMzMzMyaa2WLc0QsBz4j6evV37tK2oqIEWBsxBwnn3xyG90zMzMzM2us1WMxR8TtwOsl/SOwtM22zMzMzMza1PqPmABExHeB705GW2ZmZmZmbWhrjrOZmZmZ2SrFA2czMzMzswJ69McAp51p2zEzMzMzW6WoZCFvcTYzMzMzKzApOwcOanh4uHGm0+k0znU6nYHaWx1yk7UOnHNuVczNpNf6ZOem+7pzzrlBcqOjo41zQ0NDjXNDQ0MAzvXIDboOSnmLs5mZmZlZAQ+czczMzMwKTMpUDUm7ADsDnYhYOBltmpmZmZlNpFa2OEu6suv8W4BTgPWBoyQd0UabZmZmZmZtamuqxlpd5+cD/ysijgH2APZrqU0zMzMzs9a0NVVjDUkbkwbmioi7ASLiz5KW1YUkzScNtDn11FNb6pqZmZmZWXNtDZw3BK4mHUw6JD0lIn4naT0yB5iOiBFgZOziySef3FL3zMzMzMyaaWXgHBFzam5aDrymjTbNzMzMzNo0qT+AEhEPAL+ZzDbNzMzMzCaCj+NsZmZmZlbAA2czMzMzswKKiKnuQ51p2zEzMzMzW6XUHryi26TOcW5qeHi4cabT6TTOdTodAEZHRxvlhoaGVio3aD8n8/E1zaxsbrLW+Vhush/fZD8vM6Wfq2puZV+zq3Juuq8755ybzrmZ9Fqf7Nyg66CUp2qYmZmZmRXwwNnMzMzMrIAHzmZmZmZmBVoZOEt6vqQNqvPrSjpG0nckfVLShm20aWZmZmbWpra2OH8ZeKA6/1nST3B/srru9JbaNDMzMzNrTVtH1VgjIpZV5+dFxHOq85dIuralNs3MzMzMWtPWFueOpAOr89dJmgcgaTvg4bqQpPmSFklaNDIy0lLXzMzMzMyaa2vgfAjwEkm3ADsCl0n6NfDF6raeImIkIuZFxLz58+e31DUzMzMzs+ZamaoREfcCB0haH9i6auf2iLirjfbMzMzMzNrW6i8HRsR9wHVttmFmZmZmNhl8HGczMzMzswIeOJuZmZmZFVBETHUf6kzbjpmZmZnZKkUlC03nLc6qO0l6a+5255xbXXIzoY/OOefczMvNhD4659wE58pExIw7AYucc865mdFH55xzbublZkIfnXNusnLdp+m8xdnMzMzMbNrwwNnMzMzMrMBMHTgP+nvczjm3quVmQh+dc865mZebCX10zrnJyv3NdD6qhpmZmZnZtDFTtzibmZmZmU2uld27cLJPwMuBXwK/Ao4ozHwZ+D3QadjWlsD/B24EbgDeVZh7HHAl6efGbwCOadDmmsDPgPMa9nUJcD1wLYV7jQIbAd8AflE9xhcWZLav2hg7LQXeXdjee6rnowMsAB5XmHtXlbkh11av9QwMAecDN1d/Ny7Mvb5qbzkwr0F7J1TP52LgW8BGhbmPVplrgYXAZk3qGDicdOzzJxa2dzTw2671uGdpe8A7qtfgDcCnCtv7r662lgDXFubmApeP1TWwc2HuWcBl1WviO8AGPXI9X9+5mslksvWSyWXrJZPL1ktdrl+9ZNrL1kuuvVy9ZNrL1ksml62XTC5bL9S8pwN/B1xR1cp/AWsX5g4jfYbVvWbrcl+tnssOqe7XKsx9qbpuMek9f72SXNft/wrc36CfZwC/6VqHcwtzAo4DbqrW0TsLcxd3tXUHcE5hbnfgmip3CbBNYe6lVa4DnAnM6vHcPOazvF+tZHLZWsnksrWSyWVrpS7Xr1Yy7WVrJZPL1koml62VklOjhaf6VD0BtwBbA2tXK3fHgtyuwHNoPnDeFHhOdX79agWVtKexYgPWql4sLyhs873Af44vxoLcktyLqiZzJnBIdX5tegzyCtbH74CtCpbdvHpxrFtd/hpwQEFuuHrhzwZmAT8Cti1dz8CnqL5gAUcAnyzMPZ30JeFC6gfOvXJ7UL2JAp9s0N4GXeffCXyhtI5Jg4EfArf2qoGa9o4GDu/z3PfK/UO1DtapLj+5tJ9dt38a+EhhewuBV1Tn9wQuLMxdBbykOn8Q8NEeuZ6v71zNZDLZesnksvWSyWXrpS7Xr14y7WXrJZPL1kuun7l6ybSXrZdMLlsv1Lynk97H/qm6/gvAoYW5ZwNzqHnfzuT25NFjzi5o0F53vfwL4zY61eWqy/OA/6D3wLmuvTOAfTL1Upc7EPh3YI2aeun72QqcDexf2N5NwNOr698OnFGQ+3vgNmC76vpjgYN7PMbHfJb3q5VMLlsrmVy2VjK5bK3U5frVSqa9bK1kctlayfUzVyslp5k2VWNn4FcR8euI+CtwFrB3v1BE/AQYbdpYRNwZEddU5+8jfavZvCAXEXF/dXGt6hT9cpK2AP4ROK1pX5uStAFpwPElgIj4a0T8qeHd7A7cEhG3Fi4/C1hX0izSQPiOgszTgcsj4oGIWAZcBLym14I163lv0hcEqr+vLslFxI0R8ctcx2pyC6t+QtrytUVhbmnXxcfTo14ydfwZ4AO9Mn1yWTW5Q4HjI+KhapnfN2lPkoD/Q3ojL8kFsEF1fkN61ExNbnvgJ9X584HX9cjVvb5ra6Yu069eMrlsvWRy2Xrp895VWy8r8Z5Xl8vWS7/26uolk8vWSyaXrZfMe/pLSVvloMf7S10uIn4WEUuokcl9r7otSFtCx9dLXW4p/O35XJcV66VnTtKapP+KfKBJP+seV0HuUODYiFheLTe+XrLtSVqftE7OKcz1q5deuUeAhyLipur6Fepl/Gd59bxna6VXrupDtlYyuWytZHLZWqnL9auVulyJmly2Vvq1V1crJWbawHlz0je9MbdT8KY+ESTNIX3zu6Jw+TUlXUv6F/L5EVGSO4lUdMsH6GIACyVdLWl+wfJbA3cDp0v6maTTJD2+YZv/RI8BUM/ORfwWOBH4b+BO4N6IWFgQ7QC7SnqCpNmkb9FbNujjJhFxZ9WHO4EnN8iurIOA75cuLOk4SbcB+wEfKcy8CvhtRFw3QP8Ok7RY0pclbVyY2Q54saQrJF0k6XkN23wxcFdE3Fy4/LuBE6rn5UTgQ4W5DvCq6vzr6VMz417fRTXT9D2hIJetl/G50nrpzjWplx79LKqXcbnieql5XvrWy7hccb2My/Wtl/Hv6aT/fv6p64tPz8+jAT8LsjlJawFvAn5QmpN0Ouk/hDuQ/p1ekjsMOHfs9dCwn8dV9fIZSesU5p4G7CtpkaTvS9q2yfNC2qhywbgvlrncIcD3JN1Oej6P75cjDULXkjSvWmQfVqyX8Z/lT6CgVnrkStXmcrVSl+tXKzW5vrWS6We2VmpyfWsl0x5kaqWfmTZw7vWTiH2/3a50o9J6pE367y59kiPikYiYS/qWt7Ok4T5t7AX8PiKuHrCbL4qI5wCvAP6fpF37LD+L9O/tz0fEs4E/k/4tXUTS2qQPmq8XLr8xaUve3wGbAY+X9MZ+uYi4kfQv7PNJL/zrgGXZ0DQg6UhSP79amomIIyNiyypzWEEbs4EjKRxkj/N50hvPXNIXmU8X5mYBG5P+zfl+4GvVlolSb6Dwy1blUOA91fPyHqr/kBQ4iPQ6uJr0L/m/1i04yOt7kEwu169eeuVK6qU7V91/Ub30aK+oXnrkiuol83xm66VHrqheeuT61sv493TSf8NWWKxfrt9nQWHuc8BPIuLi0lxEHEh6770R2LcgtyvpS0SvgVO/9j5EGnQ9j7TPwAcLc+sAD0bEPOCLpLm5TZ6X2nqpyb2HNF9/C+B00tSEbA7YibTR6DOSrgTuo+szqeazvO/YZdAxQEGuZ63kcrla6ZWTtBl9aiXTXrZWMrlsrRQ8L00/ix4VDed2TOUJeCHww67LHwI+VJidQ8M5zlVuLdJ8wPeuRL+Pov980k+QvoUuIX3TewD4yoDtHV3Q3lOAJV2XXwx8t0EbewMLGyz/euBLXZf3Bz43wGP7OPD20vVM2jli0+r8psAvm9QHmTnOdTngzaQdjWYPUo/AVpnb/pYDnkHaCrKkOi0jbdF/SsP2im8jfXnZrevyLcCTCp+XWcBdwBYN1t+9PHrYTAFLB3gM2wFX1ty2wuu7X830ypTUS12uX73k2svVy/hcab0UtNfzua55LvvWS+Z5ydZLTXt966Xg8dXWS9cyR5G+CPyBR+eoP+bzKZM7vOvyEgr2TenOVefPoZrb2aS96rqX0Gcfmip3FOmzaKxelpOmSjZtb7fC9g4n7Sg7p2v93dvgeXkCcA8FO513rb9buq57KvDzAR7fHsDXui73+iz/ar9aqcl9pev2nrWSy+VqpV97dbVSk/tjv1opbG+FWqnL9auVPs9Lca30rINBQlN1Ir2R/pq01XJs58CdCrNzaL5zoEiTz09qmHsS1Y52pDlCFwN7NcivUDx9ln88sH7X+UuBlxfkLga2r84fDZzQoM2zgAMbLP980l7Js6vn9UzgHYXZJ1d/n1q9WFY4MkbdeibNuere0WuFo0Dk6oOGA2fSUV9+To/BZJ/ctl3n3wF8o2kdk/kQ7tHepl3n3wOcVZh7G2leGaQBxm1UA5V+/ayem4saPi83Ug28SHPqry7MjdXMGqTX8EE9Mj1f37maqcv0q5dMW9l6yeSy9dKvn3X1kmkvWy+ZXLZecv3M1UumvWy9ZHLZeqHmPZ30H7fuHb7eXpLLrYM+7R1Cen9ft+Z56ZV7JdXRIqrHfyJwYpN+Vtf32jmwrp+bdrV3Emmee0nu+LHnnvQ5eFVpP6taO7PB87IXaTA7tpPfwcDZhbmxelkHuAB4aU27u/HozmzZWqnL9auVTHvZWumVq9ZXtlb69bOuVjL9zNZKJpetlVw/c7VSchooNJUn0hzXm0hbLo4szCwg/XvxYdI3kBX2gK3J7UL6d8rYYZ96HrarR+6ZpMOfLCbNnVvhCAKlRVW4/NakLxFjh8wpfV7mkg7ZtJj0rbR2QDouN5v0bW3Dho/rGNLAt0Pa83adwtzFpMHFdcDuTdYz6ZvlBaRDAF0ADBXmXlOdf4i01WuFLUk1uV+RBgdj9dLr6Bi9cmdXz8ti0uGwNm9ax9R/CPdq7z9Ih95aDJxL18CoT25t0rf9DulwTCt8YNT1k7T39Nsarr9dgKurdX8F8NzC3LtI7xM3kd5gew3ue76+czWTyWTrJZPL1ksml62Xuly/esm0l62XTC5bL7l+5uol0162XjK5bL1Q855Oeu+9slqPX2fce1om986qXpaRdkg7rTC3jPTZN9b38UcbWSFH+jLw02r9dUhbP8cfbq/vZxa9B851/fxxV3tfYcXD39XlNgK+W2UvA55V2k/Sl9aeG4wy7b2mauu6Kr91Ye4E0pe0X5I/ROpuPDrQy9ZKJpetlUwuWyu9ciW1Utdev1rJ9DNbK5lctlZy/czVSsnJvxxoZmZmZlZgpu0caGZmZmY2JTxwNjMzMzMr4IGzmZmZmVkBD5zNzMzMzAp44GxmZmZmVsADZzOzKSJpC0nflnSzpF9LOqXmJ2cHvf9XS9qx6/Kxkl42UfdvZra68cDZzGwKVD89/U3gnIjYFtiW9CMLn5rAZl4N/G3gHBEfiYgfTeD9m5mtVjxwNjObGi8FHoyI0wEi4hHSr/LtL+kwSaeMLSjpPEm7Vef3kHSZpGskfV3SetX1x0v6uaTFkk6U9PfAq4ATJF0r6WmSzpC0T7X87pJ+Jul6SV8e29ItaYmkY6r7v17SDpP5pJiZTWceOJuZTY2dSL9y9zcRsZT0i36zegUkPRH4MPCyiHgO6Zc/3ytpiPQraDtFxDOBj0XEpaRf+Xt/RMyNiFu67udxpF/m2zcinlG1d2hXU3+o7v/zwOET8FjNzFYJHjibmU0NkX7+udf1dV5AmnrxU0nXAm8GtgKWAg8Cp0l6LfBAn7a3B34TETdVl88Edu26/ZvV36uBOX3uy8xsteGBs5nZ1LgBmNd9haQNgE2Ae3js+/PjxhYBzq+2IM+NiB0j4uCIWAbsDJxNmtf8gz5t5wbnAA9Vfx+hZuu3mdnqyANnM7OpcQEwW9L+AJLWBD4NnAL8BpgraQ1JW5IGxQCXAy+StE2VmS1pu2qe84YR8T3g3cDcavn7gPV7tP0LYM7Y/QBvAi6a8EdoZraK8cDZzGwKRESQ5iXvI+lm0lbm5RFxHPBT0uD5euBE4JoqczdwALBA0mLSQHoH0uD4vOq6i0g7GQKcBby/2gnwaV1tPwgcCHxd0vXAcuAL7T5iM7OZT+m928zMplJ1FIwFwGsj4up+y5uZ2eTzwNnMzMzMrICnapiZmZmZFfDA2czMzMysgAfOZmZmZmYFPHA2MzMzMyvggbOZmZmZWQEPnM3MzMzMCnjgbGZmZmZW4H8AVF/J7pcRLswAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Adapted from http://alanpryorjr.com/visualizations/seaborn/heatmap/heatmap/\n",
    "f, ax = plt.subplots(figsize=(12,6))\n",
    "ax = sns.heatmap(skills_questions_mask.T, linewidths=0.25, ax=ax, square=True, cbar=False, cmap=sns.color_palette('Greys_r'))\n",
    "ax.set(title='Skills for questions', xlabel='Question', ylabel='Skill')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"Each row corresponds to a skill and each column to a question. White squares show which skills are needed for each question\" ([see Figure 2.14](http://mbmlbook.com/LearningSkills_Moving_to_real_data.html))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Responses data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "responses = pd.read_csv('http://mbmlbook.com/Downloads/LearningSkills_Real_Data_Experiments-Original-Inputs-IsCorrect.csv', header=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "#responses.T.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs4AAAFwCAYAAACoxP20AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XuYLFdd7//3BwiBEG4dINyJAgFh0IAB5RwPIqgBRUAF4cQjCNEteBBFOVwO/oSoIAiKREUzgoDKHRQwKIeLJCLIZQMhNHJTDFej4A43ESTJ9/dH1cAwe7qmuvdePdN7v1/PM8/u6apPr9Vdq2q+Xbt6daoKSZIkScMut9sdkCRJklaBhbMkSZI0goWzJEmSNIKFsyRJkjSChbMkSZI0goWzJEmSNIKFsyRJkjSChbMkHWZJLkzyn0m+mOSiJM9Ncvxu90uSdGgsnCWpjR+qquOBU4DbAo/d5f5Ikg6RhbMkNVRVFwH/j66AJsmxSZ6W5GNJ/jXJHya5cr/sWknOSfLZJAeSvCnJ5fplFyZ5bJJ/SHJxkuckudJGO0l+Osk/9rlXJbn+pmWV5CFJPtxnfz9J+mU3S3Jeks8l+UySF2/K3TLJ6/rH/GCSH9u07Af6vnwhySeTPLL1aylJu83CWZIaSnJD4O7AP/Z3PQU4ma6QvhlwA+BX+mW/BHwCuDZwIvB/gdr0cD8OnAbctH+MX+7buAvwG8CPAdcDPgq8aEtX7gHcHvi2fr3T+vt/DXgtcE3ghsDv9o95FeB1wAuA6wD/E3hmklv3uWcDP1NVVwXWgL+Z86WRpJVj4SxJbbwiyReAjwP/Bjy+P8v708AjqupAVX0BeBJw/z7zVbrC9yZV9dWqelNVbS6cf6+qPl5VB4An0hWz0BXUf1xV76qqr9BdFnLHJCdtyj65qj5bVR8D3kh/Brxv8ybA9avqy1X1d/399wAurKrnVNUlVfUu4OXAfTblbpXkalV1cb9cko5oFs6S1Ma9+7OxdwZuCVyL7kzyccA7+8sxPgu8pr8f4Kl0Z6Zfm+QjSR6z5TE/vun2R4GNyzGu3/8OQFV9Efh3urPZGy7adPtLwMaHFR8FBHh7kvcleXB//02A79joZ9/XHweu2y//UeAHgI/2l3rcccyLIkmr7Aq73QFJOpJV1XlJngs8DfgR4D+BW1fVJ7dZ9wt0l2v8Un9JxBuTvKOq3tCvcqNNq98Y+FR/+1N0hS7wtcssTgAOamObNi+iOwtOku8CXp/kb+mK9POq6vtm5N4B3CvJMcDDgJds6Z8kHXE84yxJ7f0O8H3AtwJ/BDw9yXUAktwgyWn97Xv0H9YL8Hng0v5nw/9OcsMkE7rrnzc+yPcC4EFJTklyLN3lH2+rqgt36liS+/bXYQNcTHdN9aXAOcDJSX4iyTH9z+2TfEuSKyb58SRXr6qvbuqrJB3RLJwlqbGq+jTwJ8D/Bzya7nKMtyb5PPB64Bb9qjfvf/8i8PfAM6vq3E0P9QK6D/J9pP/59f7x39A/9suBf6H78OD9Gef2wNuSfBF4FfDzVfXP/dnv7+8f51N0l3o8BTi2z/0EcGH/HB4C/K+R7UnSyso3fu5EkrQXJbkQ+Kmqev1u90WSjlaecZYkSZJGsHCWJEmSRvBSDUmSJGkEzzhLkiRJI1g4S5IkSSPs5S9A8RoSSZIkLUPGrLSXC2cOHDgwd2Yymcydm0wmC7V3NOSWtQ3MmVuF3Nra2lyZ6XQKrMa+vuzcXt920G2/RXOrMqbNrWZulfb1ZecW3QZjeamGJEmSNIKFsyRJkjRCs0s1ktwSuBdwA7rrlT8FvKqq3t+qTUmSJKmVJmeckzwaeBHdhdZvB97R335hkse0aFOSJElqqdUZ5zOAW1fVVzffmeS3gfcBT27UriRJktREq2ucLwOuv8391+uXbSvJviT7k+xfX19v1DVJkiRpfq3OOP8C8IYkHwY+3t93Y+BmwMNmhapqHdiomGuRKUUkSZKkFpoUzlX1miQnA3eg+3BggE8A76iqS1u0KUmSJLXUbFaNqroMeGurx5ckSZKWyXmcJUmSpBEsnCVJkqQRUlW73YdZ9mzHJEmSdETJmJWaXeN8OKytrc2dmU6nzDsbx2QyATC3TW6RmU1WKbesMbbR3qo8vyP9dVnW6zmdToH5j2UbuVU4Riya2+vbDg5tTK/KPnSkvy5Hansb+96yj0mrkFt0243lpRqSJEnSCBbOkiRJ0ggWzpIkSdIIFs6SJEnSCEsvnJM8aNltSpIkSYdqN844n7kLbUqSJEmHpMl0dEkumLUIOHEgtw/YB3D22Wc36JkkSZK0mFbzOJ8InAZcvOX+AG+ZFaqqdWB949ezzjqrTe8kSZKkObUqnM8Bjq+q87cuSHJuozYlSZKkZpoUzlV1xsCy01u0KUmSJLXkdHSSJEnSCBbOkiRJ0gipqt3uwyx7tmOSJEk6omTMSq0+HHhYrK2tzZ2ZTqccOHBgrsxkMlmovel0eki5ZfdzkfbmzexWbtGxcqTnVmU7rEo/l32MWJVj0iocW5Y9NldlX9jr+9Ch5lbl9VzmvncouVU4li26DcbyUg1JkiRpBAtnSZIkaQQLZ0mSJGmEZoVzklsmuWuS47fcf7dWbUqSJEmtNCmckzwceCXwc8A0yb02LX5SizYlSZKkllrNqvHTwLdX1ReTnAS8LMlJVfUMRk73IUmSJO0lrQrny1fVFwGq6sIkd6Yrnm/CQOGcZB+wD+Dss89u1DVJkiRpfq2ucb4oySkbv/RF9D2AawG3mRWqqvWqOrWqTt23b1+jrkmSJEnza1U4PwC4aPMdVXVJVT0AuFOjNiVJkqRmmlyqUVWfGFj25hZtSpIkSS05j7MkSZI0goWzJEmSNEKqarf7MMue7ZgkSZKOKKOmS241Hd1hceDAgbkzk8lk7txkMlmovd3Kra2tzZWbTqcLt7esbbBquXm3AXTbYVWe35GeW3T7LXPfO9Jzq7LNj/R+Lmtf2Mgd6ceII71uWfYxcJH2Ft0GY3mphiRJkjSChbMkSZI0goWzJEmSNEKza5yT3AGoqnpHklsBdwM+UFV/1apNSZIkqZUmhXOSxwN3B66Q5HXAdwDnAo9JctuqemKLdiVJkqRWWp1xvg9wCnAs3Vdv37CqPp/kqcDbAAtnSZIkrZRW1zhfUlWXVtWXgH+qqs8DVNV/ApfNCiXZl2R/kv3r6+uNuiZJkiTNr9UZ5/9KclxfOH/7xp1Jrs5A4VxV68BGxVyLzMUnSZIktdCqcL5TVX0FoKo2F8rHAA9s1KYkSZLUTJPCeaNo3ub+zwCfadGmJEmS1JLzOEuSJEkjWDhLkiRJI6SqdrsPs+zZjkmSJOmIkjErNfvmwMNhkVk1JpPJ3LnJZLJQe0dDblnbwJy5Vcitra3NlZlOp8Bq7OvLzu31bQfd9lt2blX2BXO7m9utfX0VjoGL7ntjeamGJEmSNIKFsyRJkjSChbMkSZI0wtIK5yR/sqy2JEmSpMOtyYcDk7xq613A9yS5BkBV3bNFu5IkSVIrrWbVuCHwD8Cz6KaVC3Aq8FuN2pMkSZKaanWpxqnAO4HHAZ+rqnOB/6yq86rqvFmhJPuS7E+yf319vVHXJEmSpPk1OeNcVZcBT0/y0v7ffx3TVlWtAxsVcy0yH6IkSZLUQtMvQKmqTwD3TfKDwOdbtiVJkiS1tJRvDqyqVwOvXkZbkiRJUgvO4yxJkiSNYOEsSZIkjZCq2u0+zLJnOyZJkqQjSsastJRrnBe1trY2d2Y6nc6dm06nAMw7i8dkMgHm7+dGe6uQW9Y2ONTcIjOwTCaTlXl+5lYzd6j77JF8TPIYYW6enNvv4Aysxr6+7Nyi22AsL9WQJEmSRrBwliRJkkawcJYkSZJGWMo1zkm+C7gDMK2q1y6jTUmSJOlwanLGOcnbN93+aeD3gKsCj0/ymBZtSpIkSS21ulTjmE239wHfV1VnAt8P/HijNiVJkqRmWl2qcbkk16QrzFNVnwaoqv9IcsmsUJJ9dIU2Z599dqOuSZIkSfNrVThfHXgn3WTSleS6VXVRkuMZmGC6qtaB9Y1fzzrrrEbdkyRJkubTpHCuqpNmLLoM+OEWbUqSJEktLfWbA6vqS8A/L7NNSZIk6XBwHmdJkiRpBAtnSZIkaYRU1W73YZY92zFJkiQdUWZOXrHZUq9xnteBAwfmzkwmk7lzk8lkofaOhtyytsHRkltbW5s7N51OV+b5Hem5ebffdDoFVmNfX3Zur2876Lbfkb7PHunP70jNrdK+vuzcottgLC/VkCRJkkawcJYkSZJGsHCWJEmSRmhSOCf5jiRX629fOcmZSf4yyVOSXL1Fm5IkSVJLrc44/zHwpf72M+i+gvsp/X3PadSmJEmS1EyrWTUuV1WX9LdPrarb9bf/Lsn5jdqUJEmSmml1xnma5EH97fckORUgycnAV2eFkuxLsj/J/vX19UZdkyRJkubX6ozzTwHPSPLLwGeAv0/yceDj/bJtVdU6sFEx1yJz8UmSJEktNCmcq+pzwE8muSrwzX07n6iqf23RniRJktRa028OrKovAO9p2YYkSZK0DM7jLEmSJI1g4SxJkiSNkKra7T7Msmc7JkmSpCNKxqzU9BrnQ7XIrBqTyWTu3GQyAWBtbW2u3HQ6PaTcov1cZm5Z28CcuSMxtxv77Krklr3t5j1OQ3esXnZur4/pQ82tyuu517ffbu3rq1DvLLoNxvJSDUmSJGkEC2dJkiRpBAtnSZIkaYQmhXOShye5UYvHliRJknZDqzPOvwa8Lcmbkvxskms3akeSJElailaF80eAG9IV0N8O/EOS1yR5YP813JIkSdJKaVU4V1VdVlWvraozgOsDzwTuRldUbyvJviT7k+xfX19v1DVJkiRpfq3mcf6GSaSr6qvAq4BXJbnyrFBVrQMbFXMtMhefJEmS1EKrM873m7Wgqv6zUZuSJElSM00K56r6UIvHlSRJknaL8zhLkiRJI1g4S5IkSSOkqna7D7Ps2Y5JkiTpiJKdV2k3q8ZhscisGpPJZO7cZDIBYG1tba7cdDoF5u/nRnurkFvWNjBnbhVyyz5GHMnt7fVtB93zW5WxaW5v5JY1zlapjlh2btFtN5aXakiSJEkjWDhLkiRJI1g4S5IkSSM0ucY5yRWB+wOfqqrXJzkd+G/A+4H1/psEJUmSpJXR6sOBz+kf+7gkDwSOB/4cuCtwB+CBjdqVJEmSmmhVON+mqr41yRWATwLXr6pLk/wZ8J5ZoST7gH0AZ599Nve5z30adU+SJEmaT6vC+XL95RpXAY4Drg4cAI4FjpkVqqp1YH3j10WmFJEkSZJaaFU4Pxv4AHB54HHAS5N8BPhO4EWN2pQkSZKaaVI4V9XTk7y4v/2pJH8CfC/wR1X19hZtSpIkSS01++bAqvrUptufBV7Wqi1JkiSpNedxliRJkkawcJYkSZJGSFXtdh9m2bMdkyRJ0hElY1Zqdo3z4bDIdHSTyWTu3GQyWai9oyG3rG2warm1tbW5c9PpdGWe36rkFt0Oy8pNp1NgNfb1ZedWZaysyr5g7ujKrdK+Dix87Fwkt+g2GMtLNSRJkqQRLJwlSZKkESycJUmSpBGaXeOc5KbADwM3Ai4BPgy8sKo+16pNSZIkqZUmZ5yTPBz4Q+BKwO2BK9MV0H+f5M4t2pQkSZJaanWpxk8Dd6uqX6f7qu1bVdXjgLsBT58VSrIvyf4k+9fX1xt1TZIkSZpfy+norgBcChwLXBWgqj6W5JhZgapaBzYq5lpkShFJkiSphVaF87OAdyR5K3An4CkASa4NWA1LkiRp5TQpnKvqGUleD3wL8NtV9YH+/k/TFdKSJEnSSml2qUZVvQ94X6vHlyRJkpbJeZwlSZKkESycJUmSpBFSVTuvlJwIPAm4flXdPcmtgDtW1bMb9m3njkmSJEmHLqNWGlk4/zXwHOBxVfVtSa4AvLuqbnNofRxUa2trc4em0ynzTmM3mUwAViY37+synU4Xbm+RKQF3I7essbLRnrmjM7fovrfMfXZVcnt924HHiKHckf567vXnt0p1xLLbW3TbMbJwHnupxrWq6iXAZQBVdQndHM2SJEnSUWFs4fwfSU6gv3wiyXcCn2vWK0mSJGmPGTsd3S8CrwJumuTNwLWB+zTrlSRJkrTHjCqcq+pdSb4buAXdNSAfrKqvNu2ZJEmStIcMFs5JfmTGopOTUFV/3qBPkiRJ0p6z0xnnHxpYVsC2hXOSqwOPBe5Nd1kHwL8BrwSeXFWfnZHbB+wDOPvss3fomiRJkrQ8g4VzVT1owcd9CfA3wJ2r6iKAJNcFHgi8FPi+Ge2tA+sbv5511lkLNi9JkiQdXqNm1UhyQpKzkrwryTuTPKOfZWOWk6rqKRtFM0BVXVRVTwFufKidliRJkpZt7HR0LwI+Dfwo3WwanwZePLD+R5M8qv/GQaD79sEkjwY+vmhnJUmSpN0ytnCeVNWvVdU/9z+/DlxjYP37AScA5yU5kOQAcC4wAe57SD2WJEmSdsHYwvmNSe6f5HL9z48Br561clVdXFWPrqpbVtWk//mWqno03QcGJUmSpJUytnD+GeAFwH/1Py8CfjHJF5J8fs42z5xzfUmSJGnXpaoO/4MmF8xaBJxcVceOeJjD3zFJkiTpYBmz0tiv3CbJPYE79b+eW1XnDKx+InAacPE2nXrL2DbX1tbGrvo10+l07tx0Ol2ovY3cgQMH5spNJpOVyc2bWbXcssbYbuWW/bqsynZf1utyNBwjFs3t9W0Hq7Wvr8qxxdz2uWXXLatyLFukn4uO6bFGFc5JngzcHnh+f9fPJ/muqnrMjMg5wPFVdf42j3Xu6N5JkiRJe8TYM84/AJxSVZcBJHke8G5g28K5qs6Y9UBVdfq8nZQkSZJ229gPB8I3Tj939cPdEUmSJGkvG3vG+TeAdyd5I911yncCHtusV5IkSdIeM6pwrqoX9tcm356ucH705q/TliRJko50oy7VSBLgrnTXOb8SuGKSOyzSYJK/Hli2L8n+JPvX19cXeXhJkiSpibGXajwTuAy4C/CrwBeAl9OdgT5IktvNeJwAp8xqpKrWgY2Kuc4666yR3ZMkSZLaGls4f0dV3S7Ju6H7Su0kVxxY/x3AeWw/mfQ1trlPkiRJ2tPGFs5fTXJ5+m/zS3JtujPQs7wf+Jmq+vDWBUk+PncvJUmSpF02djq6s4C/AK6T5InA3wFPGlj/CQOP/XOjeydJkiTtEWNn1Xh+knfSfUAwwL2r6v0D679s4OGuOV8XJUmSpN03WDgnuRLwEOBmwHuBs6vqkkNs80zgOYf4GJIkSdJSpapmL0xeDHwVeBNwd+DCqvqFHR80uWDWIuDkqjp2RN9md0ySJEk6fLab0OIgO12qcauqug1AkmcDbx/Z+InAacDF23TqLSMfgwMHDoxd9Wsmk8ncuclkslB7u5VbW1ubKzedThdub1nb4GjJzbvtoNt+q/L8jvTcMve9Iz2317cddNvPfXb73JH+uuz157dK+/qh5BY55i66DcbaqXD+6saNqrqk+x6UUc4Bjq+q87cu6L+BUJIkSVopOxXO35bk8/3tAFfufw9QVXW17UJVdcasB6yq0xfqqSRJkrSLBgvnqrr8sjoiSZIk7WVj53GWJEmSjmoWzpIkSdIITQrnJFdL8htJ/jTJ6VuWPXMgty/J/iT719fXW3RNkiRJWsiobw5cwHOADwMvBx6c5EeB06vqK8B3zgpV1TqwUTHXIlOKSJIkSS20ulTjplX1mKp6RVXdE3gX8DdJTmjUniRJktRUqzPOxya5XFVdBlBVT0zyCeBvgeMbtSlJkiQ10+qM818Cd9l8R1U9D/gl4L8atSlJkiQ10+SMc1U9asb9r0nypBZtSpIkSS3txnR0Z+5Cm5IkSdIhSVUd/gdNLpi1CDi5qo4d8TCHv2OSJEnSwTJmpVYfDjwROA24eMv9Ad4y9kEWmY5uMpmwtrY2V2Y6nS7U3mQyOaTcov1c5vNbdBsc6bl5twF022HR3Kq8LquSW9Z22K1jxLL7uQrHlmXve+7r2+dWZTvs9e23G/ssLP+YtEh7i267sVoVzucAx1fV+VsXJDm3UZuSJElSM60+HHjGwLLTZy2TJEmS9qrd+HCgJEmStHIsnCVJkqQRLJwlSZKkEZoUzkmum+QPkvx+khOSPCHJe5O8JMn1BnL7kuxPsn99fb1F1yRJkqSFtJpV47nAq4GrAG8Eng/8IHAv4A/7fw9SVevARsVci0zrIkmSJLXQ6lKNE6vqd6vqycA1quopVfWxqvpd4CaN2pQkSZKaaVU4b37cP9my7PKN2pQkSZKaaVU4vzLJ8QBV9csbdya5GfDBRm1KkiRJzbT6ApRfmXH/PyZ5dYs2JUmSpJZ2Yzq6M3ehTUmSJOmQpKoO/4MmF8xaBJxcVceOeJjD3zFJkiTpYBmzUqvp6E4ETgMu3nJ/gLeMfZBFpqObTCZz5yaTyULt7VZubW1trtx0Ol24vWVtg0PNzfuaQPe6LJpbdj9Xpb1V2e7Lau9Q9j1Y7r6+7NyqjBX3oYOt0vM7Ul+XVatblnkMXHQbjNWqcD4HOL6qzt+6IMm5jdqUJEmSmmn14cAzBpad3qJNSZIkqaXd+HCgJEmStHIsnCVJkqQRllY4J7nOstqSJEmSDrcmhXOSyZafE4C3J7lmkpkfXUyyL8n+JPvX19dbdE2SJElaSKtZNT4DfHTLfTcA3kU3P/M3bxeqqnVgo2KuRaYUkSRJklpodanGo4APAvesqm+qqm8CPtHf3rZoliRJkvayJoVzVT0N+CngV5L8dpKr4jcBSpIkaYU1+3BgVX2iqu4LvBF4HXBcq7YkSZKk1prPqlFVfwl8D/C9AEke1LpNSZIk6XBbynR0VfWfVTXtfz1zGW1KkiRJh1OqDv+lx0kumLUIOLmqjh3xMF4TLUmSpGXImJVaTUd3InAacPGW+wO8ZeyDLDId3WQymTs3mUwWam8jt7a2NlduOp0eUnvLzC1rG5gztwq5Rff1I/kYsWhuVbb5kZ6bd2xCNz5X5fkdqa/LKu3ry84tug3GalU4nwMcX1Xnb12Q5NxGbUqSJEnNNCmcq+qMgWWnt2hTkiRJamkpHw6UJEmSVp2FsyRJkjSChbMkSZI0QpPCOcndNt2+epJnJ7kgyQuSnDiQ25dkf5L96+vrLbomSZIkLaTVrBpPAl7T3/4t4F+AHwJ+BDgbuPd2oapaBzYq5lpkShFJkiSphVaF82anVtUp/e2nJ3ngEtqUJEmSDqtWhfN1kvwi3ReeXC1J6utfUeh11ZIkSVo5rYrYPwKuChwPPA+4FkCS6wIHfSmKJEmStNe1+gKUM2fcf1GSN7ZoU5IkSWopX7+CYkkNJh+rqhuPWHW5HZMkSdLRKmNWanLGOckFsxYBM6ejkyRJkvaqVh8OPBE4Dbh4y/0B3jL2QRaZjm4ymcydm0wmC7V3NOSWtQ3MHd25tbW1uXPT6XTP93M6nQKrsa8vO7fXtx10229Vxqa5w5vb69t9lfb1Q8ktcsxddBuM1apwPgc4vqoO+iBgknMbtSlJkiQ10+rDgWcMLDu9RZuSJElSS86pLEmSJI1g4SxJkiSNsLTCOckJy2pLkiRJOtyaFM5Jnpxk49sCT03yEeBtST6a5LsHcvuS7E+yf319vUXXJEmSpIW0mlXjB6vqMf3tpwL3q6p3JDkZeAFw6nahqloHNirmWmRKEUmSJKmFVpdqHJNkoyi/clW9A6CqPgQc26hNSZIkqZlWhfPvA3+V5C7Aa5L8TpI7JTkTOGhuZ0mSJGmvazWP8+8meS/wUODkvp2TgVcAv9aiTUmSJKmlVtc4U1XnAuduvT/Jg4DntGpXkiRJaiFVtdwGk49V1Y1HrLrcjkmSJOlolTErNTnjnOSCWYuAE1u0KUmSJLXU6lKNE4HTgIu33B/gLWMfZJHp6CaTydy5yWSyUHu7lVtbW5srN51OF25vWdvAnLlVyC1z3zvSc3t920G3/RbNrcqYXpXtsCrPb1XqllU5li3Sz0XH2FitCudzgOOr6qAZNJKc26hNSZIkqZlWs2qcMbDs9BZtSpIkSS21msdZkiRJOqJYOEuSJEkjWDhLkiRJIzQpnJO8K8kvJ7npnLl9SfYn2b++vt6ia5IkSdJCWs2qcU3gGsAbk1wEvBB4cVV9aihUVevARsVci0zrIkmSJLXQ6lKNi6vqkf03BP4ScHPgXUnemGRfozYlSZKkZppf41xVb6qqnwVuADwFuGPrNiVJkqTDrdWlGh/aekdVXQq8pv+RJEmSVkqTM85Vdf9Zy5I8qEWbkiRJUkupquU2mHysv/Z5J8vtmCRJko5WGbNSk0s1klwwaxFwYos2JUmSpJZaXeN8InAacPGW+wO8ZeyDLDId3WQymTs3mUwAWFtbmys3nU6B+fu50d4q9HPetjbaW9a2O9Tckf78zB3e3KocI1Yht9e3HRwd+7qvy/a5vf667Na+vugxcJm5RbfBWK0K53OA46vq/K0LkpzbqE1JkiSpmSaFc1WdMbDs9BZtSpIkSS01n8dZkiRJOhJYOEuSJEkjWDhLkiRJIzQpnJOcmuSNSf4syY2SvC7J55K8I8ltB3L7kuxPsn99fb1F1yRJkqSFtJpV45nA44Fr0E0/94iq+r4kd+2X3XG7UFWtAxsVcy0ypYgkSZLUQqtLNY6pqr+uqhcCVVUvo7vxBuBKjdqUJEmSmmlVOH85yfcnuS9QSe4NkOS7gUsbtSlJkiQ10+pSjYcAvwlcRvcNgg9N8lzgk8BPN2pTkiRJaqbJGeeqek9VnVZVd6+qD1TVz1fVNarq1sAtWrQpSZIktZSqWm6Dyceq6sYjVl1uxyRJknS0ypiVmlyqkeSCWYuAE1u0KUmSJLXU6hrnE+mubb54y/2hm55ulEWmo5tMJnPnJpPJQu0dDbllbQNz5o7E3Crt68vOLXvbra2tzZ2bTqd7fowdas7XZTVzq7SvLzu36DYYq1XhfA5wfFWdv3VBknMbtSlJkiQ106RwrqozBpad3qJNSZIkqaVW8zhLkiRJRxQLZ0mSJGkEC2dJkiRphCaFc5Ljk/xqkvcl+VySTyd5a5Jjl1tqAAAT+klEQVSf3CG3L8n+JPvX19dbdE2SJElaSKtZNZ4P/AXdlHQ/BlwFeBHwy0lOrqr/u12oqtaBjYq5FplSRJIkSWqh1aUaJ1XVc6vqE1X128A9q+rDwIOAH2nUpiRJktRMq8L5P5J8F0CSHwIOAFTVZYz8SkNJkiRpL2l1qcZDgGclORmYAg8GSHJt4PcbtSlJkiQ10+oLUC4A7rDN/Z9O8oUWbUqSJEktpaqW22Dysaq68YhVl9sxSZIkHa1GXUrc5IxzkgtmLQJObNGmJEmS1FKra5xPpJuK7uIt9wd4y9gHWWQ6uslkMnduMpks1N6q5dbW1ubKTafTpW0Dc+ZWIbfIPgSrc4xYZm6vbzvwGGhu7+ZWaV9fdm7RbTBWq8L5HOD4qjp/64Ik5zZqU5IkSWqm1YcDzxhYdnqLNiVJkqSWWs3jLEmSJB1RLJwlSZKkEZoUzkmunuTJST6Q5N/7n/f3912jRZuSJElSS63OOL+EbkaNO1fVCVV1AvA9/X0vnRVKsi/J/iT719fXG3VNkiRJml+rWTVOqqqnbL6jqi4CnpLkwbNCVbUObFTMtciUIpIkSVILrc44fzTJo5J87ctOkpyY5NHAxxu1KUmSJDXTqnC+H3ACcF6Si5McAM4FJsCPNWpTkiRJaqbVPM4XJ3kO8DrgrVX1xY1lSe4GvKZFu5IkSVIrrWbVeDjwSuBhwDTJvTYtflKLNiVJkqSWUlWH/0GT9wJ3rKovJjkJeBnwp1X1jCTvrqrbjniYw98xSZIk6WAZs1KrWTUuv3F5RlVdmOTOwMuS3GRsxwAWmVVjMpnMnZtMJgCsra3NlZtOp8D8/dxobxVyy9oGR0tu3jEG3Tjz+W3f3l5/fhvHCI8tB+f2+raD1Rqbvi5HV26V9vVl5xbdBmO1+nDgRUlO2filL6LvAVwLuE2jNiVJkqRmWhXODwAu2nxHVV1SVQ8A7tSoTUmSJKmZVrNqfGJg2ZtbtClJkiS11OqMsyRJknREsXCWJEmSRlh64Zzkr5fdpiRJknSomlzjnOR2sxYBp8xYRpJ9wD6As88+m/vc5z4NeidJkiTNr9U8zu8AzmP7OZuvMStUVevA+savi8zFJ0mSJLXQqnB+P/AzVfXhrQuSfLxRm5IkSVIzra5xfsLAY/9cozYlSZKkZpoUzlX1MiBJ7prk+C2Lv9yiTUmSJKmlJoVzkocDr6Q7uzxNcq9Ni5/Uok1JkiSppVTV4X/Q5L3AHavqi0lOAl4G/GlVPSPJu6vqtiMe5vB3TJIkSTrYdhNaHKTVNc6Xr6ovAlTVhcCdgbsn+e2xHevX2/Ynyc8MLTdn7mjJrUIfzZkzt3q5VeijOXOHOTdOVR32H+BvgFO23HcF4E+ASw/D4+83Z87cavTRnDlzq5dbhT6aM7es3OafVmecHwBctPmOqrqkqh4A3KlRm5IkSVIzTeZxrqpPDCx7c4s2JUmSpJZanXFubX3nVcyZOypyq9BHc+bMrV5uFfpoztyycl/TZFYNSZIk6UizqmecJUmSpOU61E8XLvsHuBvwQeAfgceMzPwx8G/AdM62bgS8EXg/8D7g50fmrgS8HXhPnztzjjYvD7wbOGfOvl4IvBc4n5GfGgWuQTfH9gf653jHEZlb9G1s/Hwe+IWR7T2ifz2mwAuBK43M/Xyfed9QW9ttZ2ACvA74cP/vNUfm7tu3dxlw6hztPbV/PS8A/gK4xsjcr/WZ84HXAtefZxwDj6Sb+/xaI9t7AvDJTdvxB8a2R/fFRh/sX5/fHNneize1dSFw/sjcKcBbN8Y1cIeRuW8D/r7fJ/4SuNo2uW3376ExM5AZHC8DucHxMpAbHC+zcjuNl4H2BsfLUHtD42WgvcHxMpAbHC8DucHxwoxjOvBNwNv6sfJi4Iojcw+j+xs2a5+dlXt+/1pO6cb9MSNzz+7vu4DumH/8mNym5b8LfHGOfj4X+OdN23DrLFuzcgGeCHyo30YPH5l706a2PgW8YmTursC7+tzfATcbmbtLn5sCzwOusM1r8w1/y3caKwO5wbEykBscKwO5wbEyK7fTWBlob3CsDOQGx8pAbnCsjPmZa+Xd/ulfgH8Cvhm4Yr9xbzUidyfgdsxfOF8PuF1/+6r9BhrTXjYGG3BMv7N858g2fxF4wdbBOCJ34dBONSPzPOCn+ttXZJsib8T2uAi4yYh1b9DvHFfuf38J8JMjcmv9jn8c3YdZXw/cfOx2Bn6T/g0W8BjgKSNz30L3JuFcZhfO2+W+n/4gCjxljvautun2w4E/HDuO6YqB/wd8dLsxMKO9JwCP3OG13y73Pf02OLb//Tpj+7lp+W8BvzKyvdcCd+9v/wBw7sjcO4Dv7m8/GPi1bXLb7t9DY2YgMzheBnKD42UgNzheZuV2Gi8D7Q2Ol4Hc4HgZ6ufQeBlob3C8DOQGxwszjul0x7H79/f/IfDQkbnbAicx47g9kPsBvj7n7AvnaG/zePlttpx0mpXrfz8V+FO2L5xntfdc4D4D42VW7kF009ZebsZ42fFvK/By4AEj2/sQ8C39/T8LPHdE7r8BHwdO7u//VeCMbZ7jN/wt32msDOQGx8pAbnCsDOQGx8qs3E5jZaC9wbEykBscK0P9HBorY35W7VKNOwD/WFUfqar/Al4E3GuHDFX1t8CBeRurqn+pqnf1t79A967mBiNyVf0XwNDtbMcw4psQk9wQ+EHgWfP2dV5JrkZXcDwboKr+q6o+O+fD3BX4p6r66Mj1rwBcOckV6ArhT43IfAvw1qr6UlVdApwH/PB2K87Yzveie4NA/++9x+Sq6v1V9cGhjs3IvbbvJ3Rnvm44Mvf5Tb9ehW3Gy8A4fjrwqO0yO+QGzcg9FHhyVX2lX+ff5mkvSYAfozuQj8kVcLX+9tXZZszMyN0C+Nv+9uuAH90mN2v/njlmZmV2Gi8DucHxMpAbHC87HLtmjpdDOObNyg2Ol53amzVeBnKD42UgNzheBo7pd6E7KwfbHF9m5arq3dV9Odi2BnJ/1S8rujOhW8fLrNzn4Wuv55U5eLxsm0tyebr/FXnUPP2c9bxG5B4K/GpVXdavt3W8DLaX5Kp02+QVI3M7jZftcpcCX6mqD/X3HzRetv4t71/3wbGyXa7vw+BYGcgNjpWB3OBYmZXbaazMyo0xIzc4VnZqb9ZYGWPVCucb0L3T2/AJRhzUD4f+q8NvS/eOc8z6l09yPt1/Ib+uqsbkfodu0F22QBcLeG2SdybZN2L9bwY+DTwnybuTPCvJVeZs8/5sUwBt27mqTwJPAz4G/Avwuap67YjoFLhTkhOSHEf3LvpGc/TxxKr6l74P/wJcZ47soXow8NdjV07yxCQfB34c+JWRmXsCn6yq9yzQv4cluSDJHye55sjMycD/SPK2JOcluf2cbf4P4F+r6sMj1/8F4Kn96/I04LEjc1Pgnv3t+7LDmNmyf48aM/MeE0bkBsfL1tzY8bI5N8942aafo8bLltzo8TLjddlxvGzJjR4vW3I7jpetx3S6//387KY3Ptv+PVrwb8FgLskxwE8ArxmbS/Icuv8hvCXdf6ePyT0MeNXG/jBnP5/Yj5enJzl2ZO6mwP2S7E/y10luPs/rQndS5Q1b3lgO5X4K+Kskn6B7PZ+8U46uCD0myan9Kvfh4PGy9W/5CYwYK9vkxpqZGxors3I7jZUZuR3HykA/B8fKjNyOY2WgPRgYKztZtcJ5u69E3PHd7SE3mhxPd0r/F8a+yFV1aVWdQvcu7w5J1nZo4x7Av1XVOxfs5n+vqtsBdwf+d5KdvmjmCnT/vf0HVXVb4D/o/lt6lCRXpPtD89KR61+T7kzeNwHXB66S5H/tlKuq99P9F/br6Hb89wCXDIb2gCSPo+vn88dmqupxVXWjPvOwEW0cBzyOkUX2Fn9Ad+A5he6NzG+NzF0BuCbdf3P+H+Al/ZmJsf4nI99s9R4KPKJ/XR5B/z8kIzyYbj94J91/yf/XrBUX2b8XyQzldhov2+XGjJfNuf7xR42XbdobNV62yY0aLwOv5+B42SY3arxsk9txvGw9ptP9b9hBq+2U2+lvwcjcM4G/rao3jc1V1YPojr3vB+43IncnujcR2xVOO7X3WLqi6/Z0nxl49MjcscCXq+pU4I/ors2d53WZOV5m5B5Bd73+DYHn0F2aMJgDbk130ujpSd4OfIFNf5Nm/C3fsXZZtAYYkdt2rAzlhsbKdrkk12eHsTLQ3uBYGcgNjpURr8u8f4u+rua8tmM3f4A7Av9v0++PBR47MnsSc17j3OeOobse8BcPod+PZ+frSX+D7l3ohXTv9L4E/NmC7T1hRHvXBS7c9Pv/AF49Rxv3Al47x/r3BZ696fcHAM9c4Lk9CfjZsduZ7sMR1+tvXw/44Dzjg4FrnGflgAfSfdDouEXGI3CTgWVfywG3oTsLcmH/cwndGf3rztne6GV0b17uvOn3fwKuPfJ1uQLwr8AN59h+n+Pr02YG+PwCz+Fk4O0zlh20f+80ZrbLjBkvs3I7jZeh9obGy9bc2PEyor1tX+sZr+WO42XgdRkcLzPa23G8jHh+M8fLpnUeT/dG4DN8/Rr1b/j7NJB75KbfL2TEZ1M25/rbr6C/tnOe9vr7vpsdPkPT5x5P97doY7xcRnep5Lzt3Xlke4+k+6DsSZu23+fmeF1OAP6dER8637T9/mnTfTcG/mGB5/f9wEs2/b7d3/Ln7zRWZuT+bNPybcfKUG5orOzU3qyxMiN38U5jZWR7B42VWbmdxsoOr8vosbLtOFgktFs/dAfSj9Cdtdz4cOCtR2ZPYv4PB4bu4vPfmTN3bfoP2tFdI/Qm4B5z5A8aPDusfxXgqptuvwW424jcm4Bb9LefADx1jjZfBDxojvW/g+5Tycf1r+vzgJ8bmb1O/++N+53loJkxZm1numuuNn/Q66BZIIbGB3MWznSzvvwD2xSTO+Ruvun2zwEvm3ccM/BHeJv2rrfp9iOAF43MPYTuujLoCoyP0xcqO/Wzf23Om/N1eT994UV3Tf07R+Y2xszl6PbhB2+T2Xb/HhozszI7jZeBtgbHy0BucLzs1M9Z42WgvcHxMpAbHC9D/RwaLwPtDY6XgdzgeGHGMZ3uf9w2f+DrZ8fkhrbBDu39FN3x/cozXpftcj9EP1tE//yfBjxtnn7292/34cBZ/bzepvZ+h+469zG5J2+89nR/B98xtp/9WHveHK/LPeiK2Y0P+Z0BvHxkbmO8HAu8AbjLjHbvzNc/zDY4VmbldhorA+0NjpXtcv32GhwrO/Vz1lgZ6OfgWBnIDY6VoX4OjZUxPwuFdvOH7hrXD9GduXjcyMwL6f578at070AO+gTsjNx30f13ysa0T9tO27VN7lvppj+5gO7auYNmEBg7qEau/810byI2pswZ+7qcQjdl0wV070pnFqRbcsfRvVu7+pzP60y6wndK98nbY0fm3kRXXLwHuOs825nuneUb6KYAegMwGZn74f72V+jOeh10JmlG7h/pioON8bLd7Bjb5V7evy4X0E2HdYN5xzGz/whv196f0k29dQHwKjYVRjvkrkj3bn9KNx3TQX8wZvWT7tPTD5lz+30X8M5+278N+PaRuZ+nO058iO4Au11xv+3+PTRmBjKD42UgNzheBnKD42VWbqfxMtDe4HgZyA2Ol6F+Do2XgfYGx8tAbnC8MOOYTnfsfXu/HV/KlmPaQO7h/Xi5hO4Dac8ambuE7m/fRt+3zjZyUI7uzcCb++03pTv7uXW6vR3/ZrF94Tyrn3+zqb0/4+Dp72blrgG8us/+PfBtY/tJ96Z12xNGA+39cN/We/r8N4/MPZXuTdoHGZ4i9c58vdAbHCsDucGxMpAbHCvb5caMlVnt7TRWBvo5OFYGcoNjZaifQ2NlzI/fHChJkiSNsGofDpQkSZJ2hYWzJEmSNIKFsyRJkjSChbMkSZI0goWzJEmSNIKFsyTtkiQ3TPLKJB9O8pEkvzfjK2cXffx7J7nVpt9/Ncn3Hq7Hl6SjjYWzJO2C/qun/xx4RVXdHLg53Zcs/OZhbObewNcK56r6lap6/WF8fEk6qlg4S9LuuAvw5ap6DkBVXUr3rXwPSPKwJL+3sWKSc5Lcub/9/Un+Psm7krw0yfH9/U9O8g9JLkjytCT/Dbgn8NQk5ye5aZLnJrlPv/5dk7w7yXuT/PHGme4kFyY5s3/89ya55TJfFEnayyycJWl33JruW+6+pqo+T/eNflfYLpDkWsAvA99bVbej++bPX0wyofsWtFtX1bcCv15Vb6H7lr//U1WnVNU/bXqcK9F9M9/9quo2fXsP3dTUZ/rH/wPgkYfhuUrSEcHCWZJ2R+i+/nm7+2f5TrpLL96c5HzggcBNgM8DXwaeleRHgC/t0PYtgH+uqg/1vz8PuNOm5X/e//tO4KQdHkuSjhoWzpK0O94HnLr5jiRXA04E/p1vPD5faWMV4HX9GeRTqupWVXVGVV0C3AF4Od11za/Zoe2h4hzgK/2/lzLj7LckHY0snCVpd7wBOC7JAwCSXB74LeD3gH8GTklyuSQ3oiuKAd4K/PckN+szxyU5ub/O+epV9VfALwCn9Ot/AbjqNm1/ADhp43GAnwDOO+zPUJKOMBbOkrQLqqrorku+T5IP051lvqyqngi8ma54fi/wNOBdfebTwE8CL0xyAV0hfUu64vic/r7z6D5kCPAi4P/0HwK86aa2vww8CHhpkvcClwF/2PYZS9LqS3fsliTtpn4WjBcCP1JV79xpfUnS8lk4S5IkSSN4qYYkSZI0goWzJEmSNIKFsyRJkjSChbMkSZI0goWzJEmSNIKFsyRJkjSChbMkSZI0wv8PTE0Ww12xmQ8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f, ax = plt.subplots(figsize=(12,6))\n",
    "ax = sns.heatmap(responses.T, linewidths=0.25, ax=ax, square=True, cbar=False, cmap=sns.color_palette('Greys_r'))\n",
    "ax.set(title='Responses', xlabel='Question', ylabel='People')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"Each row corresponds to a person and again each column corresponds to a question. Here, white squares show which questions each person got correct.\" ([see Figure 2.14](http://mbmlbook.com/LearningSkills_Moving_to_real_data.html))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "skills_needed = torch.tensor(skills_questions_mask.T.values.astype(np.float32))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "people_n = 22\n",
    "skills_n = 7\n",
    "questions_n = 48\n",
    "\n",
    "guess_probs = 0.2\n",
    "mistake_probs = 0.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "metadata": {},
   "outputs": [],
   "source": [
    "def basic_model():\n",
    "    \n",
    "    # declare dimensions as conditionally independent\n",
    "    with pyro.plate('people', people_n):\n",
    "        with pyro.plate('skills', skills_n):\n",
    "            skills = pyro.sample('skills', pyro.distributions.Bernoulli(0.5))\n",
    "    \n",
    "    # setup empty arrays\n",
    "    has_skills = torch.zeros((people_n, questions_n))\n",
    "    is_correct_probs = torch.zeros((people_n, questions_n))\n",
    "    is_correct = torch.zeros((people_n, questions_n))\n",
    "    \n",
    "    # declare dimensions as conditionally independent\n",
    "    for p in pyro.plate('people', people_n):\n",
    "        for q in pyro.plate('questions', questions_n):\n",
    "            \n",
    "            # skills_needed AND skills\n",
    "            has_skills[p,q] = torch.tensor(torch.allclose(skills_needed[:,q], skills[:,p])).float()\n",
    "            \n",
    "            # probabilites\n",
    "            is_correct_probs[p,q] = has_skills[p,q] * (1 - mistake_probs) + (1 - has_skills[p,q]) * guess_probs\n",
    "            \n",
    "            # sample is correct\n",
    "            is_correct[p,q] = pyro.sample('is_correct_{}'.format(q), pyro.distributions.Bernoulli(is_correct_probs[p,q]))\n",
    "    \n",
    "    return skills, has_skills, is_correct"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 155,
   "metadata": {},
   "outputs": [],
   "source": [
    "skills, has_skills, is_correct = basic_model()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(torch.Size([7, 22]), torch.Size([22, 48]), torch.Size([22, 48]))"
      ]
     },
     "execution_count": 115,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "skills.shape, has_skills.shape, is_correct.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "metadata": {},
   "outputs": [],
   "source": [
    "is_correct_avg_pp = torch.zeros(22,48)\n",
    "n = 100\n",
    "\n",
    "for i in range(n):\n",
    "    _, _, is_correct = basic_model()\n",
    "    is_correct_avg_pp += is_correct\n",
    "    \n",
    "is_correct_avg_pp /= n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 138,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([[0.2700, 0.2100, 0.1700,  ..., 0.2200, 0.2000, 0.2200],\n",
       "        [0.3000, 0.2100, 0.1700,  ..., 0.1900, 0.1900, 0.3300],\n",
       "        [0.1600, 0.1900, 0.1400,  ..., 0.3100, 0.2200, 0.2300],\n",
       "        ...,\n",
       "        [0.1600, 0.2100, 0.2400,  ..., 0.2400, 0.2100, 0.1600],\n",
       "        [0.2000, 0.1600, 0.1500,  ..., 0.2600, 0.2500, 0.2000],\n",
       "        [0.1900, 0.2400, 0.2200,  ..., 0.1100, 0.2000, 0.1900]])"
      ]
     },
     "execution_count": 138,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "is_correct_avg_pp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs4AAAFwCAYAAACoxP20AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmcLGdZ9//PFxICISGhAxx2IkvYBg14QPTBEEEJKIIbonkUxOgo/hDZQXkeISoYBEEjImcAAWVfFDAsD4gkIpHlBMKhkR0DCRAEJxB2SHL9/qgaGOdM11T3meqZOefzfr3mdXq66qrr7qq7a66uc/ddqSokSZIkdbvSVjdAkiRJ2gksnCVJkqQeLJwlSZKkHiycJUmSpB4snCVJkqQeLJwlSZKkHiycJWmOkjwxyYvmHStJOnAWzpIOCUnukuTcJF9OspzkHUnuuNXtkiTtHIdtdQMkaWhJrgGcBTwYeAVwFeBHgW9tZbuGlOSwqrpsq9shSQcTrzhLOhScAFBVL62qy6vqG1X15qraB5DkZkn+Jcl/J/likhcnOXYlOMkFSR6dZF+SryV5XpJdSd6Y5CtJ/jnJNdt1j09SSRaTfDbJ55I8clLDkty5vRL+pSTvT3LyqmXfl+ScNsdbgGt1bOfkJBcleWySi4Hnt8/fO8n57fbPTfL9q2Iem+Qz7fY/kuTu7fNPTPKqJC9vl703yQ+sirt1krPbbX4wyX1WLXtBkr9O8vo29l1JbtYuS5JnJPmv9sr/viQL7bIjkjwtyaeTfD7Js5NcrV12rSRntfmWk7w9iX+/JM2dJx5Jh4KPApcneWGSe60UuasE+FPg+sCtgRsBT1yzzs8DP0FThP808EbgD2iK2SsBD12z/o8BtwDuATwuyY+vbVSSGwCvB/4EGAGPAl6d5NrtKi8Bzmtz/DHwwA1e53Xb7dwEWExyB+Bvgd8CjgP2AK9ri9RbAg8B7lhVRwOnABes2tZ9gVe223sJ8Jokhyc5HPgn4M3AdYDfBV7cbm/FLwOnA9cEPg48qX3+HsBJNPvwWOD+wH+3y57SPn8icHPgBsAftsseCVwEXBvYRbPfa4N9IUmbzsJZ0kGvqi4F7kJTbD0H+EKS1yXZ1S7/eFW9paq+VVVfAJ4O3HXNZv6qqj5fVZ8B3g68q6reV1XfAv4RuP2a9U+vqq9V1Qdorv7+8jpN+xXgDVX1hqq6oqreAuwFfjLJjYE7Av+3bde/0hSsXa4AntCu/w3gN4E9VfWu9kr7C2mGp9wZuBw4ArhNksOr6oKq+sSqbZ1XVa+qqu+0++OqbdydgaOAM6rq21X1LzTDYFa/vn+oqne3Q0VeTFMMA3wHOBq4FZCq+lBVfS5J2rY+vKqWq+orwJOBX1oVdz3gJlX1nap6e1VZOEuaOwtnSYeEtkj7taq6IbBAc3X5LwCSXCfJy9phC5cCL2L/YRGfX/X4G+v8ftSa9S9c9fhTbb61bgLcrx2C8KUkX6Ip8K/Xrn9JVX1tzXa6fKGqvrlm+49cs/0bAdevqo8DD6O5sv5f7etf3cbvtr+qrqC54nv99ufC9rnV7brBqt8vXvX467T7pi2ynwn8NfD5JEvt+PNrA0cC561q55va5wGeSnPl+s1JPpnkcRvsB0kahIWzpENOVX0YeAFNAQ3NMI0Cvr+qrkFzJTgHmOZGqx7fGPjsOutcCPx9VR276ufqVXUG8DngmkmuvmY7XdZehb0QeNKa7R9ZVS8FqKqXVNVdaArsohkusV/72/HEN2xfw2eBG60ZY3xj4DMbtI0255lV9YPAbWmGZjwa+CLNh4/brmrnMVW1UnB/paoeWVU3pRkm84iV8diSNE8WzpIOekluleSRSW7Y/n4jmqEF72xXORr4KvCldtzxozch7f9NcmSS2wIPAl6+zjovAn46ySlJrpzkqu2X/G5YVZ+iGbZxepKrJLkLTdE4jecAv53kh9ov5l09yU8lOTrJLZPcLckRwDdpCtfLV8X+YJKfS3IYzZXpb9Hsr3cBXwMe0455Prlt18s2akySO7ZtObzdxjeBy9ur188BnpHkOu26N0hySvv43klu3g7puLRt5+XrZ5Gk4Vg4SzoUfAX4IeBdSb5GUwCOab50Bs0X2e4AfJnmy3r/sAk5z6EZXvBW4GlV9ea1K1TVhTRfwvsD4As0V4gfzffOzae27V4GngD83TQNqKq9NGOHnwlc0rbn19rFRwBn0FztvZjmi35/sCr8tTRf3rsE+FXg59rxxd8G7gPcq419FvCA9ir+Rq5BUyBfQjO847+Bp7XLHtu2753tcJl/Bla+cHiL9vevAv8OPKuqzu65GyRp08TvV0jS5klyPPCfwOE7dR7lJE8Ebl5Vv7LVbZGk7cQrzpIkSVIPFs6SJElSDw7VkCRJknrwirMkSZLUg4WzJEmS1MNhW92ADo4hkSRJ0jz0uunVdi6cWV5enjpmNBqxsLCw8YqrjMdjgLnHTfv6RqPRAeWbJW7amJW4eR27rcp3sLdz3q/vYN0vO+m9Pu+47X7stipup7wXPMdvbdxOeq/DfM+Bsx6DvhyqIUmSJPVg4SxJkiT1MNhQjSS3ormV7A1oxit/FnhdVX1oqJySJEnSUAa54pzkscDLaAZavxt4T/v4pUkeN0ROSZIkaUhDXXE+DbhtVX1n9ZNJng58EDhjoLySJEnSIIYa43wFcP11nr9eu2xdSRaT7E2yd2lpaaCmSZIkSdMb6orzw4C3JvkYcGH73I2BmwMPmRRUVUvASsVcs0w/I0mSJA1hkMK5qt6U5ATgTjRfDgxwEfCeqrp8iJySJEnSkAabVaOqrgDeOdT2JUmSpHlyHmdJkiSpBwtnSZIkqYdU1Va3YZJt2zBJkiQdVNJnpcHGOG+GhYWFqWPG4/HUcePxeKZ8Bxo37awho9Fo7vlmmdnkQOLmdcwPNM79srlx894v88q3VeeWnXBO2u7HbqvyGbez4zy3rB8367llnnVLXw7VkCRJknqwcJYkSZJ6sHCWJEmSerBwliRJknqYe+Gc5EHzzilJkiQdqK244nz6FuSUJEmSDsgg09El2TdpEbCrI24RWATYs2fPAC2TJEmSZjPUPM67gFOAS9Y8H+DcSUFVtQQsrfx65plnDtM6SZIkaUpDFc5nAUdV1flrFyQ5e6CckiRJ0mAGKZyr6rSOZacOkVOSJEkaktPRSZIkST1YOEuSJEk9pKq2ug2TbNuGSZIk6aCSPisN9eXATbG8vDx1zGg0YmFhYaqY8XgMMHPctO0cjUYHlG+ecdPGHCpx8+qbW5XvYG/nvOK24j0Ls5+T5nku2+7HzrjuuFnf654j1o87mOsI2Bnnlr4cqiFJkiT1YOEsSZIk9WDhLEmSJPUwWOGc5FZJ7p7kqDXP33OonJIkSdJQBimckzwUeC3wu8A4yX1XLX7yEDklSZKkIQ01q8ZvAj9YVV9NcjzwqiTHV9Vf0nO6D0mSJGk7GWqoxpWr6qsAVXUBcDJwryRPp6NwTrKYZG+SvUtLSwM1TZIkSZreUFecL05yYlWdD9Beeb438LfA7SYFVdUSsFIx1yzzPUqSJElDGOqK8wOAi1c/UVWXVdUDgJMGyilJkiQNZpArzlV1UceydwyRU5IkSRqS8zhLkiRJPVg4S5IkST2kqra6DZNs24ZJkiTpoNJruuShZtXYFAsLC1PHjMfjqePG4/FM+Q6FuFmPwSwzooxGo7kdc+OMm0fcyntv2vfDaDQCdsY5YtY4zxGbGzfr/px33MG+P3dK3XIwn5NmPQZ9OVRDkiRJ6sHCWZIkSerBwlmSJEnqYbAxzknuBFRVvSfJbYB7Ah+uqjcMlVOSJEkayiCFc5InAPcCDkvyFuCHgLOBxyW5fVU9aYi8kiRJ0lCGuuL8C8CJwBE0t96+YVVdmuSpwLsAC2dJkiTtKEMVzpdV1eXA15N8oqouBaiqbyS5YlJQkkVgEWDPnj0DNU2SJEma3lBfDvx2kiPbxz+48mSSY4CJhXNVLVXV7qravbi4OFDTJEmSpOkNdcX5pKr6FkBVrS6UDwceOFBOSZIkaTCDFM4rRfM6z38R+OIQOSVJkqQhOY+zJEmS1IOFsyRJktRDqmqr2zDJtm2YJEmSDirps9Jgdw7cDMvLy1PHjEYjFhYWpooZj8cAM8dN287RaHRAcfN8fbMeg3kdOziwds6az7jNjTtYj9+Bnlt2Stws57Ltfux2WtzB+h5aiTvYX9+864+dcm6ZZ93Sl0M1JEmSpB4snCVJkqQeLJwlSZKkHuZWOCf5u3nlkiRJkjbbIF8OTPK6tU8BP5bkWICqus8QeSVJkqShDDWrxg2B/wCeSzOtXIDdwJ8PlE+SJEka1FBDNXYD5wGPB75cVWcD36iqc6rqnElBSRaT7E2yd2lpaaCmSZIkSdMb5IpzVV0BPCPJK9t/P98nV1UtASsVc80yF58kSZI0hEFvgFJVFwH3S/JTwKVD5pIkSZKGNJc7B1bV64HXzyOXJEmSNATncZYkSZJ6sHCWJEmSekhVbXUbJtm2DZMkSdJBJX1WmssY51ktLCxMHTMej6eOG4/HAEw7i8doNAKmb+dKvp0QN69jcKBxs8zAMhqNdky+nRK3U/bLrO30HLF5cTulTxtn3HaM26r3+rzPgbPkm/UY9OVQDUmSJKkHC2dJkiSpBwtnSZIkqYe5jHFOchfgTsC4qt48j5ySJEnSZhrkinOSd696/JvAM4GjgSckedwQOSVJkqQhDTVU4/BVjxeBn6iq04F7AP97oJySJEnSYIYaqnGlJNekKcxTVV8AqKqvJblsUlCSRZpCmz179gzUNEmSJGl6QxXOxwDn0UwmXUmuW1UXJzmKjgmmq2oJWFr59cwzzxyoeZIkSdJ0Bimcq+r4CYuuAH52iJySJEnSkOZ658Cq+jrwn/PMKUmSJG0G53GWJEmSerBwliRJknpIVW11GybZtg2TJEnSQWXi5BWrzXWM87QWFhamjhmPx1PHjcfjmfIdCnHzOgbGGXcwxu2k9/q845aXl6eKARiNRnM/5rZz/XwHe9x2358r771p2zkajYCdcY6A2V7frMegL4dqSJIkST1YOEuSJEk9WDhLkiRJPQxSOCf5oSTXaB9fLcnpSf4pyVOSHDNETkmSJGlIQ11x/lvg6+3jv6S5BfdT2ueeP1BOSZIkaTBDzapxpaq6rH28u6ru0D7+tyTnD5RTkiRJGsxQV5zHSR7UPn5/kt0ASU4AvjMpKMlikr1J9i4tLQ3UNEmSJGl6QxXOvwHcNckngNsA/57kk8Bz2mXrqqqlqtpdVbsXFxcHapokSZI0vUGGalTVl4FfS3I0cNM2z0VV9fkh8kmSJElDG/TOgVX1FeD9Q+aQJEmS5sF5nCVJkqQeLJwlSZKkHlJVW92GSbZtwyRJknRQSZ+VBh3jfKCWl5enjhmNRiwsLEwVMx6PAeYeN+3rG41GB5RvlrhpY3Za3Kx9bF59Ew6snTvlOBys+2UnvdfnHXew97Gdcm7ZKXHul/1jYGe81+cdN+sx6MuhGpIkSVIPFs6SJElSDxbOkiRJUg+DFM5JHprkRkNsW5IkSdoKQ11x/mPgXUnenuR3klx7oDySJEnSXAxVOH8SuCFNAf2DwH8keVOSB7a34ZYkSZJ2lKEK56qqK6rqzVV1GnB94FnAPWmK6nUlWUyyN8nepaWlgZomSZIkTW+oeZz/xyTSVfUd4HXA65JcbVJQVS0BKxVzzTJvoyRJkjSEoa4433/Sgqr6xkA5JUmSpMEMUjhX1UeH2K4kSZK0VZzHWZIkSerBwlmSJEnqIVW11W2YZNs2TJIkSQeVbLzKcLNqbIqFhYWpY8bj8dRx4/F4pnwHGjftrCGj0Wju+eZ1DA6VuFlmijmQ4zBrvp3Szu1+3Lfqvb4T4g7WY27c9oo7WPuZ55bJcbMeg74cqiFJkiT1YOEsSZIk9WDhLEmSJPUwyBjnJFcBfgn4bFX9c5JTgR8BPgQstXcSlCRJknaMob4c+Px220cmeSBwFPAPwN2BOwEPHCivJEmSNIihCufbVdX3JzkM+Axw/aq6PMmLgPdPCkqyCCwC7NmzZ6CmSZIkSdMbaozzldrhGkcDRwLHtM8fARw+Kaiqlqpqd1XtXlxcHKhpkiRJ0vSGuuL8PODDwJWBxwOvTPJJ4M7AywbKKUmSJA1mkMK5qp6R5OXt488m+Tvgx4HnVNW7h8gpSZIkDWmwOwdW1WdXPf4S8KqhckmSJElDcx5nSZIkqQcLZ0mSJKmHVNVWt2GSbdswSZIkHVTSZ6XBxjhvhoWFhaljxuPx1HHj8XimfAcat7y8PFXcaDQ6oHyzxM3rGOy0uGmPHTTHb6e8vp2yX7b7cdhJ73WYbzvnfcy3e1/ZaXHuz62N26q6ZSfEzXoM+nKohiRJktSDhbMkSZLUg4WzJEmS1MNgY5yT3Az4WeBGwGXAx4CXVtWXh8opSZIkDWWQK85JHgo8G7gqcEfgajQF9L8nOXmInJIkSdKQhrri/JvAiVV1eZKnA2+oqpOT7AFeC9x+vaAki8AiwJ49ewZqmiRJkjS9Icc4rxTlRwBHA1TVp4HDJwVU1VJV7a6q3YuLiwM2TZIkSZrOUFecnwu8J8k7gZOApwAkuTYw/eSPkiRJ0hYbpHCuqr9M8s/ArYGnV9WH2+e/QFNIS5IkSTvKYLNqVNUHgQ8OtX1JkiRpnpzHWZIkSerBwlmSJEnqIVW18UrJLuDJwPWr6l5JbgP8cFU9b8C2bdwwSZIk6cCl10o9C+c3As8HHl9VP5DkMOB9VXW7A2tjp1pYWJg6aDweM23ceDwGMG6duHkdg50Wt7w8/eQwo9Fo5rh5t3OnHIftHreT3uvA1P1lNBrNnG+7H7uVuJ3yHtop7fT1bU7cTju37IS6hZ6Fc9+hGteqqlcAVwBU1WXA5VO3TJIkSdqh+hbOX0tyHO3wiSR3Br48WKskSZKkbabvdHSPAF4H3CzJO4BrA78wWKskSZKkbaZX4VxV701yV+CWNGNAPlJV3xm0ZZIkSdI20lk4J/m5CYtOSEJV/cMAbZIkSZK2nY2uOP90x7IC1i2ckxwD/D7wMzTDOgD+C3gtcEZVfWlC3CKwCLBnz54NmiZJkiTNT2fhXFUPmnG7rwD+BTi5qi4GSHJd4IHAK4GfmJBvCVha+fXMM8+cMb0kSZK0uXrNqpHkuCRnJnlvkvOS/GU7y8Ykx1fVU1aKZoCquriqngLc+EAbLUmSJM1b3+noXgZ8Afh5mtk0vgC8vGP9TyV5THvHQaC5+2CSxwIXztpYSZIkaav0LZxHVfXHVfWf7c+fAMd2rH9/4DjgnCTLSZaBs4ERcL8DarEkSZK0BfoWzm9L8ktJrtT+/CLw+kkrV9UlVfXYqrpVVY3an1tX1WNpvjAoSZIk7Sh9C+ffAl4CfLv9eRnwiCRfSXLplDlPn3J9SZIkaculqjZ/o8m+SYuAE6rqiB6b2fyGSZIkSftLn5X63nKbJPcBTmp/PbuqzupYfRdwCnDJOo06t2/OhYWFvqt+13g8Znl5eaqY0Wg0U77xeHzQx816DOZ17KA5fvNup3HbI27e/cVzy+bF+V7f3Dj356EVt1Xv9XmfA2fJN+sx6KtX4ZzkDOCOwIvbp34vyV2q6nETQs4Cjqqq89fZ1tm9WydJkiRtE32vOP8kcGJVXQGQ5IXA+4B1C+eqOm3Shqrq1GkbKUmSJG21vl8OhP85/dwxm90QSZIkaTvre8X5T4H3JXkbzTjlk4DfH6xVkiRJ0jbTq3Cuqpe2Y5PvSFM4P3b17bQlSZKkg12voRpJAtydZpzza4GrJLnTLAmTvLFj2WKSvUn2Li0tzbJ5SZIkaRB9h2o8C7gCuBvwR8BXgFfTXIHeT5I7TNhOgBMnJamqJWClYq4zzzyzZ/MkSZKkYfUtnH+oqu6Q5H3Q3FI7yVU61n8PcA7rTyZ97DrPSZIkSdta38L5O0muTHs3vyTXprkCPcmHgN+qqo+tXZDkwqlbKUmSJG2xvtPRnQn8I3CdJE8C/g14csf6T+zY9u/2bp0kSZK0TfSdVePFSc6j+YJggJ+pqg91rP+qjs1dc7omSpIkSVuvs3BOclXgt4GbAx8A9lTVZQeY83Tg+Qe4DUmSJGmuUlWTFyYvB74DvB24F3BBVT1sw40m+yYtAk6oqiN6tG1ywyRJkqTNs96EFvvZaKjGbarqdgBJnge8u2fyXcApwCXrNOrcnttgeXm576rfNRqNWFhYmCpmPB7PlG80Gh1Q3KztnOfrmzbXSj7jtkfcrO+heb33YGftz53wnt2KuFle33Y/ditxvoeM245xB3pu2SnnsnmeW/raqHD+zsqDqrqsuQ9KL2cBR1XV+WsXtHcglCRJknaUjQrnH0hyafs4wNXa3wNUVV1jvaCqOm3SBqvq1JlaKkmSJG2hzsK5qq48r4ZIkiRJ21nfeZwlSZKkQ5qFsyRJktTDIIVzkmsk+dMkf5/k1DXLntURt5hkb5K9S0tLQzRNkiRJmkmvOwfO4PnAx4BXA7+e5OeBU6vqW8CdJwVV1RKwUjHXLNP5SJIkSUMYaqjGzarqcVX1mqq6D/Be4F+SHDdQPkmSJGlQQ11xPiLJlarqCoCqelKSi4B/BY4aKKckSZI0mKGuOP8TcLfVT1TVC4FHAt8eKKckSZI0mEGuOFfVYyY8/6YkTx4ipyRJkjSkrZiO7vQtyClJkiQdkFTV5m802TdpEXBCVR3RYzOb3zBJkiRpf+mz0lBfDtwFnAJcsub5AOf23cjCwsLUicfj8dRx4/F4pnyHQtwsUwKORqOZ4+Z1zI3rjvP4rR837X4ZjUbAznivzzvOPrY94uZ9HDzumxO3k97r846b9Rj0NVThfBZwVFWdv3ZBkrMHyilJkiQNZqgvB57WsezUScskSZKk7WorvhwoSZIk7TgWzpIkSVIPFs6SJElSD4MUzkmum+Rvkvx1kuOSPDHJB5K8Isn1OuIWk+xNsndpaWmIpkmSJEkzGeqK8wuA/wAuBN4GfAP4KeDtwLMnBVXVUlXtrqrdi4uLAzVNkiRJmt5QhfOuqvqrqjoDOLaqnlJVn66qvwJuMlBOSZIkaTBDFc6rt/t3a5ZdeaCckiRJ0mCGKpxfm+QogKr6PytPJrk58JGBckqSJEmDGeoGKH844fmPJ3n9EDklSZKkIW3FdHSnb0FOSZIk6YCkqjZ/o8m+SYuAE6rqiB6b2fyGSZIkSftLn5UGGaoB7AJOAS5Z83yAc/tuZHl5eerEo9GIhYWFqWLG4zGAcevETRtzqMTNq2+u5DNuZ8atvPem7S+j0QjYGecImO31bfdjZ5xx2zluJ9URMPs5cJ7nlr6GKpzPAo6qqvPXLkhy9kA5JUmSpMEM9eXA0zqWnTpETkmSJGlIW/HlQEmSJGnHsXCWJEmSephb4ZzkOvPKJUmSJG22QQrnJKM1P8cB705yzSSjjrjFJHuT7F1aWhqiaZIkSdJMhppV44vAp9Y8dwPgvTTzM990vaCqWgJWKuaaZcovSZIkaQhDDdV4DPAR4D5V9X1V9X3ARe3jdYtmSZIkaTsbpHCuqqcBvwH8YZKnJzka7wQoSZKkHWywLwdW1UVVdT/gbcBbgCOHyiVJkiQNbfBZNarqn4AfA34cIMmDhs4pSZIkbba5TEdXVd+oqpUbgZ8+j5ySJEnSZkrV5g89TrJv0iLghKo6osdmHBMtSZKkeUiflYaajm4XcApwyZrnA5zbdyMLCwtTJx6Px1PHjcfNxfBpp78bjZopqWfNN++4WV7fvI7BgcbNMnXhgby+WfPNO87jt36+eZ9bjNs/brsfu62KO1jfCytxB/s5cF7tPNC6Zd5xO6Fu6Wuowvks4KiqOn/tgiRnD5RTkiRJGswghXNVndax7NQhckqSJElDmsuXAyVJkqSdzsJZkiRJ6sHCWZIkSephkMI5yT1XPT4myfOS7EvykiS7OuIWk+xNsndpaWmIpkmSJEkzGeqK85NXPf5z4HPATwPvAfZMCqqqparaXVW7FxcXB2qaJEmSNL2hpqNbbXdVndg+fkaSB84hpyRJkrSphiqcr5PkETQ3PLlGktT3blHouGpJkiTtOEMVsc8BjgaOAl4IXAsgyXWB/W6KIkmSJG13Q90A5fQJz1+c5G1D5JQkSZKGlO+NoJhTwuTTVXXjHqvOt2GSJEk6VKXPSoNccU6yb9IiYOJ0dJIkSdJ2NdSXA3cBpwCXrHk+wLl9N7K8vDx14tFoxMLCwlQx4/EYwLh14uZ1DFbyGWfcPOLmfW6ZNt9oNAJ2xjli1jjPLdsjzuOwPeIO9nPEPF/frMegr6EK57OAo6pqvy8CJjl7oJySJEnSYIb6cuBpHctOHSKnJEmSNCTnVJYkSZJ6sHCWJEmSephb4ZzkuHnlkiRJkjbbIIVzkjOSrNwtcHeSTwLvSvKpJHftiFtMsjfJ3qWlpSGaJkmSJM1kqFk1fqqqHtc+fipw/6p6T5ITgJcAu9cLqqolYKVirlmmyZEkSZKGMNRQjcOTrBTlV6uq9wBU1UeBIwbKKUmSJA1mqML5r4E3JLkb8KYkf5HkpCSnA/vN7SxJkiRtd0PN4/xXST4APBg4oc1zAvAa4I+HyClJkiQNaagxzlTV2cDZa59P8iDg+UPllSRJkoaQqppvwuTTVXXjHqvOt2GSJEk6VKXPSoNccU6yb9IiYNcQOSVJkqQhDTVUYxdwCnDJmucDnNt3IwsLC1MnHo/HU8eNx2MApp3+bjQaAdO3cyXfvONmeX3zOgYrcbNMQTgajeYeN+/9crDH7ZTjPus54mA/t8wSt1PeezulnQf7fvE47B8D8z+37IRz4KzHoK+hCuezgKOqar8ZNJKcPVBOSZIkaTBDzapxWseyU4fIKUmSJA1pqHmcJUmSpIOKhbMkSZLUg4WzJEmS1MMghXOS9yb5P0luNmXcYpK9SfYuLS0N0TRJkiRpJkNdcb4mcCzwtiTvTvLwJNffKKiqlqpqd1XtXlxcHKhpkiRJ0vSGKpwvqapHtXcIfCRwC+C9Sd6WxIpYkiRJO87gY5yr6u1V9TvADYCnAD88dE5JkiRpsw11A5SPrn2iqi4H3tS/q8avAAAWcElEQVT+SJIkSTvKIFecq+qXJi1L8qAhckqSJElDSlXNN2Hy6Xbs80bm2zBJkiQdqtJnpUGGaiTZN2kRsGuInJIkSdKQhhrjvAs4BbhkzfMBzu27kYWFhakTj8fjqePG4/FM+Q6FuHkdA+MO7bjl5eWp40aj0bbPd6Dv2WnbORqNDijfPOPmfcwP9j59sLfTuP1jYP7niFnzzbOdsx6DvoYqnM8Cjqqq89cuSHL2QDklSZKkwQxSOFfVaR3LTh0ipyRJkjSkwedxliRJkg4GFs6SJElSDxbOkiRJUg+DFM5Jdid5W5IXJblRkrck+XKS9yS5fUfcYpK9SfYuLS0N0TRJkiRpJkNdcX4W8GfA62mmn9tTVccAj2uXrauqlqpqd1XtXlxcHKhpkiRJ0vSGKpwPr6o3VtVLgaqqV9E8eCtw1YFySpIkSYMZqnD+ZpJ7JLkfUEl+BiDJXYHLB8opSZIkDWaoG6D8Ns1QjSto7iD44CQvAD4D/OZAOSVJkqTBDHLFuareX1WnVNW9qurDVfV7VXVsVd0WuOUQOSVJkqQhparmmzD5dFXduMeq822YJEmSDlXps9IgQzWS7Ju0CNg1RE5JkiRpSEONcd5FM7b5kjXPh2Z6ul4WFhamTjwej6eOG4/HM+U7FOLmdQyMM+5gjNtJ73WA5eXlqeJGo9HM+bb7sTPOuO0ct5Pe6/OOm/UY9DVU4XwWcFRVnb92QZKzB8opSZIkDWaQwrmqTutYduoQOSVJkqQhDTWPsyRJknRQsXCWJEmSerBwliRJknoYpHBOclSSP0rywSRfTvKFJO9M8msbxC0m2Ztk79LS0hBNkyRJkmYy1BXnFwOfpJmS7nTgTOBXgR9L8uRJQVW1VFW7q2r34uLiQE2TJEmSpjdU4Xx8Vb2gqi6qqqcD96mqjwEPAn5uoJySJEnSYIYqnL+W5C4ASX4aWAaoqivoeUtDSZIkaTsZ6gYovw08N8kJwBj4dYAk1wb+eqCckiRJ0mCGugHKPuBO6zz/hSRfGSKnJEmSNKRU1XwTJp+uqhv3WHW+DZMkSdKhqtdQ4kGuOCfZN2kRsGuInJIkSdKQhhrjvItmKrpL1jwf4Ny+G1lYWJg68Xg8njpuPB7PlO9QiJv1GCwvL08dNxqNZo6bV185VOIO9uMwr/45Go2AnfFen3fcTukrO+U9tFPes/5tWD/Oc8vmxc3aV/oaqnA+Cziqqs5fuyDJ2QPllCRJkgYz1JcDT+tYduoQOSVJkqQhDTWPsyRJknRQsXCWJEmSehikcE5yTJIzknw4yX+3Px9qnzt2iJySJEnSkIa64vwKmhk1Tq6q46rqOODH2udeOSkoyWKSvUn2Li0tDdQ0SZIkaXpDFc7HV9VTqurilSeq6uKqegow8eYnVbVUVburavfi4uJATZMkSZKmN1Th/Kkkj0ny3ZudJNmV5LHAhQPllCRJkgYzVOF8f+A44JwklyRZBs4GRsAvDpRTkiRJGsxQ8zhfkuT5wFuAd1bVV1eWJbkn8KYh8kqSJElDGWpWjYcCrwUeAoyT3HfV4icPkVOSJEkaUqpq8zeafAD44ar6apLjgVcBf19Vf5nkfVV1+x6b2fyGSZIkSftLn5UGGaoBXHlleEZVXZDkZOBVSW7St2EACwsLUycej8dTx43H45nyHWjc8vLyVHGj0Wju+eZ1DLYqbtp9AofGfjnY9+e82rlV54h5n1tmidspfXq797Gtymfc1sbttPpjJ5xb+hrqy4EXJzlx5Ze2iL43cC3gdgPllCRJkgYzVOH8AODi1U9U1WVV9QDgpIFySpIkSYMZalaNizqWvWOInJIkSdKQhrriLEmSJB1ULJwlSZKkHuZeOCd547xzSpIkSQdqkDHOSe4waRFw4oRlJFkEFgH27NkzQMskSZKk2Qw1j/N7gHNYf87mYycFVdUSsLTy65lnnjlA0yRJkqTpDVU4fwj4rar62NoFSS4cKKckSZI0mKHGOD+xY9u/O1BOSZIkaTCDFM5V9SogSe6e5Kg1i785RE5JkiRpSIMUzkkeCryW5uryOMl9Vy1+8hA5JUmSpCGlqjZ/o8kHgB+uqq8mOR54FfD3VfWXSd5XVbfvsZnNb5gkSZK0v/UmtNjPUGOcr1xVXwWoqguAk4F7JXl634a16637k+S3upYbZ9yhErcT2micccbtvLid0EbjjNvkuH6qatN/gH8BTlzz3GHA3wGXb8L29xpnnHE7o43GGWfczovbCW00zrh5xa3+GeqK8wOAi1c/UVWXVdUDgJMGyilJkiQNZpB5nKvqoo5l7xgipyRJkjSkoa44D21p41WMM+6QiNsJbTTOOON2XtxOaKNxxs0r7rsGmVVDkiRJOtjs1CvOkiRJ0nwd6LcL5/0D3BP4CPBx4HE9Y/4W+C9gPGWuGwFvAz4EfBD4vZ5xVwXeDby/jTt9ipxXBt4HnDVlWy8APgCcT89vjQLH0syx/eH2Nf5wj5hbtjlWfi4FHtYz38Pb/TEGXgpctWfc77UxH+zKtd5xBkbAW4CPtf9es2fc/dp8VwC7p8j31HZ/7gP+ETi2Z9wftzHnA28Grj9NPwYeRTP3+bV65nsi8JlVx/En++ajubHRR9r982c98718Va4LgPN7xp0IvHOlXwN36hn3A8C/t++JfwKusU7cuu/vrj7TEdPZXzriOvtLR1xnf5kUt1F/6cjX2V+68nX1l458nf2lI66zv3TEdfYXJpzTge8D3tX2lZcDV+kZ9xCav2GT3rOT4l7c7ssxTb8/vGfc89rn9tGc84/qE7dq+V8BX52inS8A/nPVMVw7y9akuABPAj7aHqOH9ox7+6pcnwVe0zPu7sB727h/A27eM+5ubdwYeCFw2Dr75n/8Ld+or3TEdfaVjrjOvtIR19lXJsVt1Fc68nX2lY64zr7SEdfZV/r8TLXyVv+0O+ATwE2Bq7QH9zY94k4C7sD0hfP1gDu0j49uD1CffFnpbMDh7Zvlzj1zPgJ4ydrO2CPugq431YSYFwK/0T6+CusUeT2Ox8XATXqse4P2zXG19vdXAL/WI26hfeMfSfNl1n8GbtH3OAN/RvsBC3gc8JSecbem+ZBwNpML5/Xi7kF7EgWeMkW+a6x6/FDg2X37MU0x8P+AT63XBybkeyLwqA32/XpxP9YegyPa36/Tt52rlv858Ic9870ZuFf7+CeBs3vGvQe4a/v414E/Xidu3fd3V5/piOnsLx1xnf2lI66zv0yK26i/dOTr7C8dcZ39paudXf2lI19nf+mI6+wvTDin05zHfql9/tnAg3vG3R44ngnn7Y64n+R7c86+dIp8q/vL01lz0WlSXPv7buDvWb9wnpTvBcAvdPSXSXEPopm29koT+suGf1uBVwMP6Jnvo8Ct2+d/B3hBj7gfAS4ETmif/yPgtHVe4//4W75RX+mI6+wrHXGdfaUjrrOvTIrbqK905OvsKx1xnX2lq51dfaXPz04bqnEn4ONV9cmq+jbwMuC+G8RQVf8KLE+brKo+V1XvbR9/heZTzQ16xFW1N4ChebMdTo87ISa5IfBTwHOnbeu0klyDpuB4HkBVfbuqvjTlZu4OfKKqPtVz/cOAqyU5jKYQ/myPmFsD76yqr1fVZcA5wM+ut+KE43xfmg8ItP/+TJ+4qvpQVX2kq2ET4t7cthOaK1837Bl36apfr846/aWjHz8DeMx6MRvEdZoQ92DgjKr6VrvOf02TL0mAX6Q5kfeJK+Aa7eNjWKfPTIi7JfCv7eO3AD+/Ttyk9/fEPjMpZqP+0hHX2V864jr7ywbnron95QDOeZPiOvvLRvkm9ZeOuM7+0hHX2V86zul3o7kqB+ucXybFVdX7qrk52Lo64t7QLiuaK6Fr+8ukuEvhu/vzauzfX9aNS3Jlmv8Vecw07Zz0unrEPRj4o6q6ol1vbX/pzJfkaJpj8pqecRv1l/XiLge+VVUfbZ/fr7+s/Vve7vfOvrJeXNuGzr7SEdfZVzriOvvKpLiN+sqkuD4mxHX2lY3yTeorfey0wvkGNJ/0VlxEj5P6ZmhvHX57mk+cfda/cpLzaf4L+S1V1SfuL2g63RUzNLGANyc5L8lij/VvCnwBeH6S9yV5bpKrT5nzl1inAFq3cVWfAZ4GfBr4HPDlqnpzj9AxcFKS45IcSfMp+kZTtHFXVX2ubcPngOtMEXugfh14Y9+VkzwpyYXA/wb+sGfMfYDPVNX7Z2jfQ5LsS/K3Sa7ZM+YE4EeTvCvJOUnuOGXOHwU+X1Uf67n+w4CntvvlacDv94wbA/dpH9+PDfrMmvd3rz4z7TmhR1xnf1kb17e/rI6bpr+s085e/WVNXO/+MmG/bNhf1sT17i9r4jbsL2vP6TT/+/mlVR981v17NOPfgs64JIcDvwq8qW9ckufT/A/hrWj+O71P3EOA1628H6Zs55Pa/vKMJEf0jLsZcP8ke5O8McktptkvNBdV3rrmg2VX3G8Ab0hyEc3+PGOjOJoi9PAku9tVfoH9+8vav+XH0aOvrBPX18S4rr4yKW6jvjIhbsO+0tHOzr4yIW7DvtKRDzr6ykZ2WuG83i0RN/x0e8BJk6NoLuk/rO9OrqrLq+pEmk95d0qysEGOewP/VVXnzdjM/1VVdwDuBfx/STa60cxhNP+9/TdVdXvgazT/Ld1LkqvQ/KF5Zc/1r0lzJe/7gOsDV0/yKxvFVdWHaP4L+y00b/z3A5d1Bm0DSR5P084X942pqsdX1Y3amIf0yHEk8Hh6Ftlr/A3NiedEmg8yf94z7jDgmjT/zflo4BXtlYm+fpmeH7ZaDwYe3u6Xh9P+D0kPv07zPjiP5r/kvz1pxVne37PEdMVt1F/Wi+vTX1bHtdvv1V/Wyderv6wT16u/dOzPzv6yTlyv/rJO3Ib9Ze05neZ/w/ZbbaO4jf4W9Ix7FvCvVfX2vnFV9SCac++HgPv3iDuJ5kPEeoXTRvl+n6bouiPNdwYe2zPuCOCbVbUbeA7N2Nxp9svE/jIh7uE04/VvCDyfZmhCZxxwW5qLRs9I8m7gK6z6mzThb/mGtcusNUCPuHX7SldcV19ZLy7J9dmgr3Tk6+wrHXGdfaXHfpn2b9H31JRjO7byB/hh4P+t+v33gd/vGXs8U45xbuMOpxkP+IgDaPcT2Hg86Z/SfAq9gOaT3teBF82Y74k98l0XuGDV7z8KvH6KHPcF3jzF+vcDnrfq9wcAz5rhtT0Z+J2+x5nmyxHXax9fD/jINP2DjjHOk+KAB9J80ejIWfojcJOOZd+NA25HcxXkgvbnMpor+tedMl/vZTQfXk5e9fsngGv33C+HAZ8HbjjF8fsy35s2M8ClM7yGE4B3T1i23/t7oz6zXkyf/jIpbqP+0pWvq7+sjevbX3rkW3dfT9iXG/aXjv3S2V8m5Nuwv/R4fRP7y6p1nkDzQeCLfG+M+v/4+9QR96hVv19Aj++mrI5rH7+GdmznNPna5+7KBt+haeOeQPO3aKW/XEEzVHLafCf3zPcomi/KHr/q+H15iv1yHPDf9PjS+arj94lVz90Y+I8ZXt89gFes+n29v+Uv3qivTIh70arl6/aVrriuvrJRvkl9ZULcJRv1lZ759usrk+I26isb7JfefWXdfjBL0Fb90JxIP0lz1XLly4G37Rl7PNN/OTA0g8//Ysq4a9N+0Y5mjNDbgXtPEb9f59lg/asDR696fC5wzx5xbwdu2T5+IvDUKXK+DHjQFOv/EM23ko9s9+sLgd/tGXud9t8bt2+W/WbGmHScacZcrf6i136zQHT1D6YsnGlmffkP1ikmN4i7xarHvwu8atp+TMcf4XXyXW/V44cDL+sZ99s048qgKTAupC1UNmpnu2/OmXK/fIi28KIZU39ez7iVPnMlmvfwr68Ts+77u6vPTIrZqL905OrsLx1xnf1lo3ZO6i8d+Tr7S0dcZ3/pamdXf+nI19lfOuI6+wsTzuk0/+O2+gtfv9MnrusYbJDvN2jO71ebsF/Wi/tp2tki2tf/NOBp07SzfX69LwdOauf1VuX7C5px7n3izljZ9zR/B9/Tt51tX3vhFPvl3jTF7MqX/E4DXt0zbqW/HAG8FbjbhLwn870vs3X2lUlxG/WVjnydfWW9uPZ4dfaVjdo5qa90tLOzr3TEdfaVrnZ29ZU+PzMFbeUPzRjXj9JcuXh8z5iX0vz34ndoPoHs9w3YCXF3ofnvlJVpn9adtmuduO+nmf5kH83Yuf1mEOjbqXquf1OaDxErU+b03S8n0kzZtI/mU+nEgnRN3JE0n9aOmfJ1nU5T+I5pvnl7RM+4t9MUF+8H7j7Ncab5ZPlWmimA3gqMesb9bPv4WzRXvfa7kjQh7uM0xcFKf1lvdoz14l7d7pd9NNNh3WDafszkP8Lr5ft7mqm39gGvY1VhtEHcVWg+7Y9ppmPa7w/GpHbSfHv6t6c8fncBzmuP/buAH+wZ93s054mP0pxg1yvu131/d/WZjpjO/tIR19lfOuI6+8ukuI36S0e+zv7SEdfZX7ra2dVfOvJ19peOuM7+woRzOs25993tcXwla85pHXEPbfvLZTRfSHtuz7jLaP72rbR97Wwj+8XRfBh4R3v8xjRXP9dOt7fh3yzWL5wntfNfVuV7EftPfzcp7ljg9W3svwM/0LedNB9a171g1JHvZ9tc72/jb9oz7qk0H9I+QvcUqSfzvUKvs690xHX2lY64zr6yXlyfvjIp30Z9paOdnX2lI66zr3S1s6uv9PnxzoGSJElSDzvty4GSJEnSlrBwliRJknqwcJYkSZJ6sHCWJEmSerBwliRJknqwcJakLZLkhklem+RjST6Z5JkTbjk76/Z/JsltVv3+R0l+fLO2L0mHGgtnSdoC7a2n/wF4TVXdArgFzU0W/mwT0/wM8N3Cuar+sKr+eRO3L0mHFAtnSdoadwO+WVXPB6iqy2nuyveAJA9J8syVFZOcleTk9vE9kvx7kvcmeWWSo9rnz0jyH0n2JXlakh8B7gM8Ncn5SW6W5AVJfqFd/+5J3pfkA0n+duVKd5ILkpzebv8DSW41z50iSduZhbMkbY3b0tzl7ruq6lKaO/odtl5AkmsB/wf48aq6A82dPx+RZERzF7TbVtX3A39SVefS3OXv0VV1YlV9YtV2rkpzZ777V9Xt2nwPXpXqi+32/wZ41Ca8Vkk6KFg4S9LWCM3tn9d7fpI70wy9eEeS84EHAjcBLgW+CTw3yc8BX98g9y2B/6yqj7a/vxA4adXyf2j/PQ84foNtSdIhw8JZkrbGB4Hdq59Icg1gF/Df/M/z81VXVgHe0l5BPrGqblNVp1XVZcCdgFfTjGt+0wa5u4pzgG+1/17OhKvfknQosnCWpK3xVuDIJA8ASHJl4M+BZwL/CZyY5EpJbkRTFAO8E/hfSW7exhyZ5IR2nPMxVfUG4GHAie36XwGOXif3h4HjV7YD/Cpwzqa/Qkk6yFg4S9IWqKqiGZf8C0k+RnOV+YqqehLwDpri+QPA04D3tjFfAH4NeGmSfTSF9K1oiuOz2ufOofmSIcDLgEe3XwK82arc3wQeBLwyyQeAK4BnD/uKJWnnS3PuliRtpXYWjJcCP1dV5220viRp/iycJUmSpB4cqiFJkiT1YOEsSZIk9WDhLEmSJPVg4SxJkiT1YOEsSZIk9WDhLEmSJPVg4SxJkiT18P8DJW06XjKQsG0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f, ax = plt.subplots(figsize=(12,6))\n",
    "ax = sns.heatmap(is_correct.detach().numpy(), linewidths=0.25, ax=ax, square=True, cbar=False, cmap=sns.color_palette('Greys_r'))\n",
    "ax.set(title='Sampled responses', xlabel='Question', ylabel='People')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extended model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 179,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ext_model():\n",
    "    \n",
    "    alpha0 = torch.ones(questions_n) * 2.5\n",
    "    beta0 = torch.ones(questions_n) * 7.5\n",
    "    \n",
    "    guess_probs = pyro.sample('guess_probs', pyro.distributions.Beta(alpha0, beta0).independent(1))\n",
    "    \n",
    "    # setup empty arrays\n",
    "    has_skills = torch.zeros((people_n, questions_n))\n",
    "    is_correct_probs = torch.zeros((people_n, questions_n))\n",
    "    is_correct = torch.zeros((people_n, questions_n))\n",
    "    \n",
    "    # declare dimensions as conditionally independent\n",
    "    with pyro.plate('people', people_n):\n",
    "        with pyro.plate('skills', skills_n):\n",
    "            skills = pyro.sample('skill', pyro.distributions.Bernoulli(0.5))\n",
    "    \n",
    "    # declare dimensions as conditionally independent\n",
    "    for p in pyro.plate('people', people_n):\n",
    "        for q in pyro.plate('questions', questions_n):\n",
    "            \n",
    "            # skills_needed AND skills\n",
    "            has_skills[p,q] = torch.tensor(torch.allclose(skills_needed[:,q], skills[:,p])).float()\n",
    "            \n",
    "            # probabilites\n",
    "            is_correct_probs[p,q] = has_skills[p,q] * (1 - mistake_probs) + (1 - has_skills[p,q]) * guess_probs[q]\n",
    "            \n",
    "            # sample is correct\n",
    "            is_correct[p,q] = pyro.sample('is_correct_{}'.format(q), pyro.distributions.Bernoulli(is_correct_probs[p,q]))\n",
    "    \n",
    "    return skills, has_skills, is_correct"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 180,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'p' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-180-9c3d3a57c27d>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mskills\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhas_skills\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mis_correct\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mext_model\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m<ipython-input-179-723b402ef8c4>\u001b[0m in \u001b[0;36mext_model\u001b[0;34m()\u001b[0m\n\u001b[1;32m     19\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     20\u001b[0m             \u001b[0;31m# skills_needed AND skills\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 21\u001b[0;31m             \u001b[0mhas_skills\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mp\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mq\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtorch\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtensor\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtorch\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mallclose\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mskills_needed\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mq\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mskills\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mp\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfloat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     22\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     23\u001b[0m             \u001b[0;31m# probabilites\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'p' is not defined"
     ]
    }
   ],
   "source": [
    "skills, has_skills, is_correct = ext_model()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(torch.Size([7, 22]), torch.Size([22, 48]), torch.Size([22, 48]))"
      ]
     },
     "execution_count": 173,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "skills.shape, has_skills.shape, is_correct.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ext_guide():\n",
    "    alpha0 = pyro.param('alpha0', torch.ones(questions_n) * 4)\n",
    "    beta0 = pyro.param('beta0', torch.ones(questions_n) * 4)\n",
    "    \n",
    "    guess_probs = pyro.sample('guess_probs', pyro.distributions.Beta(alpha0, beta0).independent(1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 176,
   "metadata": {},
   "outputs": [],
   "source": [
    "# From https://forum.pyro.ai/t/model-based-machine-learning-book-chapter-2-skills-example-in-pyro-tensor-dimension-issue/464/3\n",
    "\n",
    "def reset_torch_seed(seed=99):\n",
    "    torch.manual_seed(seed);\n",
    "\n",
    "def SVI(vals, model, guide, lr = 0.0001, steps=1000):\n",
    "    reset_torch_seed()\n",
    "    pyro.clear_param_store()\n",
    "    svi = pyro.infer.SVI(model = model,\n",
    "                        guide = guide,\n",
    "                        optim = pyro.optim.SGD({\"lr\":lr}),\n",
    "                        #optim = pyro.optim.Adam({\"lr\":lr}), \n",
    "                        loss = pyro.infer.Trace_ELBO())\n",
    "    \n",
    "    losses = []\n",
    "    for t in range(steps):\n",
    "        losses.append(svi.step(*vals)) # vals will be passed to model and guide\n",
    "        \n",
    "    plt.plot(losses)\n",
    "    plt.title(\"ELBO\")\n",
    "    plt.xlabel(\"step\")\n",
    "    plt.ylabel(\"loss\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'data' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-177-7d5d1097983f>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mSVI\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpyro\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcondition\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mext_model\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mext_guide\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0.001\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msteps\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1000\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m: name 'data' is not defined"
     ]
    }
   ],
   "source": [
    "SVI([], pyro.condition(ext_model, data=data), ext_guide, lr = 0.001, steps=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 178,
   "metadata": {},
   "outputs": [
    {
     "ename": "RuntimeError",
     "evalue": "Multiple sample sites named 'people'\n   Trace Shapes:          \n    Param Sites:          \n   Sample Sites:          \nguess_probs dist      | 48\n           value      | 48\n     people dist      |   \n           value   22 |   \n     skills dist      |   \n           value    7 |   \n      skill dist 7 22 |   \n           value 7 22 |   ",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mRuntimeError\u001b[0m                              Traceback (most recent call last)",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pyro/poutine/trace_messenger.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m    146\u001b[0m             \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 147\u001b[0;31m                 \u001b[0mret\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    148\u001b[0m             \u001b[0;32mexcept\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mValueError\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mRuntimeError\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pyro/poutine/messenger.py\u001b[0m in \u001b[0;36m_wraps\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m     26\u001b[0m             \u001b[0;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 27\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0mfn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     28\u001b[0m         \u001b[0m_wraps\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmsngr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-171-f7dfb8fa7de7>\u001b[0m in \u001b[0;36mext_model\u001b[0;34m()\u001b[0m\n\u001b[1;32m     18\u001b[0m     \u001b[0;31m# declare dimensions as conditionally independent\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 19\u001b[0;31m     \u001b[0;32mfor\u001b[0m \u001b[0mp\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mpyro\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'people'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpeople_n\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     20\u001b[0m         \u001b[0;32mfor\u001b[0m \u001b[0mq\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mpyro\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'questions'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mquestions_n\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pyro/poutine/subsample_messenger.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, name, size, subsample_size, subsample, dim, use_cuda, device)\u001b[0m\n\u001b[1;32m     75\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msize\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msubsample_size\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 76\u001b[0;31m             self._indices, self.use_cuda, self.device)\n\u001b[0m\u001b[1;32m     77\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pyro/poutine/subsample_messenger.py\u001b[0m in \u001b[0;36m_subsample\u001b[0;34m(name, size, subsample_size, subsample, use_cuda, device)\u001b[0m\n\u001b[1;32m    104\u001b[0m             }\n\u001b[0;32m--> 105\u001b[0;31m             \u001b[0mapply_stack\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmsg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    106\u001b[0m             \u001b[0msubsample\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmsg\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"value\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pyro/poutine/runtime.py\u001b[0m in \u001b[0;36mapply_stack\u001b[0;34m(initial_msg)\u001b[0m\n\u001b[1;32m    197\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mframe\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mstack\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mpointer\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m  \u001b[0;31m# reversed(stack[0:pointer])\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 198\u001b[0;31m         \u001b[0mframe\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_postprocess_message\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmsg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    199\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pyro/poutine/messenger.py\u001b[0m in \u001b[0;36m_postprocess_message\u001b[0;34m(self, msg)\u001b[0m\n\u001b[1;32m    124\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mhasattr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmethod_name\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 125\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mgetattr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmethod_name\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmsg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    126\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pyro/poutine/trace_messenger.py\u001b[0m in \u001b[0;36m_pyro_post_sample\u001b[0;34m(self, msg)\u001b[0m\n\u001b[1;32m    108\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mparam_only\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 109\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtrace\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madd_node\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmsg\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"name\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mmsg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcopy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    110\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pyro/poutine/trace_struct.py\u001b[0m in \u001b[0;36madd_node\u001b[0;34m(self, site_name, *args, **kwargs)\u001b[0m\n\u001b[1;32m    102\u001b[0m                 \u001b[0;31m# Cannot sample after a previous sample statement.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 103\u001b[0;31m                 \u001b[0;32mraise\u001b[0m \u001b[0mRuntimeError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Multiple {} sites named '{}'\"\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'type'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msite_name\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    104\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mRuntimeError\u001b[0m: Multiple sample sites named 'people'",
      "\nDuring handling of the above exception, another exception occurred:\n",
      "\u001b[0;31mRuntimeError\u001b[0m                              Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-178-47da56238e01>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mSVI\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mext_model\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mext_guide\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0.001\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msteps\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1000\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m<ipython-input-176-867adea04ac7>\u001b[0m in \u001b[0;36mSVI\u001b[0;34m(vals, model, guide, lr, steps)\u001b[0m\n\u001b[1;32m     15\u001b[0m     \u001b[0mlosses\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     16\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mt\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msteps\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 17\u001b[0;31m         \u001b[0mlosses\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msvi\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstep\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mvals\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;31m# vals will be passed to model and guide\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     18\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     19\u001b[0m     \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlosses\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pyro/infer/svi.py\u001b[0m in \u001b[0;36mstep\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m     97\u001b[0m         \u001b[0;31m# get loss and compute gradients\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     98\u001b[0m         \u001b[0;32mwith\u001b[0m \u001b[0mpoutine\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtrace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mparam_only\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mparam_capture\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 99\u001b[0;31m             \u001b[0mloss\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mloss_and_grads\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mguide\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    100\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    101\u001b[0m         params = set(site[\"value\"].unconstrained()\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pyro/infer/trace_elbo.py\u001b[0m in \u001b[0;36mloss_and_grads\u001b[0;34m(self, model, guide, *args, **kwargs)\u001b[0m\n\u001b[1;32m    123\u001b[0m         \u001b[0mloss\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0.0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    124\u001b[0m         \u001b[0;31m# grab a trace from the generator\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 125\u001b[0;31m         \u001b[0;32mfor\u001b[0m \u001b[0mmodel_trace\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mguide_trace\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_get_traces\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mguide\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    126\u001b[0m             \u001b[0mloss_particle\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msurrogate_loss_particle\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_differentiable_loss_particle\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmodel_trace\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mguide_trace\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    127\u001b[0m             \u001b[0mloss\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0mloss_particle\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnum_particles\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pyro/infer/elbo.py\u001b[0m in \u001b[0;36m_get_traces\u001b[0;34m(self, model, guide, *args, **kwargs)\u001b[0m\n\u001b[1;32m    161\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    162\u001b[0m             \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnum_particles\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 163\u001b[0;31m                 \u001b[0;32myield\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_get_trace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mguide\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pyro/infer/trace_elbo.py\u001b[0m in \u001b[0;36m_get_trace\u001b[0;34m(self, model, guide, *args, **kwargs)\u001b[0m\n\u001b[1;32m     50\u001b[0m         \"\"\"\n\u001b[1;32m     51\u001b[0m         model_trace, guide_trace = get_importance_trace(\n\u001b[0;32m---> 52\u001b[0;31m             \"flat\", self.max_plate_nesting, model, guide, *args, **kwargs)\n\u001b[0m\u001b[1;32m     53\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mis_validation_enabled\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     54\u001b[0m             \u001b[0mcheck_if_enumerated\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mguide_trace\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pyro/infer/enum.py\u001b[0m in \u001b[0;36mget_importance_trace\u001b[0;34m(graph_type, max_plate_nesting, model, guide, *args, **kwargs)\u001b[0m\n\u001b[1;32m     42\u001b[0m     \u001b[0mguide_trace\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpoutine\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtrace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mguide\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgraph_type\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mgraph_type\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_trace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     43\u001b[0m     model_trace = poutine.trace(poutine.replay(model, trace=guide_trace),\n\u001b[0;32m---> 44\u001b[0;31m                                 graph_type=graph_type).get_trace(*args, **kwargs)\n\u001b[0m\u001b[1;32m     45\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mis_validation_enabled\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     46\u001b[0m         \u001b[0mcheck_model_guide_match\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmodel_trace\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mguide_trace\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmax_plate_nesting\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pyro/poutine/trace_messenger.py\u001b[0m in \u001b[0;36mget_trace\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m    167\u001b[0m         \u001b[0mCalls\u001b[0m \u001b[0mthis\u001b[0m \u001b[0mpoutine\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mreturns\u001b[0m \u001b[0mits\u001b[0m \u001b[0mtrace\u001b[0m \u001b[0minstead\u001b[0m \u001b[0mof\u001b[0m \u001b[0mthe\u001b[0m \u001b[0mfunction\u001b[0m\u001b[0;31m'\u001b[0m\u001b[0ms\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    168\u001b[0m         \"\"\"\n\u001b[0;32m--> 169\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    170\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmsngr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_trace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pyro/poutine/trace_messenger.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m    151\u001b[0m                 six.reraise(exc_type,\n\u001b[1;32m    152\u001b[0m                             \u001b[0mexc_type\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mu\"{}\\n{}\"\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexc_value\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mshapes\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 153\u001b[0;31m                             traceback)\n\u001b[0m\u001b[1;32m    154\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmsngr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtrace\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madd_node\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"_RETURN\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"_RETURN\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"return\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mret\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    155\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mret\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/six.py\u001b[0m in \u001b[0;36mreraise\u001b[0;34m(tp, value, tb)\u001b[0m\n\u001b[1;32m    690\u001b[0m                 \u001b[0mvalue\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtp\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    691\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__traceback__\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mtb\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 692\u001b[0;31m                 \u001b[0;32mraise\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwith_traceback\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtb\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    693\u001b[0m             \u001b[0;32mraise\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    694\u001b[0m         \u001b[0;32mfinally\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pyro/poutine/trace_messenger.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m    145\u001b[0m                                       args=args, kwargs=kwargs)\n\u001b[1;32m    146\u001b[0m             \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 147\u001b[0;31m                 \u001b[0mret\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    148\u001b[0m             \u001b[0;32mexcept\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mValueError\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mRuntimeError\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    149\u001b[0m                 \u001b[0mexc_type\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mexc_value\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtraceback\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msys\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexc_info\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pyro/poutine/messenger.py\u001b[0m in \u001b[0;36m_wraps\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m     25\u001b[0m         \u001b[0;32mdef\u001b[0m \u001b[0m_wraps\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     26\u001b[0m             \u001b[0;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 27\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0mfn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     28\u001b[0m         \u001b[0m_wraps\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmsngr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     29\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0m_wraps\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-171-f7dfb8fa7de7>\u001b[0m in \u001b[0;36mext_model\u001b[0;34m()\u001b[0m\n\u001b[1;32m     17\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     18\u001b[0m     \u001b[0;31m# declare dimensions as conditionally independent\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 19\u001b[0;31m     \u001b[0;32mfor\u001b[0m \u001b[0mp\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mpyro\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'people'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpeople_n\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     20\u001b[0m         \u001b[0;32mfor\u001b[0m \u001b[0mq\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mpyro\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'questions'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mquestions_n\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     21\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pyro/poutine/subsample_messenger.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, name, size, subsample_size, subsample, dim, use_cuda, device)\u001b[0m\n\u001b[1;32m     74\u001b[0m         self.size, self.subsample_size, self._indices = self._subsample(\n\u001b[1;32m     75\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msize\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msubsample_size\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 76\u001b[0;31m             self._indices, self.use_cuda, self.device)\n\u001b[0m\u001b[1;32m     77\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     78\u001b[0m     \u001b[0;34m@\u001b[0m\u001b[0mstaticmethod\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pyro/poutine/subsample_messenger.py\u001b[0m in \u001b[0;36m_subsample\u001b[0;34m(name, size, subsample_size, subsample, use_cuda, device)\u001b[0m\n\u001b[1;32m    103\u001b[0m                 \u001b[0;34m\"continuation\"\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    104\u001b[0m             }\n\u001b[0;32m--> 105\u001b[0;31m             \u001b[0mapply_stack\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmsg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    106\u001b[0m             \u001b[0msubsample\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmsg\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"value\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    107\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pyro/poutine/runtime.py\u001b[0m in \u001b[0;36mapply_stack\u001b[0;34m(initial_msg)\u001b[0m\n\u001b[1;32m    196\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    197\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mframe\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mstack\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mpointer\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m  \u001b[0;31m# reversed(stack[0:pointer])\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 198\u001b[0;31m         \u001b[0mframe\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_postprocess_message\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmsg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    199\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    200\u001b[0m     \u001b[0mcont\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmsg\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"continuation\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pyro/poutine/messenger.py\u001b[0m in \u001b[0;36m_postprocess_message\u001b[0;34m(self, msg)\u001b[0m\n\u001b[1;32m    123\u001b[0m         \u001b[0mmethod_name\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m\"_pyro_post_{}\"\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmsg\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"type\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    124\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mhasattr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmethod_name\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 125\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mgetattr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmethod_name\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmsg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    126\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    127\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pyro/poutine/trace_messenger.py\u001b[0m in \u001b[0;36m_pyro_post_sample\u001b[0;34m(self, msg)\u001b[0m\n\u001b[1;32m    107\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_pyro_post_sample\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmsg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    108\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mparam_only\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 109\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtrace\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madd_node\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmsg\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"name\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mmsg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcopy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    110\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    111\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_pyro_post_param\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmsg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pyro/poutine/trace_struct.py\u001b[0m in \u001b[0;36madd_node\u001b[0;34m(self, site_name, *args, **kwargs)\u001b[0m\n\u001b[1;32m    101\u001b[0m             \u001b[0;32melif\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'type'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0;34m\"param\"\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    102\u001b[0m                 \u001b[0;31m# Cannot sample after a previous sample statement.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 103\u001b[0;31m                 \u001b[0;32mraise\u001b[0m \u001b[0mRuntimeError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Multiple {} sites named '{}'\"\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'type'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msite_name\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    104\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    105\u001b[0m         \u001b[0;31m# XXX should copy in case site gets mutated, or dont bother?\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mRuntimeError\u001b[0m: Multiple sample sites named 'people'\n   Trace Shapes:          \n    Param Sites:          \n   Sample Sites:          \nguess_probs dist      | 48\n           value      | 48\n     people dist      |   \n           value   22 |   \n     skills dist      |   \n           value    7 |   \n      skill dist 7 22 |   \n           value 7 22 |   "
     ]
    }
   ],
   "source": [
    "SVI([], ext_model, ext_guide, lr = 0.001, steps=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = {'isCorrect'+str(index):torch.tensor(row.values.astype(float)).float() for index, row in responses.iterrows() }"
   ]
  },
  {
   "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}