notebooks/causal-models-exercises.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Examples and exercises for causal models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from pathlib import Path\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "\n",
    "import pytest\n",
    "\n",
    "# Only needed to generate graphs, may be safely ommitted \n",
    "# once you comment out relevant cells below\n",
    "from graphviz import Digraph\n",
    "from brent import DAG, Query"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Expanded proof of quantities unchanged between $G$ and $G' = G_{\\underline{\\textrm{days}}}$\n",
    "\n",
    "from slides ``Correlation and Causality''\n",
    "\n",
    "Graph of ``given-days'' graph $G$:\n",
    "\n",
    "![Graph of hit-rate for conditional calculations](../slides/correlation-causality/graphics/given_days.png)\n",
    "\n",
    "Graph of ``do-days'' graph $G'$:\n",
    "\n",
    "![Graph of hit-rate for do calculations](../slides/correlation-causality/graphics/do_days.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Proof of $P_{G'}(\\textrm{producttype} = p, \\textrm{rating} = r)=  P_G( \\textrm{producttype} = p, \\textrm{rating} = r)$ \n",
    "\n",
    "We swap notation for this proof a bit, writing $Pr$ for probability, and $P$ for $\\textrm{producttype}$, $R$ for $\\textrm{rating}$, etc.\n",
    "\n",
    "By the definition of graphical models and the decomposition of conditional probabilities\n",
    "\n",
    "$$\n",
    "Pr_G(P, R, D, H) = Pr(P) \\cdot Pr_G(R|P) \\cdot {\\color{red}{Pr_G(D|P)}} \\cdot Pr_G(H | D, R)\n",
    "$$\n",
    "\n",
    "Note how we have dropped the subscript from $Pr(P)$, as the marginal probabilities of each vertex variable are independent of the underlying graph. For $G'$ we have\n",
    "\n",
    "$$\n",
    "Pr_{G'}(P, R, D, H) = Pr(P) \\cdot Pr_G(R|P) \\cdot {\\color{red}{Pr(D)}} \\cdot Pr_G(H | D, R)\n",
    "$$\n",
    "\n",
    "where we've highlighted in $\\color{red}{\\textrm{red}}$ the terms that differ between the two expressions.\n",
    "\n",
    "Note further that, by definition of the graph surgery we do in going from $G$ to $G'$, all conditional probabilities are preserved that aren't affected by the ``do''-surgery.\n",
    "\n",
    "Now consider $P_{G'}(\\textrm{producttype} = p, \\textrm{rating} = r)$. Using the conditioning trick\n",
    "\n",
    "\\begin{align*}\n",
    "Pr_{G'}(P, R) & = \\sum_{d, h} Pr_{G'}(P=p, R=r, D=d, H=h) \\\\\n",
    "& = Pr(P) \\cdot Pr_{G'}(R|P) \\cdot \\sum_{d, h} {\\color{red}{Pr(D)}} \\cdot Pr_G(H | D, R) \\textrm{, by the graph property} \\\\\n",
    "& = Pr(P) \\cdot Pr_{G'}(R|P) \\cdot \\sum_{d} {\\color{red}{Pr(D)}} \\sum_h \\cdot Pr_G(H | D, R) \\textrm{, by lack of dependence on $h$ summands}\\\\\n",
    "& = Pr(P) \\cdot Pr_{G'}(R|P) \\cdot \\sum_{d} {\\color{red}{Pr(D)}} \\cdot 1 \\textrm{, by the law of total probability} \\\\\n",
    "& = Pr(P) \\cdot Pr_{G'}(R|P) \\cdot 1 \\textrm{, by the law of total probability for $D$}\n",
    "\\end{align*}\n",
    "\n",
    "Similarly\n",
    "\n",
    "\\begin{align*}\n",
    "Pr_{G}(P, R) & = \\sum_{d, h} Pr_{G'}(P=p, R=r, D=d, H=h) \\\\\n",
    "& = Pr(P) \\cdot Pr_G(R|P) \\cdot \\sum_{d, h} {\\color{red}{Pr(D | P)}} \\cdot Pr_G(H | D, R) \\textrm{, by the graph property} \\\\\n",
    "& = Pr(P) \\cdot Pr_G(R|P) \\cdot \\sum_{d} {\\color{red}{Pr(D | P)}} \\sum_h Pr_G(H | D, R) \\textrm{, by lack of dependence on $h$ summands} \\\\\n",
    "& = Pr(P) \\cdot Pr_G(R|P) \\cdot \\sum_{d} {\\color{red}{Pr(D | P)}} \\cdot 1 \\textrm{, by the law of total probability} \\\\\n",
    "& = Pr(P) \\cdot Pr_G(R|P) \\cdot 1 \\textrm{, by the law of total probability for each of $D|P$}\n",
    "\\end{align*}\n",
    "\n",
    "Since conditional probabilities unaffected by the ``do'' graph surgery are unaffected, we have $Pr_{G'}(R|P) = Pr_G(R|P)$, and the result follows."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Proof of $P_{G'}(\\textrm{hit}=1 | \\textrm{producttype} = p, \\textrm{rating} = r)=  P_G(\\textrm{hit}=1 | \\textrm{producttype} = p, \\textrm{rating} = r)$ \n",
    "\n",
    "Using the same abbreviated notation as above,\n",
    "\n",
    "\\begin{align*}\n",
    "Pr_{G'}(H=1 | P, R) &= \\sum_{d} Pr_{G'}(H=1 | P, R, D=d) \\textrm{, by the conditioning trick}\\\\\n",
    "&= \\sum_{d} \\frac{\n",
    "    Pr(P) \\cdot Pr_G(R|P) \\cdot Pr(D) \\cdot Pr_{G'}(H = 1 | D, R)\n",
    "    }{\n",
    "        \\sum_h Pr(P) \\cdot Pr_{G'}(R|P) \\cdot Pr(D) \\cdot Pr_{G'}(H = h | D, R)\n",
    "    } \\textrm{, by definition of conditional probability and graphical models}\n",
    "\\end{align*}\n",
    "\n",
    "\n",
    "by the defining propert of graphical model $G'$ and using $Pr_{G'}(D) = Pr_G(D)$ (also for $(R | P)$)\n",
    "\n",
    "\\begin{align*}\n",
    "Pr_{G'}(H=1 | P, R) &=\\sum_{d} \\frac{\n",
    "    Pr(P) \\cdot Pr_{G'}(R|P) \\cdot Pr(D) \\cdot Pr_{G'}(H = 1 | D, R)\n",
    "    }{\n",
    "        Pr(P) \\cdot Pr_{G'}(R|P) \\cdot Pr(D) \\sum_h  Pr_{G'}(H = h | D, R)\n",
    "    } \\\\\n",
    "    &=\\sum_{d} \\frac{\n",
    "        Pr(P) \\cdot Pr_{G'}(R|P) \\cdot Pr(D) \n",
    "        \\cdot Pr_{G'}(H = 1 | D, R)\n",
    "    }{\n",
    "        Pr(P) \\cdot Pr_{G'}(R|P) \\cdot Pr(D) \\sum_h  Pr_{G'}(H = h | D, R)\n",
    "    } \\\\\n",
    "    &=\\sum_{d} \\frac{\n",
    "        Pr_{G'}(H = 1 | D, R)\n",
    "    }{\n",
    "        \\sum_h  Pr_{G'}(H = h | D, R)\n",
    "    } \\textrm{, as long as none of the cancelled terms is 0} \\\\\n",
    "    &=\\sum_{d} \\frac{\n",
    "        Pr_{G'}(H = 1 | D, R)\n",
    "    }{\n",
    "        \\sum_h  Pr_{G'}(H = h | D, R)\n",
    "    } \\textrm{, as long as none of the cancelled terms is 0} \\\\\n",
    "    &=\\sum_{d} \\frac{\n",
    "        Pr_{G}(H = 1 | D, R)\n",
    "    }{\n",
    "        \\sum_h  Pr_{G}(H = h | D, R)\n",
    "    } \\textrm{, since the events $(H | D=d, R=r)$ are unaffected by the ``do'' surgery} \\\\\n",
    "\\end{align*}\n",
    "\n",
    "For $Pr_{G}(H=1 | P, R)$ note that the only different terms are the ones in red-above, namely $Pr(D)$ from above is replaced by $Pr_G(D|P)$. As $Pr_G(D|P)$ is also independent of the value of $h$, we may similarly pull it outside of the summation over $h$. Provided these probabilities are also non-zero, we obtain\n",
    "\n",
    "\\begin{equation*}\n",
    "Pr_{G}(H=1 | P, R) = =\\sum_{d} \\frac{\n",
    "        Pr_{G}(H = 1 | D, R)\n",
    "    }{\n",
    "        \\sum_h  Pr_{G}(H = h | D, R)\n",
    "    }\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Causal model example: hit rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "datadir = Path(os.getcwd()) / 'data'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(5389, 4)\n"
     ]
    },
    {
     "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>producttype</th>\n",
       "      <th>days</th>\n",
       "      <th>rating</th>\n",
       "      <th>hit</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>liability</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>liability</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>property</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>property</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>property</td>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  producttype  days  rating  hit\n",
       "0   liability     0       0    1\n",
       "1   liability     2       1    0\n",
       "2    property     3       1    1\n",
       "3    property     3       1    0\n",
       "4    property     3       0    0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "path_or_url = datadir / 'hits.csv'\n",
    "if not path_or_url.exists():\n",
    "    path_or_url = 'https://raw.githubusercontent.com/munichpavel/risk-ai-workshop/main/notebooks/data/hits.csv'\n",
    "df = pd.read_csv(path_or_url)\n",
    "print(df.shape)\n",
    "# Remove _ from column names and it messes with latex rendering\n",
    "column_names = df.columns\n",
    "column_names = [c.replace('_', '') for c in column_names]\n",
    "df.columns = column_names\n",
    "n_records = df.shape[0]\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 8.1.0 (20230707.0739)\n",
       " -->\n",
       "<!-- Pages: 1 -->\n",
       "<svg width=\"296pt\" height=\"146pt\"\n",
       " viewBox=\"0.00 0.00 295.98 146.22\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 142.22)\">\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-142.22 291.98,-142.22 291.98,4 -4,4\"/>\n",
       "<!-- producttype -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>producttype</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"54.93\" cy=\"-71.41\" rx=\"54.93\" ry=\"54.93\"/>\n",
       "<text text-anchor=\"middle\" x=\"54.93\" y=\"-66.36\" font-family=\"Times,serif\" font-size=\"14.00\">producttype</text>\n",
       "</g>\n",
       "<!-- days -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>days</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"178.28\" cy=\"-110.41\" rx=\"27.81\" ry=\"27.81\"/>\n",
       "<text text-anchor=\"middle\" x=\"178.28\" y=\"-105.36\" font-family=\"Times,serif\" font-size=\"14.00\">days</text>\n",
       "</g>\n",
       "<!-- producttype&#45;&gt;days -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>producttype&#45;&gt;days</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M107.75,-88.07C118.94,-91.66 130.53,-95.39 140.98,-98.75\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"139.79,-102.36 150.38,-102.09 141.93,-95.7 139.79,-102.36\"/>\n",
       "</g>\n",
       "<!-- rating -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>rating</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"178.28\" cy=\"-32.41\" rx=\"32.41\" ry=\"32.41\"/>\n",
       "<text text-anchor=\"middle\" x=\"178.28\" y=\"-27.36\" font-family=\"Times,serif\" font-size=\"14.00\">rating</text>\n",
       "</g>\n",
       "<!-- producttype&#45;&gt;rating -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>producttype&#45;&gt;rating</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M107.75,-54.76C117.38,-51.67 127.31,-48.47 136.55,-45.5\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"137.57,-48.53 146.02,-42.14 135.43,-41.87 137.57,-48.53\"/>\n",
       "</g>\n",
       "<!-- hit -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>hit</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"267.34\" cy=\"-71.41\" rx=\"20.64\" ry=\"20.64\"/>\n",
       "<text text-anchor=\"middle\" x=\"267.34\" y=\"-66.36\" font-family=\"Times,serif\" font-size=\"14.00\">hit</text>\n",
       "</g>\n",
       "<!-- days&#45;&gt;hit -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>days&#45;&gt;hit</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M204,-99.34C214.59,-94.6 227.02,-89.03 237.98,-84.12\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"239.24,-86.94 246.93,-79.66 236.38,-80.55 239.24,-86.94\"/>\n",
       "</g>\n",
       "<!-- rating&#45;&gt;hit -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>rating&#45;&gt;hit</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M208.27,-45.4C217.93,-49.73 228.68,-54.54 238.27,-58.84\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"236.56,-62.36 247.12,-63.25 239.42,-55.97 236.56,-62.36\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.graphs.Digraph at 0x11cb42740>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dag = DAG(df).add_edge('producttype', 'rating').add_edge('producttype', 'days')\\\n",
    "    .add_edge('rating', 'hit').add_edge('days', 'hit')\n",
    "dag.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 8.1.0 (20230707.0739)\n",
       " -->\n",
       "<!-- Pages: 1 -->\n",
       "<svg width=\"296pt\" height=\"154pt\"\n",
       " viewBox=\"0.00 0.00 295.98 154.22\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 150.22)\">\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-150.22 291.98,-150.22 291.98,4 -4,4\"/>\n",
       "<!-- producttype -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>producttype</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"54.93\" cy=\"-73.41\" rx=\"54.93\" ry=\"54.93\"/>\n",
       "<text text-anchor=\"middle\" x=\"54.93\" y=\"-68.36\" font-family=\"Times,serif\" font-size=\"14.00\">producttype</text>\n",
       "</g>\n",
       "<!-- days -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>days</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"178.28\" cy=\"-114.41\" rx=\"27.81\" ry=\"27.81\"/>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"178.28\" cy=\"-114.41\" rx=\"31.81\" ry=\"31.81\"/>\n",
       "<text text-anchor=\"middle\" x=\"178.28\" y=\"-109.36\" font-family=\"Times,serif\" font-size=\"14.00\">days</text>\n",
       "</g>\n",
       "<!-- producttype&#45;&gt;days -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>producttype&#45;&gt;days</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M107.4,-90.8C117.55,-94.23 128.05,-97.78 137.74,-101.05\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"136.27,-104.59 146.87,-104.48 138.51,-97.96 136.27,-104.59\"/>\n",
       "</g>\n",
       "<!-- rating -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>rating</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"178.28\" cy=\"-32.41\" rx=\"32.41\" ry=\"32.41\"/>\n",
       "<text text-anchor=\"middle\" x=\"178.28\" y=\"-27.36\" font-family=\"Times,serif\" font-size=\"14.00\">rating</text>\n",
       "</g>\n",
       "<!-- producttype&#45;&gt;rating -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>producttype&#45;&gt;rating</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M107.4,-56.03C117.25,-52.7 127.44,-49.25 136.89,-46.06\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"137.83,-49.1 146.18,-42.58 135.59,-42.47 137.83,-49.1\"/>\n",
       "</g>\n",
       "<!-- hit -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>hit</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"267.34\" cy=\"-73.41\" rx=\"20.64\" ry=\"20.64\"/>\n",
       "<text text-anchor=\"middle\" x=\"267.34\" y=\"-68.36\" font-family=\"Times,serif\" font-size=\"14.00\">hit</text>\n",
       "</g>\n",
       "<!-- days&#45;&gt;hit -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>days&#45;&gt;hit</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M207.31,-101.22C217.27,-96.53 228.47,-91.25 238.43,-86.56\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"239.68,-89.37 247.24,-81.94 236.7,-83.03 239.68,-89.37\"/>\n",
       "</g>\n",
       "<!-- rating&#45;&gt;hit -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>rating&#45;&gt;hit</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M207.79,-45.84C217.65,-50.48 228.69,-55.68 238.5,-60.3\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"236.64,-63.77 247.18,-64.86 239.62,-57.43 236.64,-63.77\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.graphs.Digraph at 0x11cb42c20>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q_given = Query(dag).given(days=1)\n",
    "q_given.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 8.1.0 (20230707.0739)\n",
       " -->\n",
       "<!-- Pages: 1 -->\n",
       "<svg width=\"296pt\" height=\"177pt\"\n",
       " viewBox=\"0.00 0.00 295.98 176.74\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 172.74)\">\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-172.74 291.98,-172.74 291.98,4 -4,4\"/>\n",
       "<!-- producttype -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>producttype</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"54.93\" cy=\"-54.93\" rx=\"54.93\" ry=\"54.93\"/>\n",
       "<text text-anchor=\"middle\" x=\"54.93\" y=\"-49.88\" font-family=\"Times,serif\" font-size=\"14.00\">producttype</text>\n",
       "</g>\n",
       "<!-- days -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>days</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"178.28\" cy=\"-136.93\" rx=\"27.81\" ry=\"27.81\"/>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"178.28\" cy=\"-136.93\" rx=\"31.81\" ry=\"31.81\"/>\n",
       "<text text-anchor=\"middle\" x=\"178.28\" y=\"-131.88\" font-family=\"Times,serif\" font-size=\"14.00\">days</text>\n",
       "</g>\n",
       "<!-- producttype&#45;&gt;days -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>producttype&#45;&gt;days</title>\n",
       "<path fill=\"none\" stroke=\"lightgray\" stroke-dasharray=\"5,2\" d=\"M101.12,-85.47C114.78,-94.7 129.57,-104.69 142.45,-113.4\"/>\n",
       "<polygon fill=\"lightgray\" stroke=\"lightgray\" points=\"140.2,-116.78 150.45,-119.47 144.12,-110.98 140.2,-116.78\"/>\n",
       "</g>\n",
       "<!-- rating -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>rating</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"178.28\" cy=\"-54.93\" rx=\"32.41\" ry=\"32.41\"/>\n",
       "<text text-anchor=\"middle\" x=\"178.28\" y=\"-49.88\" font-family=\"Times,serif\" font-size=\"14.00\">rating</text>\n",
       "</g>\n",
       "<!-- producttype&#45;&gt;rating -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>producttype&#45;&gt;rating</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M110.22,-54.93C118.45,-54.93 126.85,-54.93 134.79,-54.93\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"134.64,-58.43 144.64,-54.93 134.64,-51.43 134.64,-58.43\"/>\n",
       "</g>\n",
       "<!-- hit -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>hit</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"267.34\" cy=\"-95.93\" rx=\"20.64\" ry=\"20.64\"/>\n",
       "<text text-anchor=\"middle\" x=\"267.34\" y=\"-90.88\" font-family=\"Times,serif\" font-size=\"14.00\">hit</text>\n",
       "</g>\n",
       "<!-- days&#45;&gt;hit -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>days&#45;&gt;hit</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M207.31,-123.74C217.27,-119.04 228.47,-113.77 238.43,-109.08\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"239.68,-111.89 247.24,-104.46 236.7,-105.55 239.68,-111.89\"/>\n",
       "</g>\n",
       "<!--   -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title> </title>\n",
       "<text text-anchor=\"middle\" x=\"54.93\" y=\"-140.88\" font-family=\"Times,serif\" font-size=\"14.00\"> </text>\n",
       "</g>\n",
       "<!--  &#45;&gt;days -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title> &#45;&gt;days</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M81.67,-144.03C97.34,-142.86 117.73,-141.35 135.65,-140.02\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"135.71,-143.45 145.43,-139.22 135.19,-136.47 135.71,-143.45\"/>\n",
       "</g>\n",
       "<!-- rating&#45;&gt;hit -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>rating&#45;&gt;hit</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M207.79,-68.36C217.65,-73 228.69,-78.2 238.5,-82.82\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"236.64,-86.29 247.18,-87.38 239.62,-79.95 236.64,-86.29\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.graphs.Digraph at 0x11cbc8370>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q_do = Query(dag).do(days=1)\n",
    "q_do.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Probability equations\n",
    "\n",
    "for the two graphs, the original, $G$, and the ``do'' one, $G'$. If the calculation applies to both graphs, we write $P_{G, G'}$.\n",
    "\n",
    "$$\\begin{align}\n",
    "P_{G, G'}(H=1 | D=d) &= \\sum_{p, r} P_{G, G'}(H=1 | D=d, P=p, R=r) P_{G, G'}(P=p, R=r | D=d)\n",
    "\\end{align}$$\n",
    "\n",
    "by the law of total probability, aka conditioning trick.\n",
    "\n",
    "Only for $G'$ do we have $P_{G'}(P = p, R = r | D = d) = P_{G'}(P = p, R = r)$, as the pair $P, R$ are independent of $D$ in the subgraph of $G'$ (no edge).\n",
    "\n",
    "Hence we have \n",
    "\n",
    "$$\\begin{align}\n",
    "P_{G}(H=1 | D=d) &= \\sum_{p, r} P_{G}(H=1 | D=d, P=p, R=r) \\color{red}{P_{G}(P=p, R=r | D=d)}\n",
    "\\end{align}$$\n",
    "\n",
    "for $G$, but for $G'$\n",
    "\n",
    "$$\\begin{align}\n",
    "P_{G'}(H=1 | D=d) &= \\sum_{p, r} P_{G'}(H=1 | D=d, P=p, R=r) \\color{red}{P_{G'}(P=p, R=r)}\n",
    "\\end{align}$$\n",
    "\n",
    "Using equalities from lecture, have right-hand side now in terms of original graph $G$:\n",
    "\n",
    "$$\\begin{equation}\n",
    "P_{G'}(H=1 | D=d) = \\sum_{p, r} P_{G}(H=1 | D=d, P=p, R=r) P_{G}(P=p, R=r)\n",
    "\\end{equation}$$\n",
    "\n",
    "Before calculating the whole right-hand-sides, let's convince ourselves that there is a difference between the two differing terms from $G$ to $G'$.\n",
    "\n",
    "In the $G'$ expression we can use the law of total probability again to see"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "producttype  rating\n",
       "financial    0         0.040453\n",
       "             1         0.169605\n",
       "liability    0         0.339210\n",
       "             1         0.149935\n",
       "property     0         0.120059\n",
       "             1         0.180739\n",
       "dtype: float64"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# P_G(P=p, R=r)\n",
    "P_Gpr = df.groupby(['producttype', 'rating']).size() / n_records\n",
    "P_Gpr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# P_G(P=p, R=r | D=d)\n",
    "P_Gpr_given_ds = []\n",
    "for d in range(4):\n",
    "    df_mask_d = df['days'] == d\n",
    "    n_records_d = sum(df_mask_d)\n",
    "    P_Gpr_given_d = df.loc[df_mask_d].groupby(['producttype', 'rating', 'days']).size() / n_records_d\n",
    "    P_Gpr_given_ds.append(P_Gpr_given_d)\n",
    "\n",
    "# Non-do calculation, 2nd term in product\n",
    "p_pr_given_ds = pd.concat(P_Gpr_given_ds, axis=0)\n",
    "p_pr_given_ds.name = 'hit'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the red terms in the above equations were identical, then the values for each group of rows with fixed $D=d$ would be identical copies of the above variable `P_Gpr`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "producttype  rating  days\n",
       "financial    0       0       0.091294\n",
       "             1       0       0.374588\n",
       "liability    0       0       0.343059\n",
       "             1       0       0.143529\n",
       "property     0       0       0.016941\n",
       "             1       0       0.030588\n",
       "financial    0       1       0.017527\n",
       "             1       1       0.079844\n",
       "liability    0       1       0.510224\n",
       "             1       1       0.252191\n",
       "property     0       1       0.062317\n",
       "             1       1       0.077897\n",
       "financial    0       2       0.002413\n",
       "             1       2       0.015682\n",
       "liability    0       2       0.533172\n",
       "             1       2       0.225573\n",
       "property     0       2       0.095296\n",
       "             1       2       0.127865\n",
       "financial    0       3       0.002841\n",
       "             1       3       0.016335\n",
       "liability    0       3       0.094460\n",
       "             1       3       0.040483\n",
       "property     0       3       0.332386\n",
       "             1       3       0.513494\n",
       "Name: hit, dtype: float64"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_pr_given_ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The probabilities do indeed differ across different $D=d$ values, so the ``do'' probabilities from $G'$ are different from the non-do probabilities from $G$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sanity check on P_G(P=p, R=r | D=d) and P_G(P=p, R=r) calculations\n",
    "# assertion will fail if law of total probability does not hold\n",
    "idx_pr = ['financial', 0]\n",
    "\n",
    "expected_P_Gpr = P_Gpr.loc[tuple(idx_pr)]\n",
    "\n",
    "P_Gd = df.groupby(['days']).size() / n_records\n",
    "cond_prob_sum = 0\n",
    "for d, P_Gpr_given_d in enumerate(P_Gpr_given_ds):\n",
    "    idx_prd = tuple(idx_pr + [d])\n",
    "    \n",
    "    cond_prob_sum += P_Gpr_given_d.loc[idx_prd] * P_Gd.loc[d]\n",
    "assert cond_prob_sum == pytest.approx(expected_P_Gpr), f'Conditional proby sum {cond_prob_sum} does not match expected {expected_P_Gpr}'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculation of $P_G(H=1 | D=d, P=p, R=r)$\n",
    "\n",
    "Can use, e.g.\n",
    "\n",
    "$$\\begin{equation}\n",
    "P_G(H=1 | P=p, D=d, R=r) = \\frac{P_G(P=p, D=d, R=r, H=1)}{P_G(P=p, D=d, R=r)}\n",
    "\\end{equation}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "days\n",
       "0    0.565343\n",
       "1    0.397330\n",
       "2    0.240322\n",
       "3    0.215639\n",
       "Name: prob, dtype: float64"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# P(H=1 | do(D)=d) = \\sum P(H=1 | D=d, P=p, R=r) * P(P=p, R=r)\n",
    "# First term in sum\n",
    "hit_given_prd = df.groupby(['producttype', 'rating', 'days'])\n",
    "hit_given_prd = hit_given_prd['hit'].mean()\n",
    "\n",
    "# 2nd sum in P(H=1 | do(D)=1)\n",
    "pr = df.groupby(['producttype', 'rating'])\n",
    "p_pr = pr['hit'].count()\n",
    "p_pr = p_pr / df.shape[0]\n",
    "\n",
    "# Combine\n",
    "hit_given_prd = hit_given_prd.reset_index()\n",
    "p_pr = p_pr.reset_index()\n",
    "\n",
    "res_do = pd.merge(hit_given_prd, p_pr, on=['producttype', 'rating'])\n",
    "res_do['prob'] = res_do['hit_x'] * res_do['hit_y']\n",
    "\n",
    "res_do_grouped = res_do.groupby('days')\n",
    "res_do_sum = res_do_grouped['prob'].sum()\n",
    "res_do_sum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quick calculation of $P(H=1 | D=d)$\n",
    "\n",
    "but not helpful for debugging, convincing myself that the results are correct"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# P(H=1 | D=d)\n",
    "hit_given_d = df.groupby(['days'])\n",
    "hit_given_d = hit_given_d['hit'].mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In more detail for debugging / sanity-checking"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Non-causal, first term in product different\n",
    "# Non-do calculation, 2nd term in product\n",
    "\n",
    "p_pr_given_ds = p_pr_given_ds.reset_index()\n",
    "res_given = pd.merge(hit_given_prd, p_pr_given_ds, on=['producttype', 'rating', 'days'])\n",
    "res_given['prob'] = res_given['hit_x'] * res_given['hit_y']\n",
    "res_given_grouped = res_given.groupby('days')\n",
    "res_given_sum = res_given_grouped['prob'].sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate average treatment effect (ATE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "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>from-d</th>\n",
       "      <th>to-d</th>\n",
       "      <th>ate-given</th>\n",
       "      <th>ate-do</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0.170153</td>\n",
       "      <td>0.297187</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>0.252329</td>\n",
       "      <td>0.395158</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>0.473538</td>\n",
       "      <td>0.102707</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   from-d  to-d  ate-given    ate-do\n",
       "0       0     1   0.170153  0.297187\n",
       "1       1     2   0.252329  0.395158\n",
       "2       2     3   0.473538  0.102707"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ates = []\n",
    "ds = range(3)\n",
    "for d in ds:\n",
    "    rel_ate_do = (res_do_sum[d] - res_do_sum[d + 1]) / res_do_sum[d]\n",
    "    rel_ate_given = (res_given_sum[d] - res_given_sum[d + 1]) / res_given_sum[d]\n",
    "    ates.append({'from-d': d, 'to-d': d+1, 'ate-given': rel_ate_given, 'ate-do': rel_ate_do})\n",
    "ates_df = pd.DataFrame(ates)\n",
    "ates_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Causal models exercise: correlation\n",
    "\n",
    "Reproduce and try to break the spurious correlation between deaths by poisonous spider bites and the lenghts of winning words in the Scripps national spelling bee.\n",
    "\n",
    "The fatality data from the CDC can be found here: `notebooks/data/cdc-underlying-cause-of-death-1998-2018.txt`, and the spelling bee data can be found below.\n",
    "\n",
    "Difficulty: **"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# From https://en.wiktionary.org/wiki/Appendix:Scripps_winning_words\n",
    "scripps_winners_raw = '''\n",
    "    1925: gladiolus\n",
    "    1926: abrogate\n",
    "    1927: luxuriance\n",
    "    1928: albumen\n",
    "    1929: asceticism\n",
    "    1930: fracas\n",
    "    1931: foulard\n",
    "    1932: knack\n",
    "    1933: torsion\n",
    "    1934: deteriorating\n",
    "    1935: intelligible\n",
    "    1936: interning\n",
    "    1937: promiscuous\n",
    "    1938: sanitarium\n",
    "    1939: canonical\n",
    "    1940: therapy\n",
    "    1941: initials\n",
    "    1942: sacrilegious\n",
    "\n",
    "The Bee was suspended during the WWII years of 1943–1945.\n",
    "\n",
    "    1946: semaphore\n",
    "    1947: chlorophyll\n",
    "    1948: psychiatry\n",
    "    1949: dulcimer\n",
    "    1950: meerschaum [1] / meticulosity\n",
    "    1951: insouciant\n",
    "    1952: vignette\n",
    "    1953: soubrette\n",
    "    1954: transept\n",
    "    1955: crustaceology\n",
    "    1956: condominium\n",
    "    1957: n/a [2]\n",
    "    1958: syllepsis\n",
    "    1959: catamaran\n",
    "    1960: eudaemonic\n",
    "    1961: smaragdine\n",
    "    1962: n/a [3]\n",
    "    1963: equipage\n",
    "    1964: sycophant\n",
    "    1965: eczema\n",
    "    1966: ratoon\n",
    "    1967: chihuahua\n",
    "    1968: abalone\n",
    "    1969: interlocutory\n",
    "    1970: croissant\n",
    "    1971: shalloon\n",
    "    1972: macerate\n",
    "\n",
    "    1973: vouchsafe\n",
    "    1974: hydrophyte\n",
    "    1975: incisor\n",
    "    1976: narcolepsy\n",
    "    1977: cambist\n",
    "    1978: deification\n",
    "    1979: maculature\n",
    "    1980: elucubrate\n",
    "    1981: sarcophagus\n",
    "    1982: psoriasis\n",
    "    1983: Purim\n",
    "    1984: luge\n",
    "    1985: milieu\n",
    "    1986: odontalgia\n",
    "    1987: staphylococci\n",
    "    1988: elegiacal\n",
    "    1989: spoliator\n",
    "    1990: fibranne\n",
    "    1991: antipyretic\n",
    "    1992: lyceum\n",
    "    1993: kamikaze\n",
    "    1994: antediluvian\n",
    "    1995: xanthosis\n",
    "    1996: vivisepulture\n",
    "    1997: euonym\n",
    "    1998: chiaroscurist\n",
    "    1999: logorrhea\n",
    "    2000: demarche\n",
    "    2001: succedaneum\n",
    "    2002: prospicience\n",
    "    2003: pococurante\n",
    "    2004: autochthonous\n",
    "    2005: appoggiatura\n",
    "    2006: Ursprache\n",
    "    2007: serrefine\n",
    "    2008: guerdon\n",
    "    2009: Laodicean\n",
    "    2010: stromuhr\n",
    "    2011: cymotrichous\n",
    "    2012: guetapens\n",
    "    2013: knaidel\n",
    "    2014: stichomythia / feuilleton\n",
    "    2015: scherenschnitte / nunatak\n",
    "    2016: Feldenkrais / gesellschaft\n",
    "    2017: marocain\n",
    "    2018: koinonia\n",
    "    2019: auslaut / erysipelas / bougainvillea [4] / aiguillette / pendeloque / palama / cernuous / odylic\n",
    "'''"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Causal models exercise: do-calculus\n",
    "\n",
    "As before, take K to be your Karma, H to be the hours you spend in the gym lifting weight, and then W be the weight you can bench press. \n",
    "\n",
    "You are the parent of a very young child, so Karma will punish you for devoting too much time to your triceps and neglecting your partner and baby. Let $G$ be this causal graph, as shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 8.1.0 (20230707.0739)\n",
       " -->\n",
       "<!-- Pages: 1 -->\n",
       "<svg width=\"134pt\" height=\"106pt\"\n",
       " viewBox=\"0.00 0.00 134.21 106.41\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 102.41)\">\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-102.41 130.21,-102.41 130.21,4 -4,4\"/>\n",
       "<!-- K -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>K</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-79.93\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-74.88\" font-family=\"Times,serif\" font-size=\"14.00\">K</text>\n",
       "</g>\n",
       "<!-- H -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>H</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"99.21\" cy=\"-80.41\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"99.21\" y=\"-75.36\" font-family=\"Times,serif\" font-size=\"14.00\">H</text>\n",
       "</g>\n",
       "<!-- K&#45;&gt;H -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>K&#45;&gt;H</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M54.3,-80.11C56.54,-80.13 58.84,-80.14 61.15,-80.16\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"60.85,-83.66 70.87,-80.23 60.89,-76.66 60.85,-83.66\"/>\n",
       "</g>\n",
       "<!-- W -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>W</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"63.59\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"63.59\" y=\"-12.95\" font-family=\"Times,serif\" font-size=\"14.00\">W</text>\n",
       "</g>\n",
       "<!-- K&#45;&gt;W -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>K&#45;&gt;W</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M37,-63.01C40.47,-57.14 44.44,-50.41 48.24,-43.98\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"51.53,-46.29 53.61,-35.9 45.51,-42.73 51.53,-46.29\"/>\n",
       "</g>\n",
       "<!-- H&#45;&gt;W -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>H&#45;&gt;W</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M89.47,-63.35C86.1,-57.44 82.23,-50.66 78.53,-44.18\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"81.31,-42.98 73.31,-36.03 75.23,-46.45 81.31,-42.98\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.graphs.Digraph at 0x11ccb6ef0>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dot = Digraph(engine='neato')\n",
    "dot.attr('node')\n",
    "dot.node('K')\n",
    "dot.node('H')\n",
    "dot.node('W')\n",
    "\n",
    "dot.edge('K', 'H')\n",
    "dot.edge('K', 'W')\n",
    "dot.edge('H', 'W')\n",
    "\n",
    "dot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Draw the graphs $G_\\underline{W}$ and $G_\\overline{H}$. Difficulty: *\n",
    "2. Write out formulas for $P(W=1 | H=1)$ and $P(W=1|\\, \\mathrm{do}(H) = 1)$. Difficulty: **\n",
    "\n",
    "3. Calculate $P(W=1 | H=1)$ and $P(W=1|\\, \\mathrm{do}(H) = 1)$ for a Bayesian network fitted to the sample data from $(K, H, W)$ in `notebooks/data/karma_weights.csv`. Hint: the `Query` class of [https://koaning.github.io/brent/](https://koaning.github.io/brent/) can be used. Interpret the results in a qualitative way, i.e. how do you think Karma should work in this situation? Difficulty: **"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Causal models exercise: Causal calculus\n",
    "\n",
    "Prove in gory detail that the special case of Causal rule 3 holds. Difficulty: *"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}