Note
Go to the end to download the full example code.
Estimation and simulation of a nested logit model¶
We estimate a nested logit model, and we perform simulation using the estimated model.
Michel Bierlaire, EPFL Sat Jun 28 2025, 16:08:07
from IPython.core.display_functions import display
import biogeme.biogeme_logging as blog
from biogeme.biogeme import BIOGEME
from biogeme.calculator import get_value_c
from biogeme.data.optima import read_data
from biogeme.models import lognested
from biogeme.results_processing import get_pandas_estimated_parameters
from scenarios import scenario
logger = blog.get_screen_logger(level=blog.INFO)
logger.info('Example plot_b02estimation')
Example plot_b02estimation
Obtain the specification for the default scenario. The definition of the scenarios is available in Specification of a nested logit model.
V, nests, choice, _ = scenario()
The choice model is a nested logit, with availability conditions For estimation, we need the log of the probability.
log_probability = lognested(util=V, availability=None, nests=nests, choice=choice)
Get the database
database = read_data()
Create the Biogeme object for estimation.
the_biogeme = BIOGEME(database, log_probability)
the_biogeme.model_name = 'b02estimation'
Biogeme parameters read from biogeme.toml.
Estimate the parameters. Perform bootstrapping.
results = the_biogeme.estimate(run_bootstrap=True)
*** Initial values of the parameters are obtained from the file __b02estimation.iter
Parameter values restored from __b02estimation.iter
Starting values for the algorithm: {'beta_time_fulltime': -1.5964670520847424, 'beta_time_other': -0.5526862214775381, 'beta_cost': -0.719072321557382, 'mu_no_car': 1.5237160520807354, 'asc_sm': 0.06300034044143561, 'beta_dist_male': -0.6876768914763282, 'beta_dist_female': -0.8321525469977931, 'beta_dist_unreported': -0.7047372090923671, 'asc_car': 0.2586072914689233}
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
Optimization algorithm has converged.
Relative gradient: 4.478787282617695e-06
Cause of termination: Relative gradient = 4.5e-06 <= 6.1e-06
Number of function evaluations: 1
Number of gradient evaluations: 1
Number of hessian evaluations: 0
Algorithm: BFGS with trust region for simple bound constraints
Number of iterations: 0
Optimization time: 0:00:00.466837
Calculate second derivatives and BHHH
Re-estimate the model 100 times for bootstrapping
Bootstraps: 0%| | 0/100 [00:00<?, ?it/s]
0%| | 0/100 [00:00<?, ?it/s]
1%| | 1/100 [00:03<04:57, 3.01s/it]
3%|▎ | 3/100 [00:04<02:15, 1.40s/it]
5%|▌ | 5/100 [00:06<01:48, 1.15s/it]
7%|▋ | 7/100 [15:18<5:00:22, 193.79s/it]
9%|▉ | 9/100 [15:20<3:03:42, 121.12s/it]/Users/bierlair/python_envs/venv313/lib/python3.13/site-packages/joblib/externals/loky/process_executor.py:782: UserWarning: A worker stopped while some jobs were given to the executor. This can be caused by a too short worker timeout or by a memory leak.
warnings.warn(
11%|█ | 11/100 [15:24<1:57:42, 79.35s/it]
12%|█▏ | 12/100 [15:24<1:33:13, 63.56s/it]
13%|█▎ | 13/100 [15:25<1:12:06, 49.73s/it]
14%|█▍ | 14/100 [15:26<54:11, 37.81s/it]
15%|█▌ | 15/100 [15:27<40:20, 28.48s/it]
16%|█▌ | 16/100 [15:27<29:16, 20.91s/it]
17%|█▋ | 17/100 [15:29<21:27, 15.51s/it]
18%|█▊ | 18/100 [15:29<15:15, 11.16s/it]
19%|█▉ | 19/100 [15:30<11:19, 8.39s/it]
20%|██ | 20/100 [15:31<08:05, 6.07s/it]
21%|██ | 21/100 [15:34<06:44, 5.12s/it]
22%|██▏ | 22/100 [15:34<04:49, 3.71s/it]
23%|██▎ | 23/100 [15:35<03:51, 3.01s/it]
24%|██▍ | 24/100 [15:36<02:47, 2.21s/it]
25%|██▌ | 25/100 [15:37<02:28, 1.98s/it]
26%|██▌ | 26/100 [15:37<01:47, 1.45s/it]
27%|██▋ | 27/100 [15:39<01:44, 1.43s/it]
28%|██▊ | 28/100 [15:39<01:17, 1.08s/it]
29%|██▉ | 29/100 [15:40<01:23, 1.18s/it]
30%|███ | 30/100 [15:41<01:06, 1.05it/s]
31%|███ | 31/100 [15:44<01:43, 1.49s/it]
32%|███▏ | 32/100 [15:44<01:16, 1.12s/it]
33%|███▎ | 33/100 [15:45<01:20, 1.21s/it]
34%|███▍ | 34/100 [15:45<00:59, 1.10it/s]
35%|███▌ | 35/100 [15:47<01:11, 1.09s/it]
37%|███▋ | 37/100 [15:49<01:00, 1.04it/s]
39%|███▉ | 39/100 [15:50<00:56, 1.08it/s]
41%|████ | 41/100 [15:53<01:08, 1.15s/it]
42%|████▏ | 42/100 [15:54<00:56, 1.02it/s]
43%|████▎ | 43/100 [15:55<01:00, 1.06s/it]
44%|████▍ | 44/100 [15:55<00:47, 1.17it/s]
45%|████▌ | 45/100 [15:57<00:54, 1.00it/s]
46%|████▌ | 46/100 [15:57<00:43, 1.25it/s]
47%|████▋ | 47/100 [15:58<00:50, 1.05it/s]
48%|████▊ | 48/100 [15:59<00:39, 1.33it/s]
49%|████▉ | 49/100 [16:00<00:48, 1.06it/s]
50%|█████ | 50/100 [16:00<00:39, 1.26it/s]
51%|█████ | 51/100 [16:03<01:05, 1.33s/it]
52%|█████▏ | 52/100 [16:03<00:50, 1.06s/it]
53%|█████▎ | 53/100 [16:05<00:53, 1.13s/it]
54%|█████▍ | 54/100 [16:05<00:40, 1.13it/s]
55%|█████▌ | 55/100 [16:06<00:46, 1.04s/it]
56%|█████▌ | 56/100 [16:07<00:36, 1.22it/s]
57%|█████▋ | 57/100 [16:08<00:39, 1.08it/s]
58%|█████▊ | 58/100 [16:08<00:32, 1.29it/s]
59%|█████▉ | 59/100 [16:10<00:37, 1.10it/s]
60%|██████ | 60/100 [16:10<00:30, 1.29it/s]
61%|██████ | 61/100 [16:13<00:51, 1.32s/it]
62%|██████▏ | 62/100 [16:13<00:41, 1.10s/it]
63%|██████▎ | 63/100 [16:14<00:39, 1.08s/it]
64%|██████▍ | 64/100 [16:15<00:33, 1.08it/s]
65%|██████▌ | 65/100 [16:16<00:34, 1.03it/s]
66%|██████▌ | 66/100 [16:16<00:28, 1.18it/s]
67%|██████▋ | 67/100 [16:18<00:30, 1.09it/s]
68%|██████▊ | 68/100 [16:18<00:24, 1.28it/s]
69%|██████▉ | 69/100 [16:19<00:27, 1.11it/s]
70%|███████ | 70/100 [16:20<00:23, 1.28it/s]
71%|███████ | 71/100 [16:22<00:39, 1.37s/it]
72%|███████▏ | 72/100 [16:23<00:29, 1.05s/it]
73%|███████▎ | 73/100 [16:24<00:30, 1.13s/it]
74%|███████▍ | 74/100 [16:24<00:23, 1.12it/s]
75%|███████▌ | 75/100 [16:26<00:25, 1.03s/it]
76%|███████▌ | 76/100 [16:26<00:19, 1.25it/s]
77%|███████▋ | 77/100 [16:27<00:22, 1.02it/s]
78%|███████▊ | 78/100 [16:27<00:15, 1.39it/s]
79%|███████▉ | 79/100 [16:29<00:20, 1.05it/s]
80%|████████ | 80/100 [16:29<00:15, 1.26it/s]
81%|████████ | 81/100 [16:32<00:25, 1.36s/it]
82%|████████▏ | 82/100 [16:32<00:18, 1.05s/it]
83%|████████▎ | 83/100 [16:34<00:19, 1.14s/it]
84%|████████▍ | 84/100 [16:34<00:14, 1.12it/s]
85%|████████▌ | 85/100 [16:35<00:15, 1.01s/it]
86%|████████▌ | 86/100 [16:36<00:11, 1.21it/s]
87%|████████▋ | 87/100 [16:37<00:12, 1.03it/s]
88%|████████▊ | 88/100 [16:37<00:09, 1.31it/s]
89%|████████▉ | 89/100 [16:39<00:10, 1.05it/s]
90%|█████████ | 90/100 [16:39<00:07, 1.26it/s]
91%|█████████ | 91/100 [16:42<00:11, 1.33s/it]
92%|█████████▏| 92/100 [16:42<00:08, 1.07s/it]
93%|█████████▎| 93/100 [16:43<00:07, 1.10s/it]
94%|█████████▍| 94/100 [16:44<00:05, 1.12it/s]
95%|█████████▌| 95/100 [16:45<00:04, 1.01it/s]
96%|█████████▌| 96/100 [16:46<00:03, 1.19it/s]
97%|█████████▋| 97/100 [16:47<00:02, 1.08it/s]
98%|█████████▊| 98/100 [16:47<00:01, 1.28it/s]
99%|█████████▉| 99/100 [16:48<00:00, 1.09it/s]
100%|██████████| 100/100 [16:49<00:00, 1.31it/s]
100%|██████████| 100/100 [16:49<00:00, 10.09s/it]
File b02estimation~00.html has been generated.
File b02estimation~00.yaml has been generated.
Get the results in a pandas table
pandas_results = get_pandas_estimated_parameters(estimation_results=results)
display(pandas_results)
Name Value ... Bootstrap t-stat. Bootstrap p-value
0 beta_time_fulltime -1.596467 ... -5.445201 5.174692e-08
1 beta_time_other -0.552686 ... -2.037927 4.155725e-02
2 beta_cost -0.719072 ... -5.167539 2.371969e-07
3 mu_no_car 1.523716 ... 5.342280 9.178487e-08
4 asc_sm 0.063000 ... 0.284539 7.759970e-01
5 beta_dist_male -0.687677 ... -3.925931 8.639485e-05
6 beta_dist_female -0.832153 ... -3.929043 8.528455e-05
7 beta_dist_unreported -0.704737 ... -2.987710 2.810762e-03
8 asc_car 0.258607 ... 2.881082 3.963129e-03
[9 rows x 5 columns]
Simulation
simulated_choices = get_value_c(
expression=log_probability,
betas=results.get_beta_values(),
database=database,
numerically_safe=False,
use_jit=True,
)
display(simulated_choices)
Bootstraps: 0%| | 0/100 [16:49<?, ?it/s]
[-0.65553681 -1.42271161 -0.13348367 ... -0.37368096 -0.29814134
-0.27030377]
loglikelihood = get_value_c(
expression=log_probability,
betas=results.get_beta_values(),
database=database,
aggregation=True,
numerically_safe=False,
use_jit=True,
)
print(f'Final log likelihood: {results.final_log_likelihood}')
print(f'Simulated log likelihood: {loglikelihood}')
Final log likelihood: -1295.1202531057634
Simulated log likelihood: -1295.1202531057634
Total running time of the script: (16 minutes 53.085 seconds)