Update code from previous BIOGEME version

[nata...@gmail.com]

Dear Professor Bierlaire,
I am inquiring about updating a code written for BIOGEME.

In 2014, I developed a reasonably working model using an earlier version of BIOGEME. However, with the recent changes in the software, I find myself quite perplexed about how to update the code effectively. The code is rather complex (see below) as it addresses an inherent endogeneity issue associated with stated preference surveys, particularly when a hypothetical choice is conditioned on previous choice responses.

I hope you can provide guidance or insights on effectively updating the code for the newer version of BIOGEME. Your assistance in this matter would be greatly appreciated.

Best wishes,
Natan


CODE

from biogeme import *
from headers import *
from loglikelihood import *
from statistics import *
from distributions import *

#Constants
C_P1SI             = Beta('C_P1SI',0,-10,10,0)
C_P1NO             = Beta('C_P1NO',0,-10,10,1)

C_P1SI_CONIS       = Beta('C_P1SI_CONIS',0,-10,10,0)
C_P1NO_SICONIS     = Beta('C_P1NO_SICONIS',0,-10,10,0)
C_P1NO_NOCON       = Beta('C_P1NO_NOCON',0,-10,10,1)

SIGMA_eP1NO          = Beta('SIGMA_eP1NO',1,-10,10000,0)
SIGMA_eP1SI          = Beta('SIGMA_eP1SI',1,-10,10000,0)

# Best-Worst variables
qBCC               = Beta('qBCC',0,-20,20,0)
qBCS               = Beta('qBCS',0,-20,20,0)
qBTL               = Beta('qBTL',0,-20,20,1)      #Dejo fijo este parámetro
qCON               = Beta('qCON',0,-20,20,0)
qDEN               = Beta('qDEN',0,-20,20,0)
qEST               = Beta('qEST',0,-20,20,0)
qIPU               = Beta('qIPU',0,-20,20,0)

tBCC               = Beta('tBCC',1,-10,10,0)
tBCS               = Beta('tBCS',1,-10,10,0)
tBTL               = Beta('tBTL',1,-10,10,0)
tCON               = Beta('tCON',1,-10,10,0)
tDEN               = Beta('tDEN',1,-10,10,0)
tEST               = Beta('tEST',1,-10,10,0)
tIPU               = Beta('tIPU',1,-10,10,0)

ThethaB         = Beta('ThethaB',1,0,100,0)
ThethaW            = Beta('ThethaW',1,0,100,0)
ThethaED            = Beta('ThethaED',0,0,10,1)     
ThethaEDP1SI            = Beta('ThethaEDP1SI',0,0,10000,0)    
ThethaEDP1NO            = Beta('ThethaEDP1NO',0,0,10000,0)    

# Utility functions
# Best
V1= ThethaB * (qBCC + tBCC * L_BCC)
V2= ThethaB * (qBCS + tBCS * L_BCS)
V3= ThethaB * (qBTL + tBTL * L_BTL)
V4= ThethaB * (qCON + tCON * L_CON)
V5= ThethaB * (qDEN)
V6= ThethaB * (qEST)
V7= ThethaB * (qIPU + tIPU * L_IPU)
#Worst
V8= ThethaW * (qBCC * C_BCC + tBCC * L_BCC)
V9= ThethaW * (qBCS * C_BCS + tBCS * L_BCS)
V10= ThethaW * (qBTL * C_BTL + tBTL * L_BTL)
V11= ThethaW * (qCON * C_CON + tCON * L_CON)
V12= ThethaW * (qDEN * C_DEN + tDEN * L_DEN)
V13= ThethaW * (qEST * C_EST + tEST * L_EST)
V14= ThethaW * (qIPU * C_IPU + tIPU * L_IPU)
# ED
V15 = C_P1NO
V16 = C_P1SI + tBCC * L_BCC + tBCS * L_BCS + tBTL * L_BTL + tCON * L_CON +  tDEN * L_DEN + tEST * L_EST + tIPU * L_IPU

# Utilities
V1 = {1: V1,
     2: V2,
     3: V3,
     4: V4,
     5: V5,
     6: V6,
     7: V7,
     8: V8,
     9: V9,
     10: V10,
     11: V11,
     12: V12,
     13: V13,
     14: V14,
     15: V15,
     16: V16}

V1_aux = {15: V15,
          16: V16}

#V[15]

# AV
av1 = {1: DB_BCC,
      2: DB_BCS,
      3: DB_BTL,
      4: DB_CON,
      5: DB_DEN,
      6: DB_EST,
      7: DB_IPU,
      8: DW_BCC,
      9: DW_BCS,
      10: DW_BTL,
      11: DW_CON,
      12: DW_DEN,
      13: DW_EST,
      14: DW_IPU,
      15: D_P1_NO,
      16: D_P1_SI}

av2 = {17: D_SI_P2_SINIS,
      18: D_SI_P2_CONIS,
      19: D_NO_P2_CONIS,
      20: D_NO_P2_NOCON}

av_aux = {15: 1,
          16: 1}

av = {1: DB_BCC,
      2: DB_BCS,
      3: DB_BTL,
      4: DB_CON,
      5: DB_DEN,
      6: DB_EST,
      7: DB_IPU,
      8: DW_BCC,
      9: DW_BCS,
      10: DW_BTL,
      11: DW_CON,
      12: DW_DEN,
      13: DW_EST,
      14: DW_IPU,
      15: D_P1_NO,
      16: D_P1_SI,
      17: D_SI_P2_SINIS,
      18: D_SI_P2_CONIS,
      19: D_NO_P2_CONIS,
      20: D_NO_P2_NOCON}

P1 = bioLogit(V1,av1,CHOICE)

mu_P1_i = bioUniformDraws('mu_P1_i')
P1_aux = bioLogit(V1_aux,av_aux,CHOICE_aux)

error_P1SI_i = - log(P1_aux) - log(-log(mu_P1_i))

# P1SI
V17 = ThethaEDP1SI * (tBCC * L_BCC + tBCS * L_BCS + tBTL * L_BTL + tCON * L_CON +  tDEN * L_DEN + tEST * L_EST + tIPU * L_IPU + SIGMA_eP1SI * error_P1SI_i )

V18 = ThethaEDP1SI * (C_P1SI_CONIS + tBCC * L_BCC_IS + tBCS * L_BCS_IS + tBTL * L_BTL_IS + tCON * L_CON_IS + (tDEN * (L_DEN_IS * (1 - INC_B4) + B4_DEN * INC_B4)) + tEST * L_EST_IS + tIPU * L_IPU_IS)
 
# P1NO
V19 = ThethaEDP1NO * (C_P1NO_SICONIS + tBCC * L_BCC_IS + tBCS * L_BCS_IS + tBTL * L_BTL_IS + tCON * L_CON_IS + tDEN * L_DEN_IS + tEST * L_EST_IS + tIPU * L_IPU_IS)

V20 = ThethaEDP1NO * (SIGMA_eP1NO * error_P1SI_i)

V2 = {17: V17,
      18: V18,
      19: V19,
      20: V20}

# Probabilities P1SI y P1NO
P2 = bioLogit(V2,av2,CHOICE)

drawIterator('drawIter')
l_P2  = (Sum(P2,   'drawIter'))

P_FINAL = log(P1 + l_P2)

#Model
rowIterator('obsIter')
BIOGEME_OBJECT.ESTIMATE = Sum(P_FINAL,'obsIter')

exclude = (INC_ED_M == 1) + (INC_ED_C == 1) + ( INC_W_NP == 1 ) + ( INC_W_NA == 1 ) + ( INC_B_NP == 1 ) + ( INC_B_0 == 1 )  > 0
BIOGEME_OBJECT.EXCLUDE = exclude

# Statistics
nullLoglikelihood(av,'obsIter')
choiceSet = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20]
cteLoglikelihood(choiceSet,CHOICE,'obsIter')
availabilityStatistics(av,'obsIter')

BIOGEME_OBJECT.PARAMETERS['RandomDistribution'] = "MLHS"
BIOGEME_OBJECT.PARAMETERS['NbrOfDraws'] = "1000"
BIOGEME_OBJECT.PARAMETERS['dumpDrawsOnFile'] = "1"
BIOGEME_OBJECT.PARAMETERS['optimizationAlgorithm'] = "CFSQP"
BIOGEME_OBJECT.PARAMETERS['checkDerivatives'] = "1"
BIOGEME_OBJECT.PARAMETERS['numberOfThreads'] = "4"

[Michel Bierlaire]

This is indeed quite an old version of Biogeme. But most of what you have done is still valid.
1) Install the latest version. https://biogeme.epfl.ch/#install
2) Prepare a script that imports your data: https://biogeme.epfl.ch/sphinx/auto_examples/swissmetro/swissmetro_data.html#swissmetro-data
3) Look at the online examples to adapt your code to the newest version:
https://biogeme.epfl.ch/sphinx/auto_examples/swissmetro/plot_b01logit.html#sphx-glr-auto-examples-swissmetro-plot-b01logit-py

Other examples are available at
https://biogeme.epfl.ch/sphinx/auto_examples/index.html#biogeme-examples-for-the-swissmetro-data

> --
> 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/80205bf1-39ee-4586-803d-d7a53d30c901n%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

[nata...@gmail.com]

Many thanks for your quick response. 
I managed to update a large part of the code, yet I'm still confused on how to update with a few lines: 

1) mu_P1_i = bioUniformDraws('mu_P1_i')
2) error_P1SI_i = - log(P1_aux) - log(-log(mu_P1_i)) --> the expression log is still the same in the new version?
3) drawIterator('drawIter')
4) rowIterator('obsIter')
5) BIOGEME_OBJECT.ESTIMATE = Sum(P_FINAL,'obsIter')

For the above lines of code, which is the corresponding code for the current version of Biogeme?

Best wishes,

N.

[Michel Bierlaire]

> On 8 May 2024, at 15:09, nata...@gmail.com <nata...@gmail.com> wrote:
>
> Many thanks for your quick response.
> I managed to update a large part of the code, yet I'm still confused on how to update with a few lines:
>

> 1) mu_P1_i = bioUniformDraws('mu_P1_i').

A detailed discussion about all types of draws is available here: https://transp-or.epfl.ch/documents/technicalReports/Bier19.pdf
In your case:

from biogeme.expressions import bioDraws
mu_P1_i = bioDraws(’mu_P1_i’,’UNIFORM’)

> 2) error_P1SI_i = - log(P1_aux) - log(-log(mu_P1_i)) --> the expression log is still the same in the new version?

Yes. But you need to import it:

from biogeme.expressions import log


> 3) drawIterator('drawIter')

--> Not necessary anymore. Replaced by the expression 'MonteCarlo'. See the above report.

> 4) rowIterator('obsIter')

--> Not necessary anymore.

> 5) BIOGEME_OBJECT.ESTIMATE = Sum(P_FINAL,'obsIter')

--> not necessary. The new syntax is:

the_biogeme = bio.BIOGEME(database, logprob)
the_biogeme.modelName = 'my_model'
the_biogeme.calculateNullLoglikelihood(av). [if needed]
results = the_biogeme.estimate()

> To view this discussion on the web visit https://groups.google.com/d/msgid/biogeme/438e26b6-43a4-4a13-aad6-866d9ee2b23cn%40googlegroups.com.

[nata...@gmail.com]

Many thanks for your response. I still see an error that I cannot trace, but I will check the references provided in detail.

Best wishes,

Natan