# 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)