# Program 10
# This program considers choosing between quadratic and linear regression
# models, and looks for the Ds_optimal design  The design does not depend
# values of the parameters beta0, beta1 and beta2
#
f <- function(x)
{ c(1,x, x^2)}
f2 <- function(x)
{ c(1,x)}

p <- 3
p2 <- 2
s <- 3
twos <- 2*s


detinfomat <- function(variables)
{
 infomat <- matrix(0,p,p)
 infomat2 <- matrix(0,p2,p2)
 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)
  f2vec <- f2(xvec)
  wt <- 1
  infomat <- infomat + deltavec[i]*wt*fvec%*%t(fvec)
  infomat2 <- infomat2 + deltavec[i]*wt*f2vec%*%t(f2vec)
 }
 -det(infomat)/det(infomat2)
}


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)