##### R-code for Example 9.2. It writes the related data to the file
##### "example92.dat." It also makes "fig94.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=10,sigma=Sigma,df=5),rmvt(n=10,sigma=Sigma,df=5))
x = x+rbind(matrix(rep(mu1,10),10,2,byrow=T),matrix(rep(mu2,10),10,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,...,20

rXn = rep(0,20)

for(i in 1:20){

   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 Snstar(F_0M), n=1,2,...,20

SnF0M = rep(0,20)
Snstar = rep(0,20)

SnF0M[1] = rXn[1]-0.5
Snstar[1] = SnF0M[1]/(1*sqrt((1/M+1/1)/12))

for(i in 2:20){
   SnF0M[i] = SnF0M[i-1] + (rXn[i]-0.5)
   Snstar[i] = SnF0M[i]/(i*sqrt((1/M+1/i)/12))
}

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


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

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

n = seq(1,20)
n1= seq(-99,0)

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

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

plot(n,rXn,type="o",lty=1,pch=16,xlab="n",ylab=expression(r(bold(X)[n]*";"*F[0*M])),mgp=c(2,1,0),
     xlim=c(0,21), ylim=c(0,1),cex=0.8)
lines(n,rep(0.05,20),lty=2,cex=0.8)
title(xlab="(c)",cex=0.9)

plot(n,Snstar,type="o",lty=1,pch=16,xlab="n",ylab=expression(S[n]^{"*"}*(F[0*M])),mgp=c(2,1,0),
     xlim=c(0,21), ylim=c(-2.5,1.5),cex=0.8)
lines(n,rep(-1.95,20),lty=2,cex=0.8)
title(xlab="(d)",cex=0.9)

graphics.off()