##### R-code for making Figure 3.4 to demonstrate the selef-starting CUSUM
##### chart with ARL0=200 and k=0.25

set.seed(10)

x = c(rnorm(50,0,1),rnorm(50,0.5,1))

N = length(x)

xbar = rep(0,N)
xstd = rep(0,N)
zn   = rep(0,N)
cnss = rep(0,N)
cn   = rep(0,N)

for(i in 1:N){
   xbar[i] = mean(x[1:i])
}

for(i in 2:N){
   xstd[i] = sd(x[1:i])
}

for(i in 3:N){
   zn[i] = qnorm(pt(sqrt((i-1)/i)*(x[i]-xbar[i-1])/xstd[i-1],i-2))
}

k=0.25
h=5.597

cn[1] = max(0,x[1]-k)

for(i in 2:N){
   cn[i] = max(0,cn[i-1]+x[i]-k)
}


for(i in 3:length(x)){
   cnss[i] = max(0,cnss[i-1]+zn[i]-k)
}

ii = seq(1,N)

pdf("fig34.pdf",width=4.5,height=4.5)

plot(ii,cnss,type="o",lty=1,pch=16,xlab="n",
     ylab="",mgp=c(2,1,0),xlim=c(0,N), 
     ylim=c(0,20),cex=0.8)
lines(ii,rep(h,length(x)),lty=2,cex=0.8)
lines(ii,cn,lty=3,cex=0.8)

graphics.off()