##### R-code for Example 7.5. It writes the related data to 
##### the file "example75.dat." It also makes "fig78.ps"

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

library(MASS) # The R library in which the command for generating multivariate
              # normal random vectors is contained.

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))

### Set parameters in the MEWMA chart
lambda=0.2        # smoothing parameter
h=11.81891        # control limit

n = length(x[,1])
p = 3

Yn=matrix(0,n,p)     # the MEWMA charting statistic Yn
Tn2=rep(0,n)         # Tn2=Yn' Sigma^{-1} Yn

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

for(i in 2:n){       # monitoring at the i-th time point
   Yn[i,]=lambda*(x[i,]-mu1)+(1-lambda)*Yn[i-1,]
   Tn2[i]=((2-lambda)/lambda)*(t(Yn[i,])%*%solve(Sigma)%*%Yn[i,])
}


write.table(round(cbind(x,Yn,Tn2),digits=3),"example75.dat",row.names=F)


postscript("fig78.ps",width=4.5,height=4,horizontal=F)

ii = seq(1,n)

plot(ii,Tn2,type="o",lty=1,pch=16,xlab="n",ylab=expression(V[n]^2),mgp=c(2,1,0),
     xlim=c(0,n),ylim=c(0,40),cex=0.8)
lines(ii,rep(h,n),lty=2,cex=0.8)

graphics.off()