Generate Objective Functions and Others for Factor Analysis

Description

Generate Objective Functions and Others for Factor Analysis

Usage

generate_fa_funcs(S, print = 1)

Arguments

Value

This function generates the following functions as a list:

 critfa       因子分析の目的関数：パラメタは Lambda と psi or log(psi)
 critfapsi    因子分析の目的関数：パラメタは psi or log(psi) のみ
 gradfa_a     因子分析の目的関数の解析的一次微分：パラメタは Lambda と psi or log(psi)
 gradfa       因子分析の目的関数の数値的一次微分：パラメタは Lambda と psi or log(psi)
 HessQfa_a    因子分析の目的関数の数値的２次微分の近似：パラメタは Lambda と psi or log(psi)
 gradfapsi_a  因子分析の目的関数の解析的一次微分：パラメタは psi or log(psi) のみ
 gradfapsi    因子分析の目的関数の数値的一次微分：パラメタは psi or log(psi) のみ
 residfa      LS 因子分析の残差：パラメタは Lambda と psi or log(psi)
 Jacfa_a      LS 因子分析の解析的 Jacobian：パラメタは Lambda と psi or log(psi)
 residfapsi      LS 因子分析の残差：パラメタは psi or log(psi)
 Jacfapsi_a      LS 因子分析の解析的 Jacobian：パラメタは psi or log(psi)

The arguments to the functions:

 critfa       param, ndim, Phi = diag(ndim), logpsi = 0, method = "ML", LH = 0
 critfapsi    psi, ndim, Phi = diag(ndim), logpsi = 0, method = "ML", LH = 0
 gradfa_a     param, ndim, Phi = diag(ndim), logpsi = 0, method = "ML", LH = 0
 gradfa       param, ndim, Phi = diag(ndim), logpsi = 0, method = "ML", LH = 0
 HessQfa_a    param, ndim, Phi = diag(ndim), logpsi = 0, method = "LS", LH = 0
 gradfapsi_a  psiorlogpsi, ndim, Phi = diag(ndim), logpsi = 0, method = "ML", LH = 0
 gradfapsi    psiorlogpsi, ndim, Phi = diag(ndim), logpsi = 0, method = "ML", LH = 0
 residfa      param, ndim, Phi = diag(ndim), LH = 0, logpsi = 0, switch = 0, attrib = 0
 Jacfa_a      param, ndim, Phi = diag(ndim), logpsi = 0, LH = 0, switch = 0, attrib = 0
 residfapsi      param, ndim, Phi = diag(ndim), LH = 0, logpsi = 0, switch = 0, attrib = 0
 Jacfapsi_a      param, ndim, Phi = diag(ndim), logpsi = 0, LH = 0, switch = 0, attrib = 0

where

param        c( c(Lambda),psi ) or  c( c(Lambda),log(psi) )
psiorlogpsi  psi or log(psi)
ndim         # of factors
Phi          The factor correlation matrix
logpsi       = 1 if log(psi) is used as the parameter
method       = "ML" or "LS"
LH           = 1 to use lower triangular part of S to calculate LS criterion
switch       = 1 to define residual as 'model - data' for nlsr::nlfb
attrib       = 1 to return Jacobian as the attribute

Examples

#
# generate data
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)

# # of dimensions to be used
ndim <- 2

# 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


# S or vech(S)
LH <- 0


# initial values by pca
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))


# ML 因子分析： NR を利用：微分はお任せ
# ML fa by NR:  param=c( c(Lambda),psi )
resNR <- NR( param0, critfa, ndim=ndim, logpsi=0, maxiter=50, flipd=1 )
critNR <- resNR$objective
param <- resNR$par
Lambda <- matrix(param[1:(nvar*ndim)],nvar)
psi <- param[-(1:(nvar*ndim))]
Print(Lambda,psi)
gradNRa <- gradfa_a( param, ndim=ndim )
gradNRn <- gradfa( param, ndim=ndim )
Print(critNR,gradNRa, gradNRn)


# ML 因子分析： NR を利用：解析的１次微分を利用
# ML fa by NR with analytic gradient:  param=c( c(Lambda),psi )
resNR <- NR( param0, critfa, ndim=ndim, logpsi=0, maxiter=50, flipd=1
          , gradient=gradfa_a)
critNR <- resNR$objective
param <- resNR$par
Lambda <- matrix(param[1:(nvar*ndim)],nvar)
psi <- param[-(1:(nvar*ndim))]
Print(Lambda,psi)
gradNRa <- gradfa_a( param, ndim=ndim )
gradNRn <- gradfa( param, ndim=ndim )
Print(critNR,gradNRa, gradNRn)