##### R-code for Example 6.5. It writes the related data to
##### the file "example65.dat." It also creates "fig65.ps."

set.seed(100)

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

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])^2)

### Start the phase II process monitoring

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

hn[10]=10.451
hn[11]=9.840
hn[12]=9.653
hn[13]=9.634
hn[14]=9.658
hn[15]=9.692

for(i in 10:n){
    if (i>15) {
       hn[i] = -1.38-2.241*log(0.005)+(1.61+0.691*log(0.005))/sqrt(i-9)
    }

    Wn[i] = Wn[i-1]+x[i]
    Sn2[i] = Sn2[i-1]+((i-1)*x[i]-Wn[i-1])^2/i/(i-1)

    Xbarjip = rep(0,n)
    Vji2p = rep(0,n)
    Bji = rep(0,n)
    Cji = rep(0,n)

    for(j in 9:(i-2)){
 
        Xbarjip[j] = (Wn[i]-Wn[j])/(i-j)
        Vji2p[j] = Sn2[i]-Sn2[j]-j*(i-j)/i*(Wn[j]/j-Xbarjip[j])^2
        Cji[j]=1+(1/(j-1)+1/(i-j-1)-1/(i-2))/3
        Bji[j]=((j-1)*log(j-1+(i-j-1)*(Vji2p[j]/(i-j-1))/(Sn2[j]/(j-1)))+
                (i-j-1)*log((j-1)*(Sn2[j]/(j-1))/(Vji2p[j]/(i-j-1))+i-j-1)-
                (i-2)*log(i-2))/Cji[j]
    }

    Bmax[i] = max(Bji[9:(i-2)])

}

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

### Make a figure

postscript("fig65.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,4), cex=0.8)
title(xlab="(a)",cex=0.9)

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

graphics.off()