Behavior Study for Transportation, Graduate school, University of Yamanashi

モデルのパラメータ(母数)とは• 離散型選択モデルの母数

ランダム効用モデル– 確率分布を特徴付ける数– 例えば，期待値，分散等期待値と分散が母数– 変動要因を規定する少数の要素に分解する
V1 = βx

U ~ N (V1 , σ )• パラメータ推定– 観察されたデータに基づいて母数を統計的に推定すること

パラメータ推定の理論と実践

V1 = βx

山梨大学
佐々木邦明
2

効用関数の構造母数• 選択条件と結果が異なる– 一人一人の効用を推定すればよいが大変

推定アルゴリズム• 期待値に構造を仮定
仮定する
仮定– 期待値の変動要因を観測可能変数で規定• 手段選択の例：費用，所要時間，頻度(待ち時間)，アクセ
ス・イグレス時間・・・– 期待値の関数形は• 線形モデル，非線形モデル・・・
34

確率モデルによる個人の選択記述

最尤推定法
Ln (θ | x ) = ∏ f (xi |θ )U = U1 − U 2n• 点推定量を求める最もポピュラーな方法• 右上の式をθの関数とみなしたものが尤度関数

i =1

例：データとして（1，2）が得られたとき，全体の平均が何であるかを考える• 確率的選択モデルは効用がある分布に従うと仮定• ある選択肢の選択確率は確率分布を積分する

２項モデル

平均がいくつだったら（1，2）が得られやすいか？

• 尤度関数を最大化するθの値を最尤推定量とするのが最尤推定法

選択モデルの場合，f が選択確率

MaxLikelihood = max ∏ Pi ( xi | θ )個人の選択確率を全員で掛け合わせる
θ ∈Θ

i5

• ロジットモデルやNL等のGEVモデルはこの積分した結果が積分のない形式で表現可能

P(1)MaxLikelihood = max ∏ ∫ f ( xi | θ )０∏ ∫ f (xi | θ ) をどうやって最大化するか？
θ ∈Θ

ii

最尤推定のアルゴリズムの考え方••• Newton-Raphson法– テイラー展開の２次近似を利用して進める– 解の収束が早い(ステップ数が少ない)– ヘッセ行列を計算するので時間がかかるθˆML

対数尤度関数の段階的な最大化

1. 初期値を与える
2. 初期値周りで勾配(1次微分）等を用いて次の推定値の方向を決める
3. 初期値付近のステップサイズを１次微分，２次微分を用いて適切に決めて次の推定値を決める
4. 収束基準(尤度関数の一時微分ベクトル)を判定し，収束していない場合は，現在の値を初期値として2に進む

代表的な繰り返し計算法

Ln (θ | x )最大値をいきなり求めるのは無理な場合

θ2

θ ST6

θ

θ 2 = θ1 − H −1 g1
H:尤度関数の二階微分
ヘッセ行列
g: 尤度関数の一階微分• 準Newton法（BFGS法）– ヘッセ行列を単位行列から，尤度関数と尤度関数の一階微分の差を用いて逐次近似する．最終的にヘッセ行列の逆行列に収束．
– 二階微分をする必要がないので計算が早いθ1

ロジットモデルのみ全域で単峰性あり

78

そもそも推定不能

パラメータ推定がうまくいかない• 収束するとはθ1とθ2が同じになる– g1が0になる• 収束しない– θ2が計算不能– 無限に繰り返す

H-1が存在しない
変数が完全相関
変数が効用関数に影響していない関数の近似状況
2次関数近似初期値の問題• 最大値において唯一解が求まらない可能性がある（最大値となるパラメータベクトルが無数にある）→Identification Problem
例
U1

U3

U2

そもそも推定不可能

効用に適当な数を足しても選択確率に変化はない

推定プログラムに誤り

定数項は相対的数値• 局所最適解– 見かけ上の最大化

0910

駐車場の台数や料金無料は効果が計測不能？

ちょっと考えてみる• 例：3カ所の買物目的地選択• 例えば自動車保有ダミーddを使ってみる

A:駐車場有料，B：駐車場有料，C:駐車場無料
このとき，駐車場無料は大きな選択要因
ならば駐車場無料ダミー変数を選択肢に導入

U1i = (α1 + d di × β1 ×100 )× 1 + K + ε 1i

U 2i = (α 2 + d di × β1 × 200)×1 + K + ε 2i

U1i = α1 ×1 + β1 × 0 + K + ε 1i

U 3i = (d di × β1 × 1000)× 1 + K + ε 3i

U 2 i = α 2 × 1 + β1 × 0 + K + ε 2 i
U 3 i = α 3 × 0 + β1 × 1 + K + ε 3 i
これは区別できない1112

複雑なモデル（Logit Kernel)のIdentification Problem（Walker 2004）

積分（INTEGRATION））
積分（• オーダー条件– 選択肢固有の誤差パラメータは S ≤ J (J − 1) − 12

• ランク条件– 選択肢固有の誤差パラメータはS = Rank ( Jacobian(vecu(Ω ∆ ))) − 1– ΩΔ：制約無しとした場合の誤差共分散– vecu： ΩΔの独自項をベクトル化• 等号条件Ω ∆ = Ω normalized
∆
13

数値積分

シミュレーション(乱数)による積分• 個人ごとの確率を計算するためには定積分することが必要だ
が，閉じた形で表現できないものは数値積分を行う必要がある．

• 面積を求める別の方法を考える– Gaussの求積法：近似誤差を最も小さくなるように区間と幅を定める方法．
– Gauss-Legendre法：多項式近似ができる関数においては，多項式の次元にあわせた分割を行うと非常に精度よく積分を行える
方法

98• 2次元乱数を[３×９]の範囲に散らばらせて，それが積分したい関数より上か(青丸)下か(赤丸) を判定する．

9• 数値積分：関数を多項式近似などを用いて積分可能な形にして，ある間隔の関数の和に直して積分を行う方法．

1487654123• ドットをうまく大量に発生させれば精度よく積分を近似できる．

複雑な関数を簡単な１次関数で近似した例．

157654123

この図だと点の合計は20．そのうち赤丸１１個，青丸９個なのでこれより関数の積分結果は14.85となる

16

乱数の生成

制約つきの乱数• 何を積分したいのか？

1• 例えば平均０で分散1の正規分布の-3から1までの間の乱数．

0.8– 乱数の発生ルーチンを条件に合わせて改良

0.60.4• 対応する分布関数の間の乱数を発生させてそれを密度に対応させる• 制約つきの乱数• 多変量正規分布の乱数0.2• この変換によって，密度の高いところはたくさん，密度の薄いところは少しという密度関数に応じた乱数になる．

• 乱数と同等の性質を持ち，同じ値が出ず，必要な全域がカバーされる数値の生成法があるとより効率的な積分が可能

0.6

0.4
0.2-4-22417

多変量正規分布の発生

18

準乱数

• 例えば標準正規分布μを平均b標準偏差sに変換するためにはε = b + sµ と線形変換をすればいい• 同じ事を多変量正規分布に応用する．標準正規分
布の乱数ベクトルηをε = b + Lηと変換する．ただしLはLL’=Ωとなる行列．すると はMVN(b,Ω)になる．

• 高次積分において乱数による積分の精度を上げるためには，計算速度が多大に犠牲になる

32

イメージ1

• 対処方法として系統的な数列（準乱数）がある

 ε 1   b1   l11 0  η1   =   +  
 ε 2   b2   l12 l22 η 2 
ε 1 = b1 + l11η1

例えば，2次元の例を示すと

123

ε 2 = b2 + l12η1 + l22η 21920

準乱数の例

シミュレーションによる尤度計算

 代表的例としてHalton数列がある．その計算方法は素数pに対して
ttt{st +1 = st , st + 1 p , st + 2 p L, st + ( p − 1) p}

例えばp=3ならば，初期値0として１/3, 2/3, 1/9, 4/9, 7/9, 2/9, 5/9, 8/9・・・• 多次元化– 数列の異なる素数pを決めて，それぞれに応じて数列を作り多次元化する．

• 正規分布化– 数列を制約つきの乱数発生と同様の変換をして正規分布化•手順
1. 誤差項の密度関数から選択肢の数の次元の(準)乱数を発生させる
2. この乱数を誤差の値として，各代替案の効用値を計算（積分）する
3. 代替案iの効用値とその他の代替案の効用との値を比較し，それらの大小関係を1-0の変数Gで記述する．
4. 1～3のステップを繰り返す．その反復回数をRとする．
1 R
Pi = ∑ G r
5. シミュレーションされた確率は R r =1となり，この値は不偏推定量である．

効用を確定値にする

確定的に選択を決定

比率を確率に置き換える

これを尤度として最大になるようにパラメータをアップデートする

21

連続型シミュレーション

22

シミュレーションベースのパラメータ推定法

1. 誤差項の密度関数から選択肢の数の次元の乱数を発生させる
2. この乱数を誤差の値として，各代替案の効用値を計算する
3. この効用をロジット変換を行い，連続的な0-1変数に変換する．この変換後の値をSとする
4. 1～3のステップを繰り返す．その反復回数をRとする．

• シミュレーション尤度最大化（MSL）– シミュレーションによって計算された確率を尤度として，最大化を行う．

• 特性

5. シミュレーションされた確率は Pi = R ∑ S r となり，この値r =1は不偏推定量である．
1 R

23– サンプル数と乱数発生回数に依存する．
– 乱数発生回数が十分大きいと一致性や漸近的有効性を持ち解析積分と一緒の特性を持つ．
– 乱数発生回数がサンプル数に対して小さく固定されると一致性もない．

24

統計的な視点• 推定尤度に基づくモデル構造の選択– 投入する情報量(データ)が一緒なら尤度関数の大小で判定– 情報量の異なるモデル間ではAIC等の統計量による比較まとめ• 推定パラメータに基づく変数の選択– ｔ値による取捨選択– 推定値の符号による判定• 費用と所要時間のパラメータは負2526モデル構造はどうやって選ぶ？

モデル構造の再考について• 行動に関するアプリオリな知識による（仮説）– 手段選択には費用は影響する• 線形効用関数に費用を導入する．でもタクシーは選択肢に入らない• 費用は非補償型のモデル構造とする

検証– バスとパークアンドライドの選択において，誤差項に相関があるだろう• NLやCNLの採用を考える– 誤差項の相関が無視できるような構造を考える• 推定結果からの後天的な知識による（仮説の修正）– モデルの推定結果が仮説を支持しないときに仮説が間違っている可能性を考える• 線形効用関数で費用のパラメータが負にならない→キャプティブな選択があるのでは？
• ネスティッドロジットモデルのログサムの係数が0～1にならない→ネスト構造が違っている？

• 誤差項の構造の仮定＝選択肢の属性における非観測要因の仮定• 誤差構造を複雑にするより，観測変数をきちんと取る方がいい

検証

27• 質の高いデータベースはどこに存在するか？
• いやいや，そもそもデータはどうなんだ．．．．
28

ベイズの定理• ベイズの公式

ベイズ推定/更新

P(B | A)P( A)P (B )
P(B | Ai )P ( Ai )P( Ai | B ) =∑ P(B | Ai )P( Ai )P( A | B ) =i

例

バスに乗っていたときのGPSデータに某アルゴリズムを適用すると98%の確率でバスと判定する
バスに乗っていないときのデータに同じアルゴリズムを適用すると5%の確率でバスと誤って判定する
母集団ではバスに乗っている割合は3%である．
あるデータがバスに乗っていると判定されたときに，本
当にバスに乗っている確率を求めよ

ベイズ統計• ベイズの公式– H:仮説– D:データ– Θ：母数– D:データ

P(D | H )P(H )P (H | D ) =P (D )P (H i | D ) =事前確率• ベイズ統計の基本式

P(θ | D ) =P(D | H i )P(H i )∑ P(D | H i )P(H i )P(D | θ )P(θ )P (D )π (θ | D ) = kf (D | θ )π (θ )Kは確率を１にする定数• ゼミ出席確率θの分布を考える

i

事後確率

※P(D)は定数母数を求めるベイズ推定尤度

事前分布がなければ最尤推定– ゼミ開始前 事前確率は一様分布
π (θ | D1 ) = k × θ ×1 → 2θ
– 第１回に出席する
π (θ | D2 ) = k × θ × 2θ → 3θ 2
– 第２回に出席する– 第３回に欠席する π (θ | D3 ) = k × (1 − θ )× 3θ 2 → 12 × (1 − θ )θ 2
– 第４回に出席する

出席状況と出席確率密度

自然な共役分布3• 事後分布があっさり求められるのは，尤度と事前分布の数式がまとめられるから

2.5
2

第1回出
1.5

第2回出
第3回欠

1第4回出

• 尤度を介したベイズ更新をしても事後分布が変らないと毎回同じ形式で計算可能
これを自然な共役分布という

0.5000.20.4

0.60.81

事前・事後分布

自然な共役分布
事前分布

尤度

事後分布

ベータ分布

二項分布

ベータ分布

正規分布

正規分布

正規分布

逆ガンマ分布

正規分布

逆ガンマ分布

ガンマ分布

ポアソン分布

ガンマ分布

尤度

事前・事後分布尤度

事前分布

事前分布と尤度が正規分布の例• 事前分布が正規分布，尤度が正規分布の例( x − µ )2

1 − 12eπ (µ | x1 ) = k ×2π

1 −×e2π( µ − µ 0 )22∝ex +µ  µ− 1 0 2 

−12×2

※データを追加すると母数の分散が必ず小さくなる

2MCMC法によるベイズ推定

メトロポリス法p( x') ≥ p( xt ) ⇒ xt +1 = x'

• 自然な共役分布をあきらめるとき– たとえば，先ほどの例で事前分布を正規分布とする( D −θ )21 − 12π (θ | D1 ) = k × θ ×e2πp( x') < p( xt ) ⇒→ ??

P( xt +1 = x') = r

P( xt +1 = xt ) = 1 − r

発生した乱数による次の点の発生確率に応じて採用と不採用が変化する• 数値積分の応用– 事前分布と尤度からなる事後分布に，確率密度に応じたサンプリングを行い母数を計算– サンプリングした母数に対応する頻度から母数の期待値を求める

6.0
8.1
7.2
3.9
4.5
2.5
6.6
4.1
2.05.5続く• 例– 平日一日の平均トリップ数が右の表のように与えられる– 分散の事前分布を逆ガンマ，平均値の事前分布を正規分布とする– 詳細つり合い条件 p (wt )s(wt → wt +1 ) = p (wt +1 )s(wt +1 → wt )を満たすマルコフ連鎖

メトロポリス法のMCMCで平均値と分散の分布を計算してみる

1日の平均トリップ数の例

μの分布の実現値
120• 事後分布∝（逆ガンマ分布×正規分布）カーネルに対して，乱数を発生させメトロポリスアルゴリズ
ムで繰り返し計算を行うと下の図になる

100

80
7660

54μ

3

σ

2
140

20

145
89
133
177
221
265
309
353
397
441
485
529
573
617
661
705
749
793
837
881
925
969
1013
1057
1101
1145
1189
1233
1277
1321
1365
1409
1453
1497
1541
1585
1629
1673
1717
1761
1805
1849
1893
1937
1981

0• この時μとσの平均値は単純に平均で計算可能で5.01と2.17（※1000まで廃棄）

04

4.5

55.5

66.57

まとめ
• ベイズ推定は基本的には事後情報＝事前情報＋データ情報という構造であり，更新型の推定に適している• 事前分布を適切に選ぶ（一様分布）と通常の統計学の最尤推定と同じ解が得られる• 推定結果には明確な確率分布がない場合も可能（複雑な分布形状を許容）• 参考・引用文献– 入門ベイズ統計学 松原望– 道具としてのベイズ統計 涌井良幸
– ベイズモデリングによるマーケティング分析 照井伸彦