# Program 24
# Initiates a search for a Bayesian D-optimal design.  Uses output from Program 23

m <- 2
s <- 10
lim1 <- m*s
npts <- length(wts)

fx <- function(xvec)
{
 fvec <- c(1,xvec)
 fvec
}

detinfomatrix <- function(betavec,xvec,deswts)
{#calculates the determinant of the information matrix for
#a logistic regression in m variables.
#The input consists of the values of betavec, xvec and deswts
 xmat <- t(matrix(xvec,s,m))
 fxmat <- apply(xmat,2,fx)
 etavec <- as.vector(t(betavec)%*%fxmat)
 expetavec <- exp(etavec)
 modelwtvec <- expetavec/((1+expetavec)^2)
 infomat <- fxmat %*% diag(deswts*modelwtvec) %*% t(fxmat)
 det(infomat)
}

#calculate an approximation to the utility function

objective <- function(variables)
{
 xvec <- matrix(cos(pi*variables[1:lim1]),m,s,byrow=T)
 zvec <- variables[(lim1+1):((m+1)*s)]
 deswts <- zvec^2
 deswts <- deswts/sum(deswts)
 approx <- 0
 for (i in 1:npts)
 {
  temp <- detinfomatrix(abscissae[i,],xvec,deswts)
  approx <- approx + wts[i]*log(temp)
 }
 -approx
}

nsims <- 100
mindet <- 1000000
for (i in 1:nsims)
{
 initial <- c(runif(lim1),runif(s))
 out <- optim(initial,objective,NULL,method="Nelder-Mead")
 if(out$value < mindet) {mindet <- out$value
 bestdesign <- out$par
  }
 }
answer <- bestdesign
ansa <- matrix(cos(pi*answer[1:lim1]),m,s,byrow=T)
zvec <- answer[(lim1+1):((m+1)*s)]
deswts <- zvec^2
deswts <- deswts/sum(deswts)
solution <- rbind(ansa,deswts)
solution
mindet