##### R-code for Example 9.1. It writes the related data to the file
##### "example91.dat." It also makes "fig93.ps."

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

library(mvtnorm) # The R library in which the commands for generating 
                 # multivariate t or normal random vectors is contained.

set.seed(100)

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

x = rbind(rmvt(n=500,sigma=Sigma,df=5),rmvt(n=500,sigma=Sigma,df=5))
x = x+rbind(matrix(rep(mu1,500),500,2,byrow=T),matrix(rep(mu2,500),500,2,byrow=T))

m=50             # batch size

### Compute the MNS charting statistic T_S^2

xi1 = rep(0,20)
xi2 = rep(0,20)
Vhat=array(rep(0,80),dim=c(2,2,20))

for(i in 1:20){
   xi1[i] = sum(sign(x[((i-1)*50+1):(i*50),1]))
   xi2[i] = sum(sign(x[((i-1)*50+1):(i*50),2]))
   Vhat[1,1,i]=m
   Vhat[1,2,i]=sum(sign(x[((i-1)*50+1):(i*50),1])*sign(x[((i-1)*50+1):(i*50),2]))
   Vhat[2,1,i]=Vhat[1,2,i]
   Vhat[2,2,i]=Vhat[1,1,i]
}

TS2=rep(0,20)

for(n in 1:20){
   TS2[n] = t(c(xi1[n],xi2[n]))%*%solve(Vhat[,,n])%*%c(xi1[n],xi2[n])
}

### Compute the MNSR charting statistic T_SR^2

psi1 = rep(0,20)
psi2 = rep(0,20)
What=array(rep(0,80),dim=c(2,2,20))

for(i in 1:20){
   psi1[i] = sum(sign(x[((i-1)*50+1):(i*50),1])*rank(abs(x[((i-1)*50+1):(i*50),1])))
   psi2[i] = sum(sign(x[((i-1)*50+1):(i*50),2])*rank(abs(x[((i-1)*50+1):(i*50),2])))
   What[1,1,i]=m*(m+1)*(2*m+1)/6
   What[1,2,i]=sum(sign(x[((i-1)*50+1):(i*50),1])*rank(abs(x[((i-1)*50+1):(i*50),1]))*
                 sign(x[((i-1)*50+1):(i*50),2])*rank(abs(x[((i-1)*50+1):(i*50),2])))
   What[2,1,i]=What[1,2,i]
   What[2,2,i]=What[1,1,i]
}

TSR2=rep(0,20)

for(n in 1:20){
   TSR2[n] = t(c(psi1[n],psi2[n]))%*%solve(What[,,n])%*%c(psi1[n],psi2[n])
}


write.table(round(x,digits=3),
            "example91.dat",row.names=F,col.names=F)


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

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

n = seq(1,20)
nn = rep(n,rep(50,20))

plot(nn,x[,1],type="p",pch=16,xlab="n",ylab=expression(X[nj1]),mgp=c(2,1,0),
     xlim=c(0,21), ylim=c(-8,8), cex=0.8)
title(xlab="(a)",cex=0.9)

plot(nn,x[,2],type="p",pch=16,xlab="n",ylab=expression(X[nj2]),mgp=c(2,1,0),
     xlim=c(0,21), ylim=c(-8,8), cex=0.8)
title(xlab="(b)",cex=0.9)

plot(n,TS2,type="o",lty=1,pch=16,xlab="n",ylab=expression(T[S]^2),mgp=c(2,1,0),
     xlim=c(0,21), ylim=c(0,35),cex=0.8)
lines(n,rep(10.59663,20),lty=2,cex=0.8)
title(xlab="(c)",cex=0.9)

plot(n,TSR2,type="o",lty=1,pch=16,xlab="n",ylab=expression(T[SR]^2),mgp=c(2,1,0),
     xlim=c(0,21), ylim=c(0,35),cex=0.8)
lines(n,rep(10.59663,20),lty=2,cex=0.8)
title(xlab="(d)",cex=0.9)

graphics.off()