biogeme巢式logit选择学习

导入库

import pandas as pd
import csv
from biogeme import models
import biogeme.biogeme as bio
import biogeme.database as db
from biogeme.expressions import Beta, Variable
import biogeme.results as res
import biogeme.exceptions as excep
import sys

输入选择集

df=pd.read_csv('choice.csv')
database=db.Database('choice',df)##输入选择集创建数据库


globals().update(database.variables)##定义全局变量

database.variables

{'Unnamed: 0.2': Unnamed: 0.2,
'Unnamed: 0.1': Unnamed: 0.1,
'Unnamed: 0': Unnamed: 0,
'id': id,
'hc_time_all': hc_time_all,
'hc_cishu': hc_cishu,
'hc_fare': hc_fare,
'path_length': path_length,
'car_fare': car_fare,
'car_time': car_time,
'car_cishu': car_cishu,
'bus_time': bus_time,
'metro_fare': metro_fare,
'metro_cishu': metro_cishu,
'metro_time': metro_time,
'bus_fare': bus_fare,
'bus_cishu': bus_cishu,
'choice': choice}

定义变量

#首先定义常数项
ASC_CAR = Beta('ASC_CAR',0,None,None,0)
ASC_BUS = Beta('ASC_BUS',0,None,None,0)
ASC_METRO = Beta('ASC_METRO',0,None,None,0)
ASC_HC = Beta('ASC_HC',0,None,None,0)
#定义与效用有关参数的系数
BETA_TIME_CAR = Beta('BETA_TIME_CAR',0,None,None,0)
BETA_FARE_CAR = Beta('BETA_FARE_CAR',0,None,None,0)
BETA_LENGTH = Beta('BETA_LENGTH',0,None,None,0)
BETA_TIME_BUS =Beta('BETA_TIME_BUS',0,None,None,0)
BETA_FARE_BUS = Beta('BETA_FARE_BUS',0,None,None,0)
BEAT_CISHU_BUS = Beta('BETA_CISHU_BUS',0,None,None,0)
BETA_TIME_METRO = Beta('BETA_TIME_METRO',0,None,None,0)
BETA_FARE_METRO = Beta('BETA_FARE_METRO',0,None,None,0)
BETA_CISHU_METRO = Beta('BETA_CISHU_METRO',0,None,None,0)
BETA_TIME_HC = Beta('BETA_TIME_HC',0,None,None,0)
BETA_FARE_HC = Beta('BETA_FARE_HC',0,None,None,0)
BETA_CISHU_HC = Beta('BETA_CISHU_HC',0,None,None,0)

定义效用函数

V_CAR = (
ASC_CAR
+BETA_TIME_CAR*car_time
+BETA_FARE_CAR*car_fare
+BETA_LENGTH*path_length
)
V_BUS = (
ASC_BUS
+BETA_TIME_BUS*bus_time
+BETA_FARE_BUS*bus_fare
+BEAT_CISHU_BUS*bus_cishu
)
V_METRO = (
ASC_METRO
+BETA_TIME_METRO*metro_time
+BETA_FARE_METRO*metro_fare
+BETA_CISHU_METRO*metro_cishu
)
V_HC = (
ASC_HC
+BETA_TIME_HC*hc_time_all
+BETA_FARE_HC*hc_fare
+BETA_CISHU_METRO*hc_cishu
)

效用对应方案

V = {
0:V_CAR,
1:V_BUS,
2:V_METRO,
3:V_HC
}

定义树形结构

PUBLIC = Beta('PUBLIC',1,1,None,0)

CAR_NEST = 1.0,[0]
PUBLIC_NEST = PUBLIC,[1,2,3]
nests = CAR_NEST,PUBLIC_NEST

计算logprob

logprob = models.lognested(V,None,nests,choice)

估计参数

biogeme = bio.BIOGEME(database, logprob)

biogeme.modelName = 'mode_choice_estimation'
results=biogeme.estimate()
Results=results.getEstimatedParameters()
print(Results)

[16:12:00] < Warning > Numerical problem with the second derivative matrix. Norm = inf. Replaced by the identity matrix.
[16:12:01] < Warning > Numerical problem with the second derivative matrix. Norm = inf. Replaced by the identity matrix.
[16:12:01] < Warning > Numerical problem with the second derivative matrix. Norm = inf. Replaced by the identity matrix.
[16:12:02] < Warning > Numerical problem with the second derivative matrix. Norm = inf. Replaced by the identity matrix.
[16:12:03] < Warning > Numerical problem with the second derivative matrix. Norm = inf. Replaced by the identity matrix.
[16:12:04] < Warning > Numerical problem with the second derivative matrix. Norm = inf. Replaced by the identity matrix.
[16:12:06] < Warning > Numerical problem with the second derivative matrix. Norm = inf. Replaced by the identity matrix.
[16:12:10] < Warning > Numerical problem with the second derivative matrix. Norm = inf. Replaced by the identity matrix.
[16:12:12] < Warning > Numerical problem with the second derivative matrix. Norm = inf. Replaced by the identity matrix.
[16:12:15] < Warning > Numerical problem with the second derivative matrix. Norm = inf. Replaced by the identity matrix.
[16:12:17] < Warning > Numerical problem with the second derivative matrix. Norm = inf. Replaced by the identity matrix.
[16:12:18] < Warning > Numerical problem with the second derivative matrix. Norm = inf. Replaced by the identity matrix.
[16:12:20] < Warning > Numerical problem with the second derivative matrix. Norm = inf. Replaced by the identity matrix.
[16:12:23] < Warning > Numerical problem with the second derivative matrix. Norm = inf. Replaced by the identity matrix.
[16:12:24] < Warning > Numerical problem with the second derivative matrix. Norm = inf. Replaced by the identity matrix.
[16:12:25] < Warning > Numerical problem with the second derivative matrix. Norm = inf. Replaced by the identity matrix.
[16:12:26] < Warning > Numerical problem with the second derivative matrix. Norm = inf. Replaced by the identity matrix.
[16:12:27] < Warning > Numerical problem with the second derivative matrix. Norm = inf. Replaced by the identity matrix.
[16:12:28] < Warning > Numerical problem with the second derivative matrix. Norm = inf. Replaced by the identity matrix.
[16:12:29] < Warning > Numerical problem with the second derivative matrix. Norm = inf. Replaced by the identity matrix.
[16:12:30] < Warning > Numerical problem with the second derivative matrix. Norm = inf. Replaced by the identity matrix.
[16:12:31] < Warning > Numerical problem with the second derivative matrix. Norm = inf. Replaced by the identity matrix.
[16:12:31] < Warning > Numerical problem with the second derivative matrix. Norm = inf. Replaced by the identity matrix.
[16:12:32] < Warning > Numerical problem with the second derivative matrix. Norm = inf. Replaced by the identity matrix.
[16:12:33] < Warning > Numerical problem with the second derivative matrix. Norm = inf. Replaced by the identity matrix.
Value Rob. Std err Rob. t-test Rob. p-value
ASC_BUS -7.965898e-01 0.0 1.797693e+308 0.0
ASC_CAR -1.250052e+00 0.0 1.797693e+308 0.0
ASC_HC 9.795643e-01 0.0 1.797693e+308 0.0
ASC_METRO 1.067077e+00 0.0 1.797693e+308 0.0
BETA_CISHU_BUS -1.776695e+00 0.0 1.797693e+308 0.0
BETA_CISHU_METRO -6.685987e-03 0.0 1.797693e+308 0.0
BETA_FARE_BUS 8.880739e-01 0.0 1.797693e+308 0.0
BETA_FARE_CAR 2.888911e-07 0.0 1.797693e+308 0.0
BETA_FARE_HC 3.266817e-04 0.0 1.797693e+308 0.0
BETA_FARE_METRO -3.187853e-02 0.0 1.797693e+308 0.0
BETA_LENGTH 1.094267e-04 0.0 1.797693e+308 0.0
BETA_TIME_BUS -8.687563e-07 0.0 1.797693e+308 0.0
BETA_TIME_CAR -8.526598e-04 0.0 1.797693e+308 0.0
BETA_TIME_HC 4.400077e-06 0.0 1.797693e+308 0.0
BETA_TIME_METRO -3.083993e-07 0.0 1.797693e+308 0.0
PUBLIC 6.521609e+02 0.0 1.797693e+308 0.0

获取每个备选方案的选择概率

prob_CAR = models.nested(V, None, nests, 0)
prob_BUS = models.nested(V, None, nests, 1)
prob_METRO = models.nested(V, None, nests, 2)
prob_HC = models.nested(V, None, nests, 3)

读取估计结果

# Read the estimation results from the file
try:
results = res.bioResults(pickleFile='mode_choice_estimation.pickle')
except excep.biogemeError:
sys.exit(
'Run first the script 02simulation.py '
'in order to generate the '
'file 02estimation.pickle.'
)

进行模拟，计算选择概率和效用值

simulate_formulas = {
'Utility CAR': V_CAR.getValue_c(
betas=results.getBetaValues(), database=database, prepareIds=True
),
'Utility BUS': V_BUS.getValue_c(
betas=results.getBetaValues(), database=database, prepareIds=True
),
'Utility METRO': V_METRO.getValue_c(
betas=results.getBetaValues(), database=database, prepareIds=True
),
'Utility HC': V_HC.getValue_c(
betas=results.getBetaValues(), database=database, prepareIds=True
),
'Prob. CAR': prob_CAR.getValue_c(
betas=results.getBetaValues(), database=database, prepareIds=True
),
'Prob. BUS': prob_BUS.getValue_c(
betas=results.getBetaValues(), database=database, prepareIds=True
),
'Prob. METRO': prob_METRO.getValue_c(
betas=results.getBetaValues(), database=database, prepareIds=True
),
'Prob. HC': prob_HC.getValue_c(
betas=results.getBetaValues(), database=database, prepareIds=True
),
}
simulated_values = pd.DataFrame.from_dict(
simulate_formulas,
)


simulated_values
Utility CAR Utility BUS Utility METRO Utility HC Prob. CAR Prob. BUS Prob. METRO Prob. HC
0 -1.934563 1.258228 1.318423 1.282484 0.037193 0.006809 0.907109 0.048888
1 -1.934563 1.258228 1.318434 1.282424 0.037193 0.006805 0.907380 0.048621
2 -1.934563 1.274554 1.306116 1.282484 0.037573 0.060507 0.786654 0.115265
3 -1.934563 1.274554 1.306116 1.282424 0.037574 0.060543 0.787116 0.114767
4 -1.934563 1.274554 1.306092 1.282484 0.037574 0.060606 0.786367 0.115453
... ... ... ... ... ... ... ... ...
2425 -1.909643 1.258093 1.309319 1.289848 0.038368 0.012255 0.787545 0.161833
2426 -1.909643 1.257291 1.308803 1.289848 0.038384 0.011891 0.782124 0.167601
2427 -1.909643 1.274355 1.309319 1.289848 0.038352 0.044393 0.760898 0.156357
2428 -1.909643 1.266182 1.308803 1.289848 0.038378 0.024174 0.772013 0.165435
2429 -1.909643 1.266453 1.309319 1.289848 0.038362 0.023878 0.777908 0.159852

2430 rows × 8 columns

另一种模拟方法，加入biogeme

模拟结果相同，但速度更快

simulate = {
'Utility CAR': V_CAR,
'Utility BUS': V_BUS,
'Utility METRO': V_METRO,
'Utility HC': V_HC,
'Prob. CAR': prob_CAR,
'Prob. BUS': prob_BUS,
'Prob. METRO': prob_METRO,
'Prob. HC': prob_HC,
}

biogeme = bio.BIOGEME(database, simulate)
biogeme_simulation = biogeme.simulate(results.getBetaValues())

biogeme_simulation
Utility CAR Utility BUS Utility METRO Utility HC Prob. CAR Prob. BUS Prob. METRO Prob. HC
0 -1.934563 1.258228 1.318423 1.282484 0.037193 0.006809 0.907109 0.048888
1 -1.934563 1.258228 1.318434 1.282424 0.037193 0.006805 0.907380 0.048621
2 -1.934563 1.274554 1.306116 1.282484 0.037573 0.060507 0.786654 0.115265
3 -1.934563 1.274554 1.306116 1.282424 0.037574 0.060543 0.787116 0.114767
4 -1.934563 1.274554 1.306092 1.282484 0.037574 0.060606 0.786367 0.115453
... ... ... ... ... ... ... ... ...
2425 -1.909643 1.258093 1.309319 1.289848 0.038368 0.012255 0.787545 0.161833
2426 -1.909643 1.257291 1.308803 1.289848 0.038384 0.011891 0.782124 0.167601
2427 -1.909643 1.274355 1.309319 1.289848 0.038352 0.044393 0.760898 0.156357
2428 -1.909643 1.266182 1.308803 1.289848 0.038378 0.024174 0.772013 0.165435
2429 -1.909643 1.266453 1.309319 1.289848 0.038362 0.023878 0.777908 0.159852

2430 rows × 8 columns

posted @ 2022-11-01 19:03  权BB  阅读(307)  评论(0)    收藏  举报