cor.fd <- function(evalarg1, fdobj1, evalarg2=evalarg1, fdobj2=fdobj1)
{
# compute correlation matrix / matrices for functional observations
# Last modified 16 January 2010
if (!(is.fd(fdobj1) || is.fdPar(fdobj1)))
stop("Second argument is neither an fd object or an fdPar object")
if (!(is.fd(fdobj2) || is.fdPar(fdobj2)))
stop("Fourth argument is neither an fd object or an fdPar object")
##
## 1. Compute var1 = bivariate data object for variance(fdobj1)
##
var1 <- var.fd(fdobj1)
##
## 2. Evaluate var1 at evalarg1
##
print(class(var1))
evalVar1 <- eval.bifd(evalarg1, evalarg1, var1)
##
## 3. If missing(fdobj2) convert evalVar1 to correlations
##
{
if(missing(fdobj2)){
dimV1 <- dim(evalVar1)
ndV1 <- length(dimV1)
{
if(ndV1<3){
s1 <- sqrt(diag(evalVar1))
return(evalVar1/outer(s1, s1))
}
else{
if(dimV1[3] != 1)
stop("Bug in cor.fd: Programmed only for ",
"matrices or 4-d arrays with dim[3]==1. oops.")
dim1 <- dim(fdobj1$coefs)
nVars <- dim1[3]
# The following identifies the levels of evalVar1
# containing the covariance matrix of each variable with itself
evalV.diag <- choose(2:(nVars+1), 2)
# Compute the standard deviation vector for each variable
nPts <- length(evalarg1)
s1 <- array(NA, dim=c(nPts, nVars))
for(i in 1:nVars)
s1[, i] <- sqrt(diag(evalVar1[,,1,evalV.diag[i]]))
# Now compute the correlations
cor1 <- evalVar1
m <- 0
for(i in 1:nVars)for(j in 1:i){
m <- m+1
cor1[,,1,m] <- (evalVar1[,,1,m] / outer(s1[, i], s1[, j]))
}
# end for i in 1:nVars and j in 1:i
return(cor1)
}
# end if...else nV1>2
}
}
else {
##
## 4. fdobj2 was also provided
##
# 4.1 var.df(fdobj2)
var2 <- var.fd(fdobj2)
# 4.2. Evalate at evalarg2
evalVar2 <- eval.bifd(evalarg2, evalarg2, var2)
# 4.3. var12 cross covariance
#*** If fdobj1 or fdobj2 are multivariate, var.fd will complain.
var12 <- var.fd(fdobj1, fdobj2)
# 4.4. Evaluate the cross covariances
evalVar12 <- eval.bifd(evalarg1, evalarg2, var12)
# 4.5. Convert evalVar12 to correlations
s1 <- sqrt(diag(evalVar1))
s2 <- sqrt(diag(evalVar2))
return(evalVar12 / outer(s1, s2))
}
}
##
## Done: Nothing should get this far;
## all should have returned earlier,
## either via univariate or multivariate
## with fdobj1 only
## or with both fdobj1 and fdobj2
##
}