Simulation of a logit model

Example of simulation with a logit model

author:

Michel Bierlaire, EPFL

date:

Sun Apr 9 17:13:23 2023

import biogeme.biogeme as bio
from biogeme import models
from biogeme.expressions import Beta, Derive

See the data processing script: Data preparation for Swissmetro.

from swissmetro_data import (
    database,
    SM_AV,
    CAR_AV_SP,
    TRAIN_AV_SP,
    TRAIN_TT,
    TRAIN_TT_SCALED,
    TRAIN_COST_SCALED,
    SM_TT,
    SM_TT_SCALED,
    SM_COST_SCALED,
    CAR_TT,
    CAR_TT_SCALED,
    CAR_CO_SCALED,
)

Parameters.

ASC_CAR = Beta('ASC_CAR', 0, None, None, 0)
ASC_TRAIN = Beta('ASC_TRAIN', 0, None, None, 0)
ASC_SM = Beta('ASC_SM', 0, None, None, 1)
B_TIME = Beta('B_TIME', 0, None, None, 0)
B_COST = Beta('B_COST', 0, None, None, 0)

Definition of the utility functions.

V1 = ASC_TRAIN + B_TIME * TRAIN_TT_SCALED + B_COST * TRAIN_COST_SCALED
V2 = ASC_SM + B_TIME * SM_TT_SCALED + B_COST * SM_COST_SCALED
V3 = ASC_CAR + B_TIME * CAR_TT_SCALED + B_COST * CAR_CO_SCALED

Associate utility functions with the numbering of alternatives.

V = {1: V1, 2: V2, 3: V3}

Associate the availability conditions with the alternatives.

av = {1: TRAIN_AV_SP, 2: SM_AV, 3: CAR_AV_SP}

Choice probability.

prob1 = models.logit(V, av, 1)
prob2 = models.logit(V, av, 2)
prob3 = models.logit(V, av, 3)

Elasticities.

Elasticities can be computed. We illustrate below two formulas. Check in the output file that they produce the same result.

First, the general definition of elasticities. This illustrates the use of the Derive expression, and can be used with any model, however complicated it is. Note the quotes in the Derive opertor.

genelas1 = Derive(prob1, 'TRAIN_TT') * TRAIN_TT / prob1
genelas2 = Derive(prob2, 'SM_TT') * SM_TT / prob2
genelas3 = Derive(prob3, 'CAR_TT') * CAR_TT / prob3

Second, the elasticity of logit models. See Ben-Akiva and Lerman for the formula

logitelas1 = TRAIN_AV_SP * (1.0 - prob1) * TRAIN_TT_SCALED * B_TIME
logitelas2 = SM_AV * (1.0 - prob2) * SM_TT_SCALED * B_TIME
logitelas3 = CAR_AV_SP * (1.0 - prob3) * CAR_TT_SCALED * B_TIME

Quantities to be simulated.

simulate = {
    'Prob. train': prob1,
    'Prob. Swissmetro': prob2,
    'Prob. car': prob3,
    'logit elas. 1': logitelas1,
    'generic elas. 1': genelas1,
    'logit elas. 2': logitelas2,
    'generic elas. 2': genelas2,
    'logit elas. 3': logitelas3,
    'generic elas. 3': genelas3,
}

Create the Biogeme object.

As we simulate the probability for all aternatives, even when one of them is not available, Biogeme may trigger some warnings.

biosim = bio.BIOGEME(database, simulate)
biosim.modelName = 'b01logit_simul'

Values of the parameters.

betas = {
    'ASC_TRAIN': -0.701188,
    'B_TIME': -1.27786,
    'B_COST': -1.08379,
    'ASC_CAR': -0.154633,
}

Simulation

results = biosim.simulate(the_beta_values=betas)
results.describe()

	Prob. train	Prob. Swissmetro	Prob. car	logit elas. 1	generic elas. 1	logit elas. 2	generic elas. 2	logit elas. 3	generic elas. 3
count	6768.000000	6768.000000	6768.000000	6768.000000	6768.0	6768.000000	6768.0	6768.000000	6768.0
mean	0.134161	0.604315	0.261525	-1.872612	0.0	-0.447850	0.0	-1.136701	0.0
std	0.060525	0.182309	0.203746	0.886439	0.0	0.498262	0.0	1.074337	0.0
min	0.000010	0.000281	0.000000	-13.059604	0.0	-10.168912	0.0	-19.934616	0.0
25%	0.093416	0.487012	0.095533	-2.377303	0.0	-0.553065	0.0	-1.624370	0.0
50%	0.128314	0.611886	0.247461	-1.769928	0.0	-0.334896	0.0	-0.949290	0.0
75%	0.166545	0.753873	0.384113	-1.196807	0.0	-0.203747	0.0	-0.462240	0.0
max	0.642731	0.971398	0.999710	-0.232835	0.0	-0.024123	0.0	0.000000	0.0