solved problems from Wichern - Multivariate Statistical analysis

wichern6.21.R

delta0 <- c(0,0,0,0)
xbar1 <- (c(2.287, 12.600, 0.347, 14.830))
xbar2 <- (c(2.404, 7.155, 0.524,  12.840))
S1 <- matrix(c(.459, .254, -.026, -.244, .254, 27.465, -.589, -.267, -.206, -.589, .030, .102, -.244, -.267, .102, 6.854), nc=4)
S2 <- matrix(c(.944, -.089, .002, -.719, -.089, 16.432, -.400, 19.044, .002, -.400, .024, -.094, -.719, 19.044, -.094, 61.854), nc =4)
n1 <- n2 <- 20
p <- length(delta0)
alpha <- 0.05
Spooled <- ((n1-1)*S1 + (n2-1)*S2 )/(n1+n2-2)
T2 <- (1/n1+1/n2)^{-1}*t(xbar1- xbar2 - delta0)%*%solve(Spooled)%*%(xbar1 - xbar2)
critical <- (n1+n2-2)*p/(n1+n2-p-1)*qf(alpha, p, n1+n2-p-1, lower.tail=F)

twosampleinterval <- function(xbar1=xbar1,xbar2=xbar2, S1=S1,S2=S2,  alpha=alpha, n1=n1,n2=n2, option=c("Simul", "Bonf")){
  # create (1-alpha)% simultaneous interval
  Spooled <- ((n1-1)*S1 + (n2-1)*S2 )/(n1+n2-2)
  p <- length(xbar1)
  if(option=="Simul"){
    critical <- sqrt((n1+n2-2)*p/(n1+n2-p-1)*qf(alpha, p,n1+n2-p-1, lower.tail=F))
  }else{
    critical <- qt(alpha/2/p, n1+n2-2, lower.tail=F)	
  }
  output <- cbind(xbar1-xbar2 - critical*sqrt(diag(Spooled))*sqrt(1/n1+1/n2),  xbar1-xbar2 + critical*sqrt(diag(Spooled))*sqrt(1/n1+1/n2))
  for(i in 1:p){
    print(paste("Var", i,":(", zapsmall(output[i,1],3), ",", zapsmall(output[i, 2],3),")" , sep=""))
  }
  return(output)
}

computenu <- function(p=p, S1=S1, S2=S2, n1=n1, n2=n2){
  grand <-solve( 1/n1*S1 + 1/n2*S2)
  p1 <-  1/n1*S1%*%grand
  p2 <- 1/n2*S2%*%grand
  
  nu <- (p+p^2)/(1/n1*(sum(diag(p1%*%p1 )) +  sum(diag(p1))^2 )+1/n2*(sum(diag(p2%*%p2))  +  sum(diag(p2))^2))
  
  return(nu)
}

nu <- computenu(p=p,  S1=S1, S2=S2, n1=n1, n2=n2)
T2new <-  t(xbar1-xbar2 -delta0)%*%solve(1/n1*S1 + 1/n2*S2)%*%(xbar1-xbar2 -delta0)
criticalnew<- nu*p/(nu-p+1)*qf(alpha, p, nu-p+1, lower.tail=F)