Gamma-profile MDCEV forecasting

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: '64.8', 2: '201', 3: '188', 4: '46.1'} objective 56.9, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '63.7', 2: '200', 3: '189', 4: '46.7'} objective 56.9, constraint 500, choice set {1, 2, 3, 4}
Difference between optimal utility with analytical [56.88] and brute force [56.88] algorithms.
Solution with brute force: 1: 64.8, 2: 201, 3: 188, 4: 46.1
Solution with analytical: 1: 63.7, 2: 200, 3: 189, 4: 46.7
============ Comparison ===================
Brute force: {1: '17.4', 2: '76.9', 3: '389', 4: '17.1'} objective 128, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '17.4', 2: '76.9', 3: '389', 4: '17.1'} objective 128, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '77.9', 2: '191', 3: '176', 4: '55.8'} objective 53.9, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '76.8', 2: '190', 3: '177', 4: '56.4'} objective 53.9, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '52.2', 2: '157', 3: '257', 4: '33.6'} objective 57.3, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '52.7', 2: '155', 3: '258', 4: '33.5'} objective 57.3, constraint 500, choice set {1, 2, 3, 4}
Difference between optimal utility with analytical [57.26] and brute force [57.26] algorithms.
Solution with brute force: 1: 52.2, 2: 157, 3: 257, 4: 33.6
Solution with analytical: 1: 52.7, 2: 155, 3: 258, 4: 33.5
============ Comparison ===================
Brute force: {1: '81', 2: '210', 3: '174', 4: '35.7'} objective 54, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '81', 2: '209', 3: '174', 4: '35.6'} objective 54, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '58.5', 2: '196', 3: '206', 4: '40.1'} objective 61.9, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '58.6', 2: '195', 3: '207', 4: '40.1'} objective 61.9, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '33', 2: '279', 3: '157', 4: '31.4'} objective 80.2, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '33.1', 2: '280', 3: '156', 4: '31.4'} objective 80.2, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '27.5', 2: '363', 3: '90', 4: '19.7'} objective 108, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '27.5', 2: '363', 3: '90', 4: '19.7'} objective 108, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '36.8', 2: '280', 3: '137', 4: '45.8'} objective 69, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '37.2', 2: '281', 3: '136', 4: '45.6'} objective 69, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '29.1', 2: '326', 3: '113', 4: '32.3'} objective 85.2, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '29.1', 2: '326', 3: '113', 4: '32.3'} objective 85.2, constraint 500, choice set {1, 2, 3, 4}
Forecasting observation 1 / 2 [10 draws]
============ Comparison ===================
Brute force: {1: '40', 2: '218', 3: '200', 4: '41.5'} objective 69.1, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '39.4', 2: '217', 3: '201', 4: '42.1'} objective 69.1, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '34', 2: '217', 3: '231', 4: '18.5'} objective 74.6, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '34.1', 2: '215', 3: '232', 4: '18.5'} objective 74.6, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '27', 2: '252', 3: '193', 4: '27.9'} objective 102, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '27', 2: '252', 3: '193', 4: '27.9'} objective 102, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '37.8', 2: '149', 3: '287', 4: '26.1'} objective 74.1, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '37.8', 2: '149', 3: '287', 4: '26.1'} objective 74.1, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '40.6', 2: '184', 3: '245', 4: '30.3'} objective 66, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '40.9', 2: '182', 3: '247', 4: '30.2'} objective 66, constraint 500, choice set {1, 2, 3, 4}
Difference between optimal utility with analytical [66.02] and brute force [66.02] algorithms.
Solution with brute force: 1: 40.6, 2: 184, 3: 245, 4: 30.3
Solution with analytical: 1: 40.9, 2: 182, 3: 247, 4: 30.2
============ Comparison ===================
Brute force: {1: '46.4', 2: '199', 3: '224', 4: '31.3'} objective 73.7, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '46.5', 2: '197', 3: '225', 4: '31.2'} objective 73.7, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '49.5', 2: '133', 3: '251', 4: '66.8'} objective 66, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '49.5', 2: '133', 3: '251', 4: '66.7'} objective 66, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '52.2', 2: '154', 3: '265', 4: '28.5'} objective 77.8, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '52.2', 2: '154', 3: '265', 4: '28.5'} objective 77.8, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '86.4', 2: '207', 3: '179', 4: '28.5'} objective 74.4, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '86.4', 2: '206', 3: '179', 4: '28.5'} objective 74.4, constraint 500, choice set {1, 2, 3, 4}
============ Comparison ===================
Brute force: {1: '41.7', 2: '269', 3: '162', 4: '27.3'} objective 74.7, constraint 500, choice set {1, 2, 3, 4}
Analytical:  {1: '41.8', 2: '270', 3: '161', 4: '27.3'} objective 74.7, 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: 941 seconds
Forecasting observation 0 / 2 [2000 draws]
Forecasting observation 1 / 2 [2000 draws]
Execution time for 2000 draws with analytical algorithm: 2.2 seconds
                 1            2            3            4
count  2000.000000  2000.000000  2000.000000  2000.000000
mean     50.689629   209.821465   188.804548    50.684358
std      21.672126    51.900031    50.710181    20.575598
min       3.063425    52.138017    11.358548     3.676516
25%      36.456433   175.633250   155.896950    37.475468
50%      46.918936   206.143995   183.228205    46.717334
75%      60.061013   239.497875   215.591551    59.766440
max     216.995249   481.901511   412.552424   226.084040
                 1            2            3            4
count  2000.000000  2000.000000  2000.000000  2000.000000
mean     50.539994   209.506294   189.119551    50.834161
std      21.481870    52.417949    51.240780    20.488075
min       3.063248    52.140852    11.353412     3.676602
25%      36.444822   174.238863   155.166167    37.482355
50%      46.784260   205.877118   183.852707    46.947076
75%      59.701625   239.923761   216.911875    60.153328
max     216.655047   481.906738   412.539641   226.100450
                 1            2            3            4
count  2000.000000  2000.000000  2000.000000  2000.000000
mean     53.366985   195.731492   213.482213    37.419310
std      22.041282    50.539478    51.053308    16.538686
min      10.304699    40.051978    45.006069     7.506174
25%      38.470250   161.166883   178.873448    26.688697
50%      49.210175   191.502999   209.735488    34.168739
75%      63.507554   224.600702   244.575753    44.186392
max     243.093648   430.297764   441.563183   179.465845
                 1            2            3            4
count  2000.000000  2000.000000  2000.000000  2000.000000
mean     53.393998   195.163814   213.991130    37.451058
std      21.978001    50.863387    51.276689    16.554229
min      10.305108    40.055129    45.006866     7.507159
25%      38.529667   159.873296   178.981012    26.671153
50%      49.263827   190.484703   210.787014    34.115880
75%      63.640671   224.303647   246.090367    44.246687
max     242.979945   430.296348   441.567172   178.050901

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