##### R-code for Example 7.4. It writes the related data to 
##### the file "example74a.dat", "example74b.dat", and
##### "example74c.dat". It also makes "fig77.ps"

### Case a: Process mean shift only

#### Generate 30 random vectors from N_3(mu,Sigma)

library(MASS)

set.seed(100)

mu1 = c(0,0,0)  ### IC mean vector
mu2 = c(1,0,0)  ### OC mean vector
Sigma = matrix(c(1,0.8,0.5,0.8,1,0.8,0.5,0.8,1),3,3)  ### IC covariance matrix

x = rbind(mvrnorm(10,mu1,Sigma),mvrnorm(20,mu2,Sigma))

n = length(x[,1])
p = 3
k1 = 4                    # Specify the allowance constant for CUSUM of T2
k2 = 1                    # Specify the allowance constant for CUSUM of T

SnT2=rep(0,n)             # Charting statistic of CUSUM of T2
SnT=rep(0,n)              # Charting statistic of CUSUM of T
Tn2=rep(0,n)              # Hotelling's T^2

Tn2[1]=t(x[1,]-mu1)%*%solve(Sigma)%*%(x[1,]-mu1)

SnT2[1]=max(0,Tn2[1]-k1)
SnT[1]=max(0,Tn2[1]^0.5-k2)

for(i in 2:n){       # monitoring at the i-th time point

Tn2[i]=t(x[i,]-mu1)%*%solve(Sigma)%*%(x[i,]-mu1)

SnT2[i]=max(0,SnT2[i-1]+Tn2[i]-k1)
SnT[i]=max(0,SnT[i-1]+(Tn2[i])^0.5-k2)

}

write.table(round(cbind(x,Tn2,SnT2,SnT),digits=3),"example74a.dat",row.names=F)

postscript("fig77.ps",width=6.5,height=7,horizontal=F)

par(mfrow=c(3,3), mar=c(4,4,2,2))

ii = seq(1,n)

### Plot of the Hotelling's T^2 Shewhart chart

plot(ii,Tn2,type="o",lty=1,pch=16,xlab="n",ylab=expression(T[list(0,n)]^2),
     mgp=c(2,1,0),xlim=c(0,n),ylim=c(0,15),cex=0.8)
lines(ii,rep(12.838,n),lty=2,cex=0.8)
title(xlab="(a)",cex=0.9)

### Plot of the CUSUM of T^2 chart

plot(ii,SnT2,type="o",lty=1,pch=16,xlab="n",ylab=expression(C[n]),
     mgp=c(2,1,0),xlim=c(0,n),ylim=c(0,60),cex=0.8)
lines(ii,rep(13.062,n),lty=2,cex=0.8)
title(xlab="(b)",cex=0.9)

### Plot of the CUSUM of T chart

plot(ii,SnT,type="o",lty=1,pch=16,xlab="n",ylab=expression(widetilde(C)[n]),
     mgp=c(2,1,0),xlim=c(0,n),ylim=c(0,32),cex=0.8)
lines(ii,rep(2.713,n),lty=2,cex=0.8)
title(xlab="(c)",cex=0.9)


### Case b: Covariance matrix shift only

set.seed(100)

mu1 = c(0,0,0)  ### IC mean vector
mu2 = c(0,0,0)  ### OC mean vector
Sigma1 = matrix(c(1,0.8,0.5,0.8,1,0.8,0.5,0.8,1),3,3)    ### IC covariance matrix
Sigma2 = 1.1*matrix(c(1,0.8,0.5,0.8,1,0.8,0.5,0.8,1),3,3)  ### OC covariance matrix

x = rbind(mvrnorm(10,mu1,Sigma1),mvrnorm(20,mu2,Sigma2))

n = length(x[,1])
p = 3
k1 = 4                    # Specify the allowance constant for CUSUM of T2
k2 = 1                    # Specify the allowance constant for CUSUM of T

SnT2=rep(0,n)             # Charting statistic of CUSUM of T2
SnT=rep(0,n)              # Charting statistic of CUSUM of T
Tn2=rep(0,n)              # Hotelling's T^2

Tn2[1]=t(x[1,]-mu1)%*%solve(Sigma1)%*%(x[1,]-mu1)

SnT2[1]=max(0,Tn2[1]-k1)
SnT[1]=max(0,(Tn2[1])^0.5-k2)

for(i in 2:n){       # monitoring at the i-th time point

Tn2[i]=t(x[i,]-mu1)%*%solve(Sigma1)%*%(x[i,]-mu1)

SnT2[i]=max(0,SnT2[i-1]+Tn2[i]-k1)
SnT[i]=max(0,SnT[i-1]+(Tn2[i])^0.5-k2)

}

write.table(round(cbind(x,Tn2,SnT2,SnT),digits=3),"example74b.dat",row.names=F)

### Plot of the Hotelling's T^2 Shewhart chart

plot(ii,Tn2,type="o",lty=1,pch=16,xlab="n",ylab=expression(T[list(0,n)]^2),
     mgp=c(2,1,0),xlim=c(0,n),ylim=c(0,15),cex=0.8)
lines(ii,rep(12.838,n),lty=2,cex=0.8)
title(xlab="(d)",cex=0.9)

### Plot of the CUSUM of T^2 chart

plot(ii,SnT2,type="o",lty=1,pch=16,xlab="n",ylab=expression(C[n]),
     mgp=c(2,1,0),xlim=c(0,n),ylim=c(0,16),cex=0.8)
lines(ii,rep(13.062,n),lty=2,cex=0.8)
title(xlab="(e)",cex=0.9)

### Plot of the CUSUM of T chart

plot(ii,SnT,type="o",lty=1,pch=16,xlab="n",ylab=expression(widetilde(C)[n]),
     mgp=c(2,1,0),xlim=c(0,n),ylim=c(0,20),cex=0.8)
lines(ii,rep(2.713,n),lty=2,cex=0.8)
title(xlab="(f)",cex=0.9)


### Case c: Both mean and covariance matrix shifts

set.seed(100)

mu1 = c(0,0,0)   ### IC mean vector
mu2 = c(1,0,0)   ### OC mean vector
Sigma1 = matrix(c(1,0.8,0.5,0.8,1,0.8,0.5,0.8,1),3,3)      ### IC covariance matrix
Sigma2 = 1.1*matrix(c(1,0.8,0.5,0.8,1,0.8,0.5,0.8,1),3,3)  ### OC covariance matrix

x = rbind(mvrnorm(10,mu1,Sigma1),mvrnorm(20,mu2,Sigma2))

n = length(x[,1])
p = 3
k1 = 4                    # Specify the allowance constant for CUSUM of T2
k2 = 1                    # Specify the allowance constant for CUSUM of T

SnT2=rep(0,n)             # Charting statistic of CUSUM of T2
SnT=rep(0,n)              # Charting statistic of CUSUM of T
Tn2=rep(0,n)              # Hotelling's T^2

Tn2[1]=t(x[1,]-mu1)%*%solve(Sigma1)%*%(x[1,]-mu1)

SnT2[1]=max(0,Tn2[1]-k1)
SnT[1]=max(0,(Tn2[1])^0.5-k2)

for(i in 2:n){       # monitoring at the i-th time point

Tn2[i]=t(x[i,]-mu1)%*%solve(Sigma1)%*%(x[i,]-mu1)

SnT2[i]=max(0,SnT2[i-1]+Tn2[i]-k1)
SnT[i]=max(0,SnT[i-1]+(Tn2[i])^0.5-k2)

}

write.table(round(cbind(x,Tn2,SnT2,SnT),digits=3),"example74c.dat",row.names=F)

### Plot of the Hotelling's T^2 Shewhart chart

plot(ii,Tn2,type="o",lty=1,pch=16,xlab="n",ylab=expression(T[list(0,n)]^2),
     mgp=c(2,1,0),xlim=c(0,n),ylim=c(0,15),cex=0.8)
lines(ii,rep(12.838,n),lty=2,cex=0.8)
title(xlab="(g)",cex=0.9)

### Plot of the CUSUM of T^2 chart

plot(ii,SnT2,type="o",lty=1,pch=16,xlab="n",ylab=expression(C[n]),
     mgp=c(2,1,0),xlim=c(0,n),ylim=c(0,60),cex=0.8)
lines(ii,rep(13.062,n),lty=2,cex=0.8)
title(xlab="(h)",cex=0.9)

### Plot of the CUSUM of T chart

plot(ii,SnT,type="o",lty=1,pch=16,xlab="n",ylab=expression(widetilde(C)[n]),
     mgp=c(2,1,0),xlim=c(0,n),ylim=c(0,32),cex=0.8)
lines(ii,rep(2.713,n),lty=2,cex=0.8)
title(xlab="(i)",cex=0.9)


graphics.off()