基礎ゼミ5
Rによる多項ロジットモデル・ロジットモデル推定
2 0 1 7/05/16
福 田 研 修 士 1年鈴木新目次
前回の復習多項ロジットモデルの推定ネスティッドロジットモデルの推定
他のモデルに関して
付録
参考文献

2

前回の復習
選択行動のモデル化（東京-大阪間の交通手段の例）

観測要因

𝑥13620円
2.5時間
座席間隔104cm

非観測要因

𝜀25490円
1時間
座席間隔69cm

年収
移動に費やせる時間
鉄オタかどうか𝑈（選択行動） = 𝑉(𝛽𝑥) + 𝜀

選択行動を関数で表す(効用関数)
3

前回の復習離散選択モデル
個人𝑛は、選択肢集合𝐶𝑛 の中から、最大の効用𝑈𝑖𝑛 を持つ選択肢を選ぶ

𝑈𝑖𝑛 = 𝑉𝑖𝑛 + 𝜀𝑖𝑛
確定項

ランダム項
（確率変数）

非観測要因（ランダム項）は観測できないため，代わりに確率変数を用いる．
この確率変数がガンベル分布に従うモデルがロジットモデル

個人は選択肢集合の中から効用が最大となる選択肢を選ぶ4

前回の復習最尤推定法・・・尤度を最大化するパラメーター値を探索する方法
例題
表が出る確率が θ であるようなコインがある。このコインを 100回投げたら 75 回表が出た。最尤法により θ を推定せよ。
コインを 100 回投げたら 75 回表が出る確率は𝐿 𝜃 =100 𝐶75 𝜃 75 (1 − 𝜃)25 （尤度関数）
𝐿 𝜃 を最大にする𝜃を考える．（計算を簡単にするため対数を取る）

𝑙𝑜𝑔𝐿 𝜃 = 75𝑙𝑜𝑔𝜃 + 25𝑙𝑜𝑔(1 − 𝜃) + 𝑙𝑜𝑔100 𝐶75・・・（微分とかすると）・・・

𝜃=0.75の時𝐿 𝜃 最大，尤もらしい値

𝜃=0.01でも表75回出る可能性あるけど確率低い
5

多項ロジットモデルの推定（MNL）
ボーイング社における乗客の航空路選択モデル
Boeing社の路線選択調査データ，2004年

Option1（直行便），Option2（他空港立ち寄り，同機体），Option3（他空港立ち寄り，乗り換え）

最も魅力的な路線を１つ選んでもらう

6

多項ロジットモデルの推定（MNL）●データの項目
項目

内容

項目

内容

Subj_Age

年齢

SP1_MostAttractive

選択したルート

Subj_Occupation

職業

OptX_Fare

料金（$）

Subj_Income

収入（$）

OptX_TotalTriptime

移動時間(hours)

Trip_Purpose

移動目的

Trip_PartySize

同行者数

OptX_Legroom

座席間隔

Trip_Miles

移動距離(miles)

今回，Option1, Option2, Option3の選択において，それぞれのルートの効用関数を考え，各項目のパラメータを推定する．

7

多項ロジットモデルの推定(MNL)
今回，「移動時間」，「料金」，「座席間隔」の項目を使用す
ると，Option1, Option2, Option3の各効用関数𝑉𝑜𝑝 は次のようになる．
𝑉𝑜𝑝1 = 𝐴𝐶𝑆𝑜𝑝1 + 𝛽𝑡𝑖𝑚𝑒 ∙ 𝑥𝑜𝑝1_𝑡𝑖𝑚𝑒 +𝛽𝑓𝑎𝑟𝑒 ∙ 𝑥𝑜𝑝1_𝑓𝑎𝑟𝑒 +𝛽𝑙𝑒𝑔 ∙ 𝑥𝑜𝑝1_𝑙𝑒𝑔
𝑉𝑜𝑝2 = 𝐴𝐶𝑆𝑜𝑝2 + 𝛽𝑡𝑖𝑚𝑒 ∙ 𝑥𝑜𝑝2_𝑡𝑖𝑚𝑒 +𝛽𝑓𝑎𝑟𝑒 ∙ 𝑥𝑜𝑝2_𝑓𝑎𝑟𝑒 +𝛽𝑙𝑒𝑔 ∙ 𝑥𝑜𝑝2_𝑙𝑒𝑔
𝑉𝑜𝑝3 =𝛽𝑡𝑖𝑚𝑒 ∙ 𝑥𝑜𝑝3_𝑡𝑖𝑚𝑒 +𝛽𝑓𝑎𝑟𝑒 ∙ 𝑥𝑜𝑝3_𝑓𝑎𝑟𝑒 +𝛽𝑙𝑒𝑔 ∙ 𝑥𝑜𝑝3_𝑙𝑒𝑔

option1 option2 option3
直行便 transit transit（同機） （異機）

𝐴𝐶𝑆：選択肢固有定数（定数項）𝛽𝑡𝑖𝑚𝑒 ：所要時間に対するパラメータ
𝛽𝑓𝑎𝑟𝑒 : 料金に対するパラメータ𝛽𝑙𝑒𝑔 : 座席間隔に対するパラメータ(=快適さの指標)8データの読み込み

Dataの行数を調べる

要素が0のベクトルを作成する
c(0,0,0,0,0)9

尤度関数の定義，初期値の設定
Function関数で尤度関数を定義する※function関数の例
a <- function(x){x+5}a(2)でまわすと7が出る

変数の定義
b1=x[1]=𝐴𝐶𝑆𝑜𝑝1b2=x[2]=𝐴𝐶𝑆𝑜𝑝2
t1 =x[3]=𝑥𝑡𝑖𝑚𝑒

f1=x[4]=𝑥𝑓𝑎𝑟𝑒
L1=x[5]=𝑥𝑙𝑒𝑔

尤度をLLと定義し，初期値を0に10効用関数の定義
𝑉𝑜𝑝1 = 𝐴𝐶𝑆𝑜𝑝1 + 𝛽𝑜𝑝1_𝑡𝑖𝑚𝑒 ∙ 𝑥𝑡𝑖𝑚𝑒 +𝛽𝑜𝑝1_𝑓𝑎𝑟𝑒 ∙ 𝑥𝑓𝑎𝑟𝑒 +𝛽𝑜𝑝1_𝑙𝑒𝑔 ∙ 𝑥𝑙𝑒𝑔
𝑉𝑜𝑝2 = 𝐴𝐶𝑆𝑜𝑝2 + 𝛽𝑜𝑝2_𝑡𝑖𝑚𝑒 ∙ 𝑥𝑡𝑖𝑚𝑒 +𝛽𝑜𝑝2_𝑓𝑎𝑟𝑒 ∙ 𝑥𝑓𝑎𝑟𝑒 +𝛽𝑜𝑝2_𝑙𝑒𝑔 ∙ 𝑥𝑙𝑒𝑔
𝑉𝑜𝑝3 =𝛽𝑜𝑝3_𝑡𝑖𝑚𝑒 ∙ 𝑥𝑡𝑖𝑚𝑒 +𝛽𝑜𝑝3_𝑓𝑎𝑟𝑒 ∙ 𝑥𝑓𝑎𝑟𝑒 +𝛽𝑜𝑝3_𝑙𝑒𝑔 ∙ 𝑥𝑙𝑒𝑔
移動時間

座席間隔

料金

定数項

サンプル数×1のマトリックス

料金については，説明変数のおよそのケタを揃えて推定計算を安定させるため，100で割り100$を単位として変数をとる．
11

選択確率の導出および選択結果の変換
各Optionの選択確率，例えばOption1の場合は
𝑒 𝑉𝑜𝑝1
𝑃𝑜𝑝1 = 𝑉𝑜𝑝1𝑒+ 𝑒 𝑉𝑜𝑝2 + 𝑒 𝑉𝑜𝑝3

選択結果の数値変換

データのSP1_MostAttractive（選択結果）の列が「Option3」の場合は1を，そうでない行に関しては0を出力する．
12

対数尤度関数の設定
ルート3を選択した場合，ルート3の選択確率が算出される
ex) Option1を選択したデータの場合→ 1*log( 𝑉𝑜𝑝1𝑒𝑒 𝑉𝑜𝑝1+𝑒 𝑉𝑜𝑝2 +𝑒 𝑉𝑜𝑝3)

その和をとることで尤度関数を算出

尤度関数の定義は終了，この関数を最大化する．

13対数尤度関数の最大化

Optim・・・最適化関数
対数尤度関数（fr）を最適化（尤度を最大化）する
パラメーターを求める$parパラメーターの最適解の値$value最大尤度$message
エラーのありなし，NULLなら無し$hessian
ヘッセ行列（t値の計算に必要）
14

t値の算出solve・・・ 逆行列を算出する．
ヘッセ行列を分散共分散行列に
diag ・・・行列の作成sqrt ・・・平方根の計算

標準偏差の算出

最尤推定量を標準偏差で割ることでt値を算出
t値が1.96以上なら有意水準5%を満たす
15ネスティッドロジットモデルの推定(NL)MNLの限界（IIA特性）→NL（ネスティッドロジットモデル）
今回は直行便かTransit便かで分ける．
プログラム大体MNLと同じ
パラメーターの初期設定や尤度関数など一部変更あり

option1 option2 option3
直行便 transit transit（同機） （異機）

16

ネスティッドロジットモデルの推定(NL)

スケールパラメーター分1個増やす（初期値は1に）関数にもスケールパラメーターを定義

17ネスティッドロジットモデルの推定(NL)
ツリー上部の選択確率の作成(直行便かどうか，ログサム変数)

1
𝑉𝑖 = ln(𝛽𝑖𝑃(𝑚) =𝑖∈𝐶exp(𝛽𝑖 𝑉𝑖 ))exp(𝜇𝑚 𝑉𝑙 )𝑙 exp(𝜇𝑚 𝑉𝑙 )i:op1, op2, op3m:direct, transit

ツリー下部においても同様，条件付き確率を算出（transitの場合同機か異機か）

exp(𝜇𝑚 𝑉𝑖 )𝑃 𝑖𝑚 =
𝑗 exp(𝜇𝑚 𝑉𝑖 )18

ネスティッドロジットモデルの推定(NL)

同時確率を算出
exp(𝜇𝑚 𝑉𝑖 ) exp(𝜇𝑚 𝑉𝑙 )𝑃 𝑖, 𝑚 = 𝑃 𝑖|𝑚 𝑃 𝑚 =𝑗 exp(𝜇𝑚 𝑉𝑖 ) 𝑙 exp(𝜇𝑚 𝑉𝑙 )MethodをBFGSに

19

他のモデルに関して（CNL，MXL）

下記のサイトにコードが載っているので，参考にしてください．
http://bin.t.u-tokyo.ac.jp/kaken

20

付録：データの補足説明

21