Error while conducting a MonteCarlo integration estimation
Yuan Wang
Dear Professor Bierlaire and the Biogeme Community,
I hope you are all well.
I am reaching out for help since I encountered something wrong while using Biogeme to estimate. In my hybrid choice model, I have five latent variables, I know it is recommended to utilize a Monte-Carlo integration rather than a normal integration. I have done latent-variable-only estimation and choice-only estimation to get the initial values of the coefficients, and they all ran smoothly. However, while conducting the Monte-Carlo integration, something went wrong.
The error popped out while running this command:
biogeme = bio.BIOGEME(database, logprob, numberOfDraws = 500)
The error reads as follows:
biogeme.exceptions.BiogemeError: The following random variables are defined outside the Integrate operator: set()
The argument of MonteCarlo must contain a bioDraws: MonteCarlo(exp(_bioLogLogit[choice=helmetchoice]U=(0:(((((((((((((ASC_WEARINGHELMET(init=-0.38904884960128994) + (BETA_PENALTYPRICE(init=-0.029062031105895475) * penaltyprice)) + (BETA_PENALTYCHANCE(init=-0.4246962728594313) * penaltychance)) + (BETA_HELMETPERCENT(init=-0.0016743818058945331) * helmetpercent)) + (BETA_GENDER(init=-0.5007207157800359) * gender)) + (BETA_EDUCATION(init=-0.2001399820182707) * education)) + (BETA_INCOME(init=0.042334994995635505) * income)) + (BETA_MARRIAGE(init=1.018321338420314) * marriage)) + (BETA_AGE(init=0.007584758915581489) * age)) + (BETA_PBA(init=0) * ((((((intercept_2(init=0.0) + (coef_age_2(init=0.0) * age)) + (coef_gender_2(init=0.0) * gender)) + (coef_education_2(init=0.0) * education)) + (coef_income_2(init=0.0) * income)) + (coef_marriage_2(init=0.0) * marriage)) + (sigma_s2(init=1) * omega)))) + (BETA_PBE(init=0) * ((((((intercept_1(init=0.0) + (coef_age_1(init=0.0) * age)) + (coef_gender_1(init=0.0) * gender)) + (coef_education_1(init=0.0) * education)) + (coef_income_1(init=0.0) * income)) + (coef_marriage_1(init=0.0) * marriage)) + (sigma_s1(init=1) * omega)))) + (BETA_PSU(init=0) * ((((((intercept_3(init=0.057908251263525244) + (coef_age_3(init=0.008089544421022029) * age)) + (coef_gender_3(init=-0.3795903939693467) * gender)) + (coef_education_3(init=-0.0827081649727329) * education)) + (coef_income_3(init=-0.0827081649727329) * income)) + (coef_marriage_3(init=0.12407650963342834) * marriage)) + (sigma_s3(init=1) * omega)))) + (BETA_PSE(init=0) * ((((((intercept_4(init=0.02832275041016156) + (coef_age_4(init=0.0076098099971716035) * age)) + (coef_gender_4(init=-0.3001216474429452) * gender)) + (coef_education_4(init=-0.02470114184285887) * education)) + (coef_income_4(init=-0.02470114184285887) * income)) + (coef_marriage_4(init=0.2924756362319662) * marriage)) + (sigma_s4(init=1) * omega)))) + (BETA_NML(init=0) * ((((((intercept_5(init=-0.06570067252771641) + (coef_age_5(init=0.014186313627749162) * age)) + (coef_gender_5(init=-0.10546093216847484) * gender)) + (coef_education_5(init=0.007627976573005183) * education)) + (coef_income_5(init=0.007627976573005183) * income)) + (coef_marriage_5(init=-0.12201722237771571) * marriage)) + (sigma_s5(init=1) * omega)))), 1:ASC_NON_WEARING(fixed=0))av=(0:`1.0`, 1:`1.0`)))
The error seems horrifying because in my model there are many latent variables and each one contains 3-4 indicators. For a better understanding of my problem, I posted my code below. I hope you can give me some advice on how to solve my problems.
I am looking forward for your response.
Best regards,
Yuan
import numpy as np
import biogeme.database as db
import biogeme.models as models
import biogeme.biogeme as bio
import biogeme.distributions as dist
from biogeme.expressions import *
database = db.Database("helmet", df)
df.head()
globals().update(database.variables)
### Coefficients
# Read the estimates from the structural equation estimation and use them as starting values
structResults = res.bioResults(pickleFile='Latent.pickle')
structBetas = structResults.getBetaValues()
coef_intercept_1 = Beta('intercept_1', structBetas['intercept_1'], None, None, 0)
coef_gender_1 = Beta('coef_gender_1', structBetas['coef_gender_1'], None, None, 0)
coef_age_1 = Beta('coef_age_1', structBetas['coef_age_1'], None, None, 0)
coef_marriage_1 = Beta('coef_marriage_1', structBetas['coef_marriage_1'], None, None, 0)
coef_education_1 = Beta('coef_education_1', structBetas['coef_income_1'], None, None, 0)
coef_income_1 = Beta('coef_income_1', structBetas['coef_income_1'], None, None, 0)
coef_intercept_2 = Beta('intercept_2', structBetas['intercept_2'], None, None, 0)
coef_gender_2 = Beta('coef_gender_2', structBetas['coef_gender_2'], None, None, 0)
coef_age_2 = Beta('coef_age_2', structBetas['coef_age_2'], None, None, 0)
coef_marriage_2 = Beta('coef_marriage_2', structBetas['coef_marriage_2'], None, None, 0)
coef_education_2 = Beta('coef_education_2', structBetas['coef_income_2'], None, None, 0)
coef_income_2 = Beta('coef_income_2', structBetas['coef_income_2'], None, None, 0)
coef_intercept_3 = Beta('intercept_3', structBetas['intercept_3'], None, None, 0)
coef_gender_3 = Beta('coef_gender_3', structBetas['coef_gender_3'], None, None, 0)
coef_age_3 = Beta('coef_age_3', structBetas['coef_age_3'], None, None, 0)
coef_marriage_3 = Beta('coef_marriage_3', structBetas['coef_marriage_3'], None, None, 0)
coef_education_3 = Beta('coef_education_3', structBetas['coef_income_3'], None, None, 0)
coef_income_3 = Beta('coef_income_3', structBetas['coef_income_3'], None, None, 0)
coef_intercept_4 = Beta('intercept_4', structBetas['intercept_4'], None, None, 0)
coef_gender_4 = Beta('coef_gender_4', structBetas['coef_gender_4'], None, None, 0)
coef_age_4 = Beta('coef_age_4', structBetas['coef_age_4'], None, None, 0)
coef_marriage_4 = Beta('coef_marriage_4', structBetas['coef_marriage_4'], None, None, 0)
coef_education_4 = Beta('coef_education_4', structBetas['coef_income_4'], None, None, 0)
coef_income_4 = Beta('coef_income_4', structBetas['coef_income_4'], None, None, 0)
coef_intercept_5 = Beta('intercept_5', structBetas['intercept_5'], None, None, 0)
coef_gender_5 = Beta('coef_gender_5', structBetas['coef_gender_5'], None, None, 0)
coef_age_5 = Beta('coef_age_5', structBetas['coef_age_5'], None, None, 0)
coef_marriage_5 = Beta('coef_marriage_5', structBetas['coef_marriage_5'], None, None, 0)
coef_education_5 = Beta('coef_education_5', structBetas['coef_income_5'], None, None, 0)
coef_income_5 = Beta('coef_income_5', structBetas['coef_income_5'], None, None, 0)
### LATENT VARIABLE: STRUCTURAL MODEL
# 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_s1 = Beta('sigma_s1', 1, None, None, 0)
sigma_s2 = Beta('sigma_s2', 1, None, None, 0)
sigma_s3 = Beta('sigma_s3', 1, None, None, 0)
sigma_s4 = Beta('sigma_s4', 1, None, None, 0)
sigma_s5 = Beta('sigma_s5', 1, None, None, 0)
PBE = coef_intercept_1 + coef_age_1 * age + coef_gender_1 * gender + coef_education_1 * education + coef_income_1 * income + coef_marriage_1 * marriage + sigma_s1 * omega
PBA = coef_intercept_2 + coef_age_2 * age + coef_gender_2 * gender + coef_education_2 * education + coef_income_2 * income + coef_marriage_2 * marriage + sigma_s2 * omega
PSU = coef_intercept_3 + coef_age_3 * age + coef_gender_3 * gender + coef_education_3 * education + coef_income_3 * income + coef_marriage_3 * marriage + sigma_s3 * omega
PSE = coef_intercept_4 + coef_age_4 * age + coef_gender_4 * gender + coef_education_4 * education + coef_income_4 * income + coef_marriage_4 * marriage + sigma_s4 * omega
NML = coef_intercept_5 + coef_age_5 * age + coef_gender_5 * gender + coef_education_5 * education + coef_income_5 * income + coef_marriage_5 * marriage + sigma_s5 * omega
### Measurement equations
INTER_Pbe1 = Beta('INTER_Pbe1', structBetas['INTER_Pbe1'], None, None, 0)
INTER_Pbe2 = Beta('INTER_Pbe2', structBetas['INTER_Pbe2'], None, None, 0)
INTER_Pbe3 = Beta('INTER_Pbe3', structBetas['INTER_Pbe3'], None, None, 0)
INTER_Pba1 = Beta('INTER_Pba1', structBetas['INTER_Pba1'], None, None, 0)
INTER_Pba2 = Beta('INTER_Pba2', structBetas['INTER_Pba2'], None, None, 0)
INTER_Pba3 = Beta('INTER_Pba3', structBetas['INTER_Pba3'], None, None, 0)
INTER_Pba4 = Beta('INTER_Pba4', structBetas['INTER_Pba4'], None, None, 0)
INTER_Psu1 = Beta('INTER_Psu1', structBetas['INTER_Psu1'], None, None, 0)
INTER_Psu2 = Beta('INTER_Psu2', structBetas['INTER_Psu2'], None, None, 0)
INTER_Psu3 = Beta('INTER_Psu3', structBetas['INTER_Psu3'], None, None, 0)
INTER_Pse1 = Beta('INTER_Pse1', structBetas['INTER_Pse1'], None, None, 0)
INTER_Pse2 = Beta('INTER_Pse2', structBetas['INTER_Pse2'], None, None, 0)
INTER_Pse3 = Beta('INTER_Pse3', structBetas['INTER_Pse3'], None, None, 0)
INTER_Pse4 = Beta('INTER_Pse4', structBetas['INTER_Pse4'], None, None, 0)
INTER_Nml1 = Beta('INTER_Nml1', structBetas['INTER_Nml1'], None, None, 0)
INTER_Nml2 = Beta('INTER_Nml2', structBetas['INTER_Nml2'], None, None, 0)
INTER_Nml3 = Beta('INTER_Nml3', structBetas['INTER_Nml3'], None, None, 0)
INTER_Nml4 = Beta('INTER_Nml4', structBetas['INTER_Nml4'], None, None, 0)
B_Pbe1 = Beta('B_Pbe1', structBetas['B_Pbe1'], None, None, 0)
B_Pbe2 = Beta('B_Pbe2', structBetas['B_Pbe2'], None, None, 0)
B_Pbe3 = Beta('B_Pbe3', structBetas['B_Pbe3'], None, None, 0)
B_Pba1 = Beta('B_Pba1', structBetas['B_Pba1'], None, None, 0)
B_Pba2 = Beta('B_Pba2', structBetas['B_Pba2'], None, None, 0)
B_Pba3 = Beta('B_Pba3', structBetas['B_Pba3'], None, None, 0)
B_Pba4 = Beta('B_Pba4', structBetas['B_Pba4'], None, None, 0)
B_Psu1 = Beta('B_Psu1', structBetas['B_Psu1'], None, None, 0)
B_Psu2 = Beta('B_Psu2', structBetas['B_Psu2'], None, None, 0)
B_Psu3 = Beta('B_Psu3', structBetas['B_Psu3'], None, None, 0)
B_Pse1 = Beta('B_Pse1', structBetas['B_Pse1'], None, None, 0)
B_Pse2 = Beta('B_Pse2', structBetas['B_Pse2'], None, None, 0)
B_Pse3 = Beta('B_Pse3', structBetas['B_Pse3'], None, None, 0)
B_Pse4 = Beta('B_Pse4', structBetas['B_Pse4'], None, None, 0)
B_Nml1 = Beta('B_Nml1', structBetas['B_Nml1'], None, None, 0)
B_Nml2 = Beta('B_Nml2', structBetas['B_Nml2'], None, None, 0)
B_Nml3 = Beta('B_Nml3', structBetas['B_Nml3'], None, None, 0)
B_Nml4 = Beta('B_Nml4', structBetas['B_Nml4'], None, None, 0)
MODEL_Pbe1 = INTER_Pbe1 + B_Pbe1 * PBE
MODEL_Pbe2 = INTER_Pbe2 + B_Pbe2 * PBE
MODEL_Pbe3 = INTER_Pbe3 + B_Pbe3 * PBE
MODEL_Pba1 = INTER_Pba1 + B_Pba1 * PBA
MODEL_Pba2 = INTER_Pba2 + B_Pba2 * PBA
MODEL_Pba3 = INTER_Pba3 + B_Pba3 * PBA
MODEL_Pba4 = INTER_Pba4 + B_Pba4 * PBA
MODEL_Psu1 = INTER_Psu1 + B_Psu1 * PSU
MODEL_Psu2 = INTER_Psu2 + B_Psu2 * PSU
MODEL_Psu3 = INTER_Psu3 + B_Psu3 * PSU
MODEL_Pse1 = INTER_Pse1 + B_Pse1 * PSE
MODEL_Pse2 = INTER_Pse2 + B_Pse2 * PSE
MODEL_Pse3 = INTER_Pse3 + B_Pse3 * PSE
MODEL_Pse4 = INTER_Pse4 + B_Pse4 * PSE
MODEL_Nml1 = INTER_Nml1 + B_Nml1 * NML
MODEL_Nml2 = INTER_Nml2 + B_Nml2 * NML
MODEL_Nml3 = INTER_Nml3 + B_Nml3 * NML
MODEL_Nml4 = INTER_Nml4 + B_Nml4 * NML
SIGMA_STAR_Pbe1 = Beta('SIGMA_STAR_Pbe1', structBetas['SIGMA_STAR_Pbe1'], None, None, 1)
SIGMA_STAR_Pbe2 = Beta('SIGMA_STAR_Pbe2', structBetas['SIGMA_STAR_Pbe2'], None, None, 1)
SIGMA_STAR_Pbe3 = Beta('SIGMA_STAR_Pbe3', structBetas['SIGMA_STAR_Pbe3'], None, None, 1)
SIGMA_STAR_Pba1 = Beta('SIGMA_STAR_Pba1', structBetas['SIGMA_STAR_Pba1'], None, None, 1)
SIGMA_STAR_Pba2 = Beta('SIGMA_STAR_Pba2', structBetas['SIGMA_STAR_Pba2'], None, None, 1)
SIGMA_STAR_Pba3 = Beta('SIGMA_STAR_Pba3', structBetas['SIGMA_STAR_Pba3'], None, None, 1)
SIGMA_STAR_Pba4 = Beta('SIGMA_STAR_Pba4', structBetas['SIGMA_STAR_Pba4'], None, None, 1)
SIGMA_STAR_Psu1 = Beta('SIGMA_STAR_Psu1', structBetas['SIGMA_STAR_Psu1'], None, None, 0)
SIGMA_STAR_Psu2 = Beta('SIGMA_STAR_Psu2', structBetas['SIGMA_STAR_Psu2'], None, None, 0)
SIGMA_STAR_Psu3 = Beta('SIGMA_STAR_Psu3', structBetas['SIGMA_STAR_Psu3'], None, None, 0)
SIGMA_STAR_Pse1 = Beta('SIGMA_STAR_Pse1', structBetas['SIGMA_STAR_Pse1'], None, None, 1)
SIGMA_STAR_Pse2 = Beta('SIGMA_STAR_Pse2', structBetas['SIGMA_STAR_Pse2'], None, None, 1)
SIGMA_STAR_Pse3 = Beta('SIGMA_STAR_Pse3', structBetas['SIGMA_STAR_Pse3'], None, None, 1)
SIGMA_STAR_Pse4 = Beta('SIGMA_STAR_Pse4', structBetas['SIGMA_STAR_Pse4'], None, None, 1)
SIGMA_STAR_Nml1 = Beta('SIGMA_STAR_Nml1', structBetas['SIGMA_STAR_Nml1'], None, None, 1)
SIGMA_STAR_Nml2 = Beta('SIGMA_STAR_Nml2', structBetas['SIGMA_STAR_Nml2'], None, None, 1)
SIGMA_STAR_Nml3 = Beta('SIGMA_STAR_Nml3', structBetas['SIGMA_STAR_Nml3'], None, None, 1)
SIGMA_STAR_Nml4 = Beta('SIGMA_STAR_Nml4', structBetas['SIGMA_STAR_Nml4'], None, None, 0)
### CHOICE MODEL
choiceResults = res.bioResults(pickleFile='Baseline_w_social.pickle')
choiceBetas = choiceResults.getBetaValues()
ASC_WEARINGHELMET = Beta('ASC_WEARINGHELMET', choiceBetas['ASC_WEARINGHELMET'], None, None, 0)
ASC_NON_WEARING = Beta('ASC_NON_WEARING', 0, None, None, 1)
BETA_PENALTYPRICE = Beta('BETA_PENALTYPRICE', choiceBetas['BETA_PENALTYPRICE'], None, None, 0)
BETA_PENALTYCHANCE = Beta('BETA_PENALTYCHANCE', choiceBetas['BETA_PENALTYCHANCE'], None, None, 0)
BETA_HELMETPERCENT = Beta('BETA_HELMETPERCENT', choiceBetas['BETA_HELMETPERCENT'], None, None, 0)
BETA_AGE = Beta('BETA_AGE', choiceBetas['BETA_AGE'], None, None, 0)
BETA_MARRIAGE = Beta('BETA_MARRIAGE', choiceBetas['BETA_MARRIAGE'], None, None, 0)
BETA_GENDER = Beta('BETA_GENDER', choiceBetas['BETA_GENDER'], None, None, 0)
BETA_INCOME = Beta('BETA_INCOME', choiceBetas['BETA_INCOME'], None, None, 0)
BETA_EDUCATION = Beta('BETA_EDUCATION', choiceBetas['BETA_EDUCATION'], None, None, 0)
BETA_PBE = Beta('BETA_PBE', 0, None, None, 0)
BETA_PBA = Beta('BETA_PBA', 0, None, None, 0)
BETA_PSU = Beta('BETA_PSU', 0, None, None, 0)
BETA_PSE = Beta('BETA_PSE', 0, None, None, 0)
BETA_NML = Beta('BETA_NML', 0, None, None, 0)
### DEFINITION OF UTILITY FUNCTIONS:
V1 = ASC_WEARINGHELMET + BETA_PENALTYPRICE * penaltyprice + BETA_PENALTYCHANCE * penaltychance + BETA_HELMETPERCENT * helmetpercent +\
BETA_GENDER * gender + BETA_EDUCATION * education + BETA_INCOME * income + BETA_MARRIAGE * marriage + BETA_AGE * age +\
BETA_PBA * PBA + BETA_PBE * PBE + BETA_PSU * PSU + BETA_PSE * PSE + BETA_NML * NML
V2 = ASC_NON_WEARING
# Associate utility functions with the numBering of alternatives
V = {0: V1,
1: V2}
# In this example all alternatives are available for each individual.
av = {0: 1,
1: 1}
# The choice model is a logit, conditional to the value of the latent variable
prob = models.logit(V, av, helmetchoice)
logprob = log(MonteCarlo(prob))
biogeme = bio.BIOGEME(database, logprob, numberOfDraws = 500)
biogeme.modelName = "LatentChoice_Full_MonteCarlo_Estimation"
results = biogeme.estimate()
print(f"Estimating 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}")
Michel Bierlaire
If you want to use Monte-Carlo, replace RandomVariable by bioDraws. See the documentation and examples for details.
> You received this message because you are subscribed to the Google Groups "Biogeme" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to biogeme+u...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/biogeme/ea0a7c9e-b704-4d37-a916-7a9d18ca33f7n%40googlegroups.com.
Michel Bierlaire
Transport and Mobility Laboratory
School of Architecture, Civil and Environmental Engineering
EPFL - Ecole Polytechnique Fédérale de Lausanne
http://transp-or.epfl.ch
http://people.epfl.ch/michel.bierlaire
Michel Bierlaire
If you want to use Monte-Carlo, replace RandomVariable by bioDraws. See the documentation and examples for details.
> On 9 Apr 2024, at 10:12, Yuan Wang <zmeng...@gmail.com> wrote:
>