Lecture05_04 - Understanding AIC.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Understanding AIC\n",
    "\n",
    "##### Keywords: AIC, regularization, deviance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook is based on McElreath, Rethinking Statistics, Chapter 6."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we use the empirical distribution and sample quantities here we are working with our training sample (s).\n",
    "\n",
    "Clearly we can calculate deviance on the validation and test samples as well to remedy this issue. And the results will be similar to what we found in lecture for MSE, with the training deviance decreasing with complexity and the testing deviance increasing at some point. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A trick to generate data\n",
    "\n",
    "We generate data from a gaussian with standard deviation 1 and means given by:\n",
    "\n",
    "$$\\mu_i = 0.15 x_{1,i} - 0.4 x_{2,i}, y \\sim N(\\mu, 1).$$\n",
    "\n",
    "This is a **2 parameter** model.\n",
    "\n",
    "We use an interesting trick to generate this data, directly using the regression coefficients as correlations with the response variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def generate_data(N, k, rho=[0.15, -0.4]):\n",
    "    n_dim = 1 + len(rho)\n",
    "    if n_dim < k:\n",
    "        n_dim = k\n",
    "    Rho = np.eye(n_dim)\n",
    "    for i,r in enumerate(rho):\n",
    "        Rho[0, i+1] = r\n",
    "    index_lower = np.tril_indices(n_dim, -1)\n",
    "    Rho[index_lower] = Rho.T[index_lower]\n",
    "    mean = n_dim * [0.]\n",
    "    Xtrain = np.random.multivariate_normal(mean, Rho, size=N)\n",
    "    Xtest = np.random.multivariate_normal(mean, Rho, size=N)\n",
    "    ytrain = Xtrain[:,0].copy()\n",
    "    Xtrain[:,0]=1.\n",
    "    ytest = Xtest[:,0].copy()\n",
    "    Xtest[:,0]=1.\n",
    "    return Xtrain[:,:k], ytrain, Xtest[:,:k], ytest"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to generate data for 5 different cases, a one parameter (intercept) fit, a two parameter (intercept and $x_1$), three parameters (add a $x_2), and four and five parameters. Here is what the data looks like for 2 parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[ 1.        , -0.83978695],\n",
       "        [ 1.        , -0.60882982],\n",
       "        [ 1.        ,  1.02567296],\n",
       "        [ 1.        ,  0.24801809],\n",
       "        [ 1.        , -1.08181661],\n",
       "        [ 1.        , -1.85677575],\n",
       "        [ 1.        ,  1.82835523],\n",
       "        [ 1.        ,  0.35622585],\n",
       "        [ 1.        , -0.04159412],\n",
       "        [ 1.        ,  0.58678675],\n",
       "        [ 1.        , -0.24396323],\n",
       "        [ 1.        , -0.07081137],\n",
       "        [ 1.        ,  0.46510137],\n",
       "        [ 1.        , -1.02993129],\n",
       "        [ 1.        , -2.08756332],\n",
       "        [ 1.        ,  0.60666556],\n",
       "        [ 1.        ,  0.45913243],\n",
       "        [ 1.        ,  0.60083017],\n",
       "        [ 1.        , -1.05726496],\n",
       "        [ 1.        , -0.52258973]]),\n",
       " array([-1.11513393, -0.50856507,  0.50782261, -0.09031626,  0.41992084,\n",
       "        -0.82404287,  0.27567933,  0.3626567 ,  0.99109211,  1.14742966,\n",
       "         0.53597334, -1.2959274 ,  2.12659247,  0.09595858,  0.05845798,\n",
       "         0.47581813, -1.02115871,  0.83942264,  0.33097791, -1.07482199]),\n",
       " array([[ 1.        , -0.09664154],\n",
       "        [ 1.        ,  1.68464504],\n",
       "        [ 1.        , -1.63102144],\n",
       "        [ 1.        , -0.83585358],\n",
       "        [ 1.        , -1.0022563 ],\n",
       "        [ 1.        , -0.40901251],\n",
       "        [ 1.        ,  0.52024856],\n",
       "        [ 1.        ,  0.64056776],\n",
       "        [ 1.        , -0.26979402],\n",
       "        [ 1.        ,  0.57670424],\n",
       "        [ 1.        , -0.13580787],\n",
       "        [ 1.        , -0.74665431],\n",
       "        [ 1.        , -0.34801499],\n",
       "        [ 1.        , -0.53583385],\n",
       "        [ 1.        ,  1.49971783],\n",
       "        [ 1.        ,  0.47265248],\n",
       "        [ 1.        , -0.40158879],\n",
       "        [ 1.        ,  0.59618203],\n",
       "        [ 1.        ,  0.63314497],\n",
       "        [ 1.        , -0.85947691]]),\n",
       " array([-1.2830826 , -0.0773539 , -1.22779576, -1.43955577,  0.39940223,\n",
       "        -1.52159959, -0.41312511, -0.95244826,  0.20244849, -0.32376445,\n",
       "         0.09636247,  0.03435469,  0.29289397,  0.70749769, -2.20920373,\n",
       "         0.36671712, -0.73570139, -0.381103  , -0.3126861 , -0.61196652]))"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "generate_data(20,2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And for four parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[  1.00000000e+00,  -5.64117484e-01,  -1.30408291e+00,\n",
       "          -4.06307198e-01],\n",
       "        [  1.00000000e+00,   2.45856192e-01,  -1.13160363e+00,\n",
       "           6.99099707e-01],\n",
       "        [  1.00000000e+00,  -5.92401483e-01,  -5.51929080e-01,\n",
       "           1.70288811e-01],\n",
       "        [  1.00000000e+00,   1.40350006e+00,  -7.42482462e-01,\n",
       "           6.90299071e-01],\n",
       "        [  1.00000000e+00,  -1.14026512e+00,   2.27882734e-01,\n",
       "          -2.80250494e-01],\n",
       "        [  1.00000000e+00,  -1.79114172e-01,   1.71257237e+00,\n",
       "           1.32182974e+00],\n",
       "        [  1.00000000e+00,   8.39677171e-01,  -2.07787502e-01,\n",
       "           1.20281542e+00],\n",
       "        [  1.00000000e+00,  -9.38668901e-01,  -5.87192846e-01,\n",
       "           9.91223102e-01],\n",
       "        [  1.00000000e+00,  -4.11883974e-01,  -1.31283133e+00,\n",
       "          -9.42131126e-01],\n",
       "        [  1.00000000e+00,   5.27622295e-01,   2.98370087e-01,\n",
       "          -3.13398528e-01],\n",
       "        [  1.00000000e+00,   1.75945182e+00,  -9.55446150e-01,\n",
       "          -5.65605486e-01],\n",
       "        [  1.00000000e+00,   3.66192473e-01,   1.39659262e+00,\n",
       "           5.25449367e-01],\n",
       "        [  1.00000000e+00,  -8.27065820e-01,   2.07687904e-01,\n",
       "          -4.07828527e-01],\n",
       "        [  1.00000000e+00,  -9.94963330e-01,  -3.45817822e-01,\n",
       "          -1.71339667e-01],\n",
       "        [  1.00000000e+00,   1.35510550e+00,   2.80027702e-01,\n",
       "           9.10748378e-03],\n",
       "        [  1.00000000e+00,  -5.79066026e-01,   1.67201032e+00,\n",
       "          -1.21304440e+00],\n",
       "        [  1.00000000e+00,   6.61699882e-01,   5.01830803e-01,\n",
       "           3.88336944e-01],\n",
       "        [  1.00000000e+00,   1.33639603e-01,   1.61570562e-01,\n",
       "          -1.96165163e-04],\n",
       "        [  1.00000000e+00,  -1.15492222e+00,  -1.19608702e+00,\n",
       "           7.44868766e-01],\n",
       "        [  1.00000000e+00,  -1.16104055e-03,   5.60105046e-01,\n",
       "           1.13087330e+00]]),\n",
       " array([ 0.59049849, -0.32849878,  0.17111991, -0.18214341, -1.06834661,\n",
       "        -2.01806297,  0.64395116,  1.36524519, -0.39568506,  0.50772571,\n",
       "         2.53131129, -1.54337658, -0.23082485,  0.23672394, -1.72834828,\n",
       "         0.03969115,  0.84937923, -0.04779334, -0.12796287,  0.25091162]),\n",
       " array([[ 1.        ,  0.87805634, -1.3252486 ,  0.9147951 ],\n",
       "        [ 1.        ,  0.04573015,  0.26129224, -0.39094182],\n",
       "        [ 1.        , -1.91453329,  0.26466733, -0.31760853],\n",
       "        [ 1.        ,  0.50704342, -0.91631862,  0.71741282],\n",
       "        [ 1.        , -1.12324824, -1.01977223, -0.52418926],\n",
       "        [ 1.        ,  0.33666476,  0.27513745,  0.63589532],\n",
       "        [ 1.        , -0.38934175,  1.00961086, -0.61853231],\n",
       "        [ 1.        , -0.79629895,  1.28994305, -1.54766776],\n",
       "        [ 1.        , -0.71721218,  0.15741145,  0.475622  ],\n",
       "        [ 1.        ,  1.27433083,  0.6941836 ,  0.90788968],\n",
       "        [ 1.        , -0.68170154, -0.04929347,  0.27673865],\n",
       "        [ 1.        , -0.08147767,  0.81537351, -0.55685094],\n",
       "        [ 1.        , -0.32449028, -0.25664319, -0.57435918],\n",
       "        [ 1.        ,  0.64398473,  0.08495211,  0.47839262],\n",
       "        [ 1.        , -0.91519445, -2.06795016, -0.42275664],\n",
       "        [ 1.        , -0.29770443,  1.08641811,  1.94490458],\n",
       "        [ 1.        , -0.19844572, -2.70288281,  0.69816899],\n",
       "        [ 1.        ,  0.2687054 ,  1.53633517, -0.84276465],\n",
       "        [ 1.        , -0.06377957, -0.05227578,  1.96973628],\n",
       "        [ 1.        , -2.04828043,  0.47438443, -2.28134046]]),\n",
       " array([ 0.02585828, -0.1284505 , -0.39679531, -0.22793629,  0.25824818,\n",
       "        -0.35346682, -1.30972978, -1.2250625 ,  0.97636939, -1.50891804,\n",
       "         0.84376139, -2.22258544,  0.47901623, -1.08042407,  2.35481165,\n",
       "         0.38210287,  1.27868801, -1.71473971,  0.6498696 , -0.27036068]))"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "generate_data(20,4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "//anaconda/envs/py3l/lib/python3.6/site-packages/statsmodels/compat/pandas.py:56: FutureWarning: The pandas.core.datetools module is deprecated and will be removed in a future version. Please use the pandas.tseries module instead.\n",
      "  from pandas.core import datetools\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import norm\n",
    "import statsmodels.api as sm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis, n=20\n",
    "\n",
    "Here is the main loop of our analysis. We take the 5 models we talked about. For each model we generate 10000 samples of the data, split into an equal sized (N=20 each) training and testing set. We fit the regression on the training set, and calculate the deviance on the training set. Notice how we have simply used the `logpdf` from `scipy.stats`. You can easily do this for other distributions.\n",
    "\n",
    "We then use the fit to calculate the $\\mu$ on the test set, and calculate the deviance there. We then find the average and the standard deviation across the 10000 simulations.\n",
    "\n",
    "Why do we do 10000 simulations? These are our **multiple samples from some hypothetical population**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "reps=10000\n",
    "results_20 = {}\n",
    "for k in range(1,6):\n",
    "    trdevs=np.zeros(reps)\n",
    "    tedevs=np.zeros(reps)\n",
    "    for r in range(reps):\n",
    "        Xtr, ytr, Xte, yte = generate_data(20, k)\n",
    "        ols = sm.OLS(ytr, Xtr).fit()\n",
    "        mutr = np.dot(Xtr, ols.params)\n",
    "        devtr = -2*np.sum(norm.logpdf(ytr, mutr, 1))\n",
    "        mute = np.dot(Xte, ols.params)\n",
    "        #print(mutr.shape, mute.shape)\n",
    "        devte = -2*np.sum(norm.logpdf(yte, mute, 1))\n",
    "        #print(k, r, devtr, devte)\n",
    "        trdevs[r] = devtr\n",
    "        tedevs[r] = devte\n",
    "    results_20[k] = (np.mean(trdevs), np.std(trdevs), np.mean(tedevs), np.std(tedevs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>train</th>\n",
       "      <th>train_std</th>\n",
       "      <th>test</th>\n",
       "      <th>test_std</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>55.669331</td>\n",
       "      <td>6.185150</td>\n",
       "      <td>57.688110</td>\n",
       "      <td>6.825813</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>54.301259</td>\n",
       "      <td>5.930990</td>\n",
       "      <td>58.513912</td>\n",
       "      <td>7.284745</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>50.669082</td>\n",
       "      <td>4.778731</td>\n",
       "      <td>56.111054</td>\n",
       "      <td>6.668364</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>49.744479</td>\n",
       "      <td>4.585390</td>\n",
       "      <td>57.414846</td>\n",
       "      <td>7.505377</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>49.026424</td>\n",
       "      <td>4.431733</td>\n",
       "      <td>58.718321</td>\n",
       "      <td>8.279063</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       train  train_std       test  test_std\n",
       "1  55.669331   6.185150  57.688110  6.825813\n",
       "2  54.301259   5.930990  58.513912  7.284745\n",
       "3  50.669082   4.778731  56.111054  6.668364\n",
       "4  49.744479   4.585390  57.414846  7.505377\n",
       "5  49.026424   4.431733  58.718321  8.279063"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "df = pd.DataFrame(results_20).T\n",
    "df = df.rename(columns = dict(zip(range(4), ['train', 'train_std', 'test', 'test_std'])))\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(0.12156862745098039, 0.4666666666666667, 0.7058823529411765),\n",
       " (1.0, 0.4980392156862745, 0.054901960784313725),\n",
       " (0.17254901960784313, 0.6274509803921569, 0.17254901960784313),\n",
       " (0.8392156862745098, 0.15294117647058825, 0.1568627450980392),\n",
       " (0.5803921568627451, 0.403921568627451, 0.7411764705882353),\n",
       " (0.5490196078431373, 0.33725490196078434, 0.29411764705882354),\n",
       " (0.8901960784313725, 0.4666666666666667, 0.7607843137254902),\n",
       " (0.4980392156862745, 0.4980392156862745, 0.4980392156862745),\n",
       " (0.7372549019607844, 0.7411764705882353, 0.13333333333333333),\n",
       " (0.09019607843137255, 0.7450980392156863, 0.8117647058823529)]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import seaborn.apionly as sns\n",
    "colors = sns.color_palette()\n",
    "colors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We plot the traing and testing deviances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "//anaconda/envs/py3l/lib/python3.6/site-packages/matplotlib/axes/_axes.py:2818: MatplotlibDeprecationWarning: Use of None object as fmt keyword argument to suppress plotting of data values is deprecated since 1.4; use the string \"none\" instead.\n",
      "  warnings.warn(msg, mplDeprecation, stacklevel=1)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFx5JREFUeJzt3X2wHXWd5/H3ZwJKDA8RiSyCmtGajbtOITgRhwFxlVFm\nfBiD44K1ygg7I7giPrDiyNQs4jizyy4g49TWIAFUVJaFQggqCoM8rKuFYkIwBBF1ISgRTVATAaOB\nm+/+cfrqTbzce26Svufc2+9X1anT/Tvdp7/nJHU/p7t//etUFZKk7vqdQRcgSRosg0CSOs4gkKSO\nMwgkqeMMAknqOINAkjrOIJCkjjMI1GlJ1iRZl2TemLa/SnLLFN7j6UkuS/LDJBuTfDXJi7dZ5j8k\nuT/Jo0mWJdl7J34MaYcYBBLMAd61A+vvDnwD+ANgb+AS4NokuwMkeT5wAXAcsC/wC+Cfd6RgaWcy\nCCQ4G3hvkvnbs3JV3VtVH66qB6tqpKqWAk8CFjWLvAn4XFV9uaoeAf4L8Poke+yU6qUdZBBIsBy4\nBXjvti8kWZVkwxM8xv1Vn+QgekHwvabp+cA3R1+vqv8H/Ar41zv5c0jbZZdBFyANiTOAryb5yNjG\nqjpwKm+SZE/gU8AHq2pj07w7sHGbRX8OuEegoeAegQRU1Wrg88D7t/c9kswFPgd8rar+25iXHgH2\n3GbxvYCHt3db0s5kEEi/8QHgrcD+ow1J7kryyBM8PjpmuScDy4AHgJO2ed+7gBeMWfa59A4dfafF\nzyL1zUNDUqOqvpfkcuCdwJ1N2/MnWy/JrsCVwCbgLVW1ZZtFLgVuTfIS4HbgQ8BVVeUegYaCewTS\n1v4OmDfpUlv7I+A1wCuBDWP2GF4CUFV3AW+jFwjrmvd/+84rWdox8cY0ktRt7hFIUscZBJLUcQaB\nJHWcQSBJHTcjuo/us88+tXDhwkGXIUkzyooVKx6qqgWTLTcjgmDhwoUsX7580GVI0oyS5P5+lvPQ\nkCR1nEEgSR1nEEhSxxkEktRxBoEkdZxBIEkdZxBIUscZBJLUcQaBJA2rj7+692iZQSBJHWcQSFLH\nGQSS1HEGgSR1nEEgSR1nEEhSxxkEktRxBoEkdZxBIEkdZxBIUscZBNo+03Tpu6T2GQSS1HEGgSR1\nnEEgSR3XahAkmZ/kyiTfTnJ3kkOTnJlkbZI7mser2qxBkjSxXVp+/48A11XVG5I8CXgKcBRwXlWd\n0/K2JUl9aC0IkuwFHAEcD1BVm4HNSdrapCRpO7R5aOh3gfXAx5OsTHJRknnNa6ckWZXkY0meOt7K\nSU5MsjzJ8vXr17dYpiR1W5tBsAvwQuD8qjoYeBR4P3A+8BzgIOBB4NzxVq6qpVW1uKoWL1iwoMUy\nJQ2E16IMjTaD4AHggar6ejN/JfDCqvpxVY1U1RbgQuCQFmuQJE2itSCoqh8BP0iyqGk6EvhWkv3G\nLHY0sLqtGiRJk2u719ApwKVNj6F7gROAf0pyEFDAGuCklmvYPqO7rCdcO9g6JKllrQZBVd0BLN6m\n+bg2tylJmhqvLJakjjMIJKnjDAJJ6jiDQJI6ziCQpI4zCCSp4wwCSeo4g0CSOs4gkKSOMwgkqeMM\nAknqOINAkjrOIJCkjjMIpLZ4By7NEAaBJHWcQSBJHWcQSFLHGQSS1HEGgSR1nEEgSR1nEEhSxxkE\nktRxszoIjr3gVo694NZBlyFJQ21WB4EkaXIGgSR1nEEgSR1nEEhSxxkEktRxBoEkdZxBIEkdZxBI\nUscZBJLUcQaBJHWcQSBJHWcQSFLHtRoESeYnuTLJt5PcneTQJHsnuSHJd5vnp7ZZgyRpYm3vEXwE\nuK6qnge8ALgbeD9wY1X9HnBjM68BcYRWSa0FQZK9gCOAiwGqanNVbQBeB1zSLHYJsKStGiRJk2tz\nj+B3gfXAx5OsTHJRknnAvlX1YLPMj4B9x1s5yYlJlidZvn79+hbLlKRuazMIdgFeCJxfVQcDj7LN\nYaCqKqDGW7mqllbV4qpavGDBghbLlKRuazMIHgAeqKqvN/NX0guGHyfZD6B5XtdiDZKkSbQWBFX1\nI+AHSRY1TUcC3wI+C7ylaXsLcE1bNUiSJrdLy+9/CnBpkicB9wIn0AufK5L8JXA/cEzLNUiSJtBq\nEFTVHcDicV46ss3tSpL655XFktRxBoEkdZxBIEnDaNUV8MA34P6vwHm/35tviUEgScNm1RXwuXfC\nyK968xt/0JtvKQwMAk3dNP5SkTrpxr+DxzZt3fbYpl57CwwCTc00/1KROmnjA1Nr30EGgaZmmn+p\nSJ201wFTa99BBoGmZpp/qUiddOQZsOvcrdt2ndtrb4FBoKmZ5l8qUicdeAy89p9gzpN783s9szd/\nYDsDMRgEmppp/qUiddaBx8ABL4JnHw7vWd1aCIBBoKma5l8qktrX9qBzM9No98iRX/W6Rx55hn/o\nxjrwGFjR3GTuhGsHW4ukHdbXHkGSfZNcnOSLzfy/bUYPnX3sHimpY/o9NPQJ4HrgGc38d4B3t1HQ\nwNk9UlLH9BsE+1TVFcAWgKp6HBhprapBsnuk1D6vTh8q/QbBo0meRnN/4SR/CGxsrapBsnuk1C4P\nvw6dfoPgVHq3mHxukq8Cn6R397HZx+6RUrs8/Dp0+uo1VFW3J3kpsAgIcE9VPdZqZYMy2jvomnf0\nfrHs9Ux7DUk7k4dfh06/vYZOBnavqruqajWwe5K3t1vaAE3jhRxS53j4dej0e2jorVW1YXSmqn4G\nvLWdkiTNah5+HTr9BsGcJBmdSTIHeFI7JUma1bw6fej0e2XxdcDlSS5o5k9q2iRp6rw6faj0GwR/\nTe+P/39q5m8ALmqlIknStOq319AW4PzmIUmaRfoKgiSHAWcCz27WCVBV9Zz2SpMkTYd+Dw1dDLwH\nWMFsHVpC2pkcwVYzSL9BsLGqvthqJdJs8URDKIBhoKHUb/fRm5OcneTQJC8cfbRamTRTOYSCZph+\n9whe3DwvHtNWwMt3bjnSLOAQCpph+u019LK2C9nZlq1cy8rvb2DzyBYOO+smTjtqEUsO3n/QZakL\n9jqgdzhovHZpCPV9q8okrwaeD+w22lZVQ7mvu2zlWk6/6k42j2wBYO2GTZx+1Z0AhoHad+QZvXMC\nYw8POYSChli/g859FDiW3tDTAf49va6kQ+ns6+9h02Nbd27a9NgIZ19/z4AqUqc4hIJmmH73CP6o\nqg5MsqqqPpjkXGBoexH9cMOmKbVLO51DKGgG6bfX0Ohf0F8keQbwGLBfOyXtuGfMnzuldknqsn6D\n4PNJ5gNnA7cDa4DL2ipqR5121CLm7jpnq7a5u87htKMWDagizVTHXnArx15w66DLkFrVb6+hDzWT\nn0nyeWC3qpr0nsVJ1gAP07sa+fGqWpzkTHr3MljfLPY3VfWFqRY+kdETwu+7chWbR7aw//y59hqS\npCcwYRAkeXlV3ZTk9eO8RlVd1cc2XlZVD23Tdl5VnTOVQqdqycH7c9lt3wfg8pMObXNTkjSjTbZH\n8FLgJuC147xWQD9BIEkaYhMGQVV9oJn8q6ransHmCvhSkhHggqpa2rSfkuQvgOXAf25ufbmVJCcC\nJwI861nP2o5NS5L60e/J4vuSLE1y5NhbVvbh8Ko6CPhT4OQkR9C7p8FzgIOAB4Fzx1uxqpZW1eKq\nWrxgwYIpbFKSNBX9BsHzgC8BJ9MLhf+Z5PDJVqqqtc3zOuBq4JCq+nFVjTQ3u7kQOGT7Spck7Qx9\nBUFV/aKqrqiq1wMHA3sC/2eidZLMS7LH6DTwSmB1krHXHxwNrN6uyrXDRsdj+vp9P+Wws25i2cq1\ngy5J0gBMZayhl9IbZuJP6B3bn+x6+X2Bq5sjSbsA/6uqrkvyqSQH0Tt/sIbevZA1zRyPSdKofm9V\nuQZYCVwBnFZVj062TlXdC7xgnPbjplijWjDReEwGgdQt/e4RHFhVP2+1Ek0rx2PSzjJ65bXX68xc\n/Z4s/ldJbkyyGiDJgUn+tsW61DLHY5I0qt8guBA4nd5gc1TVKuCNbRWl9jkek6RR/R4aekpV3bbN\nJQSPt1CPponjMUka1W8QPJTkufR6+pDkDfQuBtMM5nhMkqD/IDgZWAo8L8la4D7gTa1VJUmaNpON\nPnrqmNkvADfTO6/wKPDnwIfbK02SNB0m2yPYo3leBLwIuIbePYuPA25rsS5J0jSZbPTRDwIk+TLw\nwqp6uJk/E/BGrJI0iZlwnUW/3Uf3BTaPmd/ctEmSZrh+TxZ/ErgtydXN/BLgE61UJEmaVv3es/gf\nknwReEnTdEJVrWyvLEnSdOl79NGquh24vcVaJEkD0O85AknSLGUQSFLHGQSS1HEGgSR1nEEgSR1n\nEEhSxxkEktRxBoEkdZxBIEkdZxBIUscZBJLUcQaBJHWcQSBJHWcQSFLH9T0MtbSVE7xTqTRbuEcg\nSR1nEEhPYNnKtaz8/ga+ft9POeysm1i2cu2gS5JaYRBI41i2ci2nX3Unm0e2ALB2wyZOv+pOw0Cz\nkucInojHwDvt7OvvYdNjI1u1bXpshLOvv4clB+8/oKqkdrhHII3jhxs2TaldmskMAmkcz5g/d0rt\n0kxmEEjjOO2oRczddc5WbXN3ncNpRy0aUEXqpBOunZbD1K2eI0iyBngYGAEer6rFSfYGLgcWAmuA\nY6rqZ23WIU3V6HmA9125is0jW9h//lxOO2qR5wc0K03HyeKXVdVDY+bfD9xYVWcleX8z/9fTUIc0\nJUsO3p/Lbvs+AJefdOiAq5HaM4hDQ68DLmmmLwGWDKAGSTuB11rMDm0HQQFfSrIiyYlN275V9WAz\n/SNg3/FWTHJikuVJlq9fv77lMiVNlddazB5tHxo6vKrWJnk6cEOSb499saoqSY23YlUtBZYCLF68\neNxlJA3ODl9r4bU6Q6PVPYKqWts8rwOuBg4BfpxkP4DmeV2bNUhqh9dazB6tBUGSeUn2GJ0GXgms\nBj4LvKVZ7C3ANW3VIKk9Xmsxe7S5R7Av8JUk3wRuA66tquuAs4BXJPku8MfNvKQZxmstZo/WzhFU\n1b3AC8Zp/wlwZFvbHcsufxqoWX4M3GstZg8HnZO03bzWYmKj3Ws3j2zhsLNuGtqgdIgJSWrBTOpe\naxBIUgsm6l47bAwCSWrBTOpeaxBIUgtmUvdag0CSWjCTutfaa0iSWjCTutcaBJLUkpnSvdZDQ5LU\nce4RdNww/0qRND3cI5CkjjMIJKnjDAJJ6jiDQJI6ziCQpI4zCCSp4wwCSeo4g0CSOs4gkKSOMwgk\nqeMMAknqOINAkjrOIJCkjjMIJKnjDAJJ6jiDQJI6ziCQpI4zCCSp4wwCSeo4g0CSOs6b10sTuPyk\nQwddgtQ69wgkqeMMAknqOINAkjrOIJCkjms9CJLMSbIyyeeb+TOTrE1yR/N4Vds1SJKe2HT0GnoX\ncDew55i286rqnGnYtiRpEq0GQZIDgFcD/wCc2ua2JA2GXWxnvrYPDf0j8D5gyzbtpyRZleRjSZ46\n3opJTkyyPMny9evXt1ymJHVXa0GQ5DXAuqpasc1L5wPPAQ4CHgTOHW/9qlpaVYuravGCBQvaKlOS\nOq/NQ0OHAX/WnAzeDdgzyaer6s2jCyS5EPh8izVIkibR2h5BVZ1eVQdU1ULgjcBNVfXmJPuNWexo\nYHVbNUiSJjeIsYb+R5KDgALWACcNoAZJmhYz4WT6tARBVd0C3NJMHzcd25Qk9ccriyWp4wwCSeo4\ng0CSOs4gkKSOMwgkqeMMAknqOINAkjrOIJCkjktVDbqGSSVZD9y/navvAzy0E8uZjfyOJub3Mzm/\no4kN6vt5dlVNOmrnjAiCHZFkeVUtHnQdw8zvaGJ+P5PzO5rYsH8/HhqSpI4zCCSp47oQBEsHXcAM\n4Hc0Mb+fyfkdTWyov59Zf45AkjSxLuwRSJImYBBIUsfN2iBI8rEk65J4K8xxJHlmkpuTfCvJXUne\nNeiahk2S3ZLcluSbzXf0wUHXNIySzEmyMon3Hx9HkjVJ7kxyR5Llg65nPLP2HEGSI4BHgE9W1e8P\nup5h09w7er+quj3JHsAKYElVfWvApQ2NJAHmVdUjSXYFvgK8q6q+NuDShkqSU4HFwJ5V9ZpB1zNs\nkqwBFlfV0F5wN2v3CKrqy8BPB13HsKqqB6vq9mb6YeBuYP/BVjVcqueRZnbX5jE7fzltpyQHAK8G\nLhp0Ldp+szYI1L8kC4GDga8PtpLh0xz2uANYB9xQVX5HW/tH4H3AlkEXMsQK+FKSFUlOHHQx4zEI\nOi7J7sBngHdX1c8HXc+wqaqRqjoIOAA4JImHGRtJXgOsq6oVg65lyB3e/B/6U+Dk5rD1UDEIOqw5\n7v0Z4NKqumrQ9QyzqtoA3Az8yaBrGSKHAX/WHAP/38DLk3x6sCUNn6pa2zyvA64GDhlsRb/NIOio\n5kToxcDdVfXhQdczjJIsSDK/mZ4LvAL49mCrGh5VdXpVHVBVC4E3AjdV1ZsHXNZQSTKv6YxBknnA\nK4Gh68k4a4MgyWXArcCiJA8k+ctB1zRkDgOOo/cr7o7m8apBFzVk9gNuTrIK+Aa9cwR2kdRU7At8\nJck3gduAa6vqugHX9FtmbfdRSVJ/Zu0egSSpPwaBJHWcQSBJHWcQSFLHGQSS1HEGgWasJLckaf2G\n4EnemeTuJJe2va2dJcnfDLoGzRwGgTopyS5TWPztwCuq6k07uYY5O/P9tjHlIGi5Hg0xg0CtSrKw\n+TV9YTOm/780V+lu9Ys+yT7NUAUkOT7JsiQ3NGO5vyPJqc2Y919LsveYTRzXXAy3OskhzfrzmvtR\n3Nas87ox7/vZJDcBN45T66nN+6xO8u6m7aPAc4AvJnnPNssfn+Sa5nN8N8kHxry2rBlk7K6xA40l\neSTJuc0FRocmOSPJN5ptLm2u+B79bs5Lsrz5/l6U5KpmO38/5v3e3HzOO5Jc0AySdxYwt2m79ImW\ne4J6zkrvHhWrkpyzff/qmnGqyoeP1h7AQuBx4KBm/grgzc30LfTGaQfYB1jTTB8PfA/YA1gAbATe\n1rx2Hr0B8kbXv7CZPgJY3Uz/1zHbmA98B5jXvO8DwN7j1PkHwJ3NcrsDdwEHN6+tAfYZZ53jgQeB\npwFz6Q0dMPp59m6eR9uf1swXcMyY99h7zPSngNeO+Wz/vZl+F/BDelc6P7n5DE8D/g3wOWDXZrl/\nBv6imX5kzPtOtNyv62ne8x5+c6Hp/EH///ExPY+p7B5L2+u+qrqjmV5BLxwmc3P17pPwcJKN9P6Q\nQe+P9YFjlrsMevefSLJnMzbQK+kNhvbeZpndgGc10zdU1Xj3qTgcuLqqHgVIchXwEmDlJHXeUFU/\nGbPO4cBy4J1Jjm6WeSbwe8BPgBF6A/2NelmS9wFPAfamF0Cjn/WzYz7zXVX1YLOde5v3PJxegH2j\n2ZGYS2+47G0dOcFyY+vZCPwSuDi9u405nEZHGASaDr8aMz1C7w8R9PYURg9P7jbBOlvGzG9h6/+3\n246RUkCAP6+qe8a+kOTFwKNTqnxyv7X9JP8O+GPg0Kr6RZJb+M3n+2VVjTT17Ebv1/niqvpBkjPZ\n+nsY+5m3/T52ofc5L6mq0yepcaLlfl1PVT3eHF47EngD8A7g5ZO8t2YBzxFokNbQ+6UKvT882+NY\ngCSHAxuraiNwPXDKmOPtB/fxPv8XWJLkKc0okUc3bZN5RZK9m/MeS4CvAnsBP2tC4HnAHz7BuqN/\n9B9K774QU/0ObgTekOTpAE0dz25eeyy9YcYnW+7Xmhr2qqovAO8BXjDFejRDuUegQToHuKI5mXrt\ndr7HL5OspHcbyf/YtH2I3p2zViX5HeA+YMJ76Vbv3s2foDdCJMBFVTXZYSGa5T9D78Y1n66q5Unu\nBN6W5G56x9zHvcdxVW1IciG9cwg/ojfCad+q6ltJ/hb4l+ZzPgacDNwPLKX3+W+vqjdNsNxYewDX\nNHsqAU6dSj2auRx9VNpOSY6nd1jnHYOuRdoRHhqSpI5zj0CSOs49AknqOINAkjrOIJCkjjMIJKnj\nDAJJ6rj/D2AflEruiNekAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x111b3a438>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(df.index, df.train, 'o', color = colors[0])\n",
    "plt.errorbar(df.index, df.train, yerr=df.train_std, fmt=None, color=colors[0])\n",
    "plt.plot(df.index+0.2, df.test, 'o', color = colors[1])\n",
    "plt.errorbar(df.index+0.2, df.test, yerr=df.test_std, fmt=None, color=colors[1])\n",
    "plt.xlabel(\"number of parameters\")\n",
    "plt.ylabel(\"deviance\")\n",
    "plt.title(\"N=20\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice:\n",
    "\n",
    "- the best fit model may not be the original generating model. Remember that the choice of fit depends on the amount of data you have and the less data you have, the less parameters you should use\n",
    "- on average, out of sample deviance must be larger than in-sample deviance, through an individual pair may have that order reversed because of sample peculiarity.\n",
    "\n",
    "## AIC, or the difference in deviances\n",
    "\n",
    "Let us see the difference between the mean testing and training deviances. This is the difference in *bias* between the two sets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1    2.018779\n",
       "2    4.212654\n",
       "3    5.441971\n",
       "4    7.670368\n",
       "5    9.691897\n",
       "dtype: float64"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.test - df.train"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voila, this seems to be roughly twice the number of parameters. In other words we might be able to get away without a test set if we \"correct\" the bias on the traing set by $2n_p$. This is the observation that motivates the AIC."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analysis N=100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "reps=10000\n",
    "results_100 = {}\n",
    "for k in range(1,6):\n",
    "    trdevs=np.zeros(reps)\n",
    "    tedevs=np.zeros(reps)\n",
    "    for r in range(reps):\n",
    "        Xtr, ytr, Xte, yte = generate_data(100, k)\n",
    "        ols = sm.OLS(ytr, Xtr).fit()\n",
    "        mutr = np.dot(Xtr, ols.params)\n",
    "        devtr = -2*np.sum(norm.logpdf(ytr, mutr, 1))\n",
    "        mute = np.dot(Xte, ols.params)\n",
    "        devte = -2*np.sum(norm.logpdf(yte, mute, 1))\n",
    "        #print(k, r, devtr, devte)\n",
    "        trdevs[r] = devtr\n",
    "        tedevs[r] = devte\n",
    "    results_100[k] = (np.mean(trdevs), np.std(trdevs), np.mean(tedevs), np.std(tedevs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "df100 = pd.DataFrame(results_100).T\n",
    "df100 = df100.rename(columns = dict(zip(range(4), ['train', 'train_std', 'test', 'test_std'])))\n",
    "df100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.plot(df100.index, df100.train, 'o', color = colors[0])\n",
    "plt.errorbar(df100.index, df100.train, yerr=df100.train_std, fmt=None, color=colors[0])\n",
    "plt.plot(df100.index+0.2, df100.test, 'o', color = colors[1])\n",
    "plt.errorbar(df100.index+0.2, df100.test, yerr=df100.test_std, fmt=None, color=colors[1])\n",
    "plt.xlabel(\"number of parameters\")\n",
    "plt.ylabel(\"deviance\")\n",
    "plt.title(\"N=100\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 234,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1    1.900383\n",
       "2    3.908291\n",
       "3    5.193728\n",
       "4    7.003606\n",
       "5    8.649499\n",
       "dtype: float64"
      ]
     },
     "execution_count": 234,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df100.test - df100.train"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We get pretty much the same result at N=100."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Assumptions for AIC\n",
    "\n",
    "This observation leads to an estimate of the out-of-sample deviance by what is called an **information criterion**, the Akaike Information Criterion, or AIC:\n",
    "\n",
    "$$AIC = D_{train} + 2n_p$$\n",
    "\n",
    "which does carry as assumptions that\n",
    "\n",
    "1. the likelihood is approximately multivariate gaussian\n",
    "2. the sample size is much larger than the number of parameters\n",
    "3. priors are flat \n",
    "4. The AIC does not assume that the true data generating process $p$ is in the set of models being fitted. The overarching goal of the AIC approach to model selection is to select the \"best\" model for our given data set without assuming that the \"true\" model is in the family of models from which we're selecting. The true model \"cancels out\" except in the expectation.\n",
    "\n",
    "We wont derive the AIC here, but if you are interested, see  http://www.stat.cmu.edu/~larry/=stat705/Lecture16.pdf\n",
    "\n",
    "Why would we want to use such information criteria? Cross validation can be expensive, especially with multiple hyper-parameters.\n",
    "\n",
    "## AIC for Linear Regression\n",
    "\n",
    "The AIC for a model is the training deviance plus twice the number of parameters:\n",
    "\n",
    "$$AIC = D_{train} + 2n_p.$$\n",
    "\n",
    "That is, -2 times the log likelihood of the model.\n",
    "\n",
    "So, one we find the MLE solution for the linear regression, we plugin the values we get, which are\n",
    "\n",
    "$$\\sigma_{MLE}^2 =  \\frac{1}{N} RSS $$\n",
    "\n",
    "where RSS is the sum of the  squares of the errors.\n",
    "\n",
    "$$AIC = -2(-\\frac{N}{2}(log(2\\pi) + log(\\sigma^2)) -2(-\\frac{1}{2\\sigma_{MLE}^2} \\times RSS) + 2p$$\n",
    "\n",
    "Thus:\n",
    "\n",
    "$$D = Nlog(RSS/N) $$\n",
    "\n",
    "$$AIC = Nlog(RSS/N) + 2p + constant$$\n",
    "\n",
    "Since the deviance for a OLS model is just proportional to the log(MSE) upto a proportionality, we'll use the MSE to derive this split.\n",
    "\n",
    "The fact that the (log-likelihood) and thus the deviance  carries an expectation over the true distribution as estimated on the sample means that the **Deviance is a stochastic quantity, varying from sample to sample**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A complete understanding of the comparison diagram\n",
    "\n",
    "(taken from McElreath, but see upstairs as well)\n",
    "\n",
    "![](../images/inoutdeviance.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are equipped to understand this diagram completely. Lets focus on the training (in) set first: blue points. \n",
    "\n",
    "1. There is some irreducible noise which contributes to the deviance no matter the number of parameters.\n",
    "2. If we could capture the true model exactly there would be no **bias**, and the deviance would go to that which comes from the irreducible noise.\n",
    "3. But we cant, so the positions of the circles tells us how much bias plus irreducible noise we have\n",
    "4. The error bars now tell us our **variance**, since they tell us how much our deviance, or MSE varies around our \"mean\" model. In real life our sample will lie somewhere along this error bar.\n",
    "5. The training set deviances go down as the number of parameters increase. The test set deviances go down and then go up\n",
    "6. Notice that testing deviance is higher on a 2 parameter model than on a 1, even though our generating \"true\" model is a 2 parameter one. Deviance and the AIC do not pick the true model, but rather the one with the highest predictive accuracy."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}