.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/swissmetro/plot_b13panel_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_b13panel_simul.py: Simulation of panel model ========================= Calculates each contribution to the log likelihood function using simulation. We also calculate the individual parameters. Michel Bierlaire, EPFL Sat Jun 21 2025, 17:06:31 .. GENERATED FROM PYTHON SOURCE LINES 13-26 .. code-block:: Python import sys from IPython.core.display_functions import display import biogeme.biogeme_logging as blog from biogeme.biogeme import BIOGEME from biogeme.calculator.single_formula import calculate_single_formula_from_expression from biogeme.expressions import Beta, Draws, MonteCarlo, PanelLikelihoodTrajectory, log from biogeme.models import logit from biogeme.results_processing import EstimationResults from biogeme.second_derivatives import SecondDerivativesMode .. GENERATED FROM PYTHON SOURCE LINES 27-28 See the data processing script: :ref:`swissmetro_panel`. .. GENERATED FROM PYTHON SOURCE LINES 28-42 .. code-block:: Python from swissmetro_panel import ( CAR_AV_SP, CAR_CO_SCALED, CAR_TT_SCALED, CHOICE, SM_AV, SM_COST_SCALED, SM_TT_SCALED, TRAIN_AV_SP, TRAIN_COST_SCALED, TRAIN_TT_SCALED, database, ) .. GENERATED FROM PYTHON SOURCE LINES 43-44 Define the number of draws to be used for Monte-Carlo integration. .. GENERATED FROM PYTHON SOURCE LINES 44-49 .. code-block:: Python NUMBER_OF_DRAWS = 100_000 logger = blog.get_screen_logger(level=blog.INFO) logger.info('Example b13panel_simul.py') .. rst-class:: sphx-glr-script-out .. code-block:: none Example b13panel_simul.py .. GENERATED FROM PYTHON SOURCE LINES 50-51 Parameters to be estimated. .. GENERATED FROM PYTHON SOURCE LINES 51-53 .. code-block:: Python b_cost = Beta('b_cost', 0, None, None, 0) .. GENERATED FROM PYTHON SOURCE LINES 54-56 Define a random parameter, normally distributed across individuals, designed to be used for Monte-Carlo simulation. .. GENERATED FROM PYTHON SOURCE LINES 56-58 .. code-block:: Python b_time = Beta('b_time', 0, None, None, 0) .. GENERATED FROM PYTHON SOURCE LINES 59-60 It is advised not to use 0 as starting value for the following parameter. .. GENERATED FROM PYTHON SOURCE LINES 60-63 .. code-block:: Python b_time_s = Beta('b_time_s', 1, None, None, 0) b_time_rnd = b_time + b_time_s * Draws('b_time_rnd', 'NORMAL_ANTI') .. GENERATED FROM PYTHON SOURCE LINES 64-65 We do the same for the constants, to address serial correlation. .. GENERATED FROM PYTHON SOURCE LINES 65-77 .. code-block:: Python asc_car = Beta('asc_car', 0, None, None, 0) asc_car_s = Beta('asc_car_s', 1, None, None, 0) asc_car_rnd = asc_car + asc_car_s * Draws('asc_car_rnd', 'NORMAL_ANTI') asc_train = Beta('asc_train', 0, None, None, 0) asc_train_s = Beta('asc_train_s', 1, None, None, 0) asc_train_rnd = asc_train + asc_train_s * Draws('asc_train_rnd', 'NORMAL_ANTI') asc_sm = Beta('asc_sm', 0, None, None, 0) asc_sm_s = Beta('asc_sm_s', 1, None, None, 0) asc_sm_rnd = asc_sm + asc_sm_s * Draws('asc_sm_rnd', 'NORMAL_ANTI') .. GENERATED FROM PYTHON SOURCE LINES 78-79 Definition of the utility functions. .. GENERATED FROM PYTHON SOURCE LINES 79-83 .. code-block:: Python v_train = asc_train_rnd + b_time_rnd * TRAIN_TT_SCALED + b_cost * TRAIN_COST_SCALED v_swissmetro = asc_sm_rnd + b_time_rnd * SM_TT_SCALED + b_cost * SM_COST_SCALED v_car = asc_car_rnd + b_time_rnd * CAR_TT_SCALED + b_cost * CAR_CO_SCALED .. GENERATED FROM PYTHON SOURCE LINES 84-85 Associate utility functions with the numbering of alternatives. .. GENERATED FROM PYTHON SOURCE LINES 85-87 .. code-block:: Python v = {1: v_train, 2: v_swissmetro, 3: v_car} .. GENERATED FROM PYTHON SOURCE LINES 88-89 Associate the availability conditions with the alternatives. .. GENERATED FROM PYTHON SOURCE LINES 89-91 .. code-block:: Python av = {1: TRAIN_AV_SP, 2: SM_AV, 3: CAR_AV_SP} .. GENERATED FROM PYTHON SOURCE LINES 92-94 Conditional on the random parameters, the likelihood of one observation is given by the logit model (called the kernel). .. GENERATED FROM PYTHON SOURCE LINES 94-96 .. code-block:: Python choice_probability_one_observation = logit(v, av, CHOICE) .. GENERATED FROM PYTHON SOURCE LINES 97-100 Conditional on the random parameters, the likelihood of all observations for one individual (the trajectory) is the product of the likelihood of each observation. .. GENERATED FROM PYTHON SOURCE LINES 100-104 .. code-block:: Python conditional_trajectory_probability = PanelLikelihoodTrajectory( choice_probability_one_observation ) .. GENERATED FROM PYTHON SOURCE LINES 105-106 We integrate over the random parameters using Monte-Carlo .. GENERATED FROM PYTHON SOURCE LINES 106-108 .. code-block:: Python log_probability = log(MonteCarlo(conditional_trajectory_probability)) .. GENERATED FROM PYTHON SOURCE LINES 109-110 We retrieve the parameters estimates. .. GENERATED FROM PYTHON SOURCE LINES 110-119 .. code-block:: Python try: results = EstimationResults.from_yaml_file(filename='saved_results/b12panel.yaml') except FileNotFoundError: sys.exit( 'Run first the script b12panel.py ' 'in order to generate the ' 'file b12panel.yaml and move it to the directory saved_results.' ) .. GENERATED FROM PYTHON SOURCE LINES 120-122 Simulate to recalculate the log likelihood directly from the formula, without the Biogeme object .. GENERATED FROM PYTHON SOURCE LINES 122-132 .. code-block:: Python simulated_loglike = calculate_single_formula_from_expression( expression=log_probability, database=database, number_of_draws=NUMBER_OF_DRAWS, the_betas=results.get_beta_values(), numerically_safe=False, second_derivatives_mode=SecondDerivativesMode.NEVER, use_jit=True, ) .. rst-class:: sphx-glr-script-out .. code-block:: none Flattening database [(6768, 38)]. Database flattened [(752, 362)] Flattening database [(6768, 38)]. Database flattened [(752, 362)] Flattening database [(6768, 38)]. Database flattened [(752, 362)] Flattening database [(6768, 38)]. Database flattened [(752, 362)] .. GENERATED FROM PYTHON SOURCE LINES 133-135 .. code-block:: Python print(f'Simulated log likelihood: {simulated_loglike}') .. rst-class:: sphx-glr-script-out .. code-block:: none Simulated log likelihood: -3574.7230792117853 .. GENERATED FROM PYTHON SOURCE LINES 136-137 We also calculate the individual parameters for the time coefficient. .. GENERATED FROM PYTHON SOURCE LINES 137-145 .. code-block:: Python numerator = MonteCarlo(b_time_rnd * conditional_trajectory_probability) denominator = MonteCarlo(conditional_trajectory_probability) simulate = { 'Numerator': numerator, 'Denominator': denominator, } .. GENERATED FROM PYTHON SOURCE LINES 146-147 Creation of the Biogeme object. .. GENERATED FROM PYTHON SOURCE LINES 147-149 .. code-block:: Python biosim = BIOGEME(database, simulate, number_of_draws=NUMBER_OF_DRAWS, seed=1223) .. rst-class:: sphx-glr-script-out .. code-block:: none Biogeme parameters read from biogeme.toml. Flattening database [(6768, 38)]. Database flattened [(752, 362)] Flattening database [(6768, 38)]. Database flattened [(752, 362)] Flattening database [(6768, 38)]. Database flattened [(752, 362)] Flattening database [(6768, 38)]. Database flattened [(752, 362)] .. GENERATED FROM PYTHON SOURCE LINES 150-151 Simulation. .. GENERATED FROM PYTHON SOURCE LINES 151-153 .. code-block:: Python sim = biosim.simulate(results.get_beta_values()) .. GENERATED FROM PYTHON SOURCE LINES 154-156 .. code-block:: Python sim['Individual-level parameters'] = sim['Numerator'] / sim['Denominator'] .. GENERATED FROM PYTHON SOURCE LINES 157-159 .. code-block:: Python print(f'{sim.shape=}') .. rst-class:: sphx-glr-script-out .. code-block:: none sim.shape=(752, 3) .. GENERATED FROM PYTHON SOURCE LINES 160-161 .. code-block:: Python display(sim) .. rst-class:: sphx-glr-script-out .. code-block:: none Numerator Denominator Individual-level parameters 0 -0.071506 0.010421 -6.861496 1 -5.400202 0.762721 -7.080176 2 -2.116223 0.269996 -7.837982 3 -3.576361 0.463302 -7.719288 4 -3.351181 0.476229 -7.036911 .. ... ... ... 747 -0.008062 0.001617 -4.987226 748 -0.000451 0.000129 -3.498607 749 -0.000016 0.000007 -2.303046 750 -5.102442 0.697186 -7.318624 751 -0.034959 0.024362 -1.434957 [752 rows x 3 columns] .. rst-class:: sphx-glr-timing **Total running time of the script:** (1 minutes 54.029 seconds) .. _sphx_glr_download_auto_examples_swissmetro_plot_b13panel_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_b13panel_simul.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_b13panel_simul.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_b13panel_simul.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_