.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/swissmetro/plot_b11b_cnl_simul.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_swissmetro_plot_b11b_cnl_simul.py: 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 .. GENERATED FROM PYTHON SOURCE LINES 11-24 .. code-block:: Python 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, ) .. GENERATED FROM PYTHON SOURCE LINES 25-26 See the data processing script: :ref:`swissmetro_data`. .. GENERATED FROM PYTHON SOURCE LINES 26-42 .. code-block:: Python 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, ) .. GENERATED FROM PYTHON SOURCE LINES 43-44 Parameters to be estimated. .. GENERATED FROM PYTHON SOURCE LINES 44-56 .. code-block:: Python 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) .. GENERATED FROM PYTHON SOURCE LINES 57-60 .. code-block:: Python existing_nest_parameter = Beta('existing_nest_parameter', 1, 1, 5, 0) public_nest_parameter = Beta('public_nest_parameter', 1, 1, 5, 0) .. GENERATED FROM PYTHON SOURCE LINES 61-62 Nest membership parameters. .. GENERATED FROM PYTHON SOURCE LINES 62-65 .. code-block:: Python alpha_existing = Beta('alpha_existing', 0.5, 0, 1, 0) alpha_public = 1 - alpha_existing .. GENERATED FROM PYTHON SOURCE LINES 66-69 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. .. GENERATED FROM PYTHON SOURCE LINES 69-85 .. code-block:: Python 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 .. GENERATED FROM PYTHON SOURCE LINES 86-87 Associate utility functions with the numbering of alternatives .. GENERATED FROM PYTHON SOURCE LINES 87-89 .. code-block:: Python v = {1: v_train, 2: v_swissmetro, 3: v_car} .. GENERATED FROM PYTHON SOURCE LINES 90-91 Associate the availability conditions with the alternatives .. GENERATED FROM PYTHON SOURCE LINES 91-93 .. code-block:: Python av = {1: TRAIN_AV_SP, 2: SM_AV, 3: CAR_AV_SP} .. GENERATED FROM PYTHON SOURCE LINES 94-95 Definition of nests. .. GENERATED FROM PYTHON SOURCE LINES 95-112 .. code-block:: Python 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) ) .. GENERATED FROM PYTHON SOURCE LINES 113-114 Read the estimation results from the pickle file. .. GENERATED FROM PYTHON SOURCE LINES 114-124 .. code-block:: Python 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() .. GENERATED FROM PYTHON SOURCE LINES 125-129 .. code-block:: Python print( 'Estimation results: ', get_pandas_estimated_parameters(estimation_results=results) ) .. rst-class:: sphx-glr-script-out .. code-block:: none 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] .. GENERATED FROM PYTHON SOURCE LINES 130-137 .. code-block:: Python 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) .. rst-class:: sphx-glr-script-out .. code-block:: none 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 .. GENERATED FROM PYTHON SOURCE LINES 138-139 The choice model is a cross-nested logit, with availability conditions. .. GENERATED FROM PYTHON SOURCE LINES 139-143 .. code-block:: Python probability_train = cnl(v, av, nests, 1) probability_swissmetro = cnl(v, av, nests, 2) probability_car = cnl(v, av, nests, 3) .. GENERATED FROM PYTHON SOURCE LINES 144-147 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. .. GENERATED FROM PYTHON SOURCE LINES 147-157 .. code-block:: Python 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 ) .. GENERATED FROM PYTHON SOURCE LINES 158-159 We report the probability of each alternative and the elasticities. .. GENERATED FROM PYTHON SOURCE LINES 159-168 .. code-block:: Python 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, } .. GENERATED FROM PYTHON SOURCE LINES 169-170 Create the Biogeme object. .. GENERATED FROM PYTHON SOURCE LINES 170-173 .. code-block:: Python biosim = BIOGEME(database, simulate) biosim.model_name = 'b11b_cnl_simul' .. GENERATED FROM PYTHON SOURCE LINES 174-175 Perform the simulation. .. GENERATED FROM PYTHON SOURCE LINES 175-177 .. code-block:: Python simulation_results = biosim.simulate(results.get_beta_values()) .. GENERATED FROM PYTHON SOURCE LINES 178-181 .. code-block:: Python print('Simulation results') display(simulation_results) .. rst-class:: sphx-glr-script-out .. code-block:: none 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] .. GENERATED FROM PYTHON SOURCE LINES 182-185 .. code-block:: Python print( f'Aggregate share of train: {100 * simulation_results["Prob. train"].mean():.1f}%' ) .. rst-class:: sphx-glr-script-out .. code-block:: none Aggregate share of train: 13.3% .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 32.069 seconds) .. _sphx_glr_download_auto_examples_swissmetro_plot_b11b_cnl_simul.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_b11b_cnl_simul.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_b11b_cnl_simul.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_b11b_cnl_simul.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_