derivation_validation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Imports\n",
    "\n",
    "import warnings\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from datetime import datetime\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as scs\n",
    "import scipy.integrate as si\n",
    "import seaborn as sns\n",
    "import numpy.random as npr\n",
    "from sklearn.preprocessing import OneHotEncoder\n",
    "from sklearn.linear_model import LinearRegression, LogisticRegression\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "\n",
    "\n",
    "# Settings\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "plt.rcParams.update({'font.size': 16})\n",
    "plt.rcParams.update({'figure.figsize': (14, 7)})\n",
    "\n",
    "# Suppress deprecation warnings.\n",
    "\n",
    "\n",
    "def fxn():\n",
    "    warnings.warn(\"deprecated\", DeprecationWarning)\n",
    "\n",
    "\n",
    "with warnings.catch_warnings():\n",
    "    warnings.simplefilter(\"ignore\")\n",
    "    fxn()\n",
    "\n",
    "\n",
    "def sigmoid(x):\n",
    "    '''Return value of sigmoid function (inverse of logit) at x.'''\n",
    "\n",
    "    return 1 / (1 + np.exp(-x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook we try to validate the derivation of the following equation:\n",
    "\n",
    "\\begin{equation}\\label{int}\n",
    "P(Y=0|do(R=r))=\\int_x P(Y=0|T=1, X=x)P(T=1|R=r, X=x)f_x(x)~dx\n",
    "\\end{equation}\n",
    "\n",
    "where $f_x(x)$ is the pdf of x. (Now $X \\sim U(0, 1)$, so $f_x(x)=1$.)\n",
    "\n",
    "First we generate data and estimate the causal effect from a causal model without unobservables Z."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0 % 10.0 % 20.0 % 30.0 % 40.0 % 50.0 % 60.0 % 70.0 % 80.0 % 90.0 % \n",
      "Analytical: 0.03709585053394618\n",
      "Estimated:  0.035425483336914504\n",
      "Difference: 0.0016703671970316747\n",
      "\n",
      "Values for P(y=0|do(r=1)) and P(y=0|do(r=0))\n",
      "\n",
      "Analytical: 0.14973849934787756 0.11264264881393138\n",
      "Estimated:  0.15289661989960157 0.11747113656268707\n",
      "\n",
      "Average difference: 0.0003938032875558321\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x1a1d351438>"
      ]
     },
     "execution_count": 80,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAasAAAEICAYAAADhmdstAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xe4FOXZ+PHvfc6h93LowqErRUVQUaQoolixR42KLSZR85qoSczrm8S0nyYxMUWj0aABuxILYhdFhChIR6TDQXrv/Zy9f388c3RZ9vTdfWZ378917bW7U++ZnZl755lnnhFVxRhjjAmzHN8BGGOMMeWxZGWMMSb0LFkZY4wJPUtWxhhjQs+SlTHGmNCzZGWMMSb0EpKsRGS+iAxJxLTSlYhcLCKrRGS3iPSpwPATReTm4PO3ReS9qH4DRGRJMK2LRKSliEwSkV0i8qdkLocJl2Ab6JSE6Q4RkdXVGP8xEfl5ImOq4HynVGT/SsB8CkXkzCqO62XdVJSItA+2q9wQxHKhiLxQkWHLTVbxfjQRuV5EJpd8V9WeqjqxnOkUiIiKSF5FAktDDwK3q2p9VZ1VmRFV9VlVPSuq06+Bh4NpvQbcAmwGGqrqXYkLObuJyHkiMllEtovIehF5QkQaRPWvJSJPisjOoP+dMeMPFZGFIrJXRD4SkQ4VHbeigm1gedWXsvpi93cAVf2eqv4mxXFcAOyq7P6VTGFZN2WJPYar6lfBdlWchHndJyLPVHR4VR0H9BKRY8sbNmOKAUOQBDsA85M0rQ7Al1qFO7hDsF5CoZT10Aj4LdAGOAZoB/wxqv99QFfc+j8d+ImIDA+m1xx4Bfg50BSYDrxYkXFNlX0PeNp3ECbhnsf9IS+bqpb5AgqBM2O6XQ9MjjcMcBJux90JbAD+HHT/ClBgd/A6BZcs/w9YCWwExgCNoqZ7XdBvC+6gED2f+4CxwDPBvG4O5v0psB1YBzwM1IyangK3AkuAXcBvgM7BODuBl6KHj1nmuLECtYLlUWAPsKyU8YcBC4EdQVwfAzfHrk9gGRAB9gXTfR44BBwMvp8ZxHJPMOyWIO6mwfgFQSw3Bet8UtC9P/DfYN3MAYZExTYxWBdTgvXyHtA8qv9pUeOuAq4PutfCnVF+FfzWjwF1gn7NgfHBOFuBT4CcUtbNqcDnwbr5HDg16H4lMD1m2B8B4yow/yHAauCnwHrg6Qps65cA86K+rwHOivr+G+CF4PMtwH+j+tULfrOjyxs3zny7BNvDDtwZ9Isx22yX4PO/gX8AbwfbwhSgFfAXYBtu++oTb9yo8X8bvX6i+pVsT7uAL4GLg+7HAPuB4mCe22OnFXz/DrA0+K3HAW1i4vgebr/bBjwCSHnLHrOOagbrt11Ut4rs76XNtzPwIW7/2Qw8CzSOPaYF63cv0CyqX19gE9C7gutmBDAbd4xZBgwvb1ssZR1Uen/DJffo48lP+OYYkRe1//8Wt4/vBt4AmgXrZCdunyyIiuOvuOPATmAGMDDoPhx3nDoUTGdO0L0RMCr4jdYE88qNmt4AYEW5y1+BFVRI5ZLVp8C1wef6QP+Yg2he1Hg34jbwTsGwrxAcVIAewQKfhttQHwxWQnSyOgRcFPwodYKNqD+QF8xvAfDDmI13HNAQ6AkcACYE82+E20lHlrIeSo013oEhZtzmwQ97GVADd8AtIk6yirfOOXLj/yHwGe5MoBbwT+D5mPU8BncArQO0xe2U5wbraljwPT9qY10GdAuGnwg8EPRrjzuAXRXE3gw4Puj3l2B9NgUa4Dby+4N+9+N2phrBayDBgSJm3TTFHUiuDX63q4LvzYC6wby7Rg3/OXBlBeY/JFjHvw/WUZ0KbOt/4Ztk1CRYjy2j+l9GkMxwO+yjMeN/AVxa3rhx5vs8cG/w29QGTou3XQXbwWbcdl4bd7BdgftTl4s7CHxU2jZJ2cnqctwZZg7wLdwfr9bxts840zojiOuEYF3/neBPUlQc44HGwfa0ieCAXdayx8yvJ7AnpltF9vfS5tsFtx/UAvKBScBfSjmmvQV8P6rfQ8DfK7huTsIl4mHBMrYl+ENT2RdV3N848nhSwJHJaikugZccBxfjknUe7ljyVNT41+D2zzzgLtyfwdpRx+VnYuJ+DXeMqge0AKYB3405BijuMke1k9VuXMYuee2l9GQ1CfgVUf/M462goNsE4Nao791xCSgP+AXBATjoVxeXtaOT1aRyYv8h8GrMxjsg6vsM4KdR3/9E1AYbM61SY413YIgZ9zrgs6jvgvvXX9VktQAYGvW9ddR6K1nPnaL6/5SYMwvgXYLEjNtY/y+q363AO8Hnn0Wvw5hl2AN0jup2CsE/JNx1t9dLWydR41wLTIvp9infnL09A/wi+NwVl7zqVmD+Q4LtpXYFDwTDcEmyW/D9qGA91o4ZpjD4PIogoUf1nxL8lmWOG2feY4DHiTpriNlmo5PVE1H9fgAsiPrem+DffbxtkjKSVZz5zgZGxNs+40xrFPCHqH71g+2xICqO6AT8EnBPecseM78BwPpyhom3v8edb5xxLwJmxdsHccl7SvA5F3dwPqmC6+afwEMV2QbLWbYq729ULFndG9X/T8DbUd8vAGaXEds24Ljg831EJSugJe6koE5Ut6s4/E9VjSCe9mWtg4pes7pIVRuXvHAHs9LchPuHvlBEPheR88sYtg2uWK3EStwBt2XQb1VJD1XdizsbiLYq+ouIdBOR8cEF7Z3A/8Od1UTbEPV5X5zv9asQa3lil0VjY6+kDsCrQcWA7bjkVRwTy6qY4S8vGT4Y5zRckiuxPurzXr5ZD0fhzrpi5eOSxoyoab4TdAd37Wcp8J6ILBeRe0pZltj1SvC9bfD5OdzGDXA18FqwLZQ3f4BNqrq/lPl+TUT6B/O5TFUXB513B+8NowZtiEuWJf2j+0X3L2/cWD/BHYymBTVrbywj3Kpuv2USketEZHbUuuzFkftOaQ77DVV1N25fbRs1TGnbV0WXfRvubCI65ors73HnKyItROQFEVkTjPtMnHFLvA70CGplDgN2qOq0UoaNVdr+c5igRvDu4PV2nEEStb+VpsLblYjcJSILRGRHEEcjSl93HXDJaF1U3P/EnWGVKPldt5cVYMIrWKjqElW9Kgjm98BYEamHy5yx1uIWpkR7XNHNBlz5ZruSHiJSB3fqedjsYr4/iiu376qqDYH/xe0IiVBWrOVZh9toARARif5eBauAc6L/QKhqbVVdEzWMxgz/dMzw9VT1gQrOq3Oc7ptxG3HPqGk2UtX6AKq6S1XvUtVOuH9md4rI0DjTiV2v4NZtybK8BzQXkeNxSeu5isw/zjqIK6gGPQ64UVUnfD2i6jbc73Zc1ODH8U3Fl/nR/YJtvDMwvwLjHkZV16vqd1S1DfBd4B8i0qW82CtgL+4AV6JVvIGCWoxPALfjrs00xhVpluw75a3Hw37DYF0045vfsFSVWPYlbtISnQCrs7/fj1uuY4Nxrylt3OAPz0vAt3ElAdGVPMpbN6XtP7HzeFZdDb36qnpOnEGqs7+Vux9UlIgMxJXUXAE0CbaVHZS+razCnVk1j4q7oar2jBrmGFypw86y5p3wZCUi14hIvqpG+CZTFuPKiyO4az4lngd+JCIdRaQ+7p/Ri6pahKs8cYGInCoiNXFFi+VtiA1w14Z2i8jRwPcTtmBlx1qeN4GeInJJUCvtfyjlwFFBjwG/K6kqLSL5IjKijOGfwa3Ls0UkV0Rqi7vPpl0Z45R4FjhTRK4QkTwRaSYixwe/7xPAQyLSIoijrYicHXw+X0S6BIl5J24biFdV9i2gm4hcHUz/W7jrleMBoraFP+LKtt8Pupc5/4oQkV64f6c/UNU34gwyBvg/EWkSbE/fwRXxALyKq3J7qYjUxhVbz1XVhRUYNzaOy6N+i224HT4R1YpnA1cHv/lwYHApw5X8mdwUxHMD7syqxAagXbAfxvMccIOIHC8itXD7xlRVLSwvwIouu6oeAj6IWYbq7O8NCC5vBAnwx+UMPwZX5Hchbn8qUd66GYVbN0NFJCfYRo+uRJxA+dt7OfvbBg4/7lZHA9yf9E1Anoj8gsNLEDYABSKSE8S9DveH808i0jBYB51FJPp3HIyrNFSmZFRdHw7MF5HduIvQV6rq/qDo5nfAlOB0sD/wJO5fyiTcheL9uHJ4VHV+8PkF3L/UXbhaeAfKmPfduKKiXbgf9sUyhq2sUmMtj6puxl3AfgBXPNIVd32jqv6KOxt4T0R24SpbnFzG/FfhaiT9L24jW4XbOcv9/VX1K1zFjLtwtYxm880Zw09xRQ+fBUUpH+Cu5YFbxg9wB4RPgX9onHvxVHULcH4w/S24YqHzg3VW4jncxd6XY/4clDX/irgLV4wyKqoIJvrs55e4IpyVuBprf1TVd4K4N+EqU/wOd5A9GVd7sdxx4zgRmBrsM+OAO1R1RSWWozR34P5lb8edFbwWbyBV/RJ3neJT3MGmN4dvnx/izgrXi8jmOONPwNXW/Q9uX+3M4euiLJVZ9n/izmxKVGd//xWuQsgO3J/JV8oaWFWn4P5sz4xJwuWtm2nADbhKGTtw20JsSUJFVXV/ux/3x2m7iNxdxXmXeBeXWBbjtu39HH7J4eXgfYuIzAw+X4erJPclbl8Zy+GXIK7C/bZlKqktEnrB2cx23Cl/InZkY0yaEXcD7g/Uw43BIvIh8Jyq/ivV885U4m70vlZVryh32DAnq2BBJuCK//6E+/d6goY5aGNMxhGRE3FF0EepamkVZUwShb0FixG4i7drcae5V1qiMsakkoiMxhWx/dASlT+hPrMyxhhjIPxnVsYYYwxp28hp8+bNtaCgwHcYxlTNokXuvXtlKi8aU30zZszYrKr55Q8ZLmmbrAoKCpg+fbrvMIypmiFD3PvEiT6jMFlIRGJbjEkLVgxojDEm9CxZGWOMCT1LVsYYY0LPkpUxxpjQs2RljDEm9CxZGWOMCT1LVsYYY0LPkpUxYaAKkYjvKIwJrbS9KdiYjLBnM8x4Cj5/EvZuhoZtoXk3OPV26DjId3TGhIYlK2N82bcd/tIbDu2FzmdAq96wYw2snAKjL4CCgXDO76Flz/KnZUyGs2RljA8HdsHGL6HpiXDpKGgR9aTzQ/vd2dakB+GJM+DcB+GEa0ufljFZwK5ZGZNqGxfAhvmQWxOueeXwRAVQozb0/z7c+ikcdRKMux1evw2KD/mJ15gQsGRlTCpFiuGV74AItOoFDVqWPmz9FnDtazDwbpj1DDx/JRzYnbpYjQkRS1bGpNLM0bB+HjTtBHm1yx8+JxeG/hwu+Bss+xDGXAh7tyY/TmNCxpKVMamybxtM+A10GAD1Kvk4ob4j4VvPwPov4OmLXeUMY7KIJStjUuWj+2H/dlfDryqOPs8lrA3z4dnLYP/OxMZnTIhZsjImFbavgs//BX2vd1XUq6rbWXDFaFg7C1642ipdmKxhVdeNSYWpj7n30+48rHPBPW8eMWjhA+eVPa2jz4MRj8Cr34U373TXs0QSFakxoWTJyphk278DZoyGnhdD46MSM83jroTNS+CTB4MWL36QmOkaE1JJKQYUkSdFZKOIfBHV7Y8islBE5orIqyLSOKrfz0RkqYgsEpGzkxGTMd7MHAMHd7kmlBLp9Huhxwh47+dQOCWx0zYmZJJ1zerfwPCYbu8DvVT1WGAx8DMAEekBXAn0DMb5h4jkJikuY1Kr+BB89qhrOqlNn8ROOycHRvwDmnZ0RYJWQ9BksKQkK1WdBGyN6faeqhYFXz8D2gWfRwAvqOoBVV0BLAVOSkZcxqTcl6/DzjXJK6arVR8ueQJ2roW37k7OPIwJAV/XrG4EXgw+t8UlrxKrg25HEJFbgFsA2rdvn8z4jEmMGf+GJgXQZdhhlSleWL7FfeifgHm06wdD7oGPfgfdz4VelyRgosaES8qrrovIvUAR8GxJpziDabxxVfVxVe2nqv3y8yt5U6UxqbZ5KRR+Aidc54rskum0O6H18fDuvXBwT3LnZYwHKU1WIjISOB/4tqqWJKTVQHQVqXbA2lTGZUxSzBwNkgvHfzv588rNg3P+ALvWwid/Tv78jEmxlCUrERkO/BS4UFX3RvUaB1wpIrVEpCPQFZiWqriMSYqigzD7Oeh+DjRolZp5tj8Zel8B//07bCtMzTyNSZFkVV1/HvgU6C4iq0XkJuBhoAHwvojMFpHHAFR1PvAS8CXwDnCbqhYnIy5jUmbRm+7Jv32vT+18h/0KcvJcdXZjMkhSKlio6lVxOo8qY/jfAb9LRizGeDFjNDQ6yj0BOJUatoEB/wMT73eN3rbqldr5G5Mk1jagMYm2rRCWf+SuVeV4uGXw5O9CzQbwyZ9SP29jksSSlTGJNusZQKDPNX7mX6cJnHQzzH/VNclkTAawZGVMIhUXwaxnocuZiWsHsCr63+Ye7jj5IX8xGJNAlqyMSaRlE1z18ROu8xtH/XxXuWPOC7D9K7+xGJMAlqyMSaQZo91TgLvFNo3pwam3AwrTn/IdiTHVZsnKmETZtR4WvwPHXw15NX1HA43aQbdzYNbT7r4vY9KYJStjEmX2c6DF0MdzEWC0fjfCnk2w8A3fkRhTLZasjEmESMQ9t6rDAGjexXc03+h8BjRub0WBJu1ZsjImEVZOhm0r4ISRviM5XE4O9L3BNai7abHvaIypMktWxiTCzDFQqxH0uNB3JEfqcy3k1IAZdnZl0pclK2Oqa+9W+HIcHHsF1KjjO5oj1c+H7sNh3svuPjBj0pAlK2Oqa+5LUHwA+oasCDBa7ytcRYsVH/uOxJgqsWRlTHVEIvD5E9C2L7Tq7Tua0nU9yxVTznvZdyTGVIklK2OqY+kHsGUpnPx935GUrUZt6HEBLHgDDu3zHY0xlWbJypjqmPoo1G8FPUb4jqR8va+Ag7th0du+IzGm0ixZGVNVmxbBsg/hxJvD0WJFeQpOgwatYd5Y35EYU2lJefiiMVlh6mOQWwv63eBl9gX3vHlEt8IHzit9hJxc6HUpTP0n7NsOdRonMTpjEsvOrIypij2bXYvmvS+Hes19R1NxPS6CyCFY8p7vSIypFDuzMqYqPn3YVVQYcEdKZhfvLKpK2vZ119gWvOHuCzMmTdiZlTGVtXcrTHsCel0C+d18R1M5OTlw9HmuFqPVCjRpxJKVMZX12T9crbqBd/uOpGqOPg8O7YVlH/mOxJgKs2JAYypj33ZXQeGYC6Flj6TMImFFfqXOYKC7QXjheDj63OTOy5gEsTMrYypj6j/hwE4Y9GPfkVRdXk3odjYsesvaCjRpI2nJSkSeFJGNIvJFVLemIvK+iCwJ3psE3UVE/iYiS0VkroickKy4jKmy/Tvhs0eg+3nQ+ljf0VTPMefDvm2wcorvSIypkGSeWf0bGB7T7R5ggqp2BSYE3wHOAboGr1uAR5MYlzFVM+1x2L8DBqfxWVWJLme6e8SsNQuTJpKWrFR1ErA1pvMIYHTweTRwUVT3Mep8BjQWkdbJis2YSjuwGz59xDUI26aP72iqr2Y96DgQlr7vOxJjKiTV16xaquo6gOC9RdC9LbAqarjVQbfDiMgtIjJdRKZv2rQp6cEa87Xpo2DfVhj0E9+RJE7Xs1wjvFuW+Y7EmHKFpYKFxOmmR3RQfVxV+6lqv/z8/BSEZQxwcC9M+Rt0PgOOOtF3NInTdZh7X2JnVyb8Ul11fYOItFbVdUEx38ag+2rgqKjh2gFrUxybMfHNeAr2bk6Ls6pKtRfYtBM06+qaXur/vSRHZkz1pPrMahxQ8jjVkcDrUd2vC2oF9gd2lBQXGuPVoX0w5a/u3qQOp/iOJvG6ngWFk+HgHt+RGFOmZFZdfx74FOguIqtF5CbgAWCYiCwBhgXfAd4ClgNLgSeAW5MVlzGVMvNp2L0BBv/UdyTJ0XUYFB+AFZ/4jsSYMiWtGFBVryql19A4wypwW7JiMaZKig7A5Ieg/anuWVClqPSjOlKszPg6nAo168OSd6F77J0mxoRHWCpYGBM+c16AXWth0N0g8eoAZYC8WtBxMCz5APSIOk3GhIYlK2PiiRS7a1Wtj3O1ADNZ59Nhx1ewdbnvSIwplSUrY+JZ8AZsXQan/Shzz6pKlCTjZR/6jcOYMliyMiaWqrtW1bSTa1090zXtBI3aw/KJviMxplT2iBBjYi2fCOtmwwV/hZxc39EkRWyli/vzOnPe9g9pWFwEuXZYMOFjW6UxsT59GOq3pPvLjTjw8uEH9TDV8kukKZFeXJX3EaydlVmtdJiMYcWAxkTbvNQ98r3fTRygpu9oUmZKpCcRFSsKNKFlycqYaNMeh5wa0Pd635Gk1DYaMl87wHJ71L0JJ0tWxpTYvxNmPwe9LoEGLX1Hk3KTI71h1TT3OBRjQsaSlTEl5jwPB3fBSd/1HYkXkyO9IHLInh5sQskqWBgDrrr6tMehbT9o17fak4vXxFHYTY90h7za7rpVt7N9h2PMYezMyhhwZxNblsKJN/uOxJsD1IT2p8Ayu25lwseSlTEAM8dArUbQY4TvSPzqfDpsWgA77Qk9JlwsWRmzbxt8+ToceznUrOs7Gr86DXHvVoXdhIwlK2PmjYWi/XDCdb4j8a9lb6jb3JKVCR1LVsbMHAOtjnUtrGe7nBzoNNglK3tkiAkRS1Ymu62dDevn2llVtE6nw+71sHGB70iM+ZolK5Pd5jwPubWg9+W+IwmPTkPcuxUFmhCxZGWyV/Ehd72q+3Co09h3NOHR+Cho1sWaXjKhYsnKZK9lH8HezXDslb4jCZ+Og2Hlf11CNyYELFmZ7DX3RajTFLqc6TuS8Ok4CA7udo8MMSYELFmZ7HRgFyx80zVam5c9jwKpsIKB7n3Fx37jMCZgycpkpwVvQNE+OPZbviMJp3rNoFVvWDHJdyTGAB6SlYj8SETmi8gXIvK8iNQWkY4iMlVElojIiyJif3VNcs19EZoUQDt7Km6pOg6Gr6bCoX2+IzEmtclKRNoC/wP0U9VeQC5wJfB74CFV7QpsA25KZVwmy+ze6M4Yel8OIr6jCa+Og6D4gHvGlTGe+SgGzAPqiEgeUBdYB5wBjA36jwYu8hCXyRYLxoFGoOfFviMJtw6nguRaUaAJhZQmK1VdAzwIfIVLUjuAGcB2VS0KBlsNtI03vojcIiLTRWT6pk2bUhGyyURfvArNu0OLHr4jCbdaDaBtX6tkYUIh1cWATYARQEegDVAPOCfOoHEbJVPVx1W1n6r2y8/PT16gJnPtXOeeXdXrEisCrIiOg2DNTNi/03ckJsuluhjwTGCFqm5S1UPAK8CpQOOgWBCgHbA2xXGZbPHl64BCz0t8R5IeOg0GLYavPvUdiclyqU5WXwH9RaSuiAgwFPgS+Ai4LBhmJPB6iuMy2WL+K9CyF+R38x1Jemh3kms7cbkVBRq/8sofJHFUdaqIjAVmAkXALOBx4E3gBRH5bdBtVCrjMllix2pYNRXO+L8qT6LgnjcTGFAaqFEb2p9slSyMdylNVgCq+kvglzGdlwMnpToWk2UWjHfvPawWYKV0HAwf/gb2bHE3CxvjgbVgYbLHgnGuBmDzLr4jSS8dB7v3Qju7Mv5YsjLZYfdG14r4MRf4jiT9tOkDNRtYUaDxypKVyQ4L3wQUjrnQdyTpJzfP3SBsycp4lPJrVsZ4sWAcNOkILXse0SvrKk1URafBsORd2LEGGsW9Z9+YpLIzK5P59m1zZwU9LrQbgauq4yD3bmdXxhNLVibzLXoHIkVWBFgdLXpC3WaWrIw3lqxM5ls4Hhq0gTYn+I4kfeXkuAcyrvgYNG5raMYklSUrk9kO7YdlH0G3s90B11Rdx0Gwcw1sXe47EpOFbO81ma1wMhzaA93jtZdsKqXkfitrhd14YLUBTWZb/Dbk1fmmgoApU7yakYUPnOc+NOsMDdu661b9bkxxZCbb2ZmVyVyqrnJF59OhRh3f0aQ/EZf0V0yCSMR3NCbL2JmVyVwbvoCdq2HwT77uZPdUVVPHQTDnedj4JbTq5Tsak0XszMpkrkXvuPduZ/uNI5PY/VbGE0tWJnMtfttVV2/QynckmaNRO2ja2SpZmJSzZGUy064NsGaG1QJMho6DoHAKFBf5jsRkEUtWJjMtede9dxvuN45M1HEQHNwF62b7jsRkEUtWJjMtegcatoNWvX1HknlKrlstn+g1DJNdLFmZzHNoPywPWq2whmsTr15zaNnLKlmYlLJkZTLPiklwaK9dr0qmjoNg1VT3x8CYFLD7rExGiL5/6rd5o7g4txb1CgZ6jCjDdRwMn/0DVk+z1kFMStiZlckwyhm5s5gc6Q01avsOJnN1OBUk14oCTcpYsjIZpYespI1s5YOIPQ4kqWo3hDZ9LFmZlLFkZTLKkBxXnXpi8fGeI8kCnQa7e9kO7PIdickClqxMRhmcO5d5kQI20dh3KJmv4yD3BObCKb4jMVkg5clKRBqLyFgRWSgiC0TkFBFpKiLvi8iS4L1JquMy6a8Be+kri/k4cpzvULLDUf3d41eWTfAdickCPs6s/gq8o6pHA8cBC4B7gAmq2hWYEHw3plIG5HxBnkT4uNiSVUrUqA0Fp8GyD31HYrJASpOViDQEBgGjAFT1oKpuB0YAo4PBRgMXpTIukxkG58xhp9ZhlnbxHUr26DIUtiyFbSt9R2IyXKrvs+oEbAKeEpHjgBnAHUBLVV0HoKrrRKRFvJFF5BbgFoD27dunJmKTJpTBuXOYEulFkd0+mFBlPj2481D3vmyCPT3YJFWqiwHzgBOAR1W1D7CHShT5qerjqtpPVfvl5+cnK0aThrrKGtrIViZGrBZgSjXvCo2OgqV23cokV6qT1WpgtapODb6PxSWvDSLSGiB435jiuEyaG5wzB4BJxcd6jiTLiEDn0939VsWHfEdjMlhKk5WqrgdWiUj3oNNQ4EtgHDAy6DYSeD2VcZn0NzhnDosi7VhHM9+hZJ/OQ+HATnfPlTFJ4qNw/wfAsyJSE1gO3IBLmi+JyE3AV8DlHuIy6ergHk7KWcjoYnt8vRedBoPkuKLA9v19R2MyVMqTlarOBvrF6TU01bGYDFE4mVpSxMcRKwL0ok4TaNvPVbI4417f0ZgMZS1YmPS39AP2ai0+jxyNKnlWAAAUNklEQVTtO5Ls1WUorJkJe7f6jsRkKEtWJv0teZ//RnpwkBq+I8lenYcC6h56aUwS2A0pJr1tWQbbVvBx5PojesW7P8gkSZs+ULuRa82i16W+ozEZyJKVSW/B/T3WHmBqxb1R+IQhsPRDUHVV2o1JICsGNOlt6QfQtBNfaUvfkZjOQ2HXWti00HckJgNZsjLpq+gAFH4CXc70HYkBV8kCrDULkxSWrEz6WjkFDu21ZBUWjdpB8+72yBCTFJasTPpa8j7k1YaCgb4jMSU6nwEr/wsH9/qOxGQYS1YmfS15zz1PqWZd35GYEt3OgqL9rq1AYxLIkpVJT1uWuecodT3LdyQmWocBULM+LH7bdyQmw1iyMulp6Qfu3a5XhUteLVcUuPhdV4XdmASxZGXS05L3oVkXaNbZdyQmVvdzYNc6WDfHdyQmg1iyMunn4F5XZd2KAMOpyzBA3NmVMQliycqkn8LJ7iJ+12G+IzHx1M+HdifadSuTUJasTPpZ8h7UqOsu5ptw6nY2rJ0Fu9b7jsRkCEtWJr2owpJ3oeNgdzHfhFP3c9z74nf8xmEyhjVka9LL5iWw/Svu3XQmz1qr6uHVogc0KYAF46Hv9b6jMRnAzqxMelnyHgATi62V9VATgWMugOUTYf8O39GYDGDJyqSXJe+xKNKONeT7jsSU55gREDlktQJNQliyMunjwG5Y+V8+ihzvOxJTEW37QoPWsGCc70hMBrBkZdLHio8hcoiJlqzSQ04OHH0+LPnAGrY11WbJyqSPJe9BzQZMj3TzHYmpqB4XQtG+b5rHMqaKLFmZ9KDqmljqPIQiq8SaPtqfCnWaWlGgqTYvyUpEckVkloiMD753FJGpIrJERF4UkZo+4jIhtnYW7FwD3c/1HYmpjNw8OOZ8WPQ2HNrnOxqTxnydWd0BLIj6/nvgIVXtCmwDbvISlQmvBW+A5EK34b4jMZXV+3I4uNtuEDbVkvJkJSLtgPOAfwXfBTgDGBsMMhq4KNVxmZBbON49aLFuU9+RmMrqMMDVCpz7su9ITBrzcWb1F+AnQCT43gzYrqpFwffVQNt4I4rILSIyXUSmb9q0KfmRmnDYtAg2L3Y3mZr0k5MLvS51FWT2bfMdjUlTKU1WInI+sFFVZ0R3jjNo3Ke2qerjqtpPVfvl59tNoVljwRvu/ejz/MZhqq73Ze4G4S+tooWpmlRXqxoAXCgi5wK1gYa4M63GIpIXnF21A9amOC4TZgvegLb9oGEb35GYMhTEaaux8IHgD0br493DMue9DH1HpjgykwlSemalqj9T1XaqWgBcCXyoqt8GPgIuCwYbCbyeyrhMiG1fBetmWxFguhNxFS0KJ8OONb6jMWkoLPdZ/RS4U0SW4q5hjfIcjwmLL4P/LZas0t+xVwAKs5/zHYlJQ96SlapOVNXzg8/LVfUkVe2iqper6gFfcZmQmfcytOkDzTr7jsRUV9NO0HEQzBoDkUj5wxsTJSxnVsYcafMSVwTY+3LfkZhEOWEkbP8KVkz0HYlJM9ZujQmveWMBgZ6X+I7EVFFspYta5LKoSROYOQY6n+EpKpOO7MzKhJOqKwIsOA0atvYdjUmQA9SEY690TxDes9l3OCaNWLIy4bRuNmxdZkWAmeiE69w9V3Oe9x2JSSOWrEw4zX0Zcmq4R0yYzNKyBxzVH6Y9AZFi39GYNGHJyoRP0UGY9xJ0OxvqNPEdjUmGU26F7Sth4ZE3EhsTjyUrEz6L34Y9m1zNMZOZjj4fGneATx/xHYlJE5asTPjMHAMN20KXob4jMcmSkwv9vw+rPoPVM8of3mQ9S1YmXLavgqUToM817oBmMlefa6BWQ/jMzq5M+SxZmXCZ9Yx773ON3zhM8tVq4Bq1nf8abF3uOxoTcpasTHhEil2y6nwGNG7vOxqTCqfcDrk1YNKDviMxIWfJyoTH4ndg52p7hEQ2adAK+t0Ec16ALct8R2NCzJKVCY/PHoVGR0F3e8hiVhlwB+TWhI//4DsSE2KWrEw4rP8CCj+Bk74DudZkZVZp0BJOvMndW7d5ie9oTEhZsjLhMPVRqFHXNcVjss+AOyCvDnxwn+9ITEhZsjL+7dnsmlc67iprsSJb1W8BA++EheNh+ce+ozEhZMnK+Pf5KCg+ACd/z3ckxqdTboNG7eHd/7U2A80RLFkZvw7sdkWA3YZDfjff0RifatSBs34NG76AWU/7jsaEjF3JNn7NHA37tsHAu494UB9A4QNWMzCr9LgI2p8KH/zK1Qqtn+87IhMSdmZl/Ck6AP/9OxQMhKNO9B2NCQMROP8hOLAL3v2Z72hMiFiyMv7MeR52rYOBd/mOxIRJi6Nh0N3uSdGL3/UdjQkJS1bGj+JDMPkv0OYE6DTEdzQmbE67E/KPgfF3wv6dvqMxIWDXrIwfc56HbStg+P2u6MdkjQpdm8yrCSMehlHD4O2fwMWPpSg6E1YpTVYichQwBmgFRIDHVfWvItIUeBEoAAqBK1R1WypjMylUdIDVr/+KzdqZi54qBuxpsdmu1AQ26Cfw8QPQdRj0utRDZCYsUl0MWATcparHAP2B20SkB3APMEFVuwITgu8mU80cQzvZzINFVwB2VmXKMOjH0O5EGP8j96wzk7VSemalquuAdcHnXSKyAGgLjACGBIONBiYCP01lbCZFDu2DSQ8yNXI0kyO9yh083j9uk0Vy8+CSx+GxgTD2Brj+Tcir5Tsq44G3ChYiUgD0AaYCLYNEVpLQWpQyzi0iMl1Epm/atClVoZpEmvoY7F7Pnw9djp1VmQpp2sldv1r9Obxj1dmzlZdkJSL1gf8AP1TVClf1UdXHVbWfqvbLz7ebBdPOni3wyZ+h23Cm6jG+ozHppOfFrrHb6aO+eZq0ySoprw0oIjVwiepZVX0l6LxBRFqr6joRaQ1sTHVcJgU+/j0c3APDfg1zl/qOxoRcbBFwLv1YdsxgV529RQ9oe4KnyIwPKT2zEhEBRgELVPXPUb3GASWPhx0JvJ7KuEwKbF7q/hX3HQn53X1HY9JQMblw2VOuhfYXr3Wt9ZuskepiwAHAtcAZIjI7eJ0LPAAME5ElwLDgu8kk7/8c8mrDELvmYKqhXjP41tOwZ5OrcFFc5DsikyKprg04mdKvqg9NZSwmhRa/C4vegjPvc/+KjamONn3ggr/Aa993jxM59w++IzIpYC1YmOQ6tA/e+jE07w79b/MdjckUx18NG+bDpw+72oL97Vlomc6SlUmuyQ/B9pUwcrxrQseYRBn2a9hW6Fpnb1IA3Yf7jsgkkTVka5Jn0yLXWG3vy6HjQN/RmEyTk+tuGG51LIy9EdbN8R2RSSJLViY5iovg1e9CzXpw9v/zHY3JVDXrceKK77DmYG3WPzaC/veMsVZPMpQVA5rkmPwQrJ0Fl4+m4Lef+47GZLBNNOHGgz9mbM1fMarmg1xx8Be+QzJJYGdWJvHWzXEtZfe6FHpe5DsakwUWaXtuO/Q/dJdV/K3Gw1alPQNZsjKJtW87vDQS6uXDuQ/6jsZkkUmR4/hl0fUMzZ3lqrSbjGLFgCZxIhF45RbYsdq1jl23qe+ITIYp73rUs8Vn0kE2cMu0f7oagqfcmprATNLZmZVJnI9/D0vedU//bX+y72hMlrq/6Co45gJXpX3eWN/hmASxZGUSY/pT7jrVcVfDiTf7jsZkMSUHLvkXdBgAr34Pln3oOySTAJasTPXNfdk9ybXr2XDBX0HsOVXGsxq14crnoHk31+jt2lm+IzLVZMnKVM+cF9z9VAWnwRWjrZUKEx51GsM1/4E6TeGZy2DLMt8RmWqwZGWqRhU+ut8lqg6nwlXPQ406vqMy5nANW8O1rwIKT18Mu9b7jshUkSUrU3l7t8LL17trVMd/G655BWo18B2VMfE17wJXv+yef/Xv8y1hpSlLVqZyln4Aj54KC8fDmb+CEY9Y0Z8Jv3Z9XZHgrnXw7/Ng5zrfEZlKsvusTMVsXgLv/xIWvQn5R8PVL0Lr444YzNplM2FQ2nZYeOt/4JlL4anh8O2x0LxriiMzVWVnVqZsGxfAa7fCIyfDikkw9Bdwy8dxE5Uxode+P1w3Dg7shlHDYOWnviMyFWRnVuZIxUXu5t7pT7pivxp14aRbYOBdUD//68HsLMqkpXZ94eYP4NnLYMwIOPt37t5Au+Ui1CxZmW9smO+qos99CXavhwat4fR73Y5sTSeZTNK0I9z0vqvN+tbdsOwjuPDvUK+Z78hMKSxZZbtd62HeyzDnRdgwD3LyoMswOOE66HoW5ObZGZTJGLHbsnAdK0YMcddj/36C+3PW70bItUNj2Ngvko0O7oEF45k09u8MyPmCXFFmRzrzSvFIxhefwta5DWGuAu/6jtSYpFJy4JTboPMZ8PZP4e0fw7TH4dTb4dhv2b2DIWLJKlsUHYTCSa5ppAVvwKE9dJR8HikewWvFp7Fc2/iO0Bh/WhwD170Oi96CiffDG3fAhF9Dr8ug1yXQ7iTIsfpoPlmyymQH98KyCbDgDXbOeYOGspedWpc3i0/mleKBTNdu7p+lMcZVsDj6POh+Llfe+yDXFb/H0KlPUmvaP9mkDcnveYZrVqzgNHf7hlXISClLVpmk6ABs+AIKp0DhJ7DiEyjaB3Wa8G5xP96JnMjkSG8OYDfxGlMqET6L9OCzSA/qs5czc2YyMHcul66eDl++5oap2xza9oXWx0Kr3tDqWPf8LEtgSROqZCUiw4G/ArnAv1T1Ac8hhUtxEezdDLs3uteejbBtJWxdBpsWunuiIsHjvJt1hT7fds/16TCAH9/7nt/YjUlDu6nLa5HTeC1yGpf+6FzYvhIKJ7s/hOtmU7T4ffIkAsBOrcMC7cDJpwxxxYpNO0Ozzq5WrSWxagtNshKRXOARYBiwGvhcRMap6pdJn7kqRIpBi4P3yOGfI8UQOQTFB6H4kDuDKS75frDUz794dRY1KKIWReRRTK4Uc8fpnVxCiQTTjxQF8ypyT9ot2gcH9/Dpgq+oK/upx37qBe8NZR+gh4UeUWEtzVgeac0Xei5fRAr4xz23ugY8jTHlqmht14KfvRV8agKcD5xPLQ7STVbTM6eQnlJIj5yVMHMMHNr7zYg16kHTTtCoLdRvAfVbBe8toFZD165mzXruVaOeq4mYE/WSHEt2hChZAScBS1V1OYCIvACMABKbrDYuhMeHuARRkohiEkCi/LpGnI5T8kBygw0xN3iVdMuFvNpQsx45EmGLNmQVLdgTqc1eanHD0OOhXj7Ub+k29Hr5HPPHuUcW61miMiYlDlCTedqJecWdvu5WeN9w2LHalXhsCV5bl8HONbBmJuzZRKWPOSXHjJKk1e1suGJM4hYkDYhqcg7UlSUilwHDVfXm4Pu1wMmqenvUMLcAtwRfuwOLUh5o6jUHNvsOIkWyaVnBljfThXV5O6hqfvmDhUuYzqzinecelklV9XHg8dSEEw4iMl1V+/mOIxWyaVnBljfTZdvyJluY6i2vBo6K+t4OWOspFmOMMSESpmT1OdBVRDqKSE3gSmCc55iMMcaEQGiKAVW1SERux7Xxkws8qarzPYcVBtlU7JlNywq2vJku25Y3qUJTwcIYY4wpTZiKAY0xxpi4LFkZY4wJPUtWHohIUxF5X0SWBO9NShluZDDMEhEZGdX9dyKySkR2xwxfS0ReFJGlIjJVRAqSuyQVk4Dl7Ssi84Ll+puIuzNSRO4TkTUiMjt4nZuqZYpHRIaLyKIgznvi9C/19xGRnwXdF4nI2RWdpi9JWtbC4HeeLSLTU7MkFVPV5RWRZiLykYjsFpGHY8aJu12bUqiqvVL8Av4A3BN8vgf4fZxhmgLLg/cmwecmQb/+QGtgd8w4twKPBZ+vBF70vawJWt5pwCm4e/HeBs4Jut8H3O17+YJYcoFlQCegJjAH6FGR3wfoEQxfC+gYTCe3ItPMlGUN+hUCzX0vX4KXtx5wGvA94OGYceJu1/aK/7IzKz9GAKODz6OBi+IMczbwvqpuVdVtwPvAcABV/UxV15Uz3bHA0JD8W6vy8opIa6Chqn6qbg8fU8r4vn3dXJiqHgRKmguLVtrvMwJ4QVUPqOoKYGkwvYpM04dkLGuYVXl5VXWPqk4G9kcPnEbbdWhYsvKjZUmyCd5bxBmmLbAq6vvqoFtZvh5HVYuAHUCzakdbfdVZ3rbB59juJW4Xkbki8mRpxYspUpHfq7Tfp6xlr+w2kArJWFZwLda8JyIzgqbVwqI6y1vWNMvark2M0NxnlWlE5AOgVZxe91Z0EnG6lXefQVXGSYgkLm9Zy/Qo8Jvg+2+APwE3VnB+iVaRdV/ZZYz3ZzIM95okY1kBBqjqWhFpAbwvIgtVdVI14kyU6ixvdaZpoliyShJVPbO0fiKyQURaq+q6oDhgY5zBVgNDor63AyaWM9uSJqtWi0ge0AjYWpm4qyqJy7s6+BzdfW0wzw1R83gCGF/V+BOgIs2Flfb7lDVuGJsgS8qyqmrJ+0YReRVX/BaGZFWd5S1rmnG3axOfFQP6MQ4oqe02Eng9zjDvAmeJSJOgeOusoFtFp3sZ8GFQHu5blZc3KDbcJSL9g2se15WMHyS+EhcDXyRrASqgIs2Flfb7jAOuDGqUdQS64i6+h7UJsoQvq4jUE5EGACJSD/f7+/w9o1VneeMqa7s2pfBdwyMbX7iy7AnAkuC9adC9H+4JySXD3Yi7AL0UuCGq+x9w/8wiwft9QffawMvB8NOATr6XNUHL2w934FoGPMw3La88DcwD5uIOFq09L+e5wOIgznuDbr8GLizv98EVly7DPfbmnLKmGYZXopcVV9NuTvCaH6ZlTcDyFuLOsnYH+2uPsrZre8V/WXNLxhhjQs+KAY0xxoSeJStjjDGhZ8nKGGNM6FmyMsYYE3qWrIwxxoSeJStjjDGhZ8nKGGNM6P1/3qZgYTEgCwsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "###\n",
    "# Synthetic data 1 - without unobservables Z.\n",
    "\n",
    "# R ~ U(0,1)\n",
    "# X ~ U(0,1)\n",
    "# T ~ Bin(1, sigmoid(a * r + b * x))\n",
    "# Y ~ Bin(1, sigmoid(c * t + d * x)), but if T = 0 then Y = 1\n",
    "\n",
    "# Weights for edges:\n",
    "# R -> T: a\n",
    "# X -> T: b\n",
    "# T -> Y: c\n",
    "# X -> Y: d\n",
    "\n",
    "N = 10000\n",
    "a = 1\n",
    "b = 1\n",
    "c = 1\n",
    "d = 1\n",
    "\n",
    "\n",
    "def generateData(N=N, a=a, b=b, c=c, d=d):\n",
    "\n",
    "    r = npr.uniform(size=N)\n",
    "    x = npr.uniform(size=N)\n",
    "\n",
    "    t = npr.binomial(n=1, p=sigmoid(a * r + b * x), size=N)\n",
    "\n",
    "    y = npr.binomial(n=1, p=sigmoid(c * t + d * x), size=N)\n",
    "\n",
    "    # If decision is negative, we set Y to positive.\n",
    "    y[t == 0] = 1\n",
    "\n",
    "    return r, x, t, y\n",
    "\n",
    "\n",
    "def analytic_R_on_Y(r, a, b, c, d):\n",
    "    '''Return probability of Y=0 with given parameters.'''\n",
    "\n",
    "    # P(y=0|t=1, x) * p(t=1|r, x) * p(x)\n",
    "    t_1 = si.quad(\n",
    "        lambda x: (1 - sigmoid(c * 1 + d * x)) * sigmoid(a * r + b * x) * 1, 0,\n",
    "        1)\n",
    "\n",
    "    # P(y=0|t=0, x) * p(t=0|r, x) * p(x)\n",
    "    # Now equal to 0, hence commented out. Y=0 and T=0 mutually exclusive.\n",
    "    t_0 = si.quad(\n",
    "        lambda x: (1 - sigmoid(c * 0 + d * x)) * (1 - sigmoid(a * r + b * x)) *\n",
    "        1, 0, 1)\n",
    "\n",
    "    return t_1[0] #+ t_0[0]\n",
    "\n",
    "# Calculate the difference between the analytic value and an estimate obtained \n",
    "# from synthetic data.\n",
    "\n",
    "diffs = np.zeros(0)\n",
    "\n",
    "for i in range(2000):\n",
    "\n",
    "    r, x, t, y = generateData()\n",
    "\n",
    "    # Fit predictive models.\n",
    "    lr_t = LogisticRegression(solver='lbfgs')\n",
    "    lr_y = LogisticRegression(solver='lbfgs')\n",
    "\n",
    "    lr_t = lr_t.fit(np.array([r, x]).T, t)\n",
    "    lr_y = lr_y.fit(np.array([t, x]).T, y)\n",
    "\n",
    "    def causal_effect_of_R_on_Y(r):\n",
    "\n",
    "        # Integrate P(Y=0|T=1, X=x) * P(T=1|R=r, X=x) * f(x) from 0 to 1\n",
    "        return si.quad(\n",
    "            lambda x: lr_y.predict_proba(np.array([[1, x]]))[0, 0] * lr_t.\n",
    "            predict_proba(np.array([[r, x]]))[0, 1], 0, 1)\n",
    "    \n",
    "    r0 = causal_effect_of_R_on_Y(0)\n",
    "    r1 = causal_effect_of_R_on_Y(1)\n",
    "\n",
    "    analytical = analytic_R_on_Y(1, a, b, c, d) - \\\n",
    "        analytic_R_on_Y(0, a, b, c, d)\n",
    "\n",
    "    diffs = np.append(diffs, analytical - r1[0] + r0[0])\n",
    "\n",
    "    if i % 200 == 0:\n",
    "        print(i / 20, \"%\", end=\" \")\n",
    "\n",
    "print()\n",
    "print(\"Analytical:\", analytical)\n",
    "print(\"Estimated: \", r1[0] - r0[0])\n",
    "print(\"Difference:\", analytical - r1[0] + r0[0])\n",
    "print()\n",
    "print(\"Values for P(y=0|do(r=1)) and P(y=0|do(r=0))\\n\")\n",
    "print(\"Analytical:\", analytic_R_on_Y(1, a, b, c, d),\n",
    "      analytic_R_on_Y(0, a, b, c, d))\n",
    "print(\"Estimated: \", r1[0], r0[0])\n",
    "print()\n",
    "print(\"Average difference:\", np.mean(diffs))\n",
    "plt.hist(diffs, bins=50)\n",
    "plt.title(\"Histogram of differences over 2000 simulations (analytic - estimate)\")\n",
    "sns.kdeplot(diffs)\n",
    "plt.axvline(x=0, c='r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It can be observed that when there are no unobservables, the causal effect $$P(Y=0|do(R=r))=\\int_x P(Y=0|T=1, X=x)P(T=1|R=r, X=x)f_x(x)~dx$$ can be estimated without bias.\n",
    "\n",
    "Next we generate data and estimate the causal effect from a causal model with unobservables Z. \n",
    "\n",
    "Z influences both T and Y."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0 % 10.0 % 20.0 % 30.0 % 40.0 % 50.0 % 60.0 % 70.0 % 80.0 % 90.0 % \n",
      "Analytical: 0.02047548901524536\n",
      "Estimated:  0.019144988392514814\n",
      "\n",
      "Values for P(y=0|do(r=1)) and P(y=0|do(r=0))\n",
      "\n",
      "Analytical: 0.10799301700601992 0.08751752799077456\n",
      "Estimated:  0.10764148848083047 0.08849650008831565\n",
      "\n",
      "Average difference: 0.0012700827692348328\n",
      "Std of differences: 0.0017989932120400834\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x109619a90>"
      ]
     },
     "execution_count": 81,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAasAAAEICAYAAADhmdstAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XecVPX1//HXmdldem8ioEgVUUREQRFFsGAFY4ktMZZooiaamBhjvsnPJBpJjNGosUaNBXvFgnSkqCDSkV4EpCNVyu7OPb8/PndlGGd3h92dvXdmzvPxmMfM3DL3PXfuvWfuvZ+5I6qKMcYYE2aRoAMYY4wx5bFiZYwxJvSsWBljjAk9K1bGGGNCz4qVMcaY0LNiZYwxJvSqpFiJyDwR6VcVr5WpROQCEVklIjtF5JgUhh8vItf5j68QkZFx/fqIyGL/tQaLSAsRmSAiO0Tk/nS+DxMu/jLQLg2v209EVldi/MdF5I9VmSnF6U5OZf2qgumsEJHTKjhuIPMmVSJyiL9cRUOQ5XwReSWVYcstVsk+NBH5iYhMKnmuql1VdXw5r9NWRFRE8lIJloH+CdysqnVVdcaBjKiqQ1X1jLhOfwEe8V/rHeB6YBNQX1Vvq7rIuU1EzhGRSSKyVUTWichTIlIvrn8NEXlGRLb7/X+dMP4AEVkgIrtEZJyIHJrquKnyl4FlFX+XlZe4vgOo6s9U9a/VnOM8YMeBrl/pFJZ5U5bEbbiqrvSXq1gapnWXiLyY6vCqOgw4UkS6lTds1hwGDEERPBSYl6bXOhT4UivwC+4QzJdQKGU+NADuBg4GugCtgfvi+t8FdMTN/1OB20VkoP96TYG3gD8CjYFpwKupjGsq7GfAC0GHMFXuZdwX8rKpapk3YAVwWkK3nwCTkg0DHI9bcbcD64F/+d1XAgrs9G8n4Irl/wFfARuA54EGca/7Y7/fZtxGIX46dwFvAC/607rOn/anwFZgLfAIUBD3egrcCCwGdgB/Bdr742wHXosfPuE9J80K1PDfjwLfAktLGf90YAGwzc/1MXBd4vwElgIesNt/3ZeBIqDQf36an+UOf9jNfu7G/vht/SzX+vN8gt+9N/CJP29mAf3iso3358Vkf76MBJrG9T8pbtxVwE/87jVwe5Qr/c/6caCW368p8L4/zjfARCBSyrw5EfjcnzefAyf63S8FpiUM+ytgWArT7wesBn4HrANeSGFZ/wEwJ+7518AZcc//CrziP74e+CSuXx3/Mzu8vHGTTLeDvzxsw+1Bv5qwzHbwH/8PeBQY7i8Lk4GDgAeBLbjl65hk48aNf3f8/InrV7I87QC+BC7wu3cB9gAxf5pbE1/Lf/5TYIn/WQ8DDk7I8TPcercF+A8g5b33hHlU4M/f1nHdUlnfS5tue2Asbv3ZBAwFGiZu0/z5uwtoEtfvWGAjcFSK82YQMBO3jVkKDCxvWSxlHhzw+oYr7vHbk9vZt43Ii1v/78at4zuB94Am/jzZjlsn28bl+DduO7Ad+ALo63cfiNtOFfmvM8vv3gB42v+MvvanFY17vT7A8nLffwozaAUHVqw+BX7kP64L9E7YiObFjXcNbgFv5w/7Fv5GBTjCf8Mn4RbUf/ozIb5YFQGD/Q+llr8Q9Qby/OnNB25NWHiHAfWBrsBeYIw//Qa4lfSqUuZDqVmTbRgSxm3qf7AXAfm4DW4xSYpVsnnO9xf+W4HPcHsCNYAngJcT5vPzuA1oLaAVbqU8259Xp/vPm8UtrEuBTv7w44Ehfr9DcBuwy/zsTYDufr8H/fnZGKiHW8jv9fvdi1uZ8v1bX/wNRcK8aYzbkPzI/9wu8583AWr70+4YN/znwKUpTL+fP4//7s+jWiks6w+yrxg18udji7j+F+EXM9wK+1jC+HOBC8sbN8l0Xwb+4H82NYGTki1X/nKwCbec18RtbJfjvtRFcRuBcaUtk5RdrC7G7WFGgB/ivni1TLZ8Jnmt/n6uHv68fhj/S1JcjveBhv7ytBF/g13We0+YXlfg24RuqazvpU23A249qAE0AyYAD5ayTfsQ+HlcvweAh1OcN8fjCvHp/ntshf+F5kBvVHB94/vbk7Z8v1gtwRXwku3gIlyxzsNtS56NG/9K3PqZB9yG+zJYM267/GJC7ndw26g6QHNgKnBDwjZAcac5Kl2sduIqdsltF6UXqwnAn4n7Zp5sBvndxgA3xj3vjCtAecCf8DfAfr/auKodX6wmlJP9VuDthIW3T9zzL4DfxT2/n7gFNuG1Ss2abMOQMO6Pgc/ingvuW39Fi9V8YEDc85Zx861kPreL6/87EvYsgBH4hRm3sP5fXL8bgY/8x7+Pn4cJ7+FboH1ctxPwvyHhzru9W9o8iRvnR8DUhG6fsm/v7UXgT/7jjrjiVTuF6ffzl5eaKW4ITscVyU7+8zb+fKyZMMwK//HT+AU9rv9k/7Msc9wk034eeJK4vYaEZTa+WD0V1+8XwPy450fhf7tPtkxSRrFKMt2ZwKBky2eS13oa+Edcv7r+8tg2Lkd8AX4NuKO8954wvT7AunKGSba+J51uknEHAzOSrYO44j3ZfxzFbZyPT3HePAE8kMoyWM57q/D6RmrF6g9x/e8Hhsc9Pw+YWUa2LcDR/uO7iCtWQAvcTkGtuG6Xsf+Xqnw/zyFlzYNUz1kNVtWGJTfcxqw01+K+oS8Qkc9F5Nwyhj0Yd1itxFe4DW4Lv9+qkh6qugu3NxBvVfwTEekkIu/7J7S3A3/D7dXEWx/3eHeS53UrkLU8ie9FE7MfoEOBt/2GAVtxxSuWkGVVwvAXlwzvj3MSrsiVWBf3eBf75kMb3F5Xoma4ovFF3Gt+5HcHd+5nCTBSRJaJyB2lvJfE+Yr/vJX/+CXcwg1wOfCOvyyUN32Ajaq6p5TpfkdEevvTuUhVF/mdd/r39eMGrY8rliX94/vF9y9v3ES34zZGU/2WtdeUEbeiy2+ZROTHIjIzbl4eyffXndLs9xmq6k7cutoqbpjSlq9U3/sW3N5EfOZU1vek0xWR5iLyioh87Y/7YpJxS7wLHOG3yjwd2KaqU0sZNlFp689+/BbBO/3b8CSDVNX6VpqUlysRuU1E5ovINj9HA0qfd4fiitHauNxP4PawSpR8rlvLCljlDSxUdbGqXuaH+TvwhojUwVXORGtwb6bEIbhDN+txxzdbl/QQkVq4Xc/9Jpfw/DHccfuOqlofuBO3IlSFsrKWZy1uoQVARCT+eQWsAs6K/wKhqjVV9eu4YTRh+BcShq+jqkNSnFb7JN034RbirnGv2UBV6wKo6g5VvU1V2+G+mf1aRAYkeZ3E+Qpu3pa8l5FAUxHpjitaL6Uy/STzICm/GfQw4BpVHfPdiKpbcJ/b0XGDH82+hi/z4vv5y3h7YF4K4+5HVdep6k9V9WDgBuBREelQXvYU7MJt4EoclGwgvxXjU8DNuHMzDXGHNEvWnfLm436foT8vmrDvMyzVAbz3xe6lJb4AVmZ9vxf3vrr5415Z2rj+F57XgCtwRwLiG3mUN29KW38SpzFUXQu9uqp6VpJBKrO+lbsepEpE+uKO1FwCNPKXlW2Uvqyswu1ZNY3LXV9Vu8YN0wV31GF7WdOu8mIlIleKSDNV9dhXKWO448Ue7pxPiZeBX4nIYSJSF/fN6FVVLcY1njhPRE4UkQLcocXyFsR6uHNDO0XkcODnVfbGys5ang+AriLyA79V2i8pZcORoseBe0qaSotIMxEZVMbwL+Lm5ZkiEhWRmuJ+Z9O6jHFKDAVOE5FLRCRPRJqISHf/830KeEBEmvs5WonImf7jc0Wkg1+Yt+OWgWRNZT8EOonI5f7r/xB3vvJ9gLhl4T7cse1Rfvcyp58KETkS9+30F6r6XpJBngf+T0Qa+cvTT3GHeADexjW5vVBEauIOW89W1QUpjJuY4+K4z2ILboWvimbFM4HL/c98IHBKKcOVfJnc6Oe5GrdnVWI90NpfD5N5CbhaRLqLSA3cujFFVVeUFzDV966qRcDohPdQmfW9Hv7pDb8A/rac4Z/HHfI7H7c+lShv3jyNmzcDRCTiL6OHH0BOoPzlvZz1bT37b3crox7uS/pGIE9E/sT+RxDWA21FJOLnXov7wnm/iNT350F7EYn/HE/BNRoqUzqarg8E5onITtxJ6EtVdY9/6OYeYLK/O9gbeAb3LWUC7kTxHtxxeFR1nv/4Fdy31B24Vnh7y5j2b3CHinbgPthXyxj2QJWatTyqugl3AnsI7vBIR9z5jYr6N25vYKSI7MA1tuhVxvRX4Vok3YlbyFbhVs5yP39VXYlrmHEbrpXRTPbtMfwOd+jhM/9QymjcuTxw73E0boPwKfCoJvktnqpuBs71X38z7rDQuf48K/ES7mTv6wlfDsqafipuwx1GeTruEEz83s//wx3C+QrXYu0+Vf3Iz70R15jiHtxGtheu9WK54yZxHDDFX2eGAbeo6vIDeB+luQX3LXsrbq/gnWQDqeqXuPMUn+I2Nkex//I5FrdXuE5ENiUZfwyute6buHW1PfvPi7IcyHt/ArdnU6Iy6/ufcQ1CtuG+TL5V1sCqOhn3ZXt6QhEub95MBa7GNcrYhlsWEo8kpKqi69u9uC9OW0XkNxWcdokRuMKyCLds72H/Uw6v+/ebRWS6//jHuEZyX+LWlTfY/xTEZbjPtkwlrUVCz9+b2Yrb5a+KFdkYk2HE/QD3FxrAD4NFZCzwkqr+t7qnna3E/dD7R6p6SbnDhrlY+W9kDO7w3/24b689NMyhjTFZR0SOwx2CbqOqpTWUMWkU9itYDMKdvF2D28291AqVMaY6ichzuENst1qhCk6o96yMMcYYCP+elTHGGEPGXuS0adOm2rZt26BjmEyxcKG773wgjQWNyT5ffPHFJlVtVv6Q4ZKxxapt27ZMmzYt6BgmU/Tr5+7Hjw8yhTGBE5HEK8ZkBDsMaIwxJvSsWBljjAk9K1bGGGNCz4qVMcaY0LNiZYwxJvSsWBljjAk9K1bGGGNCz4qVMcaY0LNiZYwxJvQy9goWxqRb2zs+SGm4FUPOSXMSY4ztWRljjAk9K1bGGGNCz4qVMcaY0LNiZYwxJvSsWBljjAk9K1bGGGNCz4qVMcaY0LPfWRlTKqWHLKatrKOFbOUrbc5H3vF49h3PmGpnxcrklFR/6FtAEX/Pe5Yf5o3fr/sy7yAeKR7MW15fQKo+oDEmKStWxiRoxlYeL3iAYyOLebh4MK/HTmGTNqBvZDa/yHuHfxU8Tu/i+dxZfC3FtgoZUy1sTTMmjuDxSMFDdJGV/LzwFoZ7vb7rN8I7nhGFx/GrvDe5Je8tDpJvuLHolgDTGpM77OC7MXEujY6jV2QBdxX/eL9CtY/wQPFF/Lboek6MzOP5giGwd0e15zQm11ixMsbXjC38Pu9lPokdwWuxfmUO+3qsHzcV3UI3WQZDL4a9O6snpDE5yoqVMb678p+jBkXcWXwtqTSeGOEdxy1FN8OqKfDypVC4K/0hjclRVqyMATrLSs6JTuXR4vNZoS1THu8Drzdc8ASsmASvXAZFe9KY0pjcZcXKGODq6Efs1gL+FzvzwEfudgkMfhSWfQyvXmEFy5g0sGJlcl4TtnFBdDJvxvqyjboVe5Hul8N5/4Ylo+H582HnxqoNaUyOs2Jlct7l0THUkCKejQ2s3AsdexVc/BysnQ1P9Yd1c6smoDHGipXJbQUU8aO80YyLHc1SbVX5F+w6GK7+EGKF8NSpMP7vUFxY+dc1JselpViJSBsRGSci80Vknojc4ndvLCKjRGSxf9/I7y4i8pCILBGR2SLSIx25jEk0MPI5zWVr5feq4rXqAT+bBF3Oh/F/gyf6wpfvgudV3TSMyTHpuoJFMXCbqk4XkXrAFyIyCvgJMEZVh4jIHcAdwO+As4CO/q0X8Jh/b0xanRf9hDXamIneUVX7wnWbwUVPu8YXI+6E137MQq81DxdfwIderzIvhrtiyDlVm8WYLJCWPStVXauq0/3HO4D5QCtgEPCcP9hzwGD/8SDgeXU+AxqKSOrth42pgAbs5JTILN6PnYCm64h4pzPhpqlw4dNEUB4peJiRBbdzXuQTQNMzTWOyUNrPWYlIW+AYYArQQlXXgitoQHN/sFbAqrjRVvvdjEmbM6OfUyAx3oudkN4JRaJw1EWcUfh3biz8JUVEebjgEV4r+Atd5Kv0TtuYLJHWYiUidYE3gVtVdXtZgybp9r2vnSJyvYhME5FpGzda02BTOedFPmWF14I5eli1TE+J8KHXm7ML7+X2op/SXtbwfsGdXB99D9vLMqZsaStWIpKPK1RDVfUtv/P6ksN7/v0Gv/tqoE3c6K2BNYmvqapPqmpPVe3ZrFmzdEU3OaAp2zgxMo9h3glU9/9SKRFei53KqXvvZ7h3PHfmv8z9+Y9RA2s1aExp0tUaUICngfmq+q+4XsOAq/zHVwHvxnX/sd8qsDewreRwoTHpcFZ0ClFR3oudGFiG7dTl5qJf8s+ii7kwOoln8u+jgKLA8hgTZulqDdgH+BEwR0Rm+t3uBIYAr4nItcBK4GK/34fA2cASYBdwdZpyGQPAudHPWOi1ZrG2DjiJ8EjsAtZoE/5V8Dj/5HHwzoeI/QTSmHhpKVaqOonSj60MSDK8AjelI4sxiRqwk56ykP/EBlXJ67W944NKv8Zb3sk0K9rG7/NfhlF/hDPvqYJkxmQP+/pmcs4pkdlERRkbC9dvz5+IncuzxWfCp4+4HxEbY75jxcrknP7R6WzS+szSdkFHSSDcU3wFtOwO790KO9YHHciY0EjXOStjQilKjH6RWYz2jk3fD4EroZg8Bqy4gg8L7mTCPy7lp0W3UdoRdbvShckl4VtbjUmjHrKYhvItY2LHBB2lVEu1Ff8o/iGnR6fzg8jEoOMYEwpWrExOGRCdQZFGq/5agFXsmdhAZngduD3/VWphf+ZojBUrk1P6R6YzxTucndQOOkqZlAh3F13BQbKF66IfBh3HmMBZsTK5o3gPnSJfM84L7yHAeF9oZz6MHc/P8t6jGVuCjmNMoKxYmdyxeysA472jAw6Sur8XX0o+xfw6742goxgTKCtWJnfs2co6bcRSPTjoJCn7Sg9iaOw0Lo5+TCvs4s0md1mxMrlj91Yme0dS3Reurawnis/FQ7guz85dmdxlxcrkhsKd4BUzKXZk0EkO2Dqa8E7sJC6NjqMxZf3TjjHZy4qVyQ3++Sq3Z5V5noidSw2K+EneR0FHMSYQVqxMbtizFfJrs4FGQSepkKXaihFeT66KjqQOu4OOY0y1s2Jlsl/RHtizHWo1DDpJpTxefB4NZBcXRicEHcWYamfFymS/VVNAPaiZ2cVqlnZglteOK6JjAA06jjHVyoqVyX7LxgMCNRsEnaTShsYG0DmymuNkYdBRjKlWVqxM9ls2HmrUg0g06CSV9l7sBLZrba7MGx10FGOqlRUrk912b4E1MzL+fFWJ3dTkzVhfzopMgZ32I2GTO6xYmey2fCKgGX++Kt7Q2AAKJAYzXww6ijHVxoqVyW7LxkNBXXcYMEss0dZM8Q6HGS+CWkMLkxusWJnstmw8tD0JJLMusVSeN2N9YfMS+Hp60FGMqRZWrEz22roSvlkK7foFnaTKDY/1gryaMOvloKMYUy2sWJnstWy8u2/XL8AQ6bGD2tD5bJj7BhQXBh3HmLSzYmWy17LxULcFNDs86CTpcfRlrrXj4pFBJzEm7axYmezkebDsYzjslKw7X/Wd9v2hTjOY/UrQSYxJOytWJjttmAe7NkH7U4NOkj7RPDjqYlj4kdvDMiaLWbEy2SmLz1ft58iLwCuCBfbHjCa7WbEy2WnpOGjaGepnzl/YV0irHtDgEJj3dtBJjEmrvKADGFMV2t7xwXePCyhiVo2JvBI7lT/73V9Zttn17B1EujQSga6D4LPH3KHAWpn5f13GlMf2rEzW6RFZTC0pzNh/BT5gXS8ArxgWfFD+sMZkKCtWJuv0icylWCN85nUJOkr1OLgHNDwE5r0TdBJj0saKlck6fSNzmKkd2EntoKNUDxE4YjAsGwe7vgk6jTFpYcXKZJX67OQoWZY7hwBLdB3sDgUutFaBJjtZsTJZ5YTIfKKiTIzlWLE6uAc0aAPz3w86iTFpYcXKZJU+kbns1JrM1A5BR6leIu5agcvGQeG3QacxpspZ03WTVU6KzGGK14XiHFi045vrA5wQacrLBXu44a77GOEd9133FUPOqe5oxlQ527MyWaMVG2kXWZd756t8n3ud2ap1OD36RdBRjKlyaSlWIvKMiGwQkblx3e4Ska9FZKZ/Ozuu3+9FZImILBSRM9ORyWS/E6PzAJjoHRVwkmAUk8dY7xgGRKYTJRZ0HGOqVLr2rP4HDEzS/QFV7e7fPgQQkSOAS4Gu/jiPikg0TblMFjspMpcN2pDF2iroKIEZGetJI9nJcZGFQUcxpkqlpVip6gQg1R98DAJeUdW9qrocWAIcn45cJnsJHn0ic5nkHQlk6V+CpGCC1429ms/pETsUaLJLdZ+zullEZvuHCUsuYtYKWBU3zGq/2/eIyPUiMk1Epm3cuDHdWU0GOVxW0VS2MznXmqwn2EVNJnlHcnpkGqBBxzGmylRnsXoMaA90B9YC9/vdk30NTrqWqeqTqtpTVXs2a9YsPSlNRuoTcadHJ+Vo44p4Y71jOCSykfayJugoxlSZaitWqrpeVWOq6gFPse9Q32qgTdygrQFby8wB6RuZwxLvYNbTOOgogRsfOxqAfpFZAScxpupUW7ESkZZxTy8ASloKDgMuFZEaInIY0BGYWl25TBYo3EXvyHw+9o4OOkkofE0zFnmt6BeZGXQUY6pMWn45KSIvA/2ApiKyGvh/QD8R6Y47xLcCuAFAVeeJyGvAl0AxcJOqWrtbk7qvJlNDivjY6xZ0ktAY53Xn6uhH1GZP0FGMqRJpKVaqelmSzk+XMfw9wD3pyGJywOJR7NYCpuTKX4KkYLzXnRvyPvDP5V0YdBxjKs2uYGEy35JRfOodwV4Kgk4SGtO8zuzQWpxqhwJNlrBiZTLb5qXwzTI7X5WgiDwme0dySnQWqDVhN5nPipXJbEvGADDeitX3jPO600o2w4b5QUcxptKsWJnMtmQUNG7HV3pQ0ElCp6QJO4tHBhvEmCpgxcpkrqI9sHwidDgt6CShtJ7GzPcOgcWjgo5iTKVZsTKZa/kEKN4NHc8IOklojfO6w6rPYM+2oKMYUylWrEzmWjQc8utA275BJwmtcbHu4BXDsvFBRzGmUqxYmcykCguHQ4f+kF8z6DShNV07Qo0Gdt7KZDwrViYzrZ0JO9ZC57PLHzaHxYhC+1Nh8Whrwm4ymhUrk5kWDgeJ2PmqVHQ8A3aug3Vzgk5iTIVZsTKZaeGH0KYX1GkadJLwK2ktaYcCTQazYmUyz9ZVbi+h81lBJ8kM9VpAy6NhyeigkxhTYVasTOZZ9JG7t/NVqWs/AFZNhT3bg05iTIVYsTKZZ+GH0KQDNO0YdJLM0b4/aMz9Ns2YDGTFymSWPdvdVStsr+rAtOkFBXVh6digkxhTIVasTGZZOga8IitWByqvwP14eumYoJMYUyFWrExmWTgcajWGNscHnSTztO8PW1bAN8uCTmLMAbNiZTJHrBgWjYBOAyESDTpN5ukwwN0vsb0rk3msWJnMseoz2LPVmqxXVON20PAQWDou6CTGHDArViZzLBwO0QJ3OMscOBHXhH35BIgVBZ3GmANixcpkBlVY8AEcdgrUqBt0mszVvj8U7oDVnwedxJgDYsXKZIZNi2DLcjsEWFmHnQwStfNWJuNYsTKZYeGH7r7TwGBzZLpaDaF1T/u9lck4eUEHMKYsbe/4AIA3Cl6igMM4/96ZwMxgQ2W69gNg/L3w7Wao0yToNMakxPasTOg1YRs9ZDGjY8cGHSU7tO8PKCwfH3QSY1Jme1Ym9PpHZxARZbTXI+goGalk77REBI/pNeow8pXnuL14378srxhyTnVHMyZltmdlQu+0yHS+1iZ8qYcGHSUreESY5B1J3+gcwP492GQGK1Ym1GpQSN/IHEbHegASdJysMdHrRkv5ho7yddBRjEmJFSsTaidG5lFb9jLas/NVVWli7CgATo7MDjiJMamxYmVC7fTIF+zQWkzxugQdJausoSlLvIOtWJmMYcXKhJfnMSA6nY+9bhSSH3SarDPB60avyHxqUBh0FGPKZcXKhNfaGbSQrdZkPU0meEdRU4o4LrIw6CjGlMuKlQmvhcMp1gjjvO5BJ8lKU7wu7NU8OxRoMoIVKxNeC4czTTuzDbtwbTrspibTvM70tWJlMoAVKxNOW76C9XMZFbMfAqfTBK8bXSKraM6WoKMYUyYrViacFn0EwBi7akVaTfRcE/a+kTkBJzGmbFasTDgt/BCadmKFtgw6SVabr4ewURtwctQOBZpwS1uxEpFnRGSDiMyN69ZYREaJyGL/vpHfXUTkIRFZIiKzRcS+TueyvTtgxWT7O5BqoESY4B3FSZE54HlBxzGmVOncs/ofkLi1uQMYo6odgTH+c4CzgI7+7XrgsTTmMmG3dBx4RVasqsnEWDeayA5YNyvoKMaUKm3FSlUnAN8kdB4EPOc/fg4YHNf9eXU+AxqKiB3/yVWLRkDNBtCmV9BJcsIk/7yV/SGjCbPqPmfVQlXXAvj3zf3urYBVccOt9rvtR0SuF5FpIjJt48aNaQ9rAuB5sHgEdDgNovYPNtVhEw2Y5x0KS6xYmfAKSwOLZJfT/t5/F6jqk6raU1V7NmvWrBpimWq3ZgZ8u9EOAVazCV43WDXFnS80JoSq+6vrehFpqapr/cN8G/zuq4E2ccO1BtZUczZTjRL/ELDEr/Le4OaocOxLHltfSj6MqXoTvG783HsPVkyCzmcFHceY76nuPathwFX+46uAd+O6/9hvFdgb2FZyuNDklv6R6UzXjmylXtBRcsoXXifIr23nrUxopbPp+svAp0BnEVktItcCQ4DTRWQxcLr/HOBDYBmwBHgKuDFduUx4teAbjoqsYKxdtaLaFZIPbU+CJWOCjmJMUmk7DKiql5XSa0CSYRW4KV1ZTGY4NToTgDHeMQEnyVHtB8DikbBlBTRqG3QaY/YTlgYWxjAgMoPV2pRF2jru6wPjAAAU+UlEQVToKLmpw2nufvGoYHMYk4QVKxMKNSikT2QuY2PHkLxxqEm7ph2gcXv3OzdjQsaKlQmF3pH51Ja9jLVDgMHqNBCWT4DCb4NOYsx+rFiZUOgfmc4urcGn3hFBR8ltnc6E2F5Y9nHQSYzZjxUrEwLKgOgMJntd2UtB0GFy2yEnQI363/1FizFhYcXKBK6jfE1r2WSHAMMgrwDa93fnrfR7F5ExJjBWrEzgBkSmAzAu1j3gJAZw5612roO1dhV2Ex5WrEzg+kdnMNdryzqaBB3FAHQ8HRA7FGhCxYqVCVRDdnCsLLIfAodJnabQ5nhYYNdmNOFhxcoE6pTILKKi/u+rTGh0OQ/WzYZvlgedxBjAipUJ2IDoDDZqfWZru6CjmHiHn+vuF7wfbA5jfFasTGDyKKZfZBbjYsegtiiGS+PD4KCjYP57QScxBrBiZQLUM7KI+rLLzleFVZdB7g8Zd6wLOokxVqxMcPpHZrBX85jkHRV0FJNMl/PcvR0KNCFgxcoEZkBkOlO8LnxLraCjmGSadYYmHeHLYUEnMcaKlQlGW1lL+8haxnj2R4uhJeL2rlZMgm83BZ3G5DgrViYQ/SPujxbHenbVilA78kLQGMx7O+gkJsdZsTKB6B+ZziKvFau0RdBRTFkOOhKaHwGzXws6iclxVqxM9duznV6RBYy1Q4CZodslsHqq/UDYBMqKlal+S8eSLzHG2FUrMsORF7n7OW8Em8PkNCtWpvotGsFWrcN07Rh0EpOKhm3g0D4w+1X72xATGCtWpnp5MVg8kvHe0cSIBp3GpKrbJbB5MaydGXQSk6OsWJnq9fUXsGsTY2J2viqjHDEIojVgxtCgk5gclRd0AJNjFrwPkTw+9roFncQkaHtH2X8J8q/84zht6ov0mngC84dcWE2pjHFsz8pUH1WY/z4cdjLbqRt0GnOAXiruT33ZzbnRz4KOYnKQFStTfTYugG+W7vv7CZNRpmlnFnmtuCI6OugoJgdZsTLVZ/77gMDh5wSdxFSIMDR2Gt0jy2DtrKDDmBxjxcpUn/nDoPVxUO+goJOYCno7dhK7tQCmPRN0FJNjrFiZ6rHlK/c36V3sEGAm204dhsVOhFmvwrebg45jcogVK1M9Fvgtzex8VcZ7KnY2FO+Gz58KOorJIVasTPWYPwyad4Um7YNOYippibaGTmfBlCegcFfQcUyOsGJl0m/balj5KRx5QdBJTFXpcwvs/gZmvBh0EpMjrFiZ9Cv5L6SuPwg2h6k6h54AbXrBpw9DrDjoNCYHWLEy6Tf3TTj4GDsEmG363ApbV8LsV4JOYnKAFSuTXpuXwpoZ7h9nTXbpfJb7EjJ+CBTtCTqNyXJ2bUCTXnPfcvdd7XxVNim5juCJkbN4qeBv/OWu3/BM7KzvDbdiiP0A3FQN27My6TX3TTjkRGjQOugkJg0+8Y5kYuxIbsp7h7pYy0CTPoEUKxFZISJzRGSmiEzzuzUWkVEisti/bxRENlOF1s2BjfPhSGtYkc3+UXwpTWQH1+e9H3QUk8WC3LM6VVW7q2pP//kdwBhV7QiM8Z+bTDb9BYgW2PmqLDdH2/Fu7ERuiH7AIbI+6DgmS4XpMOAg4Dn/8XPA4ACzmMoq2uP+Bv3wc6F246DTmDS7p+gKCsnjz3n/AzToOCYLBVWsFBgpIl+IyPV+txaquhbAv2+eOJKIXC8i00Rk2saNG6sxrjlgCz+EPVvhmCuDTmKqwQYa8WDxhZwancUZkWlBxzFZKKhi1UdVewBnATeJyMmpjKSqT6pqT1Xt2axZs/QmNJUz4wVo0Aba9Qs6iakm/4udyXyvDX/Kf4E67A46jskygRQrVV3j328A3gaOB9aLSEsA/35DENlMFdi6CpaOg+6XQyQadBpTTWJE+UPRtbRkM3/IGxp0HJNlqr1YiUgdEalX8hg4A5gLDAOu8ge7Cni3urOZKjLjRUBdsTI5Zbp24qnYOVyeN5ZTIvYHjabqBLFn1QKYJCKzgKnAB6r6ETAEOF1EFgOn+89NpinaA9Oeho5nQKO2QacxAXig+CIWea0Ykv8U7N4SdByTJar9Chaqugw4Okn3zcCA6s5jqtjcN+HbjdD7xqCTmIDspYBfF/2cdwr+BMN+AZe8ACJBxzIZzi63ZKqOKvPfHoLQhoFPfQt8EHQiE5C52o4hxZfxf/OHwmePwQn25cVUTph+Z2Uy3YqJdIms5JnYQMC+See6/8bOdr+zG/VHWDU16Dgmw1mxMlXn00fZrPV4N9Yn6CQmFAQG/cddF/L1n8C3m4MOZDKYFStTNdbMhEXDeb74DPZSEHQaExa1GsLFz8G3m+Ctn4LnBZ3IZCgrVqZqjPsb1GyY9G8iTI47uDucNQSWjoGJ/ww6jclQ1sDCVN6qz2HxCBjwJ3Z8UDvoNCZESv73CprzQH4fBo39G1ePKOZjb/8Gwfa/V6Y8tmdlKm/c3VC7KRx/Q9BJTGgJdxZdy0Jtw8P5D9NO1gQdyGQYK1amcpaOg2Xjoe+voUbdoNOYENtNTa4rvI1C8ngq/37qszPoSCaDWLEyFVe8Fz78rbtSRc9rgk5jMsDXNOPnhbfSRjbwcP4jRIkFHclkCCtWpuI+eQg2L4az74f8WkGnMRnicz2cPxZfwynR2fw+76Wg45gMYQ0sTMV8sxwm/BOOGAQdTws6jckwr8ZO5XBZyXV5w1mobQBrYGHKZntW5sB5Hrx/K0Ty4Mx7g05jMtTdxVcyKdaVu/OegZVTgo5jQs6KlTlwkx9wjSrOuBsatAo6jclQMaLcVHQLa7UJvHqF+x80Y0phhwFNSkp+L9NDFvFawd185PXm5jeawRt2sVpTcduoy3VFtzG66K/wyuVwzQgosN/qme+zPSuTsqZs46GCR/ham/L7ouuwi9WaqrBEW8NFT8O6OfDujaAadCQTQlasTEpqs4enC+6jCdu5ueiX7MC+/Zoq1OlMOO0umPe2a7hjTAI7DGjKFyviP/n/5khZzk+LbmOOtgs6kclGfW6BDV+6K6I0Pxy6nBd0IhMitmdlylZcCG9ex6nRWfxf8TWM9XoEnchkKxE47yFodSy8+VNYPS3oRCZErFiZ0hXuglcugy/f4e6iK3g5NiDoRCbb5deEy16Bei3gpUtg05KgE5mQsGJlktuxHl4YDEvGwHkP8d+Y/WjTVJO6zeHKt9zjF38AOzcEm8eEghUr832rpsKTp8Da2XDx/+DYq4JOZHJNk/Zw+evw7UYYehHs3RF0IhMwK1ZmHy8Gkx6EZ8+GvBpw3WjoOjjoVCZXtT7WfVlaNxdeuwpiRUEnMgGy1oDG2bwU3vk5rJriWmGd9xDUbhx0KpMj9v1J4/ddEr2Gfyx9infvOp9Bf34fItFqTGbCwopVrivcBZP+BZP/7a6c/oP/wlEXuZZZxoTAa7FTacwO7sh/Bd7+GVzwuBWsHGTFKld5Hsx5Hcb+Fbatgm4/hNP/AvUOCjqZMd/zeOx8InjcPuc190Vq0KMQtc1XLrFPO9eowqIR7oeX6+ZAy6PhgiegbZ+gkxlTpkdjg7n9zM7uC9bene4STfY/ajnDilWu8DyYP8xdymb9HGh4KFz4NIcNrYE+vhWwC9KaDHDyb6BGfRh+O7x4IVz6EtRqGHQqUw2sWGWh+JPVeRRzdmQKN+e9Q6fI1yz1WvJo8c94d92JFA+1j99kFrdst+L8yE38c8VjrLq3N9cW/YYV2nK/4VYMsd8FZhvbWmWpRmznsug4rswbxcHyDQu91vyi8GY+8Hrj2S8WTIYb5p3IusJGPFbwIO8U/Imbi37JJO+ooGOZNLJilW3WzWVI3pMMjk6mphQxKdaVP8auZqx3DGpFymSRqdqFQYV/5b/59/N8/hAeiQ3i38UXEsNaCmYjK1bZYPdWmPsmzHgR1kxnULSAt2J9eTY2kMXaOuh0xqTNam3OBYV/5s95z/HLvHfoFVnAb4puCDqWSQMrVpkqVgwrJsDMl13DieI90LwrnPk3er/blG3UDTqhMdViNzW5vfgGPvG68tf8ZxlRcAd8theOv95+j5VFrFhlklgxfDXZ/UHd/GGwazPUbADHXOluLbuDCNvetZZ9Jve8453ElL1duCf/afp/dAfMegUG3guHnhh0NFMFrFiFXeEuWP4xLBzOxi/eppls51utwWjvWD6I9eLjPUezd2IBTFwDrAk6rTGBWksTrin6LSsu2QOj/gTPnuUuH3by7dCyW9DxTCVYsQqBxOuiNWMLA6IzGBCZzkmRudSSQnZoLaZ4R/N+rDfjvaPZQ42A0hoTduIuGdb5bPjkYfj0EZj/HnQ8E3pdD+36Q8QaG2UaK1YhkE8x3WUJfaJz6ReZSffIMgBWa1NejfVjtHcsU7wuFNnHZUzqCmpDv99Brxvg86fgs8fdD4kbtIFul8Dh58DBPew6mBnCtn5B2LsT1s+FlZ/B8gnMqjGJ2rIXT4VZ2p77ii5htNeDhdoGsBXJmEqp1RBO/i2c+EtY+CFMf979Fc7E+6HuQdD2JDj0BGh5DDTtCDXrB53YJBGqYiUiA4F/A1Hgv6o6JOBIlbN3B2z5CrauhM1LYO0sd9u8BFA3TLPDeS12Cp94XfnM68J2a8VnTHrk1YCuF7jbrm/cNTIXj4QVk2DuG/uGq98KmnZytwatXEGr18Ld123uGjVZK8NqF5piJSJR4D/A6cBq4HMRGaaqX6Z94qrujwc15t97+z/2iqFol2vsULQLCnfu/3j3Fvh2s2udt2uT+3fTrStd93j1W7sLxx51sbs/+Bio14K7yvgvH2PMgSvr/7H2qQ9cBFxIG9lAF1lJB1lD+y1f02HrV7Rf+il1ZU/yUfPruD2wGvX33deoBwV1IL+2u8Bu4n1B7VL61YZogTuPJlGI5LliKFH/3o6uQIiKFXA8sERVlwGIyCvAIKBqi9WGBfBkP1eM1HMFqWQvpzLy60CdJlC7CdRpDq2OZcinu1ilzVmtTflKW7B1Tz3YAMwCiAHTKj9dY0wlCau0Bau0BSPBrZoAKHXZTXPZ6m5spalsox67qFe8i3q7d1NPdrnn8g312EX7hlH3JbZol/vtYxXl+66AlZwW6HQGXPJ8Fb1+ZhDVKthQVwERuQgYqKrX+c9/BPRS1ZvjhrkeuN5/2hlYWEWTbwpsqqLXSpewZwx7PrCMVSXsGcOeD4LNeKiqNgto2hUWpj2rZPu6+1VSVX0SeLLKJywyTVV7VvXrVqWwZwx7PrCMVSXsGcOeDzIjY9iE6ccGq4E2cc9bY79yNcYYQ7iK1edARxE5TEQKgEuBYQFnMsYYEwKhOQyoqsUicjMwAtd0/RlVnVdNk6/yQ4tpEPaMYc8HlrGqhD1j2PNBZmQMldA0sDDGGGNKE6bDgMYYY0xSVqyMMcaEXlYXKxFpLCKjRGSxf9+olOGu8odZLCJXxXU/VkTmiMgSEXlIZN9PyUXkFyKyUETmicg/wpbP7/8bEVERaVqRfOnMKCL3icgCEZktIm+LSMMKZBvofwZLROSOJP1riMirfv8pItI2rt/v/e4LReTMVF8zyHwi0kZExonIfH+5u6Uy+dKRMa5fVERmiMj7YcwoIg1F5A1/GZwvIieEMOOv/M95roi8LCI1K5Mx46lq1t6AfwB3+I/vAP6eZJjGwDL/vpH/uJHfbypwAu43YMOBs/zupwKjgRr+8+Zhyuf3a4NrrPIV0DSE8/AMIM9//Pdkr1tOriiwFGgHFOCuC3JEwjA3Ao/7jy8FXvUfH+EPXwM4zH+daCqvGXC+lkAPf5h6wKKK5ktXxrjxfg28BLxfyXU4LRmB54Dr/McFQMMwZQRaAcuBWv5wrwE/qcy8zPRbVu9Z4S7X9Jz/+DlgcJJhzgRGqeo3qroFGAUMFJGWQH1V/VTd0vJ83Pg/B4ao6l4AVd0QsnwADwC3U/lrSaUlo6qOVNVif/zPcL+rOxDfXZ5LVQuBkstzlZb9DWCAv2c3CHhFVfeq6nJgif96qbxmYPlUda2qTgdQ1R3AfNxGraLSMQ8RkdbAOcB/K5EtbRlFpD5wMvA0gKoWqurWMGX0h8sDaolIHlCbHP/dabYXqxaquhbAv2+eZJhWwKq456v9bq38x4ndAToBff3d+Y9F5Lgw5ROR84GvVXVWBXOlPWOCa3B7XQeitGkmHcYvjNuAJuXkLe81g8z3Hf8w0jHAlArmS2fGB3FflLxKZEtnxnbARuBZ/1Dlf0WkTpgyqurXwD+BlcBaYJuqjqxExowXmt9ZVZSIjAYOStLrD6m+RJJuWkZ3cPOtEdAbOA54TUTa+XsPgeYTkdr+a5+R4usHNQ9Lpv0HoBgYmuK0yptmZXIl+/JW0T3TdORzI4nUBd4EblXV7RXMV9b0UxmmtOXvXGCDqn4hIv0qka286acyTGnd84AewC9UdYqI/Bt3iPuPYcko7tzwINyhwa3A6yJypaq+WMGMGS/ji5WqnlZaPxFZLyItVXWtf0gq2eG61UC/uOetgfF+99YJ3dfEjfOWX5ymioiHuzDlxhDka49bwGe5owy0BqaLyPGqui5ZjoDmIeIaYpwLDEhW6MuRyuW5SoZZ7R9KaQB8U864VXXJr7TkE5F8XKEaqqpvVTBbOjOeD5wvImcDNYH6IvKiql4ZooyrgdWqWrJX+gauWFVUOjKeBixX1Y0AIvIWcCKQs8Uq8JNm6bwB97F/44B/JBmmMe5EZiP/thxo7Pf7HLf3VNI44Gy/+8+Av/iPO+F24yUs+RLGX0HlGlikax4OxP39S7MK5srDNeQ4jH0ntbsmDHMT+5/Ufs1/3JX9T2ovw53ULvc1A84nuPN+D1bR+lHlGRPG7UflG1ikJSMwEejsP74LuC9MGYFewDzcuSrBne/6RVV87pl6CzxAWt+cOyY8Bljs35dsQHvi/om4ZLhrcCc2lwBXx3XvCczFtdB5hH1X/CjAfcOZC0wH+ocpX8I0VlC5YpWuebgEV+Rn+rfHK5DtbFyLuKXAH/xufwHO9x/XBF73pzUVaBc37h/88RayfyvK771mJeZdleYDTsIdOpodN9++9wUl6HkY178flSxWafycu+P+UG428A5+69WQZfwzsAC3/ryA3/o4V292uSVjjDGhl+2tAY0xxmQBK1bGGGNCz4qVMcaY0LNiZYwxJvSsWBljjAk9K1bGGGNCz4qVMcaY0Pv/mkm+J7Z3YWQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "###\n",
    "# Synthetic data 2 - with unobservable Z\n",
    "\n",
    "# R ~ U(0,1)\n",
    "# X ~ U(0,1)\n",
    "# Z ~ U(0,1)\n",
    "# T ~ Bin(1, sigmoid(a * r + b * x + e * z))\n",
    "# Y ~ Bin(1, sigmoid(c * t + d * x + f * z))\n",
    "\n",
    "# Weights for edges:\n",
    "# R -> T: a\n",
    "# X -> T: b\n",
    "# T -> Y: c\n",
    "# X -> Y: d\n",
    "# Z -> T: e\n",
    "# Z -> Y: f\n",
    "\n",
    "N = 10000\n",
    "a = 1\n",
    "b = 1\n",
    "c = 1\n",
    "d = 1\n",
    "e = 1\n",
    "f = 1\n",
    "\n",
    "\n",
    "def generateDataWithZ(N=N, a=a, b=b, c=c, d=d, e=e, f=f):\n",
    "\n",
    "    r = npr.uniform(size=N)\n",
    "    x = npr.uniform(size=N)\n",
    "    z = npr.uniform(size=N)\n",
    "\n",
    "    t = npr.binomial(n=1, p=sigmoid(a * r + b * x + e * z), size=N)\n",
    "\n",
    "    y = npr.binomial(n=1, p=sigmoid(c * t + d * x + f * z), size=N)\n",
    "\n",
    "    # If decision is negative, we set Y to positive.\n",
    "    y[t == 0] = 1\n",
    "\n",
    "    return r, x, z, t, y\n",
    "\n",
    "\n",
    "def analytic_R_on_Y(r, a, b, c, d, e, f):\n",
    "    '''Return probability of Y=0 with given parameters.'''\n",
    "\n",
    "    # P(y=0|T=t, x, z) * p(T=t|r, x, z) * p(x) * p(z)\n",
    "    def integrand(t, x, z):\n",
    "        p_y0 = 1 - sigmoid(c * t + d * x + f * z)\n",
    "\n",
    "        if t == 1:\n",
    "            p_t = sigmoid(a * r + b * x + e * z)\n",
    "        elif t == 0:\n",
    "            p_t = 1 - sigmoid(a * r + b * x + e * z)\n",
    "\n",
    "        p_x = 1\n",
    "        \n",
    "        p_z = 1\n",
    "\n",
    "        return p_y0 * p_t * p_x * p_z\n",
    "\n",
    "    t_1 = si.dblquad(lambda x, z: integrand(1, x, z), 0,\n",
    "                     1, lambda x: 0, lambda x: 1)\n",
    "\n",
    "    # P(y=0|t=0, x) * p(t=0|r, x) * p(x)\n",
    "    # Now equal to 0, hence commented out. See data generation.\n",
    "    t_0 = si.dblquad(lambda x, z: integrand(0, x, z), 0,\n",
    "                     1, lambda x: 0, lambda x: 1)\n",
    "\n",
    "    return t_1[0]  # + t_0[0]\n",
    "\n",
    "\n",
    "diffs = np.zeros(0)\n",
    "\n",
    "for i in range(2000):\n",
    "\n",
    "    r, x, z, t, y = generateDataWithZ()\n",
    "\n",
    "    # Fit predictive models.\n",
    "    lr_t = LogisticRegression(solver='lbfgs')\n",
    "    lr_y = LogisticRegression(solver='lbfgs')\n",
    "\n",
    "    lr_t = lr_t.fit(np.array([r, x]).T, t)\n",
    "    lr_y = lr_y.fit(np.array([t, x]).T, y)\n",
    "\n",
    "    def causal_effect_of_R_on_Y(r):\n",
    "\n",
    "        # Integrate P(Y=0|T=1, X=x) * P(T=1|R=r, X=x) * f(x) from 0 to 1\n",
    "        return (si.quad(\n",
    "            lambda x: lr_y.predict_proba(np.array([[1, x]]))[0, 0] * lr_t.\n",
    "            predict_proba(np.array([[r, x]]))[0, 1], 0, 1))\n",
    "\n",
    "    r0 = causal_effect_of_R_on_Y(0)\n",
    "    r1 = causal_effect_of_R_on_Y(1)\n",
    "\n",
    "    analytical = analytic_R_on_Y(1, a, b, c, d, e, f) - analytic_R_on_Y(\n",
    "        0, a, b, c, d, e, f)\n",
    "\n",
    "    diffs = np.append(diffs, analytical - r1[0] + r0[0])\n",
    "\n",
    "    if i % 200 == 0:\n",
    "        print(i / 20, \"%\", end=\" \")\n",
    "\n",
    "print()\n",
    "print(\"Analytical:\", analytical)\n",
    "print(\"Estimated: \", r1[0] - r0[0])\n",
    "print()\n",
    "print(\"Values for P(y=0|do(r=1)) and P(y=0|do(r=0))\\n\")\n",
    "print(\"Analytical:\", analytic_R_on_Y(1, a, b, c, d, e, f),\n",
    "      analytic_R_on_Y(0, a, b, c, d, e, f))\n",
    "print(\"Estimated: \", r1[0], r0[0])\n",
    "print()\n",
    "print(\"Average difference:\", np.mean(diffs))\n",
    "print(\"Std of differences:\", np.std(diffs))\n",
    "plt.hist(diffs, bins=25)\n",
    "plt.title(\"Histogram of differences over 2000 simulations (analytic - estimate)\")\n",
    "sns.kdeplot(diffs)\n",
    "plt.axvline(x=0, c='r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now it is visible that the quantity $P(Y=0|do(R=r))$ is slightly biased when estimated with equation \\ref{int}.\n",
    "\n",
    "\n",
    "According to [Pearl's Theorem 3](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2836213/): \"A *sufficient* condition for identifying the causal effect P(y|do(x)) is that every path between X and any of its children traces at least one arrow emanating from a measured variable.\""
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}