23b. Binary probit model

Bayesian estimation of a binary probit model. Two alternatives: Train and Car. All observations such that the Swissmetro was chosen haven been removed from the sample.

Michel Bierlaire, EPFL Sat Jun 28 2025, 12:43:40

from pathlib import Path

from IPython.core.display_functions import display

See the data processing script: Data preparation for Swissmetro (binary choice).

from swissmetro_binary import (
    CAR_CO_SCALED,
    CAR_TT_SCALED,
    CHOICE,
    TRAIN_COST_SCALED,
    TRAIN_TT_SCALED,
    database,
)

from biogeme.bayesian_estimation import (
    BayesianResults,
    BayesianResultsSummary,
    get_pandas_estimated_parameters,
)
from biogeme.biogeme import BIOGEME
from biogeme.expressions import Beta, Elem, NormalCdf, log

Parameters to be estimated.

asc_car = Beta('asc_car', 0, None, None, 0)
b_time_car = Beta('b_time_car', 0, None, None, 0)
b_time_train = Beta('b_time_train', 0, None, None, 0)
b_cost_car = Beta('b_cost_car', 0, None, None, 0)
b_cost_train = Beta('b_cost_train', 0, None, None, 0)

Definition of the utility functions. We estimate a binary probit model. There are only two alternatives.

v_train = b_time_train * TRAIN_TT_SCALED + b_cost_train * TRAIN_COST_SCALED
v_car = asc_car + b_time_car * CAR_TT_SCALED + b_cost_car * CAR_CO_SCALED

Associate choice probability with the numbering of alternatives.

log_probability_dict = {
    1: log(NormalCdf(v_train - v_car)),
    3: log(NormalCdf(v_car - v_train)),
}

Definition of the model. This is the contribution of each observation to the log likelihood function.

log_probability = Elem(log_probability_dict, CHOICE)

Create the Biogeme object

the_biogeme = BIOGEME(database, log_probability)
the_biogeme.model_name = 'b23b_binary_probit'

Estimate the posterior distribution of the parameters, or read the results if already available.

yaml_file = Path('saved_results') / f'{the_biogeme.model_name}.yaml'
try:
    summary_results = BayesianResultsSummary.from_yaml_file(filename=yaml_file)
except FileNotFoundError:
    results: BayesianResults = the_biogeme.bayesian_estimation()
    summary_results = results.to_summary()

print(summary_results.short_summary())

Sample size                                              2232
Sampler                                                  NUTS
Number of chains                                         4
Number of draws per chain                                2000
Total number of draws                                    8000
Acceptance rate target                                   0.9
Run time                                                 0:00:23.751901
Posterior predictive log-likelihood (sum of log mean p)  -903.90
Expected log-likelihood E[log L(Y|θ)]                    -909.44
Best-draw log-likelihood (posterior upper bound)         -906.96
LOO (Leave-One-Out Cross-Validation)                     -916.41
LOO Standard Error                                       35.25
Effective number of parameters (p_LOO)                   12.52

Present the parameter estimates in a pandas table.

pandas_results = get_pandas_estimated_parameters(
    estimation_results=summary_results,
)
display(pandas_results)

           Name  Value (mean)  ...   ESS (bulk)   ESS (tail)
0  b_time_train     -0.651620  ...  4967.410324  4821.833483
1  b_cost_train     -0.981907  ...  6312.568778  4962.868173
2       asc_car     -0.353930  ...  5964.584157  4990.678539
3    b_time_car     -0.186066  ...  5762.288745  5084.061069
4    b_cost_car     -0.530560  ...  5789.071745  5178.484471

[5 rows x 12 columns]

Report the variables stored in the Bayesian estimation results.

display(summary_results.report_stored_variables())

             group           variable                dims            shape
0    constant_data      CAR_CO_SCALED               [obs]           [2232]
1    constant_data      CAR_TT_SCALED               [obs]           [2232]
2    constant_data             CHOICE               [obs]           [2232]
3    constant_data  TRAIN_COST_SCALED               [obs]           [2232]
4    constant_data    TRAIN_TT_SCALED               [obs]           [2232]
5   log_likelihood            _choice  [chain, draw, obs]  [4, 2000, 2232]
6        posterior            asc_car       [chain, draw]        [4, 2000]
7        posterior         b_cost_car       [chain, draw]        [4, 2000]
8        posterior       b_cost_train       [chain, draw]        [4, 2000]
9        posterior         b_time_car       [chain, draw]        [4, 2000]
10       posterior       b_time_train       [chain, draw]        [4, 2000]
11       posterior           log_like  [chain, draw, obs]  [4, 2000, 2232]
12           prior            asc_car       [chain, draw]        [1, 2000]
13           prior         b_cost_car       [chain, draw]        [1, 2000]
14           prior       b_cost_train       [chain, draw]        [1, 2000]
15           prior         b_time_car       [chain, draw]        [1, 2000]
16           prior       b_time_train       [chain, draw]        [1, 2000]
17           prior           log_like  [chain, draw, obs]  [1, 2000, 2232]
18    sample_stats    acceptance_rate       [chain, draw]        [4, 2000]
19    sample_stats          diverging       [chain, draw]        [4, 2000]
20    sample_stats             energy       [chain, draw]        [4, 2000]
21    sample_stats                 lp       [chain, draw]        [4, 2000]
22    sample_stats            n_steps       [chain, draw]        [4, 2000]
23    sample_stats          step_size       [chain, draw]        [4, 2000]
24    sample_stats         tree_depth       [chain, draw]        [4, 2000]