"""

Ordinal probit model
====================

Example of an ordinal probit model.  This is just to illustrate the
syntax, as the data are not ordered.  But the example assume, for the
sake of it, that the alternatives are ordered as 1->2->3

Michel Bierlaire, EPFL
Thu Jun 26 2025, 15:54:37
"""

from IPython.core.display_functions import display

# %%
# See the data processing script: :ref:`swissmetro_data`.
from swissmetro_data import CHOICE, TRAIN_COST_SCALED, TRAIN_TT_SCALED, database

import biogeme.biogeme_logging as blog
from biogeme.biogeme import BIOGEME
from biogeme.expressions import Beta, Elem, log
from biogeme.models import ordered_probit
from biogeme.results_processing import get_pandas_estimated_parameters

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

# %%
# Parameters to be estimated
b_time = Beta('b_time', 0, None, None, 0)
b_cost = Beta('b_cost', 0, None, None, 0)

# %%
# Threshold parameters for the ordered probit.

# %%
# :math:`\tau_1 \leq 0`.
tau1 = Beta('tau1', -1, None, 0, 0)

# %%
# :math:`\delta_2 \geq 0`.
delta2 = Beta('delta2', 2, 0, None, 0)

# %%
# :math:`\tau_2 = \tau_1 + \delta_2`
tau2 = tau1 + delta2

# %%
# Utility
utility = b_time * TRAIN_TT_SCALED + b_cost * TRAIN_COST_SCALED

# %%
# Associate each discrete indicator with an interval.
#
#   1. :math:`-\infty \to \tau_1`,
#   2. :math:`\tau_1 \to \tau_2`,
#   3. :math:`\tau_2 \to +\infty`.
the_probability = ordered_probit(
    continuous_value=utility,
    list_of_discrete_values=[1, 2, 3],
    reference_threshold_parameter=tau1,
    scale_parameter=1.0,
)

# %%
# Extract from the dict the formula associated with the observed choice.
the_chosen_proba = Elem(the_probability, CHOICE)

# %%
# Definition of the model. This is the contribution of each
# observation to the log likelihood function.
log_probability = log(the_chosen_proba)

# %%
# Create the Biogeme object.
the_biogeme = BIOGEME(database, log_probability)
the_biogeme.model_name = 'b18ordinal_probit'

# %%
# Estimate the parameters.
results = the_biogeme.estimate()

# %%
print(results.short_summary())

# %%
pandas_results = get_pandas_estimated_parameters(estimation_results=results)
display(pandas_results)