########################################################## ############### Handout #3 for ST476/576 ############### ########## T^2 hyp test and Simultaneous CI's ############ ########################################################## ## We will use the 'ellipse' R-package in this handout ## To install this package use: install.packages('ellipse') ## You only have to install it once. However once installed it needs to be ## loaded each time you start R in order to use any of its functions. To load ## it use: library(ellipse) library(MASS) ################################################# ####################### Example 1 ############### #### Test of hyp, and CI's on Generated Data #### ################################################# #### Generate a random sample of n=100 from N(mu, Sigma) #### set.seed(24) mu <- c(0,.5,-1) Sigma <- rbind(c(5, 4, 5), c(4, 10, 3), c(5, 3, 15)) n <- 100 X <- mvrnorm(n, mu, Sigma) ## T2 test for mu = c(0,0,0) or not x_bar <- colMeans(X) S <- cov(X) n <- nrow(X) mu_0 <- c(0,0,0) T2 <- t(x_bar-mu_0)%*%solve(1/n*S)%*%(x_bar-mu_0) a <- (n-3)/((n-1)*3)*T2 pval <- 1-pf(a,3,n-3) pval ## 95% T2 intervals for individual means ## mu_1 a <- c(1,0,0) n <- nrow(X) p <- ncol(X) lower <- t(a)%*%x_bar - sqrt(p*(n-1)/(n*(n-p))*qf(.95,p,n-p)*t(a)%*%S%*%a) upper <- t(a)%*%x_bar + sqrt(p*(n-1)/(n*(n-p))*qf(.95,p,n-p)*t(a)%*%S%*%a) mu_1_CI <- c(lower,upper) mu_1_CI ## mu_2 a <- c(0,1,0) lower <- t(a)%*%x_bar - sqrt(p*(n-1)/(n*(n-p))*qf(.95,p,n-p)*t(a)%*%S%*%a) upper <- t(a)%*%x_bar + sqrt(p*(n-1)/(n*(n-p))*qf(.95,p,n-p)*t(a)%*%S%*%a) mu_2_CI <- c(lower,upper) mu_2_CI ## mu_3 a <- c(0,0,1) lower <- t(a)%*%x_bar - sqrt(p*(n-1)/(n*(n-p))*qf(.95,p,n-p)*t(a)%*%S%*%a) upper <- t(a)%*%x_bar + sqrt(p*(n-1)/(n*(n-p))*qf(.95,p,n-p)*t(a)%*%S%*%a) mu_3_CI <- c(lower,upper) mu_3_CI ## Plot T^2 Confidence Ellipse for (mu_1, mu_2)' crit <- (n-1)*p/(n*(n-p))*qf(.95,p,n-p) plot(ellipse(S[1:2,1:2], t=sqrt(crit), centre=x_bar[1:2]), type='l', main='95% T^2 Confidence Ellipse for (mu_1, mu_2)') ## Note: ellipse(S, t=t, centre=x) generates a 100 points to be plotted. ## These points are on the ellipse goverened by (mu-x)'S^(-1)(mu-x) < t^2. ## Hence plot(ellipse(S, t=t, centre=x), lty='l') says to make a plot and ## connect the points with a line. ## Add the point (mu_1, mu_2)=0 to the conf ellipse points(0,0) ## Add T^2 CI lines lines(c(mu_1_CI[1], mu_1_CI[1]), c(-100,100), lty=2) lines(c(mu_1_CI[2], mu_1_CI[2]), c(-100,100), lty=2) lines(c(-100,100), c(mu_2_CI[1], mu_2_CI[1]), lty=2) lines(c(-100,100), c(mu_2_CI[2], mu_2_CI[2]), lty=2) ## Plot T^2 Confidence Ellipse for (mu_2, mu_3)' plot(ellipse(S[2:3,2:3], t=sqrt(crit), centre=x_bar[2:3]), type='l', main='95% T^2 Confidence Ellipse for (mu_2, mu_3)', xlim=c(-1,1.2), ylim=c(-2.8,0)) ## Add the point (mu_2, mu_3)=0 to the conf ellipse points(0,0) ## Add T^2 CI lines lines(c(mu_2_CI[1], mu_2_CI[1]), c(-100,100), lty=2) lines(c(mu_2_CI[2], mu_2_CI[2]), c(-100,100), lty=2) lines(c(-100,100), c(mu_3_CI[1], mu_3_CI[1]), lty=2) lines(c(-100,100), c(mu_3_CI[2], mu_3_CI[2]), lty=2) ## 95% T2 interval for mu_1 - mu_2 a <- c(1,-1,0) lower <- t(a)%*%x_bar - sqrt(p*(n-1)/(n*(n-p))*qf(.95,p,n-p)*t(a)%*%S%*%a) upper <- t(a)%*%x_bar + sqrt(p*(n-1)/(n*(n-p))*qf(.95,p,n-p)*t(a)%*%S%*%a) mu_diff_CI <- c(lower,upper) mu_diff_CI ## 95% Bonferroni CI's for individual means ## mu_1 a <- c(1,0,0) n <- nrow(X) p <- ncol(X) lower <- t(a)%*%x_bar - qt(1-.05/(2*p), n-1)*sqrt(t(a)%*%S%*%a/n) upper <- t(a)%*%x_bar + qt(1-.05/(2*p), n-1)*sqrt(t(a)%*%S%*%a/n) mu_1_Bon_CI <- c(lower,upper) mu_1_Bon_CI ## mu_2 a <- c(0,1,0) lower <- t(a)%*%x_bar - qt(1-.05/(2*p), n-1)*sqrt(t(a)%*%S%*%a/n) upper <- t(a)%*%x_bar + qt(1-.05/(2*p), n-1)*sqrt(t(a)%*%S%*%a/n) mu_2_Bon_CI <- c(lower,upper) mu_2_Bon_CI ## mu_3 a <- c(0,0,1) lower <- t(a)%*%x_bar - qt(1-.05/(2*p), n-1)*sqrt(t(a)%*%S%*%a/n) upper <- t(a)%*%x_bar + qt(1-.05/(2*p), n-1)*sqrt(t(a)%*%S%*%a/n) mu_3_Bon_CI <- c(lower,upper) mu_3_Bon_CI ######################################### ############## Example 2 ################ ### Test of hyp, CR and T2 intervals #### ######################################### ## Read in data ex.data <- read.table(file="handout2_data.txt", header=T) ## Only the first 3 columns contain the X's of interest X <- as.matrix(ex.data[,1:3]) ## T2 test for mu = c(11,9,11) or not x_bar <- colMeans(X) S <- cov(X) n <- nrow(X) mu_0 <- c(11,9,11) T2 <- t(x_bar-mu_0)%*%solve(1/n*S)%*%(x_bar-mu_0) a <- (n-3)/((n-1)*3)*T2 pval <- 1-pf(a,3,n-3) pval ## 95% T2 intervals for individual means ## mu_1 a <- c(1,0,0) n <- nrow(X) p <- ncol(X) lower <- t(a)%*%x_bar - sqrt(p*(n-1)/(n*(n-p))*qf(.95,p,n-p)*t(a)%*%S%*%a) upper <- t(a)%*%x_bar + sqrt(p*(n-1)/(n*(n-p))*qf(.95,p,n-p)*t(a)%*%S%*%a) mu_1_CI <- c(lower,upper) mu_1_CI ## mu_2 a <- c(0,1,0) lower <- t(a)%*%x_bar - sqrt(p*(n-1)/(n*(n-p))*qf(.95,p,n-p)*t(a)%*%S%*%a) upper <- t(a)%*%x_bar + sqrt(p*(n-1)/(n*(n-p))*qf(.95,p,n-p)*t(a)%*%S%*%a) mu_2_CI <- c(lower,upper) mu_2_CI ## mu_3 a <- c(0,0,1) lower <- t(a)%*%x_bar - sqrt(p*(n-1)/(n*(n-p))*qf(.95,p,n-p)*t(a)%*%S%*%a) upper <- t(a)%*%x_bar + sqrt(p*(n-1)/(n*(n-p))*qf(.95,p,n-p)*t(a)%*%S%*%a) mu_3_CI <- c(lower,upper) mu_3_CI ## Plot T^2 Confidence Ellipse for (mu_1, mu_2) crit <- (n-1)*p/(n*(n-p))*qf(.95,p,n-p) plot(ellipse(S, t=sqrt(crit), centre=x_bar[1:2]), type='l', main='95% T^2 Confidence Ellipse for (mu_1, mu_2)') ## Add T^2 CI lines lines(c(mu_1_CI[1], mu_1_CI[1]), c(-100,100), lty=2) lines(c(mu_1_CI[2], mu_1_CI[2]), c(-100,100), lty=2) lines(c(-100,100), c(mu_2_CI[1], mu_2_CI[1]), lty=2) lines(c(-100,100), c(mu_2_CI[2], mu_2_CI[2]), lty=2) ## 95% T2 interval for mu_1 - mu_2 a <- c(1,-1,0) lower <- t(a)%*%x_bar - sqrt(p*(n-1)/(n*(n-p))*qf(.95,p,n-p)*t(a)%*%S%*%a) upper <- t(a)%*%x_bar + sqrt(p*(n-1)/(n*(n-p))*qf(.95,p,n-p)*t(a)%*%S%*%a) mu_diff_CI <- c(lower,upper) mu_diff_CI ## 95% Bonferroni CI's for individual means ## mu_1 a <- c(1,0,0) n <- nrow(X) p <- ncol(X) lower <- t(a)%*%x_bar - qt(1-.05/(2*p), n-1)*sqrt(t(a)%*%S%*%a/n) upper <- t(a)%*%x_bar + qt(1-.05/(2*p), n-1)*sqrt(t(a)%*%S%*%a/n) mu_1_Bon_CI <- c(lower,upper) mu_1_Bon_CI ## mu_2 a <- c(0,1,0) lower <- t(a)%*%x_bar - qt(1-.05/(2*p), n-1)*sqrt(t(a)%*%S%*%a/n) upper <- t(a)%*%x_bar + qt(1-.05/(2*p), n-1)*sqrt(t(a)%*%S%*%a/n) mu_2_Bon_CI <- c(lower,upper) mu_2_Bon_CI ## mu_3 a <- c(0,0,1) lower <- t(a)%*%x_bar - qt(1-.05/(2*p), n-1)*sqrt(t(a)%*%S%*%a/n) upper <- t(a)%*%x_bar + qt(1-.05/(2*p), n-1)*sqrt(t(a)%*%S%*%a/n) mu_3_Bon_CI <- c(lower,upper) mu_3_Bon_CI ######################################### ############## Example 3 ################ ####### Multivariate Control Chart ###### ######################################### #### Generate a random sample of n=150 from N(mu.1, Sigma) #### set.seed(22) mu.1<- c(0,0,0) Sigma <- rbind(c(5, 4, 5), c(4, 10, 3), c(5, 3, 15)) X1 <- mvrnorm(150, mu.1, Sigma) #### Generate a random sample of n=50 from N(mu.2, Sigma) #### mu.2 <- c(3,5,0) X2 <- mvrnorm(50, mu.2, Sigma) X <- rbind(X1, X2) #### Estimate mu and Sigma from historic data (first 150 observations) xbar <- colMeans(X1) S <- cov(X1) #### Create T^2 statistics and critical values n <- nrow(X) T2 <- rep(0, n) for(i in 1:n) T2[i] <- t(X[i,] - xbar)%*%solve(S)%*%(X[i,] - xbar) crit1 <- qchisq(.997,3) crit2 <- qchisq(.95,3) #### Plot control chart plot(1:n, T2, xlab="time", ylim=c(min(T2), max(T2,crit2))) lines(1:n, T2) abline(crit1, 0) ## 99.7% line abline(crit2, 0, lty=2) ## 95% line #### Plot univariate control charts lower <- xbar[1]+3*sqrt(S[1,1]) upper <- xbar[1]-3*sqrt(S[1,1]) plot(1:n, X[,1], xlab="time", ylab="X1", ylim=c(lower, upper)) lines(1:n, X[,1]) abline(xbar[1]+3*sqrt(S[1,1]), 0) abline(xbar[1]-3*sqrt(S[1,1]), 0) abline(xbar[1]+2*sqrt(S[1,1]), 0, lty=2) abline(xbar[1]-2*sqrt(S[1,1]), 0, lty=2) lower <- xbar[2]+3*sqrt(S[2,2]) upper <- xbar[2]-3*sqrt(S[2,2]) plot(1:n, X[,2], xlab="time", ylab="X2", ylim=c(lower, upper)) lines(1:n, X[,2]) abline(xbar[2]+3*sqrt(S[2,2]), 0) abline(xbar[2]-3*sqrt(S[2,2]), 0) abline(xbar[2]+2*sqrt(S[2,2]), 0, lty=2) abline(xbar[2]-2*sqrt(S[2,2]), 0, lty=2)