"contrast" <- function(da,cmtx,alpha=0.05){
# The data matrix is "da".
# Written by Ruey S. Tsay for Bus 41912 class.
#
 nr = dim(da)[1]
 nc = dim(da)[2]
 ave=matrix(colMeans(da),1,nc)
 S=cov(da)
 cm = cmtx%*%t(ave)
 co=cmtx%*%S%*%t(cmtx)
 Si=solve(co)
 Tsq = nr*(t(cm)%*%Si%*%cm)
 qm1=dim(cmtx)[1]
 tmp = Tsq*(nr-qm1)/((nr-1)*qm1)
 deg2=nr-qm1
 pv=1-pf(tmp,qm1,deg2)
 print("Hotelling Tsq statistics & p-value")
 print(c(Tsq,pv))
 print("Simultaneous C.I. for each contrast")
 cri = sqrt((nr-1)*qm1*qf(1-alpha,qm1,deg2)/deg2)
 simucr=matrix(0,qm1,2)
 for (i in 1:qm1){
 c=cmtx[i,]
 me=t(c)%*%t(ave)
 s = t(c)%*%S%*%c
 se=sqrt(s)/sqrt(nr)
 simucr[i,1]=me-cri*se
 simucr[i,2]=me+cri*se
}
 print(simucr)
}