##### R-code for Example 7.7. It writes the related data to 
##### the file "example77.dat." It also makes "fig710.ps"

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

library(MASS)   ### R package for generating multivariate normal random vectors.
library(Matrix) ### R package for certain matrix manipulation (e.g., Cholesky
                ### decomposition of a square positive definite matrix).

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.5,0.2,0.5,1,0.5,0.2,0.5,1),3,3)  ### OC covariance matrix

### Note: In this example, correlations decreased and variances increased

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

n = length(x[,1])
p = 3
lambda = 0.2              # Specify the weighting parameter
h = 1.884                 # Specify the control limit for ARL_0=200.

U = chol(Sigma1)          # Cholesky decomposition of Sigma1: Sigma1=U'U
L = solve(t(U))           # L Sigma1 L' = I_{pxp}

y=matrix(0,n,p)
for(i in 1:n){
   y[i,] = L %*% (x[i,]-mu1)
}

En = array(0,c(n,3,3))

En[1,,]=lambda*(y[1,]%*%t(y[1,]))+(1-lambda)*diag(c(1,1,1))
for(i in 2:n){
   En[i,,]=lambda*(y[i,]%*%t(y[i,]))+(1-lambda)*En[i-1,,]
}

Cn=rep(0,n)             # Charting statistic of the MEWMA chart

for(i in 1:n){
   Cn[i] = sum(diag(En[i,,]))-log(det(En[i,,]))-p
}


write.table(round(cbind(x,y,Cn),digits=3),"example77.dat",row.names=F)

postscript("fig710.ps",width=4.5,height=5,horizontal=F)

ii = seq(1,n)

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

graphics.off()