File gamma_forecasting.py

Michel Bierlaire, EPFL Fri Jul 25 2025, 16:38:12 Forecasting with a MDCEV model and the “gamma_profile” specification.

Example: gamma profile utility
Forecasting observation 0 / 2 [10 draws]
============ Comparison ===================
Brute force: {1: '19.3', 2: '191', 3: '184', 4: '106'} objective 75.9, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '19.5', 2: '190', 3: '185', 4: '105'} objective 75.9, constraint 500, choice set {1, 2, 3, 4}
Difference between optimal utility with analytical [75.9] and brute force [75.89] algorithms.
Solution with brute force: 1: 19.3, 2: 191, 3: 184, 4: 106
Solution with analytical: 1: 19.5, 2: 190, 3: 185, 4: 105
============ Comparison ===================
Brute force: {1: '96.5', 2: '184', 3: '148', 4: '71.4'} objective 48.1, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '95.6', 2: '184', 3: '148', 4: '71.9'} objective 48.1, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '98.8', 2: '224', 3: '128', 4: '48.6'} objective 55.2, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '97.5', 2: '227', 3: '126', 4: '49.4'} objective 55.2, constraint 500, choice set {1, 2, 3, 4}
Difference between optimal utility with analytical [55.22] and brute force [55.22] algorithms.
Solution with brute force: 1: 98.8, 2: 224, 3: 128, 4: 48.6
Solution with analytical: 1: 97.5, 2: 227, 3: 126, 4: 49.4
============ Comparison ===================
Brute force: {1: '36.5', 2: '160', 3: '271', 4: '32.7'} objective 70.3, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '36.5', 2: '160', 3: '271', 4: '32.7'} objective 70.3, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '69.3', 2: '227', 3: '169', 4: '35'} objective 65, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '69.3', 2: '227', 3: '169', 4: '35'} objective 65, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '43.4', 2: '262', 3: '150', 4: '44.3'} objective 54.9, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '43', 2: '264', 3: '148', 4: '45'} objective 54.9, constraint 500, choice set {1, 2, 3, 4}
Difference between optimal utility with analytical [54.86] and brute force [54.86] algorithms.
Solution with brute force: 1: 43.4, 2: 262, 3: 150, 4: 44.3
Solution with analytical: 1: 43, 2: 264, 3: 148, 4: 45
============ Comparison ===================
Brute force: {1: '38.6', 2: '190', 3: '208', 4: '63.7'} objective 52.2, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '38.6', 2: '189', 3: '209', 4: '63.7'} objective 52.2, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '49.2', 2: '247', 3: '143', 4: '61.2'} objective 47.2, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '49.1', 2: '248', 3: '141', 4: '61.6'} objective 47.2, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '42.5', 2: '197', 3: '197', 4: '63.6'} objective 56.9, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '42.4', 2: '195', 3: '198', 4: '63.7'} objective 56.9, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '48.4', 2: '256', 3: '163', 4: '32.7'} objective 68.2, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '48.5', 2: '257', 3: '162', 4: '32.7'} objective 68.2, constraint 500, choice set {1, 2, 3, 4}
Forecasting observation 1 / 2 [10 draws]
============ Comparison ===================
Brute force: {1: '41.4', 2: '196', 3: '226', 4: '37.2'} objective 66.5, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '41.6', 2: '194', 3: '227', 4: '37.1'} objective 66.5, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '58.4', 2: '138', 3: '260', 4: '44'} objective 63, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '59.2', 2: '136', 3: '261', 4: '43.8'} objective 63, constraint 500, choice set {1, 2, 3, 4}
Difference between optimal utility with analytical [62.98] and brute force [62.98] algorithms.
Solution with brute force: 1: 58.4, 2: 138, 3: 260, 4: 44
Solution with analytical: 1: 59.2, 2: 136, 3: 261, 4: 43.8
============ Comparison ===================
Brute force: {1: '51.4', 2: '189', 3: '204', 4: '55.8'} objective 75.6, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '51', 2: '187', 3: '206', 4: '56.3'} objective 75.6, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '33.5', 2: '233', 3: '213', 4: '21'} objective 90, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '33.5', 2: '232', 3: '214', 4: '20.9'} objective 90, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '46', 2: '237', 3: '181', 4: '36'} objective 59.4, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '46', 2: '237', 3: '181', 4: '36'} objective 59.4, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '53.9', 2: '249', 3: '168', 4: '29.7'} objective 79.9, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '54', 2: '249', 3: '167', 4: '29.7'} objective 79.9, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '55.7', 2: '129', 3: '281', 4: '33.9'} objective 65, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '55.7', 2: '129', 3: '281', 4: '33.9'} objective 65, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '75.9', 2: '182', 3: '180', 4: '61.8'} objective 59.4, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '75.1', 2: '182', 3: '181', 4: '62.4'} objective 59.4, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '54.6', 2: '204', 3: '192', 4: '49.1'} objective 103, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '54.6', 2: '203', 3: '193', 4: '49'} objective 103, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '63.6', 2: '224', 3: '187', 4: '24.7'} objective 77.1, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '63.7', 2: '224', 3: '188', 4: '24.7'} objective 77.1, constraint 500, choice set {1, 2, 3, 4}
Forecasting observation 0 / 2 [2000 draws]
Forecasting observation 1 / 2 [2000 draws]
Execution time for 2000 draws with brute force algorithm: 12.3 seconds
Forecasting observation 0 / 2 [2000 draws]
Forecasting observation 1 / 2 [2000 draws]
Execution time for 2000 draws with analytical algorithm: 0.608 seconds
                 1            2            3            4
count  2000.000000  2000.000000  2000.000000  2000.000000
mean     51.518231   211.054094   187.229841    50.197834
std      23.130843    53.947558    53.337971    21.550578
min       8.728227    43.901284    44.890577     8.076306
25%      36.143132   175.047718   153.027933    36.446719
50%      47.029913   206.956391   182.614542    46.386606
75%      60.692928   242.079979   215.805571    59.587119
max     234.492435   408.033021   419.227802   366.715204
                 1            2            3            4
count  2000.000000  2000.000000  2000.000000  2000.000000
mean     51.379089   210.774331   187.507669    50.338911
std      22.941211    54.448191    53.867760    21.476516
min       8.728471    43.905352    44.893016     8.076337
25%      36.194367   173.900701   152.041808    36.585668
50%      47.048983   206.302898   183.120942    46.715967
75%      60.471824   242.938136   217.143337    59.804237
max     234.449458   408.017802   419.227170   366.717718
                 1            2            3            4
count  2000.000000  2000.000000  2000.000000  2000.000000
mean     53.343868   195.614404   213.815285    37.226442
std      22.464387    50.835646    52.409488    16.870743
min      10.628288    47.110956    52.865721     9.782055
25%      38.782480   161.272774   177.984059    26.807845
50%      48.828279   191.583525   209.828357    33.587426
75%      62.290773   225.470817   245.341326    43.262908
max     247.160308   408.160130   426.049392   177.230517
                 1            2            3            4
count  2000.000000  2000.000000  2000.000000  2000.000000
mean     53.363183   195.091525   214.290462    37.254830
std      22.388807    51.122401    52.624925    16.869329
min      10.626606    47.111895    52.871675     9.781980
25%      38.782628   159.984587   178.026109    26.777725
50%      48.899298   190.640290   211.003794    33.561971
75%      62.348376   225.067649   246.438383    43.377101
max     247.026101   408.156280   426.049082   176.091517

import sys
import time

import numpy as np
import pandas as pd
from IPython.core.display_functions import display

import biogeme.biogeme_logging as blog
from biogeme.database import Database
from biogeme.results_processing import EstimationResults
from gamma_specification import the_gamma_profile
from process_data import database

logger = blog.get_screen_logger(level=blog.INFO)
logger.info('Example: gamma profile utility')

result_file = 'saved_results/gamma_profile.yaml'
try:
    results = EstimationResults.from_yaml_file(filename=result_file)
except FileNotFoundError as e:
    print(e)
    print(f'File {result_file} is missing.')
    sys.exit()

the_gamma_profile.estimation_results = results

# %
# We apply the model only on the first two rows of the database.
two_rows_of_database: Database = database.extract_rows([0, 1])

# %
budget_in_hours = 500

# %
# # Validation

# %
# As the implementation is still experimental, we compare the result obtained by the bruteforce algorithm and
# the analytical algorithm for a few draws.

# Note that minor discrepancies between the outcome of the two algorithms are likely to occur, due to numerical
# imprecision, inevitable in finite arithmetic.

# However, if there are major differences, it should be reported.

# %
number_of_draws = 10

# %
# We generate the draws
epsilons = [
    np.random.gumbel(
        loc=0, scale=1, size=(number_of_draws, the_gamma_profile.number_of_alternatives)
    )
    for _ in range(two_rows_of_database.num_rows())
]

# %
# We first compare the results obtained from the brute force and the analytical algorithms, for each draw.
the_gamma_profile.validate_forecast(
    database=two_rows_of_database, total_budget=budget_in_hours, epsilons=epsilons
)

# %
# # Forecasting
# We use a larger number of draws to obtain the forecast.

# %
number_of_draws = 2000

# %
# We generate the draws
epsilons = [
    np.random.gumbel(
        loc=0, scale=1, size=(number_of_draws, the_gamma_profile.number_of_alternatives)
    )
    for _ in range(two_rows_of_database.num_rows())
]

# %
# First, the brute force algorithm.
start_time = time.time()
optimal_consumptions_brute_force: list[pd.DataFrame] = the_gamma_profile.forecast(
    database=two_rows_of_database,
    total_budget=budget_in_hours,
    epsilons=epsilons,
    brute_force=True,
)
end_time = time.time()

# %
print(
    f'Execution time for {number_of_draws} draws with brute force algorithm: {end_time-start_time:.3g} seconds'
)

# %
# Then, the analytical algorithm.
start_time = time.time()
optimal_consumptions_analytical: list[pd.DataFrame] = the_gamma_profile.forecast(
    database=two_rows_of_database,
    total_budget=budget_in_hours,
    epsilons=epsilons,
    brute_force=False,
)
end_time = time.time()

# %
print(
    f'Execution time for {number_of_draws} draws with analytical algorithm: {end_time-start_time:.3g} seconds'
)

# %
# Results for the first observation, brute force method
display(optimal_consumptions_brute_force[0].describe())

# %
# Results for the first observation, analytical method
display(optimal_consumptions_analytical[0].describe())

# %
# Results for the second observation, brute force method
display(optimal_consumptions_brute_force[1].describe())

# %
# Results for the second observation, analytical method
display(optimal_consumptions_analytical[1].describe())