##### R-code for Example 7.6. It writes the related data to 
##### the file "example76.dat." It also makes "fig79.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(20,mu1,Sigma),mvrnorm(10,mu2,Sigma))

n = length(x[,1])
p = 3
lambda = 0.2               # Specify the weighting parameter

xbar=matrix(0,n,p)        # Xbar_n
Sn2=array(0,c(n,p,p))     # S_n^2
un=matrix(0,n,p)          # u_n
En=matrix(0,n,p)          # Self-starting charting statistic E_{n,SS}
Tn2=rep(0,n)              # T_{n,SS}^2

xbar[2,]=(x[1,]+x[2,])/2
Sn2[2,,]=(x[2,]-x[1,])%*%t(x[2,]-x[1,])/2

for(i in 3:4){
   xbar[i,]=((i-1)*xbar[i-1,]+x[i,])/i
   un[i,]=(x[i,]-xbar[i-1,])*sqrt((i-1)/i)
   Sn2[i,,]=Sn2[i-1,,]*(i-2)/(i-1)+(un[i,]%*%t(un[i,]))/(i-1)
}

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

xbar[i,]=((i-1)*xbar[i-1,]+x[i,])/i
un[i,]=(x[i,]-xbar[i-1,])*sqrt((i-1)/i)
Sn2[i,,]=Sn2[i-1,,]*(i-2)/(i-1)+(un[i,]%*%t(un[i,]))/(i-1)

En[i,]=lambda*un[i,]+(1-lambda)*En[i-1,]
Tn2[i]=(t(En[i,]) %*% solve(Sn2[i-1,,]) %*% En[i,])*
       (2-lambda)/(lambda*(1-(1-lambda)^(2*i)))

}

write.table(round(cbind(x,un,En,Tn2),digits=3),"example76.dat",row.names=F)

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

ii = seq(1,n)

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

graphics.off()