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

.. code-block:: default


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








.. GENERATED FROM PYTHON SOURCE LINES 21-22

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

.. GENERATED FROM PYTHON SOURCE LINES 22-38

.. code-block:: default

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








.. GENERATED FROM PYTHON SOURCE LINES 39-42

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 42-48

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








.. GENERATED FROM PYTHON SOURCE LINES 49-50

Nest parameters.

.. GENERATED FROM PYTHON SOURCE LINES 50-53

.. code-block:: default

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








.. GENERATED FROM PYTHON SOURCE LINES 54-55

Nest membership parameters.

.. GENERATED FROM PYTHON SOURCE LINES 55-58

.. code-block:: default

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








.. GENERATED FROM PYTHON SOURCE LINES 59-62

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 62-66

.. code-block:: default

    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 67-68

Associate utility functions with the numbering of alternatives.

.. GENERATED FROM PYTHON SOURCE LINES 68-70

.. code-block:: default

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








.. GENERATED FROM PYTHON SOURCE LINES 71-72

Associate the availability conditions with the alternatives.

.. GENERATED FROM PYTHON SOURCE LINES 72-74

.. code-block:: default

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








.. GENERATED FROM PYTHON SOURCE LINES 75-76

Definition of nests.

.. GENERATED FROM PYTHON SOURCE LINES 76-91

.. code-block:: default


    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 92-93

Read the estimation results from the pickle file.

.. GENERATED FROM PYTHON SOURCE LINES 93-99

.. code-block:: default

    try:
        results = res.bioResults(pickleFile='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 100-102

.. code-block:: default

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





.. 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 103-110

.. code-block:: default

    print('Calculating correlation matrix. It may generate numerical warnings from scipy.')
    corr = nests.correlation(
        parameters=results.getBetaValues(),
        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.


.. 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.000000e+00</td>
          <td>8.240277e-12</td>
          <td>8.245306e-12</td>
        </tr>
        <tr>
          <th>Swissmetro</th>
          <td>8.240277e-12</td>
          <td>1.000000e+00</td>
          <td>8.237105e-12</td>
        </tr>
        <tr>
          <th>Car</th>
          <td>8.245306e-12</td>
          <td>8.237105e-12</td>
          <td>1.000000e+00</td>
        </tr>
      </tbody>
    </table>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 111-112

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

.. GENERATED FROM PYTHON SOURCE LINES 112-116

.. code-block:: default

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





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

 .. code-block:: none

    /Users/bierlair/OnlineFiles/FilesOnGoogleDrive/github/biogeme/docs/examples/swissmetro/plot_b11cnl_simul.py:112: DeprecationWarning: The function cnl_avail is deprecated. It has been replaced by the function cnl
      prob1 = models.cnl_avail(V, av, nests, 1)
    /Users/bierlair/OnlineFiles/FilesOnGoogleDrive/github/biogeme/docs/examples/swissmetro/plot_b11cnl_simul.py:113: DeprecationWarning: The function cnl_avail is deprecated. It has been replaced by the function cnl
      prob2 = models.cnl_avail(V, av, nests, 2)
    /Users/bierlair/OnlineFiles/FilesOnGoogleDrive/github/biogeme/docs/examples/swissmetro/plot_b11cnl_simul.py:114: DeprecationWarning: The function cnl_avail is deprecated. It has been replaced by the function cnl
      prob3 = models.cnl_avail(V, av, nests, 3)




.. GENERATED FROM PYTHON SOURCE LINES 117-120

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 120-124

.. code-block:: default

    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 125-126

We report the probability of each alternative and the elasticities.

.. GENERATED FROM PYTHON SOURCE LINES 126-135

.. code-block:: default

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








.. GENERATED FROM PYTHON SOURCE LINES 136-137

Create the Biogeme object.

.. GENERATED FROM PYTHON SOURCE LINES 137-140

.. code-block:: default

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








.. GENERATED FROM PYTHON SOURCE LINES 141-142

Perform the simulation.

.. GENERATED FROM PYTHON SOURCE LINES 142-144

.. code-block:: default

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








.. GENERATED FROM PYTHON SOURCE LINES 145-148

.. code-block:: default

    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 149-150

.. code-block:: default

    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:** (1 minutes 25.901 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-python

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

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

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


.. only:: html

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

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