# # 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 ##################