| ExamplesOfGNNR {lazy.mat} | R Documentation |
Examples of GN and NR in the context of Logistic Regression. Multidimensional Scaling, Factor Analysis and PCA.
## Not run:
#
# For comparison reason, the following package is required:
library(nlsr)
#
#####################################################################
#
# ML and Bayes MAP Estimation of Logistic Regression Parameters
# by Newton-Raphson Method
#
#####################################################################
#
# probability (logistic function with matrix X, no 1.7)
# X is the regressor matrix, beta is the regression coefficient vector
# Returns a vector of length nrow(X).
Prob <- function( beta, X ){
if( is.matrix(X) | is.data.frame(X) ) z=X%*%beta
else z=X*beta
P=1/(1+exp(-z))
return( P )
} # end of P
#
#
#
# This is the objective function to be minimized:
# Minus log likelihood
# parameter beta is the first argument and the second argument data
# will be passed to this via ... parameter.
# Returns a scalar.
#
mllh <- function( beta, data ){
# minus log likelihood: multiple regressor matrix X
locx=which(regexpr("x",colnames(data)) > 0)
X=as.matrix(data[,locx]); f=data$f; r=data$r
p=Prob(beta,X)
llh=sum( r*log(p)+(f-r)*log(1-p) )
return( -llh )
} # end of mllh
#
#
# First derivative of mllh wrt beta.
# parameter beta is the first argument and the second argument data
# will be passed to this via ... parameter.
# Returns a vector of size length(beta).
#
dmllh_a <- function( beta, data ){
# analytic first derivative of mllh: multiple regressor matrix X
locx=which(regexpr("x",colnames(data)) > 0)
X=as.matrix(data[,locx]); f=data$f; r=data$r
p=Prob(beta,X)
dab=colSums(c(r-f*p)*X)
# cat(".da.")
return( -dab )
} # end of dmllh_a
#
#
dmllh_n <- function( beta, data, mllh=NULL ){
# numeric first derivative of mllh: multiple regressor matrix X
# Not Used.
res=JacobianMat( beta, mllh, data=data)
return( res )
} # end of dmllh_n
#
#
#
# Second derivative of mllh wrt beta.
# parameter beta is the first argument and the second argument data
# will be passed to this via ... parameter.
# Returns a square matrix of size length(beta).
#
d2mllh_nn <- function( param, data, mllh=NULL ){
# numeric second derivative matrix using numeric first derivative of mllh
# dmllh_n comes from the parent environment
# Not Used.
res=JacobianMat( param, dmllh_n, data=data, mllh=mllh )
# cat(".d2an.")
return( res )
} # end of d2mllh_nn
#
#
d2mllh_an <- function( param, data ){
# numeric second derivative matrix using analytic first derivative of mllh
# dmllh_a comes from the parent environment
# Not Used.
res=JacobianMat( param, dmllh_a, data=data )
# cat(".d2an.")
return( res )
} # end of d2mllh_an
#
#
d2mllh_aa <- function( param, data ){
# analytic second derivative matrix of mllh
# not yet available
locx=which(regexpr("x",colnames(data)) > 0)
X=as.matrix(data[,locx]); f=data$f; r=data$r
p=Prob(param,X)
dPdbetaa=c(p*(1-p))*X
H=t(f*X)%*%dPdbetaa
return(H)
} # end of d2mllh_aa
#
#
#
#
#
#
# generate data
#
#
# seed for random numbers
seed=1701
set.seed(seed)
#
# true value of beta
beta=c(-1,2,3)
# # of parameters = # of regressors = # of columns of X including intercept.
nq=length(beta)
#
# # of trials per obs
ff=c(1,2,3) # will be recycled.
# ff=1
#
# generate regressor matrix X
n=50
#
# regressor matrix as 1 and normal random variable
X=cbind(1, matrix(rnorm(n*(nq-1)),n,nq-1) )
colnames(X)=paste("x",(1:nq)-1,sep="")
#
# true response probability
p=Prob(beta,X)
#
# generate f=# of trials, r=# of successes, y=sample proportion
f=sample( ff, n, rep=1 )
r=mapply( function(size,prob) rbinom(1,size,prob) , f, p )
y=r/f
#
# data frame
data=data.frame( id=1:n, X, f, r, y, p )
print(data)
#
# initial value
# beta0=c(0,1,1)
beta0=rev(beta)
#
#
#
# Estimation Starts Here.
#
#
#
# Arguments to nlminb or NR
#
# First "beta0" argument is the starting value of the parameter.
# Second argument "mllh" is the objective function to be minimized
# which must have "beta0" as its first argument.
# The additional argument to this function can be passed through ... .
# or, the objects in the environment in which the objective function
# "mllh" is defined can be accessed from within the objective function.
# In the example below, data, which contains f and r,
# is passed to "mllh" via ... argument of nlminb or NR.
#
# Note that, unlike nlm or nlfb which require the gradient to be returned
# as an attribute "gradient" of the function value,
# nlminb or NR uses the gradient or hessian functions which return
# the gradient or hessian as their values.
#
#
# mle by native nlminb
# no analytic gradient and hessian
# Ones like dmllh_n or d2mllh_nn must be used in nlminb.
resnlminb=nlminb( beta0, mllh, data=data, control=list(trace=0) )
# with analytic gradient
resnlminb2=nlminb( beta0, mllh, data=data, control=list(trace=0)
, gradient=dmllh_a )
# with analytic gradient and Hessian
resnlminb3=nlminb( beta0, mllh, data=data, control=list(trace=0)
, gradient=dmllh_a, hessian=d2mllh_aa )
#
#
# solution by NR
# no analytic gradient and hessian
# Functions similar to dmllh_n or d2mllh_nn will be created by in NR.
resNR=NR( x=beta0, mllh, data=data, print=2 )
# with analytic gradient
resNR2=NR( x=beta0, mllh, data=data, gradient=dmllh_a, print=2 )
# with analytic gradient and Hessian
resNR3=NR( x=beta0, mllh, data=data, gradient=dmllh_a, hessian=d2mllh_aa
, print=2 )
#
#
# comparison of results
p1=resnlminb$par; p2=resnlminb2$par; p3=resnlminb3$par
pNR1=resNR$par; pNR2=resNR2$par; pNR3=resNR3$par
cat("\n\nComparizon of the Results: Parameters\n")
Print(p1,p2,p3)
Print(pNR1,pNR2,pNR3)
#
# comparison of gradient
g1=dmllh_a( p1, data ); g2=dmllh_a( p2, data ); g3=dmllh_a( p3, data )
gNR1=dmllh_a( pNR1, data ); gNR2=dmllh_a( pNR2, data )
gNR3=dmllh_a( pNR3, data )
cat("\n\nComparison of the Results: Gradients\n")
Print(g1,g2,g3)
Print(gNR1,gNR2,gNR3)
#
#
#
#
#
# Bayes MAP
#
#
# prior
#
# The prior distribution of beta is N(mu, inv(invSigma)), where
# mu is the prior mean vector and invSigma is the prior dispersion matrix.
#
# In addition to data which contains (r,f), mu and invSigma will be
# passed to the objective function "mlogpostpdf" via ... argument
# of NR or nlminb.
#
#
logpriorpdf <- function( beta, mu, invSigma ){
# logarithm of multivariate normal pdf (constant terms omitted)
lpdf=-0.5*c(beta-mu)%*%invSigma%*%c(beta-mu) + 0.5*det(invSigma)
return(lpdf)
} # logpriorpdf
#
# popsterior pdf is proportional to the likelihood x prior pdf
mlogpostpdf <- function( beta, data, mu, invSigma ){
# minus log posterior density
# log post pdf of beta is log likelihood + log prior pdf
lppdf= -mllh( beta, data ) + logpriorpdf( beta, mu, invSigma )
return( -lppdf )
} # end of mlogpostpdf
#
#
#
# prior constant (hyper parameters) vague prior N(0,100 I)
mu=rep(0,length(beta0)); invSigma=0.01*diag(length(beta0))
#
# MAP estimates
resNR4=NR( x=beta0, mlogpostpdf, data=data, mu=mu, invSigma=invSigma )
pNR4=resNR4$par;
resnlminb4=nlminb( beta0, mlogpostpdf, data=data, mu=mu, invSigma=invSigma
, control=list(trace=0) )
p4=resnlminb4$par
#
# evaluate the final gradient numerically
gNR4=JacobianMat( pNR4, mlogpostpdf, data=data, mu=mu, invSigma=invSigma )
g4=JacobianMat( p4, mlogpostpdf, data=data, mu=mu, invSigma=invSigma )
#
#
# prior constant (hyper parameters) sharp prior N(0,I)
mu=rep(0,length(beta0)); invSigma=diag(length(beta0))
#
# MAP estimates
resNR5=NR( x=beta0, mlogpostpdf, data=data, mu=mu, invSigma=invSigma )
pNR5=resNR5$par
resnlminb5=nlminb( beta0, mlogpostpdf, data=data, mu=mu, invSigma=invSigma
, control=list(trace=0) )
p5=resnlminb5$par
#
# evaluate the final gradient numerically
gNR5=JacobianMat( pNR5, mlogpostpdf, data=data, mu=mu, invSigma=invSigma )
g5=JacobianMat( p5, mlogpostpdf, data=data, mu=mu, invSigma=invSigma )
#
cat("\n\nComparison\n")
Print(p1,p4,p5)
Print(pNR1,pNR4,pNR5)
Print(gNR4,g4,gNR5,g5)
#
#
#####################################################################
#
# End of
# ML and Bayes MAP Estimation of Logistic Regression Parameters
# by Newton-Raphson Method
#
#####################################################################
#
#
#
#
#
#
#
#
#
#####################################################################
#
# Multidimensional Scaling
#
#####################################################################
#
#
residmds <- function( vecX ){
# residual of metric mds: vech(O-D)
# Shin-ichi Mayeakawa
# 20210713,20220713,14,22
#
# Note that ns, O, and locLH come from the environment where
# this function was defined.
#
# convert vecX to matrix X
X=matrix(vecX,ns)
#
# Euclidean distance
XX=rowSums(X*X)
XX11=XX%*%matrix(1,1,ns)
D=sqrt( abs( XX11-2*X%*%t(X)+t(XX11) ) )
#
# lower half of O-D
resid=c( (O-D)[locLH] )
#
return( resid )
#
} # end of residmds
#
#
residmds2 <- function( vecX ){
# residual of metric mds: vech(D-O) to be used by nlfb
# Shin-ichi Mayeakawa
# 20210713,20220713,14,22
#
return( -residmds( vecX ) )
} # end of residmds2
#
#
rssmds <- function( vecX ){
# residual of metric mds: ssq(vech(O-D))
# Shin-ichi Mayeakawa
# 20210713,20220713,14,22
#
return( sum( residmds(vecX)^2 ) )
} # end of rssmds
#
#
Jacmds <- function( vecX ){
# Numeric Jacobian of vech(O-D), i.e. -vech(D)
# Shin-ichi Mayeakawa
# 20210713,20220713,14,22
F=JacobianMat( vecX, residmds )
return( F )
} # end of gradmds
#
#
Jacmds_a <- function( vecX ){
# Analytic Jacobian of vech(D)
# Shin-ichi Mayeakawa
# 20220723
# convert vecX to matrix X
X=matrix(vecX,ns)
# Euclidean distance
XX=rowSums(X*X)
XX11=XX%*%matrix(1,1,ns)
D=sqrt( abs( XX11-2*X%*%t(X)+t(XX11) ) )
vecD=D[locLH]
#
# This is the Jacobian of vech(D) wrt vec(X)
F=matrix(0,npair,ns*ndim)
#
if(0){
k=0
for( i in 2:ns ){
for( j in 1:(i-1) ){
k=k+1
for( a in 1:ndim ){
F[k,(a-1)*ns+i]=(X[i,a]-X[j,a])/D[i,j]
F[k,(a-1)*ns+j]=-(X[i,a]-X[j,a])/D[i,j]
}
}
}
} # en of skip
#
GX=G%*%X
for( k in 1:npair ){
for( a in 1:ndim ){
ximxja=GX[k,a]
F[k,(a-1)*ns+locvLH[k,1]]=ximxja/vecD[k]
F[k,(a-1)*ns+locvLH[k,2]]=-ximxja/vecD[k]
}
}
#
return( F )
#
} # end of Jacmds_a
#
#
HessQmds_a <- function( vecX ){
# approximate Hessian calculated from analytic Jacobian of vech(D)
# Shin-ichi Mayeakawa
# 20220722,23
F=Jacmds_a( vecX )
H=t(F)%*%F
H=H + 0.01*min(abs(diag(H)))*diag(nrow(H))
return( 2*H )
} # end of HessQmds_a
#
#
gradmds <- function( vecX ){
# numeric gradient of rssmds
# Shin-ichi Mayeakawa
# 20210713,20220713,14,22,23
res=JacobianMat( vecX, rssmds )
return( res )
} # end of gradmds
#
#
gradmds_a <- function( vecX ){
# analytic gradient of rssmds
# Shin-ichi Mayeakawa
# 20220723
res=c( -2*residmds(vecX)%*%Jacmds_a(vecX) )
return( res )
} # end of gradmds_a
#
#
normalizemds <- function( vecX ){
# normalize X in each of GN/NR iteration
# Shin-ichi Mayekawa
# 20220723
#
# convert vecX to matrix X
X=matrix(vecX,ns)
# column center X
return( c( normalize_config(X,center=1) ))
} # end of normalizemds
#
#
#
#
# generate test data
#
# X=as.matrix( expand.grid(-3:3, -3:3) )
X=as.matrix( expand.grid(-2:2, -2:2) )
# X=as.matrix( expand.grid(-1:1,-1:1) )
ns=nrow(X); ndim=ncol(X)
sname=paste("s",1:ns,sep="")
#
set.seed(1701)
sigmaE=0.5
XX=rowSums(X*X)
XX11=XX%*%matrix(1,1,ns)
D=sqrt( abs( XX11-2*X%*%t(X)+t(XX11) ) )
O=D+sigmaE*rnorm(ns*ns)
O=round( 0.5*(O+t(O))-diag(diag(O)), 2 )
locLH=lower.tri(O, diag=FALSE)
npair=sum(locLH)
# locvLH=vechindex(ns,1,0,type=2)
locvLH=cbind(row(O)[locLH],col(O)[locLH])
# G matrix for paired comparison
G=matrix(0,npair,ns)
for( k in 1:npair ){
G[k,locvLH[k,1]]=1; G[k,locvLH[k,2]]=-1
}
#
#
#
# # of dimensions
ndim=2
#
#
# initial by classic mds
J=diag(ns)-matrix(1,ns,ns)/ns
B=-0.5*J%*%O^2%*%J
ev=eigen(B)
P=ev$vectors[,1:ndim]
L=diag(ev$values[1:ndim])
Xcmds=P%*%sqrt(L)
#
X0=Xcmds
rss0=rssmds( c(X0) )
g0=gradmds( c(X0) )
title=paste("CMDS: rss =", round(rss0,3))
plot(X0, main=title)
text( X0, sname, pos=4 )
#
#
#
# Gauss-Newton
#
# by nlfb
resnlfb=nlfb(c(X0),residmds2, control=list(femax=1000,japprox="jacentral"))
X1=matrix(resnlfb$coefficients,ns)
rss1=rssmds( c(X1) )
g1=gradmds( c(X1) )
title=paste("GN by nlfb: rss =", round(rss1,3))
plot(X1, main=title)
text( X1, sname, pos=4 )
#
#
# by GN
resGN=GN( c(X0), residmds, maxiter=100, normalize=normalizemds )
X2=matrix(resGN$par,ns)
rss2=rssmds( c(X2) )
g2=gradmds( c(X2) )
title=paste("GN by GN: rss =", round(rss2,3))
plot(X2, main=title)
text( X2, sname, pos=4 )
#
#
# by GN with analytic gradient
resGN2=GN( c(X0), residmds, maxiter=100, normalize=normalizemds
, jacobian=Jacmds_a )
X21=matrix(resGN2$par,ns)
rss21=rssmds( c(X21) )
g21=gradmds( c(X21) )
title=paste("GN by GN with analytic gradient: rss =", round(rss2,3))
plot(X21, main=title)
text( X21, sname, pos=4 )
#
Print(rss0,rss1,rss2, rss21)
#
#
#
#
# Newton-Raphson
#
# by nlminb
resnlminb=nlminb( c(X0), rssmds, control=list(trace=0) )
X3=matrix(resnlminb$par,ns)
rss3=rssmds( c(X3) )
g3=gradmds( c(X3) )
title=paste("NR by nlminb: rss =", round(rss3,3))
plot(X3, main=title)
text( X3, sname, pos=4 )
#
#
# by NR
resNR=NR( c(X0), rssmds, maxiter=1000, flipd=1, normalize=normalizemds )
X4=matrix(resNR$par,ns)
rss4=rssmds( c(X4) )
g4=gradmds( c(X4) )
title=paste("NR by NR: rss =", round(rss4,3))
plot(X4, main=title)
text( X4, sname, pos=4 )
#
Print(rss0, rss1, rss2, rss3, rss4)
#
#
# by NR by supplying analytic gradient and quasi Hessian.
resNR2=NR( c(X0), rssmds, maxiter=1000, flipd=1
, normalize=normalizemds, grad=gradmds_a, hessian=HessQmds_a )
#
#
X41=matrix(resNR2$par,ns)
rss41=rssmds( c(X41) )
g41=gradmds( c(X41) )
title=paste("NR by NR with quasi H: rss =", round(rss41,3))
plot(X41, main=title)
text( X41, sname, pos=4 )
#
Print(rss1, rss2, rss3, rss4, rss41)
#
#
#
#
if(0){
#
library(lazy.mds)
#
# wmdsGN
resmds=wmdsGN( O, ndim=2, mlevel=9, maxiter=1000, print=1 )
Xw=resmds$X
rssw=rssmds( c(Xw) )
gw=gradmds( c(Xw) )
title=paste("GV by wmdsGN: rss =", round(rssw,3))
plot(Xw, main=title)
text( Xw, sname, pos=4 )
#
Print(rss0, rss1, rss2, rss3, rss4, rssw)
#
Print(max(abs(g0)), max(abs(g1)),max(abs(g2))
,max(abs(g3)),max(abs(g4)),max(abs(gw)))
#
#
}
#
#
#
#####################################################################
#
# End of Multidimensional Scaling
#
#####################################################################
#
#
#
#
#
#
#
#
#
#####################################################################
#
# Factor Analysis
#
#####################################################################
#
#
#
critfa <- function( param, logpsi=0, method="ML", LH=0 ){
# evaluation of ML or LS criterion function: param=c( c(Lambda), psi )
# Shin-ichi Mayekawa
# 20210705,20220723
# uselogpsiasparam -> logpsi: 20220728
#
# S, logdets, ndim and nvar are from the parent env.
#
#
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=param[-(1:(nvar*ndim))]
if( logpsi ) psi=exp(psi)
Sigma=Lambda%*%Phi%*%t(Lambda)+diag(psi)
#
if( method == "ML" ){
invSigma=solve(Sigma)
crit=tr(S%*%invSigma)+log(det(Sigma)) - logdetS - nvar
}
else{
if( LH == 1 ) crit=sum((S-Sigma)[locLH]^2)
else crit=sum((S-Sigma)^2)
}
#
return( crit )
#
} # end of critfa
#
#
#
#
critfapsi <- function( psi, logpsi=0, method="ML", LH=0 ){
# evaluation of ML or LS criterion function: param=psi
# Shin-ichi Mayekawa
# 20210707
# uselogpsiasparam -> logpsi: 20220728
#
# S, logdets, ndim and nvar are from the parent env.
#
#
# psi or log(psi)
# If the parameter is log(psi), recover original as exp(log(psi)).
if( logpsi ) psi=exp(psi)
#
psi=c(psi)
#
if( method == "ML" ){
# conditional MLE of Lambda given psi
Psimh=diag(1/sqrt(psi))
ev=eigen( Psimh%*%S%*%Psimh )
P=ev$vectors[,1:ndim]
D=ev$values[1:ndim]
Lambda=sqrt(psi)*P%*%diag(sqrt(D-1), nrow=ndim)%*%invPhimh
# log likelihood ratio stat
Sigma=Lambda%*%t(Lambda)+diag(psi)
invSigma=solve(Sigma)
crit=tr(S%*%invSigma)+log(det(Sigma)) - logdetS - nvar
}
else{
# conditional LSE of Lambda given psi
ev=eigen( S-diag(psi) )
P=ev$vectors[,1:ndim]
Gamma=ev$values[1:ndim]
Lambda=P%*%diag(sqrt(Gamma))
# rss
Sigma=Lambda%*%t(Lambda)+diag(psi)
if( LH == 1 ) crit=sum((S-Sigma)[locLH]^2)
else crit=sum((S-Sigma)^2)
}
#
return( crit )
#
} # end of critfapsi
#
#
#
#
gradfa_a <- function( param, logpsi=0, method="ML", LH=0 ){
# analytic gradient of LS or ML critfa wrt param=c(c(Lambda,psi)
# Shin-ichi Mayekawa
# 20220724,25
# avoid Jacfa_a: 20220726dnc
# correction for rss_LH: 20220726dnc,0728
# uselogpsiasparam -> logpsi: 20220728
#
#
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=param[-(1:(nvar*ndim))]
if( logpsi ) psi=exp(psi)
Sigma=Lambda%*%Phi%*%t(Lambda)+diag(psi)
#
if( method == "LS" ){
SS=S-Sigma
# dLambda=0.5*(-4*diag(diag(SS))%*%Lambda -4*SS%*%Lambda )
if( LH == 0 ) dLambda=-4*SS%*%Lambda%*%Phi
else dLambda=-2*((diag(diag(SS))+SS)%*%Lambda%*%Phi )
dpsi=-2*diag(SS)
if( logpsi == 1 ) dpsi=dpsi*psi
g=c( c(dLambda), dpsi )
}
else{
invSigma=matSwp(Sigma)
iSSmSiS=invSigma%*%(Sigma-S)%*%invSigma
dLambda=2*iSSmSiS%*%Lambda%*%Phi
dpsi=diag(iSSmSiS)
if( logpsi == 1 ) dpsi=dpsi*psi
g=c(c(dLambda),dpsi)
}
#
return( g )
#
} # gradfa_a
#
#
gradfa <- function( param, logpsi=0, method="ML" ){
# numerical gradient of LS or ML critfa wrt param=c(c(Lambda,psi)
# Shin-ichi Mayekawa
# 20220724,25
# uselogpsiasparam -> logpsi: 20220728
#
return( JacobianMat( param, critfa, method=method
, logpsi=logpsi ) )
} # end of gradfa
#
#
#
HessQfa_a <- function( param, logpsi=0, method="LS", LH=0 ){
# analytic quasi Hessian of LS or ML critfa wrt param=c(c(Lambda,psi)
# Shin-ichi Mayekawa
# 20220725
# LH: 20220728
# uselogpsiasparam -> logpsi: 20220728
#
if( method == "ML" ){
cat("\nerror(HessQfa) This is for method=LS only.\n")
stop()
}
F=Jacfa_a( param, logpsi=logpsi, LH=LH )
H=t(F)%*%F
H=H + 0.01*min(abs(diag(H)))*diag(nrow(H))
return( 2*H )
} # end of HessQfa_a
#
#
#
gradfapsi <- function( psiorlogpsi, logpsi=0, method="ML" ){
# numerical gradient of LS or ML critfa wrt param=c(c(Lambda,psi)
# Shin-ichi Mayekawa
# 20220726
# uselogpsiasparam -> logpsi: 20220728
#
# This works with LH=0 and LH=1
#
return( JacobianMat( psiorlogpsi, critfapsi, method=method
, logpsi=logpsi ) )
} # end of gradfapsi
#
#
gradfapsi_a <- function( psiorlogpsi, logpsi=0, method="ML" ){
# analytic gradient of LS or ML critfapsi wrt param=psi or log(psi)
# Shin-ichi Mayekawa
# 20220726
# uselogpsiasparam -> logpsi: 20220728
#
# This is for LH=0
#
#
# psi or log(psi)
# If the parameter is log(psi), recover original as exp(log(psi)).
if( logpsi ) psi=exp(psiorlogpsi)
else psi=psiorlogpsi
psi=c(psi)
#
if( method == "ML" ){
# conditional MLE of Lambda given psi
Psimh=diag(1/sqrt(psi))
ev=eigen( Psimh%*%S%*%Psimh )
Omega2=ev$vectors[,(ndim+1):nvar]^2
theta=ev$values[(ndim+1):nvar]
o2t=Omega2%*%diag(1-theta)
g=c( rowSums(o2t)/psi )
if( logpsi == 1 ) g=g*psi
}
else{
# conditional LSE of Lambda given psi
ev=eigen( S-diag(psi) )
P2=ev$vectors[,(ndim+1):nvar]^2
gamma=ev$values[(ndim+1):nvar]
g=c(-2*rowSums(P2%*%diag(gamma)) )
# above is correct if rss=ssq(S-Sigma), not if rss=ssq((S-Sigma)[locLH])
if( logpsi == 1 ) g=g*psi
}
#
return(g)
#
} # end of gradfapsi_a
#
#
#
#
#
#
residfa <- function( param, logpsi=0 ){
# residual of LS fa vec(S-Sigma)
# Shin-ichi Mayekawa
# 20220728
# uselogpsiasparam -> logpsi: 20220728
#
#
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=param[-(1:(nvar*ndim))]
if( logpsi ) psi=exp(psi)
Sigma=Lambda%*%Phi%*%t(Lambda)+diag(psi)
#
return( c(S-Sigma) )
#
} # end of residfa
#
#
residfa2 <- function( param, logpsi=0 ){
# residual of LS fa vec(Sigma-S)
# Shin-ichi Mayekawa
# 20220728
# uselogpsiasparam -> logpsi: 20220728
#
#
return( -residfa( param, logpsi=logpsi ) )
#
} # end of residfa2
#
#
#
Jacfa_a <- function( param, logpsi=0, LH=0 ){
# analytic Jacobian of vech(Sigma) or vec(Sigma) wrt param=c(c(Lambda,psi)
# Shin-ichi Mayekawa
# 20220724,25
# LH=0: 20220728
# uselogpsiasparam -> logpsi: 20220728
#
#
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=param[-(1:(nvar*ndim))]
if( logpsi ) psi=exp(psi)
Sigma=Lambda%*%Phi%*%t(Lambda)+diag(psi)
#
#
if( LH == 1 ){
npair=nvar*(nvar-1)/2
# First, calculate d vec(Sigma) / d vec(Lambda')
F=matrix(0,npair,nvar*ndim+nvar)
#
for( k in 1:npair ){
i=locvLH[k,1]; j=locvLH[k,2]
if( i == j ){
F[k,ndim*(i-1)+(1:ndim)]=2*Lambda[i,]
if( logpsi == 0 ) F[k,ndim*nvar+i]=1
else F[k,ndim*nvar+i]=psi[i]
}
else{
F[k,ndim*(i-1)+(1:ndim)]=Lambda[j,]
F[k,ndim*(j-1)+(1:ndim)]=Lambda[i,]
}
}
} # end of LH
else{
npair=nvar*nvar
# First, calculate d vec(Sigma) / d vec(Lambda')
F=matrix(0,npair,nvar*ndim+nvar)
#
for( k in 1:npair ){
# col-major rollout
j=trunc((k-1)/nvar)+1; i=((k-1)%%nvar)+1
if( i == j ){
F[k,ndim*(i-1)+(1:ndim)]=2*Lambda[i,]
if( logpsi == 0 ) F[k,ndim*nvar+i]=1
else F[k,ndim*nvar+i]=psi[i]
}
else{
F[k,ndim*(i-1)+(1:ndim)]=Lambda[j,]
F[k,ndim*(j-1)+(1:ndim)]=Lambda[i,]
}
}
} # end of LH!=0
#
# F0=matrix(as.numeric(abs(F)>0),npair)
#
# swap the columns
# d vec(S) / d vec(Lambda') -> d vec(S) / d vec(Lambda)
# od=c(t(matrix(1:(nvar*ndim),ndim,nvar)))]
od=matindex(ndim,nvar, type=1, byrow=1)
F[,1:(nvar*ndim)]=F[,od]
#
# F1=matrix(as.numeric(abs(F)>0),npair)
#
return( F )
#
} # end of Jacfa_a
#
#
#
Jacfa_aa <- function( param, logpsi=0, LH=0 ){
# returns analytic Jacobian of vec(S-Sigma) as gradient attribute
# Shin-ichi Mayekawa
# 20220728
# uselogpsiasparam -> logpsi: 20220728
#
res=Jacfa_a( param, logpsi=logpsi, LH=LH )
res0=0
attributes(res0)$gradient=res
return( res0 )
} # end of Jacfa_aa
#
#
gradfa_Ja <- function( param, logpsi=0, method="ML", LH=0 ){
# analytic gradient of LS or ML critfa wrt param=c(c(Lambda,psi)
# If method="LS" Jacfa_a is used.
# Shin-ichi Mayekawa
# 20220724,25
# LH: 20220728
# uselogpsiasparam -> logpsi: 20220728
# Should be the same as gradfa_a
#
#
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=param[-(1:(nvar*ndim))]
if( logpsi ) psi=exp(psi)
Sigma=Lambda%*%Phi%*%t(Lambda)+diag(psi)
#
if( method == "LS" ){
# by Jacobian: applicable for vec(S) and vech(S)
F=Jacfa_a( param, logpsi=logpsi, LH=LH )
# gradient
g=c( -2*t(F)%*%(S-Sigma)[locLH] )
}
else{
invSigma=matSwp(Sigma)
iSSmSiS=invSigma%*%(Sigma-S)%*%invSigma
dLambda=2*iSSmSiS%*%Lambda%*%Phi
dpsi=diag(iSSmSiS)
if( logpsi == 1 ) dpsi=dpsi*psi
g=c(c(dLambda),dpsi)
}
#
return( g )
#
} # gradfa_Ja
#
#
#
#
#
#
#
#
# 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)
#
# const to be used later
#
# log determinant of S
logdetS=log(det(S))
#
# lower half of S: NOT used.
locLH=lower.tri(S,diag=TRUE)
# vecS=S[locLH]
npair=sum(locLH)
locvLH=cbind(row(S)[locLH],col(S)[locLH])
#
#
#
# # of dimensions to be used
ndim=2
#
#
#
# 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)) )
#
# factor correlation matrix (constant)
Phi=diag(ndim)
invPhi=diag(ndim)
invPhimh=diag(ndim)
#
#
# initial value
param0=c(c(Lambda0),psi0)
param00=c(c(Lambda0),log(psi0))
#
#
#
#
#
# Maximum Likelihood Solution
#
#
#
#
#
#
# ML fa by nlminb: param=c( c(Lambda),psi )
resnlminb=nlminb( param0, critfa, control=list(trace=0)
, lower=1e-3, upper=1-1e-3 )
critnlminb=resnlminb$objective
param=resnlminb$par
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=param[-(1:(nvar*ndim))]
Print(critnlminb)
#
# ML fa by nlminb with analytic gradient: param=c( c(Lambda),psi )
resnlminb=nlminb( param0, critfa, control=list(trace=0)
, lower=1e-3, upper=1-1e-3, gradient=gradfa_a )
critnlminb=resnlminb$objective
param=resnlminb$par
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=param[-(1:(nvar*ndim))]
Print(critnlminb)
#
#
# ML fa by NR: param=c( c(Lambda),psi )
resNR=NR( param0, critfa, 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(critnlminb,critNR)
#
# ML fa by NR with analytic gradient: param=c( c(Lambda),psi )
resNR=NR( param0, critfa, 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(critnlminb,critNR)
#
#
#
#
# ML fa by nlminb: param=c( c(Lambda),log(psi) )
resnlminb=nlminb( param00, critfa, control=list(trace=0)
, logpsi=1 )
critnlminb=resnlminb$objective
param=resnlminb$par
LambdaML=matrix(param[1:(nvar*ndim)],nvar)
psiML=exp( param[-(1:(nvar*ndim))] )
Print(critnlminb)
#
# ML fa by nlminb with analytic gradient: param=c(c(Lambda),log(psi))
resnlminb=nlminb( param00, critfa, control=list(trace=0)
, logpsi=1, gradient=gradfa_a )
critnlminb=resnlminb$objective
param=resnlminb$par
LambdaML=matrix(param[1:(nvar*ndim)],nvar)
psiML=exp( param[-(1:(nvar*ndim))] )
Print(critnlminb)
#
#
# ML fa by NR: param=c( c(Lambda),log(psi) )
resNR=NR( param00, critfa, logpsi=1, maxiter=50, flipd=1 )
critNR=resNR$objective
param=resNR$par
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=exp( param[-(1:(nvar*ndim))] )
Print(critnlminb,critNR)
#
#
# ML fa by NR with analytic gradient: param=c( c(Lambda),log(psi) )
resNR=NR( param00, critfa, logpsi=1, maxiter=50, flipd=1
, gradient=gradfa_a )
critNR=resNR$objective
param=resNR$par
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=exp( param[-(1:(nvar*ndim))] )
Print(critnlminb,critNR)
#
#
#
#
#
#
#
# ML fa by nlminb: param=psi
# with constraints on psi. (no effect)
resnlminb=nlminb( psi0, critfapsi, control=list(trace=0)
, lower=0.02, upper=0.98 )
critnlminb=resnlminb$objective
psinlminb=resnlminb$par
Print(critnlminb,psinlminb)
#
g1=gradfapsi( resnlminb$par )
g2=gradfapsi_a( resnlminb$par )
Print(g1,g2)
#
#
# ML fa by nlminb: param=log(psi)
resnlminb=nlminb( log(psi0), critfapsi, control=list(trace=0)
, logpsi=1 )
critnlminb=resnlminb$objective
psinlminb=exp( resnlminb$par )
Print(critnlminb,psinlminb)
#
g1=gradfapsi( resnlminb$par, logpsi=1 )
g2=gradfapsi_a( resnlminb$par, logpsi=1 )
Print(g1,g2)
#
#
# ML fa by NR: param=psi
resNR=NR( psi0, critfapsi, maxiter=50, flipd=1 )
critNR=resNR$objective
psiNR=resNR$par
Print(critNR,psiNR)
#
g1=gradfapsi( resNR$par )
g2=gradfapsi_a( resNR$par )
Print(g1,g2)
#
#
# ML fa by NR with analytic gradient: param=psi
resNR=NR( psi0, critfapsi, maxiter=50, flipd=1, gradient=gradfapsi_a )
critNR=resNR$objective
psiNR=resNR$par
Print(critNR,psiNR)
#
g1=gradfapsi( resNR$par )
g2=gradfapsi_a( resNR$par )
Print(g1,g2)
#
#
#
# ML fa by NR: param=log(psi)
resNR=NR( log(psi0), critfapsi, logpsi=1, maxiter=50, flipd=1 )
critNR=resNR$objective
psiNR=exp( resNR$par )
Print(critNR,psiNR)
#
g1=gradfapsi( resNR$par, logpsi=1 )
g2=gradfapsi_a( resNR$par, logpsi=1 )
Print(g1,g2)
#
#
# ML fa by NR with analytic gradient: param=log(psi)
resNR=NR( log(psi0), critfapsi, logpsi=1, maxiter=50, flipd=1
, gradient=gradfapsi_a )
critNR=resNR$objective
psiNR=exp( resNR$par )
Print(critNR,psiNR)
#
g1=gradfapsi( resNR$par, logpsi=1 )
g2=gradfapsi_a( resNR$par, logpsi=1 )
Print(g1,g2)
#
#
Print(critnlminb,critNR)
#
#
#
#
#
#
# Least Squares Solution by Newton-Raphson Method
#
#
#
#
#
#
# LS fa by nlminb: param=c(c(Lambda),psi)
resnlminb=nlminb( param0, critfa, control=list(trace=0)
, logpsi=0, method="LS" )
critnlminb=resnlminb$objective
param=resnlminb$par
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=param[-(1:(nvar*ndim))]
Print(critnlminb)
#
g1=gradfa( param, method="LS" )
g2=gradfa_a( param, method="LS" )
Print(g1,g2)
#
#
# LS fa by nlminb: param=psi
resnlminb=nlminb( psi0, critfapsi, control=list(trace=0)
, method="LS" )
critnlminb=resnlminb$objective
Print(critnlminb,psinlminb)
g1=gradfapsi( resnlminb$par, method="LS" )
g2=gradfapsi_a( resnlminb$par, method="LS" )
Print(g1,g2)
#
#
#
# LS fa by nlminb with reparametrization: param=log(psi)
resnlminb=nlminb( log(psi0), critfapsi, control=list(trace=0)
, logpsi=1, method="LS" )
critnlminb=resnlminb$objective
psinlminb=exp( resnlminb$par )
Print(critnlminb,psinlminb)
g1=gradfapsi( resnlminb$par, logpsi=1, method="LS" )
g2=gradfapsi_a( resnlminb$par, logpsi=1, method="LS" )
Print(g1,g2)
#
#
#
#
# LS fa by NR with analytic g and H: param=c(c(Lambda),psi)
resNR=NR( param0, critfa, logpsi=0, maxiter=500
, flipd=1, method="LS", gradient=gradfa_a, hessian=HessQfa_a )
critNR=resNR$objective
param=resNR$par
Print(critNR,critnlminb)
#
g1=gradfa( resNR$par, logpsi=0, method="LS" )
g2=gradfa_a( resNR$par, logpsi=0, method="LS" )
Print(g1,g2)
#
#
#
# LS by NR with analytic g and H:
# param=c(c(Lambda),log(psi))
resNR=NR( param00, critfa, logpsi=1, maxiter=500
, flipd=1, method="LS", gradient=gradfa_a, hessian=HessQfa_a )
critNR=resNR$objective
param=resNR$par
Print(critNR,critnlminb)
#
g1=gradfa( resNR$par, logpsi=1, method="LS" )
g2=gradfa_a( resNR$par, logpsi=1, method="LS" )
Print(g1,g2)
#
#
#
#
# LS fa by NR: param=log(psi)
resNR=NR( log(psi0), critfapsi, logpsi=1, maxiter=50
, flipd=1, method="LS" )
critNR=resNR$objective
psiNR=exp( resNR$par )
Print(critNR,critnlminb,psiNR)
#
g1=gradfapsi( resNR$par, logpsi=1, method="LS" )
g2=gradfapsi_a( resNR$par, logpsi=1, method="LS" )
Print(g1,g2)
#
#
#
#
#
#
# Least Squares Solution by Gauss-Newton Method
#
#
#
#
# LS fa by nlfb param=c(c(Lambda),psi)
resnlfb=nlfb( param0,residfa2,control=list(femax=1000,japprox="jacentral"))
param=resnlfb$coefficients
critnlfb=resnlfb$ssquares
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=param[-(1:(nvar*ndim))]
Print(critnlfb)
#
g1=gradfa( param, logpsi=0, method="LS" )
g2=gradfa_a( param, logpsi=0, method="LS" )
Print(g1,g2)
#
#
# LS fa by nlfb with analytic jacobian: param=c(c(Lambda),psi)
resnlfb=nlfb( param0, residfa2, jacfn=Jacfa_aa, control=list(femax=1000) )
param=resnlfb$coefficients
critnlfb=resnlfb$ssquares
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=param[-(1:(nvar*ndim))]
Print(critnlfb)
#
g1=gradfa( param, logpsi=0, method="LS" )
g2=gradfa_a( param, logpsi=0, method="LS" )
Print(g1,g2)
#
#
#
# LS fa by nlfb param=c(c(Lambda),log(psi))
resnlfb=nlfb( param00, residfa2
, control=list(femax=1000,japprox="jacentral"), logpsi=1 )
param=resnlfb$coefficients
critnlfb=resnlfb$ssquares
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=exp( param[-(1:(nvar*ndim))] )
Print(critnlfb)
#
g1=gradfa( param, logpsi=1, method="LS" )
g2=gradfa_a( param, logpsi=1, method="LS" )
Print(g1,g2)
#
#
# LS fa by nlfb with analytic jacobian: param=c(c(Lambda),log(psi))
resnlfb=nlfb( param00, residfa2, control=list(femax=1000)
, logpsi=1, jacfn=Jacfa_aa )
param=resnlfb$coefficients
critnlfb=resnlfb$ssquares
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=exp( param[-(1:(nvar*ndim))] )
Print(critnlfb)
#
g1=gradfa( param, logpsi=1, method="LS" )
g2=gradfa_a( param, logpsi=1, method="LS" )
Print(g1,g2)
#
# LS fa by GN: param=c(c(Lambda),psi)
resGN=GN( param0, residfa, logpsi=0, maxiter=50 )
param=resGN$par
critGN=resGN$objective
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=param[-(1:(nvar*ndim))]
Print(critGN, critnlfb)
#
g1=gradfa( resGN$par, logpsi=0, method="LS" )
g2=gradfa_a( resGN$par, logpsi=0, method="LS" )
Print(g1,g2)
#
#
# LS fa by GN with analytic jacobian: param=c(c(Lambda),psi)
resGN=GN( param0, residfa, logpsi=0, maxiter=50
, jacobian=Jacfa_a )
param=resGN$par
critGN=resGN$objective
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=param[-(1:(nvar*ndim))]
Print(critGN, critnlfb)
#
g1=gradfa( resGN$par, logpsi=0, method="LS" )
g2=gradfa_a( resGN$par, logpsi=0, method="LS" )
Print(g1,g2)
#
#
# LS fa by GN: param=c(c(Lambda),log(psi))
resGN=GN( param00, residfa, logpsi=1, maxiter=50 )
param=resGN$par
critGN=resGN$objective
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=exp( param[-(1:(nvar*ndim))] )
Print(critGN, critnlfb)
#
g1=gradfa( resGN$par, logpsi=1, method="LS" )
g2=gradfa_a( resGN$par, logpsi=1, method="LS" )
Print(g1,g2)
#
#
#
#
# LS fa by GN with analytic jacobian: param=c(c(Lambda),log(psi))
resGN=GN( param00, residfa, logpsi=1, maxiter=50, jacobian=Jacfa_a )
param=resGN$par
critGN=resGN$objective
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=exp( param[-(1:(nvar*ndim))] )
Print(critGN, critnlfb)
#
g1=gradfa( resGN$par, logpsi=1, method="LS" )
g2=gradfa_a( resGN$par, logpsi=1, method="LS" )
Print(g1,g2)
#
#
#
#####################################################################
#
# End of Factor Analysis
#
#####################################################################
#
#
#
#
#
#
#
#
#
#
#####################################################################
#
# PCA by Newton-Raphson and Gauss-Newton method
#
#####################################################################
#
#
#
# Given n x nvar Y matrix, find n x ndim F and nvar x ndim A such that
# rss = tr( t(Y - F \
# is minimized where ndim << min(n,nvar).
#
# The solution using svd or eigen decomposition is well known
# but here we try to minimize rss directly w.r.t c( c(F), c(A) )
# by Newton-Raphson or Gauss-Newton method.
#
# It's not efficient at all but works.
#
#
#
# We use nlfb function of nlsr package for comparison.
#
library(nlsr)
#
#
rssPCA <- function( vecPCA, ndim=2, Y=0 ){
# rss of PCA: ssq(Y-Yhat)=ssq(Y-F%*%t(A))
# recover F and A from vecPCA and evaluate rss
n=nrow(Y); nvar=ncol(Y)
F=matrix(vecPCA[1:(n*ndim)],n,ndim)
A=matrix(vecPCA[-(1:(n*ndim))],nvar,ndim)
rss=ssq(Y-F%*%t(A))
return( rss )
} # end of rssPCA
#
#
drssPCA <- function( vecPCA, ndim=2, Y=0 ){
# gradient of rss of PCA
# recover F and A from vecPCA and evaluate rss
n=nrow(Y); nvar=ncol(Y)
F=matrix(vecPCA[1:(n*ndim)],n,ndim)
A=matrix(vecPCA[-(1:(n*ndim))],nvar,ndim)
dF=2*(F%*%t(A)%*%A-Y%*%A)
dA=2*(A%*%t(F)%*%F-t(Y)%*%F)
res=c(c(dF),c(dA))
return( res )
} # end of drssPCA
#
#
ymyhatPCA <- function( vecPCA, ndim=2, Y=0 ){
# residual as c(Y-Yhat)
# recover F and A from vecPCA and evaluate vec(Y-F%*%t(A))
n=nrow(Y); nvar=ncol(Y)
F=matrix(vecPCA[1:(n*ndim)],n,ndim)
A=matrix(vecPCA[-(1:(n*ndim))],nvar,ndim)
resid=c(Y-F%*%t(A))
return( resid )
} # end of ymyhatPCA
#
#
yhatmyPCA <- function( vecPCA, ndim=2, Y=0 ){
# residual as c(Yhat-Y) to be used by nlfb
# recover F and A from vecPCA and evaluate vec(F%*%t(A)-Y)
return( -ymyhatPCA( vecPCA, ndim=ndim, Y=Y ) )
} # end of yhatmyPCA
#
#
JacPCA <- function( vecPCA, ndim=2, Y=0 ){
# Jacobian: c( vec( d c(Yhat) / d F ), vec( d c(Yhat) / d A ) )
# recover F and A from vecPCA and evaluate vec(Y-F%*%t(A))
n=nrow(Y); nvar=ncol(Y)
F=matrix(vecPCA[1:(n*ndim)],n,ndim)
A=matrix(vecPCA[-(1:(n*ndim))],nvar,ndim)
#
# Jac
temp=NULL
for( i in 1: n ) temp=append(temp,list(A))
temp=Reduce( b_diag, temp )
# dF=temp[,c(t(matrix(1:(ndim*n),ndim,n)))]
# dF=dF[c(t(matrix(1:(nvar*n),nvar,n))),]
dF=temp[,matindex(ndim,n, type=1, byrow=1)]
dF=dF[matindex(nvar,n, type=1, byrow=1),]
temp=NULL
for( j in 1: nvar ) temp=append(temp,list(F))
temp=Reduce( b_diag, temp )
# dA=temp[,c(t(matrix(1:(ndim*nvar),ndim,nvar)))]
dA=temp[,matindex(ndim,nvar, type=1, byrow=1)]
#
res=cbind( dF, dA )
#
return( res )
} # end of JacPCA
#
#
#
normalizePCA <- function( vecPCA, ndim=2, Y=0 ){
# normalize X in each of GN/NR iteration
# so that F'F=n I and A'A=diag and A[1,] > 0
# 20220723,30
#
# convert vecX to matrix X
n=nrow(Y); nvar=ncol(Y)
F=matrix(vecPCA[1:(n*ndim)],n,ndim)
A=matrix(vecPCA[-(1:(n*ndim))],nvar,ndim)
#
# normalize
res=normalize_config(F,center=0, rotate=1, A=A)
#
F=res$F; A=res$A
for( a in 1:ndim ){
if( A[1,a] < 0 ){
A[,a]=-A[,a]
F[,a]=-F[,a]
}
}
res=list(F=F,A=A)
#
return( unlist(res) )
} # end of normalizePCA
#
#
#
#
# generate data
seed <- 1701
set.seed(seed)
n <- 50; nvar <- 9
Y <- matrix(rnorm(n*nvar),n)%*%(diag(1/(1:nvar)))
# centering
Y=scale( Y, center=TRUE, scale=FALSE )
#
#
#
# initial value
ndim=2
#
set.seed(1701)
F0=matrix(rnorm(n*ndim),n,ndim)
A0=t( solve(t(F0)%*%F0)%*%t(F0)%*%Y )
vecPCA0=c(c(F0),c(A0))
rss0=ssq(Y-F0%*%t(A0))
#
#
#
#
# Estimation starts here.
#
# Using all of F and A as parameters
# normalizing the solution after each NR iteration.
# with or without analytic gradient or Jacobian.
#
#
#
# Newton-Raphson
#
#
#
# by nlminb
cat("\n\n**** Starting PCA by NR by nlminb with analytic gradient ****\n\n")
resnlm=nlminb( vecPCA0, rssPCA, control=list(trace=0), Y=Y, ndim=ndim
, gradient=drssPCA )
vecPCAnlm=resnlm$par
F=matrix(vecPCAnlm[1:(n*ndim)],n,ndim)
A=matrix(vecPCAnlm[-(1:(n*ndim))],nvar,ndim)
resnorm=normalize_config(F,center=0, rotate=1, A=A)
Fnlm=resnorm$F; Anlm=resnorm$A
rssnlm=ssq(Y-Fnlm%*%t(Anlm))
gradnlm=drssPCA( vecPCAnlm, ndim=ndim, Y=Y )
Print(rss0,rssnlm, min(abs(gradnlm)))
#Print(A1,Anlm)
#
#
# by NR
cat("\n\n**** Starting PCA by NR by NR with analytic gradient **** \n\n")
resNR=NR( vecPCA0, rssPCA, Y=Y, ndim=ndim, normalize=normalizePCA
, gradient=drssPCA, flipd=1, maxiter2=3, eps=-9, print=2 )
vecPCANR=resNR$par
F=matrix(vecPCANR[1:(n*ndim)],n,ndim)
A=matrix(vecPCANR[-(1:(n*ndim))],nvar,ndim)
resnorm=normalize_config(F,center=0, rotate=1, A=A)
FNR=resnorm$F; ANR=resnorm$A
rssNR=ssq(Y-FNR%*%t(ANR))
gradNR=drssPCA( vecPCANR, ndim=ndim, Y=Y )
Print(rss0,rssNR, min(abs(gradNR)))
# Print(A1,ANR)
#
#
#
# Gauss-Newton
#
#
#
# by nlsr
cat("\n\n**** Startin PCA by GN by nlfb ****\n\n")
resnlsr=nlfb( vecPCA0, yhatmyPCA, Y=Y, ndim=ndim
, control=list(femax=500,japprox="jacentral") )
vecPCAnlsr=resnlsr$coefficients
F=matrix(vecPCAnlsr[1:(n*ndim)],n,ndim)
A=matrix(vecPCAnlsr[-(1:(n*ndim))],nvar,ndim)
resnorm=normalize_config(F,center=0, rotate=1, A=A)
Fnlsr=resnorm$F; Anlsr=resnorm$A
rssnlsr=ssq(Y-Fnlsr%*%t(Anlsr))
gradnlsr=drssPCA( vecPCAnlsr, ndim=ndim, Y=Y )
Print(rss0,rssnlsr, min(abs(gradnlsr)))
#
#
#
# by GN
cat("\n\n**** Starting PCA by GN w/o analytic Jacobian ****\n\n")
resGN=GN( vecPCA0, ymyhatPCA, Y=Y, ndim=ndim, normalize=normalizePCA
, print=2 )
vecPCAGN=resGN$par
F=matrix(vecPCAGN[1:(n*ndim)],n,ndim)
A=matrix(vecPCAGN[-(1:(n*ndim))],nvar,ndim)
resnorm=normalize_config(F,center=0, rotate=1, A=A)
FGN=resnorm$F; AGN=resnorm$A
rssGN=ssq(Y-FGN%*%t(AGN))
gradGN=drssPCA( vecPCAGN, ndim=ndim, Y=Y )
Print(rss0,rssGN, min(abs(gradGN)))
#
#
# by GN
cat("\n\n**** Starting PCA by GN with analytic Jacobian ****\n\n")
resGN=GN( vecPCA0, ymyhatPCA, Y=Y, ndim=ndim, normalize=normalizePCA
, jacobian=JacPCA, print=2 )
vecPCAGN=resGN$par
F=matrix(vecPCAGN[1:(n*ndim)],n,ndim)
A=matrix(vecPCAGN[-(1:(n*ndim))],nvar,ndim)
resnorm=normalize_config(F,center=0, rotate=1, A=A)
FGN=resnorm$F; AGN=resnorm$A
rssGN2=ssq(Y-FGN%*%t(AGN))
gradGN2=drssPCA( vecPCAGN, ndim=ndim, Y=Y )
Print(rss0,rssGN2, min(abs(gradGN2)))
#
#
#
# final results
# rss
rss=cbind(rssNR,rssnlm, rssGN,rssGN2, rssnlsr)
cat("\n\n Comparison of final rss\n")
Print(rss)
grad=cbind( min(abs(gradnlm)),min(abs(gradNR)),min(abs(gradGN))
, min(abs(gradGN2)), min(abs(gradnlsr)) )
colnames(grad)=c("NR", "nlm", "GN", "GN2", "nlsr")
cat("\n\n Comparison of final gradient\n")
Print(grad)
#
#
#####################################################################
#
# End of PCA by Newton-Raphson and Gauss-Newton method
#
#####################################################################
#
## End(Not run)