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.
No YAML file found at b01model_000000.yaml. Estimation is performed.
*** Initial values of the parameters are obtained from the file __b01model_000000.iter
Cannot read file __b01model_000000.iter. Statement is ignored.
Starting values for the algorithm: {}
Analytical Hessian method: full
As the model is rather complex, we cancel the calculation of second derivatives. If you want to control the parameters, change the algorithm from "automatic" to "simple_bounds" in the TOML file.
Optimization algorithm: hybrid Newton/BFGS with simple bounds [simple_bounds]
** Optimization: BFGS with trust region for simple bounds
Iter. asc_train b_time b_cost asc_car mu_existing Function Relgrad Radius Rho
0 0 0 0 0 1 1.1e+04 0.26 0.5 0.0044 -
1 -0.5 -0.5 -0.5 0.5 1.5 9.2e+03 0.14 0.5 0.51 +
2 -0.5 -1 1.1e-16 0 2 8.9e+03 0.094 0.5 0.25 +
3 -0.5 -1 1.1e-16 0 2 8.9e+03 0.094 0.25 -0.53 -
4 -0.25 -0.75 -0.25 0.009 1.9 8.7e+03 0.056 0.25 0.5 +
5 -0.5 -1 -0.5 0.21 2 8.7e+03 0.074 0.25 0.12 +
6 -0.32 -1 -0.6 -0.045 2 8.5e+03 0.027 0.25 0.65 +
7 -0.32 -1 -0.6 -0.045 2 8.5e+03 0.027 0.12 -0.64 -
8 -0.38 -0.97 -0.73 0.037 2 8.5e+03 0.011 0.12 0.32 +
9 -0.38 -0.97 -0.73 0.037 2 8.5e+03 0.011 0.062 -0.45 -
10 -0.39 -1 -0.67 -0.0094 2 8.5e+03 0.0041 0.062 0.74 +
11 -0.39 -1 -0.67 -0.0094 2 8.5e+03 0.0041 0.031 -0.33 -
12 -0.38 -0.98 -0.65 0.022 2 8.5e+03 0.0052 0.031 0.34 +
13 -0.38 -0.98 -0.65 0.022 2 8.5e+03 0.0052 0.016 -0.61 -
14 -0.37 -0.99 -0.64 0.0063 2 8.5e+03 0.0021 0.016 0.54 +
15 -0.37 -0.98 -0.64 0.0077 2 8.5e+03 0.0011 0.016 0.9 +
16 -0.37 -0.98 -0.64 0.0077 2 8.5e+03 0.0011 0.0078 -0.23 -
17 -0.38 -0.97 -0.63 -9.1e-05 2 8.5e+03 0.00085 0.0078 0.35 +
18 -0.38 -0.97 -0.64 0.0015 2 8.5e+03 0.0014 0.0078 0.29 +
19 -0.38 -0.97 -0.64 0.0019 2 8.5e+03 0.00048 0.0078 0.78 +
20 -0.37 -0.97 -0.63 0.00056 2 8.5e+03 0.00085 0.0078 0.69 +
21 -0.38 -0.97 -0.63 0.0021 2 8.5e+03 0.00051 0.0078 0.61 +
22 -0.37 -0.96 -0.63 -0.0007 2 8.5e+03 0.00038 0.0078 0.43 +
23 -0.37 -0.96 -0.63 -0.0007 2 8.5e+03 0.00038 0.0039 -0.039 -
24 -0.37 -0.96 -0.63 -0.00043 2 8.5e+03 0.00022 0.0039 0.74 +
25 -0.37 -0.96 -0.63 -0.00043 2 8.5e+03 0.00022 0.002 -0.053 -
26 -0.37 -0.96 -0.63 8.8e-05 2 8.5e+03 0.00036 0.002 0.43 +
27 -0.37 -0.96 -0.63 -0.00023 2 8.5e+03 0.00014 0.002 0.86 +
28 -0.37 -0.96 -0.63 -0.0013 2 8.5e+03 0.00032 0.002 0.18 +
29 -0.37 -0.96 -0.63 -0.00068 2 8.5e+03 0.00011 0.002 0.9 +
30 -0.37 -0.96 -0.63 -0.00083 2 8.5e+03 0.00017 0.002 0.69 +
31 -0.37 -0.96 -0.63 -0.0014 2 8.5e+03 0.00012 0.002 0.29 +
32 -0.37 -0.96 -0.63 -0.0014 2 8.5e+03 0.00011 0.002 0.31 +
33 -0.37 -0.96 -0.63 -0.0014 2 8.5e+03 0.00011 0.00098 0.048 -
34 -0.37 -0.96 -0.63 -0.0011 2 8.5e+03 4.9e-05 0.00098 0.79 +
35 -0.37 -0.96 -0.63 -0.0011 2 8.5e+03 4.9e-05 0.00049 -3.5 -
36 -0.37 -0.96 -0.63 -0.0011 2 8.5e+03 4.9e-05 0.00024 -0.69 -
37 -0.37 -0.96 -0.63 -0.0012 2 8.5e+03 4.9e-05 0.00024 0.16 +
38 -0.37 -0.96 -0.63 -0.0014 2.1 8.5e+03 3.9e-05 0.00024 0.22 +
39 -0.37 -0.96 -0.63 -0.0011 2.1 8.5e+03 5.1e-05 0.00024 0.26 +
40 -0.37 -0.96 -0.63 -0.0013 2.1 8.5e+03 1.4e-05 0.00024 0.87 +
41 -0.37 -0.96 -0.63 -0.0013 2.1 8.5e+03 4.9e-06 0.00024 0.87 +
Optimization algorithm has converged.
Relative gradient: 4.8911388627737665e-06
Cause of termination: Relative gradient = 4.9e-06 <= 6.1e-06
Number of function evaluations: 103
Number of gradient evaluations: 61
Number of hessian evaluations: 0
Algorithm: BFGS with trust region for simple bound constraints
Number of iterations: 42
Proportion of Hessian calculation: 0/30 = 0.0%
Optimization time: 0:00:00.214547
Calculate final gradient and BHHH
Calculate second derivatives
Biogeme parameters provided by the user.
No YAML file found at b01model_000001.yaml. Estimation is performed.
*** Initial values of the parameters are obtained from the file __b01model_000001.iter
Cannot read file __b01model_000001.iter. Statement is ignored.
Starting values for the algorithm: {}
Analytical Hessian method: full
As the model is rather complex, we cancel the calculation of second derivatives. If you want to control the parameters, change the algorithm from "automatic" to "simple_bounds" in the TOML file.
Optimization algorithm: hybrid Newton/BFGS with simple bounds [simple_bounds]
** Optimization: BFGS with trust region for simple bounds
Iter. asc_train b_time b_cost asc_car mu_public Function Relgrad Radius Rho
0 0 0 0 0 1 1.1e+04 0.26 0.5 -0.063 -
1 -0.5 -0.5 -0.5 0.5 1.5 9.6e+03 0.15 0.5 0.53 +
2 -0.65 -1 0 0 1 9e+03 0.068 0.5 0.36 +
3 -0.65 -1 0 0 1 9e+03 0.068 0.25 -1.5 -
4 -0.65 -1 0 0 1 9e+03 0.068 0.12 -0.28 -
5 -0.77 -1.1 -0.12 0.076 1.1 8.9e+03 0.059 0.12 0.37 +
6 -0.65 -1.1 -0.25 0.068 1 8.8e+03 0.043 0.12 0.82 +
7 -0.52 -1.2 -0.38 0.034 1.1 8.7e+03 0.034 0.12 0.77 +
8 -0.45 -1.2 -0.5 0.16 1.2 8.7e+03 0.026 0.12 0.36 +
9 -0.4 -1.3 -0.62 0.034 1.1 8.7e+03 0.012 1.2 0.92 ++
10 -0.4 -1.3 -0.62 0.034 1.1 8.7e+03 0.012 0.62 -3.9 -
11 -0.4 -1.3 -0.62 0.034 1.1 8.7e+03 0.012 0.31 -1.7 -
12 -0.4 -1.3 -0.62 0.034 1.1 8.7e+03 0.012 0.16 0.01 -
13 -0.49 -1.3 -0.78 0.017 1.1 8.7e+03 0.0043 0.16 0.75 +
14 -0.49 -1.3 -0.78 0.017 1.1 8.7e+03 0.0043 0.078 -1 -
15 -0.49 -1.3 -0.78 0.017 1.1 8.7e+03 0.0043 0.039 -0.04 -
16 -0.5 -1.3 -0.77 -0.022 1.1 8.7e+03 0.0019 0.039 0.64 +
17 -0.54 -1.2 -0.78 -0.017 1.1 8.7e+03 0.0022 0.039 0.25 +
18 -0.54 -1.2 -0.78 -0.017 1.1 8.7e+03 0.0022 0.02 -0.15 -
19 -0.56 -1.2 -0.77 -0.031 1.1 8.7e+03 0.0018 0.02 0.14 +
20 -0.56 -1.2 -0.78 -0.011 1.1 8.7e+03 0.0004 0.02 0.8 +
21 -0.56 -1.2 -0.78 -0.011 1.1 8.7e+03 0.0004 0.0098 -0.47 -
22 -0.56 -1.3 -0.78 -0.0094 1.1 8.7e+03 0.00055 0.0098 0.26 +
23 -0.57 -1.2 -0.78 -0.012 1.1 8.7e+03 0.0003 0.0098 0.12 +
24 -0.57 -1.2 -0.78 -0.012 1.1 8.7e+03 0.0003 0.0049 -1.5 -
25 -0.57 -1.3 -0.78 -0.0095 1.1 8.7e+03 8.4e-05 0.0049 0.88 +
26 -0.57 -1.3 -0.78 -0.0095 1.1 8.7e+03 8.4e-05 0.0024 -5.2 -
27 -0.57 -1.3 -0.78 -0.0095 1.1 8.7e+03 8.4e-05 0.0012 -3.5 -
28 -0.57 -1.3 -0.78 -0.0095 1.1 8.7e+03 8.4e-05 0.00061 -0.71 -
29 -0.57 -1.3 -0.78 -0.01 1.1 8.7e+03 7.5e-05 0.00061 0.46 +
30 -0.57 -1.3 -0.78 -0.0095 1.1 8.7e+03 6e-05 0.00061 0.53 +
31 -0.57 -1.3 -0.78 -0.0096 1.1 8.7e+03 1.2e-05 0.00061 0.84 +
32 -0.57 -1.3 -0.78 -0.0096 1.1 8.7e+03 1.2e-05 0.00031 -0.85 -
33 -0.57 -1.3 -0.78 -0.0096 1.1 8.7e+03 1.2e-05 0.00015 0.099 -
34 -0.57 -1.3 -0.78 -0.0096 1.1 8.7e+03 6e-06 0.00015 0.81 -
Optimization algorithm has converged.
Relative gradient: 6.018995530655885e-06
Cause of termination: Relative gradient = 6e-06 <= 6.1e-06
Number of function evaluations: 74
Number of gradient evaluations: 39
Number of hessian evaluations: 0
Algorithm: BFGS with trust region for simple bound constraints
Number of iterations: 35
Proportion of Hessian calculation: 0/19 = 0.0%
Optimization time: 0:00:00.182564
Calculate final gradient and BHHH
Calculate second derivatives
Biogeme parameters provided by the user.
No YAML file found at b01model_000002.yaml. Estimation is performed.
*** Initial values of the parameters are obtained from the file __b01model_000002.iter
Cannot read file __b01model_000002.iter. Statement is ignored.
Starting values for the algorithm: {}
Analytical Hessian method: full
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.76 -0.77 -0.7 -0.29 8.8e+03 0.04 10 1.1 ++
1 -0.66 -1.2 -0.77 -0.0015 8.7e+03 0.0064 1e+02 1.1 ++
2 -0.65 -1.3 -0.79 0.016 8.7e+03 0.00012 1e+03 1 ++
3 -0.65 -1.3 -0.79 0.016 8.7e+03 4e-08 1e+03 1 ++
Optimization algorithm has converged.
Relative gradient: 3.954408090874478e-08
Cause of termination: Relative gradient = 4e-08 <= 6.1e-06
Number of function evaluations: 13
Number of gradient evaluations: 9
Number of hessian evaluations: 4
Algorithm: Newton with trust region for simple bound constraints
Number of iterations: 4
Proportion of Hessian calculation: 4/4 = 100.0%
Optimization time: 0:00:00.302128
Calculate final gradient and BHHH
Calculate second derivatives
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 5 ... 4
Sample size 10719 ... 10719
Final log likelihood -8526.89 ... -8670.163
Akaike Information Criterion 17063.78 ... 17348.33
Bayesian Information Criterion 17100.18 ... 17377.45
asc_train (t-test) -0.373 (-7.17) ... -0.652 (-12)
b_time (t-test) -0.958 (-14.7) ... -1.28 (-19.5)
b_cost (t-test) -0.629 (-14.8) ... -0.79 (-15.5)
asc_car (t-test) -0.00128 (-0.0375) ... 0.0162 (0.438)
mu_existing (t-test) 2.05 (15.8) ...
mu_public (t-test) ...
[11 rows x 3 columns]
Glossary
for short_name, spec in specs.items():
print(f'{short_name}\t{spec}')
Model_000000 model_catalog:nested existing
Model_000001 model_catalog:nested public
Model_000002 model_catalog:logit
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.163 -8526.89
Akaike Information Criterion 17348.33 17063.78
Bayesian Information Criterion 17377.45 17100.18
asc_train (t-test) -0.652 (-12) -0.373 (-7.17)
b_time (t-test) -1.28 (-19.5) -0.958 (-14.7)
b_cost (t-test) -0.79 (-15.5) -0.629 (-14.8)
asc_car (t-test) 0.0162 (0.438) -0.00128 (-0.0375)
mu_existing (t-test) 2.05 (15.8)
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 existing
Total running time of the script: (0 minutes 2.341 seconds)