Note
Go to the end to download the full example code.
Investigation of several choice modelsΒΆ
Investigate several choice models:
logit
nested logit with two nests: public and private transportation
nested logit with two nests existing and future modes
for a total of 3 specifications. See Bierlaire and Ortelli (2023).
Michel Bierlaire, EPFL Sun Apr 27 2025, 15:46:15
from IPython.core.display_functions import display
import biogeme.biogeme_logging as blog
from biogeme.biogeme import BIOGEME
from biogeme.catalog import Catalog
from biogeme.data.swissmetro 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,
read_data,
)
from biogeme.expressions import Beta
from biogeme.models import loglogit, lognested
from biogeme.nests import NestsForNestedLogit, OneNestForNestedLogit
from biogeme.results_processing import compile_estimation_results, pareto_optimal
logger = blog.get_screen_logger(level=blog.INFO)
Parameters to be estimated
asc_car = Beta('asc_car', 0, None, None, 0)
asc_train = Beta('asc_train', 0, None, None, 0)
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 = 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 logit model. This is the contribution of each observation to the log likelihood function.
log_probability_logit = loglogit(v, av, CHOICE)
Nested logit model: nest with existing alternatives.
mu_existing = Beta('mu_existing', 1, 1, 10, 0)
existing = OneNestForNestedLogit(
nest_param=mu_existing, list_of_alternatives=[1, 3], name='Existing'
)
nests_existing = NestsForNestedLogit(choice_set=list(v), tuple_of_nests=(existing,))
log_probability_nested_existing = lognested(v, av, nests_existing, CHOICE)
The following elements do not appear in any nest and are assumed each to be alone in a separate nest: {2}. If it is not the intention, check the assignment of alternatives to nests.
Nested logit model: nest with public transportation alternatives.
mu_public = Beta('mu_public', 1, 1, 10, 0)
public = OneNestForNestedLogit(
nest_param=mu_public, list_of_alternatives=[1, 2], name='Public'
)
nests_public = NestsForNestedLogit(choice_set=list(v), tuple_of_nests=(public,))
log_probability_nested_public = lognested(v, av, nests_public, CHOICE)
The following elements do not appear in any nest and are assumed each to be alone in a separate nest: {3}. If it is not the intention, check the assignment of alternatives to nests.
Catalog.
model_catalog = Catalog.from_dict(
catalog_name='model_catalog',
dict_of_expressions={
'logit': log_probability_logit,
'nested existing': log_probability_nested_existing,
'nested public': log_probability_nested_public,
},
)
Read the data
database = read_data()
Create the Biogeme object.
the_biogeme = BIOGEME(database, model_catalog, generate_html=False, generate_yaml=False)
the_biogeme.model_name = 'b01model'
Biogeme parameters read from biogeme.toml.
Estimate the parameters.
dict_of_results = the_biogeme.estimate_catalog()
Estimating 3 models.
Biogeme parameters provided by the user.
*** Initial values of the parameters are obtained from the file __b01model_000000.iter
Parameter values restored from __b01model_000000.iter
Starting values for the algorithm: {'asc_train': -0.37295726674884866, 'b_time': -0.9579800087154183, 'b_cost': -0.6286908253191936, 'mu_existing': 2.0510735367894326, 'asc_car': -0.0013088324487795164}
As the model is not too complex, we activate the calculation of second derivatives. To change this behavior, modify the algorithm to "simple_bounds" in the TOML file.
Optimization algorithm: hybrid Newton/BFGS with simple bounds [simple_bounds]
** Optimization: Newton with trust region for simple bounds
Iter. asc_train b_time b_cost asc_car Function Relgrad Radius Rho
0 -0.62 -1.2 -0.77 -0.0043 8.7e+03 0.0099 10 1.1 ++
1 -0.65 -1.3 -0.79 0.013 8.7e+03 0.0034 1e+02 0.93 ++
2 -0.65 -1.3 -0.79 0.013 8.7e+03 0.0034 0.0088 -0.0093 -
3 -0.66 -1.3 -0.8 0.022 8.7e+03 0.0023 0.0088 0.46 +
4 -0.67 -1.3 -0.79 0.03 8.7e+03 0.002 0.0088 0.13 +
5 -0.67 -1.3 -0.79 0.03 8.7e+03 0.002 0.0044 2.3e-07 -
6 -0.67 -1.3 -0.79 0.03 8.7e+03 0.002 0.0022 3.3e-07 -
7 -0.67 -1.3 -0.79 0.03 8.7e+03 0.002 0.0011 5.7e-07 -
8 -0.67 -1.3 -0.79 0.03 8.7e+03 0.002 0.00055 1.1e-06 -
9 -0.67 -1.3 -0.79 0.03 8.7e+03 0.002 0.00027 2.1e-06 -
10 -0.67 -1.3 -0.79 0.03 8.7e+03 0.002 0.00014 4.2e-06 -
11 -0.67 -1.3 -0.79 0.03 8.7e+03 0.002 6.9e-05 8.3e-06 -
12 -0.67 -1.3 -0.79 0.03 8.7e+03 0.002 3.4e-05 1.6e-05 -
13 -0.67 -1.3 -0.79 0.03 8.7e+03 0.002 1.7e-05 3.3e-05 -
14 -0.67 -1.3 -0.79 0.03 8.7e+03 0.002 8.6e-06 6.6e-05 -
15 -0.67 -1.3 -0.79 0.03 8.7e+03 0.002 4.3e-06 0.00013 -
16 -0.67 -1.3 -0.79 0.03 8.7e+03 0.002 2.1e-06 0.00026 -
17 -0.67 -1.3 -0.79 0.03 8.7e+03 0.002 1.1e-06 0.00042 -
18 -0.67 -1.3 -0.79 0.03 8.7e+03 0.002 5.4e-07 0.00053 -
19 -0.67 -1.3 -0.79 0.03 8.7e+03 0.002 2.7e-07 0.00059 -
20 -0.67 -1.3 -0.79 0.03 8.7e+03 0.002 1.3e-07 0.00062 -
21 -0.67 -1.3 -0.79 0.03 8.7e+03 0.002 6.7e-08 0.00064 -
22 -0.67 -1.3 -0.79 0.03 8.7e+03 0.002 3.3e-08 0.00064 -
23 -0.67 -1.3 -0.79 0.03 8.7e+03 0.002 1.7e-08 0.00065 -
24 -0.67 -1.3 -0.79 0.03 8.7e+03 0.002 8.4e-09 0.00065 -
Optimization algorithm has *not* converged.
Algorithm: Newton with trust region for simple bound constraints
Cause of termination: Trust region is too small: 8.374050017648238e-09
Number of iterations: 25
Proportion of Hessian calculation: 5/5 = 100.0%
Optimization time: 0:00:00.313746
Calculate second derivatives and BHHH
It seems that the optimization algorithm did not converge. Therefore, the results may not correspond to the maximum likelihood estimator. Check the specification of the model, or the criteria for convergence of the algorithm.
Biogeme parameters provided by the user.
*** Initial values of the parameters are obtained from the file __b01model_000001.iter
Parameter values restored from __b01model_000001.iter
Starting values for the algorithm: {'asc_train': -0.652238664271019, 'b_time': -1.2789413398819158, 'b_cost': -0.7897904566401142, 'asc_car': 0.01622793815045202}
Optimization algorithm: hybrid Newton/BFGS with simple bounds [simple_bounds]
** Optimization: Newton with trust region for simple bounds
Iter. asc_train b_time b_cost mu_existing asc_car Function Relgrad Radius Rho
0 -1.2 -2.3 -0.65 1 -0.66 2.2e+04 0.042 10 3.2 ++
1 0.91 -0.29 -1 1 -0.47 1.1e+04 0.0018 10 0.76 +
2 0.91 -0.29 -1 1 -0.47 1.1e+04 0.0018 0.36 -0.81 -
3 0.91 -0.29 -1 1 -0.47 1.1e+04 0.0018 0.18 -1.7 -
4 0.91 -0.29 -1 1 -0.47 1.1e+04 0.0018 0.089 -0.33 -
5 0.99 -0.22 -0.94 1 -0.56 1e+04 0.0036 0.089 0.22 +
6 0.99 -0.22 -0.94 1 -0.56 1e+04 0.0036 0.045 -1.2 -
7 0.99 -0.22 -0.94 1 -0.56 1e+04 0.0036 0.022 -0.6 -
8 0.99 -0.22 -0.94 1 -0.56 1e+04 0.0036 0.011 -0.47 -
9 0.99 -0.22 -0.94 1 -0.56 1e+04 0.0036 0.0056 -0.42 -
10 0.99 -0.22 -0.94 1 -0.56 1e+04 0.0036 0.0028 -0.39 -
11 0.99 -0.22 -0.94 1 -0.56 1e+04 0.0036 0.0014 -0.38 -
12 0.99 -0.22 -0.94 1 -0.56 1e+04 0.0036 0.0007 -0.37 -
13 0.99 -0.22 -0.94 1 -0.56 1e+04 0.0036 0.00035 -0.37 -
14 0.99 -0.22 -0.94 1 -0.56 1e+04 0.0036 0.00017 -0.37 -
15 0.99 -0.22 -0.94 1 -0.56 1e+04 0.0036 8.7e-05 -0.37 -
16 0.99 -0.22 -0.94 1 -0.56 1e+04 0.0036 4.4e-05 -0.37 -
17 0.99 -0.22 -0.94 1 -0.56 1e+04 0.0036 2.2e-05 -0.37 -
18 0.99 -0.22 -0.94 1 -0.56 1e+04 0.0036 1.1e-05 -0.37 -
19 0.99 -0.22 -0.94 1 -0.56 1e+04 0.0036 5.4e-06 -0.37 -
20 0.99 -0.22 -0.94 1 -0.56 1e+04 0.0036 2.7e-06 -0.37 -
21 0.99 -0.22 -0.94 1 -0.56 1e+04 0.0036 1.4e-06 -0.37 -
22 0.99 -0.22 -0.94 1 -0.56 1e+04 0.0036 6.8e-07 -0.37 -
23 0.99 -0.22 -0.94 1 -0.56 1e+04 0.0036 3.4e-07 -0.37 -
24 0.99 -0.22 -0.94 1 -0.56 1e+04 0.0036 1.7e-07 -0.37 -
25 0.99 -0.22 -0.94 1 -0.56 1e+04 0.0036 8.5e-08 -0.37 -
26 0.99 -0.22 -0.94 1 -0.56 1e+04 0.0036 4.3e-08 -0.37 -
27 0.99 -0.22 -0.94 1 -0.56 1e+04 0.0036 2.1e-08 -0.37 -
28 0.99 -0.22 -0.94 1 -0.56 1e+04 0.0036 1.1e-08 -0.37 -
Optimization algorithm has *not* converged.
Algorithm: Newton with trust region for simple bound constraints
Cause of termination: Trust region is too small: 1.0625215485227102e-08
Number of iterations: 29
Proportion of Hessian calculation: 4/4 = 100.0%
Optimization time: 0:00:01.224185
Calculate second derivatives and BHHH
It seems that the optimization algorithm did not converge. Therefore, the results may not correspond to the maximum likelihood estimator. Check the specification of the model, or the criteria for convergence of the algorithm.
Biogeme parameters provided by the user.
*** Initial values of the parameters are obtained from the file __b01model_000002.iter
Parameter values restored from __b01model_000002.iter
Starting values for the algorithm: {'asc_train': -0.5667948673966152, 'b_time': -1.2528443783265466, 'b_cost': -0.7817835231149908, 'mu_public': 1.0817400108564705, 'asc_car': -0.009746105754462168}
Optimization algorithm: hybrid Newton/BFGS with simple bounds [simple_bounds]
** Optimization: Newton with trust region for simple bounds
Optimization algorithm has converged.
Relative gradient: 2.7150579439738833e-08
Cause of termination: Relative gradient = 2.7e-08 <= 6.1e-06
Number of function evaluations: 1
Number of gradient evaluations: 1
Number of hessian evaluations: 0
Algorithm: Newton with trust region for simple bound constraints
Number of iterations: 0
Optimization time: 0:00:00.338178
Calculate second derivatives and BHHH
Number of estimated models.
print(f'A total of {len(dict_of_results)} models have been estimated')
A total of 3 models have been estimated
All estimation results
compiled_results, specs = compile_estimation_results(
dict_of_results, use_short_names=True
)
display('All estimated models')
display(compiled_results)
All estimated models
Model_000000 ... Model_000002
Number of estimated parameters 4 ... 5
Sample size 10719 ... 10719
Final log likelihood -8670.636 ... -8669.497
Akaike Information Criterion 17349.27 ... 17348.99
Bayesian Information Criterion 17378.39 ... 17385.39
asc_train (t-test) -0.665 (-12.3) ... -0.567 (-6.87)
b_time (t-test) -1.27 (-19.5) ... -1.25 (-16.2)
b_cost (t-test) -0.793 (-15.5) ... -0.782 (-15)
asc_car (t-test) 0.0304 (0.821) ... -0.00975 (-0.195)
mu_existing (t-test) ...
mu_public (t-test) ... 1.08 (12.5)
[11 rows x 3 columns]
Glossary
for short_name, spec in specs.items():
print(f'{short_name}\t{spec}')
Model_000000 model_catalog:logit
Model_000001 model_catalog:nested existing
Model_000002 model_catalog:nested public
Estimation results of the Pareto optimal models.
pareto_results = pareto_optimal(dict_of_results)
compiled_pareto_results, pareto_specs = compile_estimation_results(
pareto_results, use_short_names=True
)
No Pareto file has been provided
display('Non dominated models')
display(compiled_pareto_results)
Non dominated models
Model_000000 Model_000001
Number of estimated parameters 4 5
Sample size 10719 10719
Final log likelihood -8670.636 -8669.497
Akaike Information Criterion 17349.27 17348.99
Bayesian Information Criterion 17378.39 17385.39
asc_train (t-test) -0.665 (-12.3) -0.567 (-6.87)
b_time (t-test) -1.27 (-19.5) -1.25 (-16.2)
b_cost (t-test) -0.793 (-15.5) -0.782 (-15)
asc_car (t-test) 0.0304 (0.821) -0.00975 (-0.195)
mu_public (t-test) 1.08 (12.5)
Glossary.
for short_name, spec in pareto_specs.items():
print(f'{short_name}\t{spec}')
Model_000000 model_catalog:logit
Model_000001 model_catalog:nested public
Total running time of the script: (0 minutes 5.361 seconds)