Examples of ML Factor Anlysis

Description

Examples of Maximum Likelihood Factor Analysis using general optimization functions

Examples

## Not run: 

#
# データの生成
#
# 単純構造を持つ Lambda から Sigma を生成し
# それを母数として持つ Wishart 乱数で S を生成
#

seed <- 1701
set.seed(seed)

nvar <- 20;  ndim0 <- 3
ps <- 0.2
df <- 500

A <- matrix(runif(nvar*ndim0),nvar)
Sigma <- A%*%t(A)
Sigma <- Sigma+ps*diag(nvar)
dS <- sqrt(diag(Sigma))
Sigma <- diag(1/dS)%*%Sigma%*%diag(1/dS)

# correlation matrix from Wishart random matrix
S <- rWishart( 1, df, Sigma )
S <- S[,,1]/df
dS <- sqrt(diag(S))
S <- diag(1/dS)%*%S%*%diag(1/dS)


#
# 最適化のための関数 nlminb や lazy.mat::GN もしくは nlsr::nlfb で用いる
# 目的関数や一次微分を与える関数を生成する。
#

# generate objective functions and others
funcs <- generate_fa_funcs( S )
critfa <- funcs$critfa
critfapsi <- funcs$critfapsi
gradfa_a <- funcs$gradfa_a
gradfa <- funcs$gradfa
HessQfa_a <- funcs$HessQfa_a
gradfapsi_a <- funcs$gradfapsi_a
gradfapsi <- funcs$gradfapsi
residfa <- funcs$residfa
Jacfa_a <- funcs$Jacfa_a
residfapsi <- funcs$residfapsi
Jacfapsi_a <- funcs$Jacfapsi_a



# 初期値の設定
# # of dimensions to be used
ndim=2

# S or vech(S)
LH=0


# 主成分分析による初期値
ev=eigen( S )
Lambda0=ev$vectors[,1:ndim]%*%diag(sqrt(ev$values[1:ndim]), nrow=ndim)
psi0=pmax(0.001, diag(S-Lambda0%*%t(Lambda0)) )

# initial value
param0=c(c(Lambda0),psi0)
param00=c(c(Lambda0),log(psi0))

printm(param0,log(psi0))

#
# 以下で最尤法を用いて直交解を求める。
#


#
# 最尤解
#
#
# 直交解
#

#
# Lawley の方法: パラメタは Lambda と psi
#

res0=efa( S, ndim=ndim, method="ML", print=0 )
Lambda0=res0$Lambda
psi0=res0$psi
crit0=res0$crit
printm(crit0, Lambda0,psi0, fmt="12.8 10.5")
param=c(c(Lambda0),psi0)
grada=gradfa_a( param, ndim=ndim, method="ML", logpsi=0, LH=LH )
gradn=gradfa( param, ndim=ndim, method="ML", logpsi=0, LH=LH )
printm(crit0,grada,gradn, fmt="15.10")


#
# EM アルゴリズムによる解法
#

res00=efa( S, ndim=ndim, method="EM", print=0 )
Lambda00=res0$Lambda
psi00=res0$psi
crit00=res0$crit
printm(crit00, Lambda00,psi00, fmt="12.8 10.5")
param=c(c(Lambda00),psi00)
grada=gradfa_a( param, ndim=ndim, method="ML", logpsi=0, LH=LH )
gradn=gradfa( param, ndim=ndim, method="ML", logpsi=0, LH=LH )
printm(crit00,grada,gradn, fmt="15.10")



#
# Newton-Raphson 法による解法
#


# ML 因子分析： nlminb を利用：数値的一次微分 パラメタは Lambda, log(psi)
# ML fa by nlminb:  param=c(c(Lambda),log(psi))
res1=nlminb( param00, critfa, ndim=ndim
                  , control=list(trace=0), logpsi=1, method="ML", LH=LH )
crit1=res1$objective
param=res1$par
Lambda1=matrix(param[1:(nvar*ndim)],nvar)
Lambda1=rotate2canon(Lambda1)
psi1=exp(param[-(1:(nvar*ndim))])
printm(crit1, Lambda1,psi1, fmt="12.8 10.5")
grada=gradfa_a( param, ndim=ndim, method="ML", logpsi=1, LH=LH )
gradn=gradfa( param, ndim=ndim, method="ML", logpsi=1, LH=LH )
printm(crit1,grada,gradn, fmt="15.10")



# ML 因子分析： nlminb を利用：解析的一次微分 パラメタは Lambda, log(psi)
# ML fa by nlminb:  param=c(c(Lambda),log(psi))
res2=nlminb( param00, critfa, ndim=ndim, control=list(trace=0)
                  , logpsi=1, method="ML", LH=LH, gradient=gradfa_a )
crit2=res2$objective
param=res2$par
Lambda2=matrix(param[1:(nvar*ndim)],nvar)
Lambda2=rotate2canon(Lambda2)
psi2=exp(param[-(1:(nvar*ndim))])
printm(crit2, Lambda2,psi2, fmt="12.8 10.5")
grada=gradfa_a( param, ndim=ndim, method="ML", logpsi=1, LH=LH )
gradn=gradfa( param, ndim=ndim, method="ML", logpsi=1, LH=LH )
printm(crit2,grada,gradn, fmt="15.10")



# ML 因子分析： nlminb を利用：数値的一次微分  パラメタは log(psi) のみ
# ML fa by nlminb with reparametrization:  param=log(psi)
res3=nlminb( log(psi0), critfapsi, ndim=ndim, control=list(trace=0)
                  , logpsi=1, method="ML", LH=LH )
crit3=res3$objective
psi3=exp( res3$par )
Lambda3=recoverL( psi3, ndim=ndim, S=S, method="ML" )
printm(crit3,Lambda3,psi3, fmt="12.8 10.5")
grada=gradfapsi_a( res3$par, ndim=ndim, method="ML", logpsi=1 )
gradn=gradfapsi( res3$par, ndim=ndim, method="ML", logpsi=1, LH=LH )
printm(crit3,grada,gradn, fmt="15.10")



# ML 因子分析： nlminb を利用：解析的一次微分  パラメタは log(psi) のみ
# ML fa by nlminb with reparametrization:  param=log(psi)
res4=nlminb( log(psi0), critfapsi, ndim=ndim, control=list(trace=0)
                  , logpsi=1, method="ML", LH=LH, gradient=gradfapsi_a )
crit4=res4$objective
psi4=exp( res4$par )
Lambda4=recoverL( psi4, ndim=ndim, S=S, method="ML" )
printm(crit4,Lambda4,psi4, fmt="12.8 10.5")
grada=gradfapsi_a( res4$par, ndim=ndim, method="ML", logpsi=1 )
gradn=gradfapsi( res4$par, ndim=ndim, method="ML"
                 , logpsi=1, LH=LH )
printm(crit4,grada,gradn, fmt="15.10")


#
# 結果の比較：Lawley と解析的一次微分を使った解のみ
#
printm(crit0, crit2, crit4, fmt="15.10")
printm(psi0, psi2, psi4, fmt="8.5")
printm(Lambda0,Lambda2,Lambda4, fmt="8.5")




#
# 直交解を変換して斜交解を求める。
#

# factor correlation matrix (constant)
Phi=diag(ndim)
Phi[1,2]=Phi[2,1]=0.3

Lambda=orthog2obliq( Lambda0, Phi )
printm(Lambda0, Lambda, Phi, fmt="10.5")
param=c(c(Lambda),psi0)
crit=critfa(param, ndim=ndim, Phi=Phi, method="ML", LH=LH )
grada=gradfa_a( param, ndim=ndim, method="ML", logpsi=0, Phi=Phi )
gradn=gradfa( param, ndim=ndim, method="ML", logpsi=0, Phi=Phi, LH=LH )
printm(crit,grada,gradn, fmt="15.10")








## End(Not run)