# Program 12
# Calculates the standardised variance for the alleged D_s-optimal design when
# choosing between a quadratic and a linear regression
# As suggested in Program 11, a loop has been eliminated in the calculation of the 
# standardised variance for plotting purposes

#Look at quadratic model
p <- 3
s <- 3
f <- function(x)
{ c(1,x, x^2)}
f2 <- function(x)
{ c(1,x)}

infomat <- matrix(0,p,p)
xvals <- c(-1,0,1)
deltavec <- c(0.25,0.5,0.25)
for (i in 1:s)
{
 xvec <- xvals[i]
 fvec <- f(xvec)
 infomat <- infomat + deltavec[i]*fvec%*%t(fvec)
}

invinfo1 <- solve(infomat)
stdvar1 <- function(x)
{
 t(f(x))%*%invinfo1%*%f(x)
}

#Now look at linear model
p <- 2
 infomat <- matrix(0,p,p)
 xvals <- c(-1,0,1)
 deltavec <- c(0.25,0.5,0.25)
 for (i in 1:s)
 {
  xvec <- xvals[i]
  fvec <- f2(xvec)
  infomat <- infomat + deltavec[i]*fvec%*%t(fvec)
 }

invinfo2 <- solve(infomat)
stdvar2 <- function(x)
{
 t(f2(x))%*%invinfo2%*%f2(x)
}

# Standardised variance for Ds_optimality
stdvar <- function(x)
{
 stdvar1(x) - stdvar2(x)
}

#Calculate std var at support points
stdvar(-1)
stdvar(0)
stdvar(1)

#Plot stdvar over design space
x <- seq(from=-1,to=1,by=0.01)
y <- sapply(x,stdvar)
par(las=1)
plot(x,y,ty="l",xlab="z",ylab="Standardised Variance",lwd=2)#,xaxt="n")
lines(c(-1,1),c(1,1),lty=2,lwd=2)
lines(c(-1,-1),c(-0.2,1),lty=2,lwd=2)
lines(c(0,0),c(-0.2,1),lty=2,lwd=2)
lines(c(1,1),c(-0.2,1),lty=2,lwd=2)
points(c(-1,0,1),c(1,1,1),pch=16,cex=2)