docs/source/examples/causal_inference.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Causal Inference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda/lib/python2.7/site-packages/sklearn/cross_validation.py:41: DeprecationWarning: This module was deprecated in version 0.18 in favor of the model_selection module into which all the refactored classes and functions are moved. Also note that the interface of the new CV iterators are different from that of this module. This module will be removed in 0.20.\n",
      "  \"This module will be removed in 0.20.\", DeprecationWarning)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "from fklearn.training.regression import xgb_regression_learner\n",
    "from fklearn.training.classification import xgb_classification_learner, logistic_classification_learner\n",
    "from fklearn.training.causal_inference import IPTW_learner\n",
    "from fklearn.validation.evaluators import r2_evaluator, roc_auc_evaluator\n",
    "from fklearn.data.datasets import make_confounded_data\n",
    "\n",
    "from sklearn.model_selection import train_test_split\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "N=50000\n",
    "features = [\"sex\", \"age\", \"severity\"]\n",
    "treatment = [\"medication\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Synthetic Data Model\n",
    "\n",
    "To better understand counterfactual inference, we will use a synthetic data, for which the generating process is known. More precisely, the data which is generated from the following model\n",
    "\n",
    "\n",
    "$Sex \\sim  \\mathcal{B}(0.5)$  \n",
    "  \n",
    "$Age \\sim \\operatorname{Gamma}(8, 4)$\n",
    "  \n",
    "$Severity \\sim \\mathbb{1}_{\\{age < 30\\}} \\operatorname{Beta}(1,3) + \\mathbb{1}_{\\{age \\geq 30\\}} \\operatorname{Beta}(3,1.5)$  \n",
    "  \n",
    "$Treat \\sim \\mathcal{B}(-0.8 + 0.33*sex + 1.5*severity)$  \n",
    "  \n",
    "$Recovery = \\operatorname{Poisson}(2+0.5 * sex+0.03 * age+ 1*severity-1 * medication)$  \n",
    "  \n",
    "We can see that the linear factor that multiply the medication is $-1$. Since this relation is inside a Poisson model, we can say that the expected causal effect of medication on recovery time is $exp(-1)\\approx 0.36$. We can also see that there are lots of confounding in the above model. For instance, Severity directly influences both the probability of treatment and the recovery time directly. The same if true for sex, with the additional complication that sext also influences on the severity of the daisies. Finally, age directly influences the recovery time and affect treatment only indirectly through influence in severity. Moreover, there are some non trivial data, like the severity, which is a mixture of two Beta distribution.\n",
    "\n",
    "The function bellow generates such data. It returns a Data Frame for which the treatment is randomly assigned, one DataFrame that simulates a observational dataset with the a above treatment probability. It also returns a third DataFrame, just like the simulated observational one, but with the treatment inverted. This will be the counterfactual dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda/lib/python2.7/site-packages/pandas/core/computation/check.py:17: UserWarning: The installed version of numexpr 2.4.3 is not supported in pandas and will be not be used\n",
      "The minimum supported version is 2.4.6\n",
      "\n",
      "  ver=ver, min_ver=_MIN_NUMEXPR_VERSION), UserWarning)\n"
     ]
    }
   ],
   "source": [
    "df_rnd, df_obs, df_cf = make_confounded_data(N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets now see the correlations present in this dataframe. As you can see, all variables are correlated with both recovery and treatment assignment, which is the typical confounding scenario."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sex</th>\n",
       "      <th>age</th>\n",
       "      <th>severity</th>\n",
       "      <th>medication</th>\n",
       "      <th>recovery</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>sex</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>-0.005599</td>\n",
       "      <td>-0.003693</td>\n",
       "      <td>0.208767</td>\n",
       "      <td>0.264228</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>age</th>\n",
       "      <td>-0.005599</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.566635</td>\n",
       "      <td>0.463141</td>\n",
       "      <td>0.601405</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>severity</th>\n",
       "      <td>-0.003693</td>\n",
       "      <td>0.566635</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.792594</td>\n",
       "      <td>0.177204</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>medication</th>\n",
       "      <td>0.208767</td>\n",
       "      <td>0.463141</td>\n",
       "      <td>0.792594</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>-0.097037</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>recovery</th>\n",
       "      <td>0.264228</td>\n",
       "      <td>0.601405</td>\n",
       "      <td>0.177204</td>\n",
       "      <td>-0.097037</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                 sex       age  severity  medication  recovery\n",
       "sex         1.000000 -0.005599 -0.003693    0.208767  0.264228\n",
       "age        -0.005599  1.000000  0.566635    0.463141  0.601405\n",
       "severity   -0.003693  0.566635  1.000000    0.792594  0.177204\n",
       "medication  0.208767  0.463141  0.792594    1.000000 -0.097037\n",
       "recovery    0.264228  0.601405  0.177204   -0.097037  1.000000"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_obs.corr()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Naive Model\n",
    "\n",
    "Lets first try to fit a model directly to the responde variable (recovery) without much attention to confounding. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "regr_learner = xgb_regression_learner(features=features + treatment, target=\"recovery\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we fit it to the the data where the treatment is randomly assigned."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "sem_fn, df_rnd, _ = regr_learner(df_rnd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get the average causal effect, we create two dataframes, one for which everyone receives the treatment and one for which no one receives it. We then run the model for both datasets and divide their predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.369575560092926"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_neg = df_rnd.copy()\n",
    "X_pos = df_rnd.copy()\n",
    "\n",
    "X_neg[\"medication\"] = 0\n",
    "X_pos[\"medication\"] = 1\n",
    "\n",
    "preds_rnd = sem_fn(X_pos)[\"prediction\"] / sem_fn(X_neg)[\"prediction\"]\n",
    "preds_rnd.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see it is trivial to identify the causal effect if the treatment is randomly assigned. Bellow, we plot the predicted treatment effect distribution with the true average treatment effect and the predicted one"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3XucXHV9//HXZ2Znd2av2Ru7STYhmxAIBDCaEFSwpV7xUmmFtmK1au2DekHKr/p41PZXK735oy3VtlpUailoSW2tbUUqKqVEBSEQIEBCAiHJEjbZZe/321y+vz/mTBiW3ezs7sycmZ338/GYR845c86cz5zMzGe/53sz5xwiIlK6An4HICIi/lIiEBEpcUoEIiIlTolARKTEKRGIiJQ4JQIRkRKnRCAiUuKUCERESpwSgYhIiSvzO4BMNDU1uQ0bNvgdhpSiZ55hfDrG0ca2l22+YG2dTwGJZO7RRx/tc841L7RfUSSCDRs2sHfvXr/DkFJ02WU8dLSf977vxpdt3nvjO30KSCRzZvZ8Jvvp1pCISIlTIhARKXFKBCIiJa4o6ghERBYSjUbp7OxkamrK71DyLhwO09bWRigUWtLxSgQisiJ0dnZSU1PDhg0bMDO/w8kb5xz9/f10dnbS3t6+pNfQrSERWRGmpqZobGwsqSQAYGY0NjYuqySkRCAiK0apJYGU5b5vJQIRkRKnRCAyj28+2MHJoUm/wxDJOVUWi8yya89xZmIJ/vz7T7N5YMLvcERyTolAZA6HukeIxh014RCjU1G/w5FF2PCZ/875OTrmGGKko6ODyy+/nNe+9rX87Gc/46KLLuLDH/4wn/vc5+jp6eGOO+5g69atfPKTn2T//v1Eo1FuuOEGrrjiCjo6OvjABz7A+Pg4AF/+8pd5/etfz+7du7nhhhtoampi//79bN++nX/+53/Oel2IEoHIHJ7sHKYmXEZTdbkSgWTsueee49vf/ja33norF110Ebt27eL+++/nzjvv5POf/zznnXceb3zjG7n11lsZGhpi586dvPnNb+aMM87gnnvuIRwOc/jwYa6++upT46s9/vjjHDhwgDVr1nDJJZfwwAMPcOmll2Y1biUCkVmmonGefXGUne0NQGm2QpGlaW9v54ILLgBg69atvOlNb8LMuOCCC+jo6KCzs5M777yTm266CUg2eT1+/Dhr1qzh2muvZd++fQSDQZ599tlTr7lz507a2pKj327bto2Ojg4lApFcO9g1QizhuFBDTcsiVVRUnFoOBAKn1gOBALFYjGAwyHe+8x3OOeeclx13ww030NLSwhNPPEEikSAcDs/5msFgkFgslvW41WpIZJYnO4dZFQmxrqHS71BkhXnb297Gl770JZxzQPK2D8Dw8DCrV68mEAjwzW9+k3g8nte4VCIQSTM8EeVwzyiXbGoq2c5JxW6uitxC8dnPfpbrr7+eCy+8kEQiQXt7O3fddRcf//jHufLKK/nGN77B5ZdfTlVVVV7jslRmKmQ7duxwmphG8mH3Mz186J8e4bcubWdjczVv+tivcaxv7BUT0xTyj02pOnjwIOeee67fYfhmrvdvZo8653YsdKxuDYmkOdQ9CsDqushp99u15zi79hzPR0giOadEIJLmUNcIdZEQkfLgafdTk1JZSZQIRNIc6h6ltTa84H7H1eNYVhAlAhHPTCzBcz1jtNYpEUhpUSIQ8RzpHSOWcEoEUnKUCEQ8h7pHADK6NXRicJJYIpHrkETyQolAxHOoa5TyYICm6ooF940lHF1DpTc3rsxvaGiIm2++edHH3XbbbZw8efLU+oYNG+jr68tmaAvKWSIws7CZPWxmT5jZATP7Y297u5ntMbPnzOxfzaw8VzGILMbB7lE2t1QTDGTWkUy3hyTdfIlgoSEhZicCP+SyZ/E08Ebn3JiZhYD7zexu4HeBLzrnvmVmXwU+Anwlh3GIZORQ1whv2Nyc0b51kZASgbzMZz7zGY4cOcK2bdsIhUKEw2Hq6+s5dOgQP/rRj3jXu97F/v37AbjpppsYGxvj/PPPZ+/evfz6r/86kUiEBx98EIAvfelLfO973yMajfLtb3+bLVu25DT2nCUCl+yyPOathryHA94IvM/bfjtwA0oE4rP+sWl6Rqc5d3VNRvuvb6jkhUElgoJ1/fWwb192X3PbNvibv5n36RtvvJH9+/ezb98+du/ezTvf+U72799Pe3s7HR0dcx5z1VVX8eUvf5mbbrqJHTte6gDc1NTEY489xs0338xNN93E17/+9ey+l1lyWkdgZkEz2wf0APcAR4Ah51yqrNQJrM1lDCKZeMbrUbyltTaj/VtqwwxNRJmYyf5IkLIy7Ny5k/b29iUd+573vAeA7du3z5tEsimng8455+LANjNbBfwnkHH5xsyuAa4BWL9+fW4CFPEcTCWC1TUZ3fJprklWKB/rG2frGg1XXXBO85d7vqQPHFdWVkYirZXZ1NTpGxqkhp7O1bDTs+Wl1ZBzbgi4D3gdsMrMUgmoDTgxzzG3OOd2OOd2NDdndt9WZKkOdY3QVF2RUYshgKbqZBuHo73juQxLikhNTQ2jo6NzPtfS0kJPTw/9/f1MT09z1113ZXRcvuSsRGBmzUDUOTdkZhHgLcBfkEwIVwHfAj4IfDdXMYhk6lD3aMb1AwBN1RUYSgTyksbGRi655BLOP/98IpEILS0tp54LhUL80R/9ETt37mTt2rUvq/z90Ic+xEc/+tGXVRbnWy5vDa0GbjezIMmSx7855+4ys6eBb5nZnwGPA/+YwxhEFhSLJ3j2xVF+43VnZnxMKBigrjLE0b6xhXeWkrFr1655n7vuuuu47rrrXrH9yiuv5Morrzy1nl4nsGPHDnbv3p3NEOeUy1ZDTwKvnmP7UWBnrs4rslgd/RNMxxIZVxSnNFVXqEQgK4J6FkvJSw0tsWURt4YglQjGKIbJnUROR4lASt6hrlGCAeOsM6oXdVxzdTnjM3F6RqdzFJksVqkm5eW+b81ZLCUrNcPY/xx8kcaqcr7z6JwN2ObV5DUhPdI7RksGA9VJboXDYfr7+2lsbCyp+aadc/T39xMOL/0zqEQgJa97ZIr1DZWLPq7Za2p6tHec129qynZYskhtbW10dnbS29vrdyh5Fw6HaWtrW/LxSgRS0qaicYYmoly8YfF/TdVGQoRDAVUYF4hQKLTknrylTnUEUtK6h5M9PDOZjGa2gBntTdVqQipFT4lASlr3SCoRRJZ0/MbmKpUIpOgpEUhJ6x6eIhIKUhte2l3STU1VdA5OMB2LZzkykfxRIpCS1j0yRWtdeMmtTDY2V5NwcLxfQ1JL8VIikJLlnKN7ZGpZTT83tyT7HjzZOZytsETyTolAStbIVIyZWIIzajIbcXQu57bW0lob5gcHurMYmUh+KRFIyer1egQ3LyMRBALG5ee38uNnexmb1iQ1UpyUCKRk9Y55iSDDOQjm844LVjMTS/C/h3qyEZZI3ikRSMnqHZ2moixAzRJbDKVsP7Oe5poK7n6qK0uRieSXEoGUrL7R6eQEM8sclyYYMC7f2sruZ3o1h7EUJSUCKVm9Y9PLqh9I9/YLWpmMxvnxM6U3zo0UPyUCKUkTMzGGJ6MZz1G8kJ0bGmisKufu/Wo9JMVHiUBKUmpYiGyVCMqCAS7d3MTPjvSX7Jj4UryUCKQkHelNDhSXrUQA8NqNjfSNTXO0T2MPSXHRMNRSko72jmNAY1X5sl4nNbkNJCufAf7u3sP87XtfMV23SMFSiUBK0pHeMeqrygkFs/cVaKwupyZcxjGVCKTIKBFISTraO05T9fJKA7OZGe1NVRzrG1c9gRSVnCUCM1tnZveZ2dNmdsDMfsfbfoOZnTCzfd7jHbmKQWQuiYTjaN/YsnsUz6W9qYrRqZhKBVJUcllHEAM+5Zx7zMxqgEfN7B7vuS86527K4blF5tU1MsVUNHFq8vls2tiUHI10z7EBNjZXZ/31RXIhZyUC51yXc+4xb3kUOAiszdX5RDJ1pCf7LYZSmqrLqa4o46Gj/Vl/bZFcyUsdgZltAF4N7PE2XWtmT5rZrWZWn48YRFKeH0hOItNYlf1EkKoneOio+hNI8ch5IjCzauA7wPXOuRHgK8AmYBvQBfz1PMddY2Z7zWxvb6+67Uv2dA9PUhawZQ82N58NjZW8ODJ9aj5kkUKX00RgZiGSSeAO59x/ADjnXnTOxZ1zCeAfgJ1zHeucu8U5t8M5t6O5uTmXYUqJ6RpOzkoWWOZgc/Np9CqhNX2lFItcthoy4B+Bg865L6RtX5222y8D+3MVg8hcuoeT8xTnSkNlslnqC4OTOTuHSDblstXQJcAHgKfMbJ+37Q+Aq81sG+CADuC3cxiDyCt0D09x7pranL3+qsoQZvDCgEoEUhxylgicc/cDc5W9v5+rc4osxDlH1/AUb9xyRs7OURYM0FIT5oVBJQIpDupZLCVlZDLGZDSe01tDAOsaInQO6NaQFAclAikpXSPJH+fVdZGcnmddfaVKBFI0lAikpHQNJ5t05rpE0NZQSffIFNOxeE7PI5INSgRSUrq9RLA617eG6iM4ByeH1JdACp8SgZSUruEpApab4SXSrWuoBNRySIqDEoGUlO7hSZprKrI6D8FcTiUC1RNIEVAikJLSNTxFa44rigFaa8OEgsYLajkkRUCJQEpK9/AUq2tzWz8AEAwYa1ZFVCKQoqBEICUl18NLpFtXX0mn6gikCCgRSMkYnYoyOh3LeYuhlHUNEY03JEVBiUBKxosj+elDkNJWX8nA+Azj07G8nE9kqZQIpGR0nepDkPvKYlDLISkeSgRSMrry1JksZV19MuGo5ZAUOiUCKRmpXsVn1Oa2M1mKOpVJsVAikJLRNTxFU3U5FWXBvJyvsaqcirIAXcMqEUhhy+XENCIFpXt4Mm8Vxbv2HAegqqKMPccGTq2/7+L1eTm/yGKoRCAl48TQZN4qilNqwyGGJ6N5PafIYikRSElwztE5OMm6+sq8nrcuUsaIEoEUOCUCKQkD4zNMzMRZ15DnEkEkxMhUDOdcXs8rshiqI5AVb9ee46da7hztHT91vz4f6iIh4gnH+Eyc6gp93aQwqUQgJWFwYgaA+sryvJ63NhwC0O0hKWhKBFISBse9RFAVyut5ayNKBFL4lAikJAxMRKksD+atD0FKnZcIhqeUCKRw5SwRmNk6M7vPzJ42swNm9jve9gYzu8fMDnv/1ucqBpGUoYkZGqrye1sIoLqiDANGJjXwnBSuXJYIYsCnnHPnAa8FPmFm5wGfAe51zm0G7vXWRXJqYHwm7/UDkJygpiasJqRS2HKWCJxzXc65x7zlUeAgsBa4Arjd2+124JdyFYMIQMI5hiaj1Ffmt34gJdmEVIlAClde6gjMbAPwamAP0OKc6/Ke6gZa5jnmGjPba2Z7e3t78xGmrFCjUzHiCUe9D7eGQL2LpfBllAjM7BfNbElJw8yqge8A1zvnRtKfc8leNnP2tHHO3eKc2+Gc29Hc3LyUU4sAydtCkP+moykqEUihy/TH/deAw2b2l2a2JdMXN7MQySRwh3PuP7zNL5rZau/51UDPYgIWWawhrw9Bg0+JoC5cxlQ0wUws4cv5RRaSUSJwzr2f5K2dI8BtZvagd+umZr5jzMyAfwQOOue+kPbUncAHveUPAt9dUuQiGRrwEsEqH+sIQH0JpHBlfLvHu63z78C3gNXALwOPmdkn5znkEuADwBvNbJ/3eAdwI/AWMzsMvNlbF8mZwfEoteEyyoL+dJupVV8CKXAZDX5iZlcAHwLOAr4B7HTO9ZhZJfA08KXZxzjn7gdsnpd805KiFVmCwQl/mo6m1KlEIAUu01Gw3gN80Tn3k/SNzrkJM/tI9sMSyZ7B8Rk2NFX5dn6NNySFLtOycvfsJGBmfwHgnLs361GJZEk0nmB4MupriaC8LEA4FNCtISlYmSaCt8yx7e3ZDEQkF04OTeKAhjwPNjdbXSSkYSakYJ321pCZfQz4OLDJzJ5Me6oGeCCXgYlkw7G+cQAaqyp8jaM2rL4EUrgWqiPYBdwN/D9ePibQqHNuIGdRiWTJkd5kImiq8TcR1EVCdA9P+RqDyHwWujXknHMdwCeA0bQHZtaQ29BElu9o7xiRUJCq8vwOPz3bqsoQo9MxpqJxX+MQmUsmJYJ3AY+SHAoivTmoAzbmKC6RrDjSO0ZzTQXJ/o3+SVVWnxyaZGNzta+xiMx22kTgnHuX9297fsIRya4jveOsr6/0OwxWeYnghBKBFKBMB527xMyqvOX3m9kXzGx9bkMTWZ6RqSi9o9M0+1w/AJwaArtzcNLnSEReKdPmo18BJszsVcCnSI459M2cRSWSBUdTFcXV/ieCmnCIgMEJJQIpQJkmgpg3ZPQVwJedc39PsgmpSME60jMGUBAlgmDAqIuE6Byc8DsUkVfIdIiJUTP7feD9wM95cxP420NHZAFH+8YoC5gvcxXPZVVluW4NSUFazHwE08BHnHPdQBvwVzmLSiQLjvSMs76xkmDA3xZDKfWVIU4MKRFI4cmoROD9+H8hbf04yVFIRQrW0b4xNhVQC51VleU8/sIQM7EE5WX+DIktMpdMWw29x8wOm9mwmY2Y2aiZjSx8pIg/YvEEHX0TbGz2b9TR2eorQziHehhLwcn0z5K/BN7tnKtzztU652qcc7W5DExkOToHJ5mJJwquRACowlgKTqaJ4EXn3MGcRiKSRUf7ki2GNhVUicBLBKonkAKTaauhvWb2r8B/kaw0BiBtQnqRgnKkJ9mHYGNTNc90j/kcTVJtpIyAqVOZFJ5ME0EtMAG8NW2bA5QIpCAd7Rujoaqc+gJpOgpQFgjQUhvWrSEpOJm2GvpwrgMRyaZjfeNsaPR/jKHZ2uoj6l0sBSfTVkNnm9m9ZrbfW7/QzP4wt6GJLN3z/RO+zlM8n7WrIro1JAUn08rifwB+H4gCOOeeBN6bq6BElmNyJk7X8BQbGgsvEbTVV9I9MkUsnvA7FJFTMk0Elc65h2dtO+0ErGZ2q5n1pEoR3rYbzOyEme3zHu9YbMAiC3l+IFlRXIglgrb6CPGEo3tEfQmkcGSaCPrMbBPJCmLM7Cqga4FjbgMun2P7F51z27zH9zOOVCRDHd48xe0FWCJYWx8B4IUB3R6SwpFpq6FPALcAW8zsBHAM+PXTHeCc+4mZbVhWdCKLsGvPcQB+8mwvAI90DPDUiWE/Q3qFVAe353rHeN2mRp+jEUk6bSIws99NW/0+cB/JUsQ4cCVp4w8twrVm9hvAXuBTzrnBJbyGyLz6xqapKg8SDvk7T/FcVteFqQ2XcahLI7RI4Vjo1lCN99gBfAyoB1YBHwVes4TzfQXYBGwjeWvpr+fb0cyuMbO9Zra3t7d3CaeSUtU/PkNjAUxGMxczY8vqWg51j/odisgpp00Ezrk/ds79Mclhp1/jnPu0c+5TwHZg0VNVOudedM7FnXMJki2Rdp5m31ucczucczuam5sXeyopYf1j0zRVF05HstnOba3hme5REgnndygiQOaVxS3ATNr6jLdtUcxsddrqLwP759tXZClmYglGpmIFWyIA2LK6lrHpmOYmkIKRaWXxN4CHzew/vfVfItkqaF5m9i/AZUCTmXUCnwMuM7NtJFsfdQC/vfiQRebXP54cCquxgIaWmG1La3KW14NdI6xrKLzez1J6Mh1i4s/N7G7gDd6mDzvnHl/gmKvn2PyPi4xPZFH6x5IF10IuEZzdUoMZHOoe5a1bW/0ORyTjEgHOuceAx3IYi8iy9Y8lSwRNBVwiqKoo48yGSg51q+WQFAbNlycrSv/4DNUVZVQUYNPRdFtaaznUpZZDUhgyLhGIFIO+sZmCrh9IdXqLxhMc6xvntgc6KC8L8L6LF90ITyRrVCKQFaV/fLqg6wdSWuvCOKBnVGMOif+UCGTFmI7FGZ2KFXQfgpTW2jCgieylMCgRyIoxMF74LYZS6qvKKQ8G6NIopFIAlAhkxehLNR0t4DqClIAZLbUVKhFIQVAikBUj1XS0GBIBQGtdhO7hKZzTUBPiLyUCWTH6x2aoKYKmoymtdWEmo3FGpk47x5NIzikRyIrRPz5NQxFUFKesVoWxFAglAlkx+sdmaKoq/IrilNY6LxGowlh8pkQgK8LYdIzR6RiNRVQiCIeCrKoM0T2sUUjFX0oEsiKk5ikuhqaj6VprwyoRiO+UCGRFeL5/AiieFkMprXVhekenmY7F/Q5FSpgSgawIHf2pEkGRJYLaMAkHR3rG/Q5FSpgSgawIx/rGqQmXUVFWHE1HU1IVxhqSWvykRCArQkffOI1F1GIopbGqgrKAaTJ78ZUSgawIHf0TRXdbCCAYMFpqwxzsUolA/KNEIEVvdCpK39h0Qc9KdjqttWGVCMRXSgRS9E61GCqypqMpqZZD6mEsflEikKJXrC2GUjY2VwFw/3N9PkcipUqJQIreCwPJnrn1lcWZCFpqwzRVV/DTw71+hyIlKmeJwMxuNbMeM9uftq3BzO4xs8Pev/W5Or+UjhNDE9RFQoSLZNTR2QJmXHpWIw8810cioSGpJf9yWSK4Dbh81rbPAPc65zYD93rrIstyYnCStvqI32Esyxs2N9M3NsNB9ScQH+QsETjnfgIMzNp8BXC7t3w78Eu5Or+Ujs7BSdauKu5EcOnmJgDuP6x6Asm/fNcRtDjnurzlbqAlz+eXFcY5x4mhSdrqK/0OZVlaasOc01LDT5UIxAe+VRa75Px8894QNbNrzGyvme3t7VUlmsxtcCLKxEyctUV+awiSpYKHOwaYimoAOsmvfCeCF81sNYD3b898OzrnbnHO7XDO7Whubs5bgFJcTgwmWwwV+60hSCaCmViCh4/NvqMqklv5TgR3Ah/0lj8IfDfP55cVpnMw2Zms2CuLAS5ubyASCvL9p7oW3lkki3LZfPRfgAeBc8ys08w+AtwIvMXMDgNv9tZFluzEULJEsBISQWV5Ge9+1Rq+u+8kI1NRv8ORElKWqxd2zl09z1NvytU5pfR0Dk5SXVFGXSTkdyjLsmvPcQCaqiuYjMb5g/94itdvauJ9F6/3OTIpBepZLEUt1XTUzPwOJSvW1kdoq4+w59gAyfYUIrmnRCBFLdl0tPhvC6W7uL2R3tFpjvVp1jLJDyUCKWqdgxMroulougvb6oiEguxR6yHJEyUCKVojU1FGp2IrrkQQCga4sK2OQ90jROMJv8OREqBEIEXrpT4Exd2reC7tTVVE444DJzX2kOSeEoEUrc5UIlhhJQKADU3JOQoe0e0hyQMlAilaJ1ZQZ7LZasMhGqrKebhDiUByT4lAilbn4CThUIDGIp2reCEbGqvY2zGgOQok53LWoUwkV1Kdrx482k9NRYh/efgFnyPKjfamSh47PsiR3jE2t9T4HY6sYCoRSNEaGJ+hvqq4exSfzobGZD2Bbg9JrikRSFFKOEff2DTN1RV+h5IzDVXlNNdUqMJYck6JQIrS0ESUaNxxRk3Y71ByxszYuaGBRzoG/Q5FVjglAilKvaNTADTXrNwSAcBFG+o5MTR5apRVkVxQIpCi1Ds6DcAZKzwR7NjQAMDjx1UqkNxRIpCi1DM6TVV5kMqKld3wbXNLNaGgqYex5JQSgRSl3tHpFX9bCKCiLMjZLTXsPzHsdyiygikRSNFxztEzOk3zCq4oTrd1TS1PnxzR/ASSM0oEUnTGZ+JMRuMrvn4g5fy1dfSPz9A9MuV3KLJCrewbrLIi9ZRIiyFI9qI+6Q2ud/N9Rzh3dS2AprCUrFKJQIpOqbQYSmmti2DASTUhlRxRIpCi0zs6TSho1Bb5hPWZKi8L0FRToUQgOaNEIEUn1WIosEImrM/EmrowJ4dVRyC5oUQgRadndHpFDy0xlzWrIgxPRhmbjvkdiqxAvlQWm1kHMArEgZhzbocfcUjxGZ+OMTwZLYmK4nRrViUn3+kamtSQ1JJ1frYa+gXnXJ+P55cilOph21pbYiWCumQiODk8pUQgWadbQ1JUHniuD+OlsfpLRaQ8SENVOcf7x/0ORVYgvxKBA35kZo+a2TVz7WBm15jZXjPb29vbm+fwpFA9eKSfNasiRMqDfoeSd+e01nC4Z4ypaNzvUGSF8SsRXOqcew3wduATZvZzs3dwzt3inNvhnNvR3Nyc/wil4EzMxHj8hUE2NZdWaSDlwrV1xBKOg10agE6yy5dE4Jw74f3bA/wnsNOPOKS4PNIxSDTu2NRc7XcovljXUMmqSIgnOzUAnWRX3hOBmVWZWU1qGXgrsD/fcUjx+dmRPkJB48wSqx9ICZhxwdo6DveMMjQx43c4soL4USJoAe43syeAh4H/ds79wIc4pMg8eKSfV6+rp7ysdNs4XNi2ioSDHx7o9jsUWUHy/o1yzh11zr3Ke2x1zv15vmOQ4jM8EeWpE8O8/qxGv0Px1ZpVYRqqyvneE11+hyIrSOn+aSVF5aFj/TgHr9/U5HcovjIzLlxbx8+O9DE4rttDkh1KBFIUdj/TQ2V5kG3rVvkdiu+2rK4l4eD+59QfU7JDiUAK3uRMnO890cXbz19d0vUDKW31EeoiIX78rPrXSHboWyUF7wcHuhibjvErO9r8DqUgBMx4w+Ymfvxsr6avlKzQDGVSsHbtOQ7A1396lIaqco70jHG0V0MsAPzc2c3c9WQXB7tGOW9Nrd/hSJFTiUAK2sD4DEf7xnnN+nqshOYfWMjPn53sbf+Tw7o9JMunRCAF7bHjgxjwmvWqJE7XUhtmS2sNP35GiUCWT4lACtZ0LM7ejgHOOqOaVZXlfodTcH7+7Gb2Pj/AuCarkWVSIpCC9YP93YxOxbjsnDP8DqXg7NpznFjCEY07brz70Kn6FJGlUCKQgvSTZ3vZc2yAS85qor2pNMcWWsiZDZVEQkEeOz7odyhS5JQIpOAMTczwe995kuaaCt5yXovf4RSssmCAizbUc7BrRIPQybIoEUhBOTE0ya989UH6xqb5le1thIL6iJ7OxRsbcQ72HBvwOxQpYupHIAXj6ZMjfOifHmZyJs7tv7mTjr4Jv0MqePWV5Zy7upZHOgaYisYJh0pv5jZZPv25JQXhaz8+wq997UGmYwk+fEm7ksAivG5TIxMzce584qTfoUiRUiIQ38XBnnOiAAAL3ElEQVQTjn/d+wKT0TgffN0GWuvCfodUVDY2VdFSW8HN9z3HgEYklSVQIhDf/d29hznaO84V29YoCSyBmfGLr1rDyeEpfuPWPQxPRv0OSYqMEoH46psPPc/f/e9hXrN+FdvPbPA7nKK1samar71/O890j/Lhf3pYrYhkUZQIxBfOOb54z7N89r/286YtZ3DFtrV+h1T0fmHLGXzp6tfw1Ilh3v63P+Who/1+hyRFQolA8momluCuJ0/y3lse4m/vPcxV29v46vu3q5loFuzac5yB8RmuecMmZmIJrr7lIT5y2yMkEhqqWk5PzUclL14YmOBbjxzn9p89z9h0jPrKEO+8YDWvXreKf9vb6Xd4K8ra+gjXvvEs7tx3knsP9fCxOx7lr391G9UV+rrL3PTJkJzqHp7ir3/0DP/+WCcGnNNSw872Rja3VBPQsNI5U1EW5KrtbaxZFeHu/V1c8eX7+fRbz+FtW1sJBHTd5eWUCCSrovEEB06O8NSJYZ7qHOJ7T3QRTzh+85J2fvPSdg2bnEdmxiVnNXH1zvV89rv7+dgdj7H5jGretrWVC9rq2H5mPU3VFX6HKQXAl0RgZpcDfwsEga875270Iw5Zulg8wZHecQ52jXBiaJLe0WmO9I7x6PODTMzEAYiEgpzdUs1bzmuloapcScAnxwcm+Mil7Tx1YpgHnuvj5t3PkXAQMNjZ3sA7L1zD5Vtbaa5RUihVeU8EZhYE/h54C9AJPGJmdzrnns53LPKSqWicwy+Ocah7hMGJGaajCaZjCaZjcaZjCSZm4kzOxBkYn6FreJKTQ1PMxBOnjg+HAqyKlHNhWx3tTdW01UdYFQlpVrECETDjVW2reFXbKqLxBF1Dk5SHgvz3kyf57H/t53Pf3c/F7Y1cclYjm1tq2NRcRVN1BbXhEADTsQRT0TiT0TjReILqijJqIyFV8q8QfpQIdgLPOeeOApjZt4ArAF8TgXMO58ABiVPL3r/py2nP421POO94vH3Tl9Of9xpvJNLO5VzyeXDEEu7UD+7ETJyJmRjRuCMYSH6RU49gIFnsjyccfWPT9I1OMzETf9kPd/KHPLk8MD7DiyNTjEzFiISCVJYHqaooo7I8eOo1BsZnmN24xICyoFEWCFBeFiAUDFBZHqQuEuK1GytprQvTWhehsapcPwhFJBQMsL4xObT3b17Szouj0zzVOcyBk8M8OKvJacB4xeciXZX3eairLKcuUpZcnvWojYQoCwQIBiAYCFAWMMqCRiQUJFIepLK8jEgoSCho8/7hMN+fE/P9neEcxJ0jkXDEnSOecJgZoYBRFgxQFjRCgYD3+X7pvKnvY8I5DO97VwJ1Kn4kgrXAC2nrncDFuTjRn971NP/y8HGc9x/r/d6eWk7/wS5mBoTKkl+wUPClL1qZ96WrLA9ydksNkVCQaDz5l/5MPMFMLIEB7U1VXNi2ipbaMKtrw9REyrwv7sr/ApQ6M6O1NkzreWHecl4L07E4PSPT9I9PMz4dZ3wmRsCSn6tQMPlv0IypWLJ0MDWT/HdiJk7X8BRHe8fTSg7F88VaKOEFA0bQ7OUZadb+bvaGWSx1sKVv8/61V+6X2va1D2znDZubT/vay2Uuz7+CZnYVcLlz7re89Q8AFzvnrp213zXANd7qOUA/0JfPWAtYE7oW6XQ9XqJr8RJdCzjTObdgFvGjRHACWJe23uZtexnn3C3ALal1M9vrnNuR+/AKn67Fy+l6vETX4iW6Fpnz48buI8BmM2s3s3LgvcCdPsQhIiL4UCJwzsXM7FrghySbj97qnDuQ7zhERCTJl34EzrnvA99f5GG3LLxLydC1eDldj5foWrxE1yJDea8sFhGRwqLG3yIiJa4gEoGZXW5mz5jZc2b2mTme/10ze9rMnjSze83szLTn4ma2z3sUfaVzBtfio2b2lPd+7zez89Ke+33vuGfM7G35jTz7lnotzGyDmU2mfS6+mv/os2uha5G235Vm5sxsR9q2kvpcpO33smuxEj8XWZPs8erfg2SF8RFgI1AOPAGcN2ufXwAqveWPAf+a9tyY3+8hz9eiNm353cAPvOXzvP0rgHbvdYJ+vyefrsUGYL/f7yGf18Lbrwb4CfAQsKNUPxenuRYr6nORzUchlAhODTnhnJsBUkNOnOKcu885N+GtPkSy78FKlMm1GElbreKl/o1XAN9yzk07544Bz3mvV6yWcy1WmgWvhedPgb8AptK2ldznwjPXtZB5FEIimGvIidPNW/gR4O609bCZ7TWzh8zsl3IRYB5ldC3M7BNmdgT4S+C6xRxbRJZzLQDazexxM/uxmb0ht6Hm3ILXwsxeA6xzzv33Yo8tMsu5FrCyPhdZUwiJIGNm9n5gB/BXaZvPdMneg+8D/sbMNvkSXB455/7eObcJ+D3gD/2Ox0/zXIsuYL1z7tXA7wK7zKzWrxhzzcwCwBeAT/kdi98WuBYl9blYjEJIBBkNOWFmbwb+L/Bu59x0artz7oT371FgN/DqXAabYxldizTfAlKloMUeW+iWfC282yD93vKjJO8pn52jOPNhoWtRA5wP7DazDuC1wJ1eJWmpfS7mvRYr8HORPX5XUpDs1HaUZEVWqvJn66x9Xk3yP23zrO31QIW33AQcZo6Ko2J5ZHgtNqct/yKw11veyssrBY9S3JWCy7kWzan3TrJS8QTQ4Pd7yuW1mLX/bl6qIC25z8VprsWK+lxk8+H7VJVuniEnzOxPSH6x7yR5K6ga+LY3bvhx59y7gXOBr5lZgmTp5kZXxBPcZHgtrvVKR1FgEPigd+wBM/s3kvM6xIBPOOfivryRLFjOtQB+DvgTM4sCCeCjzrmB/L+L7MjwWsx3bCl+Luazoj4X2aSexSIiJa4Q6ghERMRHSgQiIiVOiUBEpMQpEYiIlDglAhGREqdEIAXLzFaZ2cez/JrXm1llNl8zW+cxsy3eqJiPm9kmM7vOzA6a2R35OL+ULjUflYJlZhuAu5xz58/xXJlzLraE1+wg2cGob9kBZvk83pDKZc65P/PWDwFvds515uP8UrpUIpBCdiOwyfsr+a/M7DIz+6k378TTkBx/yswe9vb5mpkFve1f8QYjPGBmf+xtuw5YA9xnZvd528a81z5gZv9jZjvNbLeZHTWzd3v7BL19HrHknBi/7W2/zNv3383skJndYUmvOE86M9vuDXr2qJn90MxWm9k7gOuBj5nZfd5Y+RuBu83s/5hZlZnd6r3Xx83sirTYbjKz/V5sn1zo/CKv4HfXZj30mO/BrPHjgcuAcaDdWz8X+B4Q8tZvBn7DW27w/g2SHGbgQm+9A2hKe00HvN1b/k/gR0AIeBWwz9t+DfCH3nIFsJfkEAeXAcMkx7sJAA8Cl851nrTzhYCfAc3e+q+R7B0LcAPw6bR9T70G8Hng/d7yKuBZkkNvfwz4d5IlifT3Pef59dBjrofvQ0yILNLDLjmuPsCbgO3AI97QIxGgx3vuV83sGpJj06wmOUHLk3O83gzwA2/5KWDaORc1s6dIJiKAtwIXmtlV3nodsNk79mHn3boxs33eMfefJv5zSA6Kdo8Xc5DkqJgLeSvwbjP7tLceBtYDbwa+6rzbZE5DJsgSKBFIsRlPWzbgdufc76fvYGbtwKeBi5xzg2Z2G8kfzrlEnXOpirIEMA3gnEuYWer7YcAnnXM/nHWey1L7e+Is/J0y4IBz7nUL7DfXcVc6556ZFcMiX0bklVRHIIVslOSwwvO5F7jKzM4AMLMGS85nXUsyYQybWQvw9kW85lx+SPLefcg7z9lmVrXE2J8Bms3sdd5rhcxsa4YxfNK8X34zSw23fg/w26mkZWYNC5xf5BWUCKRgueTY8Q94FaF/NcfzT5OcjOZHZvYkyR/F1c65J4DHgUPALuCBtMNuAX6wyErUr5OsnH7MzPYDX2Phv/znPI9LTq94FfAXZvYEsA94fQYx/CnJ+oUnzeyAt56K7bi3/QmSEzTNe36Ruaj5qIhIiVOJQESkxCkRiIiUOCUCEZESp0QgIlLilAhEREqcEoGISIlTIhARKXFKBCIiJe7/A0HDkvsRy3nMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f788cdd5990>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = sns.distplot(preds_rnd)\n",
    "ax.set_xlabel('treatment effect')\n",
    "ax.set_ylabel('density')\n",
    "plt.axvline(np.mean(preds_rnd), label='mean', lw=5)\n",
    "plt.axvline(np.exp(-1), color='r', label='truth')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now for the hard part. Lets try to use the same naive model for the data where there is confounding.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "sem_fn, df_obs, _ = regr_learner(df_obs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.42235738039016724"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_neg = df_obs.copy()\n",
    "X_pos = df_obs.copy()\n",
    "\n",
    "X_neg[\"medication\"] = 0\n",
    "X_pos[\"medication\"] = 1\n",
    "\n",
    "preds_rnd = sem_fn(X_pos)[\"prediction\"] / sem_fn(X_neg)[\"prediction\"]\n",
    "preds_rnd.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see now that the treatment effect is now underestimated ($log(0.41) = -0.89$), so we predict that the treatment has less impact than it actually does. But why is that? Well just to give one reason, if we recall the data generating process, we will notice that severity is causes a longer recovery period and a greater probability of treatment. Thus, some of the negative effect of severity goes to the correlation of treatment with recovery time, which appears to be smaller than the actual causal effect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xt8XHWd//HXZ3K/p23SW9KScG2hLQVKRUDtAiqKwgqIqLjgsrLiLiy77v5++HNVdH/rovBTdwVWWZZFkIqClwUUhUUqcitNL0AvXFqall6gae63SSaZ7++POSlpSJtJMnPOycz7+XjMg5kzJ+d8OqXzzvf7Pef7NeccIiKSvSJBFyAiIsFSEIiIZDkFgYhIllMQiIhkOQWBiEiWUxCIiGQ5BYGISJZTEIiIZDkFgYhIlssNuoBkVFVVubq6uqDLSK9XXkn897jjUnK4l3a3v2Pb4pqKlBxbRKaGtWvX7nfOVY+135QIgrq6OhoaGoIuI71WrEj8d9WqlByu7vpfv2Nbw43npeTYIjI1mNmOZPZT15CISJZTEIiIZDkFgYhIlpsSYwQiIocTi8XYtWsX0Wg06FICUVhYSG1tLXl5eRP6eQWBiEx5u3btoqysjLq6Osws6HJ85ZyjubmZXbt2UV9fP6FjqGtIRKa8aDTKjBkzsi4EAMyMGTNmTKo1pCAQkYyQjSEwZLJ/dgWBiEiWUxCIiGQ5DRaHxFsdfQA8vnrnQds/9a75QZQjIllEQSAiGWO0qVVSrfEQU7U0NjZy7rnnctppp/HMM89w6qmn8tnPfpavfe1r7Nu3j3vvvZcTTjiBa665ho0bNxKLxbjhhhu44IILaGxs5DOf+Qzd3d0A3HLLLZx++umsWrWKG264gaqqKjZu3Mgpp5zCj3/845SPh6QtCMzsTuAjwD7n3CJv23Tgp0Ad0Ahc4pxrTVcNU93O5m627O1g4ZzyoEsRkSRs3bqV+++/nzvvvJNTTz2VlStX8tRTT/Hggw/yzW9+k+OPP56zzjqLO++8k7a2NpYvX84555zDzJkzeeyxxygsLOS1117jk5/85IH51davX8+mTZuYO3cuZ5xxBk8//TRnnnlmSutO5xjBXcC5I7ZdDzzunDsGeNx7LaPYsreD2//4Ov/wwAtBlyIiSaqvr2fx4sVEIhFOOOEEzj77bMyMxYsX09jYyKOPPsqNN97I0qVLWbFiBdFolJ07dxKLxfjc5z7H4sWL+fjHP87mzZsPHHP58uXU1tYSiURYunQpjY2NKa87bS0C59yTZlY3YvMFwArv+Y+AVcD/TlcNU9W2pi5+8vxODGPzng46ozHKCid2x6CI+KegoODA80gkcuB1JBJhYGCAnJwcfv7zn3PciOnmb7jhBmbNmsULL7xAPB6nsLBw1GPm5OQwMDCQ8rr9vmpolnNur/f8TWCWz+cPva6+Ae55bgfTS/K55NR5xB2s3aHeM5FM8MEPfpDvf//7OOeARLcPQHt7O3PmzCESiXDPPfcwODjoa12BDRY755yZuUO9b2ZXAVcBzJ+fPVfO7G7tpX8gzgVLa6ipLCI3Yjy/vYUVx80MujSR0DvUQG5YfOUrX+G6665jyZIlxONx6uvrefjhh/nCF77ARRddxN133825555LSUmJr3XZUDKl5eCJrqGHhw0WvwKscM7tNbM5wCrn3JhLci1btsxl+sI0b538bgC+8cVb+fVLe/k/H15IaUEuP2t4g9yI8cDVp4/reKNdPRH2fyQiE7VlyxYWLlwYdBmBGu0zMLO1zrllY/2s311DDwKXe88vB/7b5/OH3v6uPgrzIpTk5wDwrvrpvLirnWjM36aiiGSPtAWBmf0EeBY4zsx2mdmVwI3A+83sNeAc77UMs7+rj+rSggPXCZ9aN53+wTgb3mgLuDIRyVTpvGrok4d46+x0nTMT7O/q58iqt/sHT62bjhk8v72F046cEWBlIpKpNNdQiDjnaO+NUVX29uViFcV5HDerjDWNLQFWJiKZTEEQIrHBOABVpQUHbV9eP521O1oZ8N4XEUklBUGIxAYTV3BVleYftP3k+dPo6R9kW1N3EGWJSIZTEITIUItgRsnBLYIjqxNjBtv3KwhEwqqtrY3bbrtt3D931113sWfPngOv6+rq2L9/fypLG5OCIERig3EqivLIzz34r6WuSkEgEnaHCoKxpoQYGQRB0DTUIRIbdFSPGB8AKC/Mo6q0gEYFgUhoXX/99Wzbto2lS5eSl5dHYWEh06ZN4+WXX+bRRx/lIx/5CBs3bgTg5ptvpquri0WLFtHQ0MCnP/1pioqKePbZZwH4/ve/z0MPPUQsFuP+++9nwYIFaa1dQRAajthgnBkjxgdWegvVlBbksHp784HXWrBG5BCuuw42bEjtMZcuhe9977C73HjjjWzcuJENGzawatUqzjvvPDZu3Eh9ff0hZwy9+OKLueWWW7j55ptZtuztG4CrqqpYt24dt912GzfffDN33HFHKv8076CuoZAYdBB37h1XDA2ZUVrA/q5+n6sSkYlavnw59fX1E/rZCy+8EIBTTjklLdNOj6QWQUjEBhIDxdVlowdBVWkBa/taicYGKczL8bM0kalljN/c/TJ84rjc3Fzi8bcv/45Go4f92aGpp9M17fRIahGExKHuIRgydElps1oFIqFUVlZGZ2fnqO/NmjWLffv20dzcTF9fHw8//HBSP+cXtQhCIjbowIzK4tEXoBkKiP1dfdRMK/KzNBFJwowZMzjjjDNYtGgRRUVFzJr19nIreXl5fPWrX2X58uXU1NQcNPh7xRVX8PnPf/6gwWK/KQhCYjAeJzdiRA6xKPX0knyMRBCISDitXLnykO9de+21XHvtte/YftFFF3HRRRcdeD18TGDZsmWsWrUqlSWOSl1DITHoOGQIAOTlRKgozqO5W11DIpJaCoKQiMcdOWP8bVSVFqhFICIppyAIiUHnyDlMiwASA8b7u/pI56pyIlNVNv+7mOyfXUEQEoNxRyQyVhAUEI3F6e7XamUiwxUWFtLc3JyVYeCco7m5mcLCwgkfQ4PFIdA/EMc5R84YQTA0GV2zuodEDlJbW8uuXbtoamoKupRAFBYWUltbO+GfVxCEQFtPYgD4cIPF8Pa9BBonEDlYXl7ehO/iFXUNhUKLFwRjtQgqi/PJMaOpU1cOiUjqKAhCoMW7W3isweKciDGjNJ8mtQhEJIUUBCEw1CIYa7AYEnMRNXUefp4SEZHxUBCEQGv3UItg7H2rywpo6e6nf0DrF4tIaigIQqClOwYk1yKYWVZA3MGOZi1SIyKpoSAIgdaefiJmGMl0DSWuFd66ryvdZYlIllAQhEBLd39SrQHgwFKWCgIRSRUFQQi09vSPeenokPzcCJVFeWxtUhCISGooCEKguat/zEtHh6suK2CbgkBEUkRBEAKtPcl3DYEXBPu6icezb14VEUk9BUHAnHO0dPcndenokOqyAnpjg+xp701fYSKSNRQEAeuNDdI3EE96jABgpq4cEpEUUhAErKU7uQnnhqsuS1w5tK1J9xKIyOQFEgRm9rdmtsnMNprZT8xs4hNpT3Gt3s1k42kRlOTnUFmcpxaBiKSE70FgZjXAtcAy59wiIAe41O86wmI88wwNMTOOri5lm4JARFIgqK6hXKDIzHKBYmBPQHUErqU7MZPoeC4fBTh6Zimv7uvMyhWZRCS1fA8C59xu4GZgJ7AXaHfOPep3HWHRMoGuIYAltZW09cTY0dyTjrJEJIsE0TU0DbgAqAfmAiVmdtko+11lZg1m1pDJy8+1difuKh5nDnDS/EoA1r/RmoaqRCSbBNE1dA6w3TnX5JyLAb8ATh+5k3PudufcMufcsurqat+L9EtLTz/TivMgiQnnhjt2VhnF+Tms39mWnsJEJGsEEQQ7gdPMrNjMDDgb2BJAHaHQ2t3PtOL8cf9cTsQ4sbZSQSAikxbEGMFq4AFgHfCSV8PtftcRFi3d/UwrGX8QQKJ7aMveDnr7B1NclYhkk0CuGnLOfc05t8A5t8g59xnnXNYuwtva08/0CbQIAE6aP42BuGPjnvYUVyUi2UR3FgespTs2oRbBytU7ed2bgfTOp7azcvXOVJcmIllCQRAg5xxtBwaLx6+sMI9pxXnsbNElpCIycQqCAPXGBhmIO8qLJhYEAPOnF/OGgkBEJkFBEKCO3gEAygsnHgTzphfTER2gvTeWqrJEJMsoCALUEU18eZcX5U74GPOnFwNo3iERmTAFQYA6vSAom0SLoKayiOkl+brDWEQmTEEQoLe7hibeIjAzTp5fyetN3exq1ViBiIyfgiBAb3cNTbxFAIn7CRzwy3W7U1CViGQbBUGAOqKJFkHZJFoEANOK8zmyuoQH1u3StNQiMm4KggB1eFf6TOaqoSEnz5/GjuYeGnZorEBExkdBEKCOaIz83AiFeTmTPtaiuRWU5OfwszVvpKAyEckmCoIAdUYHJjVQPFx+boTzl9bw0It7aO3uT8kxRSQ7KAgC1NEbS0m30JDLTz+CaCzOTxvUKhCR5CkIAtQRHaBsklcMDbdgdjmnHTmde57dkbJjikjmUxAEqDMaS1nX0JArTq9jd1tvSo8pIplNQRCgVHcNAZyzcBZzKwpTekwRyWyp/XVUxqUjOjCpeYZGGlqTYFFNBXvaoyk7rohkNrUIAtQZjU1qnqFDOXZWWcqPKSKZS0EQkL6BQaKxeMrHCABmq2tIRMZBQRCQTm96icnOMzSaiFnKjykimUtBEJDOFM0zJCIyWQqCgKRyniERkclQEAQkVVNQi4hMloIgIOoaEpGwUBAERF1DIhIWCoKAqGtIRMJCQRCQzugAEYOS/MmvRSAiMhnqoPbZ0DQQaxpbKcjN4SfPJ6aMPjvIokQkq6lFEJBobJDCPH38IhI8fRMFJBEE6hYSkeApCAKiIBCRsEgqCMzso2am0EihaCxOkYJAREIg2S/3TwCvmdm3zWzBZE9qZpVm9oCZvWxmW8zs3ZM95lSjMQIRCYukvomcc5cBJwHbgLvM7Fkzu8rMJjrx/b8Cv3XOLQBOBLZM8DhTVq+6hkQkJJL+ldQ51wE8ANwHzAE+Bqwzs2vGc0IzqwDeC/ynd9x+51zbeI4x1cWdo28griAQkVBIdozgAjP7JbAKyAOWO+c+ROK3+S+O85z1QBPwX2a23szuMLOSUc55lZk1mFlDU1PTOE8Rbn2xOICCQERCIdkWwYXAd51zi51zNznn9gE453qAK8d5zlzgZODfnXMnAd3A9SN3cs7d7pxb5pxbVl1dPc5ThFs0NghAkc9jBL39g76eT0SmhmS/id50zj05fIOZfQvAOff4OM+5C9jlnFvtvX6ARDBkjV4vCPxuEWzf3+3r+URkakg2CN4/yrYPTeSEzrk3gTfM7Dhv09nA5okca6qKDgQTBNuaunw9n4hMDYeda8jMrga+ABxlZi8Oe6sMeHoS570GuNfM8oHXgc9O4lhTTrQ/mDECBYGIjGasSedWAo8A/8LB/fidzrmWiZ7UObcBWDbRn5/qDrQIcv0dI9jWpK4hEXmnsYLAOecazeyvRr5hZtMnEwbZLBrQGMG2fWoRiMg7JdMi+AiwFnCADXvPAUemqa6M1jeQ6Boq8LlF8Pr+LuJxRyRiY+8sIlnjsEHgnPuI9996f8rJDn2xQXIjRm6Ov0EQjcXZ3dbLvOnFvp5XRMIt2RvKzhi66cvMLjOz75jZ/PSWlrmiA3EKArqZTAPGIjJSsr+S/jvQY2ZDdxJvA+5JW1UZri826PtA8RANGIvISMl+Gw045xxwAXCLc+5WEpeQygREY3EKAph5tLI4Ty0CEXmHZNcs7jSzLwGXAe/11ibIS19Zma1vYJCCXP+7ho6qLtWVQyLyDuNZj6APuNK7M7gWuCltVWW4IGce3bSng5Wrd7Jy9c5Azi8i4ZPsegRvOue+45z7o/d6p3Pu7vSWlrmiAY0RVJcW0NU3oMnnROQgyV41dKGZvWZm7WbWYWadZtaR7uIyVVBjBNVlBQA0dfX5fm4RCa9kv42+DZzvnKtwzpU758qcc+XpLCxTOecCGyM4EASdCgIReVuyQfCWcy7rlpNMh4G4I+78n2cIYFpxPjlmCgIROUiyVw01mNlPgV+RGDQGwDn3i7RUlcGG5hkK4oaynIgxozRfXUMicpBkg6Ac6AE+MGybAxQE4zQ0z1BhAGMEkOgeeqtDQSAib0sqCJxzWbVeQDodaBEEMEYAiSDYsreDgXg8kPOLSPgke9XQsWb2uJlt9F4vMbN/TG9pmenAzKNBtQhKC4g7aOnqD+T8IhI+yX4b/QfwJSAG4Jx7Ebg0XUVlsr6htQgCbBEA7NOAsYh4kg2CYufc8yO2DaS6mGwQHQhmmcohM8sKMeDNjmgg5xeR8Ek2CPab2VEkBogxs4uBvWmrKoP1HRgjCKZrKD83wvSSfN5SEIiIJ9mrhv4KuB1YYGa7ge3Ap9NWVQaLBrQ62XCzKwp5s11BICIJhw0CM/u7YS9/AzxBohXRDVwEfCd9pWWmoFYnG252RSGb93TQ0z9AcX6yvwuISKYa69uozHssA64GpgGVwOeBk9NbWmYKcnWyIbPLC3HAq29pSmoRGXvN4q8DmNmTwMnOuU7v9Q3Ar9NeXQYKcnWyIbPLCwF45c0Ols6rDLQWEQlest9Is4DhF573e9tknIKaeXS4aSX55OdE2LK3M9A6RCQcku0gvht43sx+6b3+U+CutFSU4YKaeXS4iBmzygt45U0FgYgkP8XEP5vZI8B7vE2fdc6tT19ZmatvIE5lUfCrfM4qL+TlNztwzmFmQZcjIgFK+pIR59w6YF0aa8kK0dgghV4ffZBmVxTSsKOVps4+ZoagHhEJTrCd1VmobyD4MQJIBAHAFnUPiWS94L+Rsohzjmgs+DECOPjKIRHJbgoCH/UNxANbnWyk4vxcZpcX6sohEVEQ+KkzmpinL+gbyoYcP7ecTXvagy5DRAIWWBCYWY6ZrTezh4OqwW9dfYkgCGp1spEW1VSwdV8XPf2aSFYkmwX5jfQ3wJYAz++7zmgMCG51spEW11QQd7Blr8YJRLJZIEFgZrXAecAdQZw/KF0HuobC0iIoB2DjbgWBSDYL6hvpe8D/Ag65cK6ZXWVmDWbW0NTU5F9ladQ51DUUkhbB7PJCqkrzeWm3xglEspnvQWBmHwH2OefWHm4/59ztzrllzrll1dXVPlWXXkMtgqBWJxvJzFhUU8FGBYFIVguiRXAGcL6ZNQL3AWeZ2Y8DqMN3Q4PFQS5KM9Limgpe29dF1Fs5TUSyj+/fSM65Lznnap1zdcClwO+dc5f5XUcQ3h4sDk8QLKqpYDDu2KwBY5GsFZ5vpCzQ2TcQ+OpkIy2uqQBQ95BIFgt0nULn3CpgVZA1+KkrOhCam8mGzKkoZEZJPi/tUhCIZKvw/GqaBbr6BkIxvcRwQwPGunJIJHuF61spw3VGB0JzD8FwGjAWyW6Bdg1lm67oQGjuKgZYuXonkBjEHow7vv3bVzh6Zimfetf8gCsTET+F79fTDNbeG6MoZGMEAPVVpeSY8do+zUQqko0UBD5q6+2nOD98QZCfG2H+jGK27usKuhQRCYCCwEdhbREAHDOzlL3t0QM3vYlI9lAQ+CQaGyQai1MUwhYBwNEzSwHYplaBSNZREPikvTdxV3FYg2BuZRFFeTm8piAQyToKAp8cCIKQdg1FzDiquoSt+zpxzgVdjoj4SEHgk7aeRBAU54f3it2jZ5bRER1gW5NaBSLZREHgk7C3CCAxYAzw+JZ9AVciIn5SEPikracfCO8YAcC0knzqZhRz1zON9A3oLmORbKEg8MlUaBEA/MlxM9nbHuUX63YHXYqI+ERB4JP23hgRC896xYdy9MxSTqyt4LZVWxkYPORKoiKSQcL9rZRB2npilBflETELupTDMjP++qxjeKOllwdf2BN0OSLiAwWBT9p7Y1QW5QVdRlLOXjCTBbPL+O7/vEprd3/Q5YhImikIfNLWG6NiigTBfWve4H3HVrOnLcrHbnuGu59pPDBTqYhkHgWBT9p7Y1QU5wddRtKOmFHCxSfX0tjczS/W79ZNZiIZLLx3N2WY9p5+jpheHHQZ43LivEpaevp5bPNb5OVEuHT5fHIi4R7jEJHxU4vAJ1Opa2i4FcdW875jq1nT2MLf/WwDMV1JJJJx1CLwQTzu6OiNUVk89YLAzPjgCbMpyI3w3xv20NM/yC2fOilUK62JyOSoReCDzr4B4o4p2SIYsuK4mXz9/BN4bPNb/MWPGujp17oFIplCQeCDDu+u4qkcBACXn17HTRcv4emt+7n8zue12L1IhlAQ+GBo5tGpHgQrV+8kNui4ZNk8GhpbufT253Q1kUgGUBD4oK03cVNW5RS6fPRwltRWcvbCmWx4o407n24MuhwRmSQFgQ+GJpybioPFh7LiuJkcP6ecb/5mC2saW4IuR0QmQUHgg0zpGhouYsbHT6llZlkBN//ulaDLEZFJUBD4oD1DBotHKsjL4coz61m9vYX1O1uDLkdEJkhB4IP23hgFuREKQ74WwURcunw+5YW53P7k60GXIiITpCDwQVtPf0aNDwxXWpDLZ959BL/d9Cbb93cHXY6ITICCwAftU3R6iWRdfnodeTkRtQpEpijfg8DM5pnZE2a22cw2mdnf+F2D39p6YlQWZcalo6OZWVbIhSfV8Mv1u+iIxoIuR0TGKYi5hgaALzrn1plZGbDWzB5zzm0OoBZftPfGmDfFZh5N1tA6BdNL8onG4nzlVxt5V/0MPvWu+QFXJiLJ8r1F4Jzb65xb5z3vBLYANX7X4adM7xoCqKksYk5FIQ2NunpIZKoJdIzAzOqAk4DVQdaRblNpmcqJMjOWHTGN3W297GnrDbocERmHwILAzEqBnwPXOec6Rnn/KjNrMLOGpqYm/wtMkb6BQXr6BzO+RQCwdN40ciNGww7daSwylQQSBGaWRyIE7nXO/WK0fZxztzvnljnnllVXV/tbYApl4vQSh1KUn8Oimgo2vNGmmUlFppAgrhoy4D+BLc657/h9fr8dmF4iQyacG8upddOJxuL88A+6lFRkqgiiRXAG8BngLDPb4D0+HEAdvhjqL59TURhwJf6oryrhxNoKvv/719i0pz3ockQkCUFcNfSUc86cc0ucc0u9x2/8rsMve9ujQPYEAcBHT5zLtJJ8vvizF+gf0BrHImGnO4vTbG9bLxGDWeXZEwTF+bnceOFiXn6zk288vIl4XIvXiISZgiDNdrdFmVlWSF5Odn3UZy+cxefeU8+Pn9vJtfet1+CxSIhl17dTAPa29zKnMntaA0NWrt5J3YwSzj1hNg+/uJdzv/dH7tJqZiKhpCBIs73tUeZWFgVdRiDMjPceW80nls1jR3M3963ZycCgxgxEwkZBkEbOOfa09TI3iwaKR3PivEo+euJcXn6zk68+uEkL3ouETBCTzmWFlat30tU3QN9AnD1t0QOTs2Wr046cQXtvjJWrd3JUdSlXnlkfdEki4lEQpFE23VWcjPcfP4u8nAg3PrKFd9VPZ1FNRdAliQjqGkqr9p5+IPPWKp6oiBnfvngJ00vyufa+9fT0DwRdkoigIEirtgxdtH4yppfk891LlrJ9fzdffzBjl6AQmVLUNZRG7T0xciNGSYE+5iFDYyXvObqanza8QUFuhAVzyrWQjUiA1CJIo7beGOVFeUTMgi4ldM5ZOJPZ5YX8cv1uevrURSQSJAVBGmXDymQTlZsT4eJTaunuH+DBF/cEXY5IVlMQpFE2rEw2GXMrizhrwSxe3NXOrU9sDbockaylzus0GYw7OnpjVOjS0cN637HV7O/q46bfvUJv/yBf/MCxmLrSRHylIEiTzmgMh64YGktOxLj4lFoWzC7jlie2snlvB584dR5/ctxM8nPVYBXxg4IgTQ7cTFaUHSuTTUbEjG9+bDE1lUX86NlGfv/yPorzcziyuoQjq0q54ow6Tp4/LegyRTKWgiBNDtxDoK6hpEQixjVnH8PVK47ij6/t5w+vNrF9fzdPb93Pr1/ay1//ydFcc9bR5GbZdN4iflAQpEl7z1CLQEGQjJFzMR07q4xjZ5XxvmOreeiFPfzr46/x4At7+LN3H8Fnz9A8RSKppF+v0mR/Vx9FeTkU5uUEXcqUVpiXw8eXzePjp9Syo7mbHz3TSLfuOxBJKQVBmuxo7mH+9OKgy8gYJ82fxiXL5rGjuYfP3rWGLoWBSMooCNKguauPpq4+6mYoCFJpSW0llyybx9odrXz8B8+yt7036JJEMoKCIA3WNLYCUFdVEnAlmefEeZXcecWp7Gzu5mO3PkNDY0vQJYlMeQqCNFjT2EJuxKjJ0iUq0213ay9/fmY9vbFBLv7Bs5zz//7A1x/apGUwRSZIVw2lwZrGFmqnFetSxzSaU1HEdeccw/PbW/jja/v5r6cbeWDtLt57TDUfPXEuZy+cSZ4+f5GkKAhSrLtvgE17OnjPMVVBl5LxCnJzeM8x1Zx25AxeebOTwbjj8Zf38euX9lJVWsC5i2axpKaSE2rKqZtRounARQ5B/zJSbN3OVgbjjroZGh/wS15O5MCyl4tqKnjtrU7W7Gjl/oZd/Pi5t+9PmF6Sz9EzSzlhbjmzygvZ0dzN603dNHX10dzVz1HVJfzLhUs4bnZZUH8UkUAoCFJszfYWIoYuHQ1ITsRYMKecBXPKiTtHS1c/e9p7ae3up6Wnn71tvazf2Ups0FGcn0NVaQHlRXlUzS5g0+52Pvxvf+QfPngcV73nSCIRTX4n2UFBkGLPN7Zw/Nxy3UgWAhEzqsoKqCorOGh73Dn6B+Lv+Ds6Z+Es/nvDbm585GX2tvVyw/knaCZUyQoKghTa2dzDuh1tXH76EUGXIocRMRs1qEsLcvnU8vls39/NHU9tpzAvh+s/tEBhIBlPQZBC33h4M7k5xpVnHsnvX94XdDkyAWbGl89bSN9AnB8++Tq9sUG+fN5CCnLVwpPMpevrUmTVK/v4ny1vcc1ZxzC7ojDocmQSfvL8Gxw3u4wzjprB3c/uYMVNq9i6rzPoskTSJpAgMLNzzewVM9tqZtcHUUMq9Q/E+cZDm6mvKuHPz6wLuhxJgYgZ5y359VeZAAAJeklEQVSZy2dOO4L23hgf+O6TXH7n8/xq/W5efrOD5q4+4nEXdJkiKeF715CZ5QC3Au8HdgFrzOxB59xmv2tJhY272/nLe9ayu62Xy99dx8/X7g66JEmhhXPKufbsY+iMxvjV+j1c99MNB97LiRhVpfnMLi9kUU0FJ86r5KR5lRxVXaorjmRKCWKMYDmw1Tn3OoCZ3QdcAAQeBM45nAM39By814nt/YNxWrr62dfZx6Y97TTsaOWRl/ZSnJ/LpafO0/XnGaq8MI/ywjyuXnEUu1t7aeuN0RmN0RUdoKtvgNaefn6+bhf3emsqFORGOLG2ktppRcytTDxqphUxoySfovwcivJyKM5PTFGeMyIw3CEaGWZgJFoqZmgAW1IqiCCoAd4Y9noX8K50nOjqH6/lD682Hfgyjztg2Bf7yC/88ZpdXshlpx3BEdNLKMrXYGKmi5gxb3ox80Z5L+4c+7v62NXSyxutPbzVEeXVtzrpiMZIVw/SUDgMDwUb8f7b2w96Ma79h7YreoLx0DVncmR1aVrPEdqrhszsKuAq72WXmb0SYDlVwP6RG3cAq1N9ptPGdenpqHUdyqe/Ne5qJmNctflIdY1fWGsLa12QwtqO+qdJ/XhSXyhBBMFuOOiXqlpv20Gcc7cDt/tV1OGYWYNzblnQdYwU1rogvLWprvELa21hrQvCXdtogrhqaA1wjJnVm1k+cCnwYAB1iIgIAbQInHMDZvbXwO+AHOBO59wmv+sQEZGEQMYInHO/AX4TxLknKBRdVKMIa10Q3tpU1/iFtbaw1gXhru0dzE3kchkREckYmmJCRCTLKQg8Y017YWafN7OXzGyDmT1lZseHpbZh+11kZs7MfLlaIYnP7Aoza/I+sw1m9hd+1JVMbd4+l5jZZjPbZGYrw1CXmX132Of1qpm1+VFXkrXNN7MnzGy9mb1oZh8OSV1HmNnjXk2rzKzWp7ruNLN9ZrbxEO+bmf2bV/eLZnayH3VNSOJu2ux+kBi03gYcCeQDLwDHj9infNjz84HfhqU2b78y4EngOWBZGOoCrgBuCenf5zHAemCa93pmGOoasf81JC6mCMtndjtwtff8eKAxJHXdD1zuPT8LuMenz+y9wMnAxkO8/2HgERL34p0GrPajrok81CJIODDthXOuHxia9uIA51zHsJclJG5MDkVtnn8CvgVEQ1ZXEJKp7XPArc65VgDnnB/zho/3M/sk8BMf6oLkanNAufe8AtgTkrqOB37vPX9ilPfTwjn3JNBymF0uAO52Cc8BlWY2x4/axktBkDDatBc1I3cys78ys23At4Frw1Kb1+Sc55z7tU81JVWX5yKvWfyAmY02O0M6JFPbscCxZva0mT1nZueGpC4g0d0B1PP2F1y6JVPbDcBlZraLxFV/14SkrheAC73nHwPKzGyGD7WNJem/76ApCMbBOXerc+4o4H8D/xh0PQBmFgG+A3wx6FpG8RBQ55xbAjwG/CjgeobLJdE9tILEb97/YWaVgVZ0sEuBB5xzg0EXMswngbucc7Ukuj3u8f7/C9rfA+8zs/XA+0jMVBCmzy30wvCXGAZJTXsxzH3An6a1oreNVVsZsAhYZWaNJPoiH/RhwHjMz8w51+yc6/Ne3gGckuaakq6NxG9nDzrnYs657cCrJIIh6LqGXIp/3UKQXG1XAj8DcM49CxSSmFMn0Lqcc3uccxc6504Cvuxt822Q/TDG+70SnKAHKcLwIPHb4eskmuJDA1InjNjnmGHPPwo0hKW2Efuvwp/B4mQ+sznDnn8MeC4snxlwLvAj73kViSb8jKDr8vZbADTi3ecTos/sEeAK7/lCEmMEaa0xybqqgIj3/J+Bb/j4udVx6MHi8zh4sPh5v+oa958j6ALC8iDR1H2VxBUKX/a2fQM433v+r8AmYAOJAalDfhn7XduIfX0JgiQ/s3/xPrMXvM9sQVg+M+8f53dIrIPxEnBpGOryXt8A3OjXZzWOz+x44Gnv73MD8IGQ1HUx8Jq3zx1AgU91/QTYC8RItDCvBD4PfH7Y/2O3enW/5Ne/y4k8dGexiEiW0xiBiEiWUxCIiGQ5BYGISJZTEIiIZDkFgYhIllMQSGiZWaWZfSHFx7zOzIpTecxUncfMFnizjq43s6PM7Foz22Jm9/pxfsleunxUQsvM6oCHnXOLRnkv1zk3MIFjNpK4nnv/pAtM8Xm8KZZznXP/13v9MnCOc26XH+eX7KUWgYTZjcBR3m/JN5nZCjP7o5k9SOJGMMzsMjN73tvnh2aW423/dzNr8NYa+Lq37VpgLvCEmT3hbevyjr3JzP7HzJZ7c9q/bmbne/vkePus8SbQ+0tv+wpv3wfM7GUzu9ebg/4d5xnOzE4xsz+Y2Voz+52ZzfHm9r8OuNqb8/8HJKZefsTM/tbMSrz575/3WgwXDKvtZjPb6NV2zVjnF3mHoO9o00OPQz0Ycfs+iQniuoF67/VCEhPb5XmvbwP+zHs+3ftvDom7rZd4rxuBqmHHdMCHvOe/BB4F8oATgQ3e9quAf/SeFwANJKY8WAG0k5hDJgI8C5w52nmGnS8PeAao9l5/Am/NARJ3FP/9sH0PHAP4JnCZ97ySxF20JcDVwAMkWhLD/9yjnl8PPUZ7BLJ4vcgkPO8Sk8QBnE1iIrs1ZgZQBAytK3CJmV1FYq6aOSSmR3hxlOP1A7/1nr8E9DnnYmb2EokgAvgAsMTMLvZeV5CYoK7fq2cXgJlt8H7mqcPUfxyJSQIf82rOITFNwVg+AJxvZn/vvS4E5gPnAD9wXjeZc+5w8+OLjEpBIFNN97DnRmLiuC8N38HM6klMTXyqc67VzO4i8cU5mphzbmigLA70ATjn4mY29O/DgGucc78bcZ4VQ/t7Bhn735QBm5xz7x5jv9F+7iLn3CsjahjnYUTeSWMEEmadJKbZPpTHgYvNbCaAmU33FnQpJxEY7WY2C/jQOI45mt+R6LvP885zrJmVTLD2V4BqM3u3d6w8MzshyRquMe+b38xO8rY/BvzlUGiZ2fQxzi/yDgoCCS3nXDPwtDcQetMo728msUDQo2b2IokvxTnOuRdIrEf8MrCSxIyZQ24HfjvOQdQ7SAxOr7PEQuU/ZOzf/Ec9j0sst3gx8C0zG5rF8/QkavgnEuMLL5rZJu/1UG07ve0vAJ863PlFRqPLR0VEspxaBCIiWU5BICKS5RQEIiJZTkEgIpLlFAQiIllOQSAikuUUBCIiWU5BICKS5f4/RAxKjttOrEMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f788cd42690>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = sns.distplot(preds_rnd)\n",
    "ax.set_xlabel('treatment effect')\n",
    "ax.set_ylabel('density')\n",
    "plt.axvline(np.mean(preds_rnd), label='mean', lw=5)\n",
    "plt.axvline(np.exp(-1), color='r', label='truth')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, lets see how does out model performs in both the factual and counter factual data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Factual Evaluation: {'r2_evaluator__recovery': 0.8669406871788605}\n",
      "Counter Factual Evaluation {'r2_evaluator__recovery': 0.8341172823328981}\n"
     ]
    }
   ],
   "source": [
    "eval_fn = r2_evaluator(prediction_column=\"prediction\",target_column=\"recovery\")\n",
    "\n",
    "print \"Factual Evaluation:\", eval_fn(sem_fn(df_obs))\n",
    "print \"Counter Factual Evaluation\", eval_fn(sem_fn(df_cf))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## IPTW  Model\n",
    "\n",
    "Having played a little with the a naive benchmark, lets now use a framework that is meant to uncover causal relation, the Inverse Probability of Treatment Waiting framework. The idea is very simple. We first fit a model to predict the probability of treatment. Then, we overweight two kinds of samples: the treated with low form whom the predicted probability of treatment is low and the untreated for whom the predicted probability is high. The intuition is that this will make the model focus on the treated that looks like the untreated and on the untreated that looks like the treated.\n",
    "\n",
    "In order for this to work we have to watch out for two things:\n",
    "\n",
    "    1) Treatment assignment should not be deterministic on the covariates. In other words, everyone should have a probability of being treated. This is some times referred as the positivity assumption of causal inference.\n",
    "    2) All variables that influence the responde (recovery time) and treatment assignment should be observable. Another way to say this is that treatment is randomly assigned given once we account for the covariates. This is sometimes referred as the unconfoundedness assumption.\n",
    "    \n",
    "Trouble is, although we can check for assumption 1) wit the data, there is no way to check for assumption 2). Here, we know how the data was generated and we know that we have unconfoundedness, but, in general, we have to rely on experts intuition about the problem we try to model.\n",
    "\n",
    "OK. Back to the model. Lets split our data fist."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "train, test, _, _ = train_test_split(df_obs, df_obs, test_size=0.33, random_state=42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we will fit the IPTW learner. It takes tow extra learners as inputs: a classification model to predict the probability of treatment and a response model that, on this case, will be a regression one, since the our response variable is continuous."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "causal_learner = IPTW_learner(features=features, target=\"recovery\", treatment = \"medication\",\n",
    "                              ps_learner = xgb_classification_learner(num_estimators=25, extra_params={\"max_depth\": 2}),\n",
    "                              response_learner = xgb_regression_learner)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The IPTW predict function returns a DataFrame with the predicted response variable, the propensity scores and the weight of the each samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sex</th>\n",
       "      <th>age</th>\n",
       "      <th>severity</th>\n",
       "      <th>medication</th>\n",
       "      <th>recovery</th>\n",
       "      <th>prediction</th>\n",
       "      <th>propensity_score</th>\n",
       "      <th>weight</th>\n",
       "      <th>counterfactual</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>23990</th>\n",
       "      <td>0</td>\n",
       "      <td>22.357267</td>\n",
       "      <td>0.797063</td>\n",
       "      <td>1.0</td>\n",
       "      <td>9</td>\n",
       "      <td>12.010621</td>\n",
       "      <td>0.903162</td>\n",
       "      <td>1.107221</td>\n",
       "      <td>12.305794</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8729</th>\n",
       "      <td>0</td>\n",
       "      <td>10.764319</td>\n",
       "      <td>0.165695</td>\n",
       "      <td>0.0</td>\n",
       "      <td>12</td>\n",
       "      <td>11.791059</td>\n",
       "      <td>0.046857</td>\n",
       "      <td>1.049161</td>\n",
       "      <td>12.667633</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3451</th>\n",
       "      <td>1</td>\n",
       "      <td>24.046842</td>\n",
       "      <td>0.494215</td>\n",
       "      <td>1.0</td>\n",
       "      <td>13</td>\n",
       "      <td>15.043891</td>\n",
       "      <td>0.791081</td>\n",
       "      <td>1.264093</td>\n",
       "      <td>15.231667</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2628</th>\n",
       "      <td>1</td>\n",
       "      <td>30.389036</td>\n",
       "      <td>0.567990</td>\n",
       "      <td>1.0</td>\n",
       "      <td>24</td>\n",
       "      <td>20.715506</td>\n",
       "      <td>0.889086</td>\n",
       "      <td>1.124750</td>\n",
       "      <td>20.982327</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>38352</th>\n",
       "      <td>1</td>\n",
       "      <td>23.535528</td>\n",
       "      <td>0.050463</td>\n",
       "      <td>0.0</td>\n",
       "      <td>28</td>\n",
       "      <td>26.280563</td>\n",
       "      <td>0.137160</td>\n",
       "      <td>1.158963</td>\n",
       "      <td>25.861427</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       sex        age  severity  medication  recovery  prediction  \\\n",
       "23990    0  22.357267  0.797063         1.0         9   12.010621   \n",
       "8729     0  10.764319  0.165695         0.0        12   11.791059   \n",
       "3451     1  24.046842  0.494215         1.0        13   15.043891   \n",
       "2628     1  30.389036  0.567990         1.0        24   20.715506   \n",
       "38352    1  23.535528  0.050463         0.0        28   26.280563   \n",
       "\n",
       "       propensity_score    weight  counterfactual  \n",
       "23990          0.903162  1.107221       12.305794  \n",
       "8729           0.046857  1.049161       12.667633  \n",
       "3451           0.791081  1.264093       15.231667  \n",
       "2628           0.889086  1.124750       20.982327  \n",
       "38352          0.137160  1.158963       25.861427  "
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sem_fn, train_pred, _ = causal_learner(train)\n",
    "test_pred = sem_fn(test)\n",
    "train_pred.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we process, it a good idea to validate our propensity score classifier. It better not overffit! If iid does, it will be very confident, giving extremely low probability of treatment to the untreated and extremely hight probability of treatment for the treated. This will result in no overlap in in the propensity score for the treated and untreated, which will be a near violation of the positivity assumption. Also, if our model can precisely predict the treatment assignment, it means it is a deterministic function of the covariates, also violating the positivity assumption."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "({'roc_auc_evaluator__medication': 0.9780837845595602},\n",
       " {'roc_auc_evaluator__medication': 0.976750491771322})"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ps_eval_fn = roc_auc_evaluator(prediction_column=\"propensity_score\", target_column=\"medication\")\n",
    "ps_eval_fn(train_pred), ps_eval_fn(test_pred) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "No handlers could be found for logger \"matplotlib.legend\"\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xd85HWd+PHXe3p6z262ZJMt7C69BJYiClIED0FPUFEsyP1Qz3Iqej+9837q6RX19O48PREVsYEFRYogIh2EhQV2YXsv2U02vU6mf35/fGZ2s8kkmSTTsvN+Ph55zMx3vjOfz2Sz8/5+2vsjxhiUUkoVLkeuK6CUUiq3NBAopVSB00CglFIFTgOBUkoVOA0ESilV4DQQKKVUgdNAoJRSBU4DgVJKFTgNBEopVeBcua5AKmpra01TU1Ouq6GUUnPKSy+91GWMqZvqvDkRCJqamli3bl2uq6GUUnOKiOxL5TztGlJKqQKngUAppQqcBgKllCpwc2KMQCmlCl04HKa1tZVAIDDuOZ/Px6JFi3C73TN6bw0ESik1B7S2tlJWVkZTUxMicuS4MYbu7m5aW1tpbm6e0Xtr15BSSs0BgUCAmpqaY4IAgIhQU1OTtKWQKg0ESik1R4wNAlMdT5UGAqWUKnAaCJRSqsDpYLFSSo1y59r9SY+/e01jlmsynjEmaTeQMWZW75uxFoGI3C4iHSKycczxj4vIVhHZJCJfz1T5Sil1PPH5fHR3d4/70k/MGvL5fDN+70y2CO4AvgP8NHFARC4GrgFOM8YERaQ+g+XPDet+nPx4y43ZrYdSKq8tWrSI1tZWOjs7xz2XWEcwUxkLBMaYp0SkaczhjwD/bowJxs/pyFT5SimVC5nqWnK73TNeJzCVbA8WnwBcKCJrReRJETk7y+UrpZQaI9uDxS6gGjgXOBv4tYgsNUlGOkTkZuBmgMbG3A/SKKXU8SrbLYJW4HfGegGIAbXJTjTG3GaMaTHGtNTVTbmvglJKqRnKdiD4PXAxgIicAHiArizXQSml1CgZ6xoSkbuAi4BaEWkFvgjcDtwen1IaAt6frFtIKaVU9mRy1tD1Ezx1Q6bKVEopNX2aYkIppQqcBgKllCpwGgiUUqrAaSBQSqkCp4FAKaUKnAYCpZQqcBoIlFKqwGkgUEqpAqeBQCmlCpwGAqWUKnAaCJRSqsBpIFBKqQKngUAppQqcBgKllCpwGgiUUqrAZXvP4rlj3Y/HH2u5Mfv1UEqpDMtYi0BEbheRjvhuZGOfu0VEjIgk3a9YKaVU9mSya+gO4IqxB0VkMXA5sD+DZSullEpRxgKBMeYpoCfJU/8J/D2gexUrpVQeyOpgsYhcAxw0xmzIZrlKKaUmlrXBYhEpBv4B2y2Uyvk3AzcDNDY2ZrBmSilV2LLZIlgGNAMbRGQvsAh4WUTmJzvZGHObMabFGNNSV1eXxWoqpVRhyVqLwBjzGlCfeBwPBi3GmK5s1UEppdR4mZw+ehfwHLBSRFpF5KZMlaWUUmrmMtYiMMZcP8XzTZkqWymlVOo0xYRSShU4DQRKKVXgNBAopVSB00CglFIFTgOBUkoVOA0ESilV4DQQKKVUgdNAoJRSBU4DgVJKFTgNBEopVeA0ECilVIHTQKCUUgVOA4FSShU4DQRKKVXgNBAopVSB00CglFIFLpM7lN0uIh0isnHUsW+IyFYReVVE7hGRykyVr5RSKjWZbBHcAVwx5tgjwMnGmFOB7cDnM1i+UkqpFGQsEBhjngJ6xhz7kzEmEn/4PLAoU+UrpZRKTS7HCD4IPJTD8pVSSpGjQCAi/whEgF9Mcs7NIrJORNZ1dnZmr3JKKVVgsh4IROQDwFXAe4wxZqLzjDG3GWNajDEtdXV1WaufUkoVGlc2CxORK4C/B95gjPFns2yllFLJZXL66F3Ac8BKEWkVkZuA7wBlwCMisl5Ebs1U+UoppVKTsRaBMeb6JId/lKnylFJKzYyuLFZKqQKngUAppQqcBgKllCpwGgiUUqrAaSBQSqkCp4FAKaUKnAYCpZQqcBoIlFKqwGkgUEqpAqeBQCmlCpwGAqWUKnAaCJRSqsBpIFBKqQKngUAppQqcBgKllCpwGgiUUqrAaSBQSqkCl8mtKm8XkQ4R2TjqWLWIPCIiO+K3VZkqXymlVGoy2SK4A7hizLHPAY8aY1YAj8YfK6WUyqGMBQJjzFNAz5jD1wA/id//CfDWTJWvlFIqNdkeI5hnjGmL328H5k10oojcLCLrRGRdZ2dndmqnlFIFKGeDxcYYA5hJnr/NGNNijGmpq6vLYs2UUqqwpBQIROQtIpKOoHFYRBri79kAdKThPZVSSs1Cql/u7wR2iMjXRWTVLMq7D3h//P77gXtn8V5KKaXSIKVAYIy5ATgD2AXcISLPxfvwyyZ6jYjcBTwHrBSRVhG5Cfh34DIR2QFcGn+slFIqh1ypnmiMGRCRu4Ei4JPA24DPisi3jTH/k+T86yd4q0tmVFOllFIZkeoYwTUicg/wBOAGzjHGXAmcBtySueoppZTKtFRbBH8N/Gd8bcARxhh/vMtHKaXUHJXqYHH72CAgIl8DMMY8mvZaKaWUyppUA8FlSY5dmc6KKKWUyo1Ju4ZE5CPA3wLLROTVUU+VAc9msmJKKaWyY6oxgjuBh4B/49gEcYPGmLF5hJRSSs1BUwUCY4zZKyIfHfuEiFRrMFBKqbkvlRbBVcBL2LxAMuo5AyzNUL2UUkplyaSBwBhzVfy2OTvVUUoplW2pLii7QERK4vdvEJFviUhjZqumlFIqG1KdPvo9wC8iiZXEu4CfZaxWSimVJ17e18tDG9umPnEOSzUQROL7B1wDfMcY813sFFKllDquvbi3h2d3dtE/Es51VTIm1UAwKCKfB24A/hDfm8CduWoppVTuxYyhbSBAzMBzu7pzXZ2Mmc5+BEHgJmNMO7AI+EbGaqWUUnmgdzhEKBID4Okdx++WuSklnYt/+X9r1OP9wE8zVSmllMoHh/oDAFQWuXl6R1eOa5M5qc4a+msR2SEi/SIyICKDIjKQ6coppVQutfWP4BA4b1kN+3v87OseznWVMiLVrqGvA1cbYyqMMeXGmDJjTPlMCxWRT4nIJhHZKCJ3iYhvpu+llFKZ0tYXoK7My+oG+3X31HHaKkg1EBw2xmxJR4EishD4BNBijDkZcALvSsd7K6VUOrX1j9BQUURNiYdFVUU8vf34HCdIdWOadSLyK+D32EFjAIwxv5tFuUUiEgaKgUMzfB+llMqIoWCEgUCEhgofIsKFK+p4YMMhwtEYbmeq19BzQ6qfphzwA5cDb4n/XDWTAo0xB4H/APYDbUC/MeZPM3kvpZTKlPb4QHFDRREAFyyvYTAYYVv7YC6rlRGpzhq6MV0FikgVdmFaM9AH/EZEbjDG/HzMeTcDNwM0Nmo2C6VUdrX1jwDQUGGHMBuri+PHA5y8sCJn9cqEVGcNnSAij4rIxvjjU0XkCzMs81JgjzGm0xgTBn4HnD/2JGPMbcaYFmNMS11d3QyLUkqpmWnrD1BR5KbEa6+X55fbgHB4IJDLamVEql1DPwA+D4QBjDGvMvMB3v3AuSJSLCICXAKkZSBaKaXS5VDfyJHWAEBNqReHFHYgKDbGvDDmWGQmBRpj1gJ3Ay8Dr8XrcNtM3ksppTIhZgxdQ0HmlR8NBE6HUFfmLehA0CUiy7Cb0SAi12IHemfEGPNFY8wqY8zJxpj3GmOCU7/qOBUcgvbXcl0LpdQo/lCUmIHyomNTqs0r93F4YHpfV/5QhHvXH2QwkL9J61INBB8Fvg+sEpGDwCeBD2esVoVk79Ow7kfQszvXNVFKxSW+tMu8x86nqS/zTbtF8GprP2v39PDEtvxdgzBpIBCRT4vIp4G3Ag8C/wLcih3gfXvmq1cABuMNqx06g1apfDEUsD3fpWMCwbxyLx2D02sR7OwYAmDTofzNyjNVi6As/tMCfASoAiqxrYEzM1u1AjHYBg4ndG6F3n25ro1SCruYDKDMd2wgmF/uo2c4RDASTel9ojHDrs5EIOhPbyXTaNJAYIz5sjHmy9i002caYz5jjLkFOAvQyf2zFQ7AcBcseR24i7VVoFSeGEy0CHxjWwR28LgjxXGCg71+gpEYJR4nmw4NYPf3yj+pjhHMA0KjHofix9RsdO8ADFQtgaUXQccmGNBsG0rl2lAwgtspeF3OY47Xl3sB6BhMbZxgR8cQApy7rIae4RDteTrjKNVA8FPgBRH5koh8CVgL3JGpShWMjq32trQBFrbY+33aPaRUrg0GwpT5xm/COO/IorLUWgQ7O4ZYUFnE8rpSADYezM9xgpQCgTHmX4Abgd74z43GmH/LZMUKQucWEAeU1oGvAhAY6c11rZQqeEPByLiBYhgdCKa+sg+Eoxzo9bOivpT5FT5E8necINXsoxhjXsYuAlPp0rEVSmrBEf9n8FVoIFAqDwwGItSVeccdryp243E6UmoR7O4cJmZgeX0pXpeT5tqSvJ05dHzlUp1rOrfabqGEoioY6ctdfZRSwMQtAhGhvjy11cW7u4ZwO+VIsrqTF1Sw6WB+tgg0ECQTHIRoaOrzZiMcgN49UDb/6LGiSm0RKJVj0ZjBH4qOmzGUYFcXTx0IuodC1JZ6ccX3LjhpQTmH+gP0Dmf4u2UGNBAk86M3wZ/+CTb8EgbbM1NG13YwsTGBoAoCvfa4Uionjqwh8I4fLAa7qCyVQNA/EqZyVIqKROrqfOwe0kAwVshvp3EWVcKhl+H570Em5v52xmcMjQ0EsajNP6SUyomJVhUn1Jf5UlpH0DcSoqLYc+TxSQvsvscb83DAWAPBWN077e0JV8Lqt0CwHwIZ+Ifr2GIHiUtG7bVQVGVvAzpOoFSuDAXjeYYm6RoaDEYYDk6cgDkQjhIIx45pEVQWe6gr87Knczi9FU4DDQRjdW23t6X1ULbA3h/MwCKvzm1Qs/zojCE4Ggh0nECpnJloVXHCvPiissm6h/pGbDCpLD62e6m21EvXUP4lW9ZAMFbXdju3v6QOyuIzegZmnHF7Yv37oXLJscd8GgiUyrXEGMFEXUOpLCrr99sB4cqisYHAQ5cOFs8BXdvtF7TTDZ5i8FVmpkUw1AFlY7J0uIvA6dVAoFQODQYi+NwO3M7kX49H8g1NkmYi0SIYPUYAtkXQrS0CS0QqReRuEdkqIltE5Lxc1COprh1Qe8LRx+UL0t8iiEVhuBNKxwQCkfhaAg0ESuWKXUOQfMYQpNg15A/jkPHjDLWlHrqGgnmXfC5XLYL/Bv5ojFkFnEa+7Fkci9rB4toVR4+VNcDQYftcuvi77RTRsYEAdC2BUjk2GIhMOFAMtsuo2OOkrX+yQBCiosiNQ+SY4zWlXgLhGP5QGr9P0iDrgUBEKoDXAz8CMMaEjDH5MU2m/wBEAuNbBCZqu3LSZeiwvS2tH/+ctgiUyqmhYHjC8QGwq4sbKny0TxYIRsJUFHnGHa8tta2JfBswzkWLoBnoBH4sIq+IyA9FpCQH9Riva4e9HR0IEgPG6RwnOBIIkrUIqiA0BOGR9JWnlErZUDAy4YyhhAWVRRzqm/j/aL8/TFXx+O6lmlIbHLqG8mvAOBeBwIXd3ex7xpgzgGHgc2NPEpGbRWSdiKzr7MzSXp+JqaOjA0FpPYgzvfsEJFoXE3UNge5LoFQOJOb/j92reKyGCh+HJmgRRGOGgUCYiiSBoLZEWwQJrUCrMWZt/PHdJNn20hhzmzGmxRjTUldXN/bpzOjaDkXVUFJz9JjDZYPBYBoHjCfrGkpMIe0/kL7ylFIp6YzvRzxZ1xDYFkHXUJBQZHw6mMFAmJiBymRdQ2X2WHehtwiMMe3AARFZGT90CbA52/VIqmsH1K0cf7x8QfpbBJ4y8CTpEUssKutvTV95SqmUJK7UJxssBlhQUYQxyWcO9U+wmAyguiTRNaQtAoCPA78QkVeB04F/zVE9jtW1/dgZQwllDTbtQ7pSRA8dTt4aALtuAdFAoFQOHGkRJNmdbLSGSruW4GCScYI+f3wNQdH49/C6nJT7XHm3liDljWnSyRizHmjJRdkT8vfYuf2jxwcSEgPGndugcc3syxrqSD4+AOB0gbdMA4FSOdA5lFrXUENFEQBt/ckCQfJVxQk2zUSBdw3lrcRewVXN458riV+99+xKT1mD7RO3CAB85UfHEZRSWZP6GIFtERzqG9811DcSpsjtxOt2jnsO8jPfkAaChMH4F+/otNAJxVU2/1DP7vSUNVmLAMBbnt7BaaVUSrqGghR7nDgdMul5xR4XlcXuCVoE4aTjAwk1pR668yzfkAaChMnm9jtcdhA3HYEgPGJTW0/WIvCWHw1MSqms6RwMTjlQnNBQUURbkhZB/0g46fhAgrYI8tlkgQBsNtLuNHQNTbaGIMFXbscrohPnO1dKpV/nYHDKbqGEBRW+cYPFxhh6/aEpWwR9/jDhaP7sRKiBIGGw3a4hcI2f+wtAcS307Jn9bmWpBAJvOWBsMFBKZU3XUIiyKWYMJTRU+sblG2ofCBCMxKgv8034ukSaiZ486h7SQJAwdHjyL+eSOtul4++efTkw9WAxwFCG9ktWSo1jjJlWi6Chooj+kTD+0NGW+7b2QeBoqupkakvzby2BBoKEwfbx+wOMVlJrb2c7TjBVFxTEWwToOIFSWTQcijISjqbeNZRk5tD2w/FAUOad8HU1RxLPaYsg/wx1QGmSGUMJib2FZztOMNQByLF7FY/lq4ifqy0CpbIlMXU01cHiBUnWEmxrH6LM56J4kmCS6BrKp0VlGgjA9vsPTdEiKK5OzxTSocO2deGc5I/NW2ZvtUWgVNYkumqmyjyasKDSBoLRWUi3Hx6ctFsIRmcg1UCQX0Z6IRqavLvG4YLKxtkvKptqDUGirKJqbREolUWpLiZLmFfuQ+Ro11A0ZtjRMThptxBAmdeFx+XIq8RzGgggtX57gOql6WkRTDZQnFA2X1sESmXR0a6h1GYNeVwOaku9R7qGDvT4CYRjU7YIRITaEo+OEeSdoUlWFY9WvQy6d89uCmkqLQKw52iLQKms6RoK4nQIxZ7kqSGSWVBxdArptsNTzxhKqC3Lr0VlGgjg6JX3ZIPFYFsEwX6boG4mjNEWgVJ5qnMwSE2JZ9w+w5NZUFnE3u5hjDFsj08drS+fvGsIoKbEQ/ewBoL8krjynuoLumaZvZ1p91CgD6LBabQIDs9+AZtSKiWdg8EjM3pSddHKOg70jPDwpsNsOzzI4uoivK6pWxS1pV66BrVrKL8MHgZ38dHZOhOpXmpvZzpgPBgPOFN1QSXOiYVn3vpQSk1L11CQuikGesd6+5mLWFZXwtcf3srmtgFWzpviOySuptRL93AQkycXehoI4Oiq4qmahJVL7P7FiU3upysRCKbqgoKjrQYdJ1AqKzoHpx8IXE4H//eKVezuHGZ35zAnpBgIaks9hKOGgZH8yCeWk41p8s7Q4dSu0l0e2yro3DqzcqbbIki8Zt5J0ytn3Y+TH2+5cXrvo1SBMMbQOTT9riGAy06cx1lLqnhpXy8r55cxHIxO+ZpEOV3DwaSb3GdbzloEIuIUkVdE5IFc1eGIqTaKGa1+1SwCQXyPgVQCQaJFMKgtAqUyrX8kTDhqpt0iADsd9J+uOpHm2hLObqpO6TVHFpUN5seAcS67hv4O2JLD8o8aOpxadw1A3Wo7WBwen4c8pXK85ck3rR8rESy0a0ipjEtM5ZxJIAA4fXElj3/moiOrjadyJM1EnmQgzUkgEJFFwF8BP8xF+ccI+SE4MHl6idHqV4GJQfcMxgkG21JrDYANFp4ynUKqVBZ0xK/ME5lBMy3RIsiXfEO5GiP4L+DvgQlHVkTkZuBmgMbGxszVZCjFNQRg+94HDtn7a78PC8+y91Ptex+cItX1WGW6qEypbEisKq4v87K3y5/x8qqLPYhAZ56sLs56i0BErgI6jDEvTXaeMeY2Y0yLMaalrm6STJ2zlWp6iYSSept8biZ994NtUNaQ+vllDdoiUCoLEuke6kqnXhWcDi6ng6piT960CHLRNXQBcLWI7AV+CbxRRH6eg3pYR2bypBgInC6bPXS6gcCYqfc8GEvTTCiVFZ2DQTxOB+VF2eskqS315E2aiawHAmPM540xi4wxTcC7gMeMMTdkux5HHNk6MsW+e4DShul/QSdWFU+nRVCxCPoPQmzq6WhKqZmzq4o9yDTSS8xWTYk3bzKQ6oKywTab9rm4JvXXlM2H4S6IhqdRTmIx2TRaBFVNdnVxYlxCKZURrb1+GlKc8ZMu+ZR4LqeBwBjzhDHmqlzWwa4hmA+OafwqyuYB5uj4QqrlwPRaBFVL7G3fvtRfo5Satn3dfppqUpjWnUY1JR5tEeSNwUNQPo0vZ7BdQzC9cYLprCpOqGqyt717U3+NUmpaRkJR2gcCNNcWZ7Xc2lIPg8EIgXDuu341EAy2T+/LGaC0zs4cms44wXRWFSdULLbl9GqLQKlM2ds9DMCSLLcI8mlRmQaCgWlO6QQ7plBSN70WwXRWFSc43VC+SFsESmXQvnggaK7NctdQIt9QHqSZKOxAEBq2G81MNxAAlC+A/tbUzx9sm95AcULVEh0jUCqD9sQXkC2pyX7XEJAXG9QUdiCYyQBuQuUSOyU00J9iWSlmOB2raom2CJTKoH3dw9SWelLeqzhdjmQgzYMB4wIPBDPot0+oTMzo2Z96WTMJOFVNtlsplPll70oVoj1dw1mfMQSjMpDmwRTSAg8E8RZB+YLpv7ZioR3ITSUQzGRVcUJVs71NNeAopaZlb/dw1geKAYo9Loo9zryYQlrggWAWLQKnB8oWpNZ/P5NVxQmVupZAqUzxhyIcHghmfepoQm1pfiwqK+xAMNAW36u4fGavr2q0V+qx2OTnzWRV8ZEymuytjhMolXb7um2Xa1OWZwwl1JTmx6Kywg4EiX77meYXqVwCkQB075yinFkMSpfU2mClgUCptNvbZaeO5mKMAPKnRVDYexbPdAA3oTK+T8LBdVB3wsTn9ew+9vzpELGtAl1UptQ4d64dP3b27jWp/z/b252bqaMJtaUeXtnfl5OyR9MWwUzGBxJK54HTCwcn3VrB7nHsKbPZRGeiqklbBEplwN6u3EwdTagp8dIzHCQaMzkpP6FwA0FiJs908wyNJg57lT9VIOjYAnUrZ9cF1bfP1lkplTZ7u3MzdTShttRDzECfP7fjBIUbCEZ6bf/+bLqGwAaC9o0QHpn4nI4tdq/jmapqgtAQ+LunPrf1JXj+e7D9jxo4lJrC3u7hnA0Uw9E0E7nON1S4YwSzGcAdrXoZ7HoU9v0Fll8y/vnhLvB3Qd3qmZdRu8Letq2H5ZcmPycShAc+Det/bqe2dm2zKTROepttuSiljtExEODwQJDl9aUpnT/b8YhkEquLOwaCnDBvwi3cM65wvyGOrCGYZSCoWQYON+x+PPnzHVvsbf0sAkHjeXYsYudjE5/z3HdsEDj/E3DpP8PSi2Hv07DxtzMvV6nj2KNb7e6EF63M4J7oU1hUZTfDae3NbeaAXGxev1hEHheRzSKySUT+Ltt1AGa3mGw0lxcaz4VdTyR/vnOrvZ1NIPAUQ9MFsPPPyZ8faIOnvgmrroLLvwJuH6y+GprfAPuePTprSSl1xKNbOlhYWcTKHF6JN1T4cDmEA4UWCIAIcIsx5kTgXOCjInJi1muRrhYBwNKL4PBrR/c/Hq1jC3grZl/O8kttd0+yVBOPftluaXn5V48eE4GVb7aL5bbcp+MFSo0SCEd5Zmcnl66uz+o+xWO5nA4WVBaxv2eSMcYsyMXm9W3GmJfj9weBLcDCbNeDgTYoqrJXz7O17GJ7u/vJ8c8lBopn+8e2/DJ7u/PRY4+3roMNd8F5H4Pq5mOfc3lh5ZV26mn7q7MrX6njyF92dREIx3jj6hms9k+zxupi9vcUXovgCBFpAs4A1ma98MH29LQGABpOB18l7BrTh28MdG6BulnMGEqoXQEVjcd2D0Uj8MAn7Z7LF346+esWnWOf33I/RHK/lF2pfPDnLR2UeJycu7Q611VhcXUxB3IcCHI2a0hESoHfAp80xgwkef5m4GaAxsbZjcwn1btnZit9k3E4bffQ7sftl3/i6n+ow05TrU9Dz5eInZX02t32C93lgbXfg/bX4B0/Be8E/ZwOJ6x+C7z4A3jhNjj/Y7Ovi1Ikn0UDs59Jk2nGGB7b0sGFK+rwupy5rg6N1cX0DIcYCkYo9ebmKzknLQIRcWODwC+MMb9Ldo4x5jZjTIsxpqWuLs2j+tGIzQ9UO0laiOladrEdd+jcdvRYZ2LGUBpaBGDHCUKDtiuofSM8/q9wwpV2YHgy9SfaVsmTX0s+jqFUAXntYD/tAwEuWV2f66oAsLjazhzKZasgF7OGBPgRsMUY861slw/YVbrRkF3tmy7LLwNxwrrbjx7riM8Yms0agtGaX28Hf+//BNx6ASCw+Bx46Q5Y9+OjP2OJ2PUEYb8dWFaqQAXCUT7329co97m4NA/GB8C2CICcjhPkoh1yAfBe4DURWR8/9g/GmAezVoOu7fY2nS2CioVwxnvgpR/DBZ+wexVsfQCKa6E0TVcevnL4uw1w6GVo2wDzT4P+A6m9tnQenPsR+Mt34KwbYVGLPZ4scLTcmJ76KpVHjDF84fcb2dw2wO0faKGqxJPrKgFHA0EuWwRZDwTGmGeA3M3XglGBYEV63/f1fw8bfglPfcPmB9r7NFz1n7OfMTRacbXtIkqsME72RT5Z/TbeA79+P9z8BJTmbiGNyoy52m8/G4FwlL3dw5R4XMyvSD4L8I5n9/LUjk4e29rBxSvrae8Pcufa/Xnxe6koclPmcxVWIMgLnduhpN5OH02nysVw5vttq8AY2x1zVh5dXfvK4V0/h9uvgF+/F953X65rpNSMbTzYz4+f3cPuzmGi8XUyDoFfvXiAc5qrOWVhBYFIlEN9I9zxl30MByOcvKD8mLGBiQJnNokIi6tyO4W0MANB1/b0jg+MduEt8MrP7Irlt/x3elsD6bDgDLjmu/Dbm+APn4IFZ+VfHVOl3VoF6+FN7Xzyl+txOYXzl9VwwvwyguEoB/tGCEZi/PLF/dzxl6M7By6vL+XilfU05zD6UcbYAAAaZUlEQVTB3GQaq4vZ0TGYs/ILLxAYYwPByX+dmfcvb4AP/tH2yfsqMlPGbJ1yrU198dQ3YPAwrPqrqV8zUReUfvGqLPvRM3v46h82c+qiSt588vxj9hI4cUEF717TSCgS40Cvn1Kvi4oiN797+WAOazy1xppiHt/WQSxmcDiyf2FWeEnnhjvtZvK1GWoRgL3qLl+QufdPh4v/Ec76AOx8BHY/kevaKJWSn/xlL195YDNvOnE+v7r53Ak3lPG4HCyrK2VeuQ+fO/drBaayuLqYYCRGZ462rSy8FkEmBopT7aKYzsBuponAX30LDr4Cm++1q6wz1V2mVBr8+sUDfPG+TVx24jz+591n4HYeP9exi+NZSPf3+JlXnoa0N9NUeIEgseArnVNH5yqH0055ffqwTWF94WftgLIqWLOZdRQIRzk8EMjICtmfP7+Pf7p3IxeuqOU7x1kQgFFrCbr9nN2U/bQXhRcIunaAuwTKs5/nLi85PXDW++Hpb9lgsObDua6RmmPC0RjP7+7miW2djISj/ODp3axuKOfDb1jGVac2TDu759hg9OS2Dh7efJg3rqrnf99zZl6khUi3hVVFiJCzdNQFGAi2Qe1ycBxfVxSzUtYAJ78dXv2lzW569k25rpGaI/zBCLc9vZuOwSAr6ktpaaqmvszLw5va+fhdr3DfhkN89a0n8+iW8alNpmplRGIx/vBqG2v39HDN6Qv4j+tOO+5aAglel5NFVUVsaRuXdi0rCjAQ7LA7fuW7bE+NXLzG/m62PWi33VxyfuqvjcU0sE7mOJ3mGgxHueO5vfQMh3jfeUtYNd92K757TSMff+Nybn92D9/803Yu/daTXLp6Hi1LqlJuHQyMhLnzhf3s7/Fz4Ypa/vMdp+dkNk02nb+0lgc3thGJxnBlOeAVViAY6rApGern/n/CtBOBU6+zeZjuvgk+/AyU1CQ/NxqCvc/YKaiP/wv4e6Bikd0PoaoZ/N1QXAMltVBUY/dFOPuD2f08ikN9I8T2PMUCX4ihWAml7lGbE80yEI2Eovxs7T4O9Y3wnjVHg0CCy+ng5tcv4/IT5/O5373KPa8cZMOBPt52xsIjG7ZP5NXWPu7bcIhI1PCusxdz6qLK4z4IALz+hDp+te4A6w/00ZLlcYLCCgRb7re3J1yZ23rkK5fPjhf85X/gF2+H638FZaMSc5mY3Qhn2x8g0G+//Fe8yaaq6NsPPXtg8+9t6u3RnB5Yeys0X2iT8y2/BJzJp/2lpGc3vPBDm2+pbpWdAeaa/MtlrkhHiohNh/p5ansnB3pHADuN+YvbDG9ZHOSjq4ZZUR6dVR0HAmH+5o517Okc5tqzFrG6YeIJBk21Jdz5N+dyy6838ODGNr792A4uXT2P85aOv8jY1j7If/15Ow9tbGdRVRHXnrmI+hzMoJmOdK5Mft3yWhwCT+3o0kCQUZvvhZoVs9s/eC5LZfpqxWJ4x8/g7hvhh5fC275nxxA6tsDW+2HgkN0g54z3Qc2y5FeWz30XhrvB32WDQnDAtsZe/im8+EObiG/llfCWb0+vS2mkFx74NGy6h6g4QQTnvmeIOH3sWHQtA6VLWdOS+tsdb4wx/HnLYR7b2kFtqYerTm3g4siztAU99Eg1v9zj4/f7fVw4L8jVHOD85bWEIzF6/SFe2tfL87t72Hyon0AkhsshNNWWsLyulKV1JUfef2fHEJ/69Xq2tg3yzvjV+lQcDuHs5mpOmF/GvesP8tDGdh7b2sErB/pYXl/KYCDCxoP9PLOzC5/bweUnzuPCFXU4Z9gKyIe0ETNRUezmtMWVPLW9k09flt1ZjYUTCIa7bBK413167qZUyJaVV8CND8Kd74Q7Rq06Lqq2AWDB6SCTfIG7i6Gy2OZeGi0agc7NsO0hm4ajcyu89Xupreno2Q2/eIfddvOCv+M+z1tY2P4oZf79LGl/mJX77mRvQ3629Nbu6Rl3bFc0vQnPItEY/3jPRh7b2sGZjVW87YyFOB3Cov0hFhWFWNPs4WOrh7ljZzG/2+fjs3eP37q0ubaE8iI3tS4HwXCM9Qf6eGFPDwLcu/4Q/SNh9vf48boc/OB9LbT1B455/bL9v7F3nKOuZkddKFQUuXnvuUvY1TnMhtY+HtvawT2vHMTtFOZX+Pjsm1bynjWNPPhae9p+L3PNhSvq+M5jO+jzh6gszl521MIJBFsfsF0bJ16T65pkTLIvnDXNM2xiLjgDPvws7HsWIgF7O+9UcM7iT8bpgvmnwryTofVFGxBufR288Qt22upE3UW7Hre5kUwM3ncvNF3AyNr9GIeLgdKlbG6+keWtv2Vp2x/gj5+Hy79q10jkWiQEe5+mse1hioOHCXqq8Xvr6C1P00ZFcaFIjE/c9Qp/3NTOxSvruHJ5EUv338nyA3dTPrSbmMMNe8qoXnwOn15+Pp86sZRX5r2VzYcGKPY4KfW6OHVRJfMrfMdcTUdjhtZePzs7hxgKRKgt9fB/Xr+Uy1bPG3duqkSE5fWlLK8v5R0ti4jEDF6XY8pB5NFlHQk4wLL47a7G66Zdl3z0hhNq+fajO3hmZxdXnZq97ASFEwg232sHMuefkuuaZFwkBl0hNz1hF74eF+VuQ0NxFN90vxtL6+Ckt9r74ZH0VVAcdpbSpV+C+z8Jf/oCvPILuOLf7JafiS+FSAge+wr85ds2Jcj1d0HNMnqGQ7T2+unpK6XBG2KhD7Y1Xs+S9j8x//n/tbvPvf1HM1ocl+zLbSZX7mbno4Tv/SSewf1U42YfC2gY3kY9L9PY/ggl3Zvw13yKoqXn4w/bvDj7u/08s7OLoUCYimIPdaVeFlcXTTpvPhCO8uGfv8QT2zr5f3+1ipPbfsfpT/4XnsgQXRWn0F57LhKL0uDotTPCdjyMLLmAM0+8mjMbl0z6GZwOYUlNCUtqSqb+HRiDxKKYaQTgX69rTfncQnHaokrKfC6e3q6BIP38PbDnKTjvY3O7W2iSPv6OEQcPHq7i2d5ydg37MIktH+IZNRwYGkujLC+LsqI8gjPYz4qSEeZ5w0d+JTNuPUwiEoPfbRlm/4gXf9RByDgodUa5fKmXEzb+kdIVl9tEfZt+Dz97K3grYMVldjzgwFoIDUHLBxm+6Mvcv7mPu375LBsO9MXffREAZc4Ip1UMc0nt2/jAee9CHvws/M9ZcPE/wBnvnV0rZprM4GEO330L8/fdz4FYA1+PfIrdnhNwutz4I4I31Mu1jsd5V9fjFP/8z7wWa+bO6Bt5IbaK3aYBgwMnhmj836/IEeWUJXWcvLCcUxZWHHPl/Pzubv7p9xvZ2TnE9y71cuW2D0Hri7TXrOGVlZ+it+KkI1fPDc3VMNhu80rteRq+fTqc/Tdw5nuheuk0P6SBrh0s33cv83peoGpgG0XBDtzREcLOYjhYC5VNNmVJcAi8pWn67R4lJoojFgYgJsfP15jL6eB1y2t5dGsHA4Ew5RPkUko3McZMfVa6CxW5AvhvwAn80Bjz75Od39LSYtatWzezwmIx+O0HYdM98KGnoeHU1F4327xAWcg1ZAy80uPixzuLeajVS8QITUUBzqgYosEXosodYWldOX1hB/uHnOwYcLJj0MWeQScRY79Qyl0RTirzc1r5MB842c38otj06t5yI8YYBoMReoZCtPaOsO+VR9ja7+K1Xjdb+lwEY8mDrwPDyVURzq0LcUGNnzWh5/D1brMZUYsqiS4+j921F3N7WxP3rT/IcCjKivpS3nbmQg71Bljd/yQHRrxsGSrmxb5S/FEnS2qKec/CTq7ruZWqrpeIli1gT/kaAPy+egKeGozDZYNeks+59jffpD/sZN+Il76wi4gRlp39JlbNL2Pl/LIJr85NJMSWB/+Xxle+jjsW5Oeuv8bzhlt40+lN7Hn4u0fOixloDXh4xXs2V8aeZPW+O6kN7AUg6ixi2FVBxFNBr6uOfbKQP42s5O7+VURidm/bN66sJxIzHOgd4antnayohB81P0HjttvB5WVn3WV0V5wy7oLnmCA/2A6HN9oWgonBwrPsvtYjveztixATJyb+c7j2XMREOX+xD/pboXcP7HsOhmw//rBvPt2Vp+D3zadk+ACeyADzHIN2LCcWBocbGs9lvedM2mrPo7d89eTjSxPwBbuo7V1Pbd8GFrf/mZJAGw5jZz8ZhIGSZnrLV9JVdTodVWfRV7ZiRuVkUqoty5f29fDO7z/Pxavq+f4NZ81q6qyIvGSMmXIKRdYDgYg4sdeplwGtwIvA9caYzRO9ZsaBwBh46O/hhdtsv/H5H0/9tRkMBA9v72e338euYR8dQQ/+qINAzEFThZM6X4zFJVGWlkVZVhZhQXEM56i/g3AMdg44ebzdyz37fewYcFHmjvHOpgAnetpZ6AsdU2Syq/xwDH6/dZjtQz62DRXz2mAxvWF75bGyPMKFp67gxAXlLKsrpbzITcwY/Ovv4ZDfQduIk0N+B691RukOuTgcq2AwECY25s+o1BXjxMoIp1RFKAr30VwcoNQVxeMwDESclJTX8Gqvm+c73azvcROKCU4xNBTFKKuoIhYz7O4aIhw1+NwOrjp1Adefs5gzG+2ipDvX7j+mrzgYE57vLWOT5zTW7ukmEI5yuWMdb3U+ywWOjVSIXbofxUGb1NPrXYi/+kSCRfWMOEqJhIL09vcR7GmlwgxQTy/10oeXEAOU0GdK2csChsqXUVzfTPm8Jko9QlGgk9jeZzm1836qTR8vchKPLP0ci1aciis+I2p0PRN2NV5nvxgSadFbX4T1d9HX3Y433Ic31IPD2Hz6I+4quqvP5MXAQp7oqcbrdjDPE+LNJdtZ1f80EhmBM26AutWsPZha9spdjdfx7lVO2HCnHYPp2m4z804iJk78vvl0VZ7K4ZpzOFx9DkPFi48EncTnXNNcDdGwDRqd2+zPgO0GCjuL6Ko8nf6y5QwUL8FfNJ+At5aIs4iow4szFsQZHaFkpJ1S/wGqBrdS17ue0hGbRjoqbvy+eQwVLyLktl1/rsgIxuGkamALJYHDAIRcZXRWnUl35ckMFi9hqHgRAU8lIXc5MXFjxInDRHBGAzhjQVzRAM5oMF5+4tgIMYeboLuKEW8t/qIGO94yQ9PpYrz9mT388wOb+eybVvLRi5fPuMx8DgTnAV8yxrwp/vjzAMaYf5voNTMOBE9/Ex79Z9sl9KZ/Oeap3uEQ/SNhqoo9lPlc46NukkBgDBjsVV0sfmsMxBB7LP6rjBmInnY9/SPhI/3ZuzqG2bZtExt7XbSN2KtKwVDtjhz5goyIm46Ag+HI0SsZj8OwqCSKSyBioHXYSSh+hd1SE+JtSwJcszhIqdskHSxOhTFwIOClx1nHU4c9vNjtIxSNTXi+x2Godoep9YTxVC2k3Oem2OOk2OOistjNGf2PUO2OkOqFTDAmbB8qYvNgMZ0hN/6oA4OwuChIU1GAD57ioeL8DxzzmrGBIGHNdbcQjER5tbWftv4APUNBnt5+mMWdT1IVPsy8SBsLTDtLzCHm04VXIse8PoKTIUcZIVcZUXcJRtxU1tQR7m+naGAPntj4sZKoETYUr2HklPchg20pTXs8EghGW/fjI/+GEoviC3VROnIIRFg6stGu02DU/9eiarsL3hnvsVf1o14/o/Kfv5WX9nQiJnr0B4PBgXG4CLlKQRwTDsweEwjGeGnHASqGdlM+vAenCVMxuAtXLDDuvLH83nq6Kk+jq+o0OitPp7d8Nc0H7036eTCGkpFD1PW+TH3PSyzseIKiUHcKv43UGAS/bz6DxYsZKl7EUPzW75tH2FUa/ykh4irBiP0/7oiG8EQGcIcHecslb0i9LGP4+F2v8OBrbfzspjVcsLx2RnVONRDkonNtITB6x/VWYE1GSqpeZvuIL/vKuKfuf/UQ/+/eTYC9oHHEr2oSgdEYu5+vmen2yvf/+ZiHToewpMTJObVhKkwXS4sDNBcH8DmPDcTGwEDEyaGAh0NBD20BDx1B95F6nlQXpqk4yAklI9R7bR/pplmOuYlAY1GQ65r9fGiln1AM9g852TXoxB8RHAI+p2FBcYyGoig1XsMLexNfOAeOfbMRYJqz3rwOwynlfk4pT55wq8IzvbELr8t5TAZHj8vJsv0u7J+eTTa4pLma2BnvJzTcgys0gMPjA5ePlx64fVy3yoLrbrF3YjEYOEisr5VA9378EQgW1TO/+STOLLPbH679zTenVdeJGIeTEd88Rnzz2NV4HUvXNEJoGLp3gcMFnhK758UMF+Yt2/+bY6d5Ari8RFzFaaj9eBFXCd2Vp9BdeUr8SztGUbCL4kA7vmAPCzoexxGLEHO4iImbkLuMoKeaqPPogrLKwe1UDm6f+POM0l1xEt0VJ+GIhfGGevCG+3BFA7giI/RUnoSYGDFxEXV6iTp81PatJyauI+UnbsVEcUf9uCNDeEO9hN1llPpb40Fm4qBrEIQxF9mvbwd3UUq/LxHha28/lRKPixXz0j/GMq68HLQIrgWuMMb8Tfzxe4E1xpiPjTnvZuDm+MOVwLYZFlkLdM3wtXOVfubjX6F9XtDPPBNLTOKqdhK5aBEcBEavNFoUP3YMY8xtwG2zLUxE1qXSNDqe6Gc+/hXa5wX9zJmUi2H1F4EVItIsIh7gXcB9OaiHUkopctAiMMZERORjwMPY6aO3G2M2ZbseSimlrJysxDDGPAg8mKXiZt29NAfpZz7+FdrnBf3MGZOTBWVKKaXyR34tvVNKKZV1x00gEJErRGSbiOwUkc8led4rIr+KP79WRJqyX8v0SeHzflpENovIqyLyqIhMnmFsDpjqM4867+0iYkRkzs8wSeUzi8g74v/Wm0TkzmzXMd1S+NtuFJHHReSV+N/3m3NRz3QRkdtFpENENk7wvIjIt+O/j1dF5My0V8IYM+d/sIPOu4Cl2OVMG4ATx5zzt8Ct8fvvAn6V63pn+PNeDBTH739kLn/eVD9z/Lwy4CngeaAl1/XOwr/zCuAVoCr+uD7X9c7CZ74N+Ej8/onA3lzXe5af+fXAmcDGCZ5/M/AQIMC5wNp01+F4aRGcA+w0xuw2xoSAXwJjNx64BvhJ/P7dwCWS6k7a+WfKz2uMedwYk1iq+zyJVJ1zVyr/xgBfAb4GTJ2/IP+l8pn/D/BdY0wvgDGmI8t1TLdUPrMBEjnGK4BDWaxf2hljngImyw1yDfBTYz0PVIpIQzrrcLwEgmRpKxZOdI4xJgL0AxPszp73Uvm8o92EvaKYy6b8zPEm82JjzB+yWbEMSuXf+QTgBBF5VkSej2f2nctS+cxfAm4QkVbs7MNpZJOck6b7/33ajp9E3iopEbkBaAFSz3g1B4mIA/gW8IEcVyXbXNjuoYuwrb6nROQUY0zfpK+a264H7jDGfDOexPJnInKyMWbiTIlqUsdLiyCVtBVHzhERF7ZJmb7UhNmVUpoOEbkU+EfgamNMavmJ89dUn7kMOBl4QkT2YvtS75vjA8ap/Du3AvcZY8LGmD3YFO8pbAKdt1L5zDcBvwYwxjwH+LA5eY5XKf1/n43jJRCkkrbiPuD98fvXAo+Z+EjMHDTl5xWRM4DvY4PAXO83hik+szGm3xhTa4xpMsY0YcdFrjbGzHBHo7yQyt/177GtAUSkFttVtDublUyzVD7zfuASABFZjQ0Ek2+mMLfdB7wvPnvoXKDfGNOWzgKOi64hM0HaChH5Z2CdMeY+4EfYJuRO7MDMu3JX49lJ8fN+AygFfhMfE99vjLk6Z5WepRQ/83Elxc/8MHC5iGwGosBnjTFztaWb6me+BfiBiHwKO3D8gTl8UYeI3IUN5rXxcY8vAm4AY8yt2HGQNwM7AT+QZNerWdZhDv/+lFJKpcHx0jWklFJqhjQQKKVUgdNAoJRSBU4DgVJKFTgNBEopVeA0EKi8JSKVIvK3aX7PT4pIcTrfM13liMgqEVkfz6q5TEQ+ISJbROQX2ShfFS6dPqryVjxV+APGmJOTPOeK54ya7nvuxWYl7Zp1BdNcTjzlsssY89X4463ApcaY1myUrwqXtghUPvt3YFn8KvkbInKRiDwtIvcBm8HmUhKRF+LnfF9EnPHj3xORdfEc/V+OH/sEsAB4XEQejx8bir/3JhH5s4icIyJPiMhuEbk6fo4zfs6L8XzwH4ofvyh+7t0islVEfhFf/TmunNFE5CwReVJEXhKRh0WkQWxO/U8CHxGba/9WbCrmh0TkUyJSIjZv/QvxFsM1o+r2HyKyMV63j09VvlLj5DoXt/7oz0Q/QBOjcrRjV18OA83xx6uB+wF3/PH/Au+L36+O3zqBJ4BT44/3ArWj3tMAV8bv3wP8Cbuq8zRgffz4zcAX4ve9wDqgOV6ffmzuFwfwHPC6ZOWMKs8N/AWoiz9+J3b1LNismp8Zde6R9wD+Fbghfr8Sm1OoBLvXxN3YlsToz520fP3Rn2Q/x0WKCVVQXjA2uRrYfDNnAS/G02gUAYm8Su8QkZuxaVQasBuYvJrk/ULAH+P3XwOCxpiwiLyGDUQAlwOnisi18ccV2MRuoXh9WgFEZH38Nc9MUv+V2OR4j8Tr7ARSyRtzOXC1iHwm/tgHNAKXYjdcigAYYybLa69UUhoI1FwzPOq+AD8xxnx+9Aki0gx8BjjbGNMrIndgvziTCRtjEgNlMSAIYIyJic1Smyjn48aYh8eUc1Hi/LgoU/+fEmCTMea8Kc5L9rq3G2O2janDNN9GqfF0jEDls0FseumJPApcKyL1ACJSLXZv5nJswOgXkXnAldN4z2Qexvbdu+PlnCAiJTOs+zagTmwefUTELSInpViHj0v8m19sdlmAR4APJYKWiFRPUb5S42ggUHnL2Cyaz8YHQr+R5PnNwBeAP4nIq9gvxQZjzAbsPr5bgTuBZ0e97Dbgj9McRP0hdnD6ZbEbjH+fqa/8k5Zj7PaL1wJfE5ENwHrg/BTq8BXs+MKrIrIp/jhRt/3x4xuAd09WvlLJ6PRRpZQqcNoiUEqpAqeBQCmlCpwGAqWUKnAaCJRSqsBpIFBKqQKngUAppQqcBgKllCpwGgiUUqrA/X8oL1HKr559wAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f788cd424d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# check for positivity\n",
    "ax = sns.distplot(train_pred[train_pred[\"medication\"] == 1][\"propensity_score\"])\n",
    "sns.distplot(train_pred[train_pred[\"medication\"] == 0][\"propensity_score\"])\n",
    "ax.set_xlabel('treatment effect')\n",
    "ax.set_ylabel('density')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK. We can see that our propensity score model is not overfilling, but it is also very good. When we plot the propensity score distribution for treated and untreated we can see that there is almost no overlap between them. This might give us some trouble, since is a near violation of the positivity assumption. \n",
    "\n",
    "Now, lets validate our responde prediction on both factual and counterfactual data. We can see that it is already better than the naive model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "({'r2_evaluator__recovery': 0.8686622319366993},\n",
       " {'r2_evaluator__recovery': 0.8539786349462983})"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sem_eval_fn = r2_evaluator(prediction_column=\"counterfactual\", target_column=\"recovery\")\n",
    "sem_eval_fn(train_pred), sem_eval_fn(test_pred) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, lets see the predicted average treatment effect. We can see that it is not perfect, but is much better already."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.38716813921928406"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_neg = train_pred.copy()\n",
    "X_pos = train_pred.copy()\n",
    "\n",
    "X_neg[\"medication\"] = 0\n",
    "X_pos[\"medication\"] = 1\n",
    "\n",
    "preds_rnd = sem_fn(X_pos)[\"counterfactual\"] / sem_fn(X_neg)[\"counterfactual\"]\n",
    "preds_rnd.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xt8XWWd7/HPb+/sZCdtek9L2wAtyK1cBCkVRQUFHRQGHOWlOOCMHo+MNxxmnDmjM47gmXMcHBmPRxGVQUa5jY73ijLAIOiAWEihhd6AImlJW2gbaC5NdrIvv/lj7R3SNpedZO+99sr+vl+vvLKva/1yW988z7OeZ5m7IyIitSsWdgEiIhIuBYGISI1TEIiI1DgFgYhIjVMQiIjUOAWBiEiNUxCIiNQ4BYGISI1TEIiI1Li6sAsoxoIFC3zZsmVhl1Gdnnoq+HzccQA8uaPrkJecvHR2JSsSkSqxdu3ave7eMt7rIhEEy5Yto62tLewyqtM55wSfH3gAgGWf/sUhL2m79oLK1SMiVcPMthXzOnUNiYjUOAWBiEiNUxCIiNS4SIwRiIiMJZ1O09HRQSqVCruUUCSTSVpbW0kkEpN6v4JARCKvo6OD5uZmli1bhpmFXU5FuTudnZ10dHSwfPnySW1DXUMiEnmpVIr58+fXXAgAmBnz58+fUmtIQSAi00IthkDBVL92BYGISI1TEEwjH7/9sbBLEJEIUhBME6l0lns2vRB2GSISQTpraJp4fPs+0lkPuwyRUI20xEqptY+yZEt7ezvnn38+Z555Jr/97W8544wz+OAHP8jVV1/N7t27uf322znxxBO58sor2bBhA+l0mmuuuYaLL76Y9vZ23v/+97N//34Arr/+el7/+tfzwAMPcM0117BgwQI2bNjA6aefzm233Vby8RAFwTTxyHMvhV2CSM3bunUrP/jBD7j55ps544wzuOOOO3jwwQdZvXo1X/jCF1ixYgVvectbuPnmm9m3bx+rVq3ivPPOY+HChdx7770kk0meeeYZ3ve+9w2tr/b444+zceNGlixZwllnncVDDz3EG97whpLWrSCYJh5tVxCIhG358uWcfPLJAJx44omce+65mBknn3wy7e3tdHR0sHr1aq677jogOO11+/btLFmyhE984hOsW7eOeDzO008/PbTNVatW0draCsCpp55Ke3u7gkAOlXNYu+3lsMsQqXkNDQ1Dt2Ox2ND9WCxGJpMhHo/zox/9iOPyy8YXXHPNNSxatIj169eTy+VIJpMjbjMej5PJZEpetwaLp4H9gxn601la5zaGXYqIjOEP/uAP+NrXvoZ7MJ73+OOPA9DV1cXixYuJxWLceuutZLPZitalFsE00NOfBuCc41q47XfbQ65GJDyjDeRWi7//+7/nqquu4pRTTiGXy7F8+XLuvPNOPvaxj/Hud7+bW265hfPPP58ZM2ZUtC4rJFM1W7lypevCNKM45xy2vNDDxz78z1x6xuF84ZdbDnlJtf9xiEzV5s2bOeGEE8IuI1QjfQ/MbK27rxzvveoaijgHelJpXrt8Ho31auCJyMTpyBFxqXSWbM45/ch5xGp3qRURmQK1CCIunc0BsHh2kqb6eMjViEgUqUUQcdlcMMbz0Na99A1W9kwDEZkeFAQRVwiCZF186LaIyESoayjiMvmDf0MiRiKuH6eITJyOHBFXaAU01MWpr9OPUyQs+/bt44Ybbpjw+77zne+wc+fOofvLli1j7969pSxtXDpyRFwhCBJxUxCIhGi0IBhvSYiDgyAMGiOIuGzOiZlhZtSra0gkNJ/+9Kd59tlnOfXUU0kkEiSTSebOncuWLVu45557uPDCC9mwYQMA1113Hb29vZx00km0tbVx2WWX0djYyMMPPwzA1772NX7+85+TTqf5wQ9+wPHHH1/W2hUEEZfNObH8BAK1CESAq66CdetKu81TT4WvfGXMl1x77bVs2LCBdevW8cADD3DBBRewYcMGli9fTnt7+4jvueSSS7j++uu57rrrWLnylQnACxYs4LHHHuOGG27guuuu46abbirlV3MIHTkiLus+NJEsVsMX7xapNqtWrWL58uWTeu+73vUuAE4//fRRQ6SU1CKIuELXkIjkjfOfe6UMXziurq6OXC43dD+VSo353sLS0+VadvpgahFEXDbnJb9snYhMXHNzMz09PSM+t2jRInbv3k1nZycDAwPceeedRb2vUsrWIjCzm4ELgd3uflL+sXnA94FlQDvwHnfXFVWmIJtz4lpkSCR08+fP56yzzuKkk06isbGRRYsWDT2XSCT43Oc+x6pVq1i6dOkBg78f+MAH+MhHPnLAYHGllW0ZajN7E9AL3DIsCP4JeMndrzWzTwNz3f1vxtuWlqEe3drlryaZiLP+1p8C8Lc/efKQ12gZapnutAx1lS5D7e6/AQ6+kO7FwHfzt78LvLNc+68VmZxr1VERmZJKjxEscvdd+dsvAIvGerGMbTCTw/2V00dFRCYjtMFiD/qkRu2XMrMrzKzNzNr27NlTwcqio3cgOJtAZw2JQBSutlguU/3aKx0EL5rZYoD8592jvdDdb3T3le6+sqWlpWIFRklvqhAEIRciErJkMklnZ2dNhoG709nZSTKZnPQ2Kj2PYDXwp8C1+c8/q/D+p5XuVHDRerUIpNa1trbS0dFBrfYeJJNJWltbJ/3+cp4++m/AOcACM+sAriYIgH83sw8B24D3lGv/taDQNaR5BFLrEonEpGfxShmDwN3fN8pT55Zrn7WmN5VhJuoaEpGp0cziCOsZUNeQiEydgiDCNFgsIqWgIIiwnsLpo0oCEZkCBUGE9aYymBmKARGZCgVBhPWkMvkF5xQFIjJ5CoII6x3IaOVREZkyBUGE9aQyxHXGkIhMkYIgwnpSabUIRGTKFAQRpq4hESkFBUGEKQhEpBQUBBHWm1IQiMjUKQgirCeVoU5BICJTpCCIqIFMlsFsTi0CEZkyBUFEFdYZUhCIyFQpCCKqR0EgIiWiIIiowkVpNKFMRKZKQRBRahGISKkoCCKqJ3+9YgWBiEyVgiCiCl1DOn1URKZKQRBRQ2MECgIRmSIFQURpjEBESkVBEFE9qQyJuGE6a0hEpkhBEFG9A2lmNtTp2mQiMmUKgojqG8zSVF8XdhkiMg0oCCIqlc7SWB8PuwwRmQYUBBEVtAgUBCIydQqCiOobzJJMKAhEZOpCCQIz+wsz22hmG8zs38wsGUYdUZZKq0UgIqVR8SAws6XAJ4GV7n4SEAcurXQdUdc3mKVRLQIRKYGwuobqgEYzqwOagJ0h1RFZ/YMaLBaR0qh4ELj7DuA6YDuwC+hy93sqXUfU9atrSERKJIyuobnAxcByYAkww8wuH+F1V5hZm5m17dmzp9JlVr1+dQ2JSImE0TV0HvCcu+9x9zTwY+D1B7/I3W9095XuvrKlpaXiRVazXM7pT2dp1IQyESmBMIJgO3CmmTVZsFDOucDmEOqIrIFMDkAtAhEpiTDGCNYAPwQeA57M13BjpeuIsr7BYOVRjRGISCmE0rfg7lcDV4ex7+mgbzALqEUgIqWhmcURlErng0AtAhEpAQVBBBVaBOoaEpFSUBBEUH9aXUMiUjoKggjqH1TXkIiUjoIggvo1RiAiJaQgiKChMYKEJpSJyNQpCCKo0CJI1uvHJyJTpyNJBPUPTShTi0BEpk5BEEGaUCYipaQgiKD+dJb6uhjxmIVdiohMAwqCCOrXhetFpIQUBBGkaxGISCkpCCKoL63LVIpI6SgIIiilFoGIlJDOP4yQO9ZsB+C5vfvJunPHmu2c2z0QclUiEnVqEURQOpujPq4fnYiURlFHEzP7QzPTkadKDGZzJBQEIlIixR5N3gs8Y2b/ZGbHl7MgGd9gJkd9nYJAREqjqKOJu18OnAY8C3zHzB42syvMrLms1cmI0llXi0BESqboo4m7dxNcdP57wGLgj4DHzOzKMtUmoxjM5qiPa1axiJRGsWMEF5vZT4AHgASwyt3fDrwa+FT5ypODuTvpTI6EuoZEpESKPX30XcD/c/ffDH/Q3fvM7EOlL0tGk8k5DjprSERKptijyQsHh4CZfRHA3e8reVUyqnQ2B6DBYhEpmWKPJm8d4bG3l7IQKc5gJgiCiQwWD2Sy5SpHRKaBMbuGzOyjwMeAo83siWFPNQMPlbMwGVk668DEuoZ6UxkaZmpJChEZ2XhjBHcAdwH/CHx62OM97v5S2aqSUQ1mJ94i6B3IMH9mQ7lKEpGIGy8I3N3bzezjBz9hZvMUBpVX6BqayBhB70CmXOWIyDRQTIvgQmAt4MDwk9cdOKpMdcko0kMtguLnEfSmFAQiMroxg8DdL8x/Xl7KnZrZHOAm4CSCQPkf7v5wKfcxXalFICKlVuyEsrPMbEb+9uVm9mUzO2IK+/3/wH+4+/EEk9I2T2FbNWXo9NEJjhGIiIym2KPJN4A+MyvMJH4WuHUyOzSz2cCbgG8DuPugu++bzLZq0WQGi3vUNSQiYyj2aJJxdwcuBq53968TnEI6GcuBPcC/mtnjZnZTobUxXH5RuzYza9uzZ88kdzX9DJ0+qq4hESmRYo8mPWb2GeBy4Bf5axMkJrnPOuA1wDfc/TRgPweemgqAu9/o7ivdfWVLS8skdzX9TGZCWVd/ulzliMg0MJHrEQwAH3L3F4BW4EuT3GcH0OHua/L3f0gQDFKEdDZH3Ix4rPizhl7qHSxjRSISdUUtOpc/+H952P3twC2T2aG7v2Bmz5vZce7+FHAusGky26pFg9kcibqJLUHduV/XNRaR0RUVBGb2LuCLwEKCuQRGMNls1iT3eyVwu5nVA78HPjjJ7dScdGbi1yveqxaBiIyh2GWo/wn4Q3cvyWme7r4OWFmKbdWayVyvWC0CERlLsUeUF0sVAjI16Ulcr7hTLQIRGUOxLYI2M/s+8FOCQWMA3P3HZalKRjWZ6xX3DWbpG8zQVF/sj1tEakmxR4ZZQB/wtmGPOaAgqLDgesUTvyhNZ+8gTfMUBCJyqGLPGtJgbpVIZ3PMaJj4Ab1z/yCHz2sqQ0UiEnXFrjV0rJndZ2Yb8vdPMbPPlrc0GclgJjehlUcLOns1YCwiIyu2j+FfgM8AaQB3fwK4tFxFyejSU+gaEhEZSbFHlCZ3f+Sgx7SATQiCCWWTCIL9CgIRGVmxR5S9ZnY0wQAxZnYJsKtsVcmI3J3BSUwoa6qPq2tIREZV7Kjjx4EbgePNbAfwHHBZ2aqSEaWzTs4hmZjYhejnz6xXi0BERjVmEJjZXw67+0vgfoJWxH7g3Qxbf0jKrz+dBaBxokEwo4G9ahGIyCjGaxEUrjlwHHAG8DOCdYbeDxw8ZiBlNhQE9RMLggUz69m5L1WOkkRkGhjvmsWfBzCz3wCvcfee/P1rgF+UvTo5QGowCIJkYmJjBPNnNPDkjq5ylCQi00CxR5RFwPBO5sH8Y1JBk+4amllPZ+8gwUXmREQOVOxg8S3AI2b2k/z9dwLfKUtFMqrJBsG2zj4yOefmB9tprI/zx689ohzliUhEFbvExP81s7uAN+Yf+qC7P16+smQkqUkGQWFJit6BzITHF0Rk+it60Rp3fwx4rIy1yDj6C2MEEzyYzxwWBC3NDSWvS0SibeJTVCU0/eksDXUxYjaxtYZmNATBsX9Ak8FF5FAKgghJpbMT7haCA1sEIiIHUxBESP9gdlJ9/IUL0qhFICIjURBESH86N+HlJQDiMaOpPq4WgYiMSEEQIZPtGoLgzCG1CERkJAqCCOmfQhDMbKijdyBb4opEZDpQEERIfzo74eUlCtQiEJHRKAgiIp3NMZjJTXpCWHNDHd2pdImrEpHpQEEQEd39wUF8sl1DsxsTDGRyDKTVPSQiB1IQRER3KujWmcxZQxAEAcC+frUKRORACoKI6Cq0CCbZNTQrHwTdCgIROUhoQWBmcTN73MzuDKuGKOmaYtfQnHwQdCkIROQgYbYI/hzYHOL+I6Xwn/xku4ZmNSYw1DUkIocKJQjMrBW4ALgpjP1H0VS7huIxY2ayTi0CETlEWC2CrwD/C8iFtP/ImWrXEAQDxhojEJGDVTwIzOxCYLe7rx3ndVeYWZuZte3Zs6dC1VWv7lSaupiRiE/+Rza7MaGuIRE5RBgtgrOAi8ysHfge8BYzu+3gF7n7je6+0t1XtrS0VLrGqtPdn55SawCCIOjqT+vaxSJygIoHgbt/xt1b3X0ZcCnwK3e/vNJ1RE13f2bSA8UFsxsTDGZyQ3MSRERA8wgio6s/PeXrDRcmlb3QlSpFSSIyTYQaBO7+gLtfGGYNUdFVoq4hgJ1d/aUoSUSmCbUIIqI7lZ70yqMFhSDYtU8tAhF5hYIgIkrRNdScDCaV7VKLQESGURBEQC7nJTlrKB4zmpN17NIYgYgMoyCIgP2DGXI++eUlhpvdmFCLQEQOoCCIgFLMKi6Y3ZjQGIGIHEBBEAFdU1xwbrjZjQl2dvVrUpmIDFEQREB3fzABbKqDxQCzm+pJpXNafE5EhigIIqDUXUMAO9U9JCJ5CoII2Nc3CJSmRTCvqR6Ah7bunfK2RGR6UBBEwPMv91EXM2YlE1Pe1pI5Sd58XAv/fO9T/H5PbwmqE5GoUxBEQHtnH0vnNhKP2ZS3ZWZc++5TqI/H+OsfPkE2p0FjkVqnIIiA7Z19HDGvqWTbWzQryecvPpG1217mXx96rmTbFZFoUhBEwLbO/SybP6Ok23znqUs57Yg5rF6/s6TbFZHoURBUuX19g3SnMhw5v3QtAgi6iFYtn8fmXd0MZLIl3baIREtd2AXI2LZ19gFwxLwm9vYOlmSbd6zZDgTzE9JZ5yv3PsPh85r449ceUZLti0i0qEVQ5ba9FATBkSXuGgI4fG4jAB37tPaQSC1TEFS57Z37AUo6WFwwuzHBjIY6drzcV/Jti0h0KAiqXHtnHwubG0oymexgZkbrnEaef1ktApFapiCocts7+0o+UDxc67xG9vYMMJDWgLFIrVIQVLltL+3niHmlHx8oaJ3ThAM7NE4gUrMUBFUslc7yYvdAeVsEhQFjdQ+J1CwFQRXbPnTGUPmCYEZDHXObEnRowFikZmkeQRUqnOe/eVc3AFt29bB/YHvZ9tc6t4nnFQQiNUstgirWuT+YQDZ/Rn1Z99M6t5F9fWle7NY1CkRqkYKginX2DpBMxMpy6uhwxy5qBuBXW3aXdT8iUp0UBFWss3eQeTPqMZv68tNjWdjcwNymBP+56cWy7kdEqpOCoErt6urn2T29HL1gZtn3ZWacsHgWD27dS99gpuz7E5HqoiCoQu7OXRteIJmIc/ZxLRXZ5wmLZzGQyfHgM7qEpUitqXgQmNnhZna/mW0ys41m9ueVrqHaPf1iL1t39/Lm4xfSVF+ZE7uWzZ9Bc7KO/9ys7iGRWhNGiyADfMrdVwBnAh83sxUh1FEVtnXup3/wleUdMtkcd23YxbwZ9Zx51LyK1RGPGW8+biH3bd6ty1eK1JiKB4G773L3x/K3e4DNwNJK11EN0tkcF371Qa69a/PQY//1zF529wzwthWLqItV9sdz3opFdO4fZN3z+yq6XxEJV6hjBGa2DDgNWDPCc1eYWZuZte3Zs6fSpVXEMy/20jOQYfX6naSzOQB+/sROkokYKxbPqng9Zx/bQiJu3L3xhYrvW0TCE1oQmNlM4EfAVe7effDz7n6ju69095UtLZUZMK20DTu7AHi5L82Dz+wllc5y78YXOXHxbOrilf/RzG5M8KZjWvj5+p3k1D0kUjNCCQIzSxCEwO3u/uMwaqgGm3Z201QfZ3ZjgtXrd/Lrp/fQM5Dh5NbZodV00alL2NWV4tH2l0KrQUQqq+JrDVkwO+rbwGZ3/3Kl919N7t+ym5aZDSyc1cAvntjFlhd6aKqPc3RL+ecOjOatKxbRmIjzs/U7ee1R80OrQ0QqJ4xF584C3g88aWbr8o/9rbv/MoRaQpPLObu6Upx+5FxWLJnFo+0vs3lXN6uWzSMeK+9M4tEUFrs7ZtFMfvLYDo4/rJm6WEwXtReZ5ioeBO7+IBDOka6KPNe5n8FsjiVzkixfMINZyTq6U+F2CxWc2jqHJzq62Lq7l+MPq/ygtYhUlmYWh2TDjmCgeMmcRmJmrFw2jwUz61m+oHxXIyvWqxbNpDERZ71OIxWpCboeQUg27ewmHjMWNicBOPf4hZx7/MKyLzBXjLpYjJOWzmL9811Dp7WKyPSlFkFINuzs4rBZyaHxADOrihAoOHHJbAazOZ7d0xt2KSJSZgqCELg7G3Z0s2ROMuxSRnVUywwa6mJs2nnIFA8RmWYUBCHYsa+frv40S+Y0hl3KqOpiMY47rJnNu7q19pDINKcgCMH65/MDxbOrNwgg6B7aP5ilTZPLRKY1BUEIVq/fwfwZ9VXdIgA4duFM6mLG3Ru1NLXIdKYgqLC9vQPct3k373rN0tAmjhWrIRHMcr574wu4q3tIZLpSEFTYTx7bQSbnvGfl4WGXUpQTl8xix75+PvPjJ3lo614yOp1UZNpREFSQu/Pvbc9z2hFzOGZRc9jlFOWU1jm889Ql/GzdTi67aQ1/dutatQ5EphkFQQWte34fz+zu5b0RaQ0A1NfF+Mqlp/H4597Kp956LPdt2c0dj2wPuywRKSHNLK6gf2/roDER54JTFoddyoQUFqObO6OeV7XM5POrN/FS7yBXnntMyJWJSCmoRVBBD23dy9nHttCcTIRdyqTEzHjXa5YSi8H3256nrf0lXcBGZBpQi6AC7liznVQ6y/aX+jj+sOah/7CjaE5TPe88dSk/XNvBJd98mEWzGvjqpafp2gUiEaYWQYXs7OoHqPq5A8U4pXUOf/eOE/jq+04jbsYXfrlZA8giEaYgqJCd+1IALJ5dvesLTURDIs5Fr17CR9/8KtZ3dPHIc5p9LBJVCoIK2bWvn+ZkXWTHB0Zyx5rt5HJOU32cq1dv5I412yPd7SVSqxQEFbKzq7/q1xaajEQ8xuuOms+WF3p4sTtFJpdj6+4edRWJRIgGiysgnc2xp2eAFYun52UfX3vUfH799B5u/d02egcyDGZyfPaCE/ifbzwq7NJEpAhqEVTAC10pcj49BopHMrOhjjceswB359TWOZx9bAv/eNcWfvvs3rBLE5EiqEVQAUNnDE3DrqGCt644jLeuOAyAgXSWTTu7+fB327h01RHUxYzmZIJPvOVVIVcpIiNREFTAzn0pGhNx5jRNn4HisTQk4lx+5pHc8MBWvv3gcwAYcHLrbM4+tiXc4kTkEOoaqoBdXf0snpOsqmsSl1tLcwOfPPcY/vR1R/KB1y9j0awkV97xGM/t3R92aSJyEAVBmaWzOV7oSk3rbqHRzG2q57jDZnHsomYuP/NI4jHjw7e00Z1Kh12aiAyjICizwvUHWufWXhAMN29GPV+/7DU8t3c/f/LtR9jXNxh2SSKSpyAoow07uvjszzZwdMsMTlwyO+xyQte+t49LzzicJ3d0cf5X/otv/fpZTUATqQIKgjLZ1zfIR25by/wZ9bz3jCOq/rKUlXLiktn86euW8dL+Qa6/fytrt72sFUxFQhZKEJjZ+Wb2lJltNbNPh1FDOe3uTnH5t9ewu3uAb1x+OjMbdHLWcK9aOJMPv+koZjcm+NFjHVzwtQe5fc02Xt6v7iKRMFT8CGVmceDrwFuBDuBRM1vt7psqXUup3fa7bezal+K2NdvoH8xy6arD2bSzO+yyqtLSOY185OyjebKji7ZtL/F3P9nA1T/byIWnLOavzz+epdN08p1INQrjX9VVwFZ3/z2AmX0PuBgoWxC4O+4QG9Y9k8s5OXfq4rGh+33pLO7OjPo6YjFjIJOluz9DfTxGc7KOwWyOp1/sYcuuHnZ1pdjdk2J3z0Dw0Z3ixe5gBnFzso4r3nTUtJ1JXCoxM159+BxOaZ3Nrq4Uj29/mTuf2MWdT+zi9CPn8sZjFpBMxGmqr6OpPk5zso4FMxuYN6OehroYdfEY8ZhRF7MDPgPkPPi55xyc4OfvDmbBR8wMy9cQPKauO6k+uZxX5PczjCBYCjw/7H4H8Npy7Oijt63lV1t2M5jNBUFgwTV4czkYzOYAiMeMRNwYyASvKaiPx4ZeA/kDCMEBpqBwcGpOJjhsVpLjDmtmdmOCFYtnTatVRsvNzFgyp5Elcxo561ULuGfTi7S1v8yaCi5tXQiHWP6PLpa/7wcFieNDIePAVNbWG+tve6w/+9EOCmO/Z6x9jfLkBOo7ePsHb3Os/Q//Hgbf1dGeG9vwXQzf33i1FPZR+DkfsC8/tKYDtmmv7NfslecKt4f/zuCQG/q9CbY59I8INvRPimE4TjrrZHPOfZ86m6NbZo7z1U9N1XZem9kVwBX5u71m9lSY9YxgAVA9i+mceeRozyy47ItVVOfYqut7Orao1BqVOkG1juhVX5zS20c9MAwXRhDsAA4fdr81/9gB3P1G4MZKFTVRZtbm7ivDrmM8UakTVGs5RKVOUK1hCuOsoUeBY8xsuZnVA5cCq0OoQ0RECKFF4O4ZM/sEcDcQB252942VrkNERAKhjBG4+y+BX4ax7xKq2m6rg0SlTlCt5RCVOkG1hsZ0SUERkdqmJSZERGqcgmAM4y2FYWZ/aWabzOwJM7vPzIo6Vascil22w8zebWZuZqGd8VBMrWb2nvz3dqOZ3VHpGvM1jPfzP8LM7jezx/O/A+8Iqc6bzWy3mW0Y5Xkzs6/mv44nzOw1la5xWC3j1XpZvsYnzey3ZvbqStc4rJYxax32ujPMLGNml1SqtpILZt3q4+APgoHsZ4GjgHpgPbDioNe8GWjK3/4o8P1qrTX/umbgN8DvgJXVWitwDPA4MDd/f2GV1nkj8NH87RVAe0jf0zcBrwE2jPL8O4C7COY+nQmsCaPOImt9/bCf+9urudZhvye/IhjzvCSsWqf6oRbB6IaWwnD3QaCwFMYQd7/f3fvyd39HMCciDOPWmvcPwBeBVCWLO0gxtX4Y+Lq7vwzg7rsrXCMUV6cDs/K3ZwM7K1jfK0W4/wYYaxr2xcAtHvgdMMfMFlemugONV6u7/7bwcyfcv6livq8AVwI/AsL4HS0ZBcHoRloKY+kYr/8QwX9dYRi31nx3wOHu/otKFjaCYr6vxwIc7ChqAAAFaElEQVTHmtlDZvY7Mzu/YtW9opg6rwEuN7MOgv8Ir6xMaRM20d/lahHm39S4zGwp8EfAN8KuZaqqdomJKDGzy4GVwNlh1zISM4sBXwY+EHIpxaoj6B46h+A/wt+Y2cnuvi/Uqg71PuA77v7PZvY64FYzO8ndc+O9UcZmZm8mCII3hF3LGL4C/I2756K+aKGCYHRFLYVhZucBfwec7e4DFartYOPV2gycBDyQ/4U9DFhtZhe5e1vFqgwU833tIOgbTgPPmdnTBMHwaGVKBIqr80PA+QDu/rCZJQnWoKm2boKifperhZmdAtwEvN3dO8OuZwwrge/l/6YWAO8ws4y7/zTcsiZOXUOjG3cpDDM7DfgWcFFI/dgFY9bq7l3uvsDdl7n7MoK+1zBCYNxa835K0BrAzBYQdBX9vpJFUlyd24FzAczsBCAJ7KlolcVZDfxJ/uyhM4Eud98VdlEjMbMjgB8D73f3p8OuZyzuvnzY39QPgY9FMQRALYJR+ShLYZjZ/wba3H018CVgJvCD/H8F2939oiqttSoUWevdwNvMbBOQBf660v8ZFlnnp4B/MbO/IBg4/oDnTyWpJDP7N4LgXJAfr7gaSOS/jm8SjF+8A9gK9AEfrHSNBUXU+jlgPnBD/m8q4yEt7lZErdOGZhaLiNQ4dQ2JiNQ4BYGISI1TEIiI1DgFgYhIjVMQiIjUOAWBVC0zm2NmHyvxNq8ys6ZSbrNU+zGz481sXX4106PN7JNmttnMbq/E/qV26fRRqVpmtgy4091PGuG5OnfPTGKb7QQrr+6dcoEl3k9+qes6d/8/+ftbgPPcvaMS+5fapRaBVLNrgaPz/yV/yczOMbP/MrPVwCYI1nkys0fyr/mWmcXzj3/DzNosuJ7B5/OPfRJYAtxvZvfnH+vNb3ujmf2nma0yswfM7PdmdlH+NfH8ax7Nr5X/Z/nHz8m/9odmtsXMbs/P3j1kP8OZ2elm9mszW2tmd5vZYguuZXAV8FELrnHwTYIlsO8ys78wsxkWrI//SL7FcPGw2q4zsw352q4cb/8ihwh7HWx96GO0D2AZw9aCJ5jluR9Ynr9/AvBzIJG/fwPwJ/nb8/Kf48ADwCn5++3AgmHbdII1bQB+AtxDMHv01cC6/ONXAJ/N324A2oDl+Xq6CNbuiQEPA28YaT/D9pcAfgu05O+/l2DWMgSrmf7VsNcObQP4AnB5/vYc4GlgBsF1MH5I0JIY/nWPuH996GOkDy0xIVHziLs/l799LnA68Gh+OYJGXlnw7T1mdgXBMiqLCS4c88QI2xsE/iN/+0lgwN3TZvYkQRABvA04xV65AtVsgkXwBvP1dACY2br8ex4co/7jCBYAvDdfcxwoZt2ftwEXmdlf5e8ngSOA84Bver6bzN3HWz9f5BAKAoma/cNuG/Bdd//M8BeY2XLgr4Az3P1lM/sOwYFzJGl3LwyU5YABAA+WFi78fRhwpbvffdB+zim8Pi/L+H9TBmx099eN87qR3vdud3/qoBomuBmRQ2mMQKpZD8ES2qO5D7jEzBYCmNk8C64bPYsgMLrMbBHBJQ+L3eZI7ibou0/k93Osmc2YZO1PAS0WXL8AM0uY2YlF1nCl5Y/8Fqx8C3Av8GeF0DKzeePsX+QQCgKpWh6sOPpQfiD0SyM8vwn4LHCPmT1BcFBc7O7rCa55vAW4A3ho2NtuBP5jgoOoNxEMTj9mwYXMv8X4//mPuB8PLnt5CfBFM1sPrCO4Tu94/oFgfOEJM9uYv1+obXv+8fXAH4+1f5GR6PRREZEapxaBiEiNUxCIiNQ4BYGISI1TEIiI1DgFgYhIjVMQiIjUOAWBiEiNUxCIiNS4/wbfnDIDGx8SGgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f788ca8c390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = sns.distplot(preds_rnd)\n",
    "ax.set_xlabel('treatment effect')\n",
    "ax.set_ylabel('density')\n",
    "plt.axvline(np.mean(preds_rnd), label='mean', lw=5)\n",
    "plt.axvline(np.exp(-1), color='r', label='truth')\n",
    "plt.legend();"
   ]
  }
 ],
 "metadata": {
  "nbsphinx": {
   "orphan": true
  },
  "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}