biogeme-nest_logit(交通方式选择)

biogeme-nest_logit

基础数据:

optima.dat

image-20220930090135788

 

变量的描述:出处

  • OccupStat:职业

  • TimePT:公共交通通行时间

  • TimeCar:小汽车通行时间

  • MarginalCostPT:公共交通总成本

  • CostCarCHF:小汽车的总汽油成本

  • distance_km:距离

  • choice:Choice variable: 0 = public transports (train, bus, tram, etc.); 1 = private modes (car, motorbike, etc.); 2 = soft modes (bike, walk, etc.)

  • Gender:性别

 

 

image-20220930091519387

  • public transportation 和slow modes 被分为了一个nest

  • car单独组成一个nest

    不同主枝下面可以有不同数量的分支,也可以只有一个分支

     

 

定义效用函数

三种交通方式的效用模型。这个效用函数,实质上 = 常数+时间成本+经济成本

image-20220930132302195

image-20220930130955627

  • ASC CAR, ASC SM, BETA TIME FULLTIME, BETA TIME OTHER, BETA DIST MALE, BETA DIST FEMALE, BETA DIST UNREPORTED, BETA COST, are parameters to be estimated(是我们需要进行标定的参数)

  • TimePT scale, MarginalCostPT scaled, TimeCar scale, CostCarCHF scale, distance km scale are attributes and fulltime , notfulltime, male, female, unreportedGender are socio-economic characteristics. (其他的是社会经济变量:时间、成本、距离、性别)

  • ASC CAR, ASC SM,ASC_PT :常数项

    需要注意的是:并不是所以的常数项都可以标定得到。(参考教程

 


代码

(整体的代码结构参照:夜间清风-交通出行预测

整体代码见 biogeme -02estimation.ipynb

1.导入库+读取数据库文件并创建数据库对象

 

 # 导入库
 from biogeme import models
 import biogeme.biogeme as bio
 from scenarios import scenario, database
 import pandas as pd
 import biogeme.database as db
 from biogeme.expressions import Beta, Variable
 
 #读取数据
 df = pd.read_csv('optima.dat', sep='\t')  #dataframe
 #创建数据库对象
 database = db.Database('optima', df)

 

2. 根据数据设置变量

数据库中的每一个列名对应的模型中的一个变量。用以下语句设置变量:

 Choice = Variable('Choice')
 TimePT = Variable('TimePT')
 TimeCar = Variable('TimeCar')
 MarginalCostPT = Variable('MarginalCostPT')    # 边际成本
 CostCarCHF = Variable('CostCarCHF')  #
 distance_km = Variable('distance_km')   #距离
 Gender = Variable('Gender')
 OccupStat = Variable('OccupStat')
 Weight = Variable('Weight')

 

3. 设置待估计变量

 # List of parameters to be estimated   #待估计的参数列表
 ASC_CAR = Beta('ASC_CAR', 0, None, None, 0)  #Beta:从数据中估计未知参数
 ASC_PT = Beta('ASC_PT', 0, None, None, 1)
 ASC_SM = Beta('ASC_SM', 0, None, None, 0)
 BETA_TIME_FULLTIME = Beta('BETA_TIME_FULLTIME', 0, None, None, 0)
 BETA_TIME_OTHER = Beta('BETA_TIME_OTHER', 0, None, None, 0)
 BETA_DIST_MALE = Beta('BETA_DIST_MALE', 0, None, None, 0)
 BETA_DIST_FEMALE = Beta('BETA_DIST_FEMALE', 0, None, None, 0)
 BETA_DIST_UNREPORTED = Beta('BETA_DIST_UNREPORTED', 0, None, None, 0)
 BETA_COST = Beta('BETA_COST', 0, None, None, 0)

Biogeme用Beta库来定义待估参数(模型系数)。ASC CAR, ASC SM,ASC_PT :常数项。其他参数是各变量对应的系数。仔细看一下效用函数就能明白。

在定义待估计系数时,需要提供5个信息:

  • 参数名称

  • 参数的默认值

  • 参数下限,没有用None表示

  • 参数上限,没有用None表示

  • 还要给参数一个标记:0表示要估计该参数;1表示保持它的默认值

4. 变量数值缩放

由于数值原因,最好将数据缩放到1左右

 TimePT_scaled = TimePT / 200
 TimeCar_scaled = TimeCar / 200
 CostCarCHF_scaled = CostCarCHF / 10
 distance_km_scaled = distance_km / 5
 male = Gender == 1
 female = Gender == 2
 unreportedGender = Gender == -1
 
 fulltime = OccupStat == 1
 notfulltime = OccupStat != 1
 #OccupStat 1:表示全职 !1 :非全职

 

5. 定义效用函数

 MarginalCostScenario = MarginalCostPT * factor
 MarginalCostPT_scaled = MarginalCostScenario / 10
 #定义效用函数
 V_PT = (
     ASC_PT
     + BETA_TIME_FULLTIME * TimePT_scaled * fulltime
     + BETA_TIME_OTHER * TimePT_scaled * notfulltime
     + BETA_COST * MarginalCostPT_scaled
 )
 V_CAR = (
     ASC_CAR
     + BETA_TIME_FULLTIME * TimeCar_scaled * fulltime
     + BETA_TIME_OTHER * TimeCar_scaled * notfulltime
     + BETA_COST * CostCarCHF_scaled
 )
 V_SM = (
     ASC_SM
     + BETA_DIST_MALE * distance_km_scaled * male
     + BETA_DIST_FEMALE * distance_km_scaled * female
     + BETA_DIST_UNREPORTED * distance_km_scaled * unreportedGender
 )

 

6. 对应效用函数与方案

 V = {
  0: V_PT,
    1: V_CAR,
    2: V_SM
    }

V: 代表各备选方案的效用函数的字典

7. 定义nest结构(分层结构)

 MU_NOCAR = Beta('MU_NOCAR', 1.0, 1.0, None, 0)
 
 CAR_NEST = 1.0, [1]  # car单独分为一类 nest:CAR_NEST
 NO_CAR_NEST = MU_NOCAR, [0, 2]  # PT、SM分为一类 nest:NO_CAR_NEST
 nests = CAR_NEST, NO_CAR_NEST

 

8. 计算logprob(概率的对数)

 logprob = models.lognested(V, None, nests, Choice)  #V是效用函数 
 #需要输入1. 各备选方案的效用函数字典 2. 各方案可用性字典(可以为None)
 #3. nests:tuple
 """
    Example::
        nesta = MUA, [1, 2, 3]
        nestb = MUB, [4, 5, 6]
        nests = nesta, nestb
 """
 #logprob:基于MEV,对嵌套logit模型的选择概率取log

9. estimation

# Create the Biogeme object for estimation  #估计
biogeme = bio.BIOGEME(database, logprob)
biogeme.modelName = '02estimation'

print('Estimation...')
# Estimate the parameters. Perform bootstrapping.
#估计参数
results = biogeme.estimate(bootstrap=100)

# Get the results in a pandas table
pandasResults = results.getEstimatedParameters()
print(pandasResults)

10. Simulation

print('Simulation...')

simulated_choices = logprob.getValue_c(
betas=results.getBetaValues(), database=database
)
#getValue_c ,如果提供了一整个数据库,则返回的是一个list,是将表达式应用到每一行的一个结果(概率值,浮点数)
#getValue_c()函数中需要输入上面标定的参数+数据库

print(simulated_choices)

loglikelihood = logprob.getValue_c(
betas=results.getBetaValues(), #标定系数值
database=database, #数据库
aggregation=True, #如果一共了一个database,且参数为true,则表达式应用于每个条目(行),并聚合所有值,返回和。
)

print(f'Final log likelihood: {results.data.logLike}')
print(f'Simulated log likelihood: {loglikelihood}')

结果

image-20220930151342493

 

整体代码

#%%
"""File scenarios.py

:author: Michel Bierlaire, EPFL
:date: Sun Oct 31 09:40:59 2021

Specification of a nested logit model, that will be estimated, and
used for simulation. Three alternatives: public transporation, car
and slow modes. RP data.
三种选择方式:public car 慢速模式
Based on the Optima data.
It contains a function that generates scenarios where the current
cost of public transportation is multiplied by a factor.
包含一个函数,生成当前公共交通成本*一个因素

"""

import pandas as pd
import biogeme.database as db
from biogeme.expressions import Beta, Variable
from biogeme import models
import biogeme.biogeme as bio
from scenarios import scenario, database


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

# Variables from the data #设置变量
Choice = Variable('Choice')
TimePT = Variable('TimePT')
TimeCar = Variable('TimeCar')
MarginalCostPT = Variable('MarginalCostPT') # 边际成本
CostCarCHF = Variable('CostCarCHF') #
distance_km = Variable('distance_km') #距离
Gender = Variable('Gender')
OccupStat = Variable('OccupStat')
Weight = Variable('Weight')

# Exclude observations such that the chosen alternative is -1
database.remove(Choice == -1.0)

# Normalize the weights (规范化权重)
sumWeight = database.data['Weight'].sum()
numberOfRows = database.data.shape[0] #行数
normalizedWeight = Weight * numberOfRows / sumWeight

# List of parameters to be estimated #待估计的参数列表
ASC_CAR = Beta('ASC_CAR', 0, None, None, 0) #Beta:从数据中估计未知参数
ASC_PT = Beta('ASC_PT', 0, None, None, 1)
ASC_SM = Beta('ASC_SM', 0, None, None, 0)
BETA_TIME_FULLTIME = Beta('BETA_TIME_FULLTIME', 0, None, None, 0)
BETA_TIME_OTHER = Beta('BETA_TIME_OTHER', 0, None, None, 0)
BETA_DIST_MALE = Beta('BETA_DIST_MALE', 0, None, None, 0)
BETA_DIST_FEMALE = Beta('BETA_DIST_FEMALE', 0, None, None, 0)
BETA_DIST_UNREPORTED = Beta('BETA_DIST_UNREPORTED', 0, None, None, 0)
BETA_COST = Beta('BETA_COST', 0, None, None, 0)

# Definition of variables: 定义变量
# For numerical reasons, it is good practice to scale the data to
# that the values of the parameters are around 1.0.
# 由于数值原因,最好将数据缩放到1左右

TimePT_scaled = TimePT / 200
TimeCar_scaled = TimeCar / 200
CostCarCHF_scaled = CostCarCHF / 10
distance_km_scaled = distance_km / 5
male = Gender == 1
female = Gender == 2
unreportedGender = Gender == -1

fulltime = OccupStat == 1
notfulltime = OccupStat != 1
#OccupStat 1:表示全职 !1 :非全职

# def scenario(factor=1.0):
"""Provide the model specification for a scenario with the price of
public transportation is multiplied by a factor

:param factor: factor that multiples the price of public transportation.
:type factor: float

:return: a dict with the utility functions, the nesting structure,
and the choice expression.
返回的是一个字典,包含效用函数,嵌套结构,选择表达式

:rtype: dict(int: biogeme.expression), tuple(biogeme.expression,
list(int)), biogeme.expression
"""
factor = 1.0
MarginalCostScenario = MarginalCostPT * factor
MarginalCostPT_scaled = MarginalCostScenario / 10
# Definition of utility functions:
V_PT = (
ASC_PT
+ BETA_TIME_FULLTIME * TimePT_scaled * fulltime
+ BETA_TIME_OTHER * TimePT_scaled * notfulltime
+ BETA_COST * MarginalCostPT_scaled
)
V_CAR = (
ASC_CAR
+ BETA_TIME_FULLTIME * TimeCar_scaled * fulltime
+ BETA_TIME_OTHER * TimeCar_scaled * notfulltime
+ BETA_COST * CostCarCHF_scaled
)
V_SM = (
ASC_SM
+ BETA_DIST_MALE * distance_km_scaled * male
+ BETA_DIST_FEMALE * distance_km_scaled * female
+ BETA_DIST_UNREPORTED * distance_km_scaled * unreportedGender
)

# Associate utility functions with the numbering of alte1rnatives
#将效用函数和方案编号联系起来
V = {0: V_PT, 1: V_CAR, 2: V_SM}

# Definition of the nests:
# 1: nests parameter 嵌套系数
# 2: list of alternatives #备选列表
MU_NOCAR = Beta('MU_NOCAR', 1.0, 1.0, None, 0)

CAR_NEST = 1.0, [1] # car单独分为一类 nest:CAR_NEST
NO_CAR_NEST = MU_NOCAR, [0, 2] # PT、SM分为一类 nest:NO_CAR_NEST
nests = CAR_NEST, NO_CAR_NEST
# return V, nests, Choice, MarginalCostScenario
# The choice model is a nested logit, with availability conditions
# For estimation, we need the log of the probability
# log(prob) :概率的对数
logprob = models.lognested(V, None, nests, Choice) #V是效用函数

# Create the Biogeme object for estimation #估计
biogeme = bio.BIOGEME(database, logprob)
biogeme.modelName = '02estimation'

print('Estimation...')
# Estimate the parameters. Perform bootstrapping.
results = biogeme.estimate(bootstrap=100)

# Get the results in a pandas table
pandasResults = results.getEstimatedParameters()
print(pandasResults)

print('Simulation...')

simulated_choices = logprob.getValue_c(
betas=results.getBetaValues(), database=database
)

print(simulated_choices)

loglikelihood = logprob.getValue_c(
betas=results.getBetaValues(),
database=database,
aggregation=True,
)

print(f'Final log likelihood: {results.data.logLike}')
print(f'Simulated log likelihood: {loglikelihood}')

 

附:变量说明

image-20220930125231560

image-20220930125301937

image-20220930125430334

 

posted @ 2022-09-30 15:24  sheyueyu  阅读(824)  评论(2编辑  收藏  举报