{ "cells": [ { "cell_type": "code", "metadata": {}, "source": [ "#%%\n", "\"\"\"File 11cnl_simul.py\n", "\n", ":author: Michel Bierlaire, EPFL\n", ":date: Sun Sep 8 11:13:22 2019\n", "\n", " Example of simulation with a cross-nested logit model.\n", " Three alternatives: Train, Car and Swissmetro\n", " Train and car are in the same nest.\n", " SP data\n", "\"\"\"\n", "\n", "import sys\n", "import pandas as pd\n", "import biogeme.database as db\n", "import biogeme.biogeme as bio\n", "from biogeme import models\n", "import biogeme.results as res\n", "from biogeme.expressions import Beta, Derive\n", "\n", "# Read the data\n", "df = pd.read_csv('swissmetro.dat', sep='\\t')\n", "database = db.Database('swissmetro', df)\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", "# Simulation should be done with estimated value of the\n", "# parameters. You can include them manually. Here, we prefer to set\n", "# them to zero, and read the values from the file created after\n", "# estimation.\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", "MU_EXISTING = Beta('MU_EXISTING', 1, 1, None, 0)\n", "MU_PUBLIC = Beta('MU_PUBLIC', 1, 1, None, 0)\n", "ALPHA_EXISTING = Beta('ALPHA_EXISTING', 0.5, 0, 1, 0)\n", "ALPHA_PUBLIC = 1 - ALPHA_EXISTING\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: in simulation, do not use the\n", "# DefineVariable operator, as it hides the functional\n", "# relationships. In particular, derivatives cannot be calculated.\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 + B_TIME * TRAIN_TT_SCALED + B_COST * TRAIN_COST_SCALED\n", "V2 = ASC_SM + B_TIME * SM_TT_SCALED + B_COST * SM_COST_SCALED\n", "V3 = ASC_CAR + B_TIME * CAR_TT_SCALED + B_COST * CAR_CO_SCALED\n", "\n", "# Associate utility functions with the numbering of alternatives\n", "V = {1: V1, 2: V2, 3: V3}\n", "\n", "# Associate the availability conditions with the alternatives\n", "av = {1: TRAIN_AV_SP, 2: SM_AV, 3: CAR_AV_SP}\n", "# Definition of nests\n", "# Nest membership parameters\n", "alpha_existing = {1: ALPHA_EXISTING, 2: 0.0, 3: 1.0}\n", "\n", "alpha_public = {1: ALPHA_PUBLIC, 2: 1.0, 3: 0.0}\n", "\n", "nest_existing = MU_EXISTING, alpha_existing\n", "nest_public = MU_PUBLIC, alpha_public\n", "nests = nest_existing, nest_public\n", "\n", "# Instead of estimating the parameters, read the estimation\n", "# results from the pickle file.\n", "try:\n", " results = res.bioResults(pickleFile='11cnl.pickle')\n", "except FileNotFoundError:\n", " print(\n", " 'Run first the script 11cnl.py in order to generate the file '\n", " '11cnl.pickle.'\n", " )\n", " sys.exit()\n", "\n", "print('Estimation results: ', results.getEstimatedParameters())\n", "\n", "# The choice model is a cross-nested logit, with availability conditions\n", "prob1 = models.cnl_avail(V, av, nests, 1)\n", "prob2 = models.cnl_avail(V, av, nests, 2)\n", "prob3 = models.cnl_avail(V, av, nests, 3)\n", "\n", "# We calculate elasticities\n", "genelas1 = Derive(prob1, 'TRAIN_TT') * TRAIN_TT / prob1\n", "genelas2 = Derive(prob2, 'SM_TT') * SM_TT / prob2\n", "genelas3 = Derive(prob3, 'CAR_TT') * CAR_TT / prob3\n", "\n", "# We report the probability of each alternative and the elasticities\n", "simulate = {\n", " 'Prob. train': prob1,\n", " 'Prob. Swissmetro': prob2,\n", " 'Prob. car': prob3,\n", " 'Elas. 1': genelas1,\n", " 'Elas. 2': genelas2,\n", " 'Elas. 3': genelas3,\n", "}\n", "\n", "# Create the Biogeme object\n", "biosim = bio.BIOGEME(database, simulate)\n", "biosim.modelName = '11cnl_simul'\n", "\n", "# Perform the simulation\n", "simresults = biosim.simulate(results.getBetaValues())\n", "print('Simulation results')\n", "print(simresults.describe())\n", "\n", "print(f'Aggregate share of train: {100*simresults[\"Prob. train\"].mean():.1f}%')\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": 4 }