.. 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 <sphx_glr_download_auto_examples_swissmetro_plot_b13panel_simul.py>`
        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.

:author: Michel Bierlaire, EPFL
:date: Sun Apr  9 18:16:29 2023

.. GENERATED FROM PYTHON SOURCE LINES 13-29

.. code-block:: default


    import sys
    import pickle
    import biogeme.biogeme_logging as blog
    import biogeme.biogeme as bio
    from biogeme import models
    import biogeme.exceptions as excep
    import biogeme.results as res
    from biogeme.expressions import (
        Beta,
        bioDraws,
        PanelLikelihoodTrajectory,
        MonteCarlo,
        log,
    )








.. GENERATED FROM PYTHON SOURCE LINES 30-31

See the data processing script: :ref:`swissmetro_panel`.

.. GENERATED FROM PYTHON SOURCE LINES 31-47

.. code-block:: default

    from swissmetro_panel 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,
    )

    print(f'Samples size = {database.getSampleSize()}')





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Samples size = 752




.. GENERATED FROM PYTHON SOURCE LINES 48-51

We use a low number of draws, as the objective is to illustrate the
syntax. In practice, this value is insufficient to have a good
approximation of the integral.

.. GENERATED FROM PYTHON SOURCE LINES 51-56

.. code-block:: default

    NUMBER_OF_DRAWS = 50

    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 57-58

Parameters.

.. GENERATED FROM PYTHON SOURCE LINES 58-60

.. code-block:: default

    B_COST = Beta('B_COST', 0, None, None, 0)








.. GENERATED FROM PYTHON SOURCE LINES 61-63

Define a random parameter, normally distributed across individuals,
designed to be used for Monte-Carlo simulation.

.. GENERATED FROM PYTHON SOURCE LINES 63-67

.. code-block:: default

    B_TIME = Beta('B_TIME', 0, None, None, 0)
    B_TIME_S = Beta('B_TIME_S', 1, None, None, 0)
    B_TIME_RND = B_TIME + B_TIME_S * bioDraws('B_TIME_RND', 'NORMAL_ANTI')








.. GENERATED FROM PYTHON SOURCE LINES 68-69

We do the same for the constants, to address serial correlation.

.. GENERATED FROM PYTHON SOURCE LINES 69-81

.. code-block:: default

    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 * bioDraws('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 * bioDraws('ASC_TRAIN_RND', 'NORMAL_ANTI')

    ASC_SM = Beta('ASC_SM', 0, None, None, 1)
    ASC_SM_S = Beta('ASC_SM_S', 1, None, None, 0)
    ASC_SM_RND = ASC_SM + ASC_SM_S * bioDraws('ASC_SM_RND', 'NORMAL_ANTI')








.. GENERATED FROM PYTHON SOURCE LINES 82-83

Definition of the utility functions.

.. GENERATED FROM PYTHON SOURCE LINES 83-87

.. code-block:: default

    V1 = ASC_TRAIN_RND + B_TIME_RND * TRAIN_TT_SCALED + B_COST * TRAIN_COST_SCALED
    V2 = ASC_SM_RND + B_TIME_RND * SM_TT_SCALED + B_COST * SM_COST_SCALED
    V3 = ASC_CAR_RND + B_TIME_RND * CAR_TT_SCALED + B_COST * CAR_CO_SCALED








.. GENERATED FROM PYTHON SOURCE LINES 88-89

Associate utility functions with the numbering of alternatives.

.. GENERATED FROM PYTHON SOURCE LINES 89-91

.. code-block:: default

    V = {1: V1, 2: V2, 3: V3}








.. GENERATED FROM PYTHON SOURCE LINES 92-93

Associate the availability conditions with the alternatives.

.. GENERATED FROM PYTHON SOURCE LINES 93-95

.. code-block:: default

    av = {1: TRAIN_AV_SP, 2: SM_AV, 3: CAR_AV_SP}








.. GENERATED FROM PYTHON SOURCE LINES 96-98

Conditional on the random parameters, the likelihood of one observation is
given by the logit model (called the kernel).

.. GENERATED FROM PYTHON SOURCE LINES 98-100

.. code-block:: default

    obsprob = models.logit(V, av, CHOICE)








.. GENERATED FROM PYTHON SOURCE LINES 101-104

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 104-106

.. code-block:: default

    condprobIndiv = PanelLikelihoodTrajectory(obsprob)








.. GENERATED FROM PYTHON SOURCE LINES 107-108

We integrate over the random parameters using Monte-Carlo

.. GENERATED FROM PYTHON SOURCE LINES 108-110

.. code-block:: default

    logprob = log(MonteCarlo(condprobIndiv))








.. GENERATED FROM PYTHON SOURCE LINES 111-112

We retrieve the parameters estimates.

.. GENERATED FROM PYTHON SOURCE LINES 112-121

.. code-block:: default

    try:
        results = res.bioResults(pickleFile='saved_results/b12panel.pickle')
    except excep.BiogemeError:
        sys.exit(
            'Run first the script b12panel.py '
            'in order to generate the '
            'file b12panel.pickle.'
        )








.. GENERATED FROM PYTHON SOURCE LINES 122-124

Simulate to recalculate the log likelihood directly from the
formula, without the Biogeme object

.. GENERATED FROM PYTHON SOURCE LINES 124-132

.. code-block:: default

    simulated_loglike = logprob.getValue_c(
        database=database,
        betas=results.getBetaValues(),
        numberOfDraws=NUMBER_OF_DRAWS,
        aggregation=True,
        prepareIds=True,
    )








.. GENERATED FROM PYTHON SOURCE LINES 133-135

.. code-block:: default

    print(f'Simulated log likelihood: {simulated_loglike}')





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Simulated log likelihood: -3746.058773770821




.. 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:: default

    numerator = MonteCarlo(B_TIME_RND * condprobIndiv)
    denominator = MonteCarlo(condprobIndiv)

    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:: default

    biosim = bio.BIOGEME(database, simulate, parameter_file='few_draws.toml')





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    File few_draws.toml has been parsed. 




.. GENERATED FROM PYTHON SOURCE LINES 150-151

Suimulation.

.. GENERATED FROM PYTHON SOURCE LINES 151-153

.. code-block:: default

    sim = biosim.simulate(results.getBetaValues())








.. GENERATED FROM PYTHON SOURCE LINES 154-156

.. code-block:: default

    sim['Individual-level parameters'] = sim['Numerator'] / sim['Denominator']








.. GENERATED FROM PYTHON SOURCE LINES 157-159

.. code-block:: default

    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:: default

    sim





.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>Numerator</th>
          <th>Denominator</th>
          <th>Individual-level parameters</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>1</th>
          <td>-0.061842</td>
          <td>0.010676</td>
          <td>-5.792511</td>
        </tr>
        <tr>
          <th>2</th>
          <td>-5.461292</td>
          <td>0.745673</td>
          <td>-7.323982</td>
        </tr>
        <tr>
          <th>3</th>
          <td>-2.165100</td>
          <td>0.284843</td>
          <td>-7.601041</td>
        </tr>
        <tr>
          <th>4</th>
          <td>-3.931312</td>
          <td>0.472809</td>
          <td>-8.314807</td>
        </tr>
        <tr>
          <th>5</th>
          <td>-3.746690</td>
          <td>0.444783</td>
          <td>-8.423627</td>
        </tr>
        <tr>
          <th>...</th>
          <td>...</td>
          <td>...</td>
          <td>...</td>
        </tr>
        <tr>
          <th>935</th>
          <td>-0.010185</td>
          <td>0.002201</td>
          <td>-4.628451</td>
        </tr>
        <tr>
          <th>936</th>
          <td>-0.000394</td>
          <td>0.000110</td>
          <td>-3.596566</td>
        </tr>
        <tr>
          <th>937</th>
          <td>-0.000049</td>
          <td>0.000016</td>
          <td>-3.106041</td>
        </tr>
        <tr>
          <th>938</th>
          <td>-5.311541</td>
          <td>0.684250</td>
          <td>-7.762576</td>
        </tr>
        <tr>
          <th>939</th>
          <td>-0.000230</td>
          <td>0.000540</td>
          <td>-0.425370</td>
        </tr>
      </tbody>
    </table>
    <p>752 rows × 3 columns</p>
    </div>
    </div>
    <br />
    <br />


.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 1.431 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-python

      :download:`Download Python source code: plot_b13panel_simul.py <plot_b13panel_simul.py>`

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_b13panel_simul.ipynb <plot_b13panel_simul.ipynb>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_