#### Base functions for calculating UCLs. ####
#These functions are not intended to be used directly. They are the distilled
#essence of the UCL calculations included in rucl. They assume that all
#necessary inputs have been calculated previously. This approach was taken in an
#attempt to maximize efficiency. Who knows if it really worked though... If you
#are interested in directly accessing the functionality provided by these
#functions please refer to the ucl function.
### Naming Convention ###
# Each of these functions begin with a single character followed by a period.
# This initial character indicated which type of dataset each is intended for as
# follows:
# b - function will work with either a censored or uncensored dataset
# u - function will only work with an uncensored dataset
# c - function will only work with a censored data
### Arguement Convention ### Most UCL calculations require similar inputs. Short
#(usually 1-2 letters) argument names were used for efficiency. The following is a
#key to those short names:
# m - the arithmetic mean of the dataset (log-transformed for b.hucl)
# se - the standard error of the mean of the dataset
# sd - the standard deviation (of the log transformed data)
# n - the sample size
# con - the desired confidence for the UCL (1 - alpha) (0.95 by default)
# u3 - the unbiased third moment estimate
# skew - the skewness of the dataset
# k - the shape parameter for the gamma distribution
# v - the variance of the dataset
# B - a bootstrap dataset containing N bootstrap estimates of a population parameter
# N - the number of bootstrap iteration that need to be performed
# x - the dataset being analyzed
# d - a logical vector indicating whether the values in 'x' are detections or
# non-detects: TRUE = detection, FALSE = nondetect
# Student's t UCL
b.tucl <- function(m, se, n, con, ...) m + qt(con, n - 1) * se
# UCL based on central limit theorem
b.zucl <- function(m, se, con, ...) m + qnorm(con) * se
# Modified Student's t UCL
b.modt <- function(m, u3, se, n, con) {
m + u3 / (6 * (se * sqrt(n))^2 * n) + qt(con, n - 1) * se
}
# Modified Central Limit Theorem (CLT) UCL
b.modz <- function(m, se, skew, n, con) {
m + (qnorm(con) + (skew / (6 * sqrt(n))) * (1 + 2 * qnorm(con)^2)) * se
}
# Land's H UCL !!PROVIDE WITH LOG TRANSFORMED DATA!!
b.hucl <- function(m, sd, n, con, ...) {
H <- mns.interpolate(mns.h95, new.x = c(n = n, sd = sd))
exp(m + 0.5 * sd^2 + sd * H/sqrt(n - 1))
}
# Chebyshev Inequality
b.cheb <- function(m, se, con) m + sqrt((1 / (1 - con)) - 1) * se
# Approximate Gamma
b.appgamma <- function(m, k, n, con, ...) 2 * n * k * m /
qchisq(1 - con, 2 * n * k)
# Adjusted Gamma
b.adjgamma <- function(m, k, n, con, ...) {
b <- max(0, mns.interpolate(mns.adja, new.x = c(n=n, a=1-con)))
2 * n * k * m / qchisq(b, 2*n*k)
}
# Chebyshev MVUE
b.mvue <- function(m, v, n, con) {
g <- function(t, n) {
w <- 1
new.w <- 1
p <- 1
while(new.w > 5e-10 & p < 100) {
new.num <- ((n - 1)^(p*2 - 1) * t^p)
new.den <- (factorial(p) * n^p * prod(seq(n+1, by=2, length.out = p-1)))
new.w <- new.num/new.den
w <- w + new.w
p = p + 1
}
if(p >= 100) stop("Too many iterations.")
w
}
mean.mvue <- exp(m) * g(v/2, n)
sd.mvue <- sqrt(exp(2 * m) * (g(v/2, n)^2 - g((n - 2) * v / (n - 1), n)))
mean.mvue + (sqrt((1 / (1 - con)) - 1) * sd.mvue)
}
### Bootstrap UCL functions ###
# Standard Bootstrap
b.zboot <- function(B, con) {
N <- length(B)
m <- sum(B)/N
sigma.B <- sqrt((1/(N-1))*sum((B - m)^2))
b.zucl(m, sigma.B, con)
}
# Percentile Bootstrap
b.pboot <- function(B, con) as.vector(quantile(B, con)[[1]])
# Uncensored Student's-t bootstrap
u.tboot <- function(x, m, se, con, N, ...) {
ti <- function(data, m) {
n <- length(data)
m.i <- sum(data)/n
s.i <- sqrt(sum((data-m.i)^2)/(n-1))
sqrt(n)*((m.i - m)/s.i)
}
B <- u.bootstrap(x, N, ti, m = m)
m - as.vector(quantile(B, 1-con)) * se
}
# Censored Student's-t bootstrap
c.tboot <- function(x, d, m, se, con, N, ...) {
ti <- function(x, d, m) {
n <- length(x)
p <- ple(x, d)
sqrt(n)*((p$mean - m)/(p$se*sqrt(n)))
}
B <- c.bootstrap(x, d, N, ti, m = m)
m - as.vector(quantile(B, 1-con)) * se
}
# Hall's Bootstrap
b.hallboot <- function(x, m, n, skew, se, con, N, ...) {
Wi <- function(data, m) {
n.i <- length(data)
m.i <- sum(data)/n.i
s.i <- sqrt(sum((data-m.i)^2)/(n.i-1))
w.i <- (m.i - m)/s.i
k3.i <- skew(data)
w.i + k3.i * w.i^2 / 3 + k3.i^2 * w.i^3 / 27 + k3.i / (6*n.i)
}
m - 3 * ((1 + skew * (as.vector(quantile(u.bootstrap(x, N, Wi, m = m), 1-con))
- skew / (6 * n)))^(1/3) - 1) / skew * (se * sqrt(n))
}
# Uncensored BCA (bias corrected accelerated) bootstrap
u.bcaboot <- function(B, x, m, n, con, N, ...) {
m.B <- sum(B)/N
z.0 <- qnorm(length(B[B < m.B])/N)
m.i <- sapply(seq(n), function(i, d, n) sum(d[-i])/n, d = x, n = n-1)
a.hat <- sum((m - m.i)^3) / (6 * (sum((m - m.i)^2))^1.5)
a.2 <- pnorm(z.0 + ((z.0 + qnorm(con)) / (1 - a.hat*(z.0 + qnorm(con)))))
as.vector(quantile(B, a.2)[[1]])
}
# Censored BCA (bias corrected accelerated) bootstrap
c.bcaboot <- function(B, x, d, m, n, con, N, ...) {
m.B <- sum(B)/N
z.0 <- qnorm(length(B[B < m.B])/N)
m.i <- sapply(seq(n), function(i, x, d) ple.lite(x[-i], d[-i]), x = x, d = d)
a.hat <- sum((m - m.i)^3) / (6 * (sum((m - m.i)^2))^1.5)
qcon <- qnorm(con)
a.2 <- pnorm(z.0 + (z.0 + qcon) / (1 - a.hat*(z.0 + qcon)))
as.vector(quantile(B, a.2)[[1]])
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.