# cvm.R
# Power for test of zero correlation for an entire matrix
# (Large Sample LR Test)
# Execute with   source("cvm.R")
#
M <- 10000
sim <- numeric(M)
set.seed(32448)
n <- 50 ; v1 <- .42 ;       v2 <- .18
          s1 <- sqrt(v1) ; s2 <- sqrt(v2)
G <- function(datamat)
   {
   nn <- dim(datamat)[1] ; kk <- dim(datamat)[2] ; df <- kk*(kk-1)/2
   G <- numeric(3)
   names(G) <- c("Chisq","df","P-value")
   S <- var(datamat)
   G[1] <- nn * ( sum(log(diag(S))) - sum(log(eigen(S)$values)) ) #$
   G[2] <- df
   G[3] <- 1 - pchisq(G[1],df)
   G
   } # End function G
merror <- function(phat,m,alpha=0.01) # (1-alpha)*100% merror for a proportion
    {
    z <- qnorm(1-alpha/2)
    merror <- z * sqrt(phat*(1-phat)/m)  # m is (Monte Carlo) sample size
    merror
    }

for(j in 1:M)
   {

e1 <- rnorm(n,0,s1) ; e2 <- rnorm(n,0,s2)
x1 <- rnorm(n)+e1 ; x2 <- rnorm(n)+e1 ;  x3 <- rnorm(n)+e1 ;
x4 <- rnorm(n)+e2 ;  x5 <- rnorm(n)+e2
dat <- cbind(x1,x2,x3,x4,x5)
# print(G(dat))
cat("Simulation ",j,"\n")
sim[j] <- G(dat)[3] < .05

    }
poww <- length(sim[sim==1])/M
cat("Power = ", poww ,"\n")
cat("Plus or Minus 99% Margin of error: ",merror(poww,M),"\n")