##### R-code for Example 8.2. It writes the related data to 
##### the file "example82.dat" and makes "fig83.ps"

set.seed(100)

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

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

L=y[order(y)[53]]
U=y[order(y)[948]]

#### 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)+2

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

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

# Define the charting statistic psi_i

Xmedian=rep(0,n)   
for(i in 1:n){
   Xmedian[i]=median(x[i,])
}

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

postscript("fig83.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)

plot(ii,xx,type="p",pch=16,xlab="i",ylab=expression(X[ij]),mgp=c(2,1,0),
     xlim=c(0,31), ylim=c(-1.2,10), cex=0.8)
title(xlab="(a)",cex=0.9)

i1=seq(1,n)

plot(i1,Xmedian,type="o",lty=1,pch=16,xlab="i",
     ylab=expression(X[i(3)]),mgp=c(2,1,0),xlim=c(0,31), 
     ylim=c(-1.1,3),cex=0.8)
lines(i1,rep(U,n),lty=2,cex=0.8)
lines(i1,rep(L,n),lty=2,cex=0.8)
title(xlab="(b)",cex=0.9)

graphics.off()