# Program 29
# This "program" includes the various segments of R code referred to in Examples 5.4.5 and 5.4.6

# Finding values of y for which P(Y = y) >= 0.001
lambda <- 10
ylo <- qpois(0.0005,lambda)
yhi <- qpois(0.9995,lambda)
y1values <- ylo:yhi
prob1 <- dpois(y1values,lambda)
indic <- prob1 >= 0.001
y1values <- y1values[indic]
prob1 <- prob1[indic]
ylo
yhi
y1values

########################################
# Program for Example 5.4.6
# First a program to calculate MSE(x) vs x
fx <- function(x)
{
 c(1,x)
}

s <- 2
betavec <- c(0,1)
nvec <- c(4,4)
xlo <- -2
xhi <- 8
nsubintervals <- 100

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)

p <- length(betavec)
nsimuls <- 20

mse <- function(x,probs1,probs2,mat0,mat1)
{
 lambdastarmat <- exp(mat0 + mat1*x)
 lambdax <- exp(sum(betavec*fx(x)))
 mse <- t(probs1)%*%((lambdastarmat - lambdax)^2)%*%probs2
mse
}

#Plotting the MSE when the support points are at 6.5 and 8
x1 <- 6.5
x2 <- 8.0
distx <- x1 - x2
xvec <- c(x1,x2)
etavec <- t(betavec)%*%sapply(xvec,fx)
lambdavec <- exp(etavec)
muvec <- nvec*lambdavec
y1values <- qpois(0.0001,muvec[1]):qpois(0.9999,muvec[1])
y2values <- qpois(0.0001,muvec[2]):qpois(0.9999,muvec[2])
probs1 <- dpois(y1values,muvec[1])
probs2 <- dpois(y2values,muvec[2])
incid1 <- probs1 >= 0.0001
incid2 <- probs2 >= 0.0001
y1values <- y1values[incid1]
y2values <- y2values[incid2]
probs1 <- probs1[incid1]
probs2 <- probs2[incid2]
a1star <- log((y1values + 0.5)/nvec[1])
a2star <- log((y2values + 0.5)/nvec[2])
mat0 <- outer((-a1star*x2),a2star*x1,"+")/distx
mat1 <- outer(a1star,a2star,"-")/distx
 
par(las=1)
xvalues <- seq(from=-2,to=8,by=0.01)
yvalues <- sapply(xvalues,mse,probs1=probs1,probs2=probs2,mat0=mat0,mat1=mat1)
plot(xvalues,yvalues,ty="l",xlab="z",ylab="MSE")


################################################################################
# Calculate the IMSE when the input is a vector of z-values that will be transformed
# to x-values in the appropriate domain by the first line of the function
imse <- function(zvec)
{
 xvec <- 3 + 5*cos(pi*zvec)
 x1 <- xvec[1]
 x2 <- xvec[2]
 distx <- x1 - x2
 etavec <- t(betavec)%*%sapply(xvec,fx)
 lambdavec <- exp(etavec)
 muvec <- nvec*lambdavec
 y1values <- qpois(0.0001,muvec[1]):qpois(0.9999,muvec[1])
 y2values <- qpois(0.0001,muvec[2]):qpois(0.9999,muvec[2])
 probs1 <- dpois(y1values,muvec[1])
 probs2 <- dpois(y2values,muvec[2])
 incid1 <- probs1 >= 0.0001
 incid2 <- probs2 >= 0.0001
 y1values <- y1values[incid1]
 y2values <- y2values[incid2]
 probs1 <- probs1[incid1]
 probs2 <- probs2[incid2]
 a1star <- log((y1values + 0.5)/nvec[1])
 a2star <- log((y2values + 0.5)/nvec[2])
 mat0 <- outer((-a1star*x2),a2star*x1,"+")/distx
 mat1 <- outer(a1star,a2star,"-")/distx
 
 yvalues <- sapply(xvalues,mse,probs1=probs1,probs2=probs2,mat0=mat0,mat1=mat1)
 integral <- sum(multiplier*yvalues)*h/3
 integral
}


#Try to find the "IMSE-optimal" design using xvec = c(6.272,8) as the initial vector
# of x-values (which must first be back-transformed to z-values)

xguess <- c(6.272,8)
zguess <- acos((xguess-3)/5)/pi

#Try to find "IMSE-optimal" design. "tweaking" the initial guess on each search
nsims <- 20
minimse <- 1.0e08
for(i in 1:nsims)
{
 initial <- zguess + 0.06*(runif(s) - 0.5)
 out <- optim(initial,imse,NULL,method="Nelder-Mead")
 if(out$value < minimse) {minimse <- out$value
    bestdesign <- out$par}
}
answer <-  bestdesign
solution <- 3 + 5*cos(pi*answer)
minimse
solution

#################################################################
# This lets you specify the values of the probabilities qlo and qhi, from which R will
# calculate the quantiles xlo and xhi

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

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