#Program 25

#Calculates the standardised variance for a suspected Bayesian D-optimal design.
#Then evaluates it at the design's support points, and produces a contour plot.

#The design used here is Design 7.2 in Chapter 7.

infomat <- function(betavec,xvec,deswts)
{
 xmat <- t(matrix(xvec,s,m))
 fxmat <- apply(xmat,2,fx)
 etavec <- as.vector(t(betavec)%*%fxmat)
 expetavec <- exp(etavec)
 modelwtvec <- expetavec/((1+expetavec)^2)
 infomat <- fxmat %*% diag(deswts*modelwtvec) %*% t(fxmat)
 infomat
}

invinfomatrices <- rep(0,p*p*npts)
dim(invinfomatrices) <- c(p,p,npts)

optxvec <- c(-0.042,1,-1,-1,0.052,1,1,-1,1,-1,0.078,-1,-1,-0.057,1,1)
deltavec <- c(0.074,0.180,0.071,0.183,0.071,0.067,0.183,0.171)

for (h in 1:npts)
{
 betavec <- abscissae[h,]
 invinfomatrices[,,h] <- solve(infomat(betavec,optxvec,deltavec))
}

stdvar <- function(xvec)
{
 sum <- 0
 for (h in 1:npts)
 {
  fvec <- fx(xvec)
  betavec <- abscissae[h,]
  expeta <- exp(sum(betavec*fvec))
  deswt <- expeta/(1+expeta)^2
  temp <- wts[h]*deswt*t(fvec)%*%invinfomatrices[,,h]%*%fvec
  sum <- sum + temp
 }
 sum
}

#Evaluate stdvar at each support point
xmat <- matrix(optxvec,m,s,byrow=T)
for (i in 1:s)
{
 point <- xmat[,i]
 sv <- stdvar(point)
 cat("support point = ",point,"  stdvar = ",sv,"\n")
}

#Plot the standardised variance
exp1 <- expression(x[1])
exp2 <- expression(x[2])
x1 <- x2 <- seq(from=-1.1,to=1.1,by=0.05)
lx1 <- length(x1)
y <- rep(0,lx1^2)
kount <- 1
for (i in 1:lx1)
{
 for (j in 1:lx1)
 {
  y[kount] <- stdvar(c(x1[i],x2[j]))
  kount <- kount + 1
 }
}
ymat <- matrix(y,lx1,lx1,byrow=T)
par(las=1)
contour(x1,x2,ymat,xlab=exp1,ylab=exp2, levels = seq(from=2.2,to=3.2,by=0.2),lwd = c(1,1,1,1,2,1),
   labcex=1)
lines(c(-1.05,1.05),c(1,1),lty=2)
lines(c(-1.05,1.05),c(-1,-1),lty=2)
lines(c(-1,-1),c(-1.05,1.05),lty=2)
lines(c(1,1),c(-1.05,1.05),lty=2)
points(c(-0.042,1,-1,-1,0.052,1,1,-1),c(1,-1,0.078,-1,-1,-0.057,1,1),pch=16,cex=2)