<NME> Python Biogeme Benchmark--01Logit.ipynb
<BEF> {
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from collections import OrderedDict    # For recording the model specification \n",
    "\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": 2,
   "metadata": {
    "collapsed": true
   },
   "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')"
   ]
  },
  {
   "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": 3,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The cleaned number of observations is 6,768.\n"
     ]
    }
   ],
   "source": [
    "# Select obervations whose choice is known (i.e. CHOICE != 0)\n",
    "# **AND** whose PURPOSE is either 1 or 3\n",
    "include_criteria = (swissmetro_wide.PURPOSE.isin([1, 3]) &\n",
    "                    (swissmetro_wide.CHOICE != 0))\n",
    "\n",
    "# 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.loc[include_criteria].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": 4,
   "metadata": {
    "collapsed": false
   },
   "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": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index([u'GROUP', u'SURVEY', u'SP', u'ID', u'PURPOSE', u'FIRST', u'TICKET',\n",
       "       u'WHO', u'LUGGAGE', u'AGE', u'MALE', u'INCOME', u'GA', u'ORIGIN',\n",
       "       u'DEST', u'TRAIN_AV', u'CAR_AV', u'SM_AV', u'TRAIN_TT', u'TRAIN_CO',\n",
       "       u'TRAIN_HE', u'SM_TT', u'SM_CO', u'SM_HE', u'SM_SEATS', u'CAR_TT',\n",
       "       u'CAR_CO', u'CHOICE', u'custom_id'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Look at the columns of the swissmetro data\n",
    "clean_sm_wide.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "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": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\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>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",
       "    <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>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>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",
       "  </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",
       "seat_configuration    0.0   0.0    0.0    0.0   0.0    0.0    0.0   0.0    0.0\n",
       "travel_time         112.0  63.0  117.0  103.0  60.0  117.0  130.0  67.0  117.0\n",
       "headway             120.0  20.0    0.0   30.0  10.0    0.0   60.0  30.0    0.0\n",
       "travel_cost          48.0  52.0   65.0   48.0  49.0   84.0   48.0  58.0   52.0"
      ]
     },
     "execution_count": 7,
     "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": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/timothyb0912/anaconda/lib/python2.7/site-packages/pandas/computation/expressions.py:181: UserWarning: evaluating in Python space because the '*' operator is not supported by numexpr for the bool dtype, use '&' instead\n",
      "  unsupported[op_str]))\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": 9,
   "metadata": {
    "collapsed": true
   },
   "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": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Log-likelihood at zero: -6,964.6630\n",
      "Initial Log-likelihood: -6,964.6630\n",
      "Estimation Time: 0.03 seconds.\n",
      "Final log-likelihood: -5,331.2520\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/timothyb0912/anaconda/lib/python2.7/site-packages/scipy/optimize/_minimize.py:385: RuntimeWarning: Method BFGS does not use Hessian information (hess).\n",
      "  RuntimeWarning)\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": 11,
   "metadata": {
    "collapsed": false
   },
   "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>6,768</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>         <td>Multinomial Logit Model</td> <th>  Df Residuals:      </th>    <td>6,764</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, 30 Aug 2016</td>     <th>  Pseudo R-squ.:     </th>    <td>0.235</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>13:09:30</td>         <th>  Pseudo R-bar-squ.: </th>    <td>0.234</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>converged:</th>              <td>True</td>           <th>  Log-Likelihood:    </th> <td>-5,331.252</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th> </th>                         <td> </td>            <th>  LL-Null:           </th> <td>-6,964.663</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>[95.0% Conf. Int.]</th> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ASC Car</th>   <td>   -0.1546</td> <td>    0.043</td> <td>   -3.577</td> <td> 0.000</td> <td>   -0.239    -0.070</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ASC Train</th> <td>   -0.7012</td> <td>    0.055</td> <td>  -12.778</td> <td> 0.000</td> <td>   -0.809    -0.594</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>B_COST</th>    <td>   -1.0838</td> <td>    0.052</td> <td>  -20.910</td> <td> 0.000</td> <td>   -1.185    -0.982</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>B_TIME</th>    <td>   -1.2779</td> <td>    0.057</td> <td>  -22.465</td> <td> 0.000</td> <td>   -1.389    -1.166</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:                6,768\n",
       "Model:             Multinomial Logit Model   Df Residuals:                    6,764\n",
       "Method:                                MLE   Df Model:                            4\n",
       "Date:                     Tue, 30 Aug 2016   Pseudo R-squ.:                   0.235\n",
       "Time:                             13:09:30   Pseudo R-bar-squ.:               0.234\n",
       "converged:                            True   Log-Likelihood:             -5,331.252\n",
       "                                             LL-Null:                    -6,964.663\n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [95.0% Conf. Int.]\n",
       "------------------------------------------------------------------------------\n",
       "ASC Car       -0.1546      0.043     -3.577      0.000        -0.239    -0.070\n",
       "ASC Train     -0.7012      0.055    -12.778      0.000        -0.809    -0.594\n",
       "B_COST        -1.0838      0.052    -20.910      0.000        -1.185    -0.982\n",
       "B_TIME        -1.2779      0.057    -22.465      0.000        -1.389    -1.166\n",
       "==============================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Look at the estimated coefficients and goodness-of-fit statistics\n",
    "example_mnl.get_statsmodels_summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\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.154632</td>\n",
       "      <td>0.043235</td>\n",
       "      <td>-3.576518</td>\n",
       "      <td>3.482020e-04</td>\n",
       "      <td>0.066932</td>\n",
       "      <td>-2.310308</td>\n",
       "      <td>2.087112e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ASC Train</th>\n",
       "      <td>-0.701187</td>\n",
       "      <td>0.054874</td>\n",
       "      <td>-12.778138</td>\n",
       "      <td>2.172063e-37</td>\n",
       "      <td>0.130980</td>\n",
       "      <td>-5.353404</td>\n",
       "      <td>8.631471e-08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>B_COST</th>\n",
       "      <td>-1.083791</td>\n",
       "      <td>0.051830</td>\n",
       "      <td>-20.910412</td>\n",
       "      <td>4.305015e-97</td>\n",
       "      <td>0.204736</td>\n",
       "      <td>-5.293592</td>\n",
       "      <td>1.199367e-07</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>B_TIME</th>\n",
       "      <td>-1.277860</td>\n",
       "      <td>0.056883</td>\n",
       "      <td>-22.464576</td>\n",
       "      <td>9.218679e-112</td>\n",
       "      <td>0.229761</td>\n",
       "      <td>-5.561684</td>\n",
       "      <td>2.671842e-08</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.154632  0.043235  -3.576518   3.482020e-04        0.066932   \n",
       "ASC Train   -0.701187  0.054874 -12.778138   2.172063e-37        0.130980   \n",
       "B_COST      -1.083791  0.051830 -20.910412   4.305015e-97        0.204736   \n",
       "B_TIME      -1.277860  0.056883 -22.464576  9.218679e-112        0.229761   \n",
       "\n",
       "           robust_t_stats  robust_p_values  \n",
       "ASC Car         -2.310308     2.087112e-02  \n",
       "ASC Train       -5.353404     8.631471e-08  \n",
       "B_COST          -5.293592     1.199367e-07  \n",
       "B_TIME          -5.561684     2.671842e-08  "
      ]
     },
     "execution_count": 12,
     "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>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}

<MSG> Added example showing how to estimate a nested logit model.

<DFF> @@ -889,8 +889,9 @@
   }
  ],
  "metadata": {
+  "anaconda-cloud": {},
   "kernelspec": {
-   "display_name": "Python 2",
+   "display_name": "Python [default]",
    "language": "python",
    "name": "python2"
   },

Added example showing how to estimate a nested logit model.