《數據分析實戰-托馬茲.卓巴斯》讀書筆記第10章--離散選擇模型理論
python學習筆記-目錄索引html
第10章解釋了選擇模型理論以及一些流行的模型：多項式Logit模型、嵌套Logit模型以及混合Logit模型。python
本章中，會學習如下技巧：
·準備數據集以估算離散選擇模型
·估算知名的多項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上的指南安裝。裝好後就能夠經過下面的命令調用了：ios
>>pythonbiogeme <model_file> <dataset>
注意：這段話來自原書做者git
/*
咱們使用軟件包時遇到了問題。這個軟件包要求Python 3.2之上的版本。儘管本書使用的是Python 3.5（安裝在~/anaconda文件夾中），
可是Python Biogeme會到/usr/local/Frameworks/Python.framework/Versions/3.5路徑下尋找Python，這個路徑在咱們的機器上根本不存在。
咱們只好手動建立這個路徑，而後從新安裝軟件。 咱們也在Windows上試了，能夠正常運行。不過，有些朋友可能會在安裝和運行時遇到問題，這也是要你們注意一下的。 若是你遇到了任何問題，你能夠在Biogeme用戶羣中詢問；通常很快就有回覆了。若是你實在無法子了，也能夠給原做者發郵件，mail@tomdrabas.com，能幫的原做者必定幫。 */
如下來自邀月的安裝步驟：web
/* 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文件格式以下：json
/* 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章完。
- 1. 《數據分析實戰-托馬茲.卓巴斯》讀書筆記第6章-迴歸模型
- 2. 《數據分析實戰-托馬茲.卓巴斯》讀書筆記第1章-數據格式與數據交互
- 3. 《數據分析實戰-托馬茲.卓巴斯》讀書筆記第4章-聚類技巧（K均值、BIRCH、DBSCAN）
- 4. 《數據分析實戰-托馬茲.卓巴斯》讀書筆記第8章--圖(NetworkX、Gephi)修訂版
- 5. 《數據分析實戰-托馬茲.卓巴斯》讀書筆記第5章-降維技巧
- 6. 《數據分析實戰-托馬茲.卓巴斯》讀書筆記第7章-時間序列技術（ARMA模型、ARIMA模型）
- 7. 《數據分析實戰-托馬茲.卓巴斯》讀書筆記第2章-變量分佈與相關性、圖表
- 8. 大數據導論托馬斯—第一章理解大數據
- 9. jQuery 實戰讀書筆記之第二章：選擇元素
- 10. 大數據之路、阿里巴巴大數據實踐讀書筆記 --- 第四章、離線數據開發
- 更多相關文章...
- • SQLite 分離數據庫 - SQLite教程
- • MySQL選擇數據庫（MySQL USE語句） - MySQL教程
- • 委託模式
- • Tomcat學習筆記（史上最全tomcat學習筆記）
-
每一个你不满意的现在，都有一个你没有努力的曾经。
- 1. 《數據分析實戰-托馬茲.卓巴斯》讀書筆記第6章-迴歸模型
- 2. 《數據分析實戰-托馬茲.卓巴斯》讀書筆記第1章-數據格式與數據交互
- 3. 《數據分析實戰-托馬茲.卓巴斯》讀書筆記第4章-聚類技巧（K均值、BIRCH、DBSCAN）
- 4. 《數據分析實戰-托馬茲.卓巴斯》讀書筆記第8章--圖(NetworkX、Gephi)修訂版
- 5. 《數據分析實戰-托馬茲.卓巴斯》讀書筆記第5章-降維技巧
- 6. 《數據分析實戰-托馬茲.卓巴斯》讀書筆記第7章-時間序列技術（ARMA模型、ARIMA模型）
- 7. 《數據分析實戰-托馬茲.卓巴斯》讀書筆記第2章-變量分佈與相關性、圖表
- 8. 大數據導論托馬斯—第一章理解大數據
- 9. jQuery 實戰讀書筆記之第二章：選擇元素
- 10. 大數據之路、阿里巴巴大數據實踐讀書筆記 --- 第四章、離線數據開發