23rd Sep 2022
行動モデル夏の学校

MNL, NL, RL コードの説明
小林里瑳

Multinomial Logit ModelMNLモデルのパラメータ推定

1. まずはデータを読み込む
2. データ数を数える
3. パラメータ数の初期値(=0)を代入したベクトル を作る
今回はパラメータ数を5にするので、要素が5 の0ベクトルを作る

2

MNLモデルのパラメータ推定
ベクトルxを入れると尤度を出す関数frを作る
今回の効用関数
𝑉!"#$% = 𝑏1 + 𝑑1 ∗ 鉄道所要時間/100
𝑉&'( = 𝑏2 + 𝑑1 ∗ バス所要時間/100
𝑉)#" = 𝑏3 + 𝑑1 ∗ 自動車所要時間/100
𝑉&$*+ = 𝑏4 + 𝑑1 ∗ 自転車所要時間/100
𝑉,#-* =+𝑑1 ∗ 徒歩所要時間/100

MNLの確率exp(𝜇𝑉! )𝑃 𝑖 =∑"∈$ exp(𝜇𝑉! )3

MNLモデルのパラメータ推定
exp(𝜇𝑉! )𝑃 𝑖 =∑"∈$ exp(𝜇𝑉! )?

?

exp(𝜇𝑉! )𝑃 𝑖 =∑"∈$ exp(𝜇𝑉! )選択確率が0だと尤度も0になる
するとln 0 = −∞ になってしまいエ
ラーが出る
選択確率0の時，1にする処理をする(𝐿𝐿 𝜽 = - - 𝑦!% log 𝑃% (𝑖|𝜽)
?

%&' !
4

MNLモデルのパラメータ推定

l 最適化関数optim
l optim(par, fn, gr = NULL, method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B","SANN"), lower = -Inf, upper = Inf, control = list(), hessian = FALSE, ...)l Par というパラメータを初期値としてfnを最小化するようにパラメータを動かしながら反復して探索する（最大化の場合はfnscale =-1）
l 探索方法をmethod から選ぶ-“Nelder-Mead”法：関数値だけを用い頑健(初期値の選択に敏感)だが遅い-“BFGS”法：準ニュートン法．関数値と勾配関数を関数の曲面近似に使う-“CG”法：共役勾配法．破綻しやすいがメモリ使用量が少ない-“L-BFGS-B”法：変数の上限・下限を設定した準ニュートン法-“SANN”法：確率的手法である焼きなまし法．関数値だけを用い遅い．
5

MNLモデルのパラメータ推定
l t値
l 母集団の平均値の仮説検定に利用
l |t| >1.96 で5%有意
l ヘッセ行列の逆行列の対角成分の平方根が分母

l 修正済み尤度比l)*+())+パラメータ数))*

6

Nested Logit Model𝑚𝑖

NLモデルのパラメータ推定
MNL同様にベクトルxを入れると尤度を出す関数frを作る
x[1]~x[4] は定数項（MNLと同じ）

x[5] は時間：上位ネストと下位ネストの組
み合わせで決まる効用x[6] は運賃：上位ネストと下位ネストの組
み合わせで決まる効用（今回は考慮しない）
x[7] はスケールパラメータ

NLではスケールパラメータを推定
スケールパラメータが1ならばMNLと同等
8

NLモデルのパラメータ推定
l 各選択肢の効用の定義

l 今回はf1で公共交通の運賃を定義しています9

NLモデルのパラメータ推定
ネストの作成：上位ネストの計算

公共／私有の選択についてログサム変数を設定

ログサム変数を使って公共／私有の選択確率の分子を計算

分母を計算し(deno)，公共／私有の選択確率を計算
10

NLモデルのパラメータ推定
ネストの作成：上位ネスト内の条件付確率公共ー鉄道／バス というネスト構造
公共を条件とした鉄道, 公共を条件としたバスの選択確率をそれぞれ計算

私有ー車／自転車／徒歩 というネスト構造
私有を条件とした車，自転車，徒歩の選択確率をそれぞれ計算
11

NLモデルのパラメータ推定
同時確率の計算ネストの作り方を含めたNLモデルの詳細については2日目の倉内先生の資料を参照ください

12

NLモデルのパラメータ推定
l 最適化関数optim による最尤推定l ここでは7番目のパラメータがスケールパラメータなので，変数に下限０上限１を設定して，準ニュートン法でパラメータを探索
l スケールパラメータが(0,1)の範囲に収まらない場合はネスト構造が不適切
l 同様に１に貼り付いている場合は不適切であると同時にMNLと同義

13

Recursive Logit Model
Fosgerau, Mogens, Emma Frejinger, and Anders Karlstrom. "A link based
network route choice model with unrestricted choice set." Transportation
Research Part B: Methodological 56 (2013): 70-80.

はじめに
l 自宅から最寄り駅に行くルートをランダム効用最大化理論で考える
l MNL モデルでは経路の羅列が不可能（組み合わせ爆発を起こすため）出発

l 解決方法

到着l 代表的なルートのみ抽出する：明示的な経路集合を定める
l 逐次的な行動モデルを考えることで明示的な経路列挙を必要としない

15

Recursive Logit モデルとは
l 即時効用，誤差項，目的地までの期待最大効用に従って各ノードで選択
l 選択そのものはMNL モデル

l 逐次的かつ再帰的な方法で経路が選択される
l 逐次的：リンク選択の積み重ねが経路全体の選択になると考える
l 再帰的：次の状態だけでなく将来の状態まで考えた効用の式...Bellman 方程式の利用
l マルコフ過程に従うと仮定する

l 選択肢は無限集合であり，ループを含む経路を選択する可能性がある
l 道路ネットワークが与えられている場合，経路集合を設定しなくとも良い-アクティビティパス選択などネットワークを限定する必要がある場合はまた別

16

定式化
l 有向連結グラフ（ネットワーク）𝒢 = 𝒮, 𝒜 を考える
l 𝒮: ノードの集合，𝒜: リンクの集合
l リンク 𝑘, 𝑎 ∈ 𝒜

l リンク𝑘 の終点から出るリンクの集合 𝒜 𝑘 , 𝑎 ∈ 𝒜 𝑘 ∈ 𝒜

l OからD の経路はリンクの列 𝑘* … 𝑘. で表される ∀𝑖 < 𝐼, 𝑘!/' ∈ 𝒜 𝑘! ,17

定式化
l 効用関数l 即時効用𝑢% 𝑎 𝑘 = 𝑣% 𝑎 𝑘 + 𝜇𝜀% 𝑎l 期待効用

𝑉%0 𝑎

l 終点のダミーリンク 𝑑l 目的地ノード(D) からダミーリンク 𝑑 を追加し，吸収状態を定義
l 全リンクの集合を𝒜8 . = 𝒜 ∪ 𝑑l 目的地ノードを終点とする全てのリンク𝑘に対して𝑣% 𝑑 𝑘 = 0 を仮定する

l 主体はマルコフ過程に従って次の状態を選択する
l 各リンク𝑘 でランダム効用𝜀% 𝑎 を観測
l 即時効用と期待効用の和で表される効用を最大化する

𝑣% 𝑎 𝑘 + 𝜇𝜀% 𝑎 + 𝑉%. 𝑎18

定式化
l 効用関数l 即時効用𝑢% 𝑎 𝑘 = 𝑣% 𝑎 𝑘 + 𝜇𝜀% 𝑎- 旅行者 𝑛 がリンク𝑘 にいる状態で選択肢集合𝒜 𝑘 からリンク𝑎 を選択する- 確定項 𝑣% 𝑎 𝑘 = 𝑣% (𝑥%,#|* ; 𝛽)
- 誤差項 𝜇𝜀% 𝑎 ：平均値0のガンベル分布を仮定l 期待効用

𝑉%0 𝑎- リンク𝑎 を選択した場合の目的地Dから伸びるダミーリンク𝑑までの期待効用

19

定式化
l 効用関数l 即時効用𝑢% 𝑎 𝑘 = 𝑣% 𝑎 𝑘 + 𝜇𝜀% 𝑎l 期待効用𝑉%0 𝑎l 価値関数 𝑉%0 𝑘 : Bellman 方程式と呼ぶ
𝑉%0 𝑘 = E
リンク𝑘 を選んだ場合の目的地𝑑 までの期待効用

max

1∈𝒜 3

リンク 𝑘 において次にリンク𝑎を選んだ場合の効用

𝑣% 𝑎 𝑘 + 𝑉%0 𝑎 + 𝜇𝜀% 𝑎∀𝑘 ∈ 𝒜

𝒜 𝑘 の中から効用が最大のものを選ぶ

20

選択確率
l リンク 𝑘 における次のリンク 𝑎 の選択確率はMNL モデルで表される
exp
𝑃%. 𝑎 𝑘 =1𝑣% (𝑎|𝑘) + 𝑉%. 𝑎𝜇(RL 3)

1∑#! ∈𝒜 * exp
𝑣% (𝑎′|𝑘) + 𝑉%. 𝑎′𝜇l 価値関数は次のように書くことができる
1. 𝑎 )
𝜇lnJ
𝛿𝑎𝑘exp(𝑣𝑎𝑘+𝑉%%
𝜇𝑉%. 𝑘 =∀𝑘 ∈ 𝒜#∈𝒜(*)

0𝑘=𝑑𝛿 𝑎𝑘 = M(RL 4)

1 𝑖𝑓 𝑎 ∈ 𝒜(𝑘)0 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒21価値関数の求め方
l Bellman 方程式を解く 𝑉%0 𝑘 = Emax1∈𝒜 3

𝑣% 𝑎 𝑘 + 𝑉%0 𝑎 + 𝜇𝜀% 𝑎l 式(RL 4) の両辺に対数をとって変換する' %5 3𝑒4=-'% 1 )(6𝑎𝑘/5𝛿 𝑎 𝑘 𝑒4

∀𝑘 ∈ 𝒜1∈𝒜(3)

1(RL 5)𝑘=𝑑l さらにこれを行列式で定義する
l 以下に示す要素𝑀 から構成される行列 𝑴 ( 𝒜 × 𝒜 ) を定義する
5(7 𝑎 𝑘 )
𝛿 𝑎 𝑘 𝑒6

𝑀*# = Z0𝑎 ∈ 𝒜(𝑘)𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒(RL 6)

22

価値関数の求め方l Bellman 方程式を解く 𝑉%0 𝑘 = Emax

1∈𝒜 3𝑣% 𝑎 𝑘 + 𝑉%0 𝑎 + 𝜇𝜀% 𝑎l さらに以下の要素からなるベクトル𝒃 ( 𝒜 ×1)と𝒛 ( 𝒜 ×1)を定義する

1 𝑘≠𝑑𝑏3 = K
0 𝑘=𝑑𝑧3' %5 3= 𝑒4

l すると式(RL 5) は以下のように簡潔な線形方程式で書き下すことができる

𝐳 = 𝐌𝐳 + 𝐛 ⟺ 𝐈 − 𝐌 𝐳 = 𝐛

𝐈 は単位行列(RL 7)

l 𝐳 が解を持つ条件は𝐈 − 𝐌 が逆行列を持つこと-逆行列を持たない可能性もある．持つかどうかはパスの数や瞬間効用のサイズに依存する

" !
$ %
#𝑒=$"
(+ 𝑎 𝑘 ,$ ! & )#𝛿 𝑎𝑘 𝑒∀𝑘 ∈ 𝒜&∈𝒜(%)

1(RL 5)𝑘=𝑑23

遷移確率の求め方
l 選択確率に対応する行列𝑷 を作る
l 目的地に対するリンクの選択確率は出発点に依存しないことに注意
l 式(3) を𝒜 内のリンクについて定義した行列𝑷 に整理する
l 状態 𝑘 に対応する行は以下の通り． ⊙は要素ごとの積， 𝑴* は𝐌の𝑘 列目を表す

𝑴 b ⊙ 𝒛c𝑷b =𝑴b 𝒛

(RL 8)

行列 𝑴 の要素'(6 𝑎 𝑘 )4
𝛿 𝑎𝑘 𝑒𝑀31 = Q0𝑎 ∈ 𝒜(𝑘)𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒(RL 6)

24

経路選択確率l 経路選択確率
l 経路𝜎 = 𝑘$ :$89 = {𝑘9 , 𝑘5 , … , 𝑘: } について考える

𝑘9 は出発リンク， 𝑘; = 𝑑とする

l 経路𝜎 の選択確率は:<55:<5 𝑒 6 7 𝑘$=5 𝑘$ => *-./$89$89

𝑃 𝜎 = h𝑃 𝑘$=5 𝑘$ = h∑#∈? *57 𝑎 𝑘$ => #𝑒6:<557 𝑘$=5 𝑘$ => *-./ <> *= h 𝑒6$89
5:<5 5<6> *07 𝑘$=5 𝑘$
6=𝑒h 𝑒$89(RL 9)25

経路選択確率l 経路選択確率
l 𝑣 𝜎 = ∑:<5$89 𝑣(𝑘$=5 |𝑘$ ) と置くと𝑃 𝜎

5:<5 5<6> *07 𝑘𝑘=𝑒h 𝑒 6 $=5 $$89=

57 @𝑒65> *𝑒6 0=57 @𝑒6∑@A∈B

57 @A𝑒6(RL 10)l ここでΩ は全選択肢の集合
l 式(10) の分子は経路𝜎の確定項，分母は可能な経路の確定項の和であるため，このモデルは無限の選択肢を持つMNL モデルと同等
l 従ってIIA 特性を持つ

26

Discounted Recursive Logit Model (Oyama, Hato, 2017)l Bellman方程式に時間割引率𝛽を導入
𝑉%0 𝑘 =max

7&'( ∈$ 3=𝐸max

9𝐸 -8&"

7&'( ∈$ 3𝛽 8+" 𝑢 𝑠"/'𝑠" ; 𝜃𝑣% 𝑠"/' 𝑠" + 𝛽 𝑉%0 𝑠"/' + 𝜇𝜀% (𝑠"/')l 時間割引率を導入し，将来価値の低減，未来の不確実性や近視眼的な選択を表現
𝑃+, 𝑠()* 𝑠(!!" "#𝑒 #=1(𝑣+ 𝑠()* 𝑠( + 𝛽 𝑉+, 𝑠()* )𝜇=1∑"#$%-∈.("#) exp (𝑣+ 𝑠′()* 𝑠( 𝛽 𝑉+, 𝑠()* ′ )𝜇exp<"#$%∈𝒜("# )1,𝛿 𝑠()* 𝑠( exp1(𝑣+ 𝑠()* 𝑠( + 𝛽 𝑉+, 𝑠()* ) , 𝑠( ∈ 𝒜
𝜇𝑠( = 𝑑27

Discounted Recursive Logit Model (Oyama, Hato, 2017)l 時間割引率を含んだ価値関数の計算方法
行列 𝑴 の要素
𝑀(1 (1./ = Z0ベクトル 𝒛 の要素
𝛿 𝑠C=5 𝑠C

5(7 𝑠C=5 𝑠C )6𝑒

∀𝑠C ∈ 𝒜𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒𝑧7&' %
5 7&= 𝑒4

RLモデルと異なり，価値関数の解は非線形方程式系に従うことになる𝑧"# =< 𝑀"#"#$% (𝑧"#$% )/ ,∀𝑠( ∈ 𝒜"#$%∈𝒜

1,式(dRL 1)を行列表記で書くと以下の通り
𝐳 = 𝐌𝐗(𝐳) + 𝐛(dRL 1)

𝑠( = 𝑑(dRL 2)

28

Discounted Recursive Logit Model (Oyama, Hato, 2017)l 時間割引率を含んだ価値関数の計算方法
行列 𝐌 の要素
𝑀(1 (1./ = Z0

ベクトル 𝐳 の要素
𝛿 𝑠C=5 𝑠C

5(7 𝑠C=5 𝑠C )6𝑒

∀𝑠C ∈ 𝒜𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒𝑧7&' %
5 7&= 𝑒4

𝐳 = 𝐌𝐗(𝐳) + 𝐛(dRL 2)

𝐗(𝐳) は 𝒜 × 𝒜 の行列で要素𝑋(𝑧)(1 = (𝑧(1 ) D式(dRL 2)を解くには価値関数が不動点に達するまで反復計算する必要がある

29

Discounted Recursive Logit Model (Oyama, Hato, 2017)l 時間割引率を含んだ価値関数の計算方法→不動点までの反復計算

𝐳 = 𝐌𝐗(𝐳) + 𝐛(dRL 2)

𝐗(𝐳) は 𝒜 × 𝒜 の行列で要素𝑋(𝑧)(1 = (𝑧(1 ) D
反復計算の手順（Mai, Fosgerau, and Frejinger, 2015 ）1. ベクトル 𝐳 (𝟎) として初期化する
2.

𝐳 (𝟏) = 𝐌𝐗 𝐳 𝟎

3.

𝐳 (𝒏=𝟏) − 𝐳 (𝒏) < 𝛾 となるまで2を繰り返す+ 𝐛 として更新する

30

dRL モデルの推定：コードの概要
データのインプット

価値関数（V, z）の初期値の設定
対数尤度関数の定義• 選択確率の計算

𝑃% =𝑴"⊙𝒛#𝑴"𝒛• 対数尤度の計算
価値関数の定義• 時間割引率を導入しているため価値関数z を閾値まで求解
最尤推定
収束するまで繰り返し価値関数の計算

31

dRL モデルの推定：データのインプット
l インプットデータとして最低限必要なもの
l 経路データ... トリップID／発リンク／着リンク／終点リンクの列からなる経路の集合

l ネットワークデータ... リンクの接続関係を表す

l リンクデータ... リンクid ごとの説明変数を格納モデルの説明変数となる

32

dRL モデルの推定：対数尤度関数の定義
配布コード：RL_Estimation_Assignment_Node.R を参照ください！

←終点リンクdごとに計算している!
p# b← 𝑧b = 𝑒 "を計算

選択確率の計算← 𝑃b =

𝑴$ ⊙𝒛%
𝑴$ 𝒛を計算33

dRL モデルの推定：価値関数の計算（newV）
関数名：newV の中身←終点リンクd ごとに計算している← 𝐌と𝐛を定義← 即時効用 𝐌 を更新

34

dRL モデルの推定：価値関数の計算（newV）
𝐳 = 𝐌𝐗(𝐳) + 𝐛(dRL 2)𝐗(𝐳) は 𝒜 × 𝒜 の行列で要素𝑋(𝑧)7& = (𝑧7& ) P

反復計算の手順（Mai, Fosgerau, and Frejinger,
2015 ）

1. ベクトル 𝐳 (𝟎) として初期化する
2.

𝐳 (𝟏) = 𝐌𝐗 𝐳 𝟎3.

𝐳 (𝒏=𝟏) − 𝐳 (𝒏) < 𝛾 となるまで2を繰り返す+ 𝐛 として更新する

35

dRL モデルの推定：最尤推定と価値関数の更新←最尤推定の実施←パラメータの更新←価値関数の求解（newVの実行）←価値関数の更新

36

RLモデル, DRL モデルを用いた研究の例
l 経路選択モデル（RLモデル）
l Fosgerau, Mogens, Emma Frejinger, and Anders Karlstrom. "A link based network route choicemodel with unrestricted choice set." Transportation Research Part B: Methodological 56 (2013): 7080.

l 災害時におけるドライバーの経路選択モデル（DRLモデル）
l Oyama, Yuki, and Eiji Hato. "A discounted recursive logit model for dynamic gridlock network
analysis." Transportation Research Part C: Emerging Technologies 85 (2017): 509-527.

l Activity-travel スケジューリングの時空間経路選択モデル（mixed RLモデル）
l Zimmermann, Maëlle, et al. "Capturing correlation with a mixed recursive logit model for activitytravel scheduling." Transportation Research Part C: Emerging Technologies 93 (2018): 273-291.

l 歩行者環境におけるアクセシビリティ分析のための歩行者経路選択モデル（RLモデル）
l Zheng LIANG*, Ka Fai NG, Hong K. LO “A Link-Based Approach For Measuring Walking Accessibility”,2022
37