##### R-code for Example 10.1. It writes the related data to the file
##### "example101.dat." It also makes "fig101.ps."

x = c(0.2,0.4,0.6,0.8,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)

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

### Compute the values of the charting statistic T_j^2

Tj2 = rep(0,20)
Sigma =  matrix(0,2,2)

for(i in 1:20){

   n=length(x)
   xbar=mean(x)
   sxx=var(x)*(n-1)

   Sigma[1,1]=1/n+xbar^2/sxx
   Sigma[1,2]=-xbar/sxx
   Sigma[2,1]=-xbar/sxx
   Sigma[2,2]=1/sxx
   Sigma=0.5^2*Sigma
   
   Tj2[i]=t(c(a[i],b[i])-c(2,1))%*%solve(Sigma)%*%
           (c(a[i],b[i])-c(2,1))

}

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

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

j=seq(1,20)

plot(j,a,type="p",pch=16,xlab="j",xlim=c(0,21),ylab=expression(hat(a)[j]),
     ylim=c(0,4),mgp=c(2,1,0),cex=0.9,xaxt="n",yaxt="n")
axis(1,at=c(0,5,10,15,20),labels=c("0","5","10","15","20"),cex=.7)
axis(2,at=c(0,1,2,3,4),labels=c("0","1","2","3","4"),cex=.7)
title(xlab="(a)",cex=0.9)

plot(j,b,type="p",pch=16,xlab="j",xlim=c(0,21),ylab=expression(hat(b)[j]),
     ylim=c(-1,2.5),mgp=c(2,1,0),cex=0.9,xaxt="n",yaxt="n")
axis(1,at=c(0,5,10,15,20),labels=c("0","5","10","15","20"),cex=.7)
axis(2,at=c(-1,0,1,2),labels=c("-1","0","1","2"),cex=.7)
title(xlab="(b)",cex=0.9)

plot(j,mse,type="p",pch=16,xlab="j",xlim=c(0,21),ylab=expression(hat(sigma)[j]^2),
     ylim=c(0,1),mgp=c(2,1,0),cex=0.9,xaxt="n",yaxt="n")
axis(1,at=c(0,5,10,15,20),labels=c("0","5","10","15","20"),cex=.7)
axis(2,at=c(0,0.25,0.5,0.75,1),labels=c("0","0.25","0.5","0.75","1"),cex=.7)
title(xlab="(c)",cex=0.9)

h=qchisq(0.995,2)  ### control limit

plot(j,Tj2,type="o",lty=1,pch=16,xlab="n",ylab=expression(T[j]^2),
     mgp=c(2,1,0),xlim=c(0,21), ylim=c(0,22),cex=0.9,xaxt="n",yaxt="n")
lines(j,rep(h,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,5,10,15,20),labels=c("0","5","10","15","20"),cex=.7)
title(xlab="(d)",cex=0.9)

graphics.off()