a LC-NL model

166 views
Skip to first unread message

Dexue kong

unread,
Oct 14, 2021, 4:52:38 AM10/14/21
to Biogeme
Dear Professor , 
      I have written a LC-NL model using the " swissmetro" data, the model calculation results are shown in the following text,the estimated value of each unknown parameter is 0I don't know what's wrong with my code?
————
Traceback (most recent call last):
  File "D:/python/Workspace/PycharmProjects/test/LCMbiogeme/nestchuli/LCNL_Try.py", line 145, in <module>
    results = biogeme.estimate()
  File "D:\python\Anaconda3\lib\site-packages\biogeme\biogeme.py", line 701, in estimate
    self.calculateInitLikelihood()
  File "D:\python\Anaconda3\lib\site-packages\biogeme\biogeme.py", line 409, in calculateInitLikelihood
    scaled=False)
  File "D:\python\Anaconda3\lib\site-packages\biogeme\biogeme.py", line 445, in calculateLikelihood
    f = self.theC.calculateLikelihood(x, self.fixedBetaValues)
  File "src\cbiogeme.pyx", line 97, in biogeme.cbiogeme.pyBiogeme.calculateLikelihood
RuntimeError: src/bioExprLog.cc:61: Biogeme exception: Current values of the literals:
ASC_CAR_class0 = 0
ASC_CAR_class1 = 0
ASC_CAR_class2 = 0
ASC_SM_class0 = 0
ASC_SM_class1 = 0
ASC_SM_class2 = 0
ASC_TRAIN_class0 = 0
ASC_TRAIN_class1 = 0
ASC_TRAIN_class2 = 0
B_COST_class0 = 0
B_COST_class1 = 0
B_COST_class2 = 0
B_TIME_class0 = 0
B_TIME_class1 = 0
B_TIME_class2 = 0
CAR_AV_SP = 0
CAR_CO_SCALED = 0
CAR_TT_SCALED = 0
CHOICE = 1
MU_class0 = 1
MU_class1 = 1
MU_class2 = 1
PROB_class0 = 0.5
PROB_class1 = 0.5
SM_AV = 1
SM_COST_SCALED = 0
SM_TT_SCALED = 0.56
TRAIN_AV_SP = 1
TRAIN_COST_SCALED = 0
TRAIN_TT_SCALED = 1.32
Cannot take the log of a non positive number [-0.693147]

——————
The code is as follows :

import pandas as pd
import biogeme.database as db
import biogeme.biogeme as bio
import biogeme.models as models
import biogeme.messaging as msg
from biogeme.expressions import Beta, DefineVariable, bioDraws, \
PanelLikelihoodTrajectory, MonteCarlo, log

# Read the data
df = pd.read_csv('swissmetro.dat', '\t')
database = db.Database('swissmetro', df)

# They are organized as panel data. The variable ID identifies each individual.
database.panel("ID")

# The following statement allows you to use the names of the variable
# as Python variable.
globals().update(database.variables)

# Here we use the "biogeme" way for backward compatibility
exclude = ((PURPOSE != 1) * (PURPOSE != 3) + (CHOICE == 0)) > 0
database.remove(exclude)

# Parameters to be estimated. One version for each latent class.
numberOfClasses = 3

B_COST = [Beta(f'B_COST_class{i}', 0, None, None, 0) for i in range(numberOfClasses)]
B_TIME = [Beta(f'B_TIME_class{i}', 0, None, None, 0) for i in range(numberOfClasses)]
ASC_CAR = [Beta(f'ASC_CAR_class{i}', 0, None, None, 0) for i in range(numberOfClasses)]
ASC_TRAIN = [Beta(f'ASC_TRAIN_class{i}', 0, None, None, 0) for i in range(numberOfClasses)]
ASC_SM = [Beta(f'ASC_SM_class{i}', 0, None, None, 1) for i in range(numberOfClasses)]

# Definition of new variables
SM_COST = SM_CO * (GA == 0)
TRAIN_COST = TRAIN_CO * (GA == 0)

# Definition of new variables: adding columns to the database
CAR_AV_SP = DefineVariable('CAR_AV_SP', CAR_AV * (SP != 0), database)
TRAIN_AV_SP = DefineVariable('TRAIN_AV_SP', TRAIN_AV * (SP != 0), database)
TRAIN_TT_SCALED = DefineVariable('TRAIN_TT_SCALED', TRAIN_TT / 100.0, database)
TRAIN_COST_SCALED = DefineVariable('TRAIN_COST_SCALED', TRAIN_COST / 100, database)
SM_TT_SCALED = DefineVariable('SM_TT_SCALED', SM_TT / 100.0, database)
SM_COST_SCALED = DefineVariable('SM_COST_SCALED', SM_COST / 100, database)
CAR_TT_SCALED = DefineVariable('CAR_TT_SCALED', CAR_TT / 100, database)
CAR_CO_SCALED = DefineVariable('CAR_CO_SCALED', CAR_CO / 100, database)

# Utility functions
V1 = [ASC_TRAIN[i] + B_TIME[i] * TRAIN_TT_SCALED + B_COST[i] * TRAIN_COST_SCALED
for i in range(numberOfClasses)]
V2 = [ASC_SM[i] + B_TIME[i] * SM_TT_SCALED + B_COST[i] * SM_COST_SCALED
for i in range(numberOfClasses)]
V3 = [ASC_CAR[i] + B_TIME[i] * CAR_TT_SCALED + B_COST[i] * CAR_CO_SCALED
for i in range(numberOfClasses)]
V = [{1: V1[i],
2: V2[i],
3: V3[i]} for i in range(numberOfClasses)]

# Associate the availability conditions with the alternatives
av = {1: TRAIN_AV_SP,
2: SM_AV,
3: CAR_AV_SP}

# Definition of nests:
# 1: nests parameter
# 2: list of alternatives
MU = [Beta(f'MU_class{i}', 1, 1, 10, 0) for i in range(numberOfClasses)]

existing = [[MU[i], [1, 3] for i in range(numberOfClasses)]
future = 1.0, [2]
nests = [[existing[i], future] for i in range(numberOfClasses)]

# The choice model is a discrete mixture of logit, with availability conditions
# We calculate the conditional probability for each class
prob = [models.lognested(V[i], av, nests[i], CHOICE) for i in range(numberOfClasses)]

# Class membership model
# Parameters for the class membership model
CLASS_CTE = [Beta(f'CLASS_CTE_class{i}', 0, None, None, 0) for i in range(numberOfClasses)]
CLASS_INC = [Beta(f'CLASS_INC_class{i}', 0, None, None, 0) for i in range(numberOfClasses)]

CLASS_CTE[2] = 0
CLASS_INC[2] = 0
W = [CLASS_CTE[i] + CLASS_INC[i] * INCOME for i in range(numberOfClasses)]
PROB_class = [models.logit({0: W[0], 1: W[1], 2: W[2]}, None, 0) for i in range(numberOfClasses)]

# Conditional to the random variables, likelihood for the individual.
probIndiv = PROB_class[0] * prob[0] + PROB_class[1] * prob[1] + PROB_class[2] * prob[2]

# We integrate over the random variables using Monte-Carlo
logprob = log(probIndiv)

# Define level of verbosity
logger = msg.bioMessage()
#logger.setSilent()
#logger.setWarning()
logger.setGeneral()
#logger.setDetailed()

# Create the Biogeme object
biogeme = bio.BIOGEME(database, logprob)
biogeme.modelName = '2_LCNL_kdx'

# Estimate the parameters
results = biogeme.estimate()
pandasResults = results.getEstimatedParameters()
print(pandasResults)

Bierlaire Michel

unread,
Oct 14, 2021, 4:55:52 AM10/14/21
to cqj...@gmail.com, Bierlaire Michel, Biogeme
The error message is clear: you are trying to take the log of a negative number. 
In order to discover where the problem comes from, use the simulate function of biogeme to report the various quantities involved in your expression.


--
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/82137407-4e34-481f-b83b-792f35251124n%40googlegroups.com.

Reply all
Reply to author
Forward
0 new messages