# Program 17
# Constructing the D-optimal design for a Poisson distribution with m variables when
# eta = beta0 + beta1 x1 + ... + beta_{p-1} x_{p-1}

# Test whether the conditions for the design to exist are satisfied
m <- 2
betavec <- c(1,2,-2)
lvec <- c(-1,-1)
uvec <- c(1,1)

construct_Poisson_Dopt <- function()
{
 beta <- betavec[-1]
 product <- abs(beta*(uvec - lvec))
 test <- product < 2
 if(sum(test) > 0){
  stop("Variable(s)",c(1:m)[test]," have domains that are not wide enough")
 }
 
# Conditions satisfied, so proceed
 test <- beta > 0
 cvec <- uvec*test + lvec*(!test)
 xjstar <- matrix(cvec,m,m) - diag(2/beta,m)%*%diag(1,m)
 supports <- t(cbind(xjstar,cvec))
 weights <- rep(1/(m+1),m+1)
 design <- cbind(supports, weights)
 rownames(design) <- seq(1:(m+1))
 design
}

construct_Poisson_Dopt()

###########################################################
# This segment considers the second example in Chapter 5 of finding a D-optimal design
m <- 4
betavec <- c(1,2,1,-1,-2)
lvec <- c(-1,-1,-1,-1)
uvec <- c(1,1,1,1)
construct_Poisson_Dopt()