Commit 7ae85a8f by Eric Coissac

More efficient version of the permutation simulation

parent 924f9e93
......@@ -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)
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
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment