#
# This file has functions for computing various quantities
# associated with the exact test of the hypotheses
# H: Q(p.0) >= delta.0 versus K: Q(p.0) < delta.0, or 
# equivalently, H: F(delta.0) <= p.0 versus K: F(delta.0) > p.0. 
# In particular,
# (a) "critval.exact" computes the critical value, and
# (b) "ucb.exact" computes the UCB of Q(p.0).
#
# See the "example.R" file for an illustration.
#
# Date: June 24, 2005
#
# Note: These functions rely heavily on the root-finding
# function "uniroot". This function requires the specification
# of an interval that contains the root, and the current
# default settings are not infallible. So if you get an error
# message, tweak these settings a little. 
# 
# Update 1: March 13, 2008. "up" and "low" arguments in "uniroot"
# 			function changed to "upper" and "lower".


################## BEGIN CODE ##############

#
# Exact UCB for Q(p.0)
#

ucb.exact <- function(n, p.0, mu.mle, sigma.mle, alpha=0.05)
{
ncp <- (mu.mle/sigma.mle)^2
cval <- critval.exact(n, p.0, alpha)
ucb <- sqrt(qchisq(cval, df=1, ncp=ncp))*sigma.mle
return(ucb)
}


#
# Exact critical value for testing (H, K)
#

critval.exact <- function (n, p.0, alpha=0.05) 
{
delta.0 <- 1.0
ans <- uniroot (function(x) cvalue.equation (x, n, p.0, delta.0, alpha), 
        lower=p.0, upper=0.9999) $ root
return(ans)
}


#
# maximize the type-I error prob on the boundary of H
# and set up the equation
#

cvalue.equation <- function(cvalue, n, p.0, delta.0, alpha) 
{
ans <- optimize(function(x) obj.function(x, cvalue, n, p.0, delta.0), 
                interval=c(0, 1), maximum = TRUE)$objective - alpha
return(ans)
}

#
# type-I error error prob on the boundary
#

obj.function <- function(mu, cvalue, n, p.0, delta.0)
{
sigma <- uniroot (function(x) sigma.equation(x, mu, p.0, delta.0), 
            lower=0, upper = (delta.0/qnorm( (1.0+p.0)/2.0) ) )$root
ans <- power.function (cvalue, mu, sigma, n, delta.0)
return(ans)
}


#
# solve for sigma for a given mu on the boundary
#

sigma.equation <- function(sigma, mu, p.0, delta.0)
{
z.r <- (delta.0-mu)/sigma
z.l <- (-delta.0-mu)/sigma
ans <- pnorm(z.r) - pnorm(z.l) - p.0
return(ans)
}


#
# the power function 
#


power.function <- function(cvalue, mu, sigma, n, delta.0) 
{
q <- qnorm((1.0+cvalue)/2.0)
up <-  n*(delta.0/(sigma*q))^2
ans <- integrate(function(y) power.integrand(y, cvalue=cvalue, 
                mu = mu, sigma = sigma, n =n, delta.0 = delta.0), 
                lower = 0, upper = up)$value
return(ans)
}


#
# the integrand for the power function
# 

power.integrand <- function(w, cvalue, mu, sigma, n, delta.0) 
{
m <- length(w)
ans <- numeric(m)
for (i in 1:m)
{
t<- uniroot (function(x) t.equation (x, w[i], cvalue, n, sigma, delta.0),
                 lower=0, upper = delta.0) $root
z.r <- sqrt(n)*(t- mu)/sigma
z.l <- sqrt(n)*(-t - mu)/sigma
ans[i] <- (pnorm(z.r) - pnorm(z.l))*dchisq(w[i], df=(n-1))
}
return(ans)
}
    

#
# the t-equation in the integrand
#

t.equation <- function(t, w, cvalue, n, sigma, delta.0) 
{ 
s <- sigma*sqrt(w/n)
z.r <- (delta.0-t)/s
z.l <- (-delta.0-t)/s
ans <- pnorm(z.r)-pnorm(z.l)-cvalue
return(ans)
}



#################  END OF CODE ##################