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

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

.. _sphx_glr_auto_examples_latentbis_plot_m03_simultaneous_estimation.py:


Choice model with the latent variable: maximum likelihood estimation
====================================================================

Mixture of logit.
Measurement equation for the indicators.
Maximum likelihood (full information) estimation.

:author: Michel Bierlaire, EPFL
:date: Fri Apr 14 10:07:43 2023

.. GENERATED FROM PYTHON SOURCE LINES 13-58

.. code-block:: default


    import sys
    from functools import reduce
    import biogeme.biogeme_logging as blog
    import biogeme.biogeme as bio
    import biogeme.exceptions as excep
    import biogeme.distributions as dist
    import biogeme.results as res
    from biogeme import models
    from biogeme.expressions import (
        Beta,
        Variable,
        log,
        RandomVariable,
        Integrate,
        Elem,
        bioNormalCdf,
        exp,
    )

    from read_or_estimate import read_or_estimate

    from optima import (
        database,
        male,
        age,
        haveChildren,
        highEducation,
        childCenter,
        childSuburb,
        SocioProfCat,
        WaitingTimePT,
        Choice,
        TimePT_scaled,
        TimeCar_scaled,
        MarginalCostPT_scaled,
        CostCarCHF_scaled,
        distance_km_scaled,
        PurpHWH,
        PurpOther,
    )

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





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

 .. code-block:: none

    Example m03_simultaneous_estimation.py 




.. GENERATED FROM PYTHON SOURCE LINES 59-60

Read the estimates from the structural equation estimation.

.. GENERATED FROM PYTHON SOURCE LINES 60-72

.. code-block:: default

    MODELNAME = 'm01_latent_variable'
    try:
        struct_results = res.bioResults(pickleFile=f'saved_results/{MODELNAME}.pickle')
    except excep.BiogemeError:
        print(
            f'Run first the script {MODELNAME}.py in order to generate the '
            f'file {MODELNAME}.pickle, and move it to the directory saved_results'
        )
        sys.exit()
    struct_betas = struct_results.getBetaValues()









.. GENERATED FROM PYTHON SOURCE LINES 73-74

Coefficients

.. GENERATED FROM PYTHON SOURCE LINES 74-94

.. code-block:: default

    coef_intercept = Beta('coef_intercept', struct_betas['coef_intercept'], None, None, 0)
    coef_age_30_less = Beta(
        'coef_age_30_less', struct_betas['coef_age_30_less'], None, None, 0
    )
    coef_haveChildren = Beta(
        'coef_haveChildren', struct_betas['coef_haveChildren'], None, None, 0
    )
    coef_highEducation = Beta(
        'coef_highEducation', struct_betas['coef_highEducation'], None, None, 0
    )
    coef_artisans = Beta('coef_artisans', struct_betas['coef_artisans'], None, None, 0)
    coef_employees = Beta('coef_employees', struct_betas['coef_employees'], None, None, 0)
    coef_male = Beta('coef_male', struct_betas['coef_male'], None, None, 0)
    coef_child_center = Beta(
        'coef_child_center', struct_betas['coef_child_center'], None, None, 0
    )
    coef_child_suburb = Beta(
        'coef_child_suburb', struct_betas['coef_child_suburb'], None, None, 0
    )








.. GENERATED FROM PYTHON SOURCE LINES 95-96

Latent variable: structural equation

.. GENERATED FROM PYTHON SOURCE LINES 98-100

Define a random parameter, normally distributed, designed to be used
for numerical integration

.. GENERATED FROM PYTHON SOURCE LINES 100-104

.. code-block:: default

    omega = RandomVariable('omega')
    density = dist.normalpdf(omega)
    sigma_s = Beta('sigma_s', 1, None, None, 0)








.. GENERATED FROM PYTHON SOURCE LINES 105-118

.. code-block:: default

    ACTIVELIFE = (
        coef_intercept
        + coef_child_center * childCenter
        + coef_child_suburb * childSuburb
        + coef_highEducation * highEducation
        + coef_artisans * (SocioProfCat == 5)
        + coef_employees * (SocioProfCat == 6)
        + coef_age_30_less * (age <= 30)
        + coef_male * male
        + coef_haveChildren * haveChildren
        + sigma_s * omega
    )








.. GENERATED FROM PYTHON SOURCE LINES 119-120

Measurement equations

.. GENERATED FROM PYTHON SOURCE LINES 120-129

.. code-block:: default

    indicators = [
        'ResidCh01',
        'ResidCh04',
        'ResidCh05',
        'ResidCh06',
        'LifSty07',
        'LifSty10',
    ]








.. GENERATED FROM PYTHON SOURCE LINES 130-131

We define the intercept parameters. The first one is normalized to 0.

.. GENERATED FROM PYTHON SOURCE LINES 131-134

.. code-block:: default

    inter = {k: Beta(f'inter_{k}', 0, None, None, 0) for k in indicators[1:]}
    inter[indicators[0]] = Beta(f'INTER_{indicators[0]}', 0, None, None, 1)








.. GENERATED FROM PYTHON SOURCE LINES 135-136

We define the coefficients. The first one is normalized to 1.

.. GENERATED FROM PYTHON SOURCE LINES 136-139

.. code-block:: default

    coefficients = {k: Beta(f'coeff_{k}', 0, None, None, 0) for k in indicators[1:]}
    coefficients[indicators[0]] = Beta(f'B_{indicators[0]}', 1, None, None, 1)








.. GENERATED FROM PYTHON SOURCE LINES 140-141

We define the measurement equations for each indicator

.. GENERATED FROM PYTHON SOURCE LINES 141-143

.. code-block:: default

    linear_models = {k: inter[k] + coefficients[k] * ACTIVELIFE for k in indicators}








.. GENERATED FROM PYTHON SOURCE LINES 144-145

We define the scale parameters of the error terms.

.. GENERATED FROM PYTHON SOURCE LINES 145-148

.. code-block:: default

    sigma_star = {k: Beta(f'sigma_star_{k}', 1, 1.0e-5, None, 0) for k in indicators[1:]}
    sigma_star[indicators[0]] = Beta(f'sigma_star_{indicators[0]}', 1, None, None, 1)








.. GENERATED FROM PYTHON SOURCE LINES 149-150

Symmetric threshold.

.. GENERATED FROM PYTHON SOURCE LINES 150-157

.. 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 158-159

Ordered probit models.

.. GENERATED FROM PYTHON SOURCE LINES 159-177

.. code-block:: default

    tau_1_residual = {k: (tau_1 - linear_models[k]) / sigma_star[k] for k in indicators}
    tau_2_residual = {k: (tau_2 - linear_models[k]) / sigma_star[k] for k in indicators}
    tau_3_residual = {k: (tau_3 - linear_models[k]) / sigma_star[k] for k in indicators}
    tau_4_residual = {k: (tau_4 - linear_models[k]) / sigma_star[k] for k in indicators}
    dict_prob_indicators = {
        k: {
            1: bioNormalCdf(tau_1_residual[k]),
            2: bioNormalCdf(tau_2_residual[k]) - bioNormalCdf(tau_1_residual[k]),
            3: bioNormalCdf(tau_3_residual[k]) - bioNormalCdf(tau_2_residual[k]),
            4: bioNormalCdf(tau_4_residual[k]) - bioNormalCdf(tau_3_residual[k]),
            5: 1 - bioNormalCdf(tau_4_residual[k]),
            6: 1.0,
            -1: 1.0,
            -2: 1.0,
        }
        for k in indicators
    }








.. GENERATED FROM PYTHON SOURCE LINES 178-179

Product of the likelihood of the indicators.

.. GENERATED FROM PYTHON SOURCE LINES 179-186

.. code-block:: default

    prob_indicators = reduce(
        lambda x, y: x * Elem(dict_prob_indicators[y], Variable(y)),
        indicators,
        Elem(dict_prob_indicators[indicators[0]], Variable(indicators[0])),
    )









.. GENERATED FROM PYTHON SOURCE LINES 187-190

Choice model
Read the estimates from the sequential estimation, and use
them as starting values

.. GENERATED FROM PYTHON SOURCE LINES 190-201

.. code-block:: default

    MODELNAME = 'm02_sequential_estimation'
    try:
        choice_results = res.bioResults(pickleFile=f'saved_results/{MODELNAME}.pickle')
    except excep.BiogemeError:
        print(
            f'Run first the script {MODELNAME}.py in order to generate the '
            f'file {MODELNAME}.pickle, and move it to the directory saved_results'
        )
        sys.exit()
    choice_betas = choice_results.getBetaValues()








.. GENERATED FROM PYTHON SOURCE LINES 202-228

.. code-block:: default

    ASC_CAR = Beta('ASC_CAR', choice_betas['ASC_CAR'], None, None, 0)
    ASC_PT = Beta('ASC_PT', 0, None, None, 1)
    ASC_SM = Beta('ASC_SM', choice_betas['ASC_SM'], None, None, 0)
    BETA_COST_HWH = Beta('BETA_COST_HWH', choice_betas['BETA_COST_HWH'], None, None, 0)
    BETA_COST_OTHER = Beta(
        'BETA_COST_OTHER', choice_betas['BETA_COST_OTHER'], None, None, 0
    )
    BETA_DIST = Beta('BETA_DIST', choice_betas['BETA_DIST'], None, None, 0)
    BETA_TIME_CAR_REF = Beta(
        'BETA_TIME_CAR_REF', choice_betas['BETA_TIME_CAR_REF'], None, 0, 0
    )
    BETA_TIME_CAR_AL = Beta(
        'BETA_TIME_CAR_AL', choice_betas['BETA_TIME_CAR_AL'], None, None, 0
    )
    BETA_TIME_PT_REF = Beta(
        'BETA_TIME_PT_REF', choice_betas['BETA_TIME_PT_REF'], None, 0, 0
    )
    BETA_TIME_CAR = BETA_TIME_CAR_REF * exp(BETA_TIME_CAR_AL * ACTIVELIFE)
    BETA_TIME_PT_AL = Beta(
        'BETA_TIME_PT_AL', choice_betas['BETA_TIME_PT_AL'], None, None, 0
    )
    BETA_TIME_PT = BETA_TIME_PT_REF * exp(BETA_TIME_PT_AL * ACTIVELIFE)
    BETA_WAITING_TIME = Beta(
        'BETA_WAITING_TIME', choice_betas['BETA_WAITING_TIME'], None, None, 0
    )








.. GENERATED FROM PYTHON SOURCE LINES 229-230

Definition of utility functions:

.. GENERATED FROM PYTHON SOURCE LINES 230-247

.. code-block:: default

    V0 = (
        ASC_PT
        + BETA_TIME_PT * TimePT_scaled
        + BETA_WAITING_TIME * WaitingTimePT
        + BETA_COST_HWH * MarginalCostPT_scaled * PurpHWH
        + BETA_COST_OTHER * MarginalCostPT_scaled * PurpOther
    )

    V1 = (
        ASC_CAR
        + BETA_TIME_CAR * TimeCar_scaled
        + BETA_COST_HWH * CostCarCHF_scaled * PurpHWH
        + BETA_COST_OTHER * CostCarCHF_scaled * PurpOther
    )

    V2 = ASC_SM + BETA_DIST * distance_km_scaled








.. GENERATED FROM PYTHON SOURCE LINES 248-249

Associate utility functions with the numbering of alternatives

.. GENERATED FROM PYTHON SOURCE LINES 249-251

.. code-block:: default

    V = {0: V0, 1: V1, 2: V2}








.. GENERATED FROM PYTHON SOURCE LINES 252-254

Conditional on omega, we have a logit model (called the kernel) for
the choice

.. GENERATED FROM PYTHON SOURCE LINES 254-256

.. code-block:: default

    condprob = models.logit(V, None, Choice)








.. GENERATED FROM PYTHON SOURCE LINES 257-259

Conditional on omega, we have the product of ordered probit for the
indicators.

.. GENERATED FROM PYTHON SOURCE LINES 259-261

.. code-block:: default

    condlike = prob_indicators * condprob








.. GENERATED FROM PYTHON SOURCE LINES 262-263

We integrate over omega using numerical integration

.. GENERATED FROM PYTHON SOURCE LINES 263-265

.. code-block:: default

    loglike = log(Integrate(condlike * density, 'omega'))








.. GENERATED FROM PYTHON SOURCE LINES 266-267

Create the Biogeme object

.. GENERATED FROM PYTHON SOURCE LINES 267-270

.. code-block:: default

    the_biogeme = bio.BIOGEME(database, loglike)
    the_biogeme.modelName = 'm03_simultaneous_estimation'





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

 .. code-block:: none

    File biogeme.toml has been parsed. 




.. GENERATED FROM PYTHON SOURCE LINES 271-273

If estimation results are saved on file, we read them to speed up the process.
If not, we estimate the parameters.

.. GENERATED FROM PYTHON SOURCE LINES 273-275

.. code-block:: default

    results = read_or_estimate(the_biogeme=the_biogeme, directory='saved_results')








.. GENERATED FROM PYTHON SOURCE LINES 276-278

.. code-block:: default

    print(results.short_summary())





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

 .. code-block:: none

    Results for model m03_simultaneous_estimation
    Nbr of parameters:              37
    Sample size:                    1906
    Excluded data:                  359
    Final log likelihood:           -15178.83
    Akaike Information Criterion:   30431.66
    Bayesian Information Criterion: 30637.11





.. GENERATED FROM PYTHON SOURCE LINES 279-282

.. code-block:: default

    print(f'Final log likelihood: {results.data.logLike:.3f}')
    print(f'Output file: {results.data.htmlFileName}')





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

 .. code-block:: none

    Final log likelihood: -15178.829
    Output file: m03_simultaneous_estimation.html




.. GENERATED FROM PYTHON SOURCE LINES 283-284

.. code-block:: default

    results.getEstimatedParameters()





.. 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>Value</th>
          <th>Rob. Std err</th>
          <th>Rob. t-test</th>
          <th>Rob. p-value</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>ASC_CAR</th>
          <td>0.664627</td>
          <td>0.119378</td>
          <td>5.567438</td>
          <td>2.585113e-08</td>
        </tr>
        <tr>
          <th>ASC_SM</th>
          <td>0.169755</td>
          <td>0.306409</td>
          <td>0.554016</td>
          <td>5.795680e-01</td>
        </tr>
        <tr>
          <th>BETA_COST_HWH</th>
          <td>-1.180142</td>
          <td>0.267148</td>
          <td>-4.417565</td>
          <td>9.981925e-06</td>
        </tr>
        <tr>
          <th>BETA_COST_OTHER</th>
          <td>-0.556328</td>
          <td>0.114508</td>
          <td>-4.858404</td>
          <td>1.183359e-06</td>
        </tr>
        <tr>
          <th>BETA_DIST</th>
          <td>-1.151594</td>
          <td>0.262733</td>
          <td>-4.383143</td>
          <td>1.169792e-05</td>
        </tr>
        <tr>
          <th>BETA_TIME_CAR_AL</th>
          <td>0.001630</td>
          <td>0.004710</td>
          <td>0.346112</td>
          <td>7.292588e-01</td>
        </tr>
        <tr>
          <th>BETA_TIME_CAR_REF</th>
          <td>-6.414696</td>
          <td>1.292676</td>
          <td>-4.962337</td>
          <td>6.965005e-07</td>
        </tr>
        <tr>
          <th>BETA_TIME_PT_AL</th>
          <td>-0.000477</td>
          <td>0.005487</td>
          <td>-0.086948</td>
          <td>9.307131e-01</td>
        </tr>
        <tr>
          <th>BETA_TIME_PT_REF</th>
          <td>-1.777138</td>
          <td>0.607850</td>
          <td>-2.923644</td>
          <td>3.459605e-03</td>
        </tr>
        <tr>
          <th>BETA_WAITING_TIME</th>
          <td>-0.020978</td>
          <td>0.007567</td>
          <td>-2.772418</td>
          <td>5.564161e-03</td>
        </tr>
        <tr>
          <th>coef_age_30_less</th>
          <td>47.924453</td>
          <td>3.586199</td>
          <td>13.363577</td>
          <td>0.000000e+00</td>
        </tr>
        <tr>
          <th>coef_artisans</th>
          <td>-0.347921</td>
          <td>0.381494</td>
          <td>-0.911997</td>
          <td>3.617704e-01</td>
        </tr>
        <tr>
          <th>coef_child_center</th>
          <td>2.401157</td>
          <td>2.303511</td>
          <td>1.042390</td>
          <td>2.972308e-01</td>
        </tr>
        <tr>
          <th>coef_child_suburb</th>
          <td>1.586959</td>
          <td>2.300117</td>
          <td>0.689947</td>
          <td>4.902274e-01</td>
        </tr>
        <tr>
          <th>coef_employees</th>
          <td>-0.007955</td>
          <td>0.250222</td>
          <td>-0.031793</td>
          <td>9.746373e-01</td>
        </tr>
        <tr>
          <th>coef_haveChildren</th>
          <td>13.727760</td>
          <td>2.233246</td>
          <td>6.146998</td>
          <td>7.896321e-10</td>
        </tr>
        <tr>
          <th>coef_highEducation</th>
          <td>-15.357200</td>
          <td>1.412346</td>
          <td>-10.873536</td>
          <td>0.000000e+00</td>
        </tr>
        <tr>
          <th>coef_intercept</th>
          <td>-46.875602</td>
          <td>3.602223</td>
          <td>-13.012964</td>
          <td>0.000000e+00</td>
        </tr>
        <tr>
          <th>coef_male</th>
          <td>15.747680</td>
          <td>1.296598</td>
          <td>12.145379</td>
          <td>0.000000e+00</td>
        </tr>
        <tr>
          <th>coeff_LifSty07</th>
          <td>0.219754</td>
          <td>0.032807</td>
          <td>6.698373</td>
          <td>2.107514e-11</td>
        </tr>
        <tr>
          <th>coeff_LifSty10</th>
          <td>0.220491</td>
          <td>0.027794</td>
          <td>7.932942</td>
          <td>2.220446e-15</td>
        </tr>
        <tr>
          <th>coeff_ResidCh04</th>
          <td>0.251910</td>
          <td>0.027454</td>
          <td>9.175749</td>
          <td>0.000000e+00</td>
        </tr>
        <tr>
          <th>coeff_ResidCh05</th>
          <td>0.602666</td>
          <td>0.047296</td>
          <td>12.742425</td>
          <td>0.000000e+00</td>
        </tr>
        <tr>
          <th>coeff_ResidCh06</th>
          <td>0.383646</td>
          <td>0.041095</td>
          <td>9.335647</td>
          <td>0.000000e+00</td>
        </tr>
        <tr>
          <th>delta_1</th>
          <td>33.783060</td>
          <td>3.120931</td>
          <td>10.824674</td>
          <td>0.000000e+00</td>
        </tr>
        <tr>
          <th>delta_2</th>
          <td>71.483870</td>
          <td>6.216996</td>
          <td>11.498137</td>
          <td>0.000000e+00</td>
        </tr>
        <tr>
          <th>inter_LifSty07</th>
          <td>-58.539864</td>
          <td>5.548443</td>
          <td>-10.550683</td>
          <td>0.000000e+00</td>
        </tr>
        <tr>
          <th>inter_LifSty10</th>
          <td>-13.721599</td>
          <td>2.182355</td>
          <td>-6.287520</td>
          <td>3.225777e-10</td>
        </tr>
        <tr>
          <th>inter_ResidCh04</th>
          <td>13.152788</td>
          <td>2.031354</td>
          <td>6.474887</td>
          <td>9.488255e-11</td>
        </tr>
        <tr>
          <th>inter_ResidCh05</th>
          <td>-111.806138</td>
          <td>9.999910</td>
          <td>-11.180714</td>
          <td>0.000000e+00</td>
        </tr>
        <tr>
          <th>inter_ResidCh06</th>
          <td>-29.530546</td>
          <td>3.713242</td>
          <td>-7.952767</td>
          <td>1.776357e-15</td>
        </tr>
        <tr>
          <th>sigma_s</th>
          <td>69.519678</td>
          <td>6.177397</td>
          <td>11.253878</td>
          <td>0.000000e+00</td>
        </tr>
        <tr>
          <th>sigma_star_LifSty07</th>
          <td>78.148260</td>
          <td>7.047972</td>
          <td>11.088049</td>
          <td>0.000000e+00</td>
        </tr>
        <tr>
          <th>sigma_star_LifSty10</th>
          <td>66.904224</td>
          <td>5.969069</td>
          <td>11.208486</td>
          <td>0.000000e+00</td>
        </tr>
        <tr>
          <th>sigma_star_ResidCh04</th>
          <td>63.302189</td>
          <td>5.767499</td>
          <td>10.975674</td>
          <td>0.000000e+00</td>
        </tr>
        <tr>
          <th>sigma_star_ResidCh05</th>
          <td>91.806284</td>
          <td>8.382677</td>
          <td>10.951905</td>
          <td>0.000000e+00</td>
        </tr>
        <tr>
          <th>sigma_star_ResidCh06</th>
          <td>93.782658</td>
          <td>8.467819</td>
          <td>11.075185</td>
          <td>0.000000e+00</td>
        </tr>
      </tbody>
    </table>
    </div>
    </div>
    <br />
    <br />


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

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


.. _sphx_glr_download_auto_examples_latentbis_plot_m03_simultaneous_estimation.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_m03_simultaneous_estimation.py <plot_m03_simultaneous_estimation.py>`

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

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


.. only:: html

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

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