##### R-code for making Figure 8.7

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 En

muWn=5*(5+100+1)/2
lambda=0.1

L=219.5
U=310.5

En=rep(0,n)

En[1]=lambda*Wn[1]+(1-lambda)*muWn

for(i in 2:n){
   En[i]=lambda*Wn[i]+(1-lambda)*En[i-1]
}

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

postscript("fig87.ps",width=4,height=4,horizontal=F)

i1=seq(1,n)

plot(i1,En,type="o",lty=1,pch=16,xlab="n",ylab=expression(E[n]),
     mgp=c(2,1,0),xlim=c(0,31), ylim=c(200,400),cex=0.8)
lines(i1,rep(L,n),lty=2,cex=0.8)
lines(i1,rep(U,n),lty=2,cex=0.8)

graphics.off()