スライド 1

クラシックな機械学習の入門
10. マルコフ連鎖モンテカルロ法
Imputation Posterior
次元の呪い
ＭＣＭＣ
Metropolis-Hastings アルゴリズム
Ｇｉｂｂｓサンプリング
by 中川裕志（東京大学）
Sampling法の必要性
 ＥＭやＶＢで解析的に処理して数値計算に持ち込
める場合ばかりではない。
 ＶＢの場合、因子化仮定による近似も問題。
 シミュレーションで解く方法は対象に係わらず使え
る。
 ただし、一般にシミュレーションで良い解に収束させ
るには時間が膨大にかかる。
 次元の呪いも影響してくるので、事態は良くない。
 そこで、効率の良いシミュレーション法が必要にな
る。
Samplingによるシミュレーションを用いる局面
 EM,VBでは以下の更新式を解析解として求めておき、
繰り返し計算を行うようになる。
 繰り返し
 P(Z,X| θold,M)argmaxθexp(E[logP(Z,X| θold,M)])でθを推定
 この途中で行う期待値計算E[logP(Z,X | θold,M)])の中で
周辺化などの積分計算の解析解を求めることが困難
である場合は、その部分を期待値計算で用いる確率
分布p(Z)にかなうようにsamplingしたデータでシミュ
レーションによって計算する。
EMのQ計算をsamplingで置き換える方法
Q関数
Q( , old )   p( Z |  old , X ) log p(Z , X |  )d
 この積分を以下のようにしてsamplingで求める
事後分布 p( Z |  old , X ) を用いて観測データか
らsamplingして｛ Z(l) ｝(l=1,2,..,L)を得る。
こうして得た｛ Z(l) ｝を用いて、以下のようにθ
L
を求める.
1
Q( , old )   log p( Z ( l ) , X |  )
L
l 1
M-stepは、これらの結果を使って通常のEMと
同じように行う。
観測データ
観測データ
の値Ｙ
これらを
Σ して
分布
Ｐ（Ｘ）に
沿う積分
を計算
Ｘ
Ｐ（Ｘ）
この部分の面積に比例するように
観測データをサンプリングする
Ｘ
Imputation Posterior: IP Algorithm
I-step
 p( Z | X )   p( Z |  , X ) p( | X )d の計算にあたって、sample
θ(l) l=1,2,..,L をp(θ|X）にしたがって、取り出し、これを用い
てp(Z|θ(l),X)によってZ(l)を取り出してp(Z|X)を計算
P-step
 p( | X )   p( | Z , X ) p(Z | X )dZ の計算にあたって、 I-stepで得たsample ｛ Z(l) ｝を用いて、以下のよう
にθを求める.
1 L
p( | X )   p( | Z ( l ) , X )
L l 1
しかし、sampleはうまくやらないと効率が悪く
なってしまう. （次元の呪い）
次元の呪い(Curse of dimensinality)
(x1,x2)の赤、
青、黄の教師
データが図のよ
うに配置。
新規のデータ●
が何色かを決
める方法を教師
データから学習
x1
次元が高いと、周辺部分ばかりに
サンプル点が集まり、中心部分に
サンプルデータが非常に少ない。
x2
 データを記述するための情報の種類(feature)を表
す次元Dが増加するとどうなるか？
 各次元のfeatureの間の関係を示す係数の数がfeature
のとりうる値N(f)の積で増加。
 各featureの値:rによって以下のように識別するとしよう。
 If 1-e<r<1 then A else B
 1-e<r<1という条件を満たすデータ空間の
体積：K（１－(1-e)D ）のデータ空間全体:K 1Dに対する割
合は 1-(1-e) D
 この割合はDが大きくなると急激に１に近づく。
 つまり、Dが大きい次元の空間では、縁の部分の体積ば
かり大きくなる。よって、実測された教師データは縁の部
分にばかり集まる。
 よって、rの値の閾値が比較的小さいところになる場合は、
教師データが極端に少なくなり、学習の精度が落ちる。
Markov Chain Monte Carlo:
MCMC の基本アイデア
次のデータ点が指定した領域から外に出ていたら、
元のデータ点に戻って、次のデータ点を選びなおす。
近似分布の導入
 sampling でp(Z)の正確な評価が困難でも、 p(Z)と比例
する分布 ~p ( Z ) なら乱数によってsamplingすることがで
きるとしよう. 両者の関係は正規化定数Zpで下のように
1 ~
書ける。
p( Z ) 
p( Z )
本当に知りたいのは
これだが、
Zp
~
p(Z )
 ただし、 Zpは未知。だから、
のままで計算できる
ことが望ましい。
 さらに近似分布(proposal分布）q(z(t+1)|z(t))にしたがい、
z(t)からsamplingしてz(t+1)を求める。
この条件付き確率（遷移確率）
は、問題から求めておく
Markov連鎖
 Markov連鎖の定常性
when p(z (t 1) )  Z( t ) p(z(t 1) | z(t ) ) p(z(t ) )　, then p(z(t 1) )  p(z (t ) )
 ただし、定常的であっても、一つに状態に収束しているとは言い
切れない。循環的になっているかもしれない。
 下の図のような状況はまずい。つまり、a,b,cばかりにシミュレー
ションでサンプルする点が集中してしまい、p(z)の示す領域全体
の体積（＝確率）が正確に求まらない。
循環的になると、
サンプル点はこの
a,b,cばかりにな
る！
p(z)
a
Ｚ
１
c
１
ｚ
Ｚ１
b
bｚ
Markov連鎖
 Markov連鎖の定常性
when p(z (t 1) )  Z( t ) p(z(t 1) | z(t ) ) p(z(t ) )　, then p(z(t 1) )  p(z (t ) )
 ただし、定常的であっても、一つに状態に収束しているとは言い
切れない。循環的になっているかもしれない。
Ｚ１
b
１ｚ
a
Ｚ１
c
１ｚ
循環的でない定常状態になるためには次の詳細釣り合い条件
が必要。この条件から定常性もいえる。
詳細釣り合い条件：
ある p*が存在し , すべてのz, z 'に対して　p* ( z ' ) p ( z '  z )  p* ( z ) p ( z  z ' )
詳細釣り合い条件を書き換えると
~
~
p( z ) / Z p
p( z '  z )
p( z )
p( z )

~
~
p( z  z ' ) p ( z ' ) p ( z ' ) / Z p p ( z ' )
つまり、zz’の遷移確率を決めるには、正規
化定数：Zpは知らなくてもよい。
Metropolis アルゴリズム
Step 1: 近似分布(proposal分布）として、q(z(t+1)|z(t))を設定しておく。
Step 2: q(z(t+1)|z(t)) を用いてsample z(1), z(2),…. z(t)を求めておき、
~
pˆ  が計算できる状態にした上で、以下のstep,3,4,5で生成する。
Step 3: 新規のsample znew を生成する。
~ new

Step 4: A( z new , z ( t ) )  min1, pˆ ( z )  とする。
 


~ (t ) 
pˆ ( z ) 
また u  0,1
をあらかじめ決めた閾値とする。
Step 5: if A(znew,z(t)) > u then z(t+1)=znew
else z(t)からznewを生成してstep 3へ戻る
z(t)を捨てずにとっておき、 q(z(t+1)|z(t))を用いて再度のsample znew
の生成で役立てるところがポイント
Metropolisアルゴリズムの味噌
Step 4: A( z new , z ( t ) )  min1, ~p ( z new )  とする。
 ~ (t ) 

p( z ) 
 
また。 u  0,1
をあらかじめ決めた閾値とする。
Step 5: if A(znew,z(t)) > u then z(t+1)=znew
else z(t)からznewを生成してstep 3へ戻る
new
(t )
~
~



 u であること、すなわち新たな znew が
p
z

p
z
 つまり、
前の値 z(t)よりある程度以上、狙った分布に近い場合のみ採
用されることになる。
~
p ( z new )
 ~ (t ) だから、正規化定数Zpを知らなくてもよいところが味噌
p( z )
 対称性 q(za|zb)=q(zb|za) の条件は t∞のとき、samplingされ
た z の分布がp(z)に近づくことを保証する。
Metropolis-Hastings アルゴリズム
Step 1: 近似分布(proposal分布）として、q(z(t+1)|z(t))を設定しておく。
Step 2: q(z(t+1)|z(t)) を用いてsample z(1), z(2),…. z(t)を求めて、
~
pˆ  が計算できる状態にしたうえで、以下のstep3,4,5で生成する。
Step 3: 新規のsample znew を生成する。
~ new
(t )
new
Step 4:
とする。


ˆ
p
(
z
)
q
z
|
z
new
(t )
k

Ak ( z , z )  min1, ~ ( t )
new
(t ) 
 pˆ ( z )qk z | z 




ただしkは検討すべき状態遷移の集合上のインデックス。Step5の条件を満た
すものを選ぶ。（Znewに行ったきり戻ってこられないようなもの、つまりmin
の第２項の分子が小さいものは使わない）
 
また u  0,1 をあらかじめ決めた閾値とする。
Step 5: if A(znew,z(t)) > u then z(t+1)=znew
else z(t)からznewを生成してstep 3へ戻る
Gibbs Sampling
z=(z1z2,…zM)からなるベクトルの場合, zの全
要素を一度に生成して、条件A(znew,z(t))>uを
満たすかをチェックするのは効率が悪い。
そこで、zjnewをz1(t),.., zj-1(t)が与えられた条件で
samplingする方法が有効。
Gibbs Samplingでは、zjのsampling と評価
を、条件A(znew,z(t))>uのチェックの代わりに、
それ以前の結果を利用して行う。
Gibbs Sampling
1. z1z2,…zMの初期化
2. 以下をt=1からTまで繰り返す
1. Sample z1(t+1) z1new :p(z1new|z2(t),z3(t),…,zM(t))
2. Sample z2(t+1) z2new :p(z2new|z1(t+1),z3(t),…,zM(t))
3. Sample z3(t+1) z3new :p(z3new|..,z2(t+1),z4(t),…,zM(t))
…
j. Sample zj(t+1) zjnew :p(zjnew|..,zj-1(t+1),zj+1(t),…,zM(t))
…
M. Sample zM(t+1) zMnew :p(zMnew|..,zM-1(t+1))
水平の線はｘ２の値を固定
して次の点に進むイメージ
垂直の線はｘ１の値を固定
して次の点に進むイメージ
Ｘ２
Ｘ１
zkを更新するときには z \ kは固定されていること

に注意
z ( t ) \ k  z new \ k
qk z new | z ( t )   p ( z newk | z ( t ) \ k )
qk z ( t ) | z new   p ( z ( t ) k | z new \ k )
p( z new )  p( z newk | z ( t ) \ k ) p z ( t ) \ k 
p( z ( t ) )  p( z ( t ) k | z new \ k ) p z new \ k 
これを利用して Metroplis Hasting法を書き換える
~ new

pˆ ( z ) qk z ( t ) | z new  
new
(t )

Ak ( z , z )  min1, ~ ( t )
new
(t ) 
 pˆ ( z ) qk z | z  
 p( z newk | z ( t ) \ k ) p z ( t ) \ k qk z ( t ) | z new  

 min1,
(t )
new
new
new
(t ) 
 p( z k | z \ k ) p z \ k qk z | z  
 p( z newk | z ( t ) \ k ) p z ( t ) \ k  p( z ( t ) k | z new \ k ) 
  1
 min1,
(t )
new
new
new
(t )
 p( z k | z \ k ) p z \ k  p ( z k | z \ k ) 
よって、新規のデータは常に採用される  Gibbs Sampling
条件付正規分布
 Gibbs Sampling で行った方法（他の変数の条件
は同一にした条件付き分布で、1変数だけを予測す
る一例を正規分布で示す。
 変数ベクトルｚをｘとｙに分割すると
P(y|x=a)
y
p(y)
Ｘ＝a
 D次元変数ベクトル z をスカラー x とD-1次元ベクト
ル y に分割する。
N ( z | μ,  )
多次元正規分布
 x
z   
 y
 μx 
μ   
 μy 
where   
精度行列：   
1
 x 2


 yx
T
とすると
 xy 

 yy 
 xy   yx
T
 xx


 yx
 xy 

 yy 
 xy   yx
T
 ここで多次元正規分布の指数の肩の項は次式
1
 ( z  μ)T  1 ( z  μ) 
2
1
1
 ( x  μ x )xx ( x  μ x )  ( x  μ x ) xy ( y  μ y )
2
2
1
1
T
 ( y  μ y )  yx ( x  μ x ) 
( y  μ y )T  yy ( y  μ y )
2
2
-(G-10)
一般に正規分布 N ( z | μ, ) の指数の肩は次式で書け、右
辺の第１項、第２項の係数の部分に着目すれば期待値、共分
散が求まる。
1
1 T 1
T 1
 ( z  μ)  ( z  μ)   z  z  zT  1μ  const
2
2
-（G-20)
条件付正規分布p(x|y)の期待値μx|yと共分散Σ x|yをこの方法を
（G-10)式に適用して求めよう。ー問題
yを定数とみなしてxの分布を求めれば、条件付分布になるか
ら（G-10)の第１項のxの２次の項の係数が共分散。すなわち
1
 xxx x
2
により
 x| y 2  xx 1
一方、（G-10)においてxの１次の項がΣ -1 μ
xxx  x   xy ( y  μ y )
これは次式
これにより
μ x| y   x| y xx  x   xy ( y  μ y )
2
1
 x| y  xx より
2
1
 μ x  xx  xy ( y  μ y )
次に、これらの結果を共分散行列を用いて書き直す
1
 xy 
  2 x| y  xy 


   xx
 において ( Matrix 1 )を使えば


 




yx
yy
yx
yy




xx  ( x2   xy  yy1  yx ) 1  xy  ( x2   xy  yy1  yx ) 1  xy  yy1

μ x| y  μ x   xy  yy1 ( y  μ y )
 x| y 2   x2   xy  yy1  yx
粒子フィルタ
• 内部状態𝑥𝑡が 𝑥𝑡+1 = 𝑔 𝑥𝑡 + 𝑣𝑡 (ガウス雑
音）という更新式で変化していく
• 観測されるのは𝑦𝑡 = ℎ 𝑥𝑡 + 𝑢𝑡 𝑣𝑡 (ガウス雑
音）
• 𝑡毎の観測値を用いて𝑥𝑡を推定、𝐸𝑡 𝑓 𝑥𝑡 を
計算するメカニズム
– カルマンフィルタによく似ている。
2-３
)
)
)
f(
)

Download Report