sfericity.test<-function(n,s1,s2=NULL,estsigma=TRUE){
+ #### Performs a hypothesis test that a covariance matrix is of specified
+ #### form. Test is of the form H0: S1=sigma^2*S2. n is the number of
+ #### observations on which the sample covariance matrix is based.
+ #### If the input parameter estsigma is TRUE:
+ #### Perform test of the hypothesis that S1=sigma^2 S2, for unknown sigma.
+ #### If S2 not specified, assumed that S2=I. Reference is Basilevsky,
+ #### Statistical Factor Analysis and Related Methods, page 191.
+ #### If the input parameter estsigma is FALSE:
+ #### Perform test of the hypothesis that S1=S2. If S2 not specified,
+ #### assumed that S2=I. Reference is Seber, Multivariate Observations,
+ #### sec 3.5.4
+ #### Only the lower triangle+diagonal is required at entry, and the upper
+ #### triangle is ignored.
+ #### DAW July 2000
+ dname <- paste(substitute(s1))
+ p<-nrow(s1)
+ for (i in 1:(p-1)){for (j in ((i+1):p)){
+     s1[i,j]<-s1[j,i]
+     s2[i,j]<-s2[j,i] }}
+ if (!is.null(s2)){
+     b<-eigen(s2,symmetric=T,only.values=F)
+     r<-b$vectors %*% diag(1/sqrt(b$values))
+     s<-t(r) %*% s1 %*% r }
+ else { s<-s1 }
+
+ d<-eigen(s,symmetric=T,only.values=T)$values
+ ldet<-sum(log(d))
+ tr<-sum(d)
+
+ if (estsigma==TRUE){
+     sighat<-tr/p
+     cc<--(n-(2*p^2+p+2)/(6*p))*(ldet-p*log(tr/p))
+     statistic <- cc
+     sighat<-sighat
+     names(statistic) <- "L statistic"
+     parameter <- 0.5*(p+2)*(p-1)
+     names(parameter) <- "df"
+     rval<-list(data.name=dname,sighat=sighat,statistic=statistic,parameter=parameter,p.value=1-pchisq(statistic,parameter),method="Sphericity test") }
+ else {
+     cc<--n*(p+ldet-tr)
+     statistic <- cc
+     names(statistic) <- "L statistic"
+     parameter <- 0.5*(p+1)*p
+     names(parameter) <- "df"
+     rval<-list(data.name="",statistic=statistic,parameter=parameter,p.value=1-pchisq(statistic,parameter),method="Covariance equality test statistic")
+ }
+ class(rval) <- "htest"
+ return(rval)
+ }
> pw <- function(q,n) {
+    pdf <- function(w) { 1/2 * (n-2) * w^((n-3)/2) }
+    integrate(pdf,0,q)
+ }
>
> varcomp <- function(covmat,n) {
+    if (is.list(covmat)) {
+     if (length(covmat) < 2)
+         stop("covmat must be a list with at least 2 elements")
+     ps <- as.vector(sapply(covmat,dim))
+     if (sum(ps[1] == ps) != length(ps))
+         stop("all covariance matrices must have the same dimension")
+     p <- ps[1]
+         q <- length(covmat)
+         if (length(n) == 1)
+         Ng <- rep(n,q)
+     else if (length(n) == q)
+         Ng <- n
+     else
+         stop("n must be equal length(covmat) or 1")
+
+     DNAME <- deparse(substitute(covmat))
+    }
+
+    else
+     stop("covmat must be a list")
+
+    ng <- Ng - 1
+    Ag <- lapply(1:length(covmat),function(i,mat,n) { n[i] * mat[[i]] },mat=covmat,n=ng)
+    A <- matrix(colSums(matrix(unlist(Ag),ncol=p^2,byrow=T)),ncol=p)
+    detAg <- sapply(Ag,det)
+    detA <- det(A)
+    V1 <- prod(detAg^(ng/2))/(detA^(sum(ng)/2))
+    kg <- ng/sum(ng)
+    l1 <- prod((1/kg)^kg)^(p*sum(ng)/2) * V1
+    rho <- 1 - (sum(1/ng) - 1/sum(ng))*(2*p^2+3*p-1)/(6*(p+1)*(q-1))
+    w2 <- p*(p+1) * ((p-1)*(p+2) * (sum(1/ng^2) - 1/(sum(ng)^2)) - 6*(q-1)*(1-rho)^2) / (48*rho^2)
+    f <- 0.5 * (q-1)*p*(p+1)
+    STATISTIC <- -2*rho*log(l1)
+    PVAL <- 1 - (pchisq(STATISTIC,f) + w2*(pchisq(STATISTIC,f+4) - pchisq(STATISTIC,f)))
+    names(STATISTIC) <- "corrected lambda*"
+    names(f) <- "df"
+    RVAL <- structure(list(statistic = STATISTIC, parameter = f,p.value = PVAL, data.name = DNAME, method = "Equality of Covariances Matrices Test"),class="htest")
+    return(RVAL)
+ }
> sphericity.test<-function(n,s1,s2=NULL,estsigma=TRUE){