{ "cells": [ { "cell_type": "code", "metadata": {}, "source": [ "#%%\n", "\"\"\"File 25triangularMixture.py\n", "\n", ":author: Michel Bierlaire, EPFL\n", ":date: Mon Sep 9 10:19:24 2019\n", "\n", " Example of a mixture of logit models, using Monte-Carlo integration.\n", " The mixing distirbution is specified by the user. Here, a triangular\n", " distribution.\n", " Three alternatives: Train, Car and Swissmetro\n", " SP data\n", "\n", "\"\"\"\n", "\n", "import numpy as np\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, log, MonteCarlo\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_COST = Beta('B_COST', 0, None, None, 0)\n", "\n", "# Define a random parameter, normally distributed, designed to be used\n", "# for Monte-Carlo simulation\n", "B_TIME = Beta('B_TIME', 0, 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, None, None, 0)\n", "\n", "def theTriangularGenerator(sampleSize, numberOfDraws):\n", " \"\"\"\n", " User-defined random number generator to the database.\n", " See the numpy.random documentation to obtain a list of other distributions.\n", " \"\"\"\n", " return np.random.triangular(-1, 0, 1, (sampleSize, numberOfDraws))\n", "\n", "myRandomNumberGenerators = {'TRIANGULAR': (theTriangularGenerator,\n", " 'Draws from a triangular distribution')}\n", "\n", "database.setRandomNumberGenerators(myRandomNumberGenerators)\n", "\n", "# Define a random parameter with a triangular distribution, designed to be used\n", "# for Monte-Carlo simulation\n", "B_TIME_RND = B_TIME + B_TIME_S * bioDraws('B_TIME_RND', 'TRIANGULAR')\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_RND * TRAIN_TT_SCALED + \\\n", " B_COST * TRAIN_COST_SCALED\n", "V2 = ASC_SM + \\\n", " B_TIME_RND * SM_TT_SCALED + \\\n", " B_COST * SM_COST_SCALED\n", "V3 = ASC_CAR + \\\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 B_TIME_RND, we have a logit model (called the kernel)\n", "prob = models.logit(V, av, CHOICE)\n", "\n", "# We integrate over B_TIME_RND using Monte-Carlo\n", "logprob = log(MonteCarlo(prob))\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, numberOfDraws=100000)\n", "biogeme.modelName = '25triangularMixture'\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 }