Highest Density Region-Based Designs

• Highest Density Region-Based Designs are set construction methods that identify minimal-volume regions defined by density superlevel sets with prescribed probability content.
• They are applied in diverse fields such as environmental contour design, surrogate-model sampling, and conformal prediction to isolate the most probable data subsets.
• Methodologies optimize boundary recovery and bandwidth selection by addressing estimation challenges near density thresholds to balance bias and variance.

In the literature surveyed here, the common object underlying highest density region-based designs is a density superlevel set with prescribed probability content. In the univariate formulation, the $100(1-\tau)\%$ HDR is $R_\tau=\{x\in\mathbb R:f(x)\ge f_\tau\}$, where $f_\tau$ is the unique level satisfying $\int f(x)\mathbf 1_{\{f(x)\ge f_\tau\}}\,dx=1-\tau$; in multivariate form, the same construction appears as $R(f_m)=\{x:f(x)\ge f_m\}$ with $\int_{R(f_m)}f(x)\,dx=1-\alpha$ and boundary $C(f_m)=\{x:f(x)=f_m\}$ (Samworth et al., 2010, Haselsteiner et al., 2017). This superlevel-set viewpoint recurs across density estimation, environmental contour construction, conformal prediction, surrogate-model sampling, and graph-based region search, because it isolates the most probable part of a distribution while allowing disconnected or nonconvex geometry when the underlying density is multimodal (Sampson et al., 2024, Minini et al., 12 Sep 2025).

1. Formal structure and geometric meaning

An HDR is a probability-content-constrained level set. In the statistical formulation used for density estimation, the inferential target is not the entire density $f$ but the set where $f$ exceeds the threshold $f_\tau$, and this threshold is determined by coverage rather than by a user-specified density value (Samworth et al., 2010). In environmental design, the same structure is written as $R_\tau=\{x\in\mathbb R:f(x)\ge f_\tau\}$0, with $R_\tau=\{x\in\mathbb R:f(x)\ge f_\tau\}$1 chosen so that the enclosed set contains probability $R_\tau=\{x\in\mathbb R:f(x)\ge f_\tau\}$2; the corresponding contour $R_\tau=\{x\in\mathbb R:f(x)\ge f_\tau\}$3 is the highest density contour (HDC) (Haselsteiner et al., 2017). In conditional prediction, the analogous object is $R_\tau=\{x\in\mathbb R:f(x)\ge f_\tau\}$4, or, for standardized residuals, $R_\tau=\{x\in\mathbb R:f(x)\ge f_\tau\}$5 (Sampson et al., 2024, Sampson et al., 2024).

A central property is minimum volume at fixed probability content. The environmental-contour formulation states explicitly that the HDR occupies the smallest possible volume in variable space among all regions with the same included probability, which gives the HDC a direct interpretation as the tightest region of likely design conditions (Haselsteiner et al., 2017). The uncertainty-quantification framework QUEST makes the same point in measure-theoretic language, defining $R_\tau=\{x\in\mathbb R:f(x)\ge f_\tau\}$6 and $R_\tau=\{x\in\mathbb R:f(x)\ge f_\tau\}$7, so that uncertainty is quantified by the Lebesgue volume needed to capture $R_\tau=\{x\in\mathbb R:f(x)\ge f_\tau\}$8 of the mass (Goring et al., 17 Jun 2026).

The geometry is boundary-driven. In the univariate asymptotic theory, $R_\tau=\{x\in\mathbb R:f(x)\ge f_\tau\}$9 is determined by finitely many crossing points $f_\tau$0 satisfying $f_\tau$1, with alternating derivative signs, so the estimation problem reduces to recovering where $f_\tau$2 changes sign near the boundary (Samworth et al., 2010). In multivariate settings, the same principle appears through boundary manifolds $f_\tau$3, and risk expansions are written as Hausdorff integrals over that boundary (Doss et al., 2018).

Discrete spaces require a modified formulation. For countable $f_\tau$4, the relevant object is the canonical smallest covering region, defined by minimizing $f_\tau$5 subject to $f_\tau$6, then maximizing $f_\tau$7 among all sets of that optimal size. Its weak HDR characterization is $f_\tau$8, $f_\tau$9 for all $\int f(x)\mathbf 1_{\{f(x)\ge f_\tau\}}\,dx=1-\tau$0, and $\int f(x)\mathbf 1_{\{f(x)\ge f_\tau\}}\,dx=1-\tau$1, meaning every included atom is at least as probable as every excluded atom (O'Neill, 2022). This replaces the naive continuous threshold rule in the presence of ties and countable support.

2. Statistical estimation and bandwidth selection

The classical estimator is plug-in thresholding of a kernel density estimator,

$\int f(x)\mathbf 1_{\{f(x)\ge f_\tau\}}\,dx=1-\tau$2

with estimated level $\int f(x)\mathbf 1_{\{f(x)\ge f_\tau\}}\,dx=1-\tau$3 and estimated region $\int f(x)\mathbf 1_{\{f(x)\ge f_\tau\}}\,dx=1-\tau$4 (Samworth et al., 2010). What distinguishes HDR estimation from ordinary density estimation is the loss: the paper on univariate asymptotics uses

$\int f(x)\mathbf 1_{\{f(x)\ge f_\tau\}}\,dx=1-\tau$5

so risk is the $\int f(x)\mathbf 1_{\{f(x)\ge f_\tau\}}\,dx=1-\tau$6-probability content of the symmetric difference rather than an integrated squared error for the density itself (Samworth et al., 2010). This directly measures the probability that a fresh draw falls in exactly one of the true and estimated HDRs.

The main asymptotic message is that bandwidths optimal for global density estimation are generally not optimal for HDR recovery. The univariate theory shows that the relevant terms are local to the boundary points $\int f(x)\mathbf 1_{\{f(x)\ge f_\tau\}}\,dx=1-\tau$7, involving $\int f(x)\mathbf 1_{\{f(x)\ge f_\tau\}}\,dx=1-\tau$8, $\int f(x)\mathbf 1_{\{f(x)\ge f_\tau\}}\,dx=1-\tau$9, $R(f_m)=\{x:f(x)\ge f_m\}$0, and the interaction between pointwise density estimation and threshold estimation (Samworth et al., 2010). The resulting bias-variance balance is nonstandard but still yields the familiar order $R(f_m)=\{x:f(x)\ge f_m\}$1; what changes is the optimal constant, because the threshold $R(f_m)=\{x:f(x)\ge f_m\}$2 is itself random and materially contributes to risk (Samworth et al., 2010).

Uniform-in-bandwidth expansions are crucial because bandwidth selection requires optimization over $R(f_m)=\{x:f(x)\ge f_m\}$3. The remainder term in the risk approximation is uniform over $R(f_m)=\{x:f(x)\ge f_m\}$4, which is precisely what makes minimization of the asymptotic risk meaningful (Samworth et al., 2010). The same principle is extended to $R(f_m)=\{x:f(x)\ge f_m\}$5 through bandwidth matrix selection for KDE level sets and HDRs, using the symmetric-difference loss $R(f_m)=\{x:f(x)\ge f_m\}$6 and deriving a multivariate asymptotic approximation $R(f_m)=\{x:f(x)\ge f_m\}$7 that includes an HDR-specific threshold-estimation term $R(f_m)=\{x:f(x)\ge f_m\}$8 (Doss et al., 2018). That paper turns the expansion into a feasible plug-in selector, with pilot estimators for $R(f_m)=\{x:f(x)\ge f_m\}$9, $\int_{R(f_m)}f(x)\,dx=1-\alpha$0, and $\int_{R(f_m)}f(x)\,dx=1-\alpha$1, and implements the method in the R package lsbs (Doss et al., 2018).

The statistical implication is direct: if the design target is a level-set-type region, smoothing should be tuned for set recovery, not for global density fit. This principle is explicit in both the univariate and multivariate bandwidth-selection papers (Samworth et al., 2010, Doss et al., 2018).

3. Computational variants and generalized HDR estimators

When the distribution is known analytically, HDR computation can be reformulated as a nonlinear optimization problem over candidate regions. For continuous univariate distributions with known density and quantile functions, monotone densities yield one-sided intervals, strictly quasi-concave densities yield a single interval $\int_{R(f_m)}f(x)\,dx=1-\alpha$2 obtained by minimizing width $\int_{R(f_m)}f(x)\,dx=1-\alpha$3, strictly quasi-convex densities yield unions of boundary intervals, and general multimodal cases reduce to allocation across segments separated by local minima (O'Neill, 2022). This provides exact computational infrastructure for model-based HDR constructions.

When the density is unknown, the paper "Alternative Approaches for Estimating Highest-Density Regions" generalizes Hyndman’s density-quantile method by replacing the density estimate with any neighbourhood measure $\int_{R(f_m)}f(x)\,dx=1-\alpha$4 that preserves the density-induced ordering asymptotically (Deliu et al., 2023). Theorem 1 there shows that if $\int_{R(f_m)}f(x)\,dx=1-\alpha$5 is a sparsity measure, the estimated HDR is $\int_{R(f_m)}f(x)\,dx=1-\alpha$6, whereas for a concentration measure it is $\int_{R(f_m)}f(x)\,dx=1-\alpha$7 (Deliu et al., 2023). The paper studies direct KDE, $\int_{R(f_m)}f(x)\,dx=1-\alpha$8-nearest-neighbour distances, local $\int_{R(f_m)}f(x)\,dx=1-\alpha$9-neighbourhood probability masses, and copula-based density or local-mass scores, and reports that the copula-based strategy is especially advantageous in scenarios with multimodality or particular dependence structures (Deliu et al., 2023).

Directional domains require their own HDR machinery. For data on $C(f_m)=\{x:f(x)=f_m\}$0, the directional HDR is $C(f_m)=\{x:f(x)=f_m\}$1, estimated by a von Mises–Fisher kernel density estimator and a plug-in threshold based on the empirical quantile of $C(f_m)=\{x:f(x)=f_m\}$2 (Saavedra-Nieves et al., 2020). The paper proposes a bootstrap bandwidth selector

$C(f_m)=\{x:f(x)=f_m\}$3

thereby targeting boundary reconstruction rather than global MISE (Saavedra-Nieves et al., 2020).

Geometric regularization can also be imposed after smoothing. In the COVID-19 hotspot paper, the classical estimator $C(f_m)=\{x:f(x)=f_m\}$4 is compared with a hybrid estimator based on $C(f_m)=\{x:f(x)=f_m\}$5-convex hulls, where the shape parameter is estimated from high-density and low-density subsamples separated by a bootstrap-calibrated uncertainty band (Saavedra-Nieves, 2020). The paper reports that the hybrid method performs best for large sample sizes and large $C(f_m)=\{x:f(x)=f_m\}$6, that is, when the inferential target is the most concentrated modal core rather than a broad effective support (Saavedra-Nieves, 2020).

4. HDRs as design objects across application domains

In probabilistic engineering design, the HDC is the boundary of the HDR of the environmental joint density. The contour is defined directly in the original environmental-variable space by $C(f_m)=\{x:f(x)=f_m\}$7, with $C(f_m)=\{x:f(x)=f_m\}$8 and $C(f_m)=\{x:f(x)=f_m\}$9 (Haselsteiner et al., 2017). The paper argues that this contour is particularly attractive because the outside region has total probability $f$0, giving a direct return-period interpretation, and because the HDR is the smallest-volume region with that probability content (Haselsteiner et al., 2017). The same paper also stresses a common misconception: minimum volume at fixed probability does not imply that the resulting contour is least conservative in engineering practice; compared with IFORM, the HDC can be larger in original variable space and may be “overly conservative” when structural failure surfaces are known to be convex (Haselsteiner et al., 2017).

For surrogate-model training under uncertainty propagation, HDRs are used as sampling domains rather than as acceptance contours. The proposed method in "A sampling method based on highest density regions" first computes a compact HDR $f$1 with probability $f$2, then draws samples uniformly within that region rather than from the original input law (Minini et al., 12 Sep 2025). The rationale is that natural sampling overconcentrates near the mode, whereas uniform sampling inside the HDR spreads design points more evenly over the probabilistically relevant domain. The paper reports that these HDR-based designs globally outperform random-vector-based designs in relative mean square error and in estimating probability of failure, especially for low-dimensional and moderately correlated inputs (Minini et al., 12 Sep 2025).

In weakly supervised object localization, the same superlevel-set intuition is transferred to graphs of superpixels. The level set maximum-weight connected subgraph (LS-MWCS) solves

$f$3

so that increasing the level parameter $f$4 makes the selected connected region smaller while its average weight per node increases monotonically (Zhao et al., 2014). This is not a classical probabilistic HDR, but it is explicitly framed as a search for the densest connected score region rather than the region with maximum total score, which is the relevant correction when projected classifier scores are nearly all positive (Zhao et al., 2014).

These applications share a common design logic: the operational region is chosen by ranking locations, states, or graph nodes by a density-like score and then retaining the most concentrated subset subject to a content or size criterion. The details differ, but the boundary-first geometry is the same.

5. Predictive regions, conformalization, and uncertainty quantification

A major recent development is the use of HDRs as predictive sets. In heteroscedastic regression, the paper "Highest Probability Density Conformal Regions" shows that the oracle predictive region can be written as $f$5, where $f$6 is the density of the standardized residual $f$7 (Sampson et al., 2024). The practical method estimates the one-dimensional residual density by KDE, extracts the smallest estimated $f$8 upper density set, and then conformalizes its interval endpoints using signed residual scores; the resulting predictive region can be a single interval under unimodality or a union of intervals under multimodality, while preserving finite-sample marginal coverage (Sampson et al., 2024).

A complementary strategy, conformal highest conditional density sets (CHCDS), starts from any conditional density estimator $f$9 and its estimated highest-density cutoff $f$0, defines calibration scores

$f$1

and returns

$f$2

where $f$3 is the conformal quantile of the $f$4 (Sampson et al., 2024). The conceptual contribution is that conformalization is achieved by shifting the density threshold rather than expanding an interval in $f$5-space. Under regularity conditions, the conformal adjustment is $f$6, and under correct specification it vanishes asymptotically (Sampson et al., 2024).

Prototype-based discretizations provide a further variant. Conformalized High-Density Quantile Regression (CHDQR) approximates the conditional density by a piecewise-constant model over adaptive Voronoi cells in output space, orders cells by estimated density, and takes the union of the highest-density cells until a conformally calibrated cumulative-mass threshold is reached (Cengiz et al., 2024). This yields non-convex, potentially disconnected prediction regions with marginal coverage guarantees, while using fewer prototypes than fixed grids in the reported two-dimensional experiments (Cengiz et al., 2024).

QUEST turns HDRs into scalar uncertainty measures. For a density $f$7, it defines the $f$8-HDR $f$9 and its volume $f_\tau$0, then uses $f_\tau$1 or the integrated volume $f_\tau$2 as aleatoric or epistemic uncertainty scores (Goring et al., 17 Jun 2026). The paper proves that these measures satisfy monotonicity-under-spread axioms unavailable to variance-based summaries, connects $f_\tau$3 to cross-entropy inside the HDR, and reports favorable selective-prediction performance against variance and differential entropy (Goring et al., 17 Jun 2026).

Finally, conformal prediction regions themselves can be reinterpreted as imprecise highest density regions. Under consonance, the conformal transducer $f_\tau$4 induces a credal set $f_\tau$5, and the paper proves

$f_\tau$6

so the conformal region is exactly the smallest set with lower probability $f_\tau$7 for that credal set (Caprio et al., 10 Feb 2025). This gives conformal threshold sets a robust HDR-style interpretation.

6. Adjacent methodologies, misconceptions, and limits

Several nearby research strands are useful but should not be conflated with HDR estimation proper. Minimum density hyperplanes minimize the integral of the empirical density along a separating hyperplane,

$f_\tau$8

so that the separator avoids cutting through contiguous high-density clusters (Pavlidis et al., 2015). This is a low-density-separation surrogate for HDR structure, not a direct reconstruction of an HDR boundary. Likewise, density-based interpretable hypercube region partitioning for mixed numeric and categorical data constructs a partition of the observed feature space into hyper-rectangles ranked by a density proxy and explicitly carves out empty regions; it is presented as adjacent to HDR methodology rather than a direct HDR estimator (Ackerman et al., 2021). Lattice-based designs with quasi-optimal separation on all projections concern packing density and projection-wise space-filling quality, not probability-density superlevel sets (He, 2017).

The main limitations of HDR-based designs follow from the underlying estimation problem. Core bandwidth theory is one-dimensional in (Samworth et al., 2010), while multivariate HDR bandwidth selection is practically developed mainly for $f_\tau$9 in (Doss et al., 2018). Grid-based HDC computation is straightforward in low dimensions but subject to the curse of dimensionality, and the environmental-contour paper explicitly recommends Monte Carlo methods for higher $R_\tau=\{x\in\mathbb R:f(x)\ge f_\tau\}$00 (Haselsteiner et al., 2017). Hybrid geometric estimators can underperform on complex multimodal structures at moderate sample sizes (Saavedra-Nieves, 2020). Surrogate-model HDR sampling is reported to be most useful for low-dimensional and moderately correlated inputs (Minini et al., 12 Sep 2025). Conformal predictive HDRs retain exact marginal coverage, but conditional coverage remains impossible in the usual distribution-free sense for continuous responses without pathological set size (Sampson et al., 2024).

A persistent misconception is that HDRs are intrinsically connected or convex. The surveyed literature consistently rejects that view. Multimodal environmental distributions can produce disconnected HDCs (Haselsteiner et al., 2017), multimodal residual densities yield unions of predictive intervals (Sampson et al., 2024), CHCDS explicitly treats disconnected conditional sets as a feature rather than a bug (Sampson et al., 2024), and the discrete theory shows that exact HDR analogues may be unions of disjoint atoms or intervals with boundary ties (O'Neill, 2022). Another misconception is that density estimation and HDR estimation have the same smoothing target; the bandwidth-selection papers show that they do not (Samworth et al., 2010, Doss et al., 2018).

Taken together, these results support a unified interpretation. Highest density region-based designs are not a single algorithmic family but a design principle: specify probability content first, rank the space by a density or density-order-preserving score, and optimize the construction for boundary recovery or minimum-volume set recovery rather than for global fit. The versatility of the principle explains its appearance in environmental contours, predictive inference, surrogate training, graph search, and uncertainty quantification, while the limitations show that practical success depends on how accurately the relevant density-induced ordering can be estimated in the domain of interest.