_01_DAN_Python_Biogeme_Benchmark.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import OrderedDict    # For recording the model specification \n",
    "import pickle\n",
    "import pandas as pd                    # For file input/output\n",
    "import numpy as np                     # For vectorized math operations\n",
    "\n",
    "import pylogit as pl                   # For MNL model estimation\n",
    "                                       # To convert from wide to long format\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Load the Swissmetro Dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the raw swiss metro data\n",
    "# Note the .dat files are tab delimited text files\n",
    "# swissmetro_wide = pd.read_table(\"../data/swissmetro.dat\", sep='\\t')\n",
    "swissmetro_wide = pickle.load(open('../data/swissmetro_clean.pkl', \"rb\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Clean the dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the 01Logit.py file provided is an example from Python Biogeme (see: <a href=\"http://biogeme.epfl.ch/examples_swissmetro.html\">http://biogeme.epfl.ch/examples_swissmetro.html</a>). See http://www.strc.ch/conferences/2001/bierlaire1.pdf for a detailed explanation of the variables. The 01Logit.py file excludes observations meeting the following critera:\n",
    "<pre>\n",
    "exclude = (( PURPOSE != 1 ) * (  PURPOSE   !=  3  ) + ( CHOICE == 0 )) > 0\n",
    "</pre>\n",
    "As a result, their dataset has 6,768 observations. Below, I make the same exclusions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The cleaned number of observations is 10,692.\n"
     ]
    }
   ],
   "source": [
    "# Use \".copy()\" so that later on, we avoid performing operations \n",
    "# on a view of a dataframe as opposed to on an actual dataframe\n",
    "clean_sm_wide = swissmetro_wide.copy()\n",
    "\n",
    "# Look at how many observations we have after removing unwanted\n",
    "# observations\n",
    "final_num_obs = clean_sm_wide.shape[0]\n",
    "num_obs_statement = \"The cleaned number of observations is {:,.0f}.\"\n",
    "print (num_obs_statement.format(final_num_obs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Create an id column that ignores the repeat observations per individual\n",
    "\n",
    "In the simple example given on the Python Biogeme website for 01Logit.py, the repeated observations per individual are treated as separate and independent observations We will do the same"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a custom id column that ignores the fact that this is a \n",
    "# panel/repeated-observations dataset, and start the \"custom_id\" from 1\n",
    "clean_sm_wide[\"custom_id\"] = np.arange(clean_sm_wide.shape[0], dtype=int) + 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Convert the data from 'wide' format to 'long' format"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4a. Determine the 'type' of each column in the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['GROUP', 'SURVEY', 'SP', 'ID', 'PURPOSE', 'FIRST', 'TICKET', 'WHO',\n",
       "       'LUGGAGE', 'AGE', 'MALE', 'INCOME', 'GA', 'ORIGIN', 'DEST', 'TRAIN_AV',\n",
       "       'CAR_AV', 'SM_AV', 'TRAIN_TT', 'TRAIN_CO', 'TRAIN_HE', 'SM_TT', 'SM_CO',\n",
       "       'SM_HE', 'SM_SEATS', 'CAR_TT', 'CAR_CO', 'CHOICE', 'custom_id'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Look at the columns of the swissmetro data\n",
    "clean_sm_wide.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the list of individual specific variables\n",
    "ind_variables = clean_sm_wide.columns.tolist()[:15]\n",
    "\n",
    "# Specify the variables that vary across individuals **AND** \n",
    "# across some or all alternatives\n",
    "alt_varying_variables = {u'travel_time': dict([(1, 'TRAIN_TT'),\n",
    "                                               (2, 'SM_TT'),\n",
    "                                               (3, 'CAR_TT')]),\n",
    "                          u'travel_cost': dict([(1, 'TRAIN_CO'),\n",
    "                                                (2, 'SM_CO'),\n",
    "                                                (3, 'CAR_CO')]),\n",
    "                          u'headway': dict([(1, 'TRAIN_HE'),\n",
    "                                            (2, 'SM_HE')]),\n",
    "                          u'seat_configuration': dict([(2, \"SM_SEATS\")])}\n",
    "\n",
    "# Specify the availability variables\n",
    "availability_variables = dict(zip(range(1, 4), ['TRAIN_AV', 'SM_AV', 'CAR_AV']))\n",
    "\n",
    "# Determine the columns that will denote the\n",
    "# new column of alternative ids, and the columns\n",
    "# that denote the custom observation ids and the \n",
    "# choice column\n",
    "new_alt_id = \"mode_id\"\n",
    "obs_id_column = \"custom_id\"\n",
    "choice_column = \"CHOICE\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4b. Actually perform the conversion from wide to long formats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "      <th>5</th>\n",
       "      <th>6</th>\n",
       "      <th>7</th>\n",
       "      <th>8</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>custom_id</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mode_id</th>\n",
       "      <td>1.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CHOICE</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GROUP</th>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SURVEY</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SP</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ID</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PURPOSE</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>FIRST</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TICKET</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>WHO</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>LUGGAGE</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AGE</th>\n",
       "      <td>3.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MALE</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>INCOME</th>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GA</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ORIGIN</th>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>DEST</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>travel_time</th>\n",
       "      <td>112.0</td>\n",
       "      <td>63.0</td>\n",
       "      <td>117.0</td>\n",
       "      <td>103.0</td>\n",
       "      <td>60.0</td>\n",
       "      <td>117.0</td>\n",
       "      <td>130.0</td>\n",
       "      <td>67.0</td>\n",
       "      <td>117.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>travel_cost</th>\n",
       "      <td>48.0</td>\n",
       "      <td>52.0</td>\n",
       "      <td>65.0</td>\n",
       "      <td>48.0</td>\n",
       "      <td>49.0</td>\n",
       "      <td>84.0</td>\n",
       "      <td>48.0</td>\n",
       "      <td>58.0</td>\n",
       "      <td>52.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>headway</th>\n",
       "      <td>120.0</td>\n",
       "      <td>20.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>30.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>60.0</td>\n",
       "      <td>30.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>seat_configuration</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                        0     1      2      3     4      5      6     7      8\n",
       "custom_id             1.0   1.0    1.0    2.0   2.0    2.0    3.0   3.0    3.0\n",
       "mode_id               1.0   2.0    3.0    1.0   2.0    3.0    1.0   2.0    3.0\n",
       "CHOICE                0.0   1.0    0.0    0.0   1.0    0.0    0.0   1.0    0.0\n",
       "GROUP                 2.0   2.0    2.0    2.0   2.0    2.0    2.0   2.0    2.0\n",
       "SURVEY                0.0   0.0    0.0    0.0   0.0    0.0    0.0   0.0    0.0\n",
       "SP                    1.0   1.0    1.0    1.0   1.0    1.0    1.0   1.0    1.0\n",
       "ID                    1.0   1.0    1.0    1.0   1.0    1.0    1.0   1.0    1.0\n",
       "PURPOSE               1.0   1.0    1.0    1.0   1.0    1.0    1.0   1.0    1.0\n",
       "FIRST                 0.0   0.0    0.0    0.0   0.0    0.0    0.0   0.0    0.0\n",
       "TICKET                1.0   1.0    1.0    1.0   1.0    1.0    1.0   1.0    1.0\n",
       "WHO                   1.0   1.0    1.0    1.0   1.0    1.0    1.0   1.0    1.0\n",
       "LUGGAGE               0.0   0.0    0.0    0.0   0.0    0.0    0.0   0.0    0.0\n",
       "AGE                   3.0   3.0    3.0    3.0   3.0    3.0    3.0   3.0    3.0\n",
       "MALE                  0.0   0.0    0.0    0.0   0.0    0.0    0.0   0.0    0.0\n",
       "INCOME                2.0   2.0    2.0    2.0   2.0    2.0    2.0   2.0    2.0\n",
       "GA                    0.0   0.0    0.0    0.0   0.0    0.0    0.0   0.0    0.0\n",
       "ORIGIN                2.0   2.0    2.0    2.0   2.0    2.0    2.0   2.0    2.0\n",
       "DEST                  1.0   1.0    1.0    1.0   1.0    1.0    1.0   1.0    1.0\n",
       "travel_time         112.0  63.0  117.0  103.0  60.0  117.0  130.0  67.0  117.0\n",
       "travel_cost          48.0  52.0   65.0   48.0  49.0   84.0   48.0  58.0   52.0\n",
       "headway             120.0  20.0    0.0   30.0  10.0    0.0   60.0  30.0    0.0\n",
       "seat_configuration    0.0   0.0    0.0    0.0   0.0    0.0    0.0   0.0    0.0"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Perform the desired conversion\n",
    "long_swiss_metro = pl.convert_wide_to_long(clean_sm_wide, \n",
    "                                           ind_variables, \n",
    "                                           alt_varying_variables, \n",
    "                                           availability_variables, \n",
    "                                           obs_id_column, \n",
    "                                           choice_column,\n",
    "                                           new_alt_id_name=new_alt_id)\n",
    "\n",
    "# Look at the first 9 rows of the long-format dataframe\n",
    "long_swiss_metro.head(9).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. Create the variables used in the Python Biogeme Logit Model Example\n",
    "\n",
    "In 01Logit.py, the travel time and travel cost variables are scaled for ease of numeric optimization. We will do the same such that our estimated coefficients are comparable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/danph/anaconda3/envs/dev/lib/python3.8/site-packages/pandas/core/computation/expressions.py:203: UserWarning: evaluating in Python space because the '*' operator is not supported by numexpr for the bool dtype, use '&' instead\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "# Scale both the travel time and travel cost by 100\n",
    "long_swiss_metro[\"travel_time_hundredth\"] = (long_swiss_metro[\"travel_time\"] /\n",
    "                                             100.0)\n",
    "\n",
    "# Figure out which rows correspond to train or swiss metro \n",
    "# alternatives for individuals with GA passes. These individuals face no \n",
    "# marginal costs for a trip\n",
    "train_pass_train_alt = ((long_swiss_metro[\"GA\"] == 1) *\n",
    "                        (long_swiss_metro[\"mode_id\"].isin([1, 2]))).astype(int)\n",
    "# Note that the (train_pass_train_alt == 0) term accounts for the\n",
    "# fact that those with a GA pass have no marginal cost for the trip\n",
    "long_swiss_metro[\"travel_cost_hundredth\"] = (long_swiss_metro[\"travel_cost\"] *\n",
    "                                             (train_pass_train_alt == 0) /\n",
    "                                             100.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6. Specify and Estimate the Python Biogeme Logit Model Example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6a. Specify the Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the model's specification dictionary and variable names dictionary\n",
    "# NOTE: - Keys should be variables within the long format dataframe.\n",
    "#         The sole exception to this is the \"intercept\" key.\n",
    "#       - For the specification dictionary, the values should be lists\n",
    "#         or lists of lists. Within a list, or within the inner-most\n",
    "#         list should be the alternative ID's of the alternative whose\n",
    "#         utility specification the explanatory variable is entering.\n",
    "\n",
    "example_specification = OrderedDict()\n",
    "example_names = OrderedDict()\n",
    "\n",
    "# Note that 1 is the id for the Train and 3 is the id for the Car.\n",
    "# The next two lines are placing alternative specific constants in\n",
    "# the utility equations for the Train and for the Car. The order\n",
    "# in which these variables are placed is chosen so the summary\n",
    "# dataframe which is returned will match that shown in the HTML\n",
    "# file of the python biogeme example.\n",
    "example_specification[\"intercept\"] = [3, 1]\n",
    "example_names[\"intercept\"] = ['ASC Car', 'ASC Train']\n",
    "\n",
    "# Note that the names used below are simply for consistency with\n",
    "# the coefficient names given in the Python Biogeme example.\n",
    "# example_specification[\"travel_cost_hundredth\"] = [[1, 2, 3]]\n",
    "# example_names[\"travel_cost_hundredth\"] = ['B_COST']\n",
    "\n",
    "example_specification[\"travel_cost_hundredth\"] = [[1, 2, 3]]\n",
    "example_names[\"travel_cost_hundredth\"] = ['B_COST']\n",
    "\n",
    "example_specification[\"travel_time_hundredth\"] = [[1, 2, 3]]\n",
    "example_names[\"travel_time_hundredth\"] = ['B_TIME']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6b. Estimate the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Log-likelihood at zero: -11,071.2632\n",
      "Initial Log-likelihood: -11,071.2632\n",
      "Estimation Time for Point Estimation: 0.16 seconds.\n",
      "Final log-likelihood: -8,647.8792\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/danph/anaconda3/envs/dev/lib/python3.8/site-packages/scipy/optimize/_minimize.py:522: RuntimeWarning: Method BFGS does not use Hessian information (hess).\n",
      "  warn('Method %s does not use Hessian information (hess).' % method,\n"
     ]
    }
   ],
   "source": [
    "# Provide the module with the needed input arguments to create\n",
    "# an instance of the MNL model class\n",
    "example_mnl = pl.create_choice_model(data=long_swiss_metro,\n",
    "                                     alt_id_col=new_alt_id,\n",
    "                                     obs_id_col=obs_id_column,\n",
    "                                     choice_col=choice_column,\n",
    "                                     specification=example_specification,\n",
    "                                     model_type=\"MNL\",\n",
    "                                     names=example_names)\n",
    "\n",
    "# Start the model estimation from initial values of all zeros\n",
    "# i.e. 4 zeros for the 4 coefficients being estimated\n",
    "example_mnl.fit_mle(np.zeros(4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.c Compare the model output with that of Python Biogeme"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Multinomial Logit Model Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>         <td>CHOICE</td>          <th>  No. Observations:  </th>   <td>10,692</td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>         <td>Multinomial Logit Model</td> <th>  Df Residuals:      </th>   <td>10,688</td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>                  <td>MLE</td>           <th>  Df Model:          </th>      <td>4</td>     \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Tue, 02 Feb 2021</td>     <th>  Pseudo R-squ.:     </th>    <td>0.219</td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>09:30:25</td>         <th>  Pseudo R-bar-squ.: </th>    <td>0.219</td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>AIC:</th>                 <td>17,303.758</td>        <th>  Log-Likelihood:    </th> <td>-8,647.879</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>BIC:</th>                 <td>17,332.867</td>        <th>  LL-Null:           </th> <td>-11,071.263</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ASC Car</th>   <td>    0.0169</td> <td>    0.031</td> <td>    0.537</td> <td> 0.591</td> <td>   -0.045</td> <td>    0.078</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ASC Train</th> <td>   -0.6564</td> <td>    0.042</td> <td>  -15.659</td> <td> 0.000</td> <td>   -0.739</td> <td>   -0.574</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>B_COST</th>    <td>   -0.7891</td> <td>    0.036</td> <td>  -21.726</td> <td> 0.000</td> <td>   -0.860</td> <td>   -0.718</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>B_TIME</th>    <td>   -1.2773</td> <td>    0.043</td> <td>  -29.962</td> <td> 0.000</td> <td>   -1.361</td> <td>   -1.194</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                     Multinomial Logit Model Regression Results                    \n",
       "===================================================================================\n",
       "Dep. Variable:                      CHOICE   No. Observations:               10,692\n",
       "Model:             Multinomial Logit Model   Df Residuals:                   10,688\n",
       "Method:                                MLE   Df Model:                            4\n",
       "Date:                     Tue, 02 Feb 2021   Pseudo R-squ.:                   0.219\n",
       "Time:                             09:30:25   Pseudo R-bar-squ.:               0.219\n",
       "AIC:                            17,303.758   Log-Likelihood:             -8,647.879\n",
       "BIC:                            17,332.867   LL-Null:                   -11,071.263\n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "ASC Car        0.0169      0.031      0.537      0.591      -0.045       0.078\n",
       "ASC Train     -0.6564      0.042    -15.659      0.000      -0.739      -0.574\n",
       "B_COST        -0.7891      0.036    -21.726      0.000      -0.860      -0.718\n",
       "B_TIME        -1.2773      0.043    -29.962      0.000      -1.361      -1.194\n",
       "==============================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Look at the estimated coefficients and goodness-of-fit statistics\n",
    "example_mnl.get_statsmodels_summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Get accuracy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "              precision    recall  f1-score   support\n",
      "\n",
      "           1       0.68      0.01      0.02      1413\n",
      "           2       0.66      0.90      0.76      6199\n",
      "           3       0.66      0.47      0.55      3080\n",
      "\n",
      "    accuracy                           0.66     10692\n",
      "   macro avg       0.67      0.46      0.44     10692\n",
      "weighted avg       0.66      0.66      0.60     10692\n",
      "\n"
     ]
    }
   ],
   "source": [
    "def model_pred(data, model, alt_id_col, obs_id_col, choice_col):\n",
    "    data['predicted'] = model.predict(data)  \n",
    "    is_chosen = data.groupby([obs_id_col])['predicted'].idxmax()\n",
    "    data['predicted_choice'] = 0\n",
    "    data.loc[is_chosen.values,'predicted_choice'] = 1\n",
    "    \n",
    "    actual = data.loc[data[choice_column] ==1,alt_id_column]\n",
    "    pred = data.loc[data['predicted_choice'] ==1,alt_id_column]\n",
    "    return actual, pred   \n",
    "\n",
    "\n",
    "alt_id_column = \"mode_id\"\n",
    "obs_id_column = \"custom_id\"\n",
    "choice_column = \"CHOICE\"\n",
    "\n",
    "actual, predict = model_pred(long_swiss_metro, example_mnl, new_alt_id, obs_id_column, choice_column)\n",
    "\n",
    "from sklearn.metrics import classification_report\n",
    "print(classification_report(actual, predict))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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>parameters</th>\n",
       "      <th>std_err</th>\n",
       "      <th>t_stats</th>\n",
       "      <th>p_values</th>\n",
       "      <th>robust_std_err</th>\n",
       "      <th>robust_t_stats</th>\n",
       "      <th>robust_p_values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>ASC Car</th>\n",
       "      <td>0.016871</td>\n",
       "      <td>0.031410</td>\n",
       "      <td>0.537117</td>\n",
       "      <td>5.911867e-01</td>\n",
       "      <td>0.037121</td>\n",
       "      <td>0.454493</td>\n",
       "      <td>6.494738e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ASC Train</th>\n",
       "      <td>-0.656446</td>\n",
       "      <td>0.041921</td>\n",
       "      <td>-15.659217</td>\n",
       "      <td>2.874124e-55</td>\n",
       "      <td>0.054510</td>\n",
       "      <td>-12.042623</td>\n",
       "      <td>2.121026e-33</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>B_COST</th>\n",
       "      <td>-0.789140</td>\n",
       "      <td>0.036322</td>\n",
       "      <td>-21.726386</td>\n",
       "      <td>1.155465e-104</td>\n",
       "      <td>0.050933</td>\n",
       "      <td>-15.493837</td>\n",
       "      <td>3.818078e-54</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>B_TIME</th>\n",
       "      <td>-1.277273</td>\n",
       "      <td>0.042630</td>\n",
       "      <td>-29.961676</td>\n",
       "      <td>3.100099e-197</td>\n",
       "      <td>0.065605</td>\n",
       "      <td>-19.469057</td>\n",
       "      <td>2.009356e-84</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           parameters   std_err    t_stats       p_values  robust_std_err  \\\n",
       "ASC Car      0.016871  0.031410   0.537117   5.911867e-01        0.037121   \n",
       "ASC Train   -0.656446  0.041921 -15.659217   2.874124e-55        0.054510   \n",
       "B_COST      -0.789140  0.036322 -21.726386  1.155465e-104        0.050933   \n",
       "B_TIME      -1.277273  0.042630 -29.961676  3.100099e-197        0.065605   \n",
       "\n",
       "           robust_t_stats  robust_p_values  \n",
       "ASC Car          0.454493     6.494738e-01  \n",
       "ASC Train      -12.042623     2.121026e-33  \n",
       "B_COST         -15.493837     3.818078e-54  \n",
       "B_TIME         -19.469057     2.009356e-84  "
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Look at robust p-values in case one wants to see them\n",
    "example_mnl.summary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Summary\n",
    "My estimation results match those of Python Biogeme. <br>\n",
    "The Python Biogeme log-likelihood is -5,331.252 and their estimated parameters are:\n",
    "<pre>\n",
    "ASC Car:    -0.155\n",
    "ASC Train:  -0.701\n",
    "B_COST:     -1.08\n",
    "B_TIME:     -1.28\n",
    "</pre>\n",
    "\n",
    "As shown above, my log-likelihood is -5,331.252, and my estimated parameters are:\n",
    "<pre>\n",
    "ASC Car:    -0.1546\n",
    "ASC Train:  -0.7012\n",
    "B_COST:     -1.0838\t\n",
    "B_TIME:     -1.2779\n",
    "</pre>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test with Biogeme"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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>Value</th>\n",
       "      <th>Std err</th>\n",
       "      <th>t-test</th>\n",
       "      <th>p-value</th>\n",
       "      <th>Rob. Std err</th>\n",
       "      <th>Rob. t-test</th>\n",
       "      <th>Rob. p-value</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>ASC_CAR</th>\n",
       "      <td>-0.160781</td>\n",
       "      <td>0.043247</td>\n",
       "      <td>-3.717742</td>\n",
       "      <td>0.000201</td>\n",
       "      <td>0.058178</td>\n",
       "      <td>-2.763588</td>\n",
       "      <td>0.005717</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ASC_TRAIN</th>\n",
       "      <td>-0.716768</td>\n",
       "      <td>0.055016</td>\n",
       "      <td>-13.028317</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.082651</td>\n",
       "      <td>-8.672180</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>B_COST</th>\n",
       "      <td>-1.083953</td>\n",
       "      <td>0.051807</td>\n",
       "      <td>-20.922856</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.068156</td>\n",
       "      <td>-15.903978</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>B_TIME</th>\n",
       "      <td>-1.269009</td>\n",
       "      <td>0.056857</td>\n",
       "      <td>-22.319395</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.104124</td>\n",
       "      <td>-12.187436</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              Value   Std err     t-test   p-value  Rob. Std err  Rob. t-test  \\\n",
       "ASC_CAR   -0.160781  0.043247  -3.717742  0.000201      0.058178    -2.763588   \n",
       "ASC_TRAIN -0.716768  0.055016 -13.028317  0.000000      0.082651    -8.672180   \n",
       "B_COST    -1.083953  0.051807 -20.922856  0.000000      0.068156   -15.903978   \n",
       "B_TIME    -1.269009  0.056857 -22.319395  0.000000      0.104124   -12.187436   \n",
       "\n",
       "           Rob. p-value  \n",
       "ASC_CAR        0.005717  \n",
       "ASC_TRAIN      0.000000  \n",
       "B_COST         0.000000  \n",
       "B_TIME         0.000000  "
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#%%\n",
    "\"\"\"File 01logit.py\n",
    "\n",
    ":author: Michel Bierlaire, EPFL\n",
    ":date: Thu Sep  6 15:14:39 2018\n",
    "\n",
    " Example of a logit model.\n",
    " Three alternatives: Train, Car and Swissmetro\n",
    " SP data\n",
    "\"\"\"\n",
    "\n",
    "import pandas as pd\n",
    "import pickle\n",
    "import biogeme.database as db\n",
    "import biogeme.biogeme as bio\n",
    "import biogeme.models as models\n",
    "from biogeme.expressions import Beta\n",
    "\n",
    "# Read the data\n",
    "# df = pd.read_csv('swissmetro.dat', '\\t')\n",
    "# df = pd.read_table(\"../data/swissmetro.dat\", sep='\\t')\n",
    "df = pickle.load(open('../data/swissmetro_clean.pkl', \"rb\"))\n",
    "database = db.Database('swissmetro', df)\n",
    "\n",
    "# The following statement allows you to use the names of the variable\n",
    "# as Python variable.\n",
    "globals().update(database.variables)\n",
    "\n",
    "# Removing some observations\n",
    "exclude = ((PURPOSE != 1) * (PURPOSE != 3) + (CHOICE == 0)) > 0\n",
    "database.remove(exclude)\n",
    "\n",
    "# Parameters to be estimated\n",
    "ASC_CAR = Beta('ASC_CAR', 0, None, None, 0)\n",
    "ASC_TRAIN = Beta('ASC_TRAIN', 0, None, None, 0)\n",
    "ASC_SM = Beta('ASC_SM', 0, None, None, 1)\n",
    "B_TIME = Beta('B_TIME', 0, None, None, 0)\n",
    "B_COST = Beta('B_COST', 0, None, None, 0)\n",
    "\n",
    "\n",
    "# Definition of new variables\n",
    "SM_COST = SM_CO * (GA == 0)\n",
    "TRAIN_COST = TRAIN_CO * (GA == 0)\n",
    "CAR_AV_SP = CAR_AV * (SP != 0)\n",
    "TRAIN_AV_SP = TRAIN_AV * (SP != 0)\n",
    "TRAIN_TT_SCALED = TRAIN_TT / 100.0\n",
    "TRAIN_COST_SCALED = TRAIN_COST / 100\n",
    "SM_TT_SCALED = SM_TT / 100.0\n",
    "SM_COST_SCALED = SM_COST / 100\n",
    "CAR_TT_SCALED = CAR_TT / 100\n",
    "CAR_CO_SCALED = CAR_CO / 100\n",
    "\n",
    "# Definition of the utility functions\n",
    "V1 = ASC_TRAIN + \\\n",
    "     B_TIME * TRAIN_TT_SCALED + \\\n",
    "     B_COST * TRAIN_COST_SCALED\n",
    "V2 = ASC_SM + \\\n",
    "     B_TIME * SM_TT_SCALED + \\\n",
    "     B_COST * SM_COST_SCALED\n",
    "V3 = ASC_CAR + \\\n",
    "     B_TIME * CAR_TT_SCALED + \\\n",
    "     B_COST * CAR_CO_SCALED\n",
    "\n",
    "# Associate utility functions with the numbering of alternatives\n",
    "V = {1: V1,\n",
    "     2: V2,\n",
    "     3: V3}\n",
    "\n",
    "# Associate the availability conditions with the alternatives\n",
    "av = {1: TRAIN_AV_SP,\n",
    "      2: SM_AV,\n",
    "      3: CAR_AV_SP}\n",
    "\n",
    "# Definition of the model. This is the contribution of each\n",
    "# observation to the log likelihood function.\n",
    "logprob = models.loglogit(V, av, CHOICE)\n",
    "\n",
    "# Create the Biogeme object\n",
    "biogeme = bio.BIOGEME(database, logprob)\n",
    "biogeme.modelName = '01logit'\n",
    "\n",
    "# Estimate the parameters\n",
    "results = biogeme.estimate()\n",
    "\n",
    "# Get the results in a pandas table\n",
    "pandasResults = results.getEstimatedParameters()\n",
    "pandasResults"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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>parameters</th>\n",
       "      <th>std_err</th>\n",
       "      <th>t_stats</th>\n",
       "      <th>p_values</th>\n",
       "      <th>robust_std_err</th>\n",
       "      <th>robust_t_stats</th>\n",
       "      <th>robust_p_values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>ASC Car</th>\n",
       "      <td>0.016871</td>\n",
       "      <td>0.031410</td>\n",
       "      <td>0.537117</td>\n",
       "      <td>5.911867e-01</td>\n",
       "      <td>0.037121</td>\n",
       "      <td>0.454493</td>\n",
       "      <td>6.494738e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ASC Train</th>\n",
       "      <td>-0.656446</td>\n",
       "      <td>0.041921</td>\n",
       "      <td>-15.659217</td>\n",
       "      <td>2.874124e-55</td>\n",
       "      <td>0.054510</td>\n",
       "      <td>-12.042623</td>\n",
       "      <td>2.121026e-33</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>B_COST</th>\n",
       "      <td>-0.789140</td>\n",
       "      <td>0.036322</td>\n",
       "      <td>-21.726386</td>\n",
       "      <td>1.155465e-104</td>\n",
       "      <td>0.050933</td>\n",
       "      <td>-15.493837</td>\n",
       "      <td>3.818078e-54</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>B_TIME</th>\n",
       "      <td>-1.277273</td>\n",
       "      <td>0.042630</td>\n",
       "      <td>-29.961676</td>\n",
       "      <td>3.100099e-197</td>\n",
       "      <td>0.065605</td>\n",
       "      <td>-19.469057</td>\n",
       "      <td>2.009356e-84</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           parameters   std_err    t_stats       p_values  robust_std_err  \\\n",
       "ASC Car      0.016871  0.031410   0.537117   5.911867e-01        0.037121   \n",
       "ASC Train   -0.656446  0.041921 -15.659217   2.874124e-55        0.054510   \n",
       "B_COST      -0.789140  0.036322 -21.726386  1.155465e-104        0.050933   \n",
       "B_TIME      -1.277273  0.042630 -29.961676  3.100099e-197        0.065605   \n",
       "\n",
       "           robust_t_stats  robust_p_values  \n",
       "ASC Car          0.454493     6.494738e-01  \n",
       "ASC Train      -12.042623     2.121026e-33  \n",
       "B_COST         -15.493837     3.818078e-54  \n",
       "B_TIME         -19.469057     2.009356e-84  "
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# From Pylogit\n",
    "example_mnl.summary"
   ]
  },
  {
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}