{ "cells": [ { "cell_type": "code", "metadata": {}, "source": [ "#%%\n", "\"\"\"File 13panel.py\n", "\n", ":author: Michel Bierlaire, EPFL\n", ":date: Sun Sep 8 18:55:38 2019\n", "\n", " Example of a mixture of logit models, using Monte-Carlo integration.\n", " The datafile is organized as panel data.\n", " Same as 12panel, where the starting values for the parameters are close to optimal.\n", " Useful for a faster estimation of the model with a large number of draws.\n", " Three alternatives: Train, Car and Swissmetro\n", " SP data\n", "\"\"\"\n", "\n", "import pandas as pd\n", "import biogeme.database as db\n", "import biogeme.biogeme as bio\n", "import biogeme.models as models\n", "import biogeme.messaging as msg\n", "from biogeme.expressions import Beta, DefineVariable, bioDraws, \\\n", " PanelLikelihoodTrajectory, MonteCarlo, log\n", "\n", "# Read the data\n", "df = pd.read_csv('swissmetro.dat', '\\t')\n", "database = db.Database('swissmetro', df)\n", "\n", "# They are organized as panel data. The variable ID identifies each individual.\n", "database.panel(\"ID\")\n", "\n", "# The Pandas data structure is available as database.data. Use all the\n", "# Pandas functions to invesigate the database\n", "#print(database.data.describe())\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 can be done directly using pandas.\n", "#remove = (((database.data.PURPOSE != 1) &\n", "# (database.data.PURPOSE != 3)) |\n", "# (database.data.CHOICE == 0))\n", "#database.data.drop(database.data[remove].index,inplace=True)\n", "\n", "# Here we use the \"biogeme\" way for backward compatibility\n", "exclude = ((PURPOSE != 1) * (PURPOSE != 3) + (CHOICE == 0)) > 0\n", "database.remove(exclude)\n", "\n", "\n", "# Parameters to be estimated\n", "B_COST = Beta('B_COST', -3.32, None, None, 0)\n", "\n", "# Define a random parameter, normally distributed across individuals,\n", "# designed to be used for Monte-Carlo simulation\n", "B_TIME = Beta('B_TIME', -5.4, None, None, 0)\n", "\n", "# It is advised not to use 0 as starting value for the following parameter.\n", "B_TIME_S = Beta('B_TIME_S', 1.55, None, None, 0)\n", "B_TIME_RND = B_TIME + B_TIME_S * bioDraws('B_TIME_RND', 'NORMAL_ANTI')\n", "\n", "# We do the same for the constants, to address serial correlation.\n", "ASC_CAR = Beta('ASC_CAR', 0.357, None, None, 0)\n", "ASC_CAR_S = Beta('ASC_CAR_S', 3.37, None, None, 0)\n", "ASC_CAR_RND = ASC_CAR + ASC_CAR_S * bioDraws('ASC_CAR_RND', 'NORMAL_ANTI')\n", "\n", "ASC_TRAIN = Beta('ASC_TRAIN', -0.617, None, None, 0)\n", "ASC_TRAIN_S = Beta('ASC_TRAIN_S', 3.13, None, None, 0)\n", "ASC_TRAIN_RND = ASC_TRAIN + ASC_TRAIN_S * bioDraws('ASC_TRAIN_RND', 'NORMAL_ANTI')\n", "\n", "ASC_SM = Beta('ASC_SM', 0, None, None, 1)\n", "ASC_SM_S = Beta('ASC_SM_S', 1.36, None, None, 0)\n", "ASC_SM_RND = ASC_SM + ASC_SM_S * bioDraws('ASC_SM_RND', 'NORMAL_ANTI')\n", "\n", "# Definition of new variables\n", "SM_COST = SM_CO * (GA == 0)\n", "TRAIN_COST = TRAIN_CO * (GA == 0)\n", "\n", "# Definition of new variables: adding columns to the database\n", "CAR_AV_SP = DefineVariable('CAR_AV_SP', CAR_AV * (SP != 0), database)\n", "TRAIN_AV_SP = DefineVariable('TRAIN_AV_SP', TRAIN_AV * (SP != 0), database)\n", "TRAIN_TT_SCALED = DefineVariable('TRAIN_TT_SCALED', TRAIN_TT / 100.0, database)\n", "TRAIN_COST_SCALED = DefineVariable('TRAIN_COST_SCALED', TRAIN_COST / 100, database)\n", "SM_TT_SCALED = DefineVariable('SM_TT_SCALED', SM_TT / 100.0, database)\n", "SM_COST_SCALED = DefineVariable('SM_COST_SCALED', SM_COST / 100, database)\n", "CAR_TT_SCALED = DefineVariable('CAR_TT_SCALED', CAR_TT / 100, database)\n", "CAR_CO_SCALED = DefineVariable('CAR_CO_SCALED', CAR_CO / 100, database)\n", "\n", "# Definition of the utility functions\n", "V1 = ASC_TRAIN_RND + \\\n", " B_TIME_RND * TRAIN_TT_SCALED + \\\n", " B_COST * TRAIN_COST_SCALED\n", "V2 = ASC_SM_RND + \\\n", " B_TIME_RND * SM_TT_SCALED + \\\n", " B_COST * SM_COST_SCALED\n", "V3 = ASC_CAR_RND + \\\n", " B_TIME_RND * 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", "# Conditional to the random parameters, the likelihood of one observation is\n", "# given by the logit model (called the kernel)\n", "obsprob = models.logit(V, av, CHOICE)\n", "\n", "# Conditional to the random parameters, the likelihood of all observations for\n", "# one individual (the trajectory) is the product of the likelihood of\n", "# each observation.\n", "condprobIndiv = PanelLikelihoodTrajectory(obsprob)\n", "\n", "# We integrate over the random parameters using Monte-Carlo\n", "logprob = log(MonteCarlo(condprobIndiv))\n", "\n", "# Define level of verbosity\n", "logger = msg.bioMessage()\n", "#logger.setSilent()\n", "#logger.setWarning()\n", "#logger.setGeneral()\n", "logger.setDetailed()\n", "#logger.setDebug()\n", "\n", "# Create the Biogeme object\n", "biogeme = bio.BIOGEME(database, logprob, numberOfDraws=100000)\n", "biogeme.modelName = '13panel'\n", "\n", "# Estimate the parameters.\n", "results = biogeme.estimate()\n", "pandasResults = results.getEstimatedParameters()\n", "print(pandasResults)\n" ], "outputs": [], "execution_count": null } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.1" } }, "nbformat": 4, "nbformat_minor": 1 }