Note
Go to the end to download the full example code.
Simulation of a cross-nested logit modelΒΆ
Illustration of the application of an estimated 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/b11cnl.yaml')
except FileNotFoundError:
print(
'Run first the script b11cnl.py in order to generate the file b11cnl.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.308242 ... -1.541166 1.232764e-01
1 b_time_train -1.073895 ... -7.579577 3.463896e-14
2 b_cost -0.973727 ... -14.711113 0.000000e+00
3 b_headway_train -0.004366 ... -4.491978 7.056485e-06
4 ga_train 1.142963 ... 4.935557 7.992239e-07
5 alpha_existing 0.644541 ... 3.746332 1.794391e-04
6 existing_nest_parameter 1.771197 ... 7.697990 1.376677e-14
7 b_time_swissmetro -0.991554 ... -5.574492 2.482537e-08
8 b_headway_swissmetro -0.007723 ... -2.601120 9.292006e-03
9 ga_swissmetro -0.138638 ... -0.860222 3.896665e-01
10 asc_car -0.606253 ... -4.886201 1.028002e-06
11 b_time_car -0.857033 ... -6.760819 1.372147e-11
12 public_nest_parameter 1.840637 ... 3.956373 7.609646e-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,
/Users/bierlair/python_envs/venv313/lib/python3.13/site-packages/scipy/integrate/_quadpack_py.py:1260: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
If increasing the limit yields no improvement it is advised to analyze
the integrand in order to determine the difficulties. If the position of a
local difficulty can be determined (singularity, discontinuity) one will
probably gain from splitting up the interval and calling the integrator
on the subranges. Perhaps a special-purpose integrator should be used.
quad_r = quad(f, low, high, args=args, full_output=self.full_output,
Train Swissmetro Car
Train 1.000000 3.827456e-01 5.308235e-01
Swissmetro 0.382746 1.000000e+00 8.253506e-12
Car 0.530823 8.253506e-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 = 'b11cnl_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.079453 0.683679 0.236869 -1.864135 -0.205801 -0.938887
1 0.160898 0.692937 0.146165 -1.425184 -0.197440 -1.234250
2 0.091761 0.615750 0.292490 -2.170137 -0.268820 -0.869686
3 0.127237 0.552379 0.320384 -1.659511 -0.302277 -0.534566
4 0.094348 0.685563 0.220089 -2.099597 -0.206042 -0.762397
... ... ... ... ... ... ...
6763 0.180536 0.704441 0.115023 -1.398696 -0.159526 -1.483052
6764 0.144838 0.708049 0.147113 -1.531130 -0.164446 -0.829455
6765 0.108647 0.696163 0.195190 -1.682743 -0.158975 -0.724200
6766 0.133255 0.725392 0.141353 -1.835460 -0.153694 -0.828747
6767 0.154699 0.697214 0.148087 -1.510788 -0.171445 -1.044934
[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: (1 minutes 11.873 seconds)