# Program 20
# This follows on immediately from Program 19 (i.e., do NOT quit your
# R run before running Program_20), and produces the standardised variance, 
# calculates it at the design's support points, and plots it against the value of x

xvec <- c(-1.5529,1.4004)
deltavec <- c(0.5,0.5)
invinfomatrices <- rep(0,p*p*h)
dim(invinfomatrices) <- c(p,p,h)
for (m in 1:h)
{
 betavec <- betamat[m,]
 invinfomatrices[,,m] <- solve(infomat(xvec,deltavec,betavec))
}
invinfomatrices

stdvar <- function(x)
{
 sv <- 0
 fx <- f(x)
 for (m in 1:h)
 {
  betavec <- betamat[m,]
  eta <- sum(f(x)*betavec)
  expeta <- exp(eta)
  wt <- expeta/((1+expeta)^2)
  matinv <- invinfomatrices[,,m]
  sv <- sv + psivec[m]*wt*t(fx)%*%matinv%*%fx
 }
sv
}

stdvar(-1.5529)
stdvar(1.4004)
xvec <- seq(from=-5,to=5,by=0.01)
yvec <- sapply(xvec,stdvar)
par(las=1)
plot(xvec,yvec,ty="l",xlab="x",ylab="Standardised variance")
lines(c(-5,5),c(2,2),lty=2)