##### R-code for Example 7.8. It writes the related data to 
##### the file "example78.dat." It also makes "fig711.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(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 = Sigma1                                         ### OC covariance matrix   

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

n = length(x[,1])
p = 3
hn=rep(0,n)
for(i in 26:n){
   hn[i]=exp(2.706+0.230*p-(p+3)/50*log(i-25))
}

xbar=matrix(0,n,p)
WnInv=array(0,c(n,p,p))

### Compute the initial values of xbar and WnInv from the first 25 IC obs

xbar[25,] = colMeans(x[1:25,])
WnInv[25,,] = solve(var(x[1:25,])*(25-1))

### Start the phase II process monitoring

Tmax2 = rep(0,n)

for(i in 26:n){

    xbar[i,]=xbar[i-1,]+(x[i,]-xbar[i-1,])/i
    temp=WnInv[i-1,,]%*%(x[i,]-xbar[i-1,])%*%t(x[i,]-xbar[i-1,])%*%WnInv[i-1,,]
    temp1=c(t(x[i,]-xbar[i-1,])%*%WnInv[i-1,,]%*%(x[i,]-xbar[i-1,]))

    WnInv[i,,]=WnInv[i-1,,]-(i-1)*temp/(i*(1+(i-1)/i*temp1))

    Tji2 = rep(0,n)

    for(j in 25:(i-1)){
        cj=j*i/(i-j)
        Tji2[j] = (i-2)*cj*(t(xbar[j,]-xbar[i,])%*%WnInv[i,,]%*%(xbar[j,]-xbar[i,]))/
                  (1-cj*(t(xbar[j,]-xbar[i,])%*%WnInv[i,,]%*%(xbar[j,]-xbar[i,])))
    }

    Tmax2[i] = max(Tji2[25:(i-1)])

}

write.table(round(cbind(x,xbar,Tmax2,hn),digits=3), 
            file="example78.dat",col.names=T, row.names=F)

### Make a figure

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

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

ii = seq(1,n)

### Plots of X1, X2, X3

plot(ii,x[,1],type="o",lty=1,pch=16,xlab="n",
     ylab=expression(X[n1]),mgp=c(2,1,0),xlim=c(0,n), 
     ylim=c(-2.5,3.5), cex=0.8)
title(xlab="(a)",cex=0.9)

plot(ii,x[,2],type="o",lty=1,pch=16,xlab="n",
     ylab=expression(X[n2]),mgp=c(2,1,0),xlim=c(0,n), 
     ylim=c(-2.5,3.5), cex=0.8)
title(xlab="(b)",cex=0.9)

plot(ii,x[,3],type="o",lty=1,pch=16,xlab="n",
     ylab=expression(X[n3]),mgp=c(2,1,0),xlim=c(0,n), 
     ylim=c(-2.5,3.5), cex=0.8)
title(xlab="(c)",cex=0.9)

### Plot of the CPD chart

plot(ii[26:50],Tmax2[26:50],type="o",lty=1,pch=16, xlab="n",
     ylab=expression(T[list(max,n)]^2),mgp=c(2,1,0),xlim=c(26,n), 
     ylim=c(0,50),cex=0.8)
lines(ii[26:50],hn[26:50],lty=2,cex=0.8)

graphics.off()