32_MoreModelEvaluation_PyMC3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=float:right><img src=\"assets/img/appliedai-logo.png\" width=100 style=\"margin: 0px 20px\"></img></div>\n",
    "\n",
    "\n",
    "##### Jonathan Sedar Personal Project\n",
    "## PyMC3 vs PyStan Comparison\n",
    "_Spring 2016_\n",
    "\n",
    "This set of Notebooks and scripts comprise the **pymc3_vs_pystan** personal project by Jonathan Sedar of Applied AI Ltd, written primarily for presentation at the PyData London 2016 Conference.\n",
    "\n",
    "The project demonstrates hierarchical linear regression using two Bayesian inference frameworks: PyMC3 and PyStan. The project borrows heavily from code written for Applied AI Ltd and is supplied here for educational purposes only. No copyright or license is extended to users.\n",
    "\n",
    "\n",
    "    \n",
    "# 32_MoreModelEvaluation_PyMC3\n",
    "\n",
    "#### Demonstrate some model evaluation built into PyMC3\n",
    "\n",
    "      \n",
    "+ [Setup](#Setup)\n",
    "    + [Local Functions](#Local-Functions)\n",
    "    + [Load Data](#Load-Data)\n",
    "    + [Describe Dataset](#Describe-Dataset)\n",
    "\n",
    "\n",
    "+ [Reload Models and Traces](#Reload-Models-and-Traces)\n",
    "    + [Reload PyMC3 model and traces](#Reload-PyMC3-model-and-traces)\n",
    "    + [Reload PyStan fitted model](#Reload-PyStan-fitted-model)\n",
    "\n",
    "\n",
    "+ [Information Criteria](#Information-Criteria)\n",
    "    + [Deviance Information Criterion](Deviance-Information-Criterion)  \n",
    "    + [Widely-Accepted Information Criterion](#Widely-Accepted-Information-Criterion)  \n",
    "\n",
    "\n",
    "+ [Posterior Predictive Checks on Validation Data](#Posterior-Predictive-Checks-on-Validation-Data)\n",
    "    + [Held-out test set](#Held-out-test-set)\n",
    "    + [Cross-Validation](#Cross-Validation)\n",
    "    + [Leave-One-Out Cross-Val](#Leave-One-Out-Cross-Val)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Note:\n",
    "\n",
    "For more information on Model Selection in PyMC3, and about DIC and WAIC, you could start with:\n",
    "+ Thomas Wiecki's [detailed response](https://stats.stackexchange.com/questions/161082/bayesian-model-selection-in-pymc3/166383#166383) to a question on Cross Validated\n",
    "+ The Deviance Information Criterion: 12 Years On ([Speigelhalter et al 2014](http://onlinelibrary.wiley.com/doi/10.1111/rssb.12062/abstract))\n",
    "+ A Widely Applicable Bayesian Information Criterion ([Watanabe 2013](http://www.jmlr.org/papers/volume14/watanabe13a/watanabe13a.pdf))\n",
    "+ Efficient Implementation of Leave-One-Out Cross-Validation and WAIC for Evaluating Fitted Bayesian Models ([Gelman et al 2015](http://arxiv.org/abs/1507.04544))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## Interactive magics\n",
    "%matplotlib inline\n",
    "%qtconsole --colors=linux"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# filter warnings for presentation's sake\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "# general packages\n",
    "import sqlite3\n",
    "from convenience_functions import *\n",
    "from ipywidgets import interactive, fixed\n",
    "from itertools import combinations\n",
    "#from io import StringIO\n",
    "#from collections import OrderedDict\n",
    "\n",
    "\n",
    "# scientific packages\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import patsy as pt\n",
    "from scipy import optimize\n",
    "from scipy.stats import ks_2samp\n",
    "\n",
    "#from sklearn.neighbors.kde import KernelDensity\n",
    "import statsmodels.api as sm\n",
    "\n",
    "# pymc3 libraries\n",
    "import pymc3 as pm\n",
    "import theano as thno\n",
    "import theano.tensor as T \n",
    "import pystan\n",
    "\n",
    "sns.set(style=\"darkgrid\", palette=\"muted\")\n",
    "pd.set_option('display.mpl_style', 'default')\n",
    "plt.rcParams['figure.figsize'] = 12, 4\n",
    "np.random.seed(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Local Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def strip_derived_rvs(rvs):\n",
    "    '''Convenience fn: remove PyMC3-generated RVs from a list'''\n",
    "    ret_rvs = []\n",
    "    for rv in rvs:\n",
    "        if not (re.search('_log',rv.name) or re.search('_interval',rv.name)):\n",
    "            ret_rvs.append(rv)     \n",
    "    return ret_rvs\n",
    "\n",
    "\n",
    "def trace_median(x):\n",
    "    return pd.Series(np.median(x,0), name='median')\n",
    "\n",
    "\n",
    "def plot_traces_pymc(trcs, varnames=None):\n",
    "    ''' Convenience fn: plot traces with overlaid means and values '''\n",
    "\n",
    "    nrows = len(trcs.varnames)\n",
    "    if varnames is not None:\n",
    "        nrows = len(varnames)\n",
    "    ax = pm.traceplot(trcs, varnames=varnames, figsize=(12,nrows*1.4)\n",
    "        ,lines={k: v['mean'] for k, v in \n",
    "            pm.df_summary(trcs,varnames=varnames).iterrows()})\n",
    "\n",
    "    for i, mn in enumerate(pm.df_summary(trcs, varnames=varnames)['mean']):\n",
    "        ax[i,0].annotate('{:.2f}'.format(mn), xy=(mn,0), xycoords='data'\n",
    "                    ,xytext=(5,10), textcoords='offset points', rotation=90\n",
    "                    ,va='bottom', fontsize='large', color='#AA0022')    \n",
    "\n",
    "        \n",
    "def plot_stan_trc(dftrc):\n",
    "    \"\"\"\n",
    "       Create simple plots of parameter distributions and traces from \n",
    "       output of pystan sampling. Emulates pymc traceplots.\n",
    "    \"\"\"\n",
    "\n",
    "    fig, ax2d = plt.subplots(nrows=dftrc.shape[1], ncols=2, figsize=(14, 1.8*dftrc.shape[1]),\n",
    "                                facecolor='0.99', edgecolor='k')\n",
    "    fig.suptitle('Distributions and traceplots for {} samples'.format(\n",
    "                                dftrc.shape[0]),fontsize=14)\n",
    "    fig.subplots_adjust(wspace=0.2, hspace=0.5)\n",
    "\n",
    "    k = 0\n",
    "    \n",
    "    # create density and traceplot, per parameter coeff\n",
    "    for i, (ax1d, col) in enumerate(zip(ax2d, dftrc.columns)):\n",
    "\n",
    "        samples = dftrc[col].values\n",
    "        scale = (10**np.round(np.log10(samples.max() - samples.min()))) / 20\n",
    "        kde = KernelDensity(bandwidth=scale).fit(samples.reshape(-1, 1))\n",
    "        x = np.linspace(samples.min(), samples.max(), 100).reshape(-1, 1)\n",
    "        y = np.exp(kde.score_samples(x))\n",
    "        clr = sns.color_palette()[0]\n",
    "\n",
    "        # density plot\n",
    "        ax1d[0].plot(x, y, color=clr, linewidth=1.4)\n",
    "        ax1d[0].vlines(np.percentile(samples, [2.5, 97.5]), ymin=0, ymax=y.max()*1.1,\n",
    "                       alpha=1, linestyles='dotted', colors=clr, linewidth=1.2)\n",
    "        mn = np.mean(samples)\n",
    "        ax1d[0].vlines(mn, ymin=0, ymax=y.max()*1.1,\n",
    "                       alpha=1, colors='r', linewidth=1.2)\n",
    "        ax1d[0].annotate('{:.2f}'.format(mn), xy=(mn,0), xycoords='data'\n",
    "                    ,xytext=(5,10), textcoords='offset points', rotation=90\n",
    "                    ,va='bottom', fontsize='large', color='#AA0022')    \n",
    "        ax1d[0].set_title('{}'.format(col), fontdict={'fontsize':10})\n",
    "\n",
    "\n",
    "        # traceplot\n",
    "        ax1d[1].plot(np.arange(len(samples)),samples, alpha=0.2, color=clr, linestyle='solid'\n",
    "                              ,marker=',', markerfacecolor=clr, markersize=10)\n",
    "        ax1d[1].hlines(np.percentile(samples,[2.5, 97.5]), xmin=0, xmax=len(samples),\n",
    "                       alpha=1, linestyles='dotted', colors=clr)\n",
    "        ax1d[1].hlines(np.mean(samples), xmin=0, xmax=len(samples), alpha=1, colors='r')\n",
    "\n",
    "        k += 1\n",
    "                \n",
    "        ax1d[0].set_title('{}'.format(col), fontdict={'fontsize':14})#,'fontweight':'bold'})\n",
    "        #ax1d[0].legend(loc='best', shadow=True)\n",
    "        \n",
    "        _ = [ax1d[j].axes.grid(True, linestyle='-', color='lightgrey') for j in range(2)]\n",
    "            \n",
    "    plt.subplots_adjust(top=0.94)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "cnxsql = sqlite3.connect('data/car_emissions.db')\n",
    "dfs = pd.read_sql('select * from cars_post_exclusions_2sd', cnxsql, index_col=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## convert sqlite bool storage (as ints) back to bools\n",
    "for ft in ['parent_is_vw', 'mfr_is_vw', 'is_tdi']:\n",
    "    dfs[ft] = dfs[ft].astype(bool)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2593, 13)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>1653</th>\n",
       "      <th>835</th>\n",
       "      <th>763</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",
       "      <th>dtype</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>emissions_nox_mgkm</th>\n",
       "      <td>12</td>\n",
       "      <td>36</td>\n",
       "      <td>32</td>\n",
       "      <td>2593</td>\n",
       "      <td>37.32</td>\n",
       "      <td>17.9</td>\n",
       "      <td>1</td>\n",
       "      <td>23.000000</td>\n",
       "      <td>35.000000</td>\n",
       "      <td>51.000000</td>\n",
       "      <td>76</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>parent_is_vw</th>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>2593</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>False</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>True</td>\n",
       "      <td>bool</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mfr_is_vw</th>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>2593</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>False</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>True</td>\n",
       "      <td>bool</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>parent</th>\n",
       "      <td>daimler-ag</td>\n",
       "      <td>bmw</td>\n",
       "      <td>bmw</td>\n",
       "      <td>2593</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>aston</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>volksw</td>\n",
       "      <td>object</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mfr</th>\n",
       "      <td>mercedes-benz</td>\n",
       "      <td>bmw</td>\n",
       "      <td>bmw</td>\n",
       "      <td>2593</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>abarth</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>volvo</td>\n",
       "      <td>object</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>trans</th>\n",
       "      <td>auto</td>\n",
       "      <td>semiauto</td>\n",
       "      <td>auto</td>\n",
       "      <td>2593</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>auto</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>semiau</td>\n",
       "      <td>object</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fuel_type</th>\n",
       "      <td>petrol</td>\n",
       "      <td>petrol</td>\n",
       "      <td>petrol</td>\n",
       "      <td>2593</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>diesel</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>petrol</td>\n",
       "      <td>object</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>is_tdi</th>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>2593</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>False</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>True</td>\n",
       "      <td>bool</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>metric_combined</th>\n",
       "      <td>-0.0728208</td>\n",
       "      <td>0.80692</td>\n",
       "      <td>0.220426</td>\n",
       "      <td>2593</td>\n",
       "      <td>-0.00</td>\n",
       "      <td>0.5</td>\n",
       "      <td>-0.685973</td>\n",
       "      <td>-0.339409</td>\n",
       "      <td>-0.152797</td>\n",
       "      <td>0.167108</td>\n",
       "      <td>2.75301</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>metric_extra_urban</th>\n",
       "      <td>-0.0751821</td>\n",
       "      <td>0.47462</td>\n",
       "      <td>0.0728415</td>\n",
       "      <td>2593</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.5</td>\n",
       "      <td>-0.47696</td>\n",
       "      <td>-0.180913</td>\n",
       "      <td>-0.075182</td>\n",
       "      <td>0.093988</td>\n",
       "      <td>21.5997</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>metric_urban_cold</th>\n",
       "      <td>-0.0249258</td>\n",
       "      <td>0.724487</td>\n",
       "      <td>0.291493</td>\n",
       "      <td>2593</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.5</td>\n",
       "      <td>-0.657763</td>\n",
       "      <td>-0.308037</td>\n",
       "      <td>-0.158155</td>\n",
       "      <td>0.158264</td>\n",
       "      <td>2.92276</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>engine_capacity</th>\n",
       "      <td>-0.259591</td>\n",
       "      <td>0.502214</td>\n",
       "      <td>-0.0383156</td>\n",
       "      <td>2593</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.5</td>\n",
       "      <td>-0.920116</td>\n",
       "      <td>-0.278857</td>\n",
       "      <td>-0.045471</td>\n",
       "      <td>-0.037765</td>\n",
       "      <td>2.57901</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>emissions_co_mgkm</th>\n",
       "      <td>-0.325803</td>\n",
       "      <td>-0.456874</td>\n",
       "      <td>0.0567836</td>\n",
       "      <td>2593</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.5</td>\n",
       "      <td>-0.846545</td>\n",
       "      <td>-0.382482</td>\n",
       "      <td>-0.088458</td>\n",
       "      <td>0.272874</td>\n",
       "      <td>2.28145</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                             1653       835        763  count   mean   std  \\\n",
       "emissions_nox_mgkm             12        36         32   2593  37.32  17.9   \n",
       "parent_is_vw                False     False      False   2593    NaN   NaN   \n",
       "mfr_is_vw                   False     False      False   2593    NaN   NaN   \n",
       "parent                 daimler-ag       bmw        bmw   2593    NaN   NaN   \n",
       "mfr                 mercedes-benz       bmw        bmw   2593    NaN   NaN   \n",
       "trans                        auto  semiauto       auto   2593    NaN   NaN   \n",
       "fuel_type                  petrol    petrol     petrol   2593    NaN   NaN   \n",
       "is_tdi                      False     False      False   2593    NaN   NaN   \n",
       "metric_combined        -0.0728208   0.80692   0.220426   2593  -0.00   0.5   \n",
       "metric_extra_urban     -0.0751821   0.47462  0.0728415   2593   0.00   0.5   \n",
       "metric_urban_cold      -0.0249258  0.724487   0.291493   2593   0.00   0.5   \n",
       "engine_capacity         -0.259591  0.502214 -0.0383156   2593   0.00   0.5   \n",
       "emissions_co_mgkm       -0.325803 -0.456874  0.0567836   2593   0.00   0.5   \n",
       "\n",
       "                         min        25%        50%        75%      max  \\\n",
       "emissions_nox_mgkm         1  23.000000  35.000000  51.000000       76   \n",
       "parent_is_vw           False        NaN        NaN        NaN     True   \n",
       "mfr_is_vw              False        NaN        NaN        NaN     True   \n",
       "parent                aston         NaN        NaN        NaN   volksw   \n",
       "mfr                   abarth        NaN        NaN        NaN    volvo   \n",
       "trans                   auto        NaN        NaN        NaN   semiau   \n",
       "fuel_type             diesel        NaN        NaN        NaN   petrol   \n",
       "is_tdi                 False        NaN        NaN        NaN     True   \n",
       "metric_combined    -0.685973  -0.339409  -0.152797   0.167108  2.75301   \n",
       "metric_extra_urban  -0.47696  -0.180913  -0.075182   0.093988  21.5997   \n",
       "metric_urban_cold  -0.657763  -0.308037  -0.158155   0.158264  2.92276   \n",
       "engine_capacity    -0.920116  -0.278857  -0.045471  -0.037765  2.57901   \n",
       "emissions_co_mgkm  -0.846545  -0.382482  -0.088458   0.272874  2.28145   \n",
       "\n",
       "                      dtype  \n",
       "emissions_nox_mgkm  float64  \n",
       "parent_is_vw           bool  \n",
       "mfr_is_vw              bool  \n",
       "parent               object  \n",
       "mfr                  object  \n",
       "trans                object  \n",
       "fuel_type            object  \n",
       "is_tdi                 bool  \n",
       "metric_combined     float64  \n",
       "metric_extra_urban  float64  \n",
       "metric_urban_cold   float64  \n",
       "engine_capacity     float64  \n",
       "emissions_co_mgkm   float64  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "custom_describe(dfs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Declare feats for use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fts_cat = ['parent_is_vw', 'mfr_is_vw', 'parent', 'mfr', 'trans', 'fuel_type', 'is_tdi']\n",
    "fts_cat_smp = ['mfr_is_vw','trans','fuel_type','is_tdi']\n",
    "fts_num = ['metric_combined', 'metric_extra_urban', 'metric_urban_cold'\n",
    "           ,'engine_capacity', 'emissions_co_mgkm']\n",
    "fts_num_smp = ['metric_combined', 'engine_capacity', 'emissions_co_mgkm']\n",
    "ft_endog = 'emissions_nox_mgkm'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Describe dataset\n",
    "\n",
    "+ The dataset is 2593 rows, with 12 exog features, 1 endog feature.\n",
    "+ These are observations of car emissions tests, one row per car.\n",
    "+ You can read off the basic distributional statistics of the features in the table above. Numeric features have been standardized according to [Gelman's 2sd principle](http://www.stat.columbia.edu/~gelman/research/published/standardizing7.pdf).\n",
    "+ I have selected these particular 12 features to work with. Some are derivatives of original features.\n",
    "\n",
    "We have the following features to choose from:\n",
    "\n",
    "```\n",
    "+ Categoricals:\n",
    "    + `trans`     - the car transmission, simplified to 'auto', 'semiauto', 'manual'\n",
    "    + `fuel_type` - the car power supply, simplified to 'petrol', 'diesel'\n",
    "    + `parent`    - the parent company of the car manufacturer, 20 values\n",
    "    + `mfr`       - the car manufacturer, 38 values\n",
    "\n",
    "+ Booleans:\n",
    "    + `parent_is_vw` - if the parent company of the car manufacturer is Volkswagen\n",
    "    + `mfr_is_vw`    - if the car manufacturer is Volkswagen\n",
    "    + `is_tdi`       - (processed feature) if the car engine type is a turbo diesel\n",
    "    \n",
    "+ Numerics:\n",
    "    + `metric_combined`    - a score for fuel efficiency in combined driving\n",
    "    + `metric_extra_urban` - a score for fuel efficiency in an extra-urban driving\n",
    "    + `metric_urban_cold`  - a score for fuel efficiency in an urban setting, cold start\n",
    "    + `emissions_co_mgkm`  - a count of CO particulates emitted mg/km\n",
    "    \n",
    "+ Numeric endogenous feature:\n",
    "    + `emissions_nox_mgkm` - a count of NOx particulates emitted mg/km    \n",
    "```\n",
    "\n",
    "For the purposes of this Notebook, the final feature mentioned `emissions_nox_mgkm` will be used as the _endogenous_ / _dependent_ / _output_ feature of the linear models. All other features may be used as _exogenous_ / _independent_ / _input_ features."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Create design matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'emissions_nox_mgkm ~ metric_combined + engine_capacity + emissions_co_mgkm + mfr_is_vw + trans + fuel_type + is_tdi'"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fml_all = '{} ~ '.format(ft_endog) + ' + '.join(fts_num_smp + fts_cat_smp)\n",
    "fml_all"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2593, 9)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>1731</th>\n",
       "      <th>1033</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",
       "      <th>dtype</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Intercept</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>2593</td>\n",
       "      <td>1.00</td>\n",
       "      <td>0.00</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mfr_is_vw[T.True]</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>2593</td>\n",
       "      <td>0.04</td>\n",
       "      <td>0.20</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>trans[T.manual]</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>2593</td>\n",
       "      <td>0.49</td>\n",
       "      <td>0.50</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>trans[T.semiauto]</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>2593</td>\n",
       "      <td>0.12</td>\n",
       "      <td>0.32</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fuel_type[T.petrol]</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>2593</td>\n",
       "      <td>0.50</td>\n",
       "      <td>0.50</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>is_tdi[T.True]</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>2593</td>\n",
       "      <td>0.12</td>\n",
       "      <td>0.32</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>metric_combined</th>\n",
       "      <td>-0.312750</td>\n",
       "      <td>-0.472703</td>\n",
       "      <td>2593</td>\n",
       "      <td>-0.00</td>\n",
       "      <td>0.50</td>\n",
       "      <td>-0.685973</td>\n",
       "      <td>-0.339409</td>\n",
       "      <td>-0.152797</td>\n",
       "      <td>0.167108</td>\n",
       "      <td>2.753013</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>engine_capacity</th>\n",
       "      <td>0.042048</td>\n",
       "      <td>-0.278857</td>\n",
       "      <td>2593</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.50</td>\n",
       "      <td>-0.920116</td>\n",
       "      <td>-0.278857</td>\n",
       "      <td>-0.045471</td>\n",
       "      <td>-0.037765</td>\n",
       "      <td>2.579014</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>emissions_co_mgkm</th>\n",
       "      <td>0.099293</td>\n",
       "      <td>-0.421449</td>\n",
       "      <td>2593</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.50</td>\n",
       "      <td>-0.846545</td>\n",
       "      <td>-0.382482</td>\n",
       "      <td>-0.088458</td>\n",
       "      <td>0.272874</td>\n",
       "      <td>2.281452</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                         1731      1033  count  mean   std       min  \\\n",
       "Intercept            1.000000  1.000000   2593  1.00  0.00  1.000000   \n",
       "mfr_is_vw[T.True]    0.000000  0.000000   2593  0.04  0.20  0.000000   \n",
       "trans[T.manual]      0.000000  1.000000   2593  0.49  0.50  0.000000   \n",
       "trans[T.semiauto]    0.000000  0.000000   2593  0.12  0.32  0.000000   \n",
       "fuel_type[T.petrol]  0.000000  0.000000   2593  0.50  0.50  0.000000   \n",
       "is_tdi[T.True]       0.000000  0.000000   2593  0.12  0.32  0.000000   \n",
       "metric_combined     -0.312750 -0.472703   2593 -0.00  0.50 -0.685973   \n",
       "engine_capacity      0.042048 -0.278857   2593  0.00  0.50 -0.920116   \n",
       "emissions_co_mgkm    0.099293 -0.421449   2593  0.00  0.50 -0.846545   \n",
       "\n",
       "                          25%       50%       75%       max    dtype  \n",
       "Intercept            1.000000  1.000000  1.000000  1.000000  float64  \n",
       "mfr_is_vw[T.True]    0.000000  0.000000  0.000000  1.000000  float64  \n",
       "trans[T.manual]      0.000000  0.000000  1.000000  1.000000  float64  \n",
       "trans[T.semiauto]    0.000000  0.000000  0.000000  1.000000  float64  \n",
       "fuel_type[T.petrol]  0.000000  1.000000  1.000000  1.000000  float64  \n",
       "is_tdi[T.True]       0.000000  0.000000  0.000000  1.000000  float64  \n",
       "metric_combined     -0.339409 -0.152797  0.167108  2.753013  float64  \n",
       "engine_capacity     -0.278857 -0.045471 -0.037765  2.579014  float64  \n",
       "emissions_co_mgkm   -0.382482 -0.088458  0.272874  2.281452  float64  "
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(mx_en, mx_ex) = pt.dmatrices(fml_all, dfs, return_type='dataframe', NA_action='raise')\n",
    "custom_describe(mx_ex, 2, )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reload Models and Traces"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reload PyMC3 model and traces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loaded object from disk at models/mdl_robust_pymc.pkl\n"
     ]
    }
   ],
   "source": [
    "with pm.Model():\n",
    "    mdl_robust_pymc = read_pickle(relnm='models/mdl_robust_pymc.pkl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## load traces within model context\n",
    "with mdl_robust_pymc:\n",
    "    trc_robust_pymc = pm.backends.text.load('traces/trc_robust_pymc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# alternatively, load with a ref to reloaded model\n",
    "trc_robust_pymc = pm.backends.text.load('traces/trc_robust_pymc', model=mdl_robust_pymc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reload PyStan fitted model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Fit object contains traces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loaded object from disk at models/mdl_robust_stan_model.pkl\n",
      "Loaded object from disk at models/mdl_robust_stan_fit.pkl\n"
     ]
    }
   ],
   "source": [
    "## load the model first, and then the fit:\n",
    "\n",
    "mdl_robust_pystan_model = read_pickle(relnm='models/mdl_robust_stan_model.pkl')\n",
    "\n",
    "mdl_robust_pystan_fit = read_pickle(relnm='models/mdl_robust_stan_fit.pkl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# More model comparison techniques available in PyMC3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Lets create another pymc3 model which ought to yield worse results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with pm.Model() as mdl_ols_pymc:\n",
    "      \n",
    "    pm.glm.glm(fml_all,\n",
    "               dfs,\n",
    "               intercept_prior=pm.Uniform.dist(lower=-1e3, upper=1e3),\n",
    "               regressor_prior=pm.Uniform.dist(lower=-1e3, upper=1e3),\n",
    "               family=pm.glm.families.Normal())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with mdl_robust_pymc:\n",
    "\n",
    "    ## find MAP using Powell, seems to be more robust\n",
    "    start_MAP = pm.find_MAP(fmin=optimize.fmin_powell)\n",
    "\n",
    "    ## take samples\n",
    "    trc_ols_pymc = pm.sample(2000, start=start_MAP, njobs=1, step=pm.NUTS()\n",
    "                        ,trace=pm.backends.Text('traces/trc_ols_pymc'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Deviance Information Criterion (DIC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "datasets by using the Deviance Information Criterion (DIC) and Watanabe - Akaike (or Widest Available) Information Criterion (WAIC).\n",
    "+ DIC (`stats.dic`) and WAIC (`stats.waic`) are new additions to PyMC3, so this example shows their usage in a more concrete fashion, also usage of the new `glm` submodule.\n",
    "+ The example was inspired by Jake Vanderplas' [recent blogpost](https://jakevdp.github.io/blog/2015/08/07/frequentism-and-bayesianism-5-model-selection/) on model selection, although in this first iteration, Cross-Validation and Bayes Factor comparison are not implemented.\n",
    "+ The datasets are tiny and generated within this Notebook. They contain errors in the measured value (y) only.\n",
    "\n",
    "\n",
    "For more information on Model Selection in PyMC3, and about DIC and WAIC, you could start with:\n",
    "\n",
    "+ Thomas Wiecki's [detailed response](https://stats.stackexchange.com/questions/161082/bayesian-model-selection-in-pymc3/166383#166383) to a question on Cross Validated\n",
    "+ The Deviance Information Criterion: 12 Years On [(Speigelhalter et al 2014)](http://onlinelibrary.wiley.com/doi/10.1111/rssb.12062/abstract)\n",
    "+ A Widely Applicable Bayesian Information Criterion [(Watanabe 2013)](http://www.jmlr.org/papers/volume14/watanabe13a/watanabe13a.pdf)\n",
    "+ Efficient Implementation of Leave-One-Out Cross-Validation and WAIC for Evaluating Fitted Bayesian Models [(Gelman et al 2015)](http://arxiv.org/abs/1507.04544)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "The [Deviance Information Criterion](https://en.wikipedia.org/wiki/Deviance_information_criterion) (DIC) is a fairly unsophisticated method for comparing the deviance of likelihood across the sample traces of a model run. \n",
    "\n",
    "However, this simplicity apparently yields quite good results in a variety of cases, see the discussion worth reading in ([Speigelhalter et al 2014](http://onlinelibrary.wiley.com/doi/10.1111/rssb.12062/abstract))\n",
    "\n",
    "DIC has recently been added to PyMC3, so lets see what it tells us about our model fits. Lower numbers are better"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Deviance Information Criterion (DIC) is a fairly unsophisticated method for comparing the deviance of likelhood across the the sample traces of a model run. However, this simplicity apparently yields quite good results in a variety of cases, see the discussion worth reading in [(Speigelhalter et al 2014)](http://onlinelibrary.wiley.com/doi/10.1111/rssb.12062/abstract)\n",
    "\n",
    "DIC has recently been added to PyMC3, so lets see what it tells us about our model fits for both datasets."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Manual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dftrc_lin = pm.trace_to_dataframe(traces_lin['k1'])\n",
    "trc_lin_logp = dftrc_lin.apply(lambda x: models_lin['k1'].logp(x.to_dict()), axis=1)\n",
    "mean_deviance = -2 * trc_lin_logp.mean(0)\n",
    "mean_deviance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "deviance_at_mean = -2 * models_lin['k1'].logp(dftrc_lin.mean(0).to_dict())\n",
    "deviance_at_mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dic_k1 = 2 * mean_deviance - deviance_at_mean\n",
    "dic_k1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Or use stats.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pm.stats.dic(model=models_lin['k1'], trace=traces_lin['k1'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Observe**\n",
    "\n",
    "+ We should prefer the model(s) with lower DIC, which (happily) directly opposes the increasing likelihood we see above.\n",
    "\n",
    "\n",
    "+ Linear-generated data (lhs):\n",
    "    + The DIC increases monotonically with model complexity, this is great too see!\n",
    "    + The more complicated the model, the more it would appear we are overfitting.\n",
    "\n",
    "\n",
    "+ Quadratic-generated data (rhs):\n",
    "    + The DIC dips slightly for the correct model k2\n",
    "    + The difference is slight though!\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# OLS\n",
    "pm.stats.dic(model=mdl_ols, trace=trc_ols[-1000:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Lasso\n",
    "pm.stats.dic(model=mdl_lasso, trace=trc_lasso[-7000::7])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Ridge\n",
    "pm.stats.dic(model=mdl_ridge, trace=trc_ridge[-1000:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Observe:**\n",
    "\n",
    "+ All 3 values are quite similar! \n",
    "+ Interestingly, the Lasso has the highest (worst) DIC value, probably because I deliberately chose a suboptimal regularization parameter $\\lambda$ in order to favour the dropping of feature coefficients. This likely resulted in an _underfitted_ model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dics = [pm.stats.dic(model=mdl, trace=trc) for mdl, trc in zip(mdls, trcs)]\n",
    "dfdics = pd.DataFrame({'model':['ols','lasso','ridge','ridge_hn','student'],'dic':dics})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "waics = [pm.stats.waic(model=mdl, trace=trc) for mdl, trc in zip(mdls, trcs)]\n",
    "dfwaics = pd.DataFrame({'model':['ols','lasso','ridge','ridge_hn','student'],'waic':waics})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "g = sns.barplot(y='model',x='waic', data=dfwaics, orient='h')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "See my example in pymc3 docs: https://github.com/pymc-devs/pymc3/blob/master/pymc3/examples/GLM-model-selection.ipynb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Compare Watanabe - Akaike Information Criterion [WAIC]\n",
    "\n",
    "The Widely Applicable Bayesian Information Criterion (WBIC), a.k.a the Watanabe - Akaike Information Criterion (WAIC) is another simple option for calculating the goodness-of-fit of amodel using numerical techniques. See [(Watanabe 2013)](http://www.jmlr.org/papers/volume14/watanabe13a/watanabe13a.pdf) for details.\n",
    "\n",
    "WAIC has also recently been added to PyMC3, so lets see what it tells us about our model fits for both datasets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "##### This time go straight for the implementation in pymc3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pm.stats.waic(model=models_lin['k1'], trace=traces_lin['k1'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "**Observe**\n",
    "\n",
    "+ We should prefer the model(s) with lower WAIC\n",
    "\n",
    "\n",
    "+ Linear-generated data (lhs):\n",
    "    + The WAIC seems quite flat across models\n",
    "    + The WAIC seems best (lowest) for simpler models, but **k1** doesn't stand out as much as it did when using DIC\n",
    "\n",
    "\n",
    "+ Quadratic-generated data (rhs):\n",
    "    + The WAIC is certainly wrong for **k1**, but otherwise also quite flat across the models\n",
    "    + There does appear to be a slight dip in the right place at **k2**\n",
    "    \n",
    "    \n",
    "For these particular models and data, I would prefer to use the DIC scores in order to choose models.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hold out set (Theano shared vars)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "##### Create split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "kfold_shuffle = ShuffleSplit(mx_ex.shape[0], n_iter=1, test_size=0.2, random_state=0)\n",
    "\n",
    "for train_idx, test_idx in kfold_shuffle:\n",
    "\n",
    "    mx_ex_train = mx_ex.iloc[train_idx]\n",
    "    mx_en_train = mx_en.iloc[train_idx]\n",
    "    mx_ex_test = mx_ex.iloc[test_idx]\n",
    "    mx_en_test = mx_en.iloc[test_idx]\n",
    "\n",
    "print('ntrain: {}, ntest: {}'.format(mx_ex_train.shape, mx_ex_test.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Use theano shared variables to allow us to switch train for test and re-use the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "shd_b0 = thno.shared(mx_ex_train['Intercept'].values)\n",
    "shd_b1 = thno.shared(mx_ex_train['is_tdi[T.True]'].values)\n",
    "shd_b2 = thno.shared(mx_ex_train['engine_capacity'].values)\n",
    "shd_b3 = thno.shared(mx_ex_train['metric_combined'].values)\n",
    "\n",
    "shd_obs = thno.shared(mx_en_train['emissions_nox_mgkm'].values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Train model "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "with pm.Model() as mdl_holdout:\n",
    "\n",
    "    # define priors, use Normal for Ridge (sd=100, weakly informative)\n",
    "    b0 = pm.Normal('b0_intercept', mu=0, sd=100)\n",
    "    b1 = pm.Normal('b1_is_tdi', mu=0, sd=100)\n",
    "    b2 = pm.Normal('b2_engine_capacity', mu=0, sd=100)\n",
    "    b3 = pm.Normal('b3_metric_combined', mu=0, sd=100)    \n",
    " \n",
    "    # define linear model\n",
    "    yest = (b0 * shd_b0 +\n",
    "            b1 * shd_b1 +\n",
    "            b2 * shd_b2 + \n",
    "            b3 * shd_b3)\n",
    "\n",
    "    # create MLE with stdev chosen from HalfNormal dist and free param on nu             \n",
    "    epsilon = pm.HalfNormal('epsilon', sd=10)\n",
    "    nu = pm.DiscreteUniform('nu', lower=1, upper=100)\n",
    "   \n",
    "    ## Student T likelihood with variable degress of freedom nu\n",
    "    likelihood = pm.StudentT('likelihood', nu=nu, mu=yest, sd=epsilon\n",
    "                            ,observed=shd_obs)\n",
    "    \n",
    "    start_MAP = pm.find_MAP(fmin=optimize.fmin_powell, disp=True)\n",
    "      \n",
    "    trc_holdout = pm.sample(2000, step=pm.NUTS(), start=start_MAP, progressbar=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "### Test model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "##### Now switch values in the shared variables to the test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "shd_b0.set_value(mx_ex_test['Intercept'].values)\n",
    "shd_b1.set_value(mx_ex_test['is_tdi[T.True]'].values)\n",
    "shd_b2.set_value(mx_ex_test['engine_capacity'].values)\n",
    "shd_b3.set_value(mx_ex_test['metric_combined'].values)\n",
    "\n",
    "shd_obs.set_value(mx_en_test['emissions_nox_mgkm'].values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "##### And run the posterior check which will now create estimates for the test data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ppc = pm.sample_ppc(trc_holdout[-1000:], samples=500, model=mdl_holdout, size=50)\n",
    "ppc['likelihood'].shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluate pediction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PPC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ax = plt.subplot()\n",
    "sns.distplot([n.mean() for n in ppc['likelihood']], kde=False, ax=ax)\n",
    "ax.axvline(mx_en_test[ft_endog].mean())\n",
    "# ax.set(title='Posterior predictive of the mean', xlabel='mean(x)', ylabel='Frequency');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Coverage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MSE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## KS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Leave One Out LOO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## K-Fold Cross Validation and/or Leave-One-Out (LOO)\n",
    "\n",
    "http://www.stat.columbia.edu/~gelman/research/unpublished/waic_stan.pdf\n",
    "\n",
    "\n",
    "\n",
    "http://arxiv.org/abs/1507.04544\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Bayes Factor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Following text lifted directly from [JakeVDP blogpost](https://jakevdp.github.io/blog/2015/08/07/frequentism-and-bayesianism-5-model-selection/)\n",
    "\n",
    "The Bayesian approach proceeds very differently. Recall that the Bayesian model involves computing the odds ratio between two models:\n",
    "\n",
    "$$O_{21}=\\frac{P(M_{2} \\;|\\; D)}{P(M_{1} \\;|\\; D)}=\\frac{P(D \\;|\\; M_{2})}{P(D \\;|\\; M_{1})}\\frac{P(M_{2})}{P(M_{1})}$$\n",
    "\n",
    "Here the ratio $\\frac{P(M2)}{P(M1)}$ is the prior odds ratio, and is often assumed to be equal to 1 if no compelling prior evidence favors one model over another. The ratio $\\frac{P(D \\;|\\; M2)}{P(D \\;|\\; M1)}$ is the **Bayes factor**, and is the key to Bayesian model selection.\n",
    "\n",
    "\n",
    "The Bayes factor can be computed by evaluating the integral over the parameter likelihood:\n",
    "\n",
    "$$P(D \\;|\\; M)=\\int_{\\Omega}P(D \\;|\\; \\theta,M) \\; P(\\theta \\;|\\; M) \\;d\\theta$$\n",
    "\n",
    "This integral is over the entire parameter space of the model, and thus can be extremely computationally intensive, especially as the dimension of the model grows beyond a few. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "**&copy; Applied AI Ltd 2016**  \n",
    "<a href='http://www.applied.ai'>applied.ai</a>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}