##### R-code for Example 8.5. It writes the related data to the file
##### "example85.dat" and makes "fig86.ps"

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

set.seed(10)

x1=matrix(rt(200,3)/sqrt(3),nrow=20,ncol=10)
x2=matrix(rt(100,3)/sqrt(3)+1,nrow=10,ncol=10)

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

# Define the Wilcoxon singed-rank statistic psi_i

psi=rep(0,n)   
for(i in 1:n){
   psi[i]=sum(sign(x[i,])*rank(abs(x[i,])))
}

# Define the charting statistic Cn+

lambda=0.1
h=2.684*sqrt((10*11*21/6)*(lambda/(2-lambda)))

En=rep(0,n)
En[1]=lambda*psi[1]

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

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

postscript("fig86.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),10)

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

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(-12.1,31),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()