##### R-code for Example 10.3 (continued). It also makes
##### "fig105.ps"

x = seq(0.1,1,0.15)
n = length(x)
y = matrix(0,20,n)

set.seed(100)

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

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

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

for(i in 11:20){
   y[i,] = 1/(1+x^2)+rnorm(n,0,0.1)
}

### Fit NLS function for each profile

d = rep(0,20)
b = rep(0,20)

for(i in 1:20){
   yy=y[i,]
   fit = nls(yy ~ d1/(1+x^b1),start=list(d1=1,b1=1))
   d[i] = coef(fit)[1]
   b[i] = coef(fit)[2]
}

thetabar=c(mean(b),mean(d))
Stheta=var(cbind(b,d))

### Compute the values of the Shewhart charting statistic T_j^2 in (10.31)

Tj2 = rep(0,20)

for(i in 1:20){
   
   Tj2[i]=t(c(b[i],d[i])-thetabar)%*%solve(Stheta)%*%
           (c(b[i],d[i])-thetabar)

}

postscript("fig105.ps",width=4.5,height=4.5,horizontal=F)

j=seq(1,20)

h=(20-1)^2/20*qbeta(0.995,1,8.5)  ### control limit of the Shewhart chart

plot(j,Tj2,type="o",lty=1,pch=16,xlab="j",ylab=expression(tilde(T)[j]^2),
     mgp=c(2,1,0),xlim=c(0,21), ylim=c(0,10),cex=0.8,xaxt="n",yaxt="n")
lines(j,rep(h,20),lty=2,cex=0.8)
axis(1,at=c(0,5,10,15,20),labels=c("0","5","10","15","20"),cex=.6)
axis(2,at=c(0,2,4,6,8,10),labels=c("0","2","4","6","8","10"),cex=.6)

graphics.off()