"parallel" <- function(subject=100,var=10,rep=100,cent=0.05) { r <- subject c <- var y <- matrix(c(1:r*c), nrow=r, ncol=c) ycor <- matrix(c(1:c*c), nrow=c, ncol=c) evpea <- NULL leg.txt <- "Pearson" # Simulation of k samples to obtain k random eigenvalues vectors # for Pearson correlation coefficients for (k in c(1:rep)) { y <- rnorm(y) y <- matrix(y, nrow=r, ncol=c) evpea <- rbind(evpea, eigen(cor(y))[[1]]) } # Temporay function to compute the standard error of a centile SEcentile <- function(sd, n = 100, p = 0.95) {return(sd/sqrt(n) * sqrt(p*(1-p))/dnorm(qnorm(p))) } # Summary statistics sprob <- c(cent) mevpea <- sapply(as.data.frame(evpea), mean) # Eigenvalues means sevpea <- sapply(as.data.frame(evpea), sd ) # Eigenvalues Standard deviations quant <- function(x, sprobs = sprobs) {return(as.vector(quantile(x, probs = sprob))) } qevpea <- sapply(as.data.frame(evpea), quant) # Eigenvalues centiles sqevpea <- sevpea sqevpea <- sapply(as.data.frame(sqevpea), SEcentile, n = rep, p = cent) # Standard error of the centiles # List of results return result <- list(eigen = data.frame(mevpea, sevpea, qevpea, sqevpea), subject = r, variables = c, centile = cent ) class(result) <- 'parallel' # For future use return(result) }