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}
として計算できます．条件付き期待値については，こちらの記事(条件付き多変量正規分布)をご覧ください．
アルゴリズムは以下のような反復ステップで構成されています．
- A,Dの初期値を発生させる．
- $E[f|x_i],E[ff^T|x_i]$を計算する
- (1),(2)のようにA,Dを更新する
- 対数尤度を計算し，収束していれば終了．していなければ2に戻る．
これをプログラムにしたのが次章のコードです．
4. サンプルコード
今回も論文に書いてある通り(=数式をそのままプログラムにしたandどこかを工夫して計算スピードをあげたりしてない)に実装しています．データはBig5と呼ばれる心理学のデータで，ここからダウンロードしています．
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アルゴリズム徹底解説