Z,X

モデル推定
潜在変数がある場合のモデル推論
EMアルゴリズム
EMの混合正規分布への適用例
変分ベイズ法
EP法
潜在変数を考慮する推論
 prior:w(=潜在変数+分布のパラメター)と観測データ
{x,t}(ただし、xは入力、tは結果)の分布にハイパーパ
ラメターα、βを導入する。
 まず、観測データから事後確率を最大化するpriorをベ
イズで求める。このためにいろいろな技法あり。(例え
ば経験ベイズ法)
p(w | x, t,,  )  p(t | x, w,  ) p(w |  )
 新規データxに対する予測 t を計算するのは以下の式

p(t | x, x, t, ,  )  p(t | x, w,  ) p( w | x, t, ,  )dw
モデル推定のための学習法
 事前知識のない場合はK-meansなど
 EM(Expectation and Maximization)アルゴリズム
 モデル学習法の古典
 変分ベイズ法
 予測モデルqを既知のモデル(prior)とのKL divergence
を最小化する関数として変分法で繰り返し計算で求める。
速度は速い。
 KL(p||q)=Σp log(p/q)
 MCMC(Markov Chain Monte Carlo)
 モデルのパラメター推定をシミュレーションで解いてしまう
方法。速度は遅い。
経験ベイズ法 :
事前分布のパラメターの初期値の推定
 ベイズモデルのパラメターを恣意的に決めると気持ち悪い
 ベイズモデルのパラメターを観測データに基づいて求める
方法
 π( θ | α ):事前分布 ただし、 α は事前分布πを決めるパ
  arg max
ラメター
 観測データx から尤度p (x|θ)が与えられたとき、これを用
いて事前分布のパラメターαmaxを求める。
max

 max  arg max  p x |    |  d

 ( EB10)
経験ベイズ法 :
事前分布のパラメターの初期値の推定例
 max  arg max  p x |    |  d

 ( EB10)
 多項分布、ディリクレ分布の例
観測データ iの出現回数
Dir (  | X ,  )  Mult  X |  Dir (  |  )   ii mi 1
i 1
 max  arg max Dir (  | X ,  )

 arg max Dir (  |   X )

K
i 1
i 1
m:i X  (m1 ,..., mK ),  mi  M , i  0
K

K
K
 0  M 
 arg max
ii mi 1

1  m1   K  mK  i 1

EMアルゴリズム
観測変数 observed variable(顕在変数
manifest variable)
隠れ変数 hidden variable(潜在変数 latent
variable)
 例:混合正規分布のとき、どの正規分布から生成された
データかを示す変数
観測変数のデータを用いて母集団の統計モ
デルの潜在変数を推定する(最尤推定値)
 パラメータの値を点で推定
このための数値解法:EMアルゴリズム
Expectation and Maximization Algorithm
 観測されたデータ(=観測変数の実測データ):X
 潜在変数のとる値:Z
 統計モデルのおける未知のパラメター:θ
logの中のΣが現
 対数尤度関数(log likelihood func.)
れると嫌だ!
L( X |  )  log p( X |  )  log  p( X , Z |  )
Z
 未知パラメターθの最尤推定値は
ˆ  arg max L( X |  )

EMアルゴリズムの導出 その1
P( X , Z |  )  P( Z | X , ) P( X |  )
を用いると
log P( X , Z |  )  log P( Z | X ,  )  log P( X |  )
ここで Zに関する確率分布 q( Z )を用いると上の式は以 下のようになる
P( X , Z |  )
P( Z | X , )
q
(
Z
)
log

q
(
Z
)
log
 log P( X |  )
Z

q( Z )
q(Z )
Z
この式は以下のように みなせる
P( X , Z |  )
q
(
Z
)
log
  KL(q( Z ) || P( Z | X ,  ))  log P( X |  )
Z
q( Z )
すなわち、
log P( X |  )  KL(q( Z ) || P( Z | X ,  ))   q ( Z ) log
Z
P( X , Z |  )
q( Z )
EMアルゴリズムの導出 その2
log P( X |  )  KL( q( Z ) || P( Z | X ,  ))   q( Z ) log
Z
P( X , Z |  )
q( Z )
ここで Zに対する分布の推定は 、 の現在の推定値  oldを用いると
P( Z | X ,  old )なので、これを q( Z )とする。すなわち
log P( X |  )  KL( P( Z | X ,  old ) || P( Z | X ,  ))   P( Z | X ,  old ) log
Z
さて、  pldを更新した  newを、log P( X |  )を改善(より大きくす
選びたい。そこで、 log P( X | 
log P( X | 
new
new
P( X , Z |  )
P( Z | X ,  old )
る)ように
)  log P( X |  ) を評価してみよう。
old
)  log P( X |  )
old
=0
 KL( P( Z | X ,  old ) || P( Z | X ,  new ))  KL( P( Z | X ,  old ) || P( Z | X ,  old ))
new
old
P
(
X
,
Z
|

)
P
(
X
,
Z
|

)
old
  P( Z | X ,  old ) log

P
(
Z
|
X
,

)
log

old
old
P
(
Z
|
X
,

)
P
(
Z
|
X
,

)
Z
Z
EMアルゴリズムの導出 その3
log P( X | 
new
)  log P( X |  )
old
 KL( P( Z | X ,  old ) || P( Z | X ,  new ))
  P( Z | X ,  old ) log P( X , Z |  new )   P( Z | X , old ) log P( X , Z |  old )
Z
Z
KL( P( Z | X , old ) || P( Z | X ,  new ))は定義より非負。
また、 P( Z | X , old ) log P( X , Z |  old )は newには関係ないので、
Z
結局、できるだけ真の に近い  newは

 P( Z | X ,
old
) log P( X , Z |  )を最大化するような  である
Z

 new  arg max  P( Z | X , old ) log P( X , Z |  )

Z
なお、以後 Q( |  old )   P( Z | X ,  old ) log P( X , Z |  ) と書く。
Z
 E Z   oldに固定 log P( X , Z |  )
EMアルゴリズムの詳細
Q( |  old )   P( Z | X , old ) log P( X , Z |  )  E Z, old log P( X , Z |  )
Z
P( Z | X ,
old
P( Z , X , old )
)
 P(Z , X , old )
Z
 初期化:θoldを適当に決める
 以下のEstep, Mstepを収束するまで繰り返す
 E step: P(Z|X, θold)を計算
 M step: θnew =argmax Q(θ| θold) とし、 θold をθnewに更新
θ
しかし、上記の導出から分かるように、log P( X |  ) をθを動か
しながら最大化しているので、局所解に陥る可能性あり。
EMとQ関数の再考
 QをEZ( θold)[logP(Z,X| θ)]で書き直すと
 以下のEstep, Mstepを収束するまで繰り返す
 E step: P(Z|X, θold)を計算
 M step: θnew =argmax θ EZ( θold)[logP(Z,X| θ)]とし、 θold を
θnewに更新
 つまり、 P(Z,X| θ)を θold を固定してZで期待値をとることで
θに関する情報を教師データZから集めてにθ反映させること
を繰り返しての良い推定値を求めている。
K-meansに似ている
θはベクトル。よって、 Mstepでは、θの全要素
を一度に更新する式を求めている点に注意。
楽屋裏の話:なぜQ(θ| θ(t))か?
 なぜ
L( X |  )  log p( X |  )  log  p( X , Z |  )
を直接最適化せずにQ(θ| θold)か?
Z
 Q(θ| θold)=ΣP(Z|X ,θold)logP(Z,X |θold)
である。すなわち、 LはlogΣの形で解析的に扱いにくい。QはΣlogの形で
解析的に扱いやすい
 では、そもそもなぜ尤度ではなく対数尤度なのか?
 多くの分布は exp(f(x)) だったり(ex 正規分布)、べき乗の形であるから、
logをとると扱いやすい。
 なぜ、expやべき乗なのか?
 複数の確率変数の共起は、各々の確率の積だから、という説明も可能
 理論的な背景から見れば「指数分布族:exponential family」であることが効
果を発揮している。
EMの適用例:1変数正規分布
このような観測データから、混合
正規分布の各パラメタを推定する
混合正規分布
 1 N ( x | 1 ,  1 )   2 N ( x | 2 ,  2 )
要素の正規分布 N ( x | 1 ,  1 )
N ( x | 2 ,  2 )
where
1   2  1
問題の定式化
 1 , 1 ,  1 ,  2 , 2 ,  2
p( x )   1 N ( x | 1 ,  1 )   2 N ( x | 2 ,  2 )
1   2  1
潜在変数zk  0,1
p( z1  1)   1
p( z2  1)   2
推定するパラメター:
p( z)   1 1   2
z
z2
(GM 10)
p( x | z )  N ( x | 1 ,  1 ) z1 N ( x | 2 ,  2 ) z2
 p( x ) 
(GM 11)
 p( z ) p( x | z )   N ( x |  ,  )  
k 1, 2
k
k
ここで次のように定義
  zk   p ( zk  1 | x ) 
1
1
1
2
N ( x | 2 ,  2 )
される  zk を導入し、ベイズの定
p( zk  1) p( x | zk  1)
 p( zk  1) p( x | zk  1)
k 1, 2

 k N ( x | k ,  k )
 j N ( x |  j ,  j )
j 1, 2
(GM 30)
(GM 20)
理で書きなおす
いよいよEMアルゴリズムの適用
 次のパラメターに適当な初期値を設定:  ,  ,  2 k
k
k
 E step: P(Z|X, θ(t))を計算
 ただし、観測されたデータはN個あるとする。
 実際には、P(Z|X, θ(t))ではなく Zの期待値 を求めておくことにする
。
N
z
z
2
(GM 10)と (GM 11)より p( Z | X ,  ,  ,  )    1 N ( xn | 1 ,  12 )  2 N ( xn | 2 ,  22 )
n1
n 1
ここで Z(すなわち zn1 , zn 2 まとめて znk where k  1,2)を評価する。
Eznk  

 znk  k N ( xn | k , k )
2
znk
  N ( x
j
|  j , j )
2
n

znk
 k N ( xn |  k ,  k 2 )

  znk 
2
  j N ( xn |  j ,  j )

znj
(GM 40)
j 1, 2
znj
これはデータ xnの znk (  どちらの正規分布が選 ばれるか)への寄与を 表す
N
note!
定義より
2
N
2
  z    Ez  N
nk
n 1 k 1
N
n 1 k 1
N
  z    Ez  N
n 1
nk
nk
n 1
nk
k
 k , k ,  k の現在の値
の
から計算した更新値
n2
 Mstep
(GM 10)( GM 11)より N
p( X , Z |  ,  ,  )  p( X | Z ) p( Z )    1 N ( xn | 1 ,  1 ) n1  2 N ( xn | 2 ,  2 ) n 2
z
z
n 1
N
2
 log p( X , Z |  ,  ,  )   znk log  k  log N ( xn | k ,  k )
n 1 k 1
(GM 40)より 

N

2

E Z log p( X , Z |  ,  ,  )   Eznk  log  k  log N ( xn | k ,  k )
2
2
n 1 k 1
N
2


   znk  log  k  log N ( xn | k ,  k )
n 1 k 1
次に  znk を固定した上で (GM 50)において
 k , k ,  k 2を最大化してこの stepでの最適値を求める。
2
(GM 50)
kの最適化し k に更新する
(GM 50)を kで微分してゼロとおく
new

k
N
2
  z log 
nk
n 1 k 1
k
。
 log N ( xn | k ,  k )


( xn   k ) 2
 znk log  k 
 const 

2
2 k
n 1 k 1


N
( xn   k )
   znk 
0
2


k
N
2
k
n 1
N
k
new

 x  z 
n 1
N
n
nk
  znk 
n 1
1

Nk
N
 x  z 
n 1
n
nk
(GM 70)
2 new
 k の最適化し、  k に更新する
(GM 50)を  kで微分してゼロとおく
2

N
 k
2

。

2

2



z
log


log
N
(
x
|

,

 nk
k
n
k
k )
n 1 k 1

 k
 1

( xn   k ) 2
2
 znk  log  k 
 const 

2
2 k
n 1 k 1
 2

N
2
2

1
( xn   k ) 2 


   znk 

0
2
2 2
n 1


2 k
 2 k

N
 
N
N
n 1
n 1
(GM 801)
(GM 801)より  k    znk     znk  ( xn   k ) 2
2
N
 k
2 new

2



z
(
x


)
 nk n k
n 1
Nk
(GM 80)
 kの最適化


 2

L    znk  log  k  log N ( xn | k ,  k )      k  1
n 1 k 1
 k 1

N
N
L
 znk 
new

     znk    k  N k   k  0
 k n 1  k
n 1
N
2
2
N

一方    znk    k  
k 1  n 1

2
以上の2式より
k 
N  0
Nk
N
(GM 90)
ここでは、  ( z k ) が古い  k , k ,  k を反映して計算された値であった。
それを固定して、loglikelihoodを最大化する新たな  k , k ,  k 2 を求めているわ
けだ。
2
以上の(GM70)(GM80)(GM90)によって  k , k ,  k
log likelihood (GM50)が収束しなければEstepに戻る
2
の更新式が求められた。
EMの適用例:混合多変数正規分布
1変数の場合と似ているが、少し難しいところあり
このような観測データから、混合
正規分布の各パラメタを推定する
この例では1変数の正規分布だ
が以下の導出は多変数の正規分
布を仮定している。
 x, 1 , 2はベクトル ,
1 ( 11 ),  2 ( 21 ),は行列
混合正規分布
 1 N ( x | 1 , 1 )   2 N ( x | 2 , 2 )
要素の正規分布 N ( x | 1 , 1 )
N ( x | 2 ,  2 )
where
1   2  1
Σは共分散行列、Λは精度行列で
あることに注意
下の式はk変数の正規分布においてN個のデータがある場合
x



i
1
1
N
1


     xi1  1 ,..., xik  k 
N i 1
 xik  k 
  x21   x1, x 2 


1
 

 
 x1, x 2   xk2 


問題の定式化
 1 , 1 , 1 ,  2 ,  2 ,  2
p ( x)   1 N ( x | 1 , 1 )   2 N ( x |  2 ,  2 )
1   2  1
潜在変数z k  0,1
p ( z1  1)   1
p ( z 2  1)   2
推定するパラメター:
p(z )   1 1   2
z
z2
(GM 10)
p ( x | z )  N ( x | 1 , 1 ) z1 N ( x |  2 ,  2 ) z2
 p ( x) 
 p( z
k 1, 2
k
(GM 11)
) p ( x | z k )   1 N ( x | 1 , 1 )   2 N ( x |  2 ,  2 )
ここで次のように定義
 z k   p( z k  1 | x) 
される   z k を導入し、ベイズの定
p ( z k  1) p( x | z k  1)
 p( zk  1) p( x | zk  1)
k 1, 2

 k N ( x | k ,  k )
 jN (x |  j ,  j )
j 1, 2
(GM 30)
(GM 20)
理で書きなおす
いよいよEMアルゴリズムの適用
 次のパラメターに適当な初期値を設定:  k , k ,  k
 E step: P(Z|X, θ(t))を計算
 ただし、観測されたデータはN個あるとする。
 実際には、P(Z|X, θ(t))ではなく Zの期待値 を求めておくことにす
る。
N
z
z
(GM 10)と (GM 11)より p ( Z | X ,  , ,  )    1 N ( xn | 1 , 1 )  2 N ( xn |  2 ,  2 )
n1
n 1
ここで Z(すなわち z n1 , z n 2 まとめて z nk where k  1,2)を評価する。
 znk  k N ( xn | k ,  k ) nk
z
Ez nk  
z nk
  N ( x
j
n |  j, j )

z nj
 k N ( xn |  k ,  k )
   z nk 

N
(
x
|

,

)
 j n j j

(GM 40)
j 1, 2
z nj
これはデータ xnの z nk (どちらの正規分布が選 ばれるか)への寄与を 表す
N
note!
定義より
2
N
2
  z    Ez  N
nk
n 1 k 1
N
n 1 k 1
nk
n 1
nk
のk , k ,  kの現在の値か
ら計算した更新値
N
  z    Ez  N
n 1
nk
k
n2
 Mstep
(GM 10)(GM 11)より N
p ( X , Z |  , ,  )  p ( X | Z ) p ( Z )    1 N ( xn | 1 , 1 ) n1  2 N ( xn | 2 ,  2 ) n 2
z
z
n 1
N
2
 log p ( X , Z |  , ,  )   znk log  k  log N ( xn | k ,  k )
n 1 k 1
(GM 40)より N
2
E Z log p ( X , Z |  , ,  )   Eznk log  k  log N ( xn | k ,  k )
n 1 k 1
N
2
    znk log  k  log N ( xn | k ,  k )
n 1 k 1
N
2


    znk  log  k  log N ( xn | k , k1 )
n 1 k 1
次に   znk を固定した上で (GM 50)において
 k , k ,  kを最適化してこの stepでの最適値を求める。
(GM 50)
kの最適化し k newに更新する
(GM 50)を kで微分してゼロとおく

k
N
2
  z log 
nk
n 1 k 1


k
k
。
 log N ( xn | k ,  k )
1


1
T



z
log


(
x


)

(
x


)

const


nk 
k
n
k
k
n
k
2


n 1 k 1
N
2
N
   znk  k
n 1
1
 xn   k   0
N
k
new

 x  z 
n 1
N
n
nk
  z 
n 1
nk
1

Nk
N
 x  z 
n 1
n
nk
(GM 70)
 kの最適化し、  k に更新する
new
ここが多変数だと難しくなる部分
new
1
  k (  k1 )の最適化し、 new
k (  k
(GM 50)を  k で微分してゼロとおく

 k
  z log 
2
N
nk
n 1 k 1


 k
N
2
1
2
k
)に更新する
。

 log N ( xn | k , 1k )


  znk  log  k  ( xn  k )T  k ( xn  k )  const 
n 1 k 1
1
2
 1  log  k 1  ( xn  k )T  k ( xn  k ) 

    znk 
0
 k
2
n 1

 2  k
N
(GM 801)
 1  log  k 1  ( xn  k )T  k ( xn  k ) 
 znk 

0

2


2


n 1
k
k


(GM 801)のおのおのの項の微分 を計算する
N

 log  k
1
 k
 k

(GM 801)
T
 k
1
(GM 802)
trace( AB)
 BT より
A
 log  k


1

( xn   k ) T  k ( xn   k )   k 
trace k ( xn  k ) ( xn   k ) T 
 k
 k
 k
公式 x T Ax  trace( Axx T )   k  ( xn  k ) ( xn   k ) T    k  ( xn   k ) ( xn   k ) T  0
1
1
T
N
(GM 801)(GM 802)(GM 803)より   znk  k
1
n 1
N
  znk  k  k
1
n 1
1
new
 k
1new
   znk ( xn  k ) ( xn   k ) T
n 1
N
1



z


 nk k N k
n 1
N
 k
N

  z ( x
n 1
nk
n
 k ) ( xn   k ) T
Nk
(GM 803)
(GM 80)
 kの最適化
 2

L    z nk log  k  log N ( xn |  k ,  k )     k  1
n 1 k 1
 k 1

N
N
 z nk 
L
new

      z nk    k  N k   k  0
 k n 1  k
n 1
N
2
N

一方    z nk    k   N    0
k 1  n 1

N
以上の2式より
k  k
(GM 90)
N
2
ここでは、  ( z k ) が古い  k ,  k ,  k を反映して計算された値であった。
それを固定して、loglikelihoodを最大化する新たな  k ,  k ,  k を求めているわ
けだ。
以上の(GM70)(GM80)(GM90)によって  k ,  k ,  k の更新式が求められた。
log likelihood (GM50)が収束しなければEstepに戻る
EM法の応用:不完全観測データの場合
のモデル推定:多項分布の場合
 観測値が不完全な場合としては、複数の確率変数があるの
に観測されるのは、K個の和だけなどという場合。
 例題:母集団が次の多項分布である場合にN個の観測値か
らパラメタ-を推定する問題を考える。
N!
p  x1 , x2 , x3 , x4 , x5  
1x1 2x23x3 4x45x5
x1! x2! x3! x4! x5!
 1  1 1  
ただし 1 ,  2 ,  3 ,  4 , 5    , ,
,
,  (1)
4 4
2 4 4
 観測値としては、x1,x2,x3,x4,x5ではなく、x1+x2に対応するyと
x3,x4,x5が得られていたとする。
 このため、観測値から直接にパラメタ- θを求められない。
 そこで、以下のstep1, step2を繰り返して θ を近似する。ただし
、 θの初期値を θ(0)とする。また、以下はk+1回目の繰り返し
とする。
 k を用いて x1,x2 ,x3,x4 ,x5の推定値を求める。
step1 既に求まっている
の近似値として  k が与えられていたとき
log
の対数尤度は次式
N!
 x1 log 2   x2  x3  x4  x5  log 4   x2  x5  log    x3  x4  log1   
x1! x2! x3! x4! x5!
  x2  x5  log    x3  x4  log1     const
x1 , x2 , x3 , x4 , x5のなす多項分布は次式
1 


4


 ただし、 y  x1  x2なので、
N!
1  
  
x1! x2! x3! x4! x5!  2 4 

y! 

x1! x2!  1  old

4
2
1
2
y
  old
x1





 4
 1  old
4
2
1   
  

4
 4

x1 , x2の分布は次式
x3
x4
x5
x2


y!  2

 
x1! x2!  2  old





x1
 old

 2  old



x2
step2 step1の結果を用いて Q関数


Q  |  old  Ex , x , x , x , x ,  old [ p ( x1 , x2 , x3 , x4 , x5 , )]
1
2
3
4
5
の新しい近似値  k  1を次式で求める
 k  1  arg max Ex , x , x , x , x ,  [ p( x1 , x2 , x3 , x4 , x5 , )]
old

1
2
3
4
5
具体的には


Q  |  old  Ex , x , x , x , x ,  old [ p ( x1 , x2 , x3 , x4 , x5 , )]
1
2
3
4
5
 E old  x2  x5  log    x3  x4  log1     const 
yは x1, x 2が各々確率
1
2
(*)
 old
,
4
1  old 1  old


2
4 2
4
である 2項分布で y  x1  x 2
n k
nk
2項分布  p 1  p  の期待値は np
k 
 old
 old
この場合には n  y, p 
 E x2   y
2   old
2   old


 old
(*)   y
 x5  log    x3  x4  log1     const
 2   old



 old
 y
 x5 
old
2   old
Q  | 

  x3  x4  0



1




  new  arg max Q  |  old 

y
 old
 x5
2   old
 old
y
 x3  x4  x5
2   old
準備:KL divergence
 相対エントロピー or Kullback-Leibler divergence or
KL divergence: KL(P||Q):分布PとQの類似性を測る尺度
KL( P || Q ) 

i
P ( xi ) log
P ( xi )
Q ( xi )
 KL(P||P)=0
 KL(P||Q)≠KL(Q||P)
 非対称なので擬距離
 対称性を持たせるために
SymmetricKL(P||Q)=(KL(P||Q)+KL(Q||P))/2 という尺度もある
 相互情報量:
I  x, y   KLP x, y  || P x P y    P x, y  log
x, y
この部分をpointwise mutual informationとして使うこと
もある
P  x, y 
P x P y 
変分ベイズ法(Variational Bayse:VB)
 観測データからのベイズ推定
 観測データ:X、未知パラメター:θ、
モデル構造:M、潜在変数集合:Z
M
Z
θ
X
p( X , Z , , M )  p( X , Z |  , M ) p( | M ) p( M )
p( X , Z |  , M ) p( | M ) p( M ) p( X , Z |  , M ) p( | M ) p( M )
p( , Z , M | X ) 

p( X )
  p( X , Z , , M )dθ
M
Z
p( | X )   p( , Z , M | X )
M
Z
 新たなデータxの事後予測分布
p ( x | X )   p ( x |  ) p ( | X )d
この積分の計算が
困難
計算の困難さの問題点を詳しくいうと
 ベイズ推定は、最尤推定と異なり、未知データの予測値では
なく予測分布を求める
 教師データが少ない場合でも、汎化能力の高い予測器が作
れる
 ただし、P(Z|x,θ,M)、xは1個の観測データ(L次元)で、ベイズ
の定理で次のように変形するが
P ( Z | x,  , M ) 
 右辺 P(x,Z|θ,M)は
P ( x, Z |  , M )
 P ( x, Z |  , M )
Z
 xはZの成分(z1,z2,…,zK)に組み合わせで依存しているため、次式のよ
うに分解できない。
K
L
P(x, Z |  , M )   P( xl , zk |  , M )
k 1 l 1
 よって、Kが大きくなるとZの成分の組み合わせの数が膨大になり計算
が困難
変分ベイズ法の考え方
問題であった期待値の計算を近似的に確率的シミュレーションで
解くMCMC法が有力。
ただし、MCMCは計算が膨大。数理モデルを工夫し計算を効率
化する方法として変分ベイズ法
EM法では、Q(θ)を最大にするθを求めた。
VB法では、θの値を最大化の対象にするのではなく、θの分布の
形そのものを求める変分法
変分ベイズ法のトリック
L( X )  log P( X )  log   p( X , Z ,  , M )d
M
Z
p( X , Z , , M )
 log   q( Z ,  , M )
d
q( Z ,  , M )
M Z
Jensen の不等式
log( E[ x ])  E[log( x )]
より
p( X , Z , , M )
   q( Z ,  , M ) log
d  F ( q)
q( Z ,  , M )
M Z
  q( Z , , M )d  1
M
に注意。
Z
すなわち、 log P( X )はM ,  , Zに対して
周辺化しているので、
M ,  , Zに依存しない。
L( X )  F(q)
   q( Z ,  , M ) log p( X )d    q( Z ,  , M ) log
M
Z
M
Z
p( X , Z , , M )
d  (1)
q( Z ,  , M )
p( X , Z , , M )
 p( Z ,  , M | X ) だから
p( X )
p( Z , , M | X )
(1)     q( Z ,  , M ) log
d  KL( q || p )
q( Z ,  , M )
M Z
よって
L( X )  F(q)  KL( q || p )
 L(X)はqに依存しないから、F(q)を最大化すること
はKL(q||p)を最小化、すなわちpに最も似たqを求め
ることになる
 EMでは、qをP( X , Z |  (t ) ) と決め打ちしていたが、VB
では、より柔軟な決め方を行っている。
上記のqを求めるプロセスをいきなり計算することは困難なので、
パラメターの事前分布とにq関して以下の仮定をおく。
因子化の仮定
I
p(θ | M ) 
 p(
i
|M)
i 1
I
 q(
q( Z , θ , M )  q( M ) q( Z | M )
i
|M)
i 1
この仮定の下で、変分法を使えば、次に述べる変分ベイズ
法のアルゴリズムが導ける
しかし、この仮定が成立しないこともあるので、適用する問
題によっては注意が必要。
因子化の仮定下でのVBの導出 その1
I
F(q)    q( M ) q( Z | M ) qi | M  log
M
i 1
Z
d  d1  d I
p( X , Z |  , M ) p( | M ) p( M )
d
I


q( Z | M )  qi | M q( M )
 i 1

I
p( | M )   p i | M  に注意すると
i 1
I

p( X , Z |  , M ) 


q
(
Z
|
M
)
q

|
M
log
d 

i
 
q( Z | M )
i 1
Z

I


p i | M 
F(q)   q( M )    qi | M  log
di



q

|
M
M
i
 i 1



p( M )
 log

q
(
M
)


I
 q( Z | M )  q
Z
i 1
i
| M d1  d I  1
 q( Z | M )  1
 q | M d  1
Z
j
j
因子化の仮定下でのVBの導出 その2
モデル Mが与えられたときの
 q( Z | M )  1
Zの最適な分布 q( Z | M )は
という条件下で F(q) を q( Z | M )に対して最大化
Z
すなわち、次の J [ q( Z | M )]の極値問題を解く。
I
J [ q( Z | M )]   q( M )  q( Z | M ) qi | M  log
M
Z
i 1
p( X , Z |  , M )
d
q( Z | M )


    q ( Z | M )  1
 Z

J [ q( Z | M )]が q( Z | M )の微分を含まないので
J [ q( Z | M )]
0
q( Z | M )
J [ q( Z | M )]
0

ただし、 q( Z | M )は と無関係なので、 J [ q( Z | M )]の右辺  の前に出、
次のようになる。
因子化の仮定下でのVBの導出 その3
J [q ( Z | M )]

q ( Z | M )
I




 q ( Z | M )   q  i | M  log p ( X , Z |  , M )d 
i 1
 q ( M )  Z

I

M
   q ( Z | M ) log q ( Z | M )  q i | M d 
 i 1


q ( Z | M ) 
Z






    q ( Z | M )  1


 Z

 I

  q  i | M  log p ( X , Z |  , M )d

 i 1

  q( M )
0
I
I
M
 log q ( Z | M )
q i | M d    q  i | M d   




i 1
i 1
 なぜなら  q ( M )  1だから   q ( M )  とおけるので
M
M
因子化の仮定下でのVBの導出 その4
I
さらに
  q
i
| M d  1だから、結局
i 1
I
  q
i
| M  log p( X , Z |  , M )d  1    log q( Z | M )
i 1
 I

 q( Z | M )  C exp   q i | M  log p ( X , Z |  , M )d 
 i 1

もし、 が確率変数でなく確定 していると、
q i | M     i   
はある iに対する  iの値(スカラー)
よって、 q( Z | M )  Cp( X , Z |  , M )
C
1
 p( X , Z |  , M )
Z
 (*)
この式の は(*)の からなるベクトル
因子化の仮定下でのVBの導出 その5
モデル Mが与えられたときの Z の最適な分布 q(i | M ) は
次の J [ q(i | M )](  下の式の  の内部)の極値問題を 解く。
I






q
(
Z
|
M
)
q

|
M
log
p
(
X
,
Z
|

,
M
)

log
q
(
Z
,
M
)
d


i

 

Z
i

1
  q( M ) 


I


p

|
M
i
M


   qi | M  log
di
q(i | M ) 


qi | M 
i 1






q

|
M
d


1
i


 i
 
 
   q( Z | M )   q j | M log p( X , Z |  , M )d i di 

i j
   0
  q( M )   Z
p i | M 
M


  log
 1di
  di


qi | M 


 f ( )d
i
よって
i
 0 だから f (i )  0
 q( Z | M )  q j | M log p( X , Z |  , M )d i  log
Z


i j
qi | M 
C
p i | M 


qi | M   C p i | M  exp  q( Z | M )   q j | M log p( X , Z |  , M )d i 
i j
Z

変分ベイズ法のアルゴリズム
 初期化として、以下の初期分布を設定
{q( i | M ) old },{ p( i | M ) old }
i  1,..., I
 反復計算 以下を収束するまで繰り返す。
 VB-E step
q( Z | M ) new  C exp(   q( | M ) old log p( X , Z |  , M )d )  C exp( EM , Z , [log p( X , Z |  , M )])
M
Z
 VB-M step
q( i | M ) new
 C ' p( i | M ) exp(   q( Z | M ) newq( \ { i } | M ) old log p( D, Z |  , M )dZd \ { i })
M
Z
 C ' p( i | M ) exp( EM new , Z , \{ }old [log p( X , Z |  , M )])
i
変数 newを変数 oldとする。
 \ {i } はθの構成要素からθiを除いた残りを意味する
変分ベイズ法再考
EMの再考を思い出して比較してみる。
VB  Estep :
q( Z | M )new  C exp( EM ,Z , [log p( X , Z |  , M )])
VB  Mstep :
q(i | M )new  C ' p(i | M ) exp( EM new ,Z , \{ }old [log p( X , Z |  , M )])
i
変数 newを変数 oldとする。
 P(Z,X| θ,M)を θold を固定してZ, θ,Mで期待値をとる
ことによって、Z,θ,Mに関する情報を教師データZか
ら集めて再度推定することを繰り返しての良い推定
値を求めている。
 ただし、因子化仮定によってθiを別々に更新してい
る。だから解析的に更新式が求まる場合もあるわけ
だ。
変分ベイズ法の例題
 1変数正規分布p(X| μ , τ )のパラメターを観測データから推
定。ただし、平均値μ(これを潜在変数Zとみなす) 精度
(=分散の逆数)τ(これをパラメターθとみなす)の事前分
布は以下のように与えられているとする。
 p(τ|a,b )を定義するガンマ関数は、パラメターa,bによっていろ
いろな分布形になるので、a,bを変化させて適した分布を得る
目的でVBの事前分布として使うことが多い。
  
p ( X |  , )   
 2 
N
2
  N
2
exp   ( xi   ) 
 2 i 1

p(  |  )  N (  | 0, (0 ) 1 )
b0 0 a0 1
p( | a0 , b0 )  Gamma( | a0 , b0 ) 
exp( b0 )
(a0 )
a
ガンマ分布
0.5
10
 factorizedな変分近似の事後分布を
q(μ , τ )=q μ(μ)q τ(τ)とすると、以下のようにVB-Eステップ、VB-Mステッ
プの計算ができる。左辺のqは更新した結果とする。
 VB-E:ここでは、内部変数μ,λを更新する。
q   (  )  C0 exp E [log p ( X |  , )  log p (  |  )]
 E[ ]  N
2
2 
 C1 exp 
  ( xi   )  0 (    0 ) 
 2  i 1

E[log(τ/2π)N/2]
N
N


2
2 
などはμには



x



x


0 0
i
0 0
i 
関係しないの
 E[ ]
i 1
i 1
  2  2




C
exp



N



1
0
で定数とみな


2


N


N
0
0

す





よって
q (  )  N (  |  N , N1 )
N 
0  0  Nx
0  N
N  0  N E[ ]


 VB-Mステップ
q ( )  C0 p( ) exp E [log p( X |  , ( ) 1 )  log p(  |  )]

N 1

( a0  1) log   b0  2 log  
 C1 exp 

N

2
2
 E [ ( xi   )  0 (   0 ) ]

 2
i 1
 Gamma分布が指数分布族(事前分布と事後分布が同じタイプ
だから、以下のように推論できる。
この結果を Gamma( | a N , bN )に対応させると以下の 更新式が得られる
N 1
a N  a0 
2
N
1
bN  b0  E [ ( xi   ) 2  0 (   0 ) 2 ]
2
i 1
 こうして p ( ) を定義するGamma分布のパラメターが更新さ
れた。
 以下、同様にVB-E,VB-Mを収束するまで繰り返すことになる。
Expectation Propagation
 変分ベイズ法ではKL(q||p)を最小化した。しかし、p
に近いqを求めればよいのだからKL(p||q)を最小化
する方法もありうる。
 これがExpectation Propagation:EP法。
EP法の背景
 KL(p||q)を最小化するqをpから求める。
 qはexponential family
q(z )  h(z ) exp( ηTu(z)  a( η))
KL( p || q)  a( η)  ηT E p ( z ) [u(z) ]  const
then to minimize
KL( p || q)

a( η)
 E p ( z ) [u(z) ]
η
 一方、exponential familyであるqにおいて
qの積分が1だから
T

h
(
z
)
exp
η
u(z)  a ( η) dz=1

a ( η)
h( x ) exp( ηT u( x )  a ( η))dx   h( x ) exp( ηT u( x )  a ( η))u( x )dx  0

η
a ( η)
ゆえに
 Eq ( z ) [u(z) ]
η
より

あわせると
Eq ( z ) [u(z) ]= E p ( z ) [u(z) ] となり
momentが保存されている。
EP法
 確率変数θをパラメターとするとき、第i番目の
データの出現確率がfi(θ)とする。
 このとき、観測データDの結合確率は p( D, )   f i ( )
i
 そこで、狙いは事後分布p(θ|D)を近似する分
布qを以下のように求めること。
1
q( ) 
Z

i
~
f i ( )
EP法の処理手順
1.
~
fi ( )
を全て初期化。
q( ) 
2. 事後分布の近似を以下ように設定

~
f i ( )
i
3. 収束するまで以下を繰り返す
(a)改善すべき ~f j ( ) を決める
(b)これを次のようにして取り除く
q( )
q ( )  ~
f j ( )
\j
(c)良い十分統計量(exモメント)が保存されるような新し
い事後分布 q new ( )  q \ j ( ) f ( ) を求め、正規化定数も
j
次のように求める。
Z  q \ j ( ) f ( )d
j

j
(この計算を解析的に行うので、けっこう面倒である。)
(d)以下の更新を行う。
~
q new ( )
f j ( )  Z j \ j
q ( )
4. 得られたモデルの近似度を評価する。
p( D ) 

i
~
f i ( )d
EP法の例:グラフィカルモデル
 左側の構造のグラフィカルモデルを、右側のように
分解する。
~
~
~
f a1
x1
x2
fa
x3
x1
~
fa2
fb2
x2
fc
x4
x3
~
fc2
fb
~
fc4
~
f a1
x4
fb3
p(x)  f a ( x1 , x2 ) f b ( x2 , x3 ) f c ( x3 , x4 )
~
~
~
q(x)  f a ( x1 , x2 ) f b ( x2 , x3 ) f c ( x3 , x4 )
ここで、前頁の右のグラフのように分解してみると
~
~
~
~
~
~
q( x )  f a1 ( x1 ) f a 2 ( x2 ) f b 2 ( x2 ) f b 3 ( x3 ) f c 2 ( x2 ) f c 4 ( x4 )
~
~
~
ここでEPアルゴリズ
ムにおいて f b ( x2 , x3 )=f b 2 ( x2 ) f b 3 ( x3 )を選ぶ。
すると
~
~
~
~
q ( x )  f a1 ( x1 ) f a 2 ( x2 ) f c 2 ( x2 ) f c 4 ( x4 )
\b
そして、更新した
f bをかけて
~
~
~
~
pˆ ( x )  q \b ( x ) f b ( x2 , x3 )  f a1 ( x1 ) f a 2 ( x2 ) f b ( x2 , x3 ) f c 2 ( x2 ) f c 4 ( x4 )
ここで、 q new ( x )は、 KL( pˆ || q new )を最小化するものとし
て得る。
KL( pˆ || q new )を最小化するものは以 下のようにして求める
。
KLを最小化するには
M

KL( p || q)    p( Z)  log qi ( Z i )dZ
 i 1

M

 const (じつは  p( Z)  log pi ( Z i )dZすなわち、 qに無関係な pのエントロピー )
 i 1

これを qiについて最小化するの だが、 Lagrangue未定乗数法により
M
h  KL( p || q)   (  qi  1)
i 1
0
h
1
   p( Z) dZi

q j
qi ( Z i )
i j
 p ( Z )  dZ
i j
ここで
i
 p( Z
 p(Z j )  qi ( Z i )
j
)  1   qi ( Z i )かつ  qi ( Z i )  1より   1
ゆえに最適な qすなわち q* ( Z i )   p( Z) dZi
i j
周辺化した p
よってqnewは更新したpを周辺化すればよい
関連する周辺化された p は以下の通り
~
pˆ ( x1 )  f a1 ( x1 )
~
~
ˆp ( x2 )  f a 2 ( x2 ) f c 2 ( x2 ) f b (x2 , x3 )
pˆ ( x3 )  
x3

~
~
f a 2 ( x2 ) f c 2 ( x2 ) f b ( x2 , x3 )
x2
~
pˆ ( x4 )  f c 4 ( x4 )
~ ~
よって、 q において更新された f b 2 , f b 3は以下のとおり。
~
f b 2 ( x2 )   f b (x2 , x3 )
new
x3

~
~
~
f b 3 ( x3 )   f a 2 ( x2 ) f c 2 ( x2 ) f b ( x2 , x3 )
x2
以下、同様にfa , fb , fcについてこれらの操作を収束するまで繰り返す。