Approximating a Bayesian Posterior Prediction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "99604b98",
   "metadata": {},
   "source": [
    "### From Gelman & Hill: Data Analysis Using Regression and Multilevel Models (Ch. 7)\n",
    "\n",
    "The idea here is that, before going full Bayesian, we can get a sense for how Bayesian Posterior Predictions work by creating a normal linear regression model and then approximating a Bayesian Posterior Predictive Distribution by sampling from our predictions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6cfcd30f",
   "metadata": {},
   "outputs": [],
   "source": [
    "## import libraries\n",
    "\n",
    "from sklearn import datasets\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import statsmodels.api as smf\n",
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "16504594",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sepal_length</th>\n",
       "      <th>sepal_width</th>\n",
       "      <th>petal_length</th>\n",
       "      <th>petal_width</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>5.1</td>\n",
       "      <td>3.5</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4.9</td>\n",
       "      <td>3.0</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>4.7</td>\n",
       "      <td>3.2</td>\n",
       "      <td>1.3</td>\n",
       "      <td>0.2</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   sepal_length  sepal_width  petal_length  petal_width\n",
       "0           5.1          3.5           1.4          0.2\n",
       "1           4.9          3.0           1.4          0.2\n",
       "2           4.7          3.2           1.3          0.2"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## iris data set\n",
    "iris = datasets.load_iris()\n",
    "\n",
    "## convert to pandas sata frame\n",
    "data = iris['data']\n",
    "target = iris['target']\n",
    "\n",
    "iris_df = pd.DataFrame(data)\n",
    "iris_df.columns = ['sepal_length', 'sepal_width', 'petal_length', 'petal_width']\n",
    "iris_df.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b17b6b87",
   "metadata": {},
   "source": [
    "### Build a linear regression (traditional ordinary least squares)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "9af500d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Set up our X and y variables\n",
    "X = iris_df['petal_length']\n",
    "y = iris_df['petal_width']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "3b4dbe65",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda3/lib/python3.7/site-packages/statsmodels/tsa/tsatools.py:142: FutureWarning: In a future version of pandas all arguments of concat except for the argument 'objs' will be keyword-only\n",
      "  x = pd.concat(x[::order], 1)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>       <td>petal_width</td>   <th>  R-squared:         </th> <td>   0.927</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.927</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   1882.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Sun, 16 Oct 2022</td> <th>  Prob (F-statistic):</th> <td>4.68e-86</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>18:49:58</td>     <th>  Log-Likelihood:    </th> <td>  24.796</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>   150</td>      <th>  AIC:               </th> <td>  -45.59</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>   148</td>      <th>  BIC:               </th> <td>  -39.57</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<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>const</th>        <td>   -0.3631</td> <td>    0.040</td> <td>   -9.131</td> <td> 0.000</td> <td>   -0.442</td> <td>   -0.285</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>petal_length</th> <td>    0.4158</td> <td>    0.010</td> <td>   43.387</td> <td> 0.000</td> <td>    0.397</td> <td>    0.435</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 5.765</td> <th>  Durbin-Watson:     </th> <td>   1.455</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.056</td> <th>  Jarque-Bera (JB):  </th> <td>   5.555</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 0.359</td> <th>  Prob(JB):          </th> <td>  0.0622</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 3.611</td> <th>  Cond. No.          </th> <td>    10.3</td>\n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:            petal_width   R-squared:                       0.927\n",
       "Model:                            OLS   Adj. R-squared:                  0.927\n",
       "Method:                 Least Squares   F-statistic:                     1882.\n",
       "Date:                Sun, 16 Oct 2022   Prob (F-statistic):           4.68e-86\n",
       "Time:                        18:49:58   Log-Likelihood:                 24.796\n",
       "No. Observations:                 150   AIC:                            -45.59\n",
       "Df Residuals:                     148   BIC:                            -39.57\n",
       "Df Model:                           1                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "================================================================================\n",
       "                   coef    std err          t      P>|t|      [0.025      0.975]\n",
       "--------------------------------------------------------------------------------\n",
       "const           -0.3631      0.040     -9.131      0.000      -0.442      -0.285\n",
       "petal_length     0.4158      0.010     43.387      0.000       0.397       0.435\n",
       "==============================================================================\n",
       "Omnibus:                        5.765   Durbin-Watson:                   1.455\n",
       "Prob(Omnibus):                  0.056   Jarque-Bera (JB):                5.555\n",
       "Skew:                           0.359   Prob(JB):                       0.0622\n",
       "Kurtosis:                       3.611   Cond. No.                         10.3\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# NOTE: statsmodels does not automatically add an intercept, so you need to do that manually\n",
    "\n",
    "X = smf.add_constant(X)\n",
    "\n",
    "# Build regression model\n",
    "# NOTE: the X and y variables are reversed in the function compared to sklearn\n",
    "\n",
    "fit_lm = smf.OLS(y, X).fit()\n",
    "\n",
    "# Get an R-like output of the model\n",
    "\n",
    "fit_lm.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c21ca36",
   "metadata": {},
   "source": [
    "### First simulate a distirbution around the Slope Coefficient\n",
    "\n",
    "* The slope coefficient for `petal_length` is 0.4158 with a standard error of 0.01.\n",
    "* We can use this information to build a posterior simulation of the slope coefficient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "86d6ad71",
   "metadata": {},
   "outputs": [],
   "source": [
    "## get summary stats\n",
    "mu_slope = 0.4158\n",
    "se_slope = 0.01\n",
    "\n",
    "## create simulation\n",
    "n_samples = 10000\n",
    "coef_sim = np.random.normal(loc = mu_slope, scale = se_slope, size = n_samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "5de45ea3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQnUlEQVR4nO3dbYxcV33H8e8PBxLKg0iUTWRspw6tVepUDaCtoU1VhaZqXIJwUIlqJKjbRrIqpVFoKxGHN1VVWUqlCoUXjZAFFFc8WFaAxgqIYhkiVClNcCA8OI4bN06TxW5skCggVaE2/76YG5jYs57ZnZ19OP5+pNXcOXPuzH9v1r89Offes6kqJElteclSFyBJWniGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSg0YK9ySvSXJfkieSHE7ym0kuS7I/yZPd46V9/e9KcjTJkSQ3Tq58SdIgo47cPwR8sapeD1wLHAZ2AAeqagNwoHtOko3AVuAaYDNwb5JVC124JGl2Q8M9yauB3wE+ClBVP6mqHwBbgN1dt93Azd32FmBPVT1fVceAo8CmhS1bknQ+F43Q53XAKeCfklwLPArcAVxZVScAqupEkiu6/muAf+/bf6Zrm9Xll19e69evn2PpknRhe/TRR79XVVODXhsl3C8C3gTcXlUPJ/kQ3RTMLDKg7Zw1DpJsB7YDXHXVVRw8eHCEUiRJL0jyX7O9Nsqc+wwwU1UPd8/voxf2zyVZ3X3AauBkX/91ffuvBY6f/aZVtauqpqtqempq4C8eSdI8DQ33qvpv4Nkkv9I13QA8DuwDtnVt24D7u+19wNYkFye5GtgAPLKgVUuSzmuUaRmA24FPJnkZ8BTwp/R+MexNcivwDHALQFUdSrKX3i+A08BtVXVmwSuXJM1qpHCvqseA6QEv3TBL/53AzvmXJUkah3eoSlKDDHdJapDhLkkNMtwlqUGGuyQ1aNRLISUNsX7H589pe/rum5agEsmRuyQ1yXCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalB3sQkzcOgG5ak5cSRuyQ1yHCXpAY5LSP1cX0YtcKRuyQ1yHCXpAYZ7pLUIMNdkhrkCVVpkXnSVovBkbskNchwl6QGGe6S1CDn3KUJcg0aLZWRRu5Jnk7y7SSPJTnYtV2WZH+SJ7vHS/v635XkaJIjSW6cVPGSpMHmMi3z1qp6Q1VNd893AAeqagNwoHtOko3AVuAaYDNwb5JVC1izJGmIcaZltgDXd9u7gQeBO7v2PVX1PHAsyVFgE/DQGJ8lLRmnVrQSjTpyL+BLSR5Nsr1ru7KqTgB0j1d07WuAZ/v2nenaJEmLZNSR+3VVdTzJFcD+JE+cp28GtNU5nXq/JLYDXHXVVSOWIUkaxUgj96o63j2eBD5Hb5rluSSrAbrHk133GWBd3+5rgeMD3nNXVU1X1fTU1NT8vwNJ0jmGhnuSVyR51QvbwO8D3wH2Adu6btuA+7vtfcDWJBcnuRrYADyy0IVLkmY3yrTMlcDnkrzQ/1NV9cUkXwP2JrkVeAa4BaCqDiXZCzwOnAZuq6ozE6lekjTQ0HCvqqeAawe0fx+4YZZ9dgI7x65OWiAu1qULjcsPSFKDDHdJapDhLkkNMtwlqUGuCiktA57w1UJz5C5JDTLcJalBhrskNchwl6QGeUJVFyzXaVfLHLlLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIq2WkZcolCTQOR+6S1CDDXZIaZLhLUoOcc1dTvOtU6nHkLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWrQyOGeZFWSbyR5oHt+WZL9SZ7sHi/t63tXkqNJjiS5cRKFS5JmN5eR+x3A4b7nO4ADVbUBONA9J8lGYCtwDbAZuDfJqoUpV5I0ipHCPcla4CbgI33NW4Dd3fZu4Oa+9j1V9XxVHQOOApsWpFpJ0khGHbnfA7wf+Glf25VVdQKge7yia18DPNvXb6Zre5Ek25McTHLw1KlTc61bknQeQ8M9yduBk1X16IjvmQFtdU5D1a6qmq6q6ampqRHfWpI0ilFWhbwOeEeStwGXAK9O8gnguSSrq+pEktXAya7/DLCub/+1wPGFLFqSdH5DR+5VdVdVra2q9fROlH65qt4D7AO2dd22Afd32/uArUkuTnI1sAF4ZMErlyTNapz13O8G9ia5FXgGuAWgqg4l2Qs8DpwGbquqM2NXKp3Ftdul2c0p3KvqQeDBbvv7wA2z9NsJ7ByzNknSPPmXmKQVZND/rTx9901LUImWO5cfkKQGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ3yUkitCN6wJM2NI3dJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGuSqk1CD/kLYcuUtSgwx3SWqQ0zLSCucfMtEgjtwlqUGGuyQ1aGi4J7kkySNJvpnkUJK/7dovS7I/yZPd46V9+9yV5GiSI0lunOQ3IEk61ygj9+eB362qa4E3AJuTvAXYARyoqg3Age45STYCW4FrgM3AvUlWTaB2SdIshoZ79fy4e/rS7quALcDurn03cHO3vQXYU1XPV9Ux4CiwaSGLliSd30hz7klWJXkMOAnsr6qHgSur6gRA93hF130N8Gzf7jNd29nvuT3JwSQHT506Nca3IEk620jhXlVnquoNwFpgU5JfO0/3DHqLAe+5q6qmq2p6ampqpGIlSaOZ09UyVfUD4EF6c+nPJVkN0D2e7LrNAOv6dlsLHB+3UEnS6Ea5WmYqyWu67ZcDvwc8AewDtnXdtgH3d9v7gK1JLk5yNbABeGSB65Yknccod6iuBnZ3V7y8BNhbVQ8keQjYm+RW4BngFoCqOpRkL/A4cBq4rarOTKZ8SdIgQ8O9qr4FvHFA+/eBG2bZZyewc+zqdEHydnppfN6hKkkNMtwlqUGGuyQ1yCV/pQvEbOcy/AtNbXLkLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDfIOVS0pV4CUJsORuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNcg7VLVovBtVWjyGu3SBG/RL1z+avfI5LSNJDTLcJalBQ8M9ybokX0lyOMmhJHd07Zcl2Z/kye7x0r597kpyNMmRJDdO8huQJJ1rlJH7aeCvq+pXgbcAtyXZCOwADlTVBuBA95zuta3ANcBm4N4kqyZRvCRpsKEnVKvqBHCi2/5RksPAGmALcH3XbTfwIHBn176nqp4HjiU5CmwCHlro4iVNhidZV745zbknWQ+8EXgYuLIL/hd+AVzRdVsDPNu320zXdvZ7bU9yMMnBU6dOzaN0SdJsRg73JK8EPgO8r6p+eL6uA9rqnIaqXVU1XVXTU1NTo5YhSRrBSOGe5KX0gv2TVfXZrvm5JKu711cDJ7v2GWBd3+5rgeMLU64kaRSjXC0T4KPA4ar6YN9L+4Bt3fY24P6+9q1JLk5yNbABeGThSpYkDTPKHarXAe8Fvp3ksa7tA8DdwN4ktwLPALcAVNWhJHuBx+ldaXNbVZ1Z6MIlSbMb5WqZf2PwPDrADbPssxPYOUZdkqQxeIeqJDXIcJekBrkqpCbC5X2lpeXIXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNcuEwSSMZtBjc03fftASVaBSO3CWpQYa7JDXIaRmNzbXbpeXHkbskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkJdCSpo371pdvhy5S1KDhoZ7ko8lOZnkO31tlyXZn+TJ7vHSvtfuSnI0yZEkN06qcEnS7EYZuX8c2HxW2w7gQFVtAA50z0myEdgKXNPtc2+SVQtWrSRpJEPn3Kvqq0nWn9W8Bbi+294NPAjc2bXvqarngWNJjgKbgIcWqF4tMZcakFaG+c65X1lVJwC6xyu69jXAs339Zro2SdIiWugTqhnQVgM7JtuTHExy8NSpUwtchiRd2OYb7s8lWQ3QPZ7s2meAdX391gLHB71BVe2qqumqmp6amppnGZKkQeZ7nfs+YBtwd/d4f1/7p5J8EHgtsAF4ZNwitTScX5dWrqHhnuTT9E6eXp5kBvgbeqG+N8mtwDPALQBVdSjJXuBx4DRwW1WdmVDtkqRZjHK1zLtneemGWfrvBHaOU5QkaTzeoSpJDTLcJalBLhwmaUG5mNjy4MhdkhpkuEtSgwx3SWqQc+4CvGFJao0jd0lqkCN3SRPnFTSLz5G7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkN8jr3C5B3o0rtM9wbZ5BrufLGpslyWkaSGuTIXdKy4Wh+4Thyl6QGOXKXtKw5mp8fR+6S1CBH7o3wqhhdSBzND+fIXZIa5Mh9mXOEImk+DHdJTXAg9GITC/ckm4EPAauAj1TV3ZP6rAuN8+uShpnInHuSVcA/An8AbATenWTjJD5LknSuSY3cNwFHq+opgCR7gC3A4xP6vCY4IpcW1qhTNbP921vJ0zqTCvc1wLN9z2eAN0/os8b+DzjKvuPUImn5mMu/0XHm8Rcjb85nUuGeAW31og7JdmB79/THSY4saAF/P7F9Lwe+N/93v2B4nIbzGI1mWR2ncbJlAu/3i7O9MKlwnwHW9T1fCxzv71BVu4BdE/r8iUlysKqml7qO5c7jNJzHaDQep/mZ1E1MXwM2JLk6ycuArcC+CX2WJOksExm5V9XpJH8B/Cu9SyE/VlWHJvFZkqRzTew696r6AvCFSb3/ElpxU0lLxOM0nMdoNB6neUhVDe8lSVpRXDhMkhpkuHeSbE5yJMnRJDsGvL4lybeSPJbkYJLf7nvtL5McSvKdJJ9OcsniVr94hh2nvn6/keRMknfNdd8WzPc4JVmX5CtJDnc/U3csXtWLb5yfp659VZJvJHlg8tWuMFV1wX/RO+n7n8DrgJcB3wQ2ntXnlfx8GuvXgSe67TXAMeDl3fO9wJ8s9fe0VMepr9+X6Z1zeddc9m3ha8zjtBp4U7f9KuA/PE7nHqe+1/4K+BTwwFJ/P8vty5F7z8+WS6iqnwAvLJfwM1X14+p+moBX8OKbsi4CXp7kIuAXOOua/oYMPU6d24HPACfnsW8L5n2cqupEVX292/4RcJjeAKJF4/w8kWQtcBPwkUkXuhIZ7j2Dlks45x9UkncmeQL4PPBnAFX1XeAfgGeAE8D/VNWXJl7x0hh6nJKsAd4JfHiu+zZknOPU32c98Ebg4YUvcVkY9zjdA7wf+OmE6lvRDPeeocslAFTV56rq9cDNwN8BJLmU3mjjauC1wCuSvGdypS6pUY7TPcCdVXVmHvu2Ypzj1HuD5JX0Rqvvq6ofLmx5y8a8j1OStwMnq+rRCdW24vnHOnqGLpfQr6q+muSXklwOvBU4VlWnAJJ8Fvgt4BMTrHepjHKcpoE9SaC3Jsjbkpwecd9WzPs4VdW/JHkpvWD/ZFV9djEKXiLj/Dy9GXhHkrcBlwCvTvKJqmp1YDV3Sz3pvxy+6P2Se4re6PuFEzvXnNXnl/n5CdU3Ad+lN/J4M3CI3lx7gN3A7Uv9PS3VcTqr/8f5+YnCOe27kr/GPE4B/hm4Z6m/j+V8nM5qvx5PqJ7z5cid2ZdLSPLn3esfBv4Q+OMk/wf8L/BH1fvJejjJfcDXgdPAN2j0jroRj9Oc9l2MuhfbOMcJuA54L/DtJI91bR+o3h3fTRnzOGkI71CVpAZ5QlWSGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUoP8HuPkxaS2uqQ4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "## plot simulated distribution\n",
    "\n",
    "plt.hist(coef_sim, bins = 60)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "dc08cc14",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'Mean': 0.4162082668715149,\n",
       " 'Low90': 0.3669853338376307,\n",
       " 'High90': 0.46503374384504764}"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## get summary stats from our simulation\n",
    "summary_stats = {\n",
    "    'Mean': coef_sim.mean(),\n",
    "    'Low90': np.quantile(coef_sim, 0.05),\n",
    "    'High90': np.quantile(coef_sim, 0.95)\n",
    "}\n",
    "\n",
    "summary_stats"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adcd54cd",
   "metadata": {},
   "source": [
    "### Make a prediction on a new obseration and build a posterior predictive simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "id": "0102f316",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1. , 1.7]])"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "new_dat = np.array([[1, 1.7]])    # add a 1 upfront to indicate the model intercept\n",
    "new_dat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "id": "141ae117",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<statsmodels.regression._prediction.PredictionResults at 0x1234971d0>"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prediction = fit_lm.get_prediction(new_dat)\n",
    "prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "id": "0e729754",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>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",
       "      <th>0</th>\n",
       "      <td>0.343709</td>\n",
       "      <td>0.025945</td>\n",
       "      <td>0.292438</td>\n",
       "      <td>0.394979</td>\n",
       "      <td>-0.067538</td>\n",
       "      <td>0.754956</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       mean   mean_se  mean_ci_lower  mean_ci_upper  obs_ci_lower  \\\n",
       "0  0.343709  0.025945       0.292438       0.394979     -0.067538   \n",
       "\n",
       "   obs_ci_upper  \n",
       "0      0.754956  "
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prediction.summary_frame()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da723438",
   "metadata": {},
   "source": [
    "* We want to create a posterior predictive distribution for this individual new observation (not a population average). Thus, we will use the prediction interval, which is called `obs_ci_lower` and `obs_ci_upper` here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "id": "f00f5d76",
   "metadata": {},
   "outputs": [],
   "source": [
    "mu_pred = 0.343709\n",
    "se_pred = 0.754956 - 0.343709     # subtract the upper prediction interval from the mean to get the variability\n",
    "n_sims = 10000\n",
    "\n",
    "pred_obs = np.random.normal(loc = mu_pred, scale = se_pred, size = n_sims)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "id": "228b3b96",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAP0klEQVR4nO3dXahc13nG8f8TNXVK4hCrOlIUSY5cUEvsQj4QqluX4tZpbOxSuRcGFdrqwiACDiRQaOQWWnohUHoR2kJ9IZJQlSYxgiS1SJMmihoTCokdObVjy4oqJVbtg4SkODQfN2rtvr0422UizWhG58ycmVn6/+Cw96xZM/Pu2TqP1lmz955UFZKktrxu2gVIksbPcJekBhnuktQgw12SGmS4S1KDfmbaBQCsW7eutm7dOu0yJGmuPPXUU9+vqoV+981EuG/dupVjx45NuwxJmitJ/nPQfU7LSFKDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSg2biDFVpXLbu/ee+7Wf237fKlUjT5chdkhpkuEtSgwx3SWqQ4S5JDTLcJalBI4V7kjNJnk3ydJJjXdvaJEeSnOqWN/X0fzjJ6SQnk9w9qeIlSf1dy8j9N6vqXVW1vbu9FzhaVduAo91tktwK7AJuA+4BHkmyZow1S5KGWMlx7juBO7v1g8DjwIe79ker6hLwQpLTwA7g6yt4LWlF+h3/7rHvatmoI/cCvpzkqSR7urYNVXUOoFuu79o3AS/1PHaxa5MkrZJRR+53VNXZJOuBI0m+c5W+6dNWV3Ra+k9iD8DNN988YhmSpFGMNHKvqrPd8gLwOZamWc4n2QjQLS903ReBLT0P3wyc7fOcB6pqe1VtX1jo++XdkqRlGhruSd6Y5MbX1oH3Ac8Bh4HdXbfdwGPd+mFgV5IbktwCbAOeHHfhkqTBRpmW2QB8Lslr/T9VVf+S5JvAoSQPAi8CDwBU1fEkh4DngVeAh6rq1YlUL0nqa2i4V9X3gHf2aX8ZuGvAY/YB+1ZcnbTKPKpGrfAMVUlqkOEuSQ3yyzp03Rr0xR5SCxy5S1KDHLlrbjnylgYz3DVzPGJFWjmnZSSpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDPBRSc8Fj2qVr48hdkhrkyF0awpOqNI8cuUtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGeSikNCYeMqlZYrhrqjzzVJoMp2UkqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDRo53JOsSfLvST7f3V6b5EiSU93ypp6+Dyc5neRkkrsnUbgkabBrOUP1g8AJ4M3d7b3A0aran2Rvd/vDSW4FdgG3AW8DvpLkF6vq1THWLU2VZ9Zq1o00ck+yGbgP+FhP807gYLd+ELi/p/3RqrpUVS8Ap4EdY6lWkjSSUUfufw38CXBjT9uGqjoHUFXnkqzv2jcB3+jpt9i1/ZQke4A9ADfffPO1Va255GhXWj1DR+5Jfge4UFVPjfic6dNWVzRUHaiq7VW1fWFhYcSnliSNYpSR+x3A7ya5F3gD8OYk/wicT7KxG7VvBC50/ReBLT2P3wycHWfRkqSrGzpyr6qHq2pzVW1l6YPSf62qPwAOA7u7bruBx7r1w8CuJDckuQXYBjw59solSQOt5Hru+4FDSR4EXgQeAKiq40kOAc8DrwAPeaSMJK2uawr3qnoceLxbfxm4a0C/fcC+FdYmSVomz1CVpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNWglZ6hKGqLflTDP7L9vCpXoemO4S6vMwNdqcFpGkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ3ywmHSDPBiYho3R+6S1CDDXZIaZLhLUoOcc9dE9JtDlrR6HLlLUoMMd0lqkOEuSQ0aGu5J3pDkySTPJDme5C+79rVJjiQ51S1v6nnMw0lOJzmZ5O5JboAk6UqjjNwvAb9VVe8E3gXck+R2YC9wtKq2AUe72yS5FdgF3AbcAzySZM0EapckDTA03GvJT7qbr+9+CtgJHOzaDwL3d+s7gUer6lJVvQCcBnaMs2hJ0tWNNOeeZE2Sp4ELwJGqegLYUFXnALrl+q77JuClnocvdm2XP+eeJMeSHLt48eIKNkGSdLmRwr2qXq2qdwGbgR1Jfvkq3dPvKfo854Gq2l5V2xcWFkYqVpI0mms6Wqaq/gt4nKW59PNJNgJ0ywtdt0VgS8/DNgNnV1qoJGl0oxwts5DkLd36zwHvBb4DHAZ2d912A49164eBXUluSHILsA14csx1S5KuYpTLD2wEDnZHvLwOOFRVn0/ydeBQkgeBF4EHAKrqeJJDwPPAK8BDVfXqZMrXavOyAtJ8GBruVfVt4N192l8G7hrwmH3AvhVXJ0laFs9QlaQGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIL8iWZlS/s4HP7L9vCpVoHjlyl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIL9mTwP1+5o3SfPBkbskNchwl6QGOS0jwCkYqTVDR+5JtiT5apITSY4n+WDXvjbJkSSnuuVNPY95OMnpJCeT3D3JDZAkXWmUaZlXgD+uqncAtwMPJbkV2AscraptwNHuNt19u4DbgHuAR5KsmUTxkqT+hoZ7VZ2rqm916z8GTgCbgJ3Awa7bQeD+bn0n8GhVXaqqF4DTwI4x1y1Juopr+kA1yVbg3cATwIaqOgdL/wEA67tum4CXeh622LVd/lx7khxLcuzixYvLKF2SNMjI4Z7kTcBngA9V1Y+u1rVPW13RUHWgqrZX1faFhYVRy5AkjWCko2WSvJ6lYP9kVX22az6fZGNVnUuyEbjQtS8CW3oevhk4O66CpetZv6Oazuy/bwqVaNaNcrRMgI8DJ6rqoz13HQZ2d+u7gcd62ncluSHJLcA24MnxlSxJGmaUkfsdwB8CzyZ5umv7U2A/cCjJg8CLwAMAVXU8ySHgeZaOtHmoql4dd+GSpMGGhntV/Rv959EB7hrwmH3AvhXUJUlaAS8/IEkNMtwlqUGGuyQ1yHCXpAYZ7pLUIC/5K805T2xSP47cJalBjtyvQ34xh9Q+R+6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgz1CVGuT1ZmS4N85LDUjXJ6dlJKlBhrskNchpGek6MWiKzrn4Njlyl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ3yUEjpOuelCtrkyF2SGjQ03JN8IsmFJM/1tK1NciTJqW55U899Dyc5neRkkrsnVbgkabBRRu5/D9xzWdte4GhVbQOOdrdJciuwC7ite8wjSdaMrVpJ0kiGhntVfQ34wWXNO4GD3fpB4P6e9ker6lJVvQCcBnaMp1RJ0qiWO+e+oarOAXTL9V37JuClnn6LXdsVkuxJcizJsYsXLy6zDElSP+P+QDV92qpfx6o6UFXbq2r7wsLCmMuQpOvbcsP9fJKNAN3yQte+CGzp6bcZOLv88iRJy7Hc49wPA7uB/d3ysZ72TyX5KPA2YBvw5EqL1HB+45KkXkPDPcmngTuBdUkWgb9gKdQPJXkQeBF4AKCqjic5BDwPvAI8VFWvTqh2SdIAQ8O9qn5/wF13Dei/D9i3kqIkSSvj5QfmkFMwkobx8gOS1CDDXZIa5LSMpCt4pcj558hdkhrkyH3G+eGppOVw5C5JDXLkLmkkzsPPF0fuktQgw12SGmS4S1KDDHdJapAfqEpaNj9knV2O3CWpQYa7JDXIcJekBjnnPkO81ICkcXHkLkkNcuQuaaw8gmY2OHKXpAYZ7pLUIKdlxmzUP0n98FTSJDlyl6QGGe6S1CCnZVaBUzCSVpsjd0lqkOEuSQ1yWmZETq1Iyzfq748nO42PI3dJapAjd0kzwxH++Dhyl6QGTWzknuQe4G+ANcDHqmr/pF7LCxVJ15eVfAZ2vWTDRMI9yRrg74DfBhaBbyY5XFXPT+L1+lnJn3d+eCpp3k1q5L4DOF1V3wNI8iiwE1i1cB+VQS5dX1brd37UgeOk/pKYVLhvAl7qub0I/EpvhyR7gD3dzZ8kOTmBOtYB35/A864mt2E2zPs2zHv9MGfbkI/0bb5iGwb0G9XbB90xqXBPn7b6qRtVB4ADE3r9pSKSY1W1fZKvMWluw2yY922Y9/rBbbhWkzpaZhHY0nN7M3B2Qq8lSbrMpML9m8C2JLck+VlgF3B4Qq8lSbrMRKZlquqVJB8AvsTSoZCfqKrjk3itISY67bNK3IbZMO/bMO/1g9twTVJVw3tJkuaKZ6hKUoMMd0lqUFPhnuSBJMeT/G+SgYcbJbknyckkp5PsXc0ah0myNsmRJKe65U0D+p1J8mySp5McW+06+9Rz1fc0S/62u//bSd4zjTqvZoRtuDPJD7v3/Okkfz6NOgdJ8okkF5I8N+D+edgHw7Zh1vfBliRfTXKiy6IP9umzOvuhqpr5Ad4B/BLwOLB9QJ81wHeBXwB+FngGuHXatffU91fA3m59L/CRAf3OAOumXe+o7ylwL/BFls6BuB14Ytp1L2Mb7gQ+P+1ar7INvwG8B3huwP0zvQ9G3IZZ3wcbgfd06zcC/zGt34WmRu5VdaKqhp3p+v+XRqiq/wZeuzTCrNgJHOzWDwL3T6+UkY3ynu4E/qGWfAN4S5KNq13oVcz6v4uhquprwA+u0mXW98Eo2zDTqupcVX2rW/8xcIKlM/Z7rcp+aCrcR9Tv0giXv/nTtKGqzsHSPxRg/YB+BXw5yVPdpRymaZT3dNbf91Hr+9UkzyT5YpLbVqe0sZn1fTCqudgHSbYC7waeuOyuVdkPc/dlHUm+Ary1z11/VlWPjfIUfdpW9XjQq23DNTzNHVV1Nsl64EiS73SjnmkY5T2d+vs+xCj1fQt4e1X9JMm9wD8B2yZd2BjN+j4YxVzsgyRvAj4DfKiqfnT53X0eMvb9MHfhXlXvXeFTTP3SCFfbhiTnk2ysqnPdn2oXBjzH2W55IcnnWJpWmFa4j/KeTv19H2Jofb2/pFX1hSSPJFlXVfNyMatZ3wdDzcM+SPJ6loL9k1X12T5dVmU/XI/TMrN+aYTDwO5ufTdwxV8jSd6Y5MbX1oH3AX2PLlglo7ynh4E/6o4UuB344WvTTzNi6DYkeWuSdOs7WPr9eXnVK12+Wd8HQ836Puhq+zhwoqo+OqDb6uyHaX+6PM4f4PdY+l/xEnAe+FLX/jbgCz397mXpU+zvsjSdM/Xae2r7eeAocKpbrr18G1g6ouOZ7uf4LGxDv/cUeD/w/m49LH2By3eBZxlwNNOMb8MHuvf7GeAbwK9Nu+bL6v80cA74n+734ME53AfDtmHW98GvszTF8m3g6e7n3mnsBy8/IEkNuh6nZSSpeYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJatD/AX2Pl9UqoOtNAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(pred_obs, bins = 60)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "id": "284be89b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'Mean': 0.35544530885021625,\n",
       " 'Low90': -0.3178877847515099,\n",
       " 'High90': 1.044449387252352}"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## get summary stats from our simulation\n",
    "summary_stats = {\n",
    "    'Mean': pred_obs.mean(),\n",
    "    'Low90': np.quantile(pred_obs, 0.05),\n",
    "    'High90': np.quantile(pred_obs, 0.95)\n",
    "}\n",
    "\n",
    "summary_stats"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df871173",
   "metadata": {},
   "source": [
    "### How was the prediction interval calculated?\n",
    "\n",
    "* It may help to see how a prediciton interval (and confidence interval) are calculated by hand, so that we can have an apprecication for how they work and how we can use them in other scenarios.\n",
    "\n",
    "**First, to calculate either the the confidence interval or prediction interval there are a few pieces of information that we need:**\n",
    "\n",
    "1) The model degrees of freedom, which is called `DF Residual` in our model output above.\n",
    "\n",
    "\n",
    "2) The number of observations in our data.\n",
    "\n",
    "\n",
    "3) The t-critical value that represents how certain we want to be (here we will build a 95% Interval).\n",
    "\n",
    "\n",
    "4) The average and standard deviation of the independent variable (`petal_length`) in our data.\n",
    "\n",
    "5) The residual standard error for the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "id": "46b1d117",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_df = 148    # this makes sense given that we have 2 parameters in our model (slope & intercept) so degrees of freedom is 150 - 2\n",
    "n_obs = len(iris_df)\n",
    "petal_length_mean = iris_df['petal_length'].mean()\n",
    "petal_length_sd = iris_df['petal_length'].std()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ed16308",
   "metadata": {},
   "source": [
    "Next we need the `t-critical value`. For this you need:\n",
    "\n",
    "1) Model degrees of freedom (which we have above)\n",
    "\n",
    "2) The alpha level corresponding to the percent of confidence you want for your interval\n",
    "\n",
    "We are interested in 95% confidence, meaning that the interval should encompass 95% of the distribution. Thus, we are left with 5% (1 - 0.95) of the data in the tails. Because this is a two sided test we end up shaving 2.5% (0.05 / 2) off both the right and left sides of the distribution. Thus, we end up with `0.025` in the left tail and `0.975` in the right tail. Let's use that upper tail percentage to get our t-critical value.\n",
    "\n",
    "**NOTE: General rule of thumb, if you have a good amount of data is that the critical value for 90% is 1.65, 95% is 1.96 (often times people round this up to 2, to be safe), and for 99% is 2.58** \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "id": "dfa0f10b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.9761224936033632"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t_crit = stats.t(df = model_df).ppf(0.975)\n",
    "t_crit"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9d047cd",
   "metadata": {},
   "source": [
    "Now that we have the t-critical value, the final piece is the model's residual standard error. This is simply calculated as:\n",
    "\n",
    "$rse = sqrt((sum.of.squared.residuals)^2 / model.df)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "id": "5c887c4a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.20648434891360867"
      ]
     },
     "execution_count": 110,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "predicted_values = fit_lm.predict()\n",
    "squared_resid = (y - predicted_values)**2\n",
    "sum_squared_resid = sum_squared_resid.sum()\n",
    "\n",
    "rse = np.sqrt(sum_squared_resid / model_df)\n",
    "rse"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0c8c222",
   "metadata": {},
   "source": [
    "**Prediction Interval Calculation**\n",
    "\n",
    "Now that we have the ingredients, we calculate the prediction interval:\n",
    "\n",
    "$pred.interval95 = t.crit * rse * sqrt(1 + 1 / N.obs + ((petal.length.obs - petal.length.avg)^2) / ((N - 1) * petal.length.sd^2))$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "id": "77380f33",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.4112468496052437"
      ]
     },
     "execution_count": 130,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pred95 = t_crit * rse * np.sqrt(1 + 1/n_obs + ((new_dat[0][1] - petal_length_mean)**2) / ((n_obs - 1) * petal_length_sd**2))\n",
    "pred95"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "id": "ea2469f0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.411247"
      ]
     },
     "execution_count": 132,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## Confirm that we got the same result as the function above from statsmodels.api gave us.\n",
    "se_pred"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc085cc4",
   "metadata": {},
   "source": [
    "**It worked!**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9f238eb",
   "metadata": {},
   "source": [
    "### How is a Confidence Interval Created?\n",
    "\n",
    "The calculation for our confidence interval is a similar to that of the prediction interval but we are missing one of the terms (we are estimating at a population level here, so we don't have to bake in additional uncertainty). The good news is that we already have the ingredients. The confidence level can be calculated as such:\n",
    "\n",
    "$CL95 = t.crit * rse * sqrt(1/N.obs + ((petal.length.obs - petal.length.avg)^2) / ((N.obs - 1) * petal.length.sd^2))$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "id": "5540f54c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.05127048666266827"
      ]
     },
     "execution_count": 133,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cl95 = t_crit * rse * np.sqrt(1/n_obs + ((new_dat[0][1] - petal_length_mean)**2) / ((n_obs - 1) * petal_length_sd**2))\n",
    "cl95"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "id": "ac4d644b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.05127000000000004"
      ]
     },
     "execution_count": 134,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## Subtract the mean prediction for the new observation from the upper confidence interval, from the above\n",
    "# table to get the Confidence Level from statsmodels.api\n",
    "\n",
    "0.394979 - 0.343709"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dabefd3c",
   "metadata": {},
   "source": [
    "**Our math checks out and we get the same result! Notice that the Confidence Interval is smaller than the Prediciton Interval because, again, we are making population inference with the Confidence Interval and thus we don't have as much uncertainty as we would if we were estimating for a single observation.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e28e2cf",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}