import pandas as pd import numpy as np import biogeme.database as db import biogeme.biogeme as bio import biogeme.models as models import biogeme.distributions as dist pandas = pd.read_table("optima.dat") database = db.Database("optima",pandas) from headers import * exclude = (Choice == -1.0) database.remove(exclude) ### Variables ScaledIncome = DefineVariable('ScaledIncome',\ CalculatedIncome / 1000,database) ContIncome_0_4000 = DefineVariable('ContIncome_0_4000',\ bioMin(ScaledIncome,4),database) ContIncome_4000_6000 = DefineVariable('ContIncome_4000_6000',\ bioMax(0,bioMin(ScaledIncome-4,2)),database) ContIncome_6000_8000 = DefineVariable('ContIncome_6000_8000',\ bioMax(0,bioMin(ScaledIncome-6,2)),database) ContIncome_8000_10000 = DefineVariable('ContIncome_8000_10000',\ bioMax(0,bioMin(ScaledIncome-8,2)),database) ContIncome_10000_more = DefineVariable('ContIncome_10000_more',\ bioMax(0,ScaledIncome-10),database) age_65_more = DefineVariable('age_65_more',age >= Numeric(65),database) moreThanOneCar = DefineVariable('moreThanOneCar',NbCar > 1,database) moreThanOneBike = DefineVariable('moreThanOneBike',NbBicy > 1,database) individualHouse = DefineVariable('individualHouse',\ HouseType == 1,database) male = DefineVariable('male',Gender == 1,database) haveChildren = DefineVariable('haveChildren',\ ((FamilSitu == 3)+(FamilSitu == 4)) > 0,database) haveGA = DefineVariable('haveGA',GenAbST == 1,database) highEducation = DefineVariable('highEducation', Education >= 6,database) ### Coefficients coef_intercept = Beta('coef_intercept',0.0,None,None,1) coef_age_65_more = Beta('coef_age_65_more',0.0,None,None,0) coef_haveGA = Beta('coef_haveGA',0.0,None,None,0) coef_ContIncome_0_4000 = \ Beta('coef_ContIncome_0_4000',0.0,None,None,0) coef_ContIncome_4000_6000 = \ Beta('coef_ContIncome_4000_6000',0.0,None,None,0) coef_ContIncome_6000_8000 = \ Beta('coef_ContIncome_6000_8000',0.0,None,None,0) coef_ContIncome_8000_10000 = \ Beta('coef_ContIncome_8000_10000',0.0,None,None,0) coef_ContIncome_10000_more = \ Beta('coef_ContIncome_10000_more',0.0,None,None,0) coef_moreThanOneCar = \ Beta('coef_moreThanOneCar',0.0,None,None,0) coef_moreThanOneBike = \ Beta('coef_moreThanOneBike',0.0,None,None,0) coef_individualHouse = \ Beta('coef_individualHouse',0.0,None,None,0) coef_male = Beta('coef_male',0.0,None,None,0) coef_haveChildren = Beta('coef_haveChildren',0.0,None,None,0) coef_highEducation = Beta('coef_highEducation',0.0,None,None,0) ### Latent variable: structural equation # Note that the expression must be on a single line. In order to # write it across several lines, each line must terminate with # the \ symbol omega = RandomVariable('omega') density = dist.normalpdf(omega) sigma_s = Beta('sigma_s',1,None,None,1) CARLOVERS = \ coef_intercept +\ coef_age_65_more * age_65_more +\ coef_ContIncome_0_4000 * ContIncome_0_4000 +\ coef_ContIncome_4000_6000 * ContIncome_4000_6000 +\ coef_ContIncome_6000_8000 * ContIncome_6000_8000 +\ coef_ContIncome_8000_10000 * ContIncome_8000_10000 +\ coef_ContIncome_10000_more * ContIncome_10000_more +\ coef_moreThanOneCar * moreThanOneCar +\ coef_moreThanOneBike * moreThanOneBike +\ coef_individualHouse * individualHouse +\ coef_male * male +\ coef_haveChildren * haveChildren +\ coef_haveGA * haveGA +\ coef_highEducation * highEducation +\ sigma_s * omega # Choice model ASC_CAR = Beta('ASC_CAR',0.0,None,None,0) ASC_PT = Beta('ASC_PT',0.0,None,None,1) ASC_SM = Beta('ASC_SM',0.0,None,None,0) BETA_COST_HWH = Beta('BETA_COST_HWH',0.0,None,None,0) BETA_COST_OTHER = Beta('BETA_COST_OTHER',0.0,None,None,0) BETA_DIST = Beta('BETA_DIST',0.0,None,None,0) BETA_TIME_CAR_REF = Beta('BETA_TIME_CAR_REF',0.0,None,0,0) BETA_TIME_CAR_CL = Beta('BETA_TIME_CAR_CL',0.0,None,None,0) BETA_TIME_PT_REF = Beta('BETA_TIME_PT_REF',0.0,None,0,0) BETA_TIME_PT_CL = Beta('BETA_TIME_PT_CL',-1.0,None,None,0) BETA_WAITING_TIME = Beta('BETA_WAITING_TIME',0.0,None,None,0) TimePT_scaled = DefineVariable('TimePT_scaled', TimePT / 200 ,database) TimeCar_scaled = DefineVariable('TimeCar_scaled', TimeCar / 200 ,database) MarginalCostPT_scaled = \ DefineVariable('MarginalCostPT_scaled', MarginalCostPT / 10 ,database) CostCarCHF_scaled = \ DefineVariable('CostCarCHF_scaled', CostCarCHF / 10 ,database) distance_km_scaled = \ DefineVariable('distance_km_scaled', distance_km / 5 ,database) PurpHWH = DefineVariable('PurpHWH', TripPurpose == 1,database) PurpOther = DefineVariable('PurpOther', TripPurpose != 1,database) ### DEFINITION OF UTILITY FUNCTIONS: BETA_TIME_PT = BETA_TIME_PT_REF * \ exp(BETA_TIME_PT_CL * CARLOVERS) V0 = ASC_PT + \ BETA_TIME_PT * TimePT_scaled + \ BETA_WAITING_TIME * WaitingTimePT + \ BETA_COST_HWH * MarginalCostPT_scaled * PurpHWH +\ BETA_COST_OTHER * MarginalCostPT_scaled * PurpOther BETA_TIME_CAR = BETA_TIME_CAR_REF * \ exp(BETA_TIME_CAR_CL * CARLOVERS) V1 = ASC_CAR + \ BETA_TIME_CAR * TimeCar_scaled + \ BETA_COST_HWH * CostCarCHF_scaled * PurpHWH + \ BETA_COST_OTHER * CostCarCHF_scaled * PurpOther V2 = ASC_SM + BETA_DIST * distance_km_scaled # Associate utility functions with the numbering of alternatives V = {0: V0, 1: V1, 2: V2} # Associate the availability conditions with the alternatives. # In this example all alternatives are available # for each individual. av = {0: 1, 1: 1, 2: 1} # The choice model is a logit, conditional to # the value of the latent variable condprob = models.logit(V,av,Choice) prob = Integrate(condprob * density,'omega') loglike = log(prob) biogeme = bio.BIOGEME(database,loglike) biogeme.modelName = "03choiceOnly" results = biogeme.estimate() print(f"Estimated betas: {len(results.data.betaValues)}") print(f"Final log likelihood: {results.data.logLike:.3f}") print(f"Output file: {results.data.htmlFileName}") results.writeLaTeX() print(f"LaTeX file: {results.data.latexFileName}")