[mathjax]
因子分析はマーケティングや心理学・教育学において広く用いられる多変数解析の手法で、大量のデータの背後にある潜在的な共通因子を抽出することができます。SEEDATAではこの手法を様々なデータで適用し、定性リサーチのヒントとしています。本ページでは、因子分析において潜在変数と因子負荷量などの各種パラメタを推定する手法である、EMアルゴリズムについて説明します。
1. 因子分析における最尤推定
因子分析では、アンケートなどの観測データ(x)の背後に潜在変数(z)が存在し、以下のような関係があることを仮定しました。(因子分析：基本表式とパラメタ・変数の定性的意味、データ内の異質性に基づいたモデル)
ここでは実際に観測データが与えられた時に、因子負荷量や潜在変数などを求めることを考えます。個人ごとの差異については、潜在変数が従う分布で違いを組み込むような、以下のモデルを想定します。
begin{eqnarray}
boldsymbol{x}_{i} &=& boldsymbol{alpha}+Lambdaboldsymbol{z}_{i}+boldsymbol{epsilon} tag{2}\[0.5em]
boldsymbol{z}_{i} &=& N(mu_{z_i|x_i},Phi_{z_i|x_i}) tag{3}\[0.5em]
boldsymbol{epsilon} &=& N(0,Theta)tag{4}\[0.5em]
boldsymbol{alpha_i} &=& rm {const.}tag{5}
end{eqnarray}
(mu_{z_i|x_i})と(Phi_{z_i|x_i})は、それぞれ潜在変数の事後分布の平均値と分散であり、これらが個人ごとに異なるということです。以降は、簡単のために因子負荷量(Lambda), 測定誤差の共分散行列(Theta), データの切片(boldsymbol{alpha})をまとめてパラメタ(boldsymbol{theta})とおきます。また(N)人分の全体データを(boldsymbol{rm X}={x_i}(i=1, 2, ..,N ))とし、同じく潜在変数の全データを(boldsymbol{ rm Z})とします。分析のためには、以下の観測データ(boldsymbol{rm X})に対する以下の尤度関数を最大化するようなパラメタ(hat{boldsymbol{theta}})を求める必要があります。begin{eqnarray}
hat{boldsymbol{theta}}=
underset{boldsymbol{theta}}{rm argmax} p(boldsymbol{rm X}|boldsymbol{theta})=
underset{boldsymbol{theta}}{rm argmax} int dboldsymbol{rm Z} p(boldsymbol{rm X},boldsymbol{rm Z}|boldsymbol{theta})tag{6}
end{eqnarray}
しかし、式(6)にある(p(boldsymbol{rm X}|boldsymbol{theta}))を最大化するのは一般に困難です。EMアルゴリズムは、(p(boldsymbol{rm X}|boldsymbol{theta}))の最大化の代わりに、観測データと潜在変数の結合分布(p(boldsymbol{rm X},boldsymbol{rm Z}|boldsymbol{theta}))を用いた最尤推定を考えることで、このような問題に対処します。
2. EMアルゴリズムの基礎
2.1 基本的な枠組み
ここでは因子分析は一旦忘れて、観測データ(boldsymbol{rm X})と潜在変数(boldsymbol{rm Z})、パラメタ(boldsymbol{theta})があるような一般的な確率モデルについて考えましょう。以下では、観測データと潜在変数を合わせたデータ({boldsymbol{rm X},boldsymbol{rm Z}})を完全データ、観測データ({boldsymbol{rm X}})のみを不完全データと呼びます。
EMアルゴリズムは、先述の通り観測データが潜在変数に依存したモデルにおける、パラメタの最尤推定に用いられる手法です[1]。枠組みとしては、式(6)の尤度関数を直接最大化するのは難しいため、代わりに完全データの対数尤度(ln{p(boldsymbol{rm X},boldsymbol{rm Z}|boldsymbol{theta})})を考えることにします。潜在変数は観測できないため(p(boldsymbol{rm X},boldsymbol{rm Z}|boldsymbol{theta}))を取り扱うことはできませんが、事後分布(p(boldsymbol{rm Z}|boldsymbol{rm X},boldsymbol{theta}))が計算可能であるため、この元での完全データの対数尤度の期待値を最大化することでパラメタ推定を進めるというスタンスです。
完全データの対数尤度を考えること自体は、問題なく受け入れられるはずです。もともと手持ちの観測データ(boldsymbol{rm X})で進めようとしていた式(6)ではなく、潜在変数も含めた完全データでパラメタ推定をする方が容易であり、理にかなっていることが分かるでしょう。しかし、当然ながら完全データもまた観測不可能であるため、代わりにその期待値(E[ln{p(boldsymbol{rm X},boldsymbol{rm Z}|theta)}])を求めようとするところがEMアルゴリズムの特徴となります。計算においては、期待値の計算(E-step)とその最大化(M-step)を、交互にステップ発展させていきます。
2.2 表式による理解
上記の枠組みを式から理解してみましょう。不完全データの対数尤度を、潜在変数(boldsymbol{rm Z})の任意の確率密度関数(q(boldsymbol{rm Z}))を用いて次のように変形します。
begin{eqnarray}
ln{p(boldsymbol{rm X}|boldsymbol{theta})}&=&int dboldsymbol{rm Z} q(boldsymbol{rm Z})ln{p(boldsymbol{rm X}|boldsymbol{theta})}tag{7}\[0.5em]
&=& int dboldsymbol{rm Z} q(boldsymbol{rm Z})ln{frac{p(boldsymbol{rm X},boldsymbol{theta})}{p(boldsymbol{rm X},boldsymbol{rm Z},boldsymbol{theta})}frac{p(boldsymbol{rm X},boldsymbol{rm Z},boldsymbol{theta})}{p(boldsymbol{theta})}}tag{8}\[0.5em]
&=& int dboldsymbol{rm Z} q(boldsymbol{rm Z})ln{frac{p(boldsymbol{rm X},boldsymbol{rm Z}|boldsymbol{theta})}{p(boldsymbol{rm Z}|boldsymbol{rm X},boldsymbol{theta})}}tag{9}\[0.5em]
&=& int dboldsymbol{rm Z} q(boldsymbol{rm Z})ln{frac{p(boldsymbol{rm X},boldsymbol{rm Z}|boldsymbol{theta})}{q(boldsymbol{rm Z})} frac{q(boldsymbol{rm Z})}{p(boldsymbol{rm Z}|boldsymbol{rm X},boldsymbol{theta})}}tag{10}\[0.5em]
&=& int dboldsymbol{rm Z} q(boldsymbol{rm Z})ln{frac{p(boldsymbol{rm X},boldsymbol{rm Z}|boldsymbol{theta})}{q(boldsymbol{rm Z})}} + int dboldsymbol{rm Z} q(boldsymbol{rm Z})ln{frac{q(boldsymbol{rm Z})}{p(boldsymbol{rm Z}|boldsymbol{rm X},boldsymbol{theta})}}tag{11}\[0.5em]
&=& int dboldsymbol{rm Z} q(boldsymbol{rm Z})ln{p(boldsymbol{rm X},boldsymbol{rm Z}|boldsymbol{theta})}
-int dboldsymbol{rm Z} q(boldsymbol{rm Z}) ln{q(boldsymbol{rm Z})} + rm{KL}(q(boldsymbol{rm Z})||p(boldsymbol{rm Z}|boldsymbol{rm X},boldsymbol{theta}))tag{12}\[0.5em]
&=& L(q,boldsymbol{theta}) +rm{KL}(q||p))tag{13}\[0.5em]
end{eqnarray}
これを見ると、対数尤度は最終的に式(12)で表される3つの項に分解できることが分かります。まず第1項は分布(q(boldsymbol{rm Z}))に対する完全データの対数尤度の期待値であり、これがEMアルゴリズムで最大化したい対象となります。第2項は(q(boldsymbol{rm Z}))についてのシャノン・エントロピー、そして第3項は(q(boldsymbol{rm Z}))と潜在変数(boldsymbol{rm Z})の事後分布のKL情報量となっています。式(13)では第1,2項をまとめて(L)と表記しており[2]、KL情報量は負値を取らないため(L)は不完全データの対数尤度(ln{p(boldsymbol{rm X}|boldsymbol{theta})})の下限(lower bound)となり、ELBO(evidence lower bound)とも表現されます。今後の説明のために、式を(12)と(13)で分けて書いています。
式(12)から、完全データの対数尤度を最大化するためには、KL情報量をできるだけ小さくすれば良い事がわかります。対数尤度の最大化とKL情報量を0にすることを同時に行うのは難しいので、KL情報量を0にし対数尤度の期待値を計算するE-stepと、計算した期待値をパラメタ(boldsymbol{theta})の関数として最大化するM-stepに分けて反復計算していきましょう。このような手順は式(13)を元にすればイメージがつきやすく、様々なサイトや書籍において次のような図で説明されています。
まずE-stepについて考えます。KL情報量が0となる条件は、分布(q(boldsymbol{rm Z}))と事後分布(p(boldsymbol{rm Z}|boldsymbol{rm X},boldsymbol{theta}))が等しくなることです。これらの表式を得るために、パラメタ(boldsymbol{theta})を定数(boldsymbol{theta^{old}})として与え、M-stepへと移ります。いまKL情報量は0であるため、式(13),(12)から対数尤度は以下のようになります。
begin{eqnarray}
ln{p(boldsymbol{rm X}|boldsymbol{theta})}&=&L(q,boldsymbol{theta})\[0.5em]
&=&int dboldsymbol{rm Z} p(boldsymbol{rm Z}|boldsymbol{rm X},boldsymbol{theta}^{old})ln{p(boldsymbol{rm X},boldsymbol{rm Z}|boldsymbol{theta})}
-int dboldsymbol{rm Z} p(boldsymbol{rm Z}|boldsymbol{rm X},boldsymbol{theta}^{old}) ln{ p(boldsymbol{rm Z}|boldsymbol{rm X},boldsymbol{theta}^{old})}tag{14}\[0.5em]
&=&Q(boldsymbol{theta},boldsymbol{theta}^{old})+rm{const.}tag{15}
end{eqnarray}
式(14)の第2項はパラメタの推定に関係がないため、第1項の完全データ対数尤度の事後分布の元での期待値(Q)を最大化するような(boldsymbol{theta}^{new})を求めればよい事が分かります。書籍やサイトによっては、この期待値計算までをE-stepとするものもあります。
begin{eqnarray}
boldsymbol{theta}^{new}&=&underset{boldsymbol{theta}}{rm argmax} Q(boldsymbol{theta},boldsymbol{theta}^{old})tag{16}\[0.5em]
&=&underset{boldsymbol{theta}}{rm argmax} int dboldsymbol{rm Z} p(boldsymbol{rm Z}|boldsymbol{rm X},boldsymbol{theta}^{old})ln{p(boldsymbol{rm X},boldsymbol{rm Z}|boldsymbol{theta})}tag{17}\[1.0em]
boldsymbol{theta}^{old}&=&boldsymbol{theta}^{new}tag{18}
end{eqnarray}
このように(boldsymbol{theta}^{old})を更新する事で、反復計算により対数尤度の最大化を図る事ができます。
まとめ
このページでは、因子分析における最尤推定の導入と、EMアルゴリズムの基本的な手順について説明しました。ポイントとしては、不完全データの対数尤度を直接扱うのではなく、完全データの対数尤度の期待値を最大化することで推定を行うということであり、その具体的な計算もお分かりいただけたかと思います。このEMアルゴリズムは因子分析のみならず、混合ガウスモデルなど様々なモデルで用いられており、実際にはM-stepやE-stepではモデルごとにさらなる計算をする必要があるため、次回はそれらを確認します。
参考文献
[1] A. P. Dempster, et al, “Maximum Likelihood from Incomplete Data via the EM Algorithm”, J. Royal Stat. Soc.,39,1-38 (1977).
[2] C. M. Bishop, “Pattern Recognition and Machine Learning”, Springer, 450-453 (2006).
[3] Y. Chen, M. R. Gupta, “EM Demystified: An Expectation-Maximization Tutorial”,Department of Electrical Engineering University of Washington Seattle, WA 98195 (2010).