#
# This file has functions for computing various quantities
# associated with the 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, 
# using the modified NUT approximation. In particular,
# (a) "critval.approx" computes the critical value,
# (b) "pval.approx" computes the p-value,
# (c) "ucb.approx" computes the UCB of Q(p.0), and
# (d) "lcb.approx" computes the LCB of F(delta.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, particularly for the
# "lcb.approx" function. So if you get an error message, tweak
# these settings a little. 
# 

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

#
# UCB for Q(p.0) based on the approx test
# modified NUT approximation
# recommended only for n <= 30
#

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


#
# LCB for F(delta.0) based on the approx test
# modified NUT approximation
# recommended only for n <= 30
# change llim if "uniroot" does not like it 
# 

lcb.approx <- function(n, delta.0, mu.mle, sigma.mle, llim=0.5, 
                        alpha=0.05)
{
F.hat <- pnorm((delta.0-mu.mle)/sigma.mle) - 
                pnorm((-delta.0-mu.mle)/sigma.mle)
f <- function(x,n,alpha) critval.approx(x,n,alpha)-F.hat
lcb <- uniroot(f, low=llim, up=F.hat, n=n, alpha=alpha)$root
return(lcb) 
}


#
# pvalue for testing H vs K using the approx test
# modified NUT approximation
# recommended only for n <= 30
#

pval.approx <- function(n, delta.0, p.0, mu.mle, sigma.mle)
{
F.hat <- pnorm((delta.0-mu.mle)/sigma.mle) - 
                pnorm((-delta.0-mu.mle)/sigma.mle)
pval <- pt(-qnorm(F.hat)*sqrt(n-1), df=(n-1), 
                ncp=(-sqrt(n)*qnorm(p.0)))
return(pval)
}


#
# Critical point for the approximate test
# (modified NUT approximation)
#

critval.approx <- function(p.0, n, alpha=0.05)
{
ncp <- -sqrt(n)*qnorm(p.0)
qnct <- uniroot(function(x) pt(x, df=(n-1), ncp=ncp)-alpha, 
            c(-1000, qnorm(alpha)))$root
cval <- pnorm(-qnct/sqrt(n-1))
return(cval)
}



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