hts_2018_simple_mnl_numpy_handestimate.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hand estimate the multinomial logit (mnl) model for dc metro hts 2017/2018 survey using pandas and numpy\n",
    "Written by ajit makhija (03/21/2021)\n",
    "\n",
    "first import libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\makhijaa\\Anaconda3\\lib\\site-packages\\IPython\\core\\interactiveshell.py:3146: DtypeWarning: Columns (24,25,27,29,52,55) have mixed types.Specify dtype option on import or set low_memory=False.\n",
      "  has_raised = await self.run_ast_nodes(code_ast.body, cell_name,\n"
     ]
    }
   ],
   "source": [
    "# read the data and filter hbw trips excluding outliers\n",
    "df = pd.read_csv('hts_data/trip_processed.csv')\n",
    "df = df[(df[\"PURP\"]==1) & (df[\"CHOICE\"]>0) & (df[\"SPEED\"]>3) & (df[\"SPEED\"]<100) & (df[\"DISTANCE\"]<=120) & (df[\"INC\"]==1)]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Remove unwanted columns "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.drop(columns=['DEPARTURE_TIME_HHMM','ARRIVAL_TIME_HHMM','MPO_TRAVEL_MODE','MPO_TRAVEL_MODE_DETAIL','MPO_TRANSIT_ACCESS_MODE','MPO_TRANSIT_EGRESS_MODE','MPO_SUBWAY_ACCESS_MODE','MPO_SUBWAY_EGRESS_MODE'], axis=1, inplace=True) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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>HOUSEHOLD_ID</th>\n",
       "      <th>PERSON_ID</th>\n",
       "      <th>PERSNO</th>\n",
       "      <th>TRIPID</th>\n",
       "      <th>TRIPNO</th>\n",
       "      <th>O_PURPOSE</th>\n",
       "      <th>O_ACTIVITY</th>\n",
       "      <th>MPO_O_ACTIVITY</th>\n",
       "      <th>O_STATE_FIPS</th>\n",
       "      <th>O_STATE_COUNTY_FIPS</th>\n",
       "      <th>...</th>\n",
       "      <th>ExBusCost</th>\n",
       "      <th>PURP</th>\n",
       "      <th>SovAV</th>\n",
       "      <th>HovAV</th>\n",
       "      <th>MetroAV</th>\n",
       "      <th>CrAV</th>\n",
       "      <th>BusAV</th>\n",
       "      <th>ExBusAV</th>\n",
       "      <th>SPEED</th>\n",
       "      <th>INC</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>12814</th>\n",
       "      <td>170001127</td>\n",
       "      <td>17000112701</td>\n",
       "      <td>1</td>\n",
       "      <td>1.700000e+12</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>24</td>\n",
       "      <td>24021</td>\n",
       "      <td>...</td>\n",
       "      <td>425</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>12.106953</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12815</th>\n",
       "      <td>170001127</td>\n",
       "      <td>17000112701</td>\n",
       "      <td>1</td>\n",
       "      <td>1.700000e+12</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>24</td>\n",
       "      <td>24021</td>\n",
       "      <td>...</td>\n",
       "      <td>425</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>12.106953</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12889</th>\n",
       "      <td>170001414</td>\n",
       "      <td>17000141401</td>\n",
       "      <td>1</td>\n",
       "      <td>1.700000e+12</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>51</td>\n",
       "      <td>51059</td>\n",
       "      <td>...</td>\n",
       "      <td>625</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>26.090401</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12890</th>\n",
       "      <td>170001414</td>\n",
       "      <td>17000141401</td>\n",
       "      <td>1</td>\n",
       "      <td>1.700000e+12</td>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>2</td>\n",
       "      <td>4</td>\n",
       "      <td>51</td>\n",
       "      <td>51510</td>\n",
       "      <td>...</td>\n",
       "      <td>625</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>43.484002</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12891</th>\n",
       "      <td>170001414</td>\n",
       "      <td>17000141401</td>\n",
       "      <td>1</td>\n",
       "      <td>1.700000e+12</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>51</td>\n",
       "      <td>51059</td>\n",
       "      <td>...</td>\n",
       "      <td>625</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>19.834972</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 76 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "       HOUSEHOLD_ID    PERSON_ID  PERSNO        TRIPID  TRIPNO  O_PURPOSE  \\\n",
       "12814     170001127  17000112701       1  1.700000e+12       1          1   \n",
       "12815     170001127  17000112701       1  1.700000e+12       2          2   \n",
       "12889     170001414  17000141401       1  1.700000e+12       1          1   \n",
       "12890     170001414  17000141401       1  1.700000e+12       2          3   \n",
       "12891     170001414  17000141401       1  1.700000e+12       3          1   \n",
       "\n",
       "       O_ACTIVITY  MPO_O_ACTIVITY  O_STATE_FIPS  O_STATE_COUNTY_FIPS  ...  \\\n",
       "12814           1               1            24                24021  ...   \n",
       "12815           2               3            24                24021  ...   \n",
       "12889           1               1            51                51059  ...   \n",
       "12890           2               4            51                51510  ...   \n",
       "12891           1               2            51                51059  ...   \n",
       "\n",
       "       ExBusCost  PURP  SovAV  HovAV  MetroAV  CrAV  BusAV  ExBusAV  \\\n",
       "12814        425     1      1      1        1     1      1        1   \n",
       "12815        425     1      1      1        1     1      1        1   \n",
       "12889        625     1      1      1        1     1      1        1   \n",
       "12890        625     1      1      1        1     1      1        1   \n",
       "12891        625     1      1      1        1     1      1        1   \n",
       "\n",
       "           SPEED  INC  \n",
       "12814  12.106953    1  \n",
       "12815  12.106953    1  \n",
       "12889  26.090401    1  \n",
       "12890  43.484002    1  \n",
       "12891  19.834972    1  \n",
       "\n",
       "[5 rows x 76 columns]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calcute the neccessary variables to estimate the maximum likelihood parameters using the Newton-Raphson method. "
   ]
  },
  {
   "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>count</th>\n",
       "      <th>mean</th>\n",
       "      <th>std</th>\n",
       "      <th>min</th>\n",
       "      <th>25%</th>\n",
       "      <th>50%</th>\n",
       "      <th>75%</th>\n",
       "      <th>max</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>HOUSEHOLD_ID</th>\n",
       "      <td>1562.0</td>\n",
       "      <td>1.771980e+08</td>\n",
       "      <td>4.574191e+06</td>\n",
       "      <td>1.700011e+08</td>\n",
       "      <td>1.700964e+08</td>\n",
       "      <td>1.800737e+08</td>\n",
       "      <td>1.801585e+08</td>\n",
       "      <td>1.802553e+08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PERSON_ID</th>\n",
       "      <td>1562.0</td>\n",
       "      <td>1.771980e+10</td>\n",
       "      <td>4.574191e+08</td>\n",
       "      <td>1.700011e+10</td>\n",
       "      <td>1.700964e+10</td>\n",
       "      <td>1.800737e+10</td>\n",
       "      <td>1.801585e+10</td>\n",
       "      <td>1.802553e+10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PERSNO</th>\n",
       "      <td>1562.0</td>\n",
       "      <td>1.405250e+00</td>\n",
       "      <td>7.070963e-01</td>\n",
       "      <td>1.000000e+00</td>\n",
       "      <td>1.000000e+00</td>\n",
       "      <td>1.000000e+00</td>\n",
       "      <td>2.000000e+00</td>\n",
       "      <td>6.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TRIPID</th>\n",
       "      <td>1562.0</td>\n",
       "      <td>1.770935e+12</td>\n",
       "      <td>4.542091e+10</td>\n",
       "      <td>1.700000e+12</td>\n",
       "      <td>1.700000e+12</td>\n",
       "      <td>1.800000e+12</td>\n",
       "      <td>1.800000e+12</td>\n",
       "      <td>1.800000e+12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TRIPNO</th>\n",
       "      <td>1562.0</td>\n",
       "      <td>2.111396e+00</td>\n",
       "      <td>1.704390e+00</td>\n",
       "      <td>1.000000e+00</td>\n",
       "      <td>1.000000e+00</td>\n",
       "      <td>2.000000e+00</td>\n",
       "      <td>2.000000e+00</td>\n",
       "      <td>2.300000e+01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>BusTTBusCost</th>\n",
       "      <td>1562.0</td>\n",
       "      <td>3.302125e+04</td>\n",
       "      <td>3.087361e+04</td>\n",
       "      <td>4.512000e+02</td>\n",
       "      <td>8.500000e+03</td>\n",
       "      <td>1.974325e+04</td>\n",
       "      <td>5.379375e+04</td>\n",
       "      <td>1.555125e+05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>BusCostBusCost</th>\n",
       "      <td>1562.0</td>\n",
       "      <td>2.602177e+05</td>\n",
       "      <td>1.560365e+05</td>\n",
       "      <td>4.000000e+04</td>\n",
       "      <td>4.000000e+04</td>\n",
       "      <td>3.906250e+05</td>\n",
       "      <td>3.906250e+05</td>\n",
       "      <td>3.906250e+05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ExBusTTExBusTT</th>\n",
       "      <td>1562.0</td>\n",
       "      <td>6.121018e+03</td>\n",
       "      <td>7.914975e+03</td>\n",
       "      <td>5.089536e+00</td>\n",
       "      <td>5.849143e+02</td>\n",
       "      <td>2.887451e+03</td>\n",
       "      <td>8.436606e+03</td>\n",
       "      <td>6.191139e+04</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ExBusTTExBusCost</th>\n",
       "      <td>1562.0</td>\n",
       "      <td>3.694001e+04</td>\n",
       "      <td>2.998353e+04</td>\n",
       "      <td>4.512000e+02</td>\n",
       "      <td>1.148250e+04</td>\n",
       "      <td>2.861938e+04</td>\n",
       "      <td>5.688531e+04</td>\n",
       "      <td>1.555125e+05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ExBusCostExBusCost</th>\n",
       "      <td>1562.0</td>\n",
       "      <td>3.145459e+05</td>\n",
       "      <td>1.266437e+05</td>\n",
       "      <td>4.000000e+04</td>\n",
       "      <td>1.806250e+05</td>\n",
       "      <td>3.906250e+05</td>\n",
       "      <td>3.906250e+05</td>\n",
       "      <td>3.906250e+05</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>100 rows × 8 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                     count          mean           std           min  \\\n",
       "HOUSEHOLD_ID        1562.0  1.771980e+08  4.574191e+06  1.700011e+08   \n",
       "PERSON_ID           1562.0  1.771980e+10  4.574191e+08  1.700011e+10   \n",
       "PERSNO              1562.0  1.405250e+00  7.070963e-01  1.000000e+00   \n",
       "TRIPID              1562.0  1.770935e+12  4.542091e+10  1.700000e+12   \n",
       "TRIPNO              1562.0  2.111396e+00  1.704390e+00  1.000000e+00   \n",
       "...                    ...           ...           ...           ...   \n",
       "BusTTBusCost        1562.0  3.302125e+04  3.087361e+04  4.512000e+02   \n",
       "BusCostBusCost      1562.0  2.602177e+05  1.560365e+05  4.000000e+04   \n",
       "ExBusTTExBusTT      1562.0  6.121018e+03  7.914975e+03  5.089536e+00   \n",
       "ExBusTTExBusCost    1562.0  3.694001e+04  2.998353e+04  4.512000e+02   \n",
       "ExBusCostExBusCost  1562.0  3.145459e+05  1.266437e+05  4.000000e+04   \n",
       "\n",
       "                             25%           50%           75%           max  \n",
       "HOUSEHOLD_ID        1.700964e+08  1.800737e+08  1.801585e+08  1.802553e+08  \n",
       "PERSON_ID           1.700964e+10  1.800737e+10  1.801585e+10  1.802553e+10  \n",
       "PERSNO              1.000000e+00  1.000000e+00  2.000000e+00  6.000000e+00  \n",
       "TRIPID              1.700000e+12  1.800000e+12  1.800000e+12  1.800000e+12  \n",
       "TRIPNO              1.000000e+00  2.000000e+00  2.000000e+00  2.300000e+01  \n",
       "...                          ...           ...           ...           ...  \n",
       "BusTTBusCost        8.500000e+03  1.974325e+04  5.379375e+04  1.555125e+05  \n",
       "BusCostBusCost      4.000000e+04  3.906250e+05  3.906250e+05  3.906250e+05  \n",
       "ExBusTTExBusTT      5.849143e+02  2.887451e+03  8.436606e+03  6.191139e+04  \n",
       "ExBusTTExBusCost    1.148250e+04  2.861938e+04  5.688531e+04  1.555125e+05  \n",
       "ExBusCostExBusCost  1.806250e+05  3.906250e+05  3.906250e+05  3.906250e+05  \n",
       "\n",
       "[100 rows x 8 columns]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df['SovChoice'] = np.where(df['CHOICE']==1,1,0)\n",
    "df['HovChoice'] = np.where(df['CHOICE']==2,1,0)\n",
    "df['MetroChoice'] = np.where(df['CHOICE']==3,1,0)\n",
    "df['CrChoice'] = np.where(df['CHOICE']==4,1,0)\n",
    "df['BusChoice'] = np.where(df['CHOICE']==5,1,0)\n",
    "df['ExBusChoice'] = np.where(df['CHOICE']==6,1,0)\n",
    "\n",
    "df['SovTTSovTT'] = df['SovTT']*df['SovTT']\n",
    "df['SovTTSovCost'] = df['SovTT']*df['SovCost']\n",
    "df['SovCostSovCost'] =  df['SovCost']*df['SovCost']\n",
    "\n",
    "df['HovTTHovTT'] = df['HovTT']*df['HovTT']\n",
    "df['HovTTHovCost'] = df['HovTT']*df['HovCost']\n",
    "df['HovCostHovCost'] = df['HovCost']*df['HovCost']\n",
    "\n",
    "df['MetroTTMetroTT'] = df['MetroTT']*df['MetroTT']\n",
    "df['MetroTTMetroCost'] = \t\tdf['MetroTT']*df['MetroCost']\n",
    "df['MetroCostMetroCost'] =   \tdf['MetroCost']*df['MetroCost']\n",
    "\n",
    "df['CrTTCrTT'] = df['CrTT']*df['CrTT']\n",
    "df['CrTTCrCost'] =df['CrTT']*df['CrCost']\n",
    "df['CrCostCrCost'] = df['CrCost']*df['CrCost']\n",
    "\n",
    "df['BusTTBusTT'] =df['BusTT']*df['BusTT']\n",
    "df['BusTTBusCost'] =df['BusTT']*df['BusCost']\n",
    "df['BusCostBusCost'] =df['BusCost']*df['BusCost']\n",
    "\n",
    "\n",
    "df['ExBusTTExBusTT'] = df['ExBusTT']*df['ExBusTT']\n",
    "df['ExBusTTExBusCost'] =df['ExBusTT']*df['ExBusCost']\n",
    "df['ExBusCostExBusCost'] =df['ExBusCost']*df['ExBusCost']\n",
    "\n",
    "df.describe().T\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#    This function hand estimates the multinomial logit model using the Newton Raphson method. The goal is to estimate only the coefficients for generic variables i.e., cost and time coefficients\n",
    "Computational steps are learned from Discrete Choice Analysis text book Ben-Akiva Moshe pg 82-84 and my course and thesis work with Dr Anthony Chen (USU)\n",
    "These include \n",
    "1)  Choose an initial guess betas for B0=0 (Time) and B1=0 (Cost)\n",
    "2)  Calculate utility functions and probabilites with initial betas\n",
    "3)  Calculate JacobianMatrix (JM), Hessian Matrix (HM), Inverse Hessian Matrix (invHM), Dot Product of Inverse Hessian & Jacobian\n",
    "4)  Calculate the updated betas; check the convergence if satisfied stop, if not go back to  step 1 \n",
    "\n",
    "Results match with Biogeme estimates; caveat this function only works when all modes are available for each of survey records; can't figure out how to estimate parameters if not all modes are available  \n",
    "\n",
    "![title](hts_data/solution.jpg)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def multinomiallogitestimate():\n",
    "\n",
    "    B0=0\n",
    "    B1=0\n",
    "    B0_prev=B0\n",
    "    B1_prev=B1\n",
    "    \n",
    "    iter=0\n",
    "    error=0\n",
    "    maxerror=10**(-10)\n",
    "    df1 = pd.DataFrame(columns=['iter','Beta_Time','Beta_Cost','LL','error','Beta_Time_adjust','Beta_Cost_adjust','rho_sq'])\n",
    "    for iter in range(0,100):\n",
    "        df[\"expUSov\"]= np.exp(B0*df[\"SovTT\"]+B1*df[\"SovCost\"])\n",
    "        df[\"expUHov\"]= np.exp(B0*df[\"HovTT\"]+B1*df[\"HovCost\"])\n",
    "        df[\"expUMetro\"]= np.exp(B0*df[\"MetroTT\"]+B1*df[\"MetroCost\"])\n",
    "        df[\"expUCr\"]= np.exp(B0*df[\"CrTT\"]+B1*df[\"CrCost\"])\n",
    "        df[\"expUBus\"]= np.exp(B0*df[\"BusTT\"]+B1*df[\"BusCost\"])\n",
    "        df[\"expUExBus\"]= np.exp(B0*df[\"ExBusTT\"]+B1*df[\"ExBusCost\"])\n",
    "\n",
    "        df[\"Sumexp\"] = df[\"expUSov\"]+df[\"expUHov\"]+df[\"expUMetro\"]+df[\"expUCr\"]+df[\"expUBus\"]+df[\"expUExBus\"]\n",
    "\n",
    "        df[\"ProbSov\"]= df[\"expUSov\"]/df[\"Sumexp\"]\n",
    "        df[\"ProbHov\"]= df[\"expUHov\"]/df[\"Sumexp\"]\n",
    "        df[\"ProbMetro\"]= df[\"expUMetro\"]/df[\"Sumexp\"]\n",
    "        df[\"ProbCr\"]= df[\"expUCr\"]/df[\"Sumexp\"]\n",
    "        df[\"ProbBus\"]= df[\"expUBus\"]/df[\"Sumexp\"]\n",
    "        df[\"ProbExBus\"]= df[\"expUExBus\"]/df[\"Sumexp\"]\n",
    "\n",
    "        df[\"ProbTT\"]= df[\"ProbSov\"]*df['SovTT'] + \\\n",
    "                      df[\"ProbHov\"]*df['HovTT'] + \\\n",
    "                      df[\"ProbMetro\"]*df['MetroTT'] + \\\n",
    "                      df[\"ProbCr\"]*df['CrTT'] + \\\n",
    "                      df[\"ProbBus\"]*df['BusTT'] + \\\n",
    "                      df[\"ProbExBus\"]*df['ExBusTT'] \n",
    "\n",
    "        df[\"ProbCost\"]=df[\"ProbSov\"] * df['SovCost'] + \\\n",
    "                       df[\"ProbHov\"] * df['HovCost'] + \\\n",
    "                       df[\"ProbMetro\"] * df['MetroCost'] + \\\n",
    "                       df[\"ProbCr\"] * df['CrCost'] + \\\n",
    "                       df[\"ProbBus\"] * df['BusCost'] + \\\n",
    "                       df[\"ProbExBus\"] * df['ExBusCost']  \n",
    "        # Calculate the functions to compute the Jacobian matrix Hessian matrix and Log Likelihood                       \n",
    "        df[\"F1\"]=(df[\"SovChoice\"] - df[\"ProbSov\"]) * df[\"SovTT\"] + \\\n",
    "                 (df[\"HovChoice\"] - df[\"ProbHov\"]) * df[\"HovTT\"] + \\\n",
    "                 (df[\"MetroChoice\"] - df[\"ProbMetro\"]) * df[\"MetroTT\"] + \\\n",
    "                 (df[\"CrChoice\"] - df[\"ProbCr\"]) * df[\"CrTT\"] + \\\n",
    "                 (df[\"BusChoice\"] - df[\"ProbBus\"]) * df[\"BusTT\"] + \\\n",
    "                 (df[\"ExBusChoice\"] - df[\"ProbExBus\"]) * df[\"ExBusTT\"]\n",
    "\n",
    "        df[\"F2\"]=(df[\"SovChoice\"]-df[\"ProbSov\"])*df[\"SovCost\"] + \\\n",
    "                 (df[\"HovChoice\"]-df[\"ProbHov\"])*df[\"HovCost\"] + \\\n",
    "                 (df[\"MetroChoice\"]-df[\"ProbMetro\"])*df[\"MetroCost\"] + \\\n",
    "                 (df[\"CrChoice\"]-df[\"ProbCr\"])*df[\"CrCost\"] + \\\n",
    "                 (df[\"BusChoice\"]-df[\"ProbBus\"])*df[\"BusCost\"] + \\\n",
    "                 (df[\"ExBusChoice\"]-df[\"ProbExBus\"])*df[\"ExBusCost\"]       \n",
    "\n",
    "        df[\"S11\"]=df[\"ProbTT\"]*df[\"ProbTT\"]- \\\n",
    "                  df[\"ProbSov\"]*df['SovTTSovTT'] - \\\n",
    "                  df[\"ProbHov\"]*df['HovTTHovTT'] - \\\n",
    "                  df[\"ProbMetro\"]*df['MetroTTMetroTT'] - \\\n",
    "                  df[\"ProbCr\"]*df['CrTTCrTT'] - \\\n",
    "                  df[\"ProbBus\"]*df['BusTTBusTT'] - \\\n",
    "                  df[\"ProbExBus\"]*df['ExBusTTExBusTT']\n",
    "\n",
    "        df[\"S12\"]=df[\"ProbTT\"]*df[\"ProbCost\"]- \\\n",
    "                  df[\"ProbSov\"]*df['SovTTSovCost'] - \\\n",
    "                  df[\"ProbHov\"]*df['HovTTHovCost'] - \\\n",
    "                  df[\"ProbMetro\"]*df['MetroTTMetroCost'] - \\\n",
    "                  df[\"ProbCr\"]*df['CrTTCrCost'] - \\\n",
    "                  df[\"ProbBus\"]*df['BusTTBusCost'] - \\\n",
    "                  df[\"ProbExBus\"]*df['ExBusTTExBusCost']\n",
    "\n",
    "        df[\"S22\"]=df[\"ProbCost\"]*df[\"ProbCost\"] - \\\n",
    "                 (df[\"ProbSov\"]*df['SovCostSovCost'] + \\\n",
    "                  df[\"ProbHov\"]*df['HovCostHovCost'] + \\\n",
    "                  df[\"ProbMetro\"]*df['MetroCostMetroCost'] + \\\n",
    "                  df[\"ProbCr\"]*df['CrCostCrCost'] + \\\n",
    "                  df[\"ProbBus\"]*df['BusCostBusCost'] + \\\n",
    "                  df[\"ProbExBus\"]*df['ExBusCostExBusCost'])      \n",
    "        df[\"LL\"]= df[\"SovChoice\"]*np.log(df[\"ProbSov\"]) + \\\n",
    "                  df[\"HovChoice\"]*np.log(df[\"ProbHov\"]) + \\\n",
    "                  df[\"MetroChoice\"]*np.log(df[\"ProbMetro\"]) + \\\n",
    "                  df[\"CrChoice\"]*np.log(df[\"ProbCr\"]) + \\\n",
    "                  df[\"BusChoice\"]*np.log(df[\"ProbBus\"]) + \\\n",
    "                  df[\"ExBusChoice\"]*np.log(df[\"ProbExBus\"]) \n",
    "\n",
    "\n",
    "        #Calculate  JacobianMatrix (JM), Hessian Matrix (HM), Inverse Hessian Matrix (invHM), Dot Product of Inverse Hessian Jacobian (invHMbyJM)\n",
    "        JM = np.array([[df[\"F1\"].sum()],\n",
    "                       [df[\"F2\"].sum()]])\n",
    "\n",
    "        HM = np.array([[df[\"S11\"].sum(),df[\"S12\"].sum()],\n",
    "                       [df[\"S12\"].sum(),df[\"S22\"].sum()]])\n",
    "\n",
    "        if iter==0: LL0 = df[\"LL\"].sum()\n",
    "        LL = df[\"LL\"].sum()\n",
    "        invHM = np.linalg.inv(HM)  \n",
    "        invHMbyJM = np.dot(invHM,JM)\n",
    "        B0_adjust = invHMbyJM[0,0]\n",
    "        B1_adjust = invHMbyJM[1,0]\n",
    "        # populate the dataframe with the outputs for each iter \n",
    "        df1.loc[iter, ['iter']]            = iter\n",
    "        df1.loc[iter, ['Beta_Time']]       = B0\n",
    "        df1.loc[iter, ['Beta_Cost']]       = B1\n",
    "        df1.loc[iter, ['LL']]              = LL       \n",
    "        df1.loc[iter, ['error']]           = error\n",
    "        df1.loc[iter, ['Beta_Time_adjust']]= B0_adjust       \n",
    "        df1.loc[iter, ['Beta_Cost_adjust']]=B1_adjust\n",
    "        df1.loc[iter, ['rho_sq']]=1-LL/LL0\n",
    "            \n",
    "        B0_prev=B0\n",
    "        B1_prev=B1\n",
    "        #new parameters equal diff of previous parameter minus dotproduct of invHessian and Jacobian\n",
    "        B0 = B0-invHMbyJM[0,0]\n",
    "        B1 = B1-invHMbyJM[1,0]\n",
    "        # check convergence (max error at 10^-10) and return results in df1\n",
    "\n",
    "        if iter==0:\n",
    "            sumerror = (B0)**2+(B1)**2\n",
    "        else:    \n",
    "            sumerror = (B0-B0_prev)**2+(B1-B1_prev)**2\n",
    "        error= (sumerror/2)**0.5\n",
    "         \n",
    "        if error < maxerror:\n",
    "            return df1\n",
    "            break   \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_by_iter=multinomiallogitestimate()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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>iter</th>\n",
       "      <th>Beta_Time</th>\n",
       "      <th>Beta_Cost</th>\n",
       "      <th>LL</th>\n",
       "      <th>error</th>\n",
       "      <th>Beta_Time_adjust</th>\n",
       "      <th>Beta_Cost_adjust</th>\n",
       "      <th>rho_sq</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>-2798.73</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0156825</td>\n",
       "      <td>0.00277204</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>-0.0156825</td>\n",
       "      <td>-0.00277204</td>\n",
       "      <td>-2267.14</td>\n",
       "      <td>0.0112611</td>\n",
       "      <td>0.014406</td>\n",
       "      <td>0.00115144</td>\n",
       "      <td>0.18994</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>-0.0300885</td>\n",
       "      <td>-0.00392348</td>\n",
       "      <td>-2173.03</td>\n",
       "      <td>0.0102191</td>\n",
       "      <td>0.00866328</td>\n",
       "      <td>0.000550069</td>\n",
       "      <td>0.223567</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3</td>\n",
       "      <td>-0.0387518</td>\n",
       "      <td>-0.00447355</td>\n",
       "      <td>-2159.26</td>\n",
       "      <td>0.0061382</td>\n",
       "      <td>0.00172036</td>\n",
       "      <td>0.000111842</td>\n",
       "      <td>0.228484</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4</td>\n",
       "      <td>-0.0404722</td>\n",
       "      <td>-0.00458539</td>\n",
       "      <td>-2158.89</td>\n",
       "      <td>0.00121904</td>\n",
       "      <td>5.1051e-05</td>\n",
       "      <td>3.39567e-06</td>\n",
       "      <td>0.228617</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>5</td>\n",
       "      <td>-0.0405232</td>\n",
       "      <td>-0.00458879</td>\n",
       "      <td>-2158.89</td>\n",
       "      <td>3.61783e-05</td>\n",
       "      <td>4.28226e-08</td>\n",
       "      <td>2.88489e-09</td>\n",
       "      <td>0.228617</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>6</td>\n",
       "      <td>-0.0405233</td>\n",
       "      <td>-0.00458879</td>\n",
       "      <td>-2158.89</td>\n",
       "      <td>3.03488e-08</td>\n",
       "      <td>3.02414e-14</td>\n",
       "      <td>2.04988e-15</td>\n",
       "      <td>0.228617</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  iter  Beta_Time   Beta_Cost       LL        error Beta_Time_adjust  \\\n",
       "0    0          0           0 -2798.73            0        0.0156825   \n",
       "1    1 -0.0156825 -0.00277204 -2267.14    0.0112611         0.014406   \n",
       "2    2 -0.0300885 -0.00392348 -2173.03    0.0102191       0.00866328   \n",
       "3    3 -0.0387518 -0.00447355 -2159.26    0.0061382       0.00172036   \n",
       "4    4 -0.0404722 -0.00458539 -2158.89   0.00121904       5.1051e-05   \n",
       "5    5 -0.0405232 -0.00458879 -2158.89  3.61783e-05      4.28226e-08   \n",
       "6    6 -0.0405233 -0.00458879 -2158.89  3.03488e-08      3.02414e-14   \n",
       "\n",
       "  Beta_Cost_adjust    rho_sq  \n",
       "0       0.00277204         0  \n",
       "1       0.00115144   0.18994  \n",
       "2      0.000550069  0.223567  \n",
       "3      0.000111842  0.228484  \n",
       "4      3.39567e-06  0.228617  \n",
       "5      2.88489e-09  0.228617  \n",
       "6      2.04988e-15  0.228617  "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results_by_iter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "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>iter</th>\n",
       "      <th>Beta_Time</th>\n",
       "      <th>Beta_Cost</th>\n",
       "      <th>LL</th>\n",
       "      <th>error</th>\n",
       "      <th>Beta_Time_adjust</th>\n",
       "      <th>Beta_Cost_adjust</th>\n",
       "      <th>rho_sq</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>6</td>\n",
       "      <td>-0.0405233</td>\n",
       "      <td>-0.00458879</td>\n",
       "      <td>-2158.89</td>\n",
       "      <td>3.03488e-08</td>\n",
       "      <td>3.02414e-14</td>\n",
       "      <td>2.04988e-15</td>\n",
       "      <td>0.228617</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  iter  Beta_Time   Beta_Cost       LL        error Beta_Time_adjust  \\\n",
       "6    6 -0.0405233 -0.00458879 -2158.89  3.03488e-08      3.02414e-14   \n",
       "\n",
       "  Beta_Cost_adjust    rho_sq  \n",
       "6      2.04988e-15  0.228617  "
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results_by_iter.tail(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEJCAYAAABYCmo+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAniklEQVR4nO3de5xV1Xn/8c8zM8wAA8MdRC6igoIimjBRE68BsaQata1W0zaSaEIxyS9N0zTG2rQxaZOq/cXU2JrQmESNv0ZrajSJKIcQbykYQUXlAGZAlOFyuAz3YWaYmef3x16jx3EY5rLn7DNzvu/X67xm77X2PufZ6Mxz1mWvbe6OiIhInIqSDkBERPoeJRcREYmdkouIiMROyUVERGKn5CIiIrErSTqAfDBy5EifNGlS0mGIiPQqK1eu3Onuo9qqU3IBJk2axIoVK5IOQ0SkVzGzN49Up24xERGJnZKLiIjETslFRERip+QiIiKxU3IREZHYKbmIiEjslFxERCR2us9FJE+4O/vqGtl1oJ6agw3sPNDA3kMNNDu4g+PhuHB81o6//R7R+7Qua9lv/YiNI73vO9vvrmuJsyPH+zsndOrfQXJr9rQxnD5haOzvq+Qi0kPcnf31jdQcaGDXwXp2HWhg18GGkDiiBNKSRGoORvuHm/rmH2KzpCOQIxld0V/JRSRJ7s7BhiZ2Hahn18EGdoWkECWHhrfLa96ua6ChqbnN9xpUVsLw8lKGl5cybmh/ThtXwYhBZYwoL2XEoFKGl0fbQwf2o6Qo6r02g7f/RlvLD3v7D3dLnYUCg6y6dw5q6/h3tt/7vmTVtRzZofdQRiloSi5SsNyd2oamd7UkWloXb3dNHQytigPRdkNj28liYGnx20lhTEV/ThlbwfBBpYwsL2N4SBgjysvCMaX071ec46sVyS0lF+lzag42sKmm9ohdUS2til0H66k73Hay6N+viBHlZYwcVMqoQWVMPaaCEaGl8e4WRpQ0BpQqWYhkU3KRPmNf3WHuWlrFj377xnvGLspKihg56J1WxJQxg97eH15eysisrqgRg0oZWKpfDZHu0G+Q9HpNzc5DKzbxr0+uo6a2gSvfP54/OPWYd3VFDSwt1hiASA4puUivtmz9Lr7+yzRrtu7jA5OG8eNLz+S08UOSDkuk4Cm5SK/01q5avvn4Gp5YvY1xQwdw15+9j0tOG6vWiUieUHKRXuVAfSN3La3ih8+9QUmx8aWLT+JT552g2VcieUbJRXqFpmbnZyurue3Jdew8UM8fv38cN86dypiK/kmHJiJtSDy5mNntwEeBBmA98El332Nmc4B/AUpD3d+6+9Jwzkzgx8AA4HHgr9zdzawMuA+YCewCrnb3jbm9Ionb796o4ZZfrGb1ln3MPG4Y98yr7JE7ikUkPvmwcGUKmO7uM4DXgZtC+U7go+5+GjAPuD/rnLuB+cCU8Jobyq8Hdrv7ZOAO4NaeD196yqaaWj77wIv86feXUXOwgX+75gweXvBBJRaRXiDxlou7L87aXQ5cGcpfyipfDfQPLZPhQIW7LwMws/uAK4BFwOXA18I5DwN3mZl569X6JK8drG/kP56q4j+ffYMig7++6CTmn3+CblQU6UUSTy6tXAc82Eb5nwAvuXu9mY0DqrPqqoFxYXscsAnA3RvNbC8wgqgV9C5mNp+o9cPEiRNjuwDpuuZm539e2sxtT6xl+/56rjjjWG78yFTGDhmQdGgi0kk5SS5mtgQ4po2qm9390XDMzUAj8ECrc08l6t66uKWojffxDtS9u9B9IbAQoLKyUi2bhK3YWMPXf5nmleq9nD5hKN/7+EzeP3FY0mGJSBflJLm4+0Xt1ZvZPOBSYHZ2F5aZjQceAa519/WhuBoYn3X6eGBLVt0EoNrMSoAhQE0sFyE9YvOeQ/zLorX8YtUWjqnozx1Xn87lp4+jqEj3q4j0Zol3i5nZXOBG4AJ3r80qHwr8CrjJ3X/bUu7uW81sv5mdDTwPXAt8N1Q/RjT4v4xo7GapxlvyU21DI997aj3ff2YDAJ+fNZkFF56oNb1E+oh8+E2+CygDUuHu6uXuvgD4HDAZ+KqZfTUce7G7bwdu4J2pyIvCC+Ae4H4zqyJqsVyTq4uQjmludh5dtZlbF61j2746Pnr6sXzlI1MZN1TjKiJ9iemLfTTmsmLFiqTD6PNefGs3X/9Fmpc37WHG+CH8w6WnUDlpeNJhiUgXmdlKd69sqy4fWi7Sx23de4hbF63l5y9vYfTgMv71qtP54/dpXEWkL1NykR5zqKGJhc9s4HtPr6fJnc9++EQ+c+Fkysv0v51IX6ffcomdu/PYqi3cumgtW/bWcclpY/nKR6YyYfjApEMTkRxRcpFYrdq0h6//Ms3KN3dz6rEV3HH1GZx1woikwxKRHFNykVhk9tVx6xNr+Z8XNzNyUBm3/ckM/mTmeIo1riJSkJRcpFvqDjfxg2c38B9PraexyVlwwYl89sMnMrh/v6RDE5EEKblIl7g7v3p1K996fC2b9xxi7qnHcNMfTuW4EeVJhyYieUDJRTrttc17+fov0vxuYw3TxlZw+1Uz+NCJI5MOS0TyiJKLdNj2/XXc/sQ6Hn6xmuEDS/nmH53G1R+YoHEVEXkPJRc5qrrDTfzwt2/w70uraGhq5tPnncDnZk2mQuMqInIESi5yRO7OE69t45uL1rCp5hBzThnD3/3hNI4fqXEVEWmfkou0afWWaFzl+TdqOHnMYB741FmcM1njKiLSMUou8i479tfz7dQ6fvrCJoYO6Mc3rpjOxz4wgZLioqRDE5FeRMlFAKhvbOLHv93Id5dWUXe4ievOOZ7Pz57CkAEaVxGRzlNyEQD+5qFV/PKVrcyaOpqbL5nGiaMGJR2SiPRiSi7CgfpGFq/O8PGzj+MbV0xPOhwR6QPUkS488/oOGpqauXTG2KRDEZE+QslFSKUzDBvYj5nHDUs6FBHpI5RcCtzhpmaWrt3OrKljNCNMRGKjvyYF7oWNNew9dJg5p4xJOhQR6UOUXApcKp2hrKSI80/SDZIiEh8llwLm7qTSGc6dPJKBpZo4KCLxSTy5mNntZrbWzF4xs0fMbGgoP9PMXg6vVWb2R1nnzDSzV82syszuNDML5WVm9mAof97MJiVzVb3D2m37qd59SF1iIhK7xJMLkAKmu/sM4HXgplD+GlDp7mcAc4Hvm1nL1+u7gfnAlPCaG8qvB3a7+2TgDuDWnFxBL5VKZzCD2dOUXEQkXoknF3df7O6NYXc5MD6U12aV9wccwMzGAhXuvszdHbgPuCIcdzlwb9h+GJjd0qqR90qlM7xvwlBGDS5LOhQR6WMSTy6tXAcsatkxs7PMbDXwKrAgJJtxQHXWOdWhjPBzE0A4di8woq0PMrP5ZrbCzFbs2LEj9gvJd1v2HOLVzXuZc8oxSYciIn1QTpKLmS0xs9faeF2edczNQCPwQEuZuz/v7qcCHwBuMrP+QFstEW95m3bq3l3ovtDdK929ctSoUV29tF5ryZoMgMZbRKRH5GSKkLtf1F69mc0DLgVmh66u1uevMbODwHSilsr4rOrxwJawXQ1MAKrD+MwQoKb7V9D3pNIZThhZzuTRWqBSROKXeLeYmc0FbgQuc/farPLjWwbwzew44GRgo7tvBfab2dlhPOVa4NFw2mPAvLB9JbC0rWRV6PbVHWb5hl1qtYhIj8mHmxvuAsqAVBh7X+7uC4Bzga+Y2WGgGfiMu+8M59wA/BgYQDRG0zJOcw9wv5lVEbVYrsnVRfQmT63bweEmV3IRkR6TeHIJ04bbKr8fuP8IdSuIushal9cBV8UaYB+USmcYUV7K+yZqoUoR6RmJd4tJbjU0NvPU2u1cNG0MxUWapS0iPUPJpcA8/8Yu9tc3qktMRHqUkkuBSaUzDOhXzLlTtFCliPQcJZcC4u4sSWc4b8pI+vcrTjocEenDlFwKyOot+9iyt05dYiLS45RcCsjidIYiLVQpIjmg5FJAUukMlccNZ3h5adKhiEgfp+RSIDbV1LJm6z51iYlITii5FIhUWgtVikjuKLkUiFQ6w5TRg5g0sjzpUESkACi5FIA9tQ38bmONWi0ikjNKLgXgN+u209SshSpFJHeUXApAKp1h9OAyTh8/NOlQRKRAtLsqspndzxGe5JjN3a+NLSKJVX1jE0+v28FlZ4yjSAtVikiOHK3lUgWsD6+9wBVAMdETH4uAy4E9PReedNf/rt/FwYYmLj5VXWIikjvttlzc/ZaWbTN7ErjE3Z/NKjsX+GrPhSfdlUpnKC8t5kMnjkg6FBEpIJ0ZczkbWN6q7Hngg/GFI3Fqbo4Wqrzg5FGUlWihShHJnc4kl5eAb5rZAIDw85+Bl3sgLonBK5v3sn1/vWaJiUjOdSa5fAI4B9hrZhmiMZhzAQ3m56lUehvFRcaHTx6ddCgiUmDaHXPJ5u4bgQ+Z2QTgWGCru7/VU4FJ96XSGc6cNJyhA7VQpYjkVqfuczGzYcCHgVnAhWFf8tDGnQd5PXNAXWIikogOJxcz+yDRlOQFwAzgL4H1oVzyjBaqFJEkdabl8h3gM+7+IXf/mLufA9wA3NmdAMzsdjNba2avmNkjZja0Vf1EMztgZl/KKptpZq+aWZWZ3WlmFsrLzOzBUP68mU3qTmy9WSqdYeoxg5kwfGDSoYhIAepMcjkJeKhV2cPA5G7GkAKmu/sM4HXgplb1dwCLWpXdDcwHpoTX3FB+PbDb3SeH827tZmy9Us3BBla8WcPFarWISEI6k1x+D1zTquwqoq6yLnP3xe7eGHaXA+Nb6szsCmADsDqrbCxQ4e7L3N2B+4hWDoBoxYB7w/bDwOyWVk0h+fWaDM0Oc045JulQRKRAdXi2GPAF4Jdm9nngTWASUavh0hjjuQ54EMDMyoEbgTnAl7KOGUe0/EyL6lDWUrcJwN0bzWwvMALY2fqDzGw+UeuHiRMnxngJyUulM4wd0p/p4yqSDkVEClRnpiL/r5mdCFxCNBX5F8Dj7l5ztHPNbAnQ1tfom9390XDMzUAj8ECouwW4w90PtGp8tNUS8Q7UvbvQfSGwEKCysvKoi3P2FnWHm3j29zu5cuZ4CrDRJiJ5ojMtF9x9t5k9Q9RC2NyRxBLOu6i9ejObR9QCmh26ugDOAq40s9uAoUCzmdUBPyOr6yxsbwnb1cAEoNrMSoAhQIdi7Cue+/1ODh1u0iwxEUlUh5NLGOv4KdEaYzXACDNbDlzj7lvaPbn9951L1P11gbvXtpS7+3lZx3wNOODud4X9/WZ2NtHaZtcC3w2HPgbMA5YBVwJLs5JVQUilMwwuK+HsE7RQpYgkpzMD+ncDq4Dh7j4WGEa03tj3uhnDXcBgIGVmL5tZR97vBuAHvPNIgJbZZPcQJb0q4IvAV7oZW6/S1Oz8em2GC6eOprREz4ETkeR0plvsXGCsux8GcPeDZvZlYHN3AgjTho92zNda7a8AprdxXB3RDLaC9PKm3ew80KAuMRFJXGe+3u4GTmlVdjJ6WFjeWJzO0K/YuPDkUUmHIiIFrjMtl9uAJWZ2D9FU5OOAT6KHheWNVDrD2SeMoKJ/v6RDEZEC1+GWi7v/J3A1MBL4aPj5sTClVxJWtf0AG3YcVJeYiOSFzk5FXgos7aFYpBtaFqq8aJqSi4gkrzNTkUuJHhh2BjAou87d9cCwhKXS25g+roJjhw5IOhQRkU61XO4FTie6Mz/TM+FIV+zYX89Lm/bwhdknJR2KiAjQueQyFzje3ff0UCzSRb9ek8Fdz24RkfzRmanIbwFlPRWIdF0qnWHc0AFMGzs46VBERICjtFzMbFbW7n3Ao2b2b7TqFgsD/ZKA2oZGnqvaycfOnKiFKkUkbxytW+yeNsq+2WrfgRPiCUc665nXd1Lf2KwHg4lIXmk3ubj78bkKRLomlc5Q0b+EDxw/POlQRETeptUNe7HGpmaWrs0wa+po+hXrP6WI5I+jjbmscfdpYXsTR37wVt96lGMvsfLN3eyuPczFp+pxxiKSX4425vLprO2/6MlApPNS6QylxUWcf5IWqhSR/HK0MZfnsraf7vlwpKPcndSaDB+aPIJBZZ1axUdEpMcdrVvs6x15E3f/h3jCkY76/fYDvLmrlvnna6KeiOSfo33lnZCTKKTTFq/eBmihShHJT0frFvtkrgKRzkmlM5w+YShjKvonHYqIyHt0av6qmU0zs6+a2V1h/2Qzm9EzocmRZPbVsap6r26cFJG81eHkYmZXAc8A44CWJfYHA9/ugbikHS3PbtFClSKSrzrTcvk6MMfdFwBNoWwV0TL8kkOpdIbjRgxkyuhBRz9YRCQBnUkuo4mSCbxzM6VzhBsrO8rMbjeztWb2ipk9YmZDQ/kkMztkZi+H1/eyzplpZq+aWZWZ3WlhxUYzKzOzB0P582Y2qTux5aMD9Y0sW7+LOdPGaKFKEclbnUkuK4GPtyq7BvhdN2NIAdPdfQbwOnBTVt16dz8jvBZkld8NzAemhNfcUH49sNvdJwN3ALd2M7a88/S6HTQ0NatLTETyWmeSy+eBfzKzp4FyM3sS+Abw190JwN0Xu3tj2F0OjG/veDMbC1S4+zJ3d6JHAVwRqi8nemImwMPAbOtjX+9T6W0MG9iPmccNSzoUEZEj6vCt3e6+1symApcCvwQ2hZ8NMcZzHfBg1v7xZvYSsA/4e3d/lmhCQXXWMdWhjPBzU4i30cz2AiOAnTHGmJjDTc0sXbudOaccQ4kWqhSRPNbh5GJm33H3LwAPZZWVAY/xTrfUkc5dArS1uuLN7v5oOOZmoBF4INRtBSa6+y4zmwn83MxOBdpqibSM+7RX1zqm+URda0yc2DvW3XzhjRr21TWqS0xE8l5nFqWaYWa3uPs/ApjZAOBXhJZCe9z9ovbqzWweUYtodujqwt3rgfqwvdLM1gMnEbVUsrvOxgNbwnY10aoC1WZWAgwBao4Q00JgIUBlZWW3JiXkyuJ0hrKSIs4/aWTSoYiItKszfSuXA3PN7ItmVg48CawHPtGdAMxsLnAjcJm712aVjzKz4rB9AtHA/QZ33wrsN7Ozw3jKtcCj4bTHgHlh+0pgaUuy6u3cnVQ6w3lTRjKwVAtVikh+63Bycff9wEeATwIvA6vc/dMx/PG+i+hmzFSrKcfnA6+Y2SqiwfkF7t7SCrkB+AFQRZTgFoXye4ARZlYFfBH4Sjdjyxtrtu5n855D6hITkV6hK6si/w64BNjdUt+dVZHDtOG2yn8G/OwIdSuA6W2U1wFXdTWWfJZKZzCDWVOVXEQk/3V1VeRF7dRJD1ic3sb7Jw5j1OCypEMRETkqrYrcC2zec4jVW/bxlY9MTToUEZEOOVq32CR33xi2j/hUKnffEHNckmWJFqoUkV7maN1irxINtkM0eO68914SB4pjjkuypNIZThhVzomjtFCliPQO7c4Wc/fBWdtF7l4cfma/lFh60N5Dh1m+YZdaLSLSq3RrDREz62dmz8QVjLzXU+u209jsejCYiPQq3V2gqgg4J45ApG2pdIaRg0o5Y4IWqhSR3kOrH+axhsZmnl63g9lTx1Bc1KcWdxaRPk7JJY8t37CL/fVaqFJEep+jLlJ1hLv0O3y+dF0qnWFAv2LOnaKFKkWkd+lIcjjanfj3xRGIvJu7s2RNhvNPGkn/fpqQJyK9y1GTi+7ST8Zrm/exdW8df3PxyUmHIiLSaZ15WNiR7tCvB7a6e3M8IQlEjzMuMpg1dXTSoYiIdFpnxkxa7tCH6C797KX2m83sMeAz7p6JK7hCtjidoXLScIaXlyYdiohIp3VmttiniR5BfBLQHzgZ+AnwGeA0okT173EHWIg21dSydtt+3TgpIr1WZ1outwCTwzNTAKrM7AbgdXf/vpl9Avh93AEWosVaqFJEernOtFyKgEmtyibyzqKVB9DU5Fik0ts4acwgjhtRnnQoIiJd0plk8B1gqZn9CNgEjCd65PF3Qv0lwLI4gytEe2obeGHjbhZccMQnHIiI5L0OJxd3v83MXiF6jPD7ga3A9e7+RKj/OfDzHoixoCxdu52mZmfOKcckHYqISJd1qhsrJJIneigWIborf/TgMmaMG5J0KCIiXdbhMZewvP4tZrbBzOrCz1vMTHNlY1J3uImnX9/BRaeMoUgLVYpIL9aZlsttwJnAAuBN4Djgq0AF8Nfxh1Z4lq3fRW1Dk2aJiUiv15nZYlcBl7n7Yndf5+6LgT8C/rQ7AZjZ7Wa21sxeMbNHzGxoVt0MM1tmZqvN7FUz6x/KZ4b9KjO708wslJeZ2YOh/Hkzm9Sd2HJtcTpDeWkxHzpxRNKhiIh0S2eSy5H6abrbf5MCprv7DOB14CYAMyshuklzgbufClwIHA7n3A3MB6aE19xQfj2w290nA3cAt3Yztpxpbo4Wqrzg5FGUlWihShHp3TqTXP4b+IWZ/YGZTTOzuUSzwx7qTgChJdQYdpcTTXEGuBh4xd1XheN2uXuTmY0FKtx9mbs70arMV4RzLgfuDdsPA7NbWjX5blX1Hnbsr+dizRITkT6gM8nly8ASoiVeVgLfBX4DNMQYz3XAorB9EuBm9qSZvWhmXw7l44DqrHOqQ1lL3SaAkLD2Ar2ijymVzlBcZHz4ZC1UKSK9X2fuc2kA/iG8AAhjIAeJEs8RmdkSoK2v5De7+6PhmJuBRqL1y1piOxf4AFAL/NrMVgL72gqv5aPaqWsd03yirjUmTpzYXvg5sTid4azjhzNkYL+kQxER6bbuLtfidGDMxd0vaq/ezOYBlwKzQ1cXRC2Sp919ZzjmcaKbN3/CO11nhO0tWedMAKrDmM0QoOYIMS0EFgJUVla2mYBy5Y2dB6nafoA/Pyv5JCciEofOdIsdSbf+MIexmxuJZqLVZlU9Ccwws4EhUVwApN19K7DfzM4O4ynXAo+Gcx4D5oXtK4GlWckqb6XS2wAtVCkifcdRWy5mNqud6jhuoLwLKANSYex9ubsvcPfdZvZt4AWiBPa4u/8qnHMD8GNgANEYTcs4zT3A/WZWRdRiuSaG+HpcKp1h2tgKxg8bmHQoIiKx6Ei32D1HqX+rOwGEacNHqvsJUTdY6/IVwPQ2yuuI7sfpNXYdqGflm7v53KwpSYciIhKboyYXdz8+F4EUql+v3U6zoweDiUifEseYi3RDKp3h2CH9OfXYiqRDERGJjZJLgg41NPHs76OFKnvJvZ4iIh2i5JKg56p2Une4WbPERKTPUXJJUCq9jcFlJZx1fK9YREBEpMOUXBLS1Oz8es12Lpw6mtIS/WcQkb5Ff9US8tJbu9l1sEFdYiLSJym5JCSVztCv2Ljw5FFJhyIiEjsllwS4O4vTGc4+YQQV/bVQpYj0PUouCVi/4wBv7DyoGydFpM9ScknA4nQGgIuUXESkj1JySUAqneG0cUMYO2RA0qGIiPQIJZcc276/jpc37dEsMRHp05RccuzXa7bjrme3iEjfpuSSY6l0hvHDBjD1mMFJhyIi0mOUXHLoYH0jz1XtZI4WqhSRPk7JJYee/f0OGhq1UKWI9H1KLjm0OJ1hyIB+nDlpeNKhiIj0KCWXHGlsambp2u3MmjqakmL9s4tI36a/cjmy4s3d7Kk9rC4xESkISi45kkpnKC0u4vyTtFCliPR9Si454O6k0hk+NHkEg8pKkg5HRKTHJZ5czOx2M1trZq+Y2SNmNjSU/7mZvZz1ajazM0LdTDN71cyqzOxOC/N6zazMzB4M5c+b2aTELizLusx+3qqp5eJTjkk6FBGRnEg8uQApYLq7zwBeB24CcPcH3P0Mdz8D+Diw0d1fDufcDcwHpoTX3FB+PbDb3ScDdwC35uoi2pNaHRaqnDY64UhERHIj8eTi7ovdvTHsLgfGt3HYx4D/AjCzsUCFuy9zdwfuA64Ix10O3Bu2HwZmWx7crZhak+GMCUMZXdE/6VBERHIi8eTSynXAojbKryYkF2AcUJ1VVx3KWuo2AYSEtRcY0dYHmdl8M1thZit27NgRQ+ht27a3jleq92qWmIgUlJyMLpvZEqCtAYeb3f3RcMzNQCPwQKtzzwJq3f21lqI23sc7UPfuQveFwEKAysrKNo+JQ2pN1CWmB4OJSCHJSXJx94vaqzezecClwOzQ1ZXtGt5ptUDUUsnuOhsPbMmqmwBUm1kJMASo6Ubo3ZZKZ5g0YiCTRw9KMgwRkZxKvFvMzOYCNwKXuXttq7oi4Crgpy1l7r4V2G9mZ4fxlGuBR0P1Y8C8sH0lsLSNZJUz++sOs2y9FqoUkcKTDzdd3AWUAanwB3i5uy8IdecD1e6+odU5NwA/BgYQjdG0jNPcA9xvZlVELZZrejb09j39+g4ONzlzNAVZRApM4sklTBs+Ut1TwNltlK8AprdRXkfU0skLqXSG4eWlzDxuWNKhiIjkVOLdYn3V4aZmfhMWqiwuUpeYiBQWJZce8rs3athX16gpyCJSkJRcekgqnaGspIjzpoxMOhQRkZxTcukBLQtVnjdlJANLEx/WEhHJOSWXHpDeuo/New6pS0xECpaSSw9YvDqDGcyepuQiIoVJyaUHpNIZZk4cxshBZUmHIiKSCCWXmFXvriW9dZ+6xESkoCm5xGxJOlqoUslFRAqZkkvMUmsynDiqnBNGaaFKESlcSi4x2nvoMM9vqNFaYiJS8JRcYvTUuu00Nru6xESk4Cm5xGhxOsPIQWW8b8LQpEMREUmUkktM6hubeHrdDi6aNpoiLVQpIgVOySUmyzfUcKBeC1WKiICSS2xS6W0M6FfMOZO1UKWIiJJLDNydJentnH/SSPr3K046HBGRxCm5xODVzXvZtq9OU5BFRAIllxgsXp2hyGDW1NFJhyIikheUXGKQSmf4wKThDC8vTToUEZG8oOTSTW/tqmVdZr9miYmIZFFy6abF6W0AXKzxFhGRtyWeXMzsdjNba2avmNkjZjY0lPczs3vN7FUzW2NmN2WdMzOUV5nZnWZmobzMzB4M5c+b2aSejj+VznDymMFMHDGwpz9KRKTXSDy5AClgurvPAF4HWpLIVUCZu58GzAT+MitZ3A3MB6aE19xQfj2w290nA3cAt/Zk4LsPNvDCxhp1iYmItJJ4cnH3xe7eGHaXA+NbqoByMysBBgANwD4zGwtUuPsyd3fgPuCKcM7lwL1h+2FgdkurpicsXbudZtezW0REWks8ubRyHbAobD8MHAS2Am8B/+ruNcA4oDrrnOpQRvi5CSAkrL3AiLY+yMzmm9kKM1uxY8eOLgVbMaAfc04Zw2njhnTpfBGRvqokFx9iZkuAtka8b3b3R8MxNwONwAOh7kygCTgWGAY8G96nrZaIt3xUO3XvLnRfCCwEqKysbPOYo5lzyhi1WkRE2pCT5OLuF7VXb2bzgEuB2aGrC+DPgCfc/TCw3cx+C1QCz/JO1xlhe0vYrgYmANWhO20IUBPbhYiISIck3i1mZnOBG4HL3L02q+otYJZFyoGzgbXuvhXYb2Znh/GUa4FHwzmPAfPC9pXA0qxkJSIiOZKTlstR3AWUAakw9r7c3RcA/w78CHiNqLvrR+7+SjjnBuDHRAP9i3hnnOYe4H4zqyJqsVyTo2sQEZEsiSeXMG24rfIDRNOR26pbAUxvo7zuSOeIiEjuJN4tJiIifY+Si4iIxE7JRUREYqfkIiIisTPN1AUz2wG82cXTRwI7YwwnSbqW/NNXrgN0LfmqO9dynLuPaqtCyaWbzGyFu1cmHUccdC35p69cB+ha8lVPXYu6xUREJHZKLiIiEjsll+5bmHQAMdK15J++ch2ga8lXPXItGnMREZHYqeUiIiKxU3IREZHYKbl0g5nNNbN1ZlZlZl9JOp6uMrMfmtl2M3st6Vi6w8wmmNlvzGyNma02s79KOqauMrP+ZvY7M1sVruWWpGPqDjMrNrOXzOyXScfSXWa20cxeNbOXzWxF0vF0lZkNNbOHzWxt+J35YKzvrzGXrjGzYuB1YA7RQ8peAD7m7ulEA+sCMzsfOADc5+7vWW26tzCzscBYd3/RzAYDK4Ereul/EwPK3f2AmfUDngP+yt2XJxxal5jZF4ke9lfh7pcmHU93mNlGoNLde/VNlGZ2L/Csu//AzEqBge6+J673V8ul684Eqtx9g7s3AD8FLk84pi5x92foA0/sdPet7v5i2N4PrAHGJRtV13jkQNjtF1698pugmY0HLgF+kHQsEjGzCuB8omdg4e4NcSYWUHLpjnHApqz9anrpH7K+yMwmAe8Dnk84lC4LXUkvA9uBlLv31mv5DvBloDnhOOLiwGIzW2lm85MOpotOAHYAPwrdlT8IT/yNjZJL11kbZb3ym2VfY2aDgJ8BX3D3fUnH01Xu3uTuZwDjgTPNrNd1WZrZpcB2d1+ZdCwxOsfd3w98BPhs6FbubUqA9wN3u/v7gINArOPGSi5dVw1MyNofD2xJKBYJwvjEz4AH3P1/ko4nDqG74ilgbrKRdMk5wGVhnOKnwCwz+0myIXWPu28JP7cDjxB1kfc21UB1Vmv4YaJkExsll657AZhiZseHwbBrgMcSjqmghUHwe4A17v7tpOPpDjMbZWZDw/YA4CJgbaJBdYG73+Tu4919EtHvyFJ3/4uEw+oyMysPk0UI3UgXA71ulqW7bwM2mdnJoWg2EOvEl5I436yQuHujmX0OeBIoBn7o7qsTDqtLzOy/gAuBkWZWDfyju9+TbFRdcg7wceDVMFYB8Hfu/nhyIXXZWODeMCuxCHjI3Xv9NN4+YAzwSPQ9hhLg/7n7E8mG1GX/B3ggfDneAHwyzjfXVGQREYmdusVERCR2Si4iIhI7JRcREYmdkouIiMROyUVERGKn5CLSDWHF4gsT+uyJZnYgTFcWySuaiiwSAzP7GjC5J28QDHe5f8rdl/TUZ4jERS0XkTxgZrqhWfoUJReRbggPjroU+Dvg6tBNtSrUDTGze8xsq5ltNrN/aunCMrNPmNlvzewOM6sBvmZmJ5rZUjPbZWY7zeyBrCVg7gcmAr8In/FlM5tkZt6SmMzsWDN7zMxqwgPsPp0V59fM7CEzu8/M9ofuvMrc/mtJIVFyEem+OuCbwIPuPsjdTw/l9wKNwGSi5f8vBj6Vdd5ZRMtujAb+mWil7W8BxwLTiBZG/RqAu38ceAv4aPiM29qI47+IFiQ8FrgS+KaZzc6qv4xo8cihROvg3dWdixZpj5KLSA8wszFES7J/wd0PhhV07yBavLHFFnf/rrs3uvshd69y95S717v7DuDbwAUd/LwJwLnAje5e5+4vEz2c6+NZhz3n7o+7exNwP3D6e99JJB7q5xXpGccRPT1ya1jkEKIvc9kPmMvexsxGA3cC5wGDw/G7O/h5xwI14QmcLd4kerRwi21Z27VAfzMrcffGDn6GSIep5SISj9bTLjcB9cBIdx8aXhXufmo753wrlM1w9wrgL3j3Q+nam9q5BRjeshx8MBHY3JmLEImLkotIPDLAJDMrAnD3rcBi4P+aWYWZFYUB+/a6uQYDB4A9ZjYO+Ns2PuOEtk50903A/wLfMrP+ZjYDuB54oFtXJdJFSi4i8fjv8HOXmb0Ytq8FSokewrSb6Gl/Y9t5j1uInga4F/gV0PpJmt8C/t7M9pjZl9o4/2PAJKJWzCNEz+VJdf5SRLpPN1GKiEjs1HIREZHYKbmIiEjslFxERCR2Si4iIhI7JRcREYmdkouIiMROyUVERGKn5CIiIrH7/7k7FXOcT7kDAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# red dashes, blue squares and green triangles\n",
    "plt.plot(results_by_iter['iter'], results_by_iter['LL'])\n",
    "plt.xlabel('iteration', fontsize=12)\n",
    "plt.ylabel('LogLikelihood', fontsize=12)\n",
    "plt.show()"
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}