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

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

.. _sphx_glr_auto_examples_latent_plot_b07problem_simul.py:


Investigation of the estimation problem
=======================================

 This file is an updated version of 07problem.py, where
 the probabilities are simulated in order to
 investigate the numerical issue.

:author: Michel Bierlaire, EPFL
:date: Thu Apr 13 18:21:15 2023

.. GENERATED FROM PYTHON SOURCE LINES 14-44

.. code-block:: default


    import biogeme.biogeme as bio
    from biogeme.models import piecewiseFormula
    import biogeme.biogeme_logging as blog
    from biogeme.expressions import Beta, Elem, bioNormalCdf

    from optima import (
        database,
        age_65_more,
        ScaledIncome,
        moreThanOneCar,
        moreThanOneBike,
        individualHouse,
        male,
        haveChildren,
        haveGA,
        highEducation,
        Envir01,
        Envir02,
        Envir03,
        Mobil11,
        Mobil14,
        Mobil16,
        Mobil17,
    )


    logger = blog.get_screen_logger(level=blog.INFO)
    logger.info('Example b07problem_simul.py')





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

 .. code-block:: none

    Example b07problem_simul.py 




.. GENERATED FROM PYTHON SOURCE LINES 45-46

Parameters to be estimated

.. GENERATED FROM PYTHON SOURCE LINES 46-59

.. code-block:: default

    coef_intercept = Beta('coef_intercept', 0.0, None, None, 0)
    coef_age_65_more = Beta('coef_age_65_more', 0.0, None, None, 0)
    coef_haveGA = Beta('coef_haveGA', 0.0, None, None, 0)
    coef_moreThanOneCar = Beta('coef_moreThanOneCar', 0.0, None, None, 0)
    coef_moreThanOneBike = Beta('coef_moreThanOneBike', 0.0, None, None, 0)
    coef_individualHouse = Beta('coef_individualHouse', 0.0, None, None, 0)
    coef_male = Beta('coef_male', 0.0, None, None, 0)
    coef_haveChildren = Beta('coef_haveChildren', 0.0, None, None, 0)
    coef_highEducation = Beta('coef_highEducation', 0.0, None, None, 0)

    thresholds = [None, 4, 6, 8, 10, None]
    formula_income = piecewiseFormula(variable=ScaledIncome, thresholds=thresholds)








.. GENERATED FROM PYTHON SOURCE LINES 60-61

Latent variable: structural equation.

.. GENERATED FROM PYTHON SOURCE LINES 61-74

.. code-block:: default

    CARLOVERS = (
        coef_intercept
        + coef_age_65_more * age_65_more
        + formula_income
        + coef_moreThanOneCar * moreThanOneCar
        + coef_moreThanOneBike * moreThanOneBike
        + coef_individualHouse * individualHouse
        + coef_male * male
        + coef_haveChildren * haveChildren
        + coef_haveGA * haveGA
        + coef_highEducation * highEducation
    )








.. GENERATED FROM PYTHON SOURCE LINES 75-76

Measurement equations

.. GENERATED FROM PYTHON SOURCE LINES 78-79

Intercepts.

.. GENERATED FROM PYTHON SOURCE LINES 79-87

.. code-block:: default

    INTER_Envir01 = Beta('INTER_Envir01', 0, None, None, 1)
    INTER_Envir02 = Beta('INTER_Envir02', 0, None, None, 0)
    INTER_Envir03 = Beta('INTER_Envir03', 0, None, None, 0)
    INTER_Mobil11 = Beta('INTER_Mobil11', 0, None, None, 0)
    INTER_Mobil14 = Beta('INTER_Mobil14', 0, None, None, 0)
    INTER_Mobil16 = Beta('INTER_Mobil16', 0, None, None, 0)
    INTER_Mobil17 = Beta('INTER_Mobil17', 0, None, None, 0)








.. GENERATED FROM PYTHON SOURCE LINES 88-89

Coefficients.

.. GENERATED FROM PYTHON SOURCE LINES 89-97

.. code-block:: default

    B_Envir01_F1 = Beta('B_Envir01_F1', -1, None, None, 1)
    B_Envir02_F1 = Beta('B_Envir02_F1', -1, None, None, 0)
    B_Envir03_F1 = Beta('B_Envir03_F1', 1, None, None, 0)
    B_Mobil11_F1 = Beta('B_Mobil11_F1', 1, None, None, 0)
    B_Mobil14_F1 = Beta('B_Mobil14_F1', 1, None, None, 0)
    B_Mobil16_F1 = Beta('B_Mobil16_F1', 1, None, None, 0)
    B_Mobil17_F1 = Beta('B_Mobil17_F1', 1, None, None, 0)








.. GENERATED FROM PYTHON SOURCE LINES 98-99

Linear models.

.. GENERATED FROM PYTHON SOURCE LINES 99-107

.. code-block:: default

    MODEL_Envir01 = INTER_Envir01 + B_Envir01_F1 * CARLOVERS
    MODEL_Envir02 = INTER_Envir02 + B_Envir02_F1 * CARLOVERS
    MODEL_Envir03 = INTER_Envir03 + B_Envir03_F1 * CARLOVERS
    MODEL_Mobil11 = INTER_Mobil11 + B_Mobil11_F1 * CARLOVERS
    MODEL_Mobil14 = INTER_Mobil14 + B_Mobil14_F1 * CARLOVERS
    MODEL_Mobil16 = INTER_Mobil16 + B_Mobil16_F1 * CARLOVERS
    MODEL_Mobil17 = INTER_Mobil17 + B_Mobil17_F1 * CARLOVERS








.. GENERATED FROM PYTHON SOURCE LINES 108-109

Scale parameters.

.. GENERATED FROM PYTHON SOURCE LINES 109-117

.. code-block:: default

    SIGMA_STAR_Envir01 = Beta('SIGMA_STAR_Envir01', 1, 1.0e-5, None, 1)
    SIGMA_STAR_Envir02 = Beta('SIGMA_STAR_Envir02', 0.01, 1.0e-5, None, 0)
    SIGMA_STAR_Envir03 = Beta('SIGMA_STAR_Envir03', 1, 1.0e-5, None, 0)
    SIGMA_STAR_Mobil11 = Beta('SIGMA_STAR_Mobil11', 1, 1.0e-5, None, 0)
    SIGMA_STAR_Mobil14 = Beta('SIGMA_STAR_Mobil14', 1, 1.0e-5, None, 0)
    SIGMA_STAR_Mobil16 = Beta('SIGMA_STAR_Mobil16', 1, 1.0e-5, None, 0)
    SIGMA_STAR_Mobil17 = Beta('SIGMA_STAR_Mobil17', 1, 1.0e-5, None, 0)








.. GENERATED FROM PYTHON SOURCE LINES 118-119

Symmetric thresholds.

.. GENERATED FROM PYTHON SOURCE LINES 119-126

.. code-block:: default

    delta_1 = Beta('delta_1', 0.1, 1.0e-5, None, 0)
    delta_2 = Beta('delta_2', 0.2, 1.0e-5, None, 0)
    tau_1 = -delta_1 - delta_2
    tau_2 = -delta_1
    tau_3 = delta_1
    tau_4 = delta_1 + delta_2








.. GENERATED FROM PYTHON SOURCE LINES 127-128

Ordered probit models.

.. GENERATED FROM PYTHON SOURCE LINES 128-247

.. code-block:: default

    Envir01_tau_1 = (tau_1 - MODEL_Envir01) / SIGMA_STAR_Envir01
    Envir01_tau_2 = (tau_2 - MODEL_Envir01) / SIGMA_STAR_Envir01
    Envir01_tau_3 = (tau_3 - MODEL_Envir01) / SIGMA_STAR_Envir01
    Envir01_tau_4 = (tau_4 - MODEL_Envir01) / SIGMA_STAR_Envir01
    IndEnvir01 = {
        1: bioNormalCdf(Envir01_tau_1),
        2: bioNormalCdf(Envir01_tau_2) - bioNormalCdf(Envir01_tau_1),
        3: bioNormalCdf(Envir01_tau_3) - bioNormalCdf(Envir01_tau_2),
        4: bioNormalCdf(Envir01_tau_4) - bioNormalCdf(Envir01_tau_3),
        5: 1 - bioNormalCdf(Envir01_tau_4),
        6: 1.0,
        -1: 1.0,
        -2: 1.0,
    }

    P_Envir01 = Elem(IndEnvir01, Envir01)

    Envir02_tau_1 = (tau_1 - MODEL_Envir02) / SIGMA_STAR_Envir02
    Envir02_tau_2 = (tau_2 - MODEL_Envir02) / SIGMA_STAR_Envir02
    Envir02_tau_3 = (tau_3 - MODEL_Envir02) / SIGMA_STAR_Envir02
    Envir02_tau_4 = (tau_4 - MODEL_Envir02) / SIGMA_STAR_Envir02
    IndEnvir02 = {
        1: bioNormalCdf(Envir02_tau_1),
        2: bioNormalCdf(Envir02_tau_2) - bioNormalCdf(Envir02_tau_1),
        3: bioNormalCdf(Envir02_tau_3) - bioNormalCdf(Envir02_tau_2),
        4: bioNormalCdf(Envir02_tau_4) - bioNormalCdf(Envir02_tau_3),
        5: 1 - bioNormalCdf(Envir02_tau_4),
        6: 1.0,
        -1: 1.0,
        -2: 1.0,
    }

    P_Envir02 = Elem(IndEnvir02, Envir02)

    Envir03_tau_1 = (tau_1 - MODEL_Envir03) / SIGMA_STAR_Envir03
    Envir03_tau_2 = (tau_2 - MODEL_Envir03) / SIGMA_STAR_Envir03
    Envir03_tau_3 = (tau_3 - MODEL_Envir03) / SIGMA_STAR_Envir03
    Envir03_tau_4 = (tau_4 - MODEL_Envir03) / SIGMA_STAR_Envir03
    IndEnvir03 = {
        1: bioNormalCdf(Envir03_tau_1),
        2: bioNormalCdf(Envir03_tau_2) - bioNormalCdf(Envir03_tau_1),
        3: bioNormalCdf(Envir03_tau_3) - bioNormalCdf(Envir03_tau_2),
        4: bioNormalCdf(Envir03_tau_4) - bioNormalCdf(Envir03_tau_3),
        5: 1 - bioNormalCdf(Envir03_tau_4),
        6: 1.0,
        -1: 1.0,
        -2: 1.0,
    }

    P_Envir03 = Elem(IndEnvir03, Envir03)

    Mobil11_tau_1 = (tau_1 - MODEL_Mobil11) / SIGMA_STAR_Mobil11
    Mobil11_tau_2 = (tau_2 - MODEL_Mobil11) / SIGMA_STAR_Mobil11
    Mobil11_tau_3 = (tau_3 - MODEL_Mobil11) / SIGMA_STAR_Mobil11
    Mobil11_tau_4 = (tau_4 - MODEL_Mobil11) / SIGMA_STAR_Mobil11
    IndMobil11 = {
        1: bioNormalCdf(Mobil11_tau_1),
        2: bioNormalCdf(Mobil11_tau_2) - bioNormalCdf(Mobil11_tau_1),
        3: bioNormalCdf(Mobil11_tau_3) - bioNormalCdf(Mobil11_tau_2),
        4: bioNormalCdf(Mobil11_tau_4) - bioNormalCdf(Mobil11_tau_3),
        5: 1 - bioNormalCdf(Mobil11_tau_4),
        6: 1.0,
        -1: 1.0,
        -2: 1.0,
    }

    P_Mobil11 = Elem(IndMobil11, Mobil11)

    Mobil14_tau_1 = (tau_1 - MODEL_Mobil14) / SIGMA_STAR_Mobil14
    Mobil14_tau_2 = (tau_2 - MODEL_Mobil14) / SIGMA_STAR_Mobil14
    Mobil14_tau_3 = (tau_3 - MODEL_Mobil14) / SIGMA_STAR_Mobil14
    Mobil14_tau_4 = (tau_4 - MODEL_Mobil14) / SIGMA_STAR_Mobil14
    IndMobil14 = {
        1: bioNormalCdf(Mobil14_tau_1),
        2: bioNormalCdf(Mobil14_tau_2) - bioNormalCdf(Mobil14_tau_1),
        3: bioNormalCdf(Mobil14_tau_3) - bioNormalCdf(Mobil14_tau_2),
        4: bioNormalCdf(Mobil14_tau_4) - bioNormalCdf(Mobil14_tau_3),
        5: 1 - bioNormalCdf(Mobil14_tau_4),
        6: 1.0,
        -1: 1.0,
        -2: 1.0,
    }

    P_Mobil14 = Elem(IndMobil14, Mobil14)

    Mobil16_tau_1 = (tau_1 - MODEL_Mobil16) / SIGMA_STAR_Mobil16
    Mobil16_tau_2 = (tau_2 - MODEL_Mobil16) / SIGMA_STAR_Mobil16
    Mobil16_tau_3 = (tau_3 - MODEL_Mobil16) / SIGMA_STAR_Mobil16
    Mobil16_tau_4 = (tau_4 - MODEL_Mobil16) / SIGMA_STAR_Mobil16
    IndMobil16 = {
        1: bioNormalCdf(Mobil16_tau_1),
        2: bioNormalCdf(Mobil16_tau_2) - bioNormalCdf(Mobil16_tau_1),
        3: bioNormalCdf(Mobil16_tau_3) - bioNormalCdf(Mobil16_tau_2),
        4: bioNormalCdf(Mobil16_tau_4) - bioNormalCdf(Mobil16_tau_3),
        5: 1 - bioNormalCdf(Mobil16_tau_4),
        6: 1.0,
        -1: 1.0,
        -2: 1.0,
    }

    P_Mobil16 = Elem(IndMobil16, Mobil16)

    Mobil17_tau_1 = (tau_1 - MODEL_Mobil17) / SIGMA_STAR_Mobil17
    Mobil17_tau_2 = (tau_2 - MODEL_Mobil17) / SIGMA_STAR_Mobil17
    Mobil17_tau_3 = (tau_3 - MODEL_Mobil17) / SIGMA_STAR_Mobil17
    Mobil17_tau_4 = (tau_4 - MODEL_Mobil17) / SIGMA_STAR_Mobil17
    IndMobil17 = {
        1: bioNormalCdf(Mobil17_tau_1),
        2: bioNormalCdf(Mobil17_tau_2) - bioNormalCdf(Mobil17_tau_1),
        3: bioNormalCdf(Mobil17_tau_3) - bioNormalCdf(Mobil17_tau_2),
        4: bioNormalCdf(Mobil17_tau_4) - bioNormalCdf(Mobil17_tau_3),
        5: 1 - bioNormalCdf(Mobil17_tau_4),
        6: 1.0,
        -1: 1.0,
        -2: 1.0,
    }

    P_Mobil17 = Elem(IndMobil17, Mobil17)








.. GENERATED FROM PYTHON SOURCE LINES 248-249

Expressions ot simulate.

.. GENERATED FROM PYTHON SOURCE LINES 249-258

.. code-block:: default

    simulate = {
        'P_Envir01': P_Envir01,
        'P_Envir02': P_Envir02,
        'P_Envir03': P_Envir03,
        'P_Mobil11': P_Mobil11,
        'P_Mobil14': P_Mobil14,
        'P_Mobil16': P_Mobil16,
        'P_Mobil17': P_Mobil17,
    }







.. GENERATED FROM PYTHON SOURCE LINES 259-260

Extract the list of parameters in the simulated expressions.

.. GENERATED FROM PYTHON SOURCE LINES 260-264

.. code-block:: default

    beta_values = {}
    for expr in simulate.values():
        beta_values = beta_values | expr.get_beta_values()








.. GENERATED FROM PYTHON SOURCE LINES 265-266

Create the Biogeme object

.. GENERATED FROM PYTHON SOURCE LINES 266-269

.. code-block:: default

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





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

 .. code-block:: none

    File biogeme.toml has been parsed. 




.. GENERATED FROM PYTHON SOURCE LINES 271-274

.. code-block:: default

    simulated_values = biosim.simulate(theBetaValues=beta_values)
    simulated_values






.. 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>P_Envir01</th>
          <th>P_Envir02</th>
          <th>P_Envir03</th>
          <th>P_Mobil11</th>
          <th>P_Mobil14</th>
          <th>P_Mobil16</th>
          <th>P_Mobil17</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>0.079656</td>
          <td>0.000000e+00</td>
          <td>0.382089</td>
          <td>0.382089</td>
          <td>0.382089</td>
          <td>0.382089</td>
          <td>0.382089</td>
        </tr>
        <tr>
          <th>2</th>
          <td>1.000000</td>
          <td>1.000000e+00</td>
          <td>1.000000</td>
          <td>1.000000</td>
          <td>1.000000</td>
          <td>1.000000</td>
          <td>1.000000</td>
        </tr>
        <tr>
          <th>3</th>
          <td>0.079656</td>
          <td>7.619853e-24</td>
          <td>0.079656</td>
          <td>0.078084</td>
          <td>0.382089</td>
          <td>0.079656</td>
          <td>0.078084</td>
        </tr>
        <tr>
          <th>4</th>
          <td>0.382089</td>
          <td>1.000000e+00</td>
          <td>0.079656</td>
          <td>0.382089</td>
          <td>0.382089</td>
          <td>0.382089</td>
          <td>0.078084</td>
        </tr>
        <tr>
          <th>5</th>
          <td>0.078084</td>
          <td>0.000000e+00</td>
          <td>0.079656</td>
          <td>0.382089</td>
          <td>0.078084</td>
          <td>0.382089</td>
          <td>0.079656</td>
        </tr>
        <tr>
          <th>...</th>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
        </tr>
        <tr>
          <th>2259</th>
          <td>0.079656</td>
          <td>0.000000e+00</td>
          <td>0.078084</td>
          <td>0.078084</td>
          <td>0.382089</td>
          <td>0.079656</td>
          <td>0.079656</td>
        </tr>
        <tr>
          <th>2261</th>
          <td>0.078084</td>
          <td>1.000000e+00</td>
          <td>0.078084</td>
          <td>0.382089</td>
          <td>1.000000</td>
          <td>0.078084</td>
          <td>0.078084</td>
        </tr>
        <tr>
          <th>2262</th>
          <td>0.078084</td>
          <td>1.000000e+00</td>
          <td>0.078084</td>
          <td>0.382089</td>
          <td>1.000000</td>
          <td>0.078084</td>
          <td>0.078084</td>
        </tr>
        <tr>
          <th>2263</th>
          <td>0.382089</td>
          <td>4.906714e-198</td>
          <td>0.382089</td>
          <td>0.382089</td>
          <td>0.382089</td>
          <td>0.382089</td>
          <td>1.000000</td>
        </tr>
        <tr>
          <th>2264</th>
          <td>0.382089</td>
          <td>4.906714e-198</td>
          <td>0.382089</td>
          <td>0.382089</td>
          <td>0.382089</td>
          <td>0.382089</td>
          <td>1.000000</td>
        </tr>
      </tbody>
    </table>
    <p>1906 rows × 7 columns</p>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 275-277

We identify the entries for which the likelihood is zero, so that
the log likelihood cannot be computed

.. GENERATED FROM PYTHON SOURCE LINES 277-280

.. code-block:: default

    zero_values = simulated_values.where(simulated_values == 0, other='')
    zero_values






.. 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>P_Envir01</th>
          <th>P_Envir02</th>
          <th>P_Envir03</th>
          <th>P_Mobil11</th>
          <th>P_Mobil14</th>
          <th>P_Mobil16</th>
          <th>P_Mobil17</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td></td>
          <td>0.0</td>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
        </tr>
        <tr>
          <th>2</th>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
        </tr>
        <tr>
          <th>3</th>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
        </tr>
        <tr>
          <th>4</th>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
        </tr>
        <tr>
          <th>5</th>
          <td></td>
          <td>0.0</td>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
        </tr>
        <tr>
          <th>...</th>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
        </tr>
        <tr>
          <th>2259</th>
          <td></td>
          <td>0.0</td>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
        </tr>
        <tr>
          <th>2261</th>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
        </tr>
        <tr>
          <th>2262</th>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
        </tr>
        <tr>
          <th>2263</th>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
        </tr>
        <tr>
          <th>2264</th>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
          <td></td>
        </tr>
      </tbody>
    </table>
    <p>1906 rows × 7 columns</p>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 281-282

Diagnostic.

.. GENERATED FROM PYTHON SOURCE LINES 282-283

.. code-block:: default

    print('The problematic model is Envir02')




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

 .. code-block:: none

    The problematic model is Envir02





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

   **Total running time of the script:** (0 minutes 0.100 seconds)


.. _sphx_glr_download_auto_examples_latent_plot_b07problem_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_b07problem_simul.py <plot_b07problem_simul.py>`

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

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


.. only:: html

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

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