.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/swissmetro/plot_b05normal_mixture_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_b05normal_mixture_simul.py: Simulation of a mixture model ============================= Simulation of the mixture model, with estimation of the integration error. Michel Bierlaire, EPFL Fri Jun 20 2025, 10:29:35 .. GENERATED FROM PYTHON SOURCE LINES 11-20 .. code-block:: Python import sys import numpy as np from biogeme.biogeme import BIOGEME from biogeme.expressions import Beta, Draws, MonteCarlo from biogeme.models import logit from biogeme.results_processing import EstimationResults .. GENERATED FROM PYTHON SOURCE LINES 21-22 See the data processing script: :ref:`swissmetro_data`. .. GENERATED FROM PYTHON SOURCE LINES 22-45 .. code-block:: Python from swissmetro_data 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, ) try: import matplotlib.pyplot as plt PLOT = True except ModuleNotFoundError: print('Install matplotlib to see the distribution of integration errors.') print('pip install matplotlib') PLOT = False .. GENERATED FROM PYTHON SOURCE LINES 46-47 Parameters to be estimated. .. GENERATED FROM PYTHON SOURCE LINES 47-52 .. 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_cost = Beta('b_cost', 0, None, None, 0) .. GENERATED FROM PYTHON SOURCE LINES 53-55 Define a random parameter, normally distributed, designed to be used for Monte-Carlo simulation. .. GENERATED FROM PYTHON SOURCE LINES 55-57 .. code-block:: Python b_time = Beta('b_time', 0, None, None, 0) .. GENERATED FROM PYTHON SOURCE LINES 58-59 It is advised not to use 0 as starting value for the following parameter. .. GENERATED FROM PYTHON SOURCE LINES 59-62 .. 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') .. GENERATED FROM PYTHON SOURCE LINES 63-64 Definition of the utility functions. .. GENERATED FROM PYTHON SOURCE LINES 64-68 .. code-block:: Python v_train = asc_train + b_time_rnd * TRAIN_TT_SCALED + b_cost * TRAIN_COST_SCALED v_swissmetro = asc_sm + b_time_rnd * SM_TT_SCALED + b_cost * SM_COST_SCALED v_car = asc_car + b_time_rnd * CAR_TT_SCALED + b_cost * CAR_CO_SCALED .. GENERATED FROM PYTHON SOURCE LINES 69-70 Associate utility functions with the numbering of alternatives. .. GENERATED FROM PYTHON SOURCE LINES 70-72 .. code-block:: Python v = {1: v_train, 2: v_swissmetro, 3: v_car} .. GENERATED FROM PYTHON SOURCE LINES 73-74 Associate the availability conditions with the alternatives. .. GENERATED FROM PYTHON SOURCE LINES 74-76 .. code-block:: Python av = {1: TRAIN_AV_SP, 2: SM_AV, 3: CAR_AV_SP} .. GENERATED FROM PYTHON SOURCE LINES 77-78 The estimation results are read from the pickle file. .. GENERATED FROM PYTHON SOURCE LINES 78-88 .. code-block:: Python try: results = EstimationResults.from_yaml_file( filename='saved_results/b05normal_mixture.yaml' ) except FileNotFoundError: print( 'Run first the script plot_b05normal_mixture.py in order to generate the ' 'file b05normal_mixture.yaml.' ) sys.exit() .. GENERATED FROM PYTHON SOURCE LINES 89-90 Conditional to b_time_rnd, we have a logit model (called the kernel) .. GENERATED FROM PYTHON SOURCE LINES 90-92 .. code-block:: Python conditional_probability = logit(v, av, CHOICE) .. GENERATED FROM PYTHON SOURCE LINES 93-95 We calculate the integration error. Note that this formula assumes independent draws, and is not valid for Halton or antithetic draws. .. GENERATED FROM PYTHON SOURCE LINES 95-100 .. code-block:: Python integral = MonteCarlo(conditional_probability) integral_square = MonteCarlo(conditional_probability * conditional_probability) variance = integral_square - integral * integral error = (variance / 2.0) ** 0.5 .. GENERATED FROM PYTHON SOURCE LINES 101-102 And the value of the individual parameters. .. GENERATED FROM PYTHON SOURCE LINES 102-105 .. code-block:: Python numerator = MonteCarlo(b_time_rnd * conditional_probability) denominator = integral .. GENERATED FROM PYTHON SOURCE LINES 106-113 .. code-block:: Python simulate = { 'Numerator': numerator, 'Denominator': denominator, 'Integral': integral, 'Integration error': error, } .. GENERATED FROM PYTHON SOURCE LINES 114-115 Create the Biogeme object. .. GENERATED FROM PYTHON SOURCE LINES 115-118 .. code-block:: Python biosim = BIOGEME(database, simulate, number_of_draws=10000) biosim.model_name = 'b05normal_mixture_simul' .. GENERATED FROM PYTHON SOURCE LINES 119-120 NUmber of draws .. GENERATED FROM PYTHON SOURCE LINES 120-122 .. code-block:: Python print(f'Number of draws: {biosim.number_of_draws}') .. rst-class:: sphx-glr-script-out .. code-block:: none Number of draws: 10000 .. GENERATED FROM PYTHON SOURCE LINES 123-124 Simulate the requested quantities. The output is a Pandas data frame. .. GENERATED FROM PYTHON SOURCE LINES 124-126 .. code-block:: Python simulation_results = biosim.simulate(results.get_beta_values()) .. GENERATED FROM PYTHON SOURCE LINES 127-128 95% confidence interval on the log likelihood. .. GENERATED FROM PYTHON SOURCE LINES 128-135 .. code-block:: Python simulation_results['left'] = np.log( simulation_results['Integral'] - 1.96 * simulation_results['Integration error'] ) simulation_results['right'] = np.log( simulation_results['Integral'] + 1.96 * simulation_results['Integration error'] ) .. rst-class:: sphx-glr-script-out .. code-block:: none /Users/bierlair/python_envs/venv313/lib/python3.13/site-packages/pandas/core/arraylike.py:399: RuntimeWarning: invalid value encountered in log result = getattr(ufunc, method)(*inputs, **kwargs) .. GENERATED FROM PYTHON SOURCE LINES 136-138 .. code-block:: Python print(f'Log likelihood: {np.log(simulation_results["Integral"]).sum()}') .. rst-class:: sphx-glr-script-out .. code-block:: none Log likelihood: -5214.689681629532 .. GENERATED FROM PYTHON SOURCE LINES 139-144 .. code-block:: Python print( f'Integration error for {biosim.number_of_draws} draws: ' f'{simulation_results["Integration error"].sum()}' ) .. rst-class:: sphx-glr-script-out .. code-block:: none Integration error for 10000 draws: 716.8094872624056 .. GENERATED FROM PYTHON SOURCE LINES 145-147 .. code-block:: Python print(f'In average {simulation_results["Integration error"].mean()} per observation.') .. rst-class:: sphx-glr-script-out .. code-block:: none In average 0.1059115672669039 per observation. .. GENERATED FROM PYTHON SOURCE LINES 148-149 95% confidence interval .. GENERATED FROM PYTHON SOURCE LINES 149-154 .. code-block:: Python print( f'95% confidence interval: [{simulation_results["left"].sum()} - ' f'{simulation_results["right"].sum()}]' ) .. rst-class:: sphx-glr-script-out .. code-block:: none 95% confidence interval: [-7583.2407468481015 - -2719.046943930449] .. GENERATED FROM PYTHON SOURCE LINES 155-156 Post processing to obtain the individual parameters. .. GENERATED FROM PYTHON SOURCE LINES 156-160 .. code-block:: Python simulation_results['Beta'] = ( simulation_results['Numerator'] / simulation_results['Denominator'] ) .. GENERATED FROM PYTHON SOURCE LINES 161-162 Plot the histogram of individual parameters .. GENERATED FROM PYTHON SOURCE LINES 162-166 .. code-block:: Python if PLOT: simulation_results['Beta'].plot(kind='hist', density=True, bins=20) .. image-sg:: /auto_examples/swissmetro/images/sphx_glr_plot_b05normal_mixture_simul_001.png :alt: plot b05normal mixture simul :srcset: /auto_examples/swissmetro/images/sphx_glr_plot_b05normal_mixture_simul_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 167-168 Plot the general distribution of Beta .. GENERATED FROM PYTHON SOURCE LINES 168-182 .. code-block:: Python def normalpdf(val: float, mu: float = 0.0, std: float = 1.0) -> float: """ Calculate the pdf of the normal distribution, for plotting purposes. """ d = -(val - mu) * (val - mu) n = 2.0 * std * std a = d / n num = np.exp(a) den = std * 2.506628275 p = num / den return p .. GENERATED FROM PYTHON SOURCE LINES 183-186 .. code-block:: Python betas = results.get_beta_values(['b_time', 'b_time_s']) x = np.arange(simulation_results['Beta'].min(), simulation_results['Beta'].max(), 0.01) .. GENERATED FROM PYTHON SOURCE LINES 187-190 .. code-block:: Python if PLOT: plt.plot(x, normalpdf(x, betas['b_time'], betas['b_time_s']), '-') plt.show() .. image-sg:: /auto_examples/swissmetro/images/sphx_glr_plot_b05normal_mixture_simul_002.png :alt: plot b05normal mixture simul :srcset: /auto_examples/swissmetro/images/sphx_glr_plot_b05normal_mixture_simul_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 24.630 seconds) .. _sphx_glr_download_auto_examples_swissmetro_plot_b05normal_mixture_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_b05normal_mixture_simul.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_b05normal_mixture_simul.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_b05normal_mixture_simul.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_