How to simulate Mixed Logit Model in Pandasbiogeme

[CASA K]

Dear All,

Anyone knows how to simulate a normal mixed logit model in Pandasbiogeme. It seems the Swissmetro example only gives the logit model simulation.

Thanks,

Kevin

[Michel Bierlaire]

This example is a normal mixture of logit: http://biogeme.epfl.ch/examples/swissmetro/pandas/05normalMixture.py

--
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 post to this group, send email to bio...@googlegroups.com.
Visit this group at https://groups.google.com/group/biogeme.
For more options, visit https://groups.google.com/d/optout.

[Michel Bierlaire]

Well, this is exactly the same, irrespectively of the model. You need to define a “simulate” dictionary, with all the quantities that you want to simulate, and give it as an argument when you create an instance of Biogeme. 

You will find other examples here http://biogeme.epfl.ch/examples_indicators.html

On 13 Mar 2019, at 19:18, Kevin Xie <kevin.zho...@gmail.com> wrote:

Dear Professor,

Thanks for replaying. I might not explain it very clearly. I mean how to simulate the modelling result to make predictions. For example, there is a 01logit_simul.py to allow me simulate the result of 01logit.py. Would there be something similar to simulate the result from the parameters got from 05logit.py?

Thanks,

Kevin

[Kevin Xie]

Dear Professor,

Thanks for replaying. I might not explain it very clearly. I mean how to simulate the modelling result to make predictions. For example, there is a 01logit_simul.py to allow me simulate the result of 01logit.py. Would there be something similar to simulate the result from the parameters got from 05logit.py?

Thanks,

Kevin

On 13 Mar 2019, at 17:07, Michel Bierlaire <michel.b...@epfl.ch> wrote:

[Kevin Xie]

Hi Professor,

I have tried to simulate the mixed result. However, feels a little bit confused on the bioDraws part.

I have tried the below script and got a result with all predictions NaN (No error occured). I think the problem should be around the bioDraws and MonteCarlo simulation. Could you please help me address the issue?

Thanks,

Kevin

%%time

import biogeme.version as bv

import sys

import biogeme.database as db

import biogeme.biogeme as bio

import biogeme.models as models

import biogeme.distributions as dist

from headers import *

import io

import json

from biogeme.expressions import *

#from headers import *

# script

for i in range(testing_max):

    exec("av"+str(i)+" = Variable('av'+str(i))")

ID=Variable('ID')

choice_l=Variable('choice_l')

B_ACT_FARE = Beta('B_ACT_FARE',-1.79819206e-03,None,None,0,'actual fare')

B_BUS_LEG = Beta('B_BUS_LEG',-1.95789632e+01,None,None,0,'bus leg length')

B_OFP_FARE = Beta('B_OFP_FARE',-1.79819206e-03,None,None,0,'out of pocket fare')

B_SPK_COV_END = Beta('B_SPK_COV_END',1.15987730e+01,None,None,0,'spatial knowledge coverage ends')

B_TRANSFER = Beta('B_TRANSFER',-3.98613806e-01,None,None,0,'transfer')  

B_WEIGHT=Beta('B_WEIGHT',-2.36498917e-01,None,None,0,'weight')

B_WEIGHT_S=Beta('B_WEIGHT_S',8.04988829e-02,None,None,0,'weight_s')

B_WEIGHT_RND = B_WEIGHT + B_WEIGHT_S * bioDraws('B_WEIGHT_RND','NORMAL') 

#B_TRAVELCARD = Beta('B_TRAVELCARD',0,None,None,0,'travelcard')

#B_FREEPASS = Beta('B_FREEPASS',0,None,None,0,'freepass') 

# Utility functions

# attribute

weight={}

for i in range(testing_max):

    weight[i]=Variable('weight' + '_' + str( i ))

   

    

ivt={}

for i in range(testing_max):

    ivt[i]=Variable('ivt' + '_' + str( i ))

waiting_time={}

for i in range(testing_max):

    waiting_time[i]=Variable('waiting_time' + '_' + str( i ))

walking_time={}

for i in range(testing_max):

    walking_time[i]=Variable('walking_time' + '_' + str( i ))

intermediate_walking_time={}

for i in range(testing_max):

    intermediate_walking_time[i]=Variable('intermediate_walking_time' + '_' + str( i ))

transfer={}

for i in range(testing_max):

    transfer[i]=Variable('transfer' + '_' + str( i ))

bus_leg_length={}

for i in range(testing_max):

    bus_leg_length[i]=Variable('bus_leg_length' + '_' + str( i ))

    

bus_leg_ratio={}

for i in range(testing_max):

    bus_leg_ratio[i]=Variable('bus_leg_ratio' + '_' + str( i ))

first_leg_length={}

for i in range(testing_max):

    first_leg_length[i]=Variable('first_leg_length' + '_' + str( i ))

first_leg_ratio={}

for i in range(testing_max):

    first_leg_ratio[i]=Variable('first_leg_ratio' + '_' + str( i ))

spk_cov_end={}

for i in range(testing_max):

    spk_cov_end[i]=Variable('spk_cov_end' + '_' + str( i ))

    

spk_cov_pass={}

for i in range(testing_max):

    spk_cov_pass[i]=Variable('spk_cov_pass' + '_' + str( i ))

out_of_pocket_fare={}

for i in range(testing_max):

    out_of_pocket_fare[i]=Variable('out_of_pocket_fare' + '_' + str( i ))

actual_fare={}

for i in range(testing_max):

    actual_fare[i]=Variable('actual_fare' + '_' + str( i ))

travelcard={}

for i in range(testing_max):

    travelcard[i]=Variable('travelcard' + '_' + str( i ))

freepass={}

for i in range(testing_max):

    freepass[i]=Variable('freepass' + '_' + str( i ))

    

V = {}

for i in range(testing_max):

    V[i] = B_WEIGHT_RND * weight[i] +\

        B_TRANSFER * transfer[i] +B_BUS_LEG * bus_leg_length[i] +  B_ACT_FARE * actual_fare[i]  +  B_SPK_COV_END * spk_cov_end[i] + B_OFP_FARE * out_of_pocket_fare[i]

# Associate the availability conditions with the alternatives

av={}

for i in range(testing_max):

    av[i]=Variable('av'+str(i))

# The choice model is a logit, with availability conditions

# define prob

for i in range(testing_max):

    exec("logprob"+str(i)+" = log(MonteCarlo(models.logit(V, av, i)))")

    

simulate={}

for i in range(testing_max):

    simulate['loProb. '+ str(i)]=log(MonteCarlo(models.logit(V,av,i)))

    

    

#simulate={}

#for i in range(testing_max):

    #prob=None

    #logprob=None

    

    #prob = models.logit(V,av,i)

    #logprob = log(MonteCarlo(prob))

    #simulate['Prob. '+ str(i)]=logprob

    

biogeme = bio.BIOGEME(testdatabase,simulate)

biogeme.modelName = "20140213 Mixed Logit simulation"

results = biogeme.simulate()

#print("Results=",results.describe())

[Kevin Xie]

Hi Professor,

I got it. Please don’t waste time on the code.

Thanks,

Kevin

[Michel Bierlaire]

You may want to share the solution with everybody else...

[Kevin Xie]

Dear All,

Please see the simulation of a Mixed Logit Model. This may be helpful if you are going to make predictions with the estimated(simulated) parameters.

%%time

import biogeme.distributions as dist

#from headers import *

# script

for i in range(testing_max):

    exec("av"+str(i)+" = Variable('av'+str(i))")

ID=Variable('ID')

choice_l=Variable('choice_l')

B_ACT_FARE = Beta('B_ACT_FARE',-2.38448769e+013,None,None,0,'actual fare')

B_BUS_LEG = Beta('B_BUS_LEG',-2.38448769e+01,None,None,0,'bus leg length')

#B_OFP_FARE = Beta('B_OFP_FARE',-1.84693251e-03,None,None,0,'out of pocket fare')

B_SPK_COV_END = Beta('B_SPK_COV_END',1.16201109e+01,None,None,0,'spatial knowledge coverage ends')

B_TRANSFER = Beta('B_TRANSFER',-3.97498787e-01,None,None,0,'transfer')  

#B_WEIGHT=Beta('B_WEIGHT',-2.40347558e-01,None,None,0,'weight')

#B_WEIGHT_S=Beta('B_WEIGHT_S',8.59643814e-02,None,None,0,'weight_s')

#B_WEIGHT_RND = B_WEIGHT + B_WEIGHT_S #* bioDraws('B_WEIGHT_RND','NORMAL') 

B_WEIGHT_RND = np.random.normal(-2.40347558e-01,8.59643814e-02)

# Associate the availability conditions with the alternatives

av={}

for i in range(testing_max):

    av[i]=Variable('av'+str(i))

# The choice model is a logit, with availability conditions

# define prob

for i in range(testing_max):

    exec("prob"+str(i)+" = models.logit(V, av, i)")

    

simulate={}

for i in range(testing_max):

    simulate['Prob. '+ str(i)]=models.logit(V,av,i)

    

    

#simulate={}

#for i in range(testing_max):

    #prob=None

    #logprob=None

    

    #prob = models.logit(V,av,i)

    #logprob = log(MonteCarlo(prob))

    #simulate['Prob. '+ str(i)]=logprob

    

biogeme = bio.BIOGEME(testdatabase,simulate)

biogeme.modelName = "20140213 Mixed Logit simulation"

simresults = biogeme.simulate()

#print("Results=",results.describe())

All the best,

Kevin