##### R-code for Example 9.3. It writes the related data to the file
##### "example93.dat." It also makes "fig95.ps."

#### Generate 100 2-D random vectors from t_5(mu,Sigma)

library(mvtnorm) # The R library in which the commands for generating 
                 # multivariate t or normal random vectors is contained.

set.seed(10)

mu1 = c(0,0)     # IC mean vector
mu2 = c(2,0)     # OC mean vector
Sigma = matrix(c(1,0.5,0.5,1),2,2)  # IC covariance matrix

M=100
y = rmvt(n=M,sigma=Sigma,df=5)
Ybar=colMeans(y)
S2Y=var(y)

x = rbind(rmvt(n=100,sigma=Sigma,df=5),rmvt(n=100,sigma=Sigma,df=5))
x = x+rbind(matrix(rep(mu1,100),100,2,byrow=T),matrix(rep(mu2,100),100,2,byrow=T))

### Compute D_M(Y_i,F_0M), for i=1,2,...,M

DMY = rep(0,M)

for(i in 1:M){
   DMY[i] = 1/(1+t(y[i,]-Ybar)%*%solve(S2Y)%*%(y[i,]-Ybar))
}

### Compute r(Xn,F_0M), n=1,2,...,200

rXn = rep(0,200)

for(i in 1:200){

   temp = 1/(1+t(x[i,]-Ybar)%*%solve(S2Y)%*%(x[i,]-Ybar)) # D_M(X_i,F_0M)

   for(j in 1:M){
      if(DMY[j] <= temp) rXn[i]=rXn[i]+1
   }

   rXn[i]=rXn[i]/M

}

### Compute Q(F_{0M},F_{1m,n}), n=1,2,...,20

m=10    # batch size

QFFn=rep(0,20)

for(i in 1:20){
   for(j in ((i-1)*m+1):(i*m)){
      QFFn[i]=QFFn[i]+rXn[j]
   }
   QFFn[i]=QFFn[i]/m
}

write.table(round(cbind(x,rXn),digits=3),
            "example93.dat",row.names=F)


postscript("fig95.ps",width=6.5,height=6.5,horizontal=F)

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

n = seq(1,20)
nn = rep(seq(1,20),rep(10,20))
n1= seq(-99,0)/5

plot(nn,x[,1],type="p",pch=16,xlab="n",ylab=expression(X[n1]),mgp=c(2,1,0),
     xlim=c(-20,21), ylim=c(-4,11), cex=0.8)
points(n1,y[,1],pch=18, cex=0.8)
lines(rep(0,2),c(-3.9,10.9),lty=3,cex=0.8)
title(xlab="(a)",cex=0.9)

plot(nn,x[,2],type="p",pch=16,xlab="n",ylab=expression(X[n2]),mgp=c(2,1,0),
     xlim=c(-20,21), ylim=c(-4,11), cex=0.8)
points(n1,y[,2],pch=18, cex=0.8)
lines(rep(0,2),c(-3.9,10.9),lty=3,cex=0.8)
title(xlab="(b)",cex=0.9)

h=0.5-1.96*sqrt((1/100+1/10)/12)

plot(n,QFFn,type="o",lty=1,pch=16,xlab="n",ylab=expression(Q(F[0*M]*","*F[1*m*","*n])),
     mgp=c(2,1,0),xlim=c(0,21), ylim=c(0,1),cex=0.8)
lines(n,rep(h,20),lty=2,cex=0.8)
title(xlab="(c)",cex=0.9)

graphics.off()