146 lines
5.4 KiB
R
146 lines
5.4 KiB
R
#' Generic infimum‑type quantile function
|
||
#'
|
||
#' @param F A function that evaluates the distribution function
|
||
#' (CDF). It must be monotone non‑decreasing in its first
|
||
#' argument and return values in $0,1$. Additional
|
||
#' arguments for the CDF can be supplied via `...`.
|
||
#' @param p Numeric vector of probabilities. Values must lie in
|
||
#' $0,1$; values outside this range give `NaN`.
|
||
#' @param lower Left end of the search interval. By default
|
||
#' `-Inf` is used; if the CDF is known to be zero below a
|
||
#' finite value you can give that value to speed up the
|
||
#' search.
|
||
#' @param upper Right end of the search interval. By default
|
||
#' `Inf` is used; if the CDF is known to be one above a
|
||
#' finite value you can give that value.
|
||
#' @param tol Desired absolute tolerance for the quantile. The
|
||
#' algorithm stops when the interval width is ≤ `tol`.
|
||
#' @param max.iter Maximum number of bisection iterations (safety
|
||
#' guard). The default is `1000`.
|
||
#' @param ... Additional arguments passed to `F`.
|
||
#'
|
||
#' @return A numeric vector of the same length as `p` containing the
|
||
#' infimum‑type quantiles.
|
||
#' @examples
|
||
#' ## 1. Continuous normal distribution (compare with qnorm)
|
||
#' qinf(pnorm, p = c(0.025, 0.5, 0.975), lower = -10, upper = 10)
|
||
#' qnorm(c(0.025, 0.5, 0.975))
|
||
#'
|
||
#' ## 2. Discrete distribution: Binomial(10, 0.3)
|
||
#' Fbinom <- function(x, size, prob) pbinom(floor(x), size, prob)
|
||
#' qinf(Fbinom, p = seq(0, 1, 0.1), size = 10, prob = 0.3)
|
||
#' qbinom(seq(0, 1, 0.1), size = 10, prob = 0.3) # built‑in quantile
|
||
#'
|
||
#' ## 3. Mixed distribution (continuous + point mass at 0)
|
||
#' Fmix <- function(x, sigma = 1) {
|
||
#' 0.3 * (x >= 0) + 0.7 * pnorm(x, sd = sigma)
|
||
#' }
|
||
#' qinf(Fmix, p = c(0.1, 0.3, 0.5, 0.9), lower = -5, upper = 5)
|
||
#' @export
|
||
qinf <- function(F,
|
||
p,
|
||
lower = -Inf,
|
||
upper = Inf,
|
||
tol = .Machine$double.eps^0.5,
|
||
max.iter = 1000L,
|
||
...) {
|
||
|
||
# 1. Input checks -----------------------------------------------------------
|
||
if (!is.function(F))
|
||
stop("'F' must be a function that evaluates a CDF")
|
||
if (!is.numeric(p))
|
||
stop("'p' must be numeric")
|
||
if (any(is.na(p))) # keep NA positions in the output
|
||
warning("NA probabilities supplied – corresponding quantiles will be NA")
|
||
if (any(p < 0 | p > 1, na.rm = TRUE))
|
||
stop("probabilities must lie in [0,1]")
|
||
|
||
# 2. Helper function for a single probability ---------------------------------
|
||
|
||
one_quantile <- function(p_i) {
|
||
## Edge cases: 0 → leftmost point where F(x) >= 0,
|
||
## 1 → rightmost point where F(x) >= 1
|
||
if (is.na(p_i)) return(NA_real_)
|
||
if (p_i == 0) {
|
||
## The infimum of the set {x : F(x) >= 0} is the lower bound of the
|
||
## support. If the user gave a finite lower bound we return it,
|
||
## otherwise we try to locate it by moving left until F changes.
|
||
if (is.finite(lower)) return(lower)
|
||
## Search leftwards until we see a change in the CDF (or hit -Inf)
|
||
x <- 0
|
||
step <- 1
|
||
while (is.finite(x) && F(x, ...) == 0) {
|
||
x <- x - step
|
||
step <- step * 2
|
||
if (x < -1e12) break # give up – treat as -Inf
|
||
}
|
||
return(x)
|
||
}
|
||
if (p_i == 1) {
|
||
if (is.finite(upper)) return(upper)
|
||
## Search rightwards until F reaches 1
|
||
x <- 0
|
||
step <- 1
|
||
while (is.finite(x) && F(x, ...) < 1) {
|
||
x <- x + step
|
||
step <- step * 2
|
||
if (x > 1e12) break
|
||
}
|
||
return(x)
|
||
}
|
||
|
||
## 2.1 Initialise the bracketing interval =================================
|
||
lo <- lower
|
||
hi <- upper
|
||
|
||
## If the interval is infinite we first try to find a finite bracket.
|
||
## This is done by exponential expansion from 0 (or from the sign of p)
|
||
## until the CDF straddles p.
|
||
if (!is.finite(lo) || !is.finite(hi)) {
|
||
# start from 0 (or any convenient point)
|
||
centre <- 0
|
||
# expand left if needed
|
||
if (!is.finite(lo)) {
|
||
step <- 1
|
||
while (F(centre - step, ...) >= p_i) step <- step * 2
|
||
lo <- centre - step
|
||
}
|
||
# expand right if needed
|
||
if (!is.finite(hi)) {
|
||
step <- 1
|
||
while (F(centre + step, ...) < p_i) step <- step * 2
|
||
hi <- centre + step
|
||
}
|
||
}
|
||
|
||
## 2.2 Bisection loop =====================================================
|
||
## we keep the leftmost point where F(q) >= p
|
||
iter <- 0L
|
||
while ((hi - lo) > tol && iter < max.iter) {
|
||
mid <- (lo + hi) / 2
|
||
fmid <- F(mid, ...)
|
||
if (is.na(fmid)) {
|
||
## If the CDF returns NA (e.g. because mid is outside the domain)
|
||
## we shrink the interval towards the side that is known to be
|
||
## finite.
|
||
if (is.finite(lo)) hi <- mid else lo <- mid
|
||
} else if (fmid >= p_i) {
|
||
## We have reached (or passed) the target probability – move the
|
||
## upper bound leftwards to keep the *infimum*.
|
||
hi <- mid
|
||
} else {
|
||
## Still below the target – move the lower bound rightwards.
|
||
lo <- mid
|
||
}
|
||
iter <- iter + 1L
|
||
}
|
||
|
||
## 2.3 Return the upper bound (the smallest x with F(x) >= p) =============
|
||
hi
|
||
}
|
||
|
||
# 3. Vectorised call – preserve names / attributes of p ---------------------
|
||
out <- vapply(p, one_quantile, numeric(1), USE.NAMES = FALSE)
|
||
if (!is.null(names(p))) names(out) <- names(p)
|
||
out
|
||
} |