R/harrell-davis-quantile.R
harrell_davis_quantile.RdEstimates quantiles of equally weighted draws with the
Harrell-Davis estimator (Harrell and Davis 1982)
, a weighted
average of all order statistics. Applied pointwise to the columns of a
draws-by-grid matrix, it gives smooth credible bands: the kinks that the
empirical quantile (stats::quantile()) shows when the draws that define the
quantile change from one grid point to the next are attenuated. The
estimand is the same as that of the empirical quantile; in simulations of
smooth curves the accuracy was similar (see the details).
harrell_davis_quantile(x, probs, names = FALSE)a numeric vector of draws, or a numeric matrix with the draws in rows and one quantity (e.g., a grid point) per column. Missing values are not allowed.
probabilities in [0, 1], in any order; duplicates are
allowed.
whether the result is named by the probabilities (as
stats::quantile() names them). Defaults to FALSE.
a numeric vector of the same length as probs for a vector x,
and a matrix with length(probs) rows and one column per column of x
(keeping its column names) for a matrix x.
For m draws with order statistics y(1), ..., y(m) and a
probability p in (0, 1), the estimate is the sum of w_i * y(i) with
weights w_i = I_(i/m)(a, b) - I_((i-1)/m)(a, b), where I is the
regularized incomplete beta function, a = p * (m + 1), and
b = (1 - p) * (m + 1). The weights depend only on the number of draws and
p and are computed once for all columns that need them. The estimate is
exact for constant columns and does not overflow for finite draws. The
probabilities 0 and 1 return the minimum and the maximum, and a single
draw is returned as it is.
Caveats of the estimator:
The quantiles of the columns are pointwise: the result is not a simultaneous band of the underlying curve.
Smoothness gives no general accuracy guarantee relative to the empirical
quantile, including at extreme probabilities (e.g., 0.0005). Increasing
the number of draws reduces Monte Carlo error; smoothing alone does not
remove it. In simulations
of smooth curves (lines with normal coefficients, probability 0.025,
1,000 and 10,000 draws), the root mean squared error was about as large as
that of the empirical quantile (median ratios of 0.92 and 0.97) and the
roughness of the error along the grid was several times smaller; this is
no guarantee for other draws or settings.
The input draws have equal weights (there are no importance weights), but the estimator gives every draw a positive and unequal weight that depends on its rank and is largest for the draws closest to the target rank. A small number of draws (below roughly 200) and heavy tails can therefore move the estimate far from the empirical quantile.
Draws with a point mass (an atom, e.g., a parameter fixed at zero in part of a model-averaged posterior) are not treated separately: within about two Monte Carlo standard errors (in probability) of the edge of a point mass, the estimate blends the point mass and the continuous draws.
The weighted sum is computed in double precision. Its rounding error is
of the order of the machine precision times the largest absolute draw.
Where the sorted draws jump close to the target rank (e.g., at the edge of
a point mass or in a sparse tail), the rounding of the cell boundaries
i/m can make it larger by a factor that grows roughly in proportion to
sqrt(m * p / (1 - p)). Measured binary-step examples with zeros and
the top r draws equal to one give about 40 with m = 10000, r = 249,
p = 0.975, and about 1,400 with m = 100000, r = 50, p = 0.9995.
Differences between draws that are much smaller than
this error (e.g., the middle draw of c(-1e16, 1, 1e16)) are not
resolved.
An infinite draw has a positive weight at every probability strictly
between 0 and 1, which makes the weighted average infinite (or undefined for
infinite draws of both signs). A column that contains an infinite draw is
therefore summarized by the empirical quantile of that column
(stats::quantile(), type 7). If that quantile is NaN (an interpolation
between draws of -Inf and Inf), or if the weights of a probability
cannot be computed numerically (at tiny probabilities, with a limit that
depends on the draw count and distribution-function implementation), the
function stops with an error of class BayesTools_harrell_davis_undefined
(also of the family class BayesTools_harrell_davis). The weights are
computed only when a column with finite draws needs them, so columns of
infinite draws never raise the second error.
Harrell FE, Davis CE (1982). “A new distribution-free quantile estimator.” Biometrika, 69(3), 635–640. doi:10.1093/biomet/69.3.635 .
set.seed(1)
draws <- rnorm(500)
harrell_davis_quantile(draws, c(.025, .5, .975))
#> [1] -2.03441898 -0.03369132 2.05635874
quantile(draws, c(.025, .5, .975), names = FALSE)
#> [1] -1.9950287 -0.0367783 2.0301936
# pointwise credible band of a line from draws of its intercept and slope
intercept <- rnorm(500)
slope <- rnorm(500, sd = .5)
x_grid <- seq(-1, 1, length.out = 5)
line <- outer(intercept, rep(1, 5)) + outer(slope, x_grid)
harrell_davis_quantile(line, c(.025, .5, .975), names = TRUE)
#> [,1] [,2] [,3] [,4] [,5]
#> 2.5% -2.4734338 -2.38852983 -2.27520798 -2.19579923 -2.28972607
#> 50% -0.0555556 -0.04339342 -0.02831084 -0.02491286 -0.03882316
#> 97.5% 2.3452433 2.13871506 2.01052356 2.00197820 2.09298250