selection_bias_example_2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now, let's do the estimation for the article translation problem.\n",
    "\n",
    "We'll create some data with the same structure as our causal graph. Notice that $U$ is the empty set. We'll add another variable, $d$, which is the article's author's skill level, $d$. Higher skill means better performance. Skill, $d$, causes article quality, $x$. article quality drives both pre- and post-translation performance ($y$ and $yt$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "N=5000\n",
    "\n",
    "x = np.random.choice(range(5), size=N)\n",
    "\n",
    "u = np.random.normal(size=N)\n",
    "\n",
    "y = 100* np.random.normal(x + u)\n",
    "\n",
    "p_s = 1. / (1. + np.exp(-(y - 200)/200 + 5.*(x - 2.5) / 2.5))\n",
    "s = np.random.binomial(1, p=p_s)\n",
    "\n",
    "yt = 100*np.random.normal(x + u)\n",
    "\n",
    "X = pd.DataFrame({'X': x, 'Y': y, 'Yt': yt, 'S': s})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "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>S</th>\n",
       "      <th>X</th>\n",
       "      <th>Y</th>\n",
       "      <th>Yt</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>258.777688</td>\n",
       "      <td>-4.492901</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>3</td>\n",
       "      <td>504.236178</td>\n",
       "      <td>397.600212</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>270.548991</td>\n",
       "      <td>202.441660</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>-12.801746</td>\n",
       "      <td>-18.454125</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>321.775170</td>\n",
       "      <td>91.687487</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   S  X           Y          Yt\n",
       "0  1  2  258.777688   -4.492901\n",
       "1  1  3  504.236178  397.600212\n",
       "2  1  2  270.548991  202.441660\n",
       "3  1  1  -12.801746  -18.454125\n",
       "4  1  2  321.775170   91.687487"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_sampled = X[X['S'] == 1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now that we have our data, let's compare the sample E[Yt|D] with the population E[Yt|D]. We wouldn't normally be able to make this comparison, but since we generated the data we have the population E[Yt|D]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7f5c3e8cf470>"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEGCAYAAABrQF4qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAADxhJREFUeJzt3XtsXvV5wPHvE3LrgHExFsripI4EErB0S4NFmSI0BmIjMC5/0AKrmgixRhOgwTJtzYgmqLRJQUxjqwRFUUMbplEKbBJRi1ohLp26qQyHZtzSjcBI4yhAmkK4NS2XZ3+8vzTGc/Br+7WP/fP3I0U+53eO/T68gm8Ox+9rR2YiSarXrKYHkCRNLEMvSZUz9JJUOUMvSZUz9JJUOUMvSZUz9JJUOUMvSZUz9JJUudlNDwBwwgknZG9vb9NjSNK0snXr1p9mZvdI502J0Pf29tLf39/0GJI0rUTEznbO89aNJFXO0EtS5Qy9JFVuStyjH857773HwMAABw4caHqUts2fP5+enh7mzJnT9CiS9CtTNvQDAwMcffTR9Pb2EhFNjzOizGTfvn0MDAywZMmSpseRpF+ZsrduDhw4QFdX17SIPEBE0NXVNa3+D0TSzDBlQw9Mm8gfNN3mlTQzTOnQS5LGb8reox+qd913Ovr1Xt5w4ccez0zOOuss1q9fz8qVKwG4//772bRpExdffDHXXHNNR+eRNHlu/5NHmx6Ba+88Z9Ieyyv6w4gI7rzzTtauXcuBAwd4++23ufHGG7n99tu54447mh5Pkto2ba7om7B06VIuuugibrnlFt555x1WrVrF+vXrefHFF1m2bBnnnXcet956a9NjStLHMvQjuOmmm1i+fDlz586lv7+fPXv28Oyzz7Jt27amR5Okthj6ERx55JFcfvnlHHXUUcybN6/pcSRp1LxH34ZZs2Yxa5ZPlaTpyXqN0tFHH81bb73V9BiS1LZpc+tmpJdDTpauri5WrFjB0qVLWblypd+MlTTlTZvQN+nmm2/+yP4999zTzCCSNAbeupGkyhl6SarclA59ZjY9wqhMt3klzQxTNvTz589n37590yaeB38e/fz585seRZI+Ysp+M7anp4eBgQH27t3b9ChtO/gbpiRpKpmyoZ8zZ46/qUmSOmDK3rqRJHWGoZekyrUd+og4IiJ+FBHfLvtLIuKJiNgREd+KiLllfV7Z31GO907M6JKkdozmiv56YPug/VuA2zLzJOB14OqyfjXwelm/rZwnSWpIW6GPiB7gQuBrZT+Ac4AHyimbgUvL9iVln3L83PC3ZktSY9q9ov8H4C+BD8t+F/BGZr5f9geAhWV7IbALoBzfX87/iIhYExH9EdE/nV5CKUnTzYihj4g/BF7LzK2dfODM3JiZfZnZ193d3ckvLUkapJ3X0a8ALo6IC4D5wK8D/wgcGxGzy1V7D7C7nL8bWAQMRMRs4BhgX8cnlyS1ZcQr+sz8q8zsycxe4Arg0cz8PPAYcFk5bTXwYNneUvYpxx/N6fJzDCSpQuN5Hf2XgLURsYPWPfhNZX0T0FXW1wLrxjeiJGk8RvUjEDLzceDxsv0ScMYw5xwAPtuB2SRJHeA7YyWpcoZekipn6CWpcoZekipn6CWpcoZekipn6CWpcoZekipn6CWpcoZekipn6CWpcoZekipn6CWpcoZekipn6CWpcoZekipn6CWpcoZekipn6CWpcoZekipn6CWpcoZekipn6CWpcoZekipn6CWpcoZekipn6CWpcoZekio3u+kBJE2O7aec2vQInPrj7U2PMCN5RS9JlTP0klQ5Qy9JlTP0klQ5Qy9JlTP0klS5EUMfEfMj4j8j4r8i4rmI+HJZXxIRT0TEjoj4VkTMLevzyv6Ocrx3Yv8RJEkfp50r+l8A52TmbwPLgPMj4kzgFuC2zDwJeB24upx/NfB6Wb+tnCdJasiIoc+Wt8vunPIngXOAB8r6ZuDSsn1J2accPzciomMTS5JGpa179BFxRERsA14DHgZeBN7IzPfLKQPAwrK9ENgFUI7vB7o6ObQkqX1thT4zP8jMZUAPcAZwyngfOCLWRER/RPTv3bt3vF9OknQYo3rVTWa+ATwG/A5wbEQc/Fk5PcDusr0bWARQjh8D7Bvma23MzL7M7Ovu7h7j+JKkkbTzqpvuiDi2bH8COA/YTiv4l5XTVgMPlu0tZZ9y/NHMzE4OLUlqXzs/vXIBsDkijqD1F8N9mfntiHgeuDci/gb4EbCpnL8J+KeI2AH8DLhiAuaWJLVpxNBn5tPAp4dZf4nW/fqh6weAz3ZkOknSuPnOWEmqnKGXpMoZekmqnKGXpMoZekmqnKGXpMoZekmqnKGXpMoZekmqnKGXpMoZekmqnKGXpMoZekmqnKGXpMoZekmqnKGXpMoZekmqnKGXpMoZekmqnKGXpMoZekmqnKGXpMoZekmqnKGXpMoZekmqnKGXpMoZekmqnKGXpMoZekmqnKGXpMoZekmqnKGXpMoZekmq3OymB5Am0qc2f6rpEXhm9TNNj6AZzit6SaqcoZekyo0Y+ohYFBGPRcTzEfFcRFxf1o+PiIcj4oXy8biyHhHxlYjYERFPR8Tyif6HkCQdXjtX9O8Df56ZpwFnAtdGxGnAOuCRzDwZeKTsA6wETi5/1gBf7fjUkqS2jRj6zNyTmU+V7beA7cBC4BJgczltM3Bp2b4EuDtbfggcGxELOj65JKkto7pHHxG9wKeBJ4ATM3NPOfQKcGLZXgjsGvRpA2Vt6NdaExH9EdG/d+/eUY4tSWpX26GPiKOAfwFuyMw3Bx/LzARyNA+cmRszsy8z+7q7u0fzqZKkUWgr9BExh1bk/zkz/7Usv3rwlkz5+FpZ3w0sGvTpPWVNktSAdl51E8AmYHtm/v2gQ1uA1WV7NfDgoPVV5dU3ZwL7B93ikSRNsnbeGbsC+ALwTERsK2s3AhuA+yLiamAn8Lly7CHgAmAH8C5wVUcnliSNyoihz8wfAHGYw+cOc34C145zLklSh/jOWEmqnKGXpMoZekmqnKGXpMoZekmqnKGXpMoZekmqnKGXpMoZekmqnKGXpMoZekmqnKGXpMoZekmqnKGXpMoZekmqnKGXpMoZekmqnKGXpMoZekmqnKGXpMoZekmqnKGXpMoZekmqnKGXpMoZekmqnKGXpMoZekmqnKGXpMoZekmqnKGXpMoZekmqnKGXpMoZekmqnKGXpMoZekmq3Iihj4i7IuK1iHh20NrxEfFwRLxQPh5X1iMivhIROyLi6YhYPpHDS5JG1s4V/TeA84esrQMeycyTgUfKPsBK4OTyZw3w1c6MKUkaqxFDn5n/BvxsyPIlwOayvRm4dND63dnyQ+DYiFjQqWElSaM3e4yfd2Jm7inbrwAnlu2FwK5B5w2UtT0MERFraF31s3jx4jGOoWHdfEzTE8DN+5ueQFIx7m/GZmYCOYbP25iZfZnZ193dPd4xJEmHMdbQv3rwlkz5+FpZ3w0sGnReT1mTJDVkrKHfAqwu26uBBwetryqvvjkT2D/oFo8kqQEj3qOPiG8CZwMnRMQAcBOwAbgvIq4GdgKfK6c/BFwA7ADeBa6agJklSaMwYugz88rDHDp3mHMTuHa8Q0mSOsd3xkpS5Qy9JFXO0EtS5Qy9JFXO0EtS5Qy9JFXO0EtS5Qy9JFXO0EtS5Qy9JFXO0EtS5Qy9JFXO0EtS5Qy9JFXO0EtS5Qy9JFXO0EtS5Qy9JFXO0EtS5Qy9JFXO0EtS5Qy9JFXO0EtS5Qy9JFXO0EtS5Qy9JFXO0EtS5Qy9JFXO0EtS5Qy9JFXO0EtS5Qy9JFXO0EtS5Qy9JFXO0EtS5SYk9BFxfkT8d0TsiIh1E/EYkqT2dDz0EXEEcDuwEjgNuDIiTuv040iS2jMRV/RnADsy86XM/CVwL3DJBDyOJKkNkZmd/YIRlwHnZ+Yfl/0vAJ/JzOuGnLcGWAOwePHi03fu3Dmux+1d951xfX4nvLzhwqZHkDSDRMTWzOwb6bzZkzHMcDJzI7ARoK+vb9x/2xhZSRreRNy62Q0sGrTfU9YkSQ2YiNA/CZwcEUsiYi5wBbBlAh5HktSGjt+6ycz3I+I64HvAEcBdmflcpx9HktSeCblHn5kPAQ9NxNeWJI2O74yVpMoZekmqnKGXpMoZekmqXMffGTumISL2AuN7a2xnnAD8tOkhpgifixafh0N8Lg6ZKs/FJzOze6STpkTop4qI6G/n7cQzgc9Fi8/DIT4Xh0y358JbN5JUOUMvSZUz9B+1sekBphCfixafh0N8Lg6ZVs+F9+glqXJe0UtS5Qy9JFXO0EtS5Rr7DVNNi4hTaP0u24VlaTewJTO3NzeVNHVExBlAZuaTEXEacD7w4/LTaWe0iLg7M1c1PUe7ZuQ3YyPiS8CVtH5x+UBZ7qH1S1LuzcwNTc2mZpULgIXAE5n59qD18zPzu81NNrki4iZgJa2LwYeBzwCPAecB38vMv21wvEkVEUN/cVIAvwc8CpCZF0/6UKM0U0P/P8BvZuZ7Q9bnAs9l5snNTDb1RMRVmfn1pueYDBHxp8C1wHZgGXB9Zj5Yjj2VmcubnG8yRcQztJ6DecArQE9mvhkRn6D1l+BvNTrgJIqIp4Dnga8BSSv036R1YUhmfr+56dozU+/Rfwj8xjDrC8oxHfLlpgeYRF8ETs/MS4Gzgb+OiOvLsWhsqma8n5kfZOa7wIuZ+SZAZv6cmfffSB+wFVgP7M/Mx4GfZ+b3p0PkYebeo78BeCQiXgB2lbXFwEnAdY1N1ZCIePpwh4ATJ3OWhs06eLsmM1+OiLOBByLik8y80P8yIn6thP70g4sRcQwzLPSZ+SFwW0TcXz6+yjRr54y8dQMQEbOAM/joN2OfzMwPmpuqGeVf3D8AXh96CPiPzBzu/36qExGPAmszc9ugtdnAXcDnM/OIxoabZBExLzN/Mcz6CcCCzHymgbGmhIi4EFiRmTc2PUu7ZmzodUhEbAK+npk/GObYPZn5Rw2MNekioofWLYtXhjm2IjP/vYGxpHEz9JJUuZn6zVhJmjEMvSRVztBLQ0TEooj434g4vuwfV/Z7m51MGhtDLw2RmbuArwIH3yG9AdiYmS83NpQ0Dn4zVhpGRMyh9SaZu2i9kWrZ0HdSS9PFtHrRvzRZMvO9iPgL4LvA7xt5TWfeupEObyWwB1ja9CDSeBh6aRgRsYzWT2o8E/iziFjQ8EjSmBl6aYiICFrfjL0hM38C3Ar8XbNTSWNn6KX/74vATzLz4bJ/B3BqRPxugzNJY+arbiSpcl7RS1LlDL0kVc7QS1LlDL0kVc7QS1LlDL0kVc7QS1Ll/g8Xa7KszKoUWQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# first the sample means\n",
    "\n",
    "X_sampled.groupby('X').mean().reset_index().plot(x='X', y='Yt', kind='bar')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "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>Yt</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>X</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-0.247299</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>100.466617</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>212.718946</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>323.913250</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>446.183699</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           Yt\n",
       "X            \n",
       "0   -0.247299\n",
       "1  100.466617\n",
       "2  212.718946\n",
       "3  323.913250\n",
       "4  446.183699"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_sampled.groupby('X').mean()[['Yt']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7f5c36b66320>"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEGCAYAAABrQF4qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAE4BJREFUeJzt3X+MXeV95/H3B2JwNqAAZha5HqdGjVcJoVuHzFJWbLQUlAaTBlOpTWCrYEWobhWiJZuKLQGtIFGRQGnLbiQCcmsSU5VQkjbCatl0WSCtsquQDNTldzYmMfVYBqaEEFBqwo/v/jHH8eCMPXfmzsydefx+SaN7zvM8597vXMFnjp97zn1SVUiS2nXEoAuQJM0vg16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUuDcNugCAE088sdasWTPoMiRpSXnggQf+uaqGphu3KIJ+zZo1jI6ODroMSVpSkjzVyzinbiSpcQa9JDXOoJekxi2KOfqpvPLKK4yNjbF3795Bl9Kz5cuXMzw8zLJlywZdiiT9VM9Bn+RIYBTYXVW/luRk4HZgBfAA8JGq+kmSo4FbgfcAzwEfrqqdMy1sbGyMY489ljVr1pBkpocvuKriueeeY2xsjJNPPnnQ5UjST81k6uYy4PFJ+9cDN1TV24HngUu69kuA57v2G7pxM7Z3715WrFixJEIeIAkrVqxYUv8CkXR46CnokwwDHwD+tNsPcDbwlW7IVuCCbntDt0/Xf05mmdZLJeT3WWr1Sjo89HpG/9+B/wq83u2vAH5YVa92+2PAqm57FbALoOt/oRv/Bkk2JRlNMjo+Pj7L8iVJ05l2jj7JrwHPVtUDSc6aqxeuqs3AZoCRkZFpF65dc8XfzNVLA7Dzug8csr+qeO9738tVV13F+vXrAfjyl7/Mli1bOP/88/nYxz42p/VIWjg3/u69gy6BS28+e8Feq5cz+jOB85PsZOLD17OB/wEcl2TfH4phYHe3vRtYDdD1v5WJD2WXlCTcfPPNfPKTn2Tv3r289NJLXHnlldx44418/vOfH3R5ktSzaYO+qj5VVcNVtQa4ELi3qn4LuA/4jW7YRuDObntbt0/Xf29VTXvGvhideuqpfPCDH+T666/nM5/5DBdffDFXXXUVTz75JOvWrePyyy8fdImSNK1+rqP/feD2JH8A/AOwpWvfAvxZkh3AD5j447BkXX311Zx22mkcddRRjI6OsmfPHh555BG2b98+6NIkqSczCvqq+jrw9W77e8DpU4zZC/zmHNS2KLzlLW/hwx/+MMcccwxHH330oMuRpBnzKxB6cMQRR3DEEb5VkpYm02uGjj32WF588cVBlyFJPVu033VzoOkuh1woK1as4Mwzz+TUU09l/fr1fPaznx10SZJ0SEsm6AfpmmuuecP+bbfdNphCJGkWnLqRpMYZ9JLUuEUd9EvtPqulVq+kw8OiDfrly5fz3HPPLZnw3Pd99MuXLx90KZL0Bov2w9jh4WHGxsZYSt9suW+FKUlaTBZt0C9btsyVmiRpDizaqRtJ0tww6CWpcQa9JDXOoJekxhn0ktS4aYM+yfIk30ryj0keTfLprv2LSb6fZHv3s65rT5LPJdmR5KEkp833LyFJOrheLq98GTi7ql5Ksgz4RpL/2fVdXlVfOWD8emBt9/PLwE3doyRpAHpZM7aq6qVud1n3c6jbVTcAt3bHfZOJRcRX9l+qJGk2epqjT3Jkku3As8DdVXV/13VtNz1zQ5J96+ytAnZNOnysa5MkDUBPQV9Vr1XVOmAYOD3JqcCngHcA/w44gYnFwnuWZFOS0SSjS+lrDiRpqZnRVTdV9UPgPuDcqtrTTc+8DHyB/QuF7wZWTzpsuGs78Lk2V9VIVY0MDQ3NrnpJ0rR6uepmKMlx3fabgfcBT+ybd08S4ALgke6QbcDF3dU3ZwAvVNWeealekjStXq66WQlsTXIkE38Y7qiqv05yb5IhIMB24He78XcB5wE7gB8DH537siVJvZo26KvqIeDdU7SffZDxBVzaf2mSpLngnbGS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMb1spTg8iTfSvKPSR5N8umu/eQk9yfZkeQvkhzVtR/d7e/o+tfM768gSTqUXs7oXwbOrqpfAtYB53ZrwV4P3FBVbweeBy7pxl8CPN+139CNkyQNyLRBXxNe6naXdT8FnA18pWvfysQC4QAbun26/nO6BcQlSQPQ0xx9kiOTbAeeBe4GngR+WFWvdkPGgFXd9ipgF0DX/wKwYorn3JRkNMno+Ph4f7+FJOmgpl0cHKCqXgPWJTkO+Crwjn5fuKo2A5sBRkZGqt/nk3Roj7/jnYMugXc+8figSzgszeiqm6r6IXAf8O+B45Ls+0MxDOzutncDqwG6/rcCz81JtZKkGevlqpuh7kyeJG8G3gc8zkTg/0Y3bCNwZ7e9rdun67+3qjxjl6QB6WXqZiWwNcmRTPxhuKOq/jrJY8DtSf4A+AdgSzd+C/BnSXYAPwAunIe6JUk9mjboq+oh4N1TtH8POH2K9r3Ab85JdZKkvnlnrCQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcb0sJbg6yX1JHkvyaJLLuvZrkuxOsr37OW/SMZ9KsiPJd5K8fz5/AUnSofWylOCrwO9V1YNJjgUeSHJ313dDVf3h5MFJTmFi+cB3AT8H/O8k/6aqXpvLwiVJvZn2jL6q9lTVg932i0wsDL7qEIdsAG6vqper6vvADqZYclCStDBmNEefZA0T68fe3zV9PMlDSW5JcnzXtgrYNemwMab4w5BkU5LRJKPj4+MzLlyS1Juegz7JMcBfAp+oqh8BNwG/AKwD9gB/NJMXrqrNVTVSVSNDQ0MzOVSSNAM9BX2SZUyE/J9X1V8BVNUzVfVaVb0O/An7p2d2A6snHT7ctUmSBqCXq24CbAEer6o/ntS+ctKwXwce6ba3ARcmOTrJycBa4FtzV7IkaSZ6uermTOAjwMNJtndtVwIXJVkHFLAT+B2Aqno0yR3AY0xcsXOpV9xI0uBMG/RV9Q0gU3TddYhjrgWu7aMuSdIc8c5YSWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGtfLClOrk9yX5LEkjya5rGs/IcndSb7bPR7ftSfJ55Ls6BYOP22+fwlJ0sH1ckb/KvB7VXUKcAZwaZJTgCuAe6pqLXBPtw+wnonlA9cCm5hYRFySNCDTBn1V7amqB7vtF4HHgVXABmBrN2wrcEG3vQG4tSZ8EzjugPVlJUkLqJc1Y38qyRrg3cD9wElVtafreho4qdteBeyadNhY17YHaYH94tZfHHQJPLzx4UGXoMNczx/GJjkG+EvgE1X1o8l9VVVMLBLesySbkowmGR0fH5/JoZKkGegp6JMsYyLk/7yq/qprfmbflEz3+GzXvhtYPenw4a7tDapqc1WNVNXI0NDQbOuXJE2jl6tuAmwBHq+qP57UtQ3Y2G1vBO6c1H5xd/XNGcALk6Z4JEkLrJc5+jOBjwAPJ9netV0JXAfckeQS4CngQ13fXcB5wA7gx8BH57RiSdKMTBv0VfUNIAfpPmeK8QVc2mddkqQ54p2xktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TG9bKU4C1Jnk3yyKS2a5LsTrK9+zlvUt+nkuxI8p0k75+vwiVJvenljP6LwLlTtN9QVeu6n7sAkpwCXAi8qzvm80mOnKtiJUkzN23QV9XfAz/o8fk2ALdX1ctV9X0m1o09vY/6JEl96meO/uNJHuqmdo7v2lYBuyaNGevafkaSTUlGk4yOj4/3UYYk6VBmG/Q3Ab8ArAP2AH800yeoqs1VNVJVI0NDQ7MsQ5I0nVkFfVU9U1WvVdXrwJ+wf3pmN7B60tDhrk2SNCCzCvokKyft/jqw74qcbcCFSY5OcjKwFvhWfyVKkvrxpukGJPkScBZwYpIx4GrgrCTrgAJ2Ar8DUFWPJrkDeAx4Fbi0ql6bn9IlSb2YNuir6qIpmrccYvy1wLX9FCVJmjveGStJjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJaty0QZ/kliTPJnlkUtsJSe5O8t3u8fiuPUk+l2RHkoeSnDafxUuSptfLGf0XgXMPaLsCuKeq1gL3dPsA65lYJ3YtsAm4aW7KlCTN1rRBX1V/D/zggOYNwNZueytwwaT2W2vCN4HjDlhIXJK0wGY7R39SVe3ptp8GTuq2VwG7Jo0b69p+RpJNSUaTjI6Pj8+yDEnSdPr+MLaqCqhZHLe5qkaqamRoaKjfMiRJB/GmWR73TJKVVbWnm5p5tmvfDayeNG64a9NCuuatg64Arnlh0BVI6sz2jH4bsLHb3gjcOan94u7qmzOAFyZN8UiSBmDaM/okXwLOAk5MMgZcDVwH3JHkEuAp4EPd8LuA84AdwI+Bj85DzZKkGZg26KvqooN0nTPF2AIu7bcoSdLc8c5YSWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjZrtmLABJdgIvAq8Br1bVSJITgL8A1gA7gQ9V1fP9lSlJmq25OKP/lapaV1Uj3f4VwD1VtRa4p9uXJA3IfEzdbAC2dttbgQvm4TUkST3qN+gL+F9JHkiyqWs7qar2dNtPAydNdWCSTUlGk4yOj4/3WYYk6WD6mqMH/kNV7U7yr4G7kzwxubOqKklNdWBVbQY2A4yMjEw5RpLUv77O6Ktqd/f4LPBV4HTgmSQrAbrHZ/stUpI0e7MO+iRvSXLsvm3gV4FHgG3Axm7YRuDOfouUJM1eP1M3JwFfTbLveW6rqq8l+TZwR5JLgKeAD/VfpiRptmYd9FX1PeCXpmh/Djinn6IkSXPHO2MlqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklq3LwFfZJzk3wnyY4kV8zX60iSDm1egj7JkcCNwHrgFOCiJKfMx2tJkg5tvs7oTwd2VNX3quonwO3Ahnl6LUnSIcxX0K8Cdk3aH+vaJEkLrJ/FwfuSZBOwCeBtb3tb38+35oq/6fs5+rXzug8MuoQJ17ww6AoWjYc3PjzoEhaNdz7x+KBLWDQuvfnsQZewoOYr6HcDqyftD3dtP1VVm4HNACMjI9XvCy6akJWkRWa+pm6+DaxNcnKSo4ALgW3z9FqSpEOYlzP6qno1yceBvwWOBG6pqkfn47UkSYc2b3P0VXUXcNd8Pb8kqTfeGStJjTPoJalxBr0kNc6gl6TGGfSS1LhU9X2vUv9FJOPAU4OuAzgR+OdBF7FI+F7s53uxn+/Ffovhvfj5qhqabtCiCPrFIsloVY0Muo7FwPdiP9+L/Xwv9ltK74VTN5LUOINekhpn0L/R5kEXsIj4Xuzne7Gf78V+S+a9cI5ekhrnGb0kNc6gl6TGGfSS1LiBLSU4aEnewcSC5fvWst0NbKsq11uTgCSnA1VV305yCnAu8ET3FeSHtSS3VtXFg66jV4flh7FJfh+4CLidiYXLYWK5wwuB26vqukHVpsHqTgBWAfdX1UuT2s+tqq8NrrKFleRqYD0TJ4N3A78M3Ae8D/jbqrp2gOUtqCQHro4X4FeAewGq6vwFL2qGDteg/3/Au6rqlQPajwIeraq1g6ls8Uny0ar6wqDrWAhJ/jNwKfA4sA64rKru7PoerKrTBlnfQkryMBPvwdHA08BwVf0oyZuZ+CP4bwda4AJK8iDwGPCnQDER9F9i4sSQqvq7wVXXm8N1jv514OemaF/Z9Wm/Tw+6gAX028B7quoC4CzgvyW5rOvLwKoajFer6rWq+jHwZFX9CKCq/oXD7/+REeAB4Crghar6OvAvVfV3SyHk4fCdo/8EcE+S7wK7ura3AW8HPj6wqgYkyUMH6wJOWshaBuyIfdM1VbUzyVnAV5L8PIdf0P8kyb/qgv49+xqTvJXDLOir6nXghiRf7h6fYYll52E5dQOQ5AjgdN74Yey3q+q1wVU1GN1/uO8Hnj+wC/i/VTXVv36ak+Re4JNVtX1S25uAW4DfqqojB1bcAktydFW9PEX7icDKqnp4AGUtCkk+AJxZVVcOupZeHbZBr/2SbAG+UFXfmKLvtqr6TwMoa8ElGWZiyuLpKfrOrKr/M4CypL4Z9JLUuMP1w1hJOmwY9JLUOINeOkCS1Um+n+SEbv/4bn/NYCuTZseglw5QVbuAm4B9d0hfB2yuqp0DK0rqgx/GSlNIsoyJm2RuYeJGqnUH3kktLRVL6qJ/aaFU1StJLge+BvyqIa+lzKkb6eDWA3uAUwddiNQPg16aQpJ1THxT4xnAf0mycsAlSbNm0EsHSBImPoz9RFX9E/BZ4A8HW5U0ewa99LN+G/inqrq72/888M4k/3GANUmz5lU3ktQ4z+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWrc/wdpKWInQic9owAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# and now the (usually unmeasured) population means.\n",
    "\n",
    "X.groupby('X').mean().reset_index().plot(x='X', y='Yt', kind='bar')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "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>Yt</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>X</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-2.025832</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>95.434589</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>201.557770</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>297.889978</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>398.047538</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           Yt\n",
       "X            \n",
       "0   -2.025832\n",
       "1   95.434589\n",
       "2  201.557770\n",
       "3  297.889978\n",
       "4  398.047538"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X.groupby('X').mean()[['Yt']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## We can see there's a good bit of bias in the expected values! The sample means, $E[Yt|D=d, S=1]$,  aren't good approximations to the population means, $E[Yt|D=d]$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now, let's use the estimator to get corrected expectations. \n",
    "We want to estimate\n",
    "\n",
    "$E_{population}[Yt|X] = \\sum_{Y}E[Yt|Y, X, S=1]P(Y|X)$, \n",
    "\n",
    "or, changing notation slightly, \n",
    "\n",
    "$E_{population}[Yt|X] = \\sum_{Y}E_{sample}[Yt|Y, X]P_{population}(Y|X)$\n",
    "\n",
    "But $Y$ is continuous! We actually need to take the integral over $Y$,\n",
    "\n",
    "$E_{population}[Yt|X] = \\int E_{sample}[Yt|Y, X]P_{population}(Y|X) dY$\n",
    "\n",
    "The problem is actually pretty tricky, since we really need to estimate $P_{population}(Y|X)$ as well. We've only done the binary case until now. Let's tackle the whole problem.\n",
    "\n",
    "We'll estimate $P_{population}(Y|X)$ with a kernel density estimator, and then actually perform the integral numerically using scipy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.00232888, 0.0009541 , 0.002253  , ..., 0.00250237, 0.00089457,\n",
       "       0.00260236])"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# estimate the density, $P_{population}(Y|D)$\n",
    "\n",
    "from statsmodels.nonparametric.kernel_density import KDEMultivariateConditional\n",
    "\n",
    "conditional_model = KDEMultivariateConditional(X['Y'], X['X'], 'c', 'o', 'cv_ml')\n",
    "conditional_model.pdf(X['Y'], X['X'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_y_given_x = conditional_model.pdf(X['Y'], [0.]*len(X))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Conditional Density Plot for P(Y|X=0)')"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAEWCAYAAABSaiGHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xu8XFV99/HPl5MLBBBCEikESEKJlxO1SiPVV62lJhak1ujroW0QMCqWkogPVvvUINZa1Fa8IW0NSAWN5MjlQR+NrUoJ0KqvyuVwEUgwGrkmBAgh4SJISPJ7/tjrwGQyM3ufy+y5fd+v17zOzF5r9vzmcuY3a+2111JEYGZm1mx7tDoAMzPrDU44ZmZWCiccMzMrhROOmZmVwgnHzMxK4YRjZmalcMKxppH0bkk/qbj9lKTDG9RfLenoJsc0U1JIGtfMxxkuSX8gaW1JjxWSjhijfR0o6UeSnpT0hbHYZ87jTZS0RtJBBevfK2lmc6Pa7TE/IOmcMh+zUzjh9CBJ75Q0mBLARkk/kPSGZj9uROwTEXenGL4u6VNV5XMi4r+aHUcj6QvqmfQFulXS/0g6TVJT/1ci4scR8dKqOOaPZF8VSfWpdLlX0tIR7GeXHwx1nAo8CrwoIj48knhrPOaOFPcTkm6T9Naqx/tRRGyU9D5Jd0maWHH/KZIekXRsjX3/qaSHJB1QsW2BpA2S9htmnDMlXSfpaUk/r3qv/g04UdKLh7PPXuCE02MkfQj4EvCPwIHAYcAyYEEr42ozfxoR+wIzgM8AHwEuam1II7J/ROwDnAB8vNaX8BiYAayJEZxB3qCV+dMU9/5kr/sVkianstOASwAi4qvABuDjFff9EvD9iPhh9U4j4nvAtcC56fH3B84HFkfE48MM/1LgVmAKcBZwpaRp6XF+A/wAeNcw99n9IsKXHrkA+wFPAX/WoM5Esn/aB9PlS8DEVHY0sB74MPAIsBF4T8V9pwArgSeAG4FPAj+pKA/gCLJfqc8B21I830vl9wLzxyCOPyH7MngCeAD4REXZzBTHuDrP//kYKrYdBewEXlER2+eB+4GHgQuAvQrGdhywBniS7Mvybyrvl65fkh7vmfT6/C3wH8AHquK6HXhHjeew23MEbqp4rACOqPhMfAPYBNwHfIzsh+jLgd8AO1IMW2s8zter3sf5Bd+3jwAPAZfU2Oe72fUzs3eKdy7Zj6Nnqp7XTGAL8GrgmPSYk6vez5kVt6em9+UY4GvApSP4P3oJ8Cywb8W2HwOnVdw+Ebiu1f/z7XZpeQC+lPhmw7HAdup82aY6ZwPXAy8GpgH/A3wylR2d7n82MD59eT499A8OXAZckb4kXpG+UHdLOOn614FPVT32vbyQcEYTx9HAK9MX56vIksLbU9lMhplw0vb7yX4JQ/YLeSVwALAv8D3gnwrGthH4g3R9MnBkxf3W14sD+HPghorbvwNsBibUiPX55wgI+P0Uw7wa78M3gO+m5zET+AVwSip7d+X7V+f12uV9LPi+nUOWmPaqsb/nHzPFfwZZct6P7IfE6hr3+QBwC3DP0Ptc9TrOrNp2Alk34CZgWlXZ7cDWOpdlqc47gLuq7vevwL9U3D4SeKzV//PtdnGXWm+ZAjwaEdsb1DkRODsiHomITcA/ACdXlD+Xyp+LiO+T/bJ9qaQ+4H8BH4+IX0fEncDyUcQ6ojgAIuK/IuKOiNgZEbeTdX/84ShigeyX8wGSRNZC++uIeCwiniTrnlxYJLZU1i/pRRGxJSJuKfj4K4GXSJqdbp8MXB4R2xrc51HgMeCrwNKIuKayML1nC4EzI+LJiLgX+AK7vs7Dlfe+7QT+PiKejYhn6uzjdZK2krWCTiBrxT1O1sX2ZI36/0r2ut4WEd8pEOP1ZAnsP1OMz4uIV0XE/nUuS1K1fYDqLrjHyZL2kKEkaRWccHrLZmBqzgitg8m6Vobcl7Y9v4+qhPU02T/gNLJfpA9U3XekRhoHkn4vHdDdJOlxsn7/qaOIBWA62Zf3NGAScHMaVLAV+GHanhsbWVI+DrhP0n9Len2RB4/suMDlwElpAMMJpGMZDUyNiMkR8fKI+Oda5WStsOrXeXqRmOrIe982pefSyPXpC35qRLwuIlal7VvY9UsdgMiaFHcBqwvGeCFZy+64oq9/laeAF1VtexG7JsN92T0p9TwnnN7yU7K+57c3qPMg2YHgIYelbXk2kXWXHFp133ryDjKPNA6Ab5K1CA6NiP3IjrGo4H13I+m1ZF/CPyFrNTwDzKn45btfZAe5c0XETRGxgKzL6TtkXZA1q9bYtpysBTEPeDoifjrMp1LtUbKWQfXrvKFBDHny3rfRTE9/OzBrNEPaJZ1C9hldAnwU+KqkCRXlqytG91VfLkjVVgOHS6pMfr/Drgnv5cDPRhpnt3LC6SGpW+LjwJclvV3SJEnjJb1F0mdTtUuBj0maJmlqqr+iwL53AN8GPpH22w8sanCXh4G65+SMNI5kX7L+899IOgp4Z8H77ULSi9KQ3MuAFUPddGTDXs8dGvYqabqkYwrsb4KkEyXtFxHPkQ1q2Fmn+m6vT0owO8m6vfJaN7nSe3YF8GlJ+0qaAXyIF17nh4FDKr+QCxjN+5YX73pgHdkgjmGTdDDwOeAvI+JZsh8im8lGmQ09xpzIhu/XupyW6vwCuA34e0l7SnoH2bHCb1U83B+SjVSzCk44PSYivkD2pfIxslbJA8DpZL+2AT4FDJL9mryD7GDsp3bfU02nk3UdPUR2MPlrDepeRHYsY6ukWv3uo4ljCXC2pCfJvvDqtSLq+V667wNkX0ZfBN5TUf4Rsi++6yU9AazihWM0eU4G7k33O42sxVLLP5F9cW+V9DcV279BNiBiTL7EyQ64/xq4m6wF903g4lR2Ldmv9ockPVpwf6N534r4CiM/xrQMuCwifgzPd8X9JfBBSXOGua+FZCPntpANnT9+6HiQpD3Juk1HcwyzKyl7zc2sE0h6F3BqRDT9RN12lE7yvJVsxN3GAvXvBY5OAyJKIekDZN25f1vWY3aKtprew8zqkzSJrPW2rNWxtErqCutvdRyNRMS/tDqGduUuNbMOkI4RbSI7rvLNFofTSb5Edg6NtQF3qZmZWSncwjEzs1L4GE6FqVOnxsyZM1sdhplZR7n55psfjYhpefWccCrMnDmTwcHBVodhZtZRJBWaVcRdamZmVgonHDMzK4UTjpmZlcIJx8zMSuGEY2ZmpfAoNbNEOQsY+Bxps9FxC8d62qRJWaLJSzZQvJ6Z1eYWjvWskSYPCcaPh22NFnc2s924hWM9abQtleeec2vHbLiccKznjGWimD9/7PZl1u2ccKynjHWr5JprxnZ/Zt3MCcd6RrO6wNy1ZlaME471hL6+5u5/0qTm7t+sGzjhWNcbGICdO4vXj3jhsnhxsfs888zIYjPrJU441vVOOql43eqTO5ctK37Cp7vWzBorNeFIOlbSWknrJC2tUT5R0uWp/AZJMyvKzkzb16b13ZF0qKTrJK2RtFrSGRX1PyFpg6Tb0uW4Mp6jtZfhJIFGicWzDJiNXmknfkrqA74MvBlYD9wkaWVErKmodgqwJSKOkLQQOAf4C0n9wEJgDnAwsErSS4DtwIcj4hZJ+wI3S7q6Yp/nRsTny3mG1snGKqFITk5m9ZTZwjkKWBcRd0fENuAyYEFVnQXA8nT9SmCeJKXtl0XEsxFxD7AOOCoiNkbELQAR8SRwFzC9hOdiHaBo66ZognAiMRudMhPOdOCBitvr2T05PF8nIrYDjwNTitw3db+9BrihYvPpkm6XdLGkyaN/CtZthptEVqzIr+NjOWa1dcWgAUn7AN8CPhgRT6TN5wO/Dbwa2Ah8oc59T5U0KGlw06ZNpcRrzVfkS//gg4e/3xNPHP59zCxTZsLZABxacfuQtK1mHUnjgP2AzY3uK2k8WbIZiIhvD1WIiIcjYkdE7AT+jaxLbzcRcWFEzI2IudOmTRvF07NOs6H601dQkVaRWzlmuysz4dwEzJY0S9IEskEAK6vqrAQWpevHA9dGRKTtC9MotlnAbODGdHznIuCuiPhi5Y4kHVRx8x3AnWP+jKwtFfmyL3p+jZmNndJGqUXEdkmnA1cBfcDFEbFa0tnAYESsJEsel0haBzxGlpRI9a4A1pCNTHt/ROyQ9AbgZOAOSbelh/poRHwf+KykVwMB3Av8VVnP1drfsmWju39EfmLr64MdO0b3OGbdROGhN8+bO3duDA4OtjoMG4X58/Mn1BzLIdB5/O9lvUDSzRExN69eVwwaMBtS5uzN/f3lPZZZN3DCsZ5SZFhzUatX59fx4AGzFzjhWNco8uU+1sOa99prbPdn1s2ccKxnNKML7Omn8+vMmTP2j2vWiZxwrCssWZJfp0gXWDOsWZNfx6wXOOFYVzj//MblezTxk+6RaGbFOOFYT2j1+TAePGDmhGNdYP78VkcA48e3OgKz9ueEYx0v79ybsRwKXc+2bc1/DLNO54RjXa9dZnh2t5r1Oicc62jt0J02ZCTLHZj1Eicc62h53Wllzgo90uUOzHqFE451tdHOCj3W3K1mvcwJxzpWO57B73V2zOpzwrGOlXcG/7x55cRRqd1aVGbtxAnHutaqVa2OoDZ3q1mvcsKxjjQw0OoI6vNUN2a1OeFYRzrppMblPpZi1n6ccKwrtfuxlL6+VkdgVj4nHLMmyJtOZ+fOcuIwaydOONZx8oZDt8MZ/+0ynY5ZO3HCsY6TNxy6U874nzCh1RGYlcsJx6xJ8karPfdcOXGYtQsnHDMzK4UTjnWUvPNvWjG7wGj4JFDrJU441lEWLWpc3m6zC/gkULMXOOFYR9mxo9URmNlIOeFYx8hbbK2/v5w4xtqkSa2OwKwcTjjWMfIWW1u9upw4hiuvW+2ZZ8qJw6zVSk04ko6VtFbSOklLa5RPlHR5Kr9B0syKsjPT9rWSjknbDpV0naQ1klZLOqOi/gGSrpb0y/R3chnP0czMaist4UjqA74MvAXoB06QVN0JcgqwJSKOAM4Fzkn37QcWAnOAY4FlaX/bgQ9HRD/wOuD9FftcClwTEbOBa9Jt61KdPjeZu9WsF5TZwjkKWBcRd0fENuAyYEFVnQXA8nT9SmCeJKXtl0XEsxFxD7AOOCoiNkbELQAR8SRwFzC9xr6WA29v0vOyEuRNZ7N8eePyVsubW83datYLykw404EHKm6v54XksFudiNgOPA5MKXLf1P32GuCGtOnAiNiYrj8EHFgrKEmnShqUNLhp06bhPSMrTd50Nu0+d1m7x2dWhq4YNCBpH+BbwAcj4onq8ogIoOah24i4MCLmRsTcadOmNTlSs/ryWnFmna7MhLMBOLTi9iFpW806ksYB+wGbG91X0niyZDMQEd+uqPOwpINSnYOAR8bsmZiNQF63Wl4rzqzTlZlwbgJmS5olaQLZIICVVXVWAkPnkh8PXJtaJyuBhWkU2yxgNnBjOr5zEXBXRHyxwb4WAd8d82dkbSHvi7xduFvNet24sh4oIrZLOh24CugDLo6I1ZLOBgYjYiVZ8rhE0jrgMbKkRKp3BbCGbGTa+yNih6Q3ACcDd0i6LT3URyPi+8BngCsknQLcB/x5Wc/VxtbknAHt3fRFLnk6HOteCn+6nzd37twYHBxsdRhWJW+Cy076CE+YkL8sQSc9HzMASTdHxNy8el0xaMB618SJrY5geLZta3UEZq3jhGNtbcmSxuUXXVROHGXKmzPOrFO5S62Cu9Tazx57NO5i6sSP75w5+SPSOvF5We9yl5p1hUZfvJ06O3S7TjJq1mxOONaxuvmLO29knlkncsKxtpV3/KaT7b9/4/KtW8uJw6xMTjjWtr7ylVZH0DxbtrQ6ArPyOeFY29q5s9URmNlYcsKxjrT33q2OYPTGj29cnnfCq1mnccKxtjS9euGKKt3Q3eaTQK3XOOFYW3rwwcbl3TR/mlmvcMIxa2NeI8e6iROOdZxOmz+tkcWLG5d7jRzrJk441nbyftV30/xpy5a1OgKz8jjhWNvJ+1Xfa8dv3K1m3cIJx6zFvPS09QonHOsoecc8OlGvtdisdznhWFvJWwvGxzzMOpcTjrWVa65pdQTtybMOWDdwwjFrA15wzXqBE451jHnzWh2BmY2GE461jbxFx1atKicOM2sOJxxrG72+6Fjektl9feXEYdYsw044kvaW5I++laobliPIk7dkttcHsk6Xm3Ak7SHpnZL+Q9IjwM+BjZLWSPqcpCOaH6Z1u4GBxuXdsByBWa8r0sK5Dvht4EzgtyLi0Ih4MfAG4HrgHEknNTFG6wFnnNG4vFdOjtwj5z8yb50gs3Y2rkCd+RHxXPXGiHgM+BbwLUk5axeaNbZ5c6sjaA87djQ+5yZvnSCzdpabcCLiOUkvAxYAQ7+vNgArI+KuoTrNC9F6XTdOZ2PWi4ocw/kIcBkg4MZ0EXCppKXDeTBJx0paK2ldrftKmijp8lR+g6SZFWVnpu1rJR1Tsf1iSY9IurNqX5+QtEHSbely3HBitfLkHb/xdDZm3UGRc4qzpF8Ac6pbMZImAKsjYnahB8pGtv0CeDOwHrgJOCEi1lTUWQK8KiJOk7QQeEdE/IWkfuBS4CjgYGAV8JKI2CHpjcBTwDci4hUV+/oE8FREfL5IfABz586NwcHBotVtjIwbl3Ul1bLHHvXLutWECfBcgz6D8eNh27by4jHLI+nmiJibV6/IoIGdZF/y1Q5KZUUdBayLiLsjYhtZq2lBVZ0FwPJ0/UpgniSl7ZdFxLMRcQ+wLu2PiPgR8Ngw4rA20yih9OJQ4Lxk0igZmbWzIoMGPghcI+mXwANp22HAEcDpw3is6RX3h6yV83v16kTEdkmPA1PS9uur7ltkvM7pkt4FDAIfjogtw4jX2sCMGa2OwMzGSpFBAz+U9BKyFkXloIGbIqKdOzvOBz4JRPr7BeC91ZUknQqcCnDYYYeVGZ8V8OlPtzoCMxsrRVo4RMROdm1hjMQG4NCK24ekbbXqrJc0DtgP2FzwvtUxPzx0XdK/Af9ep96FwIWQHcMp8kRs7CxZ0ri8V86/qTZ+fOOus+nTYUPD/wCz9lPmXGo3AbMlzUoDDhYCK6vqrAQWpevHA9dGNqphJbAwjWKbBcwmGy1Xl6SDKm6+A7izXl1rnfPPr1/Wy3OH5R3H8fk41omGlXAkvany73BExHayYz5XAXcBV0TEaklnS3pbqnYRMEXSOuBDwNJ039XAFcAa4IfA+4e68yRdCvwUeKmk9ZJOSfv6rKQ7JN0O/BHw18ON2Zorbzj08uWNy82ss+QOi96lsnRLRBw59LeJcbWEh0WXa+LExr/ke31RsrxVPnv99bH2MZbDomvuf4T3M3uezyVpbMWKxuVedto6jdfDsba0116tjqD1enXAhHUvJxxribzjN08/XU4cZlYeJxxrife9r9URmFnZhptwnkp/nxzrQKy3/OY3rY6gM8yb17h8zpxy4jAbC8NKOBHxxsq/ZtZcq1Y1Ll+zpnG5WTtxl5qVLu9Xude/MetOTjhWurxf5V7/xqw7FVmAbXmaisbMWsDn41i3KNLCeQD4aeXqmwCSXiXp4mYEZb2rv7/VEbQfn49j3aLI8gQfk3Q9sErSGcB4sjVy9gXOa3J81mXmz29cvnp1OXGYWfkKLU8A/Ihs0szvAY8Af55W2jQblmuuqV+2zz7lxdFp8pYrGBhwS8jaX5FjOMuAO8jOwXk5cC3wvyVNanJs1mMuuKDVEbSvvHnnTjqpnDjMRqPIMZyfAS+LiKURsTYi3km2HMD1aSVQs0LyutP8C92suxU5hvOVGtu+IOlW4PvAEc0IzLpPo+40M+t+uQlH0mF1itYB76ko3xoRT4xZZNZTpkxpdQRm1mxFBg0sB4Laa+AMbQ/g68A3xiwy6ypLljQuP8/jHXNFND7nRvKibNbehrXiZ7fzip/N49Urx4ZfR2tHY7bipySPf7Gm8nBos95QZJTayZLOk9TX9GisJ3k4dHHjxzcuz1vYzqyViiSctwDPANdKmtbkeKwL5R2/8XDo4nw+jnWy3IQTETsjYinZNDY/lnSqpKN84qcVdf75rY7AzNpBoeUJJL0VeB+wDTgS+DzwgKR1TYzNukBeF0/eipZm1j2KnIdzD7AGODcirq4qO6RZgVl3ePe7G5fnrWhpwzd9OmzY0OoozHZX5DycP4mImktmRcR6AEkKj6+2GrZvb3UE3SfvfJwHHywvFrPhKNKltkzSB6pnHJA0QdKbJC0HFjUnPOtked1pnl3ArLcUaeEcC7wXuFTS4cAWYE+gD/hP4EsRcWvzQrROddZZjcs9u4BZbxnWTAOSxgNTgWciYmvTomoRzzQwtnxWfPNMngxbG/wHLl4My5aVF4/1trGcaWCRpEclPQZ8FXiqG5ONja28c28OPricOLrVli2Nyz0U3dpRkWM4fwe8GXgZcD/wjyN9MEnHSloraZ2kpTXKJ0q6PJXfIGlmRdmZaftaScdUbL9Y0iOS7qza1wGSrpb0y/R38kjjtuHL+8LzKCqz3lMk4TwREbdGxCMR8XfAUSN5oDQ1zpfJZi7oB06Q1F9V7RRgS0QcAZwLnJPu2w8sBOaQHVNaVjHVztfTtmpLgWsiYjZwTbptJcgbLNDnSZLMelKRhHNQml3gjWlqm5zZnOo6ClgXEXdHxDbgMmBBVZ0FZMshAFwJzJOktP2yiHg2Iu4hW4vnKICI+BHwWI3Hq9zXcuDtI4zbhilvsMCpp5YTR7frr/65VmWy2/TWZoqMUvt74JXAienvPpK+T7b09O0RcWnBx5oOPFBxez3we/XqRMR2SY8DU9L266vuOz3n8Q6MiI3p+kPAgbUqSToVOBXgsMPqrTVnw3HffY3LfTB7bKxe3XhgRqNBBWatUGSJ6Qsrb6fZBV4JvAo4DiiacFomIkJSzTFR6fldCNkotVID60Lz5zcu33vvcuIws/ZTaC61ShGxPiJ+EBHnRMTJw7jrBuDQituHpG0160gaB+wHbC5432oPSzoo7esg4JFhxGojdM01jcu/8pVy4rBM3mhBszINO+GMwk3AbEmzJE0gGwSwsqrOSl6YteB44No0Zc5KYGEaxTYLmA3cmPN4lftaBHx3DJ6DjZKXIhhbK1Y0LvfwaGsnpSWciNgOnA5cBdwFXBERqyWdLeltqdpFwJQ0C/WHSCPLImI1cAXZJKI/BN4fETsAJF0K/BR4qaT1kk5J+/oM8GZJvwTmp9vWRHPmNC6fMaOcOHqJE7h1kmHNNNDtPNPA6OTNLLBihb8gm8EzOlirjdlMA2ZFFFna2MmmNfISkllZnHBsTOSde+OF1prHLRjrFE44Nibyzr3xQmtm5oRjo+ZlpNufZx2wduCEY6N2coOzsfr73bopQ97s2551wNqBE46Nypw5jY8hrF5dXiy9zLNvWydwwrFRWbOm1RFYUXnnSZk1mxOOjVjesRsPxy3X+Jx53P3jwFrNCcdG7K/+qnH5aaeVE4dltm1rdQRmjTnh2IgMDMCvf924jpchaD9udVorOeHYiLz3vY3Lp0wpJw7bVd6ibGat5IRjI5LXfXPeeeXEYbvyqEBrZ044NmxFFlnzvGnta3reWrlmTeKEY8MyMOBF1jrdgw+2OgLrVU44NiyLFjUu7+9366bVikzmWWR2b7Ox5oRjhc2fDzt2NK7jYwid4aSTWh2B9SInHCssrytt773LicPyLV7c6gjMdueEY4UU6YLxsZv2UeQcqAkTmh+HWSUnHCskb1aBffbxsZt2s9dejcufe66cOMyGOOFYriKzClxwQTmxWHFPP51fx+vkWJmccCxXo/VuwCPTOpnXybEyOeFYQ/Pn5w+z9ci09pXXrQb5J/KajRUnHGsob2TaHv4EtbUi3Wp577HZWPHXhdVVZGRa3mAC6ww+EdTKoChyWnKPmDt3bgwODrY6jLax777w1FP1y/ffH7ZsKS8eG7kiyxL4q8BGStLNETE3r55bOFbTkiWNkw042XQbr5VjzeaEYzWdf37j8nnzyonDxoZbL9YOnHBsN0uW5NdZtar5cVj53MqxZio14Ug6VtJaSeskLa1RPlHS5an8BkkzK8rOTNvXSjomb5+Svi7pHkm3pcurm/38ukVe62bGjHLisLHlVo61WmkJR1If8GXgLUA/cIKk6gVxTwG2RMQRwLnAOem+/cBCYA5wLLBMUl+Bff6fiHh1utzWxKfXNYqMVvr0p5sfh7WOWznWLGW2cI4C1kXE3RGxDbgMWFBVZwGwPF2/EpgnSWn7ZRHxbETcA6xL+yuyTxuGM85oXL54sWcV6GRFWzleFdSaocyEMx14oOL2+rStZp2I2A48DkxpcN+8fX5a0u2SzpU0sVZQkk6VNChpcNOmTcN/Vl1kYAA2b25cp8gsxNbeipys61VBrRm6edDAmcDLgNcCBwAfqVUpIi6MiLkRMXfatGllxtd2TjutcbnXWOkOeYvoDXHXmo21MhPOBuDQituHpG0160gaB+wHbG5w37r7jIiNkXkW+BpZ95vVkXfezbx5bt10k6I/Hpx0bCyVmXBuAmZLmiVpAtkggJVVdVYCi9L144FrI5sKYSWwMI1imwXMBm5stE9JB6W/At4O3NnUZ9fhGi2eNmWKh0F3G/94sFYYV9YDRcR2SacDVwF9wMURsVrS2cBgRKwELgIukbQOeIwsgZDqXQGsAbYD74+IHQC19pkeckDSNEDAbUBOh1HvGhiAnTvrl593XnmxWHkiirVgJA+ptrHhudQq9OpcalOnNh4s4I9I95ozB9asKVbXnwOrx3OpWWF5I9Ose3ktIyuTE06PyzvR07MKdL+iLRcPILDRcsLpcWed1bjcswr0hhUritVz0rHRcMLpYUuWwH331S/3rAK9Yzjvc5HJXc1qccLpUUuWNJ6kc8oUD53tNUW71vImdzWrxwmnRzX60pg0yUOhe5WP51gzOeHYbi680F1plm/OnFZHYJ3GCcd242TT24q2coqev2M2xAmnxwwMwMyZ9cvdVWLgrjVrjtKmtrHWGxiAU0+Fp5+uXydvxmjrHfvvD1u35tfz1DdWlFs4PeSss+onm76+bBi0R6bZkC1bitedPLl5cVj3cAunRzQ650blb3gVAAALRklEQVSC7dvLjcc6Q9EJPou0hMzcwukBeefcHHZYebFY5+nvL1bPx3MsjxNOD8g758bT11gjw5ng00nHGnHC6XJ5k3P6nBsrYjiDAqZPb14c1tmccLpc3uScTjZWVNGk8+CDzY3DOpcTTpe7//5WR2Dd5OCDi9Vz15rV4oTThYZO7txjj+xSz+LFpYVkXWLDhuJ1nXSsmodFd5nqkzt37Khdb948n3NjI1N0qDT4pFDblVs4XWRgAN71rtond/b1Zf/8M2Zki22tWlV+fNY9hpNE3NKxIW7hdImBAXjPe2DnztrlO3fWLzMbicWLi6+NM5R03NrpbW7hdIGBAVi0CJ57rn4dn9xpY20kXbJu7fQ2J5wON3TMpt6xmiE+udOaYSQtFsmJp1c54XSooZFoJ53UePZnyJaL9vk21iwj7SYbSjxOPr3Dx3A6zMAAnHEGbN5crP6ECV4u2ppvOCPXaqm8r4/zdC+3cDrIUPdZ0WQzZQpcfLFbN1aOiOInhjZS2fJx66e7OOG0ocoTN2fOfGE+tEbr2VSaNCkb+vzoo042Vq4NG8a+hVKdgJyEOpe71NpM9Ymb992X3YZi09TMmJENEHCisVYabRdbnlr7dldc+3MLZ5TqtUZGqlYr5umns+2NhjYPtWruvdfJxtpDRLlJoFZLqPIyZ055sVhtpSYcScdKWitpnaSlNconSro8ld8gaWZF2Zlp+1pJx+TtU9KstI91aZ8Txvr5DLVG7rsv+8caao2MJunUa8Xcf3/Wcpk0afeyKVO8zIC1r4j2mLdvzZr8pNQul/nzYerU+uV77pn9yJ06NbtU/+Ct90N4yRIYNy7bx7hx2e2x/tHcUESUcgH6gF8BhwMTgJ8B/VV1lgAXpOsLgcvT9f5UfyIwK+2nr9E+gSuAhen6BcDivBh/93d/N4Zjxoyh33C7XmbMGNZuhrXPFSuy61L2d8WKkT+WWSvU+nz7MjaXSZMiFi/O/lZvnzev9n36+navO9zvFWAwcr5f01tfWsJ5PXBVxe0zgTOr6lwFvD5dHwc8Cqi67lC9evtM93kUGFfrsetdhptwpNpvoDSs3exixYraHxYnFutGBx/c+i/pbrtUJ5CRXIb7o7lowimzS2068EDF7fVpW806EbEdeByY0uC+9bZPAbamfdR7LAAknSppUNLgpk2bhvWE6h1TGc00MieemHWPzZjxwmSb7i6zbjU0qq3yYqOTN+tIEc1aR6vnBw1ExIURMTci5k6bNm1Y9611TGXSpNFPI3PiidnB/507PQjAeo8T0Oj09Y1+H82ae7HMhLMBOLTi9iFpW806ksYB+wGbG9y33vbNwP5pH/Uea9TcGjFrvlqdPv39rY6qPU2alA1cqvVDeN682vepTlBj8aO5njITzk3A7DR6bALZoICVVXVWAovS9eOBa1P/4EpgYRrFNguYDdxYb5/pPtelfZD2+d1mPCm3RszKt3p1/pGIRqvdtrt587LRp/VMnJj9yJ0yJbtU/uBdtqz2D+FVq7LRgkMJpq8vu718eXk/mhUltlklHQd8iWx02cUR8WlJZ5MdcFopaU/gEuA1wGNko8zuTvc9C3gvsB34YET8oN4+0/bDgcuAA4BbgZMi4tlG8c2dOzcGBwfH+mmbmXU1STdHxNzcemUmnHbnhGNmNnxFE04HNzrNzKyTOOGYmVkpnHDMzKwUTjhmZlYKDxqoIGkT8GuyaXE6xVQcbzM53ubppFjB8TYyIyJyz5x3wqkiabDIaIt24Xiby/E2TyfFCo53LLhLzczMSuGEY2ZmpXDC2d2FrQ5gmBxvczne5umkWMHxjpqP4ZiZWSncwjEzs1I44ZiZWSl6PuFI+oCkn0taLemzFdvPlLRO0lpJx1RsPzZtWydpaYti/rCkkDQ13Zakf04x3S7pyIq6iyT9Ml0W1d9rU+L8XHptb5f0/yTtX1HWtq9vO8VRSdKhkq6TtCZ9Xs9I2w+QdHV6j6+WNDltr/u5KDnuPkm3Svr3dHuWpBtSXJenpUVIy49cnrbfIGlmC2LdX9KV6XN7l6TXt+vrK+mv0+fgTkmXStqznV9bgNw1qLv5AvwRsAqYmG6/OP3tB34GTARmAb8iW/6gL10/HJiQ6vSXHPOhwFXAfcDUtO044AeAgNcBN6TtBwB3p7+T0/XJJcb6x8C4dP0c4Jx2f31TfG0RR424DgKOTNf3BX6RXsvPAkvT9qUVr3PNz0UL4v4Q8E3g39PtK8iWHgG4AFicri8BLkjXFwKXtyDW5cD70vUJwP7t+PoC04F7gL0qXtN3t/NrGxE938JZDHwm0jo5EfFI2r4AuCwino2Ie4B1wFHpsi4i7o6IbWTr7SwoOeZzgb8FKkd7LAC+EZnryVY7PQg4Brg6Ih6LiC3A1cCxZQUaEf8ZEdvTzevJVl4dirddX1/aKI5dRMTGiLglXX8SuIvsi2cB2Rcl6e/b0/V6n4vSSDoE+BPgq+m2gDcBV9aJd+h5XAnMS/XLinU/4I3ARQARsS0ittK+r+84YC9lKxtPAjbSpq/tkF5POC8B/iA1Mf9b0mvT9unAAxX11qdt9baXQtICYENE/KyqqC3jrfJesl+D0P7xtkscdaUukdcANwAHRsTGVPQQcGC63g7P40tkP5B2pttTgK0VP0QqY3o+3lT+eKpfllnAJuBrqQvwq5L2pg1f34jYAHweuJ8s0TwO3Ez7vrZAliG7mqRVwG/VKDqL7PkfQNYcfi1whbKVQlsmJ96PknVTtY1G8UbEd1Ods8hWah0oM7ZuJWkf4FtkK98+UflDNSJCUluc6yDprcAjEXGzpKNbHU8B44AjgQ9ExA2SziPrQnteu7y+6TjSArIkuRX4v5TYezFSXZ9wImJ+vTJJi4FvR9axeaOknWQT3m0gO1Yy5JC0jQbbmxqvpFeSfbh+lr5gDgFukXRUg3g3AEdXbf+vMuIdIundwFuBeel1hha+vgU1iq+lJI0nSzYDEfHttPlhSQdFxMbUpTPUNdzq5/H7wNuULQO/J/Ai4Dyyrqdx6Zd2ZUxD8a5P3UT7AZtLjHc9sD4ibki3ryRLOO34+s4H7omITQCSvk32erfrawu4S+07ZAMHkPQSsoOEjwIrgYVpZMcsYDZwI3ATMDuNBJlAdvBtZRmBRsQdEfHiiJgZETPJ/jmOjIiHUgzvSqNmXgc8nroArgL+WNLk9Ivoj9O2Ukg6lqw75W0R8XRFUdu9vlXaJY5dpD73i4C7IuKLFUUrgaERiIuA71Zsr/W5KEVEnBkRh6TP60Lg2og4EbgOOL5OvEPP4/hUv7TWRPpfekDSS9OmecAa2vP1vR94naRJ6XMxFGtbvrbPa8VIhXa5kCWYFcCdwC3AmyrKziIbqbQWeEvF9uPIRgf9iqzbqFWx38sLo9QEfDnFdAcwt6Lee8kOyq8D3lNyjOvI+o1vS5cLOuj1bYs4qmJ6A9lgkdsrXtPjyPrirwF+STbq8oC8z0ULYj+aF0apHU72A2MdWVfQ0CjRPdPtdan88BbE+WpgML3G3yEb3dmWry/wD8DP0/fXJWSjPtv2tY0IT21jZmbl6PUuNTMzK4kTjpmZlcIJx8zMSuGEY2ZmpXDCMTOzUjjhmLWpdH7HTyS9pWLbn0n6YSvjMhspD4s2a2OSXkF2/sRryGYGuRU4NiJ+1dLAzEbACceszSlbp+nXwN7AkxHxyRaHZDYiTjhmbS7NWHwLsI3sbPZnWxyS2Yh0/eSdZp0uIn4t6XLgKScb62QeNGDWGXbywpoyZh3JCcfMzErhhGNmZqXwoAEzMyuFWzhmZlYKJxwzMyuFE46ZmZXCCcfMzErhhGNmZqVwwjEzs1I44ZiZWSn+P5sexAct2hUwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as pp\n",
    "pp.plot(X['Y'], p_y_given_x, 'bo'); pp.xlabel('Y'); pp.ylabel('$P(Y|X=0)$'); pp.title(\"Conditional Density Plot for P(Y|X=0)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Conditional Density Plot for P(Y|X=4)')"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAEWCAYAAABSaiGHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X28XFV97/HPlzyRACWQRIoBSSzxIbG9SI9cfdW21GCJaA19XdoGg6JiwSAUr95bg7GtF7UVqyLtNQgKigQFLrUSW5USQK2vSuDwIJjQ1MiDJAYID+EpQEjO7/6x1yGTyZzZ+5zM7Nkz832/XvM6M3utvec3e+bMb/baa6+liMDMzKzd9up0AGZm1h+ccMzMrBROOGZmVgonHDMzK4UTjpmZlcIJx8zMSuGEY20j6d2Sflzz+GlJL29Sf42ko9sc0yxJIWl8O59ntCT9rqR1JT1XSDq8Rds6SNKPJD0l6XOt2GbO802StFbSwQXr3ydpVnuj2u05z5R0bpnP2S2ccPqQpHdIGkwJYJOk70l6Y7ufNyL2jYh7Ugxfk/TJuvJ5EfGDdsfRTPqCejZ9gW6R9B+S3i+prf8rEfHvEfHKujiOGcu2apLq0+l2n6SlY9jOLj8YRnAq8AjwaxHx4bHE2+A5d6S4n5R0h6S31T3fjyJik6T3Sbpb0qSa9adJeljSggbb/iNJD0o6sGbZQkkbJe0/xnh/P+3r2s/yl4HFkl4ylm32MiecPiPpQ8AXgL8FDgJeBiwHFnYyror5o4jYDzgM+DTwEeDizoY0JlMjYl/gROCvG30Jt8BhwNoYwxXkTY4yf5Linkq236+SdEAqez9wGUBEfAXYCPx1zbpfAL4bEd+v32hEfAe4ATgvPf9U4AJgSUQ8MYb4JwDnA6vrnuc54HvAu0a7zZ4XEb71yQ3YH3ga+JMmdSaR/dP+Kt2+AExKZUcDG4APAw8Dm4D31Kw7DVgJPAncDHwC+HFNeQCHk/1KfQHYluL5Tiq/DzimBXG8Fbg9xfEA8PGaslkpjvEjvP4XY6hZdhQwBLymJrbPAr8EHgK+BEwuGNtxwFrgKbIvy/9Vu166f1l6vmfT/vlL4F+BM+viuhP44wavYbfXCNxS81wBHF7zmfg6sBm4H/gY2Q/RVwPPATtSDFsaPM/X6t7HYwq+bx8BHgQua7DNd7PrZ2afFO8A2Y+jZ+te1yzgceAI4Nj0nAfUvZ+zah5PT+/LscBXgW/uwf/TUuAzaT98sq5sMXBjp//nq3breAC+lfhmwwJgOyN82aY65wA3AS8BZgD/AXwilR2d1j8HmJC+PLcO/4MDVwBXpS+J16Qv1N0STrrf6J/0PnYmnD2J42jgN9MX52+RJYXjU9ksRplw0vJfkv0ShuwX8krgQGA/4DvA3xWMbRPwu+n+AcCRNettGCkO4E+B1TWP/xvwKDCxQawvvkZAwO+kGOY3eB++DlyTXscs4L+AU1LZu2vfvxH21y7vY8H37VyyxDS5wfZefM4U/1lkyXl/sh8SaxqscyZwG3Dv8Ptctx9n1S07kawZcDMwo67sTmDLCLflNfUOS/tq3/p9kMqPBB7r9P981W4dD8C3Et/s7FfXgzl1fgEcV/P4WOC+dP9odv+F+TDwemAc2a/dV9WU/S1jTzhjimOE1/QF4Lx0fxZjSzg3AcvIvsCfAX6jpuwNwL1FYiNLXKeRnfOo3f7RNE84e5P9kp+THn+29guwblvDr3FLWudu4C/q34f0nm0D5taUnQb8IN1/N6NPOHnv2zZg7ybbezdZUtpClhRuqvlMLAZuarCOyJq1/nmE93NW3bLZZJ/Vy/fgf+ka4M+afJbnADvGuv1evfkcTn95FJie00PrpWRNK8PuT8te3EZEbK95vJXsV94Msl+kD9StO1ZjjQNJ/13SjZI2S3qCrN1/+h7EAjATeIzsdU4Bbk2dCrYA30/Lc2MD/gfZUc/9kn4o6Q1Fnjyy8wJXAielDgwnks5lNDE9Ig6IiFdHxD80Kic7CqvfzzOLxDSCvPdtc3otzdwUEVMjYnpEvD4iVqXlj5Mdie0ism/4u4E1BWO8iOzI7rii+7+WpD8C9ouIK5tU2w8Y9XmhXueE019+AjwPHN+kzq/ImguGvSwty7OZ7JfpoXXrjiRytjfWOAC+QdbkdWhE7E92jkUF192NpNeRfQn/mOxX97PAvPSlODUi9o/sJHeuiLglIhaSNTl9m6wJsmHVBssuJfuVPx/YGhE/GeVLqfcI2S/9+v28sUkMefLet7Fsc9idwOw96dIu6RSyz+jpwEeBr0iaWFO+pqZ3X/3tS6nafGAg9Xh7EPgz4IOSrql5qlcDPx1rnL3KCaePRNYT56+BL0o6XtIUSRMkvUXSZ1K1bwIfkzRD0vRUf0WBbe8AvgV8PG13LnByk1UeAka8JmescST7kbWfPyfpKOAdBdfbhaRfS11yrwBWRMRdETFE1u31vOFur5JmSjq2wPYmSlosaf+IeIGsU8PQCNV32z8pwQwBnyP/6CZXes+uAj4laT9JhwEfYud+fgg4pPYLuYA9ed/y4t0ArCfrxDFqkl4K/D3w5xHxPNkPkUfJmkqHn2NeZN33G93en6r9FfAKso4KR5D9uPky8J6ap/t9sp5qVsMJp89ExOfIvlQ+RnZU8gBwBtmvbYBPAoNkvybvIjsZ+8ndt9TQGWRNRw+StWt/tUndi4G5qVnq2w3K9ySO04FzJD1F9oU30lHESL6T1n2A7Mvo8+z6ZfIRsi++myQ9CawCXrnbVhp7J3BfWu/9ZEcsjfwd2Rf3Fkn/q2b518k6RLTkS5zshPszwD1kR3DfAC5JZTeQNVM9KOmRgtvbk/etiAvJ9uFYLAeuiIh/hxeb4v6c7OhkXtGNRMRTEfHg8I3siPeZiHgMQNLeZM2ml44xzp6ldILLzLqApHcBp0ZE2y/UraJ0keftZD3uNhWofx9wdETc1+bQap/zTLLm3L8s6zm7RaWG9zCzkUmaQnb0trzTsXRKagqb2+k4momIf+x0DFXlJjWzLpDOEW0mO6/yjQ6H002+QNbF2irATWpmZlYKH+GYmVkpfA6nxvTp02PWrFmdDsPMrKvceuutj0TEjLx6Tjg1Zs2axeDgYKfDMDPrKpIKjSriJjUzMyuFE46ZmZXCCcfMzErhhGNmZqVwwjEzs1I44ZhZy51+Okgj36w/uVu0mbVUkYQiwV57wY4d7Y/HqsNHOGbWMqM5ehka8tFOv3HCMbOWGGvycNLpH044ZrbH9jRpOOn0ByccM6sEJ53e54RjZnuklYliXuGJnq0bOeGY2Zhdfnlrt7d2bWu3Z9XihGNmY3bSScXqjWaex3HjxhaLVZ8Tjpm11XCyKZp0hoZaf+Rk1eCEY2ZjUuTcTX2SKZp0ih45WXdxwjGzUhVNOgcc0N44rHxOOGY2akWSQbPEMmFC/vpbthSPx7qDE46ZjdqeJoNt24rV81FOb3HCMbOWW7Eiv06RpjUf5fSWUhOOpAWS1klaL2lpg/JJkq5M5aslzaopOzstXyfp2LTsUEk3SloraY2ks2rqf1zSRkl3pNtxZbxGs143ZUp+ncWLi21r8uT8OhMnFtuWVV9p0xNIGgd8EXgzsAG4RdLKiKi91OsU4PGIOFzSIuBc4M8kzQUWAfOAlwKrJL0C2A58OCJuk7QfcKuk62q2eV5EfLacV2jWH559tnXb2ro1v7fbCy+07vmss8o8wjkKWB8R90TENuAKYGFdnYXApen+1cB8SUrLr4iI5yPiXmA9cFREbIqI2wAi4ingbmBmCa/FzEYwmos8AebOza/jcdZ6Q5kJZybwQM3jDeyeHF6sExHbgSeAaUXWTc1vrwVW1yw+Q9Kdki6R1PD0o6RTJQ1KGty8efNoX5NZXznmmNZvc82a1m/TqqknOg1I2hf4J+CDEfFkWnwB8BvAEcAm4HON1o2IiyJiICIGZsyYUUq8Zt3q+uubl7/0pWPbbpFOBjPddtH1ykw4G4FDax4fkpY1rCNpPLA/8GizdSVNIEs2l0fEt4YrRMRDEbEjIoaAL5M16ZlZG22s/48uqEgng1/9amzbtuooM+HcAsyRNFvSRLJOACvr6qwETk73TwBuiIhIyxelXmyzgTnAzen8zsXA3RHx+doNSTq45uEfAz9r+Ssys5Ypcu7H53K6W2m91CJiu6QzgGuBccAlEbFG0jnAYESsJEsel0laDzxGlpRI9a4C1pL1TPtAROyQ9EbgncBdku5IT/XRiPgu8BlJRwAB3AecVtZrNetFbtKyPaUYbZeSHjYwMBCDg4OdDsOskvKOLlr1VZL3PJMnZ92prTok3RoRA3n1eqLTgJn1j1ZeB2TlcsIxs1xlzk9T5EjJzXvdyQnHzHLlzU8zf345cQxzj7Xu5IRjZnts1arWbs/X5fQmJxwzqxxfl9ObnHDMrKl585qXt6s5rchRjufL6S7uFl3D3aLNdldWd+ixPHe7n9+KcbdoM+t6RebLaceAotYeTjhmVllFLvDMG1DUqsMJx8zGbMKE3ngOK4cTjpmNKO8cyrZt7Y+hyHN4UM/u4IRjZpU31nl2rFqccMys8orMszNxYvvjsD3jhGNmYzJ3bqcj2NULL3Q6AsvjhGNmDeUdMaxZU04cwzyoZ/dzwjGzhrrxiMHD3VSbE46ZdQ0P6tndnHDMbNQ6dW2MB/Xsbk44ZrabKlx/M5IinRU83E01OeGYWVcp0lnBw91UkxOOmXWdqVPz63jqgupxwjGzUdmrAt8ajz+eX2fLlvbHYaNTgY+OmVVJ3oRrO3aUE0eeIudyTj+9/XFYcU44ZraLtWs7HUExRc7lXHBB++Ow4pxwzKxrFRnU8/LL2x+HFeOEY2Zdq8ignied1P44rBgnHDN7Ud71K0WmfC6bz+V0D0WREfH6xMDAQAwODnY6DLOOybvgs6pfF0UmYKtq7L1A0q0RMZBXr9QjHEkLJK2TtF7S0gblkyRdmcpXS5pVU3Z2Wr5O0rFp2aGSbpS0VtIaSWfV1D9Q0nWSfp7+ule+WY+aP7/TEVgRpSUcSeOALwJvAeYCJ0qqPxg+BXg8Ig4HzgPOTevOBRYB84AFwPK0ve3AhyNiLvB64AM121wKXB8Rc4Dr02Mz60GrVuXXGTeu/XFYc2Ue4RwFrI+IeyJiG3AFsLCuzkLg0nT/amC+JKXlV0TE8xFxL7AeOCoiNkXEbQAR8RRwNzCzwbYuBY5v0+sy6wtVuOCzmbyRpIeGyonDRlbmR2gm8EDN4w3sTA671YmI7cATwLQi66bmt9cCq9OigyJiU7r/IHBQo6AknSppUNLg5s2bR/eKzHpIt1zwOZIiI0l7uJvOqvhvlmIk7Qv8E/DBiHiyvjyynhENTxlGxEURMRARAzNmzGhzpGbV1S0XfDaT12PNw910VpkJZyNwaM3jQ9KyhnUkjQf2Bx5ttq6kCWTJ5vKI+FZNnYckHZzqHAw83LJXYmaVVGT0AU/Q1jllJpxbgDmSZkuaSNYJYGVdnZXAyen+CcAN6ehkJbAo9WKbDcwBbk7ndy4G7o6IzzfZ1snANS1/RWZWOXnXCnmCts4pLeGkczJnANeSndy/KiLWSDpH0ttTtYuBaZLWAx8i9SyLiDXAVcBa4PvAByJiB/A7wDuBN0m6I92OS9v6NPBmST8HjkmPzayBvAs+u6nb8dat+XU83E1n+MLPGr7w0/pVt17wOZK81yO511orVfLCTzOzMuR1ke62BNornHDMrOcU6SI9cWL747BdOeGY9bm8gS2LTOdcRUuWNC9/4YVy4rCdfA6nhs/hWD/qtfM3tfJe24QJsG1bObH0Mp/DMbO+l9e7zkc55XLCMbOeVWRQT3eRLo8TjpmNqOoDdhaRd5Rz2mnlxGFOOGZ9La+nVtUH7Cwi7yjnmWfKicOccMz6ms9hZDyKdDmccMys5+WNr+ZRpMvhhGNmPa/I+GrWfk44Zn2ql6+/aeSlL21e7mkL2s8Jx8z6wsb62bfqeNqC9nPCMTOzUjjhmPWhvIsdx40rJ46y5TUjurdaeznhmPWhk05qXn7ppeXEUbb3v795uXurtZcTjpntpsjw/t1o+fJOR9DfnHDMrK/Mndu83GOrtY8Tjlmfyev+O2FCOXF0ypo1zcvzmt1s7EadcCTtI6lHTyma9b687r/9Pj/M0093OoLelZtwJO0l6R2S/lXSw8B/ApskrZX095IOb3+YZmatk9cLL28WVBubIkc4NwK/AZwN/HpEHBoRLwHeCNwEnCspp8+LmVl15PXCu+CCcuLoN7lTTEuaEBFNx5QtUqcbeIpp63UzZ+Y3qfXakDYj6behfdqpZVNMN0okkr6eV8fMqsfJxjppfF4FSSvrFwF/IGkqQES8vR2BmZm105IlzZvOLr+8d69H6pQi53AOAZ4EPg98Lt2eqrlvZl3AJ8J3lXcRqKeebr0iCWcAuBVYBjwRET8Ano2IH0bED9sZnJm1Tt6J8Lzh+3tRs/M4nnq69YqcwxmKiPOA9wDLJP1fCjTFNSJpgaR1ktZLWtqgfJKkK1P5akmzasrOTsvXSTq2Zvklkh6W9LO6bX1c0kZJd6TbcWOJ2axf5A3f34suu6x5+THHlBNHvyh84WdEbIiIPwG+B6wY7ROli0W/CLwFmAucKKl+kIlTgMcj4nDgPODctO5cYBEwD1gALK+5+PRraVkj50XEEen23dHGbGa9Le8czfXXlxNHvxj1SAMR8a8R8dExPNdRwPqIuCcitgFXAAvr6iwEhnvIXw3Ml6S0/IqIeD4i7gXWp+0RET8CHhtDPGZ9I68LcK8PZ2PVUOZYajOBB2oeb0jLGtaJiO3AE8C0gus2coakO1OzW8OZLiSdKmlQ0uDmzZuLvRKzHtPvw9lYOUaVcCS9qfZvxV1ANkLCEcAmRuhRFxEXRcRARAzMmDGjzPjMrALmz29e7t59rTPaI5zP1v0djY3AoTWPD0nLGtaRNB7YH3i04Lq7iIiHImJHRAwBXyY1wZn1m16dvbNVVq1qXu5hblpnrE1qOS3CDd0CzJE0W9JEsk4A9ReVrgROTvdPAG6IbOydlcCi1IttNjAHuLlpgNLBNQ//GPjZSHXNetnQUPPyJUvKicNsTN2bxyIitks6A7gWGAdcEhFrJJ0DDEbESuBi4DJJ68k6AixK666RdBWwFtgOfCAidgBI+iZwNDBd0gbgbyLiYuAzko4AArgP8GVcZg14FkyYNg0efbTTUfS+3ME7d6ks3RYRR0q6PSJe28a4OsKDd1qvmTcP1q5tXsfjp2XD2JzUZMz7JUucmJtp2eCdZta98pLN5MnlxFF1edfjXHRROXH0utEmnOG58J5qdSBmVr6tWzsdQXfYsaPTEfSGUSWciPi92r9mVl3uzjs6Eyd2OoLe5yY1sx6V1513bv3AUn3ukkualzuB77lRdRrode40YL3EM1qOnvfZ2LSs04CkS9N1M2bWJfxr3KqoSJPaA8BPaqcKAJD0W5JyDkLNrBN8dfzY9OOcQGUqMh/Ox4C/AVZJequk4yX9APgq8IP2hmdm7eCmocby5gS6/PJy4uhVRUca+BHwfeA7wMPAn6ZpAcysYvLOQ9jYLVuWf82OjazIOZzlwF1k1+C8GrgB+AtJU9ocm5m1gZuNxu7++zsdQXcrcg7np8CrImJpRKyLiHcAPwFukvSK9oZnZq3Wj1NJj4YHM22fIudwLoyIZ+uWfQ74IOBpm80qZGaRaQmtqbwx03weZ+xyr8OR9LImxYcBwweZWyLiyVYF1gm+Dse6XZHzN+4wkK/Zfpw2DR55pLxYukHR63CKdBq4lGyI/0ZvwfDyAL4GfH0UMZpZyZxs9pynMRi73IQTEX9QRiBmtmfcO6115s+H66/vdBS9p0gvtSazRJhZt/DYacXlTTvt8zhjU6SX2jslnS/JM6ObVdSUAhcprFnT/jj6xbJlnY6gOxVJOG8BngVukDSjzfGY2Rg8+2x+HWsdX48zNkW6RQ9FxFLgfODfJZ0q6Shf+GnWPebP73QE3cdNkK1XaD4cSW8D3gdsA44EPgs8IGl9G2MzswLGFWjszjsnYbvLa4L0eZzRy+2lJuleYC1wXkRcV1d2SLsCM7NihoY6HUF/et/7PK7aaBW5DuetEbG2UUFEbACQpPBMbmaV5P/MsdtnH3jmmcZlzz1Xbiy9oEiT2nJJZ9aPOCBpoqQ3SboUOLk94ZlZM772pr0uvLB5uSe6G50iQ9vsDbwXWAy8HHgc2BsYB/wbsDwibm9znKXw0DbWbfISztSp8Pjj5cTSq5rtY8lNmtDCoW0i4jlgOdmRzgRgOvBsRGzZ8zDNbKyK/Lp2stlze+89cvOZmytHp8hIAydLekTSY8BXgKedbMw6z9NIl+MrX+l0BL2jyDmcvwLeDLwK+CXwt22NyMxawvO6tEZeTzR3jy6uSMJ5MiJuj4iHI+KvgKPaHZSZNXfAAfl18uZ1sdY466xOR9A9iiScg9PoAr+XhraZMNYnk7RA0jpJ6yUtbVA+SdKVqXy1pFk1ZWen5eskHVuz/BJJD0v6Wd22DpR0naSfp78F/kXNusOWnEbtCWP+L7VGmo064OkKiiuScP4G+E3gE8A64DWSvivp7ySdWPSJ0uCfXyQbm20ucKKk+rfxFODxiDgcOA84N607F1gEzAMWkHVgGL6++mtpWb2lwPURMQe4Pj026wvbtnU6gt7igU9bo8hYahdFxJkR8fsRcSBZ1+h/BLYAx43iuY4C1kfEPRGxDbgCWFhXZyHZhG8AVwPzJSktvyIino+Ie4H1aXtExI+Axxo8X+22LgWOH0WsZpXla286Y9q0kct8HqeYQmOp1YqIDRHxvYg4NyLeOYpVZwIP1DzekJY1rBMR24EngGkF1613UERsSvcfBA5qVCk1Fw5KGty8eXOR12FWaXuN+r/aijj//JHLTjutvDi6WV98NNOwOw17zKcjuIGIGJgxw7MvWLUVufZmx472x9GPmvVWG2n4G9tVmQlnI3BozeND0rKGdSSNB/YHHi24br2HJB2ctnUw8PCYIzerCF97Y92szIRzCzBH0mxJE8k6Aaysq7OSneOynQDckI5OVgKLUi+22cAc4Oac56vd1snANS14DWaVtmJFpyPoXz6Pk6+0hJPOyZwBXAvcDVwVEWsknSPp7anaxcC0NM/Oh0g9yyJiDXAV2TQJ3wc+EBE7ACR9E/gJ8EpJGySdkrb1aeDNkn4OHJMem3WtIvPeeLj8zvH1OPlyB+/sJx6806qsSO80/zu31777Nj9f06/7v+jgnX3RacCsH/Trl12ZPF3BnnHCMesCEyd2OgKD/CbLvITU75xwzLrACy80L588uZw4rPkFoJ4bpzknHLOKK9L7aevW9sdhmWYXgFpzTjhmFXfSSZ2OwGrlNasdc0w5cXQjJxyzLjd/fqcj6D/Nhg+6/vry4ug2TjhmFTZvXn6dVavaH4ftKm/sNF8E2pgTjlmFrV3bvNwDdXZG3uR2y5aVE0e38cfVrIt5oM5quv/+TkdQTU44ZhXla2+qrVn3aGvMCcesovKuvbHOyuse7fM4u3PCMetSHsqms/K6R598cvPyfuSEY1ZBnka6O+yzz8hlO3b4KKeeE45ZF5o7t9MRGOSPnebeartywjHrQmvWdDoCg/xmNfdW25UTjlnFuDmtu+SN9OBmtZ2ccMy6jIeyqZZVq5o3cb7vfeXFUnVOOGYVUmQCLw9lUz3Nmjife668OKrOCcesQi64oNMRWDt4JtCME45ZF/G1N93JPyQyTjhmFTFlSqcjsD2RNxSR58lxwjGrjGefbV7ukaGr7ZJLmpd7nhwnHLOu4ZGhqy3vmhxzwjGrBF970xuWLGle3u/X5DjhmJm1SN7EbGedVU4cVeWEY9ZhRX71unda92g2T86jj/b3UY4TjlmHnXRSpyOwVsqbJ+e008qJo4pKTTiSFkhaJ2m9pKUNyidJujKVr5Y0q6bs7LR8naRj87Yp6WuS7pV0R7od0e7XZ9YOkyd3OgIbjcWLm5/LeeaZ/j3KKS3hSBoHfBF4CzAXOFFS/QhEpwCPR8ThwHnAuWnducAiYB6wAFguaVyBbf7viDgi3e5o48szG5Nx4/LrbN3a/jistfLO5fTr5GxlHuEcBayPiHsiYhtwBbCwrs5C4NJ0/2pgviSl5VdExPMRcS+wPm2vyDbNKmtoqNMRWLs0O5ezY0d/DndTZsKZCTxQ83hDWtawTkRsB54ApjVZN2+bn5J0p6TzJE1qFJSkUyUNShrcvHnz6F+VWRu5s0D3yjuX04/D3fRyp4GzgVcBrwMOBD7SqFJEXBQRAxExMGPGjDLjsz7na296W5ELQfvtKKfMhLMROLTm8SFpWcM6ksYD+wOPNll3xG1GxKbIPA98laz5zaxreCib7pc3d1HeFNW9psyP9C3AHEmzJU0k6wSwsq7OSmD4dNoJwA0REWn5otSLbTYwB7i52TYlHZz+Cjge+FlbX53ZKBQ5uvFQNt1v1arm7/XQUH/1WCst4aRzMmcA1wJ3A1dFxBpJ50h6e6p2MTBN0nrgQ8DStO4a4CpgLfB94AMRsWOkbaZtXS7pLuAuYDrwyTJep5lZrcsua17+3veWE0cVKHxW8kUDAwMxODjY6TCsxx1zTP7Iwf637C157/mSJfldqatM0q0RMZBbzwlnJyccK0OR5jT/W/aeZu+71N1d5IsmHJ+WNKuYqVM7HYG1Q7PrcvrlB4YTjlmJihzdPP54++Ow8uVdl9MPnQeccMzMSpB3Xc6yZeXE0UlOOGYl8bkbazao5/339/5RjhOOmVlJli+HffYZufykk2DevPLiKZsTjlkJihzdrFjR/jis8y68EKZMGbl87dqsG3UvGt/pAMwsU2TsLet+w+9zs4n38q7T6lY+wjFrsyJHN3ljbllvWbwYDjuseZ1eHNjTCcesjYp+aaxa1d44rHo+9anm5RddVE4cZXLCMWujInOe+OimPy1eDHPr5zyu0YuDtzrhmLVJ0flufHTTv9asGbmsyPTj3cYJx6wNivYyanZdhvWHkT4De++dzYk0a1bvXJ/jXmpmbVC0l1E3jxBsrTH8GbjooqwZba+9sqPjZ57Jlt/zDuguAAAK30lEQVR/P5x6ana/23sy+gjHrMWKNqV5VAEbtnw5bN+efSYOPXT38zdbt8K73tX9RzpOOGYtVDTZeERoG8kvf9l4+dBQdu1ON3eXdsIxa5GiyQY8IrSN7GUva15+wQXdOxKBE45ZC4wm2bgpzZrJuz4HsnOE3Xik44Rjtgek0SUbN6VZnsWLm0/WNuzCC9sfS6s54ZiNwZQpo0s0w9yUZkWcfz5MmNC8ztBQ93WZdrdos1EYS5IZ5qY0K2q4+/Npp+3sHt1It3WZ9hGOWY7hZjMnGyvT4sXw9NMwaVLzelu3Zr3XuuFoxwnHrM64ca1JMsOcbGxPXHxxsWFuho92qpx0nHCs7w2fjxm+DQ21bttONranFi+GSy/Nn84AsqOdZcuypDNrVvWGxnHCsb5TfwTz7LOtf47Jk51srHUWL4b77stmhW02WyjsPNK5//7sM1h75HP66TB+fPa5Hz++/K7VTjjW82qTS6uPYBqJyH5pmrXa4sXZmGvNjnbGjdv987d1a9YB4YILdg6bs2PHzotIyzoacsKxnlHfNNaqczBFRfioxtqv2dHOlCkjz6MzUm+3669vfDTUDk441hUaJZL6WzuaxopworFOqD3akbK/eUc/RQyfB2qHUhOOpAWS1klaL2lpg/JJkq5M5aslzaopOzstXyfp2LxtSpqdtrE+bXNiO15TO07ONdvmcNlwG6wE06dnN2n38xPTp+/ablt7a7Ssqrcq2WuvnUnGicY6afhoZ2go+7t4cTY0TqMjn71G8W0/0gCieywiSrkB44BfAC8HJgI/BebW1Tkd+FK6vwi4Mt2fm+pPAman7Yxrtk3gKmBRuv8lYElejL/9278do7FiRcSUKbVfPdnjFStGtZnC22xU5ls5N7NusmJFxGGHRUjZ3xUrIpYsKf55P+yw0T0fMBg536+R/pVKuQFvAK6teXw2cHZdnWuBN6T744FHANXXHa430jbTOo8A4xs990i30Sacww5rzZtVdJsjlfnW+ptZL1qyJGLcuOwzPm5cxPz5rfnRXDThlNmkNhN4oObxhrSsYZ2I2A48AUxrsu5Iy6cBW9I2RnouACSdKmlQ0uDmzZtH9YJGOuzck8PRZtts22Gu7ZZyzHpR7URv27fDqlWNzwO1a5icvu80EBEXRcRARAzMmDFjVOuONG9F3nwWY93mnmzXdlqyxAnGbFij80DtUmbC2QgcWvP4kLSsYR1J44H9gUebrDvS8keBqWkbIz3XHhvp5FyR+SzGss1GZdZcowaz4TnkzaxcZSacW4A5qffYRLJOASvr6qwETk73TwBuSO2DK4FFqRfbbGAOcPNI20zr3Ji2QdrmNa1+QSN1S9yTXwjNtll/0dfw+ErTpu2cP6O+J8q0adkv+kZjMRUZn6mqip6NMbPqUJT4XynpOOALZL3LLomIT0k6h+yE00pJewOXAa8FHiPrZXZPWncZ8F5gO/DBiPjeSNtMy18OXAEcCNwOnBQRzzeLb2BgIAYHB1v9ss3MepqkWyNiILdemQmn6pxwzMxGr2jC6ftOA2ZmVg4nHDMzK4UTjpmZlcIJx8zMSuFOAzUkbQaeIRsWp1tMx/G2k+Ntn26KFRxvM4dFRO6V8044dSQNFultURWOt70cb/t0U6zgeFvBTWpmZlYKJxwzMyuFE87uLup0AKPkeNvL8bZPN8UKjneP+RyOmZmVwkc4ZmZWCiccMzMrRd8nHElnSvpPSWskfaZm+dmS1ktaJ+nYmuUL0rL1kpZ2KOYPSwpJ09NjSfqHFNOdko6sqXuypJ+n28kjb7Utcf592rd3SvpnSVNryiq7f6sURy1Jh0q6UdLa9Hk9Ky0/UNJ16T2+TtIBafmIn4uS4x4n6XZJ/5Iez5a0OsV1ZZpahDT9yJVp+WpJszoQ61RJV6fP7d2S3lDV/Svpf6bPwc8kfVPS3lXetwC5c1D38g34A2AVMCk9fkn6Oxf4KTAJmA38gmz6g3Hp/suBianO3JJjPhS4FrgfmJ6WHQd8DxDwemB1Wn4gcE/6e0C6f0CJsf4hMD7dPxc4t+r7N8VXiTgaxHUwcGS6vx/wX2lffgZYmpYvrdnPDT8XHYj7Q8A3gH9Jj68im3oE4EvAknT/dOBL6f4i4MoOxHop8L50fyIwtYr7F5gJ3AtMrtmn767yvo2Ivj/CWQJ8OtI8ORHxcFq+ELgiIp6PiHuB9cBR6bY+Iu6JiG1k8+0sLDnm84C/BGp7eywEvh6Zm8hmOz0YOBa4LiIei4jHgeuABWUFGhH/FhHb08ObyGZeHY63qvuXCsWxi4jYFBG3pftPAXeTffEsJPuiJP09Pt0f6XNRGkmHAG8FvpIeC3gTcPUI8Q6/jquB+al+WbHuD/wecDFARGyLiC1Ud/+OByYrm9l4CrCJiu7bYf2ecF4B/G46xPyhpNel5TOBB2rqbUjLRlpeCkkLgY0R8dO6okrGW+e9ZL8GofrxViWOEaUmkdcCq4GDImJTKnoQOCjdr8Lr+ALZD6Sh9HgasKXmh0htTC/Gm8qfSPXLMhvYDHw1NQF+RdI+VHD/RsRG4LPAL8kSzRPArVR33wJZhuxpklYBv96gaBnZ6z+Q7HD4dcBVymYK7ZiceD9K1kxVGc3ijYhrUp1lZDO1Xl5mbL1K0r7AP5HNfPtk7Q/ViAhJlbjWQdLbgIcj4lZJR3c6ngLGA0cCZ0bEaknnkzWhvagq+zedR1pIliS3AP+PElsvxqrnE05EHDNSmaQlwLcia9i8WdIQ2YB3G8nOlQw7JC2jyfK2xivpN8k+XD9NXzCHALdJOqpJvBuBo+uW/6CMeIdJejfwNmB+2s/Qwf1bULP4OkrSBLJkc3lEfCstfkjSwRGxKTXpDDcNd/p1/A7wdmXTwO8N/BpwPlnT0/j0S7s2puF4N6Rmov2BR0uMdwOwISJWp8dXkyWcKu7fY4B7I2IzgKRvke3vqu5bwE1q3ybrOICkV5CdJHwEWAksSj07ZgNzgJuBW4A5qSfIRLKTbyvLCDQi7oqIl0TErIiYRfbPcWREPJhieFfqNfN64InUBHAt8IeSDki/iP4wLSuFpAVkzSlvj4itNUWV2791qhLHLlKb+8XA3RHx+ZqilcBwD8STgWtqljf6XJQiIs6OiEPS53URcENELAZuBE4YId7h13FCql/a0UT6X3pA0ivTovnAWqq5f38JvF7SlPS5GI61kvv2RZ3oqVCVG1mCWQH8DLgNeFNN2TKynkrrgLfULD+OrHfQL8iajToV+33s7KUm4IsppruAgZp67yU7Kb8eeE/JMa4naze+I92+1EX7txJx1MX0RrLOInfW7NPjyNrirwd+Ttbr8sC8z0UHYj+anb3UXk72A2M9WVPQcC/RvdPj9an85R2I8whgMO3jb5P17qzk/gX+D/Cf6fvrMrJen5XdtxHhoW3MzKwc/d6kZmZmJXHCMTOzUjjhmJlZKZxwzMysFE44ZmZWCiccs4pK13f8WNJbapb9iaTvdzIus7Fyt2izCpP0GrLrJ15LNjLI7cCCiPhFRwMzGwMnHLOKUzZP0zPAPsBTEfGJDodkNiZOOGYVl0Ysvg3YRnY1+/MdDslsTHp+8E6zbhcRz0i6Enjayca6mTsNmHWHIXbOKWPWlZxwzMysFE44ZmZWCncaMDOzUvgIx8zMSuGEY2ZmpXDCMTOzUjjhmJlZKZxwzMysFE44ZmZWCiccMzMrxf8H5ewbfxpuE60AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "p_y_given_x = conditional_model.pdf(X['Y'], [4.]*len(X))\n",
    "pp.plot(X['Y'], p_y_given_x, 'bo'); pp.xlabel('Y'); pp.ylabel('$P(Y|X=4)$'); pp.title(\"Conditional Density Plot for P(Y|X=4)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The conditional plots look pretty good!\n",
    "\n",
    "Now that we have $P_{population}(Y|D)$, we need the expected value, $E_{sample}[Yt|Y, D]$. It's linear, so let's just use a linear model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "expectation_model = LinearRegression()\n",
    "expectation_model = expectation_model.fit(X_sampled[['Y', 'X']], X_sampled['Yt'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f5c33a4beb8>]"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGg9JREFUeJzt3X+MXfV55/H3Z/wDPIXNDPas1x4bD1GsVrZWEGuEiFJVWcZJwRvFrJRWRPbiUlZOPGxFN5W6ZvljVWlZJburpiAFE7fQdfBtE5Ymi4XYsnhCtNo/QjNuEhJMWCYUYxsDA8U0wUmo4dk/zvfi62HGc++5v8499/OSRnPO95wz88z1+PHj53zP9yoiMDOz8hrodgBmZtZeTvRmZiXnRG9mVnJO9GZmJedEb2ZWck70ZmYl50RvZlZyTvRmZiXnRG9mVnJLux0AwKpVq2JsbKzbYZiZ9ZQjR468FhEji51XiEQ/NjbG9PR0t8MwM+spko7Vc55bN2ZmJedEb2ZWck70ZmYl50RvZlZyTvRmZiXnRG9m1gWVCoyNwcBA9rlSad/3KsT0SjOzflKpwO7dcOZMtn/sWLYPsGNH67+fK3ozsw67445zSb7qzJlsvB2c6M3M2mxum+bYAo85vfhie76/WzdmZm00X5tGgoj3n3v55e2JwRW9mVkbzdemiciSfa3BQbjzzvbE4ERvZtZGC7VjImDDhizhb9gA+/e350Ys1JnoJQ1JekjSjyU9I+kjki6T9Lik59Ln4XSuJN0taUbSU5K2tCd0M7NimW/K5ELtmA0b4IUX4N13s8/tSvJQf0V/F/DXEfFrwJXAM8BeYCoiNgJTaR/gemBj+tgN7GtpxGZmBVTtxR87llXr1SmT27ZlbZla7WzTzGfRRC/pA8BvAPcBRMTbEXEa2A4cSKcdAG5I29uBr0bmO8CQpDUtj9zMrEAWmjL56KNZW6ZTbZr51DPr5gpgFvhzSVcCR4DbgNURcSqd8zKwOm2PAsdrrj+Rxk7VjCFpN1nFz+XtutVsZtYhC/XiX3wxS+qdTOxz1dO6WQpsAfZFxIeBtzjXpgEgIgKYZ7LQwiJif0SMR8T4yMiib5BiZlZoC9WrRahj60n0J4ATEfFk2n+ILPG/Um3JpM+vpuMngfU1169LY2ZmpXXnnd3vxS9k0UQfES8DxyX9ahqaAI4Ch4BdaWwX8HDaPgTclGbfXAO8WdPiMTMrpR07ut+LX0i9s25+D6hIegq4CvjPwBeAj0t6Dtia9gEeBZ4HZoA/BSZbGrGZWYfVu9Lkjh2dmzLZiLqWQIiI7wPj8xyamOfcAG5tMi4zs0Lo9EqT7eAnY83M5qhUYNWqrAWzc2dnV5psBy9qZmZWY3IS9tXxmGe7VppsB1f0ZmZJvUkeijFtsl5O9GbW96qtmnqTfFGmTdbLid7M+lr1Zuvrr9d3fpGmTdbLPXoz6zuVSnYz9cUXsymT77yz+DXLl8P99/dWgq9yojezvlGpwGc/C2+9dW6sniR/ySVw7729meTBid7M+kSlAjfdlD3MVC8JPvc5uOee9sXVCe7Rm1lfuO22xpL8ypXwwAO9n+TBid7MSqx26YLFbrYuWXJujZqDB+G113q3VTOXWzdmVkpzly64EAkOHChPYp/Lid7MSqdSgV276rvRClkfvqxJHty6MbOSqFTg0kvPrU9Tb5Lfs6ccffgLcaI3s543OZkl95/9bPFz5/biy57kwa0bM+th882Lv5DBwd57qrUVnOjNrCdt3gxHj9Z//pIl/Znkwa0bM+tBW7c2luQHB8s9q2YxTvRm1jMqlWw5gqmp+q+55JL+reSr3Loxs8JrtBdfNTEBhw+3J6Ze4orezAqrOmVy587Gkvwll2QzapzkM67ozayQGnm3p1qu4t+vrope0guSfijp+5Km09hlkh6X9Fz6PJzGJeluSTOSnpK0pZ0/gJmVT94kv2ePk/x8Gmnd/IuIuCoixtP+XmAqIjYCU2kf4HpgY/rYDeT44zKzflSpwIoVjSf5PXsgoj8efsqjmR79duBA2j4A3FAz/tXIfAcYkrSmie9jZiVXTfA7d8IvflH/dVJ/LGHQrHp79AH8b0kBfCUi9gOrI+JUOv4ysDptjwLHa649kcZOYWY2R6MPPlU5wdev3or+1yNiC1lb5lZJv1F7MCKC7B+DuknaLWla0vTs7Gwjl5pZCUxOZhV5o0m+OqPGSb5+dSX6iDiZPr8KfBO4Gnil2pJJn19Np58E1tdcvi6Nzf2a+yNiPCLGR0ZG8v8EZtZztm7Nf7P1pz/t74ef8lg00Uv6FUmXVreBTwA/Ag4Bu9Jpu4CH0/Yh4KY0++Ya4M2aFo+Z9bFKBZYubezJVoCLL3YV34x6evSrgW9Kqp7/FxHx15K+Czwo6RbgGPDb6fxHgW3ADHAGuLnlUZtZz9m6tfEED+7Ft8KiiT4ingeunGf8dWBinvEAbm1JdGbW8/LebN20CZ5+uvXx9CMvgWBmbVGp5LvZCtnTrU7yreMlEMys5UZH4aWXGr9Oggce8M3WVnOiN7OWyduHB1i7Fk6+b36etYJbN2bWEqOj+ZP8wYNO8u3kRG9mTalUYGAgX6tm7dpsjRq3atrLrRszy214GE6fbvy6FSvgzJnWx2Pzc0VvZg2rzqjJk+T37HGS7zRX9GbWkLzz4oeG4I03Wh+PLc4VvZnVbXQ0/7x4J/nucUVvZovKVkBpnJ9uLQZX9Ga2oOpSwo0aGspm0zjJF4MrejObl6v48nBFb2bnGR3Nl+SXLMkefHKSLx4nejMDzk2ZbPTBJylL8GfP+sGnonLrxsxyL0Lm9Wl6gyt6sz62dWu+Kh6yKZNO8r3BFb1ZnxochJ//vPHrli2Dt99ufTzWPq7ozfpMtRefJ8lPTDjJ9yJX9GZ9JO+UScjmxVtvckVv1gfyPvgE55YStt7lit6s5JYsgXffzXetE3w51F3RS1oi6XuSHkn7V0h6UtKMpK9LWp7GL0r7M+n4WHtCN7MLqc6oyZPkXcWXSyOtm9uAZ2r2vwh8KSI+BLwB3JLGbwHeSONfSueZWQdJ+d7Wr7pGjadNlktdiV7SOuBfAn+W9gVcCzyUTjkA3JC2t6d90vGJdL6Ztdnmzfl78RFeSris6u3R/wnwh8ClaX8lcDoizqb9E8Bo2h4FjgNExFlJb6bzX2tJxGY2r7wJ3vPiy2/Ril7SJ4FXI+JIK7+xpN2SpiVNz87OtvJLm/WV4eH8SX7PHif5flBP6+ajwKckvQB8jaxlcxcwJKn6P4J1QLWrdxJYD5COfwB4fe4XjYj9ETEeEeMjIyNN/RBm/ag6ZTLvm3NHwD33tD4uK55FE31E3B4R6yJiDLgR+FZE7ACeAD6dTtsFPJy2D6V90vFvRfj+vVkrSbBvX75rJyb85tz9ppkHpv498HlJM2Q9+PvS+H3AyjT+eWBvcyGaWVUzDz4NDGRV/OHDrY3Jiq+hB6Yi4tvAt9P288DV85zzC+C3WhCbmdXwg0+Wl5dAMCu4Zh58qs6Lt/7mJRDMCsyLkFkruKI3K6DBwfxJ3lW8zeWK3qxgXMVbq7miNyuIZh582rTJSd4W5orerABcxVs7uaI366Lly13FW/u5ojfrElfx1ilO9GYd5gRvnebWjVmHNLN8ATjJW36u6M06oJkEPzHh9WmsOU70Zm20eTMcPZrvWr8hiLWKWzdmbSLlT/IHDzrJW+u4ojdrsWaq+IEBeOed1sZj5kRv1kK+2WpF5NaNWQs0swjZsmVO8tZerujNmuQq3orOFb1ZTqOjXkrYeoMrerMcXMVbL3FFb9aAZqr4iQkneesOV/RmdXIVb73Kid5sEU7w1usWbd1IuljS30j6gaSnJf1RGr9C0pOSZiR9XdLyNH5R2p9Jx8fa+yOYtcfWrU7yVg719Oh/CVwbEVcCVwHXSboG+CLwpYj4EPAGcEs6/xbgjTT+pXSeWU+RYGoq37URTvJWLIsm+sj8LO0uSx8BXAs8lMYPADek7e1pn3R8QmqmLjLrnM2bXcVb+dTVo5e0BDgCfAj4MvAT4HREnE2nnABG0/YocBwgIs5KehNYCbzWwrjNWs4J3sqqrumVEfFORFwFrAOuBn6t2W8sabekaUnTs7OzzX45s9yamTIJTvJWfA3No4+I08ATwEeAIUnV/xGsA06m7ZPAeoB0/APA6/N8rf0RMR4R4yMjIznDN2uOBC+9lO9a9+KtV9Qz62ZE0lDaXgF8HHiGLOF/Op22C3g4bR9K+6Tj34rwXwcrluHh/FX8wIATvPWWenr0a4ADqU8/ADwYEY9IOgp8TdJ/Ar4H3JfOvw94QNIM8PfAjW2I2yw3t2ms3yya6CPiKeDD84w/T9avnzv+C+C3WhKdWQs1k+D9tn7Wy/xkrPUFV/HWz7yomZVaM28I4kXIrCyc6K2UKpUswf/85/muj4DDh1sbk1m3uHVjpdNMm+bgQdixo3WxmBWBK3orjcnJ5nvxTvJWRk70VgoS7NuX79o9e9yLt3Jz68Z62uho/idbwQne+oMTvfUsT5k0q49bN9Zzmlm+AJzkrf+4oree4gRv1jhX9NYTli/Pn+TXrnWSt/7mit4Kz1W8WXNc0VthSfmT/NCQk7xZlSt6KyRX8Wat44reCqWZKn7TJid5s/k40VshbN3afBX/9NOti8esTNy6sa5rJsGvXQsnTy5+nlk/c0VvXdOKKt5J3mxxTvTWFRJMTeW71m8IYtYYt26sowYH878ZyIoVcOZMa+Mx6wdO9NYxnjJp1h2Ltm4krZf0hKSjkp6WdFsav0zS45KeS5+H07gk3S1pRtJTkra0+4ewYluyxEnerJvq6dGfBf4gIjYB1wC3StoE7AWmImIjMJX2Aa4HNqaP3UDOt4OwMpDg3XfzXRvhJG/WCosm+og4FRF/m7Z/CjwDjALbgQPptAPADWl7O/DVyHwHGJK0puWRW6ENDrqKNyuKhmbdSBoDPgw8CayOiFPp0MvA6rQ9ChyvuexEGrM+IeW/4eoq3qz16k70ki4B/gr4/Yj4h9pjERFAQ389Je2WNC1penZ2tpFLraCaWb5g2TIneLN2qSvRS1pGluQrEfGNNPxKtSWTPr+axk8C62suX5fGzhMR+yNiPCLGR0ZG8sZvBVCpNN+mefvt1sVjZuerZ9aNgPuAZyLij2sOHQJ2pe1dwMM14zel2TfXAG/WtHisZCTYuTPftStWuIo364R6KvqPAv8auFbS99PHNuALwMclPQdsTfsAjwLPAzPAnwKTrQ/bum1ysvkq3g8/mXXGog9MRcT/BRb6Kz0xz/kB3NpkXFZgzST4iQk4fLh1sZjZ4rzWjdWtFVW8k7xZ53kJBKuLq3iz3uVEbxc0OgovvZT/et9sNes+J3pbkJ9sNSsH9+jtfZpZhMwPPpkVjyt6O4+reLPycUVvQHNv6+cq3qzYXNGbq3izknNF38eamRe/dq2TvFmvcEXfp1zFm/UPV/R9ppk3BNmzx0nerBe5ou8jruLN+pMr+j7QTBV/8KCTvFmvc6IvserN1jxv61edMrljR+vjMrPOcuumpIaH4fTpfNe6gjcrF1f0JbN5c1bFO8mbWZUr+hLxzVYzm48r+hJo5sEn32w1Kz9X9D2umfXineDN+oMr+h41PJxV8XmSvB98Musvruh7kHvxZtYIV/Q9pFrF5zEx4SRv1q8WTfSS7pf0qqQf1YxdJulxSc+lz8NpXJLuljQj6SlJW9oZfD/JO2Wy+uCT35zbrH/VU9H/d+C6OWN7gamI2AhMpX2A64GN6WM3sK81YfavZt4Q5OBBePvt1sZjZr1n0UQfEf8H+Ps5w9uBA2n7AHBDzfhXI/MdYEjSmlYF208qlSzBT001fu3QkJcvMLNz8vboV0fEqbT9MrA6bY8Cx2vOO5HG3kfSbknTkqZnZ2dzhlFOo6Owc2e+aycm4I03WhuPmfW2pm/GRkQADd/mi4j9ETEeEeMjIyPNhlEK1Qef8kyZ3LTJvXgzm1/e6ZWvSFoTEadSa+bVNH4SWF9z3ro0Zovwg09m1i55K/pDwK60vQt4uGb8pjT75hrgzZoWj80jbxU/MJAleCd5M1vMohW9pL8EPgasknQC+I/AF4AHJd0CHAN+O53+KLANmAHOADe3IebSmJyEfTnmJQ0NuQ9vZvVbNNFHxGcWODQxz7kB3NpsUGVXqeS/2XrwoGfTmFljvARCh+Wt4lesgDNnWh+PmZWfE30H5bnhunYtnPTtbDNrgte66YDq062NJvmDB53kzax5rujbLM/yBcuWeekCM2sdV/RtUJ0ymSfJ79njJG9mreWKvsXyPvjkKt7M2sUVfYtUFyHL+45PTvJm1i6u6Ftg69Z8q0x6Ro2ZdYIr+iZNTuZL8p5RY2ad4oq+CZVK4w8/+cEnM+s0V/Q5VCpw8cWNLWOwYkW2AJmTvJl1mhN9gyYnswT/y1/Wf83EhBO8mXWPWzd1ynPDVYIHHvAiZGbWXU70dciT5Ccm/G5PZlYMbt1cQKUCY2ONJfmVK7MZNU7yZlYUrujnUanAbbfB6683dp2reDMrIlf0NSoVuOii7GZrI0leyp5udZI3syJyRZ9UKnDTTfDuu/Vfs3w53H+/b7aaWbE50Sd33NFYknebxsx6RV+3bqo3WwcG4Nixxc8fHMxutEY4yZtZ7+jLir5Sgc9+Ft56q/5rVq6Eu+5ym8bMek9bKnpJ10l6VtKMpL3t+B55VSpw882NJfmJCXjtNSd5M+tNLU/0kpYAXwauBzYBn5G0qdXfJ49KBXbtgn/8x/rOHxjwbBoz633taN1cDcxExPMAkr4GbAeOtuF71a1Sgd274Z13Lnzehg3wwgsdCcnMrCPa0boZBY7X7J9IY11RveG6c+fiC4tJcOedHQnLzKxjunYzVtJuYDfA5Zdf3vKvn+fp1s99zn14MyufdlT0J4H1Nfvr0th5ImJ/RIxHxPjIyEhLA6i2aepN8tX1ae65p6VhmJkVQjsS/XeBjZKukLQcuBE41Ibvs6A77qhv/ffqvHjPqDGzMmt5oo+Is8C/BR4DngEejIinW/19LuTFFxc/Z8MG2L/fCd7Myq8tPfqIeBR4tB1fux6XX77wk66Dg07wZtZfenYJhNrlC8bGsv2qO+/MEvpcK1c6yZtZ/+nJRF+92XrsWLbuzLFj2X412e/YkSX0DRuyKZMbNrgXb2b9SxHR7RgYHx+P6enpus8fG5u/NeOHncysn0g6EhHji53XkxX9Qjdb67kJa2bWb3oy0S/0fFUbnrsyM+t5PZno57vZOjjo5QvMzObTk4l+vputnk1jZja/nn3jkR07nNjNzOrRkxW9mZnVz4nezKzknOjNzErOid7MrOSc6M3MSq4QSyBImgUWWG+yZVYBr7X5e7SLY+8Ox94djr1+GyJi0XduKkSi7wRJ0/WsCVFEjr07HHt3OPbWc+vGzKzknOjNzEqunxL9/m4H0ATH3h2OvTsce4v1TY/ezKxf9VNFb2bWl0qb6CX9gaSQtCrtS9LdkmYkPSVpS825uyQ9lz52dTHm/yrpxym+b0oaqjl2e4r9WUm/WTN+XRqbkbS3O5HPr+CxrZf0hKSjkp6WdFsav0zS4+l34XFJw2l8wd+fbpG0RNL3JD2S9q+Q9GSK8euSlqfxi9L+TDo+1s24U0xDkh5Kv+/PSPpIr7z2kv5d+p35kaS/lHRx4V/7iCjdB7AeeIxsbv6qNLYN+F+AgGuAJ9P4ZcDz6fNw2h7uUtyfAJam7S8CX0zbm4AfABcBVwA/AZakj58AHwSWp3M2dfv1TzEXNrYU3xpgS9q+FPh/6XX+L8DeNL635s9g3t+fLv8Mnwf+Angk7T8I3Ji27wX2pO1J4N60fSPw9QLEfgD4N2l7OTDUC689MAr8HbCi5jX/naK/9l39w27jH8ZDwJXACzWJ/ivAZ2rOeTb9Zf8M8JWa8fPO6+LP8K+AStq+Hbi95thjwEfSx2M14+ed1+X4CxvbAvE+DHy8+nuRxtYAz17o96eL8a4DpoBrgUdSEnyNc4XCe69/9fclbS9N56mLsX8gJUvNGS/8a58S/XGywnBpeu1/s+ivfelaN5K2Aycj4gdzDlX/gKpOpLGFxrvtd8mqGOi92KHYsZ0n/Xf6w8CTwOqIOJUOvQysTttF+3n+BPhD4N20vxI4HRFn035tfO/Fno6/mc7vliuAWeDPU+vpzyT9Cj3w2kfESeC/AS8Cp8heyyMU/LXvyTcekXQY+GfzHLoD+A9kLZBCulDsEfFwOucO4CxQ6WRs/UjSJcBfAb8fEf8g6b1jERGSCjctTdIngVcj4oikj3U7nhyWAluA34uIJyXdRdaqeU+BX/thYDvZP1angf8BXNfVoOrQk4k+IrbONy7pn5P9Afwg/YVdB/ytpKuBk2S9+6p1aewk8LE5499uedDJQrFXSfod4JPARKT/77Fw7FxgvNsuFHMhSFpGluQrEfGNNPyKpDURcUrSGuDVNF6kn+ejwKckbQMuBv4JcBcwJGlpqhxr46vGfkLSUrLWyeudD/s9J4ATEfFk2n+ILNH3wmu/Ffi7iJgFkPQNsj+PQr/2pWrdRMQPI+KfRsRYRIyR/UJtiYiXgUPATekO/jXAm+m/iY8Bn5A0nP61/kQa6zhJ15H9d/xTEXGm5tAh4MZ0B/8KYCPwN8B3gY3pjv9ysps9hzod9wKKHBvKKoH7gGci4o9rDh0CqjOvdpH17qvj8/3+dFxE3B4R69Lv+I3AtyJiB/AE8Ol02tzYqz/Tp9P5XauW09/H45J+NQ1NAEfpgdeerGVzjaTB9DtUjb3Yr303bmh06oPzb8YK+DLZTJAfAuM15/0uMJM+bu5ivDNk/bzvp497a47dkWJ/Fri+Znwb2YyRn5C1f7r+uvdIbL8OBPBUzeu9jax/OgU8BxwGLlvs96fLP8fHODfr5oNkBcAMWUvhojR+cdqfScc/WIC4rwKm0+v/P8lmvPXEaw/8EfBj4EfAA2Sz4Qr92vvJWDOzkitV68bMzN7Pid7MrOSc6M3MSs6J3sys5JzozcxKzonezKzknOjNzErOid7MrOT+P/LZiBNrE/cYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df = X_sampled.copy()\n",
    "df['X'] = 4\n",
    "E_yt_given_d_x = expectation_model.predict(df[['Y', 'X']])\n",
    "pp.plot(df['Y'], E_yt_given_d_x, 'bo')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time to integrate!\n",
    "\n",
    "Again, we're estimating \n",
    "$E_{population}[Yt|X] = \\int E_{sample}[Yt|Y, X]P_{population}(Y|X) dY$\n",
    "\n",
    "For more reference on the integrator, check out https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html#scipy.integrate.quad\n",
    "\n",
    "The arguments are the integrand (a function), and the upper and lower bounds for the integration. We'll integrate from the min to the max of the observed Y values.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import quad\n",
    "\n",
    "y_min = X['Y'].min()\n",
    "y_max = X['Y'].max()\n",
    "\n",
    "def integrand(Y, X):\n",
    "    return (expectation_model.predict([[Y, X]]) * conditional_model.pdf([Y], [X]))[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.10752792283294876"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integrand(100, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Having set up the integrand (the argument of the integral), we can estimate $E_{population}[Yt|X]$ at different X values!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(395.37120799067134, 8.243781150341232e-09)"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = 4\n",
    "quad(integrand, y_min, y_max, args=(x,))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That means we can plot the expected article performance when translated as a function of the author's skill level"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'Performance')"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEiCAYAAADjxEWuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmYXFWd//H3JxAIhkBY2kzIQhhFMaAGbEFFRxYZAcXgKJvKNoyRmaD4DCIo4xAfYUbGhdGfsoRBFhHCokhEFpF1kDVgDISABgiTDgFCCLtEAt/fH+cUuSlvd1cnXXWr05/X89TTdc+5y/ferrrfuts5igjMzMzqDak6ADMza09OEGZmVsoJwszMSjlBmJlZKScIMzMr5QRhZmalnCCslKTDJN3awuV9UtJCSS9K2r5Vy21XkkZJukXSC5K+V3U8Ra38bEiaICkkrdvEZZwh6RsNjHeTpH9qVhztyAmin0haIOnPeQdXe/2owngG2of5u8BREbFhRPy+6mDawBTgaWCjiDimvlLSuZJOan1YfSNpmqQLqo6jpiy5RcSREfGtqmJqZ03LyoPUPhHx26qDGEgkrRsRK4AtgbmrOY91IuK1/o2sclsCD8RqPsla2K6WNfMoZK0VEX71wwtYAHykm7rTgZ8Xhk8BrgcE7AJ0AV8n/WJcAHy2MO76pF/X/wc8CZwBbFConwzMBp4HHgb2BE4GXgNeAV4EfpTH3Qa4DngGeAjYvzCfzYCZeT53Ad8Cbu1mfSYAQfqV+ziwGPhKoX4IcHyOZylwCbBp3bRH5HW6PccYwEvAw3m8dwA3Ac+SEscnCvM/N2/Tq/I0H8llpwFX5/n9Dvgb4L+BZcCDwPaFedTiewF4APhkoe4w4Na83ZcBjwJ7Feo3Bc7J674M+GWh7uP5//EscBvwrh4+Mx8A7gaey38/UFi/V4G/5HX5SN10U+rqf1X4DB4HzAGWk34Arsl6HgY8kqd9lPy5rE1XGO8HwELSZ+ce4EO5fM8c46s5zj/k8o2Bs0mfm0XAScA6uW6dHM/TedlTSZ+NdbvZhr2t3++AU0mfw5+TvhOv5XieLWzvk3r6TuXym4B/Koz3j8C8vO2uBbbM5crLfCrP4z5gu6r3Uau1X6s6gLXlRc8J4k3AH/MH9kP5wz821+0CrAC+T0oGHybt9N6e608l7bg3BUYAvwL+M9ftSNq57EHaKY8Btsl19R/m4flLfDhpx7F9jmNirp9B2pEPB7bLX9zeEsRFefx3Aktq6w8cDdwBjM3rdCZwUd205+dpN8jlAbw1vx8KzCclzfWA3fIOoLZNzs3rvXNe72G57GngPXn4BtJO7RDSTuck4MbCOuwHbJGnPyBv89G57jDSTu3zedp/JiUD5fpfAxcDm+RYP5zLtyftFHbK0x1K+lysX7INNyXtWA7O/4+D8vBmhXU8qWz7d1eflzUbGFfYrqu1nvl/83xhm48Gti1MV0wQnyP9wFgXOAZ4AhiW66YBF9TFeTnpMzEceDPpB8kXct2RpGQ+Lm+jG+k5QfS2fiuAL+bYNqiPvX5b0uB3ipRE5pN+yKwL/BtwW677KClRjszb8h21mAbaq/IA1pZX/nK+SPrlWHt9vlC/E+mX+2PAQYXyXfKHeHih7BLgG/nD9RLwlkLd+4FH8/szgVO7ieeND3MePgD437pxzgROzDuHV2tfhFz3H/VfpELdhPylLY7/X8DZ+f08YPdC3eg8/3UL0/5t3TyLCeJDpJ3MkEL9RcC0/P5c4Py66c8FzioMfxGYVxh+J/kXYzfrNBuYnN8fBswv1L0px/c3eV1eBzYpmcfpwLfqyh4iJ5C68oOBu+rKbgcOK6zP6iSIf+zlc9roeg4nfYY/ReGItTBd6Wcj1y8D3p3fT6OQIIBRpKOb4lHwQeTkTUrsRxbq/p4eEkQD6/d/vcXOqgmioe8U6Uj1iELdEOBl0qnB3Ug/CN9X/AwPxJcvUvevfSNiZOF1Vq0iIu4kHTKLlACKlkXES4Xhx0i/ijpIX9p7JD0r6VngmlwO6VfWww3GtiWwU20+eV6fJe0MOkg774V1MfSmfvwtCsu6vLCceaTD+lHdTFtvC2BhRLxeN/8xvUz/ZOH9n0uGN6wNSDpE0uxCjNsBmxfGf6L2JiJezm83JG3zZyJiWcnytwSOqdvG41i5XerXsX4b16/j6lhlu6zueubP4wGkX/SLJf1a0jZlC5T0FUnzJD2Xl7Fx3TKKtiQddS0uxHQm6UgC8v++MH6Pn8MG1q+nz1mZRr9TWwI/KCz3GdJ3e0xE3AD8CPgx8JSk6ZI26mMcbcEJokUkTSWdbnkc+Gpd9SaShheGx+fxnibt2LYtJJ2NI6K2o1sIvKWbRUbd8ELg5roEtmFE/DPp9NAK0pejGENv6sd/vLCsveqWNSwiFvUQX9HjwDhJxc/neNJpr0am75GkLYGzgKNIp3RGAveTvuC9WQhsKmlkN3Un1633myLiopJxHyftZIrq17En3a3/G+VruJ5ExLURsQfpqOnBPK9VSPoQ6fO8P+moaiTpFE1tGWWfw+XA5oVttFFEbJvrF9Pg57DB9atffm+fm56+U/XjfaHuf71BRNwGEBE/jIj3ABOBtwHHNjDPtuME0QKS3kY6B/450qmFr0qaVDfaNyWtl79wHwcuzb+gzwJOlfTmPK8xkj6apzkbOFzS7pKG5Lrar7wngb8tzP9K4G2SDpY0NL/eK+kdke4A+gUwTdKbJE0knT/vzTfy+NuSrm1cnMvPAE7OX2AkdUia3ODmAriTdLj+1RznLsA+pOsk/WE4aUexJMd3OOmXZ68iYjHp9MJpkjbJ8f1drj4LOFLSTkqGS/qYpBEls7qK9P/4jKR1JR1A2plc2eA61P9/y6z2eubnMCbnHy7LSadPXy8ZdQTpx8USYF1J/w4Ufy0/CUyoJfu8/X4DfE/SRvlz+xZJH87jXwJ8SdJYSZuQLkL35/o9CYyVtF439T19p4rOAL6WP/tI2ljSfvn9e/NnYCjpFPErlG+7tucE0b9+VfccxOX51roLgFMi4g8R8SfSxdefSlo/T/cE6bzt48DPSOdgH8x1x5Euht0h6Xngt8DbASLiLtKO+VTSr7abWfmr9AfApyUtk/TDiHiBdD73wLycJ0h3U9ViOIp0CuUJ0jnZcxpY35tzbNcD342I3xSWPRP4jaQXSBesd2pgfuT1+gspIexFOoo6DTiksE3WSEQ8AHyPdM7/SdL1id/1YRYHk66pPEi6KP3lPN9ZpAu+PyL9P+eTznmXxbCU9EPgGNIdNl8FPh4RTzcYw9nAxHyK45fdLGNN1nMI8K+kz8ozpJsn/rlkvGtJpz3/SDod9Aqrnta5NP9dKune/P4Q0s0HD5C202WkoxRISfZa4A/AvaQfLqVWc/1uIN0V94Skv9rWvXyniuNdTvr+zMjfy/tJn1dICfKsvG6Pkf6/3+klrrZUuyvDKpJ/HV8QEWOrjqVRkiaQ7hAaGr7X3myt5SMIMzMr5QRhZmalfIrJzMxK+QjCzMxKOUHYgKPUcu5Hqo6jnlrcRHphubtI6uqm7kOSHioMv7Ht1GYtrVr7cYIwWw1qQT8F/SEi/jci3l51HDYwOUGYtaF2Tzw2ODhB2ED1XkkP5AcBz5E0DEDS/ZL2qY2Un3R+WiW91OUnoa+UtCTP50pJYwv1q5zKqjslc0v++2x+KPL9hfG+m+f3qKS9CuVbSJop6RlJ8yV9vm7el0m6ID94dVhJvHvndX5B0iJJXynbMJK+lMcb29PpJ7PeOEHYQPVZUrPKbyG1dfNvufx8UpMmNXsDi6O8l7ohpCfGtyS1+fNn0lPQjag1r1Fr0+r2PLwTqQXXzckt3EqqtQ00g9T3xxbAp4H/kLRbYZ6TSU8VjyQ9UV/vbFL7PyNITUrcUD9CburiMFILsk4MtkacIGyg+lFELIyIZ0gdJB2Uyy8A9tbK1jMPBn5aNoOIWBoRP4+Il3NTJCeTmpRYE49FxFm5favzSE1IjJI0jtR/xXER8UpEzAb+h9TsRM3tEfHLiHg9Iv5cMu9XSc1rbBQRyyLi3kKdJH2f1JzKrhGxZA3Xw8wJwgas0qbGI+JxUns8n1JqcXUvyn+No9TQ4JmSHsundW4BRkpaZw3i6q6Z8C1IzYS/UBd3b02YF32KdET0mKSbi6e1SEcdU0idST23usGbFTlB2EDVXVPjkH65f47U29jtdc2MFx1Davhwp4jYiJWnjWqnhF4i9cdR8zeF9319wvRxUjPhxZZd+9SEeUTcHRGTSX0n/JJV+xVZRmr87xxJO/cxNrNSThA2UE3NF2E3BU5gZVPjkHaeO5C6Pj2/h3mMIF13eDbP58S6+tnAgflCdyfpukHNElITzr01uQ1ARCwk9VH9n5KGSXoXqV/uhp5DUGoK/rOSNo6IV0ndga7ShHRE3ES6NvMLSTs2Ml+znjhB2EB1IalfgUdIPYCdVKvI5+9/DmxFD81FA/9N6qf4aVKT5NfU1X+DdBF8GfDNvMzaMl4mXbP4XW5y+30NxHwQqcvVx0n9Mp8YEb9tYLqag4EF+XTYkaRksIqIuA74R1LT8zv0Yd5mf8VtMdlaKd/N87aI+FyvI5tZKT+MY2udfLroCNIvbjNbTT7FZGuV/PDZQuDqiLilt/HNrHs+xWRmZqV8BGFmZqWcIMzMrNSAvki9+eabx4QJE6oOw8xsQLnnnnuejoiO3sYb0AliwoQJzJo1q+owzMwGFEmPNTKeTzGZmVkpJwgzMyvV9AQhaR1Jv5d0ZR7eStKducOUiyWtl8vXz8Pzc/2EZsdmZmbda8U1iKOBeUCtff5TgFMjYoakM0hPvJ6e/y6LiLdKOjCPd0BfF/bqq6/S1dXFK6+80j/Rt8CwYcMYO3YsQ4cOrToUM7M3NDVB5O4bP0Zq1Oxfc89auwGfyaOcB0wjJYjJ+T2kXrV+JEnRxyf5urq6GDFiBBMmTGBlR17tKyJYunQpXV1dbLXVVlWHY2b2hmafYvpv4KusbJZ4M+DZiFiRh7tY2WHKGHKHKbn+uTx+n7zyyitsttlmAyI5AEhis802G1BHPGY2ODQtQUj6OPBURNzTz/OdImmWpFlLlpT3qjhQkkPNQIvXzAaHZh5B7Ax8QtICUmftuwE/IHXpWDu1NZaVPWotIvcSlus3BpbWzzQipkdEZ0R0dnT0+pxHy0UEH/zgB7n66qvfKLv00kvZc889Oe200yqMzMysb5p2DSIivgZ8DUDSLsBXIuKzki4l9cw1AzgUuCJPMjMP357rb+jr9YcyE47/9ZrOYhULvv2xHuslccYZZ7Dffvux6667smLFCr7+9a9zzTXXMHnyZP7lX/6lX+Mxs9b58ZE3VB0CU8/YrWXLquJJ6uOAGZJOAn4PnJ3LzwZ+Kmk+8AxwYAWx9YvtttuOffbZh1NOOYWXXnqJQw45hBNOOIGHH36YSZMmsccee/Cd73yn6jDNzHrUkgSR+8q9Kb9/BPir/nIj4hVSJ/NrhRNPPJEddtiB9dZbj1mzZrF48WLuv/9+Zs+eXXVoZmYNGdBtMbWz4cOHc8ABB7Dhhhuy/vrrVx2OmVmfuamNJhoyZAhDhngTm9nA5L1Xi4wYMYIXXnih6jDMzBrmBNEim222GTvvvDPbbbcdxx57bNXhmJn1aq2/BtHbbanNNG3atFWGL7zwwmoCMTNbDT6CMDOzUk4QZmZWygnCzMxKrZUJoh9a6GipgRavmQ0Oa12CGDZsGEuXLh0wO91afxDDhg2rOhQzs1WsdXcxjR07lq6uLrprCrwd1XqUMzNrJ2tdghg6dKh7ZjMz6wdr3SkmMzPrH04QZmZWygnCzMxKOUGYmVmppiUIScMk3SXpD5LmSvpmLj9X0qOSZufXpFwuST+UNF/SHEk7NCs2MzPrXTPvYloO7BYRL0oaCtwq6epcd2xEXFY3/l7A1vm1E3B6/mtmZhVo2hFEJC/mwaH51dPTa5OB8/N0dwAjJY1uVnxmZtazpl6DkLSOpNnAU8B1EXFnrjo5n0Y6VVKtP84xwMLC5F25rH6eUyTNkjRrID0MZ2Y20DQ1QUTEaxExCRgL7ChpO+BrwDbAe4FNgeP6OM/pEdEZEZ0dHR39HrOZmSUtuYspIp4FbgT2jIjF+TTScuAcYMc82iJgXGGysbnMzMwq0My7mDokjczvNwD2AB6sXVeQJGBf4P48yUzgkHw30/uA5yJicbPiMzOznjXzLqbRwHmS1iEloksi4kpJN0jqAATMBo7M418F7A3MB14GDm9ibGZm1oumJYiImANsX1K+WzfjBzC1WfGYmVnf+ElqMzMr5QRhZmalnCDMzKyUE4SZmZVygjAzs1JOEGZmVsoJwszMSjlBmJlZKScIMzMr5QRhZmalnCDMzKyUE4SZmZVygjAzs1JOEGZmVsoJwszMSjlBmJlZqWZ2OTpM0l2S/iBprqRv5vKtJN0pab6kiyWtl8vXz8Pzc/2EZsVmZma9a+YRxHJgt4h4NzAJ2DP3NX0KcGpEvBVYBhyRxz8CWJbLT83jmZlZRZqWICJ5MQ8Oza8AdgMuy+XnAfvm95PzMLl+d0lqVnxmZtazpl6DkLSOpNnAU8B1wMPAsxGxIo/SBYzJ78cACwFy/XPAZiXznCJplqRZS5YsaWb4ZmaDWlMTRES8FhGTgLHAjsA2/TDP6RHRGRGdHR0daxyjmZmVW7cVC4mIZyXdCLwfGClp3XyUMBZYlEdbBIwDuiStC2wMLG1FfGbWvXnbvKPqEHjHg/OqDmFQauZdTB2SRub3GwB7APOAG4FP59EOBa7I72fmYXL9DRERzYrPzMx61swjiNHAeZLWISWiSyLiSkkPADMknQT8Hjg7j3828FNJ84FngAObGJuZmfWiaQkiIuYA25eUP0K6HlFf/gqwX7PiMTOzvvGT1GZmVsoJwszMSjlBmJlZKScIMzMr5QRhZmalnCDMzKyUE4SZmZVygjAzs1JOEGZmVsoJwszMSjlBmJlZKScIMzMr5QRhZmalnCDMzKyUE4SZmZVqZo9y4yTdKOkBSXMlHZ3Lp0laJGl2fu1dmOZrkuZLekjSR5sVm5mZ9a7hDoNyt6HjI+KhBidZARwTEfdKGgHcI+m6XHdqRHy3bv4TSb3IbQtsAfxW0tsi4rVGYzQzs/7T0BGEpH2A2cA1eXiSpJk9TRMRiyPi3vz+BVJ/1GN6mGQyMCMilkfEo8B8SnqeMzOz1mj0FNM00s76WYCImA1s1ehCJE0gdT96Zy46StIcST+RtEkuGwMsLEzWRc8JxczMmqjRBPFqRDxXVxaNTChpQ+DnwJcj4nngdOAtwCRgMfC9BmOozW+KpFmSZi1ZsqQvk5qZWR80miDmSvoMsI6krSX9P+C23iaSNJSUHH4WEb8AiIgnI+K1iHgdOIuVp5EWAeMKk4/NZauIiOkR0RkRnR0dHQ2Gb2ZmfdVogvgi6eLxcuBC4Dngyz1NIEnA2cC8iPh+oXx0YbRPAvfn9zOBAyWtL2krYGvgrgbjMzOzftbQXUwR8TJwQn41amfgYOA+SbNz2deBgyRNIp2iWgB8IS9jrqRLgAdId0BN9R1MZmbVaShB5NtT94uIZ/PwJqQ7jrp9ViEibgVUUnVVD9OcDJzcSExmZtZcjZ5i2ryWHAAiYhnw5uaEZGZm7aDRBPG6pPG1AUlb0uBdTGZmNjA1+iT1CcCtkm4mnTb6EDClaVGZmVnlGr1IfY2kHYD35aIvR8TTzQvLzMyq1nBbTMD6wDN5momSiIhbmhOWmZlVrdG7mE4BDgDmAq/n4gCcIMzM1lKNHkHsC7w9IpY3MxgzM2sfjd7F9AgwtJmBmJlZe2n0COJlYLak60nNbQAQEV9qSlRmZla5RhPEzPwyM7NBotHbXM9rdiBmZtZeGr2LaWvgP4GJwLBaeUT8bZPiMjOzijV6kfocUkc/K4BdgfOBC5oVlJmZVa/RBLFBRFwPKCIei4hpwMeaF5aZmVWt0YvUyyUNAf4k6ShST28bNi8sMzOrWqNHEEcDbwK+BLyH1BHQoc0KyszMqtdQgoiIuyPixYjoiojDI+IfIuKOnqaRNE7SjZIekDRX0tG5fFNJ10n6U/67SS6XpB9Kmi9pTm4c0MzMKtJQgpDUKelySffmnfccSXN6mWwFcExETCS1AjtV0kTgeOD6iNgauD4PA+xF6od6a1JT4qevxvqYmVk/afQaxM+AY4H7WNlYX48iYjGwOL9/QdI8YAwwGdglj3YecBNwXC4/PyICuEPSSEmj83zMzKzFGk0QSyJitZ+kljQB2B64ExhV2Ok/AYzK78cACwuTdeUyJwgzswo0miBOlPQ/pFNCxbaYftHbhJI2BH5O6mToeUlv1EVESOpT16WSppB7sxs/fnwvY5uZ2epqNEEcDmxDatG12B9EjwlC0lBScvhZIZk8WTt1JGk08FQuXwSMK0w+NpetIiKmA9MBOjs73S+2NcU7z3tn1SFw36H3VR2CDXKNJoj3RsTb+zJjpUOFs4F5EfH9QtVM0i2y385/ryiUHyVpBrAT8JyvP5iZVafRBHGbpIkR8UAf5r0z6XmJ+yTNzmVfJyWGSyQdATwG7J/rrgL2BuaTmhc/vA/LMjOzftZogngfqT+IR0nXIES6hPCu7iaIiFvzeGV2Lxk/gKkNxmNmZk3WaILYs6lRmJlZ2+k1QUhaB7g2IrZpQTxmZtYmen2SOiJeAx6S5HtKzcwGkUZPMW0CzJV0F/BSrTAiPtGUqMzMrHKNJohvNDUKMzNrO432SX2zpFHAe3PRXRHxVE/TmJnZwNZoa677A3cB+5GeW7hT0qebGZiZmVWr0VNMJ5Cepn4KQFIH8FvgsmYFZmZm1Wq0R7khdaeUlvZhWjMzG4AaPYK4RtK1wEV5+ABS0xhmZraW6jFBSFo/IpZHxLGS/gH4YK6aHhGXNz88MzOrSm9HELcDO0j6aUQcTC/Ne5uZ2dqjtwSxnqTPAB/IRxCraKTDIDMzG5h6SxBHAp8FRgL71NX12mGQmZkNXD0miIi4VdJtQFdEnNyimMzMrA000ljf64AfijMzG2QafZbhekmfyt2INkTSTyQ9Jen+Qtk0SYskzc6vvQt1X5M0X9JDkj7ah3UwM7MmaDRBfAG4FPiLpOclvSDp+V6mOZfyjoZOjYhJ+XUVgKSJwIHAtnma03I/FGZmVpGGEkREjIiIIRExNCI2ysMb9TLNLcAzDcYxGZiRn7l4lNQv9Y4NTmtmZk3QaGN9kvQ5Sd/Iw+Mkre4O/ChJc/IpqE1y2RhgYWGcrlxmZmYVafQU02nA+4HP5OEXgR+vxvJOB94CTAIWA9/r6wwkTZE0S9KsJUuWrEYIZmbWiEYTxE4RMRV4BSAilgHr9XVhEfFkRLyW74w6i5WnkRYB4wqjjs1lZfOYHhGdEdHZ0dHR1xDMzKxBjSaIV/NF44A3mvt+va8LkzS6MPhJoHaH00zgQEnrS9oK2JrU/4SZmVWk0dZcfwhcDrxZ0smk5yL+racJJF0E7AJsLqkLOBHYRdIkUqJZQLo7ioiYK+kS4AFgBTA1Il7r89qYmVm/abTL0Z9JugfYHRCwb0TM62Wag0qKz+5h/JMBP61tZtYmemvuexipPaa3AvcBZ0bEilYEZmZm1ertGsR5QCcpOewFfLfpEZmZWVvo7RTTxIh4J4Cks/GFYzOzQaO3I4hXa298asnMbHDp7Qji3YU2lwRskIcFRG/NbZiZ2cDVW38QbjDPzGyQavRBOTMzG2ScIMzMrJQThJmZlXKCMDOzUk4QZmZWygnCzMxKOUGYmVkpJwgzMyvlBGFmZqWcIMzMrFTTEoSkn0h6StL9hbJNJV0n6U/57ya5XJJ+KGm+pDmSdmhWXGZm1phmHkGcC+xZV3Y8cH1EbA1cn4ch9TWxdX5NAU5vYlxmZtaApiWIiLgFeKaueDKpEyLy330L5edHcgcwUtLoZsVmZma9a/U1iFERsTi/fwIYld+PARYWxuvKZWZmVpHKLlJHRADR1+kkTZE0S9KsJUuWNCEyMzOD1ieIJ2unjvLfp3L5ImBcYbyxueyvRMT0iOiMiM6Ojo6mBmtmNpj11qNcf5sJHAp8O/+9olB+lKQZwE7Ac4VTUdYq0zauOgKY9lzVEZhZ1rQEIekiYBdgc0ldwImkxHCJpCOAx4D98+hXAXsD84GXgcObFZeZmTWmaQkiIg7qpmr3knEDmNqsWMzMrO/8JLWZmZVygjAzs1JOEGZmVsoJwszMSjlBmJlZKScIMzMr5QRhZmalnCDMzKyUE4SZmZVygjAzs1JOEGZmVsoJwszMSjlBmJlZKScIMzMr5QRhZmalnCDMzKxUq7scBUDSAuAF4DVgRUR0StoUuBiYACwA9o+IZVXEZ2Zm1R5B7BoRkyKiMw8fD1wfEVsD1+dhMzOrSDudYpoMnJffnwfsW2EsZmaDXlUJIoDfSLpH0pRcNioiFuf3TwCjqgnNzMygomsQwAcjYpGkNwPXSXqwWBkRISnKJswJZQrA+PHjmx+pmdkgVckRREQsyn+fAi4HdgSelDQaIP99qptpp0dEZ0R0dnR0tCpkM7NBp+UJQtJwSSNq74G/B+4HZgKH5tEOBa5odWxmZrZSFaeYRgGXS6ot/8KIuEbS3cAlko4AHgP2ryA2MzPLWp4gIuIR4N0l5UuB3Vsdj5mZlWun21zNzKyNOEGYmVkpJwgzMyvlBGFmZqWcIMzMrJQThJmZlXKCMDOzUk4QZmZWygnCzMxKOUGYmVkpJwgzMyvlBGFmZqWcIMzMrJQThJmZlXKCMDOzUk4QZmZWqu0ShKQ9JT0kab6k46uOx8xssGqrBCFpHeDHwF7AROAgSROrjcrMbHBqqwQB7AjMj4hHIuIvwAxgcsUxmZkNSu2WIMYACwvDXbnMzMxabN2qA+grSVOAKQDjx49f4/lNOP7XazyPNbXg2x+rOoRk2nNVR9A27jv0vqpDaBvveHBe1SG0jaln7FZ1CC3VbgliETCuMDw2l70hIqYD0wE6OztjTRfYNjtnM7M2026nmO4Gtpa0laT1gAOBmRXHZGY2KLXVEURErJB0FHAtsA7wk4iYW3FYZmaDUlslCICYiNAeAAAEfElEQVSIuAq4quo4zMwGu3Y7xWRmZm3CCcLMzEo5QZiZWSknCDMzK+UEYWZmpRSxxs+aVUbSEuCxquMANgeerjqINuFtsZK3xUreFiu1w7bYMiI6ehtpQCeIdiFpVkR0Vh1HO/C2WMnbYiVvi5UG0rbwKSYzMyvlBGFmZqWcIPrH9KoDaCPeFit5W6zkbbHSgNkWvgZhZmalfARhZmalnCDMzKyUE4SZmZVqu+a+252kbYDJrOwrexEwMyLcL+Mglj8XY4A7I+LFQvmeEXFNdZG1nqQdgYiIuyVNBPYEHsxN+Q9qks6PiEOqjqNRvkjdB5KOAw4CZgBduXgsqee7GRHx7apiazeSDo+Ic6qOoxUkfQmYCswDJgFHR8QVue7eiNihyvhaSdKJwF6kH5/XATsBNwJ7ANdGxMkVhtdSkup7wxSwK3ADQER8ouVB9ZETRB9I+iOwbUS8Wle+HjA3IrauJrL2I+n/ImJ81XG0gqT7gPdHxIuSJgCXAT+NiB9I+n1EbF9pgC2Ut8UkYH3gCWBsRDwvaQPS0dW7Kg2whSTdCzwA/A8QpARxEekHJRFxc3XRNcanmPrmdWAL/rr9p9G5blCRNKe7KmBUK2Op2JDaaaWIWCBpF+AySVuStsVgsiIiXgNelvRwRDwPEBF/ljTYviOdwNHACcCxETFb0p8HQmKocYLomy8D10v6E7Awl40H3gocVVlU1RkFfBRYVlcu4LbWh1OZJyVNiojZAPlI4uPAT4B3Vhtay/1F0psi4mXgPbVCSRszyH5ERcTrwKmSLs1/n2SA7XMHVLBVi4hrJL0N2JFVL1LfnX81DTZXAhvWdoxFkm5qfTiVOQRYUSyIiBXAIZLOrCakyvxdRCyHN3aQNUOBQ6sJqVoR0QXsJ+ljwPNVx9MXvgZhZmal/ByEmZmVcoIwM7NSThA2qEk6QdJcSXMkzZa0Uy5fIGnzkvFvy38nSLo/v99F0pUl45aWr2G8/T5Ps+74IrUNWpLeD3wc2CEilueEsF5P00TEB1oSnFkb8BGEDWajgacLd908HRGPF0eQtIGkqyV9Pg+/WDKfXkkaLuknku6S9HtJk3P5HZK2LYx3k6TO7sY3ayUnCBvMfgOMk/RHSadJ+nBd/YbAr4CLIuKsNVzWCcANEbEjqbmF70gaDlwM7A8gaTQwOiJm9TC+Wcs4QdiglZ9+fg8wBVgCXCzpsMIoVwDnRMT5/bC4vweOlzQbuAkYRnrI8hLg03mc/UnNdPQ0vlnL+BqEDWr5AcebgJtyO0KHAufm6t8Be0q6MNb8gSEBn4qIh/6qQloq6V3AAcCRPY0vaTA1YWIV8xGEDVqS3i6p2MDiJFZtZ+vfSc2I/LgfFnct8EVJyssuNuB3MfBVYOOImNPA+GYt4QRhg9mGwHmSHsgND04EptWNczSwgaT/WsNlfYvU3MQcSXPzcM1lpBY+L2lwfLOWcFMbZmZWykcQZmZWygnCzMxKOUGYmVkpJwgzMyvlBGFmZqWcIMzMrJQThJmZlXKCMDOzUv8fFjjj9gvWsb8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xs  = []\n",
    "yts = []\n",
    "for xi in X['X'].unique():\n",
    "    xs.append(xi)\n",
    "    yts.append(quad(integrand, y_min, y_max, args=(xi,))[0])\n",
    "pd.DataFrame({'X': xs, 'Yt': yts}).sort_values('X').plot(kind='bar', x='X', y='Yt')\n",
    "pp.title('Expected performance of translated articles\\nby author skill'); \n",
    "pp.xlabel('Skill level'); \n",
    "pp.ylabel('Performance')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see it matches the population level plot above!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}