Nothing
## Confidence Intervals for quantiles
quantileCI <- function(x, prob = 0.5, conf.level = 0.95, method = "exact",
R = 9999, bootci.type = c("norm", "basic", "perc", "bca"),
na.rm = FALSE, alternative = c("two.sided", "less", "greater"), ...){
if(!is.na(pmatch(method, "exact")))
method <- "exact"
METHODS <- c("exact", "asymptotic", "boot")
method <- pmatch(method, METHODS)
alternative <- match.arg(alternative)
if(is.na(method))
stop("invalid method")
if(method == -1)
stop("ambiguous method")
stopifnot(is.numeric(x), is.numeric(prob), is.numeric(conf.level))
if(length(x) <= 1)
stop("'x' has to be of at least length 2")
if(length(prob) != 1)
stop("'prob' has to be of length 1 (quantile)")
if(prob <= 0 | prob >= 1)
stop("'prob' has to be in (0, 1)")
if(length(conf.level) != 1)
stop("'conf.level' has to be of length 1 (confidence level)")
if(conf.level < 0.5 | conf.level > 1)
stop("'conf.level' has to be in [0.5, 1]")
alpha <- 1 - conf.level
if(alternative != "two.sided") alpha <- 2*alpha
z <- qnorm(1-alpha/2)
if(na.rm) x <- x[!is.na(x)]
n <- length(x)
est <- quantile(x, prob = prob)
xs <- sort(x)
if(method == 1){ # exact
if(prob^n + (1-prob)^n > alpha){
fun.n <- function(n, prob, alpha){
prob^n + (1-prob)^n - alpha
}
n.min <- ceiling(uniroot(f = fun.n, lower = 1, upper = 5e4, prob = prob,
alpha = alpha, extendInt = "yes")$root)
warning("There is no exact confidence interval that achieves the requested confidence level of at least ", conf.level,
"\n", "A sample size of at least ", n.min, " would be required!")
}
fun.cover <- function(l, u, n, prob){
pbinom(u-1, size = n, prob = prob) - pbinom(l-1, size = n, prob = prob)
}
res <- outer(1:n, 1:n, fun.cover, n = n, prob = prob)
res[lower.tri(res, diag = TRUE)] <- NA
res <- as.vector(res)
names.res <- outer(1:n, 1:n, function(x, y) paste(x, y, sep = " | "))
names(res) <- names.res
if(any(res > conf.level, na.rm = TRUE)){
res <- res[res > 1-alpha]
ind.min <- which.min(res - (1-alpha))
}else{
ind.min <- which.min(abs(res-(1-alpha)))
}
conf.level.exact <- res[ind.min]
names(conf.level.exact) <- NULL
res.sel <- res[which(res == conf.level.exact)]
if(length(res.sel) > 1){
bounds.all <- do.call("rbind", strsplit(names(res.sel), " \\| "))
CI.length <- numeric(nrow(bounds.all))
for(i in 1:nrow(bounds.all)){
ind <- as.numeric(bounds.all[i,])
CI.length[i] <- diff(xs[ind])
}
ind.min <- which.min(CI.length)
bounds <- as.numeric(bounds.all[ind.min,])
}else{
bounds <- as.numeric(strsplit(names(res.sel), " \\| ")[[1]])
}
lower <- bounds[1]
upper <- bounds[2]
# plev <- 0
# upper <- ceiling((n+1)*prob)
# lower <- floor((n+1)*prob)
# while(plev < conf.level){
# upper <- upper + 1
# lower <- lower - 1
# plev <- pbinom(upper-1, size = n, prob = prob) - pbinom(lower-1, size = n, prob = prob)
# }
# lower <- max(lower, 1)
# upper <- min(upper, n)
CI.lower <- ifelse(alternative == "less", -Inf, xs[lower])
CI.upper <- ifelse(alternative == "greater", Inf, xs[upper])
CI <- matrix(c(CI.lower, CI.upper), nrow = 1)
if(alternative != "two.sided") conf.level.exact <- 1-(1-conf.level.exact)/2
alpha.exact <- 1 - round(conf.level.exact, 3)
attr(CI, "conf.level") <- round(conf.level.exact, 3)
rownames(CI) <- paste(100*prob, "% quantile")
if(alternative == "two.sided")
colnames(CI) <- c(paste(alpha.exact/2*100, "%"), paste((1-alpha.exact/2)*100, "%"))
if(alternative == "less")
colnames(CI) <- c("0 %", paste((1-alpha.exact)*100, "%"))
if(alternative == "greater")
colnames(CI) <- c(paste(alpha.exact*100, "%"), "100 %")
meth <- "exact confidence interval"
}
if(method == 2){ # approx
prob.sd <- sqrt(n*prob*(1-prob))
k.lo <- max(1, floor(n*prob - z*prob.sd))
k.up <- min(n, ceiling(n*prob + z*prob.sd))
CI.lower <- ifelse(alternative == "less", -Inf, xs[k.lo])
CI.upper <- ifelse(alternative == "greater", Inf, xs[k.up])
CI <- matrix(c(CI.lower, CI.upper), nrow = 1)
attr(CI, "conf.level") <- conf.level
rownames(CI) <- rep(paste(100*prob, "% quantile"), nrow(CI))
if(alternative == "two.sided")
colnames(CI) <- c(paste(alpha/2*100, "%"), paste((1-alpha/2)*100, "%"))
if(alternative == "less")
colnames(CI) <- c("0 %", paste((1-alpha/2)*100, "%"))
if(alternative == "greater")
colnames(CI) <- c(paste(alpha/2*100, "%"), "100 %")
meth <- "asymptotic confidence interval"
}
if(method == 3){ # boot
boot.quant <- function(x, i){
quantile(x[i], probs = prob)
}
boot.out <- boot(x, statistic = boot.quant, R = R, ...)
CI <- try(boot.ci(boot.out, type = bootci.type, conf = 1 - alpha),
silent = TRUE)
if(inherits(CI, "try-error"))
stop("Function 'boot.ci' returned an error. Please try a different 'bootci.type'.")
if(alternative == "less"){
if("normal" %in% names(CI)){
CI$normal[1,1] <- conf.level
CI$normal[1,2] <- -Inf
}
if("basic" %in% names(CI)){
CI$basic[1,1] <- conf.level
CI$basic[1,4] <- -Inf
}
if("student" %in% names(CI)){
CI$student[1,1] <- conf.level
CI$student[1,4] <- -Inf
}
if("percent" %in% names(CI)){
CI$percent[1,1] <- conf.level
CI$percent[1,4] <- -Inf
}
if("bca" %in% names(CI)){
CI$bca[1,1] <- conf.level
CI$bca[1,4] <- -Inf
}
}
if(alternative == "greater"){
if("normal" %in% names(CI)){
CI$normal[1,1] <- conf.level
CI$normal[1,3] <- Inf
}
if("basic" %in% names(CI)){
CI$basic[1,1] <- conf.level
CI$basic[1,5] <- Inf
}
if("student" %in% names(CI)){
CI$student[1,1] <- conf.level
CI$student[1,5] <- Inf
}
if("percent" %in% names(CI)){
CI$percent[1,1] <- conf.level
CI$percent[1,5] <- Inf
}
if("bca" %in% names(CI)){
CI$bca[1,1] <- conf.level
CI$bca[1,5] <- Inf
}
}
meth <- "bootstrap confidence interval"
}
names(est) <- paste(100*prob, "% quantile")
return(structure(list("estimate" = est, "conf.int" = CI,
"method" = meth),
class = "confint"))
}
medianCI <- function(x, conf.level = 0.95, method = "exact", R = 9999,
bootci.type = c("norm", "basic", "perc", "bca"), na.rm = FALSE,
alternative = c("two.sided", "less", "greater"), ...){
res <- quantileCI(x, prob = 0.5, conf.level = conf.level, method = method,
R = R, bootci.type = bootci.type,
na.rm = na.rm, alternative = alternative, ...)
if(method != "boot"){
rownames(res$conf.int) <- rep("median", nrow(res$conf.int))
}
names(res$estimate) <- "median"
res
}
madCI <- function(x, conf.level = 0.95, method = "exact",
R = 9999, bootci.type = c("norm", "basic", "perc", "bca"),
na.rm = FALSE, constant = 1.4826,
alternative = c("two.sided", "less", "greater"), ...){
alternative <- match.arg(alternative)
if(method %in% c("exact", "asymptotic")){
M <- median(x, na.rm = na.rm)
res <- medianCI(constant*abs(x-M), conf.level = conf.level,
method = method, R = R, bootci.type = bootci.type,
na.rm = na.rm, alternative = alternative, ...)
if(method != "boot"){
rownames(res$conf.int) <- rep("MAD", nrow(res$conf.int))
}
names(res$estimate) <- "MAD"
}else{
if(method == "boot"){ # boot
if(na.rm) x <- x[!is.na(x)]
alpha <- 1 - conf.level
if(alternative != "two.sided") alpha <- 2*alpha
est <- mad(x, constant = constant)
boot.mad <- function(x, i, constant){
mad(x[i], constant = constant)
}
boot.out <- boot(x, statistic = boot.mad, R = R, constant = constant, ...)
CI <- try(boot.ci(boot.out, type = bootci.type, conf = 1 - alpha),
silent = TRUE)
if(inherits(CI, "try-error"))
stop("Function 'boot.ci' returned an error. Please try a different 'bootci.type'.")
if(alternative == "less"){
if("normal" %in% names(CI)){
CI$normal[1,1] <- conf.level
CI$normal[1,2] <- -Inf
}
if("basic" %in% names(CI)){
CI$basic[1,1] <- conf.level
CI$basic[1,4] <- -Inf
}
if("student" %in% names(CI)){
CI$student[1,1] <- conf.level
CI$student[1,4] <- -Inf
}
if("percent" %in% names(CI)){
CI$percent[1,1] <- conf.level
CI$percent[1,4] <- -Inf
}
if("bca" %in% names(CI)){
CI$bca[1,1] <- conf.level
CI$bca[1,4] <- -Inf
}
}
if(alternative == "greater"){
if("normal" %in% names(CI)){
CI$normal[1,1] <- conf.level
CI$normal[1,3] <- Inf
}
if("basic" %in% names(CI)){
CI$basic[1,1] <- conf.level
CI$basic[1,5] <- Inf
}
if("student" %in% names(CI)){
CI$student[1,1] <- conf.level
CI$student[1,5] <- Inf
}
if("percent" %in% names(CI)){
CI$percent[1,1] <- conf.level
CI$percent[1,5] <- Inf
}
if("bca" %in% names(CI)){
CI$bca[1,1] <- conf.level
CI$bca[1,5] <- Inf
}
}
meth <- "bootstrap confidence interval"
}
names(est) <- "MAD"
res <- structure(list("estimate" = est, "conf.int" = CI,
"method" = meth),
class = "confint")
}
res
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.