.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/swissmetro/plot_b11cnl_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_b11cnl_simul.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_swissmetro_plot_b11cnl_simul.py:


Simulation of a cross-nested logit model
========================================

Illustration of the application of an estimated model.

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

.. GENERATED FROM PYTHON SOURCE LINES 11-21

.. code-block:: Python


    import sys

    import biogeme.biogeme as bio
    import biogeme.exceptions as excep
    import biogeme.results as res
    from biogeme import models
    from biogeme.expressions import Beta, Derive
    from biogeme.nests import OneNestForCrossNestedLogit, NestsForCrossNestedLogit








.. GENERATED FROM PYTHON SOURCE LINES 22-23

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

.. GENERATED FROM PYTHON SOURCE LINES 23-36

.. code-block:: Python

    from swissmetro_data import (
        database,
        SM_AV,
        CAR_AV_SP,
        TRAIN_AV_SP,
        TRAIN_TT,
        TRAIN_COST_SCALED,
        SM_TT,
        SM_COST_SCALED,
        CAR_TT,
        CAR_CO_SCALED,
    )








.. GENERATED FROM PYTHON SOURCE LINES 37-40

Simulation must be done with estimated value of the
parameters. We set them to zero here , and read the values from the
file created after estimation.

.. GENERATED FROM PYTHON SOURCE LINES 40-46

.. 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_TIME = Beta('B_TIME', 0, None, None, 0)
    B_COST = Beta('B_COST', 0, None, None, 0)








.. GENERATED FROM PYTHON SOURCE LINES 47-48

Nest parameters.

.. GENERATED FROM PYTHON SOURCE LINES 48-51

.. code-block:: Python

    MU_EXISTING = Beta('MU_EXISTING', 1, 1, None, 0)
    MU_PUBLIC = Beta('MU_PUBLIC', 1, 1, None, 0)








.. GENERATED FROM PYTHON SOURCE LINES 52-53

Nest membership parameters.

.. GENERATED FROM PYTHON SOURCE LINES 53-56

.. code-block:: Python

    ALPHA_EXISTING = Beta('ALPHA_EXISTING', 0.5, 0, 1, 0)
    ALPHA_PUBLIC = 1 - ALPHA_EXISTING








.. GENERATED FROM PYTHON SOURCE LINES 57-60

Definition of the utility functions. Note that we need to have
explicitly the unscaled variables if we want to calculate the
elasticities.

.. GENERATED FROM PYTHON SOURCE LINES 60-64

.. code-block:: Python

    V1 = ASC_TRAIN + B_TIME * TRAIN_TT / 100 + B_COST * TRAIN_COST_SCALED
    V2 = ASC_SM + B_TIME * SM_TT / 100 + B_COST * SM_COST_SCALED
    V3 = ASC_CAR + B_TIME * CAR_TT / 100 + B_COST * CAR_CO_SCALED








.. GENERATED FROM PYTHON SOURCE LINES 65-66

Associate utility functions with the numbering of alternatives.

.. GENERATED FROM PYTHON SOURCE LINES 66-68

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 69-70

Associate the availability conditions with the alternatives.

.. GENERATED FROM PYTHON SOURCE LINES 70-72

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 73-74

Definition of nests.

.. GENERATED FROM PYTHON SOURCE LINES 74-89

.. code-block:: Python


    nest_existing = OneNestForCrossNestedLogit(
        nest_param=MU_EXISTING,
        dict_of_alpha={1: ALPHA_EXISTING, 2: 0.0, 3: 1.0},
        name='existing',
    )

    nest_public = OneNestForCrossNestedLogit(
        nest_param=MU_PUBLIC, dict_of_alpha={1: ALPHA_PUBLIC, 2: 1.0, 3: 0.0}, name='public'
    )

    nests = NestsForCrossNestedLogit(
        choice_set=[1, 2, 3], tuple_of_nests=(nest_existing, nest_public)
    )








.. GENERATED FROM PYTHON SOURCE LINES 90-91

Read the estimation results from the pickle file.

.. GENERATED FROM PYTHON SOURCE LINES 91-97

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 98-100

.. code-block:: Python

    print('Estimation results: ', results.get_estimated_parameters())





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

 .. code-block:: none

    Estimation results:                     Value  Rob. Std err  Rob. t-test  Rob. p-value
    ALPHA_EXISTING  0.495072      0.034752    14.246000  0.000000e+00
    ASC_CAR        -0.240459      0.053450    -4.498764  6.834958e-06
    ASC_TRAIN       0.098277      0.069978     1.404403  1.601989e-01
    B_COST         -0.818885      0.058971   -13.886130  0.000000e+00
    B_TIME         -0.776846      0.102380    -7.587854  3.241851e-14
    MU_EXISTING     2.514876      0.248325    10.127363  0.000000e+00
    MU_PUBLIC       4.113595      0.496722     8.281479  2.220446e-16




.. GENERATED FROM PYTHON SOURCE LINES 101-108

.. code-block:: Python

    print('Calculating correlation matrix. It may generate numerical warnings from scipy.')
    corr = nests.correlation(
        parameters=results.get_beta_values(),
        alternatives_names={1: 'Train', 2: 'Swissmetro', 3: 'Car'},
    )
    corr





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

 .. code-block:: none

    Calculating correlation matrix. It may generate numerical warnings from scipy.
    /Users/bierlair/venv312/lib/python3.12/site-packages/scipy/integrate/_quadpack_py.py:1260: IntegrationWarning: The algorithm does not converge.  Roundoff error is detected
      in the extrapolation table.  It is assumed that the requested tolerance
      cannot be achieved, and that the returned result (if full_output = 1) is 
      the best which can be obtained.
      quad_r = quad(f, low, high, args=args, full_output=self.full_output,
    /Users/bierlair/venv312/lib/python3.12/site-packages/scipy/integrate/_quadpack_py.py:1260: IntegrationWarning: The integral is probably divergent, or slowly convergent.
      quad_r = quad(f, low, high, args=args, full_output=self.full_output,
    /Users/bierlair/venv312/lib/python3.12/site-packages/scipy/integrate/_quadpack_py.py:1260: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
      If increasing the limit yields no improvement it is advised to analyze 
      the integrand in order to determine the difficulties.  If the position of a 
      local difficulty can be determined (singularity, discontinuity) one will 
      probably gain from splitting up the interval and calling the integrator 
      on the subranges.  Perhaps a special-purpose integrator should be used.
      quad_r = quad(f, low, high, args=args, full_output=self.full_output,


.. 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>Train</th>
          <th>Swissmetro</th>
          <th>Car</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>Train</th>
          <td>1.000000</td>
          <td>6.178004e-01</td>
          <td>5.527228e-01</td>
        </tr>
        <tr>
          <th>Swissmetro</th>
          <td>0.617800</td>
          <td>1.000000e+00</td>
          <td>8.281820e-12</td>
        </tr>
        <tr>
          <th>Car</th>
          <td>0.552723</td>
          <td>8.281820e-12</td>
          <td>1.000000e+00</td>
        </tr>
      </tbody>
    </table>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 109-110

The choice model is a cross-nested logit, with availability conditions.

.. GENERATED FROM PYTHON SOURCE LINES 110-114

.. code-block:: Python

    prob1 = models.cnl(V, av, nests, 1)
    prob2 = models.cnl(V, av, nests, 2)
    prob3 = models.cnl(V, av, nests, 3)








.. GENERATED FROM PYTHON SOURCE LINES 115-118

We calculate elasticities. It is important that the variables
explicitly appear as such in the model. If not, the derivative will
be zero, as well as the elasticities.

.. GENERATED FROM PYTHON SOURCE LINES 118-122

.. code-block:: Python

    genelas1 = Derive(prob1, 'TRAIN_TT') * TRAIN_TT / prob1
    genelas2 = Derive(prob2, 'SM_TT') * SM_TT / prob2
    genelas3 = Derive(prob3, 'CAR_TT') * CAR_TT / prob3








.. GENERATED FROM PYTHON SOURCE LINES 123-124

We report the probability of each alternative and the elasticities.

.. GENERATED FROM PYTHON SOURCE LINES 124-133

.. code-block:: Python

    simulate = {
        'Prob. train': prob1,
        'Prob. Swissmetro': prob2,
        'Prob. car': prob3,
        'Elas. 1': genelas1,
        'Elas. 2': genelas2,
        'Elas. 3': genelas3,
    }








.. GENERATED FROM PYTHON SOURCE LINES 134-135

Create the Biogeme object.

.. GENERATED FROM PYTHON SOURCE LINES 135-138

.. code-block:: Python

    biosim = bio.BIOGEME(database, simulate)
    biosim.modelName = 'b11cnl_simul'








.. GENERATED FROM PYTHON SOURCE LINES 139-140

Perform the simulation.

.. GENERATED FROM PYTHON SOURCE LINES 140-142

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 143-146

.. code-block:: Python

    print('Simulation results')
    simresults





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

 .. code-block:: none

    Simulation results


.. 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>Prob. train</th>
          <th>Prob. Swissmetro</th>
          <th>Prob. car</th>
          <th>Elas. 1</th>
          <th>Elas. 2</th>
          <th>Elas. 3</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>0.151845</td>
          <td>0.627166</td>
          <td>0.220989</td>
          <td>-1.712384</td>
          <td>-0.214602</td>
          <td>-1.238134</td>
        </tr>
        <tr>
          <th>1</th>
          <td>0.191734</td>
          <td>0.648405</td>
          <td>0.159861</td>
          <td>-1.371095</td>
          <td>-0.197321</td>
          <td>-1.485995</td>
        </tr>
        <tr>
          <th>2</th>
          <td>0.107228</td>
          <td>0.604797</td>
          <td>0.287975</td>
          <td>-2.231103</td>
          <td>-0.232560</td>
          <td>-0.994183</td>
        </tr>
        <tr>
          <th>3</th>
          <td>0.117302</td>
          <td>0.534994</td>
          <td>0.347704</td>
          <td>-1.902932</td>
          <td>-0.282801</td>
          <td>-0.549944</td>
        </tr>
        <tr>
          <th>4</th>
          <td>0.122367</td>
          <td>0.645963</td>
          <td>0.231670</td>
          <td>-2.053924</td>
          <td>-0.192771</td>
          <td>-0.886473</td>
        </tr>
        <tr>
          <th>...</th>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
        </tr>
        <tr>
          <th>8446</th>
          <td>0.202324</td>
          <td>0.688887</td>
          <td>0.108788</td>
          <td>-1.246225</td>
          <td>-0.140072</td>
          <td>-1.875101</td>
        </tr>
        <tr>
          <th>8447</th>
          <td>0.170793</td>
          <td>0.659829</td>
          <td>0.169378</td>
          <td>-1.482446</td>
          <td>-0.162456</td>
          <td>-0.972154</td>
        </tr>
        <tr>
          <th>8448</th>
          <td>0.138390</td>
          <td>0.653981</td>
          <td>0.207629</td>
          <td>-1.626039</td>
          <td>-0.151241</td>
          <td>-0.853487</td>
        </tr>
        <tr>
          <th>8449</th>
          <td>0.147761</td>
          <td>0.704414</td>
          <td>0.147825</td>
          <td>-1.709362</td>
          <td>-0.132484</td>
          <td>-0.990522</td>
        </tr>
        <tr>
          <th>8450</th>
          <td>0.198138</td>
          <td>0.666110</td>
          <td>0.135751</td>
          <td>-1.349649</td>
          <td>-0.163036</td>
          <td>-1.349502</td>
        </tr>
      </tbody>
    </table>
    <p>6768 rows × 6 columns</p>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 147-148

.. code-block:: Python

    print(f'Aggregate share of train: {100*simresults["Prob. train"].mean():.1f}%')




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

 .. code-block:: none

    Aggregate share of train: 13.1%





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

   **Total running time of the script:** (3 minutes 18.801 seconds)


.. _sphx_glr_download_auto_examples_swissmetro_plot_b11cnl_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_b11cnl_simul.ipynb <plot_b11cnl_simul.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

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

    .. container:: sphx-glr-download sphx-glr-download-zip

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


.. only:: html

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

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