.. 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-30

.. code-block:: Python


    import sys

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








.. GENERATED FROM PYTHON SOURCE LINES 31-32

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

.. GENERATED FROM PYTHON SOURCE LINES 32-48

.. code-block:: Python

    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.get_sample_size()}')





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

 .. code-block:: none

    Samples size = 752




.. GENERATED FROM PYTHON SOURCE LINES 49-52

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

.. code-block:: Python

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

Parameters.

.. GENERATED FROM PYTHON SOURCE LINES 59-61

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 62-64

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

.. GENERATED FROM PYTHON SOURCE LINES 64-68

.. code-block:: Python

    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 69-70

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

.. GENERATED FROM PYTHON SOURCE LINES 70-82

.. 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 * 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 83-84

Definition of the utility functions.

.. GENERATED FROM PYTHON SOURCE LINES 84-88

.. code-block:: Python

    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 89-90

Associate utility functions with the numbering of alternatives.

.. GENERATED FROM PYTHON SOURCE LINES 90-92

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 93-94

Associate the availability conditions with the alternatives.

.. GENERATED FROM PYTHON SOURCE LINES 94-96

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 97-99

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

.. GENERATED FROM PYTHON SOURCE LINES 99-101

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 102-105

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 105-107

.. code-block:: Python

    condprobIndiv = PanelLikelihoodTrajectory(obsprob)








.. GENERATED FROM PYTHON SOURCE LINES 108-109

We integrate over the random parameters using Monte-Carlo

.. GENERATED FROM PYTHON SOURCE LINES 109-111

.. code-block:: Python

    logprob = log(MonteCarlo(condprobIndiv))








.. GENERATED FROM PYTHON SOURCE LINES 112-113

We retrieve the parameters estimates.

.. GENERATED FROM PYTHON SOURCE LINES 113-122

.. code-block:: Python

    try:
        results = res.bioResults(pickle_file='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 123-125

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

.. GENERATED FROM PYTHON SOURCE LINES 125-133

.. code-block:: Python

    simulated_loglike = logprob.get_value_c(
        database=database,
        betas=results.get_beta_values(),
        number_of_draws=NUMBER_OF_DRAWS,
        aggregation=True,
        prepare_ids=True,
    )








.. GENERATED FROM PYTHON SOURCE LINES 134-136

.. code-block:: Python

    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 137-138

We also calculate the individual parameters for the time coefficient.

.. GENERATED FROM PYTHON SOURCE LINES 138-146

.. code-block:: Python

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

    simulate = {
        'Numerator': numerator,
        'Denominator': denominator,
    }








.. GENERATED FROM PYTHON SOURCE LINES 147-148

Creation of the Biogeme object.

.. GENERATED FROM PYTHON SOURCE LINES 148-150

.. code-block:: Python

    biosim = bio.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. 




.. GENERATED FROM PYTHON SOURCE LINES 151-152

Suimulation.

.. GENERATED FROM PYTHON SOURCE LINES 152-154

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 155-157

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 158-160

.. 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 161-162

.. code-block:: Python

    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.062321</td>
          <td>0.009942</td>
          <td>-6.268342</td>
        </tr>
        <tr>
          <th>2</th>
          <td>-5.534657</td>
          <td>0.742892</td>
          <td>-7.450155</td>
        </tr>
        <tr>
          <th>3</th>
          <td>-2.157851</td>
          <td>0.241836</td>
          <td>-8.922785</td>
        </tr>
        <tr>
          <th>4</th>
          <td>-3.585943</td>
          <td>0.446293</td>
          <td>-8.034958</td>
        </tr>
        <tr>
          <th>5</th>
          <td>-3.211187</td>
          <td>0.455042</td>
          <td>-7.056910</td>
        </tr>
        <tr>
          <th>...</th>
          <td>...</td>
          <td>...</td>
          <td>...</td>
        </tr>
        <tr>
          <th>935</th>
          <td>-0.006234</td>
          <td>0.001198</td>
          <td>-5.201502</td>
        </tr>
        <tr>
          <th>936</th>
          <td>-0.000400</td>
          <td>0.000205</td>
          <td>-1.949459</td>
        </tr>
        <tr>
          <th>937</th>
          <td>-0.000012</td>
          <td>0.000006</td>
          <td>-1.942529</td>
        </tr>
        <tr>
          <th>938</th>
          <td>-5.267642</td>
          <td>0.630936</td>
          <td>-8.348929</td>
        </tr>
        <tr>
          <th>939</th>
          <td>0.017225</td>
          <td>0.045669</td>
          <td>0.377173</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 0.934 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 <plot_b13panel_simul.ipynb>`

    .. 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-zip

      :download:`Download zipped: plot_b13panel_simul.zip <plot_b13panel_simul.zip>`


.. only:: html

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

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