Edited at

EMアルゴリズムによる因子分析のパラメータ推定


1. 雑要約

 今回は因子分析の確率モデルについて説明した後,EMアルゴリズムによるパラメータ推定のpython3での実装例を紹介します.Zoubin Ghahramani and Geoffrey E. Hinton(1996).The EM algorithm for mixture of factor analyzersのAppendix Aのところを参考にしています.


2. 記号等準備

 以下に今回使用する記号を先に整理しておきます.


  • x: はデータを表します.($p\times 1$)

  • f:共通因子ベクトル($m\times 1$)を表し,$$f\sim N(0_m,I_m)$$に従うとします.ここでmは因子数です.

  • A: 因子負荷量行列($p\times m$),共通因子得点をどのように重み付けするかを表す行列.分析の際はこの重み付けパターンを読み取る.

  • e: 独自因子ベクトル($n\times 1$)を表し,$$e\sim N(0_p,D)$$に従います.ここで,Dは対角行列で,対角成分が特に独自分散と呼ばれています.

EMアルゴリズムについてはこの記事(EMアルゴリズムについて徹底解説)が詳しかったので参考にしてください.


3. 因子分析(確率モデル)

 確率モデルの因子分析のモデル式は

\begin{eqnarray}

x&=& Af + e\
&s.t&\left\{
\begin{array}
1E[f]=0_m, \ E[e]=0_p,\\
V[f]=I_m, \ V[e]=D, \ Cov[f,e]=O_{m,p}
\end{array}
\right.
\end{eqnarray}

です.ここでポイントとなるのは,制約式の$V[e_i]=D^2$(対角行列)と$Cov[f,e]=O_{m,p}$の部分です.このモデル式は「データxが少数の要素を持つf($m\times 1$の共通因子ベクトル,$m\le p$)の重み付け和(=Aによる線型結合)+e(独自因子ベクトル)で説明される」ということを表しており,制約式のポイントとなる部分は,独自因子はデータの各変数に固有の影響を持つ(←対角行列であることから.ここがPCAとの差異です),共通因子と独自因子は無相関であるということを表しています.通常因子分析で推定されるパラメータは因子負荷量行列Aと独自分散行列Dです.xの分布は$$x\sim N(0_p,AA^T+D)$$であり,そこから導出されるxとfの同時分布は

\begin{eqnarray}

\begin{pmatrix}
f\\
x
\end{pmatrix}\sim N
\begin{pmatrix}
\begin{pmatrix}
0_m\\
0_p
\end{pmatrix},
\begin{pmatrix}
I_m & A^T\\
A & AA^T+D
\end{pmatrix}
\end{pmatrix}
\end{eqnarray}

これらの情報を用いてEMアルゴリズムを実行していくのですが,$x_i(i=1,\cdots,n)$はデータとして与えられているので,モデル式を$e=x_i-Af$と変更して,$$e\sim N(0,D)$$の対数尤度関数

\begin{eqnarray}

\log p(e|x_i,A,D)&=&\log\prod_{i=1}^n\frac{1}{(2\pi)^{\frac{p}{2}}|D|^\frac{1}{2}}\exp[-\frac{1}{2}(x_i-Af)^TD^{-1}(x_i-Af)]\\
&=&\sum_{i=1}^n [\log(2\pi)^{-\frac{p}{2}}+\log|D|^{-\frac{1}{2}}-\frac{1}{2}(x_i-Af)^TD^{-1}(x_i-Af)]\\
&=&-\frac{np}{2}\log 2\pi - \frac{n}{2}\log|D|-\frac{1}{2}\sum_{i=1}^n(x_i-Af)^TD^{-1}(x_i-Af)]\\
&=& const- \frac{n}{2}\log|D|-\frac{1}{2}\sum_{i=1}^n(x_i^TD^{-1}x_i-2x_i^TD^{-1}Af+f^TA^TD^{-1}Af)
\end{eqnarray}

を利用します.今回のモデルに対するQ関数は以下の期待値で定義されます.

\begin{eqnarray}

Q&:=&E[\log\prod_{i=1}^n(2\pi)^{-\frac{p}{2}}|D|^{-\frac{1}{2}}\exp[-\frac{1}{2}(x_i-Af)^TD^{-1}(x_i-Af)]]\\
&=&const-\frac{n}{2}\log|D|-\frac{1}{2}\sum_{i=1}^nE[(x_i^TD^{-1}x_i-2x_i^TD^{-1}Af+f^TA^TD^{-1}Af)]\\
&=&const-\frac{n}{2}\log|D|-\frac{1}{2}\sum_{i=1}^nE[(x_i^TD^{-1}x_i-2x_i^TD^{-1}Af+tr(f^TA^TD^{-1}Af))]\\
&&(\because f^TA^TD^{-1}Af\mbox{はスカラー})\\
&=&const-\frac{n}{2}\log|D|-\frac{1}{2}\sum_{i=1}^nE[(x_i^TD^{-1}x_i-2x_i^TD^{-1}Af+tr(A^TD^{-1}Aff^T))]\\
&=&const-\frac{n}{2}\log|D|-\frac{1}{2}\sum_{i=1}^n(x_i^TD^{-1}x_i-2x_i^TD^{-1}AE[f|x_i]+tr(A^TD^{-1}AE[ff^T|x_i]))]
\end{eqnarray}

AとDの最尤推定値の更新式はQ関数を微分して,

\begin{eqnarray}

\frac{\partial Q}{\partial A} &=& -\sum_{i=1}^n(-2D^{-1}x_iE[f|x_i]^T+2D^{-1}AE[ff^T|x_i])\\
\frac{\partial Q}{\partial D} &=& \frac{n}{2}D-\frac{1}{2}\sum_{i=1}^n(x_ix_i^T-2AE[f|x_i]x_i^T+AE[ff^T|x_i]A^T)
\end{eqnarray}

これらを0と置くことで

\begin{eqnarray}

\widehat{A} &=& (\sum_{i=1}^nx_iE[f|x_i]^T)(\sum_{l=1}^nE[ff^T|x_l])^{-1}\tag{1}\\
\widehat{D} &=& \frac{1}{n}diag(\sum_{i=1}^nx_ix_i^T-\widehat{A}E[f|x_i]x_i^T)\tag{2}
\end{eqnarray}

となります.$E[f|x]$と$E[ff^T|x]$という条件付き期待値は,

\begin{eqnarray}

E[f|x] &=& A^T(AA^T+D)^{-1}x\\
&&(\beta:=A^T(AA^T+D)^{-1}\mbox{と置く})\\
E[ff^T|x] &=& Var[f|x] + E[f|x]E[f|x]^T\\
&=& (I_m-A^T(AA^T+D)^{-1}A) + \beta xx^T\beta^T\\
&=& I_m - \beta A+ \beta xx^T\beta^T\\
\end{eqnarray}

として計算できます.条件付き期待値については,こちらの記事(条件付き多変量正規分布)をご覧ください.

 アルゴリズムは以下のような反復ステップで構成されています.


  1. A,Dの初期値を発生させる.

  2. $E[f|x_i],E[ff^T|x_i]$を計算する

  3. (1),(2)のようにA,Dを更新する

  4. 対数尤度を計算し,収束していれば終了.していなければ2に戻る.

これをプログラムにしたのが次章のコードです.


4. サンプルコード

 今回も論文に書いてある通り(=数式をそのままプログラムにしたandどこかを工夫して計算スピードをあげたりしてない)に実装しています.データはBig5と呼ばれる心理学のデータで,ここからダウンロードしています.


python, EMforFA.py

import numpy as np

import pandas as pd
import scipy.stats as sp

def gaussian(xvec,mean,cov):
bunbo = 1/((2*np.pi)**(xvec.size/2) * (np.sqrt(np.linalg.det(cov))))
shisu = -1/2 * np.dot((xvec-mean).T,np.dot(np.linalg.inv(cov),xvec-mean))
return bunbo * np.exp(shisu)

def likelihood(X,mean,cov):
sum = 0.0
for i in range(len(X)):
sum += gaussian(X[i],mean,cov)
like = np.log(sum)
return like

def EMforFA(data,m,starts=10): # m: number of factors
X = data.astype("float64")
n,p = X.shape
#data標準化
X = sp.stats.zscore(X)

likeAll = -1000
Lhat = np.zeros((p,m))
Psihat = np.zeros((p,p))

for j in range(starts):
turn = 0
print("---------- epoc: %d -------" % (j+1))
#Lambda,Psiの初期化
L = np.random.uniform(0.1,0.3,p*m).reshape(p,m)
for k in range(p):
L[k,np.random.randint(0,m-1)] = np.random.uniform(0.6,0.9,1)
Psi = np.diag(np.diag(np.identity(p)-np.dot(L,L.T)))
like = likelihood(X,np.zeros(p),Psi)

while True:
print(turn,": ",like)
## E-step
beta = np.dot(L.T,np.linalg.inv(Psi + np.dot(L,L.T)))
Ezx = np.zeros((m,n))
for a in range(n):
Ezx[:,a] = np.dot(beta,X[a,:].T)
Ez2x = np.zeros((n,m,m))
for a in range(n):
Ez2x[a][:,:] = np.identity(m) - np.dot(beta,L) + np.outer(np.dot(beta,X[a,:].T),np.dot(X[a,:],beta.T))
sumL = np.zeros((p,m))
sumL2 = np.zeros((m,m))
sumX = np.zeros((p,p))
sumZX = np.zeros((p,p))
for i in range(n):
sumL += np.outer(X[i,:].T,Ezx[:,i].T)
sumL2 += Ez2x[i][:,:]
sumX += np.outer(X[i,:].T,X[i,:])
#M-step
L = np.dot(sumL,np.linalg.inv(sumL2))
for i in range(n):
sumZX += np.dot(L,np.outer(Ezx[:,i],X[i,:]))
Psi = np.diag(np.diag(sumX - sumZX)) / n

like_new = likelihood(X,np.zeros(p),Psi)
diff = like_new - like
if diff < 10**(-6) and diff >= 0:
if likeAll < like:
Lhat = L
Psihat = Psi
break
if diff < 0:
print("diff is negative")
break
like = like_new
if np.isnan(like):
print("the calculated likelihood is nan")
break
turn += 1
return(Lhat,Psihat)

if __name__ == "__main__":
#データを準備
X = pd.read_excel("Big5.xlsへのpathをかいてください",sheet_name="Big5",usecols=np.arange(1,26))
X = X.as_matrix()

A,D = EMforFA(X,5)
print("A\n",A)
print("\nD\n",D)


(一年以上前に書いたコードを引っ張り出してきました)


参考文献

Zoubin Ghahramani and Geoffrey E. Hinton(1996).The EM algorithm for mixture of factor analyzers

EMアルゴリズム徹底解説