17. Mixture with lognormal distribution

Bayesian estimation of a mixture of logit models. The mixing distribution is distributed as a log normal.

Michel Bierlaire, EPFL Sat Nov 15 2025, 18:20:02

from IPython.core.display_functions import display

import biogeme.biogeme_logging as blog
from biogeme.bayesian_estimation import (
    BayesianResults,
    get_pandas_estimated_parameters,
)
from biogeme.biogeme import BIOGEME
from biogeme.expressions import Beta, DistributedParameter, Draws, exp
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,
    SM_AV,
    SM_COST_SCALED,
    SM_TT_SCALED,
    TRAIN_AV_SP,
    TRAIN_COST_SCALED,
    TRAIN_TT_SCALED,
    database,
)

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

Example b17_lognormal_mixture.py

The scale parameters must stay away from zero. We define a small but positive lower bound

POSITIVE_LOWER_BOUND = 1.0e-5

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_cost = Beta('b_cost', 0, None, None, 0)

Define a random parameter, normally distributed, designed to be used for Monte-Carlo simulation.

b_time_param = Beta('b_time', 0, None, None, 0)

It is advised not to use 0 as starting value for the following parameter.

b_time_s = Beta('b_time_s', 1, POSITIVE_LOWER_BOUND, 2, 0)

Define a random parameter, log normally distributed, designed to be used for Monte-Carlo simulation.

b_time_eps = Draws('b_time_eps', 'NORMAL')
b_time_rnd = DistributedParameter(
    'b_time_rnd', -exp(b_time_param + b_time_s * b_time_eps)
)

Definition of the utility functions.

v_train = asc_train + b_time_rnd * TRAIN_TT_SCALED + b_cost * TRAIN_COST_SCALED
v_swissmetro = asc_sm + b_time_rnd * SM_TT_SCALED + b_cost * SM_COST_SCALED
v_car = asc_car + b_time_rnd * 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}

Conditional to b_time_rnd, we have a logit model (called the kernel).

conditional_log_probability = loglogit(v, av, CHOICE)

%% Create the Biogeme object.

the_biogeme = BIOGEME(database, conditional_log_probability)
the_biogeme.model_name = 'b17_lognormal_mixture'

Biogeme parameters read from biogeme.toml.

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())

Loaded NetCDF file size: 1.8 GB
load finished in 9325 ms (9.33 s)
posterior_predictive_loglike finished in 258 ms
expected_log_likelihood finished in 11 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 635 ms
waic finished in 635 ms
/Users/bierlair/python_envs/venv313/lib/python3.13/site-packages/arviz/stats/stats.py:797: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.70 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
  warnings.warn(
loo_res finished in 10440 ms (10.44 s)
loo finished in 10440 ms (10.44 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:03:48.250586
Posterior predictive log-likelihood (sum of log mean p)  -4060.54
Expected log-likelihood E[log L(Y|θ)]                    -4411.29
Best-draw log-likelihood (posterior upper bound)         -4112.89
WAIC (Widely Applicable Information Criterion)           -4934.26
WAIC Standard Error                                      47.71
Effective number of parameters (p_WAIC)                  873.71
LOO (Leave-One-Out Cross-Validation)                     -5081.00
LOO Standard Error                                       50.18
Effective number of parameters (p_LOO)                   1020.46

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

Diagnostics computation took 72.9 seconds (cached).
        Name  Value (mean)  Value (median)  ...     R hat   ESS (bulk)   ESS (tail)
0  asc_train     -0.346957       -0.346935  ...  1.000638  2334.277646  4475.210598
1     b_time      0.574283        0.575513  ...  1.001543  1682.584504  3600.557096
2     b_cost     -1.385609       -1.383394  ...  1.000515  2427.918679  4159.723448
3    asc_car      0.176598        0.176498  ...  1.001084  1810.054702  4230.176761
4   b_time_s      1.261262        1.258514  ...  1.004316   854.955689  1864.878002

[5 rows x 12 columns]