From 7ae85a8fb362b2a9ab67699b9f3f1361157b4da2 Mon Sep 17 00:00:00 2001 From: Eric Coissac Date: Thu, 18 Apr 2019 21:04:46 +0200 Subject: [PATCH] More efficient version of the permutation simulation --- R/covls.R | 42 ++++++++++++++++++++++++++++++++---------- 1 file changed, 32 insertions(+), 10 deletions(-) diff --git a/R/covls.R b/R/covls.R index b3c26ac..05cd4ed 100644 --- a/R/covls.R +++ b/R/covls.R @@ -93,19 +93,41 @@ varls = function(...,nperm = 100,rcovls=FALSE) { CovXX2s = matrix(0, nrow = nX, ncol = nX) if (! is.null(nperm)) { - rCovXXs = array(0,dim=c(nX,nX)) - for (i in seq_len(nX)) { - for (j in i:nX) { - rCovXXs[i,j] = null_varls_sim(Xs[i],Xs[j],nperm) - rCovXXs[j,i] = rCovXXs[i,j] - CovXXs[i,j] = CovXXs[i,j] - rCovXXs[i,j] - } - } + rCovXXs = array(0,dim=c(nperm,nX,nX)) + vXs=vector(mode="list",nX) + lXs=vapply(Xs,FUN = function(x) prod(dim(x)),0) + dXs=vapply(Xs,FUN = function(x) dim(ortho(x)),c(0,0)) + + for (i in seq_len(nX)) + vXs[i]=list(t(sqrt(colMeans(Xs[i]^2) - colMeans(Xs[i])^2) %*% t(rep(1,nrow(Xs[i]))))) + + for (k in seq_len(nperm)) { + # Generate the random matrix for the kth simulation + rXs=vector(mode="list",nX) + for (i in seq_len(nX)) + rXs[i]= list(tcrossprod(scale(array(rnorm(lXs[i]), + dim = dXs[,i]) + ) * vXs[[i]])) + + for (i in seq_len(nX)) + for (j in i:nX) + rCovXXs[k,i,j] = .Trace(sqrtm(rXs[[i]] %*% rXs[[j]])) + + } } + rCovXXs = apply(rCovXXs, + MARGIN = c(2,3), + FUN=function(i) MASS::fitdistr(Re(i), + "normal")\$estimate['mean'] + ) + for (i in seq_len(nX)) - for (j in i:nX) - CovXXs[j,i] = CovXXs[i,j] + for (j in i:nX) { + CovXXs[i,j] = max(CovXXs[i,j] - rCovXXs[i,j],0) + CovXXs[j,i] = CovXXs[i,j] + rCovXXs[j,i] = rCovXXs[i,j] + } CovXXs = CovXXs / N colnames(CovXXs)=Xnames -- libgit2 0.26.0