ベイズ統計データ分析＋兵藤先生のメモ
第４章 離散選択モデル

2009/6/24(水)ベイズゼミ#2M2 北川直樹
http://www2.kaiyodai.ac.jp/~hyodo/Logit_by_R.pdf

目次
1. MCMC法のアルゴリズム
ギブズサンプラー
メトロポリスへスティング

2. 離散選択モデル
二項プロビットモデル
ギブズサンプラー，MHアルゴリズム
多項プロビットモデル

3. 数値例
最尤推定
ベイズ推定

MCMC法のアルゴリズム
Markov チェイン
3蓮の葉(A,B,C)の上を蛙が飛び移る．
各葉の遷移確率が3*3行列の場合，t→∞の蛙の存在確率を計算する。
等高線は対数尤度関数．
より良い解を求めて蓮の葉を飛び回る．
蓮の葉間の飛び方を規定するのが，MCMC 法

QuickTimeý Ç²êLí£ÉvÉçÉOÉâÉÄ
Ç™Ç±ÇÃÉsÉNÉ`ÉÉÇ¾å©ÇÈÇžÇ½Ç…ÇÕïKóvÇ-ÇQuickTimeý Ç²êLí£ÉvÉçÉOÉâÉÄ
Ç™Ç±ÇÃÉsÉNÉ`ÉÉÇ¾å©ÇÈÇžÇ½Ç…ÇÕïKóvÇ-ÇÅBMCMC法のアルゴリズム
Gibbs Sampler(GS)
各パラメータが与えられたとき，他のパラメータの条件付事後分布が定義できる場合に使われる．
例）平均0,分散1 ,相関係数ρの2 変量正規分布

θ1 θ 2 ≈ N(ρθ2 ,1− ρ 2 ), θ 2 θ1 ≈ N( ρθ1,1− ρ 2 )2 つの正規乱数を入れ子で発生させれば，得るべき2 変量正規分布に従う．
蛙の例で言えば，蛙は次に跳ぶべき蓮への望ましい跳躍方向を知っており，それに従いrandomness をもって飛び回る．

MCMC法のアルゴリズム
Metropolis-Hastings アルゴリズム(MH)
条件付事後分布が書き下せない場合に使われる．
様々なモデルに対応できるが，計算効率が悪い．
アルゴリズムの手順
次に飛ぶ蓮をランダムに選ぶ．
GS と違い，その蓮に飛んでしまったことを想定し，解の”改善度”を [0,1] の範囲内に収まる基準値で計算．
[0,1] の一様乱数を発生させ，両者の大小関係で，この蓮に飛ぶか否かを決断する．
元の蓮に止まる場合，今の計算が無駄になる．

離散選択モデル
二項プロビット
個人i，選択肢集合yi={1,0}，説明変数xi
パラメータβ，誤差項εi

1, if x iβ + εi > 0yi = otherwise
0,
条件付き事後分布は，P(y i = 1 x i , β ) = Φ(x iβ )P(y i = 0 x i , β ) = 1− Φ(x iβ )Φの関数形
正規分布：プロビット
ガンベル分布：ロジット

離散選択モデル
二項プロビット
Φの数値積分が困難
潜在変数y*を用いてギブズサンプラーを適用した事後分
布推定が可能 IN(0,∞) (x iβ,1), if y i = 1∗
yi β ~  IN(−∞,0] (x iβ,1), if y i = 0
IN[a,b](μ,σ2)は，区間[a,b]で切断させた正規分布
βの事前分布をβ〜N(b0,B0)とすると，y*が与えられたときの条件付き事後分布は，

β y ∗ ~ N(b1,B1)
ただし， b1 = B1 (B0−1b0 + X T y ∗ ), B1−1 = B0−1 + X T X, X = (x1,L x n )T

二項プロビットモデル
ギブズ・サンプラー
1. 繰り返し回数を s = 0,L,ndrawとする．
β s= 0 ≈ N(b2. s = 0とし，初期値である事前分布を決定
0 ,B0 )
ss+1
3. β から β を生成
sβを生成する
1. y ∗s+1i
 N(0,∞) (x iβ,1), if y i = 1∗yi β ≈ 
 N(−∞,0] (x iβ ,1), if y i = 0
1. β s+1 y ∗s+1i を生成する

β y ∗ ≈ N(b1,B1 )−1∗

ただし， b1 = B1 (B0 b0 + X y ),T

1. s + 1 < ndrawのときに3.に戻る．

B1−1 = B0−1 + X T X, X = (x1,L x n )Tsのとき終了= ndraw

二項プロビットモデル
MHアルゴリズムβ s= 0 ~ N(b0 ,B0 ) を決定した後，提案分布の密度関数ssρ
(β, β y)q(β,βy)を与えて，受容確率を計算し，サンプリングする．
例えば，最尤法によるパラメータ推定値 と共分βˆˆˆβ とすると，次y ~ N(β s−1, τ 2 ∑)∑∑散 を使い，提案密度分布をの受容確率を計算する．
 p(β y) q(β, β s−1 y) ρ(β , β y) = min
,1
s−1
s−1
 p(β y) q(β , β y) s−1

二項プロビットモデル
MHアルゴリズム
1. 繰り返し回数を s = 1,L,ndrawとする．
2. s = 1とし，初期値である事前分布
3. β s−1から β ≈ q(β s−1, β y)を生成
4. 受容確率 ρ(β s−1, β y)を生成するを決定
β s−1 ≈ q(b0 ,B0 ) p(β y) q(β, β s−1 y) ρ(β , β y) = min,1
s−1
s−1βy)β,βy)p(q(
s−1

1. 一様乱数 u ≈ u(0,1)
を生成する
 β if u ≤ ρ(β s−1, β y)
β =  s−1otherwise
βs

2. s + 1 < ndrawのときに3.に戻る．

sのとき終了= ndraw

多項プロビットモデル
個人iがp個の選択肢集合から選択効用uijが最大となる選択肢jを選択する問題．
p個の選択肢集合に対して，j=0,‥p-1とラベル付けし，p-1個の潜在効用wを考える．
w i = X id β + εi
w ij = uij − uip , εij = ξ ij − ξ ip , εi ~ N(0,∑)
 x Ti1 − x Tip 
X id = M xT − xT ip  ip−1

wij<0のとき，個人iは選択肢pを選択する．
wijの中で1つでも正があれば，wijが最大の選択肢jを選択

多項プロビットモデル
多項プロビットを推定する場合，(β,Σ)の事前分布は，
β ~ N(b0 ,B0 ), ∑ ~ IW (v 0 ,R0 )は，逆ウィシャート分布である．
事後分布は，β ∑,w ~ N(b1,B1 ), ∑ β,w ~ IW (v1,R1 )ただし，
b1 = B1 (B b + ∑ (X ) ∑ w i ), B = B + ∑ (X id )T ∑−1 X id
n−10 0−1d Ti

i=1−11−10ni=1

v1 = v 0 + v, R = R + ∑ (w i − X id β )(w i − X id β )T−1
1−10ni=1vは逆ウィシャート分布の自由度

多項プロビットモデル
ギブズサンプラーでは，以下の手順で条件付き分布を算出

β ∑,w and ∑ β ,w → w β,∑, y, X d → y β ,∑,w, X dwiは(p-1)次元に切断された正規分布であるため，wの条件付き事後分布を直接求めるのは難しい．
wijの事後分布を以下のような切断された正規分布とする．
w ij w i− j , y i , β ,∑ ~N(m ij , τ 2jj ) × [I( j = y i )I(w ij > max(w i− j ,0)) + I( j ≠ y i )I(w ij < max(w i− j ,0))]d2
m ij = x djj ′β + (−σ jj γ j,− j )T (w i,−i − X i,−β),τjjj = 1/σ jj

γj,-iはΣ-1のj番目の要素を除いたj番目の行
σjjはΣ-1の(j,j)要素

 γ1′ −1∑ = M 
γ ′  p−1 離散選択モデル
多項プロビットのギブズ・サンプラー
繰り返し回数を s = 0,L,ndrawとする．
,R0
適当な初期値 b0 ,B0,v 0を与える．
s = 0とし，事前分布
β s, ∑ s,を生成
wssを生成w sから w s+1 β s,∑,∑ s
β sから β s+1 w s+1を生成
s+1を生成∑sから ∑ s+1 β s+1,wのときに3.に戻る．
s + 1 < ndrawのとき終了
s = ndrawパッケージと関数
二項プロビット
ギブズサンプラー
MCMCprobit(M)rbprobitGibbs(B)

多項プロビット
ギブズサンプラー
rmnpGibbs(B)※(M):MCMCpackパッケージ(B):bayesmパッケージ

二項ロジットMHアルゴリズム
MCMClogit(M)

多項ロジットMHアルゴリズム
MCMCmnl(M)rmnlIndepMetrop(B)数値例
歩行者の速度-角度選択モデル
選択肢集合
速度3分類
角度5分類
合計15分類
説明変数
目的地を目指す
距離，角度

速度を保つ
加速，減速

定数項

10°

0° -10°
-25°

25°

1.2vn
1.0vn

0.8vn

Vvd n = β ddist ddist vd n + β ddir ddird n +β acc Iv,accv n + β dec Iv,decv n +β vdn

数値例
最尤推定

ベイズ推定

尤度関数
Nj

L = ∏ ∏ Pn (i x,θ )n=1 i=1

誤差構造
多項ロジットモデル
最適化手法
準ニュートン法

ベイズの定理p(y θ ) p(θ )P(θ y) =p(y)
誤差構造
多項ロジットモデル
MCMC法
MHアルゴリズム
棄却回数:1000回
繰返回数:10000回
初期分布
平均0，分散0推定結果最尤推定

ベイズ

最尤推定

ベイズ

説明変数

推定値

T値

平均

平均/分散

説明変数

推定値

T値

平均

平均/分散

βddist-5.821-7.05**-6.968-8.11

β60.0540.17-2.622-11.85βddir-0.058-0.20-0.292

-1.03

β7

1.544

5.40**-1.051

-6.17

βacc

0.181

1.71

0.290

2.64

β8

2.713

9.66**

0.134

0.84

βdec

-1.264-10.65**-1.998-15.27β9

1.623

5.70**-0.968

-5.70

β1

0.570

1.86

-2.139-10.36

β10

0.245

0.79

-2.368-11.15β2

1.492

5.12**-1.139-6.39

β111.3524.73**-0.692

-4.08

β3

2.300

8.04**-0.315

-1.90

β12

2.263

8.44**

0.318

2.44

β4

1.292

4.41**-1.362

-7.58

β13

3.317

12.14**

1.277

11.07

β5

0.686

2.25*

-1.967-11.57

β14

2.456

9.22

0.525

4.07

サンプル

3085

サンプル

3085

初期尤度

-8354

受容率ρ01.000

適合度

0.132

最終尤度

-7251

受容率ρ

0.641

■ 赤は，推定値⇔平均値，t値⇔平均/分散，ベイズの絶対値が大きい．
■ 青は，推定値-2≒平均値，t値-10≒平均/分散
ベイズの値が小さい．

25°

Vvd n =

10°

0° -10°1413129875

4 3 2-25°
11

β ddist ddist vd + β ddir ddird +nβ acc Iv,accv n + β dec Iv,dec v n + β vd

n1.2vn610n11.0vn
0.8vn

分布図
赤の変数

青の変数

推定結果
最尤推定法との比較
MCMC 法から得られたパラメータの分散共分散と，ヘッセ行列の逆行列から得られた分散共分散の一致性は証明できるのであろうか．

解の安定性
概ね，正規分布に漸近していることが理解できるが，必ずしも単峰性が満たされているわけでもない．これは何度も乱数初期値や，繰り返し回数，MCMC法を変えたりしながら常に成立する特性か否かを判断する必要があろう．

MCMC法の比較
independent MH
θがθ(i-1)と独立
q(θ (i−1),θ y) = q(θ y)random walk Metropolis
θがθ(i-1)に依存
q(θ (i−1),θ y) = f (θ − θ (i−1) )