{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "# Ch11 Why Model?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "Data cannot speak for themselves, 所以我们需要简单非参数模型到各种参数模型。本文的内容如下：\n",
    "\n",
    "- 教材对应章节的内容摘要\n",
    "- 教材中的 Programs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false",
    "toc-hr-collapsed": true
   },
   "source": [
    "## 内容摘要"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "本章描述了本书 Part I 中使用的非参数估计量与 Part II 中使用的基于模型的参数估计量之间的区别。它还回顾了 smoothing 的概念，并简要回顾了任何建模决策中涉及的 bias-variance trade-off。本章motivates the need for models in data analysis，无论分析目标是因果推理还是预测。We will take a break from causal considerations until the next chapter. 请记住，有关建模的统计文献非常丰富；本章只能重点介绍一些关键问题。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "**Big Picturce:**\n",
    "\n",
    "非参数模型 --> parametric models --> structural mean models\n",
    "\n",
    "- Models for the marginal mean of a counterfactual outcome are referred to as **marginal structural mean models**.\n",
    "- Effect modification and marginal structural models 在上面的模型增加新的变量 for effect modification\n",
    "- Structural nested mean model 会有一个对比某个基准 potential outcome.\n",
    "\n",
    "In Part II of this book we have described two different types of models for causal inference: \n",
    "\n",
    "- propensity models and \n",
    "- structural models. \n",
    "\n",
    "Let us now compare them. \n",
    "\n",
    "- Propensity models are models for the probability of treatment $A$ given the variables $L$ used to try to achieve conditional exchangeability. \n",
    "- Structural models describe the relation between the treatment $A$ and some component of the distribution (e.g., the mean) of the counterfactual outcome"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "Table of contents\n",
    "\n",
    "- Data cannot speak for themselves\n",
    "- Parametric estimators of the conditional mean\n",
    "- Nonparametric estimators of the conditional mean\n",
    "- Smoothing\n",
    "- The bias-variance trade-off"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### 11.1 Data cannot speak for themselves\n",
    "\n",
    "考虑一个研究总体，其中有16个人感染了人类免疫缺陷病毒（HIV）。与本书第一部分不同，我们不会将这些人视为 representatives of 1 billion individuals identical to them. Rather, these are just 16 individuals randomly sampled from a large, possibly hypothetical super-population: the target population.\n",
    "\n",
    "在研究开始时，每个人都会接受一定程度的治疗 $A$ (antiretroviral therapy), which is maintained during the study. At the end of the study, a continuous outcome $Y$ (CD4 cell count, in cells/mm<sup>3</sup>) is measured in all individuals. 我们希望一致的估算\n",
    " the mean of $Y$ among individuals with treatment level $A = a$ in the population from which the 16 individuals were randomly sampled. 也就是说，估计是未知的总体参数 $E[Y^{a=1} - Y^{a=0}|A=a]$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### 11.2 Parametric estimators of the conditional mean\n",
    "\n",
    "\n",
    "当使用参数模型时，只有在模型中编码的限制正确（i.e. if the model is correctly specified）的情况下，推断才是正确的。因此，基于模型的因果推断（本书其余大部分内容都基于此因果推断）依赖于(approximately)没有模型错误指定的条件。由于很少（甚至从来没有）完美地指定参数模型，因此几乎总是会期望一定程度的模型指定错误。This can be at least partially rectified by using nonparametric estimators，我们将在下一部分中进行描述。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### 11.3 Nonparametric estimators of the conditional mean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "- In this book we define “model” as an a priori mathematical restriction on the possible states of nature (Robins, Greenland 1986). Part I was entitled “Causal inference without models” because it only described **saturated models**.\n",
    "- 本书 Part II 中描述的大多数因果推断方法都依赖于 estimators that are parametric estimators of some part of the distribution of the data. Parametric estimation and other approaches to borrow information 是我们唯一的希望，因为多数情况下数据不能自己说话。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### 11.4 Smoothing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### 11.5 The bias-variance trade-off"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "这种偏差-方差分析是数据分析的核心。Investigators using models need to decide whether some protection against bias–by, say, adding more parameters to the model–is worth the cost in terms of variance. 尽管存在一些有助于这些决定的正式程序，但实际上，许多研究人员会根据诸如传统，参数的可解释性和软件可用性等标准来决定使用哪种模型。在本书中，我们通常会假设正确指定了我们的参数模型。这是一个不现实的假设，但它使我们能够专注于因果分析所特有的问题。毕竟，Model misspecification 是在任何类型的数据分析中都可能出现的问题，无论估计值是否具有因果关系。在实践中，谨慎的研究人员将始终对模型的有效性提出质疑，并将 conduct an analysis to assess the sensitivity of their estimates to model specification."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "## Setup and imports\n",
    "import warnings\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import scipy.stats\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "warnings.filterwarnings('ignore')\n",
    "%matplotlib inline\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false",
    "toc-hr-collapsed": true
   },
   "source": [
    "## Programs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### Program 11.1\n",
    "\n",
    "Half of the individuals were treated $(A = 1)$. Figure 11.1 is a scatter plot that displays each of the 16 individuals as a dot. The height of the dot indicates the value of the individual’s outcome $Y$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "df1 = pd.DataFrame({\n",
    "    'A': [1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],\n",
    "    'Y': [\n",
    "        200, 150, 220, 110, 50, 180, 90, 170,\n",
    "        170, 30, 70, 110, 80, 50, 10, 20\n",
    "    ]\n",
    "})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfgAAAF3CAYAAABJzllyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAWYElEQVR4nO3dcYzfd33f8de7jtvdVlZDYxg2aQNbsJqBhtEJZYrUUVHVTf4gHoIKpEKK0mZiMHVrd2q8/UHVbiLqtauG1tKlKiKgQqGda6KW6dYFNqZqob3MXRxAJzIKSc4RcQuXMnHtHPPZH/dzsF3nd3c+3/f3+33yeEinu/vcz/d764vNM7/vfe/7qdZaAIC+fNukBwAArj6BB4AOCTwAdEjgAaBDAg8AHRJ4AOjQNZMe4Gq69tpr2/XXXz/pMQBgMA8++OCftdb2X7reVeCvv/76LC8vT3oMABhMVX35cutO0QNAhwQeADok8ADQIYEHgA4JPAB0SOABoEMCDwAdEngA6JDAA0CHurqTHQBMkxMnV7O4tJLTa+s5sG8uC0cO5ejhg4M8t8ADwC44cXI1x46fyvrZc0mS1bX1HDt+KkkGibxT9ACwCxaXVp6J+3nrZ89lcWllkOcXeADYBafX1re1frUJPADsggP75ra1frUJPADsgoUjhzK3d89Fa3N792ThyKFBnt9FdgCwC85fSOcqegDozNHDBwcL+qWcogeADgk8AHRI4AGgQwIPAB0SeADokMADQIcEHgA6JPAA0CGBB4AOCTwAdEjgAaBDAg8AHRJ4AOiQwANAhwQeADok8ADQIYEHgA4JPAB0SOABoEMCDwAdGizwVXVdVX2qqj5fVZ+tqp8crb+gqv6gqr4wev/80XpV1Xur6pGqeqiqXj3UrAAw64Z8Bf90kp9urX1fkpuSvLOqbkxyV5L7W2s3JLl/9HmS3JLkhtHbnUneN+CsALBjJ06u5ua7P5mX3vX7ufnuT+bEydXBnnuwwLfWnmit/a/Rx19P8vkkB5PcluTe0cPuTXJ09PFtST7YNjyQZF9VvXioeQFgJ06cXM2x46eyuraelmR1bT3Hjp8aLPIT+Rl8VV2f5HCSzyR5UWvtiWTjPwKSvHD0sINJHrvgjz0+WgOAqbe4tJL1s+cuWls/ey6LSyuDPP/gga+q70zyn5L889baX4x76GXW2mW+351VtVxVy2fOnLlaYwLAjpxeW9/W+tU2aOCram824v6brbXjo+WvnD/1Pnr/5Gj98STXXfDHX5Lk9KXfs7V2T2ttvrU2v3///t0bHgC24cC+uW2tX21DXkVfSX4jyedba//ugi/dl+T20ce3J/n4BetvG11Nf1OSp86fygeAabdw5FDm9u65aG1u754sHDk0yPNfM8izbLg5yVuTnKqqPxmt/askdyf5WFXdkeTRJG8afe0TSW5N8kiSbyR5+4CzAsCOHD28cdnY4tJKTq+t58C+uSwcOfTM+m6r1v7aj7Vn1vz8fFteXp70GAAwmKp6sLU2f+m6O9kBQIcEHgA6JPAA0CGBB4AOCTwAdEjgAaBDAg8AHRJ4AOiQwANAhwQeADok8ADQIYEHgA4JPAB0SOABoEMCDwAdEngA6JDAA0CHBB4AOiTwANAhgQeADl0z6QEAoFcnTq5mcWklp9fWc2DfXBaOHMrRwwcHeW6BB4BdcOLkao4dP5X1s+eSJKtr6zl2/FSSDBJ5p+gBYBcsLq08E/fz1s+ey+LSyiDPL/AAsAtOr61va/1qE3gA2AUH9s1ta/1qE3gA2AULRw5lbu+ei9bm9u7JwpFDgzy/i+wAYBecv5DOVfQA0Jmjhw8OFvRLOUUPAB0SeADokMADQIcEHgA6JPAA0CGBB4AO+TW5TUxyJyAAZpvd5KbUpHcCAmB2TbohTtGPMemdgACYXZNuiMCPMemdgACYXZNuiMCPMemdgACYXZNuiMCPMemdgACYXZNuiIvsxpj0TkAAzK5JN6Raa4M80RDm5+fb8vLypMcAgMFU1YOttflL152iB4AOCTwAdEjgAaBDAg8AHRJ4AOiQwANAhwQeADok8ADQIYEHgA4JPAB0SOABoEMCDwAdEngA6JDAA0CHBB4AOiTwANAhgQeADgk8AHRI4AGgQwIPAB0aLPBV9f6qerKqHr5g7WerarWq/mT0dusFXztWVY9U1UpVHRlqTgC4Wk6cXM3Nd38yL73r93Pz3Z/MiZOrgz33NYM9U/KBJP8hyQcvWf/l1tovXrhQVTcmeXOSv5/kQJL/WlUvb62dG2JQANipEydXc+z4qayf3UjX6tp6jh0/lSQ5evjgrj//YK/gW2ufTvLVLT78tiS/1Vr7q9banyZ5JMlrdm04ALjKFpdWnon7eetnz2VxaWWQ55+Gn8G/q6oeGp3Cf/5o7WCSxy54zOOjtb+mqu6squWqWj5z5sxuzwoAW3J6bX1b61fbpAP/viR/N8mrkjyR5JdG63WZx7bLfYPW2j2ttfnW2vz+/ft3Z0oA2KYD++a2tX61TTTwrbWvtNbOtda+meTX863T8I8nue6Ch74kyemh5wOAK7Vw5FDm9u65aG1u754sHDk0yPNPNPBV9eILPv3HSc5fYX9fkjdX1XdU1UuT3JDkj4aeDwCu1NHDB/OeN7wyB/fNpZIc3DeX97zhlYNcYJcMeBV9VX0kyWuTXFtVjyd5d5LXVtWrsnH6/UtJ/kmStNY+W1UfS/K5JE8neacr6AGYNUcPHxws6Jeq1i77o+2ZND8/35aXlyc9BgAMpqoebK3NX7o+6YvsAIBdIPAA0CGBB4AOCTwAdEjgAaBDAg8AHRJ4AOiQwANAhwQeADok8ADQIYEHgA4JPAB0SOABoEMCDwAdEngA6JDAA0CHBB4AOiTwANAhgQeADgk8AHRI4AGgQwIPAB0SeADokMADQIcEHgA6JPAA0CGBB4AOCTwAdEjgAaBDAg8AHRJ4AOiQwANAhwQeADok8ADQIYEHgA4JPAB0SOABoEMCDwAdEngA6JDAA0CHBB4AOjQ28FV101CDAABXz2av4D9dVT9fVdcMMg0AcFVsFvhbkrw1yR9V1Y0DzAMAXAVjA99auz/JK5OcTLJcVT81yFQAwI5sepFda+3rrbU7svFK/heq6v9W1V9c+Lb7YwIA27Gln61X1XySf5PkC0l+McnTuzkUALAzYwM/urju3Ul+JsmvJrmrtfaXQwwGAFy5zV7B/3GSFyS5ZfTz+OecEydXs7i0ktNr6zmwby4LRw7l6OGDkx4LgBkwyYZsFviHk7yrtfbUEMNMmxMnV3Ps+Kmsnz2XJFldW8+x46eSROQBGGvSDdnsKvq3PlfjniSLSyvP/A9z3vrZc1lcWpnQRADMikk3xK1qxzi9tr6tdQA4b9INEfgxDuyb29Y6AJw36YYI/BgLRw5lbu+ei9bm9u7JwpFDE5oIgFkx6Ya4x/wY5y+CcBU9ANs16YZUa22QJxrC/Px8W15envQYADCYqnqwtTZ/6bpT9ADQIYEHgA4JPAB0SOABoEMCDwAdGizwVfX+qnqyqh6+YO0FVfUHVfWF0fvnj9arqt5bVY9U1UNV9eqh5gSAHgz5Cv4DSX74krW7ktzfWrshyf2jz5PkliQ3jN7uTPK+gWYEgC4MFvjW2qeTfPWS5duS3Dv6+N4kRy9Y/2Db8ECSfVX14mEmBYDZN+mfwb+otfZEkozev3C0fjDJYxc87vHRGgCwBZMO/LOpy6xd9pZ7VXVnVS1X1fKZM2d2eSwAmA2TDvxXzp96H71/crT+eJLrLnjcS5Kcvtw3aK3d01qbb63N79+/f1eHBYBZMenA35fk9tHHtyf5+AXrbxtdTX9TkqfOn8oHADY32G5yVfWRJK9Ncm1VPZ7k3UnuTvKxqrojyaNJ3jR6+CeS3JrkkSTfSPL2oeYEgB4MFvjW2lue5Uuvu8xjW5J37u5EANCvSZ+iBwB2gcADQIcEHgA6JPAA0CGBB4AOCTwAdGiwX5MDgOeaEydXs7i0ktNr6zmwby4LRw7l6OFhtlYReADYBSdOrubY8VNZP3suSbK6tp5jx08lySCRd4oeAHbB4tLKM3E/b/3suSwurQzy/AIPALvg9Nr6ttavNoEHgF1wYN/cttavNoEHgF2wcORQ5vbuuWhtbu+eLBw5NMjzu8gOAHbB+QvpXEUPAJ05evjgYEG/lFP0ANAhgQeADgk8AHRI4AGgQwIPAB0SeADokF+T28QkdwICgCsl8GNMeicgALhSTtGPMemdgADgSgn8GJPeCQgArpTAjzHpnYAA4EoJ/BiT3gkIAK6Ui+zGmPROQABwpQR+E5PcCQgArpRT9ADQIYEHgA4JPAB0SOABoEMCDwAdEngA6JBfk9uE3eQAmEUCP4bd5ACYVU7Rj2E3OQBmlcCPYTc5AGaVwI9hNzkAZpXAj2E3OQBmlYvsxrCbHACzSuA3YTc5AGaRU/QA0CGBB4AOCTwAdEjgAaBDAg8AHRJ4AOiQwANAhwQeADok8ADQIYEHgA4JPAB0SOABoEMCDwAdEngA6JDAA0CHBB4AOiTwANAhgQeADgk8AHRI4AGgQ9dMeoAkqaovJfl6knNJnm6tzVfVC5J8NMn1Sb6U5Edaa18berYTJ1ezuLSS02vrObBvLgtHDuXo4YNDjwHADJpkQ6bpFfwPtNZe1VqbH31+V5L7W2s3JLl/9PmgTpxczbHjp7K6tp6WZHVtPceOn8qJk6tDjwLAjJl0Q6Yp8Je6Lcm9o4/vTXJ06AEWl1ayfvbcRWvrZ89lcWll6FEAmDGTbsi0BL4l+S9V9WBV3Tlae1Fr7YkkGb1/4eX+YFXdWVXLVbV85syZqzrU6bX1ba0DwHmTbsi0BP7m1tqrk9yS5J1V9f1b/YOttXtaa/Ottfn9+/df1aEO7Jvb1joAnDfphkxF4Ftrp0fvn0zyu0lek+QrVfXiJBm9f3LouRaOHMrc3j0Xrc3t3ZOFI4eGHgWAGTPphkw88FX1t6rqeec/TvJDSR5Ocl+S20cPuz3Jx4ee7ejhg3nPG16Zg/vmUkkO7pvLe97wSlfRA7CpSTekWmuDPNGzDlD1smy8ak82fm3vw621f1tV353kY0m+J8mjSd7UWvvquO81Pz/flpeXd3VeAJgmVfXgBb+B9oyJ/x58a+2LSf7BZdb/PMnrhp8IAGbfxE/RAwBXn8ADQIcEHgA6JPAA0CGBB4AOCTwAdEjgAaBDAg8AHRJ4AOiQwANAhwQeADok8ADQIYEHgA4JPAB0SOABoEMCDwAdEngA6JDAA0CHBB4AOiTwANChayY9wLQ7cXI1i0srOb22ngP75rJw5FCOHj446bEAYCyBH+PEydUcO34q62fPJUlW19Zz7PipJBF5AKaaU/RjLC6tPBP389bPnsvi0sqEJgKArRH4MU6vrW9rHQCmhcCPcWDf3LbWAWBaCPwYC0cOZW7vnovW5vbuycKRQxOaCAC2xkV2Y5y/kM5V9ADMGoHfxNHDBwUdgJnjFD0AdEjgAaBDAg8AHRJ4AOiQwANAhwQeADrk1+Q2YTc5AGaRwI9hNzkAZpVT9GPYTQ6AWSXwY9hNDoBZJfBj2E0OgFkl8GPYTQ6AWeUiuzHsJgfArBL4TdhNDoBZ5BQ9AHRI4AGgQwIPAB0SeADokMADQIcEHgA65NfkNmE3OQBmkcCPYTc5AGaVU/Rj2E0OgFkl8GPYTQ6AWSXwY9hNDoBZJfBj2E0OgFnlIrsx7CYHwKwS+E3YTQ6AWeQUPQB0SOABoEMCDwAdEngA6JDAA0CHBB4AOjT1ga+qH66qlap6pKrumvQ8ADALpjrwVbUnya8kuSXJjUneUlU3TnYqAJh+Ux34JK9J8khr7Yuttf+X5LeS3DbhmQBg6k174A8meeyCzx8frT2jqu6squWqWj5z5sygwwHAtJr2W9XWZdbaRZ+0dk+Se5Kkqs5U1Zd3aZZrk/zZLn3v5wLHb2ccv51x/HbOMdyZ3Tx+33u5xWkP/ONJrrvg85ckOf1sD26t7d+tQapqubU2v1vfv3eO3844fjvj+O2cY7gzkzh+036K/o+T3FBVL62qb0/y5iT3TXgmAJh6U/0KvrX2dFW9K8lSkj1J3t9a++yExwKAqTfVgU+S1tonknxi0nNk9HN+rpjjtzOO3844fjvnGO7M4MevWmubPwoAmCnT/jN4AOAKCPwlNrs1blV9R1V9dPT1z1TV9cNPOb22cPx+qqo+V1UPVdX9VXXZX+94rtrqrZmr6o1V1arKVc0X2Mrxq6ofGf0d/GxVfXjoGafZFv79fk9VfaqqTo7+Dd86iTmnVVW9v6qerKqHn+XrVVXvHR3fh6rq1bs6UGvN2+gtGxfy/Z8kL0vy7Un+d5IbL3nMP03ya6OP35zko5Oee1retnj8fiDJ3xx9/A7Hb3vHb/S45yX5dJIHksxPeu5pedvi378bkpxM8vzR5y+c9NzT8rbF43dPkneMPr4xyZcmPfc0vSX5/iSvTvLws3z91iT/ORv3eLkpyWd2cx6v4C+2lVvj3pbk3tHHv5PkdVV1uRvyPBdtevxaa59qrX1j9OkD2bi3ARu2emvmn0/yC0n+csjhZsBWjt9PJPmV1trXkqS19uTAM06zrRy/luRvjz7+roy5L8lzUWvt00m+OuYhtyX5YNvwQJJ9VfXi3ZpH4C+26a1xL3xMa+3pJE8l+e5Bppt+Wzl+F7ojG/81y4at3Jr5cJLrWmu/N+RgM2Irf/9enuTlVfWHVfVAVf3wYNNNv60cv59N8qNV9Xg2frvpnw0zWje2+/+ROzL1vyY3sE1vjbvFxzxXbfnYVNWPJplP8o92daLZMvb4VdW3JfnlJD821EAzZit//67Jxmn612bj7NH/qKpXtNbWdnm2WbCV4/eWJB9orf1SVf3DJB8aHb9v7v54XRi0H17BX2wrt8Z95jFVdU02TlONOyXzXLKlWwtX1Q8m+ddJXt9a+6uBZpsFmx2/5yV5RZL/VlVfysbP8O5zod0ztvrv9+OttbOttT9NspKN4LO143dHko8lSWvtfyb5G9m4xzpbs63br++UwF9sK7fGvS/J7aOP35jkk2109QSbH7/RKeb/mI24+/nnxcYev9baU621a1tr17fWrs/GNQyvb60tT2bcqbOVf78nsnGhZ6rq2mycsv/ioFNOr60cv0eTvC5Jqur7shF423hu3X1J3ja6mv6mJE+11p7YrSdziv4C7VlujVtVP5dkubV2X5LfyMZpqUey8cr9zZObeLps8fgtJvnOJL89ujbx0dba6yc29BTZ4vHjWWzx+C0l+aGq+lySc0kWWmt/Prmpp8cWj99PJ/n1qvoX2Ti1/GNe4HxLVX0kGz/+uXZ0ncK7k+xNktbar2XjuoVbkzyS5BtJ3r6r8/jfBgD64xQ9AHRI4AGgQwIPAB0SeADokMADQIcEHgA6JPDAjlTV4ao6V1V/OOlZgG8ReGCnfiLJryZ5xejuZsAUcKMb4IpV1VySJ7KxD/ZPJvlaa+1fTnYqIPEKHtiZNyb5cmvtoSQfysZ9tvdOeCYgAg/szI9nI+xJ8t+zcX9tewvAFBB44IpU1d9LcnOSDyfJaNOR38xG9IEJs5sccKV+PBu7jj062hkwSSpJquq61tpjkxoMcJEdcAWq6pokjyX590l+75IvfyjJ77bWfm7wwYBnCDywbVV1W5LfSfJ3Lt1Pvap+Jsk7krystfbNScwH+Bk8cGXuSPKpS+M+8ttJvjfJDw47EnAhr+ABoENewQNAhwQeADok8ADQIYEHgA4JPAB0SOABoEMCDwAdEngA6JDAA0CH/j+E0TOwhO+vKgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "\n",
    "ax.scatter(df1['A'], df1['Y'])\n",
    "ax.set_xlabel('A', fontsize=14)\n",
    "ax.set_ylabel('Y', fontsize=14);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "Our estimate of the population mean in the treated is the sample aver- age 146.25 for those with $A = 1$, and our estimate of the population mean in the untreated is the sample average 67.50 in those with $A = 0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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 tr th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe thead tr:last-of-type th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th colspan=\"8\" halign=\"left\">Y</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th>count</th>\n",
       "      <th>mean</th>\n",
       "      <th>std</th>\n",
       "      <th>min</th>\n",
       "      <th>25%</th>\n",
       "      <th>50%</th>\n",
       "      <th>75%</th>\n",
       "      <th>max</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>0</td>\n",
       "      <td>8.0</td>\n",
       "      <td>67.50</td>\n",
       "      <td>53.117121</td>\n",
       "      <td>10.0</td>\n",
       "      <td>27.5</td>\n",
       "      <td>60.0</td>\n",
       "      <td>87.5</td>\n",
       "      <td>170.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>8.0</td>\n",
       "      <td>146.25</td>\n",
       "      <td>58.294205</td>\n",
       "      <td>50.0</td>\n",
       "      <td>105.0</td>\n",
       "      <td>160.0</td>\n",
       "      <td>185.0</td>\n",
       "      <td>220.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      Y                                                     \n",
       "  count    mean        std   min    25%    50%    75%    max\n",
       "A                                                           \n",
       "0   8.0   67.50  53.117121  10.0   27.5   60.0   87.5  170.0\n",
       "1   8.0  146.25  58.294205  50.0  105.0  160.0  185.0  220.0"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df1.groupby('A').describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "\"Now suppose treatment $A$ is a polytomous variable that can take 4 possible values\". Now suppose treatment $A$ is a polytomous variable that can take 4 possible values: no treatment $(A = 1)$, low-dose treatment $(A = 2)$, medium-dose treatment $(A = 3)$, and high-dose treatment $(A = 4)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "df2 = pd.DataFrame({\n",
    "    'A': [\n",
    "        1, 1, 1, 1, 2, 2, 2, 2,\n",
    "        3, 3, 3, 3, 4, 4, 4, 4\n",
    "    ],\n",
    "    'Y': [\n",
    "        110, 80, 50, 40, 170, 30, 70, 50,\n",
    "        110, 50, 180, 130, 200, 150, 220, 210\n",
    "    ]\n",
    "})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfgAAAF3CAYAAABJzllyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAbKElEQVR4nO3df3Bld3nf8fdTWaVqIRXEArxrbxaI0QAu3QWN444H6gKtMKX20kJiT2oMMSzOwBQmVMVyZkJK+weN+JFQJmSW4rFNjWOChXCpqXAMwaVTQ7SWY5ksCjY1eLU79mIj24xv3V356R862tXK0lralc6596v3a+aO7n3OubqPj7+rj865X91vZCaSJKksf6vpBiRJ0voz4CVJKpABL0lSgQx4SZIKZMBLklQgA16SpAKd1nQD6+n000/P7du3N92GJEm12bt3788ys29pvaiA3759OxMTE023IUlSbSLiJ8vVvUQvSVKBagv4iDgrIr4dEfsi4gcR8cGqPhIRP4yIeyLiqxHRW9W3R0QrIu6ubn9SV6+SJHW6Os/gjwAfzsxXAOcB74+IVwK3Aedk5quBvwGGFz3n/szcUd2urLFXSZI6Wm0Bn5kHM/Ou6v4TwD5ga2Z+MzOPVLvdCZxZV0+SJJWqkffgI2I7sBP43pJNvwV8Y9Hjl0TEZER8JyJeV1N7kiR1vNpn0UfEc4GbgQ9l5uOL6r/L/GX8G6rSQWBbZj4SEa8FxiLiVYufUz1vN7AbYNu2bXX8J0iS1PZqPYOPiG7mw/2GzBxdVL8ceCvwm1mtX5uZT2XmI9X9vcD9wMuXfs/M3JOZA5k50Nf3jD8DlCRpU6pzFn0AXwD2ZeanFtXfDHwEuCgzn1xU74uIrur+S4GzgR/X1a8kSZ2szkv05wOXAVMRcXdVuxr4DPAc4Lb53wG4s5ox/3rgYxFxBJgDrszMR2vsV5KkjlVbwGfmd4FYZtOtK+x/M/OX8yVJ6mhjkzOMjE9zYLbFlt4ehgb72bVz64a+ZlEfVStJUrsZm5xheHSK1uE5AGZmWwyPTgFsaMj7UbWSJG2gkfHpo+G+oHV4jpHx6Q19XQNekqQNdGC2tab6ejHgJUnaQFt6e9ZUXy8GvCRJG2hosJ+e7q7jaj3dXQwN9m/o6zrJTpKkDbQwkc5Z9JIkFWbXzq0bHuhLeYlekqQCGfCSJBXIgJckqUAGvCRJBTLgJUkqkAEvSVKBDHhJkgpkwEuSVCA/6EaSpA3mevCSJBXG9eAlSSqQ68FLklQg14OXJKlArgcvSVKBXA9ekqQCuR68JEmFcj14SZK0Lgx4SZIKZMBLklQgA16SpAIZ8JIkFciAlySpQAa8JEkFqi3gI+KsiPh2ROyLiB9ExAer+gsi4raI+FH19flVPSLiMxFxX0TcExGvqatXSZI6XZ1n8EeAD2fmK4DzgPdHxCuBq4DbM/Ns4PbqMcCFwNnVbTfwuRp7lSRp3YxNznD+x7/FS67675z/8W8xNjmz4a9ZW8Bn5sHMvKu6/wSwD9gKXAxcV+12HbCrun8xcH3OuxPojYgz6upXkqT1sLAe/Mxsi+TYevAbHfKNvAcfEduBncD3gBdl5kGY/yUAeGG121bgwUVP21/VJEnqGJtmPfiIeC5wM/ChzHz8RLsuU8tlvt/uiJiIiIlDhw6tV5uSJK2LTbEefER0Mx/uN2TmaFV+aOHSe/X14aq+Hzhr0dPPBA4s/Z6ZuSczBzJzoK+vb+OalyTpJBS/HnxEBPAFYF9mfmrRpluAy6v7lwNfW1R/ZzWb/jzgsYVL+ZIkdYrNsB78+cBlwFRE3F3VrgY+Dnw5Iq4Afgq8o9p2K/AW4D7gSeDdNfYqSdK6aGo9+Mh8xtvaHWtgYCAnJiaabkOSpNpExN7MHFha95PsJEkqkAEvSVKBDHhJkgpkwEuSVCADXpKkAhnwkiQVyICXJKlABrwkSQUy4CVJKpABL0lSgQx4SZIKZMBLklQgA16SpAIZ8JIkFciAlySpQAa8JEkFMuAlSSqQAS9JUoEMeEmSCmTAS5JUIANekqQCGfCSJBXIgJckqUCnNd2AJJVibHKGkfFpDsy22NLbw9BgP7t2bm26LW1SBrwkrYOxyRmGR6doHZ4DYGa2xfDoFIAhr0Z4iV6S1sHI+PTRcF/QOjzHyPh0Qx1pszPgJWkdHJhtrakubTQDXpLWwZbenjXVpY1mwEvSOhga7Kenu+u4Wk93F0OD/Q11pM3OSXaStA4WJtI5i17toraAj4hrgLcCD2fmOVXtJmDh19teYDYzd0TEdmAfsDA75c7MvLKuXiXpZOzaudVAV9uo8wz+WuCzwPULhcz8jYX7EfFJ4LFF+9+fmTtq606SpILUFvCZeUd1Zv4MERHArwNvqKsfSZJK1i6T7F4HPJSZP1pUe0lETEbEdyLidU01JklSJ2qXSXaXAjcuenwQ2JaZj0TEa4GxiHhVZj6+9IkRsRvYDbBt27ZampUkqd01fgYfEacB/xK4aaGWmU9l5iPV/b3A/cDLl3t+Zu7JzIHMHOjr66ujZUmS2l7jAQ+8CfhhZu5fKEREX0R0VfdfCpwN/Lih/iRJ6ji1BXxE3Aj8b6A/IvZHxBXVpks4/vI8wOuBeyLir4CvAFdm5qN19SpJUqercxb9pSvU37VM7Wbg5o3uSZKkUrXLJDupY7jmt6ROYMBLa+Ca35I6RTtMspM6hmt+S+oUBry0Bq75LalTGPDSGrjmt6ROYcBLa+Ca35I6hZPspDVwzW9JncKAl9bINb8ldQIv0UuSVCADXpKkAhnwkiQVyICXJKlABrwkSQUy4CVJKpABL0lSgQx4SZIKZMBLklQgA16SpAIZ8JIkFciAlySpQAa8JEkFMuAlSSqQAS9JUoEMeEmSCmTAS5JUIANekqQCGfCSJBXIgJckqUAGvCRJBTLgJUkqUG0BHxHXRMTDEXHvotrvR8RMRNxd3d6yaNtwRNwXEdMRMVhXn5IkleC0Gl/rWuCzwPVL6p/OzE8sLkTEK4FLgFcBW4A/j4iXZ+ZcHY1KkrSexiZnGBmf5sBsiy29PQwN9rNr59YNfc3azuAz8w7g0VXufjHwp5n5VGb+H+A+4NwNa06SpA0yNjnD8OgUM7MtEpiZbTE8OsXY5MyGvm47vAf/gYi4p7qE//yqthV4cNE++6uaJEkdZWR8mtbh4y9Atw7PMTI+vaGv23TAfw54GbADOAh8sqrHMvvmct8gInZHxERETBw6dGhjupQk6SQdmG2tqb5eGg34zHwoM+cy82ng8xy7DL8fOGvRrmcCB1b4HnsycyAzB/r6+ja2YUmS1mhLb8+a6uul0YCPiDMWPXwbsDDD/hbgkoh4TkS8BDgb+H7d/UmSdKqGBvvp6e46rtbT3cXQYP+Gvm5ts+gj4kbgAuD0iNgPfBS4ICJ2MH/5/QHgfQCZ+YOI+DLw18AR4P3OoJckdaKF2fJ1z6KPzGXf2u5IAwMDOTEx0XQbkiTVJiL2ZubA0nrTk+wkSdIGMOAlSSqQAS9JUoEMeEmSCmTAS5JUIANekqQCGfCSJBXIgJckqUAGvCRJBTLgJUkqkAEvSVKBDHhJkgpkwEuSVCADXpKkAhnwkiQVyICXJKlABrwkSQUy4CVJKpABL0lSgQx4SZIKZMBLklQgA16SpAIZ8JIkFei0phuQpFKMTc4wMj7NgdkWW3p7GBrsZ9fOrU23pU3KgJekdTA2OcPw6BStw3MAzMy2GB6dAjDk1Qgv0UvSOhgZnz4a7gtah+cYGZ9uqCNtdga8JK2DA7OtNdWljWbAS9I62NLbs6a6tNEMeElaB0OD/fR0dx1X6+nuYmiwv6GOtNk5yU6S1sHCRDpn0atd1BbwEXEN8Fbg4cw8p6qNAP8C+H/A/cC7M3M2IrYD+4CF2Sl3ZuaVdfUqSSdj186tBrraxgkv0UfEeev4WtcCb15Suw04JzNfDfwNMLxo2/2ZuaO6Ge6SJK3Bs70Hf0dE/IeIOOUz/cy8A3h0Se2bmXmkengncOapvo4kSXr2gL8QuAz4fkS8coN7+S3gG4sevyQiJiPiOxHxug1+bUmSinLCgM/M24F/AEwCExHxOxvRRET8LnAEuKEqHQS2ZeZO4HeAL0XEL63w3N0RMRERE4cOHdqI9iRJ6jjP+mdymflEZl7B/Jn8H0TELyLi8cW3U2kgIi5nfvLdb2ZmVq/5VGY+Ut3fy/wEvJev0N+ezBzIzIG+vr5TaUWSpGKs6r31iBgA/iPwI+ATzJ9tn7KIeDPwEeAfZ+aTi+p9wKOZORcRLwXOBn68Hq8pSdJmcMKArybXfZT5EP5j4KrM/L8n80IRcSNwAXB6ROyvvu8w8BzgtoiAY38O93rgYxFxBJgDrszMR5f9xpIk6Rme7Qz+L4EXABdW78eftMy8dJnyF1bY92bg5lN5PUmSNrNnew/+XuDVpxrukiSpXic8g8/My+pqRJIkrR8Xm5EkqUAGvCRJBTLgJUkqkAEvSVKBDHhJkgpkwEuSVCADXpKkAhnwkiQVyICXJKlABrwkSQUy4CVJKpABL0lSgQx4SZIK9GzrwW9qY5MzjIxPc2C2xZbeHoYG+9m1c2vTbUlqU/7MUDsx4FcwNjnD8OgUrcNzAMzMthgenQLwH6ykZ/BnhtqNl+hXMDI+ffQf6oLW4TlGxqcb6khSO/NnhtqNAb+CA7OtNdUlbW7+zFC7MeBXsKW3Z011SZubPzPUbgz4FQwN9tPT3XVcrae7i6HB/oY6ktTO/JmhduMkuxUsTIpxRqyk1fBnhtpNZGbTPaybgYGBnJiYaLoNSZJqExF7M3Ngad1L9JIkFciAlySpQAa8JEkFMuAlSSqQAS9JUoEMeEmSCmTAS5JUoFoDPiKuiYiHI+LeRbUXRMRtEfGj6uvzq3pExGci4r6IuCciXlNnr5IkdbK6z+CvBd68pHYVcHtmng3cXj0GuBA4u7rtBj5XU4+SJHW8WgM+M+8AHl1Svhi4rrp/HbBrUf36nHcn0BsRZ9TTqSRJna0d3oN/UWYeBKi+vrCqbwUeXLTf/qomSZKeRTsE/EpimdozPjg/InZHxERETBw6dKiGtiRJan/tEPAPLVx6r74+XNX3A2ct2u9M4MDSJ2fmnswcyMyBvr6+DW9WkqRO0A4BfwtweXX/cuBri+rvrGbTnwc8tnApX5IknVit68FHxI3ABcDpEbEf+CjwceDLEXEF8FPgHdXutwJvAe4DngTeXWevkiR1sloDPjMvXWHTG5fZN4H3b2xHkiSVqR0u0UuSpHVmwEuSVCADXpKkAhnwkiQVyICXJKlABrwkSQUy4CVJKpABL0lSgQx4SZIKZMBLklQgA16SpAIZ8JIkFciAlySpQAa8JEkFMuAlSSqQAS9JUoEMeEmSCmTAS5JUIANekqQCGfCSJBXIgJckqUAGvCRJBTLgJUkq0GlNN9DOxiZnGBmf5sBsiy29PQwN9rNr59am25Ik6VkZ8CsYm5xheHSK1uE5AGZmWwyPTgEY8pKktucl+hWMjE8fDfcFrcNzjIxPN9SRJEmrZ8Cv4MBsa011SZLaiQG/gi29PWuqS5LUTgz4FQwN9tPT3XVcrae7i6HB/oY6kiRp9Zxkt4KFiXTOopckdaLGAz4i+oGbFpVeCvwe0Au8FzhU1a/OzFvr7G3Xzq0GuiSpIzUe8Jk5DewAiIguYAb4KvBu4NOZ+YkG25MkqSO123vwbwTuz8yfNN2IJEmdrN0C/hLgxkWPPxAR90TENRHx/KaakiSp07RNwEfE3wYuAv6sKn0OeBnzl+8PAp9c4Xm7I2IiIiYOHTq03C6SJG06bRPwwIXAXZn5EEBmPpSZc5n5NPB54NzlnpSZezJzIDMH+vr6amxXkqT21U4BfymLLs9HxBmLtr0NuLf2jiRJ6lCNz6IHiIi/C/xT4H2Lyn8QETuABB5Ysk2SJJ1AWwR8Zj4J/PKS2mUNtSNJUsdri4CXOsnY5IyfcCip7Rnw0hqMTc4wPDp1dCnhmdkWw6NTAIa8pLbSTpPspLY3Mj59NNwXtA7PMTI+3VBHkrQ8A15agwOzrTXVJakpBry0Blt6e9ZUl6SmGPDSGgwN9tPT3XVcrae7i6HB/oY6kqTlOclOWoOFiXTOopfU7gx4aY127dxqoEtqe16ilySpQAa8JEkFMuAlSSqQAS9JUoEMeEmSCmTAS5JUIANekqQCGfCSJBXIgJckqUAGvCRJBTLgJUkqkAEvSVKBDHhJkgpkwEuSVCADXpKkAhnwkiQVyICXJKlABrwkSQUy4CVJKpABL0lSgQx4SZIKZMBLklSg05puYEFEPAA8AcwBRzJzICJeANwEbAceAH49M3/eVI+SJHWKtgn4yj/JzJ8tenwVcHtmfjwirqoef6SuZsYmZxgZn+bAbIstvT0MDfaza+fWul5ebcpxoZU4NtRO2i3gl7oYuKC6fx3wF9QU8GOTMwyPTtE6PAfAzGyL4dEpAP/BbmKOC63EsaF2007vwSfwzYjYGxG7q9qLMvMgQPX1hXU1MzI+ffQf6oLW4TlGxqfrakFtyHGhlTg21G7a6Qz+/Mw8EBEvBG6LiB+u5knVLwO7AbZt27ZuzRyYba2prs3BcaGVODbUbtrmDD4zD1RfHwa+CpwLPBQRZwBUXx9e5nl7MnMgMwf6+vrWrZ8tvT1rqmtzcFxoJY4NtZu2CPiI+HsR8byF+8A/A+4FbgEur3a7HPhaXT0NDfbT0911XK2nu4uhwf66WlAbclxoJY4NtZt2uUT/IuCrEQHzPX0pM/9HRPwl8OWIuAL4KfCOuhpamBTjjFgt5rjQShwbajeRmU33sG4GBgZyYmKi6TYkSapNROzNzIGl9ba4RC9JktaXAS9JUoEMeEmSCmTAS5JUIANekqQCGfCSJBXIgJckqUDt8kE3bcmlHyVJncqAX4FLP0qSOpmX6Ffg0o+SpE5mwK/ApR8lSZ3MgF+BSz9KkjqZAb8Cl36UJHUyJ9mtwKUfJUmdzIA/gV07txrokqSO5CV6SZIKZMBLklQgA16SpAIZ8JIkFciAlySpQAa8JEkFMuAlSSqQAS9JUoH8oBtpjcYmZ/yEQ0ltz4CX1mBscobh0amjSwnPzLYYHp0CMOQltRUv0UtrMDI+fTTcF7QOzzEyPt1QR5K0PANeWoMDs6011SWpKQa8tAZbenvWVJekphjw0hoMDfbT0911XK2nu4uhwf6GOpKk5TnJTlqDhYl0zqKX1O4aD/iIOAu4Hngx8DSwJzP/KCJ+H3gvcKja9erMvLWZLqVjdu3caqBLanuNBzxwBPhwZt4VEc8D9kbEbdW2T2fmJxrsTZKkjtR4wGfmQeBgdf+JiNgHeHokSdIpaKtJdhGxHdgJfK8qfSAi7omIayLi+Y01JklSh2mbgI+I5wI3Ax/KzMeBzwEvA3Ywf4b/yRWetzsiJiJi4tChQ8vtIknSptMWAR8R3cyH+w2ZOQqQmQ9l5lxmPg18Hjh3uedm5p7MHMjMgb6+vvqaliSpjTUe8BERwBeAfZn5qUX1Mxbt9jbg3rp7kySpUzU+yQ44H7gMmIqIu6va1cClEbEDSOAB4H3NtCdJUudpPOAz87tALLPJv3mXJOkkNX6JXpIkrT8DXpKkAkVmNt3DuomIQ8BPNuBbnw78bAO+byfyWBzjsTiex+MYj8UxHovjbcTx+JXMfMafkRUV8BslIiYyc6DpPtqBx+IYj8XxPB7HeCyO8Vgcr87j4SV6SZIKZMBLklQgA3519jTdQBvxWBzjsTiex+MYj8UxHovj1XY8fA9ekqQCeQYvSVKBDPhKtSTtwxGx7Gfex7zPRMR91RK2r6m7x7qs4lhcEBGPRcTd1e336u6xLhFxVkR8OyL2RcQPIuKDy+yzKcbGKo/FZhobfycivh8Rf1Udj3+/zD7PiYibqrHxvWpJ7OKs8li8KyIOLRob72mi17pERFdETEbE15fZVsu4aPyjatvItcBngetX2H4hcHZ1+zXml7P9tVo6q9+1nPhYAPzPzHxrPe006gjw4cy8KyKeB+yNiNsy868X7bNZxsZqjgVsnrHxFPCGzPxFtSLmdyPiG5l556J9rgB+npm/GhGXAP8J+I0mmt1gqzkWADdl5gca6K8JHwT2Ab+0zLZaxoVn8JXMvAN49AS7XAxcn/PuBHqXrHhXjFUci00jMw9m5l3V/SeY/we7dclum2JsrPJYbBrV/+9fVA+7q9vSSU0XA9dV978CvLFaQbMoqzwWm0ZEnAn8c+C/rLBLLePCgF+9rcCDix7vZxP/cAP+UXU57hsR8aqmm6lDdRltJ/C9JZs23dg4wbGATTQ2qsuwdwMPA7dl5opjIzOPAI8Bv1xvl/VYxbEA+FfV21hfiYizam6xTn8I/Dvg6RW21zIuDPjVW+63q836G+pdzH804j8E/jMw1nA/Gy4ingvcDHwoMx9funmZpxQ7Np7lWGyqsZGZc5m5AzgTODcizlmyy6YZG6s4Fv8N2J6Zrwb+nGNnsEWJiLcCD2fm3hPttkxt3ceFAb96+4HFv3GeCRxoqJdGZebjC5fjMvNWoDsiTm+4rQ1Tvad4M3BDZo4us8umGRvPdiw229hYkJmzwF8Ab16y6ejYiIjTgL9P4W9/rXQsMvORzHyqevh54LU1t1aX84GLIuIB4E+BN0TEf12yTy3jwoBfvVuAd1Yzps8DHsvMg0031YSIePHC+0URcS7z4+iRZrvaGNV/5xeAfZn5qRV22xRjYzXHYpONjb6I6K3u9wBvAn64ZLdbgMur+28HvpUFfvjIao7FknkpFzE/h6M4mTmcmWdm5nbgEub/n//rJbvVMi6cRV+JiBuBC4DTI2I/8FHmJ4qQmX8C3Aq8BbgPeBJ4dzOdbrxVHIu3A78dEUeAFnBJiT+0KucDlwFT1fuLAFcD22DTjY3VHIvNNDbOAK6LiC7mf5H5cmZ+PSI+Bkxk5i3M/0L0xYi4j/kztEuaa3dDreZY/JuIuIj5v8Z4FHhXY902oIlx4SfZSZJUIC/RS5JUIANekqQCGfCSJBXIgJckqUAGvCRJBTLgJUkqkAEv6ZRExM6ImIuI/9V0L5KOMeAlnar3An8MnBMRr2i6GUnz/KAbSSet+ljSg8DrmV//+ueZ+W+b7UoSeAYv6dS8HfhJZt4DfJH5z+TvbrgnSRjwkk7Ne5gPdoDvMP9Z/Bc1146kBQa8pJMSEb/K/AI0XwKoFpW5gfnQl9QwV5OTdLLeA3QBP61WiAVYWCr2rMx8sKnGJDnJTtJJiIjTgAeBPwK+vmTzF4GvZubHam9M0lEGvKQ1i4iLga8AL87MR5Zs+wjw28BLM/PpJvqT5Hvwkk7OFcC3l4Z75c+AXwHeVG9LkhbzDF6SpAJ5Bi9JUoEMeEmSCmTAS5JUIANekqQCGfCSJBXIgJckqUAGvCRJBTLgJUkqkAEvSVKB/j9OnT5eH2MZ+AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "\n",
    "ax.scatter(df2['A'], df2['Y'])\n",
    "ax.set_xlabel('A', fontsize=14)\n",
    "ax.set_ylabel('Y', fontsize=14);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "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 tr th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe thead tr:last-of-type th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th colspan=\"8\" halign=\"left\">Y</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th>count</th>\n",
       "      <th>mean</th>\n",
       "      <th>std</th>\n",
       "      <th>min</th>\n",
       "      <th>25%</th>\n",
       "      <th>50%</th>\n",
       "      <th>75%</th>\n",
       "      <th>max</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>4.0</td>\n",
       "      <td>70.0</td>\n",
       "      <td>31.622777</td>\n",
       "      <td>40.0</td>\n",
       "      <td>47.5</td>\n",
       "      <td>65.0</td>\n",
       "      <td>87.5</td>\n",
       "      <td>110.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2</td>\n",
       "      <td>4.0</td>\n",
       "      <td>80.0</td>\n",
       "      <td>62.182527</td>\n",
       "      <td>30.0</td>\n",
       "      <td>45.0</td>\n",
       "      <td>60.0</td>\n",
       "      <td>95.0</td>\n",
       "      <td>170.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>3</td>\n",
       "      <td>4.0</td>\n",
       "      <td>117.5</td>\n",
       "      <td>53.774219</td>\n",
       "      <td>50.0</td>\n",
       "      <td>95.0</td>\n",
       "      <td>120.0</td>\n",
       "      <td>142.5</td>\n",
       "      <td>180.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4</td>\n",
       "      <td>4.0</td>\n",
       "      <td>195.0</td>\n",
       "      <td>31.091264</td>\n",
       "      <td>150.0</td>\n",
       "      <td>187.5</td>\n",
       "      <td>205.0</td>\n",
       "      <td>212.5</td>\n",
       "      <td>220.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      Y                                                     \n",
       "  count   mean        std    min    25%    50%    75%    max\n",
       "A                                                           \n",
       "1   4.0   70.0  31.622777   40.0   47.5   65.0   87.5  110.0\n",
       "2   4.0   80.0  62.182527   30.0   45.0   60.0   95.0  170.0\n",
       "3   4.0  117.5  53.774219   50.0   95.0  120.0  142.5  180.0\n",
       "4   4.0  195.0  31.091264  150.0  187.5  205.0  212.5  220.0"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df2.groupby('A').describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### Program 11.2\n",
    "\n",
    "Finally, suppose that treatment $A$ is a variable representing the dose of treatment in mg/day, and that it takes integer values from 0 to 100 mg. This creates a problem: how can we estimate the mean of the outcome $Y$ among individuals with treatment level $A = 90$ in the target population?\n",
    "\n",
    "上面的描述表明，我们不能总是让数据“为自己说话”以获得有意义的估计。相反，我们经常需要用模型来补充数据，如我们在下一节中所述。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "A, Y = zip(*(\n",
    "    (3, 21),\n",
    "    (11, 54),\n",
    "    (17, 33),\n",
    "    (23, 101),\n",
    "    (29, 85),\n",
    "    (37, 65),\n",
    "    (41, 157),\n",
    "    (53, 120),\n",
    "    (67, 111),\n",
    "    (79, 200),\n",
    "    (83, 140),\n",
    "    (97, 220),\n",
    "    (60, 230),\n",
    "    (71, 217),\n",
    "    (15, 11),\n",
    "    (45, 190),\n",
    "))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "df3 = pd.DataFrame({'A': A, 'Y': Y, 'constant': np.ones(16)})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfoAAAF3CAYAAABNO4lPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAWdElEQVR4nO3df4xlZ3kf8O/T9aadkFaDa0O9YyeG1ppC44ZFK+qWKHJD1ME0irdWUInaxEJQVxVRSZtM8eYf+lOm3TRpUAKSGygmIlBKlsVKo2yRoSWtBGGdiVgTZ4RFwPasa29CBmiZJsvy9o+5Y8/as7uenb33zLz385FG9573nrnz6OhovnPe8859qrUWAKBPf2roAgCA8RH0ANAxQQ8AHRP0ANAxQQ8AHRP0ANCxq4Yu4Eq65ppr2o033jh0GQAwMQ8++OAftNauvdDrXQX9jTfemJMnTw5dBgBMTFV9+WKvm7oHgI4JegDomKAHgI4JegDomKAHgI4JegDomKAHgI4JegDomKAHgI519cl4ABdzfGklR08s5/TqWg7MzmRxYT6HD84NXRaMlaAHpsLxpZUcOXYqa2fPJUlWVtdy5NipJBH2dM3UPTAVjp5YfjrkN6ydPZejJ5YHqggmQ9ADU+H06tq2xqEXgh6YCgdmZ7Y1Dr0Q9MBUWFyYz8z+feeNzezfl8WF+YEqgsmwGA+YChsL7qy6Z9oIemBqHD44J9iZOqbuAaBjgh4AOiboAaBjgh4AOiboAaBjgh4AOubf6wBgzIbsnCjoAWCMhu6caOoeAMZo6M6Jgh4Axmjozomm7gH2gCHv8bIzB2ZnsrJFqE+qc6IreoBdbuMe78rqWlqeucd7fGll6NJ4HobunCjoAXa5oe/xsjOHD87lnjtuztzsTCrJ3OxM7rnjZqvuAVg39D1edm7Izomu6AF2uQvdy53UPV72NkEPsMsNfY+Xvc3UPcAutzHla9U9l0PQA+wBQ97jZW8zdQ8AHRP0ANAxQQ8AHRP0ANAxQQ8AHRP0ANAxQQ8AHRP0ANAxQQ8AHRP0ANAxQQ8AHRP0ANAxQQ8AHRP0ANAxQQ8AHZtY0FfVDVX1yap6uKo+X1VvHY1fXVUfr6ovjB5fOBqvqnpnVT1SVZ+rqldOqlYALuz40kpe/Y5P5CV3/9e8+h2fyPGllaFL4iImeUX/zSQ/2Vp7WZJbkrylql6e5O4kD7TWbkrywGg7SW5LctPo664k755grQBs4fjSSo4cO5WV1bW0JCurazly7JSw38UmFvSttSdaa789ev71JA8nmUtye5L7Rrvdl+Tw6PntSd7f1n06yWxVXTepegF4rqMnlrN29tx5Y2tnz+XoieWBKuJSBrlHX1U3JjmY5DNJXtxaeyJZ/2MgyYtGu80leWzTtz0+Gnv2e91VVSer6uSZM2fGWTbA1Du9uratcYY38aCvqu9I8qtJfqK19rWL7brFWHvOQGv3ttYOtdYOXXvttVeqTAC2cGB2ZlvjDG+iQV9V+7Me8h9orR0bDT+5MSU/enxqNP54khs2ffv1SU5PqlYAnmtxYT4z+/edNzazf18WF+YHqohLmeSq+0ryniQPt9Z+dtNL9ye5c/T8ziQf2zT+Y6PV97ck+erGFD8Awzh8cC733HFz5mZnUknmZmdyzx035/DB59xZZZeo1p4zGz6eH1T1vUl+M8mpJN8aDf901u/TfzjJdyZ5NMnrW2tfGf1h8AtJXpvkG0ne2Fo7ebGfcejQoXby5EV3AYCuVNWDrbVDF3r9qkkV0lr7n9n6vnuSvGaL/VuSt4y1KADonE/GA4COCXoA6JigB4COCXoA6JigB4COTWzVPTCs40srOXpiOadX13JgdiaLC/P+9xmmgKCHKbDRcWyjGclGx7Ekwh46Z+oepoCOYzC9BD1MAR3HYHoJepgCOo7B9BL0MAV0HIPpZTEeTIGNBXdW3cP0EfQwJQ4fnBPsMIVM3QNAxwQ9AHRM0ANAxwQ9AHRM0ANAxwQ9AHRM0ANAxwQ9AHRM0ANAxwQ9AHRM0ANAxwQ9AHRM0ANAxwQ9AHRM0ANAxwQ9AHRM0ANAxwQ9AHRM0ANAxwQ9AHRM0ANAxwQ9AHRM0ANAxwQ9AHRM0ANAxwQ9AHRM0ANAxwQ9AHRM0ANAxwQ9AHRM0ANAxwQ9AHRM0ANAxwQ9AHRM0ANAxwQ9AHTsqqELAC7t+NJKjp5YzunVtRyYncniwnwOH5wbuixgDxD0sMsdX1rJkWOnsnb2XJJkZXUtR46dShJhD1ySqXvY5Y6eWH465DesnT2XoyeWB6oI2EsEPexyp1fXtjUOsJmgh13uwOzMtsYBNhP0sMstLsxnZv++88Zm9u/L4sL8QBUBe4nFeLDLbSy4s+oeuByCHvaAwwfnBDtwWUzdA0DHJhb0VfXeqnqqqh7aNPbPq2qlqn5n9PW6Ta8dqapHqmq5qhYmVScA9GSSV/TvS/LaLcZ/rrX2itHXrydJVb08yRuS/JXR97yrqvZt8b0AwEVMLOhba59K8pXnufvtST7UWvvj1trvJ3kkyavGVhwAdGo33KP/8ar63Ghq/4Wjsbkkj23a5/HRGACwDUMH/buT/MUkr0jyRJJ/PxqvLfZtW71BVd1VVSer6uSZM2fGUyUA7FGDBn1r7cnW2rnW2reS/Mc8Mz3/eJIbNu16fZLTF3iPe1trh1prh6699trxFgwAe8ygQV9V123a/DtJNlbk35/kDVX1p6vqJUluSvJbk64PAPa6iX1gTlV9MMmtSa6pqseTvD3JrVX1iqxPy38pyT9Mktba56vqw0l+N8k3k7yltXZuq/cFAC6sWtvy1veedOjQoXby5MmhywCAiamqB1trhy70+tCL8QCAMRL0ANAxQQ8AHRP0ANAxQQ8AHdOPHoCJO760kqMnlnN6dS0HZmeyuDCfwwd90vk4CHoAJur40kqOHDuVtbPrH4+ysrqWI8dOJYmwHwNT9wBM1NETy0+H/Ia1s+dy9MTyQBX1TdADMFGnV9e2Nc7OCHoAJurA7My2xtkZQQ/ARC0uzGdm/77zxmb278viwvxAFfXNYjwAJmpjwZ1V95Mh6AGYuMMH5wT7hJi6B4COCXoA6JigB4COCXoA6JigB4COCXoA6JigB4COCXoA6JigB4COCXoA6JigB4COCXoA6JigB4COCXoA6JigB4COCXoA6JigB4COCXoA6JigB4COXXWxF6vqltbapydVDNC/40srOXpiOadX13JgdiaLC/M5fHBu6LKgW5e6ov9UVf2rqrroHwQAz8fxpZUcOXYqK6traUlWVtdy5NipHF9aGbo06Nalgv62JD+a5Leq6uUTqAfo2NETy1k7e+68sbWz53L0xPJAFUH/Lhr0rbUHktycZCnJyar6pxOpCujS6dW1bY0DO3fJxXitta+31t6U9Sv7f1dV/6eqvrb5a/xlAj04MDuzrXFg557XvfeqOpTkXyf5QpKfSfLNcRYF9GlxYT5Hjp06b/p+Zv++LC7MD1gV9O1Sq+6vSvL2JG9L8q4kd7fW/t8kCgP6s7G63qp7mJxLXdF/NsnVSW4b3a8H2JHDB+cEO0zQpe7RP5Tkrwp5ANibLnpF31r70UkVAgBceT4CFwA6JugBoGOCHgA65jPsAcZMIx+GJOgBxmijkc/GhwRtNPJJIuyZCFP3AGOkkQ9DE/QAY6SRD0MT9ABjpJEPQxP0AGO0uDCfmf37zhvTyIdJshgPYIw08mFogh5gzDTyYUim7gGgY4IeADom6AGgY4IeADom6AGgY1bdM1GaewBM1sSu6KvqvVX1VFU9tGns6qr6eFV9YfT4wtF4VdU7q+qRqvpcVb1yUnUyPhvNPVZW19LyTHOP40srQ5cG0K1JTt2/L8lrnzV2d5IHWms3JXlgtJ0ktyW5afR1V5J3T6hGxkhzD4DJm1jQt9Y+leQrzxq+Pcl9o+f3JTm8afz9bd2nk8xW1XWTqZRx0dwDYPKGXoz34tbaE0kyenzRaHwuyWOb9nt8NPYcVXVXVZ2sqpNnzpwZa7HsjOYeAJM3dNBfSG0x1rbasbV2b2vtUGvt0LXXXjvmstgJzT0AJm/oVfdPVtV1rbUnRlPzT43GH09yw6b9rk9yeuLVcUVp7gEweUMH/f1J7kzyjtHjxzaN/3hVfSjJX0vy1Y0pfvY2zT0AJmtiQV9VH0xya5JrqurxJG/PesB/uKrelOTRJK8f7f7rSV6X5JEk30jyxknVCQA9mVjQt9Z+5AIvvWaLfVuSt4y3IgDo325djAcAXAGCHgA6JugBoGOCHgA6JugBoGOCHgA6JugBoGOCHgA6JugBoGOCHgA6JugBoGOCHgA6NnSbWhiL40sr+t4DRNDToeNLKzly7FTWzp5LkqysruXIsVNJIuyBqWPqnu4cPbH8dMhvWDt7LkdPLA9UEcBwBD3dOb26tq1xgJ4JerpzYHZmW+MAPRP0dGdxYT4z+/edNzazf18WF+YHqghgOBbj0Z2NBXdW3QMIejp1+OCcYAeIqXsA6JqgB4COCXoA6JigB4COCXoA6JigB4COCXoA6JigB4COCXoA6JigB4COCXoA6JigB4COCXoA6JigB4COCXoA6JigB4COCXoA6JigB4COCXoA6NhVQxcAPTu+tJKjJ5ZzenUtB2Znsrgwn8MH54YuC5gigh7G5PjSSo4cO5W1s+eSJCurazly7FSSCHtgYkzdw5gcPbH8dMhvWDt7LkdPLA9UETCNBD2MyenVtW2NA4yDoIcxOTA7s61xgHEQ9DAmiwvzmdm/77yxmf37srgwP1BFwDSyGA/GZGPBnVX3wJAEPYzR4YNzgh0YlKl7AOiYoAeAjgl6AOiYoAeAjgl6AOiYoAeAjgl6AOiY/6Pf5bQ5BWAnBP0ups0pADtl6n4X0+YUgJ0S9LuYNqcA7JSg38W0OQVgp3ZF0FfVl6rqVFX9TlWdHI1dXVUfr6ovjB5fOHSdk6bNKQA7tSuCfuRvttZe0Vo7NNq+O8kDrbWbkjww2p4qhw/O5Z47bs7c7EwqydzsTO6542YL8QB43nbzqvvbk9w6en5fkv+e5G1DFTMUbU4B2IndckXfkvy3qnqwqu4ajb24tfZEkoweX7TVN1bVXVV1sqpOnjlzZkLlAsDesFuu6F/dWjtdVS9K8vGq+r3n+42ttXuT3Jskhw4dauMqEAD2ol1xRd9aOz16fCrJR5O8KsmTVXVdkowenxquQgDYmwYP+qp6QVX92Y3nSf5WkoeS3J/kztFudyb52DAVAsDetRum7l+c5KNVlazX8yuttd+oqs8m+XBVvSnJo0leP2CNALAnDR70rbUvJvmeLcb/MMlrJl8RAPRj8Kl7AGB8BD0AdEzQA0DHBD0AdEzQA0DHBD0AdEzQA0DHBD0AdEzQA0DHBD0AdEzQA0DHBD0AdEzQA0DHBD0AdGzwNrXsDseXVnL0xHJOr67lwOxMFhfmc/jg3NBlAbBDgp4cX1rJkWOnsnb2XJJkZXUtR46dShJhD7DHmbonR08sPx3yG9bOnsvRE8sDVQTAlSLoyenVtW2NA7B3CHpyYHZmW+MA7B2CniwuzGdm/77zxmb278viwvxAFQFwpViMx9ML7qy6B+iPoCfJetgLdoD+mLoHgI4JegDomKAHgI4JegDomKAHgI4JegDomKAHgI4JegDomA/MuQD92QHogaDfgv7sAPTC1P0W9GcHoBeCfgv6swPQC0G/Bf3ZAeiFoN+C/uwA9MJivC3ozw5ALwT9BejPDkAPTN0DQMcEPQB0TNADQMcEPQB0TNADQMcEPQB0zL/XTQGd+ACml6DvnE58ANPN1H3ndOIDmG6CvnM68QFMN0HfOZ34AKaboO+cTnwA081ivM7pxAcw3QT9FNCJD2B6mboHgI4JegDomKAHgI4JegDomKAHgI4JegDomKAHgI4JegDomKAHgI4JegDoWLXWhq7hiqmqM0m+fJFdrknyBxMqZ5o4ruPhuI6PYzsejut4XOq4fldr7doLvdhV0F9KVZ1srR0auo7eOK7j4biOj2M7Ho7reOz0uJq6B4COCXoA6Ni0Bf29QxfQKcd1PBzX8XFsx8NxHY8dHdepukcPANNm2q7oAWCqTEXQV9Vrq2q5qh6pqruHrmevqqobquqTVfVwVX2+qt46Gr+6qj5eVV8YPb5w6Fr3oqraV1VLVfVro+2XVNVnRsf1P1fVtw1d415UVbNV9ZGq+r3RufvXnbM7V1X/ZPR74KGq+mBV/Rnn7OWpqvdW1VNV9dCmsS3P0Vr3zlGefa6qXnmp9+8+6KtqX5JfTHJbkpcn+ZGqevmwVe1Z30zyk621lyW5JclbRsfy7iQPtNZuSvLAaJvte2uShzdt/9skPzc6rn+U5E2DVLX3/XyS32it/eUk35P1Y+yc3YGqmkvyj5Mcaq19d5J9Sd4Q5+zlel+S1z5r7ELn6G1Jbhp93ZXk3Zd68+6DPsmrkjzSWvtia+1Pknwoye0D17QntdaeaK399uj517P+C3Mu68fzvtFu9yU5PEyFe1dVXZ/kbyf5pdF2Jfn+JB8Z7eK4Xoaq+nNJvi/Je5KktfYnrbXVOGevhKuSzFTVVUm+PckTcc5eltbap5J85VnDFzpHb0/y/rbu00lmq+q6i73/NAT9XJLHNm0/PhpjB6rqxiQHk3wmyYtba08k638MJHnRcJXtWf8hyT9L8q3R9p9Pstpa++Zo23l7eV6a5EyS/zS6LfJLVfWCOGd3pLW2kuRnkjya9YD/apIH45y9ki50jm4706Yh6GuLMf9qsANV9R1JfjXJT7TWvjZ0PXtdVf1gkqdaaw9uHt5iV+ft9l2V5JVJ3t1aO5jk/8Y0/Y6N7hffnuQlSQ4keUHWp5SfzTl75W37d8M0BP3jSW7YtH19ktMD1bLnVdX+rIf8B1prx0bDT25MHY0enxqqvj3q1Ul+qKq+lPVbS9+f9Sv82dG0aOK8vVyPJ3m8tfaZ0fZHsh78ztmd+YEkv99aO9NaO5vkWJK/EefslXShc3TbmTYNQf/ZJDeNVoN+W9YXjNw/cE170ui+8XuSPNxa+9lNL92f5M7R8zuTfGzSte1lrbUjrbXrW2s3Zv38/ERr7e8l+WSSHx7t5rhehtba/07yWFXNj4Zek+R345zdqUeT3FJV3z76vbBxXJ2zV86FztH7k/zYaPX9LUm+ujHFfyFT8YE5VfW6rF8h7Uvy3tbavxm4pD2pqr43yW8mOZVn7iX/dNbv0384yXdm/RfA61trz15YwvNQVbcm+anW2g9W1UuzfoV/dZKlJH+/tfbHQ9a3F1XVK7K+yPHbknwxyRuzfpHjnN2BqvoXSf5u1v8bZynJm7N+r9g5u01V9cEkt2a9S92TSd6e5Hi2OEdHf1j9QtZX6X8jyRtbaycv+v7TEPQAMK2mYeoeAKaWoAeAjgl6AOiYoAeAjgl6AOiYoAeAjgl6YEeq6mBVnauq/zV0LcBzCXpgp/5Bkncl+e6qetnQxQDn84E5wGWrqpmsdy/7viRvTfJHrbWfGrYqYDNX9MBO/HCSL7fWPpfkl7P+Gdz7B64J2ETQAzvx5qwHfJL8j6x/9vYPDVcO8GyCHrgsVfWXst5i91eSpK3fB/xA1sMf2CWuuvQuAFt6c9Y7Qj663lArSVJJUlU3tNYeG6ow4BkW4wHbVlVXJXksyc8n+bVnvfzLST7aWvuXEy8MeA5BD2xbVd2e5CNJ/kJr7Q+f9drbkvyjJC9trX1riPqAZ7hHD1yONyX55LNDfuS/JPmuJD8w2ZKArbiiB4COuaIHgI4JegDomKAHgI4JegDomKAHgI4JegDomKAHgI4JegDomKAHgI79f91ts94yDWGCAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "\n",
    "ax.scatter(df3.A, df3.Y)\n",
    "ax.set_xlabel('A', fontsize=14)\n",
    "ax.set_ylabel('Y', fontsize=14);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "\"The values $\\hat{\\theta}_0$ and $\\hat{\\theta}_1$ can be easily computed using linear algebra, as described in any statistics textbook.\"\n",
    "\n",
    "Briefly, make $X$ be the matrix with ones in the first column and $A$ in the second column, and make $Y$ an $n \\times 1$ matrix. Then $\\hat{\\theta}_0$ and $\\hat{\\theta}_1$ are calculated as $\\hat\\theta = (X^TX)^{-1}X^TY$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[24.54636872],\n",
       "       [ 2.13715198]])"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = np.ones((df3.shape[0], 2))        # make X a matrix of ones\n",
    "X[:, 1] = df3.A                       # replace second column with A\n",
    "Y = np.array(df3.Y).reshape((-1, 1))  # make Y an n x 1 matrix\n",
    "\n",
    "theta = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(Y)\n",
    "theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "theta0 est.: 24.55\n",
      "theta1 est.:  2.14\n"
     ]
    }
   ],
   "source": [
    "print(\"theta0 est.: {:>5.2f}\".format(theta[0][0]))\n",
    "print(\"theta1 est.: {:>5.2f}\".format(theta[1][0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E[Y|A=90] est.: 216.9\n"
     ]
    }
   ],
   "source": [
    "expectation_at_90 = theta[0][0] + 90 * theta[1][0]\n",
    "print(\"E[Y|A=90] est.: {:>5.1f}\".format(expectation_at_90))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "We can also compute this with a statistics package, like Statsmodels, which will also give us confidence intervals and other values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "ols = sm.OLS(Y, df3[['constant', 'A']])\n",
    "res = ols.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>        <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>constant</th> <td>   24.5464</td> <td>   21.330</td> <td>    1.151</td> <td> 0.269</td> <td>  -21.202</td> <td>   70.295</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>A</th>        <td>    2.1372</td> <td>    0.400</td> <td>    5.347</td> <td> 0.000</td> <td>    1.280</td> <td>    2.994</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "summary = res.summary()\n",
    "summary.tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "Initially, it wasn't clear how to get confidence intervals for expected values from Statsmodels, so the next cells calculate them from scratch."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "n = df3.shape[0]\n",
    "yvar = (res.resid * res.resid).sum() / (n - 2)  # = res.mse_resid\n",
    "xval = np.array([[1, 90]])\n",
    "X = df3[['constant', 'A']]\n",
    "XpXinv = np.linalg.inv(np.dot(X.T, X))\n",
    "se_mean = np.sqrt(yvar * np.dot(xval, np.dot(XpXinv, xval.T)))[0, 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "           estimate      95% C.I.\n",
      "E[Y|A=90]   216.89   (172.15, 261.63)\n"
     ]
    }
   ],
   "source": [
    "t = scipy.stats.t.ppf(0.975, n - 2)\n",
    "ypred = res.predict([[1, 90]])[0]\n",
    "print('           estimate      95% C.I.')\n",
    "print(\n",
    "    'E[Y|A=90]   {:>6.2f}   ({:>6.2f}, {:>6.2f})'.format(\n",
    "        ypred, ypred - t * se_mean, ypred + t * se_mean\n",
    "))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "Here's an easier way to calculate the confidence intervals, contributed by GitHub user @pettypychen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "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>mean</th>\n",
       "      <th>mean_se</th>\n",
       "      <th>mean_ci_lower</th>\n",
       "      <th>mean_ci_upper</th>\n",
       "      <th>obs_ci_lower</th>\n",
       "      <th>obs_ci_upper</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>0</td>\n",
       "      <td>216.89</td>\n",
       "      <td>20.86</td>\n",
       "      <td>172.15</td>\n",
       "      <td>261.63</td>\n",
       "      <td>112.27</td>\n",
       "      <td>321.51</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     mean  mean_se  mean_ci_lower  mean_ci_upper  obs_ci_lower  obs_ci_upper\n",
       "0  216.89    20.86         172.15         261.63        112.27        321.51"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pred = res.get_prediction(exog=[1, 90])\n",
    "pred.summary_frame().round(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "Collapsed": "false",
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfoAAAF3CAYAAABNO4lPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXhU5fnG8fslBI0IRASBBDAgGARkM2wiiICERQv+qtZWxZ1apdbWphJFQdwocbetldYNd6sxUqRERAF3DUQIWxAFgSTshjVASN7fH5mcJBgg28yZOfP9XBdX8rwzCY/jITfPnMN5jbVWAADAm+q53QAAAPAfgh4AAA8j6AEA8DCCHgAADyPoAQDwMIIeAAAPq+92A3WpWbNmNi4uzu02AAAImMWLF2+31jY/2uOeCvq4uDhlZGS43QYAAAFjjPnxWI/z1j0AAB5G0AMA4GEEPQAAHkbQAwDgYQQ9AAAeRtADAOBhBD0AAB5G0AMA4GEEPQAAHkbQAwDgYQQ9AAAeRtADABAA+fsPacOO/QH/fT21qQ0AHEtaZo5S0rOVm1+gmOgoJSXGa2zPWLfbQhiYPne1/rHge0nSuodHyRgTsN+boAcQFtIyc5ScmqWCwiJJUk5+gZJTsySJsIffbNl9QH0fmu/Uk0afFdCQlwh6AGEiJT3bCflSBYVFSknPJujhF9P+t1r/XPi9U39774WKPqlBwPsg6AGEhdz8gmqtAzWVt6tA/R/+yKnvHnWWbhrU3rV+CHoAYSEmOko5lYR6THSUC93Aqx6as0ozFv3g1EvvHa4mJ0W62BFX3QMIE0mJ8YqKjKiwFhUZoaTEeJc6gpfk5hcobuL7Tsjfc1FnrZ822vWQl5joAYSJ0vPwXHWPunb/7JV67tN1Tr108nA1iXI/4EsR9ADCxtiesQQ76kxOfoEGTCs7Fz/l4s66dkA7FzuqHEEPAEA1TZm1Qi9+vt6pl00ZrsYnBs8UXx5BDwBAFW3cuV8Dp3/s1FPHdNG4/nHuNVQFBD0AAFVw73vLNfOLH506mKf48gh6AACO4cgp/v6xXXV1v9Nd7Kh6CHoAAI5iUlqWXvlyg1NnTRmuRiEwxZdH0AMAcIQNO/ZrUErZFP/QJWfrN33b1vj7ubmhEkEPAEA5yalZev3rsil++X2JOvmEmsel2xsqEfQAAEj6ccc+nZ+ywKmn/d/ZuqJPzaf4Um5vqETQAwDC3p1vL9ObGRudesV9iWpYiym+PLc3VCLoASAEuHmO18vWbd+nCx5Z4NTTL+2myxPa1Onv4faGSgQ9AAQ5t8/xetUdby3VO0s2OfXKqYk6qUHdx2JSYnyF/39SYDdUIugBIMi5fY7Xa37YtldDHl3o1CmXdtNldTzFl+f2hkoEPQAEObfP8XrJn978VqmZOZKkiHpGy6ckKqpBxHG+qvbc3FCJoAeAIOf2OV4vWLt1r4Y9VjbFP3Z5d/1fr9YudhQ49dxuAABwbEmJ8YqKrDh1BvIcb6i7/Y1MJ+QbRNTTqqkjwibkJSZ6AAh6bp/jDVVrt+7RsMcWOfUTv+oRlq8ZQQ8AIcDNc7yhaMJrSzR7WZ4k6cTIevr23uE6MdL/5+KDEUEPAPCMNVv2aPjjZVP8U7/uqV90j3GxI/cR9AAAT7jl1cWak7VZknTyCfWVMWlY2E7x5RH0AICQlr15jxKfKJvin/51T10c5lN8eQQ9ACAkWWt18yuLlb5iiySp0Yn19c3dTPFHIugBACFnVd5ujXzyE6f+x5W9NOrsVi52FLwIegBAyLDWavzLizVvZckUH31SpL66a6hOqM8UfzQEPQAgJKzM3a1RT5VN8f+8qpdGdGWKPx6CHgAQ1Ky1uuGlDH20eqsk6dSGDfRF8lA1qM/NXauCoAcABK3lObt00dOfOvWzV5+jxC4tXewo9BD0AICgY63VtS98o4VrtkmSTmt0gj69cwhTfA0Q9ACAoHLkFP+vcQm6sHMLFzsKbQQ9ACAoWGs17vmv9cl32yVJLRufqE/uvECREUzxtUHQAwBct2xTvn7xt8+c+rlrEjT0LKb4ukDQAwBcY63VVc99pc/W7pAkxUZHaUHSYKb4OkTQAwBcsXRjvsb8vWyKf+Ha3rqg02kuduRNAQt6Y0wbSTMltZRULGmGtfZJY0xTSW9KipO0XtLl1tqfjDFG0pOSRknaL+laa+2SQPULAKhcWmaOUtKzlZtfoJjoKCUlxmtsz9gqf721VlfM+FJfrdspSWrb9CR9dMf5qs8U7xeBnOgPS7rDWrvEGNNI0mJjzDxJ10qab62dZoyZKGmipDsljZTU0ferr6RnfB8BAC5Jy8xRcmqWCgqLJEk5+QVKTs2SpCqFfeaGn3TJPz536hev663B8Uzx/hSwoLfW5knK832+xxizSlKspDGSBvue9pKkBSoJ+jGSZlprraQvjTHRxphWvu8DAHBBSnq2E/KlCgqLlJKefcygt9bq8me/0Dfrf5IktWvWUPP+OIgpPgBcOUdvjImT1FPSV5JalIa3tTbPGFP6V7tYSRvLfdkm31qFoDfGjJc0XpLatm3r174BINzl5hdUa12SFv/4k375TNkUP/P6Php0ZvM67w2VC3jQG2NOlvSOpNuttbtLTsVX/tRK1uzPFqydIWmGJCUkJPzscQBA3YmJjlJOJaEeEx31s7XiYqtL//m5lmzIlySd0byh0m9nig+0gL7axphIlYT8q9baVN/yFmNMK9/jrSRt9a1vktSm3Je3lpQbqF4BAD+XlBivqMiKW8JGRUYoKTG+wlrG+p1qf9ccJ+RfuaGv5t8xmJB3QSCvujeSnpO0ylr7WLmHZkm6RtI038f3yq1PMMa8oZKL8HZxfh4A3FV6Hv5oV90XF1td8sznWrqxJODPbHGy/veHQYqod9R3b+FnpuRatwD8RsacJ+kTSVkq+ed1knSXSs7TvyWpraQNki6z1u70/cXgb5JGqOSf111nrc041u+RkJBgMzKO+RQAgJ98vW6nLn/2C6d+7ca+OrdDMxc7Cg/GmMXW2oSjPR7Iq+4/VeXn3SVpaCXPt5Ju9WtTAIBaKy62GvP3z5SVs0uS1KllI71/20Cm+CDBnfEAADX21Q879KsZXzr16zf1U/8zTnWxIxyJoAcAVFtxsdVFT3+qlXm7JUldYhrrvxPOUz2m+KBD0AMAquWL73fo1/8qm+LfGN9P/dozxQcrgh4AUCVFxVajn/pEqzfvkSR1a91EabcMYIoPcgQ9AOC4Pl+7Xb/591dO/dZv+6tPu6YudoSqIuiBMFHbHccQnoqKrUY8sUjfbd0rSerRJlqpvzuXKT6EEPRAGKjtjmMIT59+t11XPVc2xb99c38lxDHFhxqCHggDNd1xDOGpqNjqwscX6odt+yRJ55x+iv7z2/5M8SGKoAfCQE12HEN4WrRmm8Y9/7VTv/O7c3XO6ae42BFqi6AHwkB1dhxDeDpcVKxhjy3U+h37JUl94prqzd/20zF2GEWIYBshIAxUdccxhKeFa7apw93/c0I+9ZZz9dbN/Ql5j2CiB8LA8XYcQ3g6XFSsCx5doI07S97t6duuqd4YzxTvNQQ9ECbG9owl2OH4ePVWXffiN06ddusA9WgT7WJH8BeCHgDCSGFRsQanLHCu2RjQ4VS9ckNfpngPI+gBIEzMX7VFN7yU4dSzJgxQt9ZM8V5H0AOAxxUWFWvgXz/W5t0HJEkDOzbTzOv7MMWHCYIeADxs3sotumlm2RQ/+/fnqWtsExc7QqAR9ADgQYcOF2vAXz/Stj0HJUmD45vrhWt7M8WHIYIeADzmgxWbNf7lxU79/m3nqUsMU3y4IugBwCMOHS5W/4fna8e+Q5KkYWedpn+NS2CKD3MEPQB4wNzlebr5lSVOPee2geoc09jFjhAsCHoACGEHDxep70Pzlb+/UJJ0YecWmnH1OUzxcBD0ABCi5mTl6ZZXy6b4ubcPVKeWTPGoiKAHgBBz8HCREh74UHsOHJYkjejSUs9c1YspHpUi6AEghMxelqsJr2U6dfrtgxTfspGLHSHYEfQAEAJ2FRSq+30fOPXos1vp71f2crEjhAqCHgCC3PiZGfpg5RannvfHQerYgikeVUPQA0CQ2rW/UN2nflBhbf200S51g1BF0ANAELr+xW/00eqtTv3Pq3ppRNdWLnaEUEXQA0AQyd9/SD2mzquwxhSP2iDoASBIjHv+ay1as82p/zUuQRd2buFiR/ACgh4AXLZz3yH1up8pHv5B0AOAi37zry/1+fc7nPr5axM0pBNTPOoOQQ8ALtix96DOeeDDCmtM8fAHgh4AAuzyZ7/Q1+t2OvWL1/XW4PjTXOwIXkbQA0CAbNtzUL0fZIpHYBH0ABAAv3zmcy3+8Sennnl9Hw06s7mLHSFcEPQA4Edb9xxQnwfnV1hjikcgEfQA4Cdj/vaplm7a5dSv3thXAzo0c7EjhCOCHgDq2JbdB9T3IaZ4BAeCHgDq0KgnP9HKvN1O/fpN/dT/jFNd7AjhjqAHgDqwedcB9XuYKR7Bh6AHgFpKfHyRsrfsceo3x/dT3/ZM8QgOBD0A1FBOfoEGTPuowhpTPIINQQ8ANTDkkQX6Yfs+p/7Pzf3VO66pix0BlSPoAaAaNu7cr4HTP66wxhSPYEbQA0AVDZr+sTbs3O/U7/zuXJ1z+ikudgQcH0EPAMfBFI9QRtADwDGc+/B85e464NRptw5QjzbRLnYEVA9BD4SAtMwcpaRnKze/QDHRUUpKjNfYnrFut+VpP+7Yp/NTFlRYY4pHKCLogSCXlpmj5NQsFRQWSSr5J13JqVmSRNj7Se8HP9S2PQedetaEAerWmikeoame2w0AOLaU9Gwn5EsVFBYpJT3bpY68a932fYqb+H6FkF8/bTQhj5DGRA8Eudz8gmqto2Z6TP1A+fsLnXr2789T19gmLnYE1A2CHghyMdFRyqkk1GOio1zoxnu+37ZXQx9dWGGNc/HwEoIeCHJJifEVztFLUlRkhJIS413syhvOnpyuPQcPO/Wc2waqc0xjFzsC6h5BDwS50gvuuOq+7qzdukfDHltUYY0pHl5F0AMhYGzPWIK9jnS65386UFjs1HNvH6hOLZni4V0EPYCwsH3vQSU88KFT169ntPahUS52BARGwP55nTHmeWPMVmPM8nJrU4wxOcaYb32/RpV7LNkYs9YYk22MSQxUnwC858kPv6sQ8h/8cRAhj7ARyIn+RUl/kzTziPXHrbWPlF8wxnSWdIWkLpJiJH1ojDnTWlskAKiibXsOqveDZQH/5+FnasKQji52BARewILeWrvIGBNXxaePkfSGtfagpHXGmLWS+kj6wk/tAfCYx+at0VPzv3PqJfdcqKYNG7jYEeCOYDhHP8EYM05ShqQ7rLU/SYqV9GW552zyrf2MMWa8pPGS1LZtWz+3CiDYbd1zQH0enO/UfxkRr1sGd3CxI8Bdbt8C9xlJZ0jqISlP0qO+dVPJc21l38BaO8Nam2CtTWjevLl/ugQQEh79ILtCyC+550JCHmHP1YneWrul9HNjzL8kzfaVmyS1KffU1pJyA9gagBCydfcB9XmoLOAnjuykm88/w8WOgODhatAbY1pZa/N85SWSSq/InyXpNWPMYyq5GK+jpK9daBFAkJs+d7X+seB7p86850Kdwrl4wBGwoDfGvC5psKRmxphNkiZLGmyM6aGSt+XXS/qtJFlrVxhj3pK0UtJhSbdyxT2A8jbvOqB+D5dN8XePOks3DWrvYkdAcDLWVnrqOyQlJCTYjIwMt9sA4GcP/2+Vnl34g1N/e++Fij6JKR7hyRiz2FqbcLTHg+GqewCokrxdBer/8EdOPWn0WbpxIFM8cCwEPYCQ8NCcVZqxqGyKX3rvcDU5KdLFjoDQQNADCGq5+QU6d1rZFD/54s66bkA7FzsCQgtBDyBo3T97pZ77dJ1TL508XE2imOKB6iDoAQSdnPwCDSg3xU8d00Xj+se51xAQwgh6AEFlyqwVevHz9U69bMpwNT6RKd5r0jJzlJKerdz8AsVERykpMV5je1Z6p3PUEkEPIChs3LlfA6d/7NT3j+2qq/ud7mJH8Je0zBwlp2apoLDk9ig5+QVKTs2SJMLeDwh6AK67973lmvnFj06dNWW4GjHFe1ZKerYT8qUKCouUkp5N0PsBQQ/ANUdO8Q9e0lVX9mWK97rc/IJqraN2CHoArpiUlqVXvtzg1MvvS9TJJ/AjKRzEREcpp5JQj4mOcqEb73N7m1oAYWbDjv2Km/i+E/IP/9/ZWj9tNCEfRpIS4xUVGVFhLSoyQkmJ8S515G38yQIQMMmpy/T61xudesV9iWpIwIed0vPwXHUfGPwJA+B367fv0+BHFjj19F920+W927jXEFw3tmcswR4gBD0Av7rz7WV6M4MpHnALf9oA+MW67ft0QbkpPuXSbrosgSkeCDSCHkCdu+OtpXpnySZJUj1TckX9SQ34cQO4gT95AOrMD9v2asijC5360cu665fntHaxIwAEPYA68ac3v1VqZo4kKTLCaNnkREU1iDjOVwHwN4IeQK2s3bpXwx4rm+If/1V3XdKTKR4IFgQ9gBq77fVMzVqaK0k6oX49LZ08XCdGMsUDwYSgB1Bt323ZowsfX+TUT17RQ2N68G+igWBE0AOolgmvLdHsZXmSpJMaRGjJPRcyxQNBjKAHUCVrtuzR8HJT/NO/7qmLu8e42BGAqiDoARzXLa8u1pyszZKkRifU1zeThjHFAyGCoAdwVNmb9yjxibIp/u+/6aXR3Vq52BGA6iLoAfyMtVa/fXmxPli5RZLUJCpSX989VCfUZ4oHQg1BD6CCVXm7NfLJT5z6mSt7aeTZTPFAqCLoAUgqmeJvmpmhD1dtlSQ1bdhAXyQPYYoHQhxBD0Arcndp9FOfOvU/rzpHI7q2dLEjAHWFoAfCmLVWN7yUoY9Wl0zxzU4+QZ9PHKIG9eu53BmAunLMoDfG9LPWfhmoZgAEzvKcXbro6bIpfsbV52h4F6Z4wGuON9EvMsb8VdJ91trDgWgIgH9Za3XNC99o0ZptkqQWjU/QJ39hige86nhBP1LSc5JGG2OustauDEBPAPxk/qotuuGlDKf+97gEDevcIqA9pGXmKCU9W7n5BYqJjlJSYrzG9uQ++YC/HDPorbXzjTFnS3pCUoYxZpK19rHAtAagrlhr1S55jlPHNDlRC/9ygSIjAjvFp2XmKDk1SwWFRZKknPwCJadmSRJhD/jJcf+UW2v3WGtvkHS1pOnGmL3GmN3lf/m/TQA19cGKzRVC/p6LOuvz5KEBD3lJSknPdkK+VEFhkVLSswPeCxAuqnTVvTEmQdIDkr6T9IgkztcDQe7IKV6Ssh8Y4eq/i8/NL6jWOoDaO95V9/UlTZZ0p6R/SJporT0QiMYA1Nzc5Zt18yuLnXrKxZ117YB2LnZUIiY6SjmVhHpMdJQL3QDh4XgT/TeSmkoaaa2dH4B+ANRCcbFV+7sqTvFrHhgZNFfUJyXGVzhHL0lRkRFKSox3sSvA244X9MslTbDW7gpEMwBqbk5Wnm55dYlT3z+mi67uH+deQ5UoveCOq+6BwDneVfdXB6oRADVT2RT/3YMjXbnYrirG9owl2IEA4ha4QAj779Jc/f71TKd+8JKuurLv6S52BCDYEPRACAq1KR6Aewh6IMS8922O/vDGt0497f/O1hV92rrYEYBgRtADIaKo2OqMI6b4tQ+OVH2meADHQNADISAtM0e3v1k2xU+/tJsuT2jjYkeoDu7vDzcR9EAQY4oPfdzfH27jpwUQpN5evKlCyD96WXetnzaakA8x3N8fbmOiB4LM4aJidbj7fxXWvn9olCLqGZc6Qm1wf3+4jaAHgshb32zUX95Z5tSP/6q7LunZ2sWOUFvc3x9uI+iBIMAU713c3x9uI+gBl73x9QZN9F2cJUlPXtFDY3pwkZZXcH9/uI2gB1xSWFSsjkdM8T88NEr1mOI9h/v7w00EPeCCV7/6UXe/u9yp//abnrqoW4yLHQHwKoIeCCCmeACBRtADAfLyF+t1z3srnPofV/bSqLNbudcQgLBA0AN+duhwsc6cxBQPwB0EPeBHL3y2Tvf9d6VT//OqXhrRlSkeQOAQ9AiocNnc4+DhIsVPmlthbd3Do2QMUzyAwArYTbONMc8bY7YaY5aXW2tqjJlnjPnO9/EU37oxxjxljFlrjFlmjOkVqD7hP6Wbe+TkF8iqbHOPtMwct1urU899uq5CyM+4+hytnzaakAfgikDujvGipBFHrE2UNN9a21HSfF8tSSMldfT9Gi/pmQD1CD/y+uYeBw8XKW7i+7p/dtlb9eseHqXhXVq62BWAcBewoLfWLpK084jlMZJe8n3+kqSx5dZn2hJfSoo2xnBiM8R5eXOPfy36ocIU/9w1CUzxAIKC2+foW1hr8yTJWptnjDnNtx4raWO5523yreUd+Q2MMeNVMvWrbdu2/u0WteLFzT0OFBap0z2ciwcQvIJ1Y+vKfkrayp5orZ1hrU2w1iY0b97cz22hNpIS4xUVGVFhLZQ39/jnwu8rhPwL1/ZmigcQdNye6LcYY1r5pvlWkrb61jdJalPuea0l5Qa8O9Qpr2zuwRQPIJS4HfSzJF0jaZrv43vl1icYY96Q1FfSrtK3+BHaQn1zj79/vLbCxYMvXd9H55/JO0kAglfAgt4Y87qkwZKaGWM2SZqskoB/yxhzg6QNki7zPX2OpFGS1kraL+m6QPUJVKbgUJHOupcpHkDoCVjQW2t/fZSHhlbyXCvpVv92BFTNU/O/02Pz1jj1yzf00cCOTPEAQoPbb90DQWv/ocPqfG96hTWmeAChhqAHKvH4vDV6cv53Tv3qjX01oEMzFzsCgJoh6IFy9h08rC6TmeIBeAdBD/g8+kG2nv5orVO/dlNfnXsGUzyA0EbQI+ztPXhYXctN8cZI6x4e7WJHAFB3CHqEtb/OXa1nFnzv1G+O76e+7U91sSMAqFsEPcLSngOFOnvKB04dGWH03YOjXOwIAPyDoEfYeXjOKj276Aen/s/N/dU7rqmLHQGA/xD0CBu7DxSqW7kp/qQGEVo5dYSLHQGA/xH0CAsPzF6pf3+6zqnf+V1/nXM6UzwA7yPo4Wm7CgrV/b6yKb7RifWVNSXRxY4AILAIenjWff9doRc+W+/Uqbecq15tT3GvIQBwAUEPz9m1v1Ddp5ZN8fWM9NjlPQh5AGGJoIenTH5vuV764scKa8VWSk7NkiSN7RnrRlsA4BqCHp6Qv/+Qekydd9THCwqLlJKeTdADCDsEPULe3e9m6dWvNhz3ebn5BQHoBgCCSz23GwBqaue+Q4qb+L4T8q2anKj100YrNjqq0ufHHGUdALyMoEdImvjOMvW6v+yt+tm/P09fJA+VJCUlxisqMqLC86MiI5SUGB/QHgEgGPDWPULKjr0Hdc4DHzp161Oi9OmdQyo8p/Q8fEp6tnLzCxQTHaWkxHjOzwMISwQ9QkbSf5bqP4s3OfX7t52nLjFNKn3u2J6xBDsAiKBHCNi+96ASyk3xcaeepAVJF7jYEQCEDoIeQe1Pb36r1Mwcp557+0B1atnYxY4AILQQ9AhK2/YcVO8Hy6b4M5o31Pw7BrvXEACEKIIeQee21zM1a2muU6ffPkjxLRu52BEAhC6CHkFj6+4D6vPQfKfu1LKR5t4+yMWOACD0EfQICre+tkTvL8tz6nl/HKSOLZjiAaC2CHq4asvuA+pbborvGttYs38/0MWOAMBbCHq45uaXF2vuis1O/eGfzleH0052sSMA8B6CHgGXt6tA/R/+yKm7t26i9yac52JHAOBdBD0C6saXMvThqi1O/dEd56t9c6Z4APAXgh4BkZtfoHOnlU3xvdpGK/WWAS52BADhgaCH313/4jf6aPVWp17w58GKa9bQxY4AIHwQ9PCbTT/t13l//dip+8Q11Vs393exIwAIPwQ9/GLc819r0ZptTr0wabBOP5UpHgACjaBHndq4c78GTi+b4vu3P1Wvj+/nYkcAEN4IetSZK//9pT5bu8OpP/nLBWrT9CQXOwIAEPSotQ079mtQStkUf16HZnrlxr4udgQAKEXQo1Z+9ewX+mrdTqdmigeA4ELQo0bWb9+nwY8scOrB8c314nV93GsoSKVl5iglPVu5+QWKiY5SUmK8xvaMdbstAGGEoEe1XfrM58r48Sen/mziEMVGR7nYUXBKy8xRcmqWCgqLJEk5+QVKTs2SJMIeQMDUc7sBhI512/cpbuL7TsgP7XSa1k8bTcgfRUp6thPypQoKi5SSnu1SRwDCERM9qmTs3z/TtxvznfqL5CFq1YSAP5bc/IJqrQOAPxD0OKbvt+3V0EcXOnVilxZ69uoEFzsKHTHRUcqpJNRjeAcEQAAR9Diqi57+RMtzdjv1l8lD1bLJiS52FFqSEuMrnKOXpKjICCUlxrvYFYBwQ9DjZ9Zu3aNhjy1y6pFdW+qZq85xsaPQVHrBHVfdA3ATQY8KRj75iVbllU3xX981VKc1ZoqvqbE9Ywl2AK4i6CFJWrNlj4Y/XjbFX9w9Rk//uqeLHQEA6gJBDw1/fKHWbNnr1F/fPVSnNWKKBwAvIOjD2OrNuzXiiU+cemyPGD1xBVM8AHgJQR+mhjyyQD9s3+fU39w9TM0bneBiRwAAfyDow8zK3N0a9VTZFP/LXq316OXdXewIAOBPBH0YGTT9Y23Yud+pMyYNU7OTmeIBwMsI+jCwq6BQiY8v0ubdByRJlye01vRLmeIBIBwQ9EGuttucfrR6i5JTs7Rl90HVM1LGpAvVtGEDP3YMAAgmBH0Qq802p7v2F+q+2SuUuiRH8S0a6d/jeuvs1k383jMAILgQ9EHsWNucHivoP1y5RXe9m6Ud+w7p90M6aMKQDjqhfoS/2wUABCGCPohVd5vT/P2HNPW/K5WamaNOLRvp+Wt7q2ssUzwAhDOCPohVZ5vTeb4p/qd9h3Tb0I6acEEHNahfLxBtAgCCWFAkgTFmvTEmyxjzrTEmw7fW1Bgzzxjzne/jKW73GWhJifGKiqz4lvuR25z+tO+Qbn8jUzfNzNCpDRso7dYB+tOFZxLyAABJwTXRX2Ct3V6unihpvrV2mjFmoq++07Zm4V8AAAu2SURBVJ3W3HG8bU7TV2zW3e8uV/7+Q/rD0I66lSkeAHCEYAr6I42RNNj3+UuSFijMgl6qfJvTn/Yd0uRZKzRraa7OatVYL13fW11iOBcPAPi5YAl6K+kDY4yV9Ky1doakFtbaPEmy1uYZY05ztcMgMXf5Zk1Ky1L+/kL9cdiZuuWCMxQZwRQPAKhcsAT9AGttri/M5xljVlf1C40x4yWNl6S2bdv6qz/X7fRN8f9dmqsuMY018/q+6hzT2O22AABBLiiC3lqb6/u41RjzrqQ+krYYY1r5pvlWkrYe5WtnSJohSQkJCTZQPQfS3OV5mpS2XLsKCvWnC8/U7wYzxQMAqsb1tDDGNDTGNCr9XNJwScslzZJ0je9p10h6z50O3bNj70FNeG2Jbn5liVo0PlGzJpyn24Z2JOQBAFUWDBN9C0nvGmOkkn5es9bONcZ8I+ktY8wNkjZIuszFHgNuTlae7klbrt0HCvXn4Wfqt+czxQMAqs/1oLfW/iDpZ1upWWt3SBoa+I7ctX3vQU1+b4Xez8rT2bFN9OplfdWpJefiAQA143rQo8zsZbm6970V2nvgsJIS4/XbQe1VnykeAFALBH0Q2L73oO59b7nmZG1Wt9ZNlHJpd8W3bOR2WwAADyDoXWSt1exlebr3veXad7BIfxkRr/EDmeIBAHWHoHfJtj0HdU/acs1dsVndWzfRI5d1V8cWTPEAgLpF0AeYtVazluZq8qwV2n+oSBNHdtKN57VjigcA+AVBH0Bb9xzQpHeX64OVW9SjTbQeuaybOpzGFA8A8B+CPgCOnOKTR3bSjQPbK6Kecbs1AIDHEfR+tnX3Ad2dtlzzVm5Rz7bRSrm0uzqcdrLbbQEAwgRB7yfWWqV9m6Mps1bqQGGR7h51lq4/rx1TPAAgoAh6P9i6+4Duene5Ply1ReecfoqmX9pNZzRnigcABB5BX4estXo3M0dTZq3QwcPFmjT6LF03gCkeAOAegr6ObNl9QHelZmn+6q1K8E3x7ZniAQAuI+hryVqrd5bkaOp/V+hQUbHuuaizrj03jikeABAUCPpa2LzrgJJTl+nj7G3qHXeKpl/aXe2aNXS7LQAAHAR9DVhr9fbiTZo6e6UKi4o1+eLOuqZ/nOqF8BSflpmjlPRs5eYXKCY6SkmJ8RrbM9bttgAAtUTQV1PergIlp2ZpQfY29WnXVNN/2U1xIT7Fp2XmKDk1SwWFRZKknPyS/0ZJhD0AhDiCvoqstfpPxibdP3ulDhdbTbm4s8aF+BRfKiU92wn5UgWFRUpJzyboASDEEfRVkJtfoImpWVq0Zpv6tmuq6Zd20+mnhvYUX15ufkG11gEAoYOgPwZrrd78ZqMeeH+Vioqtpo7poqv6nu6JKb68mOgo5VQS6jHRUS50AwCoS+yNehS5+QUa9/zXmpiapa6xjZV++yDPvFV/pKTEeEVFRlRYi4qMUFJivEsdAQDqChP9UeTvL9SyTbt0/9iuurJPW08GfKnS8/BcdQ8A3mOstW73UGcSEhJsRkZGnX2/fQcPq+EJ/F0IABC8jDGLrbUJR3uct+6PgZAHAIQ6gh4AAA8j6AEA8DCCHgAADyPoAQDwMIIeAAAPI+gBAPAwgh4AAA8j6AEA8DDuCHMUaZk53BIWABDyCPpKpGXmKDk1y9mjPSe/QMmpWZJE2AMAQgpv3VciJT3bCflSBYVFSknPdqkjAABqhqCvRG4le7Mfax0AgGBF0FciJjqqWusAAAQrgr4SSYnxioqMqLAWFRmhpMR4lzoCAKBmuBivEqUX3HHVPQAg1BH0RzG2ZyzBDgAIebx1DwCAhxH0AAB4GEEPAICHEfQAAHgYQQ8AgIdx1X0YYIMeAAhfBL3HsUEPAIQ33rr3ODboAYDwRtB7HBv0AEB4I+g9jg16ACC8EfQexwY9ABDeuBjP49igBwDCG0EfBtigBwDCF2/dAwDgYQQ9AAAeRtADAOBhBD0AAB5G0AMA4GEEPQAAHhb0QW+MGWGMyTbGrDXGTHS7HwAAQklQB70xJkLS3yWNlNRZ0q+NMZ3d7QoAgNAR1EEvqY+ktdbaH6y1hyS9IWmMyz0BABAygj3oYyVtLFdv8q0BAIAqCPZb4JpK1myFJxgzXtJ4X7nXGHOsjdabSdpeR72hDK+rf/C6+g+vrX/wuvrH8V7X04/1xcEe9JsktSlXt5aUW/4J1toZkmZU5ZsZYzKstQl11x4kXld/4XX1H15b/+B19Y/avq7B/tb9N5I6GmPaGWMaSLpC0iyXewIAIGQE9URvrT1sjJkgKV1ShKTnrbUrXG4LAICQEdRBL0nW2jmS5tTRt6vSW/yoNl5X/+B19R9eW//gdfWPWr2uxlp7/GcBAICQFOzn6AEAQC2ERdBzG926YYxpY4z52BizyhizwhjzB996U2PMPGPMd76Pp7jdaygyxkQYYzKNMbN9dTtjzFe+1/VN3wWpqCZjTLQx5m1jzGrfsdufY7b2jDF/9P0cWG6Med0YcyLHbM0YY543xmw1xiwvt1bpMWpKPOXLs2XGmF7H+/6eD3puo1unDku6w1p7lqR+km71vZYTJc231naUNN9Xo/r+IGlVufqvkh73va4/SbrBla5C35OS5lprO0nqrpLXmGO2FowxsZJuk5Rgre2qkoulrxDHbE29KGnEEWtHO0ZHSuro+zVe0jPH++aeD3pxG906Y63Ns9Yu8X2+RyU/MGNV8nq+5HvaS5LGutNh6DLGtJY0WtK/fbWRNETS276n8LrWgDGmsaRBkp6TJGvtIWttvjhm60J9SVHGmPqSTpKUJ47ZGrHWLpK084jlox2jYyTNtCW+lBRtjGl1rO8fDkHPbXT9wBgTJ6mnpK8ktbDW5kklfxmQdJp7nYWsJyT9RVKxrz5VUr619rCv5ritmfaStkl6wXda5N/GmIbimK0Va22OpEckbVBJwO+StFgcs3XpaMdotTMtHIL+uLfRRfUYY06W9I6k2621u93uJ9QZYy6StNVau7j8ciVP5bitvvqSekl6xlrbU9I+8TZ9rfnOF4+R1E5SjKSGKnlL+Ugcs3Wv2j8bwiHoj3sbXVSdMSZSJSH/qrU21be8pfStI9/HrW71F6IGSPqFMWa9Sk4tDVHJhB/te1tU4ritqU2SNllrv/LVb6sk+Dlma2eYpHXW2m3W2kJJqZLOFcdsXTraMVrtTAuHoOc2unXEd974OUmrrLWPlXtolqRrfJ9fI+m9QPcWyqy1ydba1tbaOJUcnx9Za6+U9LGkS31P43WtAWvtZkkbjTHxvqWhklaKY7a2NkjqZ4w5yfdzofR15ZitO0c7RmdJGue7+r6fpF2lb/EfTVjcMMcYM0olE1LpbXQfdLmlkGSMOU/SJ5KyVHYu+S6VnKd/S1JblfwAuMxae+SFJagCY8xgSX+21l5kjGmvkgm/qaRMSVdZaw+62V8oMsb0UMlFjg0k/SDpOpUMORyztWCMuU/Sr1Tyr3EyJd2oknPFHLPVZIx5XdJglexSt0XSZElpquQY9f3F6m8quUp/v6TrrLUZx/z+4RD0AACEq3B46x4AgLBF0AMA4GEEPQAAHkbQAwDgYQQ9AAAeRtADAOBhBD2AWjHG9DTGFBljPnO7FwA/R9ADqK2bJP1DUldjzFluNwOgIm6YA6DGjDFRKtm9bJCkP0j6yVr7Z3e7AlAeEz2A2rhU0o/W2mWSXlbJPbgjXe4JQDkEPYDauFElAS9JC1Vy7+1fuNcOgCMR9ABqxBjTQSVb7L4mSbbkPOCrKgl/AEGi/vGfAgCVulElO0JuKNlQS5JkJMkY08Zau9GtxgCU4WI8ANVmjKkvaaOkJyXNPuLhlyW9a62dGvDGAPwMQQ+g2owxYyS9LamltXbHEY/dKel3ktpba4vd6A9AGc7RA6iJGyR9fGTI+/xH0umShgW2JQCVYaIHAMDDmOgBAPAwgh4AAA8j6AEA8DCCHgAADyPoAQDwMIIeAAAPI+gBAPAwgh4AAA8j6AEA8LD/B5P0WwAUAtg3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "\n",
    "ax.scatter(df3.A, df3.Y)\n",
    "ax.plot(df3.A, res.predict(df3[['constant', 'A']]))\n",
    "ax.set_xlabel('A', fontsize=14)\n",
    "ax.set_ylabel('Y', fontsize=14);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "(Still Program 11.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "\"Let us return to the data in Figure 11.1.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "df1['constant'] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "ols = sm.OLS(df1['Y'], df1[['constant', 'A']])\n",
    "res = ols.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>        <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>constant</th> <td>   67.5000</td> <td>   19.716</td> <td>    3.424</td> <td> 0.004</td> <td>   25.213</td> <td>  109.787</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>A</th>        <td>   78.7500</td> <td>   27.883</td> <td>    2.824</td> <td> 0.014</td> <td>   18.947</td> <td>  138.553</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.summary().tables[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E[Y|A=0] est.:  67.50\n",
      "E[Y|A=1] est.: 146.25\n"
     ]
    }
   ],
   "source": [
    "est_at_0 = res.params[0]\n",
    "est_at_1 = res.params[0] + res.params[1]\n",
    "print(\"E[Y|A=0] est.: {:>6.2f}\".format(est_at_0))\n",
    "print(\"E[Y|A=1] est.: {:>6.2f}\".format(est_at_1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### Program 11.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "Starting from the same data as Section 11.2, Program 11.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "df3['A^2'] = df3.A * df3.A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "ols = sm.OLS(df3.Y, df3[['constant', 'A', 'A^2']])\n",
    "res = ols.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>        <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>constant</th> <td>   -7.4069</td> <td>   31.748</td> <td>   -0.233</td> <td> 0.819</td> <td>  -75.994</td> <td>   61.180</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>A</th>        <td>    4.1072</td> <td>    1.531</td> <td>    2.683</td> <td> 0.019</td> <td>    0.800</td> <td>    7.414</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>A^2</th>      <td>   -0.0204</td> <td>    0.015</td> <td>   -1.331</td> <td> 0.206</td> <td>   -0.053</td> <td>    0.013</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "summary = res.summary()\n",
    "summary.tables[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfoAAAF3CAYAAABNO4lPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU1f3/8fchBEjYYtgTtgAxyh4Ii/tuoNWCiHVFXFGrP5daWrBf2367aaVVa60LiOIOLix+3VKlLlQrEAwSFsMqWYFASFgyZJmc3x8zwYBhCcnMnbnzej4ePJI5uUw+Xi9555x77jnGWisAAOBOzZwuAAAABA5BDwCAixH0AAC4GEEPAICLEfQAALgYQQ8AgIs1d7qAptSxY0fbu3dvp8sAACBoVqxYsdNa2+lIX3dV0Pfu3VuZmZlOlwEAQNAYY7Ye7esM3QMA4GIEPQAALkbQAwDgYgQ9AAAuRtADAOBiBD0AAC5G0AMA4GIEPQAALkbQAwDgYq5aGQ8AjmZhVoFmZOSosNSjhLgYTU1P0fjURKfLAgKKoAcQERZmFWj6/Gx5qrySpIJSj6bPz5Ykwh6uxtA9gIgwIyPnYMjX8lR5NSMjx6GKgOAg6AFEhMJST4PaAbcg6AFEhIS4mAa1A25B0AOICFPTUxQTHXVIW0x0lKampzhUERAcTMYDEBFqJ9wx6x6RhqAHEDHGpyYS7Ig4DN0DAOBiBD0AAC5G0AMA4GIEPQAALkbQAwDgYgQ9AAAuxuN1AAAEmJM7JxL0AAAEkNM7JzJ0DwBAADm9cyJBDwBAADm9cyJD9wAQBpy8x4vGSYiLUUE9oR6snRPp0QNAiKu9x1tQ6pHV9/d4F2YVOF0ajoPTOycS9AAQ4py+x4vGGZ+aqIcmDFJiXIyMpMS4GD00YRCz7gEAPk7f40XjOblzIj16AAhxR7qXG6x7vAhvBD0AhDin7/EivDF0DwAhrnbIl1n3OBEEPQCEASfv8SK8MXQPAICLEfQAALgYQQ8AgIsR9AAAuBhBDwCAixH0AAC4GEEPAICLEfQAALgYQQ8AgIsR9AAAuBhBDwCAixH0AAC4GEEPAICLEfQAALhY0ILeGNPDGPOJMWadMWaNMeYef3u8MeYjY8wG/8eT/O3GGPOEMWajMWaVMWZYsGoFABzZwqwCnfHwv5U07T2d8fC/tTCrwOmScBTB7NFXS7rfWnuqpNGS7jTG9Jc0TdJia22ypMX+15I0VlKy/88USU8HsVYAQD0WZhVo+vxsFZR6ZCUVlHo0fX42YR/Cghb01toia+3X/s/3SlonKVHSOEkv+g97UdJ4/+fjJL1kfb6SFGeM6RasegEAPzQjI0eeKu8hbZ4qr2Zk5DhUEY6luRPf1BjTW1KqpKWSulhriyTfLwPGmM7+wxIl5dX5a/n+tqLD3muKfD1+9ezZM6B1A0CkKyz1NKg9UllrVbyvQnkl5dq6q1y5Jb4/RaUH9Nqto2SMCVotQQ96Y0wbSW9Lutdau+co/7H1fcH+oMHamZJmSlJaWtoPvg4AaDoJcTEqqCfUE+JiHKjGWQeqvMrf7Q/xXeXKLfEot2T/wVA/UFVz8FhjpG7tWqlHfKz2VVSrbavooNUZ1KA3xkTLF/KvWmvn+5u3G2O6+Xvz3STt8LfnS+pR5693l1QYvGoBAIebmp6i6fOzDxm+j4mO0tT0FAerCoy6vXJfmHuUW1Lu66WX7Nf2PRWHHB/bIko942PVu0NrnZ3cST07xKpHfKx6xccq8aQYtWwe5ch/R9CC3vi67rMlrbPWPlrnS+9ImizpYf/HRXXa7zLGzJU0SlJZ7RA/AMAZ41MTJfnu1ReWepQQF6Op6SkH28ONr1fu8Q+x7/f3yssPhnvdX2iMkbr6e+VnJ3dSz/jYg2HeMz5WHVq3COqQ/PEy1gZntNsYc6akJZKyJdWOZzwg3336NyT1lJQr6QprbYn/F4MnJY2RVC7pRmtt5tG+R1pams3MPOohAIAIYq3Vzn2Vh4T31l3ff75tz4FDjq/tldeGd686QZ4YF6NW0c70yo/GGLPCWpt2pK8HrUdvrf2P6r/vLkkX1HO8lXRnQIsCAIS9ur3y3Dp/al+XVx76lEDXdq3Us0Oszkzu6OuV+4O9V4fQ7ZU3hiOz7gEAOF7WWu3a/32vvO4s9jx/r7zu4HRM9Pe98jP6HRrm3U8KzV55IBH0AADHVVT7euUHe+K7yrX1WL3y+Fid3tcf5B1i1DO+tXrGx6pjG/f1yhuDoAcABNW2sgNaumWXlm0p0cYd+5RXUq6iw3rlraKb+Xvirf1hHqOeHXw98+4nxUZcr7wxCHoAQEDllZRr2ZYSLd2yS0u3lGjrrnJJUtuWzXVKt7Ya3bfDwYlvtUPsndq0pFfeRAh6IEIszCpwzSNRCF3WWm3dVe4L9c0lWrql5OACO+1jojUyKV6TRvfS6D4ddGq3dopqRpgHGkEPRIDajUhqnwmu3YhEEmGPRrHWalPxPn3lD/VlW3YdXEimQ+sWGtUnXlPO7qORSfFK6dJWzQj2oCPogQhwtI1ICHo0VEGpR4vXbddXm3332Xfuq5QkdW7bUqP6dNCopHiN7hOvvp3aMPweAgh6IAKwEQkaa1PxPn24epsy1mzTqvwySVJiXIzOTu6kUX3iNSqpg3p1iCXYQxBBD0QANiJBQ1lrtaZwjzLWbNOHq7dpw459kqQhPeL0qzGnKH1AF/Xp1MbhKnE8CHogAkTSRiQ4cTU1Vl/n7taHq7fpwzXblL/bo2ZGGpXUQdeN7qWLB3RRt/b8chhuCHogArhtIxI0nSpvjf67aZc+XLNNH63druK9FWoR1UxnJnfU3ecn64JTO6tDm5ZOl4lGIOiBCDE+NZFghyTJU+nV5xuKlbF6mz5et117DlQrtkWUzkvprPSBXXVeSqeg7peOwCLoASAC7DlQpU++3aEPV2/TpznF8lR51T4mWhf176oxA7vqrOSOrDbnUgQ9ALjUzn0V+njtdn24Zpu+2LhTVV6rzm1bauLw7hozsKtGJsUrOqqZ02UiwAh6AHCRglKPMvyT6TK/K1GNlXrGx+rGM5KUPqCrUnvEsWhNhCHoASDM1feMe0qXtrrr/GSNGdBVp3Zry/PtEYygB4AwwzPuaAiCHgDCgLVWq/LL9H/fFPKMOxqEoAeAEFbmqdLCrAK9vixX327byzPuaDCCHgBCjLVWy7/brbnLcvVedpEqqms0MLGd/jh+oH4yNEHteMYdDUDQA0CIKNlfqbdX5Gvu8lxtKt6vNi2ba+Lw7rp6ZE8NTGzvdHkIUwQ9ADiopsbqy0279PryXP1rzTZVea2G9YzTIxMH65LB3RTbgh/TaByuIABwwI49B/TminzNW56n3JJytY+J1nWje+mqET2V0rWt0+XBRQh6AAgSb43V5+uL9dqyXP372x3y1liN7hOv+y8+WekDurIELQKCoAeAACvZX6mX/7tV85bnqrDsgDq2aaFbzkrSVSN6Kqlja6fLg8sR9AAQIHkl5Zq1ZLPeyMzTgaoanZXcUQ9e0l8XnNpFLZqzxjyCg6AHgCa2uqBMMz/frPeyi9TMSOOHJuq2c/qoX2fuvSP4CHoAaALW+mbPP/PZJi3ZsFNtWjbXzWcm6aYzktS1fSuny0MEI+gBoBG8NVYfrC7Ss59tVnZBmTq2aalfjknRtaN6qX0MC9vAeQQ9AJyAA1VevbkiX7M+36zcknIldWythyYM0mWpicyeR0gh6AGgAUrLfTPo53z5nXbtr9SQHnF64Een6KL+XRXFPu8IQQQ9AByHglKPZi/ZornLc1Ve6dV5KZ102zl9NSopnr3eEdIIegA4ipxte/XsZ5v0zjeFkqSfDEnQlHP66JSu7RyuDDg+BD0AHMZaq2VbSvTMZ5v0SU6xYltE6frTeuvms5KUGMee7wgvBD0A+NXUWP1r7XY989kmrcwrVYfWLXT/RSdr0mm9FBfbwunygBNC0AOIeBXVXi34ukAzP9+szTv3q2d8rP4wfqCuGN6dGfQIewQ9gIi150CVXv0qV89/sUXFeys0MLGdnrwmVWMGdFXzKJaohTsQ9AAizrayA3rhiy16dWmu9lVU66zkjnr8yqE6vW8HZtDDdQh6ABEjf3e5nli8QQuyCuStsfrx4ATddnYfDUxs73RpQMAQ9ABcr8xTpac+2agXvvxOknT1yJ665cw+6tkh1tnCgCAg6AG4VkW1V698lat//HuDyjxVuiw1Ub+4OEUJPCKHCELQA3Ada63eXVWkRzK+VV6JR2cld9S0sadoQAJD9Ig8BD0AV1m2pUR/en+dvskr1Sld2+rFm0bqnJM7OV0W4BiCHoArbNyxTw9/8K0+XrddXdu10iMTB+vyYd3ZaAYRj6AHENaK91bo8Y/Xa+7yPMVER2lqeopuOiNJMS1Y6AaQCHoAYaq8slrPLdmiZz/bpIrqGl07qqfuviBZHdu0dLo0IKQQ9ADCirfG6s3MPD360Xrt2FuhMQO66pdjUtSnUxunSwNCEkEPICxYa/VpTrEe+mCd1m/fp2E94/TUtcOU1jve6dKAkEbQA2FgYVaBZmTkqLDUo4S4GE1NT9H41ESnywqa1QVl+vP76/Tlpl3q1SFWT107TGMHdmW5WuA4EPRAiFuYVaDp87PlqfJKkgpKPZo+P1uSXB/2+bvL9bd/rdeCrAKdFBut313aX9eM6qUWzdlwBjheBD0Q4mZk5BwM+VqeKq9mZOS4NujrLllrJN1xbl/dcW5ftWsV7XRpQNgh6IEQV1jqaVB7OGPJWqDpEfRAiEuIi1FBPaHupvBjyVogcAh6IMRNTU855B69pIMLw7gBS9YCgUXQAyGu9j6822bds2QtEBwEPRAGxqcmhn2w1/JUevXYx+s1+z9bWLIWCIKgBb0x5nlJl0jaYa0d6G/7naRbJRX7D3vAWvu+/2vTJd0sySvpbmttRrBqBRAYX2zcqenzs5VbUq6rRvTQL9JTWLIWCLBg9ujnSHpS0kuHtT9mrf1r3QZjTH9JV0kaIClB0sfGmJOttV4BCDtl5VX643tr9eaKfCV1bK25U0ZrdJ8OTpcFRISgBb219nNjTO/jPHycpLnW2gpJW4wxGyWNlPTfAJUHIACstfpg9Tb9ZtEa7S6v1B3n9tU9FySrVTTD9ECwhMI9+ruMMddLypR0v7V2t6RESV/VOSbf3wYgTGwrO6AHF63WR2u3a2BiO825cYQGJvK4HBBsTq8j+bSkvpKGSiqS9Dd/e33Tbm19b2CMmWKMyTTGZBYXF9d3CIAgqqmxem1pri569DN9vr5Y08eeooU/O4OQBxziaI/eWru99nNjzCxJ7/pf5kvqUefQ7pIKj/AeMyXNlKS0tLR6fxkAEBybi/dp+vxsLd1SotP6dNBDEwapd8fWTpcFRDRHg94Y081aW+R/eZmk1f7P35H0mjHmUfkm4yVLWuZAiQCOQ5W3RrOWbNbjH29Qy+bN9JfLB+mnaT3YXQ4IAcF8vO51SedK6miMyZf0W0nnGmOGyjcs/52k2yTJWrvGGPOGpLWSqiXdyYx7IDRl55fpV2+v0tqiPRo7sKv+9ycD1LldK6fLAuBnrHXPaHdaWprNzMx0ugwgIngqvXr84/WatWSzOrZpqd+PG6gxA7s6XRYQcYwxK6y1aUf6eijMugcQZr7cuFPTF2Rr665yXT2yh6aNPVXtY9hCFghFBD2A41ZWXqU/v79O8zLz1LtDrF67dZRO79vR6bIAHAVBD+C4fJBdpN+8s0Yl+yt1+zl9de+FLHwDhAOCHsBRbd9zQL9ZtFoZa7ZrQEI7vXADC98A4YSgB1CvmhqreZl5+vP761RZXaNpY0/RLWcmqXmU0+tsAWgIgh7AD2zZuV/T56/SV5tLNLpPvB6aMFhJLHyDJrQwq0AzMnJUWOpRQlyMpqanuGYr5lBD0AM4qMpbo+eWbNHjH69Xi+bN9PCEQbpyBAvfoGktzCrQ9PnZ8lT5lkcpKPVo+vxsSSLsA4CgByBJWl3gW/hmTeEepQ/oot+PG6guLHyDAJiRkXMw5Gt5qryakZFD0AcAQQ9EOE+lV48vXq/nlmxRfOsWeua6YRozsJvTZcHFCks9DWpH4xD0QAT7ctNOTZ/vW/jmqhE9NH3sqWofy8I3CKyEuBgV1BPqCXExDlTjfkyfBSJQmadK095epWtmLZUkvXbLKD18+WBCHkExNT1FMYetwRATHaWp6SkOVeRu9OiBCPPh6iI9uMi38M1t5/TRvRecrJgWLHyD4Km9D8+s++Ag6IEIsXt/pR5YkK0PVm9T/24sfANnjU9NJNiDhKAHIsCXG3fqvjdWqmR/paamp2jK2X0UzcI3QEQg6AEXq6yu0aMfrdezn29SUsfWmj2ZXjwQaQh6wKW27Nyve+ZmaVV+ma4e2VMPXnKqYlvwTx6INPyrB1zGWqs3V+Trd++sUXRUM56LByIcQQ+4SJmnSg8syNZ7q4o0uk+8HrtyqLq159lkIJIR9IBLLNtSovvmrdT2PQf0yzEpuu3svopqxhr1QKQj6IEwV+2t0ROLN+jJTzaqR3ys3rrjdA3tEed0WQBCBEEPhLG8knLdMzdLX+eW6vJh3fW/4waoTUv+WQP4Hj8RgDC1aGWB/mfBaslIT1ydqp8MSXC6JAAhiKAHwszeA1X6zaI1WpBVoLReJ+mxK4eqR3ys02UBCFEEPRBGvs7drXvmZqlgt0f3Xpisu87rp+ascAfgKAh6IAx4a6ye/nSjHvt4g7q2a6U3bjtNab3jnS4LQBgg6IEQt2tfhe6dt1JLNuzUpUMS9KfLBqpdK7aTBXB8CHoghK3Yult3vfa1du2v1EMTBumqET1kDM/GAzh+BD0Qgqy1mvPld/rTe+uUEBej+XeczmY0AE4IQQ+EmL0HqjTt7Wy9l12kC0/tor/9dIjaxzBUD+DEEPRACMnZtld3vLJCW0vKNW3sKbrt7D4M1QNoFIIeCBHzv87XAwuy1bZVtF69ZZRG9+ngdEkAXICgBxx2oMqr37+7Vq8tzdWopHj94+pUdW7XyumyAmZhVoFmZOSosNSjhLgYTU1P0fjURKfLAlyLoAcclFdSrp+9+rWyC8p0+zl99YuLT3b1AjgLswo0fX62PFVeSVJBqUfT52dLEmEPBAhBDzhk8brtum/eSllJs65P00X9uzhdUsDNyMg5GPK1PFVezcjIIeiBACHogSCr9tbo0Y/W66lPN2lAQjs9fe1w9ewQGWvVF5Z6GtQOoPGOGvTGmNHW2q+CVQzgdsV7K3T361n67+ZdunpkD/320gFqFR3ldFlBkxAXo4J6Qj0hLsaBaoDIcKybgZ8bY/5gjKHnDzTSsi0l+vETS5SVt1t/vWKIHpowOKJCXpKmpqco5rD/5pjoKE1NT3GoIsD9jhX0YyVNkrTMGNM/CPUArmOt1czPN+nqWV+pdcvmWvCzMzRxeHeny3LE+NREPTRhkBLjYmQkJcbF6KEJg7g/DwTQUXvq1trFxphBkh6XlGmM+R9r7aPBKQ0IfweqvLr/zW/03qoijR3YVY9MHKy2Eb4hzfjURIIdCKJjDslba/dKutkY876kecaY30uqOeyYdgGqDwhbO/Yc0K0vr9Cq/FJWuQPgmOO6926MSZP0R0kbJP1VUnUgiwLC3ZrCMt3yYqZKy6v07HXDdfGArk6XBCBCHWvWfXNJv5X0K0lPSZpmrT0QjMKAcPXR2u26Z26W2sdE683bT2PXOQCOOlaPfrmkeEljrbWLg1APELastZq1ZLMe+uBbDU5sr1nXp7l6KVsA4eFYQb9a0l3W2rJgFAOEq8rqGj24cLXmZebpx4O66a9XDFFMi8h6dA5Hxvr+cNKxZt1PClYhQLgqLa/U7a+s0FebS3T3+f1074Unq1kzJt3Bh/X94TT37p4BBMHm4n267Kkv9fXWUj1+5VD9/OIUQh6HONr6/kAwsOIdcIK+2LhTd7yyQtFRzfTaraOU1jve6ZIQgljfH06jRw+cgNeW5mry88vUtX0rLbzzDEIeR3SkdfxZ3x/BQtADDeCtsfrDu2v1wIJsnZncUW/fcbp6xEfGznM4MazvD6cxdA8cp30V1br79Sz9+9sduuH03vqfH5+q5lH8royjq51wx6x7OIWgB45D/u5y3fJipjbs2Kc/jB+oSaN7OV0Swgjr+8NJBD1wDF/n7taUlzJVUV2jOTeO0FnJnZwuCQCOG0EPHMWilQWa+tYqdWvfSnOnjFC/zm2cLgkAGoSgB+phrdVjH2/QE4s3aGRSvJ69brhOat3C6bIAoMEIeuAwB6q8+sWb3+jdVUW6Ynh3/emyQWrRnEl3AMITQQ/UsWPvAd36EnvIA3APgh5BFcqbe6wt3KNbXlyu3eVVeua64UpnD3kALhC08UhjzPPGmB3GmNV12uKNMR8ZYzb4P57kbzfGmCeMMRuNMauMMcOCVScCp3Zzj4JSj6y+39xjYVaB06Xp47XbNfGZL1VjpTdvP42QB+AawbzxOEfSmMPapklabK1NlrTY/1qSxkpK9v+ZIunpINWIAArFzT2stZr1+Wbd+nKm+nVuo0V3naGBie0dqwcAmlrQgt5a+7mkksOax0l60f/5i5LG12l/yfp8JSnOGNMtOJUiUEJtc4/K6hpNeztbf3p/ncYO7Kp5U05Tl3atHKkFAALF6Xv0Xay1RZJkrS0yxnT2tydKyqtzXL6/rSjI9aEJJcTFqKCeUHdic4+6e8j/v/P76T72kAfgUqH6zFB9P3FtvQcaM8UYk2mMySwuLg5wWWiMUNnco+4e8o9dOUT3s4c8ABdzuke/3RjTzd+b7yZph789X1KPOsd1l1RY3xtYa2dKmilJaWlp9f4ygNAQCpt7sIc8gEjjdNC/I2mypIf9HxfVab/LGDNX0ihJZbVD/AhvTm7u8drSXP1m0Wr16dRasyePYHtZABEhaEFvjHld0rmSOhpj8iX9Vr6Af8MYc7OkXElX+A9/X9KPJG2UVC7pxmDVCffx1lj9+f11mv2fLTo3pZP+cXWq2raKdrosAAiKoAW9tfbqI3zpgnqOtZLuDGxFiATsIQ8g0jk9dA8EDHvIAwBBD5daW7hH1z+/TBXVXvaQBxDRCHq4TuZ3JbpxznK1adlcc6ecrn6d2zpdEgA4hqCHq3ySs0N3vLJCCe1j9PIto5TowGI8ABBKCHq4xjvfFOrn81YqpWtbvXjTSHVs09LpkgDAcQQ9XOGVr7bqwUWrNaJ3vJ6bnKZ2PD4HAJIIeoQ5a62e+nSTZmTk6IJTOuuf1w5Tq8OW2QWASEbQI2xZ61sIZ9aSLRo/NEEzrhiiaJ6RB4BDEPQIS9XeGj2wIFtvZOZr8mm99NtLB7AxDQDUg6BH2Kmo9uqe11fqwzXbdPcFybrvwmQZQ8gDQH0IeoSV/RXVmvJypr7YuEu/uaS/bjozyemSACCkEfQIG7v3V+rGOcuVXVCmv10xRJcP737EYxdmFTi6HS4AhAqCHmFh+54DmjR7qb7bVa5nrhuui/p3OeKxC7MKNH1+tjxVXklSQalH0+dnSxJhDyDiMEUZIe+7nft1+dNfqmC3R3NuHHHUkJekGRk5B0O+lqfKqxkZOYEsEwBCEj16hLR1RXs0afYyeWtq9PqU0RrcPe6Yf6ew1NOgdgBwM3r0CFkrtpboymf/q+goozdvP+24Ql6SEo6wvv2R2gHAzQh6hKTP1hfr2ueWqkOblnrz9tMatAPd1PQUxRy2Ol5MdJSmpqc0dZkAEPIYukfIeXdVoe6bt1LJndvqpZsbvjlN7YQ7Zt0DAEGPEPPa0lz9emG2RvSK13M3nPjmNONTEwl2ABBBjxDy9Keb9JcPv9V5KZ301LXDFdOCzWkAoLEIejjOWquHP/xWz362WeOGJuivbE4DAE2GoIejvDVWv16QrbnL83T9ab30OzanAYAmRdDDMRXVXt03b6Xez96mu8/vp/suOpnNaQCgiRH0cMT+imrd/soKLdmwUw9e0l83szkNAAQEQY+gKy33bU7zTV6pZkwcrCvSejhdEgC4FkGPoNq1r0LXPrdUm4v36+nrhit9QFenSwIAVyPoETS1If/drv16/oYROjO5o9MlAYDrEfQIipL9lbr2uaXastMX8mf0I+QBIBh4WBkBV7K/UtfM+oqQBwAHEPQIqN11evKzJxPyABBsBD0CZvf+Sl3z3FJtLt6n5yancU8eABzAPXoERG1PflPxPj13fZrOSu7kdEkAEJHo0aPJlZZX6rrZS7XRH/Jnn0zIA4BTCHo0qdJyX09+w459mkXIA4DjCHo0mdqe/IYd+zRz0nCdQ8gDgOMIejSJ2pBfv80X8uemdHa6JACACHo0gbLyqoMh/+z1hDwAhBKCHo1ySMhPGq7zCHkACCk8XocTVuap0qTnlypn2149M2mYzjuFkD/cwqwCzcjIUWGpRwlxMZqanqLxqYlOlwUgghD0OCFlnipNmr1U64r26Jnrhuv8U7o4XVLIWZhVoOnzs+Wp8kqSCko9mj4/W5IIewBBw9A9GqzMU6Xr64T8BacS8vWZkZFzMORreaq8mpGR41BFACIRQY8G2XOgStc/v0xri/bo6WsJ+aMpLPU0qB0AAoGgx3Hbc6BKk2Yv09rCMj117XBd2J+QP5qEuJgGtQNAIBD0OC57DlTpen/I//OaYbqIkD+mqekpiomOOqQtJjpKU9NTHKoIQCRiMh6Oae+BKk1+fpnW+EP+4gFdnS4pLNROuGPWPQAnEfQ4qr3+e/LZ+WV66lpCvqHGpyYS7AAcxdA9jqi2J5+dX6Z/EvIAEJYIetRr74Eq3fDCcq3KL9OT1wxTOiEPAGGJoMcP7Kuo1g0vLNc3eaV68ppUjRlIyANAuOIePQ6xr6Jak59fppV5pXry6lSNGdjN6ZIAAI1Ajx4H7auo1g11Qn7sIEIeAMIdQQ9JkqfSq5teWK6svFL9g5AHANdg6D7EBWP3s4pqr6a8nKnMrSX6+1Wp+hEhDwCuQdCHsGDsflbtrdHdr2dpyYademTiYF06JKFJ3hcAEBoYug9hgd79rKbGaupbq5SxZrt+d2l//TStR5O8LwAgdBD0ISyQu59Za/XgotVakFWgqekputIPmzYAAA5QSURBVOGMpEa/JwAg9BD0ISxQu59Za/XQB9/q1aW5uuPcvrrzvH6Nej8AQOgKiaA3xnxnjMk2xqw0xmT62+KNMR8ZYzb4P57kdJ3BFqjdz55YvFEzP9+syaf10i/ZSQ0AXC0kgt7vPGvtUGttmv/1NEmLrbXJkhb7X0eU8amJemjCICXGxchISoyL0UMTBjVqIt5zSzbrsY/Xa+Lw7vrtpQNkjGm6ggEAISeUZ92Pk3Su//MXJX0q6VdOFeOUptz97LWlufrje+v0o0Fd9fCEQWrWjJAHALcLlR69lfQvY8wKY8wUf1sXa22RJPk/dq7vLxpjphhjMo0xmcXFxUEqN/y8802hfr0wW+eldNLjV6aqeVSo/K8HAARSqPToz7DWFhpjOkv6yBjz7fH+RWvtTEkzJSktLc0GqsBw9tn6Yv183kqN6B2vp68brhbNCXkAiBQh8RPfWlvo/7hD0gJJIyVtN8Z0kyT/xx3OVRi+vs7drdtfXqGTu7TVc5PT1OqwyX0AAHdzPOiNMa2NMW1rP5d0saTVkt6RNNl/2GRJi5ypMHxt2L5XN81Zrs7tWurFm0aqXatop0sCAARZKAzdd5G0wD/7u7mk16y1Hxpjlkt6wxhzs6RcSVc4WGPYKSj16Prnlyk6qplevmmUOrVt6XRJAAAHOB701trNkobU075L0gXBryj87dpXoUmzl2pfRbXeuO009ewQ63RJAACHOD50j6a1r6JaN85ZroLdHs2ePEKndmvndEkAAAc53qNH06mo9ur2l1doTeEePXvdcI1Mine6JACAw+jRu4S3xurn877Rfzbu1F8uH6wL+3dxuiQAQAgg6F3AWqvfLFqt97KL9OsfnaqJw7s7XRIAIEQQ9C7w2Mcb9OrSXN12Th/denYfp8sBAIQQgj7Mzflii55YvEE/TeuuaWNOcbocAECIIejD2KKVBfrd/63Vxf276M+XDWInOgDADxD0Yeqz9cW6/41vNCopXk9czSY1AID6kQ5hqO769bNYvx4AcBQEfZhh/XoAQEMQ9GGE9esBAA1F0IeJuuvXv3TTSNavBwAcF4I+DLB+PQDgRLHWfYir8tbojldYvx4AcGII+hDmW9p2jZZs2KlHArx+/cKsAs3IyFFhqUcJcTGamp6i8amJAft+AIDgIOhD2Kwlm/X6slz97Ny++umIHgH7PguzCjR9frY8VV5Jvkl/0+dnSxJhDwBhjnv0IeqD7CL9+f1v9ePB3fSLi1MC+r1mZOQcDPlaniqvZmTkBPT7AgACj6APQSvzSnXvvJUa1jNOf7tiiJo1C+zStoWlnga1AwDCB0EfYvJKynXLi74FcWZdH5xV7xLiYhrUDgAIHwR9CCnzVOmmOctVWV2jF24YoQ5tgrMgztT0FMUc9gtFTHSUpqYH9pYBACDwmIwXIqq8Nbrz1a+1Zed+vXTzSPXr3DZo37t2wh2z7gHAfQj6EGCt1YMLV+s/G3dqxsTBOr1vx6DXMD41kWAHABdi6D4EPPv5Zs1dnqe7zuunK9IC9xgdACDyEPQOez+7SA9/8K0uHZKgn190stPlAABchqB30Ne5u3XfvJUa3uskzZg4OOCP0QEAIg9B75C8knJNeSlTXdq10sxJw4PyGB0AIPIQ9A4o81TpxjnLVeW1euHG4D1GBwCIPAR9kFV5a/SzV1do6679eua64erbqY3TJQEAXIzH64LIWqv/WbBaX2zcpb9eMUSn9e3gdEkAAJejRx9ET3+2SfMy83T3+f00cXh3p8sBAEQAevRH0NT7s7+7qlCPfJijnwxJ0H08RgcACBKCvh5NvT/7iq279fM3vlFar5P0yMTBMobH6AAAwcHQfT2acn/23F2+x+i6tW+lmUHajQ4AgFoEfT2aan/2svIq3ThnmaprrF64YYTiW7doivIAADhuBH09mmJ/9srqGt3+ygrllpTr2UnD1YfH6AAADiDo69HY/dmttfr1gmz9d/Mu/eXywRrdh8foAADOYDJePRq7P/tTn27SmyvydfcFyZowjMfoAADOIeiP4ET3Z/+/bwo1IyNH44Ym6L4LkwNQGQAAx4+h+ya0YmuJ7n/zG43ozWN0AIDQQNA3ka279uvWl1YooX0rzZyUppbNeYwOAOA8gr4J+B6jW64aa/XCjSN1Eo/RAQBCBEHfSJXVNbrtlUzll3g0c1Kakjq2drokAAAOYjJeI1hrNX1+tr7aXKLHrxyqkUnxTpcEAMAhCPpGePLfG/X21/m698LkRm14E2hNvUEPACB8EPQnaNHKAv3to/W6LDVR91wQuo/RNfUGPQCA8MI9+hOw/LsSTX1zlUYmxevhyweF9GN0TblBDwAg/BD0DfTdzv2a8lKmEk+K0bPXDQ/5x+iaaoMeAEB4IugboLS8UjfNWS5JeuGGEWHxGF1TbNADAAhfBP1xqqj2asrLK5S/26OZ16epd5g8RtfYDXoAAOGNyXjHwVqr6W9na9mWEv39qqEa0Tt8HqNr7AY9AIDwRtAfhycWb9T8rAL9/KKTNW5o+AXkiW7QAwAIfwzdH8PCrAI99vF6TRiWqP93fj+nywEAoEEI+qNYtqVEv3xrlUYlxevhCexGBwAIPwT9EWzZuV9TXs5U95Ni9Oyk4WrRnFMFAAg/3KM/Am+NVa8OrfXEVUMVFxv6j9EBAFAfgv4I+nVuo4U/O53hegBAWGM8+igIeQBAuAv5oDfGjDHG5BhjNhpjpjldDwAA4SSkg94YEyXpn5LGSuov6WpjTH9nqwIAIHyEdNBLGilpo7V2s7W2UtJcSeMcrgkAgLAR6kGfKCmvzut8f9tBxpgpxphMY0xmcXFxUIsDACDUhXrQ1zcbzh7ywtqZ1to0a21ap06dglQWAADhIdSDPl9Sjzqvu0sqdKgWAADCTqgH/XJJycaYJGNMC0lXSXrH4ZoAAAgbIb1gjrW22hhzl6QMSVGSnrfWrnG4LAAAwkZIB70kWWvfl/S+03UAABCOQn3oHgAANAJBDwCAixH0AAC4mLHWHvuoMGGMKZa09SiHdJS0M0jlRBLOa2BwXgOHcxsYnNfAONZ57WWtPeJCMq4K+mMxxmRaa9OcrsNtOK+BwXkNHM5tYHBeA6Ox55WhewAAXIygBwDAxSIt6Gc6XYBLcV4Dg/MaOJzbwOC8BkajzmtE3aMHACDSRFqPHgCAiBIRQW+MGWOMyTHGbDTGTHO6nnBljOlhjPnEGLPOGLPGGHOPvz3eGPORMWaD/+NJTtcajowxUcaYLGPMu/7XScaYpf7zOs+/sRMayBgTZ4x5yxjzrf/aPY1rtvGMMff5fw6sNsa8boxpxTV7YowxzxtjdhhjVtdpq/caNT5P+PNslTFm2LHe3/VBb4yJkvRPSWMl9Zd0tTGmv7NVha1qSfdba0+VNFrSnf5zOU3SYmttsqTF/tdouHskravz+i+SHvOf192SbnakqvD3d0kfWmtPkTREvnPMNdsIxphESXdLSrPWDpRv07GrxDV7ouZIGnNY25Gu0bGSkv1/pkh6+lhv7vqglzRS0kZr7WZrbaWkuZLGOVxTWLLWFllrv/Z/vle+H5iJ8p3PF/2HvShpvDMVhi9jTHdJP5b0nP+1kXS+pLf8h3BeT4Axpp2ksyXNliRrbaW1tlRcs02huaQYY0xzSbGSisQ1e0KstZ9LKjms+UjX6DhJL1mfryTFGWO6He39IyHoEyXl1Xmd729DIxhjektKlbRUUhdrbZHk+2VAUmfnKgtbj0v6paQa/+sOkkqttdX+11y3J6aPpGJJL/hvizxnjGktrtlGsdYWSPqrpFz5Ar5M0gpxzTalI12jDc60SAh6U08bjxo0gjGmjaS3Jd1rrd3jdD3hzhhziaQd1toVdZvrOZTrtuGaSxom6Wlrbaqk/WKYvtH894vHSUqSlCCptXxDyofjmm16Df7ZEAlBny+pR53X3SUVOlRL2DPGRMsX8q9aa+f7m7fXDh35P+5wqr4wdYaknxhjvpPv1tL58vXw4/zDohLX7YnKl5RvrV3qf/2WfMHPNds4F0raYq0tttZWSZov6XRxzTalI12jDc60SAj65ZKS/bNBW8g3YeQdh2sKS/77xrMlrbPWPlrnS+9Imuz/fLKkRcGuLZxZa6dba7tba3vLd33+21p7raRPJE30H8Z5PQHW2m2S8owxKf6mCyStFddsY+VKGm2MifX/XKg9r1yzTedI1+g7kq73z74fLamsdoj/SCJiwRxjzI/k6yFFSXreWvsnh0sKS8aYMyUtkZSt7+8lPyDfffo3JPWU7wfAFdbawyeW4DgYY86V9Atr7SXGmD7y9fDjJWVJus5aW+FkfeHIGDNUvkmOLSRtlnSjfJ0crtlGMMb8r6Qr5XsaJ0vSLfLdK+aabSBjzOuSzpVvl7rtkn4raaHquUb9v1g9Kd8s/XJJN1prM4/6/pEQ9AAARKpIGLoHACBiEfQAALgYQQ8AgIsR9AAAuBhBDwCAixH0AAC4GEEPoFGMManGGK8x5gunawHwQwQ9gMa6VdJTkgYaY051uhgAh2LBHAAnzBgTI9/uZWdLukfSbmvtL5ytCkBd9OgBNMZESVuttaskvSzfGtzRDtcEoA6CHkBj3CJfwEvSZ/Ktvf0T58oBcDiCHsAJMcb0k2+L3dckyfruA74qX/gDCBHNj30IANTrFvl2hMz1baglSTKSZIzpYa3Nc6owAN9jMh6ABjPGNJeUJ+nvkt497MsvS1pgrf190AsD8AMEPYAGM8aMk/SWpK7W2l2Hfe1Xku6Q1MdaW+NEfQC+xz16ACfiZkmfHB7yfm9K6iXpwuCWBKA+9OgBAHAxevQAALgYQQ8AgIsR9AAAuBhBDwCAixH0AAC4GEEPAICLEfQAALgYQQ8AgIsR9AAAuNj/B0azvqcWmvf0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df3.sort_values('A', inplace=True)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "\n",
    "ax.scatter(df3.A, df3.Y)\n",
    "ax.plot(df3.A, res.predict(df3[['constant', 'A', 'A^2']]))\n",
    "ax.set_xlabel('A', fontsize=14)\n",
    "ax.set_ylabel('Y', fontsize=14);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "The expected value and confidence interval can be obtained with the method contributed by @pettypychen above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "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>mean</th>\n",
       "      <th>mean_se</th>\n",
       "      <th>mean_ci_lower</th>\n",
       "      <th>mean_ci_upper</th>\n",
       "      <th>obs_ci_lower</th>\n",
       "      <th>obs_ci_upper</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>0</td>\n",
       "      <td>197.13</td>\n",
       "      <td>25.16</td>\n",
       "      <td>142.77</td>\n",
       "      <td>251.49</td>\n",
       "      <td>89.63</td>\n",
       "      <td>304.62</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     mean  mean_se  mean_ci_lower  mean_ci_upper  obs_ci_lower  obs_ci_upper\n",
       "0  197.13    25.16         142.77         251.49         89.63        304.62"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pred = res.get_prediction(exog=[1, 90, 90**2])\n",
    "pred.summary_frame().round(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false",
    "toc-hr-collapsed": false
   },
   "source": [
    "## 解读知识点\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### 估计置信区间"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import scipy.stats\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>        <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>constant</th> <td>   -7.4069</td> <td>   31.748</td> <td>   -0.233</td> <td> 0.819</td> <td>  -75.994</td> <td>   61.180</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>A</th>        <td>    4.1072</td> <td>    1.531</td> <td>    2.683</td> <td> 0.019</td> <td>    0.800</td> <td>    7.414</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>A^2</th>      <td>   -0.0204</td> <td>    0.015</td> <td>   -1.331</td> <td> 0.206</td> <td>   -0.053</td> <td>    0.013</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A, Y = zip(*(\n",
    "    (3, 21),\n",
    "    (11, 54),\n",
    "    (17, 33),\n",
    "    (23, 101),\n",
    "    (29, 85),\n",
    "    (37, 65),\n",
    "    (41, 157),\n",
    "    (53, 120),\n",
    "    (67, 111),\n",
    "    (79, 200),\n",
    "    (83, 140),\n",
    "    (97, 220),\n",
    "    (60, 230),\n",
    "    (71, 217),\n",
    "    (15, 11),\n",
    "    (45, 190),\n",
    "))\n",
    "\n",
    "df3 = pd.DataFrame({'A': A, 'Y': Y, 'constant': np.ones(16)})\n",
    "\n",
    "df3['A^2'] = df3.A * df3.A\n",
    "\n",
    "ols = sm.OLS(df3.Y, df3[['constant', 'A', 'A^2']])\n",
    "res = ols.fit()\n",
    "\n",
    "summary = res.summary()\n",
    "summary.tables[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "Collapsed": "false",
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfoAAAF3CAYAAABNO4lPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU1f3/8fchBEjYYtgTtgAxyh4Ii/tuoNWCiHVFXFGrP5daWrBf2367aaVVa60LiOIOLix+3VKlLlQrEAwSFsMqWYFASFgyZJmc3x8zwYBhCcnMnbnzej4ePJI5uUw+Xi9555x77jnGWisAAOBOzZwuAAAABA5BDwCAixH0AAC4GEEPAICLEfQAALgYQQ8AgIs1d7qAptSxY0fbu3dvp8sAACBoVqxYsdNa2+lIX3dV0Pfu3VuZmZlOlwEAQNAYY7Ye7esM3QMA4GIEPQAALkbQAwDgYgQ9AAAuRtADAOBiBD0AAC5G0AMA4GIEPQAALkbQAwDgYq5aGQ8AjmZhVoFmZOSosNSjhLgYTU1P0fjURKfLAgKKoAcQERZmFWj6/Gx5qrySpIJSj6bPz5Ykwh6uxtA9gIgwIyPnYMjX8lR5NSMjx6GKgOAg6AFEhMJST4PaAbcg6AFEhIS4mAa1A25B0AOICFPTUxQTHXVIW0x0lKampzhUERAcTMYDEBFqJ9wx6x6RhqAHEDHGpyYS7Ig4DN0DAOBiBD0AAC5G0AMA4GIEPQAALkbQAwDgYgQ9AAAuxuN1AAAEmJM7JxL0AAAEkNM7JzJ0DwBAADm9cyJBDwBAADm9cyJD9wAQBpy8x4vGSYiLUUE9oR6snRPp0QNAiKu9x1tQ6pHV9/d4F2YVOF0ajoPTOycS9AAQ4py+x4vGGZ+aqIcmDFJiXIyMpMS4GD00YRCz7gEAPk7f40XjOblzIj16AAhxR7qXG6x7vAhvBD0AhDin7/EivDF0DwAhrnbIl1n3OBEEPQCEASfv8SK8MXQPAICLEfQAALgYQQ8AgIsR9AAAuBhBDwCAixH0AAC4GEEPAICLEfQAALgYQQ8AgIsR9AAAuBhBDwCAixH0AAC4GEEPAICLEfQAALhY0ILeGNPDGPOJMWadMWaNMeYef3u8MeYjY8wG/8eT/O3GGPOEMWajMWaVMWZYsGoFABzZwqwCnfHwv5U07T2d8fC/tTCrwOmScBTB7NFXS7rfWnuqpNGS7jTG9Jc0TdJia22ypMX+15I0VlKy/88USU8HsVYAQD0WZhVo+vxsFZR6ZCUVlHo0fX42YR/Cghb01toia+3X/s/3SlonKVHSOEkv+g97UdJ4/+fjJL1kfb6SFGeM6RasegEAPzQjI0eeKu8hbZ4qr2Zk5DhUEY6luRPf1BjTW1KqpKWSulhriyTfLwPGmM7+wxIl5dX5a/n+tqLD3muKfD1+9ezZM6B1A0CkKyz1NKg9UllrVbyvQnkl5dq6q1y5Jb4/RaUH9Nqto2SMCVotQQ96Y0wbSW9Lutdau+co/7H1fcH+oMHamZJmSlJaWtoPvg4AaDoJcTEqqCfUE+JiHKjGWQeqvMrf7Q/xXeXKLfEot2T/wVA/UFVz8FhjpG7tWqlHfKz2VVSrbavooNUZ1KA3xkTLF/KvWmvn+5u3G2O6+Xvz3STt8LfnS+pR5693l1QYvGoBAIebmp6i6fOzDxm+j4mO0tT0FAerCoy6vXJfmHuUW1Lu66WX7Nf2PRWHHB/bIko942PVu0NrnZ3cST07xKpHfKx6xccq8aQYtWwe5ch/R9CC3vi67rMlrbPWPlrnS+9ImizpYf/HRXXa7zLGzJU0SlJZ7RA/AMAZ41MTJfnu1ReWepQQF6Op6SkH28ONr1fu8Q+x7/f3yssPhnvdX2iMkbr6e+VnJ3dSz/jYg2HeMz5WHVq3COqQ/PEy1gZntNsYc6akJZKyJdWOZzwg3336NyT1lJQr6QprbYn/F4MnJY2RVC7pRmtt5tG+R1pams3MPOohAIAIYq3Vzn2Vh4T31l3ff75tz4FDjq/tldeGd686QZ4YF6NW0c70yo/GGLPCWpt2pK8HrUdvrf2P6r/vLkkX1HO8lXRnQIsCAIS9ur3y3Dp/al+XVx76lEDXdq3Us0Oszkzu6OuV+4O9V4fQ7ZU3hiOz7gEAOF7WWu3a/32vvO4s9jx/r7zu4HRM9Pe98jP6HRrm3U8KzV55IBH0AADHVVT7euUHe+K7yrX1WL3y+Fid3tcf5B1i1DO+tXrGx6pjG/f1yhuDoAcABNW2sgNaumWXlm0p0cYd+5RXUq6iw3rlraKb+Xvirf1hHqOeHXw98+4nxUZcr7wxCHoAQEDllZRr2ZYSLd2yS0u3lGjrrnJJUtuWzXVKt7Ya3bfDwYlvtUPsndq0pFfeRAh6IEIszCpwzSNRCF3WWm3dVe4L9c0lWrql5OACO+1jojUyKV6TRvfS6D4ddGq3dopqRpgHGkEPRIDajUhqnwmu3YhEEmGPRrHWalPxPn3lD/VlW3YdXEimQ+sWGtUnXlPO7qORSfFK6dJWzQj2oCPogQhwtI1ICHo0VEGpR4vXbddXm3332Xfuq5QkdW7bUqP6dNCopHiN7hOvvp3aMPweAgh6IAKwEQkaa1PxPn24epsy1mzTqvwySVJiXIzOTu6kUX3iNSqpg3p1iCXYQxBBD0QANiJBQ1lrtaZwjzLWbNOHq7dpw459kqQhPeL0qzGnKH1AF/Xp1MbhKnE8CHogAkTSRiQ4cTU1Vl/n7taHq7fpwzXblL/bo2ZGGpXUQdeN7qWLB3RRt/b8chhuCHogArhtIxI0nSpvjf67aZc+XLNNH63druK9FWoR1UxnJnfU3ecn64JTO6tDm5ZOl4lGIOiBCDE+NZFghyTJU+nV5xuKlbF6mz5et117DlQrtkWUzkvprPSBXXVeSqeg7peOwCLoASAC7DlQpU++3aEPV2/TpznF8lR51T4mWhf176oxA7vqrOSOrDbnUgQ9ALjUzn0V+njtdn24Zpu+2LhTVV6rzm1bauLw7hozsKtGJsUrOqqZ02UiwAh6AHCRglKPMvyT6TK/K1GNlXrGx+rGM5KUPqCrUnvEsWhNhCHoASDM1feMe0qXtrrr/GSNGdBVp3Zry/PtEYygB4AwwzPuaAiCHgDCgLVWq/LL9H/fFPKMOxqEoAeAEFbmqdLCrAK9vixX327byzPuaDCCHgBCjLVWy7/brbnLcvVedpEqqms0MLGd/jh+oH4yNEHteMYdDUDQA0CIKNlfqbdX5Gvu8lxtKt6vNi2ba+Lw7rp6ZE8NTGzvdHkIUwQ9ADiopsbqy0279PryXP1rzTZVea2G9YzTIxMH65LB3RTbgh/TaByuIABwwI49B/TminzNW56n3JJytY+J1nWje+mqET2V0rWt0+XBRQh6AAgSb43V5+uL9dqyXP372x3y1liN7hOv+y8+WekDurIELQKCoAeAACvZX6mX/7tV85bnqrDsgDq2aaFbzkrSVSN6Kqlja6fLg8sR9AAQIHkl5Zq1ZLPeyMzTgaoanZXcUQ9e0l8XnNpFLZqzxjyCg6AHgCa2uqBMMz/frPeyi9TMSOOHJuq2c/qoX2fuvSP4CHoAaALW+mbPP/PZJi3ZsFNtWjbXzWcm6aYzktS1fSuny0MEI+gBoBG8NVYfrC7Ss59tVnZBmTq2aalfjknRtaN6qX0MC9vAeQQ9AJyAA1VevbkiX7M+36zcknIldWythyYM0mWpicyeR0gh6AGgAUrLfTPo53z5nXbtr9SQHnF64Een6KL+XRXFPu8IQQQ9AByHglKPZi/ZornLc1Ve6dV5KZ102zl9NSopnr3eEdIIegA4ipxte/XsZ5v0zjeFkqSfDEnQlHP66JSu7RyuDDg+BD0AHMZaq2VbSvTMZ5v0SU6xYltE6frTeuvms5KUGMee7wgvBD0A+NXUWP1r7XY989kmrcwrVYfWLXT/RSdr0mm9FBfbwunygBNC0AOIeBXVXi34ukAzP9+szTv3q2d8rP4wfqCuGN6dGfQIewQ9gIi150CVXv0qV89/sUXFeys0MLGdnrwmVWMGdFXzKJaohTsQ9AAizrayA3rhiy16dWmu9lVU66zkjnr8yqE6vW8HZtDDdQh6ABEjf3e5nli8QQuyCuStsfrx4ATddnYfDUxs73RpQMAQ9ABcr8xTpac+2agXvvxOknT1yJ665cw+6tkh1tnCgCAg6AG4VkW1V698lat//HuDyjxVuiw1Ub+4OEUJPCKHCELQA3Ada63eXVWkRzK+VV6JR2cld9S0sadoQAJD9Ig8BD0AV1m2pUR/en+dvskr1Sld2+rFm0bqnJM7OV0W4BiCHoArbNyxTw9/8K0+XrddXdu10iMTB+vyYd3ZaAYRj6AHENaK91bo8Y/Xa+7yPMVER2lqeopuOiNJMS1Y6AaQCHoAYaq8slrPLdmiZz/bpIrqGl07qqfuviBZHdu0dLo0IKQQ9ADCirfG6s3MPD360Xrt2FuhMQO66pdjUtSnUxunSwNCEkEPICxYa/VpTrEe+mCd1m/fp2E94/TUtcOU1jve6dKAkEbQA2FgYVaBZmTkqLDUo4S4GE1NT9H41ESnywqa1QVl+vP76/Tlpl3q1SFWT107TGMHdmW5WuA4EPRAiFuYVaDp87PlqfJKkgpKPZo+P1uSXB/2+bvL9bd/rdeCrAKdFBut313aX9eM6qUWzdlwBjheBD0Q4mZk5BwM+VqeKq9mZOS4NujrLllrJN1xbl/dcW5ftWsV7XRpQNgh6IEQV1jqaVB7OGPJWqDpEfRAiEuIi1FBPaHupvBjyVogcAh6IMRNTU855B69pIMLw7gBS9YCgUXQAyGu9j6822bds2QtEBwEPRAGxqcmhn2w1/JUevXYx+s1+z9bWLIWCIKgBb0x5nlJl0jaYa0d6G/7naRbJRX7D3vAWvu+/2vTJd0sySvpbmttRrBqBRAYX2zcqenzs5VbUq6rRvTQL9JTWLIWCLBg9ujnSHpS0kuHtT9mrf1r3QZjTH9JV0kaIClB0sfGmJOttV4BCDtl5VX643tr9eaKfCV1bK25U0ZrdJ8OTpcFRISgBb219nNjTO/jPHycpLnW2gpJW4wxGyWNlPTfAJUHIACstfpg9Tb9ZtEa7S6v1B3n9tU9FySrVTTD9ECwhMI9+ruMMddLypR0v7V2t6RESV/VOSbf3wYgTGwrO6AHF63WR2u3a2BiO825cYQGJvK4HBBsTq8j+bSkvpKGSiqS9Dd/e33Tbm19b2CMmWKMyTTGZBYXF9d3CIAgqqmxem1pri569DN9vr5Y08eeooU/O4OQBxziaI/eWru99nNjzCxJ7/pf5kvqUefQ7pIKj/AeMyXNlKS0tLR6fxkAEBybi/dp+vxsLd1SotP6dNBDEwapd8fWTpcFRDRHg94Y081aW+R/eZmk1f7P35H0mjHmUfkm4yVLWuZAiQCOQ5W3RrOWbNbjH29Qy+bN9JfLB+mnaT3YXQ4IAcF8vO51SedK6miMyZf0W0nnGmOGyjcs/52k2yTJWrvGGPOGpLWSqiXdyYx7IDRl55fpV2+v0tqiPRo7sKv+9ycD1LldK6fLAuBnrHXPaHdaWprNzMx0ugwgIngqvXr84/WatWSzOrZpqd+PG6gxA7s6XRYQcYwxK6y1aUf6eijMugcQZr7cuFPTF2Rr665yXT2yh6aNPVXtY9hCFghFBD2A41ZWXqU/v79O8zLz1LtDrF67dZRO79vR6bIAHAVBD+C4fJBdpN+8s0Yl+yt1+zl9de+FLHwDhAOCHsBRbd9zQL9ZtFoZa7ZrQEI7vXADC98A4YSgB1CvmhqreZl5+vP761RZXaNpY0/RLWcmqXmU0+tsAWgIgh7AD2zZuV/T56/SV5tLNLpPvB6aMFhJLHyDJrQwq0AzMnJUWOpRQlyMpqanuGYr5lBD0AM4qMpbo+eWbNHjH69Xi+bN9PCEQbpyBAvfoGktzCrQ9PnZ8lT5lkcpKPVo+vxsSSLsA4CgByBJWl3gW/hmTeEepQ/oot+PG6guLHyDAJiRkXMw5Gt5qryakZFD0AcAQQ9EOE+lV48vXq/nlmxRfOsWeua6YRozsJvTZcHFCks9DWpH4xD0QAT7ctNOTZ/vW/jmqhE9NH3sqWofy8I3CKyEuBgV1BPqCXExDlTjfkyfBSJQmadK095epWtmLZUkvXbLKD18+WBCHkExNT1FMYetwRATHaWp6SkOVeRu9OiBCPPh6iI9uMi38M1t5/TRvRecrJgWLHyD4Km9D8+s++Ag6IEIsXt/pR5YkK0PVm9T/24sfANnjU9NJNiDhKAHIsCXG3fqvjdWqmR/paamp2jK2X0UzcI3QEQg6AEXq6yu0aMfrdezn29SUsfWmj2ZXjwQaQh6wKW27Nyve+ZmaVV+ma4e2VMPXnKqYlvwTx6INPyrB1zGWqs3V+Trd++sUXRUM56LByIcQQ+4SJmnSg8syNZ7q4o0uk+8HrtyqLq159lkIJIR9IBLLNtSovvmrdT2PQf0yzEpuu3svopqxhr1QKQj6IEwV+2t0ROLN+jJTzaqR3ys3rrjdA3tEed0WQBCBEEPhLG8knLdMzdLX+eW6vJh3fW/4waoTUv+WQP4Hj8RgDC1aGWB/mfBaslIT1ydqp8MSXC6JAAhiKAHwszeA1X6zaI1WpBVoLReJ+mxK4eqR3ys02UBCFEEPRBGvs7drXvmZqlgt0f3Xpisu87rp+ascAfgKAh6IAx4a6ye/nSjHvt4g7q2a6U3bjtNab3jnS4LQBgg6IEQt2tfhe6dt1JLNuzUpUMS9KfLBqpdK7aTBXB8CHoghK3Yult3vfa1du2v1EMTBumqET1kDM/GAzh+BD0Qgqy1mvPld/rTe+uUEBej+XeczmY0AE4IQQ+EmL0HqjTt7Wy9l12kC0/tor/9dIjaxzBUD+DEEPRACMnZtld3vLJCW0vKNW3sKbrt7D4M1QNoFIIeCBHzv87XAwuy1bZVtF69ZZRG9+ngdEkAXICgBxx2oMqr37+7Vq8tzdWopHj94+pUdW7XyumyAmZhVoFmZOSosNSjhLgYTU1P0fjURKfLAlyLoAcclFdSrp+9+rWyC8p0+zl99YuLT3b1AjgLswo0fX62PFVeSVJBqUfT52dLEmEPBAhBDzhk8brtum/eSllJs65P00X9uzhdUsDNyMg5GPK1PFVezcjIIeiBACHogSCr9tbo0Y/W66lPN2lAQjs9fe1w9ewQGWvVF5Z6GtQOoPGOGvTGmNHW2q+CVQzgdsV7K3T361n67+ZdunpkD/320gFqFR3ldFlBkxAXo4J6Qj0hLsaBaoDIcKybgZ8bY/5gjKHnDzTSsi0l+vETS5SVt1t/vWKIHpowOKJCXpKmpqco5rD/5pjoKE1NT3GoIsD9jhX0YyVNkrTMGNM/CPUArmOt1czPN+nqWV+pdcvmWvCzMzRxeHeny3LE+NREPTRhkBLjYmQkJcbF6KEJg7g/DwTQUXvq1trFxphBkh6XlGmM+R9r7aPBKQ0IfweqvLr/zW/03qoijR3YVY9MHKy2Eb4hzfjURIIdCKJjDslba/dKutkY876kecaY30uqOeyYdgGqDwhbO/Yc0K0vr9Cq/FJWuQPgmOO6926MSZP0R0kbJP1VUnUgiwLC3ZrCMt3yYqZKy6v07HXDdfGArk6XBCBCHWvWfXNJv5X0K0lPSZpmrT0QjMKAcPXR2u26Z26W2sdE683bT2PXOQCOOlaPfrmkeEljrbWLg1APELastZq1ZLMe+uBbDU5sr1nXp7l6KVsA4eFYQb9a0l3W2rJgFAOEq8rqGj24cLXmZebpx4O66a9XDFFMi8h6dA5Hxvr+cNKxZt1PClYhQLgqLa/U7a+s0FebS3T3+f1074Unq1kzJt3Bh/X94TT37p4BBMHm4n267Kkv9fXWUj1+5VD9/OIUQh6HONr6/kAwsOIdcIK+2LhTd7yyQtFRzfTaraOU1jve6ZIQgljfH06jRw+cgNeW5mry88vUtX0rLbzzDEIeR3SkdfxZ3x/BQtADDeCtsfrDu2v1wIJsnZncUW/fcbp6xEfGznM4MazvD6cxdA8cp30V1br79Sz9+9sduuH03vqfH5+q5lH8royjq51wx6x7OIWgB45D/u5y3fJipjbs2Kc/jB+oSaN7OV0Swgjr+8NJBD1wDF/n7taUlzJVUV2jOTeO0FnJnZwuCQCOG0EPHMWilQWa+tYqdWvfSnOnjFC/zm2cLgkAGoSgB+phrdVjH2/QE4s3aGRSvJ69brhOat3C6bIAoMEIeuAwB6q8+sWb3+jdVUW6Ynh3/emyQWrRnEl3AMITQQ/UsWPvAd36EnvIA3APgh5BFcqbe6wt3KNbXlyu3eVVeua64UpnD3kALhC08UhjzPPGmB3GmNV12uKNMR8ZYzb4P57kbzfGmCeMMRuNMauMMcOCVScCp3Zzj4JSj6y+39xjYVaB06Xp47XbNfGZL1VjpTdvP42QB+AawbzxOEfSmMPapklabK1NlrTY/1qSxkpK9v+ZIunpINWIAArFzT2stZr1+Wbd+nKm+nVuo0V3naGBie0dqwcAmlrQgt5a+7mkksOax0l60f/5i5LG12l/yfp8JSnOGNMtOJUiUEJtc4/K6hpNeztbf3p/ncYO7Kp5U05Tl3atHKkFAALF6Xv0Xay1RZJkrS0yxnT2tydKyqtzXL6/rSjI9aEJJcTFqKCeUHdic4+6e8j/v/P76T72kAfgUqH6zFB9P3FtvQcaM8UYk2mMySwuLg5wWWiMUNnco+4e8o9dOUT3s4c8ABdzuke/3RjTzd+b7yZph789X1KPOsd1l1RY3xtYa2dKmilJaWlp9f4ygNAQCpt7sIc8gEjjdNC/I2mypIf9HxfVab/LGDNX0ihJZbVD/AhvTm7u8drSXP1m0Wr16dRasyePYHtZABEhaEFvjHld0rmSOhpj8iX9Vr6Af8MYc7OkXElX+A9/X9KPJG2UVC7pxmDVCffx1lj9+f11mv2fLTo3pZP+cXWq2raKdrosAAiKoAW9tfbqI3zpgnqOtZLuDGxFiATsIQ8g0jk9dA8EDHvIAwBBD5daW7hH1z+/TBXVXvaQBxDRCHq4TuZ3JbpxznK1adlcc6ecrn6d2zpdEgA4hqCHq3ySs0N3vLJCCe1j9PIto5TowGI8ABBKCHq4xjvfFOrn81YqpWtbvXjTSHVs09LpkgDAcQQ9XOGVr7bqwUWrNaJ3vJ6bnKZ2PD4HAJIIeoQ5a62e+nSTZmTk6IJTOuuf1w5Tq8OW2QWASEbQI2xZ61sIZ9aSLRo/NEEzrhiiaJ6RB4BDEPQIS9XeGj2wIFtvZOZr8mm99NtLB7AxDQDUg6BH2Kmo9uqe11fqwzXbdPcFybrvwmQZQ8gDQH0IeoSV/RXVmvJypr7YuEu/uaS/bjozyemSACCkEfQIG7v3V+rGOcuVXVCmv10xRJcP737EYxdmFTi6HS4AhAqCHmFh+54DmjR7qb7bVa5nrhuui/p3OeKxC7MKNH1+tjxVXklSQalH0+dnSxJhDyDiMEUZIe+7nft1+dNfqmC3R3NuHHHUkJekGRk5B0O+lqfKqxkZOYEsEwBCEj16hLR1RXs0afYyeWtq9PqU0RrcPe6Yf6ew1NOgdgBwM3r0CFkrtpboymf/q+goozdvP+24Ql6SEo6wvv2R2gHAzQh6hKTP1hfr2ueWqkOblnrz9tMatAPd1PQUxRy2Ol5MdJSmpqc0dZkAEPIYukfIeXdVoe6bt1LJndvqpZsbvjlN7YQ7Zt0DAEGPEPPa0lz9emG2RvSK13M3nPjmNONTEwl2ABBBjxDy9Keb9JcPv9V5KZ301LXDFdOCzWkAoLEIejjOWquHP/xWz362WeOGJuivbE4DAE2GoIejvDVWv16QrbnL83T9ab30OzanAYAmRdDDMRXVXt03b6Xez96mu8/vp/suOpnNaQCgiRH0cMT+imrd/soKLdmwUw9e0l83szkNAAQEQY+gKy33bU7zTV6pZkwcrCvSejhdEgC4FkGPoNq1r0LXPrdUm4v36+nrhit9QFenSwIAVyPoETS1If/drv16/oYROjO5o9MlAYDrEfQIipL9lbr2uaXastMX8mf0I+QBIBh4WBkBV7K/UtfM+oqQBwAHEPQIqN11evKzJxPyABBsBD0CZvf+Sl3z3FJtLt6n5yancU8eABzAPXoERG1PflPxPj13fZrOSu7kdEkAEJHo0aPJlZZX6rrZS7XRH/Jnn0zIA4BTCHo0qdJyX09+w459mkXIA4DjCHo0mdqe/IYd+zRz0nCdQ8gDgOMIejSJ2pBfv80X8uemdHa6JACACHo0gbLyqoMh/+z1hDwAhBKCHo1ySMhPGq7zCHkACCk8XocTVuap0qTnlypn2149M2mYzjuFkD/cwqwCzcjIUWGpRwlxMZqanqLxqYlOlwUgghD0OCFlnipNmr1U64r26Jnrhuv8U7o4XVLIWZhVoOnzs+Wp8kqSCko9mj4/W5IIewBBw9A9GqzMU6Xr64T8BacS8vWZkZFzMORreaq8mpGR41BFACIRQY8G2XOgStc/v0xri/bo6WsJ+aMpLPU0qB0AAoGgx3Hbc6BKk2Yv09rCMj117XBd2J+QP5qEuJgGtQNAIBD0OC57DlTpen/I//OaYbqIkD+mqekpiomOOqQtJjpKU9NTHKoIQCRiMh6Oae+BKk1+fpnW+EP+4gFdnS4pLNROuGPWPQAnEfQ4qr3+e/LZ+WV66lpCvqHGpyYS7AAcxdA9jqi2J5+dX6Z/EvIAEJYIetRr74Eq3fDCcq3KL9OT1wxTOiEPAGGJoMcP7Kuo1g0vLNc3eaV68ppUjRlIyANAuOIePQ6xr6Jak59fppV5pXry6lSNGdjN6ZIAAI1Ajx4H7auo1g11Qn7sIEIeAMIdQQ9JkqfSq5teWK6svFL9g5AHANdg6D7EBWP3s4pqr6a8nKnMrSX6+1Wp+hEhDwCuQdCHsGDsflbtrdHdr2dpyYademTiYF06JKFJ3hcAEBoYug9hgd79rKbGaupbq5SxZrt+d2l//TStR5O8LwAgdBD0ISyQu59Za/XgotVakFWgqekputIPmzYAAA5QSURBVOGMpEa/JwAg9BD0ISxQu59Za/XQB9/q1aW5uuPcvrrzvH6Nej8AQOgKiaA3xnxnjMk2xqw0xmT62+KNMR8ZYzb4P57kdJ3BFqjdz55YvFEzP9+syaf10i/ZSQ0AXC0kgt7vPGvtUGttmv/1NEmLrbXJkhb7X0eU8amJemjCICXGxchISoyL0UMTBjVqIt5zSzbrsY/Xa+Lw7vrtpQNkjGm6ggEAISeUZ92Pk3Su//MXJX0q6VdOFeOUptz97LWlufrje+v0o0Fd9fCEQWrWjJAHALcLlR69lfQvY8wKY8wUf1sXa22RJPk/dq7vLxpjphhjMo0xmcXFxUEqN/y8802hfr0wW+eldNLjV6aqeVSo/K8HAARSqPToz7DWFhpjOkv6yBjz7fH+RWvtTEkzJSktLc0GqsBw9tn6Yv183kqN6B2vp68brhbNCXkAiBQh8RPfWlvo/7hD0gJJIyVtN8Z0kyT/xx3OVRi+vs7drdtfXqGTu7TVc5PT1OqwyX0AAHdzPOiNMa2NMW1rP5d0saTVkt6RNNl/2GRJi5ypMHxt2L5XN81Zrs7tWurFm0aqXatop0sCAARZKAzdd5G0wD/7u7mk16y1Hxpjlkt6wxhzs6RcSVc4WGPYKSj16Prnlyk6qplevmmUOrVt6XRJAAAHOB701trNkobU075L0gXBryj87dpXoUmzl2pfRbXeuO009ewQ63RJAACHOD50j6a1r6JaN85ZroLdHs2ePEKndmvndEkAAAc53qNH06mo9ur2l1doTeEePXvdcI1Mine6JACAw+jRu4S3xurn877Rfzbu1F8uH6wL+3dxuiQAQAgg6F3AWqvfLFqt97KL9OsfnaqJw7s7XRIAIEQQ9C7w2Mcb9OrSXN12Th/denYfp8sBAIQQgj7Mzflii55YvEE/TeuuaWNOcbocAECIIejD2KKVBfrd/63Vxf276M+XDWInOgDADxD0Yeqz9cW6/41vNCopXk9czSY1AID6kQ5hqO769bNYvx4AcBQEfZhh/XoAQEMQ9GGE9esBAA1F0IeJuuvXv3TTSNavBwAcF4I+DLB+PQDgRLHWfYir8tbojldYvx4AcGII+hDmW9p2jZZs2KlHArx+/cKsAs3IyFFhqUcJcTGamp6i8amJAft+AIDgIOhD2Kwlm/X6slz97Ny++umIHgH7PguzCjR9frY8VV5Jvkl/0+dnSxJhDwBhjnv0IeqD7CL9+f1v9ePB3fSLi1MC+r1mZOQcDPlaniqvZmTkBPT7AgACj6APQSvzSnXvvJUa1jNOf7tiiJo1C+zStoWlnga1AwDCB0EfYvJKynXLi74FcWZdH5xV7xLiYhrUDgAIHwR9CCnzVOmmOctVWV2jF24YoQ5tgrMgztT0FMUc9gtFTHSUpqYH9pYBACDwmIwXIqq8Nbrz1a+1Zed+vXTzSPXr3DZo37t2wh2z7gHAfQj6EGCt1YMLV+s/G3dqxsTBOr1vx6DXMD41kWAHABdi6D4EPPv5Zs1dnqe7zuunK9IC9xgdACDyEPQOez+7SA9/8K0uHZKgn190stPlAABchqB30Ne5u3XfvJUa3uskzZg4OOCP0QEAIg9B75C8knJNeSlTXdq10sxJw4PyGB0AIPIQ9A4o81TpxjnLVeW1euHG4D1GBwCIPAR9kFV5a/SzV1do6679eua64erbqY3TJQEAXIzH64LIWqv/WbBaX2zcpb9eMUSn9e3gdEkAAJejRx9ET3+2SfMy83T3+f00cXh3p8sBAEQAevRH0NT7s7+7qlCPfJijnwxJ0H08RgcACBKCvh5NvT/7iq279fM3vlFar5P0yMTBMobH6AAAwcHQfT2acn/23F2+x+i6tW+lmUHajQ4AgFoEfT2aan/2svIq3ThnmaprrF64YYTiW7doivIAADhuBH09mmJ/9srqGt3+ygrllpTr2UnD1YfH6AAADiDo69HY/dmttfr1gmz9d/Mu/eXywRrdh8foAADOYDJePRq7P/tTn27SmyvydfcFyZowjMfoAADOIeiP4ET3Z/+/bwo1IyNH44Ym6L4LkwNQGQAAx4+h+ya0YmuJ7n/zG43ozWN0AIDQQNA3ka279uvWl1YooX0rzZyUppbNeYwOAOA8gr4J+B6jW64aa/XCjSN1Eo/RAQBCBEHfSJXVNbrtlUzll3g0c1Kakjq2drokAAAOYjJeI1hrNX1+tr7aXKLHrxyqkUnxTpcEAMAhCPpGePLfG/X21/m698LkRm14E2hNvUEPACB8EPQnaNHKAv3to/W6LDVR91wQuo/RNfUGPQCA8MI9+hOw/LsSTX1zlUYmxevhyweF9GN0TblBDwAg/BD0DfTdzv2a8lKmEk+K0bPXDQ/5x+iaaoMeAEB4IugboLS8UjfNWS5JeuGGEWHxGF1TbNADAAhfBP1xqqj2asrLK5S/26OZ16epd5g8RtfYDXoAAOGNyXjHwVqr6W9na9mWEv39qqEa0Tt8HqNr7AY9AIDwRtAfhycWb9T8rAL9/KKTNW5o+AXkiW7QAwAIfwzdH8PCrAI99vF6TRiWqP93fj+nywEAoEEI+qNYtqVEv3xrlUYlxevhCexGBwAIPwT9EWzZuV9TXs5U95Ni9Oyk4WrRnFMFAAg/3KM/Am+NVa8OrfXEVUMVFxv6j9EBAFAfgv4I+nVuo4U/O53hegBAWGM8+igIeQBAuAv5oDfGjDHG5BhjNhpjpjldDwAA4SSkg94YEyXpn5LGSuov6WpjTH9nqwIAIHyEdNBLGilpo7V2s7W2UtJcSeMcrgkAgLAR6kGfKCmvzut8f9tBxpgpxphMY0xmcXFxUIsDACDUhXrQ1zcbzh7ywtqZ1to0a21ap06dglQWAADhIdSDPl9Sjzqvu0sqdKgWAADCTqgH/XJJycaYJGNMC0lXSXrH4ZoAAAgbIb1gjrW22hhzl6QMSVGSnrfWrnG4LAAAwkZIB70kWWvfl/S+03UAABCOQn3oHgAANAJBDwCAixH0AAC4mLHWHvuoMGGMKZa09SiHdJS0M0jlRBLOa2BwXgOHcxsYnNfAONZ57WWtPeJCMq4K+mMxxmRaa9OcrsNtOK+BwXkNHM5tYHBeA6Ox55WhewAAXIygBwDAxSIt6Gc6XYBLcV4Dg/MaOJzbwOC8BkajzmtE3aMHACDSRFqPHgCAiBIRQW+MGWOMyTHGbDTGTHO6nnBljOlhjPnEGLPOGLPGGHOPvz3eGPORMWaD/+NJTtcajowxUcaYLGPMu/7XScaYpf7zOs+/sRMayBgTZ4x5yxjzrf/aPY1rtvGMMff5fw6sNsa8boxpxTV7YowxzxtjdhhjVtdpq/caNT5P+PNslTFm2LHe3/VBb4yJkvRPSWMl9Zd0tTGmv7NVha1qSfdba0+VNFrSnf5zOU3SYmttsqTF/tdouHskravz+i+SHvOf192SbnakqvD3d0kfWmtPkTREvnPMNdsIxphESXdLSrPWDpRv07GrxDV7ouZIGnNY25Gu0bGSkv1/pkh6+lhv7vqglzRS0kZr7WZrbaWkuZLGOVxTWLLWFllrv/Z/vle+H5iJ8p3PF/2HvShpvDMVhi9jTHdJP5b0nP+1kXS+pLf8h3BeT4Axpp2ksyXNliRrbaW1tlRcs02huaQYY0xzSbGSisQ1e0KstZ9LKjms+UjX6DhJL1mfryTFGWO6He39IyHoEyXl1Xmd729DIxhjektKlbRUUhdrbZHk+2VAUmfnKgtbj0v6paQa/+sOkkqttdX+11y3J6aPpGJJL/hvizxnjGktrtlGsdYWSPqrpFz5Ar5M0gpxzTalI12jDc60SAh6U08bjxo0gjGmjaS3Jd1rrd3jdD3hzhhziaQd1toVdZvrOZTrtuGaSxom6Wlrbaqk/WKYvtH894vHSUqSlCCptXxDyofjmm16Df7ZEAlBny+pR53X3SUVOlRL2DPGRMsX8q9aa+f7m7fXDh35P+5wqr4wdYaknxhjvpPv1tL58vXw4/zDohLX7YnKl5RvrV3qf/2WfMHPNds4F0raYq0tttZWSZov6XRxzTalI12jDc60SAj65ZKS/bNBW8g3YeQdh2sKS/77xrMlrbPWPlrnS+9Imuz/fLKkRcGuLZxZa6dba7tba3vLd33+21p7raRPJE30H8Z5PQHW2m2S8owxKf6mCyStFddsY+VKGm2MifX/XKg9r1yzTedI1+g7kq73z74fLamsdoj/SCJiwRxjzI/k6yFFSXreWvsnh0sKS8aYMyUtkZSt7+8lPyDfffo3JPWU7wfAFdbawyeW4DgYY86V9Atr7SXGmD7y9fDjJWVJus5aW+FkfeHIGDNUvkmOLSRtlnSjfJ0crtlGMMb8r6Qr5XsaJ0vSLfLdK+aabSBjzOuSzpVvl7rtkn4raaHquUb9v1g9Kd8s/XJJN1prM4/6/pEQ9AAARKpIGLoHACBiEfQAALgYQQ8AgIsR9AAAuBhBDwCAixH0AAC4GEEPoFGMManGGK8x5gunawHwQwQ9gMa6VdJTkgYaY051uhgAh2LBHAAnzBgTI9/uZWdLukfSbmvtL5ytCkBd9OgBNMZESVuttaskvSzfGtzRDtcEoA6CHkBj3CJfwEvSZ/Ktvf0T58oBcDiCHsAJMcb0k2+L3dckyfruA74qX/gDCBHNj30IANTrFvl2hMz1baglSTKSZIzpYa3Nc6owAN9jMh6ABjPGNJeUJ+nvkt497MsvS1pgrf190AsD8AMEPYAGM8aMk/SWpK7W2l2Hfe1Xku6Q1MdaW+NEfQC+xz16ACfiZkmfHB7yfm9K6iXpwuCWBKA+9OgBAHAxevQAALgYQQ8AgIsR9AAAuBhBDwCAixH0AAC4GEEPAICLEfQAALgYQQ8AgIsR9AAAuNj/B0azvqcWmvf0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df3.sort_values('A', inplace=True)\n",
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "ax.scatter(df3.A, df3.Y)\n",
    "ax.plot(df3.A, res.predict(df3[['constant', 'A', 'A^2']]))\n",
    "ax.set_xlabel('A', fontsize=14)\n",
    "ax.set_ylabel('Y', fontsize=14);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "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>mean</th>\n",
       "      <th>mean_se</th>\n",
       "      <th>mean_ci_lower</th>\n",
       "      <th>mean_ci_upper</th>\n",
       "      <th>obs_ci_lower</th>\n",
       "      <th>obs_ci_upper</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>0</td>\n",
       "      <td>197.13</td>\n",
       "      <td>25.16</td>\n",
       "      <td>142.77</td>\n",
       "      <td>251.49</td>\n",
       "      <td>89.63</td>\n",
       "      <td>304.62</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     mean  mean_se  mean_ci_lower  mean_ci_upper  obs_ci_lower  obs_ci_upper\n",
       "0  197.13    25.16         142.77         251.49         89.63        304.62"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pred = res.get_prediction(exog=[1, 90, 90**2])\n",
    "pred.summary_frame().round(2)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "165px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  },
  "toc-autonumbering": false
 },
 "nbformat": 4,
 "nbformat_minor": 4
}