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

#1. 雑要約
　今回は因子分析の確率モデルについて説明した後，EMアルゴリズムによるパラメータ推定のpython3での実装例を紹介します．[Zoubin Ghahramani and Geoffrey E. Hinton(1996).The EM algorithm for mixture of factor analyzers](https://www.csie.ntu.edu.tw/~mhyang/course/u0030/papers/Mixture%20of%20Factor%20Analyzers.pdf)の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アルゴリズムについて徹底解説)](https://qiita.com/kenmatsu4/items/59ea3e5dfa3d4c161efb)が詳しかったので参考にしてください．
#3. 因子分析(確率モデル)
　確率モデルの因子分析のモデル式は

```math
\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の同時分布は

```math
\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)$$の対数尤度関数

```math
\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関数は以下の期待値で定義されます．

```math
\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関数を微分して，

```math
\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と置くことで

```math
\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]$という条件付き期待値は，

```math
\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}
```
として計算できます．条件付き期待値については，[こちらの記事(条件付き多変量正規分布)](https://qiita.com/m1t0/items/12b58f3583eab6e1e05c)をご覧ください．
　アルゴリズムは以下のような反復ステップで構成されています．

1. A,Dの初期値を発生させる．
2. $E[f|x_i],E[ff^T|x_i]$を計算する
3. (1),(2)のようにA,Dを更新する
4. 対数尤度を計算し，収束していれば終了．していなければ2に戻る．

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


#4. サンプルコード
　今回も論文に書いてある通り(=数式をそのままプログラムにしたandどこかを工夫して計算スピードをあげたりしてない)に実装しています．データはBig5と呼ばれる心理学のデータで，[ここから](http://bstat.jp/material/)ダウンロードしています．

```{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](https://www.csie.ntu.edu.tw/~mhyang/course/u0030/papers/Mixture%20of%20Factor%20Analyzers.pdf)
[EMアルゴリズム徹底解説](https://qiita.com/kenmatsu4/items/59ea3e5dfa3d4c161efb)