# R/theta_functions.R In npExact: Exact Nonparametric Hypothesis Tests for the Mean, Variance and Stochastic Inequality

```###--------------------------------------------------
### Functions to derive the optimal theta
###--------------------------------------------------

noValidTheta <- simpleError("It was not possible to find a valid theta (i.e., one that minimizes the type II error). Please adjust the test value under the null hypothesis to a less (extreme) value.")

## helper function, to ease the reading of the code
w <- function(x) {
as.numeric(x >= 0)
}

## g1 helper function
## @importFrom stats pbinom
g1 <- function(k, N, z, alpha = 0.05) {
## if alpha < term2 -> 0
## if alpha > term1 -> 1
## else calculate the fraction term
term2 <- pbinom(k, N, z, lower.tail = FALSE)
ifelse(alpha < term2, 0,
ifelse(alpha > pbinom(k - 1, N, z, lower.tail = FALSE),
1,
(alpha - term2)/dbinom(k, N, z)))
}

## g2
## helper function
g2 <- function(mu, N, z, alpha) {
k <- 0:(N - 1)

## term2
## if alpha >= z^N -> 1
## if alpha < z^N -> alpha/z^N
term1 <- ifelse(alpha >= z^N, 1, alpha/z^N)

## g1 is 0 until alpha >= term2 from above
## how to find k^* fast? -> could then sum from there!
res <- sum(dbinom(k, N, mu) * g1(k, N, z, alpha)) + mu^N*term1
res
}

## g2 <- function(mu, N, z, alpha)
##     {
##         k <- 0:(N - 1)
##         term2 <- w(alpha - z^N) + (1 - w(alpha - z^N))*(alpha/z^N)
##         res <- sum(dbinom(k, N, mu) * g1(k, N, z, alpha)) + mu^N*term2
##         res
##     }

## possibleTheta
## Calculates possible values of theta, which are in interval (0,1)
possibleTheta <- function(N, p, alpha) {
k <- 0:N
## j <- lapply(k, function(x)x:N)
## theta <- sapply(j,
##                 function(x)
##                 (1/alpha)*sum(choose(N,x)*p^x*(1-p)^(N-x)))
## r <- sapply(theta, function(x)if(x<1 && x>0) x else NA)
## res <- rbind(k, r)
## res1 <- res[, !is.na(res[2,])]

theta <- (1/alpha) * pbinom(k - 1, N, p,
lower.tail = FALSE)
sensibletheta <- theta < 1 & theta > 0
res <- theta[sensibletheta]
res <- rbind(k[sensibletheta], res)

row.names(res) <- c("k","theta")
res
## print(identical(round(res[2], 5), round(res1[2], 5)))
## list(res, res1)
}

minTypeIIError <- function(p.alt, p, N, alpha, alternative) {
## Calculates minimum value, for given difference d
## uses possibleTheta, g2
theta <- possibleTheta(N, p, alpha)

if(length(theta[2, ]) > 0) {
f <- function(x) {
(1 - g2(p.alt, N, p, alpha*x))/(1 - x)
}

## calculate the type II errors for the thetas
typeIIerrors <- sapply(theta[2,], f)

## if no sensible theta was found, set the type II error to 1
if(!is.numeric(typeIIerrors)) {
typeIIerrors <- 1
}

## the right theta is the one that minimizes the type II error
mintypeII <- min(typeIIerrors, 1)
righttheta <- theta[2, which(typeIIerrors == mintypeII)]
righttheta <- ifelse(length(righttheta) == 0, NA, righttheta)
} else {
stop(noValidTheta)
}

list(theta = righttheta,
typeII = mintypeII)
}

minTypeIIErrorWrapper <- function(p.alt, p, N, alpha, typeIIgoal = .5) {
minTypeIIError(p.alt, p, N, alpha)\$typeII - typeIIgoal
}

## uniroot(minTypeIIErrorWrapper, c(0, 1), p = .1, N = 50, alpha = .05)
```

## Try the npExact package in your browser

Any scripts or data that you put into this service are public.

npExact documentation built on May 2, 2019, 9:58 a.m.