第6回基礎ゼミ資料
Practice NL&MXL from R
平成30年5月18日（金）朝倉研究室 修士1年 小池卓武

使用データ～横浜プローブパーソンデータ～
主なデータの中身
・トリップID・目的・出発，到着時刻・総所要時間・移動距離
・交通機関別の時間，距離・アクセス，イグレス時間，距離・費用
・代表交通手段
・代替手段生成可否
・性別，年齢等の個人属性
etc..

1

NLモデルによる推定-NLの効用関数
時間

公共費用

2
誤差項
(選択肢間相関あり?)

𝑈𝑡𝑟𝑎𝑖𝑛 = β𝑡𝑖𝑚𝑒 𝑡𝑖𝑚𝑒𝑡𝑟𝑎𝑖𝑛 + β𝑐𝑜𝑠𝑡 𝑐𝑜𝑠𝑡𝑡𝑟𝑎𝑖𝑛 + ε𝑃𝑢𝑏𝑙𝑖𝑐_𝑡𝑟𝑎𝑖𝑛
𝑈𝑏𝑢𝑠 = β𝑡𝑖𝑚𝑒 𝑡𝑖𝑚𝑒𝑏𝑢𝑠 + β𝑐𝑜𝑠𝑡 𝑐𝑜𝑠𝑡𝑏𝑢𝑠 + ε𝑃𝑢𝑏𝑙𝑖𝑐_𝑏𝑢𝑠
𝑈𝑐𝑎𝑟 = β𝑡𝑖𝑚𝑒 𝑡𝑖𝑚𝑒𝑐𝑎𝑟 + β𝑐𝑜𝑠𝑡 𝑐𝑜𝑠𝑡𝑐𝑎𝑟 + ε𝑃𝑟𝑖𝑣𝑎𝑡𝑒_𝑐𝑎𝑟

私有𝑈𝑏𝑖𝑘𝑒 = β𝑡𝑖𝑚𝑒 𝑡𝑖𝑚𝑒𝑏𝑖𝑘𝑒 + β𝑐𝑜𝑠𝑡 𝑐𝑜𝑠𝑡𝑏𝑖𝑘𝑒 + ε𝑃𝑟𝑖𝑣𝑎𝑡𝑒_𝑏𝑖𝑘𝑒

𝑈𝑤𝑎𝑙𝑘 = β𝑡𝑖𝑚𝑒 𝑡𝑖𝑚𝑒𝑤𝑎𝑙𝑘 + β𝑐𝑜𝑠𝑡 𝑐𝑜𝑠𝑡𝑤𝑎𝑙𝑘 + ε𝑃𝑟𝑖𝑣𝑎𝑡𝑒_𝑤𝑎𝑙𝑘
μ

公共鉄道

1バス

1私有

自動車 自転車 徒歩

下位ネストのスケールパラメータは1に固定し，上位ネストのスケールパラメータμを推定．その際，0 < μ < 1が条件．

NLモデルによる推定-初期設定

3～初期設定～nl_maxLik_1.rを展開！
Session→ Set Working Directory → To Source File Locationをしてから実行
読み込むデータの1行目が・データでなく列名ならheader=T・データから始まるならheader=Fnrow(行列)：引数行列の行数を返す
numeric(引数)：引数個の0を返す→パラメータ推定値の初期値(交通手段数-1)個の定数項+
パラメータβ(時間，料金)2つ+
スケールパラメータ(μ)×1つ↓
推定値7つ…初期値として0を与えておく

NLモデルによる推定-対数尤度関数の設定

4～対数尤度関数を定義～
function(引数){算術式}を用いて対数尤度を計算するための関数をあらかじめ定義．
初期値として対数尤度LL=0を設定しておく．

NLモデルによる推定-対数尤度関数の設定～exp(V)を定義～

Data$代替手段生成可否trainには，trainが代替交通手段として生成されうるなら1，そうでなければ0が入る．
SOL変数を同程度のオーダーにして推定結果を安定させるためにそれぞれ100で除す．(100分，100円単位にする)

5

NLモデルによる推定-対数尤度関数の設定～上位ネストでの公共・私有選択確率～

ログサム変数：𝑉𝑑′ = 𝑙𝑛𝑚′ exp(𝑉𝑚′ )上位ネストでの(公共・私有)選択確率：𝑃(𝑑) =exp(μ 𝑉𝑑 +𝑉𝑑′ )′ ))exp(μ(𝑉+𝑉′′𝑑𝑑𝑑′

d：上位ネスト，μ：上位ネストのスケールパラメータ6

NLモデルによる推定-対数尤度関数の設定～下位ネストでの交通機関選択確率～

exp(𝑉𝑚 )下位ネストでの(交通機関𝑚)選択確率：𝑃(𝑚|𝑑) =𝑚′ exp(𝑉𝑚′ )exp(𝑉𝑚 )交通機関選択確率：𝑃 𝑚 𝑑 𝑃 𝑑 =𝑚′ exp(𝑉𝑚′ )exp(μ 𝑉𝑑 + 𝑉𝑑′ )′𝑑 ′ exp(μ(𝑉𝑑 ′ + 𝑉𝑑 ′ ))7

NLモデルによる推定-対数尤度関数の設定

8～対数尤度の計算式を定義～
log 𝑃 の計算時に真数条件：log 𝑃 : 𝑃 > 0に反してエラーになるのを防止
代表交通手段が𝑖ならば𝐶. 𝑖には1が入る

対数尤度LL：(𝐿𝑜𝑔𝐿𝑖𝑘𝑒𝑙𝑖ℎ𝑜𝑜𝑑)：

𝑁𝑛𝐽
𝑖 𝑑𝑖 ln 𝑃𝑛 (𝑖)N：全サンプル数，n：個人，J：全交通機関選択肢数，
i：交通機関選択肢，di：交通機関選択ダミー(1or0)，Pn(i)：交通機関選択確率↓
未知パラメータに対して最大化することでパラメータ推定：maxLik関数

NLモデルによる推定-MaxLik関数～対数尤度の最大化～
maxLik関数を用いて最尤推定(Optim関数の代替案！)
maxLik(対数尤度関数，start=パラメータ初期値，method=“～”)
ヘッセ行列(hess=)，勾配(grad=)，制約条件(constraint=)も引数として入れられるが，今回は単純なケースとして上記の式形で計算．
Method...最大化手法にも様々ある...
・NM(Nelder-Mead)←計算が粗い..
・NR(Newton-Raphson)←勾配(grad=)が必要
・BFGS(Broyden-Fletcher-Goldfarb-Shanno)←gradientが無くても大丈夫
Etc..
今回はmethodについて考える機会を設けるためにあえて→今回はBFGSを採用 関数内にmethodを設定したが，何も設定しなければ最適なmethodが自動適用される仕様になっている．

9

NLモデルによる推定-結果の出力

10～maxLik関数でパラメータ推定～
Broyden-Fletcher-Golden-Shannoメソッドを最適化メソッドとして採用．
sqrt関数：引数の平方根を計算
diag関数：引数(行列)の対角成分を計算
solve関数：引数(正方行列)の逆行列を計算

パラメータ推定値/標準偏差 = t値

Printで出力→

スケールパラメータμが0 < μ < 1を満たしていない!!
→ネスト構造の再検討が必要..

パラメータ推定値は左から定数項×4，β×2，μMXLモデルによる推定-MXLの効用関数

11～ランダム係数モデルの効用関数～
𝑈𝑖 = β𝑡𝑖𝑚𝑒 𝑡𝑖𝑚𝑒𝑖 + β𝑐𝑜𝑠𝑡 𝑐𝑜𝑠𝑡𝑖 + ε𝑖= β𝑡𝑖𝑚𝑒 + σ𝑡𝑖𝑚𝑒 η𝑡𝑖𝑚𝑒 𝑡𝑖𝑚𝑒𝑖 + β𝑐𝑜𝑠𝑡 + σ𝑐𝑜𝑠𝑡 η𝑐𝑜𝑠𝑡 𝑐𝑜𝑠𝑡𝑖 + ε𝑖
β𝑥 ：平均，σ𝑥 ：平均，η𝑥 ：標準正規分布に従う乱数， ε𝑖 ：IIDガンベル分布

個人の嗜好性によってパラメータβが確率的に変化→ランダム係数と呼ぶ
NLと同じ交通手段選択肢を用いてパラメータ推定！

MXLモデルによる推定-初期設定

12～データの読み込み～

nrow(行列)：引数行列の行数を返す→サンプル数numeric(引数)：引数個の0を返す→パラメータ推定値の格納配列(交通手段数-1)個の定数項 + パラメータβ(時間，料金)の平均，標準偏差↓
推定値8つ…初期値として0を与えておく．
繰り返し計算によりパラメータ推定を行うので，その繰り返し回数R=100を設定．

MXLモデルによる推定-乱数の生成13～ハルトン数列法で準乱数を生成～
ハルトン数列法
素数xの累乗を分母においた分数を0～1区間で複数生成し，qnorm関数を用いて標準正規分布に従うような順乱数にする手法．

𝑠𝑡+1 = 𝑠𝑡 , 𝑠𝑡 +

1
2
𝑥−1,
𝑠+,…
,
𝑠+𝑡
𝑥𝑡 𝑡 𝑥𝑡𝑥𝑡↓
Qnorm関数で標準正規分布に従わせる←ハルトン数列

MXLモデルによる推定-対数尤度関数の設定

14～対数尤度関数の設定～
パラメータβに個人嗜好性を持たせたランダム係数を表現するために，平均，標準偏差を推定．

対数尤度をR回計算し，その和をRで除すことで対数尤度の平均を出す→シミュレーション尤度．
その和の受け皿としてSimLLを定義．

MXLモデルによる推定-MXLの効用関数～ランダム係数を定義～
平均

標準偏差

ハルトン数列法で生成した乱数=確率的に変化するパラメータ：ランダム係数

ランダム係数：β𝑖 = β𝑖 + ση𝑖
β𝑖 ：平均，σ：標準偏差， η𝑖 ：乱数

15

MXLモデルによる推定-MXLの効用関数

16～選択確率，対数尤度の定義～

対数尤度の計算，パラメータ推定をR回(100回)行うため，待ち時間が長い．
暇なので結果を出力．
確認せずともよい人はprintの前に#を付けてコメント化してもOK．

MXLモデルによる推定-結果の出力

17～結果の出力～
NL同様，対数尤度の最大化にmaxLik関数を使用．

Printで各種推定結果を出力！

パラメータ推定値は左から定数項×4，時間のランダム係数の平均および標準偏差，料金のランダム係数の平均および標準偏差

参考文献

18

参考文献：
1)北村隆一・森川高行・佐々木邦明・藤井聡・山本俊行，
交通行動の分析とモデリング

2)東京大学羽藤研究室，http.//bin.t.u-tkyo.ac.jp/kaken/3)Rdocumentation maxLik
https://www.rdocumentation.org/packages/maxLik/versions/1.3-4/topics/maxLik