第10章解釋了選擇模型理論以及一些流行的模型：多項式Logit模型、嵌套Logit模型以及混合Logit模型。
本章中，會學習以下技巧：
·准備數據集以估算離散選擇模型
·估算知名的多項Logit模型
·測試來自無關選項的獨立性沖突
·用巢式Logit模型處理IIA沖突
·用混合Logit模型處理復雜的替代模式
10.1導論
DCM（Discrete choice models，離散選擇模型）的目標是預測出一個人會做怎樣的選擇。這些模型與邏輯回歸有共同之處，但在誤差項的分布假設方面又有不同。
DCM以隨機效用理論作為理論基礎，並假設一個理性的人，總是做出使效用最大化的選擇。
例如，假設你要做一個二選一的決定（這兩個選項是沒有共同點的）；我們考慮選擇騎自行車還是開車去上班的情形。騎自行車幾乎無成本（當然，這里假設你已經有一輛自行車了，也忽略了騎自行車消耗的能量），而開車的成本是3塊錢。騎自行車上班要花45min，而開車只要8min。假設價格系數和時間系數分別是-2和-0.13，可以計算出，騎自行車的效用是-5.85，開車的效用是-7.04。所以，在效用最大化的假設下，你會選擇騎自行車，而不是開車。這個簡化的計算不考慮性格偏好因素。
我們生成用於本章技巧的數據集，是西雅圖與洛杉磯之間4個航班16個等級中做出的10000個選擇；各選項的詳情如下：
下一章會介紹如何轉換這個數據集，用於估算多種DCM。
我們使用Python Biogeme估算DCM，這是由Michel Bierlaire博士在Bierlaire，M.（2003）；BIOGEME：A free package for the estimation of discrete choice models中發表的。雖然不是純Python（估算的引擎是用C++寫的），但這是估算DCM最好的開源軟件。Python Biogeme還在活躍更新中，可從http://biogeme.epfl.ch下載；徑直訪問Download部分，選擇你的操作系統對應的版本。由於軟件包是免費的，因此你需要在Yahoo的Biogeme用戶群中注冊：https://groups.yahoo.com/neo/groups/biogeme/info。
下載軟件包之后，按照http://biogeme.epfl.ch/install.html上的指南安裝。裝好后就可以通過下面的命令調用了：
>>pythonbiogeme <model_file> <dataset>
注意：這段話來自原書作者
/*
我們使用軟件包時遇到了問題。這個軟件包要求Python 3.2之上的版本。盡管本書使用的是Python 3.5（安裝在~/anaconda文件夾中），
但是Python Biogeme會到/usr/local/Frameworks/Python.framework/Versions/3.5路徑下尋找Python，這個路徑在我們的機器上根本不存在。
我們只好手動創建這個路徑，然后重新安裝軟件。 我們也在Windows上試了，可以正常運行。不過，有些朋友可能會在安裝和運行時遇到問題，這也是要大家注意一下的。 如果你遇到了任何問題，你可以在Biogeme用戶群中詢問；一般很快就有回復了。如果你實在沒法子了，也可以給原作者發郵件，mail@tomdrabas.com，能幫的原作者一定幫。 */
以下來自邀月的安裝步驟：
/* pip install biogeme Looking in indexes: https://pypi.tuna.tsinghua.edu.cn/simple Collecting biogeme Using cached https://pypi.tuna.tsinghua.edu.cn/packages/40/aa/cbad1883d2914451e0714981ec2332e443dabfe5b5d874199dfbc1070b98/biogeme-3.2.5-cp37-cp37m-win_amd64.whl (816 kB) Requirement already satisfied: scipy in d:\tools\python37\lib\site-packages (from biogeme) (1.4.0) Requirement already satisfied: numpy in d:\tools\python37\lib\site-packages (from biogeme) (1.18.2) Requirement already satisfied: cython in d:\tools\python37\lib\site-packages (from biogeme) (0.29.14) Requirement already satisfied: pandas in d:\tools\python37\lib\site-packages (from biogeme) (0.25.3) Requirement already satisfied: unidecode in d:\tools\python37\lib\site-packages (from biogeme) (1.1.1) Requirement already satisfied: python-dateutil>=2.6.1 in d:\tools\python37\lib\site-packages (from pandas->biogeme) (2.8.1) Requirement already satisfied: pytz>=2017.2 in d:\tools\python37\lib\site-packages (from pandas->biogeme) (2019.3) Requirement already satisfied: six>=1.5 in d:\tools\python37\lib\site-packages (from python-dateutil>=2.6.1->pandas->biogeme) (1.12.0) Installing collected packages: biogeme Successfully installed biogeme-3.2.5 FINISHED */
10.2准備數據集以估算離散選擇模型
為這些技巧准備的數據包括了10000條選擇情形。第一列是選擇的ID，第二列是用戶可選的選項列表；不是所有情形下都要從所有選項中選擇。
准備：需要安裝好pandas、NumPy和正則表達式模塊。
步驟：Python Biogeme要求數據以制表符分隔，每一行包括所有選項的屬性、選出的選項以及做出選擇時每個選項是否可選。值的類型只能是數字。下面的代碼可在dcm_dataPrep.py文件中看到：
1 import pandas as pd 2 import numpy as np 3 import re 4 5 # read datasets 6 choices_filename = '../../Data/Chapter10/choices.csv' 7 alternatives_filename = '../../Data/Chapter10/options.json' 8 9 choices = pd.read_csv(choices_filename) 10 alternatives = pd.read_json(alternatives_filename).transpose() 11 12 # retrieve all considered alternatives 取回所有考慮的選項 13 considered = [alt.split(';') 14 for alt in list(choices['available'])] 15 16 # create flag of all available alternatives 標出所有可能的選項 17 available = [] 18 alternatives_list = list(alternatives.index) 19 alternatives_list = [ 20 re.sub(r'\.', r'_', alt) for alt in alternatives_list] 21 no_of_alternatives = len(alternatives_list) 22 23 for cons in considered: 24 f = [0] * len(alternatives_list) 25 26 cons = [re.sub(r'\.', r'_', alt) for alt in cons] 27 28 for i, alt in enumerate(alternatives_list): 29 if alt in cons: 30 f[i] = 1 31 32 available.append(list(f)) 33 34 # append to the choices DataFrame 追加到選項DataFrame中 35 alternatives_av = [alt + '_AV' for alt in alternatives_list] 36 available = pd.DataFrame(available, columns=alternatives_av) 37 choices = choices.join(available) 38 39 # drop the available column as we don't need it anymore 40 #丟掉不需要的列 41 del choices['available'] 42 43 # encode the choice variable 編碼選擇變量 44 choice = list(choices['choice']) 45 choice = [re.sub(r'\.', r'_', alt) for alt in choice] 46 choice = [alternatives_list.index(c) + 1 for c in choice] 47 48 choices['choice'] = pd.Series(choice) 49 50 # and add the alternatives' attributes 加上選項屬性 51 # first, normalize price to be between 0 and 1 首先，將價格歸一為0到1之間的數 52 max_price = np.max(alternatives['price']) 53 alternatives['price'] = alternatives['price'] / max_price 54 55 # next, create a vector with all attributes 然后，創建所有選項的向量 56 attributes = [] 57 attributes_list = list(alternatives.values) 58 59 for attribute in attributes_list: 60 attributes += list(attribute) 61 62 # fill in to match the number of rows 填充以匹配行數 63 attributes = [attributes] * len(choices) 64 65 # and their names 名字 66 attributes_names = [] 67 68 for alternative in alternatives_list: 69 for attribute in alternatives.columns: 70 attributes_names.append(alternative + '_' + attribute) 71 72 # convert to a DataFrame 轉化為DataFrame 73 attributes = pd.DataFrame(attributes, 74 columns=attributes_names) 75 76 # and join with the main dataset 和主數據集連接 77 choices = choices.join(attributes) 78 79 # save as a TSV with .dat extension 輸出為.dat文件 80 with open('../../Data/Chapter10/choices.dat', 'w') as f: 81 f.write(choices.to_csv(sep='\t', index=False))
原理：按照慣例，我們先導入會用到的模塊。
數據集位於本書Git倉庫的Data/Chapter10/文件夾中。它已拆成兩個文件：observations.csv是那10000個選擇樣本，options.json是帶上屬性的所有選項。
observations.csv文件格式如下：
/* choice available AA777.4.V AA777.1.C;AA777.2.Z;AA777.4.V;AS666.1.C;AS666.3.Y;AS666.4.V;DL001.2.Z;DL001.3.Y;DL001.4.V;UA110.1.C;UA110.3.Y;UA110.4.V UA110.3.Y AA777.1.C;AA777.2.Z;AA777.3.Y;AA777.4.V;AS666.1.C;AS666.2.Z;AS666.4.V;DL001.2.Z;UA110.2.Z;UA110.3.Y;UA110.4.V DL001.1.C AA777.1.C;AA777.2.Z;AA777.3.Y;AA777.4.V;AS666.3.Y;AS666.4.V;DL001.1.C;DL001.3.Y;UA110.3.Y;UA110.4.V */
第一列是做出的選擇：比如選擇了AA777航班的V等級。用戶是從12個選項中選擇的：用；分隔的12個選項的列表。
options.json文件格式如下：
/* { "AA777.1.C": { "compartment": 1, "frequentFlyer": 1, "price": 410, "refund": 1 }, "AA777.2.Z": { "compartment": 1, "frequentFlyer": 1, "price": 320, "refund": 0 }, */
compartment值為1表示是商務艙，frequentFlyer和refund都為1，意味着可以累計常旅積分，並可以免費退票。
使用pandas讀取CSV與JSON格式的文件，參考本書1.2節。
讀入數據集后，我們先關注考慮的選項列表。從choices DataFrame對象提取出可選的列，並將行拆成獨立的選項。我們將用這個創建一個帶有所有可選選項的DataFrame對象。
要做到這一點，需要所有可選選項的列表。從DataFrame對象alternatives提取所有索引。不過這些值里包含了點號，比如AA777.1.C。如果直接使用Python Biogeme會有問題，所以用re包的.sub（...）方法將點號替換成下划線。
注意，正則表達式r'.'匹配的是除換行符之外的所有字符，所以要用r'\.'。
.sub（...）方法以正則表達式為第一個參數，以替代字符串為第二個參數，以要處理的字符串為第三個參數。
最后，遍歷所有記錄，標出可選的選項。首先，創建和選項一樣多的0。然后，為了兼容alternatives_list，也替換了cons里的點號。最后，遍歷所有選項，將其中可用的標為1。最后得到了這樣的列表：
創建相應的列名，以便在Python Biogeme中使用；要達到這一目標，在每個選項的名字后加上_AV后綴。
最后，創建一個DataFrame對象available，其中包含可用選項的向量，與DataFrame對象choices連接，並刪除后面用不到的available列：
與上圖合並
可以看到，選擇了AA777航班的V等級，和期望的一樣，這是決策者的考慮。
然后，需要給choice列中的選項編碼。首先，取出DataFrame中的所有記錄，並用下划線替換點號。接着，遍歷所有元素，找到alternatives_list中元素的.index（...）。這樣做之后，我們和數據集的其余部分保持一致。最后，用投射到pandas的Series對象中已編碼的choice列表替換choice列。
最后一項是往數據集中加所有選項的所有屬性。前面提過，數據集中的每一行需要包括所有選項的所有信息。
就本例而言，這是相同信息的列表，自頂而下。不過，這種處理方式一般在處理更普通的選項時會有用得多，比如，用火車還是用汽車，這個時候價格和距離可能都不同。如果數據集有不同的目的地，每一行會有不同的屬性。
首先，將價格歸一化到0和1之間。直接從數據集中取出最大價格，然后所有價格都除以這個值。也可以用Python Biogeme做到這一點，但是推薦直接處理，僅僅用Python Biogeme估算模型。
然后，提取DataFrame對象alternatives的.values屬性，以獲取所有屬性的值。attributes_list本質上就是一個NumPy數組列表（pandas的默認數據結構很多都由NumPy承襲而來）：
/* [array([1., 1., 1., 1.]), array([1. , 1. , 0.7804878, 0. ]), array([0. , 1. , 0.51219512, 1. ]), array([0. , 0. , 0.36585366, 0. ]), array([1. , 1. , 0.95121951, 1. ]), array([1. , 1. , 0.76829268, 0. ]), array([0. , 1. , 0.47560976, 1. ]), array([0. , 0. , 0.31707317, 0. ]), array([1. , 1. , 0.97560976, 1. ]), array([1. , 1. , 0.7804878, 0. ]), array([0. , 1. , 0.48780488, 1. ]), array([0. , 0. , 0.31707317, 0. ]), array([1. , 1. , 0.92682927, 1. ]), array([1. , 0. , 0.75609756, 0. ]), array([0. , 1. , 0.48780488, 1. ]), array([0. , 0. , 0.30487805, 0. ])] */
因為我們需要一個帶有所有屬性的大向量，所以遍歷attributes_list，並給每個數組擴展attributes列表。對於這10000行，都要復制一遍，所以將attributes列表乘以choices的長度，即行數。
再之后，為所有屬性創建列名。創建好attributes列表並放置好列名后，創建一個DataFrame對象，並與初始數據集連接。
最后一步，將數據集導出成一個.dat文件，以制表符分隔（以與Biogeme的傳統保持一致）。
更多:：如果你的數據集是CSV格式，你可以使用Biogeme的biopreparedata工具創建數據集。另外，biocheckdata可用於檢查數據集是否符合Biogeme的標准。
10.3
估算知名的多項Logit模型
MNL（Multinomial Logit，多項Logit）模型是由Daniel McFadden於1973年在其研討會論文中提出的，http://eml.berkeley.edu/reprints/mcfadden/zarembka.pdf。
這個模型基於十分嚴格的假設：IID（Independent and Identically Distributed，獨立同分布）的誤差項和IIA（Independence from Irrelevant Alternatives，獨立於無關選項）。
IID假設所有選項效用函數的誤差項是獨立（不相干）的，並服從同樣的分布。（就MNL而言，I型極值分布，一般稱為Gumbel分布，這是以發現並分析這個分布的E.J.Gumbel命名的。）
另一方面，IIA假設選項間可能性的比率是一個常量，即，從考慮范圍中移除一個選項，不會改變其余選項的比率。考慮這樣一個場景：你要從自行車、火車、汽車中挑一個作為上班的交通工具。假設每個選項的可能性是相同的，也就是每項被選中的概率都是1/3。在IIA假設下，如果汽車壞了，那么另兩項的比率保持恆定，即被選中的可能性是一半一半。下一技巧會解釋為什么這是個嚴格假設。
MNL中的概率計算式如下：
截圖（公式）
U（i）是每個選項的效用，定義為：
截圖（公式）
其中αi是與選項相關的常量，βi是由選項i所對應的屬性系數構成的向量，Xi是選項i的屬性向量，εi是誤差項。
准備：需要安裝好Python Biogeme包。
步驟：Python Biogeme的語法和Python本身的很像。對於Biogeme而言，Python只是一層很薄的包裝，里面還是C++代碼（MNL目錄中的dcm_mnl.py文件）。下面的代碼經過刪減了：
1 from biogeme import * 2 from headers import * 3 from loglikelihood import * 4 from statistics import * 5 6 # Specify parameters to be estimated 7 # 8 # Arguments: 9 # - 1 Name, typically, the same as the variable. 10 # - 2 Starting value. 11 # - 3 Lower bound. 12 # - 4 Upper bound. 13 # - 5 flag whether to estimate the parameter (0) 14 # or keep it fixed (1). 15 16 # add the coefficients to be estimated 17 18 ASC = Beta('ASC',0,-10,10,1,'ASC') 19 C_price = Beta('C_price',0,-10,10,0,'C price') 20 Z_price = Beta('Z_price',0,-10,10,0,'Z price') 21 Y_price = Beta('Y_price',0,-10,10,0,'Y price') 22 V_price = Beta('V_price',0,-10,10,0,'V price') 23 24 B_comp = Beta('B_comp',0,-10,10,0,'compartment') 25 B_refund = Beta('B_refund',0,-3,3,0,'refund') 26 27 # Utility functions 28 29 # compartment attributes 30 c = {} 31 32 c[1] = AA777_1_C_compartment 33 c[2] = AA777_2_Z_compartment 34 c[3] = AA777_3_Y_compartment 35 ... 36 c[13] = UA110_1_C_compartment 37 c[14] = UA110_2_Z_compartment 38 c[15] = UA110_3_Y_compartment 39 c[16] = UA110_4_V_compartment 40 41 # price attributes 42 p = {} 43 44 p[1] = AA777_1_C_price 45 p[2] = AA777_2_Z_price 46 p[3] = AA777_3_Y_price 47 ... 48 p[14] = UA110_2_Z_price 49 p[15] = UA110_3_Y_price 50 p[16] = UA110_4_V_price 51 52 # refund attributes 53 r = {} 54 55 r[1] = AA777_1_C_refund 56 r[2] = AA777_2_Z_refund 57 r[3] = AA777_3_Y_refund 58 ... 59 r[13] = UA110_1_C_refund 60 r[14] = UA110_2_Z_refund 61 r[15] = UA110_3_Y_refund 62 r[16] = UA110_4_V_refund 63 64 # The dictionary of utilities is constructed. 65 V = {} 66 67 V[1] = C_price * p[1] + B_refund * r[1] + B_comp * c[1] 68 V[2] = Z_price * p[2] + B_refund * r[2] + B_comp * c[2] + ASC 69 V[3] = Y_price * p[3] + B_refund * r[3] + B_comp * c[3] 70 ... 71 V[13] = C_price * p[13] + B_refund * r[13] + B_comp * c[13] 72 V[14] = Z_price * p[14] + B_refund * r[14] + B_comp * c[14] 73 V[15] = Y_price * p[15] + B_refund * r[15] + B_comp * c[15] 74 V[16] = Y_price * p[16] + B_refund * r[16] + B_comp * c[16] 75 76 # availability flags 77 availability = { 78 1: AA777_1_C_AV, 79 2: AA777_2_Z_AV, 80 3: AA777_3_Y_AV, 81 4: AA777_4_V_AV, 82 5: AS666_1_C_AV, 83 6: AS666_2_Z_AV, 84 7: AS666_3_Y_AV, 85 8: AS666_4_V_AV, 86 9: DL001_1_C_AV, 87 10: DL001_2_Z_AV, 88 11: DL001_3_Y_AV, 89 12: DL001_4_V_AV, 90 13: UA110_1_C_AV, 91 14: UA110_2_Z_AV, 92 15: UA110_3_Y_AV, 93 16: UA110_4_V_AV 94 } 95 96 # The choice model is a logit, with availability conditions 97 logprob = bioLogLogit(V, availability, choice) 98 99 # Defines an iterator on the data 100 rowIterator('obsIter') 101 102 # Define the likelihood function for the estimation 103 BIOGEME_OBJECT.ESTIMATE = Sum(logprob, 'obsIter') 104 105 # Statistics 106 nullLoglikelihood(availability,'obsIter') 107 choiceSet = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16] 108 cteLoglikelihood(choiceSet, choice, 'obsIter') 109 availabilityStatistics(availability, 'obsIter') 110 111 # Parameters 112 BIOGEME_OBJECT.PARAMETERS['optimizationAlgorithm'] = 'BIO' 113 BIOGEME_OBJECT.PARAMETERS['numberOfThreads'] = '8' 114 115 BIOGEME_OBJECT.FORMULAS['AA777 C utility'] = V[1] 116 BIOGEME_OBJECT.FORMULAS['AA777 Z utility'] = V[2] 117 ... 118 BIOGEME_OBJECT.FORMULAS['DL001 Y utility'] = V[11] 119 BIOGEME_OBJECT.FORMULAS['DL001 V utility'] = V[12] 120 BIOGEME_OBJECT.FORMULAS['UA110 C utility'] = V[13] 121 BIOGEME_OBJECT.FORMULAS['UA110 Z utility'] = V[14] 122 BIOGEME_OBJECT.FORMULAS['UA110 Y utility'] = V[15] 123 BIOGEME_OBJECT.FORMULAS['UA110 V utility'] = V[16]
原理：要執行代碼，可以使用本章主目錄下的run_pythonbiogeme.sh腳本，或者（假設你在MNL目錄）執行下面這行：
/* d: cd D:\Java2018\practicalDataAnalysis\Codes\Chapter10\MNL D:\tools\Python37\python D:\tools\Python37\Lib\site-packages\biogeme\biogeme.py dcm_mnl ../../../Data/Chapter10/choices.dat biogeme.py dcm_mnl ../../../Data/Chapter10/choices.dat import os >>> os.chdir('D:\Java2018\practicalDataAnalysis\Codes\Chapter10\MNL') os.getcwd() pythonbiogeme dcm_mnl ../../../Chapter10/choices.dat import os os.chdir('D:\Java2018\practicalDataAnalysis\Codes\Chapter10\MNL') os.getcwd() 'D:\\Java2018\\practicalDataAnalysis\\Codes\\Chapter10\\MNL' >>> biogeme dcm_mnl ../../../Chapter10/choices.dat SyntaxError: invalid syntax >>> biogeme dcm_mnl ../../../Chapter10/choices.dat SyntaxError: invalid syntax pip install biogeme --no-cache-dir import biogeme.version as ver print(ver.getText()) biogeme 3.2.3 [September 25, 2019] Version entirely written in Python Home page: http://biogeme.epfl.ch Submit questions to https://groups.google.com/d/forum/biogeme Michel Bierlaire, Transport and Mobility Laboratory, Ecole Polytechnique Fédérale de Lausanne (EPFL) */
run_pythonbiogeme.sh先清除pythonbiogeme之前生成的文件，再重新估算模型。
使用run_pythonbiogeme.sh腳本也很簡單（簡單一行命令）：
/* ../run_pythonbiogeme.sh dcm_mnl ../../../Data/Chapter10/choices.dat */
隨便哪種方式，你都會看到類似下面的內容：
執行不了！！后面再研究。
Init.Log-likelihood給出了log-likelihood函數的初始值，即在初始系數下log-likelihood函數的值。每個選項被選中的可能性都一樣。這個值后面會用來計算rho平方值——這個指標等價於線性回歸中的R2。估算DCM的目標是最小化log-likelihood函數的絕對值。然后，輸出的表展示了估算的進度；你可以看到f（x）（log-likelihood函數）的值越來越小，一直到第7次循環終止，因為第6次循環和第7次循環之間的差足夠小了，所以函數收斂了。表下面輸出了估算參數的值，不過我們后面會在dcm_mnl.html文件中查看。
現在仔細考察模型文件。
照例，先加載需要的模塊。biogeme、headers、log-likelihood以及statistics這幾個模塊幾乎是所有模塊都需要的。
然后用Beta（...）方法指定系數。這個方法的第一個參數是參數名（一般情況下，就是對象名）。第二個參數指定了系數的初始值。第三個參數和第四個參數決定了系數值的下界和上界，第五個參數指定系數值是要固定（不估算）還是估算。
對於可識別的模型，要歸一化一個參數，即固定為0。通常，這個參數是ASC（alternative specific constants，方案特定變量）中的一個。參考K.Train《Properties of Discrete Choice Models》一書的第2章http://eml.berkeley.edu/books/choice2nd/Ch02_p9-33.pdf。
Beta（...）方法最后一個參數是將要用於報告中的易記名字。
指定了系數后，創建屬性名的字典。
這些名字必須同.dat文件中完全相同。
為每個選項指定效用函數時，只使用compartment、price和refund屬性。注意，為了避免混淆，存儲同一選項的屬性時，字典中用的都是同一個數字。
最終，字典V里有全部的效用函數。只對選項2（AA777航班，Z艙位）加進了ASC。也要注意到，B_refund系數和B_comp系數不是選項特定的，而是通用的——所有選項都相同。不同航班之間，等價艙位的價格系數是共用的。
availability對象對每個選項，指定了可用性的標志。
最后，要估算參數，先指定logprob對象。使用bioLogLogit（...）方法估算MNL。這個方法以效用的字典作為第一個參數，可用性的字典作為第二個參數，以及存儲choice的變量名。然后在數據集上定義一個迭代器，命名為obsIter。為BIOGEME_OBJECT指定ESTIMATE變量：這是要最小化整個數據集的對數概率和。
剩下的部分就是處理估算的統計數據。nullLogLikelihood（...）方法計算空的log-likelihood並將其包括在輸出中。第一個參數是每個選項可用性的字典，第二個參數是行迭代器的名字。cteLoglikelihood（...）計算只有常數參數時log-likelihood函數的值。它以choiceSet（所有選項的列表）作為第一個參數，存儲選擇的變量名作為第二個參數，記錄的迭代器作為第三個參數。得到的統計數據是availabilityStatistics（...），將返回數據集中每個可用選項的起算時間。
最后指定BIOGEME_OBJECT.PARAMETERS和.FORMULAS。使用.PARAMETERS（...）方法，可以向Biogeme及其solver傳入各種參數；這里，使用BIO優化器（這是Biogeme自己的優化算法），並指定執行時啟動的線程數。.FORMULAS（...）允許我們給效用函數起一個易記名字，並在最終的報告中使用。
Python Biogeme給出的報告如下（刪減版；完整版可在瀏覽器中訪問dcm_mnl.html查看）：
報告給出了多項統計數據。有些已介紹過了，這里就不再重復。看一下init.模型數據的Rho-square-bar：這個數據告訴你，模型表現如何。0.20以上的值就非常好（相當於線性模型中85%的R2），0.15到0.2之間的值不錯，0.1到0.15之間的值是可接受的，而不到0.1就是很差的模型了。
重點是，Diagnostic要看Convergence reached。后面我們會看到，Run time也是一項重要數據，尤其是在處理復雜模型時。
考慮Estimated的參數，你需要保留統計上顯著（p-values<0.05）的變量。模型中，參數都是統計上顯著的。而且，所有價格系數都是負的，意味着更高的價格將降低選中的可能性。意外的是，如果提供一個商務艙的選項，它被選中的可能性就上升了。
參考：
關於MNL，可以參考：
·Kenneth Train的一本（免費）好書：http://eml.berkeley.edu/books/choice2.html；
·一份演示文稿：http://www.bauer.uh.edu/rsusmel/phd/ec1-20.pdf；
·Python Biogeme的介紹：http://biogeme.epfl.ch/documentation/pythonfirstmodel-2.4.pdf。
10.4
測試來自無關選項的獨立性沖突
MNL基於一個相當嚴格的IIA假設，選項的概率之間的比率保持不變。這其實只適用於選項集合沒有共同特征的情況，換種說法就是，選項之間不相關。
最著名的IIA反例就是紅藍大巴悖論。考慮這樣的場景，你要在汽車、火車、藍色大巴之間做選擇。簡化一下，假設每個選項的概率都是1/3。在IIA假設下，如果加入一個紅色大巴的選項，選項概率之間的比率恆定，所以各選項的概率都成了1/4。
然而，現實中，顏色因素會有這么大的影響？如果沒有，我們實際上還是在汽車、火車、大巴之間做選擇。加入一個新的大巴選項，不該影響汽車和火車的概率，其還應該是1/3。而紅色與藍色的概率對半分，兩個大巴選項的概率應該都是1/6。
這就違反了IIA假設：P（car）/P（train）沒變，P（car）/P（blue bus）和P（train）/P（blue bus）都變了。
本技巧中，我們會考察數據集中是否違反了IIA。下兩項技巧中，我們將提供能處理選項之間相關性的模型。
准備：需要安裝好Python Biogeme包。
步驟：前一技巧中估算了MNL模型，我們來計算各選項的概率。然后移除第一個選項，重新估算模型，再對新模型重新計算。最后比較概率之間的比率。本技巧用到的文件在TestingIIA目錄下。
代碼和估算的差不多，所以只討論不同的地方（TestingIIA/dcm_mnl_simul.py文件）：
1 from biogeme import * 2 from headers import * 3 from loglikelihood import * 4 from statistics import * 5 6 # Specify parameters to be estimated 7 # 8 # Arguments: 9 # - 1 Name, typically, the same as the variable. 10 # - 2 Starting value. 11 # - 3 Lower bound. 12 # - 4 Upper bound. 13 # - 5 flag whether to estimate the parameter (0) 14 # or keep it fixed (1). 15 16 # add the coefficients to be estimated 17 C_price = Beta('C_price',-7.29885,-10,10,0,'C price' ) 18 V_price = Beta('V_price',-5.07495,-10,10,0,'V price' ) 19 Y_price = Beta('Y_price',-4.40754,-10,10,0,'Y price' ) 20 Z_price = Beta('Z_price',-8.70638,-10,10,0,'Z price' ) 21 22 ASC = Beta('ASC',0,-10,10,1,'ASC' ) 23 B_comp = Beta('B_comp',3.52571,-10,10,0,'compartment' ) 24 B_refund = Beta('B_refund',-0.718748,-3,3,0,'refund' ) 25 26 # Utility functions 27 28 # The data are associated with the alternative index 29 # compartment attributes 30 c = {} 31 32 c[1] = AA777_1_C_compartment 33 c[2] = AA777_2_Z_compartment 34 c[3] = AA777_3_Y_compartment 35 c[4] = AA777_4_V_compartment 36 c[5] = AS666_1_C_compartment 37 c[6] = AS666_2_Z_compartment 38 c[7] = AS666_3_Y_compartment 39 c[8] = AS666_4_V_compartment 40 c[9] = DL001_1_C_compartment 41 c[10] = DL001_2_Z_compartment 42 c[11] = DL001_3_Y_compartment 43 c[12] = DL001_4_V_compartment 44 c[13] = UA110_1_C_compartment 45 c[14] = UA110_2_Z_compartment 46 c[15] = UA110_3_Y_compartment 47 c[16] = UA110_4_V_compartment 48 49 # price attributes 50 p = {} 51 52 p[1] = AA777_1_C_price 53 p[2] = AA777_2_Z_price 54 p[3] = AA777_3_Y_price 55 p[4] = AA777_4_V_price 56 p[5] = AS666_1_C_price 57 p[6] = AS666_2_Z_price 58 p[7] = AS666_3_Y_price 59 p[8] = AS666_4_V_price 60 p[9] = DL001_1_C_price 61 p[10] = DL001_2_Z_price 62 p[11] = DL001_3_Y_price 63 p[12] = DL001_4_V_price 64 p[13] = UA110_1_C_price 65 p[14] = UA110_2_Z_price 66 p[15] = UA110_3_Y_price 67 p[16] = UA110_4_V_price 68 69 # refund attributes 70 r = {} 71 72 r[1] = AA777_1_C_refund 73 r[2] = AA777_2_Z_refund 74 r[3] = AA777_3_Y_refund 75 r[4] = AA777_4_V_refund 76 r[5] = AS666_1_C_refund 77 r[6] = AS666_2_Z_refund 78 r[7] = AS666_3_Y_refund 79 r[8] = AS666_4_V_refund 80 r[9] = DL001_1_C_refund 81 r[10] = DL001_2_Z_refund 82 r[11] = DL001_3_Y_refund 83 r[12] = DL001_4_V_refund 84 r[13] = UA110_1_C_refund 85 r[14] = UA110_2_Z_refund 86 r[15] = UA110_3_Y_refund 87 r[16] = UA110_4_V_refund 88 89 # The dictionary of utilities is constructed. 90 V = {} 91 92 V[1] = C_price * p[1] + B_refund * r[1] + B_comp * c[1] 93 V[2] = Z_price * p[2] + B_refund * r[2] + B_comp * c[2] + ASC 94 V[3] = Y_price * p[3] + B_refund * r[3] + B_comp * c[3] 95 V[4] = V_price * p[4] + B_refund * r[4] + B_comp * c[4] 96 V[5] = C_price * p[5] + B_refund * r[5] + B_comp * c[5] 97 V[6] = Z_price * p[6] + B_refund * r[6] + B_comp * c[6] 98 V[7] = Y_price * p[7] + B_refund * r[7] + B_comp * c[7] 99 V[8] = V_price * p[8] + B_refund * r[8] + B_comp * c[8] 100 V[9] = C_price * p[9] + B_refund * r[9] + B_comp * c[9] 101 V[10] = Z_price * p[10] + B_refund * r[10] + B_comp * c[10] 102 V[11] = Y_price * p[11] + B_refund * r[11] + B_comp * c[11] 103 V[12] = V_price * p[12] + B_refund * r[12] + B_comp * c[12] 104 V[13] = C_price * p[13] + B_refund * r[13] + B_comp * c[13] 105 V[14] = Z_price * p[14] + B_refund * r[14] + B_comp * c[14] 106 V[15] = Y_price * p[15] + B_refund * r[15] + B_comp * c[15] 107 V[16] = Y_price * p[16] + B_refund * r[16] + B_comp * c[16] 108 109 # availability flags 110 availability = { 111 1: AA777_1_C_AV, 112 2: AA777_2_Z_AV, 113 3: AA777_3_Y_AV, 114 4: AA777_4_V_AV, 115 5: AS666_1_C_AV, 116 6: AS666_2_Z_AV, 117 7: AS666_3_Y_AV, 118 8: AS666_4_V_AV, 119 9: DL001_1_C_AV, 120 10: DL001_2_Z_AV, 121 11: DL001_3_Y_AV, 122 12: DL001_4_V_AV, 123 13: UA110_1_C_AV, 124 14: UA110_2_Z_AV, 125 15: UA110_3_Y_AV, 126 16: UA110_4_V_AV 127 } 128 129 # The choice model is a logit, with availability conditions 130 probAA777_C = bioLogit(V, availability, 1) 131 probAA777_Z = bioLogit(V, availability, 2) 132 probAA777_Y = bioLogit(V, availability, 3) 133 probAA777_V = bioLogit(V, availability, 4) 134 probAS666_C = bioLogit(V, availability, 5) 135 probAS666_Z = bioLogit(V, availability, 6) 136 probAS666_Y = bioLogit(V, availability, 7) 137 probAS666_V = bioLogit(V, availability, 8) 138 probDL001_C = bioLogit(V, availability, 9) 139 probDL001_Z = bioLogit(V, availability, 10) 140 probDL001_Y = bioLogit(V, availability, 11) 141 probDL001_V = bioLogit(V, availability, 12) 142 probUA110_C = bioLogit(V, availability, 13) 143 probUA110_Z = bioLogit(V, availability, 14) 144 probUA110_Y = bioLogit(V, availability, 15) 145 probUA110_V = bioLogit(V, availability, 16) 146 147 # Defines an itertor on the data 148 rowIterator('obsIter') 149 150 # exclude observations where AA777 C was selected 151 exclude = choice == 1 152 BIOGEME_OBJECT.EXCLUDE = exclude 153 154 # simulate 155 simulate = { 156 'P_AA777_C': probAA777_C, 157 'P_AA777_Z': probAA777_Z, 158 'P_AA777_Y': probAA777_Y, 159 'P_AA777_V': probAA777_V, 160 'P_AS666_C': probAS666_C, 161 'P_AS666_Z': probAS666_Z, 162 'P_AS666_Y': probAS666_Y, 163 'P_AS666_V': probAS666_V, 164 'P_DL001_C': probDL001_C, 165 'P_DL001_Z': probDL001_Z, 166 'P_DL001_Y': probDL001_Y, 167 'P_DL001_V': probDL001_V, 168 'P_UA110_C': probUA110_C, 169 'P_UA110_Z': probUA110_Z, 170 'P_UA110_Y': probUA110_Y, 171 'P_UA110_V': probUA110_V 172 } 173 174 ## Code for the sensitivity analysis 175 names = ['B_comp','B_refund','C_price','V_price','Y_price','Z_price'] 176 values = [[1.71083,-0.0398667,-1.67587,0.190499,0.209566,-2.13821],[-0.0398667,0.0188657,-0.00717013,-0.083915,-0.0941582,0.0155518],[-1.67587,-0.00717013,1.76813,0.0330621,0.0365816,2.18927],[0.190499,-0.083915,0.0330621,0.418485,0.452985,-0.0676863],[0.209566,-0.0941582,0.0365816,0.452985,0.498726,-0.0766095],[-2.13821,0.0155518,2.18927,-0.0676863,-0.0766095,2.74714]] 177 vc = bioMatrix(6, names, values) 178 BIOGEME_OBJECT.VARCOVAR = vc 179 180 BIOGEME_OBJECT.SIMULATE = Enumerate(simulate,'obsIter') 181 182 # Statistics 183 nullLoglikelihood(availability,'obsIter') 184 choiceSet = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16] 185 cteLoglikelihood(choiceSet, choice, 'obsIter') 186 availabilityStatistics(availability, 'obsIter') 187 188 # Parameters 189 BIOGEME_OBJECT.PARAMETERS['RandomDistribution'] ="MLHS" 190 BIOGEME_OBJECT.PARAMETERS['NbrOfDraws'] = "1" 191 192 BIOGEME_OBJECT.FORMULAS['AA777 C utility'] = V[1] 193 BIOGEME_OBJECT.FORMULAS['AA777 Z utility'] = V[2] 194 BIOGEME_OBJECT.FORMULAS['AA777 Y utility'] = V[3] 195 BIOGEME_OBJECT.FORMULAS['AA777 V utility'] = V[4] 196 BIOGEME_OBJECT.FORMULAS['AS666 C utility'] = V[5] 197 BIOGEME_OBJECT.FORMULAS['AS666 Z utility'] = V[6] 198 BIOGEME_OBJECT.FORMULAS['AS666 Y utility'] = V[7] 199 BIOGEME_OBJECT.FORMULAS['AS666 V utility'] = V[8] 200 BIOGEME_OBJECT.FORMULAS['DL001 C utility'] = V[9] 201 BIOGEME_OBJECT.FORMULAS['DL001 Z utility'] = V[10] 202 BIOGEME_OBJECT.FORMULAS['DL001 Y utility'] = V[11] 203 BIOGEME_OBJECT.FORMULAS['DL001 V utility'] = V[12] 204 BIOGEME_OBJECT.FORMULAS['UA110 C utility'] = V[13] 205 BIOGEME_OBJECT.FORMULAS['UA110 Z utility'] = V[14] 206 BIOGEME_OBJECT.FORMULAS['UA110 Y utility'] = V[15] 207 BIOGEME_OBJECT.FORMULAS['UA110 V utility'] = V[16]
原理：和估算模型一樣，先加載Python Biogeme使用的包。
然后定義系數。然而，使用已估算的系數。其可從MNL文件夾下的dcm_mnl_param.py文件導入。接着，指定屬性和效用的字典。
在IIA測試中，將排除AA777_C選項，移除完整的模型中選擇這個選項的所有觀測值。你只要指定條件——這里是choice==1——並設置BIOGEME_OBJECT的.EXCLUDE變量。
當然，條件也可能非常復雜。
決定了可用性和排除項后，我們用biologit（...）方法為每個概率創建一個變量，放到模擬的字典中。這個方法以效用的字典作為第一個參數，以可用性的字典作為第二個參數。最后一個參數指定了選項的序數。
最后，我們指定系數的名字和系數的協方差矩陣。這也可從MNL文件夾下的dcm_mnl_param.py文件導入。Python Biogeme將會在敏感度分析中使用。
使用Enumerate（...）方法仿真概率。這個方法以模擬的字典作為第一個參數，以行迭代器作為第二個參數。由於我們只對概率的估算感興趣（而不是蒙特卡洛分析），我們將RandomDistribution指定為MLHS（Modified Latin Hypercube Sampling，修正拉丁超立方抽樣），只畫一次。
和估算一樣運行模擬，輸入一行命令：
/* ../run_pythonbiogeme.sh dcm_mnl_simul ../../../Data/Chapter10/choices.dat */
同樣報錯！
軟件包創建了dcm_mnl_simul.html，包含計算出的概率表格，如下圖所示：
表格有10000行。對於敏感度分析，仿真過程畫出了100個均勻分布的系數參數，返回每個選項的第5個、第95個以及中位數。
最后一步，我們移除第一個選項，重復模擬的過程（TestingIIA目錄下的dcm_iia.py文件和dcm_iia_simul.py文件）：
可以看到，AA777_C選項不在了。如果滿足IIA的要求，我們看到的概率應該不變。我們分析下第二行。
原始模型的P（AA777_V）是0.17025，P（AA777_Y）是0.0555726；比率是3.0636。新模型中，P（AA777_V）=0.172479，P（AA777_Y）=0.0562227；新比率是3.0678。
變化不大。
更多：
然而，這只是一個觀測值。我們檢查一下其余選項（dcm_iia_testing.py文件）：
1 import pandas as pd 2 3 # read the html tables 4 old_model = 'dcm_mnl_simul.html' 5 new_model = 'dcm_iia_simul.html' 6 7 old_model_p = pd.read_html(old_model, header = 0)[3] 8 new_model_p = pd.read_html(new_model, header = 0)[3] 9 10 # let's look at only two columns 11 cols = ['P_AA777_Y', 'P_AA777_V'] 12 13 # make sure that there are no zeros 14 old_model_p = old_model_p[cols] 15 old_model_p = old_model_p[old_model_p[cols[0]] != 0] 16 old_model_p = old_model_p[old_model_p[cols[1]] != 0] 17 18 new_model_p = new_model_p[cols] 19 new_model_p = new_model_p[new_model_p[cols[0]] != 0] 20 new_model_p = new_model_p[new_model_p[cols[1]] != 0] 21 22 # calculate the ratios 23 old_model_p['ratios_old'] = old_model_p \ 24 .apply(lambda row: row['P_AA777_V'] / row['P_AA777_Y'], 25 axis=1) 26 27 new_model_p['ratios_new'] = new_model_p \ 28 .apply(lambda row: row['P_AA777_V'] / row['P_AA777_Y'], 29 axis=1) 30 31 # join with the old model results 32 differences = old_model_p.join(new_model_p, rsuffix='_new') 33 34 # and calculate the differences 35 differences['diff'] = differences\ 36 .apply(lambda row: row['ratios_new'] - row['ratios_old'], 37 axis=1) 38 39 # calculate the descriptive stats for the columns 40 print(differences[['ratios_old', 'ratios_new', 'diff']] \ 41 .describe())
我們使用pandas，因為它提供了從HTML文件直接讀取表格的方法（參考本書1.6節）。
讀入了文件，我們從完整模型和精簡模型中讀取之前測試過的兩列，確保列中沒有0。然后，計算比率並連接兩個文件以計算差異。最后，我們生成描述性統計數據：
差變化不大——變異系數（標准差/平均值）接近於0，最大的差是0.004247。所以，我們不能說違反了IIA。
10.5用巢式Logit模型處理IIA沖突
不過，如果你的模型違反了IIA，那么你就得使用更高級的模型。本技巧中，我們將認識一個稍微復雜點的巢式Logit模型。這個模型將相似的選項放到一個巢里（名字的由來）。篇幅有限，我們這里不討論模型的原理，但是我極力推薦之前提過的Kenneth Train的書。
准備：需要裝好Python Biogeme包。
步驟：
模型代碼的框架還是一樣的；這里，我們只展示不同之處（Nested/dcm_nested.py文件）：
1 from biogeme import * 2 from headers import * 3 from nested import * 4 from loglikelihood import * 5 from statistics import * 6 7 # Specify parameters to be estimated 8 # 9 # Arguments: 10 # - 1 Name, typically, the same as the variable. 11 # - 2 Starting value. 12 # - 3 Lower bound. 13 # - 4 Upper bound. 14 # - 5 flag whether to estimate the parameter (0) 15 # or keep it fixed (1). 16 17 # add the coefficients to be estimated 18 C_price = Beta('C_price',0,-10,10,0,'C price' ) 19 V_price = Beta('V_price',0,-10,10,0,'V price' ) 20 Y_price = Beta('Y_price',0,-10,10,0,'Y price' ) 21 Z_price = Beta('Z_price',0,-10,10,0,'Z price' ) 22 23 ASC = Beta('ASC',0,-10,10,1,'ASC' ) 24 25 B_refund = Beta('B_refund',0,-3,3,0,'refund' ) 26 B_comp = Beta('B_comp',0,-3,3,0,'compartment' ) 27 28 # nest mu parameter 29 biz_mu = Beta('biz_mu',0.5,0,1,0) 30 eco_mu = Beta('eco_mu',0.5,0,1,0) 31 32 # Utility functions 33 34 # compartment attributes 35 c = {} 36 37 c[1] = AA777_1_C_compartment 38 c[2] = AA777_2_Z_compartment 39 c[3] = AA777_3_Y_compartment 40 c[4] = AA777_4_V_compartment 41 c[5] = AS666_1_C_compartment 42 c[6] = AS666_2_Z_compartment 43 c[7] = AS666_3_Y_compartment 44 c[8] = AS666_4_V_compartment 45 c[9] = DL001_1_C_compartment 46 c[10] = DL001_2_Z_compartment 47 c[11] = DL001_3_Y_compartment 48 c[12] = DL001_4_V_compartment 49 c[13] = UA110_1_C_compartment 50 c[14] = UA110_2_Z_compartment 51 c[15] = UA110_3_Y_compartment 52 c[16] = UA110_4_V_compartment 53 54 # price attributes 55 p = {} 56 57 p[1] = AA777_1_C_price 58 p[2] = AA777_2_Z_price 59 p[3] = AA777_3_Y_price 60 p[4] = AA777_4_V_price 61 p[5] = AS666_1_C_price 62 p[6] = AS666_2_Z_price 63 p[7] = AS666_3_Y_price 64 p[8] = AS666_4_V_price 65 p[9] = DL001_1_C_price 66 p[10] = DL001_2_Z_price 67 p[11] = DL001_3_Y_price 68 p[12] = DL001_4_V_price 69 p[13] = UA110_1_C_price 70 p[14] = UA110_2_Z_price 71 p[15] = UA110_3_Y_price 72 p[16] = UA110_4_V_price 73 74 # refund attributes 75 r = {} 76 77 r[1] = AA777_1_C_refund 78 r[2] = AA777_2_Z_refund 79 r[3] = AA777_3_Y_refund 80 r[4] = AA777_4_V_refund 81 r[5] = AS666_1_C_refund 82 r[6] = AS666_2_Z_refund 83 r[7] = AS666_3_Y_refund 84 r[8] = AS666_4_V_refund 85 r[9] = DL001_1_C_refund 86 r[10] = DL001_2_Z_refund 87 r[11] = DL001_3_Y_refund 88 r[12] = DL001_4_V_refund 89 r[13] = UA110_1_C_refund 90 r[14] = UA110_2_Z_refund 91 r[15] = UA110_3_Y_refund 92 r[16] = UA110_4_V_refund 93 94 # The dictionary of utilities is constructed. 95 V = {} 96 97 V[1] = Z_price * p[1] + B_refund * r[1] + B_comp * c[1] 98 V[2] = Z_price * p[2] + B_refund * r[2] + B_comp * c[2] + ASC 99 V[3] = Y_price * p[3] + B_refund * r[3] + B_comp * c[3] 100 V[4] = V_price * p[4] + B_refund * r[4] + B_comp * c[4] 101 V[5] = C_price * p[5] + B_refund * r[5] + B_comp * c[5] 102 V[6] = Z_price * p[6] + B_refund * r[6] + B_comp * c[6] 103 V[7] = Y_price * p[7] + B_refund * r[7] + B_comp * c[7] 104 V[8] = V_price * p[8] + B_refund * r[8] + B_comp * c[8] 105 V[9] = C_price * p[9] + B_refund * r[9] + B_comp * c[9] 106 V[10] = Z_price * p[10] + B_refund * r[10] + B_comp * c[10] 107 V[11] = Y_price * p[11] + B_refund * r[11] + B_comp * c[11] 108 V[12] = V_price * p[12] + B_refund * r[12] + B_comp * c[12] 109 V[13] = C_price * p[13] + B_refund * r[13] + B_comp * c[13] 110 V[14] = Z_price * p[14] + B_refund * r[14] + B_comp * c[14] 111 V[15] = Y_price * p[15] + B_refund * r[15] + B_comp * c[15] 112 V[16] = Y_price * p[16] + B_refund * r[16] + B_comp * c[16] 113 114 # availability flags 115 availability = { 116 1: AA777_1_C_AV, 117 2: AA777_2_Z_AV, 118 3: AA777_3_Y_AV, 119 4: AA777_4_V_AV, 120 5: AS666_1_C_AV, 121 6: AS666_2_Z_AV, 122 7: AS666_3_Y_AV, 123 8: AS666_4_V_AV, 124 9: DL001_1_C_AV, 125 10: DL001_2_Z_AV, 126 11: DL001_3_Y_AV, 127 12: DL001_4_V_AV, 128 13: UA110_1_C_AV, 129 14: UA110_2_Z_AV, 130 15: UA110_3_Y_AV, 131 16: UA110_4_V_AV 132 } 133 134 # 1: nests parameter 135 # 2: list of alternatives 136 business = biz_mu, [1,2,5,6,9,10,13,14] 137 economy = eco_mu, [3,4,7,8,11,12,15,16] 138 nests = business, economy 139 140 # The choice model is a logit, with availability conditions 141 logprob = lognested(V, availability, nests, choice) 142 143 # Defines an itertor on the data 144 rowIterator('obsIter') 145 146 # DEfine the likelihood function for the estimation 147 BIOGEME_OBJECT.ESTIMATE = Sum(logprob,'obsIter') 148 149 # Statistics 150 151 nullLoglikelihood(availability,'obsIter') 152 choiceSet = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16] 153 cteLoglikelihood(choiceSet,choice,'obsIter') 154 availabilityStatistics(availability,'obsIter') 155 156 BIOGEME_OBJECT.PARAMETERS['optimizationAlgorithm'] = 'BIO' 157 BIOGEME_OBJECT.PARAMETERS['checkDerivatives'] = '1' 158 BIOGEME_OBJECT.PARAMETERS['numberOfThreads'] = '8' 159 BIOGEME_OBJECT.PARAMETERS['moreRobustToNumericalIssues'] = '0' 160 161 BIOGEME_OBJECT.FORMULAS['AA777 C utility'] = V[1] 162 BIOGEME_OBJECT.FORMULAS['AA777 Z utility'] = V[2] 163 BIOGEME_OBJECT.FORMULAS['AA777 Y utility'] = V[3] 164 BIOGEME_OBJECT.FORMULAS['AA777 V utility'] = V[4] 165 BIOGEME_OBJECT.FORMULAS['AS666 C utility'] = V[5] 166 BIOGEME_OBJECT.FORMULAS['AS666 Z utility'] = V[6] 167 BIOGEME_OBJECT.FORMULAS['AS666 Y utility'] = V[7] 168 BIOGEME_OBJECT.FORMULAS['AS666 V utility'] = V[8] 169 BIOGEME_OBJECT.FORMULAS['DL001 C utility'] = V[9] 170 BIOGEME_OBJECT.FORMULAS['DL001 Z utility'] = V[10] 171 BIOGEME_OBJECT.FORMULAS['DL001 Y utility'] = V[11] 172 BIOGEME_OBJECT.FORMULAS['DL001 V utility'] = V[12] 173 BIOGEME_OBJECT.FORMULAS['UA110 C utility'] = V[13] 174 BIOGEME_OBJECT.FORMULAS['UA110 Z utility'] = V[14] 175 BIOGEME_OBJECT.FORMULAS['UA110 Y utility'] = V[15] 176 BIOGEME_OBJECT.FORMULAS['UA110 V utility'] = V[16]
原理：前面已經指出，MNL和NL的主要差別在於，NL將有共性的選項放在一起。每個巢關聯一個lambda參數：這個參數定義了巢中選項間競爭有多激烈。lambda參數取值范圍在0（競爭最激烈）和1（沒有競爭）之間。
我們的例子中，我們假設商務選項（艙位C和Z）有更大的休息室和更好的服務，所以將他們放到商務巢中；經濟選項（艙位Y和V）放到經濟巢中。biz_mu和eco_mu是我們巢的lambda參數。
然后我們傳入帶選項數字的列表以及lambda參數，將選項分配到巢。
我們也要使用一個不同的估算器——lognested（...）方法，而不是bioLogLogit（...）方法；這個額外加入了nests對象：
果然，巢的結構並沒有給模型帶來什么好處。Rho平方值略好一點，但付出的代價是兩個額外的自由度。另外，biz_mu和biz_eco也不是統計上顯著的。
10.6用混合Logit模型處理復雜的替代模式
混合Logit模型，與前面介紹的模型都不同，允許有些系數隨機地服從一個正態分布，也就是說，有均值和標准差。這樣的效果是降低了對IIA假設的依賴，也允許靈活地給替代模式建模。不過，這也要付出計算時間的代價。
准備：需要裝好Python Biogeme包。
步驟：由於我們已經明確了，用我們的數據集估算出來的MNL模型，並不違反IIA性質，我們就只展示一下估算混合Logit模型的機制（MixedLogit/dcm_mixed.py文件）：
1 from biogeme import * 2 from headers import * 3 from loglikelihood import * 4 from statistics import * 5 6 # Specify parameters to be estimated 7 # 8 # Arguments: 9 # - 1 Name, typically, the same as the variable. 10 # - 2 Starting value. 11 # - 3 Lower bound. 12 # - 4 Upper bound. 13 # - 5 flag whether to estimate the parameter (0) 14 # or keep it fixed (1). 15 16 # add the coefficients to be estimated 17 18 C_price = Beta('C_price',0,-10,10,0,'C price' ) 19 V_price = Beta('V_price',0,-10,10,0,'V price' ) 20 Y_price = Beta('Y_price',0,-10,10,0,'Y price' ) 21 Z_price = Beta('Z_price',0,-10,10,0,'Z price' ) 22 23 ASC = Beta('ASC',0,-10,10,1,'ASC' ) 24 25 B_ref = Beta('B_ref',0,-3,3,0,'refund' ) 26 B_ref_S = Beta('B_ref_S',0,-3,3,0,'refund (std)' ) 27 B_comp = Beta('B_comp',0,-3,3,0,'compartment' ) 28 29 # Random parameters 30 B_ref_rnd = B_ref + B_ref_S * bioDraws('B_ref_rnd') 31 32 # Utility functions 33 34 # The data are associated with the alternative index 35 # compartment attributes 36 c = {} 37 38 c[1] = AA777_1_C_compartment 39 c[2] = AA777_2_Z_compartment 40 c[3] = AA777_3_Y_compartment 41 c[4] = AA777_4_V_compartment 42 c[5] = AS666_1_C_compartment 43 c[6] = AS666_2_Z_compartment 44 c[7] = AS666_3_Y_compartment 45 c[8] = AS666_4_V_compartment 46 c[9] = DL001_1_C_compartment 47 c[10] = DL001_2_Z_compartment 48 c[11] = DL001_3_Y_compartment 49 c[12] = DL001_4_V_compartment 50 c[13] = UA110_1_C_compartment 51 c[14] = UA110_2_Z_compartment 52 c[15] = UA110_3_Y_compartment 53 c[16] = UA110_4_V_compartment 54 55 # price attributes 56 p = {} 57 58 p[1] = AA777_1_C_price 59 p[2] = AA777_2_Z_price 60 p[3] = AA777_3_Y_price 61 p[4] = AA777_4_V_price 62 p[5] = AS666_1_C_price 63 p[6] = AS666_2_Z_price 64 p[7] = AS666_3_Y_price 65 p[8] = AS666_4_V_price 66 p[9] = DL001_1_C_price 67 p[10] = DL001_2_Z_price 68 p[11] = DL001_3_Y_price 69 p[12] = DL001_4_V_price 70 p[13] = UA110_1_C_price 71 p[14] = UA110_2_Z_price 72 p[15] = UA110_3_Y_price 73 p[16] = UA110_4_V_price 74 75 # refund attributes 76 r = {} 77 78 r[1] = AA777_1_C_refund 79 r[2] = AA777_2_Z_refund 80 r[3] = AA777_3_Y_refund 81 r[4] = AA777_4_V_refund 82 r[5] = AS666_1_C_refund 83 r[6] = AS666_2_Z_refund 84 r[7] = AS666_3_Y_refund 85 r[8] = AS666_4_V_refund 86 r[9] = DL001_1_C_refund 87 r[10] = DL001_2_Z_refund 88 r[11] = DL001_3_Y_refund 89 r[12] = DL001_4_V_refund 90 r[13] = UA110_1_C_refund 91 r[14] = UA110_2_Z_refund 92 r[15] = UA110_3_Y_refund 93 r[16] = UA110_4_V_refund 94 95 # The dictionary of utilities is constructed. 96 V = {} 97 98 V[1] = Z_price * p[1] + B_ref_rnd * r[1] + B_comp * c[1] 99 V[2] = Z_price * p[2] + B_ref_rnd * r[2] + B_comp * c[2] + ASC 100 V[3] = Y_price * p[3] + B_ref_rnd * r[3] + B_comp * c[3] 101 V[4] = V_price * p[4] + B_ref_rnd * r[4] + B_comp * c[4] 102 V[5] = C_price * p[5] + B_ref_rnd * r[5] + B_comp * c[5] 103 V[6] = Z_price * p[6] + B_ref_rnd * r[6] + B_comp * c[6] 104 V[7] = Y_price * p[7] + B_ref_rnd * r[7] + B_comp * c[7] 105 V[8] = V_price * p[8] + B_ref_rnd * r[8] + B_comp * c[8] 106 V[9] = C_price * p[9] + B_ref_rnd * r[9] + B_comp * c[9] 107 V[10] = Z_price * p[10] + B_ref_rnd * r[10] + B_comp * c[10] 108 V[11] = Y_price * p[11] + B_ref_rnd * r[11] + B_comp * c[11] 109 V[12] = V_price * p[12] + B_ref_rnd * r[12] + B_comp * c[12] 110 V[13] = C_price * p[13] + B_ref_rnd * r[13] + B_comp * c[13] 111 V[14] = Z_price * p[14] + B_ref_rnd * r[14] + B_comp * c[14] 112 V[15] = Y_price * p[15] + B_ref_rnd * r[15] + B_comp * c[15] 113 V[16] = Y_price * p[16] + B_ref_rnd * r[16] + B_comp * c[16] 114 115 # availability flags 116 availability = { 117 1: AA777_1_C_AV, 118 2: AA777_2_Z_AV, 119 3: AA777_3_Y_AV, 120 4: AA777_4_V_AV, 121 5: AS666_1_C_AV, 122 6: AS666_2_Z_AV, 123 7: AS666_3_Y_AV, 124 8: AS666_4_V_AV, 125 9: DL001_1_C_AV, 126 10: DL001_2_Z_AV, 127 11: DL001_3_Y_AV, 128 12: DL001_4_V_AV, 129 13: UA110_1_C_AV, 130 14: UA110_2_Z_AV, 131 15: UA110_3_Y_AV, 132 16: UA110_4_V_AV 133 } 134 135 # The choice model is a logit, with availability conditions 136 prob = bioLogit(V, availability, choice) 137 l = mixedloglikelihood(prob) 138 139 # Defines an itertor on the data 140 rowIterator('obsIter') 141 142 # Define the likelihood function for the estimation 143 BIOGEME_OBJECT.ESTIMATE = Sum(l, 'obsIter') 144 145 # Statistics 146 147 nullLoglikelihood(availability,'obsIter') 148 choiceSet = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16] 149 cteLoglikelihood(choiceSet, choice, 'obsIter') 150 availabilityStatistics(availability, 'obsIter') 151 152 BIOGEME_OBJECT.PARAMETERS['NbrOfDraws'] = '100' 153 BIOGEME_OBJECT.DRAWS = { 'B_ref_rnd': 'NORMAL' } 154 BIOGEME_OBJECT.PARAMETERS['optimizationAlgorithm'] = 'BIO' 155 BIOGEME_OBJECT.PARAMETERS['numberOfThreads'] = '8' 156 157 BIOGEME_OBJECT.FORMULAS['AA777 C utility'] = V[1] 158 BIOGEME_OBJECT.FORMULAS['AA777 Z utility'] = V[2] 159 BIOGEME_OBJECT.FORMULAS['AA777 Y utility'] = V[3] 160 BIOGEME_OBJECT.FORMULAS['AA777 V utility'] = V[4] 161 BIOGEME_OBJECT.FORMULAS['AS666 C utility'] = V[5] 162 BIOGEME_OBJECT.FORMULAS['AS666 Z utility'] = V[6] 163 BIOGEME_OBJECT.FORMULAS['AS666 Y utility'] = V[7] 164 BIOGEME_OBJECT.FORMULAS['AS666 V utility'] = V[8] 165 BIOGEME_OBJECT.FORMULAS['DL001 C utility'] = V[9] 166 BIOGEME_OBJECT.FORMULAS['DL001 Z utility'] = V[10] 167 BIOGEME_OBJECT.FORMULAS['DL001 Y utility'] = V[11] 168 BIOGEME_OBJECT.FORMULAS['DL001 V utility'] = V[12] 169 BIOGEME_OBJECT.FORMULAS['UA110 C utility'] = V[13] 170 BIOGEME_OBJECT.FORMULAS['UA110 Z utility'] = V[14] 171 BIOGEME_OBJECT.FORMULAS['UA110 Y utility'] = V[15] 172 BIOGEME_OBJECT.FORMULAS['UA110 V utility'] = V[16]
原理：隨機的B_ref_rnd參數由確定的部分（平均值）B_ref和變化的部分（隨機值）B_ref_S——即標准差——組成。
bioDraws（...）用於任意抽選。抽選是針對每一個觀測值進行的，這就帶來了性能問題：對於每個觀測值，估算器都要進行正態分布預定次數的抽樣，以估算系數。
之后每個效用函數都使用random參數取代B_ref；這么做可以處理替代模式，因為模型可以在這個維度上處理選項之間的相關性。
mixedloglikelihood（...）方法進行模型的預估。它唯一的參數是bioLogit（...）對象。
設置.DRAWS變量，指定從什么分布中抽；我們的例子里是正態分布。結果如下：
如我們所料，附加的估算系數並沒有給模型帶來價值；B_ref_S並不顯著，完全可以從模型中移除，也不會損失精確度。
同之前的模型相比，你可以發現ML估算的時間比較久：MNL要花1秒，而NL要花50秒。
第10章完。