##### R-code for Example 8.3. It writes the related data to file
##### "example83.dat" and makes "fig84.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(120,3)/sqrt(3),nrow=20,ncol=6)
x2=matrix(rt(60,3)/sqrt(3)+1,nrow=10,ncol=6)

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+

k=8
h=10

Cnp=rep(0,n)
Cnp[1]=max(0,psi[1]-k)

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

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

postscript("fig84.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),6)

plot(ii,xx,type="p",pch=16,xlab="n",ylab=expression(X[nj]),mgp=c(2,1,0),
     xlim=c(0,31), ylim=c(-3,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=expression(C[n]^{"+"}),mgp=c(2,1,0),xlim=c(0,31), 
     ylim=c(0,120),cex=0.8)
lines(i1,rep(h,n),lty=2,cex=0.8)
title(xlab="(b)",cex=0.9)

graphics.off()