Examples of LS Factor Anlysis

Description

Examples of Least Squares 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))

#
# 以下では最小２乗法を用いて直交解を求める。
#


#
# 最小２乗解
#

#
# ALS を用いる方法: パラメタは Lambda と psi
#
res0=efa( S, ndim=ndim, method="LS", 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="LS", logpsi=0, LH=LH )
gradn=gradfa( param, ndim=ndim, method="LS", logpsi=0, LH=LH )
printm(crit0,grada,gradn, fmt="15.10")


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

#
# パラメタは Lambda と log(psi)
#

# LS 因子分析： nlminb を利用：数値的一次微分 パラメタは Lambda, log(psi)
# LS fa by nlminb:  param=c(c(Lambda),log(psi))
res1=nlminb( param00, critfa, ndim=ndim
                  , control=list(trace=0), logpsi=1, method="LS", 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="LS", logpsi=1, LH=LH )
gradn=gradfa( param, ndim=ndim, method="LS", logpsi=1, LH=LH )
printm(crit1,grada,gradn, fmt="15.10")



# LS 因子分析： nlminb を利用：解析的一次微分 パラメタは Lambda, log(psi)
# LS fa by nlminb:  param=c(c(Lambda),log(psi))
res2=nlminb( param00, critfa, ndim=ndim, control=list(trace=0)
                  , logpsi=1, method="LS", 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="LS", logpsi=1, LH=LH )
gradn=gradfa( param, ndim=ndim, method="LS", logpsi=1, LH=LH )
printm(crit2,grada,gradn, fmt="15.10")



#
# パラメタは log(psi) のみ
#


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



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


#
# 結果の比較：als と解析的一次微分を使った解のみ
#
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="LS", LH=LH )
grada=gradfa_a( param, ndim=ndim, method="LS", logpsi=0, Phi=Phi )
gradn=gradfa( param, ndim=ndim, method="LS", logpsi=0, Phi=Phi, LH=LH )
printm(crit,grada,gradn, fmt="15.10")






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

#
# Lambda と log(psi) をパラメタとして用いる場合
#

if(0){
 # nlsr パッケージが必要
 library(nlsr)
 # LS 因子分析：nlsr::nlfb を利用：数値的ヤコビアン パラメタは Lambda, log(psi)
 # LS fa by nlsr::nlfb  param=c(c(Lambda),log(psi))
 res5=nlfb( param00, residfa, ndim=ndim, LH=LH, switch=1
            , control=list(femax=1000,japprox="jacentral"), logpsi=1 )
 param=res5$coefficients
 crit5=res5$ssquares
 Lambda5=matrix(param[1:(nvar*ndim)],nvar)
 psi5=exp( param[-(1:(nvar*ndim))] )
 printm(crit5,Lambda5,psi5, fmt="12.8 10.5")
 grada=gradfa_a( param, ndim=ndim, LH=LH, method="LS", logpsi=1 )
 gradn=gradfa( param, ndim=ndim, Phi=Phi, LH=LH, method="LS", logpsi=1 )
 printm(crit5,grada,gradn, fmt="15.10")


 # LS 因子分析：nlsr::nlfb を利用：解析的ヤコビアン パラメタは Lambda, log(psi)
 # LS fa by nlsr::nlfb with analytic jacobian:  param=c(c(Lambda),log(psi))
 res6=nlfb( param00, residfa, ndim=ndim, LH=LH, switch=1
            , control=list(femax=1000), logpsi=1, jacfn=Jacfa_a, attrib=1 )
 param=res6$coefficients
 crit6=res6$ssquares
 Lambda6=matrix(param[1:(nvar*ndim)],nvar)
 psi6=exp( param[-(1:(nvar*ndim))] )
 printm(crit6,Lambda6,psi6, fmt="12.8 10.5")
 grada=gradfa_a( param, ndim=ndim, LH=LH, method="LS", logpsi=1 )
 gradn=gradfa( param, ndim=ndim, LH=LH, method="LS", logpsi=1 )
 printm(crit6,grada,gradn, fmt="15.10")
} # if(0)



# LS 因子分析：GN を利用：数値的ヤコビアン パラメタは Lambda, log(psi)
# LS fa by GN: param=c(c(Lambda),log(psi))
res5=GN( param00, residfa, ndim=ndim, LH=LH, logpsi=1, maxiter=50 )
param=res5$par
crit5=res5$objective
Lambda5=matrix(param[1:(nvar*ndim)],nvar)
psi5=exp( param[-(1:(nvar*ndim))] )
printm(crit5,Lambda5,psi5, fmt="12.8 10.5")
grada=gradfa_a( param, ndim=ndim, LH=LH, method="LS", logpsi=1 )
gradn=gradfa( param, ndim=ndim, LH=LH, method="LS", logpsi=1 )
printm(crit5,grada,gradn, fmt="15.10")


# LS 因子分析：GN を利用：　解析的ヤコビアン パラメタは Lambda, log(psi)
# LS fa by GN with analytic jacobian: param=c(c(Lambda),log(psi))
res6=GN( param00, residfa, ndim=ndim, LH=LH, logpsi=1, maxiter=50
          , jacobian=Jacfa_a )
param=res6$par
crit6=res6$objective
Lambda6=matrix(param[1:(nvar*ndim)],nvar)
psi6=exp( param[-(1:(nvar*ndim))] )
printm(crit6,Lambda5,psi5, fmt="12.8 10.5")
grada=gradfa_a( param, ndim=ndim, LH=LH, method="LS", logpsi=1 )
gradn=gradfa( param, ndim=ndim, LH=LH, method="LS", logpsi=1 )
printm(crit6,grada,gradn, fmt="15.10")




#
# log(psi) のみをパラメタとして用いる場合
#


if(0){
 # nlsr パッケージが必要
 library(nlsr)
 # LS fa by nlsr::nlfb  param=c(c(Lambda),log(psi))
 res7=nlfb( log(psi0), residfapsi, ndim=ndim, LH=LH, switch=1
               , control=list(femax=1000,japprox="jacentral"), logpsi=1 )
 param=res7$coefficients
 crit7=res7$ssquares
 psi=exp( param )
 Lambda=recoverL( psi, ndim=ndim, S=S, method="LS" )
 printm(Lambda,psi)
 grada=gradfapsi_a( param, ndim=ndim, LH=LH, method="LS", logpsi=1 )
 gradn=gradfapsi( param, ndim=ndim, LH=LH, method="LS", logpsi=1 )
 printm(crit7, grada, gradn, fmt="15.10")


 # LS fa by nlsr::nlfb  param=c(c(Lambda),log(psi))
 res8=nlfb( log(psi0), residfapsi, ndim=ndim, LH=LH, switch=1, attrib=1
               , control=list(femax=1000,japprox="jacentral")
               , logpsi=1, jacfn=Jacfapsi_a  )
 param=res8$coefficients
 crit8=res8$ssquares
 psi=exp( param )
 Lambda=recoverL( psi, ndim=ndim, S=S, method="LS" )
 printm(Lambda,psi)
 grada=gradfapsi_a( param, ndim=ndim, LH=LH, method="LS", logpsi=1 )
 gradn=gradfapsi( param, ndim=ndim, LH=LH, method="LS", logpsi=1 )
 printm(crit8, grada, gradn, fmt="15.10")
} # if(0)


# LS 因子分析： lazy.mat::GN を利用：数値的ヤコビアン パラメタは log(psi) のみ
# LS fa by GN: param=log(psi)
res7=GN( log(psi0), residfapsi, ndim=ndim, LH=LH, logpsi=1, maxiter=50 )
param=res7$par
crit7=res7$objective
psi7=exp( param )
Lambda7=recoverL( psi7, ndim=ndim, S=S, method="LS" )
printm(crit7,Lambda7,psi7, fmt="12.8 10.5")
grada=gradfapsi_a( param, ndim=ndim, LH=LH, method="LS", logpsi=1 )
gradn=gradfapsi( param, ndim=ndim, LH=LH, method="LS", logpsi=1 )
printm(crit7,grada,gradn, fmt="15.10")


# LS 因子分析： lazy.mat::GN を利用：解析的ヤコビアン パラメタは log(psi) のみ
# LS fa by GN: param=log(psi)
res8=GN( log(psi0), residfapsi, ndim=ndim, LH=LH, logpsi=1, maxiter=50
         , jacobian=Jacfapsi_a )
param=res8$par
crit8=res8$objective
psi8=exp( param )
Lambda8=recoverL( psi8, ndim=ndim, S=S, method="LS" )
printm(crit8,Lambda8,psi8, fmt="12.8 10.5")
grada=gradfapsi_a( param, ndim=ndim, LH=LH, method="LS", logpsi=1 )
gradn=gradfapsi( param, ndim=ndim, LH=LH, method="LS", logpsi=1 )
printm(crit8,grada,gradn, fmt="15.10")



#
# 結果の比較：als と解析的一次微分を使った解のみ
#
printm(crit0, crit7, crit7, fmt="15.10")
printm(psi0, psi7, psi8, fmt="8.5")
printm(Lambda0,Lambda7,Lambda8, fmt="8.5")


## End(Not run)