Skip to content

Instantly share code, notes, and snippets.

@akelleh

akelleh/front_door_criterion.ipynb

Created March 22, 2018 18:16
Show Gist options
  • Save akelleh/2b034010b61773052aaa6e981dc46f9b to your computer and use it in GitHub Desktop.
Save akelleh/2b034010b61773052aaa6e981dc46f9b to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
front_door_criterion.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"from sklearn.linear_model import LogisticRegression, LinearRegression"
]
},
{
"cell_type": "code",
"execution_count": 184,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"N = 5000\n",
"u = np.random.normal(size=N)\n",
"\n",
"p_d = 1./(1. + np.exp(-u))\n",
"d = np.random.binomial(1, p=p_d)\n",
"\n",
"p_z = 1./(1. + np.exp(-(d + 0.1 * np.random.normal(size=N))))\n",
"z = np.random.binomial(1, p=p_z)\n",
"\n",
"y = 2. * u + 2.*z + 0.2* np.random.normal(size=N) \n",
"df = pd.DataFrame({'D': d, 'Y': y, 'Z': z})"
]
},
{
"cell_type": "code",
"execution_count": 185,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"D 1.000000\n",
"Y 2.279199\n",
"Z 0.724349\n",
"dtype: float64"
]
},
"execution_count": 185,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df[df['D'] == 1].mean()"
]
},
{
"cell_type": "code",
"execution_count": 186,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"D 0.000000\n",
"Y 0.136311\n",
"Z 0.501965\n",
"dtype: float64"
]
},
"execution_count": 186,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df[df['D'] == 0].mean()"
]
},
{
"cell_type": "code",
"execution_count": 187,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEICAYAAABfz4NwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAFthJREFUeJzt3X+s3Xd93/Hnq4QCclmTYHZlORaX\nDmtVurQh9ZJMdN2lKcEJHc5UxmCocVA6CzWsIPkPvFZburBKYRKghXVhhnh2KEtIaVFMEy11A1cZ\nWgNxVhrnx5gNSxR7Ji5xcHHYaL2998f5OJzd3Ot7zvW551eeD+nofs/n+z3f+z6f+/2e1/3+PKkq\nJEn6kVEXIEkaDwaCJAkwECRJjYEgSQIMBElSYyBIkgADQZLUGAgTKMmTSf5Xku8l+W6S/5LkfUn8\ne2oqJfndJP9hQdvfS/JsknWjqmva+AEyuf5+Vb0aeB1wM/Ah4LbRliStmg8AVyV5C0CSVwKfArZX\n1dGRVjZFDIQJV1Unqmov8I+ArUn+1qhrkgatqp4F/imwM8ka4Ebgm1W1e6SFTZlzRl2ABqOqvpbk\nMPB3gUdHXY80aFX1e0neBdwBvAm4eMQlTR0DYbr8T+D8URchraJfA74J/GZVPT3qYqaNu4ymy3rg\n+KiLkFZLVT0DfAd4bNS1TCMDYUok+dt0AuEro65F0mQyECZckr+W5JeAO4HfraoDo65J0mTyGMLk\n+mKSU8D/BR4HPgZ8crQlSZpk8QtyJEngLiNJUmMgSJIAA0GS1BgIkiRgzM8yWrt2bc3Ozi467vnn\nn2fNmjXDLWgM2Q8dZ+qHhx9++DtV9dohl7RiZ1ruV2JSlhHrHJyVLvNjHQizs7Ps379/0XHz8/PM\nzc0Nt6AxZD90nKkfkjw13GrOzpmW+5WYlGXEOgdnpcu8u4wkSYCBIElqDARJEjDmxxDO5MCRE1y3\n456+XvPkzW9bpWokrcRs1zq8/aJTPa3Trserxy0ESRJgIEiSmondZaTpNdvnrkCA3ZvH+7xwDc5K\nlg9wV1Mv3EKQJAEGgiSpMRAkSYCBIElqDARJEmAgSJIaA0GSBBgIkqTGQJAkAQaCJKkxECRJgIEg\nLSrJriTHkjza1XZ+kn1JDraf57X2JLklyaEkjyS5pOs1W9v0B5NsHcV7kXq1bCAk2ZDky0keT/JY\nkg+0dlcOTbPdwOYFbTuA+6tqI3B/ew5wFbCxPbYBt0JnHQFuBC4DLgVuPL2eSOOoly2EU8D2qroQ\nuBy4IcmFuHJoilXVA8DxBc1bgD1teA9wTVf77dXxIHBuknXAW4F9VXW8qp4D9vHikJHGxrK3v66q\no8DRNvy9JE8A6+msBHNtsj3APPAhulYO4MEkp1eOOdrKAZDk9MpxxwDfj7SaZtr6APBtYKYNrwee\n7prucGtbqv1Fkmyj8w8UMzMzzM/PD6zokydPDnR+g7T9olMvDM+86v9/PmiD6oNx7s+z1df3ISSZ\nBd4IfJVVWjl6XTFWsvBM4x9xGhfOlXwoDLsfqqqS1ADntxPYCbBp06aam5sb1KyZn59nkPMbpOsW\nfIXmRw+s3le0PPmeuYHMZ5z782z13PtJfgz4feCDVfUXSV4YN8iVo9cV4xOfvbvvhWdQC8Q4mcaF\ns9/vyobOF+QMoR+eSbKuqo62rd5jrf0IsKFrugta2xF+uBV9un1+tYuUVqqns4ySvJxOGHy2qv6g\nNT/TVgr6WDkWa5cmxV7g9MkQW4G7u9qvbSdUXA6caFvP9wFXJjmvHS+7srVJY6mXs4wC3AY8UVUf\n6xrlyqGpleQO4E+Av5nkcJLrgZuBtyQ5CPxiew5wL/At4BDwKeDXANrxsg8DD7XHTaePoUnjqJd9\nLm8CfgU4kOTrre036KwMd7UV5SngnW3cvcDVdFaO7wPvhc7KkeT0ygGuHBpjVfXuJUZdsci0Bdyw\nxHx2AbsGWJq0ano5y+grQJYY7cohSVPCK5UlSYCBIElqDARJEmAgSJIaA0GSBPR56wpJmlSzK7gC\n/smb37YKlYwvtxAkSYBbCJIGZCX/gWu8uIUgSQIMBElSYyBIkgADQZLUGAiSJMBAkCQ1BoIkCTAQ\nJEmNgSBJAgwESVJjIEiSAANBktQYCJIkwECQJDUGgiQJMBAkSY2BIEkCDARJUrNsICTZleRYkke7\n2s5Psi/JwfbzvNaeJLckOZTkkSSXdL1ma5v+YJKtq/N2pNWV5MkkB5J8Pcn+1tb3+iCNo162EHYD\nmxe07QDur6qNwP3tOcBVwMb22AbcCp0VBrgRuAy4FLjx9EojTaA3V9XFVbWpPe9rfZDG1bKBUFUP\nAMcXNG8B9rThPcA1Xe23V8eDwLlJ1gFvBfZV1fGqeg7Yx4tDRppU/a4P0lg6Z4Wvm6mqo23428BM\nG14PPN013eHWtlS7NGkK+KMkBfz7qtpJ/+vDURZIso3OVgQzMzPMz88PrOCTJ08OdH5L2X7RqbN6\n/cyrzn4eg7ZYvw2rP0dhpYHwgqqqtnIMRK8rxkoWnmn8I07jwrmSD4Uh9sPPVdWRJH8d2Jfkv3WP\nXOn60IJlJ8CmTZtqbm5uIMVCZ7kf5PyWct2Oe87q9dsvOsVHD5z1R9JAPfmeuRe1Das/R2Glvf9M\nknVVdbRtAh9r7UeADV3TXdDajgBzC9rnF5txryvGJz57d98Lz2J/3Ek3jQvnSj5Ydm9eM5R+qKoj\n7eexJF+gc0ys3/VBGksrPe10L3D6TKGtwN1d7de2sysuB060Ten7gCuTnNcOJl/Z2qSJkWRNklef\nHqazHD9K/+uDNJaW/Rc7yR10/rtfm+QwnbOFbgbuSnI98BTwzjb5vcDVwCHg+8B7AarqeJIPAw+1\n6W6qqoUHqqVxNwN8IQl01p3/WFX/KclD9LE+SONq2UCoqncvMeqKRaYt4IYl5rML2NVXddIYqapv\nAT+zSPuz9Lk+SOPIK5UlSYCBIElqDARJEmAgSJKa8boKRNJYmD3Li8w0mdxCkCQBBoIkqTEQJEmA\nxxAkaUmLHUvZftGpM95v68mb37aaJa0qtxAkSYCBIElqDARJEmAgSJIaA0GSBBgIkqTGQJAkAQaC\nJKkxECRJgIEgSWoMBEkSYCBIkhoDQZIEGAiSpMbbX0tTbHbHPcverlk6zS0ESRLgFoIkDdRiX6qz\nnHH5Uh23ECRJwAgCIcnmJN9IcijJjmH/fmnYXOY1KYa6yyjJy4DfAd4CHAYeSrK3qh4fZh3SsAxy\nmV/JrgipH8M+hnApcKiqvgWQ5E5gC2AgaFq5zGtZKw37QR97GHYgrAee7np+GLise4Ik24Bt7enJ\nJN9YYl5rge/088vzkX6mnhh998M0evNHztgPrxtmLQssu8xDX8t93359QpYR6+zfGT7TXpdkW1Xt\n7Gd+Y3eWUXsDy76JJPuratMQShpr9kPHpPdDr8v9SkxK31jnYCXZT5/L1LAPKh8BNnQ9v6C1SdPK\nZV4TY9iB8BCwMcnrk/wo8C5g75BrkIbJZV4TY6i7jKrqVJL3A/cBLwN2VdVjK5zdqmxeTyD7oWMs\n+2HAy/xKjWXfLMI6B6vvOlNVq1GIJGnCeKWyJAkwECRJzdgHwnKX/Sd5RZLPtfFfTTI7/CpXXw/9\ncF2SP0/y9fb41VHUuZqS7EpyLMmjS4xPkltaHz2S5JJh1ziOkvxWkiNdy8bVo66p26Tc2iPJk0kO\ntD7cP+p6TltsvUhyfpJ9SQ62n+f1Mq+xDoSuy/6vAi4E3p3kwgWTXQ88V1VvAD4OTN3lZz32A8Dn\nquri9vj0UIscjt3A5jOMvwrY2B7bgFuHUNOk+HjXsnHvqIs5rY9le1y8ufXhOF2HsJsXrxc7gPur\naiNwf3u+rLEOBLou+6+qvwROX/bfbQuwpw1/HrgiSYZY4zD00g9Tr6oeAI6fYZItwO3V8SBwbpJ1\nw6lOK+SyfZaWWC+6Pxf3ANf0Mq9xD4TFLvtfv9Q0VXUKOAG8ZijVDU8v/QDwy21XyeeTbFhk/LTr\ntZ9eit7flo1dve4+GJJJ+psV8EdJHm63GhlnM1V1tA1/G5jp5UXjHgjq3ReB2ar6aWAfP/zvQC8B\nSf44yaOLPLbQ2XX2N4CLgaPAR0da7OT6uaq6hM7urRuS/PyoC+pFda4t6On6grG7l9ECvVz2f3qa\nw0nOAX4ceHY45Q3Nsv1QVd3v+dPAvx5CXePmJXubiKr6xV6mS/Ip4A9XuZx+TMzfrKqOtJ/HknyB\nzu6uB0Zb1ZKeSbKuqo623abHennRuG8h9HLZ/15gaxt+B/Clmr6r7Rb2w/uAf57kZNfj+SSV5F8A\nbweeGGnFo7EXuLadbXQ5cKJrs/kla8FxlH8ALHqW1ohMxK09kqxJ8urTw8CVjFc/LtT9ubgVuLun\nV1XVWD+Aq4H/DnwT+M3WdhPw9jb8SuD3gEPA14CfGHXNI+qHe4G/Ah4Dvgz85KhrXoU+uIPOLo+/\norOv+Xo64fi+Nj50zlj5JnAA2DTqmsfhAXym9ccj7YNi3ahrWlDfi5btcXsAPwH8WXs8Nk51LrFe\nvIbO2UUHgT8Gzu9lXt66YgokeSPwn4Ffqqr5EZcjaUKN+y4jLSPJuXROt/2wYSDpbLiFMMHa9RZ3\n0zmD4JryjynpLIz7WUY6sw8BPwX8rGEg6Wy5hTChkszRufbg56vqT0dcjqQp4DGECdROI7wT+KBh\nIGlQDITJ9E/oXIr+bxZci3AyySdHXZykyeQuI0kS4BaCJKkxECRJgIEgSWoMBEkSMOYXpq1du7Zm\nZ2cXHff888+zZs2a4RY0huyHjjP1w8MPP/ydqnrtkEuSJs5YB8Ls7Cz79y/+Xdbz8/PMzc0Nt6Ax\nZD90nKkfkjw13GqkyeQuI0kS0EMgJHllkq8l+bMkjyX5l6399Um+muRQks+1L7cgySva80Nt/GzX\nvP5Za/9Gkreu1puSJPWvly2EHwC/UFU/Q+c7WTe3b6P6CPDxqnoD8BydL2Wg/XyutX+8TUeSC+l8\nG9JPAZuBf5fkZYN8M5KklVv2GEK7i+bJ9vTl7VHALwD/uLXvAX6Lzpd5b2nD0LlP/79tt2neAtxZ\nVT8A/keSQ3S+k/RPBvFGND1md9zT92t2b/bAunS2ejqo3P6Tfxh4Az/8isLvVtWpNslhYH0bXg88\nDVBVp5KcoPN1buuBB7tm2/2a7t+1DdgGMDMzw/z8/KI1nTx5cslxLyXT2A/bLzq1/EQLTGM/SMPW\nUyBU1f8BLm7fzvUF4CdXq6Cq2gnsBNi0aVMtdeaIZ9d0TGM/XLfCLYRp6wdp2Po6y6iqvkvnC9z/\nDnBuktOBcgFwpA0fATYAtPE/Djzb3b7IayRJI9bLWUavbVsGJHkV8BbgCTrB8I422VY6X+UIsLc9\np43/UjsOsRd4VzsL6fXARuBrg3ojkqSz08suo3XAnnYc4UeAu6rqD5M8DtyZ5F8Bfwrc1qa/DfhM\nO2h8nM6ZRVTVY0nuAh4HTgE3tF1RkqQx0MtZRo8Ab1yk/Vt0zhJa2P6/gX+4xLx+G/jt/suUJK02\nr1SWJAEGgiSpMRAkSYCBIElqDARJEmAgSJIaA0GSBBgIkqTGQJAkAQaCJKkxECRJgIEgSWoMBEkS\nYCBIkhoDQZIEGAiSpMZAkCQBBoIkqTEQJEmAgSBJagwESRLQQyAk2ZDky0keT/JYkg+09vOT7Ety\nsP08r7UnyS1JDiV5JMklXfPa2qY/mGTr6r0tSVK/etlCOAVsr6oLgcuBG5JcCOwA7q+qjcD97TnA\nVcDG9tgG3AqdAAFuBC4DLgVuPB0ikqTRWzYQqupoVf3XNvw94AlgPbAF2NMm2wNc04a3ALdXx4PA\nuUnWAW8F9lXV8ap6DtgHbB7ou5Ekrdg5/UycZBZ4I/BVYKaqjrZR3wZm2vB64Omulx1ubUu1L/wd\n2+hsWTAzM8P8/PyitZw8eXLJcS8l09gP2y861fdrprEfpGHrORCS/Bjw+8AHq+ovkrwwrqoqSQ2i\noKraCewE2LRpU83NzS063fz8PEuNeymZxn64bsc9fb9m9+Y1U9cP0rD1FAhJXk4nDD5bVX/Qmp9J\nsq6qjrZdQsda+xFgQ9fLL2htR4C5Be3zKy38wJETfX9wPHnz21b66yRp6vVyllGA24AnqupjXaP2\nAqfPFNoK3N3Vfm072+hy4ETbtXQfcGWS89rB5CtbmyRpDPSyhfAm4FeAA0m+3tp+A7gZuCvJ9cBT\nwDvbuHuBq4FDwPeB9wJU1fEkHwYeatPdVFXHB/IuJElnbdlAqKqvAFli9BWLTF/ADUvMaxewq58C\nJUnD4ZXKkiTAQJAkNQaCJAkwECRJjYEgSQIMBElSYyBIkgADQZLUGAiSJMBAkCQ1BoIkCTAQJEmN\ngSBJAgwESVJjIEiSAANBktQYCJIkwECQJDUGgiQJMBAkSc2ygZBkV5JjSR7tajs/yb4kB9vP81p7\nktyS5FCSR5Jc0vWarW36g0m2rs7bkSStVC9bCLuBzQvadgD3V9VG4P72HOAqYGN7bANuhU6AADcC\nlwGXAjeeDhFJ0nhYNhCq6gHg+ILmLcCeNrwHuKar/fbqeBA4N8k64K3Avqo6XlXPAft4cchIkkbo\nnBW+bqaqjrbhbwMzbXg98HTXdIdb21LtL5JkG52tC2ZmZpifn1+8gFfB9otO9VX0UvOaZCdPnpy6\n99Xv3xWmsx+kYVtpILygqipJDaKYNr+dwE6ATZs21dzc3KLTfeKzd/PRA/2V/+R7Fp/XJJufn2ep\nPppU1+24p+/X7N68Zur6QRq2lZ5l9EzbFUT7eay1HwE2dE13QWtbql2SNCZWGgh7gdNnCm0F7u5q\nv7adbXQ5cKLtWroPuDLJee1g8pWtTZI0Jpbd55LkDmAOWJvkMJ2zhW4G7kpyPfAU8M42+b3A1cAh\n4PvAewGq6niSDwMPteluqqqFB6olSSO0bCBU1buXGHXFItMWcMMS89kF7OqrOknS0HilsiQJMBAk\nSY2BIEkCDARJUmMgSJIAA0GS1BgIkiTAQJAkNQaCJAkwECRJjYEgSQIMBElSYyBIkgADQZLUGAiS\nJMBAkCQ1BoIkCTAQJEmNgSBJAgwESVJjIEiSgBEEQpLNSb6R5FCSHcP+/ZKkxQ01EJK8DPgd4Crg\nQuDdSS4cZg2SpMUNewvhUuBQVX2rqv4SuBPYMuQaJEmLOGfIv2898HTX88PAZd0TJNkGbGtPTyb5\nxhLzWgt8p59fno/0M/XE6LsfptGbP3LGfnjdMGuRJtWwA2FZVbUT2LncdEn2V9WmIZQ01uyHDvtB\nOnvD3mV0BNjQ9fyC1iZJGrFhB8JDwMYkr0/yo8C7gL1DrkGStIih7jKqqlNJ3g/cB7wM2FVVj61w\ndsvuVnqJsB867AfpLKWqRl2DJGkMeKWyJAkwECRJzdgHwnK3ukjyiiSfa+O/mmR2+FWuvh764bok\nf57k6+3xq6OoczUl2ZXkWJJHlxifJLe0PnokySXDrlGaZGMdCD3e6uJ64LmqegPwcWDqLj/r45Yf\nn6uqi9vj00Mtcjh2A5vPMP4qYGN7bANuHUJN0tQY60Cgt1tdbAH2tOHPA1ckyRBrHAZv+QFU1QPA\n8TNMsgW4vToeBM5Nsm441UmTb9wDYbFbXaxfapqqOgWcAF4zlOqGp5d+APjltqvk80k2LDJ+2vXa\nT5IWMe6BoN59EZitqp8G9vHDrSZJ6sm4B0Ivt7p4YZok5wA/Djw7lOqGZ9l+qKpnq+oH7emngZ8d\nUm3jxFujSGdh3AOhl1td7AW2tuF3AF+q6bvabtl+WLCv/O3AE0Osb1zsBa5tZxtdDpyoqqOjLkqa\nFGN3t9NuS93qIslNwP6q2gvcBnwmySE6BxzfNbqKV0eP/fDrSd4OnKLTD9eNrOBVkuQOYA5Ym+Qw\ncCPwcoCq+iRwL3A1cAj4PvDe0VQqTSZvXSFJAsZ/l5EkaUgMBEkSYCBIkhoDQZIEGAiSpMZAkCQB\nBoIkqfl/RCydfmJ8Iv4AAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x11707cf10>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"df.hist();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We want to estimate $E(Y|do(D=d)) = \\sum_Z P(Z|D) \\sum_{D'} E[Y|D', Z]P(D') = \\sum_Z P(Z|D) P(Y|do(Z=z))$, so we need each component. Y is real-valued, but Z and X are discrete (binary).\n",
"\n",
"* $P(Z|D)$ we can use a logistic regression\n",
"* $E[Y|D', Z]$ we can use linear regression \n",
"* $\\sum_{D'}E[Y|D', Z] P(D')$ we can do our usual trick, so we don't have to compute P(D'), or since it's discrete, we can just calculate it\n",
"\n",
"Let's start with $P(Z|D)$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### What do we do first? CHECK SUPPORT\n",
"\n",
"We're measuring two effects, so we should check both. We need P(Y, Z) >0, and P(Z, D) > 0.\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 203,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"Z\n",
"0 [[AxesSubplot(0.125,0.125;0.775x0.755)]]\n",
"1 [[AxesSubplot(0.125,0.125;0.775x0.755)]]\n",
"dtype: object"
]
},
"execution_count": 203,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEICAYAAABWJCMKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAE1dJREFUeJzt3X+QXeV93/H3J1LANrIlDO2OR1Ij\nOlHaUugPvMVkPE1XUcYRxEXMxHFgSC1cppqkxHUDbZGbP+gkkylMSjw249hVKwa5US0ISSuNIXEZ\nmR3GbUWM4hTxI443WAYpBMUG1Kyx48r59o97lG4Vybu6Z/deNs/7NbOz5zznOed5vler/dxzzr13\nU1VIktrzXeOegCRpPAwASWqUASBJjTIAJKlRBoAkNcoAkKRGGQCS1CgDQFqgJEeSfCPJHyd5Ncn/\nSPKTSfx/pGXJH1zp3PzDqnoz8D3AncDtwK7xTkkajgEgDaGqTlTVfuDHgW1JLhv3nKRzZQBIPVTV\nbwFHgb8/7rlI58oAkPr7A+Ct456EdK4MAKm/tcDL456EdK4MAKmHJH+PQQB8btxzkc6VASANIclb\nkrwb2Av8SlUdHvecpHMV/x6AtDBJjgATwEngT4FngF8BPlFV3x7j1KShGACS1CgvAUlSowwASWqU\nASBJjTIAJKlRK8c9ge/k4osvrg0bNgy9/9e//nUuuOCCxZvQMtBaza3VC9bcij41Hzp06KtV9Zfm\n6/e6DoANGzbwxBNPDL3/9PQ0U1NTizehZaC1mlurF6y5FX1qTvKVhfTzEpAkNcoAkKRGGQCS1CgD\nQJIaZQBIUqMMAElqlAEgSY0yACSpUQaAJDXqdf1O4L4OHzvBTTseGvm4R+78kZGPKUnnyjMASWrU\nvAGQ5N4kx5M8NaftF5P8bpInk/yXJGvmbPtQkpkkX0zyw3Pat3RtM0l2LH4pkqRzsZAzgPuALae1\nPQJcVlV/C/g94EMASS4Frgf+ZrfPLydZkWQF8DHgauBS4IauryRpTOYNgKp6DHj5tLb/VlUnu9WD\nwLpueSuwt6r+pKq+DMwAV3ZfM1X1XFV9C9jb9ZUkjcli3AT+x8D93fJaBoFwytGuDeCF09rfcaaD\nJdkObAeYmJhgenp66IlNvBFuu/zk/B0XWZ859zU7OzvW8UettXrBmlsxipp7BUCSnwVOAnsWZzpQ\nVTuBnQCTk5PV5zPA79mzj7sPj/6FTkdunBr5mKe09rnprdUL1tyKUdQ89G/HJDcB7wY2V1V1zceA\n9XO6reva+A7tkqQxGOploEm2AP8KuLaqXpuzaT9wfZLzk1wCbAR+C/g8sDHJJUnOY3CjeH+/qUuS\n+pj3DCDJp4Ap4OIkR4E7GLzq53zgkSQAB6vqJ6vq6SQPAM8wuDR0S1V9uzvOTwOfAVYA91bV00tQ\njyRpgeYNgKq64QzNu75D/18AfuEM7Q8DD5/T7CRJS8Z3AktSowwASWqUASBJjTIAJKlRBoAkNcoA\nkKRGGQCS1CgDQJIaZQBIUqMMAElqlAEgSY0yACSpUQaAJDXKAJCkRhkAktQoA0CSGmUASFKjDABJ\napQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhplAEhSo+YNgCT3Jjme5Kk5bW9N8kiSL3XfL+za\nk+SjSWaSPJnkijn7bOv6fynJtqUpR5K0UAs5A7gP2HJa2w7gQFVtBA506wBXAxu7r+3Ax2EQGMAd\nwDuAK4E7ToWGJGk85g2AqnoMePm05q3A7m55N3DdnPZP1sBBYE2StwE/DDxSVS9X1SvAI/z5UJEk\njdDKIfebqKoXu+U/BCa65bXAC3P6He3aztb+5yTZzuDsgYmJCaanp4ecIky8EW67/OTQ+w+rz5z7\nmp2dHev4o9ZavWDNrRhFzcMGwJ+pqkpSizGZ7ng7gZ0Ak5OTNTU1NfSx7tmzj7sP9y7xnB25cWrk\nY54yPT1Nn8dsuWmtXrDmVoyi5mFfBfRSd2mH7vvxrv0YsH5Ov3Vd29naJUljMmwA7AdOvZJnG7Bv\nTvv7ulcDXQWc6C4VfQZ4V5ILu5u/7+raJEljMu/1kSSfAqaAi5McZfBqnjuBB5LcDHwFeG/X/WHg\nGmAGeA14P0BVvZzk54HPd/1+rqpOv7EsSRqheQOgqm44y6bNZ+hbwC1nOc69wL3nNDtJ0pLxncCS\n1CgDQJIaZQBIUqMMAElqlAEgSY0yACSpUQaAJDXKAJCkRhkAktQoA0CSGmUASFKjDABJapQBIEmN\nMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhplAEhSowwASWqUASBJjTIAJKlRBoAkNcoAkKRG9QqAJD+T\n5OkkTyX5VJI3JLkkyeNJZpLcn+S8ru/53fpMt33DYhQgSRrO0AGQZC3wz4DJqroMWAFcD9wFfLiq\nvhd4Bbi52+Vm4JWu/cNdP0nSmPS9BLQSeGOSlcCbgBeBHwQe7LbvBq7rlrd263TbNydJz/ElSUMa\nOgCq6hjw74DnGfziPwEcAl6tqpNdt6PA2m55LfBCt+/Jrv9Fw44vSepn5bA7JrmQwbP6S4BXgV8F\ntvSdUJLtwHaAiYkJpqenhz7WxBvhtstPzt9xkfWZc1+zs7NjHX/UWqsXrLkVo6h56AAAfgj4clX9\nEUCSXwfeCaxJsrJ7lr8OONb1PwasB452l4xWA187/aBVtRPYCTA5OVlTU1NDT/CePfu4+3CfEodz\n5MapkY95yvT0NH0es+WmtXrBmlsxipr73AN4HrgqyZu6a/mbgWeAR4H3dH22Afu65f3dOt32z1ZV\n9RhfktRDn3sAjzO4mfvbwOHuWDuB24Fbk8wwuMa/q9tlF3BR134rsKPHvCVJPfW6PlJVdwB3nNb8\nHHDlGfp+E/ixPuNJkhaP7wSWpEYZAJLUKANAkhplAEhSowwASWqUASBJjTIAJKlRBoAkNcoAkKRG\nGQCS1CgDQJIaZQBIUqMMAElq1Oj/WookLSMbdjw0lnHv23LBko/hGYAkNcoAkKRGGQCS1CgDQJIa\nZQBIUqMMAElqlAEgSY0yACSpUQaAJDXKAJCkRhkAktSoXgGQZE2SB5P8bpJnk3x/krcmeSTJl7rv\nF3Z9k+SjSWaSPJnkisUpQZI0jL5nAB8BfrOq/jrwt4FngR3AgaraCBzo1gGuBjZ2X9uBj/ccW5LU\nw9ABkGQ18APALoCq+lZVvQpsBXZ33XYD13XLW4FP1sBBYE2Stw09c0lSL6mq4XZM/g6wE3iGwbP/\nQ8AHgWNVtabrE+CVqlqT5NPAnVX1uW7bAeD2qnritONuZ3CGwMTExNv37t071PwAjr98gpe+MfTu\nQ7t87erRD9qZnZ1l1apVYxt/1FqrF6x51A4fOzGWcS9ZvWLomjdt2nSoqibn69fn7wGsBK4APlBV\njyf5CP/vcg8AVVVJzilhqmong2BhcnKypqamhp7gPXv2cffh0f/JgyM3To18zFOmp6fp85gtN63V\nC9Y8ajeN8e8BLHXNfe4BHAWOVtXj3fqDDALhpVOXdrrvx7vtx4D1c/Zf17VJksZg6ACoqj8EXkjy\n17qmzQwuB+0HtnVt24B93fJ+4H3dq4GuAk5U1YvDji9J6qfv9ZEPAHuSnAc8B7yfQag8kORm4CvA\ne7u+DwPXADPAa11fSdKY9AqAqvod4Ew3GjafoW8Bt/QZT5K0eHwnsCQ1ygCQpEYZAJLUKANAkhpl\nAEhSowwASWqUASBJjTIAJKlRBoAkNcoAkKRGGQCS1CgDQJIaZQBIUqMMAElqlAEgSY0yACSpUQaA\nJDXKAJCkRhkAktQoA0CSGmUASFKjDABJapQBIEmNMgAkqVEGgCQ1qncAJFmR5AtJPt2tX5Lk8SQz\nSe5Pcl7Xfn63PtNt39B3bEnS8BbjDOCDwLNz1u8CPlxV3wu8Atzctd8MvNK1f7jrJ0kak14BkGQd\n8CPAf+zWA/wg8GDXZTdwXbe8tVun27656y9JGoNU1fA7Jw8C/xZ4M/AvgJuAg92zfJKsB36jqi5L\n8hSwpaqOdtt+H3hHVX31tGNuB7YDTExMvH3v3r1Dz+/4yyd46RtD7z60y9euHv2gndnZWVatWjW2\n8UettXrBmkft8LETYxn3ktUrhq5506ZNh6pqcr5+K4c6OpDk3cDxqjqUZGrY45yuqnYCOwEmJydr\namr4Q9+zZx93Hx66xKEduXFq5GOeMj09TZ/HbLlprV6w5lG7acdDYxn3vi0XLHnNfX47vhO4Nsk1\nwBuAtwAfAdYkWVlVJ4F1wLGu/zFgPXA0yUpgNfC1HuNLknoY+h5AVX2oqtZV1QbgeuCzVXUj8Cjw\nnq7bNmBft7y/W6fb/tnqc/1JktTLUrwP4Hbg1iQzwEXArq59F3BR134rsGMJxpYkLdCiXCCvqmlg\nult+DrjyDH2+CfzYYownSerPdwJLUqMMAElqlAEgSY0yACSpUQaAJDXKAJCkRhkAktQoA0CSGmUA\nSFKjDABJapQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhplAEhSowwASWqUASBJjTIAJKlRBoAk\nNcoAkKRGGQCS1CgDQJIaNXQAJFmf5NEkzyR5OskHu/a3JnkkyZe67xd27Uny0SQzSZ5McsViFSFJ\nOnd9zgBOArdV1aXAVcAtSS4FdgAHqmojcKBbB7ga2Nh9bQc+3mNsSVJPQwdAVb1YVb/dLf8x8Cyw\nFtgK7O667Qau65a3Ap+sgYPAmiRvG3rmkqReUlX9D5JsAB4DLgOer6o1XXuAV6pqTZJPA3dW1ee6\nbQeA26vqidOOtZ3BGQITExNv37t379DzOv7yCV76xtC7D+3ytatHP2hndnaWVatWjW38UWutXrDm\nUTt87MRYxr1k9Yqha960adOhqpqcr9/KoY4+R5JVwK8B/7yq/vfgd/5AVVWSc0qYqtoJ7ASYnJys\nqampoed2z5593H24d4nn7MiNUyMf85Tp6Wn6PGbLTWv1gjWP2k07HhrLuPdtuWDJa+71KqAk383g\nl/+eqvr1rvmlU5d2uu/Hu/ZjwPo5u6/r2iRJY9DnVUABdgHPVtUvzdm0H9jWLW8D9s1pf1/3aqCr\ngBNV9eKw40uS+ulzfeSdwD8CDif5na7tXwN3Ag8kuRn4CvDebtvDwDXADPAa8P4eY0uSeho6ALqb\nuTnL5s1n6F/ALcOOJ0laXL4TWJIaZQBIUqMMAElqlAEgSY0yACSpUQaAJDXKAJCkRhkAktQoA0CS\nGmUASFKjDABJapQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhplAEhSowwASWqUASBJjTIAJKlR\nBoAkNcoAkKRGGQCS1CgDQJIaNfIASLIlyReTzCTZMerxJUkDIw2AJCuAjwFXA5cCNyS5dJRzkCQN\njPoM4Epgpqqeq6pvAXuBrSOegyQJWDni8dYCL8xZPwq8Y26HJNuB7d3qbJIv9hjvYuCrPfYfSu4a\n9Yj/n7HUPEat1QvW3IRNd/Wq+XsW0mnUATCvqtoJ7FyMYyV5oqomF+NYy0VrNbdWL1hzK0ZR86gv\nAR0D1s9ZX9e1SZJGbNQB8HlgY5JLkpwHXA/sH/EcJEmM+BJQVZ1M8tPAZ4AVwL1V9fQSDrkol5KW\nmdZqbq1esOZWLHnNqaqlHkOS9DrkO4ElqVEGgCQ1atkHwHwfLZHk/CT3d9sfT7Jh9LNcXAuo+dYk\nzyR5MsmBJAt6TfDr2UI/QiTJjyapJMv+JYMLqTnJe7t/66eT/OdRz3GxLeBn+68keTTJF7qf72vG\nMc/FkuTeJMeTPHWW7Uny0e7xeDLJFYs6gapatl8MbiT/PvBXgfOA/wVcelqffwp8olu+Hrh/3PMe\nQc2bgDd1yz/VQs1dvzcDjwEHgclxz3sE/84bgS8AF3brf3nc8x5BzTuBn+qWLwWOjHvePWv+AeAK\n4KmzbL8G+A0gwFXA44s5/nI/A1jIR0tsBXZ3yw8Cm5NkhHNcbPPWXFWPVtVr3epBBu+3WM4W+hEi\nPw/cBXxzlJNbIgup+Z8AH6uqVwCq6viI57jYFlJzAW/pllcDfzDC+S26qnoMePk7dNkKfLIGDgJr\nkrxtscZf7gFwpo+WWHu2PlV1EjgBXDSS2S2NhdQ8180MnkEsZ/PW3J0ar6+qh0Y5sSW0kH/n7wO+\nL8l/T3IwyZaRzW5pLKTmfwP8RJKjwMPAB0YztbE51//v5+R191EQWjxJfgKYBP7BuOeylJJ8F/BL\nwE1jnsqorWRwGWiKwVneY0kur6pXxzqrpXUDcF9V3Z3k+4H/lOSyqvrTcU9sOVruZwAL+WiJP+uT\nZCWD08avjWR2S2NBH6eR5IeAnwWurao/GdHclsp8Nb8ZuAyYTnKEwbXS/cv8RvBC/p2PAvur6v9U\n1ZeB32MQCMvVQmq+GXgAoKr+J/AGBh8U9xfVkn58znIPgIV8tMR+YFu3/B7gs9XdXVmm5q05yd8F\n/j2DX/7L/bowzFNzVZ2oqourakNVbWBw3+PaqnpiPNNdFAv52f6vDJ79k+RiBpeEnhvlJBfZQmp+\nHtgMkORvMAiAPxrpLEdrP/C+7tVAVwEnqurFxTr4sr4EVGf5aIkkPwc8UVX7gV0MThNnGNxsuX58\nM+5vgTX/IrAK+NXufvfzVXXt2Cbd0wJr/gtlgTV/BnhXkmeAbwP/sqqW7dntAmu+DfgPSX6GwQ3h\nm5bzE7okn2IQ4hd39zXuAL4boKo+weA+xzXADPAa8P5FHX8ZP3aSpB6W+yUgSdKQDABJapQBIEmN\nMgAkqVEGgCQ1ygCQpEYZAJLUqP8LT1/czhET+R0AAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x115c1bed0>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEICAYAAAC55kg0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAFRNJREFUeJzt3X+wXGd93/H3J1JwUwlsiMkdRzaR\nmJEz9Y/WxbeOO23o1TgF4SYY2g6xh8QWMAgCdJrG00Y0mYGB8Qw0cZixoSai9sgEx8KFgjWxXeq4\n3HFoI0AG17IdHGRbBAlXKsjIueC6kf3tH3sEi5GsvburXV8979fMzj377HPO83yvrvTRec7Zvakq\nJElt+olpT0CSND2GgCQ1zBCQpIYZApLUMENAkhpmCEhSwwwBSWqYISANKMnuJE8m+esk303yP5O8\nPYl/j7Rk+cMrLc6vVNULgZ8DPgD8NnD9dKckDc8QkIZQVQerahvwq8AVSc6Z9pykYRgC0giq6kvA\nHuAXpz0XaRiGgDS6bwEvmfYkpGEYAtLoVgEHpj0JaRiGgDSCJP+AXgh8YdpzkYZhCEhDSPKiJL8M\nbAU+UVU7pz0naRjx9wlIg0myG5gBDgHPAA8CnwA+WlVPT3Fq0tAMAUlqmMtBktQwQ0CSGmYISFLD\nDAFJatjyaU/gWE499dRavXr1UPt+73vfY8WKFeOd0POcNZ/4WqsXrHmx7rnnnm9X1UsH6fu8D4HV\nq1ezY8eOofadn59nbm5uvBN6nrPmE19r9YI1L1aSbwza1+UgSWqYISBJDTMEJKlhhoAkNcwQkKSG\nGQKS1DBDQJIaZghIUsMMAUlq2PP+HcOSNE2rN902lXG3rJ/Mx2R4JiBJDTMEJKlhhoAkNcwQkKSG\nHTMEktyQZH+S+/vaPpnk3u6xO8m9XfvqJE/2vfbRvn3OT7Izya4k1yTJ8SlJkjSoQe4O2gJ8GPj4\n4Yaq+tXD20muBg729X+4qs47wnGuA94KfBG4HVgP3LH4KUuSxuWYZwJVdTdw4Eivdf+bfwNw83Md\nI8lpwIuqantVFb1Aed3ipytJGqdR3yfwi8C+qvp6X9uaJF8FngB+t6r+DFgF7Onrs6drO6IkG4GN\nADMzM8zPzw81uYWFhaH3Xaqs+cTXWr0w3ZqvPPfQVMadVM2jhsBl/OhZwGPAy6rqO0nOBz6b5OzF\nHrSqNgObAWZnZ2vYX7Hmr6RrQ2s1t1YvTLfmDVN8s9gkah46BJIsB/45cP7htqp6Cniq274nycPA\nmcBe4PS+3U/v2iRJUzTKLaK/BHytqn6wzJPkpUmWddsvB9YCj1TVY8ATSS7sriNcDtw6wtiSpDEY\n5BbRm4E/B34+yZ4kb+leupQfvyD8SuC+7pbRTwFvr6rDF5XfAfwnYBfwMN4ZJElTd8zloKq67Cjt\nG47Q9mng00fpvwM4Z5HzkyQdR75jWJIaZghIUsMMAUlqmCEgSQ0zBCSpYYaAJDXMEJCkhhkCktQw\nQ0CSGmYISFLDDAFJapghIEkNMwQkqWGGgCQ1zBCQpIYZApLUMENAkhpmCEhSwwwBSWrYIL9o/oYk\n+5Pc39f23iR7k9zbPS7ue+3dSXYleSjJq/va13dtu5JsGn8pkqTFGuRMYAuw/gjtH6qq87rH7QBJ\nzgIuBc7u9vmPSZYlWQZ8BHgNcBZwWddXkjRFy4/VoaruTrJ6wONdAmytqqeAR5PsAi7oXttVVY8A\nJNna9X1w0TOWJI3NKNcE3pXkvm656MVd2yrgm3199nRtR2uXJE3RMc8EjuI64P1AdV+vBt48rkkl\n2QhsBJiZmWF+fn6o4ywsLAy971JlzSe+1uqF6dZ85bmHpjLupGoeKgSqat/h7SQfA/6ke7oXOKOv\n6+ldG8/RfqTjbwY2A8zOztbc3Nww02R+fp5h912qrPnE11q9MN2aN2y6bSrjblm/YiI1D7UclOS0\nvqevBw7fObQNuDTJSUnWAGuBLwFfBtYmWZPkBfQuHm8bftqSpHE45plAkpuBOeDUJHuA9wBzSc6j\ntxy0G3gbQFU9kOQWehd8DwHvrKqnu+O8C/gcsAy4oaoeGHs1kqRFGeTuoMuO0Hz9c/S/CrjqCO23\nA7cvanaSpOPKdwxLUsMMAUlqmCEgSQ0zBCSpYYaAJDXMEJCkhhkCktQwQ0CSGmYISFLDDAFJapgh\nIEkNMwQkqWGGgCQ1zBCQpIYZApLUMENAkhpmCEhSwwwBSWqYISBJDTMEJKlhxwyBJDck2Z/k/r62\n30vytST3JflMklO69tVJnkxyb/f4aN8+5yfZmWRXkmuS5PiUJEka1CBnAluA9c9quxM4p6r+LvCX\nwLv7Xnu4qs7rHm/va78OeCuwtns8+5iSpAk7ZghU1d3AgWe1/beqOtQ93Q6c/lzHSHIa8KKq2l5V\nBXwceN1wU5YkjcvyMRzjzcAn+56vSfJV4Angd6vqz4BVwJ6+Pnu6tiNKshHYCDAzM8P8/PxQE9t/\n4CDX3nTrUPuO4txVJ098zMMWFhaG/n4tVa3V3Fq9MN2arzz30LE7HQeTqnmkEEjyO8Ah4Kau6THg\nZVX1nSTnA59NcvZij1tVm4HNALOzszU3NzfU/K696Vau3jmOnFuc3W+cm/iYh83PzzPs92upaq3m\n1uqF6da8YdNtUxl3y/oVE6l56H8hk2wAfhm4qFvioaqeAp7qtu9J8jBwJrCXH10yOr1rkyRN0VC3\niCZZD/w74LVV9f2+9pcmWdZtv5zeBeBHquox4IkkF3Z3BV0OTH6dRpL0I455JpDkZmAOODXJHuA9\n9O4GOgm4s7vTc3t3J9Argfcl+RvgGeDtVXX4ovI76N1p9FPAHd1DkjRFxwyBqrrsCM3XH6Xvp4FP\nH+W1HcA5i5qdJOm48h3DktQwQ0CSGmYISFLDDAFJapghIEkNMwQkqWGGgCQ1zBCQpIYZApLUMENA\nkhpmCEhSwwwBSWqYISBJDTMEJKlhhoAkNcwQkKSGGQKS1DBDQJIaZghIUsMGCoEkNyTZn+T+vraX\nJLkzyde7ry/u2pPkmiS7ktyX5BV9+1zR9f96kivGX44kaTEGPRPYAqx/Vtsm4K6qWgvc1T0HeA2w\ntntsBK6DXmgA7wF+AbgAeM/h4JAkTcdAIVBVdwMHntV8CXBjt30j8Lq+9o9Xz3bglCSnAa8G7qyq\nA1X1OHAnPx4skqQJWj7CvjNV9Vi3/b+BmW57FfDNvn57urajtf+YJBvpnUUwMzPD/Pz8cBP8Kbjy\n3END7TuKYec7DgsLC1Mdfxpaq7m1emG6NU/j3xCYXM2jhMAPVFUlqXEcqzveZmAzwOzsbM3NzQ11\nnGtvupWrd46lxEXZ/ca5iY952Pz8PMN+v5aq1mpurV6Ybs0bNt02lXG3rF8xkZpHuTtoX7fMQ/d1\nf9e+Fzijr9/pXdvR2iVJUzJKCGwDDt/hcwVwa1/75d1dQhcCB7tlo88Br0ry4u6C8Ku6NknSlAy0\nVpLkZmAOODXJHnp3+XwAuCXJW4BvAG/out8OXAzsAr4PvAmgqg4keT/w5a7f+6rq2RebJUkTNFAI\nVNVlR3npoiP0LeCdRznODcANA89OknRc+Y5hSWqYISBJDTMEJKlhhoAkNcwQkKSGGQKS1DBDQJIa\nZghIUsMMAUlqmCEgSQ0zBCSpYYaAJDXMEJCkhhkCktQwQ0CSGmYISFLDDAFJapghIEkNMwQkqWGG\ngCQ1bOgQSPLzSe7tezyR5DeTvDfJ3r72i/v2eXeSXUkeSvLq8ZQgSRrW8mF3rKqHgPMAkiwD9gKf\nAd4EfKiqfr+/f5KzgEuBs4GfBf40yZlV9fSwc5AkjWZcy0EXAQ9X1Teeo88lwNaqeqqqHgV2AReM\naXxJ0hBSVaMfJLkB+EpVfTjJe4ENwBPADuDKqno8yYeB7VX1iW6f64E7qupTRzjeRmAjwMzMzPlb\nt24dal77Dxxk35ND7TqSc1edPPlBOwsLC6xcuXJq409DazW3Vi9Mt+adew9OZdw1Jy8buuZ169bd\nU1Wzg/QdOQSSvAD4FnB2Ve1LMgN8Gyjg/cBpVfXmxYRAv9nZ2dqxY8dQc7v2plu5eufQK15D2/2B\nfzbxMQ+bn59nbm5uauNPQ2s1t1YvTLfm1Ztum8q4W9avGLrmJAOHwDiWg15D7yxgH0BV7auqp6vq\nGeBj/HDJZy9wRt9+p3dtkqQpGUcIXAbcfPhJktP6Xns9cH+3vQ24NMlJSdYAa4EvjWF8SdKQRlor\nSbIC+KfA2/qa/0OS8+gtB+0+/FpVPZDkFuBB4BDwTu8MkqTpGikEqup7wE8/q+3Xn6P/VcBVo4wp\nSRof3zEsSQ0zBCSpYYaAJDXMEJCkhhkCktQwQ0CSGmYISFLDDAFJapghIEkNMwQkqWGGgCQ1zBCQ\npIYZApLUMENAkhpmCEhSwwwBSWqYISBJDTMEJKlhhoAkNWzkEEiyO8nOJPcm2dG1vSTJnUm+3n19\ncdeeJNck2ZXkviSvGHV8SdLwxnUmsK6qzquq2e75JuCuqloL3NU9B3gNsLZ7bASuG9P4kqQhHK/l\noEuAG7vtG4HX9bV/vHq2A6ckOe04zUGSdAypqtEOkDwKPA4U8IdVtTnJd6vqlO71AI9X1SlJ/gT4\nQFV9oXvtLuC3q2rHs465kd6ZAjMzM+dv3bp1qLntP3CQfU8OW9nwzl118uQH7SwsLLBy5cqpjT8N\nrdXcWr0w3Zp37j04lXHXnLxs6JrXrVt3T9/KzHNaPtQIP+ofV9XeJD8D3Jnka/0vVlUlWVTSVNVm\nYDPA7Oxszc3NDTWxa2+6lat3jqPExdn9xrmJj3nY/Pw8w36/lqrWam6tXphuzRs23TaVcbesXzGR\nmkdeDqqqvd3X/cBngAuAfYeXebqv+7vue4Ez+nY/vWuTJE3BSCGQZEWSFx7eBl4F3A9sA67oul0B\n3NptbwMu7+4SuhA4WFWPjTIHSdLwRl0rmQE+01v2Zznwx1X1X5N8GbglyVuAbwBv6PrfDlwM7AK+\nD7xpxPElSSMYKQSq6hHg7x2h/TvARUdoL+Cdo4wpSRof3zEsSQ0zBCSpYYaAJDXMEJCkhhkCktQw\nQ0CSGmYISFLDDAFJapghIEkNMwQkqWGGgCQ1zBCQpIYZApLUMENAkhpmCEhSwwwBSWqYISBJDTME\nJKlhhoAkNWzoEEhyRpLPJ3kwyQNJ/nXX/t4ke5Pc2z0u7tvn3Ul2JXkoyavHUYAkaXij/KL5Q8CV\nVfWVJC8E7klyZ/fah6rq9/s7JzkLuBQ4G/hZ4E+TnFlVT48wB0nSCIY+E6iqx6rqK932XwN/Aax6\njl0uAbZW1VNV9SiwC7hg2PElSaNLVY1+kGQ1cDdwDvBbwAbgCWAHvbOFx5N8GNheVZ/o9rkeuKOq\nPnWE420ENgLMzMycv3Xr1qHmtf/AQfY9OdSuIzl31cmTH7SzsLDAypUrpzb+NLRWc2v1wnRr3rn3\n4FTGXXPysqFrXrdu3T1VNTtI31GWgwBIshL4NPCbVfVEkuuA9wPVfb0aePNijllVm4HNALOzszU3\nNzfU3K696Vau3jlyiYu2+41zEx/zsPn5eYb9fi1VrdXcWr0w3Zo3bLptKuNuWb9iIjWPdHdQkp+k\nFwA3VdV/AaiqfVX1dFU9A3yMHy757AXO6Nv99K5NkjQlo9wdFOB64C+q6g/62k/r6/Z64P5uextw\naZKTkqwB1gJfGnZ8SdLoRlkr+UfArwM7k9zbtf174LIk59FbDtoNvA2gqh5IcgvwIL07i97pnUGS\nNF1Dh0BVfQHIEV66/Tn2uQq4atgxJUnj5TuGJalhhoAkNcwQkKSGGQKS1DBDQJIaZghIUsMMAUlq\nmCEgSQ0zBCSpYYaAJDXMEJCkhhkCktQwQ0CSGmYISFLDDAFJapghIEkNMwQkqWGGgCQ1zBCQpIYZ\nApLUsImHQJL1SR5KsivJpkmPL0n6oYmGQJJlwEeA1wBnAZclOWuSc5Ak/dCkzwQuAHZV1SNV9f+A\nrcAlE56DJKmzfMLjrQK+2fd8D/ALz+6UZCOwsXu6kOShIcc7Ffj2kPsOLR+c9Ig/Yio1T1lrNbdW\nLzRY87oPjlTzzw3acdIhMJCq2gxsHvU4SXZU1ewYprRkWPOJr7V6wZqPp0kvB+0Fzuh7fnrXJkma\ngkmHwJeBtUnWJHkBcCmwbcJzkCR1JrocVFWHkrwL+BywDLihqh44jkOOvKS0BFnzia+1esGaj5tU\n1STGkSQ9D/mOYUlqmCEgSQ07IULgWB9FkeSkJJ/sXv9iktWTn+X4DFDvbyV5MMl9Se5KMvA9w89X\ng37cSJJ/kaSSLPnbCQepOckbuj/rB5L88aTnOG4D/Gy/LMnnk3y1+/m+eBrzHJckNyTZn+T+o7ye\nJNd034/7krxi7JOoqiX9oHeB+WHg5cALgP8FnPWsPu8APtptXwp8ctrzPs71rgP+drf9G0u53kFr\n7vq9ELgb2A7MTnveE/hzXgt8FXhx9/xnpj3vCdS8GfiNbvssYPe05z1iza8EXgHcf5TXLwbuAAJc\nCHxx3HM4Ec4EBvkoikuAG7vtTwEXJckE5zhOx6y3qj5fVd/vnm6n936MpWzQjxt5P/BB4P9OcnLH\nySA1vxX4SFU9DlBV+yc8x3EbpOYCXtRtnwx8a4LzG7uquhs48BxdLgE+Xj3bgVOSnDbOOZwIIXCk\nj6JYdbQ+VXUIOAj89ERmN36D1NvvLfT+J7GUHbPm7jT5jKq6bZITO44G+XM+Ezgzyf9Isj3J+onN\n7vgYpOb3Ar+WZA9wO/CvJjO1qVns3/dFe15+bITGI8mvAbPAP5n2XI6nJD8B/AGwYcpTmbTl9JaE\n5uid7d2d5Nyq+u5UZ3V8XQZsqaqrk/xD4I+SnFNVz0x7YkvViXAmMMhHUfygT5Ll9E4jvzOR2Y3f\nQB+9keSXgN8BXltVT01obsfLsWp+IXAOMJ9kN721021L/OLwIH/Oe4BtVfU3VfUo8Jf0QmGpGqTm\ntwC3AFTVnwN/i96Hy52ojvtH7ZwIITDIR1FsA67otv8l8N+ru+qyBB2z3iR/H/hDegGw1NeJ4Rg1\nV9XBqjq1qlZX1Wp610FeW1U7pjPdsRjk5/qz9M4CSHIqveWhRyY5yTEbpOa/Ai4CSPJ36IXA/5no\nLCdrG3B5d5fQhcDBqnpsnAMs+eWgOspHUSR5H7CjqrYB19M7bdxF7yLMpdOb8WgGrPf3gJXAf+6u\nf/9VVb12apMe0YA1n1AGrPlzwKuSPAg8DfzbqlqqZ7iD1nwl8LEk/4beReINS/g/dCS5mV6Qn9pd\n53gP8JMAVfVRetc9LgZ2Ad8H3jT2OSzh758kaUQnwnKQJGlIhoAkNcwQkKSGGQKS1DBDQJIaZghI\nUsMMAUlq2P8HBVraZst7wowAAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x11566c150>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"df.groupby('Z').hist('D')"
]
},
{
"cell_type": "code",
"execution_count": 205,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"Z\n",
"0 [[AxesSubplot(0.125,0.125;0.775x0.755)]]\n",
"1 [[AxesSubplot(0.125,0.125;0.775x0.755)]]\n",
"dtype: object"
]
},
"execution_count": 205,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAFA5JREFUeJzt3X+M5Hd93/HnqwYcy5vaISZby7a6\nRjKusC+43MZNhaC7dZIeGOEQpdRW6vqAcpASlKpXpcakxSqisggOpaIlMrFjENQLwvxwjZPiul3c\nSDFhjzg+G0Niu0fx1Tlj49gsWNCz3/1j5+hy7O7M7czczPeT50Na3cxnvvv9vvZ25rXf+cx3vpOq\nQpLUrr826QCSpPGy6CWpcRa9JDXOopekxln0ktQ4i16SGmfRS1LjLHrpKEk+muT3jhr7e0keT3L6\npHJJ2xXfMCX9sCQ/CdwHXF5Vtyf5MeAe4N9V1Y0TDSdtg0UvbSDJPwTeA5wP/CZwQVW9crKppO2x\n6KVNJLkZeC7wMtaK/hsTjiRty3MmHUCaYv8MeBB4hyWvLvPFWGkTVXUIeIy1+Xqpsyx6SWqcRS9J\njbPoJalxHnUjSY1zj16SGmfRS1LjLHpJapxFL0mNm4p3xp522mk1Nzc36Rg/4jvf+Q4nn3zypGMM\npCtZu5ITupO1KzmhO1m7knPfvn2PVdUL+i03FUU/NzfHysrKpGP8iOXlZRYWFiYdYyBdydqVnNCd\nrF3JCd3J2pWcSb4+yHJO3UhS4yx6SWqcRS9JjbPoJalxFr0kNc6il6TGWfSS1DiLXpIaZ9FLUuOm\n4p2x0jSbu/JzAy134JqLx5xE2h736CWpce7RqymD7n2De+D6q6PvHn2SG5I8muTedWMfT3J37+tA\nkrt743NJnl532++MM7wkqb9B9uhvBD4AfOTIQFX9oyOXk1wLPLlu+Qer6oJRBZQkDadv0VfVnUnm\nNrotSYDXAX9/tLEkSaMy7IuxLwcOVdWfrxs7O8mfJPlCkpcPuX5J0pBSVf0XWtujv7Wqzj9q/IPA\nA1V1be/6icBMVT2eZCfwGeC8qnpqg3XuAfYAzM7O7lxaWhryRxm91dVVZmZmJh1jIF3JOu6c+w8+\n2X+hnh1nnLLl7UeyDrrOfusbl6787qE7WbuSc3FxcV9VzfdbbttFn+Q5wEFgZ1U9vMn3LQP/sqq2\n/Pio+fn58hOmhtOVrOPOOcqjbo5knfbj6Lvyu4fuZO1KziQDFf0wUzc/B3x1fckneUGSE3qXXwic\nAzw0xDYkSUMa5PDKm4A/As5N8nCSN/ZuuhS46ajFXwHc0zvc8pPAW6rqW6MMLEk6NoMcdXPZJuO7\nNxi7Gbh5+FjS+PWbktm74zC7j2EqSJpWngJBkhpn0UtS4yx6SWqcRS9JjbPoJalxFr0kNc6il6TG\nWfSS1Dg/YUqdcCznsJH0w9yjl6TGWfSS1DiLXpIaZ9FLUuMseklqnEUvSY2z6CWpcRa9JDXOopek\nxln0ktQ4i16SGte36JPckOTRJPeuG7s6ycEkd/e+XrXutrcneSDJ15L8g3EFlyQNZpA9+huBXRuM\nv6+qLuh93QaQ5MXApcB5ve/5T0lOGFVYSdKx63v2yqq6M8ncgOu7BFiqqu8B/yvJA8CFwB9tO6HU\nEYOeYfPANRePOYn0w1JV/RdaK/pbq+r83vWrgd3AU8AKsLeqnkjyAeCuqvpob7nrgd+vqk9usM49\nwB6A2dnZnUtLSyP4cUZrdXWVmZmZSccYSFeybjfn/oNPjiHN1mZPgkNPj369O844ZaTr68rvHrqT\ntSs5FxcX91XVfL/ltns++g8C7wKq9++1wBuOZQVVdR1wHcD8/HwtLCxsM8r4LC8vM425NtKVrNvN\nuXsC56Pfu+Mw1+4f/Uc2HPiVhZGuryu/e+hO1q7kHNS2jrqpqkNV9UxVPQt8iLXpGYCDwFnrFj2z\nNyZJmpBtFX2S09ddfS1w5IicW4BLk5yY5GzgHOCPh4soSRpG3+elSW4CFoDTkjwMvBNYSHIBa1M3\nB4A3A1TVfUk+AXwFOAy8taqeGU90SdIgBjnq5rINhq/fYvl3A+8eJpQkaXR8Z6wkNW70hxRI2pLH\n2+t4c49ekhpn0UtS4yx6SWqcRS9JjbPoJalxHnWjiRr0CBRJ2+cevSQ1zqKXpMZZ9JLUOItekhpn\n0UtS4yx6SWqcRS9JjbPoJalxvmFKmlKDvpnsxl0njzmJus49eklqnEUvSY2z6CWpcX2LPskNSR5N\ncu+6sd9K8tUk9yT5dJJTe+NzSZ5Ocnfv63fGGV6S1N8ge/Q3AruOGrsdOL+qfhr4M+Dt6257sKou\n6H29ZTQxJUnb1bfoq+pO4FtHjX2+qg73rt4FnDmGbJKkEUhV9V8omQNurarzN7jtvwAfr6qP9pa7\nj7W9/KeA36yq/7nJOvcAewBmZ2d3Li0tbe8nGKPV1VVmZmYmHWMgXcl6dM79B5+cYJqtzZ4Eh56e\ndIr+zj7lhE787qG799Nptbi4uK+q5vstN1TRJ3kHMA/8UlVVkhOBmap6PMlO4DPAeVX11Fbrn5+f\nr5WVlb45jrfl5WUWFhYmHWMgXcl6dM5p/uCRvTsOc+3+6X+ryY27Tu7E7x66ez+dVkkGKvptH3WT\nZDfwauBXqvfXoqq+V1WP9y7vAx4EXrTdbUiShretok+yC/gN4DVV9d114y9IckLv8guBc4CHRhFU\nkrQ9fZ+XJrkJWABOS/Iw8E7WjrI5Ebg9CcBdvSNsXgH82yT/F3gWeEtVfWvDFatpm03J7N1xmN1T\nPF0jtahv0VfVZRsMX7/JsjcDNw8bSpI0Or4zVpIaZ9FLUuMseklqnEUvSY2z6CWpcRa9JDXOopek\nxln0ktQ4i16SGmfRS1LjLHpJapxFL0mNs+glqXEWvSQ1zqKXpMZZ9JLUOItekhpn0UtS4yx6SWrc\nQEWf5IYkjya5d93Y85PcnuTPe//+RG88Sf5DkgeS3JPkpeMKL0nqb9A9+huBXUeNXQncUVXnAHf0\nrgO8Ejin97UH+ODwMSVJ2zVQ0VfVncC3jhq+BPhw7/KHgV9cN/6RWnMXcGqS00cRVpJ07FJVgy2Y\nzAG3VtX5vet/WVWn9i4HeKKqTk1yK3BNVf1h77Y7gH9VVStHrW8Pa3v8zM7O7lxaWhrNTzRCq6ur\nzMzMTDrGQKYt6/6DT244PnsSHHr6OIfZpq5kPfuUE6bqd7+VabufbqYrORcXF/dV1Xy/5Z4zio1V\nVSUZ7C/G//+e64DrAObn52thYWEUUUZqeXmZacy1kWnLuvvKz204vnfHYa7dP5K73dh1JeuNu06e\nqt/9VqbtfrqZruQc1DBH3Rw6MiXT+/fR3vhB4Kx1y53ZG5MkTcAwRX8LcEXv8hXAZ9eN/5Pe0Tc/\nCzxZVY8MsR1J0hAGel6a5CZgATgtycPAO4FrgE8keSPwdeB1vcVvA14FPAB8F3j9iDNLko7BQEVf\nVZdtctNFGyxbwFuHCSVJGh3fGStJjbPoJalxFr0kNc6il6TGWfSS1DiLXpIaN/3v75a0pf0Hn9z0\nlBPrHbjm4uOQRtPIPXpJapxFL0mNs+glqXEWvSQ1zqKXpMZ51I2OydwAR3doOh3L784jdNriHr0k\nNc6il6TGWfSS1DiLXpIaZ9FLUuMseklqnEUvSY3b9nH0Sc4FPr5u6IXAvwFOBd4EfLM3flVV3bbt\nhJKkoWy76Kvqa8AFAElOAA4CnwZeD7yvqt47koSSpKGMaurmIuDBqvr6iNYnSRqRVNXwK0luAL5c\nVR9IcjWwG3gKWAH2VtUTG3zPHmAPwOzs7M6lpaWhc4za6uoqMzMzk44xkOOVdf/BJ4f6/tmT4NDT\nIwozZl3JOo6cO844ZbQr7OnKY6orORcXF/dV1Xy/5YYu+iTPA/4PcF5VHUoyCzwGFPAu4PSqesNW\n65ifn6+VlZWhcozD8vIyCwsLk44xkOOVddhz3ezdcZhr93fjFEtdyTqOnOM6101XHlNdyZlkoKIf\nxdTNK1nbmz8EUFWHquqZqnoW+BBw4Qi2IUnaplEU/WXATUeuJDl93W2vBe4dwTYkSds01PO9JCcD\nPw+8ed3we5JcwNrUzYGjbpMkHWdDFX1VfQf4yaPGLh8qkSRppHxnrCQ1zqKXpMZZ9JLUOItekhpn\n0UtS4yx6SWqcRS9JjbPoJalxFr0kNc6il6TGWfSS1DiLXpIaZ9FLUuMseklqnEUvSY2z6CWpcRa9\nJDXOopekxln0ktS4oT4zFiDJAeDbwDPA4aqaT/J84OPAHGsfEP66qnpi2G1Jko7dqPboF6vqgqqa\n712/Erijqs4B7uhdlyRNwLimbi4BPty7/GHgF8e0HUlSH6Mo+gI+n2Rfkj29sdmqeqR3+S+A2RFs\nR5K0Damq4VaQnFFVB5P8FHA78Dbglqo6dd0yT1TVTxz1fXuAPQCzs7M7l5aWhsoxDqurq8zMzEw6\nxkCGzbr/4JMjTLO52ZPg0NPHZVND60rWceTcccYpo11hT1ceU13Jubi4uG/dlPmmhi76H1pZcjWw\nCrwJWKiqR5KcDixX1bmbfd/8/HytrKyMLMeoLC8vs7CwMOkYAxk269yVnxtdmC3s3XGYa/cPfQzA\ncdGVrOPIeeCai0e6viO68pjqSs4kAxX9UFM3SU5O8uNHLgO/ANwL3AJc0VvsCuCzw2xHkrR9w+4G\nzAKfTnJkXf+5qv4gyZeATyR5I/B14HVDbkfScTToM7xx7flrtIYq+qp6CHjJBuOPAxcNs25J0mhM\n/wSkhnK85t4lTS9PgSBJjbPoJalxFr0kNc6il6TGWfSS1DiLXpIaZ9FLUuMseklqnEUvSY2z6CWp\ncRa9JDXOopekxln0ktQ4i16SGmfRS1LjLHpJapxFL0mNs+glqXEWvSQ1btufGZvkLOAjwCxQwHVV\n9f4kVwNvAr7ZW/Sqqrpt2KCSps+gn0l84JqLx5xEWxnmw8EPA3ur6stJfhzYl+T23m3vq6r3Dh9P\nkjSsbRd9VT0CPNK7/O0k9wNnjCqYJGk0UlXDrySZA+4Ezgf+BbAbeApYYW2v/4kNvmcPsAdgdnZ2\n59LS0tA5Rm11dZWZmZlJxxjIZln3H3xyAmk2N3sSHHp60ikG05WsXci544xTgO48prqSc3FxcV9V\nzfdbbuiiTzIDfAF4d1V9Ksks8Bhr8/bvAk6vqjdstY75+flaWVkZKsc4LC8vs7CwMOkYP2KjedG9\nOw5z7f5hZuKOj67khO5k7UpO6J91Wubyp/Wxf7QkAxX9UEfdJHkucDPwsar6FEBVHaqqZ6rqWeBD\nwIXDbEOSNJxtF32SANcD91fVb68bP33dYq8F7t1+PEnSsIZ5vvcy4HJgf5K7e2NXAZcluYC1qZsD\nwJuHSihJGsowR938IZANbvKYeUmaIt14BeeviEHffCJJx8JTIEhS4yx6SWqcRS9JjbPoJalxFr0k\nNc6il6TGeXilpKnh+e3Hwz16SWqcRS9JjXPq5jjwHa+SJsk9eklqnEUvSY2z6CWpcRa9JDXOopek\nxln0ktQ4D68cgodNSpMx6sde6++0dY9ekhpn0UtS48ZW9El2JflakgeSXDmu7UiStjaWOfokJwD/\nEfh54GHgS0luqaqvjGN7gzrWeb29Ow6z23l4qXlHd8NWj/0uzueP68XYC4EHquohgCRLwCXAWIre\nF0UlHS9dfCE4VTX6lSa/DOyqqn/au3458Heq6tfWLbMH2NO7ei7wtZEHGd5pwGOTDjGgrmTtSk7o\nTtau5ITuZO1Kzr9ZVS/ot9DEDq+squuA6ya1/UEkWamq+UnnGERXsnYlJ3Qna1dyQneydiXnoMb1\nYuxB4Kx118/sjUmSjrNxFf2XgHOSnJ3kecClwC1j2pYkaQtjmbqpqsNJfg34r8AJwA1Vdd84tjVm\nUz21dJSuZO1KTuhO1q7khO5k7UrOgYzlxVhJ0vTwnbGS1DiLXpIaZ9EPIMnbknw1yX1J3jPpPFtJ\nsjdJJTlt0lk2k+S3ev+f9yT5dJJTJ51pva6cviPJWUn+R5Kv9O6bvz7pTFtJckKSP0ly66SzbCXJ\nqUk+2buP3p/k704607As+j6SLLL2rt6XVNV5wHsnHGlTSc4CfgH435PO0sftwPlV9dPAnwFvn3Ce\nH1h3+o5XAi8GLkvy4smm2tRhYG9VvRj4WeCtU5wV4NeB+ycdYgDvB/6gqv4W8BK6kXlLFn1/vwpc\nU1XfA6iqRyecZyvvA34DmOpX2Kvq81V1uHf1LtbeZzEtfnD6jqr6PnDk9B1Tp6oeqaov9y5/m7VC\nOmOyqTaW5EzgYuB3J51lK0lOAV4BXA9QVd+vqr+cbKrhWfT9vQh4eZIvJvlCkp+ZdKCNJLkEOFhV\nfzrpLMfoDcDvTzrEOmcA31h3/WGmtDzXSzIH/G3gi5NNsql/z9pOyLOTDtLH2cA3gd/rTTP9bpKT\nJx1qWH7CFJDkvwF/Y4Ob3sHa/9HzWXtq/DPAJ5K8sCZwXGqfnFexNm0zFbbKWlWf7S3zDtamHz52\nPLO1JskMcDPwz6vqqUnnOVqSVwOPVtW+JAuTztPHc4CXAm+rqi8meT9wJfCvJxtrOBY9UFU/t9lt\nSX4V+FSv2P84ybOsnfDom8cr3xGb5Uyyg7U9kT9NAmtTIV9OcmFV/cVxjPgDW/2fAiTZDbwauGgS\nfzS30KnTdyR5Lmsl/7Gq+tSk82ziZcBrkrwK+DHgryf5aFX94wnn2sjDwMNVdeSZ0SdZK/pOc+qm\nv88AiwBJXgQ8jyk7q11V7a+qn6qquaqaY+3O+tJJlXw/SXax9jT+NVX13UnnOUpnTt+Rtb/q1wP3\nV9VvTzrPZqrq7VV1Zu++eSnw36e05Ok9Zr6R5Nze0EWM6fTqx5N79P3dANyQ5F7g+8AVU7YH2kUf\nAE4Ebu89A7mrqt4y2UhrOnb6jpcBlwP7k9zdG7uqqm6bYKYWvA34WO8P/UPA6yecZ2ieAkGSGufU\njSQ1zqKXpMZZ9JLUOItekhpn0UtS4yx6SWqcRS9Jjft/Fglp2spjL9QAAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x115f4fd10>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAE5FJREFUeJzt3X+MZfV93vH3EyA2YVLWLs6ILKsu\nUmgiwtbYjDApUTsT8gNDlHWkxMKiFGyqTVXc2u1KMXaqxlXqatsEu4mcom6MY1KIJwjbMiI4DSGM\nLP7ANosxy4+4WZu1zZYuwcbgwa6jxZ/+MWebyXZ27525d+bM/e77JY3mnu89955ntHef+d4z55yb\nqkKS1K7v6zuAJGl9WfSS1DiLXpIaZ9FLUuMseklqnEUvSY2z6CWpcRa9dIwktyX5/WPG/nGSryc5\nu69c0lrFE6akvy3J3wUeB66pqnuTvBJ4FPiPVfWRXsNJa2DRSytI8svAfwYuAP4tcGFVvbHfVNLa\nWPTScST5GHAacClLRf+1niNJa3Jq3wGkTexfAF8Cfs2S1yTzj7HScVTVYeA5lvbXSxPLopekxln0\nktQ4i16SGudRN5LUOGf0ktQ4i16SGmfRS1LjLHpJatymODP2rLPOqu3bt2/4dl966SXOOOOMDd/u\nKCYxM5h7o01i7knMDP3m3rdv33NV9ZpB622Kot++fTsPPfTQhm93YWGB2dnZDd/uKCYxM5h7o01i\n7knMDP3mTvKVYdZz140kNW5g0Sd5ZZLPJvlCkseT/Ptu/Nwkn0lyIMkfJfn+bvwV3fKB7v7t6/sj\nSJJOZJgZ/XeBn6qq1wIXApcnuQT4T8AHqupHgOeB67v1rwee78Y/0K0nSerJwKKvJYvd4mndVwE/\nBdzZjd8KvKm7vbNbprv/siQZW2JJ0qoMdQmEJKcA+4AfAX4X+E3gwW7WTpJtwKeq6oIkjwGXV9XT\n3X1fAt5QVc8d85y7gF0A09PTF83Pz4/vpxrS4uIiU1NTG77dUUxiZjD3RpvE3JOYGfrNPTc3t6+q\nZgauWFVDfwFbgPuBnwQOLBvfBjzW3X4MOGfZfV8CzjrR81500UXVh/vvv7+X7Y5iEjNXmXujTWLu\nScxc1W9u4KEaortXddRNVX2zK/qfALYkOXp45jnAoe72oa746e4/E/j6arYjSRqfYY66eU2SLd3t\n04GfAZ5kqfB/qVvtWuCT3e27umW6+/+8+80jSerBMCdMnQ3c2u2n/z7gjqq6O8kTwHyS/wB8Hril\nW/8W4L8nOQB8A7hqHXJLkoY0sOir6lHgdSuMfxm4eIXx/wP88ljSSau0/cY/Hmq9g3uuXOck0ubh\nmbGS1DiLXpIaZ9FLUuMseklqnEUvSY2z6CWpcRa9JDXOopekxln0ktQ4i16SGmfRS1LjLHpJapxF\nL0mNs+glqXEWvSQ1zqKXpMZZ9JLUuGE+SlA6qfmpVZp0zuglqXEWvSQ1zqKXpMZZ9JLUOItekhpn\n0UtS4yx6SWqcRS9JjRtY9Em2Jbk/yRNJHk/yjm78vUkOJXmk+7pi2WPeneRAki8m+bn1/AEkSSc2\nzJmxR4DdVfVwkh8E9iW5t7vvA1X1W8tXTnI+cBXw48APA3+W5O9X1cvjDC5JGs7AGX1VPVNVD3e3\nvwU8CWw9wUN2AvNV9d2qego4AFw8jrCSpNVLVQ2/crId+DRwAfBvgOuAF4GHWJr1P5/kg8CDVXVb\n95hbgE9V1Z3HPNcuYBfA9PT0RfPz86P+LKu2uLjI1NTUhm93FJOYGUbPvf/QC2NMAzu2njnUeouL\nizz1wnjfjA677VFM4utkEjNDv7nn5ub2VdXMoPWGvqhZkingY8A7q+rFJDcDvwFU9/0m4G3DPl9V\n7QX2AszMzNTs7OywDx2bhYUF+tjuKCYxM4ye+7ohLyw2rINXzw613sLCAjc98FIv2x7FJL5OJjEz\nTEbuoY66SXIaSyV/e1V9HKCqDlfVy1X1PeD3+JvdM4eAbcsefk43JknqwTBH3QS4BXiyqt6/bPzs\nZav9IvBYd/su4Kokr0hyLnAe8NnxRZYkrcYwu24uBa4B9id5pBt7D/CWJBeytOvmIPArAFX1eJI7\ngCdYOmLnBo+40WYz7DXmd+84gh/boEk38BVcVQ8AWeGue07wmPcB7xshlyRpTDwzVpIaZ9FLUuMs\neklqnEUvSY2z6CWpcRa9JDXOopekxln0ktQ4i16SGmfRS1LjLHpJapxFL0mNs+glqXEWvSQ1zqKX\npMZZ9JLUOItekhpn0UtS4yx6SWqcRS9JjfPj7dWr7Tf+cd8RpOY5o5ekxln0ktQ4i16SGmfRS1Lj\nLHpJatzAok+yLcn9SZ5I8niSd3Tjr05yb5K/7L6/qhtPkt9JciDJo0lev94/hCTp+IaZ0R8BdlfV\n+cAlwA1JzgduBO6rqvOA+7plgDcC53Vfu4Cbx55akjS0gUVfVc9U1cPd7W8BTwJbgZ3Ard1qtwJv\n6m7vBP6gljwIbEly9tiTS5KGkqoafuVkO/Bp4ALgq1W1pRsP8HxVbUlyN7Cnqh7o7rsPeFdVPXTM\nc+1iacbP9PT0RfPz86P/NKu0uLjI1NTUhm93FJOYGY6fe/+hF3pIM7zp0+Hwd/rZ9o6tZ675sZP4\nOpnEzNBv7rm5uX1VNTNovaHPjE0yBXwMeGdVvbjU7UuqqpIM/xtj6TF7gb0AMzMzNTs7u5qHj8XC\nwgJ9bHcUk5gZjp/7uk1+ZuzuHUe4aX8/J5AfvHp2zY+dxNfJJGaGycg91FE3SU5jqeRvr6qPd8OH\nj+6S6b4/240fArYte/g53ZgkqQfDHHUT4Bbgyap6/7K77gKu7W5fC3xy2fg/7Y6+uQR4oaqeGWNm\nSdIqDPOe9FLgGmB/kke6sfcAe4A7klwPfAV4c3ffPcAVwAHg28Bbx5pYkrQqA4u++6NqjnP3ZSus\nX8ANI+aSJI2JZ8ZKUuO8Hr20SQ17rf6De65c5ySadM7oJalxFr0kNc6il6TGWfSS1DiLXpIaZ9FL\nUuMseklqnEUvSY2z6CWpcRa9JDXOopekxln0ktQ4i16SGmfRS1LjLHpJapxFL0mNs+glqXEWvSQ1\nzqKXpMb5mbFaF8d+3unuHUe4bsjPQJU0Xs7oJalxzuilCXfsuydY+R3UwT1XblQkbTLO6CWpcRa9\nJDXOopekxg0s+iQfTvJskseWjb03yaEkj3RfVyy7791JDiT5YpKfW6/gkqThDDOj/whw+QrjH6iq\nC7uvewCSnA9cBfx495j/muSUcYWVJK3ewKKvqk8D3xjy+XYC81X13ap6CjgAXDxCPknSiFJVg1dK\ntgN3V9UF3fJ7geuAF4GHgN1V9XySDwIPVtVt3Xq3AJ+qqjtXeM5dwC6A6enpi+bn58fw46zO4uIi\nU1NTG77dUUxK5v2HXvhby9Onw+Hv9BRmBC3l3rH1zH7CDGlSXtvH6jP33NzcvqqaGbTeWo+jvxn4\nDaC67zcBb1vNE1TVXmAvwMzMTM3Ozq4xytotLCzQx3ZHMSmZjz2Ge/eOI9y0f/JO22gp98GrZ/sJ\nM6RJeW0faxJyr+mom6o6XFUvV9X3gN/jb3bPHAK2LVv1nG5MktSTNRV9krOXLf4icPSInLuAq5K8\nIsm5wHnAZ0eLKEkaxcD3pEk+CswCZyV5Gvh1YDbJhSztujkI/ApAVT2e5A7gCeAIcENVvbw+0SVJ\nwxhY9FX1lhWGbznB+u8D3jdKKEnS+EzeX5kkrclKFz9biRc/a4+XQJCkxln0ktQ4i16SGmfRS1Lj\nLHpJapxFL0mNs+glqXEWvSQ1zhOmtCrDnnQjafNwRi9JjbPoJalxFr0kNc6il6TGWfSS1DiLXpIa\nZ9FLUuMseklqnEUvSY2z6CWpcRa9JDXOopekxln0ktQ4i16SGmfRS1LjLHpJatzAok/y4STPJnls\n2dirk9yb5C+776/qxpPkd5IcSPJoktevZ3hJ0mDDzOg/Alx+zNiNwH1VdR5wX7cM8EbgvO5rF3Dz\neGJKktZqYNFX1aeBbxwzvBO4tbt9K/CmZeN/UEseBLYkOXtcYSVJq5eqGrxSsh24u6ou6Ja/WVVb\nutsBnq+qLUnuBvZU1QPdffcB76qqh1Z4zl0szfqZnp6+aH5+fjw/0SosLi4yNTW14dsdRd+Z9x96\nYU2Pmz4dDn9nzGE2wMmYe8fWM8cbZkh9v7bXqs/cc3Nz+6pqZtB6I384eFVVksG/Lf7/x+0F9gLM\nzMzU7OzsqFFWbWFhgT62O4q+M1+3xg8H373jCDftn7zPoj8Zcx+8ena8YYbU92t7rSYh91qPujl8\ndJdM9/3ZbvwQsG3Zeud0Y5Kknqx1qnIXcC2wp/v+yWXjb08yD7wBeKGqnhk5paQNs30V79oO7rly\nHZNoXAYWfZKPArPAWUmeBn6dpYK/I8n1wFeAN3er3wNcARwAvg28dR0yax2s5j+3pMkysOir6i3H\nueuyFdYt4IZRQ0mSxsczYyWpcRa9JDXOopekxln0ktQ4i16SGmfRS1LjLHpJapxFL0mNs+glqXEW\nvSQ1zqKXpMZZ9JLUOItekhpn0UtS4yx6SWrc5H0YpqRNY9gPrPGTqPrljF6SGueMvmF+PKAkcEYv\nSc2z6CWpcRa9JDXOopekxln0ktQ4i16SGmfRS1LjLHpJatxIJ0wlOQh8C3gZOFJVM0leDfwRsB04\nCLy5qp4fLaYkaa3GcWbsXFU9t2z5RuC+qtqT5MZu+V1j2I46nvEqaTXWY9fNTuDW7vatwJvWYRuS\npCGlqtb+4OQp4HmggP9WVXuTfLOqtnT3B3j+6PIxj90F7AKYnp6+aH5+fs051mpxcZGpqakN3+4o\nFhcXeeqFl/uOsWrTp8Ph7/SdYvXMPR47tp45cJ1J/P8I/eaem5vbV1Uzg9YbddfNT1bVoSQ/BNyb\n5C+W31lVlWTF3yRVtRfYCzAzM1Ozs7MjRlm9hYUF+tjuKBYWFrjpgZf6jrFqu3cc4ab9k3cNPXOP\nx8GrZweuM4n/H2Eyco/0SqiqQ933Z5N8ArgYOJzk7Kp6JsnZwLNjyClpgg3zd6XdO44wu/5RTkpr\n3kef5IwkP3j0NvCzwGPAXcC13WrXAp8cNaQkae1GmdFPA59Y2g3PqcAfVtWfJPkccEeS64GvAG8e\nPaYkaa3WXPRV9WXgtSuMfx24bJRQkqTx8cxYSWqcRS9JjbPoJalxFr0kNc6il6TGbZ5T5zT0SSX+\ns0laDWf0ktQ4i16SGmfRS1LjLHpJapxFL0mNs+glqXEWvSQ1zqKXpMZ55o2kTWOYkwYBDu65cp2T\ntMUZvSQ1zhm9pInjzH91nNFLUuMseklqnLtuNsCwbzMlaT04o5ekxjmjX4EzcEktcUYvSY1zRi+p\nWR6GucQZvSQ1zqKXpMZZ9JLUuHXbR5/kcuC3gVOAD1XVnvXa1rCO3V+3e8cRrvMIG+mktx778jfT\n3wfWpeiTnAL8LvAzwNPA55LcVVVPrMf2JGkjrFTekzBhXK8Z/cXAgar6MkCSeWAnMPai95h3STqx\nVNX4nzT5JeDyqvpn3fI1wBuq6u3L1tkF7OoWfxT44tiDDHYW8FwP2x3FJGYGc2+0Scw9iZmh39x/\nr6peM2il3o6jr6q9wN6+tg+Q5KGqmukzw2pNYmYw90abxNyTmBkmI/d6HXVzCNi2bPmcbkyStMHW\nq+g/B5yX5Nwk3w9cBdy1TtuSJJ3Auuy6qaojSd4O/A+WDq/8cFU9vh7bGlGvu47WaBIzg7k32iTm\nnsTMMAG51+WPsZKkzcMzYyWpcRa9JDXOogeS7E5SSc7qO8swkvxmkr9I8miSTyTZ0nemE0lyeZIv\nJjmQ5Ma+8wySZFuS+5M8keTxJO/oO9NqJDklyeeT3N13lmEl2ZLkzu51/WSSn+g70yBJ/nX3+ngs\nyUeTvLLvTMdz0hd9km3AzwJf7TvLKtwLXFBV/wD4n8C7e85zXMsuh/FG4HzgLUnO7zfVQEeA3VV1\nPnAJcMMEZF7uHcCTfYdYpd8G/qSqfgx4LZs8f5KtwL8CZqrqApYOOrmq31THd9IXPfAB4FeBifmr\ndFX9aVUd6RYfZOk8hc3q/10Oo6r+Gjh6OYxNq6qeqaqHu9vfYql0tvabajhJzgGuBD7Ud5ZhJTkT\n+EfALQBV9ddV9c1+Uw3lVOD0JKcCPwD8r57zHNdJXfRJdgKHquoLfWcZwduAT/Ud4gS2Al9btvw0\nE1KaAEm2A68DPtNvkqH9F5YmLt/rO8gqnAv8FfD73S6nDyU5o+9QJ1JVh4DfYmlPwDPAC1X1p/2m\nOr7miz7Jn3X70I792gm8B/h3fWdcyYDcR9f5NZZ2M9zeX9J2JZkCPga8s6pe7DvPIEl+Hni2qvb1\nnWWVTgVeD9xcVa8DXgI29d9ykryKpXem5wI/DJyR5J/0m+r4mv/M2Kr66ZXGk+xg6R/pC0lgaffH\nw0kurqr/vYERV3S83EcluQ74eeCy2twnQ0zk5TCSnMZSyd9eVR/vO8+QLgV+IckVwCuBv5Pktqra\ntAXUeRp4uqqOvmu6k01e9MBPA09V1V8BJPk48A+B23pNdRzNz+iPp6r2V9UPVdX2qtrO0ovt9Zuh\n5AfpPtTlV4FfqKpv951ngIm7HEaWfvPfAjxZVe/vO8+wqurdVXVO93q+CvjzCSh5uv9zX0vyo93Q\nZazDJc3H7KvAJUl+oHu9XMYm/gNy8zP6Rn0QeAVwb/du5MGq+uf9RlrZBF0OY7lLgWuA/Uke6cbe\nU1X39Jipdf8SuL2bDHwZeGvPeU6oqj6T5E7gYZZ2n36eTXwpBC+BIEmNO2l33UjSycKil6TGWfSS\n1DiLXpIaZ9FLUuMseklqnEUvSY37v2spfgcVLX/LAAAAAElFTkSuQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x115be4d90>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"df.groupby('Z').hist('Y', bins=30)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The support isn't great on the high and low side of Y, but the density there is small, and we'll plan to extrapolate with our models. That's always something to be concerned about! If the extrapolation fails, your estimator is probably biased."
]
},
{
"cell_type": "code",
"execution_count": 188,
"metadata": {
"collapsed": false
},
"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>D</th>\n",
" <th>Y</th>\n",
" <th>Z</th>\n",
" <th>$P(D|Z)$</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>0</td>\n",
" <td>1.813351</td>\n",
" <td>1</td>\n",
" <td>0.501965</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>1</td>\n",
" <td>1.207073</td>\n",
" <td>0</td>\n",
" <td>0.724348</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>1</td>\n",
" <td>0.317765</td>\n",
" <td>0</td>\n",
" <td>0.724348</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>1</td>\n",
" <td>1.855845</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>1</td>\n",
" <td>4.678077</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" D Y Z $P(D|Z)$\n",
"0 0 1.813351 1 0.501965\n",
"1 1 1.207073 0 0.724348\n",
"2 1 0.317765 0 0.724348\n",
"3 1 1.855845 1 0.724348\n",
"4 1 4.678077 1 0.724348"
]
},
"execution_count": 188,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p_z_given_d = LogisticRegression(C=1e12)\n",
"p_z_given_d = p_z_given_x.fit(df[['D']], df['Z'])\n",
"df[\"$P(D|Z)$\"] = p_z_given_d.predict_proba(df[['D']])[:,1]\n",
"df.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let's get $E[Y|D', Z]$"
]
},
{
"cell_type": "code",
"execution_count": 189,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"y_given_d_z = LinearRegression()\n",
"y_given_d_z = y_given_d_z.fit(df[['D', 'Z']], df['Y'])\n",
"df[\"$E[Y| D', Z]$\"] = y_given_d_z.predict(df[['D', 'Z']])"
]
},
{
"cell_type": "code",
"execution_count": 190,
"metadata": {
"collapsed": false,
"scrolled": true
},
"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>D</th>\n",
" <th>Y</th>\n",
" <th>Z</th>\n",
" <th>$P(D|Z)$</th>\n",
" <th>$E[Y| D', Z]$</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>0</td>\n",
" <td>1.813351</td>\n",
" <td>1</td>\n",
" <td>0.501965</td>\n",
" <td>1.105574</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>1</td>\n",
" <td>1.207073</td>\n",
" <td>0</td>\n",
" <td>0.724348</td>\n",
" <td>0.869490</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>1</td>\n",
" <td>0.317765</td>\n",
" <td>0</td>\n",
" <td>0.724348</td>\n",
" <td>0.869490</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>1</td>\n",
" <td>1.855845</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" <td>2.815666</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>1</td>\n",
" <td>4.678077</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" <td>2.815666</td>\n",
" </tr>\n",
" <tr>\n",
" <th>5</th>\n",
" <td>0</td>\n",
" <td>2.862237</td>\n",
" <td>1</td>\n",
" <td>0.501965</td>\n",
" <td>1.105574</td>\n",
" </tr>\n",
" <tr>\n",
" <th>6</th>\n",
" <td>1</td>\n",
" <td>1.914817</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" <td>2.815666</td>\n",
" </tr>\n",
" <tr>\n",
" <th>7</th>\n",
" <td>1</td>\n",
" <td>0.391702</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" <td>2.815666</td>\n",
" </tr>\n",
" <tr>\n",
" <th>8</th>\n",
" <td>0</td>\n",
" <td>-1.753199</td>\n",
" <td>0</td>\n",
" <td>0.501965</td>\n",
" <td>-0.840601</td>\n",
" </tr>\n",
" <tr>\n",
" <th>9</th>\n",
" <td>0</td>\n",
" <td>1.555124</td>\n",
" <td>1</td>\n",
" <td>0.501965</td>\n",
" <td>1.105574</td>\n",
" </tr>\n",
" <tr>\n",
" <th>10</th>\n",
" <td>1</td>\n",
" <td>2.222936</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" <td>2.815666</td>\n",
" </tr>\n",
" <tr>\n",
" <th>11</th>\n",
" <td>1</td>\n",
" <td>2.152429</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" <td>2.815666</td>\n",
" </tr>\n",
" <tr>\n",
" <th>12</th>\n",
" <td>0</td>\n",
" <td>-4.815750</td>\n",
" <td>0</td>\n",
" <td>0.501965</td>\n",
" <td>-0.840601</td>\n",
" </tr>\n",
" <tr>\n",
" <th>13</th>\n",
" <td>1</td>\n",
" <td>-1.869212</td>\n",
" <td>0</td>\n",
" <td>0.724348</td>\n",
" <td>0.869490</td>\n",
" </tr>\n",
" <tr>\n",
" <th>14</th>\n",
" <td>1</td>\n",
" <td>1.575286</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" <td>2.815666</td>\n",
" </tr>\n",
" <tr>\n",
" <th>15</th>\n",
" <td>0</td>\n",
" <td>0.775557</td>\n",
" <td>1</td>\n",
" <td>0.501965</td>\n",
" <td>1.105574</td>\n",
" </tr>\n",
" <tr>\n",
" <th>16</th>\n",
" <td>1</td>\n",
" <td>3.589573</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" <td>2.815666</td>\n",
" </tr>\n",
" <tr>\n",
" <th>17</th>\n",
" <td>1</td>\n",
" <td>3.947033</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" <td>2.815666</td>\n",
" </tr>\n",
" <tr>\n",
" <th>18</th>\n",
" <td>0</td>\n",
" <td>-2.040981</td>\n",
" <td>0</td>\n",
" <td>0.501965</td>\n",
" <td>-0.840601</td>\n",
" </tr>\n",
" <tr>\n",
" <th>19</th>\n",
" <td>0</td>\n",
" <td>-1.035452</td>\n",
" <td>0</td>\n",
" <td>0.501965</td>\n",
" <td>-0.840601</td>\n",
" </tr>\n",
" <tr>\n",
" <th>20</th>\n",
" <td>1</td>\n",
" <td>9.153940</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" <td>2.815666</td>\n",
" </tr>\n",
" <tr>\n",
" <th>21</th>\n",
" <td>1</td>\n",
" <td>3.093233</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" <td>2.815666</td>\n",
" </tr>\n",
" <tr>\n",
" <th>22</th>\n",
" <td>1</td>\n",
" <td>0.585575</td>\n",
" <td>0</td>\n",
" <td>0.724348</td>\n",
" <td>0.869490</td>\n",
" </tr>\n",
" <tr>\n",
" <th>23</th>\n",
" <td>1</td>\n",
" <td>2.369515</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" <td>2.815666</td>\n",
" </tr>\n",
" <tr>\n",
" <th>24</th>\n",
" <td>0</td>\n",
" <td>1.512505</td>\n",
" <td>0</td>\n",
" <td>0.501965</td>\n",
" <td>-0.840601</td>\n",
" </tr>\n",
" <tr>\n",
" <th>25</th>\n",
" <td>1</td>\n",
" <td>3.437692</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" <td>2.815666</td>\n",
" </tr>\n",
" <tr>\n",
" <th>26</th>\n",
" <td>1</td>\n",
" <td>0.649377</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" <td>2.815666</td>\n",
" </tr>\n",
" <tr>\n",
" <th>27</th>\n",
" <td>0</td>\n",
" <td>1.113973</td>\n",
" <td>0</td>\n",
" <td>0.501965</td>\n",
" <td>-0.840601</td>\n",
" </tr>\n",
" <tr>\n",
" <th>28</th>\n",
" <td>1</td>\n",
" <td>2.839603</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" <td>2.815666</td>\n",
" </tr>\n",
" <tr>\n",
" <th>29</th>\n",
" <td>0</td>\n",
" <td>-4.086222</td>\n",
" <td>0</td>\n",
" <td>0.501965</td>\n",
" <td>-0.840601</td>\n",
" </tr>\n",
" <tr>\n",
" <th>...</th>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4970</th>\n",
" <td>0</td>\n",
" <td>-1.334995</td>\n",
" <td>0</td>\n",
" <td>0.501965</td>\n",
" <td>-0.840601</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4971</th>\n",
" <td>0</td>\n",
" <td>-3.746823</td>\n",
" <td>0</td>\n",
" <td>0.501965</td>\n",
" <td>-0.840601</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4972</th>\n",
" <td>1</td>\n",
" <td>2.555232</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" <td>2.815666</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4973</th>\n",
" <td>1</td>\n",
" <td>2.639780</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" <td>2.815666</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4974</th>\n",
" <td>0</td>\n",
" <td>-1.421910</td>\n",
" <td>1</td>\n",
" <td>0.501965</td>\n",
" <td>1.105574</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4975</th>\n",
" <td>0</td>\n",
" <td>-0.019750</td>\n",
" <td>0</td>\n",
" <td>0.501965</td>\n",
" <td>-0.840601</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4976</th>\n",
" <td>0</td>\n",
" <td>0.228803</td>\n",
" <td>0</td>\n",
" <td>0.501965</td>\n",
" <td>-0.840601</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4977</th>\n",
" <td>0</td>\n",
" <td>-1.248309</td>\n",
" <td>1</td>\n",
" <td>0.501965</td>\n",
" <td>1.105574</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4978</th>\n",
" <td>0</td>\n",
" <td>0.938322</td>\n",
" <td>1</td>\n",
" <td>0.501965</td>\n",
" <td>1.105574</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4979</th>\n",
" <td>0</td>\n",
" <td>-0.113057</td>\n",
" <td>0</td>\n",
" <td>0.501965</td>\n",
" <td>-0.840601</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4980</th>\n",
" <td>0</td>\n",
" <td>1.210678</td>\n",
" <td>0</td>\n",
" <td>0.501965</td>\n",
" <td>-0.840601</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4981</th>\n",
" <td>1</td>\n",
" <td>-1.036083</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" <td>2.815666</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4982</th>\n",
" <td>0</td>\n",
" <td>-0.347940</td>\n",
" <td>1</td>\n",
" <td>0.501965</td>\n",
" <td>1.105574</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4983</th>\n",
" <td>0</td>\n",
" <td>-0.098543</td>\n",
" <td>0</td>\n",
" <td>0.501965</td>\n",
" <td>-0.840601</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4984</th>\n",
" <td>1</td>\n",
" <td>7.214081</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" <td>2.815666</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4985</th>\n",
" <td>1</td>\n",
" <td>8.972760</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" <td>2.815666</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4986</th>\n",
" <td>1</td>\n",
" <td>2.063666</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" <td>2.815666</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4987</th>\n",
" <td>1</td>\n",
" <td>4.867444</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" <td>2.815666</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4988</th>\n",
" <td>1</td>\n",
" <td>3.401611</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" <td>2.815666</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4989</th>\n",
" <td>0</td>\n",
" <td>-1.991439</td>\n",
" <td>0</td>\n",
" <td>0.501965</td>\n",
" <td>-0.840601</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4990</th>\n",
" <td>1</td>\n",
" <td>0.220615</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" <td>2.815666</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4991</th>\n",
" <td>0</td>\n",
" <td>-1.517893</td>\n",
" <td>0</td>\n",
" <td>0.501965</td>\n",
" <td>-0.840601</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4992</th>\n",
" <td>1</td>\n",
" <td>3.627046</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" <td>2.815666</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4993</th>\n",
" <td>1</td>\n",
" <td>1.361454</td>\n",
" <td>0</td>\n",
" <td>0.724348</td>\n",
" <td>0.869490</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4994</th>\n",
" <td>1</td>\n",
" <td>2.889220</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" <td>2.815666</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4995</th>\n",
" <td>0</td>\n",
" <td>1.348715</td>\n",
" <td>0</td>\n",
" <td>0.501965</td>\n",
" <td>-0.840601</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4996</th>\n",
" <td>1</td>\n",
" <td>2.383191</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" <td>2.815666</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4997</th>\n",
" <td>0</td>\n",
" <td>-0.533678</td>\n",
" <td>1</td>\n",
" <td>0.501965</td>\n",
" <td>1.105574</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4998</th>\n",
" <td>1</td>\n",
" <td>1.741694</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" <td>2.815666</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4999</th>\n",
" <td>1</td>\n",
" <td>2.364081</td>\n",
" <td>1</td>\n",
" <td>0.724348</td>\n",
" <td>2.815666</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"<p>5000 rows × 5 columns</p>\n",
"</div>"
],
"text/plain": [
" D Y Z $P(D|Z)$ $E[Y| D', Z]$\n",
"0 0 1.813351 1 0.501965 1.105574\n",
"1 1 1.207073 0 0.724348 0.869490\n",
"2 1 0.317765 0 0.724348 0.869490\n",
"3 1 1.855845 1 0.724348 2.815666\n",
"4 1 4.678077 1 0.724348 2.815666\n",
"5 0 2.862237 1 0.501965 1.105574\n",
"6 1 1.914817 1 0.724348 2.815666\n",
"7 1 0.391702 1 0.724348 2.815666\n",
"8 0 -1.753199 0 0.501965 -0.840601\n",
"9 0 1.555124 1 0.501965 1.105574\n",
"10 1 2.222936 1 0.724348 2.815666\n",
"11 1 2.152429 1 0.724348 2.815666\n",
"12 0 -4.815750 0 0.501965 -0.840601\n",
"13 1 -1.869212 0 0.724348 0.869490\n",
"14 1 1.575286 1 0.724348 2.815666\n",
"15 0 0.775557 1 0.501965 1.105574\n",
"16 1 3.589573 1 0.724348 2.815666\n",
"17 1 3.947033 1 0.724348 2.815666\n",
"18 0 -2.040981 0 0.501965 -0.840601\n",
"19 0 -1.035452 0 0.501965 -0.840601\n",
"20 1 9.153940 1 0.724348 2.815666\n",
"21 1 3.093233 1 0.724348 2.815666\n",
"22 1 0.585575 0 0.724348 0.869490\n",
"23 1 2.369515 1 0.724348 2.815666\n",
"24 0 1.512505 0 0.501965 -0.840601\n",
"25 1 3.437692 1 0.724348 2.815666\n",
"26 1 0.649377 1 0.724348 2.815666\n",
"27 0 1.113973 0 0.501965 -0.840601\n",
"28 1 2.839603 1 0.724348 2.815666\n",
"29 0 -4.086222 0 0.501965 -0.840601\n",
"... .. ... .. ... ...\n",
"4970 0 -1.334995 0 0.501965 -0.840601\n",
"4971 0 -3.746823 0 0.501965 -0.840601\n",
"4972 1 2.555232 1 0.724348 2.815666\n",
"4973 1 2.639780 1 0.724348 2.815666\n",
"4974 0 -1.421910 1 0.501965 1.105574\n",
"4975 0 -0.019750 0 0.501965 -0.840601\n",
"4976 0 0.228803 0 0.501965 -0.840601\n",
"4977 0 -1.248309 1 0.501965 1.105574\n",
"4978 0 0.938322 1 0.501965 1.105574\n",
"4979 0 -0.113057 0 0.501965 -0.840601\n",
"4980 0 1.210678 0 0.501965 -0.840601\n",
"4981 1 -1.036083 1 0.724348 2.815666\n",
"4982 0 -0.347940 1 0.501965 1.105574\n",
"4983 0 -0.098543 0 0.501965 -0.840601\n",
"4984 1 7.214081 1 0.724348 2.815666\n",
"4985 1 8.972760 1 0.724348 2.815666\n",
"4986 1 2.063666 1 0.724348 2.815666\n",
"4987 1 4.867444 1 0.724348 2.815666\n",
"4988 1 3.401611 1 0.724348 2.815666\n",
"4989 0 -1.991439 0 0.501965 -0.840601\n",
"4990 1 0.220615 1 0.724348 2.815666\n",
"4991 0 -1.517893 0 0.501965 -0.840601\n",
"4992 1 3.627046 1 0.724348 2.815666\n",
"4993 1 1.361454 0 0.724348 0.869490\n",
"4994 1 2.889220 1 0.724348 2.815666\n",
"4995 0 1.348715 0 0.501965 -0.840601\n",
"4996 1 2.383191 1 0.724348 2.815666\n",
"4997 0 -0.533678 1 0.501965 1.105574\n",
"4998 1 1.741694 1 0.724348 2.815666\n",
"4999 1 2.364081 1 0.724348 2.815666\n",
"\n",
"[5000 rows x 5 columns]"
]
},
"execution_count": 190,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df"
]
},
{
"cell_type": "code",
"execution_count": 191,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"D\n",
"0 0.5088\n",
"1 0.4912\n",
"Name: Y, dtype: float64"
]
},
"execution_count": 191,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p_d = df.groupby('D').count()['Y'] / len(df)\n",
"p_d"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's do D=1. We see $P(D=1) = 0.485$"
]
},
{
"cell_type": "code",
"execution_count": 192,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# compute E[Y|do(Z=z)] for each z. We need to perform the sum over D'.\n",
"from collections import defaultdict\n",
"\n",
"E_y_do_z = defaultdict(lambda : 0.)\n",
"\n",
"for zi in df['Z'].unique():\n",
" for di in df['D'].unique():\n",
" E_y_do_z[zi] += y_given_d_z.predict([[di, zi]])[0] * p_d[di]\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 193,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Now, we can evaluate the outer sum\n",
"E_y_do_d = defaultdict(lambda : 0.)\n",
"for di in df['D'].unique():\n",
" for zi in df['Z'].unique():\n",
" E_y_do_d[di] += E_y_do_z[zi] * p_z_given_d.predict_proba([[di]])[0][1]\n"
]
},
{
"cell_type": "code",
"execution_count": 194,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"defaultdict(<function __main__.<lambda>>,\n",
" {0: 0.9763059836877175, 1: 1.40883365975047})"
]
},
"execution_count": 194,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"E_y_do_d"
]
},
{
"cell_type": "code",
"execution_count": 195,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# We can check by intervening! Then, the naive estimator should work.\n",
"\n",
"N = 5000\n",
"u = np.random.normal(size=N)\n",
"\n",
"p_d = 0.5 # 1./(1. + np.exp(-u)) replace with random assignment to break confounding, keep the same marginal of p(D)\n",
"d = np.random.binomial(1, p=p_d, size=N)\n",
"\n",
"p_z = 1./(1. + np.exp(-(d + 0.1 * np.random.normal(size=N))))\n",
"z = np.random.binomial(1, p=p_z)\n",
"\n",
"y = 2. * u + 2.*z + 0.2* np.random.normal(size=N) \n",
"df = pd.DataFrame({'D': d, 'Y': y, 'Z': z})"
]
},
{
"cell_type": "code",
"execution_count": 199,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"1.433817023841926"
]
},
"execution_count": 199,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df[df['D'] == 1].mean()['Y']"
]
},
{
"cell_type": "code",
"execution_count": 198,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"0.9729122633926304"
]
},
"execution_count": 198,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df[df['D'] == 0].mean()['Y']"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.11"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment