solutions_pgmpy.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Causal Inference in Statistics: A Primer\n",
    "Bill Behrman 2023-05-31\n",
    "\n",
    "-   1 Preliminaries: Statistical and Causal Models\n",
    "    -   1.4 Graphs\n",
    "        -   Study question 1.4.1\n",
    "        -   Study question 1.5.1\n",
    "-   2 Graphical Models and Their Applications\n",
    "    -   2.3 Colliders\n",
    "        -   Study question 2.3.1\n",
    "    -   2.4 d-separation\n",
    "        -   Study question 2.4.1\n",
    "    -   2.5 Model Testing and Causal Search\n",
    "        -   Study question 2.5.1\n",
    "-   3 The Effects of Interventions\n",
    "    -   3.3 The Backdoor Criterion\n",
    "        -   Study question 3.3.1\n",
    "    -   3.4 The Front-Door Criterion\n",
    "        -   Study question 3.4.1\n",
    "    -   3.5 Conditional Interventions and Covariate-Specific Effects\n",
    "        -   Study question 3.5.1\n",
    "    -   3.8 Causal Inference in Linear Systems\n",
    "        -   Study question 3.8.1\n",
    "\n",
    "The following are solutions using the [pgmpy](https://pgmpy.org/) Python package to study questions from the book\n",
    "\n",
    "J Pearl, M Glymour, NP Jewell. Causal Inference in Statistics: A Primer. Wiley, 2016.\n",
    "\n",
    "The [book’s website](http://bayes.cs.ucla.edu/PRIMER/) contains PDFs for each chapter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-04-14T14:20:51.670773Z",
     "start_time": "2023-04-14T14:20:51.251306Z"
    }
   },
   "outputs": [],
   "source": [
    "import networkx as nx\n",
    "from pgmpy.base import DAG\n",
    "from pgmpy.independencies import IndependenceAssertion, Independencies\n",
    "from pgmpy.inference import CausalInference\n",
    "from pgmpy.models import BayesianNetwork"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code below uses pgmpy version 0.1.22."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1 Preliminaries: Statistical and Causal Models\n",
    "\n",
    "### 1.4 Graphs\n",
    "\n",
    "Figure 1.8. A directed graph used in study question 1.4.1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARgAAAEYCAYAAACHjumMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAAgUElEQVR4nO3de3hU5b0v8O/vXZOEEAUiqFHJAyhQio0o2yJSqEIRRYq31hu69ZGKVryUIiJHRUClirpF8RT0iLjVKu4eLxvd3th0c1QKCEW5KGwgICAIiDQFQgLJrPd3/piJGxCSua2ZWWt9P8/DY4BZ7/uqk2++6zqiqiAi8oLJ9QKIKLgYMETkGQYMEXmGAUNEnmHAEJFnGDBE5BkGDBF5hgFDRJ5hwBCRZxgwROSZSK4XQPlNRBwAvwAwGMDJAAoB7AOwEsAsAAuU95vQEQjfG3Q4IjLEGPOAtfYUEdEWLVq4bdq0MZFIROrq6nTHjh127969EVVVY8xSa+1oVZ2T63VTfmHA0EFEZIAx5nVVPeq0007ToUOHmpNPPvmwr7XWYvny5XjxxRfdyspKxxizzVo7QFVXZHnZlKcYMAQAkJh/V9WLzj77bB01apQUFhYmvP3u3bsxYcIEd82aNQbAFFUd4dliyTcYMAQRKTDGrHIcp8O4ceNMt27dUh5rzpw5ePrppxXAX1zXPS9zqyQ/YsCEnIiIMWZVUVFRp2effdaUlpamPWZlZSVGjRql1tr3rbWDMrBM8imepqZpxphO06ZNy0i4AEDHjh3x6KOPCoALReS3GRmUfIkBE2IicgaAm0aMGGFat26d0bE7d+6MSy+9FCLytIhkdnDyDQZMiBlj3urSpYuec845nox/ww03oLS0VETkz55MQHmPARNSItLVWttu5MiRnr4Hhg0b5qhqXxEp9nIeyk8MmJASkSeOP/5494QTTvB0nt69e6OoqMgCGO/pRJSXGDAhJSK9BwwY4GRjrp/+9KeOMeZX2ZiL8gsDJoRExLHWlvzsZz9r9HXWWlx11VXf/379+vUYPHgwtmzZAgCora3F9ddfj5qamkbHOeOMM6Cq5emvnPyGARNOvxARPemkkxp9kTEHvz0++OADHHPMMd8HykcffYQePXqgefPmjY7Ts2dPqGqhiJSkt2zyGwZMOLUvKCiwibywqKgI+/btQ21tLb788kucffbZqK2tBQB8+OGHGDhwYJNjtGjRouHLxhONAoePawinAhFJ6IXNmzdHbW0tPv30U/Tp0wfRaBQ1NTVYt24dHMfBkW6EPIKilFZLvsUGE05V0Wg0oYRpCJg5c+ZgwIABKC4uRm1tLT744ANccMEFCU1m7fdl6bsU10s+xYAJp/90Xdc0dXAWiAXMsmXL0Lp1axxzzDEoLi5GVVUVPvvsM/Tp0yehyZYvXw4AVlW3prds8hsGTAip6g4Rqf/000+bfG1xcTFmzZr1/bGW4uJizJ49Gz179kRRUWJ7PIsWLYLjODvTWjT5EgMmpETkq7lz5zZ5oLd58+ZQVTQ8wqG4uBhbtmxJePcIABYvXuy6rtt0mlHg8HENISUi14vIC6+//npSD5ZK1rZt2zBs2DAAOFVVV3o2EeUlNpiQUtUXAeyfOXOmp/M888wzaozZzHAJJwZMiKnqQ2+88Ybu3OnN4ZEvv/wSS5YsEWvtcE8moLzHXaSQcxynsqysrP2zzz6b0fuSotEorrnmGre2tnaOtTbxAzYUKGwwIWet7bt161Z96KGHErqyN8Ex8bvf/c6tra2tUdWLMzUu+Q8DhkpUdcOnn34q999/vz3goriU7Nu3D8OHD3e//vprR1VrVXV/htZJPsSACSkRcURkDIAvEPvExquXLl1aN2TIEHfVqlUpjTlv3jwMGTLEbt26dZeqDgZwnIjw4roQ4zGYEBKRLgD+L2LBEgHwB1WdICLFIvKOqv7iRz/6kR06dKjp2rVrk+MtWLAAL730krt582YHwIsAblBVFZGeABYA2Kaq3j7ZivISAyZkRGQUgAcANIv/0WoAFaoaPeA1/Y0xT1hrK5o3bx5t376907VrV+nUqRMKCwtRU1OD1atXY9WqVXbjxo1aV1dnRGSeqt6uqssOmY8hE2IMmJARkW8AlAEQALUAeqjqF0d47bEA7hSRPsaYrtbaoxBrPPXGmF2u6y4H8F8AnlDV2kbmZMiEFAMmRETEAdDQVOoQ3zVKYvuOANaqamLPejh4W4ZMCPEgb0gcEi7nAngCwMRsza+qCwGcDaCMB37Dgw0mBA4Jl7aquiXFcVJuMAeMwSYTImwwAZepcMkUNplwYcAEWL6FSwOGTHgwYAIqX8OlAUMmHBgwAZTv4dKAIRN8DJiA8Uu4NGDIBBsDJkD8Fi4NGDLBxYAJCL+GSwOGTDAxYALA7+HSgCETPAwYnwtKuDRgyAQLA8bHghYuDRgywcGA8amghksDhkwwMGB8KOjh0oAh438MGJ8JS7g0YMj4GwPGR8IWLg0YMv7FgPGJsIZLA4aMPzFgfCDs4dKAIeM/DJg8x3A5GEPGXxgweYzhcngMGf9gwOQphkvjGDL+wIDJQwyXxDBk8h8DJs8wXJLDkMlvDJg8wnBJDUMmfzFg8gTDJT0MmfzEgMkDDJfMYMjkHwZMjjFcMoshk18YMDnEcPEGQyZ/MGByhOHiLYZMfmDA5ADDJTsYMrnHgMkyhkt2MWRyiwGTRQyX3GDI5A4DJksYLrnFkMkNBkwWMFzyA0Mm+xgwHmO45BeGTHYxYDzEcMlPDJnsYcB4hOGS3xgy2cGA8QDDxR8YMt5jwGQYw8VfGDLeYsBkEMPFnxgy3mHAZAjDxd8YMt5gwGQAwyUYGDKZx4BJE8MlWBgymcWASQPDJZgYMpnDgEkRwyXYGDKZwYBJAcMlHBgy6WPAJInhEi4MmfQwYJLAcAknhkzqGDAJYriEG0MmNQyYBDBcCGDIpIIB0wSGCx2IIZMcBkwjGC50OAyZxDFgjoDhQo1hyCSGAXMYDBdKBEOmaQyYQzBcKBkMmcYxYA7AcKFUMGSOjAETx3ChdDBkDo8BA4YLZQZD5odCHzAMF8okhszBQh0wDBfyAkPmf4Q2YBgu5CWGTEwoA4bhQtnAkAlhwDBcKJvCHjKhChiGC+VCmEMmNAHDcKFcCmvIhCJgGC6UD8IYMoEPGIYL5ZOwhUygA4bhQvkoTCET2IBhuFA+C0vIBDJgGC7kB2EImcAFDMOF/CToIROogGG4kB8FOWQCEzAMF/KzoIZMIAKG4UJBEMSQ8X3AMFwoSIIWMr4OGIZL9ohIBEBF/OuKHC8n0IIUMr4NGIZL1v0zgD/Hv/5MRHrncjFBF5SQ8WXAMFxy4msAe+Nf1wPYlMO1hEIQQkZUNddrSArDJTdE5BgA2wAUAKgG0EL99ubxKRHpCWABgG2qekKu15MMXzUYhkvuqOrfAfwj/tsVDJfs8XOT8U3AMFzywpL4Pz/K6SpCyK8h44uAYbjkjYZgWZTTVYTUISGzR0Q+EZErc72uxmQlYESkVES6isipItI6yW0ZLnlARAoANPy3/1ZEJJfrCat4yPwSwFGIhc10EUn4uIyInCAiFSLSRUSO9mqdDTwJGBFpLSLPRyKRTcaYegB/B/AlgC8AfGeMiUYikS0i8oqInNTIOAyXHJGYYY7jLHUcZy+AOgAvxf96HgDrOM4uY8wCEbkodysNpb6IfV84AIoAvHykwI//UH8nEonsEBEXwDcAlgNYBWC3MWa/4zjrROQxESnO9EIzehZJRDobY/7NWnv60Ucf7fbo0cM588wzceaZZ6JZs2YAgLq6Onz++edYvHgxFi1a5FZVVTnGmNXW2mtUdckBYzFcciD+Rp0mIjeISEFFRQW6d+8uvXr1QllZ2fev27BhAxYsWIAlS5bYNWvWiIjstdY+oqoTc7f6cBCRUgDPINZkmiN2+cAtqvryAa8ZbIyZZq09qayszO3Ro4fTo0cPVFRUwJhYr6iursbChQuxZMkSLF682N2/f78Rkf+nqper6s6MrDVTASMijwAYXV5errfeeqs59dRTE9qusrISU6dOtWvXrhUAzwO4CbFmxXDJMhHpaYx533GcFldddZW57LLLEIlEmtyupqYGf/rTn/Duu+8qgE3W2nNVdYPnCw45ERmEWKtsgVjD7Ajg7yLytqoOOPPMM/W3v/2tHH/88QmNN2/ePDz33HNuVVWVqurtqvpM2mtMN2BExDHGLAfw45tuukkGDRqU0jgff/wxJk+ebK2131hr28b/mOGSJSIyEsC/nHbaaXbs2LGmoXEmo6qqCvfdd5+7adMmAXCpqr6d8YXSQeJtZhqAKwHME5FuxcXFze+//34n0R/yh3rhhRfw1ltvKYB3rbWD01pfOgETD5fVhYWFHaZNm2batGmTzlpQXV2Nm2++2e7Zs8eo6smq+lVaA1JCROQuAI/eeOONuPjii9Meb+rUqfr+++8DwCUMmewQkX8WkRfLy8vtU0895STSPBuzZs0ajB49Wq21f7HWnpfyutIJGMdx/qugoOCc6dOnm1atWqU8zoFqamrwm9/8xq2pqVnpuu5pGRmUjkhEBgD48KabbsLgwWn9sDpIPGQUwCncXfKWiIgx5rsTTzyx5R//+Een4RhLuiorK3HnnXeqtXaKqo5IZYyUVyIil1tr+z7yyCMZCxcAaN68OSZPnuyo6k9EZETGBqYfEJECEXnz7LPP1kyGCwAMHz5cysvL1RgzN6MD0+H8q+M4rSZPnpyxcAGAjh074vbbbxcAd4hISvtbKa1GRBwRebl///7asWPHVIZoVFlZGa644goRkcdFpGXGJyAAgIjMLCoqKh49erQn17RMnDjRAdBORO71YnwCRKQbgOtGjRqV0nGzpvTv3x+dO3e2xpg5qWyfatzd7ThOQTzdPHHttdeiefPmAPCwV3OEmYgUqeqlw4cPN+nurx9JaWkpBg0aJMaYuz2ZgCAi/7u8vNzt1auXZ3OMGzfOsdaWiUi/ZLdNKWCMMb/v06ePyWQdO5yBAwc6xpjrPZ0kvO4vLCzUvn37ejrJtddeC1U9WkQyuw9GEJESVf3Zdddd53g5T4sWLdCxY0drjHks2W2TTggRaW+tbTN06NBkN03a1VdfDVUtFhFvvwtCyBgztE+fPp6+MYHYMbVOnTpZEbnH67lCaERRUZHt2bOn5xNdd911xlp7RrK3iKTSjS8vKipyW7VqldCb8/3338dXX32F4cOHAwBefvllfPvtt7jzzjub3LawsBAtWrRwd+3adQkAHizMIGvtceecc05Cr/36668xfvx4PPfcczDGwHVdjB8/Hn379kW/fk235u7du5vKysqu6a6ZfqB/u3btEv6Gf/LJJ7F+/frvf19VVYVoNIqZM2c2uW23bt0AQACcBWBhonOmso/T98QTT0z8xX37YvHixaiursaiRYvwt7/9DbfddlvC23fo0MHh4xmTFz8Q309Eyg7zd50BmIqKxB6tW15ejrZt22Lx4sUAYj8k2rZtm1C4AECvXr1grW3BGySTF78x8dT4rTMHcRznJz/+8Y8T/h4eMWIEpkyZgilTpuCee+6B4zgYMWJEQtsaY3DUUUdFAVyS6HxACg3GGNOxbdu2CVfrZs2a4ec//zlefvllLFmyBA8++CCKiooSnq99+/ayYsWK8mTXSbgYwEwAKiI1AJYh9riFxQBOKSgocCORSML/Hy+++GLMmjUL0WgUq1atwsSJid9y1KFDh4Yv2wHYkPCGBMTaggIoEJG1iN1oOh/AEhFp2aVLl6QH3LVrF8aPH48rr7wSZ511VsLbHXvssaa6ujqp09VJB4yIFCYTEEDsVNfw4cNx33334YQTknviX1FREVT1WBHhE9RSVwTg3PgvAIAxxiYzQPfu3TFjxgy89NJLeOSRRxK6R+kwvmKJScup8V83x3+vyZ6a3rdvHx588EH07t0bAwcOTGrbwsJCQey9lLCk3yWqWldfX5/UNq+99hpatmwJ13WTnQ7RaBQi8nfE9v0ocZcAmAhgH2J33G5FrL18DKDMWntXsgN26dIFp5xyCkpLS1NdUz/EHh5OiVuD2DOQCxB72PpyxP4ffgrg9bq6uoS/h13XxWOPPYa2bdvi2muvTXoh9fX1ithNlQlLOmCstVu2bt16ChI8fvPWW2+hvr4ed999N1555RUke75+06ZNUNWtqlqZ7FrDTESmAdgNYCWApapafcDf/VN9ff0Yay2SudTg66+/Rv/+/ZNey7Zt2xq+nKeqyf10Crn4s3YcxB5XuuXAZyFHIpHqdevWtUr0e2ratGmIRqO4/fbbU1rLzp07LYB1yWyTykHejzZv3pxQvV62bBnmzJmDESNGoKKiArW1tQcdxU7E+vXro9ba+SmsM9RUda+q/h9VnXdguMR9BkBXr16d1JibNm1Cu3btkl7L/PnzYYzZy3BJnqr+h6rOUtXNhz5o3Vq76osvvkjo0MHMmTOxbt06jBkzBo6T2tUJu3fvdgC8k8w2qQTMGzU1NZF9+/Y1+qJvv/0WTz/9NMaMGdNwRS4GDx6MWbNmJTyRtRZVVVUOgPdSWCcdgaqqMabqk08+SXibHTt2oKSkBMXFyT/07LPPPlMAbKAZpqofb9iwocnjDtu3b8err76K3bt34+6778Ydd9yBO+64A5MmTUp4rsrKSqiqIMnLRVK6m9pxnD2DBw8+6sYbb0x622TMmjULzz//fL2qFvFjMjJLRJ4tKSn5zWuvvebpxXbRaBS/+tWv1Fo7TFWf93KusBGRYwF8O2nSJHTt6u1lRvfee6/94osv1ruu2ymZ7VK61t9aO3327NnJH7FN0htvvBFV1bcYLp4YvXfvXrN06VJPJ3nzzTehqnUAZng6UQip6g5jzIoZM2YkdUYwWXV1dVixYoVYa8clu22qNxPdV1tbK4lcAZiq2bNnN+wejfRskhBT1V3GmPlPPPGEZz8o9u3bh9dee82q6gz+kPCGtXbk6tWrTbLH05Lx+OOPWxGpVtVXk902pYBR1b0A7p45c6Zu3749lSEaVV1djalTpyqAp/jITO9Yawf94x//0Oeee86T8R944AHruu5uALd6MgFBVeeIyNxx48a51ma+yCxduhQLFiww1tpfp7J9yrdDq+rjIrJ61KhRbl1dUqfGG2WtxciRI11V3Q62F0+p6i5Vvfntt9/+/jaATJk1a1ZDrR7I9uItVR1UU1NTP378+IwmzM6dO/HAAw/Y+EPEZ6cyRlrPW7DW9ti9e/eem2++2W3qrFIi4ufo3W3bttVZa7vzjek9VZ0hItMffPBBzVTIvPPOO5g+fToA3Bn/oDDykKrWqmqfzz//HBMnTsxIyHz33Xe45ZZbrOu661T1klTHycSnCrQ2xqwuKioqHTt2bMI30B1q/fr1uO+++9zq6up9qnqqqm5Ma2GUFGPMc6p640UXXYRhw4alNEY0GsWkSZPswoULDYC7VPXxzK6SGiMivUTkoxNOOEH+8Ic/OK1bJ/Uhqt+bO3cupkyZYq21a6y1p6Vz/VJGPhcp/mzXWao6sGfPnnrbbbdJy5aJPemypqYGzzzzjM6dO1eMMfOttReo6p60F0VJE5GhIvJsq1atZOTIkc7pp5+e8LZz587F1KlT7f79+2tV9RJVTekRi5QeESk3xsxV1ZMvu+wyGTJkCAoLCxPadsuWLZg8ebJdvXq1IPbBbremuxeR6U92HGSM+VdrbZv27du7AwYMcHr16oVDk3TXrl2YP38+5syZY9esWWOMMbuttcNV9ZWMLYZSIiItjTH/Ya39WUlJie3Xr5/Tu3dvdOnS5aDbCqLRKJYvX46PPvoI8+bNc+vq6hwANQBa8Yrd3BORd0TklwC0W7du2q9fP3PWWWd9f9Frgy1btuCvf/0rZs+e7W7fvt0xxmy01l6sqssysg4vDnOISG8ReVRETrfWFhtjNBKJWABwXVdc1zXGmP0AVllr/5eqfpDxRVBa4g9bn2SMudxaWwpACgoK3PgDpyQajRoA1hjzrbX2eQD/DeBlAMep6o5crp2A+NMHlgGYbIy5R1U7qGqB4zg2EomoqqK+vt6oqsRv4/hEVX+vqv+d0XV4fRw1/oHa5wEoQ+yg8rcAPoyf6iafEJEzAJwBoBixz0Ker6prDnmNAtimqsk9k4MySkRuBzAFQIsDDzfEr/z9BYDWiN0VvRHAX1TVs2uhPA8YCg8RGQLgFbDF5FRDe1HV03O+FgYMZRJbTG4dqb3kirefO0JhdA2Asngdp+ybglh7yXm4AGww5AG2mNzIt/YCsMGQN9hiciOv2gvABkMeYYvJrnxsLwAbDHmHLSa78q69AGww5CG2mOzI1/YCsMGQt9hisiMv2wvABkMeY4vxVj63F4ANhrzHFuOtvG0vABsMZQFbjDfyvb0AbDCUHWwx3sjr9gKwwVCWsMVklh/aC8AGQ9nDFpNZed9eADYYyiK2mMzwS3sB2GAou9hiMsMX7QVgg6EsY4tJj5/aC8AGQ9nHFpMe37QXgA2GcoAtJjV+ay8AGwzlBltManzVXgA2GMoRtpjk+LG9AGwwlDtsMcnxXXsB2GAoh9hiEuPX9gKwwVBuscUkxpftBWCDoRxji2mcn9sLwAZDuccW0zjftheADYbyAFvM4fm9vQBsMJQf2GIOz9ftBWCDoTzBFnOwILQXgA2G8gdbzMF8314ANhjKI2wxMUFpLwAbDOUXtpiYQLQXgA2G8kzYW0yQ2gvABkP5J+wtJjDtBWCDoTwU1hYTtPYCsMFQfgpriwlUewHYYChPha3FBLG9AGwwlL/C1mIC114ANhjKY2FpMUFtLwAbDOW3sLSYQLYXgA2G8lzQW0yQ2wvABkP5L+gtJrDtBWCDIR8IaosJensB2GDIH4LaYgLdXgA2GPKJoLWYMLQXgA2G/CNoLSbw7QVggyEfCUqLCUt7AdhgyF+C0mJC0V4ANhjyGb+3mDC1F4ANhvzH7y0mNO0FYIMhH/JriwlbewHYYMif/NpiQtVeADYY8im/tZgwtheADYb8y28tJnTtBWCDIR/zS4sJa3sB2GDI3/zSYkLZXgA2GPK5fG8xYW4vABsM+V++t5jQtheADYYCIF9bTNjbC8AGQ8GQry0m1O0FYIOhgMi3FsP2EsMGQ0GRby0m9O0FYIOhAMmXFsP28j/YYChI8qXFsL3EscFQoOS6xbC9HIwNhoIm1y2G7eUAbDAUOLlqMWwvP8QGQ0GUqxbD9nIINhgKpGy3GLaXw2ODoaDKdothezkMNhgKrGy1GLaXI2ODoSDLVothezkCNhgKNK9bDNtL49hgKOi8bjFsL41gg6HA86rFsL00jQ2GwsCrFsP20gQ2GAqFhhYDoKOq7k1zrBIAQ8H20iQ2GAqLOwGUAdgV37VJiYj8E4DdiIXLDoZL4xgwFHgiYgA8CMACcAAMSWO4QQBcAArgWBE5I/0VBhcDhgJPVS2APwHYF/+j7iJSnOJwlwEoQCysvgSwMv0VBhcDhsLidgCVAKKIBU3vZAcQkRYAusZ/uxfAhaq6P2MrDCAGDIWCqtYBuBBANYCjEdvVSda5iO0e1QL4tapuytgCA4oBQ6GhqlsQ28UBgIsAQERKReRWEXkrEolsiEQiO+O/vhGRv4jIWBE5Ob7NLwE0AzBJVf8zB/8KvsPT1BQ6InIPgInGmC+stT8pKChwy8rK0KlTJ6dNmzYQEezZswdr1661mzdvtrW1tRFjzGZrbQmAzwAMiB/XoSYwYChURKTUGDPHWtu9Y8eO7tChQ52KiopGt9m+fTumT5+uixYtAoAqa+0AVV2SlQX7HAOGQkNELheRV0tLS+Whhx5yysvLk9q+pqYGEyZMsCtXrhQAT6nq771ZaXAwYCgURGQYgGcvvPBC3HLLLZLOWHPnzsXkyZMVwGvW2nSuqQk8BgwFnohcBODfhwwZIldffXVGxlyxYgXuvfdeVdU/qmrKVwYHHc8iUaCJSLGI/Llv377IVLgAQEVFBUaPHi0AbhORXhkbOGAYMBRoxpj3SkpKIiNGjEhrt+hwevfujdNPP92KyHsikvHxg4ABQ4ElIh2ttedOmDDBMcabt/rYsWON4zhHA7jHkwl8jgFDgSUiTx533HHRzp07ezZHYWEh+vTpY4wxIzybxMcYMBRI8V2WC6644oqI13MNHToU1to2InKm13P5DQOGgqq7qjrnnXee5xO1atUKLVu2jAK43vPJfMbzdCfKkctKSkqixpiE3uMffvgh3n33XQDAhg0b0L59ewCxs0XDhg1rcvsOHTo4y5Yt49mkQzBgKKh6nHTSSU6iLz7//PNx/vnnY+fOnbjrrrswZcqUpCbr2rWrLF++/JSkVxlw3EWiQBKRVkcddVTSp443btyIdu3aJT1faWkpABQmvWHAMWAoqFK6RH3Dhg0pBQyviD88BgwFkqpW7d69O+lHKmzcuPH74y/J2LlzJ0SET7c7BAOGguqv33zzTUoBk0qDWblypVpr1yS9YcAxYCio3qipqYlYm3jGWGvxzTffoG3btklPtmHDBldV5yW9YcAxYCiQVPVLEal/5513Et5m69ataNOmDQoKCpKaa8eOHdizZ08EwIwklxl4fFwDBZaI/FtpaellL730kqeXYzz88MO6cOHCra7rnuTlPH7EBkNBNrKqqiry+eefezZBdXU1Fi5cCGvtw55N4mNsMBRoxpi3ioqKLpo5c6aJRDJfZEaOHOmuW7duh7X2ROU30w+wwVCgqeqv9+/fX/PQQw9l/FMA3nvvPaxdu9ZYa/szXA6PAUOBpqquql64ZMkSTJ06NWMhMG/ePEybNg0Axqnql5kaN2i4i0Sh0PBc3rPOOgtjxoyRdHaX3nzzTbzwwgsA8ISq3pmpNQYRA4ZCQ0TOEZH3mzVrVjh27NgmPw/pUFVVVbjnnnvczZs3C4AxqvqYNysNDgYMhYqIFIjILFUdePzxx7tXXXWV079//0a3WblyJWbMmGFXr15tjDHrrbX9VHVjlpbsawwYCiUR6SwiT6rqABExLVu2dNu3b++UlpaKiGDv3r346quvojt37jSu6xpjzDJr7RhV/SDXa/cTBgyFmog4AM4D8EsR6WmMaQXAqOpea+1SAB8AeFtV9+Rwmb7FgCEiz/A0NRF5hgFDRJ5hwBCRZxgwROQZBgwReYYBQ0SeYcAQkWcYMETkGQYMEXmGAUNEnmHAEJFn/j8LEjxAD9R2fAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 266.457x266.457 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig_1_8 = DAG(\n",
    "    ebunch=[\n",
    "        (\"W\", \"Y\"),\n",
    "        (\"W\", \"Z\"),\n",
    "        (\"X\", \"W\"),\n",
    "        (\"X\", \"Y\"),\n",
    "        (\"Y\", \"T\"),\n",
    "        (\"Y\", \"Z\"),\n",
    "        (\"Z\", \"T\"),\n",
    "    ]\n",
    ")\n",
    "\n",
    "fig_1_8_node_pos = {\"T\": (3, 0), \"W\": (2, 4), \"X\": (0, 2), \"Y\": (2, 2), \"Z\": (4, 2)}\n",
    "\n",
    "fig_1_8.to_daft(node_pos=fig_1_8_node_pos).render();"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Study question 1.4.1\n",
    "\n",
    "Consider the graph shown in Figure 1.8:\n",
    "\n",
    "(a) Name all of the parents of Z."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['W', 'Y']"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig_1_8.get_parents(node=\"Z\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(b) Name all the ancestors of Z."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'W', 'X', 'Y'}"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nx.ancestors(fig_1_8, source=\"Z\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(c) Name all the children of W."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['Y', 'Z']"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig_1_8.get_children(node=\"W\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(d) Name all the descendants of W."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'T', 'Y', 'Z'}"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nx.descendants(fig_1_8, source=\"W\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(e) Draw all (simple) paths between X and T (i.e., no node should appear more than once)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[['X', 'W', 'Y', 'T'],\n",
       " ['X', 'W', 'Y', 'Z', 'T'],\n",
       " ['X', 'W', 'Z', 'Y', 'T'],\n",
       " ['X', 'W', 'Z', 'T'],\n",
       " ['X', 'Y', 'W', 'Z', 'T'],\n",
       " ['X', 'Y', 'T'],\n",
       " ['X', 'Y', 'Z', 'T']]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list(nx.all_simple_paths(fig_1_8.to_undirected(), source=\"X\", target=\"T\"))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(f) Draw all the directed paths between X and T."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[['X', 'W', 'Y', 'T'],\n",
       " ['X', 'W', 'Y', 'Z', 'T'],\n",
       " ['X', 'W', 'Z', 'T'],\n",
       " ['X', 'Y', 'T'],\n",
       " ['X', 'Y', 'Z', 'T']]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list(nx.all_simple_paths(fig_1_8, source=\"X\", target=\"T\"))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Study question 1.5.1\n",
    "\n",
    "(a) Draw the graph that complies with the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARgAAACnCAYAAADKZUOYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAAdIklEQVR4nO3deXgVVbY28HftykBCBqLQApIGZBBFRGlAieKQi3zSNA54HVC0FRAkiFfRlgYZvEQUFMGGBgQEB1BQPrVREZuLchWMSERQ6OYSAkQQmQ2BEDKcWuv74wRau8nJqUptuH5Zv+fx8RyovffiqXPeGk7VLhIRKKWUDeZMF6CU+v+XBoxSyhoNGKWUNRowSilrNGCUUtZowCilrNGAUUpZowGjlLJGA0YpZY0GjFLKmpjTMQgRpQFoBIAA7BWRQ6djXBUcIooF0BRAXQBHAewQvc/kF4eIGgGoD6ACwG4ROWpzPCt7MER0NhHNjYmJ2WmMqQDwI4C/AdgE4KAxJhQTE7ObiF4nonNt1KBqhsLudxxng+M4xwCUA9gKYAOAbQDYcZwiY8wXRHTDmaxVVY2I2hLR+zExMQeIyAXwA4BvAWwGcMQYU+Y4zjYieo6IEgIfP8iNEBG1Nsa8ycyXJCcnu507d3Y6duyIjh07ok6dOgCA8vJyrF+/Hrm5uVi7dq1bWFjoGGO2MPNdIrIusGKUL0REAGYS0X1EFNuuXTt06NCBMjIy0LBhw5PLFRQU4IsvvsC6des4Ly+PiOgYM08QkfFnrnp1AhH1MsbMZOZzGzZs6Hbu3Nnp3Lkz2rVrB2PC+xXFxcVYs2YN1q1bh9zcXLesrMwQ0X+LyK1BHWUEFjBENAHA4+np6TJkyBDTtm3bqNrl5+djxowZvHXrVgIwF8BA3fU+M4jocmPMMsdxUu644w7Tu3dvxMRUfxRdUlKCBQsWYOnSpQJgJzNfIyIF1gtW/4KI4onoPRHp3rFjR3nggQfonHPOiart6tWrMWfOHLewsFBEZKiIvFjjemr6XSYixxjzLYALBg4cSD179vTVz2effYYpU6YwM//AzG1E5FiNClOeENEwAM9ffPHFPHr0aHNij9OLwsJCjBo1yt25cycBuFlE3gu8UFUlImpKRBsTEhISx4wZ40S7kf9nL7/8Mt59910BsJSZe9WoppoETGW4bImLi2s+c+ZMU79+/ZrUguLiYgwaNMgtLi7+kZnTRaSsRh2qqBDRHwA8O2DAANx444017m/GjBmybNkyALhJQ+b0IKJziSg/PT099k9/+pMTzZ5nJHl5eXj88ceFmT9m5ut811WTgHEc55PY2NirX3rpJVOvXj3f/fxUSUkJ+vfv75aUlPzddd2LA+lUVYmIugP468CBA9GrV402Vj9TGTICoIUeLtlFRGSMOdi4cePU6dOnOyfOsdRUfn4+Hn30UWHmqSLysJ8+fFdCRLcy87UTJkwILFwAIDExEVOmTHFE5CIiejiwjtW/IKJYInqnS5cuEmS4AEBWVhalp6eLMWZloB2rU3nFcZx6U6ZMCSxcAKBly5YYOnQoAXiIiHwdb/mqhogcIprfrVs3admypZ8uImrYsCFuu+02IqJJRJQa+AAKAEBEC+Pj4xMef/xxstH/+PHjHQBNiegJG/0rgIjaA7jnscce83XerDrdunVD69at2Rizwk97v3E33HGc2Mp0s6Jv375ITEwEgGdsjVGbEVG8iNyclZVlanq8XpW0tDT07NmTjDHDrQygQER/Tk9PdzMyMqyNMXbsWIeZGxJRpte2vgLGGPNI165dTZC7Y6fSo0cPxxjze6uD1F5j4uLi5Nprr7U6SN++fSEiyUQU7DGYAhHVFZEr7rnnHsfmOCkpKWjZsiUbY57z2tZzQhBRM2au369fv4jLvfzyy3jnnXdOvhcR3HHHHSgtLY16rD59+kBEEojI7regFjLG9OvatWvED2YQ6zAxMRGtWrViIhrpv1pVhYfj4+P58ssvj7hQEOvxnnvuMcx8aeWFmFHzswtya3x8vFvdid2CggI0a9bs5Pt9+/YhJSUFXo4T4+LikJKS4gK4yUedKgJm/tXVV18dcZkg1iEAdOjQwRDRhX7qVBF1a9q0abVf+CDWY/v27YHwvYSXeSnQT8Bc27hx42oX+u6779C8efOT73fs2PGzf2S0mjdv7hDRlZ4bqioRUWsApl27dhGXC2odZmRkgJlTvG79VGSO41x0wQUXVPsdDmI9GmOQlJQUgseNveeAMca0bNKkScRd6+LiYriui7S0tJN/VlBQ8LN/ZLSaNWtGxph0zw1VJBmxsbFupJO7Qa7Dn7Rp6rmxqhIzp7Zp0ybiMkGuxwYNGhgAnn6u9nMOJi4+Pj7iMgUFBWjatOm//NmJ1JwwYQJc18Xnn3+OxYsXR+yrcqzTMq1ELZJU3Qn6SOswJycHS5cuBQAsXrwYX375ZbTjJvuoVVWt2p+mI63HoqIiTJ8+HdOnT8fYsWMxevToiH3FxcURgMhf/n8u0MvCACAi5RUVFRGX2b1798/uvGVmbN68GSfujcjMzMTs2bOxadMm3HrrrRH7CoVCABDyWqeK6DgzR1wg0jrs0qUL1q9fjxUrViApKQmXXRb1YbneXxYsKS8vj7hApPWYmpqKIUOGoHfv3khKSsKIESMi9lVRUSEIT9sRNc8Bw8y79+zZE/HT2aBBA2zduhUn/vGLFy9G27ZtkZKSAgBo3LgxPvvsM9x3333Vjrdz506IyH6vdaqINlRUVDiRQibSOiQidOrUCdu3b0ePHj2qHWzv3r0nXu4KoHZVyRhTvG3btojLVPdd3L17NxYsWIAHH3zwxHVnVTp06BAjPBdQ1Pwcenz6/fffX4kI4dShQwesX78egwcPRmxsLM4//3wMGTIEQPiO24ULF6Jfv35YvXo1MjMjX7uzffv2EDPn+KhTVe1rALJlyxa64IILTrlApHUIAMeOHUOXLl2iGiwnJwfGmGOu60be9VWeMPPmTZs2XY7wrzunFGk97tq1C4sWLcLQoUOj+kXpyJEjDoD3vdToJ2DeLikpGVtaWhqxqP79+6N///4/+7Pjx49j9uzZeOCBB5CcnIxnnnkG11xzDao6H8DMKCwsdAB86KNOVQUREcdxCletWnVWVQEDnHodnrB9+3Zcf/31UY339ddfC4B8P7WqqonIZwUFBZ1Qzff4VOuxpKQEo0aNQvv27TF37tyTy1X1nc7Pz4eIEABP95b5upvacZyjvXr1ShowYIDntl4sWbIEc+fOrRCReJ2EKlhENKtu3br9Fy1aZPUq0FAohFtuuUWY+X4RmWtzrNqGiBoA2D9x4kRceKHdy4yeeOIJ3rRp03bXdVt5aefrWn9mfmn58uWun7ZevP322yEReVfDxYrHjx07ZjZs2GB1kHfeeQciUg5gntWBaiEROWCM2Thv3rzIZ+xrqLy8HBs3biRmHuu1rd+biUYdP36cFi5c6LN59ZYvX37i8GiYtUFqMREpMsbkTJ482dqGorS0FIsWLWIRmacbCTuYediWLVvMli1brI0xadIkJqJiEXnDa1tfAVM5neXwhQsXyr59+/x0EVFxcTFmzJghAP4kIrsDH0ABAJi55+HDh2XOnDlW+h83bhy7rnsEwJBqF1a+iMgKIlo5duxYt7pLD/zYsGEDvvjiC8PM/+6nve/boUVkEhFteeyxx9zqfov3gpkxbNgwV0T2QfderBKRIhEZ9N577yE3NzfQvpcsWXJit7qH7r3YJSI9S0pKKp588slAE+bQoUMYN24cV04ivtxPHzWab4GZOx85cuTooEGDXC93ZlYlFAph6NCh7t69e8uZuYN+MO0TkXlE9FJ2drYEFTLvv/8+XnrpJQB4VETWBNKpqpKIHBeRruvXr8f48eMDCZmDBw9i8ODB7LruNhG5yW8/QTxV4GxjzJb4+Pi00aNHV3sDXVW2b9+OUaNGucXFxaUi0lZEvqtRYcoTY8wcERlwww034P777/fVRygUwsSJE3nNmjUGwB9EZFKwVapIiCiDiD5t1KgRPf30087ZZ5/tq5+VK1di6tSpzMx5zHyxiPi+fimQ5yJVzu26RER6XH755fLggw9Samp0M12WlJTgxRdflJUrV5IxJoeZr7f9OEt1akTUj4hm1atXj4YNG+ZccsklUbdduXIlZsyYwWVlZcdF5CYR8TXFoqoZIko3xqwUkfN69+5Nd955J+Li4qJqu3v3bkyZMoW3bNlCAF4EMKSmRxFBP9mxpzHmFWau36xZM7d79+5ORkYG/jlJi4qKkJOTgxUrVnBeXp4xxhxh5iwReT2wYpQvRJRqjPmAma+oW7cuZ2ZmOldeeSXatGnzswsiQ6EQvv32W3z66adYvXq1W15eTkT0roj0qckWTwWDiJ4golEA4tu3by+ZmZnmsssu+5fbAXbv3o3PP/8cy5cvd/ft2+cYY75j5htF5JtA6rBxmoOIriSiZ4noEmZOMMZITEyMAIDrunBd1xhjygBsZuYRIvJR4EWoGqmcbH2iMeZWZk4DQLGxsa4xhlzXRSgUMgDYGLOfmecC+E8Nlv99iOj3xpiRItJcRGIdx2HHcQBAKioqjIiQMeaYiKwSkUdE5H8CHd/2edTKB2pfB2ABwrfr3wLgr/rkxl8WIroUwP9BeBL2+wDkiEjema1KeVF55e+/AVgIIAfAfwL4WESsXQtlPWBODkS0H0CDyvsZ1C8QEbUEsFXX4S8bEQmAV0XkXttj2X0sgFKqVtOAUUpZowGjlLJGA0YpZY0GjFLKGg0YpZQ1GjBKKWs0YJRS1mjAKKWs0YBRSlmjAaOUskYDRilljQaMUsoaDRillDUaMEopazRglFLWaMAopazRgFFKWaMBo5SyRgNGKWWNBoxSyhoNGKWUNRowSilrNGCUUtZowCilrNGAUUpZowGjlLJGA0YpZY0GjFLKGg0YpZQ1GjBKKWs0YJRS1mjAKKWs0YBRSlmjAaOUskYDRilljQaMUsoaDRillDUaMEopazRglFLWaMAopazRgFFKWaMBo5SyRgNGKWWNBoxSyhoNGKWUNRowSilrNGCUUtZowCilrNGAUUpZowGjlLJGA0YpZY0GjFLKGg0YpZQ1GjBKKWs0YJRS1mjAKKWs0YBRSlmjAaOUskYDRilljQaMUsoaDRillDUaMEopazRglFLWaMAopazRgFFKWaMBo5SyRgNGKWWNBoxSyhoNGKWUNRowSilrNGCUUtZowCilrNGAUUpZowGjlLLmtAQMEc0E0KDy9VunY0wVLCLqCmBV5evNRHTeGS5JeUREcUT0WeXb3xPRcNtjnq49mA4/ef2b0zSmCta5AFIqXzcFUO/MlaJ8CgHo/JP3rWwPeLoC5jMAXPn6y9M0pgrWOgBS+ToGwKYzWIvyQUQYwJbKt8cA5Nge83QFzJcAigGUIBw26pdnG/7xeSkQkfIzWYzybVXl/xnAV7YHO10Bsw7hrV6o8rX6hanc+m2ufGt9y6es+QLhvZc6AP5ue7DTFTAFCO9e1wWw8TSNqYJ3Yuv3+RmtQtXEOoS/h9tEJGR7sNMSMCIiCKfldyJSejrGVFasqfy/7oX+cuUBKMNp2khQ+LtvoWOiBgD6Ash0HKc9gPoIB9pB13X/BuATAG+IyC4rBagaIyICcCOA3xljuhBRIxFJJaIfmXmXiHwO4B0RWXmGS1URENGFAO4AcLXjOK0R/i6WisgBZl4H4L8AvC4ixwIfO+iAIaKbjDFPM/MF8fHxbsOGDdG6dWvnrLPOAgAcPnwYW7du5R9++EFKS0sdY8x2Zn5SROYHWojyjYjOBfACEd1ERE5aWprbokWLmPT0dMTGxqK0tBQFBQVSUFDgFhUVOURUysyvARguIkVnun51cuPwqDHmD8z8q8TExFB6errTsmVLSkpKgohg//79yMvLC+3fv9+EQiEioi9F5FERCewcW2ABQ0TnGmM+EZHWF154IQ8YMMC0bNkyYptdu3Zhzpw5vGHDBiKiH5g5U0TyAilIeVb5oZwC4KGzzjrL7d27d0yvXr1gTNVH0uXl5Vi4cCGWLVvmHjt2DAAeE5EXTk/F6lSIqKsx5gMiSrriiiuoX79+dPbZZ0ds89VXX+HVV191CwoKHGPMKmbuEcQeTSABQ0QPAJjeqFEjGT9+vNOgQQNP7YuKijB69Gh3x44dBsBoERlf46KUJ0TUyBjzNRGdk5WVRd27d/fcx4IFC/DWW28JEf2NmTuKSJmFUlUExpiXReTejh078ogRI0xcXJyn9vn5+RgzZoxbXFzsisiNIvJRTeqpccAQ0R8BPHP77bejb9++Nerrgw8+wKxZswBgqoj8R406U1EjoqZEtLlhw4ZxkydPdpKSknz3tXfvXjzyyCNuSUnJj8zc3MZxvTo1x3H+CuC6UaNGUadOnXz3w8yYPHmyfPrppwBwq4i87bevGgUMEfUH8NLgwYPx29/+1nc/P5WTk4NnnnkGAEbpnox9RJRMRHt+/etf15k6daoT6XAoWiUlJbj//vvd4uLiPcz8a7H1S4I6yRiziIhumzJlCp13XjC3ic2ePVvef/99ALhaRFZVt/wp6/I7OBE1IKJZt9xyS2DhAgAZGRkYNGgQAGQTUevAOlanZIz5a3Jycp0XXnghkHABgMTERMyaNcsxxjQGMCuQTlWViKiXiNw+bty4wMIFAAYOHEi/+c1vhIiWEpHjpw/fnyhjzCf169fHvffe67eLKv3ud79D8+bN2RjzSeCdq5OI6C5m7pKdne3ExMQE2ndSUhIeeeQRA2AAEekNrpYQUTwRvdm1a1dp37594P2PHDnSxMXF1SWiRX7a+woYIurKzBeNHz/eV6pFIzs72xGRxkR0t60xajtjzIzMzEwJcqv3U1dddRVatWolxpjXrQygAGBSXFxc/LBhw8hG53FxcRg+fLgRkX+vvLbNE18BY4x5/rzzznMbNWrkp3lUUlNT0b59ezHGPGltkFqMiHowc8qgQYOsfDBPyMrKMsx8fuW1NSpgxpj7evToYYLeA/2pTp06ITk52QXwrNe2ngOGiBKYueO9995rbe/lhIEDBxpmPo+I0m2PVdsYY55u3bo1JyYmWh2nZcuWSEtL8/XhVJER0Q0iknjXXXdZH6tXr14OEd3htZ2fPZgbHMeRSy+9NKqFly1bhhkzZpx8P3/+fDz//PNRtU1PT0edOnVchC9zVsG64Lrrrotq/e/atQv9+/cHc3hKH9d1MXr0aHzySXSnyDp37uw4jnON70pVVe4855xzuE6dOlEt/MILL+Chhx46+d/dd9+NPn36RNX25ptvhojUISJPhy1+Aub6tLQ0rn6xsGuvvRa5ubkoLi7G2rVr8dVXX+HBBx+MerDKw7BMH3WqKhBRMjPHd+nSJarl09PT0aRJE+Tm5gIIbySaNGmCzMzoVkvHjh3BzOf4LlidkjGmU6tWraI+knj44YcxdepUTJ06FSNHjoTjOHj44YejalunTh3Exsa6AHp7qdHzgZsxpkPz5s2jblenTh1cddVVmD9/PtatW4fs7GzEx8dHPV7r1q2dnTt3tvNaZ21X+cvNhwjfPfvfAHIRvgv6BwDdHcfh1NTUqDcwN954I5YsWYJQKITNmzdj/PjoL1Hq0KEDRMQhorNF5JCnf0gtR0RrASQjfPdzDsLr8O8iUgHg3Isvvthzn0VFRXjyySdx++2347LLLou6XYMGDeSHH364CsD0aNt4DhgiSk5OTvbUplu3bsjKysKoUaPg9cRwvXr1ICINiUgv1vJOAPwKQBeEJxmKR/hW/f2O4+0UWocOHTBv3jy89tprmDBhArycVPzJ5eoHw7c7KR/aIHyqgAEkEFE+EcVWd4/RPystLUV2djauvPJK9OjRw1PbpKQkB0CqlzZ+Tj17/qIvWrQIqampcF3X+2AiICIG8KPnxrVbA4TXFQE4Mb2lA6AQwGH4WI9t2rRBixYtkJaW5remQ/jH3MwqOid+GnYRnhGyDsLTz+4D0NrLRdKu6+K5555DkyZNfN3WUzmWp8+N54ARkSNFRdHfkf/uu++ioqICw4cPx+uvv46MjAxP4/34448AsFdEfu2t0tqNiC4C8CbCMwiuQnjX+hsROU5EN4RCoXe99rlr1y5069bNcy2lpSfnGGslIoWeO6jFiGgpwjPQfQpgLYB1IrIXABzHOX7gwIHozvACmDlzJkKhEIYOHeqrlqNHjzLCG6ioeQ4YZs7dsWPHRdG0/eabb7BixQo899xzSExMxPHjx7F9+3Z4ubArLy/PdV13g9c6azsR2QSgbRV//V/MbA4dOgQvu9g7d+5E06ZNPdeydu1aEFGImTVcPBKRnhH+eue3337bumfPSIuELVy4ENu2bcPTTz8Nr4fHJxw8eJAAeJpczM+vSB8ePny42gr379+PadOm4Y9//CNOXGvRq1cvLFmyxNNg+/btA8Kz36mAiMhxY0xpTk708wodOHAAdevWRUJCgufx1q1bB2PMHs8NVUTMvDY/P7/aeXX37duHN954A0eOHMHw4cNP/kw9ceLEqMcqKSlBKBQyAN7xUqPnu6mJKBZA2YgRI8jr4Y5XeXl5ePTRRwHgVyJywOpgtYwx5oumTZt2mjZtmvULJu+880736NGjr4jIANtj1SZElAng44ULF6ImU2xE4+WXX8Zf/vKXY67rehrI8x6MiFQQ0ar58+d7P2Pr0dy5c9kY83cNl+CJyB8KCgocL+fT/Ni4cSOOHj3qALD+mNLaRkQ+McYUvfrqq9bH+uijj1xmfsVrO1/3IonIw99//72zfft2P82jcuDAAWzevNkw8whrg9RiIrLaGHPwz3/+s9Wf/2fOnOkaYzbo9S92MPP0jz/+mMvL7T0Hb+XKlSgpKTEAPH8X/QbMemNMzujRo90Tl48HbeTIkS4R5YnIe1YGUGDme9esWUMbN9p5VNUHH3yAXbt2GWa+1coACgDGhEKh4qeeesrKF7GkpARTp05lALNF5KjX9r7ng2Hm648ePcrTpk0LfAu4YMEC7N27l5j5mqD7Vv8gIkuJ6KPs7Gz+yU/JgThw4ADmzJkjACaKSH6gnauTRMQVkZ7r1683q1evDrz/MWPGMDP/CGCwn/a+A0ZEjorIbStWrMDrrwc33ceHH36IN998EwAGi4j+8mCZiNxQVlZWOHjwYDeo3ezCwkJkZWUxwpe06yGuZSKyGsC0Z599VjZs2BBYvxMmTJC8vDypfNqHrx2JICb97gfgpW7dumHo0KFUk2kXX3nlFbz99tsAMEJEJtSoMBU1IjrbGJOfkpKSPGnSJOecc/zfl5iXl4eRI0dyRUXFDmY+X0Ss/xigwowxrwHom5WVRddff73vfsrLyzF27FjetGkTAGSKyKd++wrqsSU3EdGbycnJTnZ2tuN1hrQ9e/bgiSeecA8ePAgRGSQic2tclPKEiJKNMWtF5Pw+ffpQtLfxn8DMmDZtmqxYsYKMMR8z83U62ffpR0RPARjZokULHjdunJOSkuKpfW5uLiZOnMjl5eXFIpIpIjV6THCQD16ra4xZzswZTZo0ce+++26nuutkvv76a7zyyivujh07HGPMRmb+N/1J+swioscATEhISED37t2dvn37ItJ8I4cPH8a8efNk1apV4rpuhYjcVZPHXKiaI6I2lQ9BbNi2bVvp379/xIcgMjOWLVuGxYsXu4cOHXKI6C0RuTOIvU8bj469lIimiMhVxhjUq1fPbdasmZOamkpEhCNHjqCgoCBUWFhoXNclIsqtfFxl8GeolC9ElADgKWPM/cycnJCQEGrSpIlp3LixiYmJQVlZGb7//nt3z549KCsrc4wxB5h5MoDn9JDofw8i6m2MGc/MbWJiYrh+/frcvHnzmISEBIgIDh06JAUFBW7ldUquiLwP4BER+S6wGmztxVZe8dsDQE9jTCciSqkcr7jygdvLAHwgIsetFKACQURtAdyC8IPT0wHEAih3XTcf4ftS/q+IFJzBElU1iCgNwM0Auhtj2hFRIgCXmQ+JyBcA3gOw0sYhrbWAUUqpYJ60pZRSp6ABo5SyRgNGKWWNBoxSyhoNGKWUNRowSilrNGCUUtZowCilrNGAUUpZowGjlLJGA0YpZc3/A7g2v6qvhq4cAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 266.457x153.071 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "graph = DAG(ebunch=[(\"U_X\", \"X\"), (\"U_Y\", \"Y\"), (\"U_Z\", \"Z\"), (\"X\", \"Y\"), (\"Y\", \"Z\")])\n",
    "\n",
    "graph_node_pos = {\n",
    "    \"U_X\": (0, 2),\n",
    "    \"U_Y\": (2, 2),\n",
    "    \"U_Z\": (4, 2),\n",
    "    \"X\": (0, 0),\n",
    "    \"Y\": (2, 0),\n",
    "    \"Z\": (4, 0),\n",
    "}\n",
    "\n",
    "graph.to_daft(node_pos=graph_node_pos).render();"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2 Graphical Models and Their Applications\n",
    "\n",
    "### 2.3 Colliders\n",
    "\n",
    "Figure 2.5. A directed graph for demonstrating conditional independence (error terms are not shown explicitly)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-04-14T16:37:28.267436Z",
     "start_time": "2023-04-14T16:37:26.179470Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAt4AAAA2CAYAAADwHJ6ZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAAfbklEQVR4nO3deXRUVfYv8O8+N5AEEsYAYYiEKSKjgCDQKEMziAjY2DZD42uEIIOKotIogyIJgzQNPBBaNIoyiEM7tDIIjfrzKSBji/ADGyKpCDIFQgghIZW6e78/UqHTKqFuUTe4VvZnLdci5J5zNlbVrX3P3fccEhEopZRSSiml3GVudABKKaWUUkqVBZp4K6WUUkopVQo08VZKKaWUUqoUaOKtlFJKKaVUKdDEWymllFJKqVKgibdSSimllFKlQBNvpZRSSimlSoEm3koppZRSSpUCTbyVUkoppZQqBZp4K6WUUkopVQo08VZKKaWUUqoUhLnVMRHVADAcQA/LsloDqOj/VZ5t2/8L4DMAb4rIMbdiUNeHiAjAQAD3GGM6EVFNABYAHzMfF5GtAN4Xkc9vaKCqRETUDMAQAF0ty0oAEA6ARSSbmfcA+CeANSJy6UbGqa6OiCoCGAaglzHmNiKqhMKJk3xmPiwiXwB4S0QO3tBAVYmIqCuAQUTUxRgTh8LvYFtEMph5O4D1AD4QEbmhgaqrIqK6KMxtuluW1QJApP9Xl2zb/haFuc0aETl9o2JUJSMiC8DvAfQ1xnQkouoozG0KmDndn9u8KyLbXBk/1J9vIrrXGDObmW8JDw+3Y2NjkZCQYFWrVg0AkJWVhSNHjvCJEyfk8uXLljHmKDPPEJFVIQ1EBc1/YllERPcSkVW1alW7UaNGYXFxcShXrhwuX74Mj8cjHo/HvnDhgkVEl5l5JYDJInLhRsevrlw0PWmMmcTMNStUqOCLi4uzGjduTFFRURARnDlzBocPH/adOXPG+Hw+IqIdIvKkWycb5RwRdSSiBSLSMSwsTGrWrMlNmjQJq1WrFogIly5dQmpqqhw7dsy+dOlSmDHmDDPPBzBfk7dfByKKBjDXGDNCRCpUqlTJ16BBAys+Pp4iIiJQUFCA48eP4/vvv/dlZmZaImKLyD8ATNSJqV8PIhpmjJnJzI0iIiLs2rVrIyEhwapSpQoAIDMzE0eOHLFPnjyJ/Px8yxjzHTNPFZH3b2zkqggRNSaiRQDuIiJTvXp1u3HjxmF16tRBuXLlkJeXh6NHj0p6erqdnZ0dZoy5yMyvAJgmInkhiyNU52YiqmuM+UxEEpo1a8aJiYmmcePGJbY5duwYXnnlFf7mm2+IiE4wcw8RORySgJRj/mRtIYAJ1apVswcNGhTWv39/GHP1iiSv14u1a9di48aN9qVLlwDgKRFZVDoRq19CRHcYY9YRUdRvfvMbGjlyJFWvXr3ENrt378Ybb7xhezweyxjzJTP31RnwG4eIIo0xG5m5a/369e0//elPVvv27Utsc/78ebz66qvy1VdfCTNfEpGBejfqxiKiR4hoUWRkJPr27WsNGzYM5cuXv+rxzIwNGzbg3XfftTMzMw2ApQAm6EXUjUNEjf25Tb3WrVtLYmKiqV+/foltUlNTkZKSwgcPHjRElOrPbfQi6gbx5zYrAPypRo0aviFDhoT17NmzxNwmNzcXa9aswebNm+38/HwWkYdE5PWQxBOKzzMRjQWwtHbt2jJr1iyrRo0ajtpfuHAB06dPt9PS0gyA6SIy67qDUo4QUW1jzF4iqjV+/Hjq3bu34z5Wr16Nd955R4jof5n5NhHJdyFUVQJjzAoRGXHbbbfxM888Y0r6kv8lqampePbZZ+2cnBzbn7h94lKo6iqIqCcRratYsWLY888/byUkJDhq7/V6MXfuXN61a5chotXM/IBLoaqrIKJwY8xOEWl5//330wMPOH8JtmzZghdffJH9ZSjtRORHF0JVJSCipwHMjo+P5+TkZKty5cqO2p89exZTp061T5w4QSi8gFrqSqDqqogowRjztWVZlZ966inTuXNnR+2ZGSkpKfLxxx+TMWY7M98hIvZ1xXS9ibf/jTln8ODBGD58+HX1tW7dOixfvhwAFovIY9fVmQoYEdUnokOxsbHlFyxYYEVFRQXd16lTpzBx4kQ7Nzc3k5kb6Kxp6bEsaxOAXtOmTaNrzY6WhJmxYMEC+eKLLwDgfhF5L1QxqpIR0b0A3u/SpQsmTZpEJc3IXMvu3buRlJQkAD63bfu3oYpRlcx/tyItMjIyZuHChVbt2rWD7is3NxdPPPGEfeLEiQIRaS4iR0MYqioBEc0H8GRiYiIGDhx4XX2tXbsWb775JgA8KyJJoYhPXRsRtSSiPQ0bNjRz5861IiIigu7L4/Fg0qRJ7PV6f2DmBBEpCDqu60m8iWgUgJRx48bh7rvvDrqf4rZt24Y5c+YAhTU1OvPtMiKKJqKTN910U8TixYut6/miL5Kbm4vRo0fbOTk5J5n5Jr1N6j5jzFtE9IeFCxdSw4YNQ9Lnyy+/LB9//DEAdBWRL0PSqboqIuoM4Kt+/fph7NixFIo+09LS8Pjjj4uIvM/Mvw9Fn+rqiIiMMZ6KFSvWTUlJsSpUqHDdfTIzHn/8cdvj8eSLSB19jsZ9RROKkydPRpcuXULS56ZNm/Diiy8CwDgReSkknaqrIqLaRORp1qyZNXfuXCsUfV64cAGjR4+28/Pzv7Ntu0Ww/QSdZRFRDSJaft9994Us6QaAzp07Y8yYMQCQRETO7rEqx4wxm6KjoyMWLVoUkqQbACpUqIDly5dbxpg6AJaHpFN1VUTUX0QGz5w5M2RJNwA89NBD1K5dOyGi9f6nwJVLiMgioo1t2rSRUCXdANCgQQMkJyeTiNznn01X7lpqjKn30ksvhSTpBgBjDBYtWmRVqlQp3BizOSSdqqsiooYAZicmJoYs6QaAPn36YPDgwQDwIhHVClnH6hcZYz6rUaOGNXv27JB9d1WuXBnLli2zADQjoplBxxZ0Q2M+i4mJwYgRI4Lt4qruueceNGjQgI0xn4W8c3UFEf2RmTslJSVZYWGhXVkyKioKEydONAASiahdSDtXVxBROBG9fccdd0jr1q1D3v+UKVNM+fLlKxLRWyHvXF1BRGvKlSsXNW3atJDvrdCyZUt069ZNiGgtEYWHun9ViIhaAxg7YcIEU6lSpZD2bYxBcnKyxcwdiGhESDtX/8UY8z/169fn6y0v+SXDhw9HzZo1YYzRh55dRERPiMjNs2fPDtmEYpGYmBiMHDmSAEzzX6Q5FlRERHQHM7eYNWuWa7NgSUlJlojUISJ9MMglxphlPXr0kFDOkhZ35513okmTJmKMWePKAAoA5pcvXz78iSeeCNksaXHly5fH5MmTjYj83r82vwoxIqoqIn+YNGmS44dhAzVx4kQKDw8vh8JVi5QLjDFrGjVqJN27d3el//j4ePTq1UuMMYtdGUCBiIaISL3k5GTXcptZs2ZZzHwLEfVwa4yyjArNHjRoENWq5c6NhYEDB6J27dpMRKuDaR9U4m2M+WvDhg3t63lo5FoqV66M1q1bizFmhmuDlGFE1JeZK40ZM8aVhK3I+PHjDTPf7F8bXIWYMebBvn37mlDfsSiuffv2iI6OtgHMc22Qsu0vUVFR3LFjR9cGMMagX79+ljHmT64NUoYRUS1mbj5+/HhXd4NOTEwkEYkmon5ujlNWGWOSWrZsKUVrc7shNjYWjRo1YmOMnk/dMZKIyl/vYh/XkpiYaIlIR//mZo44PkkQUSQz3zZixAjXaz4feughw8wNiSjO7bHKGmPM7ISEBA5VHeLVNG7cGFWrVtWkzQVENEBEKvzxj390faz+/ftbRDTE9YHKIGPMsH79+rl+Ph02bBhEJJKIBrk9Vhn0l8qVK9tOl350qkKFCmjatCkbY+a4OlAZRES1mblxYmKiqxdPADBq1CjDzG2DSdpUyYwxU9u0aQM3J6MAoEOHDoiMjGQAyU7bBvMGG2BZlrRp0yaggzdu3Ihly5Zd+XnVqlX461//GlDbuLg4RERE2Cjc7lqF1i29evUK+PXfuHEjHnjgATz66KMYPXo0Pvss8PL7Dh06WJZldQsmSFWiYbVq1WInSyRlZmZi3rx5mDhxIsaMGYMpU6YE1O53v/sdRCSCiNy7zVUGEVENZo4cNCiwXHjTpk2YMGECJkyYgAEDBlz58yuvvHLNtuXLl0ft2rVtAEOvM2z1E5ZldW/fvn3AF08rVqzA++//Z0NDEcGQIUNw+fLla7bt3bu3EZGbg4tUlWBoeHi43aBBg2seePLkSTzyyCP/9XcFBQUYNWoU0tPTr9m+ZcuWsCxLAIS+kLyME5H699xzT0B38o8dO4ZRo0aBmQEAtm1j+vTpAec3rVq1sowxfZzGGEzifVfVqlU50IO7d++OXbt2IScnBzt37sTu3bt/9oYtib+cRWuhQoiIopk5vFOnTgG3SU9Px7Bhw7BkyRJMmjQJKSkpAbe97bbbwMz6FHeIGWPaN2nSxNFM6YIFC9CpUycsXLgQy5cvx0MPPRRQu4iICJQrV84GoLOloXVfuXLl7EDvPPXp0weLFy/Gc889h5iYGCxevBiLFy/G6NGjA2qfkJAQZlnWbdcTsPo5Zo5t1y7wZ8g9Hg/i4+Ov/Hz69GlUqlQJgVxEd+zYESJSnoic7eairqVHbGxsQAfWqlULZ8+evZKwAcAnn3yCFi1a4Fq7WhapXr06A+gbTKDqlxFRYxExt956a0DHx8XFoV69eti1axeAwonhevXqoUePwFJO/4IG175S+wnHibcxpm2DBg0CnsOPiIjAnXfeiVWrVuHll1/G008/jfDwwB+sT0hIsCzLauk0zrKOiNoR0Wki+pKIkohoABHV9W+d2tuyLHayC5fH40HduoVl2rVq1UK5cuUCbtu2bVuIiEVEJe9brn6GiHYS0SEiSiGikUTUmoiK/ufXbdWqVcB92baN/fv3o0WL/yw/WvzL/1pq1KghAO4MuIECABDRGCI6RUQbiGgyEf2WiKr6f31HTEyM43Xu09PTA/6CL65169YQkdr+uMKIqBURPeh/fx0ioj2OOy3jiChaRMI6dOgQcJv09HQUn1lNS0sL+LMYFRUFy7KuJG1EVMd/fk8iov/nP+8Hv4NWGWVZVqtAJzKMMahRowbOnDkDAMjPz8eHH36IYcOGBTxew4YNw4wxbYOLtuwiou5EdIaIPieiGUTUr9jyjPeEh4fbTspMBg4ciA0bNmDr1q04dOgQRo0aFXDbTp06gZkd78rjuAiGiKKjo6MdtenZsyfGjx+PadOmwekDmVWqVIGIxBKRbsLinACoCaATgEsAwgHkAzhjWc5KSosSbxHBunXr4GQL5GIrNZwtzPtVEJqisOSKAUQSUSoRlatePfBrGcuycOutt2LChAm4/fbb0aNHDzRr1izg9lFRURaAPxDRH5wGrwAUJkq/BZAHoCIRnQNAFSpUcDwB4vF4gkq8q1evXjRbeghAY38sBsCVWlM91wYn0BVpcnJyYNs2qlateuXvPB4PAilxKBIWFia2bScT0UsAIlB4Xq8IwELhOWKnnmudcToZFRcXh+PHjyM2Nhbr169Hhw4d4GQVDf+Sk0318xYUG0A3AF1QmNtEEtElANlOV4Zq27YtXnvtNaxcuRJz5851VBterVo1AIXL+opIfqDtgqk+d/wmeeutt1C5cmXYtvPt7UUERMQAMh03LttqoPC1IgBe/99ZAM4DyIKD1zEjIwN5eXl4/vnnce7cOcTHxzu6si/mHAq/FFTgipbwswH4UPglmwPgNIAEp5uCzpgxAwcPHsTOnTsxY8YMTJw4EYGWHPnH8gLQnfOcKb4Moxf/SY7OAYgKpsP09HQE+pxNcf7XUFD4/okFEAngMgrfX0VX4xnBxFSGGQABXwH/0kWTx+NBt27dnI6bhcLXrJ7/Zy8KX88i+jo6ICJV4KAKoCjxbt68OdavX4/58+c7HQ8ACgBkO2qofno+Jf9/mQByg9kpu2nTpmjUqNF/XQy7yXHiLSLZFy4E/r37wQcfoKCgAJMnT8aaNWvQuXNnR+NlZmYCwCkRuclZpGUbEbUA8DaA/QC+BLAHwD4RySOiAT6f74NA+0pPT0eLFi0wa9Ys5OTk4OGHH8Z3332HW265JaD2xR4YaiIi5539S8o2IlqPwpmsLwDsBLBHRE4BgGVZeRkZGY5ucxERmjdvjubNmyMnJwcejyfgxPvixYsM4D0RCeqqq6wiomEApgDYCmAbCj+L34mIj4hW5eTkDIHDsr/09HTce++9jmM5c+YMiKiAmbv5Y4sF0A5AewBdAeSLyF2OOy7D/CtT5OTm5iKQWv0ff/wRxWuJmRmHDh3Cww8/HPCYPp+PAMwRkfeIKBJAKxS+jnf4/zxURL519i8p28LCwjyZmZkB30aqV68e9u3bh48++ghdu3Z1nLRlZWUBwGERCbxeUIGIOgJYAWA3gK9QeD7dLyL5RPSw1+v9v077PHbsGHr27Ok4lqJSIyez3UAQNd7MvCstLc0XyLH79u3Dli1b8Pjjj6Nly5bIy8vD0aNHHY13+PBh27btb5zGWdaJyAERaS4iQ0RkqYh8LSJ5/l//k5nNuXPnAuorLS0NRZvsREVFoWvXrlceRgjEzp07QUQ+TbqdE5F+ItJNRJ4TkfVFSbffD99+G/h36969e1FQUACg8KR/8OBBR7OmZ8+eJQC645pDIvKmiLQQkTEi8ob/s1l0Dv303Llzjs7DzIwTJ06gXr161z74J/bt2wcAx4vFdsr/vpohIt016XZORC4RUcGOHTsCOr5GjRo4cuQIvN7CG5HvvvsumjdvjkB3u8zKyoJt2wbAJ/7x80Rkh4gsE5Gh/vO+Jt0O2bb9TWpqasC35ePi4nDkyBFs2bIFga5KVFxaWpqPmQP/IlUAAH8uc4uIPCAiy0Vkd7HE9yOv12sVfbYC9cMPPwRVurdt2zYYY3KdtgtmVZMNWVlZ1ywQPnPmDJYsWYKnn376yixA//798Y9//MPRYKdPnwYA3To+hEQkzxhzedu2bQEdn56ejuK7W3bo0AF79gT+DNaePXtgjDnpOFBVImbemZqaGtBFMABs3boV48ePx6OPPoqZM2di+PDhaNq0aUBtc3Nz4fP5DID3r3mwcuIDn89ncnJyAm5w8uRJxMTEOHrAuUhqaqqPmXc6bqhKZIw5sXfv3oCObdu2LVq1aoVx48Zh7NixOHHihKPZ7u3bt8MYc1lELgUbr/pFW06dOnXto/zq1q0Lj8eDPn36ICrKecXY+fPnDYANjhuqqxKRY0RkO8lPMjIyULFiRURGRl774J/49ttvBcD3TtuR03IY/4oK+c888ww5LRtx6vDhw3jyyScBoKaIaL1aCBljttevX7/9kiVLSmPjDvvixYuvi0ii22OVJf4thz9du3ZtUCd+J1asWIEPP/zwkm3b7g5UBlmWdXHAgAFRTp6mD0Zubi4GDx4MAH1EZLOrg5UxRPRSVFRU4tq1a10/nz722GO2x+PZa9t24MuoqGvyr7p1dt68eQGXUQZrx44dSE5OFgCRTssUVMksy/p38+bNG8+ePdv1jZDuu+8+2+v1zhWRaU7aOQ5MRAqI6MtVq1Y5f1LSoVdffZWNMQc16Q49EZnk8XgsJ/X6wdi/fz8uXrxoAZjs6kBlkIh8Zoy58MYbb7g+1ieffGIz8+uuD1QGMfOrmzZtcv18unLlShhjsjXpdsXknJwc4y/lcU12djaOHj1qMfMkVwcqg0TknDFm/2uvveb6AgArV660iWirJt2hx8zPHjhwgJyWmzj1xRdfwOv1EoAkp22DuiIQkcePHz9uOa3XdiIjIwOHDh0yzPyMa4OUYSLylTHm7IsvvujqUkZ/+9vfbGPMNyISWEG5coSZl3766afs5knm888/R25urgGgn0V3TM3LyzNbtmxxbQCv14vNmzczM7/k2iBlmIhcIKI9L730kqsXUMuWLRNjTKaIfOHmOGUVM0/+7rvvjL/E1RVpaWn44YcfLBGZ6NogZZiIvA0g9+WXX3Z1nBUrVthE9HkwF0/BJt7/MsZsmz59ul1856ZQmjJlik1Eh0XkI1cGUGDmEV9//TXt37/flf7XrVuHY8eOGWa+35UBFAA86/P5cpKTk135IObm5mLx4sUM4GURuejGGGWdv1b3b0uXLuXcXMfP6QRkzpw57PP5LqFwdRXlAhEZfPz4cfPxxx+70v/BgwexdetWYuYHXRlAQUQ2GmMOTZ061ZULKGbG9OnTbWPMDhHZ7cYYChCRRzZt2gSPx+NK/6tXr0ZmZiZEZGgw7YOugWHmuy5evMhLliwJ+Yzp6tWrcerUKSpa8kq5Q0TWE9EnSUlJXGzJv5DIyMjAK6+8IgBeEJHUkHaurhARW0T6/etf/zJfffVVyPt/9tlnmZkzAYwLeeequEeY+ZwbX/jbt2/H7t27jYj0FxHXS1rKKhE5CmB2SkqKhHrG1Ov14vnnn2ci+qdORrmLmbufPn0abpTwLVu2TLKzs5mZe4W8c3WFiLxujNk1bdo02+cLeP2BgKSnp+Odd94REXkq2DLooBNvEbkoIn/YsmUL1qxZE2w3P7Nhwwa8/fbbADBORHQlDJeJyID8/Pzz48aNs0NVrnD+/HmMHz+eARwUES1PcJmIfAVgybx58+Sbb74JWb9z586Vw4cPCzP3CGZTAhU4ERFm7v79999j9uzZIbt7sW/fPsyZM0cALNPyBPf5H7I68PDDD3Ogy7Vei9frxdixY+3Lly9fEJH+IelUXZWInAYw5u9//ztCeffirbfewqZNm+Bf8lHvHrqMmXtlZ2fnTZgwIWSVGadOncLEiRPZX5+/KNh+HK9q8rMOiEYCSOnZsyceffRRMib4B0lff/11vPfeewDwjIjMva7AVMCIqLoxJrVSpUrR8+fPt5xse/tThw8fxpQpU7igoCCNmW/WGbbSY4xZCWD4+PHj6a67gl+O2ev14rnnnuMDBw4AQA9N2EoPEXUhov+55ZZbKCkpyTjd/ri4zZs3w/8Mx5vMPDx0UaqSEJFljDkUFhbWaPbs2ebmm28Ouq+MjAw8+eSTdlZW1iURaajPypQeIpoEYN69996L61lxiJmxdOlS2bx5MwCMFRF3i4/VFURUn4gO1qxZM3z+/PlWlSpVgu5r3759mDFjhjDzN8zc7nomo6478QYAIrqXiN6Ojo62kpKSrOJrPgfi5MmTmDp1qn327Fn4N5l49bqDUo4QUbQxZqeI3Dx06FAaOtRZ6RIzY8mSJbJlyxYyxnzKzL10lrT0EVEygCmNGjXimTNnWoFuylFk165deOGFF9jr9eaISA8RCXxBVBUSRNTGGPN5WFhY9J///Gdz++23O2qfk5OD6dOn26mpqQbAPBF52p1I1dUQERljNjFzrx49eshjjz3meFLqnXfewerVq4WIjjBzBxFxdwkq9TNENIKIUqpXr47k5GSrbt26jtqnpaVh+vTpdnZ2Nvtnut9zKVR1FURUyxizA8BNI0eOpIEDBzpq7/P5MG/ePNm+fTsR0Xsicv/15jYhSbyBwm1zjTGbmblzvXr17AceeMC61jrfe/fuxeuvv26npaVZxpj9zPxbXTrwxiKipwDMjYyMRO/eva3hw4cjIuLqu5JnZWXhtddeky+//FJs2y4QkT/qyeXGIqKmxpjPRCS2efPmMmrUKNO4ceOrHs/M2LhxI95991373LlzFhG9IyLD9G7FjUOF1ojI0GrVqvnuv//+sLvvvhslJW9Hjx5FSkoKHzhwgIjotP98erAUw1Y/QUSDiGiNZVnlu3TpQiNHjqSStha/fPky1qxZg02bNtl5eXkE4GkR+UvpRax+yn9H+FNmbh0fH28/+OCDVtu2bUtss337dqxatco+duyYZYz5mpl7a3nJjeWflHqmYsWK0q9fP2vw4MEo6Y5iRkYGUlJSZMeOHcLMeSIyKFRLsYYs8b7SIVEbIlooIncaY1ClShU7Pj7eqly5MhERsrOz4fF4fOfPnze2bRMR7RKRJ/11qupXgIgiASQbY0Yzc3RkZKSvXr16pk6dOiYsLAz5+fk4fvy4ffLkSeTn51vGmAxmXgDgL5qs/XoQ0SBjzCxmbhoWFsYxMTHcoEGDsMjISIgIzp07Jx6Px/avs26LyMcAJopI+o2OXRUiojgAC4loIAArOjrarl+/vhUTE0NEhLy8PKSlpdlnz54ln89njDH/ZuapevH760FEFoAnjDFPMXPN8PBwOzY2FnFxcVZ4eDh8Ph9OnjzJx48f59zc3DBjTA4zpwCYprtT/noQUWci+quI3G6MkWrVqnF8fHxYpUqVICLIzs6WtLQ0Oysry2JmENGXIvKE3jX89SCiygBeMMb8H2aOrFixoi8uLs6qVasWFeU26enp9unTp+H1ei1jzAlmngNgaSjv4Ic88b7SceEOl30B9DPGtCeiSv7xcph5D4CNANaJSJ4rAaiQIKLmAO4D0NWyrDgA5QB4bdtOBfA5gL+LiOcGhqiugYiqAvgdgN7GmJZEVAGAzcznRGQ7gI8AfK6lQb9eREQAugEYQESdjTHVAVgiksvM+wFsBvCBiJy/kXGqkhFRfRSeT3tYltUEQHkABbZt/wjgCwDviYg767uqkCCicAD3oDC3aUdEUSh8PvoiM+9C4Tbw60Wk4IYGqkpERO0A/I6I7jDG1AEQBiDftu1/A/gUhZ/FH10ZW79rlVJKKaWUcp/re9krpZRSSimlNPFWSimllFKqVGjirZRSSimlVCnQxFsppZRSSqlSoIm3UkoppZRSpUATb6WUUkoppUqBJt5KKaWUUkqVAk28lVJKKaWUKgWaeCullFJKKVUKNPFWSimllFKqFGjirZRSSimlVCn4/2LQRtzQfGLwAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 720x39.685 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig_2_5 = DAG(\n",
    "    ebunch=[\n",
    "        (\"R\", \"S\"),\n",
    "        (\"S\", \"T\"),\n",
    "        (\"U\", \"T\"),\n",
    "        (\"V\", \"U\"),\n",
    "        (\"V\", \"Y\"),\n",
    "        (\"X\", \"R\"),\n",
    "    ]\n",
    ")\n",
    "\n",
    "fig_2_5_node_pos = {\n",
    "    \"R\": (2, 0),\n",
    "    \"S\": (4, 0),\n",
    "    \"T\": (6, 0),\n",
    "    \"U\": (8, 0),\n",
    "    \"V\": (10, 0),\n",
    "    \"X\": (0, 0),\n",
    "    \"Y\": (12, 0),\n",
    "}\n",
    "\n",
    "fig_2_5.to_daft(node_pos=fig_2_5_node_pos).render();"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Figure 2.6. A directed graph in which P is a descendant of a collider."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAt4AAACnCAYAAAAmJlt/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAAoC0lEQVR4nO3deXhV1bk/8O+79kkChFFQUEDABh4oIINMAspQCiIFvVhboTi0CCLTT5DpCshUBS3KLa2AFUTFqXJlcAIplYsyKDMCIhAhkUmGGMCTAMnZ671/JOGCMux9cgb85ft5Hp4nSPZaL5K9z3evvfZaoqogIiIiIqLoMvEugIiIiIioKGDwJiIiIiKKAQZvIiIiIqIYYPAmIiIiIooBBm8iIiIiohhg8CYiIiIiigEGbyIiIiKiGGDwJiIiIiKKAQZvIiIiIqIYYPAmIiIiIoqBQCw6EZFyAK4HIAC+U9WMWPRLkSMiCQCqAUgG8AOAfaqq8a2K/BKR6wFUAJAL4KCq/hDnksgnEUkGUBVAAoDvVfVgnEsin0REkHc9LQMgC0C6qubGtyryS0TKA6gEQAEcVtXMOJdEPolIEvLOxeIATgD4NtrZJioj3iJSXkTmBAKBb40xuQC+B7ADwHYAx40xoUAgcFBE3hCRytGogQpH8vRxHGeL4zhZAHIA7AGwBcA3AKzjOCeNMWtFpFs8a6VLE5G6IvJ+IBA4JiIugEMAvgSwE8ApY8xZx3G+EZG/iEjx+FZLFyMiSSLyjOM4qcaYswCCyPv3+xLAARGx+f++H4pI/fhWS5ciIl2MMWscxzkJwALYh7zr6R4AOY7jZDmOs1VE+uUHc7rKiEhlEXk9EAgcNMaEABxHXq7ZAeB7Y0xuIBDYLyIv54dyusrkZ5shjuNsdxwnG8AZALuQdy6mIS/bnDDGfCYiHaJSQySDvYjUMsb801rbsFSpUm6zZs2cJk2aoEmTJihWrBgAICcnB5s3b8b69euxbt06NzMz0zHG7LLW/kFVN0asGApL/gV/poj8UUQS6tevj8aNG0vLli1RqVKlc9+XlpaGtWvXYuPGjXb37t0iIlnW2imq+lT8qqcCItLVGDPTWlu5UqVKbrNmzZxmzZqhfv36MCbvfjsYDOLzzz/Hxo0bsX79evfs2bNGRP5HVe/lU6n4E5FyIjJfVdsnJSXZpk2bOo0aNUKLFi1QunRpAIC1Fjt27MAXX3yBdevWuYcPH3aMMQettf1V9b04/xUIgIiMMsY8oaqlatasaZs0aWJatGiBGjVqnPueI0eOYPXq1di0aZNu27YNqpqrqq8A6Mcni/EnIo2MMW9Ya+uULVvWbd68udOkSRM0btwYiYmJAIAzZ85gw4YN2LhxI9atW+eeOnXKMcZstdbep6pfx/mvUOTlP7V/RUR+Z4wxDRs2lEaNGkmrVq1QoUIFAHnX07179+Lzzz/H+vXr7d69e40x5oS1dpyqTo9YLZE6p0VkCoARVatW1QEDBpi6det6Oi41NRUzZsywe/bsEQBzAPTlhSY+RKSFMWaJ4zil77vvPtO9e3cEAleejZSdnY3XX38dH374oQL41lrbVlXTol4w/UT+6Oh7qtqxSZMm2q9fP6lYsaKnY1etWoWXXnrJzczMVFUdpKqzolwuXYKI9BWRF8qWLSt9+/Z1Wrdu7em4I0eO4MUXX9T169eLiCxX1d+o6tkol0sXISLVjDErVLV6ly5d5P7770eJEiWueFwoFMLChQvx1ltv2VAo9IOq3qmqa2JQMv1IwUAUgL4pKSk6YMAAk5KS4unYHTt2YObMmW56eroB8JyqDo9mrXRpItJBRBYlJSUV79Wrl+nateu5AajLCQaDmDt3Lv71r3+piHxtrW2nqkcKXU9hM66IOMaYLwHU6du3r3Tp0iWsdj799FNMmzbNWmsPWWtrq2pWoQojX0RkKIDnbr75Zjt27FhT8ITCj8zMTIwZM8b99ttvBcB/cMQttkSkmohsK168eIknn3zS8Xrz+2Nz587FwoULFcCH1tquka2SrsQYs0hV77r77rvRu3fvsNrYuXMnJkyY4GZnZ59R1bqqmh7hMukyRKQLgPeqVKmiTz/9tFOuXDnfbeTk5GDSpEl2y5YtBsAIVf1LxAulSxKR4saYr40xVQYPHmzatWsXVjtLlizBrFmzFMAua+3NnMsfWyIyGcCoFi1a6MiRI8XLYOKPfffddxgzZox75MgRAGinqp8VqqbCBO/80L0rMTGxxsyZM03BcH24gsEgHnnkETcYDH5vra3KkZrYEJHhAJ59+OGHcddddxW6vRkzZuiSJUsA4G6G79jIn3uYWrVq1YS//vWvTjgXl/Pt3r0bI0aMUGvtv621v45QmXQFjuMsAdBpypQpUqdOnUK1FQqFMGTIEDc9PT1XVWup6v7IVEmXkx+63+/YsSMGDRpU6LnaH3zwAV588UUA+E9VnVLoAumKRCTJGJOenJxcYdasWU7B1K5wZWRk4NFHH3XPnj2731qboqpuhEqlyxCRZwCMGDp0KMK9cTrflClT7OrVqwGgbWHCd6GCt+M4nyQkJLSZPXu2KVu2bNjtnC87Oxu9e/d2s7Ozv3Jd9+aINEqXJCIdAXzct29fdO0aucHN/PCtAH7BaSfRJSJijDl+ww03lHnhhRccL4/QvEhNTcXjjz+u1trpqvpYRBqlSxKRqcaYoc8995x4fZx9JdZaDBo0yD1w4MAP1tprOI0vukSkGoC9nTp1koEDB0bsBcmPPvoIM2fOBIA7VXVJpNqli3McZ0vx4sXrvfzyy46X6UFenDx5Er1797a5ubmrXNdtE5FG6ZJEpCeAN4YNG4Y2bSL3v3vKlCl2zZo1VlUrhfsuVNif0CJyr7W23ZQpUyIWugGgRIkSmDZtmqOq9UTksYg1TD8hIgkisuDWW2/VSIZuAOjfv79UrVpVjTErItowXcwrjuOUnTZtWsRCNwCkpKQUjNgNFpHw5q2QJyLySwBDBwwYELHQDQDGGEybNs0JBAKlReS1iDVMF2WM+aRKlSoaydANAHfeeSdatWqlIvLf+S+JUZSIyEBVvXnatGkRC90AUKZMGUyZMsVYa28Xkfsi1jD9hIgki8grd9xxh0YydAPAqFGjzDXXXCPGmH+H20ZYn9Ii4ojIvA4dOmgkPyQKVKpUCb/73e9ERKaKSJmId0AAABF5KykpqfiIESOisnTVU0895QCoJiKjo9E+ASLSAMADw4YNC2te/pV06NABtWrVssaY5RFvnM4xxixLSUmxHTt2jHjbiYmJGD58uFHVXiJyS8Q7IACAiIxU1Rp//vOfnWi0P2zYMClWrFiSiMyPRvsEiEgpEfmve+65R66//vqIt5+SkoKOHTuqiLzKG6joMcZ8VKpUKfPoo49GLduo6s0iEtZLOOEOj410HCchEvPXLqVXr14Fb4BPjlYfRZmIJKnqf/Tv398Udj7wpZQrVw5dunQRY8zIqHRAEJG/V61a1W3ZsmXU+hg3bpxjra0kIu2j1kkRJiJtrLWVx48fH5XABgAtWrTAjTfe6IpIxJbEogsZY57o3LmzlC8fneWbA4EABg8e7Khqt/xNPyjyJhcrVgwPPvhg1DoYMGCABPI+dP8zap0UYSJS3lp7+6hRoyL6BPh8lStXRqtWrWCMCWv55LCqMsYMue2220y0/lIFOnfu7BhjoncGFG1PJiYmaiReOLicXr16QVVLiQhXx4gwEUlW1VYPPPBA1AIbAJQuXRopKSnWGMNVFaLAGDP1pptucsuUie7DvYceeshR1Vslb+dLiiAR6WytLR3NwAYArVu3RlJSkgUwIaodFVHGmIc6d+4c1eupMQZt2rQxxpj/F81+irBnS5Ys6davH929xB5++GGx1lYUkSZ+j/WdnEWkurW2wp/+9Ce/h/rWo0cPqGpxEYluOiyCjDF/uu2226J6gQHy5uzXrFnTisgT0e6rCHosKSnJtmjRIuodPfDAA8Za24g76kWWiIi19pb7778/6udi06ZNUaxYMQtgRLT7KmpEZHRKSoqN5JzgS2nTpo3jOM4fo95RESMira21yX/4wx+i3tdDDz0Ea+01IhL5ubpFnDHm3k6dOkX9elq+fHnccMMNLoBxfo8NZ8j63qSkJNfrC5VLlizBjBkzzv1+3rx5eO655zwdm5iYiNKlS7sA7vZfJl2OtfY6Py8dLFmyBPfffz8GDRqEPn364JNPPvF8bOPGjU3+y2MUWR2qVavmKwh///33ePbZZzFkyBA88sgjeOIJb/dDDRo0AAAB0Nx/mXQZjQFI48aNPX3zxx9/jMGDB2Pw4MHo1q3bua9feuklT8dXr17diEjbsKulixKRerfccovnz9O5c+diwYIF536vqrjvvvtw5syZKx7btm1buK5buLV76WLuKV26dKhgJ8rLOXz4MAYOHHjBf8vNzUXv3r2Rnn7lJfPLlCmDYsWKhQD8Ntxi6adExLHWlvr1r72tgLt//3707t0b1loAgOu6GDt2rOd8U79+fcdxHN/vzYQTvNvdcMMN3r+5XTusX78ewWAQ69atw4YNG37yA3s5NWrUcETE27ZtdE7+C7DtRaTSRf6sFgDj51FMeno6evbsib/97W8YPnw4Zs+e7fnYli1bwlpbmqOl/olIfRGpKyI/uYN3HKdenTp1fJ3Dzz//PG699VZMmzYNL774Ivr27evpOGMMSpYsGQJvgn0TkTIi0kZELrYYcPfk5OSQ12l7nTp1wvTp0zFu3DhUqFAB06dPx/Tp09GnTx9Px//yl78UY0y9i9To5P+ccQlXn/KfWpTx8+QpLS0N1atXP/f7I0eOoHTp0vDygnT+xljmYisNiUjF/Ot+1Ef8/n8jIq2qVavm6f9bxYoVcfz48XOBDQCWLl2KevXqoVq1ap76u+GGGwwAvjfjU/5qbO1F5GI3n78SEa1cubKntqpWrYoqVapg/fr1APIGhqtUqYL27b39szRt2hTW2uu81l7Ad/A2xqRUqVLF80ldrFgx3H777Zg3bx7+8Y9/YNSoUUhK8v5eSPXq1cUYU9VvnYS7ACwBkCYi34vIChEZn7+5wx0JCQmun5cq09LSUPDDXLFiRSQkeH8hu0aNGgVfersi0fk+B/AFgGwR2S4is0TkARGpa60tU7t2bc8Nua6Lbdu2oV69/8td53/4X8m1115rAHBZQf8mAlgG4LiIHBKRxSIyVETaAGhQoUIF3zek6enpnj/gz1e7du2Cm+C6+T9HM0VkG4As5P2cfeG7UboByFuxwqv09PTzr4vYt2+f53PRGIOEhAQXQCcRuVNExuVf378HkI6863537+UTABhjbqxevbqnc9EYg2uvvRZHjx4FAJw9exaLFi1Cz549Pfd34403Gsdxbgqv2iLtQeT9jB8UkeMiskxERotIJwBtihcv7mtzorvuugsfffQRVq9ejZ07d/raLbhRo0ZQVcfvoKLv5SxEJNFPcAbyliTr378/xowZA79L9CQlJUFVrxURbvwQviQAbfN/AQCMMfZS33wxBcFbVfHBBx/g/vvvD6eOfRz0LpS6+b8eyf+9+llC0HEcNGzYEIMHD0bz5s3Rvn17/PKX3mcAJSYmCoDf8FwslOsBdMv/BQBITEz0/f8zLS0trOBdvHhxIO+6v/0if5wEAPz3ja5gMAjXdXH+NvJpaWkXBPErcRxHcnNzLzdn8x1ea/0xxvgaFKxatSoOHDiASpUq4cMPP0SzZs1QsWJFz8fnZ5tf8HwLiyJv6mN5AL/O/wUAcBzHV/Bu3LgxXn75Zbz22muYMmUK/AxInjctqQTyBi888R28VTUnNzfX1zFvv/02ypQpA9f1v0tqKBRC/p0855b6czeApwCcQd4PxWEA6wF8CqCStXa414aOHTuG06dPY8KECcjIyED16tV93dmfpz0Ablvtz24AQQAJAHIBfIm8f8MvAPx3Tk6Or3N4/Pjx+Oqrr7Bu3TqMHz8eQ4YMwa233urp2NzcXAXwPwC8zU+hAs8C+A3yzsViAPYCWANgNYBeubm5twHwNTUgPT0djRo18l1ITk4OALgA7gHQFEAbADcDSASQA6A0gJq+Gy7aKgHwvH30xW6a0tLS0LZtW88dWmsVeT9XBwHcDqAZ8m7qspH3MzYWwIJLNkA/ISKf5+bmel4LsiB4161bFx9++CGmTp3qq7/c3FyIyLcAfuW31iLuIeS9IH4aednmIPKeDH8GoIG11veLx7Vr18YvfvGLC26GvQiFQgVfZvs5znfwttYePHz48C/gcZrKwoULkZubi5EjR+KNN96A3/WGv/32W6jqYVVN9VtrUSYiMwGcAvAVgC2qGjzvz27Jzc0dZa2Fl7ml6enpqFevHp566ikEg0EMGDAAX3/9NerUqeOplu+++67gy1Wq6u+urYgTkW7IC2UbARw8f8vvQCAQ/Oabb8r6OadEBHXr1kXdunURDAaRlpbmOXhnZGRYANt4LvojIkMBfIi8f8Md558DIlLv+++/9/0OS3p6Ou6++27ftezZswcikmWtXQxgcX4NAqAygFsAKP99/RGRfQBw8OBBeJlbevDgQVSq9H+v3lhrsXPnTgwYMMBTf9Za5OTkOAAWqurnAKbn15EMoCHynoq9ef41n67McZyDBw4cuAZ5I6lXVKVKFWzduhXvvfce2rRp4zu0HTp0SK21aTzf/BGR55A3gLcVwFZVPX3en3U/ffq0701t9u/fjw4dOviuZceOHQBgz/9c9iKclytXHjhwwNM0ha1bt2L58uV47LHHUL9+fZw+fRp79+711dnevXtD1to1YdRZpKlqlqr+Q1VXXeQCvAmA7tq1y1Nb+/btw0035U1FK1myJNq0aXPuZQQv1qxZA2NMFkO3f6r6gaouVtUDPz65rbU7t2/f7vmE37RpEwqeVp04cQJfffWVr1HTU6dOOQDe93wAAQBUNU1V56jqloucA4t++OEHX6Pd1locOnQIVapU8V3Ljh07VFW//lF9mv/ztVhV3/PdaBGnqq4xJrh69WpP33/ttddiz549BU8fMH/+fNStWxelS1/s3dufyv8MVfxoPn7+NX91/nWfodsna+3avXv3en4sX7VqVezZswfLly9H9+7+p9Tv37/fVdWVvg8s4lQ1U1VfVNXPzw/d+T601prjx4/7avPbb78Na+reunXr4DjO936PCyd4v5udnR240rJHR48exd/+9jeMGjWqYAdKdO3aFYsXL/bckbUWmZmZDoCPwqiTLkFV1RiT+dln3p6OpqennwveANCsWTNs3LjRc3+bNm1SALyrjzBV/TQtLc3zB8Xq1avRv39/DBo0CBMnTkSvXr3g9eXM1NRUqKoAWBFuvXRRn6qq7N692/MBhw8fRoUKFXy94Fxg3759rqqu8n0gXZaq7t68ebOnm+DGjRvj5ptvxqOPPop+/frh0KFDnke7AWDlypUwxpz0O8pGV/T+iRMnnPNXKrmcypUrIy0tDZ06dULJkiV9dZSTk4OsrKwAgIVh1EmXoKpnjTGnV6zw/jF17NgxJCcnF7z/4su2bdus67pf+j1Owjl3Hcf5oWvXriUffvhh38f6sXjxYsyZMydXVZN4kYksEXkxOTm599tvvx3VZadCoRDuuecetdb2UdU50eyrqBGRawEcfeaZZ3y9JBmO0aNH2+3bt+91XZfzfyPMcZxddevWTXn66aejuhXw7t278fjjjwNAJVU9Es2+ihoRedAYM3f+/PniZR3owujRo4cbDAZfUdXofgAXMZLnzEMPPZQYzgi2H3PnzsWiRYuyXNf1l9jpikRkYfny5bu+8sorUc022dnZ+P3vfw8AnVV1qZ9jw7rQW2tnL1u2zP+bkj69++67IVVdyNAdFSOysrLMli1botrJggULoKo5AF6OakdFkKoeM8Zse/nll32tUONXTk4Otm3bJtZa3zt00ZVZa8ds375dCqYeRMvs2bOtMWYHQ3fkqeqrqnr23XffjWo/27ZtQzAYNAA8vxxP3uRPuXp34cKFUc82S5cuda21c6PdTxE1JCMjw0lLS4tqJ6+++iqMMaf8hm4gzOANYMzp06flrbfeCvPwK1u2bFnBNJOhUeukCFPVk8aYNc8//3zULjJnzpzB22+/bVX1Zd48RYe1duiuXbuM1/n64Zg6daoVkaCqvhm1ToowVZ0vIqeeeeaZqN1ApaamYufOncZaOyxafRR1qjr7nXfesV52nwzX1KlTXWPMF6qaGbVOirbHT5w4YZYu9Z2lPHvnnXeQnZ0tALxtG0y+qGqaMearyZMnRy3bnDx5EkuXLlVrrbdt2H8krOCtqlkARr711lt65EjkB0+CwSBmzJihAP6qqgcj3gEBAKy1XU6cOKFet5v2a+LEidZ13VMAvE9gJF9UdbmIrBg3bpzrdW6iH1u2bMHatWuNtZZbG0eRtbb7unXrzObNm6PRNsaOHesaYz4NZ3SGPBtsrT0xYcKEqHzgz5kzB5mZmWqtvTMa7ROgqocBPD9r1iw9depUxNs/duwYXn/9dQXwhKr+EPEOCABgrb3j0KFDMn/+/Ki0P3bsWBfAIVWdGM7xYc8pVNWpIrJr2LBhbiQfkVprMXToUDf/cShHu6NIVU+q6iPvvfeer1VKvFi8eHHB9ITOHO2OLlXtkp2dnTt+/PiIJu+MjAxMnDjRish7qroskm3ThVT1ExFZMGnSJJuRkRHRtidNmmSzsrJC1to7ItowXUBV1Vrbafv27WbBgsguob1p0yYsWrRIVXUAR7ujS1WHqeqhoUOHRnQwIycnB48//rgrIqmq+kzEGqafUNX9AEbPmzdPd+7cGdG2582bh3379hlrrbd95S+iUC/zWGubnTp16odHHnnEjcTjtVAohEGDBrnfffddjrW2MQNb9KnqyyIye9KkSRqp8P3+++9j9uzZAPB4/jqzFEWqelpVb9u8eTOeeuqpiHxSHD9+HI8++qh1XfcbVb07Em3S5anqb13X3dOvXz977NixiLQ5efJku2HDBqjq7RdZeosiTFU3AHgs/+W5iLS5adMmjB8/XkXkVVX9R0Qapcuy1t5y9OjRswMHDnTP2yQlbDk5OejXr5974sSJLGtt0wiUSFegqlMAvD9q1KiIhe833ngD77zzjgL4o6p6X4rqR8Ja1eSCBkTKG2N2JSUllRs7dqypX79+WO3s3bsXY8aMcYPB4BlVrauq6YUqjHwxxrykqg9369YNffr0CauNUCiEZ555xn7++ecGwHBV9beVFxWKiLQUkZXXX3+9PP3000758p43YbvAihUrMH36dGut3W2tvZnrr8eOiDjGmC9FpPbgwYNN+/bhDapkZGRg9OjR7sGDBwGgLZcQjC0RGQLguWbNmunIkSNNuCudzJkzB4sWLVIReTWcHfkofCJSVUS+KlmyZPGJEyc6KSkpYbWzY8cOTJgwwZ49e/aktbamqkb2kRZdljFmIYC7evToIT169AirjTNnzmDSpEn2yy+/FAAPq2qhFosodPAGABFJEJHFqtq5RYsWOnDgQClTpoynY7OzszFr1ixdsWKFGGPWWGvv4Nyn+BCRP4nIi2XLlpWhQ4c6DRs29HzsihUrMGPGDHv27NnTqnq3qi6PXqV0KSJS1RizQlVv6t69u/Ts2RNeP/QPHjyIadOm2V27dgmAWQAG8KlT7OXvJPl3AI/WqlVLhw4darzsiAjkjaz985//xPz581VE9llr23MQIz5EpL2ILE5MTCzRv39/XzdRW7duxfPPP+9mZmaqqvZX1ei8iEOXJSLJxpgl1trb2rVrp3379hWva3afPHkSM2bM0DVr1oiIfKyqd6nq2SiXTBeRfyM89brrrtNhw4Y5XnfdttZi6dKlmD17tnVd95S1tpOqrit0PZH8XBWRLsaYV6y1FapXr+527NjRadmyJX488nby5EmsWbMGy5cvt7t37zbGmFPW2v6q+kbEiqGwiEgZY8wH1tpWycnJtn379k7r1q1Ru3btC7aXD4VC+PLLL7Fy5UqsWrXKzcnJERFZqKo9OEIafyIyWkTGAEhq0KCBtm/f3jRv3vzcZlYFDh48iNWrV2PZsmXukSNHHGNMurX2LlXdGp/KqYCI1DfGLLbW1rjuuutCnTp1CrRq1eon25JnZ2fjiy++wIoVK+zWrVtFVc+q6tOqOilOpVM+EXFE5E1V/W1iYqK2atXKuf3229GwYUMEAoFz32etxa5du7Bq1Sp88sknbjAYNMaYtdba33BOd/yJyH3GmFnW2jI1a9Z0O3To4LRs2RJly5a94PsyMjKwZs0aLFu2zE1LS3OMMRnW2j+qKnf8jTMRqWiM+cha27h06dKhDh06BG677TbcdNNNF2SbnJwcbN68GStXrtS1a9eq67qqqvMA/ClSA1ERDd7nGhVpLSLPikhDa21xY4wGAgELAK7riuu6xhhzFsBOa+1/8k37q4+IlAHwjDHmXmttOQCSkJDgGmPguq6EQiEDwBpjjlpr5wCYwMB99cnf2OMJVa2hqgmO49hAIKCqitzcXKOqYozJUtXPVHWI/mg7cYo/EaklIv8lIrdZa0uKiCYkJFgA4rouXNc1IpIrImnW2qdU9dV410wXEpEEAGOMMX1UtaKqmkAg4DqOA2stcnNzHQBqjMm01r4LYCQD99VHRDoYY6YAqGutLeY4jnUcRwEgFAqJtdYYY86o6lZVHa6q3raHppjJ33juOWNMV2ttGVyYbRAKhRwRcUXkO2vt3wH8RVUjulJRVIL3BR2IFAfwawCVkPcy51EAH+cvSUg/EyLSCEAjAMUBZAFYU5iXCyj28i84vwJQHkAOgHQA/470RYWiR0QcAO0AjAVwO4A/APiXqkbmbUyKCRFJAdAaQDKA0wC2AtjEqV0/HyKSDKATgOsAWADfIe9c5EvMPxP50/qaA6gLoBiAIICVqpoW1X55nhMR/byIyD8A9FFViXctRETkXaGWEyQiIiIiIm8YvImIiIiIYoDBm4iIiIgoBhi8iYiIiIhigMGbiIiIiCgGGLyJiIiIiGKAwZuIiIiIKAYYvImIiIiIYoDBm4iIiIgoBhi8iYiIiIhigMGbiIiIiCgGGLyJiIiIiGKAwZuIiIiIKAYYvImIiIiIYoDBm4iIiIgoBhi8iYiIiIhigMGbiIiIiCgGGLyJiIiIiGKAwZuIiIiIKAYYvImIiIiIYoDBm4iIiIgoBhi8iYiIiIhigMGbiIiIiCgGGLyJiIiIiGKAwZuIiIiIKAYYvImIiIiIYoDBm4iIiIgoBhi8iYiIiIhigMGbiIiIiCgGGLyJiIiIiGKAwZuIiIiIKAYYvImIiIiIYoDBm4iIiIgoBhi8iYiIiIhigMGbiIiIiCgGGLyJiIiIiGKAwZuIiIiIKAYYvImIiIiIYoDBm4iIiIgoBhi8iYiIiIhigMGbiIiIiCgGGLyJiIiIiGKAwZuIiIiIKAYYvImIiIiIYoDBm4iIiIgoBhi8iYiIiIhigMGbiIiIiCgGGLyJiIiIiGKAwZuIiIiIKAYYvImIiIiIYoDBm4iIiIgoBhi8iYiIiIhigMGbiIiIiCgGGLyJiIiIiGKAwZuIiIiIKAYYvImIiIiIYoDBm4iIiIgoBhi8iYiIiIhigMGbiIiIiCgGGLyJiIiIiGKAwZuIiIiIKAYYvImIiIiIYoDBm4iIiIgoBhi8iYiIiIhigMGbiIiIiCgGGLyJiIiIiGKAwZuIiIiIKAYYvImIiIiIYoDBm4iIiIgoBhi8iYiIiIhigMGbiIiIiCgGGLyJiIiIiGKAwZuIiIiIKAYYvImIiIiIYoDBm4iIiIgoBhi8iYh+RkTkcQB98r9eJSJJcS6JiIg8YvAmIvp5qXXe180AhOJVCBER+cPgTUT087IGQDD/692q6sazGCIi8o7Bm4jo52XDeV9/FrcqiIjINwZvIqKfl68BJALIBrA2zrUQEZEPDN5ERD8j+VNL9gAoAWBjnMshIiIfGLyJiH5+VgPIQd7oNxER/UyIqsa7BiIiugwRSQbQE8CvjTFNRKQCgOIAjltrd6vqSgBvq+pXcS2UiIgui8GbiOgqJSItROR5VW0RCAT0uuuuszVr1gxUrFgRIoKsrCykpqbq/v373aysrIAx5qi1diqAqcqLOxHRVYfBm4joKiMixY0xS6y1bapVq+Y++OCDTtOmTS97TGZmJubMmaOrVq1Sa22Wqt6lqitiVDIREXnA4E1EdBURkQ4i8kFycnJgwoQJTq1ata580HlycnIwZcoUu379eiMir1tr749SqURE5BODNxHRVUJE7gawoHXr1hg+fLgYE/777xs2bMCkSZMUwArXdX8VqRqJiCh8DN5ERFcBEWkJYFWXLl3Qr18/iUSb+/btw2OPPaaqusBa+9tItElEROHjcoJERHEmIo6ILGnUqJFGKnQDQI0aNfDnP/9ZVPWe/NF0IiKKIwZvIqI4E5E3EhISSo4ZMybi1+T69eujbdu2KiJviUhSpNsnIiLvGLyJiOJIRMqp6u+GDx9uEhMTo9LHkCFDJCkpKQHAtKh0QEREnjB4ExHF119KlixpW7RoEbUOjDHo0qWLY4x5MGqdEBHRFTF4ExHFkTGmZ5cuXZxo99OzZ0+oanER6R7tvoiI6OIC8S6AiKioEpFrARTv3t17Fl6yZAnefPNNlCtXDtnZ2bjvvvvQoUOHKx6XmJiI66+/3j106FAPAAvCr5qIiMLF4E1EFD/3JCQkuCVKlPA84p2eno6ePXuic+fOSE1NxZNPPukpeANArVq1AkeOHGkSdrVERFQonGpCRBQ/t1WoUMHXZgppaWmoUqUKAKBChQqw1no+tkGDBlDV6/2VSEREkcLgTUQUP2VLlCjh6zpcELxVFe+//z6aNm3q+djy5ctDVfmkk4goTngBJiKKH1+j3ceOHcOZM2cwbtw4OI6DWrVqoV+/ft47407FRERxxeBNRBQ/mcFg0MLj08f09HQ0aNAAEyZMCKuzo0ePQkRCYR1MRESFxqkmRETx8++MjAzP1+F9+/ahRo0aYXe2detWADgQdgNERFQoDN5ERPGzMBQKmWAw6Omb09PTUb169bA7S01NDVlr14XdABERFYpwzh8RUfw4jvNDt27dSvbu3Tuq/WRnZ+P3v/89AHRS1WVR7YyIiC6KI95ERHFkrZ3z8ccfu9Hu57XXXoMx5hRDNxFR/DB4ExHF1+jTp0+b5cuXR62DnJwcLFu2zFprZ0WtEyIiuiIGbyKiOFLVLAAzX3jhBZudnR2VPiZPnmxDoVAWgCei0gEREXnC4E1EFH8DrbUZo0ePjviUk7Vr12LDhg1GVbuqatSntBAR0aUxeBMRxZmqqrW23TfffIOnn37a+x7wV7B161ZMnjxZAcxQ1ZWRapeIiMLDVU2IiK4SItJaRP6nTp06MmnSJJOYmBh2W8uWLcPf//53BfCmtbZX5KokIqJwMXgTEV1FRKSRMWZFIBAoNWLECNO8eXNfxweDQYwdO9ZNTU01AJ5V1VHRqZSIiPxi8CYiuspInjdUtcc111wTuvfeewN33nknjLn07MC9e/di9uzZdvv27SIiR6y1v1LVr2JYNhERXQGDNxHRVUpEqgKYJiJ3AXBKlSrlVqtWzalQoYKICE6fPo19+/a5x48fl1AoZIwxu6y1o1X13XjXTkREP8XgTUR0lRMRAdAWQDcRaWmMKQ/AUdVsa+02AMsALFTVzHjWSUREl8fgTUREREQUA1xOkIiIiIgoBhi8iYiIiIhigMGbiIiIiCgGGLyJiIiIiGKAwZuIiIiIKAYYvImIiIiIYoDBm4iIiIgoBhi8iYiIiIhigMGbiIiIiCgGGLyJiIiIiGKAwZuIiIiIKAb+F/8sMbJtwc/XAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 720x153.071 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig_2_6 = fig_2_5.copy()\n",
    "fig_2_6.add_edges_from(ebunch=[(\"T\", \"P\")])\n",
    "\n",
    "fig_2_6_node_pos = fig_2_5_node_pos.copy()\n",
    "fig_2_6_node_pos.update({\"P\": (6, -2)})\n",
    "\n",
    "fig_2_6.to_daft(node_pos=fig_2_6_node_pos).render();"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Study question 2.3.1\n",
    "\n",
    "(a) List all pairs of variables in Figure 2.5 that are independent conditional on the set Z = {R, V}."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-04-14T13:04:02.248851Z",
     "start_time": "2023-04-14T13:04:02.245074Z"
    }
   },
   "outputs": [],
   "source": [
    "def pairs(v):\n",
    "    return sorted([(x, y) for x in v for y in v if x < y])\n",
    "\n",
    "\n",
    "def print_independent(graph, x, y, z=set()):\n",
    "    if not graph.is_dconnected(start=x, end=y, observed=z):\n",
    "        print(IndependenceAssertion(x, y, z))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(S ⟂ U | R, V)\n",
      "(S ⟂ X | R, V)\n",
      "(S ⟂ Y | R, V)\n",
      "(T ⟂ X | R, V)\n",
      "(T ⟂ Y | R, V)\n",
      "(U ⟂ X | R, V)\n",
      "(U ⟂ Y | R, V)\n",
      "(X ⟂ Y | R, V)\n"
     ]
    }
   ],
   "source": [
    "z = {\"R\", \"V\"}\n",
    "\n",
    "for x, y in pairs(set(fig_2_5.nodes()).difference(z)):\n",
    "    print_independent(graph=fig_2_5, x=x, y=y, z=z)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(b) For each pair of nonadjacent variables in Figure 2.5, give a set of variables that, when conditioned on, renders that pair independent.\n",
    "\n",
    "The function `print_implied_conditional_independencies()` below returns a minimal d-separator for each pair of nonadjacent, observed variables in a graph. This function differs from the R function `daggity::impliedConditionalIndependencies()` in that the later may return multiple independecies for each pair."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-04-14T17:04:56.715587Z",
     "start_time": "2023-04-14T17:04:56.715031Z"
    }
   },
   "outputs": [],
   "source": [
    "def print_implied_conditional_independencies(graph):\n",
    "    independencies = Independencies()\n",
    "    for x, y in pairs(set(graph.nodes()).difference(graph.latents)):\n",
    "        if (not graph.has_edge(x, y)) and (not graph.has_edge(y, x)):\n",
    "            z = graph.minimal_dseparator(start=x, end=y)\n",
    "            if z is not None:\n",
    "                independencies.add_assertions([x, y, z])\n",
    "    return independencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-04-14T17:04:59.061106Z",
     "start_time": "2023-04-14T17:04:59.040599Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(R ⟂ T | S)\n",
       "(R ⟂ U)\n",
       "(R ⟂ V)\n",
       "(R ⟂ Y)\n",
       "(S ⟂ U)\n",
       "(S ⟂ V)\n",
       "(S ⟂ X | R)\n",
       "(S ⟂ Y)\n",
       "(T ⟂ V | U)\n",
       "(T ⟂ X | S)\n",
       "(T ⟂ Y | V)\n",
       "(U ⟂ X)\n",
       "(U ⟂ Y | V)\n",
       "(V ⟂ X)\n",
       "(X ⟂ Y)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print_implied_conditional_independencies(fig_2_5)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(c) List all pairs of variables in Figure 2.6 that are independent conditional on the set Z = {R, P}."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(S ⟂ X | P, R)\n",
      "(T ⟂ X | P, R)\n",
      "(U ⟂ X | P, R)\n",
      "(V ⟂ X | P, R)\n",
      "(X ⟂ Y | P, R)\n"
     ]
    }
   ],
   "source": [
    "z = {\"P\", \"R\"}\n",
    "\n",
    "for x, y in pairs(set(fig_2_6.nodes()).difference(z)):\n",
    "    print_independent(graph=fig_2_6, x=x, y=y, z=z)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(d) For each pair of nonadjacent variables in Figure 2.6, give a set of variables that, when conditioned on, renders that pair independent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(P ⟂ R | T)\n",
       "(P ⟂ S | T)\n",
       "(P ⟂ U | T)\n",
       "(P ⟂ V | T)\n",
       "(P ⟂ X | T)\n",
       "(P ⟂ Y | V)\n",
       "(R ⟂ T | S)\n",
       "(R ⟂ U)\n",
       "(R ⟂ V)\n",
       "(R ⟂ Y)\n",
       "(S ⟂ U)\n",
       "(S ⟂ V)\n",
       "(S ⟂ X | R)\n",
       "(S ⟂ Y)\n",
       "(T ⟂ V | U)\n",
       "(T ⟂ X | S)\n",
       "(T ⟂ Y | V)\n",
       "(U ⟂ X)\n",
       "(U ⟂ Y | V)\n",
       "(V ⟂ X)\n",
       "(X ⟂ Y)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print_implied_conditional_independencies(fig_2_6)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(e) Suppose we generate data by the model described in Figure 2.5, and we fit them with the linear equation\n",
    "\n",
    "Y = a + b \\* X + c \\* Z .\n",
    "\n",
    "Which of the variables in the model may be chosen for Z so as to guarantee that the slope b would be equal to zero? \\[Hint: Recall, a non-zero slope implies that X and Y are dependent given Z.\\]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "set()\n",
      "R\n",
      "S\n",
      "U\n",
      "V\n"
     ]
    }
   ],
   "source": [
    "nodes = [set(), \"R\", \"S\", \"T\", \"U\", \"V\"]\n",
    "\n",
    "for z in nodes:\n",
    "    if not fig_2_5.is_dconnected(start=\"X\", end=\"Y\", observed=z):\n",
    "        print(z)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "X and Y are marginally independent, so the Z term could be omitted. X and Y are conditionally independent when Z is any member of the set {R, S, U, V}."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(f) Suppose we generate data by the model described in Figure 2.6, and we fit them with the linear equation\n",
    "\n",
    "Y = a + b \\* X + c \\* R + d \\* S + e \\* T + f \\* P .\n",
    "\n",
    "Which of the coefficients would be zero?\n",
    "\n",
    "The coefficients for the following variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P\n",
      "R\n",
      "X\n"
     ]
    }
   ],
   "source": [
    "nodes = [\"P\", \"R\", \"S\", \"T\", \"X\"]\n",
    "\n",
    "for x in nodes:\n",
    "    if not fig_2_6.is_dconnected(start=x, end=\"Y\", observed=set(nodes).difference(x)):\n",
    "        print(x)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.4 d-separation\n",
    "\n",
    "Figure 2.9. A causal graph used in study question 2.4.1. All error terms (not shown) are assumed to be mutually independent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-04-14T12:01:57.043428Z",
     "start_time": "2023-04-14T12:01:55.250166Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARgAAAEYCAYAAACHjumMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAAk9ElEQVR4nO3de3gU5d3/8c/3nk1IQkJMIIACVY7iAblEEEHwwANcRkUqFNGWCgpqieKF1scAAtJEqKciT2lF6qFaUAQrPpyEUioPQlFADoIWjSHkJyqHBEOAHAg79/f3RzaAQJKd3ZnN7sz39Y9J2LnvO5fZ987M7s4SM0MIIZygGnoBQgj3ksAIIRwjgRFCOEYCI4RwjARGCOEYCYwQwjESGCGEYyQwQgjHSGCEEI6RwAghHONzcnAiagpgKICbAKQDIAClAD4G8D4z73dyfiFENSLqAOAXAHoCaAzABHAQwD8B/C8zlzkyr93vRSKiFADPKqVGaq0bx8XFmc2aNUNKSooCgLKyMl1UVISqqiqDiCqZeSGA/2bmIlsXIoTHBaIyk4gymdmXkJBgZmRkICkpSWmt+ciRI/rHH39UpmkqpVSJ1no2gGeY+aRta7ArMEREAOYCGJOUlKQzMzONu+66C0lJSee9fWVlJT744AMsX77cPHr0qEFE7zPzcGY2bVmQEB5FRClKqdVa6+syMjLMYcOGGQMGDIDPd/4DlpKSEsybNw//93//p/1+PzPzZGZ+1pa12BEYIuqilFqrlEp79NFH1c0332xp+y1btuC5557TVVVVZcx8GzOvD3tRQngQEf2SiN5s0qSJysnJMdq1a2dp+0WLFmH+/PlMRHu01n2Y+WBY6wk3MER0M4B/XnrppZSTk6Nq22OpT1VVFWbMmKG3bt1KAIYx8/thLUwIjyGiSQCmZ2Zm8m9+8xtSKrTncIqLi/Hkk0+axcXFVcx8FTPnh7ymcAJDRNcA2NynTx/Kzs6mkAc6w9y5c3n58uUAcCszr7JjTCHcjogeBfA/Y8eOxa233hr2eFprjB8/3iwsLKxi5rah7smEHBgiMoio5Morr2w8Y8YMW5/ufvHFF/njjz8+yczpTp3dFsItAidz8+69914aNmyYbeNqrTFmzBjz8OHDBaZpdgpljJDDQEQL4uPjG0+bNs3219I8/vjjlJycbCilVto9thBuo5T66JJLLtF2xiUwLp5//nmDmTsQ0YSQxghlIyK6hJmHZWdnq/j4+FCGqJNSCjk5OYbWui8R9bZ9AiFcInBo1PqZZ54xnBi/WbNmGDFiBBHRM0QUZ3X7UPc+XmratKm/R48eIW5evw4dOuDiiy82iehFxyYRIsYppSb26tULqampjs1x1113wefzEYCnrG5rOTBU7fa77rrL0VcBA8CoUaMMZr6OiBo7PZcQsYaIrtZatxwzZowtT7DU5cYbb1RKqXFWtwslEv8FwLjlllvqveGsWbNQUFBw6vuSkhL4/X4sWLAgqIm6d+8On8/Hfr//lwBeDWGtQrjZ2LS0NH+zZs3qvB+Hez8EgHvvvRdr1qxJJ6LGVp54CSUwg1JSUkylVL3bjh8//tTXBw4cwIQJE/DII49Ymqx58+b6hx9+GAAJjBA/oZTq1b59+4jcD9PS0mAYhjZNczCAd4Jeo6VZABBRr0suucTSCaXS0lJMmzYNw4cPR8+ePS3N16lTJ59SqruljYTwhnZdunQJ+sbh3A8BID09XQMYaGUby4FRSqU3bdo06GO+yspK5Obmok+fPsjMzLQ6HZo1a1bzBkohxBmYOb5Vq1ZB3Tbc+yEApKamKgDNrWwTyiGSEexLkE3TxAsvvIDWrVtjxIgRp35+4MABLFq0CGVlZZg4cWLdkxkGADjyFJwQMY4C94861XY//OSTT/DZZ5+hvLwcAwYMQLdu3eocRylFsNgMy3swzFxeXl4e1G3nzJkDv9+PceN+evK5ZcuWePTRR4Ma4/jx4wBwwuIyhfAC88iRI/XeqLb7Ya9evTBu3Dg8/PDD2LBhQ73jVFRUaADHrCzQcmC01l/s3bu33ksqLFiwAHv27MGECRMQTGVrk5+fr03T/DrkAYRwKaXUoS+++KLO2wRzP1y4cGFQ718qKioCgM1W1hjKIdLq4uLiX9R1g4MHD+Kdd95B8+bNkZ2dfernrVq1+sn3wdi3b58GsC6EdQrhaqZpbsvLy7sQtZxCqO9+yMx46623cM0116BDhw51zuX3+1FZWWkAsHSVg1ACs9jv97+Wn59f66JatGiBZcuW1TrA0aNHMW/ePBQUFOC9995Dbe+hOHz4MMrLy30Agn/CXgjv+OD7778f5Pf7z3sxqfruh8uWLcOOHTtQVlaG/fv313nyd/Xq1SAiv9ba0qUbQno3tWEYX11++eUdf//73zt60fDnn3+e//3vfx8yTbOlk/MIEYsCr6qvHDlyZPzQoUMdneu+++4zi4uLVzDzYCvbhRQIrfVTX375JR09ejSUzYNSVVWFjRs3stb6eccmESKGcbW/L1682NHLzBYUFKC4uNgA8JjVbUMKDDO/T0QHpk6d6tgvNmPGDK21Pg7gJafmEMIFHj127BjefPNNRwbXWmPq1KkmEX3KzAX1b/FTIR/iaK3779mzR61aZf9F5zZv3oytW7cqZr6d7f7YAyFchJkPM/O4xYsX83fffWf7+K+88gofPXpUM7OlV/DWCDkwzPwfADNefvll3r59e6jDnCM/Px/Tp09nInpTLv4tRP2YeQ4RbRk/frwuKSmxbdwVK1Zg5cqVxMyjmNnS619qhHWSlpknA5j/9NNPczAv1KnP559/jt/+9rfMzMTMS8MeUAiP0FqvqqqqUg888IC2Y09m0aJFeOWVVwBgIjMH/ebGs9n1sSV/BDDu6quv1pMnT7Z8lTu/34+ZM2fy+vXriYgWMHMzAAMADGHmD8JeoBAuRkTTADwNIEcpdQsz9xgyZAiNGjXK8lilpaWYPHmyWVhYqACMZ+Y/hrU2Gz94rQ8RrTAMI7l///5q5MiRSE5OrnObyspKvP3221i5cqWuqqo6wczDmHlFYLzVkMgIUacz4vI7Zp4W+NlYIpqdnJyMoUOHGoMHD671Q9dqFBUV4bXXXuNPP/0UAL7TWvcL5+NKTq3PznOogU93fEYplaW1viAjI8PfoUMH31VXXYWMjAwQEQ4fPoydO3ciPz/ff/DgQYOIyrXWfwXwBDOfOGs8iYwQtThfXM74t6YA/oeIhhFR3EUXXaQvvfRS48orr0RqaipM08QPP/yAnTt3ckFBgVlSUuJTSu3XWs9g5j/ZtkannqQhon4ARimlegL4GTM3Cvz8BIDvtdabASxg5tpfagiJjBDnU1dczrodARgJ4OeGYXRj5hbMHA9AE1ElgD1a608AvMzMn9u+zkg9C0xEhwBkMLPl64dKZIQ4Ldi41LE9A3iLmUfZu7JzOfpSf7sEnoP/J4DFRHRnQ69HiIYSblwiLSYCA0hkhIi1uAAxFBhAIiO8KxbjAsRYYACJjPCeWI0LEIOBASQywjtiOS5AjAYGkMgI94v1uAAxHBhAIiPcyw1xAWI8MIBERriPW+ICuCAwgERGuIeb4gK4JDCAREbEPrfFBXBRYACJjIhdbowL4LLAABIZEXvcGhfAhYEBJDIidrg5LoBLAwNIZET0c3tcABcHBpDIiOjlhbgALg8MIJER0ccrcQE8EBhAIiOih5fiAngkMIBERjQ8r8UF8FBgAImMaDhejAvgscAAEhkReV6NC+DBwAASGRE5Xo4L4NHAABIZ4TyvxwXwcGAAiYxwjsSlmqcDA0hkhP0kLqd5PjCAREbYR+LyUxKYAImMCJfE5VwSmDNIZESoJC7nJ4E5i0RGWCVxqZ0E5jwkMiJYEpe6SWBqIZER9ZG41E8CUweJjKiNxCU4Eph6SGTE2SQuwZPABEEiI2pIXKyRwARJIiMkLtZJYCyQyHiXxCU0EhiLJDLeI3EJnQQmBBIZ75C4hEcCEyKJjPtJXMIngQmDRMa9JC72kMCESSLjPhIX+0hgbCCRcQ+Ji70kMDaRyMQ+iYv9JDA2ksjELomLMyQwNpPIxB6Ji3MkMA6QyMQOiYuzJDAOkchEP4mL8yQwDpLIRC+JS2RIYBwmkYk+EpfIkcBEgEQmekhcIksCEyESmYYncYk8CUwESWQajsSlYUhgIkwiE3kSl4YjgWkAEpnIkbg0LAlMA5HIOE/i0vAkMA1IIuMciUt0kMA0MImM/SQu0UMCEwUkMvaRuEQXCUyUkMiET+ISfSQwUUQiEzqJS3SSwEQZiYx1EpfoJYGJQhKZ4ElcopsEJkpJZOoncYl+EpgoJpGpncQlNkhgopxE5lwSl9ghgYkBEpnTJC6xRQITIyQyEpdYJIGJIV6OjMQlNklgYowXIyNxiV0SmBjkpchIXGKbBCZGeSEyEpfYJ4GJYW6OjMTFHSQwMc6NkZG4uIcExgXcFBmJi7tIYFzCDZGRuLiPBMZFYjkyEhd3ksC4TCxGRuLiXhIYF4qlyEhc3E0C41KxEBmJi/tJYFwsmiMjcfEGCYzLRWNkJC7eIYHxgGiKjMTFWyQwHnFWZOYQ0RYiSo/E3ETUmYi2E9HfIXHxFGLmyExEdAhABjNTRCYU50VEhwGkAzgJYBkzDw1yuzgAFwNoDOAYgL0cxB8PERkAtgO4EgABeJ+ZfxHi8oUNiIgBvMXMo5yeS/ZgvMcI/DcOwC1ENPh8N6JqDxiGscMwjDIAVQC+AbADwB4A2jCMUqXUJ0R0Rx3zPQGgHarjwgDa2/R7iBggezAeQ0S9ASwE0BRAIoAjANoz84+BfycAc4joPiKK69KlC7p160a9e/dGy5YtT41TWFiITz75BFu3btV5eXlERGVa62eZefoZc3UGsC0wTxmqw3QXM38doV9XnEck92AkMB5ERIkAngMwBkAjAP/LzEOJ6Dql1ErDMJrcfffdasiQIfD5fPWOV15ejvnz52PFihUM4Fut9U0A9uH0oVElgGkA/sDMpjO/lQiWBEZERGBvZhGAVgD+CuC+q666Sk+ZMkUlJCRYHq+kpASTJ082v/32WwIwH8C9AHZC9lqiigRGRExgb+afAK4fM2YMBg8+7ykZS15++WVeuXIlAXgLwGjZa4kucpJXRFJfANc/+OCDtsQFALKysigzM5MB/BpAG1sGFTFJAuNhRBRHRIt79erFgwYNsnXsrKwsatOmDSul1to6sIgpEhgPI6IFjRo1SnzyyScdOWydPn26AeBiInrKifFF9JPAeBQRNWLmO7OyslQwzxSFIi0tDbfddhsppbIdmUBEPQmMd02Nj4/nm2++2dFJRowYAWZOISJ7j8FETHDmoUtEPaXU/X379jXqu92sWbNQUFBw6vuSkhL4/X4sWLAgqHmSkpLQsWNH/c0330wCsCzkBYuYJIHxKK118xtvvLHe240fP/7U1wcOHMCECRPwyCOPWJqrW7duKj8//3KraxSxTw6RPIiIOgFQXbp0CXqb0tJSTJs2DcOHD0fPnj0tzde7d29orZsE3oYgPEQC40294+LizGBP7lZWViI3Nxd9+vRBZmam5cnatm1b8+XFljcWMU0OkbwpWangHltM08QLL7yA1q1bY8SIEad+vm/fPixduhRHjx5F165dceuttwYzXEpoyxWxSvZgvKlCax3UDefMmQO/349x48b95Odt2rTBww8/jOzsbOzevTvYecusLVPEOtmD8aYdJ0+eNLTWqGtPZsGCBdizZw9mzJgBwzj3CadNmzbhww8/RH1PdR84cKDmy31hrFnEIAmMN20DwF9//TVddtll573BwYMH8c4776B58+bIzj79OrlWrVqd+r5nz57o2bMnfve73+Gmm26qdbKNGzdCKVVmmuZJO38JEf0kMB7EzGwYRsn69evTawtMixYtsGxZ7S9b2bVrFzZu3IiTJ0+ie/fudc63bds2BpAfzppFbJLLNXgUEc1t3Ljx6HfffbfeF9uFw+/3Y+jQoay1foCZX3dyLhEcuVyDiIQny8rK1I4dOxydZPHixWDmKgBvODqRiEoSGI9i5lKl1MaZM2c6djGoyspKvPvuu5qZ3wjmEwiE+0hgPExrfduRI0f41VdfdWT8nJwcbZrmUQAPOzKBiHoSGA9j5lJmfmjp0qXYsmWLrWMvWbIEu3btIq11puy9eJcExuOY+Q0A23Jzc22LzLJly/Daa68B1ScSP7VlUBGTJDAeF/is6G7MvC0nJwfhHC75/X5Mnz5d/+UvfwGAPACjGvqzsEXDksB42FkfRH8NgNHLli3z33vvvabVZ5fWrl2Le+65R2/atKkCwABmvhSnPwtbIuNR8joYjzorLtPO+HmqUmq51vr6xo0b6379+hl9+vRB586df/K2Ar/fj507d2LdunXYsGGDWVVVRUT0ATPfw8wnzxhvNYABAIYw8wcR+wVFreRzkYSjaovLWbdJBfCcUmqY1joNAMXFxZlKKZimSX6/XwHQSqlDWuvXA2Od960AEpnoIoERjgkmLrVsdzWAq3H6c6Y3MnOehe0lMlEikoGR9yJ5SKhxAQBm3o7qz5oOCTMPDERmMRFJZDxCTvJ6RDhxsQszD4Sc+PUUCYwHRENcakhkvEUC43LRFJcaEhnvkMC4WDTGpYZExhskMC4VzXGpIZFxPwmMC8VCXGpIZNxNAuMysRSXGhIZ95LAuEgsxqWGRMadJDAuEctxqSGRcR8JjAu4IS41JDLuIoGJcW6KSw2JjHtIYGKYG+NSQyLjDhKYGOXmuNSQyMQ+CUwM8kJcakhkYpsEJsZ4KS41JDKxSwITQ7wYlxoSmdgkgYkRXo5LDYlM7JHAxACJy2kSmdgigYlyEpdzSWRihwQmiklcaieRiQ0SmCglcamfRCb6SWCikMQleBKZ6CaBiTISF+skMtFLAhNFJC6hk8hEJwlMlJC4hE8iE30kMFFA4mIfiUx0kcA0MImL/SQy0UMC04AkLs6RyEQHCUwDkbg4TyLT8CQwDUDiEjkSmYYlgYkwiUvkSWQajgQmgiQuDUci0zAkMBEicWl4EpnIk8BEgMQlekhkIksC4zCJS/SRyESOBMZBEpfoJZGJDAmMQyQu0U8i4zwJjAMkLrFDIuMsCYzNJC6xRyLjHAmMjSQusUsi4wwJjE0kLrFPImM/CYwNJC7uIZGxlwQmTBIX95HI2EcCEwaJi3tJZOwhgQmRxMX9JDLhk8CEQOLiHRKZ8EhgLJK4eI9EJnQSGAskLt4lkQmNBCZIEhchkbFOAhMEiYuoIZGxRgJTD4mLOJtEJngSmDpIXERtJDLBkcDUQuIi6iORqZ8E5jwkLiJYEpm6SWDOInERVklkaieBOYPERYRKInN+EpgAiYsIl0TmXBIYSFyEfSQyP+X5wEhchN0kMqd5OjASF+EUiUw1zwZG4iKcJpHxaGAkLiJSvB4ZzwVG4iIizcuR8VRgJC6ioXg1Mp4JjMRFNDQvRsYTgZG4iGjhtci4PjASFxFtvBQZVwdG4iKilVci49rASFxEtPNCZFwZGImLiBVuj4zrAiNxEbHGzZFxVWAkLiJWuTUyrgmMxEXEOjdGxhWBkbgIt3BbZGI+MBIX4TZuikxEAkNEcwBkBL5eZOO40yBxES7kRGSIKJ6IPg58O5KIsu0Yty6R2oPpdsbX19gxoMRFuJ0DkfEDuPaM7zvaMGadIhWYjwHowNebwh1M4iK8ws7IMLMG8HXg2zIAG8NcXr0iFZhNAI4DKEd1bEImcRFeY/OezPrAfzWAz8Icq16RCsxWAD5U76JtDXUQiYvwKhsj8wmq914SAPzHjrXVhZjZ6TlARATgGKp/qWRmrrSw7ZUALgLQGxIX4XFEtBrAAABDABgANjPztxa27wxgN4CvmPkyZ1Z5ms/pCQCAmZmI/gOgqcW4JAP4F4DmgR9JXISnMfPAQGQWB360nYi6B86vBCMPwAkA/3ZkgWdx7BCJiDKI6DEiWubz+b41DONKwzBa+Xy+74joH0SUTURt6hlmFoCUM753/KSUOI2q/ZyIXjMM40ufz1fk8/l+9Pl8h5RS24hoNhHd3NDr9KClZ3zdCcBjdd2YiC4nohwiWmcYxvdGteGGYRQQ0XtE9CARNXZiobYfIhHRz5VSM7TWlzVq1Mhs2bIlOnXqZKSnpwMAjhw5gm+++Ub/8MMPXFlZaSilCrTW05h53lnj9AOwDEASgEoApQD6M/MXti5YnIOIWgGYFYiLkZaWZrZv397Xpk0bxMXFobKyEoWFhVxYWGiWlpYaRFSptf4bgGxmLm3o9bsdEfVH9R5MAoA4ABUAujLzN2fchgD8Vin131rr5klJSf42bdoYHTp0oOTkZDAzDh06hLy8PP+hQ4eU3+8nItrEzL9lZtseyG0LDBG1Ukp9xMydLr/8cj1mzBjVoUOHOrfZt28fXn31Vb1jxw4ioh+01v2YOS9waLQH1YdGFQDeATCemY/bslhxXoE/ypcAPJqenm4OGTLEN2jQIChV+45uVVUVFixYgJUrV5plZWUA8AQzz4rMir2LiFoAeAtAH1Q/CO8A0J2ZNRH1VUotJ6Lk66+/nu6//35q2rRpneN99tlneOutt8zCwkJDKbVea53JzGVhr9OOwBDRbwD8+cILL+Tp06cbGRkZlrYvLS3FlClTzL179yoAUwC0BTAawEEA9zDz2rAXKepERBcGDntaZGVl0cCBAy2PMX/+fCxatIiJ6EutdXdmPuHAUkVA4AFhOIBXASQDeIKIrmTmUd27d9cTJ05U8fHxlsbMz8/H1KlTzePHj5vMPJiZV4W1xnADQ0QTAPx++PDhGDFiRFhjLV++HHPnzgUABvAGZK8lIojoYiLa3bJly/iZM2caycnJIY914MABPPbYY2Z5efmPWuu2djwKirrV7M0QUX8iUpMnT6YePXqEPJ7WGjNnzuR169YBwDBmfj/UscI6yUtEowH8fuzYsWHHBQBuv/12TJw4EQAIwF6Ji/OIKIWIvvzZz34W/8orr4QVFwBo2bIlXn/9dSM5OTldKfVV4FFWOIiZDxLRESJSL730UlhxAQClFJ544gkaNGgQALxHRH1DHivUDQPPEs0dOnQobr311lCHOUfv3r3x0EMPAUAuEXWybWBxXkqpf6SkpCTMmjXLqOtcixVJSUmYO3euoZS6CMBcWwYVtSKiQcw8PCcnh9q1a2fbuA8++CBdc801TEQriMgIZYyQ/6KUUh81a9YMo0aNCnWIWt1+++1o27atVkp9ZPvg4hQi+pXWuldubq7h89n7kqjk5GQ89thjCsAYIrLlDa7iXETUiIgW9u3bl7t27Wr7+JMmTVLx8fGNiejdULYPKTBE1FdrfeX06dNDqlowcnNzDWa+iIh+7dQcXqeUerlfv35s56PemW644QZ07NiRlVJvOzKBAIAX4+PjGz3++OOOHIrGx8cjOztbMfMviMjaszcIMTBKqT+0a9fOvPDCC0PZPCipqano2rUrK6WmOTaJhxFRpta6yUMPPeToOZKsrCyltb408NoaYTOl1H2ZmZnK7j3QM/Xo0QMpKSkmgOetbms5MESUqLXuPmrUKMf2Xmo8+OCDSmvdLohX/AqLlFIzOnXqpJOSkhydp0OHDkhLSwvpj1PUjYjuYOakX/3qV47PNWjQIIOI7ra6XSh7MHcYhsFXX311CJta06ZNGyQkJJgALP9iol6XDRgwICLvpr/22msNwzBuisRcHvPLFi1a6ISEBMcnuvPOO8HMCURk6bAllD+wW9LS0oJ9YxVWrlyJl19++dT38+bNwx/+8IegJwschvWzskBRNyJK0Vo36tWrV52301rj7rtPt72goACDBg3C999/DwCoqKjAyJEjUV5eXuc43bt3h9a6RfgrF2dSSvXo2LFjUEcS+/btw+jRo6F19V3XNE1MmTIFH30U3PMoCQkJiIuLM1H9Lu7g12jlxgCglOrWtm3boA/4br75ZmzZsgXHjx/H5s2b8dlnn+GRRx4Jer5OnToZhmF0sbpOryOia4joIBGtJ6JcIrqDiFoFXpcy0DAMnZqaWucYZz9tvWrVKqSnp58Kyrp163DttdeivsOsbt26gZkNIqr79eriHES0mYh2E9FrRHQ/EXUlorjAP7e66qqrghqnTZs2aN26NbZs2QKg+oG+devW6Ncv+MfujIwMBnCDlfVbPjNERCkpKSn13zAgISEBN9xwA+bNm4etW7ciNzcXjRo1Cnr7Cy64AMzckoicv3CN+zCq38/VC9UXGWqE6rfqHzKM4E6hNWrUCJWVlWBmfPnll+jVqxcqKioAAP/4xz8wbty4esc44+XqxfK6u5B1RvWpAg0gkYjyiSiuvvcYnWnw4MFYsmQJ/H4/du/ejenTp1taQHJysgGg7kels4Ry6tnyHb1///7IysrC5MmTYfWZJ2YGEWkAP1qd1+MyUP3/igBUBX5mACgBcARB/n9MSkpCRUUFNm3ahL59+8Lv96O8vBx79uyBYRiw+BT3YZy+NrMITs1TwyaqrwiZgOrLzx4E0MnKW326deuGN954A3/729/w7LPPwuozT4G5LN3/LQeGmY+Wllp7R/67776L1NRUmKZpdTr8+OOPAHCAmX9meWMPC1wJcCGAXai+DutWAJ8zcwUR3eH3+z8IZpyawKxZswaTJk3C2rVrUVFRgVWrVuGWW24Jai2VlaeuMdaRmUus/zbeRUQrADQGsA7AZgBbmfkAABiGUVFUVGTpDG/nzp3Rvn17pKWlWV7LsWPHNKofoIJm+RyM1nrL3r17/cHe/oMPPsDJkyeRnZ2NpUuX1r/BWfLy8kzTNHdY3tDjmPkLZr6Cme9m5j8z86fMXBH4539qrdXhw4frHScpKQmff/45mjZtivT0dCQmJqKkpATbtm1D377BvUVl8+bNICK/xMU6Zr6NmW9i5qeZeUVNXAK+3blzp6Xx9u3bh7Zt24a0luLiYgJg6coGoTyL9OGRI0eCOoD//PPPsWbNGowfPx5dunRBRUUFCgoKLE128OBBAJC3DNiImSuUUpUbN9Z/XaHExEQsWbIEmZmZp75fvXo1rrvuuqDPpW3duhVKqf1hLVqcQ2u9OT8/P+gHewD49ttvcfHFF1ueq7y8HH6/X+H0pTqDEkpgVmitUd8f56FDhzB79mxMmDDh1LMMgwYNwpIlS4KeKC8vDydOnDAAyEvNbcbMO1avXl3vMWtSUhKYGTXvc0lMTMT3338f9OERAGzZssU0TXN16KsVtfjroUOHfMePB3fRgaKiIjRu3BiJiYmWJ1q4cCGUUmXMXP9u7xlCuh6MUmpdq1atrp8zZ46jr+bNzs7WX3311VemaV7h5DxeRER9AKyfP38+6nu6Ohy7du3CpEmTAKCZ1T9OUT/DMI4MHDgw9eGHH3Z0nuHDh5vl5eWvMHPwrzFBiO9FYubx3333nWH1cMeKoqIi7N69W2mtJzo2iYcx8walVPGf/vQnR5/+nzNnjqmU2iFxcYbW+s//+te/dFVVVf03DtHatWtRXl6uAFi+L4YamO1KqY1Tpkwxa14ZaLdJkyaZRJTHzNbPDIugaK1Hffrpp7Rr1y5Hxl++fDn27duntNbDHJlAAMBUv99//JlnnnHkjlheXo4//vGPGsBfmPmY1e1Dfi+K1vqWY8eO6dmzZ9v+CDh//nwcOHCAtNY32T22OI2ZVxDRqtzcXH3GU8m2KCoqwquvvsoAnmPmfFsHF6cws8nMt23fvl1t2LDB9vGnTp2qtdY/AhgbyvYhB4aZjzHzXWvWrMHbb9t3DvbDDz/EwoULAWAsM8szDw5j5jtOnDhRMnbsWNOu3eySkhJkZWVpAP9hZjnEdRgzbwAw+/nnn+cdO3bYNu6zzz7LeXl5HPi0j5B2JOy46Pf9AF7r378/xo0bR+FcdvHNN9/E+++/DwATmfnZsBYmgkZETZVS+U2aNEl58cUXjRYtQn9fYl5eHiZNmqRPnjy5V2t9KTNbf3WlCIlS6m8ARmRlZZGVZ/nOVlVVhaefflp/8cUXANCPmdeFOpZdH1vycyJamJKSYuTm5hpWr5C2f/9+PPXUU2ZxcTGY+SFmfj3sRQlLiChFKbWZmS+955576J577rG0vdYas2fP5jVr1pBS6l9a6wGhPuqJ0BHRMwAmtW/fXufk5BhNmjSxtP2WLVvw3HPP6aqqquPM3I+Zt4a1Hhs/eK2xUmq11rp369atzV//+tdG796969xm27ZtePPNN829e/caSqldWuv/YuYiWxYkQkJETwB4NjExEQMHDjRGjBiBuq43cuTIEbzxxhu8fv16Nk3zJDP/KpyPuRDhI6LOgQ9BbHnFFVfw6NGj6/wQRK01Vq5ciffee888fPiwQUSLmPmXdux9OvHRsVcT0UvMfINSChdccIF5ySWXGKmpqUREOHr0KAoLC/0lJSXKNE0ioi2Bj6u0/wyVCAkRJQJ4Rin1gNY6JTEx0d+6dWt10UUXKZ/PhxMnTuC7774z9+/fjxMnThhKqSKt9UwAL8ghUfQgoiFKqela684+n083a9ZMt23b1peYmAhmxuHDh7mwsNA8duyYAcBk5mUAHmPm/2fbGpzaiw1csyITwG1KqR5E1CQw33Gt9VYAKwEsP+P9MSIKEdEVAIYCuNEwjDao/izkKtM081H9vpS/M3NhAy5R1IOI0gDcCWCgUqoLESUBMLXWh5n5EwBLAax14pDWscAIIURErskqhPAmCYwQwjESGCGEYyQwQgjHSGCEEI6RwAghHCOBEUI4RgIjhHCMBEYI4RgJjBDCMRIYIYRj/j+RJxoGYAeA6QAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 266.457x266.457 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig_2_9 = DAG(\n",
    "    ebunch=[\n",
    "        (\"W\", \"Y\"),\n",
    "        (\"X\", \"W\"),\n",
    "        (\"Z_1\", \"X\"),\n",
    "        (\"Z_1\", \"Z_3\"),\n",
    "        (\"Z_2\", \"Y\"),\n",
    "        (\"Z_2\", \"Z_3\"),\n",
    "        (\"Z_3\", \"X\"),\n",
    "        (\"Z_3\", \"Y\"),\n",
    "    ]\n",
    ")\n",
    "\n",
    "fig_2_9_node_pos = {\n",
    "    \"W\": (2, 0),\n",
    "    \"X\": (0, 0),\n",
    "    \"Y\": (4, 0),\n",
    "    \"Z_1\": (0, 4),\n",
    "    \"Z_2\": (4, 4),\n",
    "    \"Z_3\": (2, 2),\n",
    "}\n",
    "\n",
    "fig_2_9.to_daft(node_pos=fig_2_9_node_pos).render();"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Study question 2.4.1\n",
    "\n",
    "(a) For each pair of nonadjacent nodes in this graph, find a set of variables that d-separates that pair. What does this list tell us about independencies in the data?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(W ⟂ Z_1 | X)\n",
       "(W ⟂ Z_2 | X)\n",
       "(W ⟂ Z_3 | X)\n",
       "(X ⟂ Y | Z_3, W, Z_2)\n",
       "(X ⟂ Z_2 | Z_3, Z_1)\n",
       "(Y ⟂ Z_1 | Z_3, W, Z_2)\n",
       "(Z_1 ⟂ Z_2)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print_implied_conditional_independencies(fig_2_9)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the data to be consistent with the causal model represented by the graph, they should reflect the above conditional independencies."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(b) Repeat question (a) assuming that only variables in the set {Z_3, W, X, Z_1} can be measured."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(W ⟂ Z_1 | X)\n",
       "(W ⟂ Z_3 | X)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig_2_9.latents = set(fig_2_9.nodes()).difference({\"W\", \"X\", \"Z_1\", \"Z_3\"})\n",
    "\n",
    "print_implied_conditional_independencies(fig_2_9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig_2_9.latents = set()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(c) For each pair of nonadjacent nodes in the graph, determine whether they are independent conditional on all other variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(W ⟂ Z_1 | Y, Z_3, X, Z_2)\n",
      "(X ⟂ Y | Z_3, W, Z_1, Z_2)\n",
      "(X ⟂ Z_2 | Y, Z_3, W, Z_1)\n",
      "(Y ⟂ Z_1 | Z_3, W, X, Z_2)\n"
     ]
    }
   ],
   "source": [
    "for x, y in pairs(fig_2_9.nodes()):\n",
    "    print_independent(\n",
    "        graph=fig_2_9, x=x, y=y, z=set(fig_2_9.nodes()).difference({x, y})\n",
    "    )"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(d) For every variable V in the graph, find a minimal set of nodes that renders V independent of all other variables in the graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-04-14T13:04:48.776503Z",
     "start_time": "2023-04-14T13:04:48.769949Z"
    }
   },
   "outputs": [],
   "source": [
    "def print_markov_blanket(graph, node):\n",
    "    print(f\"Markov blanket of {node}:\", sorted(graph.get_markov_blanket(node)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Markov blanket of W: ['X', 'Y', 'Z_2', 'Z_3']\n",
      "Markov blanket of X: ['W', 'Z_1', 'Z_3']\n",
      "Markov blanket of Y: ['W', 'Z_2', 'Z_3']\n",
      "Markov blanket of Z_1: ['X', 'Z_2', 'Z_3']\n",
      "Markov blanket of Z_2: ['W', 'Y', 'Z_1', 'Z_3']\n",
      "Markov blanket of Z_3: ['W', 'X', 'Y', 'Z_1', 'Z_2']\n"
     ]
    }
   ],
   "source": [
    "for node in sorted(fig_2_9.nodes()):\n",
    "    print_markov_blanket(graph=fig_2_9, node=node)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(e) Suppose we wish to estimate the value of Y from measurements taken on all other variables in the model. Find the smallest set of variables that would yield as good an estimate of Y as when we measured all variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Markov blanket of Y: ['W', 'Z_2', 'Z_3']\n"
     ]
    }
   ],
   "source": [
    "print_markov_blanket(graph=fig_2_9, node=\"Y\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(f) Repeat question (e) assuming that we wish to estimate the value of Z_2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Markov blanket of Z_2: ['W', 'Y', 'Z_1', 'Z_3']\n"
     ]
    }
   ],
   "source": [
    "print_markov_blanket(graph=fig_2_9, node=\"Z_2\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(g) Suppose we wish to predict the value of Z_2 from measurements of Z_3. Would the quality of our prediction improve if we add measurement of W? Explain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig_2_9.is_dconnected(start=\"W\", end=\"Z_2\", observed=\"Z_3\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When conditioned on Z_3, the variables W and Z_2 are connected, hence they are not independent. Thus, it would likely be advantageous to add W."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.5 Model Testing and Causal Search\n",
    "\n",
    "#### Study question 2.5.1\n",
    "\n",
    "(a) Which of the arrows in Figure 2.9 can be reversed without being detected by any statistical test? \\[Hint: Use the criterion for equivalence class.\\]\n",
    "\n",
    "This question can be [answered](https://github.com/behrman/cis/blob/master/solutions_dagitty.md#study-question-251) with the R function `dagitty::equivalenceClass()`, which generates a complete partially directed graph (CPDAG) that represents all the graphs that are Markov equivalent to the input DAG. There is no comparable functionality in pgmpy as of version 0.1.22 (`DAG.to_pdag()` has not been implemented)."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(d) Write down a regression equation for Y such that, if a certain coefficient in that equation is nonzero, the model of Figure 2.9 is wrong.\n",
    "\n",
    "Here are conditional independencies for `fig_2_9`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-04-14T13:08:49.797815Z",
     "start_time": "2023-04-14T13:08:49.787605Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(W ⟂ Z_1 | X)\n",
       "(W ⟂ Z_2 | X)\n",
       "(W ⟂ Z_3 | X)\n",
       "(X ⟂ Y | Z_3, W, Z_2)\n",
       "(X ⟂ Z_2 | Z_3, Z_1)\n",
       "(Y ⟂ Z_1 | Z_3, W, Z_2)\n",
       "(Z_1 ⟂ Z_2)"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print_implied_conditional_independencies(fig_2_9)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we assume that the model is correct and linear, then the conditional independence of X and Y conditioned on {W, Z_1, Z_3} would imply that the coefficient b for this linear regression\n",
    "\n",
    "Y = a + b \\* X + c \\* W + d \\* Z_1 + e \\* Z_3\n",
    "\n",
    "should be zero."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(e) Repeat question (d) for variable Z_3.\n",
    "\n",
    "Similarly, the conditional independence of Z_3 and W conditioned on X would imply that the coefficient b for this linear regression\n",
    "\n",
    "Z_3 = a + b \\* W + c \\* X\n",
    "\n",
    "should be zero."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(f) Repeat question (e) assuming the X is not measured.\n",
    "\n",
    "Assuming that X is not measured, we will convert it to a latent variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-04-14T13:17:48.265174Z",
     "start_time": "2023-04-14T13:17:48.229907Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(W ⟂ Z_2 | Z_3, Z_1)\n",
       "(Y ⟂ Z_1 | Z_3, W, Z_2)\n",
       "(Z_1 ⟂ Z_2)"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig_2_9.latents = \"X\"\n",
    "\n",
    "print_implied_conditional_independencies(fig_2_9)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, Z_3 is not conditionally independent with any other variable, so no such regression exists."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-04-14T13:18:51.618773Z",
     "start_time": "2023-04-14T13:18:51.611246Z"
    }
   },
   "outputs": [],
   "source": [
    "fig_2_9.latents = set()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(g) How many regression equations of the type described in (d) and (e) are needed to ensure that the model is fully tested, namely, that if it passes all these tests, it cannot be refuted by additional tests of this kind. \\[Hint: Ensure that you test every vanishing partial regression coefficient that is implied by the product decomposition (1.29).\\]\n",
    "\n",
    "This question can be [answered](https://github.com/behrman/cis/blob/master/solutions_dagitty.md#study-question-251) by using the R function `dagitty::impliedConditionalIndependencies()` to generate a basis set of conditional independencies. There is no comparable functionality in pgmpy as of version 0.1.22."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "## 3 The Effects of Interventions\n",
    "\n",
    "### 3.3 The Backdoor Criterion\n",
    "\n",
    "Figure 3.8. Causal graph used to illustrate the backdoor criterion in the following study questions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-04-14T14:03:58.189303Z",
     "start_time": "2023-04-14T14:03:56.224500Z"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARgAAAEYCAYAAACHjumMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtbElEQVR4nO3deXgUVb7/8c/3VGclISxBEYiAhIjKoogoCC5cZGQUF7iIKIOM4gKKDy4/o8g2IIiOI47MsFyVcRRF9AoXhUERZRBEBNlxlBhCAiJLiGEJSQhd5/v7ozsYWZKu6qpOL9/X8/iYYNc5h8fk3VXVXdXEzBBCCDeo2l6AECJ6SWCEEK6RwAghXCOBEUK4RgIjhHCNBEYI4RoJjBDCNRIYIYRrJDBCCNdIYIQQrvG4OTgRNQTQD8B1ABoAIACHAXwJ4ENm3uvm/EIIHyLKBPDfAK4EUAeACWA/gM8A/B8zH3NlXqevRSKiVABTlFL3aK3rxMXFmenp6UhNTVUAcOzYMV1YWIiKigqDiMqZeR6A/8fMhY4uRIgY54/Ky0TUm5k9iYmJZqNGjZCcnKy01nzo0CH9yy+/KNM0lVKqWGs9DcBzzHzCsTU4FRgiIgCzAAxNTk7WvXv3Nu644w4kJyef8fHl5eVYsGABFi1aZB45csQgog+ZeQAzm44sSIgYRUSpSqmlWuurGjVqZPbv39+44YYb4PGc+YCluLgYb7/9Nv79739rr9fLzDyamac4shYnAkNE7ZRSy5VS9R999FF1/fXXW9p+3bp1eOGFF3RFRcUxZr6JmVcGvSghYhAR3UVEb9atW1dNmDDBuOCCCyxt//7772POnDlMRDu01t2YeX9Q6wk2MER0PYDPLrzwQpowYYI62x5LTSoqKjB58mS9fv16AtCfmT8MamFCxBgiGgVgUu/evfmhhx4ipey9hnPw4EE89dRT5sGDByuYuT0z59peUzCBIaLLAazt1q0bZWdnk+2Bqpg1axYvWrQIAH7PzJ84MaYQ0Y6IHgXw12HDhuH3v/990ONprTFy5EgzPz+/gplb2t2TsR0YIjKIqLht27Z1Jk+e7OjL3S+99BJ/+eWXJ5i5gVtnt4WIFv6TuTmDBw+m/v37Ozau1hpDhw41i4qK8kzTzLIzhu0wENHc+Pj4OuPHj3f8vTSPP/44paSkGEqpJU6PLUS0UUp90aJFC+1kXPzj4sUXXzSYOZOInrY1hp2NiKgFM/fPzs5W8fHxdoaollIKEyZMMLTW3Ymoq+MTCBEl/IdGzZ577jnDjfHT09MxaNAgIqLniCjO6vZ29z6mNmzY0HvFFVfY3LxmmZmZaN68uUlEL7k2iRARTin1TJcuXZCWlubaHHfccQc8Hg8BeNbqtpYDQz4333HHHa6+CxgAhgwZYjDzVURUx+25hIg0RHSZ1rrx0KFDHXmBpTrXXnutUkqNsLqdnT2Y/wJg3HjjjTY2taZTp07weDwM4C7XJxMi8gyrX7++Nz093fWJBg8eDK11A6tP9nb2QvqkpqaaSqmAt12yZAnmzp2LevXqoby8HOeffz6ys7MRF1fzId0555yjf/755xsAvGZjrUJELaVUl1atWln6HV61ahUWLFiAEydOwOv1onv37hg4cGCN29WvXx+GYWjTNG8F8G6g81kODBF1adGihaUTSgUFBRg8eDB69uwJrTUeeugh5Ofno3Xr1jVum5WV5dm3b18nq+sUIgZc0K5du4Af/Pnnn2PRokV49tlnkZ6ejvLycnz22WcBb9+gQQNdWFjYCxYCY/kQSSnVoGHDhpaO+fLz81H5luW9e/eCmdG0adOAtk1PT6+8gFIIUQUzxwf6e1RaWoo33ngD2dnZqDykSkxMRJ8+fQKeLy0tTQE4x8oa7RwiGVbfgrxr1y5MnToVXq8XRUVFGDt27FkvgjxtMsMAAFdeghMiwpH/96NGa9asQVZWFho3bmx7MqUUwWIzLAeGmUtLS0sDfnxhYSHS0tIwbdo0AMAXX3yBefPmoW3btgFtX1JSAgDHra5TiBhgHjp0KKDCFBQUwOqFj6cqKyvTAI5a2cbyIZLWetvOnTsDvqVCQUEBMjIyTn7fsmVLHDp0KOD5cnNztWma2y0tUogYoJQ6sG3btoAem5iYCK11UPMVFhYCwFor29h5mXrpwYMHAz4Hs3PnzpOBYWZ8/vnnuPTSSwOebPfu3RrACquLFCLamaa5IScnJ6An+8svvxxfffUViouLAQAnTpzAp59+GvBcXq8X5eXlBgBLdzmwcw5mvtfrfT03NxeZmZk1PrigoADbtm3Dt99+CyJCVlYW7r333oAmKioqQmlpqQfAXBvrFCLaLdizZ08fr9d71ptJVcrKysLAgQMxbtw4MDO8Xi+uu+66gCdaunQpiMirtbZ06wZbV1MbhvHDxRdf3Pr555939abhL774In/11VcHTNO0f2ZKiCjlf1d9+T333BPfr18/V+f64x//aB48eHAxM99qZTtbgdBaP/vdd9/RkSNH7GwekIqKCqxevZq11i+6NokQEYx9/nf+/Pmu3mY2Ly8PBw8eNAA8ZnVbW4Fh5g+JaN/YsWNd+4tNnjxZa61LAEx1aw4hosCjR48exZtvvunK4FprjB071iSiNcycZ3V724c4WuueO3bsUJ984vxN59auXYv169crZr6Znf7YAyGiCDMXMfOI+fPn808//eT4+DNnzuQjR45oZu5lZ3vbgWHm/wCYPH36dN64caPdYU6Tm5uLSZMmMRG9KTf/FqJmzDyDiNaNHDlSV75K5ITFixdjyZIlxMxDmNnS+18qBXWSlplHA5gzbtw4XrVqVTBDAQA2b96MJ554gpmZmPmjoAcUIkZorT+pqKhQ999/v3ZiT+b999/HzJkzAeAZZg742qNTOfWxJa8CGHHZZZfp0aNHW77Lndfrxcsvv8wrV64kIprLzOkAbgDQl5kXBL1AIaIYEY0HMA7ABKXUjcx8Rd++fWnIkCGWxzp8+DBGjx5t5ufnKwAjmfnVoNbm4AevdSOixYZhpPTs2VPdc889SElJqXab8vJyvPPOO1iyZImuqKg4zsz9mXmxf7ylkMgIUa0qcfkTM4/3/9kwIpqWkpKCfv36GbfeemuN75MpLCzE66+/zmvWrAGAn7TWPYL5uJKT63PyHKr/0x2fU0oN11rXa9SokTczM9PTvn17NGrUCESEoqIibNmyBbm5ud79+/cbRFSqtf4HgCeZ+fgp40lkhDiLM8Wlyn9rCOCvRNSfiOKaNGmiL7zwQqNt27ZIS0uDaZr4+eefsWXLFs7LyzOLi4s9Sqm9WuvJzPw3x9bo1os0RNQDwBCl1JUAzmfmBP+fHwewR2u9FsBcZv64hnEkMkKcorq4nPI4AnAPgNsMw+jIzOcyczwATUTlAHZorb8GMJ2ZNzu+zlC9CkxEBwA0YmbL9w+VyAjxq0DjUs32DOCfzDzE2ZWdztW3+jvF/xr8ZwDmE9Httb0eIWpLsHEJtYgIDCCRESLS4gJEUGAAiYyIXZEYFyDCAgNIZETsidS4ABEYGEAiI2JHJMcFiNDAABIZEf0iPS5ABAcGkMiI6BUNcQEiPDCAREZEn2iJCxAFgQEkMiJ6RFNcgCgJDCCREZEv2uICRFFgAImMiFzRGBcgygIDSGRE5InWuABRGBhAIiMiRzTHBYjSwAASGRH+oj0uQBQHBpDIiPAVC3EBojwwgERGhJ9YiQsQA4EBJDIifMRSXIAYCQwgkRG1L9biAsRQYACJjKg9sRgXIMYCA0hkROjFalyAGAwMIJERoRPLcQFiNDCAREa4L9bjAsRwYACJjHCPxMUnpgMDSGSE8yQuv4r5wAASGeEcictvSWD8JDIiWBKX00lgqpDICLskLmcmgTmFREZYJXE5OwnMGUhkRKAkLtWTwJyFREbUROJSMwlMNSQy4mwkLoGRwNRAIiNOJXEJnAQmABIZUUniYo0EJkASGSFxsU4CY4FEJnZJXOyRwFgkkYk9Ehf7JDA2SGRih8QlOBIYmyQy0U/iEjwJTBAkMtFL4uIMCUyQJDLRR+LiHAmMAyQy0UPi4iwJjEMkMpFP4uI8CYyDJDKRS+LiDgmMwyQykUfi4h4JjAskMpFD4uIuCYxLJDLhT+LiPgmMiyQy4UviEhoSGJdJZMKPxCV0JDAhIJEJHxKX0JLAhIhEpvZJXEJPAhNCEpnaI3GpHRKYEJPIhJ7EpfZIYGqBRCZ0JC61SwJTSyQy7pO41D4JTC2SyLhH4hIeJDC1TCLjPIlL+JDAhAGJjHMkLuFFAhMmJDLBk7iEHwlMGJHI2CdxCU8SmDAjkbFO4hK+JDBhSCITOIlLeJPAhCmJTM0kLuFPAhPGJDJnJ3GJDBKYMCeROZ3EJXJIYCKAROZXEpfIIoGJEBIZiUskksBEkFiOjMQlMklgIkwsRkbiErkkMBEoliIjcYlsEpgIFQuRkbhEPglMBIvmyEhcooMEJsJFY2QkLtFDAhMFoikyEpfoIoGJEtEQGYlL9JHARJFIjozEJTqFJDBEdA6ARv6vm4RizlgViZGRuIQOEbX0f9mSiNLcni9UezCfVvn6mxDNGbMiKTISl9AhokQA2/3fXgPgf9yeM1SB2V3l670hmjOmRUJkJC6hxczlAI74vzUB5Lg9Z6gCswJABQAG8GWI5ox54RwZiUut2ej/9zEA69yeLFSBWQ+gHEAJ5BAppMIxMhKXWrUCwAkA8fD9XroqVIHZCCAJgAfAtyGaU/iFU2QkLrVuHXyHR8cB/Oz2ZMTMbs/hm4joZwBpAFI4VJOK3yCipQBuADATQCcAv2PmX0IwbxsAcwHsANAPEpdaQ0TpAAoBrGTma1yfLxS/60RUH8B8AA0BXM/MRa5PKs6IiIoANIBvN/ljZu4X4HZxAJoDqAPgKICdgTxREJEB3x5sWwAE4ENm/m+byxdBIqLzAGwDMA9ANjMfdXM+Vw6RiKghEb3h8Xh2KaVOAPgFwHUA2gE4qJTyejyePUT0DhE1dWMN4qwM/7/jANxIRLee6UHkc79hGJsMwzgG30n6HwFsgm9PRBuGcVgp9TUR3VLNfE8CuAC+uDCAVg79PUQAiOgSIvrY4/EUEpEJ32FRAwDDABxRSh03DGMHEf2ZiJIcn9/JPRgiylJKzdNaX5qammp27tzZ6NSpEzp16oTExEQAQEVFBTZu3Ih169Zh7dq1ZnFxsaGU2q61vpuZXT/pFOuIqCt8z14N4TsvdghAq8pDJSIiADOI6I9EFNeuXTt07NiRunbtisaNG58cJz8/H19//TXWr1+vc3JyiIiOaa2nMPOkKnO1AbDBP88x+MJ0BzNvh3AVEfVRSs3QWjdt3Lix2blzZ6Nz585o164dlPLtV5SUlGDNmjVYv3491q1bZx4/flwR0b+Zub9TRxmOBYaIpgB4KiMjgx9++GF1ySWXBLRdbm4upk+frn/88UcC8AaAB+Qcjbv8z1QvABgKIAHA/zFzPyK6Sim1xDCMunfeeafq27cvPB5PjeOVlpZizpw5WLx4MQPYpbW+Dr73PlUeGpUDGA/gL8xsuvO3EgBARAlE9BEz9+rUqRM/9NBDdO655wa07apVq/Daa6+ZxcXFzMwjmHlm0OsJ9neZiAyl1BYAFz3wwAN000032Rrnyy+/xNSpU7XW+metdRtmPhbUwkSN/Hsz7wNoCuAfAP7Yvn17PWbMGFW5x2lFcXExRo8ebe7atYsAzAEwGMAWyF5LSBBRcyLampSUlDx27Fgj0Cf5U/3jH//AggULGMBirXWfoNYUTGD8cdkeHx/fcsaMGSo9PT2YtaCkpAQPPvigWVJS8ovWOoOZjwc1oKiRf2/mMwBXDx06FLfeesZTMpZMnz6dlyxZQgD+CeA+2WtxHxE1JaLcjIyMuL/+9a9GIHue1cnJycFTTz3FWuvPtdY32F5XMIExDOOLuLi4a19//XVVr1492+NUVVpaivvuu88sLS39j2ma7R0ZVJwVEfUC8OkDDzyAPn2CerL6DX9kGL7zO/mODSxOQ0SklDrYpEmTtL///e9G5TmWYOXm5uKJJ55grfWrzDzSzhi2V0JE/bXW10+ZMsWxuABAcnIypk6dajBzWyIa6djA4jREFEdE87t06cJOxgUAhg8fThkZGayUWu7owOJM3jQMo97UqVMdiwsAZGZmYsSIEQTgUSKydbxlazVEZBDR2z179uTMzEw7Q1SrcePGuOOOO4iIXgrFJeWxiojmJiQkJD311FPkxviTJk0yADQnomfdGF8ARNQBwOAnn3zS1nmzmvTs2RNZWVlaKbXMzvZ2c5dtGEacv26uGDRoEJKTkwHgebfmiGVElMDMtw8fPlwFe7x+NvXr18dNN91ESqlsVyYQIKK/ZWRkmF27dnVtjnHjxhla68ZE1MPqtrYCo5R6rHv37srJ3bEz6d27t6GUusfVSWLX2Pj4eL7++utdnWTQoEFg5lQicvYYTICI6jDz1YMHDzZqfrR9devWRWZmplZK/dnqtpYLQUQttNbp9957r9VNLRs4cCCYOYmI3P0tiEFKqXu7d+/u6g8m4Dun1rp1a01Eo9yeKwaNTEhI0FdddZXrEw0ePFhprS/zvxEzYHb2jfsnJCSY9erVs/zDOWvWLHzzzTeYPXt2QI+Pj49H3bp1zcOHD98GQE4WOkhrfc61114b8ONfeeUV5OXlnfy+uLgYXq8Xc+fOrXHbjh07qtzc3IttLVRUp2fz5s0D/oVfsmQJ5s6di3r16oGZ0bRpUwwZMuQ379A+mw4dOgC+yz2uBLAm0DntBOb6Jk2s31Z3//792LJlC7xeL0pLSyvPr9SoZcuWxubNm7tZnlCcFRFlAVDt2rULeJuRI0ee/Hrfvn14+umn8cgjjwS0bdeuXfHee+/VJSKSd2k7xzCMthdddFHARyEFBQW4++678bvf/Q6ALzjPP/88XnnlFdS0Y6KUQkpKirekpOQ2WAiM5UMkpVRms2bNLO+9vPvuuxgwYAAyMjKwa9eugLdr0aIFKaUyrM4nqtU1Li7OtHNy9/Dhwxg/fjwGDBiAK6+8MqBtWrasvM80mlueUJyV1jqtTZs2AT8+Pz8fzZv/+r+gd+/eKC4uRmFhYUDbN2rUSAGw9HK1nXMw8QkJCZa2KSgoQEFBAbp3746MjAwUFBQEvK1/Lnde5ohdKXZO0JeXl2PixIno1q0bevfubWfeVDsbibOy9NJ0QUEBzj///N/8WXx8PI4dC+yqnPj4eILv2rXAF2jlwQDAzBUnTpywtM2cOXNw9913g4gs78F4vV4A8FpbpahBmdba0gamaeLPf/4zmjVrhkGDBtmdV64vcxZXVFQE9MDCwkIkJib+5tSE1+tFcXExAr0Y8sSJEwzfbTsCZnnPQGu9Z+/eva0QYJy2b9+ODRs2IC8vDzNnzkRFRQVatGgR8Hy7du0CMx+wuk5RrU0nTpwwtNYIdE9mxowZ8Hq9GDFihOXJ9u3bV/nl7uoeJ6xRSpXs2LGjXiDvgSkoKPjN4REALFu2DO3btw/4fGhRUZGG75YbAbNz6LHip59+6oYAA/PWW29hzJgxuPTSSwH4Xn2oesKwJnl5eV6t9Wob6xRntwEAb9++nS666KIaHzx37lzs2LEDkydPhmFYf2V79erVUEodM03T2q6vqJbW+vtt27ZdBd+rO9XauXPnbwKzYcMGfPDBBxg/fnzA8x05csQA8LGVNdoJzIelpaXjysvLUdPx36ZNm+D1ek/GBfC9u7OsrAxHjx5Famr1h+RaaxQXFxsA/mVjneIsmJkNwyheuXJlg5oCs3//frz77rs455xzkJ396xtymzZt+pvvq7NhwwYGkBvMmsXpmPnL/Pz8KxDA73FBQQG2bt2KzZs3g5mRkZGBP/3pT2jWrFlAc+Xm5oKZCRbfLmLramrDMI726dMnZejQoZa3tWLhwoV44403TjBzgry86SwimlWnTp373nvvPVffbOf1etGvXz/WWt/PzG+4OVesIaJGAA688MILuPhid99m9Oyzz+pt27blmabZ2sp2tt7rr7V+fenSpa7f4+PDDz/0MvMCiYsrnjp27JjatGmTq5PMnz8fzFwBILB3V4qAMXOhUmrr7NmzrZ2xt6iiogJbt24lrfU4q9vavZhodFlZGQXyLk67li5dWnl49Lhrk8QwZj6slFr98ssvu/ZEUV5ejvfee08z82x5knCH1vrx7du3q+3b3bth4EsvvaSJqISZ37W6ra3A+G9nmT137lzev3+/nSGqVVJSgunTpzOAvzLzHscnEAAArfVNhw4d4tdee82V8SdMmKBN0zwC4GFXJhBg5mVEtHzcuHGm1bceBGLTpk34+uuvldba1kfN2L4cmplfIqLtTz75pBnoa/GB0Frj8ccfN5l5P2TvxVXMfJiZH/zoo4+wbp2zH1O8cOHCyt3q3rL34i5mvqm0tPTE+PHjHS1MUVERJkyYoP03EV9qZ4yg7regte585MiRow8++KBZXl4ezFAAUPk+C3Pfvn0VWuuO8oPpPmaeDWDDxIkTHYvMxx9/jNdffx0A/snMAV+3Iuxh5jJm7r5x40ZMmjTJkcgcPHgQw4YN06Zp7mDm2+yO48SnCjRUSm1PSEioP2bMGEsX0FWVl5eH0aNHmyUlJeXMfAkzB349gbCtymdFbwDQ8ZZbbsH9999vayyv14sXXnhBr1mzRgHIAZAFoC8zL3BqveLsiKgrEa0477zzaPLkyUbDhg1tjbN8+XK8+uqrWmudo7Vuz8y237/kyOci+e/tupCZe1911VX8yCOPUFpaYHe6LC0txcyZM3n58uWklFqttb7R7Y+zFD6nfhA9Ed1LRLPq1atHjz/+uFH1/Us1Wb58OaZPn66PHz9exsy3+c8NVH4WtkQmRIgoQym1nJkv6Nu3L911112Ij48PaNs9e/Zg6tSpevv27QTf55c/HOxRhNOf7HiTUupNrXV6ixYtzF69ehldu3bFqSU9fPgwVq9ejWXLlumcnByllDqitR7OzO84thhRrVPjUuXP05RSi7TWV9epU0f36NHD6NatG9q0afObywq8Xi+2bNmCFStWYNWqVWZFRQUR0QJmHlj1GU8iUzuI6FkiGg0goUOHDtyjRw915ZVXnnZZwJ49e/DVV19h6dKl5v79+w2lVIHW+lZm3uzIQpjZ8X8AdCOi1UqpUgCslNLx8fFmfHy8aRiG6f+zcqXURgA3urEG+afa/z/j4fuc6PHVPCYNwEylVBEADYDj4uK8CQkJXo/HY/q3N5VSewE8ByCumrGW+h9/e23/3WPtHwD3KKW2E1EFADYMw/T/LnqJSPt/F0uIaAmANk7P7+gezJn4P9jrBvg+6S8VQD8An7J8cmOtONueSwDbXQbgMvz6OdOrmTnHwvayJ1PL/O/8/S8AcwGsBvAnAJ+zix+M53pgTk5EdABAI/ZdzyBqgd24ODi/RCYMEBHD9wrfELfncvdjAUTYqO24AAAz94LvY2rnE9HttbEGEVoSmBgQDnGpJJGJLRKYKBdOcakkkYkdEpgoFo5xqSSRiQ0SmCgVznGpJJGJfhKYKBQJcakkkYluEpgoE0lxqSSRiV4SmCgSiXGpJJGJThKYKBHJcakkkYk+EpgoEA1xqSSRiS4SmAgXTXGpJJGJHhKYCBaNcakkkYkOEpgIFc1xqSSRiXwSmAgUC3GpJJGJbBKYCBNLcakkkYlcEpgIEotxqSSRiUwSmAgRy3GpJJGJPBKYCCBx+ZVEJrJIYMKcxOV0EpnIIYEJYxKXs5PIRAYJTJiSuNRMIhP+JDBhSOISOIlMeJPAhBmJi3USmfAlgQkjEhf7JDLhSQITJiQuwZPIhB8JTBiQuDhHIhNeJDC1TOLiPIlM+JDA1CKJi3skMuFBAlNLJC7uk8jUPglMLZC4hI5EpnZJYEJM4hJ6EpnaI4EJIYlL7ZHI1A4JTIhIXGqfRCb0JDAhIHEJHxKZ0JLAuEziEn4kMqEjgXGRxCV8SWRCQwLjEolL+JPIuE8C4wKJS+SQyLhLAuMwiUvkkci4RwLjIIlL5JLIuEMC4xCJS+STyDhPAuMAiUv0kMg4SwITJIlL9JHIOEcCEwSJS/SSyDhDAmOTxCX6SWSCJ4GxQeISOyQywZHAWCRxiT0SGfskMBZIXGKXRMYeCUyAJC5CImOdBCYAEhdRSSJjjQSmBhIXcSqJTOAkMNWQuIizkcgERgJzFhIXUROJTM0kMGcgcRGBkshUTwJzComLsEoic3YSmCokLsIuicyZSWD8JC4iWBKZ00lgIHERzpHI/FbMB0biIpwmkflVTAdG4iLcIpHxidnASFyE2yQyMRoYiYsIlViPTMwFRuIiQi2WIxNTgZG4iNoSq5GJmcBIXERti8XIxERgJC4iXMRaZKI+MBIXEW5iKTJRHRiJiwhXsRKZqA2MxEWEu1iITFQGRuIiIkW0RybqAiNxEZEmmiMTVYGRuIhIFa2RiZrASFxEpIvGyERFYCQuIlpEW2QiPjASFxFtoikyIQkMEc0A0Mj/9fsOjjseEhcRhdyIDBHFE9GX/m/vIaJsJ8atTqj2YDpW+fpyJwaUuIho50JkvAA6V/m+tQNjVitUgfkSgPZ//U2wg0lcRKxwMjLMrAFs9397DMDqIJdXo1AF5hsAJQBK4YuNbRIXEWsc3pNZ6f+3BvBtkGPVKFSBWQ/AA98u2nq7g0hcRKxyMDJfw7f3kgjgP06srTrEzG7PASIiAEfh+0ulMHO5hW3bAmgCoCskLiLGEdFSADcA6AvAALCWmXdZ2L4NgO8B/MDMF7mzyl953J4AAJiZieg/ABpajEsKgM8BnOP/I4mLiGnM3Msfmfn+P9pIRJ3851cCkQPgOICvXFngKVw7RCKiRkT0GBF97PF4dhmG0dYwjKYej+cnIvqUiLKJKKOGYV4BkFrle9dPSolfkc9tRPS6YRjfeTyeQo/H84vH4zmglNpARNOI6PraXmcM+qjK11kAHqvuwUR0MRFNIKIVhmHsMXwGGIaRR0QfENEDRFTHjYU6fohERLcppSZrrS9KSEgwGzdujKysLKNBgwYAgEOHDuHHH3/UP//8M5eXlxtKqTyt9XhmfvuUcXoA+BhAMoByAIcB9GTmbY4uWJyGiJoCeMUfF6N+/fpmq1atPBkZGYiLi0N5eTny8/M5Pz/fPHz4sEFE5VrrtwBkM/Ph2l5/tCOinvDtwSQCiANQBqADM/9Y5TEE4Aml1P/TWp+TnJzszcjIMDIzMyklJQXMjAMHDiAnJ8d74MAB5fV6iYi+YeYnmNmxJ3LHAkNETZVSXzBz1sUXX6yHDh2qMjMzq91m9+7deO211/SmTZuIiH7WWvdg5hz/odEO+A6NygC8C2AkM5c4slhxRv4fyqkAHm3QoIHZt29fT58+faDU2Xd0KyoqMHfuXCxZssQ8duwYADzJzK+EZsWxi4jOBfBPAN3gexLeBKATM2si6q6UWkREKVdffTXde++91LBhw2rH+/bbb/HPf/7TzM/PN5RSK7XWvZn5WNDrdCIwRPQQgL+fd955PGnSJKNRo0aWtj98+DDGjBlj7ty5UwEYA6AlgPsA7AcwkJmXB71IUS0iOs9/2HPu8OHDqVevXpbHmDNnDt5//30mou+01p2Y+bgLSxV+/ieEAQBeA5AC4EkiasvMQzp16qSfeeYZFR8fb2nM3NxcjB071iwpKTGZ+VZm/iSoNQYbGCJ6GsDzAwYMwKBBg4Iaa9GiRZg1axYAMIDZkL2WkCCi5kT0fePGjeNffvllIyUlxfZY+/btw2OPPWaWlpb+orVu6cSzoKhe5d4MEfUkIjV69Gi64oorbI+ntcbLL7/MK1asAID+zPyh3bGCOslLRPcBeH7YsGFBxwUAbr75ZjzzzDMAQAB2SlzcR0SpRPTd+eefHz9z5syg4gIAjRs3xhtvvGGkpKQ0UEr94H+WFS5i5v1EdIiI1NSpU4OKCwAopfDkk09Snz59AOADIupueyy7G/pfJZrVr18//P73v7c7zGm6du2KBx98EAAmElGWYwOLM1JKfZqampr4yiuvGNWda7EiOTkZs2bNMpRSTQDMcmRQcVZE1IeZB0yYMIEuuOACx8Z94IEH6PLLL2ciWkxEhp0xbP9EKaW+SE9Px5AhQ+wOcVY333wzWrZsqZVSXzg+uDiJiO7WWneZOHGi4fE4+5aolJQUPPbYYwrAUCJy5AJXcToiSiCied27d+cOHTo4Pv6oUaNUfHx8HSJ6z872tgJDRN211m0nTZpkq2qBmDhxosHMTYjoD27NEeuUUtN79OjBTj7rVXXNNdegdevWrJR6x5UJBAC8FB8fn/D444+7cigaHx+P7Oxsxcz/TUTWXr2BzcAopf5ywQUXmOedd56dzQOSlpaGDh06sFJqvGuTxDAi6q21rvvggw+6eo5k+PDhSmt9of+9NcJhSqk/9u7dWzm9B1rVFVdcgdTUVBPAi1a3tRwYIkrSWncaMmSIa3svlR544AGltb4ggHf8CouUUpOzsrJ0cnKyq/NkZmaifv36tn44RfWI6BZmTr777rtdn6tPnz4GEd1pdTs7ezC3GIbBl112mY1NrcnIyEBiYqIJwPJfTNToohtuuCEkV9N37tzZMAzjulDMFWPuOvfcc3ViYqLrE91+++1g5kQisnTYYucH7Mb69esHemEVlixZgunTp5/8/u2338Zf/vKXgCfzH4b1sLJAUT0iStVaJ3Tp0qXax2mtceedv7Y9Ly8Pffr0wZ49ewAAZWVluOeee1BaWlrtOJ06dYLW+tzgVy6qUkpd0bp164COJHbv3o377rsPWvt+dU3TxJgxY/DFF4G9jpKYmIi4uDgTvqu4A1+jlQcDgFKqY8uWLQM+4Lv++uuxbt06lJSUYO3atfj222/xyCOPBDxfVlaWYRhGO6vrjHVEdDkR7SeilUQ0kYhuIaKm/vel9DIMQ6elpVU7xqkvW3/yySdo0KDByaCsWLECnTt3Rk2HWR07dgQzG0RU/fvVxWmIaC0RfU9ErxPRvUTUgYji/P+5afv27QMaJyMjA82aNcO6desA+J7omzVrhh49An/ubtSoEQO4xsr6LZ8ZIqLU1NTUmh/ol5iYiGuuuQZvv/021q9fj4kTJyIhISHg7evVqwdmbkxE7t+4JvowfNdzdYHvJkMJ8F2qf8AwAjuFlpCQgPLycjAzvvvuO3Tp0gVlZWUAgE8//RQjRoyocYwqb1c/KO+7s60NfKcKNIAkIsoloriarjGq6tZbb8XChQvh9Xrx/fffY9KkSZYWkJKSYgCo/lnpFHZOPVv+Re/ZsyeGDx+O0aNHw+orT8wMItIAfrE6b4xrBN//KwJQ4f8zA0AxgEMI8P9jcnIyysrK8M0336B79+7wer0oLS3Fjh07YBgGLL7EXYRf780sAlP50rAJ3x0hE+G7/ex+AFlWLvXp2LEjZs+ejbfeegtTpkyB1Vee/HNZ+v23HBhmPnL4sLUr8t977z2kpaXBNE2r0+GXX34BgH3MfL7ljWOY/06A8wBshe8+rOsBbGbmMiK6xev1LghknMrALFu2DKNGjcLy5ctRVlaGTz75BDfeeGNAaykvP3mPsdbMXGz9bxO7iGgxgDoAVgBYC2A9M+8DAMMwygoLCy2d4W3Tpg1atWqF+vXrW17L0aNHNXxPUAGzfA5Ga71u586d3kAfv2DBApw4cQLZ2dn46KOPat7gFDk5OaZpmpssbxjjmHkbM1/CzHcy89+ZeQ0zl/n/82daa1VUVFTjOMnJydi8eTMaNmyIBg0aICkpCcXFxdiwYQO6dw/sEpW1a9eCiLwSF+uY+SZmvo6ZxzHz4sq4+O3asmWLpfF2796Nli1b2lrLwYMHCYClOxvYeRXpX4cOHQroAH7z5s1YtmwZRo4ciXbt2qGsrAx5eXmWJtu/fz8AyCUDDmLmMqVU+erVNd9XKCkpCQsXLkTv3r1Pfr906VJcddVVAZ9LW79+PZRSe4NatDiN1nptbm5uwE/2ALBr1y40b97c8lylpaXwer0Kv96qMyB2ArNYa42afjgPHDiAadOm4emnnz75KkOfPn2wcOHCgCfKycnB8ePHDQDyVnOHMfOmpUuX1njMmpycDGZG5XUuSUlJ2LNnT8CHRwCwbt060zTNpfZXK87iHwcOHPCUlAR204HCwkLUqVMHSUlJlieaN28elFLHmLnm3d4qbN0PRim1omnTplfPmDHD1XfzZmdn6x9++OEH0zQvcXOeWERE3QCsnDNnDmp6uToYW7duxahRowAg3eoPp6iZYRiHevXqlfbwww+7Os+AAQPM0tLSmcwc+HtMYPNaJGYe+dNPPxlWD3esKCwsxPfff6+01s+4NkkMY+ZVSqmDf/vb31x9+X/GjBmmUmqTxMUdWuu/f/7557qioqLmB9u0fPlylJaWKgCWfxftBmajUmr1mDFjzMp3Bjpt1KhRJhHlMLP1M8MiIFrrIWvWrKGtW7e6Mv6iRYuwe/dupbXu78oEAgDGer3ekueee86VX8TS0lK8+uqrGsD/MPNRq9vbvhZFa33j0aNH9bRp0xx/BpwzZw727dtHWuvrnB5b/IqZFxPRJxMnTtRVXkp2RGFhIV577TUG8AIz5zo6uDiJmU1mvmnjxo1q1apVjo8/duxYrbX+BcAwO9vbDgwzH2XmO5YtW4Z33nHuHOy//vUvzJs3DwCGMbO88uAyZr7l+PHjxcOGDTOd2s0uLi7G8OHDNYD/MLMc4rqMmVcBmPbiiy/ypk2bHBt3ypQpnJOTw/5P+7C1I+HETb/vBfB6z549MWLECArmtotvvvkmPvzwQwB4hpmnBLUwETAiaqiUyq1bt27qSy+9ZJx7rv3rEnNycjBq1Ch94sSJnVrrC5nZ+rsrhS1KqbcADBo+fDhZeZXvVBUVFRg3bpzetm0bAPRg5hV2x3LqY0tuI6J5qampxsSJEw2rd0jbu3cvnn32WfPgwYNg5geZ+Y2gFyUsIaJUpdRaZr5w4MCBNHDgQEvba60xbdo0XrZsGSmlPtda32D3WU/YR0TPARjVqlUrPWHCBKNu3bqWtl+3bh1eeOEFXVFRUcLMPZh5fVDrcfCD1+oopZZqrbs2a9bM/MMf/mB07dq12m02bNiAN99809y5c6ehlNqqtf4vZi50ZEHCFiJ6EsCUpKQk9OrVyxg0aBCqu9/IoUOHMHv2bF65ciWbpnmCme8O5mMuRPCIqI3/QxAbX3LJJXzfffdV+yGIWmssWbIEH3zwgVlUVGQQ0fvMfJcTe59ufHTsZUQ0lZmvUUqhXr16ZosWLYy0tDQiIhw5cgT5+fne4uJiZZomEdE6/8dVOn+GSthCREkAnlNK3a+1Tk1KSvI2a9ZMNWnSRHk8Hhw/fhw//fSTuXfvXhw/ftxQShVqrV8G8Gc5JAofRNRXKTVJa93G4/Ho9PR03bJlS09SUhKYGUVFRZyfn28ePXrUAGAy88cAHmPmAsfW4NZerP+eFb0B3KSUuoKI6vrnK9FarwewBMCiKtfHiDBERJcA6AfgWsMwMuD7LOQK0zRz4bsu5X+ZOb8WlyhqQET1AdwOoJdSqh0RJQMwtdZFzPw1gI8ALHfjkNa1wAghREjuySqEiE0SGCGEayQwQgjXSGCEEK6RwAghXCOBEUK4RgIjhHCNBEYI4RoJjBDCNRIYIYRrJDBCCNf8fwuLBIYQA4gtAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 266.457x266.457 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig_3_8 = DAG(\n",
    "    ebunch=[\n",
    "        (\"A\", \"X\"),\n",
    "        (\"B\", \"A\"),\n",
    "        (\"B\", \"Z\"),\n",
    "        (\"C\", \"D\"),\n",
    "        (\"C\", \"Z\"),\n",
    "        (\"D\", \"Y\"),\n",
    "        (\"W\", \"Y\"),\n",
    "        (\"X\", \"W\"),\n",
    "        (\"Z\", \"X\"),\n",
    "        (\"Z\", \"Y\"),\n",
    "    ]\n",
    ")\n",
    "\n",
    "fig_3_8_node_pos = {\n",
    "    \"A\": (0, 2),\n",
    "    \"B\": (0, 4),\n",
    "    \"C\": (4, 4),\n",
    "    \"D\": (4, 2),\n",
    "    \"W\": (2, 0),\n",
    "    \"X\": (0, 0),\n",
    "    \"Y\": (4, 0),\n",
    "    \"Z\": (2, 2),\n",
    "}\n",
    "\n",
    "fig_3_8.to_daft(node_pos=fig_3_8_node_pos).render();"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "#### Study question 3.3.1\n",
    "\n",
    "Consider the graph in Figure 3.8:\n",
    "\n",
    "(a) List all of the sets of variables that satisfy the backdoor criterion to determine the causal effect of X on Y.\n",
    "\n",
    "This question can be [answered](https://github.com/behrman/cis/blob/master/solutions_dagitty.md#study-question-331) by using the R function `dagitty::adjustmentSets()` to generate all sets of variables that satisfy the backdoor criterion. There is no comparable functionality in pgmpy as of version 0.1.22. (`get_all_backdoor_adjustment_sets()` below returns *minimal* sets)."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "(b) List all of the minimal sets of variables that satisfy the backdoor criterion to determine the causal effect of X on Y (i.e., any set of variables such that, if you removed any one of the variables from the set, it would no longer meet the criterion)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-04-14T14:22:42.905957Z",
     "start_time": "2023-04-14T14:22:42.882219Z"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "inference = CausalInference(BayesianNetwork(fig_3_8))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-04-14T14:24:10.776547Z",
     "start_time": "2023-04-14T14:24:10.105357Z"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "frozenset({frozenset({'C', 'Z'}),\n",
       "           frozenset({'A', 'Z'}),\n",
       "           frozenset({'D', 'Z'}),\n",
       "           frozenset({'B', 'Z'})})"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inference.get_all_backdoor_adjustment_sets(\"X\", \"Y\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: As of pgmpy version 0.1.22, the computation for `get_all_backdoor_adjustment_sets()` grows exponentially with the number of nodes in the graph, so it is only suitable for relatively small graphs."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "(c) List all minimal sets of variables that need be measured in order to identify the effect of D on Y."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-04-14T14:31:47.772302Z",
     "start_time": "2023-04-14T14:31:47.771464Z"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "frozenset({frozenset({'A', 'Z'}),\n",
       "           frozenset({'B', 'Z'}),\n",
       "           frozenset({'W', 'Z'}),\n",
       "           frozenset({'X', 'Z'}),\n",
       "           frozenset({'C'})})"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inference.get_all_backdoor_adjustment_sets(\"D\", \"Y\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Repeat, for the effect of {D, W} on Y.\n",
    "\n",
    "This question can be [answered](https://github.com/behrman/cis/blob/master/solutions_dagitty.md#study-question-331) by using the R function `dagitty::adjustmentSets()`, which allows the exposure and outcome to be sets of variables. There is no comparable functionality in pgmpy as of version 0.1.22."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "### 3.4 The Front-Door Criterion\n",
    "\n",
    "#### Study question 3.4.1\n",
    "\n",
    "Assume that in Figure 3.8, only X, Y, and one additional variable can be measured. Which variable would allow the identification of the effect of X on Y? What would that effect be?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "frozenset({frozenset({'W'})})"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inference.get_all_frontdoor_adjustment_sets(\"X\", \"Y\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "We saw in question 3.3.1 (b) that all backdoor adjustment sets for the effect of X on Y in Figure 3.8 require at least two variables. The single variable W, however, is a front-door adjustment set through which the effect of X on Y could be identified.\n",
    "\n",
    "Note: As of pgmpy version 0.1.22, the computation for `get_all_frontdoor_adjustment_sets()` grows exponentially with the number of nodes in the graph, so it is only suitable for relatively small graphs."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "### 3.5 Conditional Interventions and Covariate-Specific Effects\n",
    "\n",
    "#### Study question 3.5.1\n",
    "\n",
    "Consider the causal model of Figure 3.8.\n",
    "\n",
    "(a) Find an expression for the c-specific effect of X on Y.\n",
    "(b) Identify a set of four variables that need to be measured in order to estimate the z-specific effect of X on Y.\n",
    "\n",
    "These questions can be [answered](https://github.com/behrman/cis/blob/master/solutions_dagitty.md#study-question-351) by using the R functions `dagitty::adjustedNodes()` and `dagitty::adjustmentSets()`. There is no comparable functionality in pgmpy as of version 0.1.22."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "### 3.8 Causal Inference in Linear Systems\n",
    "\n",
    "Figure 3.18. Graph corresponding to Model 3.1 in study question 3.8.1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARgAAAEYCAYAAACHjumMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAAuQ0lEQVR4nO3de3hU1b3/8c937UlCQkIggEQgFTAiXpAKkVtBKwV+oOIFimiLQhW1gPDg5RhEboIgqBUryuV4KRYUwSMcBAqmWIpQQJCbYC0xhBwRucaQkExCmL2+vz8yQQQymT2zdzKX7+t5fJqE2WstSuY9e++Z2UPMDCGEcIKq7QUIISKXBEYI4RgJjBDCMRIYIYRjJDBCCMdIYIQQjpHACCEcI4ERQjhGAiOEcIwERgjhGJeTgxNRQwADAPwaQAoAAlAI4HMAHzPzESfnF0JUIKJ0AL8F0AlAXQAmgGMA/g7gf5m5xJF57X4vEhElAZihlBqita4bExNjNmrUCElJSQoASkpK9IkTJ1BeXm4QURkzLwHwX8x8wtaFCBHlvFF5lYj6MrOrTp06ZuPGjZGQkKC01nzq1Cn9448/KtM0lVKqQGs9G8ALzHzWtjXYFRgiIgDzAQxLSEjQffv2Ne69914kJCRc8vZlZWVYvnw5Vq1aZRYVFRlE9DEzD2Jm05YFCRGliChJKZWlte7cuHFjc+DAgUavXr3gcl36gKWgoAALFy7EP//5T+3xeJiZxzPzDFvWYkdgiKitUmq9UqrB6NGj1a233mpp++3bt2PmzJm6vLy8hJlvZ+aNQS9KiChERL8jogX16tVTU6ZMMVq1amVp+6VLl2LRokVMRAe01t2Y+VhQ6wk2MER0K4C/X3311TRlyhRV1R5LdcrLyzF9+nS9Y8cOAjCQmT8OamFCRBkiGgdgWt++ffmPf/wjKRXYczgnT57EM888Y548ebKcmW9g5pyA1xRMYIioA4Bt3bp1o8zMTAp4oPPMnz+fV61aBQC3MfNaO8YUItIR0WgAfx4+fDhuu+22oMfTWmPMmDFmXl5eOTO3DHRPJuDAEJFBRAXXX3993enTp9v6dPcrr7zCn3/++VlmTnHq7LYQkcJ7Mjf7wQcfpIEDB9o2rtYaw4YNM/Pz83NN02wdyBgBh4GIFsfGxtadPHmy7a+lefLJJykxMdFQSq2xe2whIo1S6h8tWrTQdsbFOy5eeuklg5nTiWhsQGMEshERtWDmgZmZmSo2NjaQIXxSSmHKlCmG1ro7EXW1fQIhIoT30Kj5Cy+8YDgxfqNGjTB48GAioheIKMbq9oHufcxq2LCh56abbgpw8+qlp6fjiiuuMInoFccmESLMKaWe7dKlC5KTkx2b495774XL5SIAz1nd1nJgqMId9957r6OvAgaAoUOHGszcmYjqOj2XEOGGiG7UWqcOGzbMlidYfLnllluUUmqU1e0CicRvABh9+vSp9oavvfYacnNzz31fUFAAj8eDxYsX+zVRRkYGXC4Xezye3wF4K4C1ChHJhjdo0MDTqFEjn/fjYO+HAPDggw9i3bp1KURU18oTL4EEpl9SUpKplKp22zFjxpz7+ujRoxg7diwef/xxS5Nddtll+ocffugFCYwQP6OU6nLllVfWyP2wQYMGMAxDm6Z5F4AP/F6jpVkAEFGXFi1aWDqhVFhYiMmTJ2PQoEHo1KmTpflat27tUkplWNpIiOjQqm3btn7fOJj7IQCkpKRoAL2tbGM5MEqplIYNG/p9zFdWVoapU6eiW7du6Nu3r9Xp0KhRo8o3UAohzsPMsc2aNfPrtsHeDwEgOTlZAbjMyjaBHCIZ/r4E2TRNvPzyy2jevDkGDx587udHjx7F0qVLUVJSgmeffdb3ZIYBAI48BSdEmCPv/cOnqu6HW7ZswZdffgm3241evXqhffv2PsdRShEsNsPyHgwzu91ut1+3nTt3LjweD0aN+vnJ59TUVIwePdqvMYqLiwHgjMVlChENzFOnTlV7o6ruh126dMGoUaMwcuRIbNq0qdpxSktLNYDTVhZoOTBa630HDx6s9pIKixcvxoEDBzB27Fj4U9mq5OTkaNM09wc8gBARSil1fN++fT5v48/9cMmSJX69f+nEiRMAsM3KGgM5RMo6efLkb33d4NixY/jggw9w2WWXITMz89zPmzVr9rPv/XHo0CENYEMA6xQiopmmuTM7O/tyVHEKobr7ITPjvffeQ4cOHZCenu5zLo/Hg7KyMgOApascBBKYZR6P5+2cnJwqF9WkSROsXLmyygGKioqwcOFC5Obm4qOPPkJV76HIz8+H2+12AfD/CXshosfyw4cP9/N4PJe8mFR198OVK1di9+7dKCkpwZEjR3ye/M3KygIRebTWli7dENC7qQ3D+M+111571YsvvujoRcNfeukl/te//nXcNM1UJ+cRIhx5X1VfNmTIkNgBAwY4Otcf/vAH8+TJk6uZ+S4r2wUUCK31c19//TUVFRUFsrlfysvLsXnzZtZav+TYJEKEMa7wP8uWLXP0MrO5ubk4efKkAeAJq9sGFBhm/piIjk6cONGxv9j06dO11roYwCyn5hAiAow+ffo0FixY4MjgWmtMnDjRJKKtzJxb/RY/F/Ahjta654EDB9TatfZfdG7btm3YsWOHYuY72O6PPRAigjBzPjOPWrZsGX///fe2jz9v3jwuKirSzGzpFbyVAg4MM/8bwPQ5c+bwrl27Ah3mIjk5OZg2bRoT0QK5+LcQ1WPmuUS0fcyYMbqgoMC2cVevXo01a9YQMw9lZkuvf6kU1ElaZh4PYNGkSZPYnxfqVGfPnj146qmnmJmJmT8JekAhooTWem15ebl65JFHtB17MkuXLsW8efMA4Flm9vvNjRey62NLXgcw6sYbb9Tjx4+3fJU7j8eDV199lTdu3EhEtJiZGwHoBaA/My8PeoFCRDAimgxgEoApSqk+zHxT//79aejQoZbHKiwsxPjx4828vDwFYAwzvx7U2mz84LVuRLTaMIzEnj17qiFDhiAxMdHnNmVlZXj//fexZs0aXV5efoaZBzLzau94WZDICOHTeXF5npkne382nIhmJyYmYsCAAcZdd91V5YeuVTpx4gTefvtt3rp1KwB8r7XuEczHlZxbn53nUL2f7viCUmqE1rp+48aNPenp6a4bbrgBjRs3BhEhPz8fX331FXJycjzHjh0ziMittf4LgKeZ+cwF40lkhKjCpeJy3p81BPBnIhpIRDFNmzbVV199tXH99dcjOTkZpmnihx9+wFdffcW5ublmQUGBSyl1RGs9nZnfsG2NTj1JQ0Q9AAxVSnUC8AtmjvP+/AyAw1rrbQAWM3PVLzWEREaIS/EVlwtuRwCGALjbMIz2zNyEmWMBaCIqA3BAa70FwBxm3mP7OmvqWWAiOg6gMTNbvn6oREaIn/gbFx/bM4D3mHmovSu7mKMv9beL9zn4vwNYRkT31PZ6hKgtwcalpoVFYACJjBDhFhcgjAIDSGRE9ArHuABhFhhAIiOiT7jGBQjDwAASGRE9wjkuQJgGBpDIiMgX7nEBwjgwgERGRK5IiAsQ5oEBJDIi8kRKXIAICAwgkRGRI5LiAkRIYACJjAh/kRYXIIICA0hkRPiKxLgAERYYQCIjwk+kxgWIwMAAEhkRPiI5LkCEBgaQyIjQF+lxASI4MIBERoSuaIgLEOGBASQyIvRES1yAKAgMIJERoSOa4gJESWAAiYyofdEWFyCKAgNIZETtica4AFEWGEAiI2petMYFiMLAABIZUXOiOS5AlAYGkMgI50V7XIAoDgwgkRHOkbhUiOrAABIZYT+Jy0+iPjCAREbYR+LycxIYL4mMCJbE5WISmPNIZESgJC6XJoG5gERGWCVxqZoE5hIkMsJfEhffJDBVkMiI6khcqieB8UEiI6oicfGPBKYaEhlxIYmL/yQwfpDIiEoSF2skMH6SyAiJi3USGAskMtFL4hIYCYxFEpnoI3EJnAQmABKZ6CFxCY4EJkASmcgncQmeBCYIEpnIJXGxhwQmSBKZyCNxsY8ExgYSmcghcbGXBMYmEpnwJ3GxnwTGRhKZ8CVxcYYExmYSmfAjcXGOBMYBEpnwIXFxlgTGIRKZ0CdxcZ4ExkESmdAlcakZEhiHSWRCj8Sl5khgaoBEJnRIXGqWBKaGSGRqn8Sl5klgapBEpvZIXGqHBKaGSWRqnsSl9khgaoFEpuZIXGqXBKaWSGScJ3GpfRKYWiSRcY7EJTRIYGqZRMZ+EpfQIYEJARIZ+0hcQosEJkRIZIIncQk9EpgQIpEJnMQlNElgQoxExjqJS+iSwIQgiYz/JC6hTQIToiQy1ZO4hD4JTAiTyFRN4hIeJDAhTiJzMYlL+JDAhAGJzE8kLuFFAhMmJDISl3AkgQkj0RwZiUt4ksCEmWiMjMQlfElgwlA0RUbiEt4kMGEqGiIjcQl/EpgwFsmRkbhEBglMmIvEyEhcIocEJgJEUmQkLpFFAhMhIiEyEpfII4GJIOEcGYlLZKqRwBDRZQAae79uWhNzRqtwjIzEpeYQUUvvly2JKNnp+WpqD+bT877+oobmjFrhFBmJS80hojoA9nu/vRnAfzs9Z00F5tB5Xx+poTmjWjhERuJSs5i5DECR91sTQLbTc9ZUYDYAKAfAAD6voTmjXihHRuJSa3Z5/7cEwHanJ6upwOwAUAagGHKIVKNCMTISl1q1AcBZALGouF86qqYCswtAPAAXgC9raE7hFUqRkbjUuu2oODw6A+AHpycjZnZ6joqJiH4AkAwgkWtqUvEzRJQFoBeAeQAyAPw/Zv6xBuZtA2AxgAMABkDiUmuIqBGAEwA2MvPNjs9XE/d1ImoAYBmAhgBuZeZ8xycVl0RE+QBSULGbvJKZB/i5XQyAKwDUBXAawEF/HiiIyEDFHuz1AAjAx8z82wCXL4JERJcD2AdgCYBMZj7t5HyOHCIRUUMiesflcn2nlDoL4EcAvwbQFsBJpZTH5XIdJqL3iaiZE2sQVTK8/xsDoA8R3XWpG1GFRwzD2G0YRgkqTtJ/C2A3KvZEtGEYhUqpLUR0p4/5ngbQChVxYQBX2vT3EH4gouuIaKXL5TpBRCYqDotSAAwHUKSUOmMYxgEiepmI4m2f3849GCJqrZRaorX+ZVJSktmxY0cjIyMDGRkZqFOnDgCgvLwcu3btwvbt27Ft2zazoKDAUErt11r/npkdP+kU7YioKyoevRqi4rzYKQBXVh4qEREBmEtEfyCimLZt26J9+/bUtWtXpKamnhsnLy8PW7ZswY4dO3R2djYRUYnWegYzTztvrjYAdnrnKUFFmO5l5v0QjiKifkqpuVrrZqmpqWbHjh2Njh07om3btlCqYr+iuLgYW7duxY4dO7B9+3bzzJkzioj+ycwD7TrKsC0wRDQDwDNpaWk8cuRIdd111/m1XU5ODubMmaO//fZbAvAOgEflHI2zvI9UMwEMAxAH4H+ZeQARdVZKrTEMo959992n+vfvD5fLVe14brcbixYtwurVqxnAd1rrX6PitU+Vh0ZlACYD+BMzm878rQQAEFEcEX3CzL0zMjL4j3/8IzVp0sSvbTdt2oS33nrLLCgoYGYexczzgl5PsPdlIjKUUl8BuObRRx+l22+/PaBxPv/8c8yaNUtrrX/QWrdh5pKgFiaq5d2bWQqgGYC/APjDDTfcoCdMmKAq9zitKCgowPjx483vvvuOACwC8CCAryB7LTWCiK4gor3x8fEJEydONPx9kL/QX/7yFyxfvpwBrNZa9wtqTcEExhuX/bGxsS3nzp2rGjVqFMxaUFxcjMcee8wsLi7+UWudxsxnghpQVMu7N/N3AL8aNmwY7rrrkqdkLJkzZw6vWbOGALwH4GHZa3EeETUjopy0tLSYP//5z4Y/e56+ZGdn45lnnmGt9Wda614BryuYwBiG8Y+YmJhb3n77bVW/fv2Axzmf2+3Gww8/bLrd7n+bpnmDLYOKKhFRbwCfPvroo+jXL6gHq5/xRoZRcX4nz7aBxUWIiJRSJ5s2bZr85ptvGpXnWIKVk5ODp556irXWrzPzmEDGCHglRDRQa33rjBkzbIsLACQkJGDWrFkGM19PRGNsG1hchIhiiGhZly5d2M64AMCIESMoLS2NlVLrbR1YXMoCwzDqz5o1y7a4AEB6ejpGjRpFAEYTUUDHWwGthogMIlrYs2dPTk9PD2QIn1JTU3HvvfcSEb1SE28pj1ZEtDguLi7+mWeeISfGnzZtmgHgCiJ6zonxBUBE7QA8+PTTTwd03qw6PXv2ROvWrbVSal0g2weau0zDMGK8dXPE4MGDkZCQAAAvOjVHNCOiOGa+Z8SIESrY4/WqNGjQALfffjsppTIdmUCAiN5IS0szu3bt6tgckyZNMrTWqUTUw+q2AQVGKfVE9+7dlZ27Y5fSt29fQyk1xNFJotfE2NhYvvXWWx2dZPDgwWDmJCKy9xhMgIjqMvOvHnzwQaP6WweuXr16SE9P10qpl61ua7kQRNRCa93ooYceuuSfa61x3333nfs+NzcX/fr1w+HDhwEApaWlGDJkCNxud7Vz3X///WDmeCJy9l4QhZRSD3Xv3r3aX8zXXnsNo0ePPvffAw88gPvvv9/veRISEnDVVVdpIhoX1ILFpYyJi4vTnTt3vuQf2nlffPDBB5XW+kbvCzH9Fsi+8cC4uDizfv36l/zlvHCvZu3atUhJSTn3l9iwYQM6duxYefjjU2xsLOrVq2cWFhbeDUBOFtpIa33ZLbfcUu3txowZc+7ro0ePYuzYsXj88cctzdW+fXuVk5NzrdU1imr1vOKKK6q8w9t5X2zXrh1Q8XaPTgC2+rvAQI5xbm3a1PdldePi4lBWVobS0lJ8/fXX6NKlC0pLSwEAn376Kfr27ev3ZC1btjSIqFsA6xRVIKLWAFTbtm393qawsBCTJ0/GoEGD0KlTJ0vzde3aFVrrelYf/YRvhmFcf8011/i8D9t1X1RKITEx0QPgbitrtBwYpVR68+bNfe5aJyQkoLS0FBs2bED37t2RmJgIt9uNAwcOwDAMtGrVyu/5WrRoQUqpNKvrFD51jYmJMf09uVtWVoapU6eiW7dulh4cKrVsWXmdaVxheWNRJa11cps2bXzexs77YuPGjRUAS09XB3IOJjYuLs7nbSr/UuvWrUPv3r0RHx+P0tJSrF27Fn369MHRo0fx+uuv48UXq3+CyDuXM09zRK9Ef0/Qm6aJl19+Gc2bN8fgwYPP/fzQoUN488038eKLL+Jvf/ubv/MmWV+q8KHap6aruy9u2bIFs2fPxsyZM7Fz506fY8XGxhIq3rvm/wKt3BgAmLn87NmzPm+TkJCAPXv2oGHDhkhJSUF8fDwKCgqwc+dOdO/eHampqRg9erRf83k8HgDwWF2n8KlUa+3XDefOnQuPx4NRo0b97OdpaWkYOXIkMjMz8c033/g7r7y/zF5cXl7u8wbV3Re7dOmCUaNGYeTIkdi0aZPPsc6ePcuouGyH3ywHRmt9+MiRIz5/O+Pj47FixYpzu9Px8fHIyspC586dUd3ez4W+++47MPNxq+sUPu0+e/asUV1kFi9ejAMHDmDs2LEwjIuPir/44gs8//zz6NChg89xjh49WvnlIV+3E9YopYoPHDjg8zb+3heXLFmC2267zedY+fn5GhWX3PB/jVZu7LXh+++/9/mbmZCQAGauPPOM+Ph4HD58GH369LE8WW5urkdrvTmAdYqq7QTA+/dX/QbnY8eO4YMPPkBRUREyMzPPPU09c+bMc7fp1KkTnn/+eWzYsMHnZJs3b4ZSqoSZfe/6Cku01t/s27fP55sJq7svMjMWLFiADh06oLpX5RcVFRkAVlpZYyDnNj52u92TysrKUNXx3/lPbQJA586dsXLlT+sqKirCwoULkZubi48++ggDBw685DhaaxQUFBgA/D7IF9VjZjYMo2Djxo0p11xzzSVv06RJk5/9m11o79692Lx5M86ePYuMjAyf8+3cuZMB5ASzZnExZv48Ly/vJvi4H1d3X1y5ciV2796NkpISHDlypMqT+Dk5OWBmgsWXiwT0bmrDME7369cvcdiwYZa3tWLFihV45513zjJznFyEyl5ENL9u3boPf/jhh46+CtTj8WDAgAGstX6Emd9xcq5oQ0SNARyfOXMmrr3W2ZcZPffcc3rfvn25pmleZWW7gF7rr7V+Oysry/FrfHz88cceZl4ucXHEMyUlJWr37t2OTrJs2TIwczmAdx2dKAox8wml1N53333XvzP2ASovL8fevXtJaz3J6raBvplofGlpKS1evDjAzauXlZVVeXj0pGOTRDFmLlRKbX711Vcde6AoKyvDhx9+qJn5XXmQcIbW+sn9+/crX+fTgvXKK69oIipm5g+sbhtQYLyXs8xcvHgxHzt2LJAhfCouLsacOXMYwJ+Z+bDtEwgAgNb69lOnTvFbb73lyPhTpkzRpmkWARjpyAQCzLyOiNZPmjTJ9PelB1bs3r0bW7ZsUVrrgD5qJuC3QzPzK0S0/+mnnzarey7eCq01nnzySZOZj0H2XhzFzIXM/Ngnn3yC7dvt/ZjiFStWVO5W95W9F2cx8+1ut/vs5MmTbS1Mfn4+pkyZor0XEc8KZIygrregte5YVFR0+rHHHjPLysqCGQoAKl/QZR49erRca91efjGdx8zvAtg5depU2yKzcuVKvP322wDwHjP7/cY4ERhmLmXm7rt27cK0adNsiczJkycxfPhwbZrmAWa+O9Bx7PhUgYZKqf1xcXENJkyYYOkNdOfLzc3F+PHjzeLi4jJmvo6Z/y+ohQm/nPdZ0TsBtL/zzjvxyCOPBDSWx+PBzJkz9datWxWAbACtAfRn5uV2rVdUjYi6EtGGyy+/nKZPn240bNgwoHHWr1+P119/XWuts7XWNwTz+iVbPhfJe23XFczct3Pnzvz4449TcrJ/V7p0u92YN28er1+/npRSm7XWfZz+OEtR4cIPoieih4hofv369enJJ580fvnLX/o91vr16zFnzhx95syZUma+23tuoPKzsCUyNYSI0pRS65m5Vf/+/el3v/sdYmNj/dr28OHDmDVrlt6/fz+h4vPLRwZ7FGH3JzverpRaoLVu1KJFC7N3795G165dcWFJCwsLsXnzZqxbt05nZ2crpVSR1noEM79v22KETxfG5byfJyulVmmtf1W3bl3do0cPo1u3bmjTps3Pri/i8Xjw1VdfYcOGDdi0aZNZXl5ORLScme8//xFPIlM7iOg5IhoPIK5du3bco0cP1alTp4uu/XL48GH861//QlZWlnns2DFDKfV/Wuu7mHmPLQthZtv/A9CNiDYrpdwAWCmlY2NjzdjYWNMwDNP7szKl1C4AfZxYg/zn899nMio+J3qyj9skA5inlMoHoAFwTEyMJy4uzuNyuUzv9qZS6giAFwDE+Bgry3v7e2r77x5t/wEYopTaT0TlANgwDNN7X/QQkfbeF4uJaA2ANnbPb+sezKV4P9irFyo+6S8JwAAAn7J8cmOtqGrPxY/tbgRwI376nOnNzJxtYXvZk6ll3lf+/gbAYgCbATwP4DN28IPxHA/MuYmIjgNozBXvZxC1INC42Di/RCYEEBGj4hm+oU7P5ezHAoiQUdtxAQBm7o2Kj6ldRkT31MYaRM2SwESBUIhLJYlMdJHARLhQiksliUz0kMBEsFCMSyWJTHSQwESoUI5LJYlM5JPARKBwiEsliUxkk8BEmHCKSyWJTOSSwESQcIxLJYlMZJLARIhwjksliUzkkcBEgEiISyWJTGSRwIS5SIpLJYlM5JDAhLFIjEsliUxkkMCEqUiOSyWJTPiTwIShaIhLJYlMeJPAhJloiksliUz4ksCEkWiMSyWJTHiSwISJaI5LJYlM+JHAhAGJy08kMuFFAhPiJC4Xk8iEDwlMCJO4VE0iEx4kMCFK4lI9iUzok8CEIImL/yQyoU0CE2IkLtZJZEKXBCaESFwCJ5EJTRKYECFxCZ5EJvRIYEKAxMU+EpnQIoGpZRIX+0lkQocEphZJXJwjkQkNEphaInFxnkSm9klgaoHEpeZIZGqXBKaGSVxqnkSm9khgapDEpfZIZGqHBKaGSFxqn0Sm5klgaoDEJXRIZGqWBMZhEpfQI5GpORIYB0lcQpdEpmZIYBwicQl9EhnnSWAcIHEJHxIZZ0lgbCZxCT8SGedIYGwkcQlfEhlnSGBsInEJfxIZ+0lgbCBxiRwSGXtJYIIkcYk8Ehn7SGCCIHGJXBIZe0hgAiRxiXwSmeBJYAIgcYkeEpngSGAskrhEH4lM4CQwFkhcopdEJjASGD9JXIRExjoJjB8kLqKSRMYaCUw1JC7iQhIZ/0lgfJC4iKpIZPwjgamCxEVURyJTPQnMJUhchL8kMr5JYC4gcRFWSWSqJoE5j8RFBEoic2kSGC+JiwiWROZiEhhIXIR9JDI/F/WBkbgIu0lkfhLVgZG4CKdIZCpEbWAkLsJpEpkoDYzERdSUaI9M1AVG4iJqWjRHJqoCI3ERtSVaIxM1gZG4iNoWjZGJisBIXESoiLbIRHxgJC4i1ERTZCI6MBIXEaqiJTIRGxiJiwh10RCZiAyMxEWEi0iPTMQFRuIiwk0kRyaiAiNxEeEqUiMTMYGRuIhwF4mRiYjASFxEpIi0yIR9YCQuItJEUmRqJDBENBdAY+/XS20cdzIkLiICOREZIoolos+93w4hokw7xvWlpvZg2p/3dQc7BpS4iEjnQGQ8ADqe9/1VNozpU00F5nMA2vv1F8EOJnER0cLOyDCzBrDf+20JgM1BLq9aNRWYLwAUA3CjIjYBk7iIaGPznsxG7/9qAF8GOVa1aiowOwC4ULGLtiPQQSQuIlrZGJktqNh7qQPg33aszRdiZqfnABERgNOo+EslMnOZhW2vB9AUQFdIXESUI6IsAL0A9AdgANjGzN9Z2L4NgG8A/IeZr3FmlT9xOT0BADAzE9G/ATS0GJdEAJ8BuMz7I4mLiGrM3NsbmWXeH+0iogzv+RV/ZAM4A+BfjizwAo4dIhFRYyJ6gohWulyu7wzDuN4wjGYul+t7IvqUiDKJKK2aYV4DkHTe946flBI/oQp3E9HbhmF87XK5Trhcrh9dLtdxpdROIppNRLfW9jqj0Cfnfd0awBO+bkxE1xLRFCLaYBjGYaPCIMMwconoIyJ6lIjqOrFQ2w+RiOhupdR0rfU1cXFxZmpqKlq3bm2kpKQAAE6dOoVvv/1W//DDD1xWVmYopXK11pOZeeEF4/QAsBJAAoAyAIUAejLzPlsXLC5CRM0AvOaNi9GgQQPzyiuvdKWlpSEmJgZlZWXIy8vjvLw8s7Cw0CCiMq31XwFkMnNhba8/0hFRT1TswdQBEAOgFEA7Zv72vNsQgKeUUv+ltb4sISHBk5aWZqSnp1NiYiKYGcePH0d2drbn+PHjyuPxEBF9wcxPMbNtD+S2BYaImiml/sHMra+99lo9bNgwlZ6e7nObQ4cO4a233tK7d+8mIvpBa92DmbO9h0YHUHFoVArgAwBjmLnYlsWKS/L+Us4CMDolJcXs37+/q1+/flCq6h3d8vJyLF68GGvWrDFLSkoA4Glmfq1mVhy9iKgJgPcAdEPFg/BuABnMrImou1JqFREl/upXv6KHHnqIGjZs6HO8L7/8Eu+9956Zl5dnKKU2aq37MnNJ0Ou0IzBE9EcAb15++eU8bdo0o3Hjxpa2LywsxIQJE8yDBw8qABMAtATwMIBjAO5n5vVBL1L4RESXew97mowYMYJ69+5teYxFixZh6dKlTERfa60zmPmMA0sVXt4HhEEA3gKQCOBpIrqemYdmZGToZ599VsXGxloaMycnBxMnTjSLi4tNZr6LmdcGtcZgA0NEYwG8OGjQIAwePDiosVatWoX58+cDAAN4F7LXUiOI6Aoi+iY1NTX21VdfNRITEwMe6+jRo3jiiSdMt9v9o9a6pR2PgsK3yr0ZIupJRGr8+PF00003BTye1hqvvvoqb9iwAQAGMvPHgY4V1EleInoYwIvDhw8POi4AcMcdd+DZZ58FAAJwUOLiPCJKIqKvf/GLX8TOmzcvqLgAQGpqKt555x0jMTExRSn1H++jrHAQMx8jolNEpGbNmhVUXABAKYWnn36a+vXrBwAfEVH3gMcKdEPvs0TzBwwYgNtuuy3QYS7StWtXPPbYYwAwlYha2zawuCSl1KdJSUl1XnvtNcPXuRYrEhISMH/+fEMp1RTAfFsGFVUion7MPGjKlCnUqlUr28Z99NFHqUOHDkxEq4nICGSMgH+jlFL/aNSoEYYOHRroEFW644470LJlS62U+oftg4tziOj3WusuU6dONVwue18SlZiYiCeeeEIBGEZEtrzBVVyMiOKIaEn37t25Xbt2to8/btw4FRsbW5eIPgxk+4ACQ0TdtdbXT5s2LaCq+WPq1KkGMzclogecmiPaKaXm9OjRg+181DvfzTffjKuuuoqVUu87MoEAgFdiY2PjnnzySUcORWNjY5GZmamY+bdEZO3ZGwQYGKXUn1q1amVefvnlgWzul+TkZLRr146VUpMdmySKEVFfrXW9xx57zNFzJCNGjFBa66u9r60RNlNK/aFv377K7j3Q8910001ISkoyAbxkdVvLgSGieK11xtChQx3be6n06KOPKq11Kz9e8SssUkpNb926tU5ISHB0nvT0dDRo0CCgX07hGxHdycwJv//97x2fq1+/fgYR3Wd1u0D2YO40DINvvPHGADa1Ji0tDXXq1DEBWP6LiWpd06tXryr//bXWuO++n/5vz83NRb9+/XD48GEAQGlpKYYMGQK3213tRB07djQMw/h18EsWF/hdkyZNdJ06dRyf6J577gEz1yEiS4ctgQSmT4MGDfx9YxXWrFmDOXPmnPt+4cKF+NOf/uT3ZN7DsB5WFih8I6IkrXVcly5dqrzNhc8orV27FikpKeeCsmHDBnTs2BH+7AFlZGRAa90kuFWLCymlbrrqqqv8OpI4dOgQHn74YWhdcdc1TRMTJkzAP/7h3/ModerUQUxMjImKd3H7v0YrNwYApVT7li1b+n3Ad+utt2L79u0oLi7Gtm3b8OWXX+Lxxx/3e77WrVsbhmG0tbrOaEdEHYjoGBFtJKKpRHQnETXzvi6lt2EYOjk52ecYcXFxKCsrQ2lpKb7++mt06dIFpaWlAIBPP/0Uffv29Wst7du3BzMbROT79eriIkS0jYi+IaK3ieghImpHRDHeP252ww03+DVOWloamjdvju3btwOoeKBv3rw5evTw/7G7cePGDOBmK+u3fGaIiJKSkpKqv6FXnTp1cPPNN2PhwoXYsWMHpk6diri4OL+3r1+/Ppg5lYicv3BN5GFUvJ+rCyouMhSHirfqHzeM6h/4EhISUFpaii+++ALdu3eHx+OB2+3GgQMHYBgG/H326byXq5+U190FrA0qThVoAPFElENEMdW9x+h8d911F1asWAGPx4NvvvkG06ZNs7SAxMREA4DvR6ULBHLq2fIdvWfPnhgxYgTGjx8Pq888MTOISAP40eq8Ua4xKv6tCEC592cGgAIAp+DHv2NlYNatW4dx48Zh/fr1KC0txdq1a9GnTx8cOnQIn3zyCYqKitCuXTt/XnCZj5+uzSz8U/nUsImKK0LWQcXlZ48BaG3lrT7t27fHu+++i7/+9a+YMWMGrD7z5J3L0v3fcmCYuaiw0No78j/88EMkJyfDNE2r0+HHH38EgKPM/AvLG0cx75UAlwDYi4rrsO4AsIeZS4noTo/Hs7y6MRISErBnzx40bNgQKSkpiI+PR0FBAXbu3Ilhw4YhLi4OI0eOhNYas2bNqjIwZWXnrjF2FTMX2PM3jA5EtBpAXQAbAGwDsIOZjwKAYRilJ06csHSGt02bNrjyyivRoEEDy2s5ffq0RsUDlN8sn4PRWm8/ePCgx9/bL1++HGfPnkVmZiY++eST6je4QHZ2tmma5m7LG0Y5Zt7HzNcx833M/CYzb2XmUu8f/11rrfLz832OER8fjxUrVpw71xIfH4+srCx07tz53GHuF198geeffx4dOlT9Yt1t27aBiDwSF+uY+XZm/jUzT2Lm1ZVx8fruq6++sjTeoUOH0LJly4DWcvLkSQJg6coGgTyL9LdTp075deZ6z549WLduHcaMGYO2bduitLQUubm5liY7duwYAMhbBmzEzKVKqbLNm31fVyghIQHMjMqXoMfHx+Pw4cPo06fPudt06tQJzz//PLzvvL2kHTt2QCl1xJ7Vi0pa6205OTl+P9gDwHfffYcrrrjC8lxutxsej0fhp0t1+iWQwKzWWqO6X87jx49j9uzZGDt27LmnMvv164cVK1b4PVF2djbOnDljAJCXmtuMmXdnZWX5PGYdM2YM5s+fj8oTs507d8bKlSuRllbxuse9e/di/vz5eOONN5CRkVHlONu3bzdN08yycfmiwl+OHz/uKi7276IDJ06cQN26dREfH295oiVLlkApVcLMvnd7LxDQ9WCUUhuaNWv2q7lz5zr6at7MzEz9n//85z+maV7n5DzRiIi6Adi4aNEiVPd0dTD27t2LcePGAUAjq7+conqGYZzq3bt38siRIx2dZ9CgQabb7Z7HzP6/xgQBvheJmcd8//33htXDHStOnDiBb775Rmmtn3VskijGzJuUUiffeOMNR5/+nzt3rqmU2i1xcYbW+s3PPvtMl5eXV3/jAK1fvx5ut1sBsHxfDDQwu5RSmydMmGBWvjLQbuPGjTOJKJuZrZ8ZFn7RWg/dunUr7d2715HxV61ahUOHDimt9UBHJhAAMNHj8RS/8MILjtwR3W43Xn/9dQ3gv5n5tNXtA74ejNa6z+nTp/Xs2bNtfwRctGgRjh49SlrrX9s9tvgJM68morVTp07V5z2VbIsTJ07grbfeYgAzmTnH1sHFOcxsMvPtu3btUps2bbJ9/IkTJ2qt9Y8AhgeyfcCBYebTzHzvunXr8P779p2D/dvf/oYlS5YAwHBmlmceHMbMd545c6Zg+PDhpl272QUFBRgxYoQG8G9mlkNchzHzJgCzX3rpJd69e7dt486YMYOzs7PZ+2kfAe1I2HHR74cAvN2zZ0+MGjWKgrns4oIFC/Dxxx8DwLPMPCOohQm/EVFDpVROvXr1kl555RWjSZPA35eYnZ2NcePG6bNnzx7UWl/NzNZfXSkCopT6K4DBI0aMoPNfSmBVeXk5Jk2apPft2wcAPZi56tcgVMOujy25m4iWJCUlGVOnTjWsXiHtyJEjeO6558yTJ0+CmR9j5neCXpSwhIiSlFLbmPnq+++/n+6//35L22utMXv2bF63bh0ppT7TWvcK9FFPBI6IXgAw7sorr9RTpkwx6tWrZ2n77du3Y+bMmbq8vLyYmXsw846g1mPjB6/VVUplaa27Nm/e3HzggQeMrl27+txm586dWLBggXnw4EFDKbVXa/0bZj5hy4JEQIjoaQAz4uPj0bt3b2Pw4MHwdb2RU6dO4d133+WNGzeyaZpnmfn3wXzMhQgeEbXxfghi6nXXXccPP/ywzw9B1FpjzZo1+Oijj8z8/HyDiJYy8+/s2Pt04qNjbySiWcx8s1IK9evXN1u0aGEkJycTEaGoqAh5eXmegoICZZomEdF278dV2n+GSgSEiOIBvKCUekRrnRQfH+9p3ry5atq0qXK5XDhz5gy+//5788iRIzhz5oyhlDqhtX4VwMtySBQ6iKi/Umqa1rqNy+XSjRo10i1btnTFx8eDmZGfn895eXnm6dOnDQAmM68E8AQz/59ta3BqL9Z7zYq+AG5XSt1ERPW88xVrrXcAWANg1XnvjxEhiIiuAzAAwC2GYaSh4rOQy03TzEHF+1L+h5nzanGJohpE1ADAPQB6K6XaElECAFNrnc/MWwB8AmC9E4e0jgVGCCHs+aQtIYS4BAmMEMIxEhghhGMkMEIIx0hghBCOkcAIIRwjgRFCOEYCI4RwjARGCOEYCYwQwjESGCGEY/4/MPDFrPEILpQAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 266.457x266.457 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig_3_18 = DAG(\n",
    "    ebunch=[\n",
    "        (\"W_1\", \"X\"),\n",
    "        (\"W_2\", \"Y\"),\n",
    "        (\"W_3\", \"Y\"),\n",
    "        (\"X\", \"W_3\"),\n",
    "        (\"Z_1\", \"W_1\"),\n",
    "        (\"Z_1\", \"Z_3\"),\n",
    "        (\"Z_2\", \"W_2\"),\n",
    "        (\"Z_2\", \"Z_3\"),\n",
    "        (\"Z_3\", \"X\"),\n",
    "        (\"Z_3\", \"Y\"),\n",
    "    ]\n",
    ")\n",
    "\n",
    "fig_3_18_node_pos = {\n",
    "    \"W_1\": (0, 2),\n",
    "    \"W_2\": (4, 2),\n",
    "    \"W_3\": (2, 0),\n",
    "    \"X\": (0, 0),\n",
    "    \"Y\": (4, 0),\n",
    "    \"Z_1\": (0, 4),\n",
    "    \"Z_2\": (4, 4),\n",
    "    \"Z_3\": (2, 2),\n",
    "}\n",
    "\n",
    "fig_3_18.to_daft(node_pos=fig_3_18_node_pos).render();"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "#### Study question 3.8.1\n",
    "\n",
    "Given the model depicted above, answer the following questions:\n",
    "\n",
    "(a) Identify three testable implications of this model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(W_1 ⟂ W_2)\n",
       "(W_1 ⟂ W_3 | X)\n",
       "(W_1 ⟂ Y | Z_3, Z_1, W_3)\n",
       "(W_1 ⟂ Z_2)\n",
       "(W_1 ⟂ Z_3 | Z_1)\n",
       "(W_2 ⟂ W_3 | Z_2)\n",
       "(W_2 ⟂ X | Z_2)\n",
       "(W_2 ⟂ Z_1)\n",
       "(W_2 ⟂ Z_3 | Z_2)\n",
       "(W_3 ⟂ Z_1 | X)\n",
       "(W_3 ⟂ Z_2 | X)\n",
       "(W_3 ⟂ Z_3 | X)\n",
       "(X ⟂ Y | W_2, Z_3, W_3)\n",
       "(X ⟂ Z_1 | W_1, Z_3)\n",
       "(X ⟂ Z_2 | W_1, Z_3)\n",
       "(Y ⟂ Z_1 | W_2, Z_3, W_3)\n",
       "(Y ⟂ Z_2 | W_2, Z_3, W_3)\n",
       "(Z_1 ⟂ Z_2)"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print_implied_conditional_independencies(fig_3_18)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "The above conditional independencies are all testable implications of this model."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "(b) Identify a testable implication assuming that only X, Y, W_3, and Z_3 are observed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(W_3 ⟂ Z_3 | X)"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig_3_18.latents = set(fig_3_18.nodes()).difference({\"W_3\", \"X\", \"Y\", \"Z_3\"})\n",
    "\n",
    "print_implied_conditional_independencies(fig_3_18)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "W_3 and Z_3 are independent conditional on X."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig_3_18.latents = set()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "(e) If we regress Z_1 on all other variables in the model, which regression coefficients will be zero?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Markov blanket of Z_1: ['W_1', 'Z_2', 'Z_3']\n"
     ]
    }
   ],
   "source": [
    "print_markov_blanket(graph=fig_3_18, node=\"Z_1\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'W_2', 'W_3', 'X', 'Y'}"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "set(fig_3_18.nodes()).difference({\"Z_1\"}.union(fig_3_18.get_markov_blanket(\"Z_1\")))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Z_1 is d-separated from the variables {W_2, W_3, X, Y} by the set {W_1, Z_2, Z_3}. This implies that in a regression of Z_1 on all other variables, the regression coefficients of {W_2, W_3, X, Y} would be zero."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "(g) Assume that variables Z_2 and W_2 cannot be measured. Find a way to estimate b using regression coefficients. \\[Hint: Find a way to turn Z_1 into an instrumental variable for b.\\]\n",
    "\n",
    "This question can be [answered](https://github.com/behrman/cis/blob/master/solutions_dagitty.md#study-question-381) by using the R function `dagitty::instrumentalVariables()`, which generates a list of instrumental variables under linearity assumptions. There is no comparable functionality in pgmpy as of version 0.1.22."
   ]
  }
 ],
 "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.11.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}