# Program 11
# This calculates the standardised variance function of the hypothesised
# D-optimal design, then evaluates the std var at the two support points,
# and draws a graph of the function vs x

# An option is presented below to eliminate a loop and make the program run faster

beta0 <- 0
beta1 <- 1
s <- 2
infomat <- function(x)
{
 info <- matrix(0,2,2)
 for (i in 1:s)
{
  pt <- x[i]
  delta <- x[i+s]
  expeta <- exp(beta0 + beta1*pt)
  wt <- expeta/(1+expeta)^2
  info <- info + delta*wt*matrix(c(1,pt,pt,pt^2),2,2)
 }
 info
}

optdesign <- c(-1.5434,1.5434,0.5,0.5)
optmat <- infomat(optdesign)
invmat <- solve(optmat)

stdvar <- function(x)
{
 expeta <- exp(x)
 wt <- expeta/(1+expeta)^2
 fx <- matrix(c(1,x),2,1)
 sv <- wt*t(fx)%*%invmat%*%fx
 sv
}

stdvar(-1.5434)
stdvar(1.5434)

x <- seq(from=-10,to=10,by=0.02)
lx <- length(x)                           #
y <- rep(0,lx)                            #  These six commands may be replaced by the one
for (i in 1:lx)                           #  command given at the bottom of this pages
{                                         #
 y[i] <- stdvar(x[i])                     #
}                                         #
par(las=1)
plot(x,y,ty="l",ylab="Standardised Variance",lwd=2,xaxt="n")
axis(1,at=c(-10,-5,-1.5434,1.5434,5,10),
   label=c("-10","-5","-1.5434","1.5434","5","10"))
lines(c(-10,10),c(2,2),lty=2,lwd=2)
lines(c(-1.5434,-1.5434),c(-0.2,2),lty=2,lwd=2)
lines(c(1.5434,1.5434),c(-0.2,2),lty=2,lwd=2)
points(c(-1.5434,1.5434),c(2,2),pch=16,cex=2)

# If you wish, you may replace the six commands indicated above by 
# the following one command:
y <- sapply(x,stdvar)