3. Moneymetric and heteroscedastic specification

Although normalizing the scale to 1 is a common practice in random utility models, it is sometimes preferable to normalize another parameter. For instance, normalizing the cost coefficient to -1 sets the units of the utility function as currency units (CHF here), and the estimated coefficients are easily interpreted as willingness to pay. In that case, the scale must be estimated.

We also illustrate here a heteroscedastic specification, where a different scale is associated with different segments of the sample.

This example illustrates how to estimate such a specification with Bayesian inference.

Michel Bierlaire, EPFL Thu Nov 20 2025, 11:10:03

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

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_time = Beta('b_time', 0, None, 0, 0)
b_cost = Beta('b_cost', -1, None, None, 1)
scale_not_group3 = Beta('scale_not_group3', 1, POSITIVE_LOWER_BOUND, None, 0)
scale_group3 = Beta('scale_group3', 1, POSITIVE_LOWER_BOUND, 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

Scale associated with group 3 is estimated.

scale = (GROUP != 3) * scale_not_group3 + (GROUP == 3) * scale_group3

Scale the utility functions, and associate them with the numbering of alternatives.

v = {1: scale * v_train, 2: scale * v_swissmetro, 3: scale * 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)

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

USER_NOTES = (
    'Illustrates a moneymetric heteroscedastic specification. A different scale is'
    ' associated with different segments of the sample.'
)

Create the Biogeme object.

the_biogeme = BIOGEME(database, logprob, user_notes=USER_NOTES)
the_biogeme.model_name = 'b03_scale'

Estimate the parameters.

try:
    results = BayesianResults.from_netcdf(
        filename=f'saved_results/{the_biogeme.model_name}.nc'
    )
except FileNotFoundError:
    results = the_biogeme.bayesian_estimation()
load finished in 4538 ms (4.54 s)
print(results.short_summary())
posterior_predictive_loglike finished in 262 ms
expected_log_likelihood finished in 11 ms
best_draw_log_likelihood finished in 11 ms
waic_res finished in 704 ms
waic finished in 704 ms
loo_res finished in 7974 ms (7.97 s)
loo finished in 7974 ms (7.97 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:00:46.865756
Posterior predictive log-likelihood (sum of log mean p)  -4975.45
Expected log-likelihood E[log L(Y|θ)]                    -4979.15
Best-draw log-likelihood (posterior upper bound)         -4976.71
WAIC (Widely Applicable Information Criterion)           -4982.88
WAIC Standard Error                                      53.94
Effective number of parameters (p_WAIC)                  7.42
LOO (Leave-One-Out Cross-Validation)                     -4982.88
LOO Standard Error                                       53.94
Effective number of parameters (p_LOO)                   7.43

Get the results in a pandas table

pandas_results = get_pandas_estimated_parameters(estimation_results=results)
display(pandas_results)
               Name  Value (mean)  ...   ESS (bulk)   ESS (tail)
0         asc_train     -1.254377  ...  3542.313400  4327.622563
1           asc_car     -0.042057  ...  3850.083378  4634.164946
2  scale_not_group3      0.357689  ...  3672.097395  4144.399754
3      scale_group3      1.490217  ...  3513.553977  4356.330711
4            b_time     -1.051782  ...  3508.996986  4312.974517

[5 rows x 12 columns]

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

Gallery generated by Sphinx-Gallery