Note
Go to the end to download the full example code.
Calculation of revenues¶
We use an estimated model to calculate revenues.
Michel Bierlaire, EPFL Sat Jun 28 2025, 18:57:49
import sys
from pathlib import Path
import numpy as np
from biogeme.biogeme import BIOGEME
from biogeme.models import nested
from biogeme.results_processing import EstimationResults
try:
import matplotlib.pyplot as plt
can_plot = True
except ModuleNotFoundError:
can_plot = False
from scenarios import scenario
from biogeme.data.optima import normalized_weight, read_data
Read the estimation results from the file.
try:
results = EstimationResults.from_yaml_file(
filename=Path('saved_results') / 'b02estimation.yaml'
)
except FileNotFoundError:
sys.exit(
'Run first the script plot_b02estimation.py '
'in order to generate the '
'file b02estimation.yaml.'
)
Read the data
database = read_data()
Function calculating the revenues
def revenues(factor: float) -> tuple[float, float, float]:
"""Calculate the total revenues generated by public transportation,
when the price is multiplied by a factor.
:param factor: factor that multiplies the current cost of public
transportation
:return: total revenues, followed by the lower and upper bound of
the confidence interval.
"""
filename = Path(f'revenue_{factor:.2f}.txt')
SEPARATOR = '%'
try:
with open(filename, 'r') as f:
line = f.read()
revenue, left, right = line.split(SEPARATOR)
return float(revenue), float(left), float(right)
except FileNotFoundError:
...
# Obtain the specification for the default scenario
utilities, nests, _, marginal_cost_scenario = scenario(factor=factor)
# Obtain the expression for the choice probability of each alternative
prob_pt = nested(utilities, None, nests, 0)
# We now simulate the choice probabilities,the weight and the
# price variable
simulate = {
'weight': normalized_weight,
'Revenue public transportation': prob_pt * marginal_cost_scenario,
}
the_biogeme = BIOGEME(database, simulate)
simulated_values = the_biogeme.simulate(results.get_beta_values())
# We also calculate confidence intervals for the calculated quantities
beta_bootstrap = results.get_betas_for_sensitivity_analysis()
left, right = the_biogeme.confidence_intervals(beta_bootstrap, 0.9)
revenues_pt = (
simulated_values['Revenue public transportation'] * simulated_values['weight']
).sum()
revenues_pt_left = (left['Revenue public transportation'] * left['weight']).sum()
revenues_pt_right = (right['Revenue public transportation'] * right['weight']).sum()
with open(filename, 'w') as f:
print(
f'{revenues_pt} {SEPARATOR} {revenues_pt_left} {SEPARATOR} {revenues_pt_right}',
file=f,
)
return revenues_pt, revenues_pt_left, revenues_pt_right
Current revenues for public transportation
r, r_left, r_right = revenues(factor=1.0)
print(
f'Total revenues for public transportation (for the sample): {r:.1f} CHF '
f'[{r_left:.1f} CHF, '
f'{r_right:.1f} CHF]'
)
Total revenues for public transportation (for the sample): 3043.3 CHF [2445.6 CHF, 3676.4 CHF]
We now investigate how the revenues vary with the multiplicative factor
factors = np.arange(0.0, 5.0, 0.1)
plot_revenues = [revenues(s) for s in factors]
zipped = zip(*plot_revenues, strict=True)
rev = next(zipped)
lower = next(zipped)
upper = next(zipped)
largest_revenue = max(rev)
max_index = rev.index(largest_revenue)
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 3177.12it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2906.37it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 3137.76it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2616.78it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2541.11it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2882.49it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2947.16it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2869.27it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2694.91it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2807.94it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2703.14it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2842.77it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2899.56it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 3040.34it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2821.52it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2805.35it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2809.54it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2899.86it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2703.80it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2970.26it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2812.70it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 1700.05it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2702.18it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2963.61it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2749.43it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2938.24it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2955.19it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2793.22it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2023.88it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2943.78it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2828.71it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2899.96it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2715.02it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 1741.73it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2880.94it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2926.98it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2956.48it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2853.89it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2778.24it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2320.99it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2542.06it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2719.32it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2934.62it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2863.18it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2930.58it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2289.17it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 1904.65it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2746.87it/s]
0%| | 0/100 [00:00<?, ?it/s]
100%|██████████| 100/100 [00:00<00:00, 2811.44it/s]
print(
f'Largest revenue: {largest_revenue:.1f} obtained with '
f'factor {factors[max_index]:.1f}'
)
Largest revenue: 3062.4 obtained with factor 1.2
if can_plot:
# We plot the results
ax = plt.gca()
ax.plot(factors, rev, label="Revenues")
ax.plot(factors, lower, label="Lower bound of the CI")
ax.plot(factors, upper, label="Upper bound of the CI")
ax.legend()
plt.show()
Total running time of the script: (0 minutes 9.506 seconds)