2. Logit and sample with weights (Bayesian)

Example of a logit model with a weighted sample

Michel Bierlaire, EPFL Thu Oct 30 2025, 10:15:17

from IPython.core.display_functions import display

from biogeme.bayesian_estimation import BayesianResults, get_pandas_estimated_parameters
from biogeme.biogeme import BIOGEME
from biogeme.expressions import Beta
from biogeme.models import loglogit

See the data processing script: Data preparation for Swissmetro.

from swissmetro_data import (
    CAR_AV_SP,
    CAR_CO_SCALED,
    CAR_TT_SCALED,
    CHOICE,
    GROUP,
    SM_AV,
    SM_COST_SCALED,
    SM_TT_SCALED,
    TRAIN_AV_SP,
    TRAIN_COST_SCALED,
    TRAIN_TT_SCALED,
    database,
)

Parameters to be estimated.

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)

Definition of the utility functions.

v_train = asc_train + b_time * TRAIN_TT_SCALED + b_cost * TRAIN_COST_SCALED
v_swissmetro = asc_sm + b_time * SM_TT_SCALED + b_cost * SM_COST_SCALED
v_car = asc_car + b_time * CAR_TT_SCALED + b_cost * CAR_CO_SCALED

Associate utility functions with the numbering of alternatives.

v = {1: v_train, 2: v_swissmetro, 3: v_car}

Associate the availability conditions with the alternatives.

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

Definition of the model.

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

logprob = loglogit(v, av, CHOICE)

Definition of the weight.

WEIGHT_GROUP_2 = 8.890991e-01
WEIGHT_GROUP_3 = 1.2
weight = WEIGHT_GROUP_2 * (GROUP == 2) + WEIGHT_GROUP_3 * (GROUP == 3)

These notes will be included as such in the report file.

USER_NOTES = (
    'Example of a logit model with three alternatives: '
    'Train, Car and Swissmetro.'
    ' Weighted Exogenous Sample Maximum Likelihood estimator (WESML)'
)

Create the Biogeme object. It is possible to control the generation of the HTML and the yaml files. Note that these parameters can also be modified in the .TOML configuration file.

formulas = {'log_like': logprob, 'weight': weight}
the_biogeme = BIOGEME(
    database,
    formulas,
    mcmc_sampling_strategy='pymc',
    user_notes=USER_NOTES,
    generate_html=True,
    generate_yaml=False,
)
the_biogeme.model_name = 'b02_weight'

## %%
# Estimate the parameters.
try:
    results = BayesianResults.from_netcdf(
        filename=f'saved_results/{the_biogeme.model_name}.nc'
    )
except FileNotFoundError:
    results = the_biogeme.bayesian_estimation()

## %%
print(results.short_summary())

load finished in 4343 ms (4.34 s)
posterior_predictive_loglike finished in 250 ms
expected_log_likelihood finished in 12 ms
best_draw_log_likelihood finished in 11 ms
/Users/bierlair/python_envs/venv313/lib/python3.13/site-packages/arviz/stats/stats.py:1667: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail.
See http://arxiv.org/abs/1507.04544 for details
  warnings.warn(
waic_res finished in 643 ms
waic finished in 643 ms
loo_res finished in 7607 ms (7.61 s)
loo finished in 7607 ms (7.61 s)
Sample size                                              6768
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:01:08.123196
Posterior predictive log-likelihood (sum of log mean p)  -5667.46
Expected log-likelihood E[log L(Y|θ)]                    -5671.10
Best-draw log-likelihood (posterior upper bound)         -5669.09
WAIC (Widely Applicable Information Criterion)           -5674.75
WAIC Standard Error                                      63.97
Effective number of parameters (p_WAIC)                  7.29
LOO (Leave-One-Out Cross-Validation)                     -5674.75
LOO Standard Error                                       63.97
Effective number of parameters (p_LOO)                   7.30

Get the results in a pandas table

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

        Name  Value (mean)  Value (median)  ...     R hat   ESS (bulk)   ESS (tail)
0  asc_train     -0.795184       -0.794658  ...  1.000254  3965.997402  4158.713468
1     b_time     -1.347826       -1.347359  ...  1.000156  3807.078668  4421.493905
2     b_cost     -1.141865       -1.141583  ...  1.000357  5604.456149  5152.311844
3    asc_car     -0.090900       -0.091456  ...  0.999939  4252.924344  4418.206056

[4 rows x 12 columns]