##### R-code for making Figure 6.6

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

par(mfcol=c(4,3),mar=c(4,4,1,2))


### First data with mean shift only

set.seed(50)

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

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

ii=seq(1,n)

plot(ii,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)

##### Construct the CPD chart for detecting mean shifts

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

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

}

plot(ii[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(ii[10:30],hn[10:30],lty=2,cex=0.8)
title(xlab="(b)",cex=0.9)


##### Construct the CPD chart for detecting variance shifts

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)

Wn[10] = Wn[9]+x[10]
Sn2[10] = Sn2[9]+(9*x[10]-Wn[9])^2/10/9

### 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 11: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)])

}

plot(ii[11:30],Bmax[11: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,12),cex=0.8)
lines(ii[11:30],hn[11:30],lty=2,cex=0.8)
title(xlab="(c)",cex=0.9)


##### Construct the CPD chart for detecting both mean and variance shifts

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)

Wn[10] = Wn[9]+x[10]
Sn2[10] = Sn2[9]+(9*x[10]-Wn[9])^2/10/9

### Start the phase II process monitoring

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

hn[10]=15.330
hn[11]=14.556
hn[12]=14.313
hn[13]=14.265
hn[14]=14.249

for(i in 11:n){

    if (i>14) {
       hn[i] = 1.58-2.52*log(0.005)+(0.094+0.33*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)
    Jji = 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+11/12*(1/j+1/(i-j)-1/i)+(1/j^2+1/(i-j)^2-1/i^2)
        Jji[j]=(j*log((Sn2[i]/(i-1))/(Sn2[j]/(j-1)))+
                (i-j)*log((Sn2[i]/(i-1))/(Vji2p[j]/(i-j-1))))/Cji[j]
    }

    Jmax[i] = max(Jji[9:(i-2)])

}

plot(ii[11:30],Jmax[11:30],type="l",lty=1,xlab="n",
     ylab=expression(J[list(max,n)]),mgp=c(2,1,0),xlim=c(10,n), 
     ylim=c(0,25),cex=0.8)
lines(ii[11:30],hn[11:30],lty=2,cex=0.8)
title(xlab="(d)",cex=0.9)


### Second data with variance shift only

set.seed(100)

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

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

ii=seq(1,n)

plot(ii,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(-4,4), cex=0.8)
title(xlab="(e)",cex=0.9)

##### Construct the CPD chart for detecting mean shifts

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

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

}

plot(ii[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,6),cex=0.8)
lines(ii[10:30],hn[10:30],lty=2,cex=0.8)
title(xlab="(f)",cex=0.9)


##### Construct the CPD chart for detecting variance shifts

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)

Wn[10] = Wn[9]+x[10]
Sn2[10] = Sn2[9]+(9*x[10]-Wn[9])^2/10/9

### 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 11: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)])

}

plot(ii[11:30],Bmax[11: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,25),cex=0.8)
lines(ii[11:30],hn[11:30],lty=2,cex=0.8)
title(xlab="(g)",cex=0.9)


##### Construct the CPD chart for detecting both mean and variance shifts

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)

Wn[10] = Wn[9]+x[10]
Sn2[10] = Sn2[9]+(9*x[10]-Wn[9])^2/10/9

### Start the phase II process monitoring

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

hn[10]=15.330
hn[11]=14.556
hn[12]=14.313
hn[13]=14.265
hn[14]=14.249

for(i in 11:n){

    if (i>14) {
       hn[i] = 1.58-2.52*log(0.005)+(0.094+0.33*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)
    Jji = 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+11/12*(1/j+1/(i-j)-1/i)+(1/j^2+1/(i-j)^2-1/i^2)
        Jji[j]=(j*log((Sn2[i]/(i-1))/(Sn2[j]/(j-1)))+
                (i-j)*log((Sn2[i]/(i-1))/(Vji2p[j]/(i-j-1))))/Cji[j]
    }

    Jmax[i] = max(Jji[9:(i-2)])

}

plot(ii[11:30],Jmax[11:30],type="l",lty=1,xlab="n",
     ylab=expression(J[list(max,n)]),mgp=c(2,1,0),xlim=c(10,n), 
     ylim=c(0,25),cex=0.8)
lines(ii[11:30],hn[11:30],lty=2,cex=0.8)
title(xlab="(h)",cex=0.9)


### Third data with both mean and variance shifts

set.seed(50)

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

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

ii=seq(1,n)

plot(ii,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(-4,6.5), cex=0.8)
title(xlab="(i)",cex=0.9)

##### Construct the CPD chart for detecting mean shifts

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

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

}

plot(ii[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(ii[10:30],hn[10:30],lty=2,cex=0.8)
title(xlab="(j)",cex=0.9)


##### Construct the CPD chart for detecting variance shifts

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)

Wn[10] = Wn[9]+x[10]
Sn2[10] = Sn2[9]+(9*x[10]-Wn[9])^2/10/9

### 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 11: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)])

}

plot(ii[11:30],Bmax[11: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,20),cex=0.8)
lines(ii[11:30],hn[11:30],lty=2,cex=0.8)
title(xlab="(k)",cex=0.9)


##### Construct the CPD chart for detecting both mean and variance shifts

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)

Wn[10] = Wn[9]+x[10]
Sn2[10] = Sn2[9]+(9*x[10]-Wn[9])^2/10/9

### Start the phase II process monitoring

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

hn[10]=15.330
hn[11]=14.556
hn[12]=14.313
hn[13]=14.265
hn[14]=14.249

for(i in 11:n){

    if (i>14) {
       hn[i] = 1.58-2.52*log(0.005)+(0.094+0.33*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)
    Jji = 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+11/12*(1/j+1/(i-j)-1/i)+(1/j^2+1/(i-j)^2-1/i^2)
        Jji[j]=(j*log((Sn2[i]/(i-1))/(Sn2[j]/(j-1)))+
                (i-j)*log((Sn2[i]/(i-1))/(Vji2p[j]/(i-j-1))))/Cji[j]
    }

    Jmax[i] = max(Jji[9:(i-2)])

}

plot(ii[11:30],Jmax[11:30],type="l",lty=1,xlab="n",
     ylab=expression(J[list(max,n)]),mgp=c(2,1,0),xlim=c(10,n), 
     ylim=c(0,25),cex=0.8)
lines(ii[11:30],hn[11:30],lty=2,cex=0.8)
title(xlab="(l)",cex=0.9)


graphics.off()