##### R-code for Example 8.4. It writes the related data to the file
##### "example84.dat" and makes "fig85.ps"

set.seed(100)

#### Generate 100 observations of the reference sample from the
#### IC process distribution (chisq_3-3)/sqrt(6)

y = (rchisq(100,3)-3)/sqrt(6)

#### Generate 20 batches of data with the batch size 5 from the
#### (chisq_3-3)/sqrt(6) distribution, and another 10 batches of data from 
#### (chisq_3-3)/sqrt(6)+1

set.seed(10)

x1=matrix((rchisq(100,3)-3)/sqrt(6),nrow=20,ncol=5)
x2=matrix((rchisq(50,3)-3)/sqrt(6)+1,nrow=10,ncol=5)

x = rbind(x1,x2)
n = length(x[,1])

# Define the Wilcoxon rank-sum statistic Wn

Wn=rep(0,n)   
for(i in 1:n){
   Wn[i]=sum(rank(c(x[i,],y))[1:5])
}

# Define the charting statistic Cn+

muWn=5*(5+100+1)/2
k=0.5*sqrt(5*100*(5+100+1)/12)
h=353

Cnp=rep(0,n)
Cnn=rep(0,n)
Cnp[1]=max(0,Wn[1]-muWn-k)
Cnn[1]=min(0,Wn[1]-muWn+k)

for(i in 2:n){
   Cnp[i]=max(0,Cnp[i-1]+Wn[i]-muWn-k)
   Cnn[i]=min(0,Cnn[i-1]+Wn[i]-muWn+k)
}

write.table(round(cbind(x,Wn,Cnp,Cnn),digits=3), 
            file="example84.dat",col.names=T, row.names=F)

postscript("fig85.ps",width=7,height=3.5,horizontal=F)

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

xx=c(x)
ii = rep(seq(1,n),5)

ii1=rep(seq(-20,-1),5)

plot(ii,xx,type="p",pch=16,xlab="n",ylab=expression(X[nj]),mgp=c(2,1,0),
     xlim=c(-20,31), ylim=c(-1.5,5), cex=0.8)
points(ii1,y,pch=18, cex=0.8)
lines(rep(0,2),c(-1.5,5),lty=3,cex=0.8)
title(xlab="(a)",cex=0.9)

i1=seq(1,n)

plot(i1,Cnp,type="o",lty=1,pch=16,xlab="n",ylab="",mgp=c(2,1,0),xlim=c(0,31), 
     ylim=c(-400,1400),cex=0.8)
lines(i1,Cnn,type="o",lty=3,pch=18,cex=0.8)
lines(i1,rep(h,n),lty=2,cex=0.8)
lines(i1,rep(-h,n),lty=2,cex=0.8)
title(xlab="(b)",cex=0.9)

graphics.off()