Note
Go to the end to download the full example code.
1c. Simulation of a logit model (traditional and Bayesian)¶
Example of simulation with a logit model
Michel Bierlaire, EPFL Thu Oct 30 2025, 14:03:15
import sys
import pandas as pd
from IPython.core.display_functions import display
See the data processing script: Data preparation for Swissmetro.
from swissmetro_data import (
CAR_AV_SP,
CAR_CO_SCALED,
CAR_TT,
SM_AV,
SM_COST_SCALED,
SM_TT,
TRAIN_AV_SP,
TRAIN_COST_SCALED,
TRAIN_TT,
database,
)
import biogeme.biogeme_logging as blog
from biogeme.bayesian_estimation import BayesianResults
from biogeme.biogeme import BIOGEME
from biogeme.expressions import Beta, Derive
from biogeme.models import logit
logger = blog.get_screen_logger(level=blog.INFO)
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. As we will calculate the derivative with respect to TRAIN_TT, SM_TT and CAR_TT, they must explicitly appear in the model. If not, the derivative will be zero. Therefore, we do not use the _SCALED version of the attributes. We explicitly include their definition.
v_train = asc_train + b_time * TRAIN_TT / 100 + b_cost * TRAIN_COST_SCALED
v_swissmetro = asc_sm + b_time * SM_TT / 100 + b_cost * SM_COST_SCALED
v_car = asc_car + b_time * CAR_TT / 100 + b_cost * CAR_CO_SCALED
Associate utility functions with the numbering of alternatives.
v = {1: v_train, 2: v_swissmetro, 3: v_car}
Associate the availability conditions with the alternatives.
av = {1: TRAIN_AV_SP, 2: SM_AV, 3: CAR_AV_SP}
Choice probability.
prob_train = logit(v, av, 1)
Elasticity.
time_elasticity_train = Derive(prob_train, 'TRAIN_TT') * TRAIN_TT / prob_train
Quantities to be simulated.
simulate = {
'Prob. train': prob_train,
'train time elasticity': time_elasticity_train,
'Value of time': b_time / b_cost,
}
Create the Biogeme object.
As we simulate the probability for all alternatives, even when one of them is not available, Biogeme may trigger some warnings.
biosim = BIOGEME(database, simulate)
biosim.model_name = 'b01c_logit_simul'
Biogeme parameters read from biogeme.toml.
Retrieve the estimated values of the parameters.
RESULTS_FILE_NAME = 'saved_results/b01a_logit.nc'
try:
estimation_results = BayesianResults.from_netcdf(filename=RESULTS_FILE_NAME)
except FileNotFoundError:
logger.error(
f'File {RESULTS_FILE_NAME} does not exist. Run the estimation script first.'
)
sys.exit()
betas = estimation_results.get_beta_values()
Loaded NetCDF file size: 930.9 MB
load finished in 50 ms
Diagnostics computation took 21.1 seconds (cached).
Simulation using the posterior mean of each parameter
print('Simulation using the posterior mean of each parameter')
results = biosim.simulate(the_beta_values=betas)
display(results)
Simulation using the posterior mean of each parameter
Prob. train train time elasticity Value of time
0 0.167836 -1.193375 1.181095
1 0.184112 -1.076014 1.181095
2 0.142833 -1.426787 1.181095
3 0.161131 -1.106322 1.181095
4 0.139621 -1.432134 1.181095
... ... ... ...
6763 0.172374 -1.144480 1.181095
6764 0.164425 -1.155472 1.181095
6765 0.149244 -1.176465 1.181095
6766 0.134525 -1.418452 1.181095
6767 0.176708 -1.138487 1.181095
[6768 rows x 3 columns]
Bayesian simulation using the posterior draws
print('Bayesian simulation')
bayesian_results = biosim.simulate_bayesian(
bayesian_estimation_results=estimation_results, percentage_of_draws_to_use=3
)
with pd.option_context('display.max_columns', None, 'display.expand_frame_repr', False):
display(bayesian_results)
Bayesian simulation
Bayesian simulation performed with 3% of the draws, that is 240/8000 draws. Adjust the parameter "percentage_of_draws_to_use" if you need a different number of draws.
0%| | 0/240 [00:00<?, ?it/s]
40%|████ | 97/240 [00:00<00:00, 965.73it/s]
85%|████████▌ | 204/240 [00:00<00:00, 1025.23it/s]
100%|██████████| 240/240 [00:00<00:00, 1026.31it/s]
Prob. train_mean Prob. train_q025 Prob. train_q975 train time elasticity_mean train time elasticity_q025 train time elasticity_q975 Value of time_mean Value of time_q025 Value of time_q975
__biogeme_internal_obs_id__
0 0.169049 0.159572 0.179921 -1.193815 -1.300611 -1.088641 1.191855 1.07219 1.333216
1 0.185445 0.174157 0.197510 -1.076168 -1.169746 -0.983543 1.191855 1.07219 1.333216
2 0.143830 0.135783 0.152670 -1.427853 -1.561303 -1.296890 1.191855 1.07219 1.333216
3 0.162177 0.153098 0.172204 -1.107015 -1.208190 -1.006578 1.191855 1.07219 1.333216
4 0.140498 0.133031 0.148981 -1.433463 -1.568961 -1.298703 1.191855 1.07219 1.333216
... ... ... ... ... ... ... ... ... ...
6763 0.173583 0.164135 0.184458 -1.144924 -1.248113 -1.042941 1.191855 1.07219 1.333216
6764 0.165470 0.156611 0.176083 -1.156199 -1.261435 -1.050582 1.191855 1.07219 1.333216
6765 0.150248 0.141944 0.159680 -1.177319 -1.286645 -1.069279 1.191855 1.07219 1.333216
6766 0.135366 0.128093 0.143513 -1.419873 -1.555795 -1.285192 1.191855 1.07219 1.333216
6767 0.177843 0.168219 0.188807 -1.139021 -1.240790 -1.036346 1.191855 1.07219 1.333216
[6768 rows x 9 columns]
Total running time of the script: (0 minutes 33.715 seconds)