#Program_14
# This program can be used to calculate an IMSE-optimal design for two support points
# when eta = beta0 + beta1 * x and the logit link is used.  It is described in Section 4.8

# This program does not follow exactly the form given at the beginning of Example 4.8.6.  Instead,
# it contains the modifications made to speed up the calculations and described in Subsection 4.8.3.

fx <- function(x)
{
 c(1,x)
}


mse <- function(x,x1,x2)
{#Calculates the MSE at a value x for support points x1 and x2
 distx <- x1 - x2
 eta1 <- sum(fx(x1)*betavec)
 eta2 <- sum(fx(x2)*betavec)
 prob1 <- 1/(1+exp(-eta1))
 prob2 <- 1/(1+exp(-eta2))
 y1probs <- dbinom(y1vals,n1,prob1)
 y2probs <- dbinom(y2vals,n2,prob2)
 pix <- 1/(1+exp(-sum(fx(x)*betavec)))
 mat1 <- outer(a1star*x1,a2star*x2,"-")
 mat2 <- outer((-a1star),a2star,"+")
 pistarxmat <- 1/(1+exp((-mat1 - mat2*x)/distx))
 mse <- t(y1probs)%*%((pistarxmat-pix)^2)%*%y2probs
 mse
}

imse <- function(xvec)
{#Calculates the IMSE when the support points are at xvec[1] and xvec[2]
 xvec2 <- 5*cos(pi*xvec)
 x1 <- xvec2[1]
 x2 <- xvec2[2]
 distx <- x1 - x2
 msevalues <- sapply(xvalues,mse,x1=x1,x2=x2)
 summse <- sum(multiplier*msevalues)
 integral <- (h/3)*summse
 integral
}

#################################################################
#Specify the values of some parameters, which will cause other parameters to be calculated
betavec <- c(0,1)
qlo <- 0.0001
qhi <- 0.9999
nsubintervals <- 50
n1 <- 6
n2 <- 6

xlo <- (log(qlo/(1-qlo)) - betavec[1])/betavec[2]
xhi <- (log(qhi/(1-qhi)) - betavec[1])/betavec[2]

nsteps <- 2*nsubintervals
h <- (xhi - xlo)/nsteps
xvalues <- seq(from=xlo,to=xhi,length=(nsteps+1))
temp <- 2 + 2*(((1:(nsteps-1)) %% 2) == 1)
multiplier <- c(1,temp,1)

y1vals <- 0:n1
y2vals <- 0:n2
a1star <- log((n1-y1vals+0.5)/(y1vals+0.5))
a2star <- log((n2-y2vals+0.5)/(y2vals+0.5))

#######################################################################################
#Try to find IMSE-opt design
nsims <- 100
minimse <- 100
initialz <- acos(c(-1.5434,1.5434)/5)/pi
for(i in 1:nsims)
{
 initialx <- initialz + 0.1*(runif(2)-0.5)
 out <- optim(initialx,imse,NULL,method="Nelder-Mead")
 if(out$value < minimse) {minimse <- out$value
    bestdesign <- out$par}
}
answer <-  bestdesign
xvals <- 5*cos(pi*answer)
minimse
xvals