""" File: 09-lecture_solution.py Michel Bierlaire Thu Aug 07 2025, 08:47:55 """ import biogeme.biogeme_logging as blog from IPython.core.display_functions import display from biogeme.biogeme import BIOGEME from biogeme.expressions import Beta, Draws, MonteCarlo, exp, log from biogeme.models import logit, loglogit from biogeme.results_processing import ( EstimationResults, compile_estimation_results, get_pandas_estimated_parameters, ) # Variables used for the specification of the Swissmetro model are defined in the file `swissmetro_variables.py`. from swissmetro_variables import ( BUSINESS, CAR_AV_SP, CAR_CO_SCALED, CAR_TT_SCALED, CHOICE, FIRST, GA, LOW_INC, MALE, SM_AV, SM_COST_SCALED, SM_HE_SCALED, SM_TT_SCALED, TRAIN_AV_SP, TRAIN_COST_SCALED, TRAIN_HE_SCALED, TRAIN_TT_SCALED, database, ) # As the estimation time may be long, we ask Biogeme to report the details of the iterations. logger = blog.get_screen_logger(level=blog.INFO) # **Tip:**
It is advised to start working with a low number of draws, # until the script is working well. Then, increase the number of draws to 10000, say. Then, execute the script # overnight.
number_of_draws = 10 # # Parameters 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) b_fr = Beta('b_fr', 0, None, None, 0) # # Availability conditions av = {1: TRAIN_AV_SP, 2: SM_AV, 3: CAR_AV_SP} # # Logit model # ## Utility functions v_train = ( asc_train + b_time * TRAIN_TT_SCALED + b_cost * TRAIN_COST_SCALED + b_fr * TRAIN_HE_SCALED ) v_swissmetro = b_time * SM_TT_SCALED + b_cost * SM_COST_SCALED + b_fr * SM_HE_SCALED v_car = asc_car + b_time * CAR_TT_SCALED + b_cost * CAR_CO_SCALED v = {1: v_train, 2: v_swissmetro, 3: v_car} # ## Model logprob = loglogit(v, av, CHOICE) # ## Estimation biogeme = BIOGEME(database, logprob) biogeme.model_name = '01logit' results_logit: EstimationResults = biogeme.estimate(recycle=True) # ## Results # General statistics print(results_logit.short_summary()) # Estimated parameters param_logit = get_pandas_estimated_parameters(estimation_results=results_logit) display(param_logit) # # Random parameter: normal distribution # Read the results from file results_normal = EstimationResults.from_yaml_file(filename='02normal.yaml') # # Random parameter: lognormal distribution # Read the results from file results_lognormal = EstimationResults.from_yaml_file(filename='03lognormal.yaml') # # Latent classes # Read the results from file results_latent = EstimationResults.from_yaml_file(filename='04latentClass.yaml') # # Latent classes with class membership model # Read the results from file results_latentsocio = EstimationResults.from_yaml_file(filename='05latentClass.yaml') # # Latent classes with random parameter # We consider again two classes in the population. The first class of individuals have considered all variables # when making their choice. For them, the specification of the utility function is the same as for the logit model, # where the time coefficient is now distributed in the population. b_time_s = Beta('b_time_s', 1, None, None, 0) b_time_rnd = b_time + b_time_s * Draws('b_time_rnd', 'NORMAL') # Utility function, class 1 v_train_class_1 = ( asc_train + b_time_rnd * TRAIN_TT_SCALED + b_cost * TRAIN_COST_SCALED + b_fr * TRAIN_HE_SCALED ) v_swissmetro_class_1 = ( b_time * SM_TT_SCALED + b_cost * SM_COST_SCALED + b_fr * SM_HE_SCALED ) v_car_class_1 = asc_car + b_time * CAR_TT_SCALED + b_cost * CAR_CO_SCALED v_class_1 = {1: v_train_class_1, 2: v_swissmetro_class_1, 3: v_car_class_1} # The second class of individuals ignored the travel time variable when making the choice. Therefore, this variable # is removed from the utility function. v_train_class_2 = asc_train + b_cost * TRAIN_COST_SCALED + b_fr * TRAIN_HE_SCALED v_swissmetro_class_2 = b_cost * SM_COST_SCALED + b_fr * SM_HE_SCALED v_car_class_2 = asc_car + b_cost * CAR_CO_SCALED v_class_2 = {1: v_train_class_2, 2: v_swissmetro_class_2, 3: v_car_class_2} # The following parameters are involved in the class membership model. g_intercept = Beta('g_intercept', 0, None, None, 0) g_male = Beta('g_male', 0, None, None, 0) g_ga = Beta('g_ga', 0, None, None, 0) g_business = Beta('g_business', 0, None, None, 0) g_low_inc = Beta('g_low_inc', 0, None, None, 0) g_first = Beta('g_first', 0, None, None, 0) # Class membership model. Note that `omega` can potentially take any real value. We have to transform it into a # probability using the transform `1 / (1 + exp(omega))` omega = ( g_intercept + g_male * MALE + g_ga * GA + g_business * BUSINESS + g_low_inc * LOW_INC + g_first * FIRST ) prob_class_1 = 1 / (1 + exp(omega)) prob_class_2 = 1 - prob_class_1 # ## Model # We first calculate the choice probability for each class. cond_choice_prob_class_1 = logit(v_class_1, av, CHOICE) choice_prob_class_1 = MonteCarlo(cond_choice_prob_class_1) choice_prob_class_2 = logit(v_class_2, av, CHOICE) # The choice probability is obtained by using the class membership model. prob = prob_class_1 * choice_prob_class_1 + prob_class_2 * choice_prob_class_2 logprob = log(prob) # ## Estimation biogeme = BIOGEME(database, logprob, number_of_draws=number_of_draws) biogeme.model_name = '06mixedLatentClass' results_latentrandom: EstimationResults = biogeme.estimate() # ## Results # General statistics print(results_latentrandom.short_summary()) # Estimated parameters param_latentrandom = get_pandas_estimated_parameters( estimation_results=results_latentrandom ) display(param_latentrandom) # # Comparison # We build a summary data frame. We first gather the parameter estimates for each model. summary = compile_estimation_results( { 'Logit': results_logit, 'Random param. (normal)': results_normal, 'Random param. (lognormal)': results_lognormal, 'Latent class': results_latent, 'Latent class with class mbship': results_latentsocio, 'Latent with random param.': results_latentrandom, } ) display(summary[0])