.. 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. :author: Michel Bierlaire, EPFL :date: Sun Apr 9 17:47:42 2023 .. GENERATED FROM PYTHON SOURCE LINES 12-21 .. code-block:: Python import sys import numpy as np import biogeme.biogeme as bio from biogeme import models import biogeme.results as res from biogeme.exceptions import BiogemeError from biogeme.expressions import Beta, bioDraws, MonteCarlo .. GENERATED FROM PYTHON SOURCE LINES 22-23 See the data processing script: :ref:`swissmetro_data`. .. GENERATED FROM PYTHON SOURCE LINES 23-46 .. code-block:: Python from swissmetro_data import ( database, CHOICE, SM_AV, CAR_AV_SP, TRAIN_AV_SP, TRAIN_TT_SCALED, TRAIN_COST_SCALED, SM_TT_SCALED, SM_COST_SCALED, CAR_TT_SCALED, CAR_CO_SCALED, ) 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 47-48 Parameters to be estimated. .. GENERATED FROM PYTHON SOURCE LINES 48-53 .. 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 54-56 Define a random parameter, normally distributed, 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 * bioDraws('b_time_rnd', 'NORMAL') .. GENERATED FROM PYTHON SOURCE LINES 64-65 Definition of the utility functions. .. GENERATED FROM PYTHON SOURCE LINES 65-69 .. code-block:: Python V1 = ASC_TRAIN + B_TIME_RND * TRAIN_TT_SCALED + B_COST * TRAIN_COST_SCALED V2 = ASC_SM + B_TIME_RND * SM_TT_SCALED + B_COST * SM_COST_SCALED V3 = ASC_CAR + B_TIME_RND * CAR_TT_SCALED + B_COST * CAR_CO_SCALED .. GENERATED FROM PYTHON SOURCE LINES 70-71 Associate utility functions with the numbering of alternatives. .. GENERATED FROM PYTHON SOURCE LINES 71-73 .. code-block:: Python V = {1: V1, 2: V2, 3: V3} .. GENERATED FROM PYTHON SOURCE LINES 74-75 Associate the availability conditions with the alternatives. .. GENERATED FROM PYTHON SOURCE LINES 75-77 .. code-block:: Python av = {1: TRAIN_AV_SP, 2: SM_AV, 3: CAR_AV_SP} .. GENERATED FROM PYTHON SOURCE LINES 78-79 The estimation results are read from the pickle file. .. GENERATED FROM PYTHON SOURCE LINES 79-87 .. code-block:: Python try: results = res.bioResults(pickle_file='saved_results/b05_estimation_results.pickle') except BiogemeError: print( 'Run first the script 05normalMixture.py in order to generate the ' 'file 05normalMixture.pickle.' ) sys.exit() .. GENERATED FROM PYTHON SOURCE LINES 88-89 Conditional to b_time_rnd, we have a logit model (called the kernel) .. GENERATED FROM PYTHON SOURCE LINES 89-91 .. code-block:: Python prob = models.logit(V, av, CHOICE) .. GENERATED FROM PYTHON SOURCE LINES 92-94 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 94-99 .. code-block:: Python integral = MonteCarlo(prob) integralSquare = MonteCarlo(prob * prob) variance = integralSquare - integral * integral error = (variance / 2.0) ** 0.5 .. GENERATED FROM PYTHON SOURCE LINES 100-101 And the value of the individual parameters. .. GENERATED FROM PYTHON SOURCE LINES 101-104 .. code-block:: Python numerator = MonteCarlo(B_TIME_RND * prob) denominator = integral .. GENERATED FROM PYTHON SOURCE LINES 105-112 .. code-block:: Python simulate = { 'Numerator': numerator, 'Denominator': denominator, 'Integral': integral, 'Integration error': error, } .. GENERATED FROM PYTHON SOURCE LINES 113-114 Create the Biogeme object. .. GENERATED FROM PYTHON SOURCE LINES 114-117 .. code-block:: Python biosim = bio.BIOGEME(database, simulate, number_or_draws=100) biosim.modelName = 'b05normal_mixture_simul' .. GENERATED FROM PYTHON SOURCE LINES 118-119 NUmber of draws .. GENERATED FROM PYTHON SOURCE LINES 119-121 .. code-block:: Python print(biosim.number_of_draws) .. rst-class:: sphx-glr-script-out .. code-block:: none 100 .. GENERATED FROM PYTHON SOURCE LINES 122-123 Simulate the requested quantities. The output is a Pandas data frame. .. GENERATED FROM PYTHON SOURCE LINES 123-125 .. code-block:: Python simresults = biosim.simulate(results.get_beta_values()) .. GENERATED FROM PYTHON SOURCE LINES 126-127 95% confidence interval on the log likelihood. .. GENERATED FROM PYTHON SOURCE LINES 127-134 .. code-block:: Python simresults['left'] = np.log( simresults['Integral'] - 1.96 * simresults['Integration error'] ) simresults['right'] = np.log( simresults['Integral'] + 1.96 * simresults['Integration error'] ) .. rst-class:: sphx-glr-script-out .. code-block:: none /Users/bierlair/venv312/lib/python3.12/site-packages/pandas/core/arraylike.py:399: RuntimeWarning: invalid value encountered in log result = getattr(ufunc, method)(*inputs, **kwargs) .. GENERATED FROM PYTHON SOURCE LINES 135-137 .. code-block:: Python print(f'Log likelihood: {np.log(simresults["Integral"]).sum()}') .. rst-class:: sphx-glr-script-out .. code-block:: none Log likelihood: -5221.0410529590245 .. GENERATED FROM PYTHON SOURCE LINES 138-143 .. code-block:: Python print( f'Integration error for {biosim.number_of_draws} draws: ' f'{simresults["Integration error"].sum()}' ) .. rst-class:: sphx-glr-script-out .. code-block:: none Integration error for 100 draws: 695.33010691944 .. GENERATED FROM PYTHON SOURCE LINES 144-146 .. code-block:: Python print(f'In average {simresults["Integration error"].mean()} per observation.') .. rst-class:: sphx-glr-script-out .. code-block:: none In average 0.1027378999585461 per observation. .. GENERATED FROM PYTHON SOURCE LINES 147-148 95% confidence interval .. GENERATED FROM PYTHON SOURCE LINES 148-153 .. code-block:: Python print( f'95% confidence interval: [{simresults["left"].sum()} - ' f'{simresults["right"].sum()}]' ) .. rst-class:: sphx-glr-script-out .. code-block:: none 95% confidence interval: [-7619.516345268195 - -2784.741503973364] .. GENERATED FROM PYTHON SOURCE LINES 154-155 Post processing to obtain the individual parameters. .. GENERATED FROM PYTHON SOURCE LINES 155-157 .. code-block:: Python simresults['Beta'] = simresults['Numerator'] / simresults['Denominator'] .. GENERATED FROM PYTHON SOURCE LINES 158-159 Plot the histogram of individual parameters .. GENERATED FROM PYTHON SOURCE LINES 159-163 .. code-block:: Python if PLOT: simresults['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 164-165 Plot the general distribution of Beta .. GENERATED FROM PYTHON SOURCE LINES 165-179 .. 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 180-183 .. code-block:: Python betas = results.get_beta_values(['B_TIME', 'B_TIME_S']) x = np.arange(simresults['Beta'].min(), simresults['Beta'].max(), 0.01) .. GENERATED FROM PYTHON SOURCE LINES 184-187 .. 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 2.267 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 `_