Note
Go to the end to download the full example code.
11b. Simulation of a cross-nested logit modelΒΆ
Illustration of the application of an estimated CNL model.
Michel Bierlaire, EPFL Sat Jun 21 2025, 16:53:57
import sys
from IPython.core.display_functions import display
from biogeme.biogeme import BIOGEME
from biogeme.expressions import Beta, Derive
from biogeme.models import cnl
from biogeme.nests import NestsForCrossNestedLogit, OneNestForCrossNestedLogit
from biogeme.results_processing import (
EstimationResults,
get_pandas_estimated_parameters,
)
See the data processing script: Data preparation for Swissmetro.
from swissmetro_data import (
CAR_AV_SP,
CAR_CO_SCALED,
CAR_TT,
GA,
SM_AV,
SM_COST_SCALED,
SM_HE,
SM_TT,
TRAIN_AV_SP,
TRAIN_COST_SCALED,
TRAIN_HE,
TRAIN_TT,
database,
)
Parameters to be estimated.
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_swissmetro = Beta('b_time_swissmetro', 0, None, None, 0)
b_time_train = Beta('b_time_train', 0, None, None, 0)
b_time_car = Beta('b_time_car', 0, None, None, 0)
b_cost = Beta('b_cost', 0, None, None, 0)
b_headway_swissmetro = Beta('b_headway_swissmetro', 0, None, None, 0)
b_headway_train = Beta('b_headway_train', 0, None, None, 0)
ga_train = Beta('ga_train', 0, None, None, 0)
ga_swissmetro = Beta('ga_swissmetro', 0, None, None, 0)
existing_nest_parameter = Beta('existing_nest_parameter', 1, 1, 5, 0)
public_nest_parameter = Beta('public_nest_parameter', 1, 1, 5, 0)
Nest membership parameters.
alpha_existing = Beta('alpha_existing', 0.5, 0, 1, 0)
alpha_public = 1 - alpha_existing
Definition of the utility functions. Note that in order to calculate the derivative with respect to the travel time variables, they need to explicitly appear in the specification. Therefore, we have replaced the scaled versions of the variables by their original definition.
v_train = (
asc_train
+ b_time_train * TRAIN_TT / 100
+ b_cost * TRAIN_COST_SCALED
+ b_headway_train * TRAIN_HE
+ ga_train * GA
)
v_swissmetro = (
asc_sm
+ b_time_swissmetro * SM_TT / 100
+ b_cost * SM_COST_SCALED
+ b_headway_swissmetro * SM_HE
+ ga_swissmetro * GA
)
v_car = asc_car + b_time_car * 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}
Definition of nests.
nest_existing = OneNestForCrossNestedLogit(
nest_param=existing_nest_parameter,
dict_of_alpha={1: alpha_existing, 2: 0.0, 3: 1.0},
name='existing',
)
nest_public = OneNestForCrossNestedLogit(
nest_param=public_nest_parameter,
dict_of_alpha={1: alpha_public, 2: 1.0, 3: 0.0},
name='public',
)
nests = NestsForCrossNestedLogit(
choice_set=[1, 2, 3], tuple_of_nests=(nest_existing, nest_public)
)
Read the estimation results from the pickle file.
try:
results = EstimationResults.from_yaml_file(filename='saved_results/b11a_cnl.yaml')
except FileNotFoundError:
print(
'Run first the script plot_b11a_cnl.py in order to generate the file b11a_cnl.yaml, and move it to the directory '
'saved_results.'
)
sys.exit()
print(
'Estimation results: ', get_pandas_estimated_parameters(estimation_results=results)
)
Estimation results: Name Value ... Robust t-stat. Robust p-value
0 asc_train -0.308465 ... -1.541500 1.231951e-01
1 b_time_train -1.073896 ... -7.579054 3.486100e-14
2 b_cost -0.973726 ... -14.711229 0.000000e+00
3 b_headway_train -0.004366 ... -4.491560 7.070336e-06
4 ga_train 1.142978 ... 4.934247 8.046038e-07
5 alpha_existing 0.644707 ... 3.744259 1.809270e-04
6 existing_nest_parameter 1.771215 ... 7.695616 1.398881e-14
7 b_time_swissmetro -0.991511 ... -5.574074 2.488505e-08
8 b_headway_swissmetro -0.007724 ... -2.601172 9.290596e-03
9 ga_swissmetro -0.138849 ... -0.861538 3.889421e-01
10 asc_car -0.606265 ... -4.886097 1.028545e-06
11 b_time_car -0.857011 ... -6.760465 1.375500e-11
12 public_nest_parameter 1.839948 ... 3.953511 7.701258e-05
[13 rows x 5 columns]
print('Calculating correlation matrix. It may generate numerical warnings from scipy.')
corr = nests.correlation(
parameters=results.get_beta_values(),
alternatives_names={1: 'Train', 2: 'Swissmetro', 3: 'Car'},
)
display(corr)
Calculating correlation matrix. It may generate numerical warnings from scipy.
/Users/bierlair/python_envs/venv313/lib/python3.13/site-packages/scipy/integrate/_quadpack_py.py:1260: IntegrationWarning: The integral is probably divergent, or slowly convergent.
quad_r = quad(f, low, high, args=args, full_output=self.full_output,
/Users/bierlair/python_envs/venv313/lib/python3.13/site-packages/scipy/integrate/_quadpack_py.py:1260: IntegrationWarning: The algorithm does not converge. Roundoff error is detected
in the extrapolation table. It is assumed that the requested tolerance
cannot be achieved, and that the returned result (if full_output = 1) is
the best which can be obtained.
quad_r = quad(f, low, high, args=args, full_output=self.full_output,
Train Swissmetro Car
Train 1.000000 3.825254e-01 5.309075e-01
Swissmetro 0.382525 1.000000e+00 8.286443e-12
Car 0.530907 8.286443e-12 1.000000e+00
The choice model is a cross-nested logit, with availability conditions.
probability_train = cnl(v, av, nests, 1)
probability_swissmetro = cnl(v, av, nests, 2)
probability_car = cnl(v, av, nests, 3)
We calculate elasticities. It is important that the variables explicitly appear as such in the model. If not, the derivative will be zero, as well as the elasticities.
general_time_eslaticity_train = (
Derive(probability_train, 'TRAIN_TT') * TRAIN_TT / probability_train
)
general_time_elasticity_swissmetro = (
Derive(probability_swissmetro, 'SM_TT') * SM_TT / probability_swissmetro
)
general_time_elasticity_car = (
Derive(probability_car, 'CAR_TT') * CAR_TT / probability_car
)
We report the probability of each alternative and the elasticities.
simulate = {
'Prob. train': probability_train,
'Prob. Swissmetro': probability_swissmetro,
'Prob. car': probability_car,
'Elas. 1': general_time_eslaticity_train,
'Elas. 2': general_time_elasticity_swissmetro,
'Elas. 3': general_time_elasticity_car,
}
Create the Biogeme object.
biosim = BIOGEME(database, simulate)
biosim.model_name = 'b11b_cnl_simul'
Perform the simulation.
simulation_results = biosim.simulate(results.get_beta_values())
print('Simulation results')
display(simulation_results)
Simulation results
Prob. train Prob. Swissmetro Prob. car Elas. 1 Elas. 2 Elas. 3
0 0.079452 0.683679 0.236869 -1.864056 -0.205787 -0.938859
1 0.160898 0.692939 0.146163 -1.425095 -0.197418 -1.234231
2 0.091760 0.615749 0.292490 -2.170005 -0.268798 -0.869664
3 0.127235 0.552385 0.320380 -1.659378 -0.302238 -0.534559
4 0.094349 0.685563 0.220088 -2.099494 -0.206026 -0.762381
... ... ... ... ... ... ...
6763 0.180538 0.704439 0.115023 -1.398603 -0.159509 -1.483026
6764 0.144840 0.708050 0.147110 -1.531040 -0.164430 -0.829444
6765 0.108649 0.696162 0.195189 -1.682657 -0.158963 -0.724187
6766 0.133259 0.725389 0.141352 -1.835360 -0.153682 -0.828735
6767 0.154699 0.697214 0.148087 -1.510699 -0.171428 -1.044914
[6768 rows x 6 columns]
print(
f'Aggregate share of train: {100 * simulation_results["Prob. train"].mean():.1f}%'
)
Aggregate share of train: 13.3%
Total running time of the script: (0 minutes 32.069 seconds)