interventions_and_causal_effects_in_pyro.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interventions and Causal Effects in Pyro"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statistics import mean\n",
    "import torch\n",
    "import numpy as np\n",
    "import pyro\n",
    "import pyro.distributions as dist\n",
    "from pyro.infer import Importance, EmpiricalMarginal\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "pyro.set_rng_seed(101)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To start, we implement the exact same survey model that we trained in the `bnlearn` tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'A': tensor(1), 'S': tensor(1), 'E': tensor(0), 'O': tensor(0), 'R': tensor(1), 'T': tensor(2)}\n"
     ]
    }
   ],
   "source": [
    "A_alias = ['young','adult','old']\n",
    "S_alias = ['M','F']\n",
    "E_alias = ['high','uni']\n",
    "O_alias = ['emp','self']\n",
    "R_alias = ['small','big']\n",
    "T_alias = ['car','train','other']\n",
    "\n",
    "A_prob = torch.tensor([0.3,0.5,0.2])\n",
    "S_prob = torch.tensor([0.6,0.4])\n",
    "E_prob = torch.tensor([[[0.75,0.25], [0.72,0.28], [0.88,0.12]],\n",
    "                     [[0.64,0.36], [0.7,0.3], [0.9,0.1]]])\n",
    "O_prob = torch.tensor([[0.96,0.04], [0.92,0.08]])\n",
    "R_prob = torch.tensor([[0.25,0.75], [0.2,0.8]])\n",
    "T_prob = torch.tensor([[[0.48,0.42,0.1], [0.56,0.36,0.08]],\n",
    "                     [[0.58,0.24,0.18], [0.7,0.21,0.09]]])\n",
    "\n",
    "def model():\n",
    "    A = pyro.sample(\"A\", dist.Categorical(probs=A_prob))\n",
    "    S = pyro.sample(\"S\", dist.Categorical(probs=S_prob))\n",
    "    E = pyro.sample(\"E\", dist.Categorical(probs=E_prob[S][A]))\n",
    "    O = pyro.sample(\"O\", dist.Categorical(probs=O_prob[E]))\n",
    "    R = pyro.sample(\"R\", dist.Categorical(probs=R_prob[E]))\n",
    "    T = pyro.sample(\"T\", dist.Categorical(probs=T_prob[R][O]))\n",
    "    return{'A': A,'S': S,'E': E,'O': O,'R': R,'T': T}\n",
    "\n",
    "print(model())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can condition on someone being self-employed (O = 1) and living in a big city."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "conditioned_model = pyro.condition(model, data={'O':torch.tensor(1),'R':torch.tensor(1)})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this we can predict the person's age.  Here, I use simple importance sampling for the inference algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'P(A | O = self, R = big) - Importance Sampling')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAam0lEQVR4nO3debhddX3v8fdH5goakEghCcRHUYu2IKYKihVBUWgttAKKA4PUtE9tnefair222qsVB1p6qVhwuCJVEaQWRQaHVrBBEVEUUoWSyBAio+MFv/eP9TuLzfEcspOcffZJ8n49z37OWr+19trfPX7W77fW3idVhSRJAA8YdwGSpLnDUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QyFGZbkj5O8Zz2uf1qSY2ewpFmVZP8kKwbmH5Xk8iR3JnnZOGsblKSSPGKaZS9I8vkht/PsJB+f2eo0VyU5IclH2vSuSe5Kstm465pJhsJaSnJtkp+2F8NN7UN827ZsS+DNwDsnXWfbtv6/z8Dtb5Xk7Un+p9VxTZLXJsn6bntEXgdcVFXbVdX71rRy+7D+cXu8ViZ592y/6arqo1V10JDrfgZ4TJLfWtfba6+ht63r9WfS4IfeXJJkvyT/meT2JD9K8h9JfnucNVXV/1TVtlV1zzjrmGmGwrp5dlVtC+wNLKELAoBDge9W1cpJ6z8H+DnwjCS/vp63/a/AgcAhwHbAi4ClwHvXc7ujshvw7bW8zp7t8X0q8FzgxTNe1cz6GN1zsEFLsvm4a5hKkgcB5wLvB3YAFgBvpXtPaYYZCuuhffj/O/DY1nQw8MUpVj0G+CfgCuCF63p7SQ4EDgKeU1VXVtXdVXVJ2+ZLpxsOWcfbSpITk9yc5I4k30ry2LZsqyTvar2Vm5L8U5JtptjGhcDTgJPanv8j16aGqloO/Aew10zcpykckuT7SW5J8s4kDwBIcmySr0yslOSgJN9re6n/mOSLSf5oYDsXA787EwUlWdx6S8cluT7JrUn+JMlvJ7kiyW1JThpY/9i213xSq++77XUysXyXJOe0vevlSV4ysOyEJJ9I8pEkdwB/ArwJeG57vr7Z1jsuyVVtCPD7Sf54YBv7J1mR5NXttXJDkuMGlm+T5O+TXNfq+8rEayXJPm3v/7Yk30yy/zQPyyMBqupjVXVPVf20qj5fVVe07Tw8yYVJVrfn8qNJ5g3UcG3rTV/ReqGnJtkpyb+3+/SFJNtPevyXJvlhuz+vWcNztXmbvzjJ/2rPx51JPp9kx4H1j26Pw+okf9nqevoaXxSzzFBYD0kW0e2xf6M1/SbwvUnr7AbsD3y0XY5ej5t8BnBpVV0/2FhVlwIr6HoQU9V5bnvjTXU5d5rbOgj4Hbo35IOBI4HVbdk7WvtewCPo9tz+avIGquoA4MvAn7Vu9tVrc2eTPBp4CrD8ftZZl/s24Q/oenp70/XyfqVH0t7UnwDeCDyE7vl90qTVrgIWp9ujnSlPBHan6ym9B/gL4OnAY4Ajkzx10rr/DewIvAX4VJId2rIz6F4buwCHA3+b5ICB6x5Kd//mAacCfwt8vD1fe7Z1bgZ+D3gQcBxwYpK9B7bx63SvkQXA8cA/THzIAu8CHk/3mO1AN5z4yyQLgH8D3tbaXwN8Msn8KR6Lq4F7kpye5OCBbU8I8PZ2H38DWAScMGmd59C9fx4JPJtuZ+5NwHy6z8HJx7ueRvf4HwS8fi0+vJ9P9xg9FNiy3S+S7AH8I/ACYGfufbzmnqryshYX4FrgLuA24Dq6J3qbtuwa4FmT1n8zcHmbXgDcAzzufrZ/GnDsNMs+AJwxzbJLgL+Ywft5AN2bcR/gAQPtAX4MPHygbV/gB216f2DFwLKLgT9ai9st4I52G0U3NLPVCJ7HGnyugD8FLmjTxwJfadNHA1+ddP+vH7xPwBZte7uuYy2nAW9r04vbthYMLF8NPHdg/pPAKwZq/SGQgeVfoxtWXNReb9sNLHs7cFqbPgH40qRaTgA+soZ6Pw28fOD5/imw+cDymydeN23ZnlNs4/XAhye1fQ44Zprb/I32OK0A7gbOAXaaZt3DgG8MzF8LvGDS43fywPyfA5+e9Pg/emD5/wZOnfz4DKy7+cBr/c2TXlPntem/Aj42sOzXgF8AT5/p1/b6XuwprJvDqmpeVe1WVX9aVT9t7bfSjfMPOpquh0B1w01fpBtOWhe30O1lTGXntnxGVNWFwEnAPwA3Jzml7QnPp3tBXzaxRw6c19pnyt7AtnR7yU8EHjiD2x402OO6jm5Pc7JdBter7h29YtI6E8/5bZOvnO5Mpruy9ica3DQw/dMp5rcdmF/Z6powcV92AX5UVXdOWja4h3qfXudU2t75JW0I6ja63vGOA6usrqq7B+Z/0urbEdiarhcz2W7AEYM9O2A/pnl9V9VVVXVsVS2kG67dha4HRRsKOiPdiQl3AB+ZVB+s3eMJw702pnLjwPTE4wC/+jr6Cff2vOcUQ2FmXUEb/wRI8iS6Lugbk9yY5Ea6D7nnZ90O6n0BeGIbtuoleSLdXuGFU12pjZ3eNc1l2g+qqnpfVT0e2KPdr9fSBc9Pgce0YJxXVQ+u7sDwjKnOmcBXmWJoan3vWzP4OO5Kt8c92Q3AwoHby+B88xvAtVV1xxT346PVDcVsW1UHr6GedbWg1TVh4r78ENghyXaTlg2eCDH5Z5LvM59kK7o963fR7ZnPAz5L12Nak1uAnwEPn2LZ9XQ9hXkDlwdW1TvWtNGq+i5dr2HiWN7ftrp/s6oeRHeMbX3PxhvmtbE2Jr+OtqEbjpxzDIWZ9Vm6M2YmHAOcT/ehule7PBbYhu6g9Fqpqi8AF9CNvT4myWZJ9qHbMzq5qq6Z5noHD3wwTb5MWUe6A5tPTLIF3VDOz4BfVtUvgX+mG1d+aFt3QZJnDnMf0h0YvXYt7vY7gJdkmrO21uW+DXhtku1byL4cmOr7Bv8G/GaSw1qQv5RuDH3QU+nGqMflocDLkmyR5Ai6kPpsdcee/hN4e5Kt0502ezzd62U6N9EdH5n4bNgS2ApYBdyd5GC6cfY1aq+VDwLvTnfAe7Mk+7ag+Qjw7CTPbO1bpztoPTlwSfLodAeyF7b5RcBRdEOm0PXU7gJub8cqXjtMfWvwl0l+Lclj6I4RrO93UT5Bd3+flO7U9RNY/+AaCUNhZn0GeHR7A2xNd3D2/VV148DlB8CHWfchpOcAF9EN2dxF9+Y6lW5cdCY9iO7D/1a67vNq7v3+xevpDv5e0rrrXwAeNeR2F9GdUTSUqvoW8CVm5o0+2dnAZcDldB/+p05x+7cAR9CNK6+mC/hl3Pd0yKOA/zOC+oZ1KV2P9Bbgb4DDq2piaOIourHvHwJnAW9pOxfT+df2d3WSr7ehp5cBZ9K9Fp5PN54/rNcA3wL+C/gR8Hd0x6iupzvI/Sa6wLme7jme6jPpTroe9qVJfkwXBlcCr27L30o35Hg73fP4qbWobzpfpHuNXwC8q6qG+jLjdKrq23Tv0TPoeg130R17mXOn1ea+Q5FaX0mWAntU1SvW8fqnARdX1WkzWddcke6bwi+vqqvGXcu6aHvQK+gOXF6U5NnAi6rqyDHVcyzdQe/9xnH7G5ski4EfAFtMOk4y07ezLd0xqN3bjuKcMSe/rLIhq6pTxl3DXFZDflN4LmlDY5fSHUt5LV23/xLov9H8mfFVpw1F24G4gO718y66HtS146xpKg4fzT2fphvO0NyxL90ZNLfQneN+2MAZZ9KwDuXeEwB2B55Xc3CoxuEjSVLPnoIkqbdBH1PYcccda/HixeMuQ5I2KJdddtktVTXlF0436FBYvHgxy5YtG3cZkrRBSXLddMscPpIk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9TbobzRr03Li+VePu4SN1iuf8cg1r6RNgj0FSVLPUJAk9QwFSVLPUJAk9QwFSVJvpKGQ5Nok30pyeZJlrW2HJOcnuab93b61J8n7kixPckWSvUdZmyTpV81GT+FpVbVXVS1p828ALqiq3YEL2jzAwXT/zHp3YClw8izUJkkaMI7ho0OB09v06cBhA+0fqs4lwLwkO4+hPknaZI06FAr4fJLLkixtbTtV1Q1t+kZgpza9ALh+4LorWpskaZaM+hvN+1XVyiQPBc5P8t3BhVVVSWptNtjCZSnArrvuOnOVSpJG21OoqpXt783AWcATgJsmhoXa35vb6iuBRQNXX9jaJm/zlKpaUlVL5s+fP8ryJWmTM7JQSPLAJNtNTAMHAVcC5wDHtNWOAc5u0+cAR7ezkPYBbh8YZpIkzYJRDh/tBJyVZOJ2/m9VnZfkv4AzkxwPXAcc2db/LHAIsBz4CXDcCGuTJE1hZKFQVd8H9pyifTVw4BTtBbx0VPVIktbMbzRLknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknojD4UkmyX5RpJz2/zDklyaZHmSjyfZsrVv1eaXt+WLR12bJOm+ZqOn8HLgqoH5vwNOrKpHALcCx7f244FbW/uJbT1J0iwaaSgkWQj8LvCBNh/gAOATbZXTgcPa9KFtnrb8wLa+JGmWjLqn8B7gdcAv2/xDgNuq6u42vwJY0KYXANcDtOW3t/XvI8nSJMuSLFu1atUoa5ekTc7IQiHJ7wE3V9VlM7ndqjqlqpZU1ZL58+fP5KYlaZO3+Qi3/WTg95McAmwNPAh4LzAvyeatN7AQWNnWXwksAlYk2Rx4MLB6hPVJkiYZWU+hqt5YVQurajHwPODCqnoBcBFweFvtGODsNn1Om6ctv7CqalT1SZJ+1Ti+p/B64FVJltMdMzi1tZ8KPKS1vwp4wxhqk6RN2iiHj3pVdTFwcZv+PvCEKdb5GXDEbNQjSZqa32iWJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSb2ShkGTrJF9L8s0k307y1tb+sCSXJlme5ONJtmztW7X55W354lHVJkma2ih7Cj8HDqiqPYG9gGcl2Qf4O+DEqnoEcCtwfFv/eODW1n5iW0+SNIuGCoUkFwzTNqg6d7XZLdqlgAOAT7T204HD2vShbZ62/MAkGaY+SdLMuN9QaENAOwA7Jtk+yQ7tshhYsKaNJ9ksyeXAzcD5wH8Dt1XV3W2VFQPbWQBcD9CW3w48ZIptLk2yLMmyVatWDXMfJUlD2nwNy/8YeAWwC3AZMLHnfgdw0po2XlX3AHslmQecBTx63Uvtt3kKcArAkiVLan23J0m61/2GQlW9F3hvkj+vqvev641U1W1JLgL2BeYl2bz1BhYCK9tqK4FFwIokmwMPBlav621KktbemnoKAFTV+5M8CVg8eJ2q+tB010kyH/h/LRC2AZ5Bd/D4IuBw4AzgGODsdpVz2vxX2/ILq8qegCTNoqFCIcmHgYcDlwP3tOYCpg0FYGfg9CSb0R27OLOqzk3yHeCMJG8DvgGc2tY/FfhwkuXAj4Dnre2dkSStn6FCAVgC7LE2e+5VdQXwuCnavw88YYr2nwFHDLt9SdLMGzYUrgR+HbhhhLXMqhPPv3rcJWy0XvmMR467BEnraNhQ2BH4TpKv0X0pDYCq+v2RVCVJGothQ+GEURYhSZobhj376IujLkSSNH7Dnn10J93ZRgBb0v1kxY+r6kGjKkySNPuG7SlsNzHdfo/oUGCfURUlSRqPtf6V1PZDd58GnjmCeiRJYzTs8NEfDsw+gO57Cz8bSUWSpLEZ9uyjZw9M3w1cSzeEJEnaiAx7TOG4URciSRq/Yf/JzsIkZyW5uV0+mWThqIuTJM2uYQ80/wvdr5ju0i6faW2SpI3IsKEwv6r+parubpfTgPkjrEuSNAbDhsLqJC9s/15zsyQvxH+AI0kbnWFD4cXAkcCNdL+Uejhw7IhqkiSNybCnpP41cExV3QqQZAfgXXRhIUnaSAzbU/itiUAAqKofMcU/0JEkbdiGDYUHJNl+Yqb1FIbtZUiSNhDDfrD/PfDVJP/a5o8A/mY0JUmSxmXYbzR/KMky4IDW9IdV9Z3RlSVJGoehh4BaCBgEkrQRW+ufzpYkbbwMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSb2ShkGRRkouSfCfJt5O8vLXvkOT8JNe0v9u39iR5X5LlSa5IsveoapMkTW2UPYW7gVdX1R7APsBLk+wBvAG4oKp2By5o8wAHA7u3y1Lg5BHWJkmawshCoapuqKqvt+k7gauABcChwOlttdOBw9r0ocCHqnMJMC/JzqOqT5L0q2blmEKSxXT/qe1SYKequqEtuhHYqU0vAK4fuNqK1jZ5W0uTLEuybNWqVSOrWZI2RSMPhSTbAp8EXlFVdwwuq6oCam22V1WnVNWSqloyf/78GaxUkjTSUEiyBV0gfLSqPtWab5oYFmp/b27tK4FFA1df2NokSbNklGcfBTgVuKqq3j2w6BzgmDZ9DHD2QPvR7SykfYDbB4aZJEmzYOj/vLYOngy8CPhWkstb25uAdwBnJjkeuA44si37LHAIsBz4CXDcCGuTJE1hZKFQVV8BMs3iA6dYv4CXjqoeSbPvxPOvHncJG61XPuORI9mu32iWJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSb2ShkOSDSW5OcuVA2w5Jzk9yTfu7fWtPkvclWZ7kiiR7j6ouSdL0RtlTOA141qS2NwAXVNXuwAVtHuBgYPd2WQqcPMK6JEnTGFkoVNWXgB9Naj4UOL1Nnw4cNtD+oepcAsxLsvOoapMkTW22jynsVFU3tOkbgZ3a9ALg+oH1VrS2X5FkaZJlSZatWrVqdJVK0iZobAeaq6qAWofrnVJVS6pqyfz580dQmSRtumY7FG6aGBZqf29u7SuBRQPrLWxtkqRZNNuhcA5wTJs+Bjh7oP3odhbSPsDtA8NMkqRZsvmoNpzkY8D+wI5JVgBvAd4BnJnkeOA64Mi2+meBQ4DlwE+A40ZVlyRpeiMLhao6appFB06xbgEvHVUtkqTh+I1mSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVJvToVCkmcl+V6S5UneMO56JGlTM2dCIclmwD8ABwN7AEcl2WO8VUnSpmXOhALwBGB5VX2/qn4BnAEcOuaaJGmTsvm4CxiwALh+YH4F8MTJKyVZCixts3cl+d4s1DYX7AjcMu4ihvGqcRcwN2wwzxf4nDWb0nO223QL5lIoDKWqTgFOGXcdsy3JsqpaMu46NByfrw2Pz1lnLg0frQQWDcwvbG2SpFkyl0Lhv4DdkzwsyZbA84BzxlyTJG1S5szwUVXdneTPgM8BmwEfrKpvj7msuWSTGzLbwPl8bXh8zoBU1bhrkCTNEXNp+EiSNGaGgiSpZyhIMyTJsUlOWsM6i5Nc2ab3SnLI7FSn+5PkrmnaT0ty+GzXM06GgjQ+ewGGguYUQ2EWJfnrJK8YmP+bJC9P8s4kVyb5VpLntmX7Jzl3YN2Tkhzbpq9N8tYkX2/XeXRrn5/k/CTfTvKBJNcl2XGW7+ZGK8mnk1zWHt+lre24JFcn+Rrw5IF177OHOXlPtJ12/dfAc5NcPvG8a/SSvKq9364cfD+2ZWnvte8l+QLw0DGVOTaGwuz6IHA0QJIH0H0XYwXdHuOewNOBdybZeYht3VJVewMnA69pbW8BLqyqxwCfAHad2fI3eS+uqscDS4CXJVkAvJUuDPaj+yHHobTf9/or4ONVtVdVfXwUBeu+kjweOI7uJ3T2AV6S5HEDq/wB8Ci65/Jo4EmzXuSYGQqzqKquBVa3F+FBwDfoPkw+VlX3VNVNwBeB3x5ic59qfy8DFrfp/eh+SJCqOg+4dcaKF3RB8E3gErpv378IuLiqVrUPeT/Y5779gLOq6sdVdRfd++gpA8t/h3vfjz8ELhxHkeNkKMy+DwDH0u2tfPB+1rub+z4/W09a/vP29x7m0JcQN1ZJ9qfrye1bVXvSBfp37+cq/fPXeoVbjrpGaSYYCrPvLOBZdL2BzwFfphtX3izJfLo9la8B1wF7JNkqyTzgwCG2/R/AkQBJDgK2H0H9m6oHA7dW1U/aMZx9gG2ApyZ5SJItgCMG1r8WeHyb/n1giym2eSew3ehK1hS+DByW5NeSPJBuuOjLA8u/xL3vx52Bp42jyHFyD3OWVdUvklwE3FZV9yQ5C9gX+CZQwOuq6kaAJGcCVwI/oNszXZO3Ah9L8iLgq8CNdB88Wn/nAX+S5Crge3RDSDcAJ9A91rcBlw+s/8/A2W246Tzgx1Ns8yLgDUkuB97ucYXRq6qvJzmNbscL4ANV9Y0kE6ucBRwAfAf4H7rndpPiz1zMsjaU8HXgiKq6Zoa3vRVwT/sdqX2Bk6tqr5m8DUkbN3sKs6j9e9Fz6Q50zWggNLsCZ7bg+QXwkhHchqSNmD0FSVLPA82SpJ6hIEnqGQqSpJ6hIM2wJIclqYnfpJI2JIaCNPOOAr7S/kobFM8+kmZQkm3pvtz2NOAzVfWoMZckrRV7CtLMOhQ4r6qupvvxw8ev6QrSXGIoSDPrKNov1ba/DiFpg+LwkTRDkuxA9/8xVtH9jtVm7e9u5RtNGwh7CtLMORz4cFXtVlWLq2oR3Y8ZPmUN15PmDENBmjlH0f3K5qBP4hCSNiAOH0mSevYUJEk9Q0GS1DMUJEk9Q0GS1DMUJEk9Q0GS1DMUJEm9/w+Lz3/7s/fXlgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "A_posterior = pyro.infer.Importance(conditioned_model, num_samples=5000).run()\n",
    "A_marginal = EmpiricalMarginal(A_posterior,\"A\")\n",
    "A_samples = [A_marginal().item() for _ in range(1000)]\n",
    "A_unique, A_counts = np.unique(A_samples, return_counts=True)\n",
    "\n",
    "plt.bar(A_unique, A_counts, align='center', alpha=0.5)\n",
    "plt.xticks(A_unique, A_alias)\n",
    "plt.ylabel('count')\n",
    "plt.xlabel('A')\n",
    "plt.title('P(A | O = self, R = big) - Importance Sampling')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interventions and causal effect estimation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So our model is working just as it did with `bnlearn`. Just as before, let's do the intervention on Residence (R) and see it's effect on Occupation (O).\n",
    "\n",
    "First we view the distribution of occupation in large cities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "big_city_model = pyro.do(model, data={'R': torch.tensor(1)})\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'P(T | do(R = big)) - Importance Sampling')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAYp0lEQVR4nO3de7hUdb3H8fcHgaOJCsgOL4CQoqbdrH0008o0syzDesosKyx7qOxup7LLKarjeep0UY9dOWpiV9EukmVFJJZmGiQZRKahJMRNhRTzEvY9f/x+WxfD7M2wZc3sze/zep559rqv75q15zNr/daaGUUEZmZWjiGdLsDMzNrLwW9mVhgHv5lZYRz8ZmaFcfCbmRXGwW9mVhgHf4dIerOkc/o573RJ07di+nmS3rQV018r6ZD+1NZfeZu+0cf4xZKOanFZN0g6uNL/FEm/3gZl2iAg6VRJ11T6N0h6QidrGmgc/DWRdLuk+/M/3WpJF0kakccNBz4CfEbSs/M0GyTdJykq/RskTWhz3ScA90bEjbl/uqR/5lrWS/q1pMPbWRNARBwcEfNanPyzwCcq894ErM/b1i9bemNqp8ZgGygkHSzpZ5Luzv8rCyQd3+m6ImJERCztdB0DiYO/XidExAjg6UA3KewBpgB/iogVEfGr/I85Aug5Sh3ZMywi/trmmt8CfL1h2CW5vjHAVcClba5pa80Gnidpj8qwbwJv7lA924ykoZ2uoQ8/BOYAewCPB94J3NPRiqwpB38bRMQK4ErgSXnQi4Cr61qfpGMl/UnS3yV9AVBl3BBJH5G0TNIaSRdL2i2PGw4c3VttEbGRFKB7S+qqofQdJV0i6V5Jv5P01Erdt0t6fu7eSdJMSeskLZH0fknLK3U+ACwAjqssex5wjKR/2xaF5jOz0yXdkuv9pKR98xnRPZJm5ecTSUdJWi7pQ5LuzNtySmVZu+X9sDbvl49IGpLHnZqb3s6WdBdwCfAV4PCes7A83Ysl3ZjXfUe1KVDSxFzvVEl/zTV8uDJ+h1zbX/K2LJA0Po87UNKcfBR/s6STenk+xgCTgP+LiIfy49qIuCaPHyXpiryN63L3uMr88yT9V37+Nkj6oaTdJX0zb9NvJU1seP7fKWlp3p7P9Dxnveyr/XL3RZK+KOlHeVuvl7RvZdoX5O38u6QvSbpaW9FMOlg4+Nsgv4iOB27Mg54M3FzTusYA3yOdXYwB/gIcUZnk1Px4HvAEYATwhTxuMvCviFhOEznIXg/cBazrZZoj82l+b48j+yh/CulsYjTwLeAHkoY1me5jwMRc/7HAa5tMswR45I0jv/n+Ezigj/VvreOAZwDPBN4PzMi1jCe9yb+6Mu0epP2xNzAVmCGpp5bzgN1I2/Nc0nP8hsq8hwFLgbF5+W8BrstnhCPzNPfl+UYCLwbeKunEhnqPJG3/McBHJT0xDz8j13o8sCvwRuAfknYmHcF/i3QEfzLwJUkHNXku7gJuBb4h6URJYxvGDwG+BuwDTADu59H/ux4nA6/Lz9G+wHV5ntGk/fmxhulfRjqTfjrpf+eNTepq5mTg48CoXPNZ8Mhr5zLgg8DupNfos1pc5uASEX7U8ABuBzYA64FlwJeAnfK4W4AXNplnIhDA0C0sezowvZdxrwd+U+kXsBx4U+6fC5xeGX8AKRCHkt4gVjVZ10N5Ox4mvcCPquH5mt5Q9xBgJfDsyvP5/Ny9FDiuMu2bgOUNyzsLuLBh2ArgOY+hvm9U+gM4otK/APhApf9zwDm5+yhgI7BzZfws4D+BHfLze1Bl3JuBebn7VOCvDbWcClyzhXrPAc5u+L8aVxl/A3By7r4ZmNJkGa8CftUw7KvAx3pZ5zhSmP8F+BfwS2ByL9M+DVhX6Z8HfLjh+buy0n8CsLDh+X9hpf90YG6z5ydPu1/uvgg4vzLueFKzK6TXznUNr507yK+d7enhI/56nRgRIyNin4g4PSLuz8PXAbvUtM69SP+sAET6D76jYfyySv8yUuiP7aOuWZGOLMcCi0hHuXWo1v0v0hvWXk2m22QbG7p77EJ6s9rSMLTpBfbFW1Hv6kr3/U36R1T610XEfZX+ZaTtGAMMY/N9snelv9n2bULSYZKuyk0pfyedFYxpmGxVpfsflfrGk8K60T7AYdUzNuAU0tnLZiJieUS8PSL2zfPeB1yc63ucpK/mpqx7SG8KIyXtUFnE1jyfsOnz0vN8tqK356HZa6fp2e9g5+DvjJuA/Wta9krSCxkASar2A38jvSh7TCAdja4mnfZKUjV0HhERdwLTgOmS9mw2TUOINns8u4/aq3UPIR1B/q2XbRxX6R/fZJonAr+vLG9vYDhNmtiicoE9Ig5uHL+NjMpNJz0mkLbtTtIZV+M+WVEtsWFZzb5S91uki9rjI2I30nUANZmumTtITSvNhl+dD156HiMi4q1bWmBE3AF8kUeva72XdHZ5WETsCjwnD2+1xmaq+73n+XwsNvm/yq+dcb1PPng5+Dvjx6S23Dr8CDhY0suV7gB5J5seoX0beI+kSUq3l/436a6djRHxEPDzvmqLiJuBn5LatJuNr4Zos8ev+qj9GZW63w08CPymyXSzgA/mC4Z7A2+vjpS0I+msZE5l8HOBX0TEg32sv24flzQ8v/m9BLg0Ih4mbc9ZknaRtA+pzb2vW0dXA+PyNZceuwB3R8QDkg4FXrMVdZ0PfFLSZCVPkbQ7cAWwv6TXSRqWH/9euTbwiLwvPi5pP6UbCMaQ2tx79t8upKP29ZJGs3l7fX+8L693PPAu0oXvx+JHwJPzNYqhwNvo5exmsHPwd8YPgQMltXpq2rJ8VP5K4FOk9vjJwLWVSS4k3a75S+A24AHgHZXxXyVdYOvLZ4Bpkh6/jcrucTmpXXldruHlEfHPJtN9gnQKfhvpjeoy0ptEjxNIbeTVI8BTSEfBnbKKtF1/I90Z9ZaI+FMe9w5Ss8hS4BrS0fuFfSzrF8BiYJWkO/Ow04FPSLoX+CjpzaRVn8/T/4x0++UFpOtR9wIvIF0M/Vvehk8Dze6Meoh0LeHneRmLSPvk1Dz+HGAn0hnOb4CfbEV9vbmcdG1lISm0L3gsC6u8dv6H9No5CJjPpv9b2wXlixjWZpKmkS7ovbsf804HiIjp27isnuVfC7w98oe4BjpJbyVdqHxu7r8eOC0iFuX+pwBfjYi2f/Asr/8o0oXh7bLZoBMkBenC8a01rmMI6QDjlIi4qq71dMJA/jDIdi0iZnS6ht5ExBFbnqpz8vWFJ5Bu95tMaj9+5NbAiDisOn2kT+52JPRtcJF0HHA9qVnqfaRrEM2aGwc1B//gNK/TBXTYcFKT1CTSXTrfId0ua/ZYHU5qahsO/JF0Z979fc8y+Lipx8ysML64a2ZWmEHR1DNmzJiYOHFip8swMxtUFixYcGdEbPa9WoMi+CdOnMj8+fM7XYaZ2aAiaVmz4W7qMTMrjIPfzKwwDn4zs8I4+M3MCuPgNzMrjIPfzKwwtQa/pJGSLlP6/dclkg6XNFrpNzxvyX9H1VmDmZltqu4j/nOBn0TEgaTfP10CnEn6ibTJpJ8BPLPmGszMrKK24Je0G+lXdi4AiIiHImI96UeRZ+bJZgKNPwhtZmY1qvOTu5OAtcDXJD2V9IMJ7wLGRsTKPM0q0u+4biZ/X/00gAkTJvS7iLPn/Lnf81rf3nNsXb8eaWZ1qrOpZyjwdODLEXEI6ReGNmnWyT9m3PTrQSNiRkR0R0R3V9dmXzVhZmb9VGfwLweWR8T1uf8y0hvB6p4f6s5/19RYg5mZNagt+CNiFXCHpAPyoGNIP2wwG5iah00l/W6mmZm1Sd3fzvkO4JuShpN+SPoNpDebWZJOA5YBJ9Vcg5mZVdQa/BGxEOhuMuqYOtdrZma98yd3zcwK4+A3MyuMg9/MrDAOfjOzwjj4zcwK4+A3MyuMg9/MrDAOfjOzwjj4zcwK4+A3MyuMg9/MrDAOfjOzwjj4zcwK4+A3MyuMg9/MrDAOfjOzwjj4zcwK4+A3MyuMg9/MrDAOfjOzwjj4zcwK4+A3MyuMg9/MrDAOfjOzwjj4zcwKM7TOhUu6HbgXeBjYGBHdkkYDlwATgduBkyJiXZ11mJnZo9pxxP+8iHhaRHTn/jOBuRExGZib+83MrE060dQzBZiZu2cCJ3agBjOzYtUd/AH8TNICSdPysLERsTJ3rwLGNptR0jRJ8yXNX7t2bc1lmpmVo9Y2fuDIiFgh6fHAHEl/qo6MiJAUzWaMiBnADIDu7u6m05iZ2dar9Yg/Ilbkv2uA7wOHAqsl7QmQ/66pswYzM9tUbcEvaWdJu/R0Ay8AFgGzgal5sqnA5XXVYGZmm6uzqWcs8H1JPev5VkT8RNJvgVmSTgOWASfVWIOZmTWoLfgjYinw1CbD7wKOqWu9ZmbWN39y18ysMA5+M7PCOPjNzArj4DczK4yD38ysMA5+M7PCOPjNzArj4DczK4yD38ysMA5+M7PCOPjNzArj4DczK4yD38ysMA5+M7PCOPjNzArj4DczK4yD38ysMA5+M7PCOPjNzArj4DczK4yD38ysMA5+M7PCOPjNzArj4DczK4yD38ysMA5+M7PC1B78knaQdKOkK3L/JEnXS7pV0iWShtddg5mZPaodR/zvApZU+j8NnB0R+wHrgNPaUIOZmWW1Br+kccCLgfNzv4CjgcvyJDOBE+uswczMNlX3Ef85wPuBf+X+3YH1EbEx9y8H9m42o6RpkuZLmr927dqayzQzK0dtwS/pJcCaiFjQn/kjYkZEdEdEd1dX1zauzsysXENrXPYRwEslHQ/sCOwKnAuMlDQ0H/WPA1bUWIOZmTWo7Yg/Ij4YEeMiYiJwMvCLiDgFuAp4RZ5sKnB5XTWYmdnmOnEf/weAMyTdSmrzv6ADNZiZFavOpp5HRMQ8YF7uXgoc2o71mpnZ5vzJXTOzwjj4zcwK4+A3MyuMg9/MrDAOfjOzwjj4zcwK4+A3MyuMg9/MrDAOfjOzwjj4zcwK4+A3MyuMg9/MrDAtBb+kua0MMzOzga/Pb+eUtCPwOGCMpFGA8qhd6eUnE83MbGDb0tcyvxl4N7AXsIBHg/8e4As11mVmZjXpM/gj4lzgXEnviIjz2lSTmZnVqKUfYomI8yQ9C5hYnSciLq6pLjMzq0lLwS/p68C+wELg4Tw4AAe/mdkg0+pPL3YDB0VE1FmMmZnVr9X7+BcBe9RZiJmZtUerR/xjgD9KugF4sGdgRLy0lqrMzKw2rQb/9DqLMDOz9mn1rp6r6y7EzMzao9W7eu4l3cUDMBwYBtwXEbvWVZiZmdWj1SP+XXq6JQmYAjyzrqLMzKw+W/3tnJH8ADiur+kk7SjpBkm/l7RY0sfz8EmSrpd0q6RLJA3vZ+1mZtYPrTb1vLzSO4R0X/8DW5jtQeDoiNggaRhwjaQrgTOAsyPiO5K+ApwGfHnrSzczs/5o9a6eEyrdG4HbSc09vcof9tqQe4flRwBHA6/Jw2eS7hhy8JuZtUmrbfxv6M/CJe1A+lbP/YAvAn8B1kfExjzJcnr5emdJ04BpABMmTOjP6s3MrIlWf4hlnKTvS1qTH9+VNG5L80XEwxHxNGAccChwYKuFRcSMiOiOiO6urq5WZzMzsy1o9eLu14DZpO/l3wv4YR7WkohYD1wFHA6MlNRzpjEOWNFytWZm9pi1GvxdEfG1iNiYHxcBfR6GS+qSNDJ37wQcCywhvQG8Ik82Fbi8X5WbmVm/tBr8d0l6raQd8uO1wF1bmGdP4CpJNwG/BeZExBXAB4AzJN0K7A5c0N/izcxs67V6V88bgfOAs0l35vwaOLWvGSLiJuCQJsOXktr7zcysA1oN/k8AUyNiHYCk0cBnSW8IZmY2iLTa1POUntAHiIi7aXI0b2ZmA1+rwT9E0qiennzE3+rZgpmZDSCthvfngOskXZr7XwmcVU9JZmZWp1Y/uXuxpPmkr1sAeHlE/LG+sszMrC4tN9fkoHfYm5kNclv9tcxmZja4OfjNzArj4DczK4yD38ysMA5+M7PCOPjNzArjT9/agHP2nD93uoTt1nuO3b/TJdgA4CN+M7PCOPjNzArj4DczK4yD38ysMA5+M7PCOPjNzArj4DczK4yD38ysMA5+M7PCOPjNzArj4DczK4yD38ysMA5+M7PC1Bb8ksZLukrSHyUtlvSuPHy0pDmSbsl/R9VVg5mZba7OI/6NwHsj4iDgmcDbJB0EnAnMjYjJwNzcb2ZmbVJb8EfEyoj4Xe6+F1gC7A1MAWbmyWYCJ9ZVg5mZba4tbfySJgKHANcDYyNiZR61ChjbyzzTJM2XNH/t2rXtKNPMrAi1B7+kEcB3gXdHxD3VcRERQDSbLyJmRER3RHR3dXXVXaaZWTFqDX5Jw0ih/82I+F4evFrSnnn8nsCaOmswM7NN1XlXj4ALgCUR8fnKqNnA1Nw9Fbi8rhrMzGxzdf7Y+hHA64A/SFqYh30I+BQwS9JpwDLgpBprMDOzBrUFf0RcA6iX0cfUtV4zM+ubP7lrZlYYB7+ZWWEc/GZmhXHwm5kVxsFvZlYYB7+ZWWEc/GZmhXHwm5kVxsFvZlYYB7+ZWWEc/GZmhanzS9rMrABnz/lzp0vYbr3n2P1rWa6P+M3MCuPgNzMrjIPfzKwwDn4zs8I4+M3MCuPgNzMrjIPfzKwwDn4zs8I4+M3MCuPgNzMrjIPfzKwwDn4zs8I4+M3MClNb8Eu6UNIaSYsqw0ZLmiPplvx3VF3rNzOz5uo84r8IeGHDsDOBuRExGZib+83MrI1qC/6I+CVwd8PgKcDM3D0TOLGu9ZuZWXPtbuMfGxErc/cqYGxvE0qaJmm+pPlr165tT3VmZgXo2MXdiAgg+hg/IyK6I6K7q6urjZWZmW3f2h38qyXtCZD/rmnz+s3Mitfu4J8NTM3dU4HL27x+M7Pi1Xk757eB64ADJC2XdBrwKeBYSbcAz8/9ZmbWRkPrWnBEvLqXUcfUtU4zM9syf3LXzKwwDn4zs8I4+M3MCuPgNzMrjIPfzKwwDn4zs8I4+M3MCuPgNzMrjIPfzKwwDn4zs8I4+M3MCuPgNzMrjIPfzKwwDn4zs8I4+M3MCuPgNzMrjIPfzKwwDn4zs8I4+M3MCuPgNzMrjIPfzKwwDn4zs8I4+M3MCuPgNzMrjIPfzKwwDn4zs8J0JPglvVDSzZJulXRmJ2owMytV24Nf0g7AF4EXAQcBr5Z0ULvrMDMrVSeO+A8Fbo2IpRHxEPAdYEoH6jAzK9LQDqxzb+COSv9y4LDGiSRNA6bl3g2Sbm5DbQPBGODOThfRijM6XcDAMGj2F3ifZYNmn22D/bVPs4GdCP6WRMQMYEan62g3SfMjorvTdVhrvL8GH++zzjT1rADGV/rH5WFmZtYGnQj+3wKTJU2SNBw4GZjdgTrMzIrU9qaeiNgo6e3AT4EdgAsjYnG76xjAimveGuS8vwaf4veZIqLTNZiZWRv5k7tmZoVx8JuZFcbBb9ZA0khJp/djvh9LGllHTda6xv0n6ShJV3SypoHGwT9ISBqwn7nYDo0ENgv+Le2DiDg+ItbXVpW1qun+66/t8bW33W3QYCDp9cB/AAHcBMwCPgIMB+4CTomI1ZKmA/sCTwD+Cry6IwWX51PAvpIWAv8EHgDWAQcC+0v6AemzKDsC5+YPGyLpdqAbGAFcCVwDPIv0OZUpEXF/m7ejCJLOAN6Ye88Hnsmj+28O8CNghKTLgCcBC4DXRkRIegbwedI+uxM4NSJWSpoHLASOBL4NfK6Nm1S/iPCjjQ/gYODPwJjcPxoYxaN3WL0J+Fzunk76J92p03WX9AAmAoty91HAfcCkyvjR+e9OwCJg99x/O+nrACYCG4Gn5eGzSEHT8W3b3h7AM4A/ADuTwnsxcEjP/qvsw7+TPiw6BLiOFOjDgF8DXXm6V5FuLweYB3yp09tX18NH/O13NHBpRNwJEBF3S3oycImkPUlH/bdVpp8dPlLstBsiorpP3inpZbl7PDCZdKZWdVtELMzdC0hvBrbtHQl8PyLuA5D0PeDZTaa7ISKW52kWkvbHetIZwBxJkD5XtLIyzyX1ld1ZDv6B4Tzg8xExW9JRpCP9Hvd1pCKremQf5P3zfODwiPhHbhLYsck8D1a6HyadHVjnNO6PoYCAxRFxeC/zbLevPV/cbb9fAK+UtDuApNHAbjz6fUVTO1WYPeJeYJdexu0GrMuhfyCpPdk651fAiZIeJ2ln4GXAtfS+/6puBrokHQ4gaZikg+srdeDwEX+bRcRiSWcBV0t6GLiRdIR/qaR1pDeGSR0ssXgRcZekayUtAu4HVldG/wR4i6QlpOD4TSdqtCQififpIuCGPOj8iFhQ2X9Xki7uNpv3IUmvAP5X0m6kPDyHdJ1gu+avbDAzK4ybeszMCuPgNzMrjIPfzKwwDn4zs8I4+M3MCuPbOc22Uv4MxtzcuwfpA0Frc/+hEfFQRwoza5Fv5zR7DPIX6W2IiM92uhazVrmpx8ysMA5+M7PCOPjNzArj4DczK4yD38ysMA5+M7PC+HZOM7PC+IjfzKwwDn4zs8I4+M3MCuPgNzMrjIPfzKwwDn4zs8I4+M3MCvP/jGFxE+jmRFYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "T_samples = [big_city_model()['T'] for _ in range(100)]\n",
    "T_unique, T_counts = np.unique(T_samples, return_counts=True)\n",
    "plt.bar(T_unique, T_counts, align='center', alpha=0.5)\n",
    "plt.xticks(T_unique, T_alias)\n",
    "plt.ylabel('count')\n",
    "plt.xlabel('T')\n",
    "plt.title('P(T | do(R = big)) - Importance Sampling')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we calculate the causal effect of city size on car use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "small_city_model = pyro.do(model, data={'R': torch.tensor(0)})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.1033"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t_samples_big = [\n",
    "    1 if big_city_model()['T'] == 0.0 else 0\n",
    "    for _ in range(50000)\n",
    "]\n",
    "t_samples_small = [\n",
    "    1 if small_city_model()['T'] == 0.0 else 0\n",
    "    for _ in range(50000)\n",
    "]\n",
    "\n",
    "causal_effect = mean(t_samples_big) - mean(t_samples_small)\n",
    "causal_effect"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "venv",
   "language": "python",
   "name": "venv"
  },
  "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}