L1,L2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# L1, L2 Regularization: Frequentist and Bayesian "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Lasso, also called the LASSO (least absolute shrinkage and selection operator) is a somewhat magical technique (to me) for constrained aka regularized regression, which adds an additional parameter to the cost function of a least squares regression. This additional parameter is a multiple of the sum of absolute values of the feature coefficients: the L1 norm.\n",
    "\n",
    "The effect is twofold:\n",
    "\n",
    "- to force coefficient values smaller, thus controlling overfitting and making the model more generalizable\n",
    "- to force weak or small coefficient values to zero, thus performing feature selection for us: the remaining non-zero coefficients are more important to the estimating power of the model\n",
    "\n",
    "(On whiteboard)\n",
    "Set parameter  λ  to constrain the sum of absolute values of  β"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr", 
     "output_type": "stream",
     "text": [
      "/anaconda3/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": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import statsmodels.api as sm\n",
    "from sklearn.model_selection import train_test_split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(\"yield_forecast.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def format_data(df,test_size):\n",
    "   \n",
    "    # Targets are final grade of student\n",
    "    labels = df['y_pred']\n",
    "    # Drop the school and the grades from features\n",
    "    df = df.drop(columns=['asd_desc','state','y_pred'],axis = 1)\n",
    "    \n",
    "    # One-Hot Encoding of Categorical Variables\n",
    "    #df = pd.get_dummies(df)\n",
    "    \n",
    "    #df['y'] = list(labels)\n",
    "    #most_correlated = df.corr().abs()['y'].sort_values(ascending=False)\n",
    "    #print(most_correlated)\n",
    "    \n",
    "    # Keep correlations greater than 0.2 in absolute value\n",
    "    #most_correlated = most_correlated[most_correlated >= 0.2][1:]\n",
    "    \n",
    "    #df = df.ix[:, most_correlated.index]\n",
    "    #df = df.drop(columns = 'y')\n",
    "    \n",
    "    # Split into training/testing sets with 25% split\n",
    "    X_train, X_test, y_train, y_test = train_test_split(df, labels, \n",
    "                                                        test_size = test_size,\n",
    "                                                       random_state=42)\n",
    "    \n",
    "    return X_train, X_test, y_train, y_test\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "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>cyield</th>\n",
       "      <th>irig_flag</th>\n",
       "      <th>days_under0</th>\n",
       "      <th>dewPoint</th>\n",
       "      <th>precipAccumulation</th>\n",
       "      <th>precip</th>\n",
       "      <th>days_under_n10</th>\n",
       "      <th>days_over42</th>\n",
       "      <th>days_over32</th>\n",
       "      <th>humidity</th>\n",
       "      <th>temp_delta</th>\n",
       "      <th>temperatureMin</th>\n",
       "      <th>apparentTemperatureMin</th>\n",
       "      <th>precipIntensity</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>167</th>\n",
       "      <td>43.7</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>38.253084</td>\n",
       "      <td>5.011</td>\n",
       "      <td>122.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>192.0</td>\n",
       "      <td>202.0</td>\n",
       "      <td>0.593458</td>\n",
       "      <td>22.783785</td>\n",
       "      <td>42.809065</td>\n",
       "      <td>39.582570</td>\n",
       "      <td>0.002351</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>90</th>\n",
       "      <td>126.0</td>\n",
       "      <td>1</td>\n",
       "      <td>0.0</td>\n",
       "      <td>36.287056</td>\n",
       "      <td>0.817</td>\n",
       "      <td>103.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>200.0</td>\n",
       "      <td>214.0</td>\n",
       "      <td>0.482383</td>\n",
       "      <td>24.715794</td>\n",
       "      <td>46.114159</td>\n",
       "      <td>43.728738</td>\n",
       "      <td>0.000774</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>145.0</td>\n",
       "      <td>1</td>\n",
       "      <td>0.0</td>\n",
       "      <td>34.378551</td>\n",
       "      <td>1.254</td>\n",
       "      <td>137.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>207.0</td>\n",
       "      <td>213.0</td>\n",
       "      <td>0.531636</td>\n",
       "      <td>26.164019</td>\n",
       "      <td>40.776402</td>\n",
       "      <td>38.332290</td>\n",
       "      <td>0.001316</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>128</th>\n",
       "      <td>105.8</td>\n",
       "      <td>1</td>\n",
       "      <td>0.0</td>\n",
       "      <td>33.207944</td>\n",
       "      <td>7.027</td>\n",
       "      <td>142.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>195.0</td>\n",
       "      <td>214.0</td>\n",
       "      <td>0.553925</td>\n",
       "      <td>26.733879</td>\n",
       "      <td>38.655467</td>\n",
       "      <td>35.487757</td>\n",
       "      <td>0.001543</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>118</th>\n",
       "      <td>66.7</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>35.930607</td>\n",
       "      <td>0.458</td>\n",
       "      <td>100.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>195.0</td>\n",
       "      <td>211.0</td>\n",
       "      <td>0.491963</td>\n",
       "      <td>24.069112</td>\n",
       "      <td>45.568785</td>\n",
       "      <td>43.678551</td>\n",
       "      <td>0.001033</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     cyield  irig_flag  days_under0   dewPoint  precipAccumulation  precip  \\\n",
       "167    43.7          0          0.0  38.253084               5.011   122.0   \n",
       "90    126.0          1          0.0  36.287056               0.817   103.0   \n",
       "5     145.0          1          0.0  34.378551               1.254   137.0   \n",
       "128   105.8          1          0.0  33.207944               7.027   142.0   \n",
       "118    66.7          0          0.0  35.930607               0.458   100.0   \n",
       "\n",
       "     days_under_n10  days_over42  days_over32  humidity  temp_delta  \\\n",
       "167             0.0        192.0        202.0  0.593458   22.783785   \n",
       "90              0.0        200.0        214.0  0.482383   24.715794   \n",
       "5               0.0        207.0        213.0  0.531636   26.164019   \n",
       "128             0.0        195.0        214.0  0.553925   26.733879   \n",
       "118             0.0        195.0        211.0  0.491963   24.069112   \n",
       "\n",
       "     temperatureMin  apparentTemperatureMin  precipIntensity  \n",
       "167       42.809065               39.582570         0.002351  \n",
       "90        46.114159               43.728738         0.000774  \n",
       "5         40.776402               38.332290         0.001316  \n",
       "128       38.655467               35.487757         0.001543  \n",
       "118       45.568785               43.678551         0.001033  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_train_yield, X_test_yield, y_train_yield, y_test_yield = format_data(df,0.25)\n",
    "X_train_yield.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>         <td>y_pred</td>      <th>  R-squared:         </th> <td>   0.995</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.994</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   2127.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Tue, 15 May 2018</td> <th>  Prob (F-statistic):</th> <td>6.37e-152</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>07:04:38</td>     <th>  Log-Likelihood:    </th> <td> -497.25</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>   152</td>      <th>  AIC:               </th> <td>   1019.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>   140</td>      <th>  BIC:               </th> <td>   1055.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>    12</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>cyield</th>                 <td>    0.5085</td> <td>    0.037</td> <td>   13.885</td> <td> 0.000</td> <td>    0.436</td> <td>    0.581</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>irig_flag</th>              <td>   21.9086</td> <td>    2.135</td> <td>   10.262</td> <td> 0.000</td> <td>   17.688</td> <td>   26.129</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>days_under0</th>            <td>  8.65e-13</td> <td> 1.06e-11</td> <td>    0.082</td> <td> 0.935</td> <td>   -2e-11</td> <td> 2.18e-11</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>dewPoint</th>               <td>   -3.3828</td> <td>    0.776</td> <td>   -4.361</td> <td> 0.000</td> <td>   -4.916</td> <td>   -1.849</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>precipAccumulation</th>     <td>    0.0498</td> <td>    0.080</td> <td>    0.624</td> <td> 0.533</td> <td>   -0.108</td> <td>    0.207</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>precip</th>                 <td>   -0.1072</td> <td>    0.035</td> <td>   -3.041</td> <td> 0.003</td> <td>   -0.177</td> <td>   -0.038</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>days_under_n10</th>         <td>-7.228e-13</td> <td> 1.15e-12</td> <td>   -0.630</td> <td> 0.530</td> <td>-2.99e-12</td> <td> 1.55e-12</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>days_over42</th>            <td>    0.3203</td> <td>    0.103</td> <td>    3.122</td> <td> 0.002</td> <td>    0.117</td> <td>    0.523</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>days_over32</th>            <td>   -0.2728</td> <td>    0.102</td> <td>   -2.675</td> <td> 0.008</td> <td>   -0.474</td> <td>   -0.071</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>humidity</th>               <td>  114.2470</td> <td>   30.994</td> <td>    3.686</td> <td> 0.000</td> <td>   52.970</td> <td>  175.524</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>temp_delta</th>             <td>    0.5845</td> <td>    0.247</td> <td>    2.364</td> <td> 0.019</td> <td>    0.096</td> <td>    1.073</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>temperatureMin</th>         <td>    2.8752</td> <td>    0.899</td> <td>    3.198</td> <td> 0.002</td> <td>    1.098</td> <td>    4.653</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>apparentTemperatureMin</th> <td>   -1.0454</td> <td>    0.932</td> <td>   -1.122</td> <td> 0.264</td> <td>   -2.888</td> <td>    0.797</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>precipIntensity</th>        <td>  -89.3043</td> <td>  550.564</td> <td>   -0.162</td> <td> 0.871</td> <td>-1177.798</td> <td>  999.190</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 0.743</td> <th>  Durbin-Watson:     </th> <td>   2.190</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.690</td> <th>  Jarque-Bera (JB):  </th> <td>   0.867</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 0.144</td> <th>  Prob(JB):          </th> <td>   0.648</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 2.767</td> <th>  Cond. No.          </th> <td>1.07e+18</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                 y_pred   R-squared:                       0.995\n",
       "Model:                            OLS   Adj. R-squared:                  0.994\n",
       "Method:                 Least Squares   F-statistic:                     2127.\n",
       "Date:                Tue, 15 May 2018   Prob (F-statistic):          6.37e-152\n",
       "Time:                        07:04:38   Log-Likelihood:                -497.25\n",
       "No. Observations:                 152   AIC:                             1019.\n",
       "Df Residuals:                     140   BIC:                             1055.\n",
       "Df Model:                          12                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==========================================================================================\n",
       "                             coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------------------\n",
       "cyield                     0.5085      0.037     13.885      0.000       0.436       0.581\n",
       "irig_flag                 21.9086      2.135     10.262      0.000      17.688      26.129\n",
       "days_under0              8.65e-13   1.06e-11      0.082      0.935      -2e-11    2.18e-11\n",
       "dewPoint                  -3.3828      0.776     -4.361      0.000      -4.916      -1.849\n",
       "precipAccumulation         0.0498      0.080      0.624      0.533      -0.108       0.207\n",
       "precip                    -0.1072      0.035     -3.041      0.003      -0.177      -0.038\n",
       "days_under_n10         -7.228e-13   1.15e-12     -0.630      0.530   -2.99e-12    1.55e-12\n",
       "days_over42                0.3203      0.103      3.122      0.002       0.117       0.523\n",
       "days_over32               -0.2728      0.102     -2.675      0.008      -0.474      -0.071\n",
       "humidity                 114.2470     30.994      3.686      0.000      52.970     175.524\n",
       "temp_delta                 0.5845      0.247      2.364      0.019       0.096       1.073\n",
       "temperatureMin             2.8752      0.899      3.198      0.002       1.098       4.653\n",
       "apparentTemperatureMin    -1.0454      0.932     -1.122      0.264      -2.888       0.797\n",
       "precipIntensity          -89.3043    550.564     -0.162      0.871   -1177.798     999.190\n",
       "==============================================================================\n",
       "Omnibus:                        0.743   Durbin-Watson:                   2.190\n",
       "Prob(Omnibus):                  0.690   Jarque-Bera (JB):                0.867\n",
       "Skew:                           0.144   Prob(JB):                        0.648\n",
       "Kurtosis:                       2.767   Cond. No.                     1.07e+18\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "[2] The smallest eigenvalue is 1.41e-29. This might indicate that there are\n",
       "strong multicollinearity problems or that the design matrix is singular.\n",
       "\"\"\""
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "smfit = sm.OLS(y_train_yield, X_train_yield).fit()\n",
    "smfit.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "predict_ols  = smfit.predict(X_test_yield)\n",
    "rmse_ols = sqrt(mean_squared_error(predict_ols, y_test_yield))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8.424805020095407"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rmse_ols"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import preprocessing\n",
    "\n",
    "continous_var_list = list(X_train_yield.columns[2:])\n",
    "\n",
    "continous_var_list.append('cyield')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# scale(center) training data \n",
    "for i in continous_var_list:\n",
    "    X_train_yield[i]=preprocessing.scale(X_train_yield[i].astype('float64'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import Lasso\n",
    "\n",
    "#sklearn uses the term alpha but it is the same in effect as lambda\n",
    "lasso = Lasso(alpha = 1)\n",
    "ols_lasso = lasso.fit(X_train_yield,y_train_yield)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Index(['cyield', 'irig_flag', 'days_under0', 'dewPoint', 'precipAccumulation',\n",
      "       'precip', 'days_under_n10', 'days_over42', 'days_over32', 'humidity',\n",
      "       'temp_delta', 'temperatureMin', 'apparentTemperatureMin',\n",
      "       'precipIntensity'],\n",
      "      dtype='object')\n"
     ]
    }
   ],
   "source": [
    "print(X_train_yield.columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([21.22362034, 10.62431805,  0.        , -1.90546923, -0.        ,\n",
       "       -0.25399776,  0.        , -0.        , -0.        , -0.66028305,\n",
       "        0.        , -0.        , -0.        ,  0.        ])"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ols_lasso.coef_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import mean_squared_error\n",
    "from math import sqrt\n",
    "predict_lasso  = ols_lasso.predict(X_test_yield)\n",
    "rmse_lasso = sqrt(mean_squared_error(predict_lasso, y_test_yield))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1864.3329546588554"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rmse_lasso"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Bayesian Lasso Regression\n",
    "\n",
    "In the Bayesian framework, the choice of regulariser is analogous to the choice of prior over the weights. If a Gaussian prior is used, then the Maximum a Posteriori (MAP) solution will be the same as if an L2 penalty was used.  We will see this in the Ridge section. \n",
    "\n",
    "It is possible to imitate the effect of the L1 regulariation in Bayesian linear regression by using a Laplacian prior for feature coefficients.In fact, when you place a Laplace prior over the parameters, the MAP solution should be identical (not merely similar) to regularization with the L1 penalty and the Laplace prior will produce an identical shrinkage effect to the L1 penalty.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What does a laplace distribution Look like?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"laplace_dist.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see above, the Laplace distribution or 'double exponential distribution' looks like two exponential distributions back to back, (-inf, inf)\n",
    "- The distribution is symmetric and centered on a location, by default at zero\n",
    "- The scale parameter b acts similarly to a variance parameter, making the distribution narrower &  taller when small, and wider & shallower when large.\n",
    "- The tails are heavier than a normal distribution\n",
    "\n",
    "- The effect of using this for a prior is to encourage the marginal distribution of the feature coefficients to be close to zero. 'Strong' coefficients will move from the zero point because of the effect of the data likelihood, but 'weak' coefficients will be overly influenced by the prior distribution and remain close to zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the varaince is forced to be tighter, there is a penlization of the sum of the absolute beta-values. Let's try it out, using a narrow Laplace distribution for each coefficient;  b=0.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda3/lib/python3.6/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  from ._conv import register_converters as _register_converters\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "y ~ cyield + irig_flag + days_under0 + dewPoint + precipAccumulation + precip + days_under_n10 + days_over42 + days_over32 + humidity + temp_delta + temperatureMin + apparentTemperatureMin + precipIntensity\n"
     ]
    }
   ],
   "source": [
    "\n",
    "import pymc3 as pm\n",
    "\n",
    "\n",
    "# We need to add our target back into the dataframe \n",
    "X_train_yield['y'] = list(y_train_yield)\n",
    "formula = 'y ~ ' + ' + '.join(['%s' % variable for variable in X_train_yield.columns[:-1]])\n",
    "print(formula)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "CompoundStep\n",
      ">Metropolis: [sd_log__]\n",
      ">Metropolis: [precipIntensity]\n",
      ">Metropolis: [apparentTemperatureMin]\n",
      ">Metropolis: [temperatureMin]\n",
      ">Metropolis: [temp_delta]\n",
      ">Metropolis: [humidity]\n",
      ">Metropolis: [days_over32]\n",
      ">Metropolis: [days_over42]\n",
      ">Metropolis: [days_under_n10]\n",
      ">Metropolis: [precip]\n",
      ">Metropolis: [precipAccumulation]\n",
      ">Metropolis: [dewPoint]\n",
      ">Metropolis: [days_under0]\n",
      ">Metropolis: [irig_flag]\n",
      ">Metropolis: [cyield]\n",
      ">Metropolis: [Intercept]\n",
      "100%|██████████| 40500/40500 [03:45<00:00, 179.44it/s]\n",
      "The number of effective samples is smaller than 10% for some parameters.\n"
     ]
    }
   ],
   "source": [
    "\n",
    "with pm.Model() as mdl_lasso:\n",
    "\n",
    "    ## Use GLM submodule for simplified model specification\n",
    "    ## Betas are Laplace (for Lasso)\n",
    "    ## Likelihood is Normal \n",
    "    \n",
    "    priors = {\"Intercept\": pm.Laplace.dist(mu=0,b=0.1),\n",
    "          \"Regressor\": pm.Laplace.dist(mu=0, b=0.1)\n",
    "             }\n",
    "    pm.GLM.from_formula(formula, X_train_yield\n",
    "               ,family=pm.glm.families.Normal(),\n",
    "               priors = priors)\n",
    "   \n",
    "    ## take samples using Metropolis\n",
    "    trc_lasso = pm.sample(40000, step=pm.Metropolis())\n",
    "    \n",
    "rvs_lasso = [rv.name for rv in mdl_lasso.free_RVs]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "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>sd</th>\n",
       "      <th>mc_error</th>\n",
       "      <th>hpd_2.5</th>\n",
       "      <th>hpd_97.5</th>\n",
       "      <th>n_eff</th>\n",
       "      <th>Rhat</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Intercept</th>\n",
       "      <td>0.033369</td>\n",
       "      <td>0.147348</td>\n",
       "      <td>0.001307</td>\n",
       "      <td>-0.243261</td>\n",
       "      <td>0.372392</td>\n",
       "      <td>12752.112726</td>\n",
       "      <td>1.000212</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>cyield</th>\n",
       "      <td>0.008753</td>\n",
       "      <td>0.144545</td>\n",
       "      <td>0.001244</td>\n",
       "      <td>-0.286242</td>\n",
       "      <td>0.317488</td>\n",
       "      <td>12215.797341</td>\n",
       "      <td>0.999988</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>irig_flag</th>\n",
       "      <td>0.023684</td>\n",
       "      <td>0.144673</td>\n",
       "      <td>0.001200</td>\n",
       "      <td>-0.261596</td>\n",
       "      <td>0.340129</td>\n",
       "      <td>13582.094278</td>\n",
       "      <td>0.999994</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>days_under0</th>\n",
       "      <td>0.001703</td>\n",
       "      <td>0.140275</td>\n",
       "      <td>0.001228</td>\n",
       "      <td>-0.296166</td>\n",
       "      <td>0.293691</td>\n",
       "      <td>12795.671660</td>\n",
       "      <td>1.000011</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>dewPoint</th>\n",
       "      <td>-0.005811</td>\n",
       "      <td>0.141128</td>\n",
       "      <td>0.001142</td>\n",
       "      <td>-0.286644</td>\n",
       "      <td>0.302634</td>\n",
       "      <td>14326.564104</td>\n",
       "      <td>0.999988</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>precipAccumulation</th>\n",
       "      <td>-0.002344</td>\n",
       "      <td>0.139962</td>\n",
       "      <td>0.001232</td>\n",
       "      <td>-0.285032</td>\n",
       "      <td>0.303106</td>\n",
       "      <td>13599.872401</td>\n",
       "      <td>1.000029</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>precip</th>\n",
       "      <td>-0.000225</td>\n",
       "      <td>0.143678</td>\n",
       "      <td>0.001154</td>\n",
       "      <td>-0.290755</td>\n",
       "      <td>0.321854</td>\n",
       "      <td>12685.048499</td>\n",
       "      <td>1.000184</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>days_under_n10</th>\n",
       "      <td>-0.000827</td>\n",
       "      <td>0.143340</td>\n",
       "      <td>0.001403</td>\n",
       "      <td>-0.308703</td>\n",
       "      <td>0.296381</td>\n",
       "      <td>12484.341763</td>\n",
       "      <td>1.000030</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>days_over42</th>\n",
       "      <td>0.002368</td>\n",
       "      <td>0.144378</td>\n",
       "      <td>0.001394</td>\n",
       "      <td>-0.297255</td>\n",
       "      <td>0.311021</td>\n",
       "      <td>12704.563450</td>\n",
       "      <td>0.999994</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>days_over32</th>\n",
       "      <td>0.004131</td>\n",
       "      <td>0.140339</td>\n",
       "      <td>0.001265</td>\n",
       "      <td>-0.285340</td>\n",
       "      <td>0.299848</td>\n",
       "      <td>13383.223540</td>\n",
       "      <td>1.000093</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>humidity</th>\n",
       "      <td>-0.007016</td>\n",
       "      <td>0.138938</td>\n",
       "      <td>0.001320</td>\n",
       "      <td>-0.319773</td>\n",
       "      <td>0.276530</td>\n",
       "      <td>11849.915732</td>\n",
       "      <td>0.999995</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>temp_delta</th>\n",
       "      <td>0.004460</td>\n",
       "      <td>0.140156</td>\n",
       "      <td>0.001186</td>\n",
       "      <td>-0.285004</td>\n",
       "      <td>0.308285</td>\n",
       "      <td>13570.098506</td>\n",
       "      <td>1.000014</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>temperatureMin</th>\n",
       "      <td>-0.002262</td>\n",
       "      <td>0.140839</td>\n",
       "      <td>0.001176</td>\n",
       "      <td>-0.296004</td>\n",
       "      <td>0.302429</td>\n",
       "      <td>13765.821793</td>\n",
       "      <td>1.000014</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>apparentTemperatureMin</th>\n",
       "      <td>-0.002385</td>\n",
       "      <td>0.142660</td>\n",
       "      <td>0.001121</td>\n",
       "      <td>-0.321288</td>\n",
       "      <td>0.284931</td>\n",
       "      <td>13862.051877</td>\n",
       "      <td>1.000016</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>precipIntensity</th>\n",
       "      <td>-0.001106</td>\n",
       "      <td>0.142198</td>\n",
       "      <td>0.001353</td>\n",
       "      <td>-0.283684</td>\n",
       "      <td>0.304242</td>\n",
       "      <td>12239.307366</td>\n",
       "      <td>1.000001</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sd</th>\n",
       "      <td>86.437560</td>\n",
       "      <td>4.910928</td>\n",
       "      <td>0.051030</td>\n",
       "      <td>77.158364</td>\n",
       "      <td>96.024441</td>\n",
       "      <td>7238.262262</td>\n",
       "      <td>1.000086</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                             mean        sd  mc_error    hpd_2.5   hpd_97.5  \\\n",
       "Intercept                0.033369  0.147348  0.001307  -0.243261   0.372392   \n",
       "cyield                   0.008753  0.144545  0.001244  -0.286242   0.317488   \n",
       "irig_flag                0.023684  0.144673  0.001200  -0.261596   0.340129   \n",
       "days_under0              0.001703  0.140275  0.001228  -0.296166   0.293691   \n",
       "dewPoint                -0.005811  0.141128  0.001142  -0.286644   0.302634   \n",
       "precipAccumulation      -0.002344  0.139962  0.001232  -0.285032   0.303106   \n",
       "precip                  -0.000225  0.143678  0.001154  -0.290755   0.321854   \n",
       "days_under_n10          -0.000827  0.143340  0.001403  -0.308703   0.296381   \n",
       "days_over42              0.002368  0.144378  0.001394  -0.297255   0.311021   \n",
       "days_over32              0.004131  0.140339  0.001265  -0.285340   0.299848   \n",
       "humidity                -0.007016  0.138938  0.001320  -0.319773   0.276530   \n",
       "temp_delta               0.004460  0.140156  0.001186  -0.285004   0.308285   \n",
       "temperatureMin          -0.002262  0.140839  0.001176  -0.296004   0.302429   \n",
       "apparentTemperatureMin  -0.002385  0.142660  0.001121  -0.321288   0.284931   \n",
       "precipIntensity         -0.001106  0.142198  0.001353  -0.283684   0.304242   \n",
       "sd                      86.437560  4.910928  0.051030  77.158364  96.024441   \n",
       "\n",
       "                               n_eff      Rhat  \n",
       "Intercept               12752.112726  1.000212  \n",
       "cyield                  12215.797341  0.999988  \n",
       "irig_flag               13582.094278  0.999994  \n",
       "days_under0             12795.671660  1.000011  \n",
       "dewPoint                14326.564104  0.999988  \n",
       "precipAccumulation      13599.872401  1.000029  \n",
       "precip                  12685.048499  1.000184  \n",
       "days_under_n10          12484.341763  1.000030  \n",
       "days_over42             12704.563450  0.999994  \n",
       "days_over32             13383.223540  1.000093  \n",
       "humidity                11849.915732  0.999995  \n",
       "temp_delta              13570.098506  1.000014  \n",
       "temperatureMin          13765.821793  1.000014  \n",
       "apparentTemperatureMin  13862.051877  1.000016  \n",
       "precipIntensity         12239.307366  1.000001  \n",
       "sd                       7238.262262  1.000086  "
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pm.summary(trc_lasso)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ridge Regression \n",
    "This Lasso model behaves similarly to Ridge regression, which has a similar constraint on the sum of squares of the feature coefficients: the L2 norm. The effect of Ridge regression is to regularize the model by inhibiting large coefficient values and thus discouraging overfitting.\n",
    "\n",
    "I've decided to squeeze an example of Ridge regression into this post because I wanted to demonstrate a linear regression using only the 'important' features as found by the Lasso, and also since the default behaviour of the PyMC3 glm submodule is to conduct a Ridge regression.\n",
    "\n",
    "Set parameter  λ  to constrain the sum of squared values of  β"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import RidgeCV\n",
    "ridge_cv = RidgeCV(cv=5, normalize=False,alpha=0.1)\n",
    "ols_ridge = ridge_cv.fit(X_train_yield, y_train_yield)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train_yield.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "ols_ridge.coef_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "predict_ridge  = ols_ridge.predict(X_test_yield)\n",
    "rmse_ridge = sqrt(mean_squared_error(predict_ridge, y_test_yield))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We notice that the same coeffs that were zeroed out in Lasso have the lowest beta values here. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's complete the train of thought by implementing a Ridge regression in PyMC3. It is possible to imitate the effect of the L2 regulariation in Bayesian linear regression by using a Normal prior for feature coefficients. This is the default setting for the glm submodule in PyMC3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with pm.Model() as mdl_ridge:\n",
    "\n",
    "    ## Use GLM submodule for simplified model specification\n",
    "    ## Betas are Normal (as per default settings (for Ridge)\n",
    "    ## Likelihood is Normal \n",
    "    \n",
    "    pm.GLM.from_formula(formula, X_train_yield)\n",
    "   \n",
    "    ## take samples using NUTS sampler\n",
    "    trc_ridge = pm.sample(2000, step=pm.NUTS())\n",
    "    \n",
    "rvs_ridge = [rv.name for rv in mdl_ridge.unobserved_RVs]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pm.summary(trc_ridge)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}