EMアルゴリズムとは

確率モデルのパラメータを最尤法に基づいて推定する手法のひとつであり、 観測不可能な潜在変数に確率モデルが依存する場合に用いられる。

  • EMアルゴリズムは反復法の一種であり、期待値(expectation, E) ステップと最大化 (maximization, M)ステップを交互に繰り替えすことで計算が進行する。 Eステップでは、現在推定されている潜在変数の分布に基づいて、モデルの尤度の期待値を計算する。 Mステップでは、E ステップで求まった尤度の期待値を最大化するようなパラメータを求める。 M ステップで求まったパラメータは、次の E ステップで使われる潜在変数の分布を決定するために用いられる。
  • EM アルゴリズムは の尤度を単調に増加させるので、安定的に解に近づける。
  • 混合分布や隠れマルコフモデルやグラフィカルモデルの学習に応用されている。音声認識や金融、移動体追跡などの逐次推定に用いられている。移動体の場合、複数対象のモデル化の問題を各対象の確率分布の混合分布の最尤推定問題としてとらえることで解決する。

例:混合正規分布でy=ax+uのモデルの場合の最尤法

同一正規分布の場合観測値yをm次元ベクトル、入力(隠れ状態) xをn(n>m)次元ベクトルとする。データセット{xt, yt t= 1,..., T} から、パラメータa を求める問題を考える。ただしu は分平均0、分散σ^2Iのガウスノイズである。このとき、y は中心位置ax、共分散σ^2Iの正規分布N(y|ax,σ^2) にしたがう。I:単位行列とする。

P(y|x,a)=N(y|ax,σ^2)=SQRT(2πσ^2)^(-m/2)EXP[-||y-ax||^2/)2σ^2)]

中心位置が入力x に依存しているため上式で与えられるy の分布はx の条件付きである。対数尤度を取ると

LogP(y|x,a)=(-1/2σ^2)||y-ax||^2 + aに無関係な定数

で表わされる。 観測データy1;t=(y1,y2,...yt) x1;t=(x1,x2,...,xt)の生起確率は、独立なガウス分布に従うことから、y1;tの正規確率は上記の確率密度関数ので表わされるので、対数尤度関数は和の形になり

P(y1,y2,...,yt|a)=P(y1|a)P(y2|a)....P(yt|a)

であることから、観測データ所与のもとでの尤度関数は

L(a)=(-1/2σ^2) Σ||yt-axt||^2 + aに無関係な定数

このことから、L(a) をa について最大化することと、二乗誤差E=Σ||yt-axt||^2を最小化することは等価である。したがってこの場合の最尤推定法は最小二乗法と等価である。

混合分布とは

さて、ここで、M個の正規分布、P(X|i)=N(x|μi,σi^2),i=1,Mがあり、この混合分布を示す。 混合モデルとは、M個の分布の重ね合わせである。

P(X|θ)=ΣP(x,i|θ)=ΣP(X|i,θ)P(i|θ)  i=1~Mの合計
  • 混合比率
    P(i|θ)=vi Σvi=1
  • 生成分布
    P(x|i,θ)=N(x|μi,σi^2)
  • 未知パラメータ
    θ=(v1,..,vm;μ1...,μm;σ1,...,σm)
    である。 この時、パラメータが所与のもとでのXの分布は、
    P(x|θ)=Σp(x,i|θ)=ΣP(x|i,θ)p(i|θ)=Σvi N(x|μi,σi^2)
    となる。
  • 注意すべき点は、混合正規分布において、クラスタの指標i は隠れ変数になっていることである。すなわち、データを観測した点では、上式に対応する生成過程は観測できないため、データがどのクラスタから出てきたかは分からない。できるのは推測することだけである。これを隠れ変数をもつ推定問題という。
    Gauss-Mixtue.JPG

対数尤度関数と最尤法

データセットx1;t=(x1,x2,...,xt)が、観測された場合、 尤度は

P(x1;t|θ)=P(xt,x1;t-1|θ)=P(xt,x1;t-1|θ)
         =P(xt|x1;t-1,θ)P(x1;t-1|θ)
         =P(xt|x1;t-1,θ)......P(x3|x1;2,θ)P(x2|x1;1,θ)
         =πP(xt|x1:t-i,θ) i=1~tの積
          =π {πP(k|x1:t-i,θ)P(xt,k|x1:t-i,θ)}

であるので、その対数尤度は

L(θ)=LogP(x1;t|θ)=Σlog P(xt|x1;t-i,θ) = ΣLog {Σp(k|x1:t-i,θ)P(xt,k|x1:t-i,θ)}
内側のΣは隠れ変数k=1~Mについて、外側のΣはt=1~tについての合計である。

である。

これをを最大化するθを求めるのが最尤法である。

  • ここで、各サンプル時点の観測値xtが、どの生成分布から出てきたかがわかれば、要素毎に推定を行えば良いので、隠れ変数iとサンプルxtをつけくわえた完全変数(x1;t,i1;t)の分布が重要である。 完全変数(x1;t,i1;t)の分布は、ベイズの定理より
    P(X1;t,i1;t|θ)=P(X1;t|i1;t,θ)P(i1;t|θ)
    で表わされることに留意する。
  • 期待値計算のEステップでは、データX1;tを観測したという条件でのP(j|x1;t,θ)が重要な役割を果たす。
    P(j|x1;t,θ)=P(x1;t,j|θ)/P(x1;t|θ)=P(X1;t|j,θ)P(j|θ)/Σp(x1;t,j|θ)

逐次表示でない場合、独立性を仮定してP(x1;t|θ)=P(Xt|θ)P(xt-1|θ)...P(x2|θ)p(x1|θ)として、L(θ)=ΣP(x1;t|θ)=Σlog Σp(xt,i|θ) から必要条件である停留条件はθで偏微分して0とおいて次のように計算できる。

[ΣΣ∂p(xt,i|θ)/∂θ]/Σp(xt,j|θ)=0
ΣΣ{p(i|Xt,θ)∂p(xt,i|θ)/∂θ}=0          a式
ただし、p(i|Xt,θ)=P(xt,i|θ)/Σp(xt,j|θ)    b式

p(i|Xt,θ)は、データx1;t を観測した際の、隠れ変数i の事後確率(posterior) と呼ばれる。

それぞれのパラメータについて、最尤推定を行ってみる。Xベクトルの次数をdとする。

  • 平均値μj,j=1,Mの最尤推定値
    ∂L(θ)/∂μj=Σ{[P(j|θ)/P(xt|θ)]∂P(xt|j,θ)/∂μj}=0
    Σ{[P(j|θ)/P(xt|j,θ)/P(xt|θ)](μj-xt)/σj^2}=0
    Σ{P(j|xt,θ)(μj-xt)/σj^2}=0
    ゆえに
    μj*=ΣP(j|xt,θ)xt/ΣP(j|xt,θ)
  • 標準偏差σj,j=1,Mの最尤推定値
    ∂L(θ)/∂σj=Σ{[P(j|θ)/P(xt|θ)]∂P(xt|j,θ)/∂σj}=0
    Σ{P(j|xt,θ)[(d/σj)-||xt-μj||^2/σj^3]}=0
    ゆえに
    σj*^2=(1/d){ΣP(j|xt,θ)||xt-μj||^2/ΣP(j|xt,θ)
  • 混合比率P(j)の最尤推定値 まずP(j)=vjを次の、Sofmax fynction で表わされるとすれば、綺麗な解を得る。
    P(j)=EXP(γj)/Σ EXP(γj)
    停留条件は
    ∂L(θ)/∂γj=Σ{∂L(θ)/∂P(k)・∂p(k)/∂γj}=0
    Σ{P(j|xt,θ)-P(j)}=0
    ゆえに
    P(j)*=(1/T)ΣP(j|xt,θ)
  • 以上より 最尤法の解は、P(j|xt,θ)を使って計算できる。そこで、これを計算する必要があるが、そのためには、P(xt|j,θ)とP(j|θ)が分からないと求められない。
    P(j|xt,θ)=P(xt|j,θ)P(j|θ)/P(xt|θ)

EMアルゴリズム

一般に、上記の方程式は、非線形であり、EMアルゴリズムでは、式の形に注目して、以下のような繰り返しアルゴリズムを考える。

  • E (Expectation) ステップ
    • 各データxt に対して、現在のパラメータの推定値 θ*t を用いて隠れ変数i の事後確率p(i|Xt,θ*t)をb式より計算する。
    • 隠れ変数を含む(完全)データセットに対する対数尤度Σlog P(xt,i|θ)の、隠れ変数の予測事後確率についての期待値Q(θ|θ*t)を計算する。
      Q(θ|θ*t)= ΣΣp(i|Xt,θ*t)logp(xt,i|θ)
  • M (Maximization) ステップ
    • 期待対数尤度Q(θ|θ*t) をパラメータ θ について最大化する。すなわち a式を近似計算する
      ΣΣp(i|Xt,θ*t)dp(xt,i|θ)/dθ=0
  • 繰り返し
    • 求められたパラメータをθ*t+1 としてEステップ に戻る。Eステップ とMステップ の繰り返しが収束すれば終る。

前記の、パラメータθ=(γ1,...,γm;μ1,...,μmσ1,...,σm)の正規混合分布の場合を計算してみる。但し、Xtベクトルの次元をd次元とする。

P(x,i|θ)=vi・N(μi,σi^2)=vi・SQRT(2πσi^2)^(-d/2)EXP[-(x-μi)^2/2σ^2]
P(x|θ)=Σp(x,i|θ)=ΣP(x|i,θ)p(i|θ)=Σvi・N(x|μi,σi^2)
Vi=EXP(γi)

である。

  • Eステップ
    P(i|xt,θ*t)=P(xt,i|θ*t)/P(xt|θ*t) ベイズの定理より
               =vi*t・N(μi*t,σi*t^2)/Σvi*t・N(μi*t,σi*t^2)
  • Mステップ
    LogP(x,i|θ)=LogVi-(d/2)log(2πσi^2)-(x-μi)^2/(2σi^2)
    であるので
    Q(θ|θ*t)=ΣΣP(i|xt,θ*t)LogP(xt,i|θ)
            =ΣΣ{P(i|xt,θ*t)・[LogVi-(d/2)log(2πσi^2)-(x-μi)^2/(2σi^2)]}
    前に求めた最尤推定値を用いて更新する。
    平均の更新:μj*t+1=ΣP(j|xt,θ*t)xt/ΣP(j|xt,θ*t)
    分散の更新:σj*^2t+1=(1/d){ΣP(j|xt,θ*t)||xt-μj||^2/ΣP(j|xt,θ*t)
    混合比率の更新:P(j)*t+1=(1/T)ΣP(j|xt,θ*t)
    • 求められたパラメータをθ*t+1 としてEステップ に戻る。Eステップ とMステップ の繰り返しが収束すれば終る。

対数尤度を安定的に減少させるEMアルゴリズム

EとMのステップを1回繰り返すことで、対数尤度が増加する。最初の節で述べたように対数尤度関数は、マイナスの誤差二乗和に比例する形であるので、誤差二乗和が安定的に減少しているかチェックできる。そこで E=-L(θ)とおいて、パラメータの更新によって、これがどの程度減少するかを調べてみよう。

Et+1-Et=L(θt)-L(θt+1)
       =-ΣLogp(xt+1|θt+1)+ΣLogp(xt|θt)
       =-ΣLog{p(xt+1|θt+1)/p(xt|θt)}
       =-ΣLog{Σ[P(j|θt+1)p(xt+1|j,θt+1)P(j|xt,θt)]/[P(xt|θt)p(j|xt,θt)]}

Jensenの不等式をここで使ってみる。

λj>o,Σλj=1 ---->log[Σλjxj] > Σλjxj

より

Et+1-Et<-ΣΣP(j|xt,θt)log{P(j|θt+1)p(xt+1|j,θt+1)/[P(xt|θt)p(j|xt,θt)]}=Q

となるので、Qを最小化するところの、Mステップによって、Qの値に応じて、誤差二乗和が減少することがわかる。EMアルゴリズムは、安定的に誤差二乗和を減少させる最適化アルゴリズムになっていることが理解できる

例:隠れマルコフモデル

隠れマルコフモデルは,観察データの列y1:T={y1,y2,...,yT} をモデル化するために用いる.なお,各時刻t で観察されるデータyt は,全部でM 種類あるシンボルw1,...,wM のうちのひとつである.つまり,モデル化したい観察データは,長さT のシンボル列である。

簡単のため、2種類のシンボル、(降雨あり、降雨なし)を考える。この場合、観測データーは、毎日の降雨の有無の観測結果である。そして、隠れ状態としては、(晴、曇り、雨)を考える。

  • 2つの仮定を置く。
    • 1. 時刻t において観察されたデータyt は,隠れ状態st によって生成されるとする.なお,どの隠れ状態も,全部でK 種類ある状態e1,..., eK のうちのひとつである.
    • 2. 隠れ状態の列s1:T は,一次マルコフ過程によって生成されるとする.つまり,時刻t における隠れ状態st が,e1,...., eK のいずれの状態となるかは,直前の隠れ状態st¡1 にのみ依存する.
    • 天気の場合は隠れ状態の数 (晴、曇り、雨)のK=3である。
    • eiは、i要素のみ1で他は0の単位ベクトルで表わす。晴は(1,0,0)'。曇りは(0,1,0)。雨は(0.0.1)

初期の天気の隠れ状態の確率をP(s1)で表わし、St-1からStに推移する確率が、1次のマルコフ連鎖で表わされとしよう。この時、

Prob(St=ej|St-1=ei)=aij
Σaij=1 i=1,nの式 Σはjに関する合計。

である。これらの確率をまとめてA = {aij} 行列であらわす。 また、隠れ状態がeiのとき、シンボルwjが表われる確率を、cijと書こう。これらの確率をまとめてC = {cij} 行列であらわす。

Prob{yt=wj|St=ei}=cij
Σcij=1 j=1,2の式 Σはiに関する合計

そして、最初の隠れ状態がek となる確率P(s1 = ek) を,πk と書き、これを初期状態π0=(π1,π2,π3)と呼ぶ。

仮定より、長さT の完全データ(つまり,観察されたデータも観察されていないデータも含めたすべてのデータ)の尤度は,次のように書くことができる.

P(s1:T , y1:T ) = P(s1)P(y1|s1)......P(yT|sT)

以上すべては,隠れマルコフモデルのパラメータである.つまり,θ = (A; C; π0) となる。 パラメータの数は、K^2 + KM + K 個(k=3,M=2なので合計)のパラメータの値を定める必要がある

参考


添付ファイル: fileGauss-Mixtue.JPG 500件 [詳細]

トップ   差分 バックアップ リロード   一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2010-10-29 (金) 08:58:00 (3302d)