Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
4
Help us understand the problem. What is going on with this article?
@m1t0

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

More than 1 year has passed since last update.

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アルゴリズム徹底解説

4
Help us understand the problem. What is going on with this article?
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
m1t0
Major: statistics/Multivariate analysis/Psychometrics/R/ (R) C++/python3

Comments

No comments
Sign up for free and join this conversation.
Sign Up
If you already have a Qiita account Login
4
Help us understand the problem. What is going on with this article?