aflaxman/2019_01_23c_swc-style_bayesian_computation-no_solutions.ipynb

## 2019_01_23c_swc-style_bayesian_computation-no_solutions.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np, matplotlib.pyplot as plt, pandas as pd\n",
    "pd.set_option('display.max_rows', 8)\n",
    "!date\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prerequsites:\n",
    "\n",
    "* SWC Intro to Python\n",
    "* SWC Software Testing\n",
    "\n",
    "* Basic Matplotlib Graphics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Minimal, familiar example for IHME workers new to Bayesian methods\n",
    "\n",
    "Logistic regression... who has coded their own logistic regression before? Who has used logistic regression? Who has heard of logistic regression?\n",
    "\n",
    "[poll results go here]\n",
    "\n",
    "The _Unifying Political Methodology_ version of logistic regression looks like this:\n",
    "\n",
    "\\begin{align}\n",
    "Y_i &\\sim \\text{Bernoulli}(\\pi_t),\\\\\n",
    "\\pi_i &= \\text{logit}(\\beta_0 + \\beta_1 \\cdot X_i).\n",
    "\\end{align}\n",
    "\n",
    "Gary King calls the first line is the _stochastic component_ and the second line the _systematic component_.\n",
    "\n",
    "## Exercise 1: Scientific Python Warm up\n",
    "I always need to look up the formula for $\\text{logit}$.  Let's do that now, and take a look at it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def logit(x):\n",
    "    # fill this in here\n",
    "    \n",
    "logit(0)  # what should this be?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# some minimal testing of our new logit function\n",
    "\n",
    "# write tests here, and if you do them first, you are doing test-driven development"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot this\n",
    "xx = np.linspace(-10, 10)\n",
    "yy = logit(xx)\n",
    "plt.plot(...)  # fill in these details\n",
    "plt.xlabel(...)  # label your plots---you might not remember what this was next time you look at it\n",
    "plt.ylabel(...)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2: \"forward\" simulation\n",
    "\n",
    "Is it helpful to have the bubble diagram of the logistic regression model?  To be fully Bayesian, we should probably go beyond Gary's formulation, and write down priors explicitly.\n",
    "\n",
    "\\begin{align}\n",
    "Y_i &\\sim \\text{Bernoulli}(\\pi_t),\\\\\n",
    "\\pi_i &= \\text{logit}(\\beta_0 + \\beta_1 \\cdot X_i),\\\\\n",
    "\\beta_0, \\beta_1 &\\sim \\text{Uninformative}.\n",
    "\\end{align}\n",
    "\n",
    "(This is an \"improper prior\"; you decide if you consider such a model fully Bayesian.)\n",
    "\n",
    "How to make a graphical representation of this?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# if necessary, from command line do\n",
    "# pip install daft\n",
    "import daft"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAr0AAADFCAYAAABZykpdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAHYhJREFUeJzt3XuU3VV5//H3c2YmN8iFS8LFBJIUMFCQQFBRwQqkahWBequKQDEgSimiVKBSJV5iBRelLYoiYRWBFLT+sFBMhSJeEGq5RCJUKEJCCAgkQLiEXGYy3+f3xzkiCMmcmXPLnPN+rcXKmbO+372fk8Mfnzyzv3tHZiJJkiS1s1KrC5AkSZIazdArSZKktmfolSRJUtsz9EqSJKntGXolSZLU9gy9kiRJanuGXkmSJLU9Q68kSZLanqFXkiRJbc/QK0mSpLZn6JUkSVLb6251AcNNRASwFTAa2AA8nZnrW1tVfUVEN7A1MAJYDzyZmUVrq5IkSRq6yMxW17DZi4idgWOANwCzgJHA85T/0TAWuBe4HbgGWJiZ/S0qdUgqQf4A4P3AfsBrgDVAHzCK8m8EfgncClyemXe1qFRJkqQhcXnDJkTE6yLiGmARsA3wTWBfYEJm7piZkyh3RD8G/Ao4E3ggIj4dESNaVXe1IqIUEX8JLAYuAh4BzgB2yMyJlc+4NbAr8BXKIfiHEfGziHhHq+qWJEkaLDu9ryAiRgFfAI4GzqLc3Xy+yntnAXOBnYG/zMxFjaqzFhExFbgY2BL4DHBjVvE/Q0T0AIdTDsG/AE7OzKcaV6kkSVLt7PT+gYjYAfgfYDrwmsy8sNrAC5CZdwCHAV+l3BU9rjGVDl1E/ClwG3A98KbM/FE1gRcgM/sy83vA3sCTwOKI2LNx1UqSJNXOTu+LRMT2wM+AS4F51QbBTYy3G+Vg+ZXM/GYdSqxZRLwNuAx4T2beVIfxPgicB8zOzLtrHU+SJKkRDL0VlV/b3wz8IDM/X8dxpwM3Acdm5vX1GneItcyo1HJEZt5cx3E/CJwD7O1SB0mStDky9FZExJnAm4G319rhfYWx/5Ty+tm9MvOZeo49iBq6gJ9TXp/89QaMfz4wPjOPrvfYkiRJtTL08pIO6KzMfKhBc1wIFJn58UaMX8X8p1B+AO2QRuy5GxFbUN7B4qTM/M96jy9JklQLQy8QEd8EfpuZX2jgHFsDS4DdMnNFo+bZyNzdwIPAoZl5ZwPn+QvgY5l5UKPmkCRJGoqO370hIsYDf0F5n9qGqax1vQr4SCPn2Yh3AcsaGXgrvg/MiIg9GjyPJEnSoHR86KX8K/8fZ+aj1VwcEV+IiLsi4r6I+Ogg57oIOGrQFdbuw1QZ6mv5fJnZC1xSmU+SJGmzYeiF11F+wGtAle2+9gFmAu8BjhjkXLcD0yJi7CDvq1VVn7EOn4/KPK8dwn2SJEkNY+iFWZTDaDUOo9zJ7AFOAv7fYCbKzD7gLsqhsikiYhLlU9ceqOLymj5fxR3ArIiIIdwrSZLUEIbe8slrv6ny2lnAWMonkR0AXDGE+e6rzNks04HfVLkNW82fLzMfA7qB8YO9V5IkqVEMvTASWDfQRRFRAiZn5iXAtpQ7mp+KiC0i4tsRcVFEHFnFfOsqczZLrZ9vekRcHBHfG8Sczf6MkiRJm2TohV5gRBXXvZpKRzgz11I+va0LeDfwvcw8nvLygIGMqMzZLL1UF0Bf8fNl5pLMnDPIOZv9GSVJkjbJ0AvLgWlVXLcPMDIiuiJiJPAh4N+ByZUxAPqrGGf6i65vhuWUH54baI3txj7foETENpT/v2rJyXOSJEmvxNBb/jX+flVcNxMYTfmBsJuBb2fmYuBhysEXBvj7rBwFPBNYNORqB+8RIIEpA1y3sc83WPsCv2zEqW+SJElD1d3qAjYDtwN/AnxtgOv2AY7KzLv/4P2rgK9FxDuB/xhgjD2BxzNz1ZAqHYLMzIi4HXgDsKkjll/x81U6t/OAfSLibzPz7weYcn+q3w1DkiSpKTr+GOLKll7/B0zfVBiNiOXAtMzcUMNc/wSszswzhzrGEOc9GvhQZr59E9fU4/OVgPuBD2TmrUMdR5Ikqd46PvQCRMQC4PbMPK+Bc2xBudO6T2ZuquPaiLlHVeZ+Y2be38B53gF8AXhtlVukSZIkNYVresvOAz5d+VV+o3wG+K9mB16AzFwHXAB8pVGHRkRED+VlEOcaeCVJ0ubGTm9FRJwHTMrMavbaHezY+wE/APauHN7QdJVu7yLg85n5nQaMfxbl444PNfRKkqTNjaG3IiLGAHcC52Tm/DqOO5HybghzM/Nf6zXuEGt5HeWH7d6SmffUcdxDKJ/etm9mPlyvcSVJkurF5Q0VmbkGOBT4QkQcU48xI2I74L+A77Y68AJUHi77NHB9ROxejzEj4iDKgfd9Bl5JkrS5MvS+SGbeBxxMOfj+Q6X7OyQR8RbgF8DVwGfrU2HtMvNS4EzgpxHxwaGu8a0cYnEq8B3g/Zn503rWKUmSVE8ub3gFEbEtcD4wC/gEcF21hy1ExI6UQ+XhwAmZ+YOGFVqDiHgtcAnl7drOqAT+au99HeWH//qAOZn5QEOKlCRJqhND7yZExJ8Dc4ExwIXAj4D/zczeP7huR+C1wIeBPwO2ALZu5iEUQ1F5uG1t5b+bgX+h3J1e+uKH0Sonye0GHAgcD2wLnANc6MlrkiRpODD0DqDy6/83AB+p/DmN8lG9q4Ee4FWVP++g/JDYlcBKYF5m/l0raq5WZT3ujcAfUz6R7f2Uu9tjgKVAL+Wjif8IWAHcBlwG/DAz+1tRsyRJ0lAYegepcsjErpSDYR/lMPjQH3RGfwi8DShtztt3RUQCZGb8wfvbA5OBEcB6YMnm3rWWJEnaFENvA0TESGAdm3G390Vd3l0beUqbJEnS5sDQ2yCbe7d3Y11eSZKkduSWZY1zeOXPL7a0ildQ6fJCeZmGJElS27PT20Cba7fXLq8kSeo0dnoba7Pr9trllSRJnchOb4Ntbt1eu7ySJKkT2eltvM2m22uXV5IkdSo7vU2wuXR77fJKkqROZae3OVre7bXLK0mSOpmd3iZpdbfXLq8kSepkdnqbp2XdXru8kiSp09npbaJWdXvt8kqSpE5np7e5mt7ttcsrSZJkp7cmEfEgsHOr61BNlmXm1FYXIUmSGsvQW4OIyMEuGYiIkcA6YF5m/l1jKnthroOAG4FdM/P+Rs41XA3lO5QkScOPobcGQw1MzVrb61regRl6JUnqDK7pbY2Gr+11La8kSdLv2emtQS1dwkZ3e+3yVsdOryRJncFOb+s0rNtrl1eSJOml7PTWoNYuYaO6vXZ5q2enV5KkzmCnt7Xq3u21yytJkvRyht4Wysz1wHXAmRHx84h401DHiojpEXEL5S3KcIsySZKk3zP0tlBElIDJlR/fBJxQw3B/DuxXef1sREyqpTZJkqR2YuhtrS5gJFBUfv6ziBjq+tL3AD2V1/2Ai7UlSZIqDL0tlJl9wNuB5ytvjQF2Gew4ETEKmFX5cS1weGaurEuRkiRJbcDQ22KZ+QDwIcphtQt46xCGOQDopRyeP5eZN9WvQkmSpOGvu9UFCDLz2og4HziN8jKFr8MLa353BWYCWwEBPAfcBfy60ikGeCewJfAfwLnNrb42EbEV5fXIb6T8GTYADwPXAP/dyGOaJUlS53Cf3hrUc4/XiOgCfkp5mcK7gY8ChwBPAr8EVlJepzsB2BvYCfgF8C3gy5T/AbNnZj5Xj3oaKSK2A84tlUrvKYpiVHd3dzFhwoRi5MiRFEXBc889x/PPP9+dmZRKpQeLovh8Zl7SoFrcp1eSpA5g6K1BvQNTRLwf+BdgKfBPwFWZ+eRGrh0L/BnwceANlPf6/fLm3BmNiNER8R+Zecj48eM3HHbYYd3veMc72HLLLV/x+vvuu48FCxYUd955Z2Tmusw8LjP/tc41GXolSeoAht4a1CswRcQ44DzKnd0TgOsHE14jYh/gYuAx4PjMfKTWmuotIo6IiCvHjBnTfeaZZ3bttddeVd+7YcMGvvGNb+T1118fpVLplqIoZmfm2jrVZeiVJKkDGHprUI/AVNlP9zpgEXDKUJcnREQPcCbwEeBtmXlPLXXVU0R8Ejj3oIMO4pRTTolSaWjPTy5ZsoTPfOYz/WvXrn26KIpdM3NVHWoz9EqS1AEMvTWoNTBVHuL6KXA15V0Xav4yIuIYYB7wJ5WdIVoqIj4OXHDcccdx+OGHD3j9QNatW8fHPvax/lWrVj1TFMVOmfn8wHdtsj5DryRJHcDQW4NaAlPlEIp/o7wk4a/ruRY3Iv4aOBZ4/Yt2eGi6iNgFuO/II4+MD3zgA3Ubd926dRx77LH9a9asubW/v/+NtYxl6JUkqTO4T2/rvA/YA/ibBjx89jXKYfqMOo87KKVS6capU6cW9Qy8AKNGjWLevHldRVG8ISKOrOvgkiSpLRl6WyAiRlPeneHYzFxX7/ErIfqjwMkRMbXe41ejso538pe+9KWuRow/ffp03v72t2dEzK/h6GZJktQhDL2t8X5gUWb+T6MmyMyHgcsoh9+mK5VKZxx44IGMHz++YXOccMIJEREjKT+8J0mStFGG3tY4EbigCfN8E5hTCYZNExH7FUUxac6cOQN2YBcuXMgFF/z+r+Kyyy7j3HOrO1Suu7ubWbNmUSqVPjf0aiVJUicw9DZZRGxP+WjhH1Zx7S4RsTIiHoyIOyPiqYh4oLKv74Ay8z7gQWD/mooevE9MnDhxw1ZbbTXghQcffDC33XYbq1ev5tZbb+X222/npJNOqnqio446Koqi2MklDpIkaVMMvc03C7gjM/sHujAz7wd+DhyVmTOBXwFHZOazg5jvtsqcTVMqlV6/2267dVdz7ahRo3jzm9/MZZddxre+9S3OOOMMRo6svjE9bdo0IiKBtwyxXEmS1AEMvc23L3DHIK7/Y+DuyusZwP8BRMQXq7z/jsqcVO6LiJhSOSHty5XtzeoqM3fae++9q75+9uzZLFy4kOOPP54ddthh0PONHTu2Hzhs0DdKkqSOUVU3TnU1Ebi/mgsruzyMysxVETEFeDIzeytLJKr97n4LzIiILwNvBl4D9AC9wJZAqXJIRj2N2Gmnnaq++Morr2T8+PH097+0+X355Zfz4Q9/eMD7x40bF88+++yOg65SkiR1DENv83UDAy5tqNgD+N1xwru/6PU+wJ1VjrEBmA6cTrmz3wuMAEa96JrPVzlW1Xp6eqq67vvf/z59fX2cfvrpLFiwgDe+sXzWxKpVq14Wgjemq6sLyp9JkiTpFbm8ofnWAFtUee2LlzasBfaNiBnATKoPvVsAt1T+fD3wSeBKYCnl8H1XZkY9/4uI/meeeWbAwhYvXswNN9zAKaecwl577cXatWtZsmQJAA888ADTp0+v6gOuX78+gdVV/n1IkqQOZOhtvv+jHGYHlJmXZuapldc3Zeb0zLwX2AX4TZXz7Qncm5nrMvPWzLwgMz+YmdMpL2940xA+wyaVSqVnfvWrX23ymhUrVnD++edzxhlnMGbMGADe9a53cfXVVwOwdOlSpk2bVtV8Tz75ZAC/qKloSZLU1lze0Hx3AJ+oZYDMnDOIy2cB/7aRcdYBdT8RriiKu++5554DgY1uIzZp0iTmz5//kvdmz57N7NmzAXj00UfZcceBl+muW7eOvr6+LuB7NRUtSZLamp3e5rsbmBwRr2r0RJVDKQ6gyV3QzLxu6dKlRS1jnHzyyZRKA//ved111xERvZn5eC3zSZKk9mbobbLM7AUWAMc3Ybp3U16z+2AT5nqxr/b29saPf/zjhk901VVX9WfmtQ2fSJIkDWuG3tb4BnB8RIwa8MohqpxQdhLNOe74JTKzLyJuXLBgQbW7VAzJ/fffz1NPPdUFfKqR80iSpOHP0NsCmXk38N/A3AZOcwwwBrimgXNsVGae8Pjjj5euvbYxTdiiKJg7d25/qVT678xc1pBJJElS2zD0ts5fAcdGxOvrPXBETAbOAf4yMzfUe/xqZOYS4JyLLrooq9m+bLAuvPDCfPbZZ4uiKN5W98ElSVLbMfS2SOXBq48D/xYRU+s1bkSMp9zd/YfMXFyvcYciM88AHjzxxBP7162r3yYR1113HQsXLozMPDYzn6vbwJIkqW0ZelsoM68CvgL8tHLoRE0iYhLwI+BnwNm1jlcPRVHstXr16lVz5syp6sCKgVx99dV87WtfA/hcZi6oeUBJktQRDL0tlpkXAGcBP4+IkyJiSN9JRLwbWAz8APhkZmYdyxyyzHy+KIqpq1evfvjoo4/OhQsXDmmcNWvWcOqpp/bPnz8/gdMy84v1rVSSJLWz2Eyy0bAUEZmZGz2AYZBjvRr4F6AAvgpcm5mb3P2gskPDgZR3L9gdODYzb6lHPY0QEV8FTt1+++2LY445puuAAw4Y8J5nn32Wb3/72/zoRz/KzFxRFMUhmfm/daypbt+hJEnafBl6a1DvwBQRXcAHgROBVwH/TvkEt0XAE0AC44GZlE9aeyfQRXlbsvmZubZetTRKROwSEd/KzLeMGDGi2GWXXUp77rlnzJw5k3HjxtHX18fy5ctZtGgR9913X/9vf/vbrlKptKooin8EvljvDrahV5KkzmDorUEjA1NEzARmA/tRDrmvAnqAR4C7gNspr929aXNZyjAYldPi/jYi/rRUKu3R398/jsqxxRGxoVQqPdLf3387cEFmNuyUC0OvJEmdwdBbg2YGpoi4F3i1Aa2+DL2SJHUGH2STJElS2zP0SpIkqe0ZeiVJktT2DL2SJElqe4ZeSZIktT1DryRJktqeoVeSJEltz9ArSZKktmfolSRJUtsz9EqSJKntGXolSZLU9gy9kiRJanuGXkmSJLU9Q68kSZLanqFXkiRJbc/QK0mSpLZn6JUkSVLbM/RKkiSp7Rl6JUmS1PYMvcNARPwT8OrK62tbXI4kSdKwY+gdHvZ60euZLatCkiRpmDL0Dg8/BTZUXt/eykIkSZKGI0Pv8HAb8DywnnIAliRJ0iAYeoeHO4BRlEPvHS2uRZIkadiJzGx1DcNWRGRmRpPmegoYD0zIzOeaMWcnaOZ3KEmSWsdO7/CxGHjUwCtJkjR4dnpr0MguYURMBI4CDunq6noNsE1mdkXEU/39/fcAPwEuz8wljZi/U9jplSSpMxh6a9CIwBQRh5VKpb8vimKPkSNH9u+www7stttuXdtuuy0AzzzzDL/5zW+Khx9+uFizZk13qVR6uCiKecCF6Zc5aIZeSZI6g6G3BvUMTBGxXalU+klRFDP22GOPYs6cOaXddtttk/c89thjzJ8/P2+77TaAJ4qiOCQz76pHPZ3C0CtJUmcw9NagXoEpIuYA39puu+1y3rx5Xdttt92g7l+zZg1nnXVWce+99wZwTmaeUWtNncLQK0lSZzD01qAegSkiPgWc+973vpdjjjmmpnpuuOEG/vmf/zmBS4qi+EhNg3UIQ68kSZ3B0FuDWgNTRBwFXHrcccdx+OGH16WmRYsWMXfu3MzMr2bm6XUZtI0ZeiVJ6gyG3hrUEpgiYquIWHHooYd2f/SjH61rXTfeeCPnnXdeAvtk5uK6Dt5mDL2SJHUG9+ltkVKpdMOECROi3oEX4OCDD2b33XfPUql0Q0QY6CRJUscz9LZARMwqimLfefPmdTVqjrlz55aAbYC/atQckiRJw4WhtwUi4rwpU6b0T5kypWFzjBkzhv33359SqeRODpIkqeMZepssIroy84Cjjz56wC7v8uXLmTNnDkVRAFAUBZ/97Ge58cYbq5rruOOOi6IoXhURM2qrWpIkaXgz9DbfWyOC/ffff8ALp0yZwuTJk6kcPsGll17K5MmTOfjgg6uaaOLEiWyxxRYbgA/XUrAkSdJwZ+htvkPHjx/fX+3Fhx9+OAsXLuTmm2/mnnvuYc6cOYOabMqUKV0RceCgq5QkSWojht4mi4j9dtppp6ofYNt333158sknufTSSzn99NPp7u4e1HwzZsyIUqm0+6ALlSRJaiOG3iaLiLFjx44d1DZiM2bM4IgjjmDrrbd+4b3LL7+8qnsnTJhAZo4aXJWSJEntxdDbfIM+DWT58uVMmzbthZ9XrVpFf391KySKoiAiPIFEkiR1NENvkxVFserpp58eVAh96KGH2HnnnV/4+YEHHmD69OlV3fvUU08BrB3MfJIkSe3G0Nt8//PQQw9V/SDbypUr2WKLLRg9evQL7y1duvQlnd9Nuffee4v+/v67B1+mJElS+zD0Nt81q1ev7vrd3rsDmThxIvPnz3/Je48++ig77rhjVfc/8sgjBfDjwRYpSZLUTgy9zfczoKj2gIlXcvLJJ1MqDfzVLVu2jLVr13YDlw55MkmSpDZg6G2yzEzguiuuuGJDo+eaP39+USqVHsjM5Y2eS5IkaXNm6G2BzPzEihUrun/96183bI5Vq1axePHiKIricw2bRJIkaZiIcuNRQxERmZmD2nP3d0ql0nWjR48+5IorruiqZqnCYJ144on9jzzyyEP9/f3VbfPQoWr5DiVJ0vBhp7dFMvOwtWvX9p599tl1/1fHd7/7XZYvXx5FURxU77ElSZKGI0Nvi2Tm+sw87JZbbuHiiy+u27jXX389l112GcAnM3NZ3QaWJEkaxlzeUIN6/Go8Iv4CuOJNb3oTp512WtSy1OHyyy/nO9/5DsBZmfmFWurqFC5vkCSpMxh6a1CvwBQRb42Iq8eMGdNz1llnde2+++6Duv+xxx7jzDPP7F+xYgXASZn5zVpr6hSGXkmSOoOhtwb1DEwRMTIirs3M2dtvv33/hz70oa6DDtr0ktxf/vKXXHLJJf1LlizpKpVKvy6K4uDMfLwe9XQKQ68kSZ3B0FuDRgSmiJgREedl5lsjIsaNG9c/derUrgkTJkREsHr1apYuXbph1apVpaIoIiJuz8xTM/OmetbRKQy9kiR1BkNvDRoZmCKiCzgEeFdEvD4iJkREqSiK5zJzEbAQuDYz1zdi/k5h6JUkqTMYemtgYBr+/A4lSeoMblkmSZKktmfolSRJUtsz9EqSJKntGXolSZLU9gy9kiRJanvdrS5guCv1jHwsN/Ru1+o6miG6Rzxe9K3fvtV1SJIkDZZbltUgIhJg59OvbXUpTbHs7ENpt+293LJMkqTO4PIGSZIktT1DryRJktqeoVeSJEltz9ArSZKktmfobZB1D9/D0zctGNK9Tyz8R5affyS/vfjEOlf1Us/d+Z8sO/tQ+p5Y/sJ7j1z0MTY883hD55UkSWo2Q2+DjJq8OxMOPHJI926512wmve/zda7o5XpXPEjPpOmseeA2AHJDH8Wap+kaN6nhc0uSJDWTobdBVv7737Pu4f8d0r2jpuxJ1+ixda7o5fpWPsj4/d/L2iXl0Nv7xDJ6tplChDt4SZKk9uLhFA3Su3IZIyZOfeHnxxacRtG79mXXbXXQHEZPndnEyn6v78nljN7l9Txz85UU65+nb+Uyeibu3JJaJEmSGsnQ2wC5oReKfkojt3jhve2PPKchcz38jY+w3fs/T882U1j1k0uI7h4mHHAkvY8vYeU15/Cq47/5ivdteHYlpdFjKfWMZNTUmaxdsojelUsZMXFaQ+qUJElqJUNvA/Q+8RA92055yXuN6vRu94F5dI/bFoBxr/vzF97v2WYKk943d+M1rnyQEduWu7qj/2g/nv/1T+hfvYoxu+4/5FokSZI2V4beBuhb+SA9f9AxrUen9/ErP8M27/wU3WO3feG9nq12eOF115jxL7yO7h56Jmy/0XvLNU4FYNSUvXjquq9TbFj/kiUZkiRJ7cIH2Rqgd+WDNYXHldecw2OX/Q19Tz3Cw18/hucWX09mQd+qRymNGvwDbq9074trjO4eeiZOJUo9lEZtOeS6JUmSNld2ehtg64OPq+n+iYed9rL3elc+yJjd3kipZ+Sgx+t74qGX3TvxXZ9+yTWT3vPZwRcqSZI0TERmtrqGYSsiEmDn069tdSlNsezsQ8nMttrPLCKy3T6TJEl6OZc3SJIkqe0ZeiVJktT2DL2SJElqez7IVptlQEcdYfa7dcxtZFmrC5AkSY1n6K1BZk5twxC4ST70JUmShiN3b6hRqWfkY7mhd7tW19EM0T3i8aJv/fYDXylJkrR5MfRKkiSp7fkgmyRJktqeoVeSJEltz9ArSZKktmfolSRJUtsz9EqSJKntGXolSZLU9gy9kiRJanuGXkmSJLU9Q68kSZLanqFXkiRJbc/QK0mSpLZn6JUkSVLb+//8trp3zm3pbQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 680.315x175.748 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Instantiate the PGM.\n",
    "pgm = daft.PGM([12, 3.1], origin=[0.3, 1.2])\n",
    "\n",
    "# Add in the nodes.\n",
    "pgm.add_node(daft.Node(\"beta0\", r\"$\\beta_0$\", 1.5, 4))\n",
    "pgm.add_node(daft.Node(\"beta1\", r\"$\\beta_1$\", 2.5, 4))\n",
    "pgm.add_node(daft.Node(\"pi_i\", r\"$\\pi_i$\", 2, 3))\n",
    "pgm.add_node(daft.Node(\"X_i\", r\"$X_i$\", 3, 3, observed=True))\n",
    "pgm.add_node(daft.Node(\"Y_i\", r\"$Y_i$\", 2, 2, observed=True))\n",
    "\n",
    "# Add in the edges.\n",
    "pgm.add_edge(\"beta0\", \"pi_i\")\n",
    "pgm.add_edge(\"beta1\", \"pi_i\")\n",
    "pgm.add_edge(\"X_i\", \"pi_i\")\n",
    "pgm.add_edge(\"pi_i\", \"Y_i\")\n",
    "\n",
    "# And a plate.\n",
    "pgm.add_plate(daft.Plate([1.5, 1.4, 2, 2.1], label=r\"$i = 1, \\cdots, N$\",\n",
    "    shift=-0.1))\n",
    "\n",
    "# Render\n",
    "pgm.render();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This sort of figure is sometimes called a \"DAG\" because it is a directed, acyclic graph, and one way to start figuring out what it means is to use it to simulate data starting from the \"source\" nodes and working your way in.  Try it!\n",
    "\n",
    "What are the source nodes here?\n",
    "\n",
    "(Answer: $\\beta_0$, $\\beta_1$ and $X_i$ for $i=1, \\ldots, N$.)\n",
    "\n",
    "We are stopped before we start in this process, because the source nodes are $\\beta_0$ and $\\beta_1$ and above we gave them an improper prior distribution.  Let's change that.  What prior do you want to use?\n",
    "\n",
    "(One answer: $\\beta_0, \\beta_1 \\sim \\text{Normal}(0,1^2)$.)\n",
    "\n",
    "Now can you sample from that?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate a sample value for beta_0 and beta_1\n",
    "beta_0 = ...\n",
    "beta_1 = ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now generate $\\pi_i$ with your $\\beta_0$ and $\\beta_1$ values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this is a \"deterministic\" variable in the PyMC nomenclature\n",
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wait, something is missing here.  We need a value for $X_i$ first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# assign something to be X_i\n",
    "X_i = ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# another solution --- helpful or distracting?\n",
    "X = np.array([3,1,4,1,5,9])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now we can generate a sample of pi_i\n",
    "pi_i = ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# another solution - helpful or distracting\n",
    "pi = logit(beta_0 + beta_1 * X)\n",
    "pi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To finish up this data generating simulation, we can now generate $Y_i$ from $\\pi_i$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y_i = ...\n",
    "Y_i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# another soluation --- helpful or distracting?\n",
    "Y = np.random.binomial(n=1, p=pi)\n",
    "Y"
   ]
  },
  {
   "attachments": {
    "image.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAr0AAADFCAYAAABZykpdAAAgAElEQVR4Ae3dC5hVZb3H8d9aM8NNbipXFUQyAo8UqSctL6XYzUytx65kZKAYmVl2xPIYVFJHO2alZSqeTOFAl2NpHk6Yt+poJYaiJRwUkJvchBEcGWBmr3We/7BmHIaZ2Wvvddlr7/19n8eHfXmvn3c/j39e3vW+EgkBBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAgdIIOKVptqxbNbODJfWW1CzpFUl7ynpEB3a+VtIhknoEY9smyTswG58ggAACCCCAAALlIUDQG26ejpQ0WdLbJR0vqaek1yRZcNhP0nJJT0q6T9JCSblw1WYml/0OTpH0UUknSHqzpF2SmiT1kuRKekrSE5LmSno2Mz2nIwgggAACCCCAAAKRBd4WBLK20vkDSedKOkJS+78s2IrviZK+IOkvkl6U9C/BKmnkDiRcgQWzn5H0TBC4XyXpnZL6d2h3sKT3SrpW0gZJf5R0Voc8vEUAAQQQQAABBBAoMwFb3bxe0iZJ0yQdVED/bSX4t0EgeVwB5dLOOkrSQ5L+Kmlih0C+u77USTpf0gvBqq9tgyAhgAACCCCAAAIIlJnAcElLJf1K0pAi+24rwRdI2iJpapF1JFns3ZK2SpoRbNEopi37i4Ctfq+TdGwxFVAGAQQQQAABBBBAoDQCwyStkPSvBax8dtfTMcF2h0u6y5Tyd7ZNwYLxU2Nq9xPBijiBb0ygVIMAAggggAACCCQpYP9sbw9qzYy5kdHBPtj3xFxvMdWNDVZ4Ty6mcDdlLPC1FV+2OnSDxFcIIIAAAggggEAWBK6WtCimFd6O47HtBGslDej4RYrvayT9WdLnE2rzJkl3JVQ31SKAAAIIIIAAAgjEINC6Ajoyhrq6quJWSbd09WUKn18u6ZHg+LEkmrM9vislvT+JyqkTAQQQQAABBBBAILrATyR9PXo13dZg//RvF1kU+3Bct5Xn+dLOE14vaUKefFG//lgQWEeth/IIIIAAAggggAACMQvYloN6SXZqQ9LpPyTZWbhppw9JeiyFRu0Gt42SjkmhLZpAAAEEEEAAAQQQKEDg05LuKSD/N4MbyeyUh4sLKGdZ7Ua3fxRYJo7s/xVcQhGmrijjs/q/I+nbYRoiDwIIIIAAAggggEB6AjdL+nLI5uy4L7t4wh4KGx9cORyyaEs2OyHCrve1q4vTTHaywtEhGow6PmviA5J+H6ItsiCAAAIIIIAAAqkJ2DW01Z7sBrUnQyKcI+lOSRa8XirJVlALSU3BKnHSe2vb98n2EPcNHjJr/3lnr6OOz+r8myQzbX9Vc2dt8RkCCCCAAAIIIJCaAEGvZOfoPh9S3II5W6XdJukUSfNDlmufzbZFWJtppdbx+SEajGN8dnWzPThXyuPZQgyVLAgggAACCCBQTQIEvVJPSbtDTLpZHRGs9A4KVjRtW4Qd1fUzSbdLmhSiHmvL2kwrRR2fBc13BNcyh+1z2mMM2y/yIYAAAggggECVChD0Snsl2akD+dKb2q0INwanIdje3g8HAeFFkmx7QL5kbVmbaSVrK0yQ3dX4VkmaUmBn0x5jgd0jOwIIIIAAAghUmwBB777rc48KMfFvDYJHC3QtiPykpN8Eq7/2oJilXPBnd3/Yymlr/u7yxfWdtWXjy7fHtqvxFdqPQ4MLMHYUWpD8CCCAAAIIIIBAUgIEvfsevDohBLA9fNY7eCDMzry1LQ1Lg0sfbNuDpXyeFjBbPUuC/Gn8sUGS7ecdkaexrsaXp9gBXx8n6SlJ3gHf8AECCCCAAAIIIIBAyQTsrN15IVq3Y7iO7SSf7en9aXDFcL49vW+R9EIndST90X9LstvSuktdjc9Wbu3GOrti+KvdVRB8d42kG0LkIwsCCCCAAAIIIIBAigJ2pJfdyHZwnjZtm4CdShAl/UDS7CgVFFnWLuD4XZ6ycYzPVrptD/Db8rTF1wgggAACCCCAAAIlELCV3i8l3K6tCNtRZyMTbqez6ntJ2hLygorOyof97KzgzON8+4fD1kc+BBBAAAEEEEAAgRgFbE/vS5Lsn/KTSrbCuyCpykPUOys4ZSKpgNQu7LC9vJ8I0ReyIIAAAggggAACCJRI4MaQe3uL6Z4F1ZslDSumcExlbLX3uRB7e4ttbqYk2zucVFBdbL8ohwACCCCAAAIIINBOoI8kuy1tarvP4ng5OKjXjjgrdbK9thZ8j4u5IxOD7ROtp1jEXD3VIYAAAggggAACCMQpMCbY5jA5pkqHSnpa0rUx1RdHNfZQmz20Flfge3oQ8L4zjs5RBwIIIIAAAggggEA6AmMlrZH0PUm2+ltsepek1ZK+kcF/8rfA1x5ss/23xW5HsDOHrwjqsbGSEEAAAQQQQAABBMpMYJCk+cG2hPeHuHSi/fAOk/Sj4NKKD7T/ImOv/1nSPyTdI8lWuAtJtk3CLuh4VNIbCilIXgQQQAABBBBAAIHsCXwouHXteUlfkWRX9fbopJsW6J4r6ZeSGoIb0PKd+9tJNal/ZA+32W1tuyTZ5RS279iuSe64+mururYdwi7yWBysYH+uwL8MpD44GkQAAQQQQAABBBAIL2AB4DskzQlWRi1AfFbSn4MzaTdKelnSIkmXSrJVYgsks7SPt6vR2n5c6+sxkuw2uXuDFert2nc9s43R9iS/GtzIZkeu2eq1BcEkBBBAAAEEEECgbAQ6ruiVTcdL2FG7ZOKNwX7fpmBP69ogeGztlt1+9t5gJdSCyqym1r51/B3Y0Wp2EoOtau8JblmzW+tICCCAAAIIIIAAAgi0CfQsg9Xe1lXeo9t6zQsEEEAAAQQQQAABBAoUsNVeW0ntuIpaYDWJZbe+ta70JtYIFSOAAAIIIIAAAghUtkCWV3tZ5a3s3x6jQwABBBBAAAEEUhXI6movq7yp/gxoDAEEEEAAAQQQqGyBLK72sspb2b85RocAAggggAACCJREIGurvazyluRnQKMIIIAAAggggEBlC2RptZdV3sr+rTE6BBBAAAEEEECgpAJZWe1llbekPwMaRwABBBBAAAEEKlsgC6u9rPJW9m+M0SGAAAIIIIAAApkQKPVqL6u8mfgZ0AkEEEAAAQQQQKCyBUq52ssqb2X/thgdAggggAACCCCQKYFSrfayypupnwGdQQABBBBAAAEEKlugFKu9rPJW9m+K0SGAAAIIIIBACAEnRB6ydC3woqQju/6ab8pAYI2kUWXQT7qIAAIIIIAAAgiUTMC2DRSa0lztZZU3/+wUM4f5ayUHAggggAACCCBQQQLFBkxp7e1lL2/+H1uxc5i/ZnIggAACCCCAAAIVIlBswJTGai+rvOF+ZMXOYbjayYUAAggggAACCFSAQJSAKenVXutblP5VwPSEGgJGoZjIhAACCCCAAALVLBAlYEpytZdV3vC/yihzGL4VciKAAAIIIIAAAmUsEDVgSmq11/oVtW9lPC0FdR2ngrjIjAACCCCAAALVKBA1YEpitZdV3sJ+iVHnsLDWyI0AAggggAACCJShQBwBU+tq7/9KOjmCwWhJjwcrvHH0K0JXyqooVmU1XXQWAQQQQAABBEohEDVgciX9vV2geleEQVwhaW9Q1w5JQyLUVU1Fo85hNVkxVgQQQAABBBCoUoGoAVOdpOcl5YJgdaukYm/Ja7/Ku13S4Cqdk0KHHXUOC22P/AgggAACCCCAQNkJxBEwvUHSziDofU3SG4tQ6CVpT1DHLkmnFlFHtRaJYw6r1Y5xI4AAAggggECVCMQVMJ0tyYLV3ZI+X4TdmZJeldQg6StFlK/mInHNYTUbMnYEEEAAAQQyL1DsP6VnfmApddACprgMr5N0paRHJJ0R9N/2/NrK7wRJBwdtWXD7rKTnJDUF+W6UdLmk30o6N1jxTYkgcjM2rg9JeoekvpKaJa2XdJ+kP6cwljjnMDIGFSCAAAIIIIAAAlkUiHOVsEaSneDQKOn9kn4dbHtYLekeSbdK+omkBZKWSbKtEA9J+piklZLWSOqXRaRO+jRU0lzXdW2sfm1tbW7QoEFNhx9+eNPw4cOb+vbt2+Q4jtn6ruva+D/TSR1xfRTnHMbVJ+pBAAEEEEAAgZgF4lqljLlbZVNd3KuEH5X0U0kW6P0gCHa3daFhAa4Fx5+T9HZJ35L07RRWRrvoTqiPezuO81vf9ycOGDCg+Zxzzqk966yz1LevLfAemFasWKF58+Z5Tz/9tOP7/m7f96dK+s8Dc0b6JO45jNQZCiOAAAIIIIBAMgIEvdFc4wqY+kuyLQoTJU2T9ECBwetbJd0haZOkiyRtiDasREqf5zjOgj59+tReffXVNePHjw/dSHNzs2655Rb/gQcecGpqah7P5XK2h9lWieNIcc1hHH2hDgQQQAABBBBISICgNxpsHAGTnae7SNKSYF+u7dktJtnxZ1dL+qyk9wZbIIqpJ4kyX5J0w+mnn67LL7/ccV3bqlx4WrVqlb72ta/lGhsbX/E8z/Y61xdeywEl4pjDAyrlAwQQQAABBBDIlgBBb7T5iBow2UNcf5B0r6SvF7i621XPJ0uaLemdwV7frvKl9bltv/jx1KlTde659oxdtLR7925dcsklufr6+h2e540M9jZHqTTqHEZpm7IIIIAAAgggkJIAQW806CgBk9n/MtiS8IWYAt7W0Vh9F0o6sd0JD63fpfnn0ZJWTJo0yfn4xz8eW7sW+F544YW24vtELpezUx+ipChzGKVdyiKAAAIIIIBAigLF/Ttzih2s4KY+IumY4FxdC7ziTDcHwfRVcVZaaF2u6z48atQoL86A1/rQq1cvzZ49uyaXy9kDfJMK7Rf5EUAAAQQQQKD6BAh6SzPnvYPTGWw11i6kiDtZEH2xpMskjYq78pD12T7eI6699lo7ii32NHr0aL3vfe/zHceZE+NZybH3kwoRQAABBBBAIBsCBL2lmQc7msweXPtrgs3bBQ93B8Fvgs10XrXruledeuqpGjBgQOcZYvh02rRpjuM4PYOH92KokSoQQAABBBBAoFIFCHpLM7PT7eGuFJq2yyymSLLAMM10gud5Q6ZMmZJ3z/jChQv14x+/TnH33XfrhhtuCNXX2tpaHX/88XJd1x4CJCGAAAIIIIAAAl0KEPR2SZPYF8OCq4V/F6IFexBsq6QXJT0taXtwIoOd6xsmrQjKnhQmc4x5vjh48ODmgw+2wym6T2eccYYWL16shoYGPfHEE3ryySd16aWXdl+o3bcXXHCBE5zikDfAbleMlwgggAACCCBQZQIEvelP+PGS/iYpF6LpF4KriS+QNEHSM5LOC64nDlG8JctiSdZmasl13RPHjBlTG6ZBeyjttNNOk63w3nbbbbrqqqvUs2f4hemjjjpKwZXF7wrTHnkQQAABBBBAoDoFCHrTn/fjgqA3bMv/JOnvQeaxkv4veG3XDodJFmBbm63JVkRHBMGzXVtsx5vFmnzfH/mWt7wldJ1nnnmmbJvDRRddpOHDh4cu15qxX79+9heIc1rf8ycCCCCAAAIIINBRINRqXMdCvI8kMFiSreCGSXbKQ6/g5jELVLdJ2ivJtkiEnbuXJFmwbAHuaZLeLMlub7N6+kqyv/jk34cQprdBHt/3e4wcafdGhEsLFixoeeAtl9t/8Xvu3Ln61Kc+lbeS/v37Ozt37jwsb0YyIIAAAggggEDVCoQNnKoWKIGBm/n+0V3Xjdg5vsuCr8e1e/3WYI9v1yVf/6ZZ0mhJM4IA14LdHkEw3ZrrG60v4vqzrs7i6vzp17/+tZqamjRjxgzNmzdP73jHvrsm6uvr1TEI7qq2mpqWU9FsTCQEEEAAAQQQQKBTAbY3dMqS6Ie7JB0UsoX2Wxsag20Ktmpr+3vtwbYwydp6PGjTbmiz83MXSFodBN/PBufc2raHWP5zHCe3Y8eOvH1bunSpHnzwQV1++eUaP368GhsbtWrVqpZyK1eulJ3FGybt2bPHziVuCJOXPAgggAACCCBQnQIEvenPu+3JtWA2TLpL0hVBxj8FK7bLJdmpDs+HqUDSsZKsjF2C8URwVNongrpse8PJIesJnc113R3PPGPP3HWdtmzZoptuuqnlwbU+ffq0ZPzgBz+oe++9t+X16tWrZQ+phUnbtm2zYP0vYfKSBwEEEEAAAQSqU4DtDenPuz1Y9sWIzdrZu2GTndzwyy4yWyAc+41wnuf9fdmyZad2d1PakCFDNGeOXab2erIH2uw/Sxs3btRhh+Xfprt7927bHmH7G371ek28QgABBBBAAAEE9hdgpXd/jzTe2UkMR0g6PIXG7OyvU9JeBfV9f9Hq1au9KOO77LLL7NKJvFUsWrTIjiyzfcqb82YmAwIIIIAAAghUrUD+qKJqaRIbuAVo8yRdlFgLr1f8YUm2Z9cut0gzfXfv3r3OI488knib99xzT873/fsTb4gGEEAAAQQQQKCsBQh6SzN9twRBrx1HllSyfa52tdnrd/wm1dKB9TY5jvPwvHnzwp5ScWANIT554YUXtH37dtva8OUQ2cmCAAIIIIAAAlUsQNBbmsm3LQ5/ljQrweYnS7InxO5LsI0uq/Z9f9rmzZvd++9PZhHW8zzNmjUrV1NTY45ruuwIXyCAAAIIIIAAAsG5rUCURuDzki6UZMeIxZ1sz/D1kj4jyc7pLUWys8euv/322/0wx5cV2sFbb73V37lzp5fL5d5baFnyI4AAAggggED1CbDSW7o5twevPhecrDAqxm4MCFZ3vydpaYz1FlPVVbafePr06Tk7ZSGuZA+vLVy40PF93/7S8Gpc9VIPAggggAACCFSuAEFvaef2Hkn/JukPwVXBUXszRNJDkv4o6bqolcVR3vO88Q0NDfVTpkwJdWFFvjbtHN+bb77Zsn09eCAwXxG+RwABBBBAAAEEWm7ggqF4AbsJzB4Yi5psG8K/B3t87cGzYo77spMafiTptqAe61tW0kGu6/5D0shp06Y5Z511VsH92rVrl6655prcihUr7C9qdqXydwuupPMCcc1h57XzKQIIIIAAAghkQiCOgC0TAylRJ+IMmN4k6adBwGsBnT0Blu/0A5s/uwTCTi8YF+wRtiuHs5psXFcMGzbMmzx5cs0pp9gRwt2nnTt36mc/+5keeugh3/f9LZ7nTZRkAXRcKc45jKtP1IMAAggggAACMQsQ9EYDjTtgsuO37Irg6cHlFb+RZDe4LZH0siRrz/bsTpBkN619QJKVsdVhu96sMdpwUil9tOM4t/m+/64ePXp4Rx99tHvsscc6EyZMUP/+/e12Na1bt05LlizRihUrci+99FKN67r1nud9X9K3AoM4Oxr3HMbZN+pCAAEEEEAAgZgECHqjQSYZMFlga3fynhAEuXaDW52kDcGFE08Ge3f/lEAgGE0lXGm7Le6rjuO823XdY3K5XP/WrSKO4zS7rrshl8vZGC2gT/KWiyTnMJwEuRBAAAEEEEAAgYwLWMCUVlpepsFtWj7FtpPmHBbbR8ohgAACCCCAQEQBTm+ICEhxBBBAAAEEEEAAgewLEPRmf47oIQIIIIAAAggggEBEAYLeiIAURwABBBBAAAEEEMi+AEFv9ueIHiKAAAIIIIAAAghEFCDojQhIcQQQQAABBBBAAIHsCxD0Zn+O6CECCCCAAAIIIIBARAGC3oiAFEcAAQQQQAABBBDIvgBBb/bniB4igAACCCCAAAIIRBQg6I0ISHEEEEAAAQQQQACB7AsQ9GZ/jughAggggAACCCCAQEQBgt6IgBRHAAEEEEAAAQQQyL4AQW/254geIoAAAggggAACCEQUIOiNCEhxBBBAAAEEEEAAgewLEPRmf47oIQIIIIAAAggggEBEAYLeiIAURwABBBBAAAEEEMi+AEFv9ueIHiKAAAIIIIAAAghEFCDojQhIcQQQQAABBBBAAIHsCxD0Zn+O6CECCCCAAAIIIIBARAGC3oiAFEcAAQQQQAABBBBAoNIF/JQG+ANJ1pb9d39KbVZLM2nNYbV4Mk4EEEAAAQQyKcBKbyan5YBOjW/3yYR2r3mJAAIIIIAAAgggEEKAoDcEUgay/EFSc9CPJzPQH7qAAAIIIIAAAgiUlQBBb3lM12JJr0naI8kCYBICCCCAAAIIIIAAAqkJpLUfdKik3ZJ2SDottdFVR0NpzWF1aDJKBBBAAAEEEKhIgTQDpu2ScpL6VaRk6QaV5hyWbpS0jAACCCCAQJULsL2hfH4ASyVtlPRq+XSZniKAAAIIIIAAAtkQcLLRjbLtha0SJmU4WNIFkibW1NS8WdKhvu/XOI6zPZfLLZP0qKS5klaVrV42Op7kHGZjhPQCAQQQQAABBBIL2KqFNomA6RzXdb/jed4xPXv2zA0fPlxjxoypGTRoUIvpjh079Pzzz3vr16/3du3aVeu67nrP82ZLujU4x7da7OMaZxJzGFffqAcBBBBAAAEEYhJIapUypu5lvpo4A6ahrus+6nne2GOOOcabMmWKO2bMmG4BNm3apDlz5viLF9vhDnrZ87yJkp7tthBfdhSIcw471s17BBBAAAEEEMiIAEFvtImIK2CaIum2oUOH+rNnz64ZOtQOawifdu3apZkzZ3rLly+3+bxe0lXhS1d9zrjmsOohAUAAAQQQQCDLAgS90WYnjoDpy5JuOP/88zV58uRIvXnwwQf1wx/+0Hcc507P8z4bqbLqKRzHHFaPFiNFAAEEEECgTAUIeqNNXNSAyR5Uu2vq1Kk699xzo/UkKL1kyRLNmjXL933/u5JmxFJpZVcSdQ4rW4fRIYAAAgggUCECBL3RJjJKwHSw4zhbzj777NqLL744Wi86lH744Yd14403Wt/eKsmOOiN1LRBlDruulW8QQAABBBBAIFMCnNNboulwXffBgQMHOnEHvDacM844Q+PGjfOtjQSPVCuRHM0igAACCCCAAAKFCxD0Fm4WR4njPc87zh5ai6OyzuqYNWuWze2hkj7f2fd8hgACCCCAAAIIVJMAQW8JZttxnBtHjBiRGzFiRGKt9+nTRyeddJJc1+Ukh8SUqRgBBBBAAAEEykWAoDf9marxff+UT3/603lXedetW6cpU6bI87yWXtqf11xzjWzPbpg0depUx/O8wyWNDZOfPAgggAACCCCAQKUKEPSmP7PvcRynZRU2X9O2EnzEEUcouHxCd911V8t727MbJg0ePFgHHXRQs6RPhclPHgQQQAABBBBAoFIFCHrTn9mzBwwYkAvbrB1ltnDhQj322GNatmxZy8pv2LKWb8SIETWO45xaSBnyIoAAAggggAAClSZA0JvyjDqOc8LIkSPzbm1o7dZxxx2nbdu2tazyzpgxQ7W1ta1fhfpz7Nixjuu640JlJhMCCCCAAAIIIFChAgS9KU+s4zj9+vXrV9D5yGPHjtV5552nQw45pK23c+fObXvd3YuBAwfK9/1e3eXhOwQQQAABBBBAoNIFCHrTn2G7DKGgZA+0HXXUUW1l6uvrlcuF2yFhD785jlNwm22N8QIBBBBAAAEEEKgAAYLelCfR87z6V155paAgdO3atTryyCPberpy5UqNHj267X13L7Zv325fN3aXh+8QQAABBBBAAIFKFyDoTX+G/7p27dpwy7SStm7daicwqHfv3m09Xb169X4rv21fdPJi+fLlXi6X+3snX/ERAggggAACCCBQNQIEvelP9X0NDQ01rWfv5mvejh2bM2fOftk2btyoww47bL/PunqzYcMGO+T3ka6+53MEEEAAAQQQQKAaBAh605/lP0rywl4w0Vn3LrvsMrtprbOv9vtszZo1amxstOMe7trvC94ggAACCCCAAAJVJpA/cqoykBSGa/t5F82fP98ujUg0zZkzx3Ndd6WkdYk2ROUIIIAAAggggEDGBQh6SzBBvu9/ccuWLbXPPfdcYq3bCQ9Lly61a4i/nlgjVIwAAggggAACCJSJQEHnxZbJmNLspq3aFmXoOM6i3r17T5w/f35NmK0KhQ5q+vTpuQ0bNqz1PC/cMQ+FNlA5+Yuew8ohYCQIIIAAAghUvgArvSWaY9/3z2lsbNx73XXXFXR8WZju/uIXv9C6detslff0MPnJgwACCCCAAAIIVLoAQW/pZniPBb6PP/647rjjjth68cADD+juu++2+r4kaU1sFVMRAggggAACCCBQxgJF/dN8GY837q7H8U/jH5M0/+STT9aVV17pRNnqYFcT//znP7cxzpT0zbgHW6H1xTGHFUrDsBBAAAEEEKgcAYLeaHMZV8D0Hsdx7u3Tp0/dzJkza8aNG1dQrzZt2qSrr746t2XLFit3qaSfFFRBdWeOaw6rW5HRI4AAAgggkHEBgt5oExRnwNTTcZz7fd8/c9iwYblPfvKTNaef3v2W3Keeekp33nlnbtWqVfYw3HOe550haXO0IVVd6TjnsOrwGDACCCCAAALlIkDQG22mkgiYxjqOc6Pv+7b66/Tv3z83atSomoEDB9pbNTQ0aPXq1c319fWu53n22ZO+718h6U/RhlK1pZOYw6rFZOAIIIAAAghkVYCgN9rMJBkw1UiaKOmDjuOc6DiORb0W6L7q+/4SSQsl3S9pT7QhVH3pJOew6nEBQAABBBBAAIHKELCAiVTeAsxhec8fvUcAAQQQQCCUAEeWhWIiEwIIIIAAAggggEA5CxD0lvPs0XcEEEAAAQQQQACBUAIEvaGYyIQAAggggAACCCBQzgIEveU8e/QdAQQQQAABBBBAIJQAQW8oJjIhgAACCCCAAAIIlLMAQW/E2XNq6zZJshMAKv6/YKwRxSiOAAIIIIAAAgikL8A5vdHMW467OnKGHZdb+WnNdWfbIF8TylYAAARuSURBVCvtN2NzWGljqvwfIyNEAAEEEECgQAFWegsEIzsCCCCAAAIIIIBA+QkQ9JbfnNFjBBBAAAEEEEAAgQIFCHoLBCM7AggggAACCCCAQPkJEPSW35zRYwQQQAABBBBAAIECBQh6CwQLm333+mV65U/zwmbfL9/LC7+vdTdN0kt3TN/v87jfvPr0/8geTmt6eV1b1Rtuv0TNOza3vecFAggggAACCCBQCQIEvQnNYq8jxmngqZOKqr3v+DM15CPfKKpsIYX2bnlRdUNGa9fKxS3F/OYmebteUU3/IYVUQ14EEEAAAQQQQCDzAgS9CU3R1t98R7vX/6Oo2nuNOFY1vfsVVbaQQk1bX9SAk85X46p9Qe/el9eo7tARchxO8CrEkbwIIIAAAgggkH2B2ux3sTx7uHfrGvUYPKqt85vmXSlvb2Pb+9YXB58+Rb1HTWh9m+qfTdvWqffRJ2rHYwvk7XlNTVvXqG7wkan2gcYQQAABBBBAAIE0BAh6E1D2m/dKXk5uz4Paah826fq213G+WH/LZzX0o99oWaGtf/ROObV1GnjKJO3dvEpb77teh1/0k06ba965VW7vfnLreqrXqAlqXLVEe7euVo/BR3Wanw8RQAABBBBAAIFyFiDoTWD29r68VnWDRuxXc1IrvUM/Plu1/Qe1tNX/bR9qa9O2KQz5yKy29x1f7N36onoM2req2/sNJ+i15x5VrqFefd54UsesvEcAAQQQQAABBMpegKA3gSm0vbJ1HVZM41jp3bzgazr0A19Wbb99Qa51ve7g4W0jqOkzoO21rfjWDRzW9r5j2X193Lf9oteI8dq+6EfymvfstyWjrTAvEEAAAQQQQACBMhfgQbYEJrBlFbXdft5Cm7BtCZvu/oqatm/Q+h9N1qtLH5Dve2qq3yi3V+EPuHVWtn0fWwLkwaPkuHVye/UttLvkRwABBBBAAAEEMi/AY/rRpsi34kfOuD9aLSFKW5Da8MzvdcjEi0Lk3j9LlLLta7IzfSVV2m/G5rDSxtR+2niNAAIIIIAAAvzPPvJvILWgN3JPY6iAoDcGRKpAAAEEEEAAgZIIsL2hJOw0igACCCCAAAIIIJCmAEFvmtq0hQACCCCAAAIIIFASAYLekrDTKAIIIIAAAggggECaAhxZFk17jT3HFq2Ksivdso+57HrddYdtDkkIIIAAAgggUOECBL3RJtgOuq20IDCfCCcd5BPiewQQQAABBBBAoNIEnNq6TUHga8FvRf8XjLXSppDxIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggUK/D/wpfiZFEhiZoAAAAASUVORK5CYII="
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bayesian inference (inverse probability)\n",
    "\n",
    "Now that we have done the forward simulation, perhaps concepts with cohere better when we do Bayesian inference.  This is sometimes called _the inverse problem_, because instead of starting from the source nodes of the DAG and working our way forward along the arrows, we will start with information on the sink nodes of the DAG (and the other grey shaded nodes) and work our way back.\n",
    "\n",
    "![image.png](attachment:image.png)\n",
    "\n",
    "Because of the semantics of Python, we still write our program from the source nodes forward, starting from the priors.  I'm using PyMC2, because I think it is cleanest, but you should try PyMC3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pymc as pm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# priors\n",
    "\n",
    "beta_0 = pm.Normal(...)\n",
    "beta_1 = pm.Normal(...)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Want to see what the priors look like?  One way to do that is with MCMC sampling (soon to be described in more detail):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pm.MCMC([beta_0, beta_1]).sample(iter=10_000, burn=5_000, thin=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(beta_0.trace(), density=True)\n",
    "# remember to label your plots for later\n",
    "plt.title('Prior Distribution of $\\\\beta_0$');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(beta_1.trace(), density=True)\n",
    "plt.title('Prior Distribution of $\\\\beta_1$');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now for the systematic part:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pi = beta_0 + beta_1 * X  # it is pretty tricky that this works"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# oops but we forgot the logit link\n",
    "pi = ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# failed solution\n",
    "# too tricky... what does that error mean?\n",
    "pi = logit(beta_0 + beta_1 * X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# can you do it a differen way, that is clearer?\n",
    "..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# one solution\n",
    "@pm.deterministic\n",
    "def pi(beta_0=beta_0, beta_1=beta_1, X=X):\n",
    "    return logit(beta_0 + beta_1*X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally the stochastic part (sometimes also called the likelihood):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y_obs = pm.Bernoulli('Y_obs', ...)  # called Y_obs to signify that it is about \"observed\" data, \n",
    "                                    # and because var name Y is already taken by the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ingredients for this model are now all defined.  Want to see some things about them?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y_obs.value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y_obs.logp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y_obs.parents"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y_obs.children"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pi.children"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pi.parents"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pi.value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "beta_0.children"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "beta_1.children"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "beta_0.value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here is how you can fit the model with MCMC:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y_obs.parents"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pi.value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y_obs.value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y_obs.logp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "beta_0.value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "beta_1.value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's refactor all of this into a function, make an `pm.MCMC` object out of it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def model(X, Y):\n",
    "    # prior\n",
    "    beta_0 = pm.Normal('beta_0', mu=0, tau=1)\n",
    "    beta_1 = pm.Normal('beta_1', mu=0, tau=1)\n",
    "    \n",
    "    # systematic component\n",
    "    @pm.deterministic\n",
    "    def pi(beta_0=beta_0, beta_1=beta_1, X=X):\n",
    "        return logit(beta_0 + beta_1*X)\n",
    "    \n",
    "    # stochastic component (aka likelihood)\n",
    "    Y_obs = pm.Bernoulli('Y_obs', pi, value=Y, observed=True)\n",
    "    \n",
    "    return locals() # cool python affordance, see https://stackoverflow.com/questions/7969949/whats-the-difference-between-globals-locals-and-vars\n",
    "model(X, Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = pm.MCMC(model(X,Y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "m.logp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.beta_0.value = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "m.logp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "m.sample(iter=20_000, burn=10_000, thin=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(m.beta_0.trace())\n",
    "plt.title('Posterior Distribution of $\\\\beta_0$');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(m.beta_1.trace())\n",
    "plt.title('Posterior Distribution of $\\\\beta_1$');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Questions?"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (dismod_env)",
   "language": "python",
   "name": "dismod_env"
  },
  "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
}