{ "cells": [ { "cell_type": "code", "metadata": {}, "source": [ "#%%\n", "\"\"\"File 09nested.py\n", "\n", ":author: Michel Bierlaire, EPFL\n", ":date: Sun Sep 8 00:36:04 2019\n", "\n", " Example of a 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 pandas as pd\n", "import biogeme.database as db\n", "import biogeme.biogeme as bio\n", "import biogeme.optimization as opt\n", "import biogeme.models as models\n", "import biogeme.messaging as msg\n", "from biogeme.expressions import Beta, DefineVariable\n", "\n", "# Read the data\n", "df = pd.read_csv('swissmetro.dat', '\\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", "# Parameters to be estimated\n", "ASC_CAR = Beta('ASC_CAR', 0, None, None, 0)\n", "ASC_TRAIN = Beta('ASC_TRAIN', 0, None, None, 0)\n", "ASC_SM = Beta('ASC_SM', 0, None, None, 1)\n", "B_TIME = Beta('B_TIME', 0, None, None, 0)\n", "B_COST = Beta('B_COST', 0, None, None, 0)\n", "MU = Beta('MU', 1, 1, 10, 0)\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 + \\\n", " B_TIME * TRAIN_TT_SCALED + \\\n", " B_COST * TRAIN_COST_SCALED\n", "V2 = ASC_SM + \\\n", " B_TIME * SM_TT_SCALED + \\\n", " B_COST * SM_COST_SCALED\n", "V3 = ASC_CAR + \\\n", " B_TIME * CAR_TT_SCALED + \\\n", " B_COST * CAR_CO_SCALED\n", "\n", "# Associate utility functions with the numbering of alternatives\n", "V = {1: V1,\n", " 2: V2,\n", " 3: V3}\n", "\n", "# Associate the availability conditions with the alternatives\n", "av = {1: TRAIN_AV_SP,\n", " 2: SM_AV,\n", " 3: CAR_AV_SP}\n", "\n", "#Definition of nests:\n", "# 1: nests parameter\n", "# 2: list of alternatives\n", "existing = MU, [1, 3]\n", "future = 1.0, [2]\n", "nests = existing, future\n", "\n", "# Definition of the model. This is the contribution of each\n", "# observation to the log likelihood function.\n", "# The choice model is a nested logit, with availability conditions\n", "logprob = models.lognested(V, av, nests, CHOICE)\n", "\n", "# Define level of verbosity\n", "logger = msg.bioMessage()\n", "#logger.setSilent()\n", "#logger.setWarning()\n", "logger.setGeneral()\n", "#logger.setDetailed()\n", "\n", "# Create the Biogeme object\n", "biogeme = bio.BIOGEME(database, logprob)\n", "\n", "algos = {'CFSQP ': None,\n", " 'scipy ': opt.scipy,\n", " 'Simple bounds Newton ': opt.simpleBoundsNewtonAlgorithmForBiogeme,\n", " 'Simple bounds BFGS ': opt.simpleBoundsNewtonAlgorithmForBiogeme,\n", " 'Simple bounds hybrid 20%': opt.simpleBoundsNewtonAlgorithmForBiogeme,\n", " 'Simple bounds hybrid 50%': opt.simpleBoundsNewtonAlgorithmForBiogeme,\n", " 'Simple bounds hybrid 80%': opt.simpleBoundsNewtonAlgorithmForBiogeme}\n", "\n", "algoParameters = {'Simple bounds Newton': {'proportionAnalyticalHessian': 1.0},\n", " 'Simple bounds BFGS': {'proportionAnalyticalHessian': 0.0},\n", " 'Simple bounds hybrid 20%': {'proportionAnalyticalHessian': 0.2},\n", " 'Simple bounds hybrid 50%': {'proportionAnalyticalHessian': 0.5},\n", " 'Simple bounds hybrid 80%': {'proportionAnalyticalHessian': 0.8}}\n", "\n", "results = {}\n", "msg = ''\n", "for name, algo in algos.items():\n", " biogeme.modelName = f'09nested_allAlgos_{name}'.strip()\n", " p = algoParameters.get(name)\n", " results[name] = biogeme.estimate(algorithm=algo, algoParameters=p)\n", " g = results[name].data.g\n", " msg += (f'{name}\\t{results[name].data.logLike:.2f}\\t'\n", " f'{results[name].data.gradientNorm:.2g}\\t'\n", " f'{results[name].data.optimizationMessages[\"Optimization time\"]}'\n", " f'\\t{results[name].data.optimizationMessages[\"Cause of termination\"]}\\n')\n", "\n", "print(\"Algorithm\\t\\t\\tloglike\\t\\tnormg\\ttime\\t\\tdiagnostic\")\n", "print(\"+++++++++\\t\\t\\t+++++++\\t\\t+++++\\t++++\\t\\t++++++++++\")\n", "print(msg)\n", "\n", "\"\"\"\n", "Here are the results.\n", "\n", "Algorithm\t\t loglike\t\tnormg\ttime\t\tdiagnostic\n", "+++++++++\t\t +++++++\t\t+++++\t++++\t\t++++++++++\n", "CFSQP \t-5236.90\t0.00036\t0:00:01.342914\tb'Normal termination. Obj: 6.05545e-06 Const: 6.05545e-06'\n", "scipy \t-5236.90\t0.00017\t0:00:00.549708\tb'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'\n", "Simple bounds Newton \t-5236.90\t0.0027\t0:00:00.412040\tRelative gradient = 3.5e-07 <= 6.1e-06\n", "Simple bounds BFGS \t-5236.90\t0.0027\t0:00:00.403998\tRelative gradient = 3.5e-07 <= 6.1e-06\n", "Simple bounds hybrid 20%\t-5236.90\t0.024\t0:00:00.661252\tRelative gradient = 2.7e-06 <= 6.1e-06\n", "Simple bounds hybrid 50%\t-5236.90\t0.0074\t0:00:00.531165\tRelative gradient = 1.4e-06 <= 6.1e-06\n", "Simple bounds hybrid 80%\t-5236.90\t0.0096\t0:00:00.479969\tRelative gradient = 1.9e-06 <= 6.1e-06\n", "\"\"\"\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 }