##### R-code for Example 6.4. It writes the related data to
##### the file "example64.dat." It also creates "fig63.ps."

set.seed(50)

x1 = rnorm(20,0,1)
x2 = rnorm(10,2,1)

x =c(x1,x2)
n=length(x)

Wn = rep(0,n)
Sn2 = rep(0,n)

### Compute the initial values of Wn and Sn^2 from the 9 IC observations

Wn[9] = sum(x[1:9])
Sn2[9] = sum((x[1:9]-Wn[9]/9)^2)

### Start the phase II process monitoring

Tmax = rep(0,n)
hn   = rep(0,n)

for(i in 10:n){

    hn[i] = 5.511*(0.677+0.019*log(0.005)+(1-0.115*log(0.005))/(i-6))
    Wn[i] = Wn[i-1]+x[i]
    Sn2[i] = Sn2[i-1]+((i-1)*x[i]-Wn[i-1])^2/i/(i-1)

    Tji2 = rep(0,n)
    Vji2 = rep(0,n)

    for(j in 9:(i-1)){
 
        Vji2[j] = (i*Wn[j]-j*Wn[i])^2/i/j/(i-j)
        Tji2[j] = (i-2)*Vji2[j]/(Sn2[i]-Vji2[j])

    }

    Tmax[i] = max(sqrt(Tji2[9:(i-1)]))

}

write.table(round(cbind(x,Wn,Sn2,Tmax,hn),digits=3), 
            file="example64.dat",col.names=T, row.names=F)

### Make a figure

postscript("fig63.ps",width=7.2,height=3.6,horizontal=F)

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

i=seq(1,n)

plot(i,x,type="o",lty=1,pch=16,xlab="n",
     ylab=expression(X[n]),mgp=c(2,1,0),xlim=c(0,length(x)), 
     ylim=c(-2,5), cex=0.8)
title(xlab="(a)",cex=0.9)

plot(i[10:30],Tmax[10:30],type="l",lty=1,xlab="n",
     ylab=expression(T[list(max,n)]),mgp=c(2,1,0),xlim=c(10,n), 
     ylim=c(0,8),cex=0.8)
lines(i[10:30],hn[10:30],lty=2,cex=0.8)
title(xlab="(b)",cex=0.9)

graphics.off()