.. 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:: default 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:: default 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:: default 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:: default 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:: default 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:: default 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:: default 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:: default 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-88 .. code-block:: default try: results = res.bioResults(pickleFile='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 89-90 Conditional to B_TIME_RND, we have a logit model (called the kernel) .. GENERATED FROM PYTHON SOURCE LINES 90-92 .. code-block:: default prob = models.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:: default integral = MonteCarlo(prob) integralSquare = MonteCarlo(prob * prob) variance = integralSquare - 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:: default numerator = MonteCarlo(B_TIME_RND * prob) denominator = integral .. GENERATED FROM PYTHON SOURCE LINES 106-113 .. code-block:: default 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:: default biosim = bio.BIOGEME(database, simulate) biosim.modelName = 'b05normal_mixture_simul' .. GENERATED FROM PYTHON SOURCE LINES 119-120 Simulate the requested quantities. The output is a Pandas data frame. .. GENERATED FROM PYTHON SOURCE LINES 120-122 .. code-block:: default simresults = biosim.simulate(results.getBetaValues()) .. GENERATED FROM PYTHON SOURCE LINES 123-124 95% confidence interval on the log likelihood. .. GENERATED FROM PYTHON SOURCE LINES 124-131 .. code-block:: default 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:396: RuntimeWarning: invalid value encountered in log result = getattr(ufunc, method)(*inputs, **kwargs) .. GENERATED FROM PYTHON SOURCE LINES 132-134 .. code-block:: default print(f'Log likelihood: {np.log(simresults["Integral"]).sum()}') .. rst-class:: sphx-glr-script-out .. code-block:: none Log likelihood: -5216.316731053104 .. GENERATED FROM PYTHON SOURCE LINES 135-140 .. code-block:: default 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 500 draws: 700.5475488837479 .. GENERATED FROM PYTHON SOURCE LINES 141-143 .. code-block:: default print(f'In average {simresults["Integration error"].mean()} per observation.') .. rst-class:: sphx-glr-script-out .. code-block:: none In average 0.10350879859393439 per observation. .. GENERATED FROM PYTHON SOURCE LINES 144-145 95% confidence interval .. GENERATED FROM PYTHON SOURCE LINES 145-150 .. code-block:: default print( f'95% confidence interval: [{simresults["left"].sum()} - ' f'{simresults["right"].sum()}]' ) .. rst-class:: sphx-glr-script-out .. code-block:: none 95% confidence interval: [-7529.657511926459 - -2766.3838983681285] .. GENERATED FROM PYTHON SOURCE LINES 151-152 Post processing to obtain the individual parameters. .. GENERATED FROM PYTHON SOURCE LINES 152-154 .. code-block:: default simresults['beta'] = simresults['Numerator'] / simresults['Denominator'] .. GENERATED FROM PYTHON SOURCE LINES 155-156 Plot the histogram of individual parameters .. GENERATED FROM PYTHON SOURCE LINES 156-160 .. code-block:: default 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 161-162 Plot the general distribution of beta .. GENERATED FROM PYTHON SOURCE LINES 162-176 .. code-block:: default def normalpdf(val, mu=0.0, std=1.0): """ 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 177-180 .. code-block:: default betas = results.getBetaValues(['B_TIME', 'B_TIME_S']) x = np.arange(simresults['beta'].min(), simresults['beta'].max(), 0.01) .. GENERATED FROM PYTHON SOURCE LINES 181-184 .. code-block:: default 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 8.663 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-python :download:`Download Python source code: plot_b05normal_mixture_simul.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_b05normal_mixture_simul.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_