# Program 13

# This calculates the standardised variance at points selected randomly in the neighbourhood 
# of each support point of the design.  The neighbourhood extends a distance " +/- distance"  
# from each support point in directions parallel to each coordinate axis.
#
# A program that generates the standardised variance function (sv) needed in Program 13 appears
# below the row of hash (#) signs below.

out4 <- matrix(c(1,-1,-1,-1,1,-1),s,m)

npoints <- 1000
total <- npoints*m
distance <- 0.04
for (i in 1:s)
{
 maxpoint <- rep(0,m)
 supportpt <- out4[i,1:m]
 max <- -1
 deviation <- distance*(runif(total)-0.5)
 j1  <- 1
 j2 <- m
 for (j in 1:npoints)
 {
  newpoint <- supportpt + deviation[j1:j2]
  newpoint[newpoint < -1]  <- -1
  newpoint[newpoint > 1] <- 1
  sv <- stdvar(newpoint)
  if(sv > max) {max <- sv
    maxpoint <- newpoint}
  j1 <- j1 + m
  j2 <- j2 + m
 }
 cat("For support point ",i,"\n", "Maximum std var is ",max,
  " at \n",maxpoint,"\n")
}

#########################
#Program that generates the output needed by Program 13 above.

m <- 2
betavec <- c(1,1,1)
p <- 3
s <- 3
lim1 <- m*s

fx <- function(xvec)
{
 fvec <- c(1,xvec)
 fvec
}
design <- c(1,-1,-1,-1,1,-1,1/3,1/3,1/3)

 xmat <- matrix(design[1:lim1],m,s,byrow=T)
 wtvals <- (design[(lim1+1):((m+1)*s)])^2
 deltavec <- wtvals/sum(wtvals)
 fxmat <- apply(xmat,2,fx)
 etavec <- as.vector(t(betavec)%*%fxmat)
 expetavec <- exp(etavec)
 modelwtvec <- expetavec/((1+expetavec)^2)
 infomat <- fxmat %*%diag(deltavec*modelwtvec)%*%t(fxmat)
 invinfo <- solve(infomat)

stdvar <- function(xvec)
{
 fvec <- fx(xvec)
 eta <- sum(fvec*betavec)
 expeta <- exp(eta)
 wt <- expeta/(1+expeta)^2
 sv <- wt*t(fvec)%*%invinfo%*%fvec
 sv
}