##### R-code for Example 10.1 (continued). It makes
##### "fig102.ps"

x = c(0.2,0.4,0.6,0.8,1)
n=length(x)
sxx=var(x)*(n-1)

y = matrix(0,20,5)

set.seed(10)

### generate 10 IC profiles from the model y=2+x+e, where e~N(0,0.5^2)

for(i in 1:10){
   y[i,] = 2+x+rnorm(5,0,0.5)
}

### generate 10 OC profiles from the model y=3+0.5x+e, where e~N(0,0.5^2)

for(i in 11:20){
   y[i,] = 3+0.5*x+rnorm(5,0,0.5)
}

write.table(round(y,digits=3),
            "example101.dat",row.names=F,col.names=F)

### Fit LS line for each profile

a = rep(0,20)
b = rep(0,20)
mse = rep(0,20)

x1=x-mean(x)  ### centered X

for(i in 1:20){
   fit = lm(y[i,] ~ x1)
   a[i] = fit$coefficients[1]
   b[i] = fit$coefficients[2]
   mse[i] = sum(fit$residuals^2)/fit$df.residual
}

postscript("fig102.ps",width=7.5,height=7.5,horizontal=F)

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

jj=seq(1,20)

### Compute the values of the charting statistic E_{j,I}

laI=0.2   ### lambda_I=0.2
c0=2+mean(x)

EjI = rep(0,20)
EjI[1] = laI*a[1]+(1-laI)*c0
for(j in 2:20){
   EjI[j] = laI*a[j]+(1-laI)*EjI[j-1]
}

U=c0+3.0156*0.5*sqrt(laI/((2-laI)*n))
L=c0-3.0156*0.5*sqrt(laI/((2-laI)*n))

plot(jj,EjI,type="p",pch=16,xlab="j",xlim=c(0,21),ylab=expression(E[list(j,I)]),
     ylim=c(2,4),mgp=c(2,1,0),cex=0.9,xaxt="n",yaxt="n")
lines(jj,rep(U,20),lty=2,cex=0.9)
lines(jj,rep(L,20),lty=2,cex=0.9)
axis(1,at=c(0,5,10,15,20),labels=c("0","5","10","15","20"),cex=.7)
axis(2,at=c(2,2.5,3,3.5,4),labels=c("2","2.5","3","3.5","4"),cex=.7)
title(xlab="(a)",cex=0.9)

### Compute the values of the charting statistic E_{j,S}

laS=0.2   ### lambda_S=0.2
b0=1

EjS = rep(0,20)
EjS[1] = laS*b[1]+(1-laS)*b0
for(j in 2:20){
   EjS[j] = laS*b[j]+(1-laS)*EjS[j-1]
}

U=b0+3.0109*0.5*sqrt(laS/((2-laS)*sxx))
L=b0-3.0109*0.5*sqrt(laS/((2-laS)*sxx))

plot(jj,EjS,type="p",pch=16,xlab="j",xlim=c(0,21),ylab=expression(E[list(j,S)]),
     ylim=c(0,2),mgp=c(2,1,0),cex=0.9,xaxt="n",yaxt="n")
lines(jj,rep(U,20),lty=2,cex=0.9)
lines(jj,rep(L,20),lty=2,cex=0.9)
axis(1,at=c(0,5,10,15,20),labels=c("0","5","10","15","20"),cex=.7)
axis(2,at=c(0,0.5,1,1.5,2),labels=c("0","0.5","1","1.5","2"),cex=.7)
title(xlab="(b)",cex=0.9)

### Compute the values of the charting statistic E_{j,V}

laV=0.2   ### lambda_V=0.2

EjV = rep(0,20)
EjV[1] = max(0,laV*log(mse[1]/0.5^2))
for(j in 2:20){
   EjV[j] = max(0,laV*log(mse[j]/0.5^2)+(1-laV)*EjV[j-1])
}

U=1.3723*sqrt(laV/(2-laV))*sqrt(2/(n-2)+2/(n-2)^2+4/(3*(n-2)^3)-16/(15*(n-2)^5))

plot(jj,EjV,type="p",pch=16,xlab="j",xlim=c(0,21),ylab=expression(E[list(j,V)]),
     ylim=c(0,0.5),mgp=c(2,1,0),cex=0.9,xaxt="n",yaxt="n")
lines(jj,rep(U,20),lty=2,cex=0.9)
axis(1,at=c(0,5,10,15,20),labels=c("0","5","10","15","20"),cex=.7)
axis(2,at=c(0,0.1,0.2,0.3,0.4,0.5),labels=c("0","0.1","0.2","0.3","0.4","0.5"),cex=.7)
title(xlab="(c)",cex=0.9)


graphics.off()