##### R-code for Example 8.6. It writes the related data to 
##### the file "example86.dat." It also makes "fig88.ps."

set.seed(10)

#### Generate 100 observations of the reference sample from the
#### IC process distribution t(3)/sqrt(3)

y = rt(100,3)/sqrt(3)

#### Generate 50 observations from the t(3)/sqrt(3) distribution, 
#### and another 50 observations from t(3)/sqrt(3)+1

x1=rt(50,3)/sqrt(3)
x2=rt(50,3)/sqrt(3)+1

x = c(x1,x2)
n = length(x)

lambda=0.05

# Define F_{n,lambda}(Xn)

Fnlambda=rep(0,n) 
  
for(i in 1:n){
   for(j in 1:i){
      Fnlambda[i]=Fnlambda[i]+(1-lambda)^(i-j)*as.numeric(x[j] <= x[i])
   }
   Fnlambda[i]=(Fnlambda[i]-0.5)/(1-(1-lambda)^i)*lambda 
}   # The number 0.5 above is the continuity correction 

# Define Fhat_0(Xn)

Fhat0=rep(0,n)

for(j in 1:100){
   Fhat0[1]=Fhat0[1]+as.numeric(y[j] <= x[1])
}
Fhat0[1]=Fhat0[1]/100

for(i in 2:n){

   for(j in 1:100){
      Fhat0[i]=Fhat0[i]+as.numeric(y[j] <= x[i])
   }

   for(j in 1:(i-1)){
      Fhat0[i]=Fhat0[i]+as.numeric(x[j] <= x[i])
   }

   if(Fhat0[i]<=0){Fhat0[i]=0.5/(100+i-1)}  # We have applied the continuity
   else{Fhat0[i]=(Fhat0[i]-0.5)/(100+i-1)}  # correction here.

}

# Define Gntilde

Gntilde=rep(0,n)
for(i in 1:n){
   Gntilde[i]=(log(Fnlambda[i]/Fhat0[i]))/(1-Fnlambda[i])+
              (log((1-Fnlambda[i])/(1-Fhat0[i])))/Fnlambda[i]
}

# Define the charting statistic EnSS

EnSS=rep(0,n)
EnSS[1]=lambda*Gntilde[1]

for(i in 2:n){
   EnSS[i]=((lambda*Gntilde[i]+(1-lambda)*EnSS[i-1]))
}

write.table(round(cbind(y,x,Fnlambda,Fhat0,Gntilde,EnSS),digits=3), 
            file="example86.dat",col.names=T, row.names=F)

postscript("fig88.ps",width=7,height=3.5,horizontal=F)

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

ii = seq(1,n)
ii1=seq(-99,0)

plot(ii,x,type="p",pch=16,xlab="n",ylab=expression(X[n]),mgp=c(2,1,0),
     xlim=c(-100,100), ylim=c(-2.5,3), cex=0.8)
points(ii1,y,pch=18, cex=0.8)
lines(rep(0,2),c(-2.4,2.9),lty=3,cex=0.8)
title(xlab="(a)",cex=0.9)

hn=c(0.75285, 0.75355, 0.74010, 0.72570, 0.71060, 0.68615, 0.66660, 
     0.64215, 0.62245, 0.60120, 0.57830, 0.55615, 0.53815, 0.52015, 
     0.50305, 0.48535, 0.46980, 0.45475, 0.43895, rep(0.42595,2),
     rep(0.40010,2), rep(0.37620,2), rep(0.35605,2), rep(0.33730,2), 
     rep(0.32050,5), rep(0.28415,5), rep(0.25630,10), rep(0.22020,10), 
     rep(0.19885,10), rep(0.18745,10), rep(0.17990,10), rep(0.17640,11))

plot(ii,EnSS,type="o",lty=1,pch=16,xlab="n",ylab=expression(E[list(n,SS)]),
     mgp=c(2,1,0),xlim=c(0,101), ylim=c(0,0.9),cex=0.8)
lines(ii,hn,lty=2,cex=0.8)
title(xlab="(b)",cex=0.9)

graphics.off()