# Program 9
# This considers a standard quadratic regression model, and tries to find the D-optimal design
# The value of thed eterminant does not depend upon the values of the parameters beta0, beta1
# or beta2
#
f <- function(x)
{ c(1,x, x^2)}
p <- 3
s <- 3
twos <- 2*s

detinfomat <- function(variables)
{
 infomat <- matrix(0,p,p)
 xvals <- cos(pi*variables[1:s])
 deltavec <- (variables[(s+1):twos])^2
 deltavec <- deltavec/sum(deltavec)
 for (i in 1:s)
 {
  xvec <- xvals[i]
  fvec <- f(xvec)
  wt <- 1
  infomat <- infomat + deltavec[i]*wt*fvec%*%t(fvec)
 }
 -det(infomat)
}

nsims <- 1000
mindet <- 10
#set.seed(123)
for (i in 1:nsims)
{
 initial <- runif(2*s)
 out <- optim(initial,detinfomat,NULL,method="Nelder-Mead")
 valuenow <- out$val
 if(valuenow < mindet){mindet <- valuenow
 design <- out$par}
}
cat("Min value of det\n",mindet,"\n") 
output <- design
out1 <- cos(pi*output[1:s])
out2 <- (output[(s+1):twos])^2
wts <- out2/sum(out2)
out4 <- cbind(matrix(out1,s,1),wts)
out4 <- out4[order(out4[,1]),]
cat("Design\n")
t(out4)