In spectrum0.mvar()'s PCA component now attempts to normalize the rot…

…ated components, discarding those that cannot be normalized (likely due to being redundant or numericanly unstable).

R/approx.hotelling.diff.test.R

@@ -313,6 +313,7 @@ spectrum0.mvar <- function(x, order.max=NULL, aic=is.null(order.max), tol=.Machi
   # TODO: What if the variable actually has a tiny magnitude?
   xscl <- apply(x, 2L, stats::sd)
   novar <- xscl < tol
+  p.var <- sum(!novar)
   x <- x[,!novar,drop=FALSE]
   x.full <- x.full[,!novar,drop=FALSE]
   xscl <- xscl[!novar]
@@ -332,8 +333,11 @@ spectrum0.mvar <- function(x, order.max=NULL, aic=is.null(order.max), tol=.Machi
   # s.min/s.max, where s.min and s.max are the smallest and the
   # biggest eigenvalues, respectively, is greater than the tolerance.
   e <- eigen(cov(x), symmetric=TRUE)
-  Q <- e$vectors[,sqrt(pmax(e$values,0)/max(e$values))>tol*2,drop=FALSE]
-  xr <- x.full%*%Q # Columns of xr are guaranteed to be linearly independent.
+  xr <- x.full %*% e$vectors %*% diag(1/sqrt(pmax(e$values,0)),p.var) # Columns of xr are guaranteed to be linearly independent with variance 1.
+  xr.sd <- apply(xr, 2, sd, na.rm=TRUE)
+  ind <- !is.na(xr.sd) & abs(xr.sd - 1) < tol
+  unmapper <- diag(xscl, p.var) %*% e$vectors[, ind, drop=FALSE] %*% diag(sqrt(e$values[ind]), sum(ind))
+  xr <- xr[, ind, drop=FALSE]
 
   ind.var <- cov(xr, use="complete.obs") # Get the sample variance of the transformed columns.
 
@@ -371,10 +375,9 @@ spectrum0.mvar <- function(x, order.max=NULL, aic=is.null(order.max), tol=.Machi
   infl <- exp((determinant(v.var)$modulus-determinant(ind.var)$modulus)/ncol(ind.var))
 
   # Reverse the mapping for the variance estimate.
-  v.var <- xAxT(xscl*Q, v.var)
+  v.var <- xAxT(unmapper, v.var)
 
   v[!novar,!novar] <- v.var
-
-  attr(v, "infl") <- infl
-  v
+
+  structure(v, infl = as.vector(infl), rank = sum(ind))
 }