官术网_书友最值得收藏!

首頁 >

數據科學工程實踐：用戶行為分析與建模、A/B實驗、SQLFlow

1.3　DCM模型的Python實踐

我們在1.1節、1.2節系統地學習了選擇行為的經濟學理論、DCM的設計原理以及相關數學知識，本節將基于Python對旅行出行方式選擇案例進行模型搭建及解讀，希望讀者能快速行動起來，將理論結合實際，盡快學以致用。

1.3.1　軟件包和數據格式

DCM在Python中有可以直接調用的軟件包，其中LR模型可以使用statsmodels軟件包，MNL模型、NL模型可以使用pylogit軟件包。這里推薦讀者先安裝好Anaconda（一款非常流行的數據分析平臺），安裝完畢后，再手動安裝statsmodels、pylogit軟件包，直接在終端輸入代碼清單1-1所示的命令即可。

代碼清單1-1　安裝相關軟件包

pip install statsmodels
pip install pylogit

有過建模經驗的讀者都知道，提供模型訓練/預測的數據要嚴格符合模型要求，否則模型會運行失敗。這一點對DCM尤為重要，這里先用一定的篇幅著重介紹DCM常用的兩種數據格式：寬格式和長格式。

本節的案例數據來自William Greene《微觀經濟學建模及離散選擇分析》課程中的旅行模式選擇數據（Travel Mode Choice Data）。為了方便使用，我們先對數據進行一定的處理，處理后的數據形式如表1-2所示。

表1-2　旅行模式選擇數據的字典

000

大多數讀者應該對寬格式比較熟悉，因為常用的機器學習模型的輸入數據大多是寬格式。如表1-3所示，一行數據代表一次選擇過程，字段表示選擇結果，其他與選擇有關的信息被平行放置在一行中的各字段。LR模型的輸入數據格式就是寬格式。

表1-3　寬格式示例

000

對于長格式，讀者可能就比較陌生了，不過它是多數Logit模型使用的數據格式。如表1-4所示，一次選擇行為包含多行數據，一行數據代表一次選擇過程中的一個選項，其中會有一個字段表示該選項是否被選中（選中為1、未選中為0，一次選擇只有一個選項被選中），對于備選項屬性數據可以不同，決策者屬性數據則相同。

表1-4　長格式示例

000

pylogit提供了長/寬數據格式相互轉換的函數，操作實例如代碼清單1-2所示。

代碼清單1-2　使用pylogit進行長/寬數據格式相互轉換

from collections import OrderedDict  # OrderedDict用于記錄模型的specification（聲明）
import pylogit as pl                 # 引入Logit模型軟件包pylogit
# 數據讀入
long_data_path = u'long_data.csv'
long_df.to_csv(long_data_path, sep=',',index = False)
#---------------------#
# “長格式”轉換為“寬格式”#
#---------------------#
# 指定決策者屬性的列表
individual_specific_variables = ["HINC","PSIZE"]
# 指定備選項屬性的列表
alternative_specific_variables = ['TTME', 'INVC', 'INVT', 'GC']
# 指定備選項屬性的特殊說明列表
subset_specific_variables = {}
# “觀測ID”，標識每次選擇
observation_id_column = "OBS_ID"
# “備選項ID”，標識備選方案
alternative_id_column = "ALT_ID"
# “選擇結果”，標識選擇結果
choice_column = "MODE"
# 可選變量，記錄與每個備選方案對應的名稱，并允許在寬格式數據中創建有意義的列名
alternative_name_dict = {0: "AIR",
                        1: "TRAIN",
                        2: "BUS",
                        3: "CAR"}
wide_df = pl.convert_long_to_wide(long_df,
                               individual_specific_variables,
                               alternative_specific_variables,
                               subset_specific_variables,
                               observation_id_column,
                               alternative_id_column,
                               choice_column,
                               alternative_name_dict)
#---------------------#
# “寬格式”轉換為“長格式”#
#---------------------#
# 創建決策者變量的列表
ind_variables =["HINC","PSIZE"]
# 指定每個備選項的屬性所對應的字段，0/1/2/3代表備選項，后面的字段名為其屬性在“寬格式數據”中的列名
alt_varying_variables = {"TTME": dict([(0, 'TTME_AIR'),
                                    (1, 'TTME_TRAIN'),
                                    (2, 'TTME_BUS'),
                                    (3, 'TTME_CAR')]),
                        "INVC": dict([(0, 'INVC_AIR'),
                                    (1, 'INVC_TRAIN'),
                                    (2, 'INVC_BUS'),
                                    (3, 'INVC_CAR')]),
                        "INVT": dict([(0, 'INVT_AIR'),
                                    (1, 'INVT_TRAIN'),
                                    (2, 'INVT_BUS'),
                                    (3, 'INVT_CAR')]),
                         "GC":   dict([(0, 'GC_AIR'),
                                       (1, 'GC_TRAIN'),
                                       (2, 'GC_BUS'),
                                       (3, 'GC_CAR')])}
# 指定可用性變量，字典的鍵為可選項的ID，值為數據集中標記可用性的列
# 由于篇幅所限，前面的寬格式數據中省略了這3列，實際在做數據轉化時需要標識每個選項的可用性
availability_variables = {0: 'availability_AIR',
                         1: 'availability_TRAIN',
                         2: 'availability_BUS',
                         3: 'availability_CAR'}
# “備選項ID”標識與每一行相關聯的備選方案
custom_alt_id = "ALT_ID"
# “觀測ID”標識每次選擇，觀測ID需要從1開始
obs_id_column = "OBS_ID"
wide_df[obs_id_column] = np.arange(wide_df.shape[0], dtype=int) + 1
# “選擇結果”標識選擇結果
choice_column = 'MODE'
# 執行長格式轉換
long_df = pl.convert_wide_to_long(wide_df,
                            ind_variables,
                               alt_varying_variables,
                               availability_variables,
                               obs_id_column,
                               choice_column,
                               new_alt_id_name=custom_alt_id)

1.3.2　使用邏輯回歸分析自駕選擇問題

基于前文的介紹，相信讀者已經迫不及待想使用MNL模型或NL模型進行建模分析了，這里先從LR模型的實操講起。LR模型是目前應用最廣泛的可解釋二分類模型之一，深入了解LR模型對我們的日常工作有很大幫助。

通過對案例數據進行一定的處理，可以得到一份滿足LR模型要求的寬格式數據。具體數據描述如表1-5所示，場景邏輯如圖1-5所示。

表1-5　LR模型訓練數據的字典

000
000

圖1-5　LR模型的場景邏輯示意圖

了解數據形式后，開始進行具體的模型搭建工作。

第1步：引入軟件包，讀取數據。重要的軟件包在代碼的備注中，如代碼清單1-3所示。

代碼清單1-3　引入軟件包及讀取數據

import numpy as np                     # 引入基礎軟件包numpy
import pandas as pd                    # 引入基礎軟件包pandas
import statsmodels.api as sm     # 引入Logistic regression軟件包statsmodels
from sklearn.model_selection import train_test_split # 引入訓練集/測試集構造工具包
from sklearn import metrics            # 引入模型評價指標AUC計算工具包
import matplotlib.pyplot as plt        # 引入繪圖軟件包
import scipy                           # 引入scipy軟件包完成卡方檢驗
# 數據讀入
data_path = 'wide_data.csv'
raw_data = pd.read_table(data_path, sep=',', header=0)

第2步：數據預處理。數據預處理工作對于任何模型搭建都是必要的，這里結合LR模型及后續將介紹的MNL模型、NL模型的特點著重講3個數據預處理的要點：①不要存在缺失值；②每一列數據均為數值型；③多枚舉值離散變量輸入模型前要進行啞變量處理，如代碼清單1-4所示。

代碼清單1-4　數據預處理

# 1. 缺失值探查&簡單處理
model_data.info()                      # 查看每一列的數據類型和數值缺失情況
# | RangeIndex: 210 entries, 0 to 209
# | Data columns (total 9 columns):
# | ...
# | HINC              210 non-null int64
# | ...
model_data = model_data.dropna()       # 缺失值處理—刪除
model_data = model_data.fillna(0)      # 缺失值處理—填充（零、均值、中位數、預測值等）

# 2. 數值型核查（連續變量應為int64或float數據類型）
# 若上一步中存在應為連續數值變量的字段為object，則執行下列代碼，這里假設'HINC'存在為字符串'null'的值
import re                               # 正則表達式工具包
float_patten = '^(-?\\d+)(\\.\\d+)?$'   # 定義浮點數正則patten
float_re = re.compile(float_patten)     # 編譯
model_data['HINC'][model_data['HINC'].apply(lambda x : 'not_float' if float_re.match(str(x)) == None else 'float') == 'not_float'] # 查看非浮點型數據
# | 2    null
# | Name: distance, dtype: object
model_data = model_data[model_data['HINC'] != 'null']
model_data['HINC'] = model_data['HINC'].astype(float)

第3步：單變量分析。在建模之前需要對每個自變量進行單變量分析，確定是否納入模型。變量分為離散變量和連續變量兩種，分析方式也有所不同。對于離散變量，我們使用k-1自由度的卡方檢驗，其中k為離散變量的值個數；對于連續變量，比較簡單的分析方法是直接對單變量進行邏輯回歸，查看回歸系數的顯著性，根據AUC分析自變量對y的解釋能力。保留顯著的自變量進入后續的操作，如代碼清單1-5所示。

代碼清單1-5　單變量分析

# 離散變量分析
crosstab = pd.crosstab( model_data['y'],model_data['PSIZE'])
p=scipy.stats.chi2_contingency(crosstab)[1]
print("PSIZE:",p)
# PSIZE: 0.0024577358937625327

# 連續變量分析
logistic = sm.Logit(model_data['y'],model_data['INVT_CAR']).fit()
p = logistic.pvalues['INVT_CAR']
y_predict = logistic.predict(model_data['INVT_CAR'])
AUC = metrics.roc_auc_score(model_data['y'],y_predict)
result = 'INVT_CAR:'+str(p)+'  AUC:'+str(AUC)
print(result)
# INVT_CAR:2.971604856310474e-09  AUC:0.6242563699629587

第4步：共線性檢驗。由于LR模型是一種廣義線性模型，變量間嚴重的共線性會對參數估計的準確性及泛化能力產生影響，因此需要對自變量間的共線性進行分析。若vif值大于10，可認為變量間具有很強的共線性，需要進行相應的處理，最簡單的處理方式就是剔除自變量，保留單變量分析中AUC最大的變量。共線性檢驗示例如代碼清單1-6所示。

代碼清單1-6　共線性檢驗

from statsmodels.stats.outliers_influence import variance_inflation_factor 
#共線性診斷包
X = raw_data[[ 'INVT_AIR', 'INVT_TRAIN','INVT_BUS', 'INVT_CAR']]
vif = pd.DataFrame()
vif['VIF Factor'] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
vif['features'] = X.columns
print('================多重共線性==============')
print(vif)
# | 0   14.229424    INVT_AIR
# | 1   72.782420  INVT_TRAIN
# | 2   80.279742    INVT_BUS
# | 3   35.003438    INVT_CAR

第5步：模型搭建。這里需要注意的是，對于3值及以上的離散變量要進行啞變量處理（需要記住去掉的枚舉值），并且增加截距項Intercept，同時進行訓練集和測試集的拆分（目的是防止模型過擬合，確定分析結論可以泛化），代碼如清單1-7所示。

代碼清單1-7　搭建LR模型

# 建模數據構造
X = model_data[[ 'HINC','PSIZE','TTME_TRAIN' , 'INVC_CAR']]
y = raw_data['y']
# 啞變量處理
dummies = pd.get_dummies(X['PSIZE'], drop_first=False)
dummies.columns = [ 'PSIZE'+'_'+str(x) for x in dummies.columns.values]
X = pd.concat([X, dummies], axis=1)
X = X.drop('PSIZE',axis=1)   # 刪去原離散變量
X = X.drop('PSIZE_4',axis=1) # 刪去過于稀疏的字段
X = X.drop('PSIZE_5',axis=1) # 刪去過于稀疏的字段
X = X.drop('PSIZE_6',axis=1) # 刪去過于稀疏的字段
X['Intercept'] = 1           # 增加截距項
# 訓練集與測試集的比例分別為80%和20%
X_train, X_test, y_train, y_test = train_test_split(X, y, train_size = 0.8, random_state=1234)
# 建模
logistic = sm.Logit(y_train,X_train).fit()
print(logistic.summary2())
# 重要返回信息
# | ------------------------------------------------------------------
# |                Coef.   Std.Err.     z     P>|z|    [0.025   0.975]
# | ------------------------------------------------------------------
# | HINC           0.0264    0.0100   2.6477  0.0081   0.0068   0.0459
# | TTME_TRAIN     0.0389    0.0195   1.9916  0.0464   0.0006   0.0772
# | INVC_CAR      -0.0512    0.0204  -2.5103  0.0121  -0.0913  -0.0112
# | PSIZE_1       -0.3077    0.7317  -0.4206  0.6741  -1.7419   1.1264
# | PSIZE_2       -1.0800    0.6417  -1.6829  0.0924  -2.3378   0.1778
# | PSIZE_3       -0.7585    0.7582  -1.0004  0.3171  -2.2444   0.7275
# | Intercept     -1.8879    1.1138  -1.6951  0.0901  -4.0708   0.2950
# | =================================================================
# 模型評價
print("========訓練集AUC========")
y_train_predict = logistic.predict(X_train)
print(metrics.roc_auc_score(y_train,y_train_predict))
print("========測試集AUC========")
y_test_predict = logistic.predict(X_test)
print(metrics.roc_auc_score(y_test,y_test_predict))
# | ========訓練集AUC========
# | 0.7533854166666667
# | ========測試集AUC========
# | 0.6510263929618768

第6步：模型修正。可以看到，由于不顯著變量的影響，模型的測試集AUC與訓練集AUC存在較大差異，我們需要對不顯著變量進行剔除。可以看到，新建模型的擬合優度尚可（AUC接近0.75），且自變量顯著（p < 0.05），可以進行后續解讀，如代碼清單1-8所示。

代碼清單1-8　修正LR模型

X = X.drop('PSIZE_1',axis=1) 
X = X.drop('PSIZE_2',axis=1) 
X = X.drop('PSIZE_3',axis=1)
# 訓練集與測試集的比例分別為80%和20%
X_train, X_test, y_train, y_test = train_test_split(X, y, train_size = 0.8, random_state=1234)
# 建模
logistic = sm.Logit(y_train,X_train).fit()
print(logistic.summary2())
# 重要返回信息
# | ------------------------------------------------------------------
# |                Coef.   Std.Err.     z     P>|z|    [0.025   0.975]
# | ------------------------------------------------------------------
# | HINC           0.0266    0.0096   2.7731  0.0056   0.0078   0.0454
# | TTME_TRAIN     0.0335    0.0161   2.0838  0.0372   0.0020   0.0650
# | INVC_CAR      -0.0450    0.0168  -2.6805  0.0074  -0.0778  -0.0121
# | Intercept     -2.3486    0.8275  -2.8384  0.0045  -3.9704  -0.7269
# | =================================================================
print("========訓練集AUC========")
y_train_predict = logistic.predict(X_train)
print(metrics.roc_auc_score(y_train,y_train_predict))
print("========測試集AUC========")
y_test_predict = logistic.predict(X_test)
print(metrics.roc_auc_score(y_test,y_test_predict))
# | ========訓練集AUC========
# | 0.7344618055555555
# | ========測試集AUC========
# | 0.7419354838709677

第7步：模型解讀。DCM模型解讀的對象可以分為概率（probability）和幾率（odds）。在本例中，概率為“選擇自駕的概率”，幾率為“選擇自駕的概率/不選擇自駕的概率”。限于模型的數學性質，我們無法直接從模型參數中快速得到概率，而是需要經過一定的計算，這部分會在介紹復雜MNL模型及NL模型時展示。

得益于LR模型的數學性質，數據分析師可以基于模型參數直接對幾率進行解讀（這一點類似于線性回歸）。模型解讀的話術為“在其他條件保持不變的情況下，某因素增長一個單位（或屬性a相對屬性b），幾率會變化（增長或降低）多少”，計算公式如下。

連續變量：odd(xi + 1)/odd(xi)-1=exp(βi)-1

離散變量：odd(xj = 1)/odd(xj=0)-1=exp(βj)-1

例如，根據模型可知：

在其他條件保持不變的情況下，家庭收入增長1個單位，選擇自駕的odds會變化，exp(βHINC)-1=exp(0.0266)-1=0.027，即增加0.027倍。

在其他條件保持不變的情況下，自駕成本上升1個單位，選擇自駕的odds會變化，exp(βINVC_CAR)-1=exp(-0.0450)-1=-0.044，即下降0.044倍。

1.3.3　使用多項Logit模型分析多種交通方式選擇問題

1.3.2節我們使用邏輯回歸模型分析了是否選擇自駕的二項選擇問題，如果要同時分析4種交通方式的選擇問題，則需要使用MNL模型或NL模型。本節將介紹基于IIA假定的MNL模型，模型的問題場景映射如圖1-6所示。

000

圖1-6　MNL模型的場景邏輯示意圖

需要注意的是，MNL模型的輸入數據為長格式。不同于LR模型，MNL模型需要更加詳細、復雜的初始化聲明，以指定每種選項的效用函數形式。為了保證信息的完整性，盡量先保留自變量，定義如下模型。

Vair = ASCair + βttmeTTME + βinvcINVC + βinvtINVT + βhinc,airHINC + βpsize,airPSIZE

Vtrain = ASCtrain + βttmeTTME + βinvcINVC + βinvtINVT + βhinc,trainHINC + βpsize,trainPSIZE

Vbus = ASCbus + βttmeTTME + βinvcINVC + βinvtINVT + βhinc,busHINC + βpsize,busPSIZE

Vcar = βinvcINVC + βinvtINVT

根據設計好的模型結構搭建模型。1.3.2節已經介紹了數據清洗，這里不再贅述，直接進入模型搭建的代碼學習，如代碼清單1-9所示。

代碼清單1-9　搭建MNL模型

# 第一步：模型初始化聲明
basic_specification = OrderedDict()
basic_names = OrderedDict()
# 注意截距項包含選項個數減1 
basic_specification["intercept"] = [0, 1, 2]
basic_names["intercept"] = ['ASC_air', 'ASC_train', 'ASC_bus']
# 可以靈活指定備選項屬性的影響方式
basic_specification["TTME"] = [[0, 1, 2]]
basic_names["TTME"] = ['TTME']
basic_specification["INVC"] = [[0, 1, 2, 3]]
basic_names["INVC"] = ['INVC']
basic_specification["INVT"] = [[0, 1, 2, 3]]
basic_names["INVT"] = ['INVT']
# 也可以靈活指定決策者的影響方式，但需要注意的是，由于每個選項的決策者屬性都一樣，因此保證可估計性只對部分選項生效
basic_specification["HINC"] = [0, 1, 2]
basic_names["HINC"] = ['HINC_air', 'HINC_train', 'HINC_bus']
basic_specification["PSIZE"] = [0, 1, 2]
basic_names["PSIZE"] = ['PSIZE_air', 'PSIZE_train', 'PSIZE_bus']
# 第二步：創建模型
mnl = pl.create_choice_model(data = model_data,
                alt_id_col="ALT_ID",
                obs_id_col="OBS_ID",
                choice_col="MODE",
                specification=basic_specification,
                model_type = "MNL",
                names=basic_names)
# 第三步：模型估計和模型結果
mnl.fit_mle(np.zeros(12)) # 需要輸入模型參數數量，根據之前的模型表達式即可得到
mnl.get_statsmodels_summary()
# | -------------------------------------------------------------
# |               coef     std.err z       P>|z|   [0.025  0.975]
# | -------------------------------------------------------------
# | ASC_air       6.0352   1.138   5.302   0.000   3.804   8.266
# | ASC_train     5.5735   0.711   7.836   0.000   4.179   6.968
# | ASC_bus       4.5047   0.796   5.661   0.000   2.945   6.064
# | TTME         -0.1012   0.011  -9.081   0.000  -0.123  -0.079
# | INVC         -0.0087   0.008  -1.101   0.271  -0.024   0.007
# | INVT         -0.0041   0.001  -4.627   0.000  -0.006  -0.002
# | HINC_air      0.0075   0.013   0.567   0.571  -0.018   0.033
# | HINC_train   -0.0592   0.015  -3.977   0.000  -0.088  -0.03
# | HINC_bus     -0.0209   0.016  -1.278   0.201  -0.053   0.011
# | PSIZE_air    -0.9224   0.259  -3.568   0.000  -1.429  -0.416
# | PSIZE_train   0.2163   0.234   0.926   0.355  -0.242   0.674
# | PSIZE_bus    -0.1479   0.343  -0.432   0.666  -0.820   0.524
# |==============================================================

模型的搭建完成后，我們會發現有些變量不顯著，此時需要進行模型的修正，如代碼清單1-10所示。這里受篇幅限制，主要使用屬性剔除及屬性影響合并的方式進行修正，修正后的模型聲明及模型效果如下。

Vair = ASCair + βttmeTTME + βinvtINVT + βpsize,airPSIZE

Vtrain = ASCtrain + βttmeTTME + βinvtINVT + βhinc,trainHINC

Vbus = ASCbus + βttmeTTME + βinvtINVT + βhinc,busHINC

Vcar = βinvtINVT

代碼清單1-10　修正MNL模型

basic_specification = OrderedDict()
basic_names = OrderedDict()
basic_specification["intercept"] = [0, 1, 2]
basic_names["intercept"] = ['ASC_air', 'ASC_train', 'ASC_bus']
basic_specification["TTME"] = [[0, 1, 2]]
basic_names["TTME"] = ['TTME']
basic_specification["INVT"] = [[0, 1, 2, 3]]
basic_names["INVT"] = ['INVT']
basic_specification["HINC"] = [[1, 2]]
basic_names["HINC"] = [ 'HINC_train_bus']
basic_specification["PSIZE"] = [0]
basic_names["PSIZE"] = ['PSIZE_air']
mnl = pl.create_choice_model(data = model_data,
                alt_id_col="ALT_ID",
                obs_id_col="OBS_ID",
                choice_col="MODE",
                specification=basic_specification,
                model_type = "MNL",
                names=basic_names)
mnl.fit_mle(np.zeros(7))
mnl.get_statsmodels_summary()
# | -----------------------------------------------------------------
# |                   coef     std.err z       P>|z|   [0.025  0.975]
# | -----------------------------------------------------------------
# | ASC_air          5.6860   0.937   6.068   0.000   3.849   7.523
# | ASC_train        5.4034   0.603   8.959   0.000   4.221   6.585
# | ASC_bus          5.0128   0.623   8.051   0.000   3.792   6.233
# | TTME            -0.0992   0.011  -9.428   0.000  -0.12   -0.079
# | INVT            -0.0039   0.001  -4.489   0.000  -0.006  -0.002
# | HINC_train_bus  -0.0500   0.011  -4.484   0.000  -0.072  -0.028
# | PSIZE_air       -0.8997   0.245  -3.680   0.000  -1.379  -0.420
# |==================================================================

根據模型系數可以初步判定模型的合理性，例如：TTME的系數為負，可以解釋為當某個備選項站點等待時間延長，其被選擇的概率會降低；HINC_train_bus的系數為負，可以解釋為隨著家庭收入增加，選擇火車或長途汽車的概率會降低。這種定性的合理性判斷有利于我們判斷模型搭建是否合理。當然，如果想發揮模型真正的價值，還需要對模型進行量化解讀。

對MNL模型的解讀需要基于其預測功能，原理已經在1.2.3節進行了闡述，這里主要進行代碼實現，如代碼清單1-11所示。假設其他條件保持不變，因為火車提速，使得行程耗時降低20%，通過計算可知：

  • 飛機的選擇概率會由27.6%變為25.6%，降低了2.0%。
  • 火車的選擇概率會由30.0%變為36.2%，提升了6.2%。
  • 長途汽車的選擇概率會由14.3%變為12.7%，降低1.6%。
  • 自駕的選擇概率會由28.1%變為25.5%，降低2.6%。

代碼清單1-11　解讀MNL模型

# 創建用于預測的df
prediction_df = model_data[['OBS_ID', 'ALT_ID', 'MODE','TTME', 'INVT','HINC','PSIZE']]
choice_column = "MODE"
# 對火車耗時進行變化
def INVT(x,y):
    if x == 1:
        return y*0.8
    else:
        return y
prediction_df['INVT'] = prediction_df.apply(lambda x: INVT(x.ALT_ID, x.INVT), axis = 1)
# 默認情況下，predict()方法返回的結果是每個備選方案的選擇概率
prediction_array = mnl.predict(prediction_df)
# 存儲預測概率
prediction_df["MNL_Predictions"] = prediction_array
# 對比變化前后的概率
raw_probability = prediction_df.groupby(['ALT_ID'])['MODE'].mean()
new_probability = prediction_df.groupby(['ALT_ID'])['MNL_Predictions'].mean()
print("--------原概率--------")
print(raw_probability)
print("--------新概率--------")
print(new_probability)
# | --------原概率--------
# | ALT_ID
# | 0    0.276190
# | 1    0.300000
# | 2    0.142857
# | 3    0.280952
# | Name: MODE, dtype: float64
# | --------新概率--------
# | ALT_ID
# | 0    0.255643
# | 1    0.362788
# | 2    0.126937
# | 3    0.254632

1.3.4　使用嵌套Logit模型分析多種交通方式選擇問題

本節我們要學習的是嵌套Logit模型，顧名思義，嵌套Logit模型就是人為設置選擇的層次。我們假定選擇的過程存在層次性，如圖1-7所示，旅行家庭先進行飛機或陸地交通的選擇，如果選擇陸地交通，再決定到底是選火車、長途汽車還是自駕。該模型主要應用于不滿足IIA假定的分析場景。

000

圖1-7　NL模型的場景邏輯示意圖

我們延續MNL模型的代碼進行闡述，如代碼清單1-12所示，NL模型與MNL模型的不同之處在于NL模型需要在聲明每個備選項的效用表達式之前，單獨聲明NL模型的層次。

代碼清單1-12　NL模型代碼

# 聲明嵌套形式
nest_membership = OrderedDict()
nest_membership["air_Modes"] = [0]
nest_membership["ground_Modes"] = [1, 2, 3]

# 聲明備選項的效用函數
basic_specification = OrderedDict()
basic_names = OrderedDict()
basic_specification["intercept"] = [0, 1, 2]
basic_names["intercept"] = ['ASC_air', 'ASC_train', 'ASC_bus']
# 可以靈活指定備選項屬性的影響方式
basic_specification["TTME"] = [[0, 1, 2]]
basic_names["TTME"] = ['TTME']
basic_specification["INVT"] = [[0, 1, 2, 3]]
basic_names["INVT"] = ['INVT']
# 也可以靈活指定決策者的影響方式，但需要注意的是，由于每個選項的決策者屬性都一樣，因此保證可估計性，只對部分選項生效
basic_specification["HINC"] = [[1, 2]]
basic_names["HINC"] = [ 'HINC_train_bus']
basic_specification["PSIZE"] = [0]
basic_names["PSIZE"] = ['PSIZE_air']

# 模型創建
nested_logit = pl.create_choice_model(data = model_data,
                alt_id_col="ALT_ID",
                obs_id_col="OBS_ID",
                choice_col="MODE",
                specification=basic_specification,
                model_type = "Nested Logit",
                names=basic_names,
                nest_spec=nest_membership)
nested_logit.fit_mle(np.zeros(9))
nested_logit.summary
# | -----------------------------------------------------------------
# |                 parameters   std_err   t_stats   p_values
# | -----------------------------------------------------------------
# | air_Modes        0.0000      NaN       NaN       NaN     
# | ground_Modes     0.8187      0.668     1.225     0.220   
# | ASC_air          4.0002      1.022     3.914     0.000   
# | ASC_train        4.2224      0.744     5.672     0.000   
# | ASC_bus          3.9471      0.747     5.285     0.000   
# | TTME            -0.0787      0.013    -6.180     0.000   
# | INVT            -0.0038      0.001    -4.718     0.000   
# | HINC_train_bus  -0.0364      0.010    -3.513     0.000   
# | PSIZE_air       -0.7535      0.244    -3.088     0.002   
# |==================================================================

從模型結果來看，ground_Modes的系數并不顯著，并且嘗試用其他層次條件進行建模亦會發現，層次的系數也不顯著，這從側面表明了MNL模型在該問題上的應用是合理的。當然，多數情況下IIA假定并不一定被滿足，NL模型會比MNL模型更科學和嚴謹。

NL模型的解讀與MNL模型一致，這里不再贅述。

主站蜘蛛池模板： 揭阳市| 桐城市| 孟村| 昔阳县| 枝江市| 旬邑县| 分宜县| 大石桥市| 呈贡县| 南部县| 阳泉市| 临沧市| 北宁市| 弋阳县| 定日县| 克拉玛依市| 驻马店市| 宁陵县| 桓台县| 墨脱县| 江川县| 湾仔区| 广东省| 绥化市| 乡城县| 青海省| 长宁县| 京山县| 佛教| 枣强县| 塔城市| 永州市| 南昌市| 安塞县| 阿拉善右旗| 从江县| 永兴县| 红河县| 舟曲县| 海林市| 汝阳县|