https://github.com/cran/fda
Tip revision: ebcbf854f75ed787efdd6fabf5a270fae0ac6474 authored by J. O. Ramsay on 03 November 2009, 00:00:00 UTC
version 2.2.2
version 2.2.2
Tip revision: ebcbf85
varmx.pca.fd.R
varmx.pca.fd <- function(pcafd, nharm=scoresd[2], nx=501)
{
#
# Apply varimax to the first NHARM components of a pca.fd object.
# Evaluates the harmonics at NX equally spaced points.
#
# Returns:
# An object of class pcafd
# Note that pcafd is an oldClass type object
# Last modified 22 October 2009 by Jim Ramsay
if (!(inherits(pcafd, "pca.fd"))) stop(
"Argument PCAFD is not a pca.fd object.")
harmfd <- pcafd$harmonics
harmcoef <- harmfd$coefs
coefd <- dim(harmcoef)
ndim <- length(coefd)
scoresd <- dim(pcafd$scores)
if (nharm > scoresd[2]) nharm <- scoresd[2]
basisobj <- harmfd$basis
rangex <- basisobj$rangeval
x <- seq(rangex[1], rangex[2], length = nx)
delta <- x[2]-x[1]
harmmat <- eval.fd(x, harmfd)
# If fdmat is a 3-D array, stack into a matrix
if (ndim == 3) {
harmmatd <- dim(harmmat)
dimnames(harmmat) <- NULL
harmmat <- aperm(harmmat, c(1, 3, 2))
dim(harmmat) <- c(harmmatd[1] * harmmatd[3], harmmatd[2])
}
# compute rotation matrix for varimax rotation of harmmat
rotm <- varmx(harmmat[,1:nharm])
# rotate coefficients and scores
rotharmcoef <- harmcoef
if (ndim == 2)
rotharmcoef[,1:nharm] <- harmcoef[,1:nharm] %*% rotm
else
for(j in (1:coefd[3]))
rotharmcoef[,1:nharm,j] <- harmcoef[,1:nharm,j] %*% rotm
# rotate principal component scores
rotharmscrs <- pcafd$scores
rotharmscrs[,1:nharm] <- rotharmscrs[,1:nharm] %*% rotm
# compute proportions of variance
rotharmvar <- apply(rotharmscrs^2,2,sum)
varsum <- sum(rotharmvar)
propvar <- sum(pcafd$varprop)*rotharmvar/varsum
# modify pcafd object
rotharmfd <- harmfd
rotharmfd$coefs <- rotharmcoef
rotpcafd <- pcafd
rotpcafd$harmonics <- rotharmfd
rotpcafd$varprop <- propvar
rotpcafd$scores <- rotharmscrs
rotpcafd$rotmat <- rotm
return(rotpcafd)
}