##### R-code for Example 10.3. It writes the related data to the file
##### "example103.dat." It also makes "fig104.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)
}

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

### Fit NLS function for each profile

d = rep(0,20)
b = rep(0,20)
mse = array(data=rep(0,80),dim=c(20,2,2))
G = matrix(0,n,2)

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]
   G = cbind(1/(1+x^b[i]),-d[i]*x^b[i]*log(x)/(1+x^b[i])^2)
   mse[i,,] = 0.25^2*solve(t(G)%*%G)
}

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

Tj2 = rep(0,20)

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

}

### Compute the values of the EWMA charting statistic E_j

lambda=0.2
Ej = matrix(0,20,2)
Sigma = array(data=rep(0,80),dim=c(20,2,2))

Ej[1,]=lambda*t(c(d[1],b[1])-c(2,1))
Sigma[1,,]=lambda^2*mse[1,,]
Wj2 = rep(0,20)
Wj2[1] = t(Ej[1,])%*%solve(Sigma[1,,])%*%Ej[1,]

for(i in 2:20){
   
   Ej[i,]=lambda*t(c(d[i],b[i])-c(2,1))+(1-lambda)*Ej[i-1,]
   Sigma[i,,]=lambda^2*mse[i,,]+(1-lambda^2)*Sigma[i-1,,]
   Wj2[i] = t(Ej[i,])%*%solve(Sigma[i,,])%*%Ej[i,]

}


postscript("fig104.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,b,type="p",pch=16,xlab="j",xlim=c(0,21),ylab=expression(hat(b)[j]),
     ylim=c(0,5.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(0,1,2,3,4,5),labels=c("0","1","2","3","4","5"),cex=.7)
title(xlab="(a)",cex=0.9)

plot(j,d,type="p",pch=16,xlab="j",xlim=c(0,21),ylab=expression(hat(d)[j]),
     ylim=c(0,3),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),labels=c("0","1","2","3"),cex=.7)
title(xlab="(b)",cex=0.9)

h=qchisq(0.995,2)  ### control limit of the Shewhart chart

plot(j,Tj2,type="o",lty=1,pch=16,xlab="j",ylab=expression(T[j]^2),
     mgp=c(2,1,0),xlim=c(0,21), ylim=c(0,105),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,25,50,75,100),labels=c("0","25","50","75","100"),cex=.7)
title(xlab="(c)",cex=0.9)

h=qchisq(0.995,2)  ### control limit of the Shewhart chart

plot(j,Wj2,type="o",lty=1,pch=16,xlab="j",ylab=expression(W[j]^2),
     mgp=c(2,1,0),xlim=c(0,21), ylim=c(0,25),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,25),labels=c("0","5","10","15","20","25"),cex=.7)
title(xlab="(d)",cex=0.9)

graphics.off()