Isotonic Distributional Regression

Updated 9 July 2026

The paper introduces IDR as a nonparametric method that estimates the full conditional distribution under order restrictions, ensuring stochastic monotonicity.
It employs CRPS-based isotonic regression, using techniques like PAVA, to achieve optimal calibration and efficient computation.
The method extends to various applications including forecasting, semiparametric modeling, and distribution-on-distribution regression with nearly parametric risk rates.

Searching arXiv for key papers on isotonic distributional regression and adjacent theory. Isotonic distributional regression (IDR) is a nonparametric method for estimating conditional distributions under order restrictions. In its canonical form, the covariate space is equipped with a partial order $(\mathcal{X},\preceq)$ , and the target is a family of conditional laws $x\mapsto F_x$ such that $x \preceq x'$ implies $F_x \le_{\mathrm{st}} F_{x'}$ , equivalently $F_x(y)\ge F_{x'}(y)$ for all $y$ , or $F_x^{-1}(\alpha)\le F_{x'}^{-1}(\alpha)$ for all $\alpha\in(0,1)$ (Henzi et al., 2019). The method generalizes classical isotonic regression from single functionals such as means, probabilities, or single quantiles to full conditional distributions, and it is formulated so that the fitted distributions are calibrated and optimal under broad classes of proper scoring rules, most prominently the continuous ranked probability score (CRPS) (Henzi et al., 2019).

1. Conceptual foundations

The central modeling assumption in IDR is stochastic monotonicity. For a real-valued response $Y$ and covariates $X\in\mathcal{X}$ , one seeks a map $x\mapsto F_x$ 0 such that larger covariate values in the partial order correspond to stochastically larger responses. In the standard stochastic order, $x\mapsto F_x$ 1 means

$x\mapsto F_x$ 2

and this is equivalent to monotonicity of all conditional quantiles: $x\mapsto F_x$ 3 This equivalence is foundational because it permits IDR to be viewed either through conditional CDFs or through conditional quantile functions (Henzi et al., 2019).

A closely related formulation predates the term IDR and treats bivariate observations $x\mapsto F_x$ 4 under the sole assumption that $x\mapsto F_x$ 5 is isotonic with respect to stochastic order. In that setting, one can estimate $x\mapsto F_x$ 6 for fixed $x\mapsto F_x$ 7 by isotonic least squares on the indicators $x\mapsto F_x$ 8, or estimate $x\mapsto F_x$ 9 for fixed $x \preceq x'$ 0 by isotonic regression quantiles (Mösching et al., 2019). This earlier formulation already establishes that isotonic distributional regression is fundamentally a regression problem for the full conditional law, not merely for isolated functionals.

In the main IDR formulation, training data $x \preceq x'$ 1 induce a class of isotonic tuples of distributions

$x \preceq x'$ 2

and the estimator is defined as the unique CRPS-based isotonic regression in this class (Henzi et al., 2019). This places IDR in the broader class of shape-constrained distributional estimators, but with stochastic order rather than parametric smoothness as the governing structure.

2. Optimization principles and equivalent representations

The defining loss in the original IDR paper is the continuous ranked probability score,

$x \preceq x'$ 3

and IDR is the unique minimizer of empirical CRPS over the isotonic class (Henzi et al., 2019). A major consequence is a universality property: the same fitted conditional distributions are simultaneously optimal under large classes of proper scoring rules built from quantile scores or threshold scores, and the estimator is threshold calibrated in the sense

$x \preceq x'$ 4

with respect to the empirical joint distribution (Henzi et al., 2019).

A complementary convex-analytic formulation treats distributional regression pointwise in the threshold $x \preceq x'$ 5. For each $x \preceq x'$ 6, let

$x \preceq x'$ 7

and let $x \preceq x'$ 8 encode the structural constraint. In the isotonic case,

$x \preceq x'$ 9

Then the estimator is the Euclidean projection

$F_x \le_{\mathrm{st}} F_{x'}$ 0

A central lemma shows that this projection estimator solves empirical CRPS minimization under the same structural constraint (Padilla et al., 14 May 2025). In the isotonic case, therefore, IDR can be understood as isotonic least squares on the indicator vectors $F_x \le_{\mathrm{st}} F_{x'}$ 1, indexed continuously by $F_x \le_{\mathrm{st}} F_{x'}$ 2.

This projection view also clarifies the relation between CDF-based and quantile-based estimation. In the one-dimensional ordered-covariate setting, isotonic least squares on $F_x \le_{\mathrm{st}} F_{x'}$ 3 yields an estimator $F_x \le_{\mathrm{st}} F_{x'}$ 4 whose plug-in quantiles $F_x \le_{\mathrm{st}} F_{x'}$ 5 form a band of isotonic quantile estimators. The direct isotonic quantile regression problem yields quantile curves lying within that band, so the distribution-function-based construction is more flexible than levelwise isotonic quantile regression (Mösching et al., 2019). This suggests that the distribution-function formulation is the primary object and the quantile formulation is a derived one.

3. Algorithms, interpolation, and computation

In totally ordered settings, the core computation at each threshold reduces to weighted isotonic least squares and is solved by the pool-adjacent-violators algorithm (PAVA). In the early one-dimensional formulation, if $F_x \le_{\mathrm{st}} F_{x'}$ 6 are distinct covariate values with multiplicities $F_x \le_{\mathrm{st}} F_{x'}$ 7, then for fixed $F_x \le_{\mathrm{st}} F_{x'}$ 8 the isotonic CDF estimator solves

$F_x \le_{\mathrm{st}} F_{x'}$ 9

where $F_x(y)\ge F_{x'}(y)$ 0, and PAVA computes the solution in linear time in $F_x(y)\ge F_{x'}(y)$ 1 for each $F_x(y)\ge F_{x'}(y)$ 2 (Mösching et al., 2019).

The original IDR paper extends the method to arbitrary partial orders. At each threshold $F_x(y)\ge F_{x'}(y)$ 3, one solves a quadratic program for the binary responses $F_x(y)\ge F_{x'}(y)$ 4 under isotonic constraints induced by the partial order, and modern QP solvers together with active-set methods are used when PAVA is not directly applicable (Henzi et al., 2019). For new covariate values not in the training set, the paper proposes an interpolation method based on direct predecessors and successors in the partial order. Any admissible predictive CDF must lie between the successor envelope and predecessor envelope, and the proposed predictor takes the midpoint of that envelope; in one-dimensional total orders, this reduces to nearest-neighbor extrapolation outside the training range and linear interpolation between neighboring fitted CDFs inside the range (Henzi et al., 2019).

For large samples, exact repeated isotonic projections across many thresholds become computationally demanding. A dedicated acceleration result exploits the fact that, in isotonic distributional regression, the sequence of response vectors changes only in one coordinate when moving from one ordered response level to the next. The paper develops a modified PAVA and, more importantly, an abridged PAVA that reuses the previous block structure and updates only locally. In the reported simulations, modified PAVA is about $F_x(y)\ge F_{x'}(y)$ 5 faster than standard PAVA, abridged PAVA is about $F_x(y)\ge F_{x'}(y)$ 6 faster than modified PAVA, and overall speed-ups over standard PAVA exceed an order of magnitude and can be much larger when the monotone association is strong (Henzi et al., 2020).

A further computational stabilization device is subsample aggregation. Multiple IDR fits are trained on subsamples and their predicted CDFs are averaged. This preserves stochastic monotonicity, yields smoother fitted distributions, and can improve computational efficiency through parallelization (Henzi et al., 2019).

4. Statistical guarantees and risk theory

The original IDR paper proves existence and uniqueness of the CRPS-based isotonic distributional regression estimator and establishes threshold calibration as well as simultaneous optimality for broad classes of proper scoring rules (Henzi et al., 2019). It also proves uniform consistency on interior subsets of $F_x(y)\ge F_{x'}(y)$ 7 under componentwise ordering, assuming stochastic monotonicity of the true conditional distributions, suitable density of the design points, and uniform continuity of $F_x(y)\ge F_{x'}(y)$ 8 in the covariates (Henzi et al., 2019).

A more recent theory places isotonic distributional regression within a unified convex-constrained framework and derives explicit non-asymptotic risk bounds. In the isotonic setting with

$F_x(y)\ge F_{x'}(y)$ 9

assuming $y$ 0 for all $y$ 1, the estimator satisfies

$y$ 2

and

$y$ 3

These rates match those of isotonic mean regression, formalizing the claim that convergence rates for isotonic distributional regression are consistent with those for mean estimation (Padilla et al., 14 May 2025).

The same theory yields adaptive fast rates when the isotonic signal has low complexity. If

$y$ 4

counts strict jumps at threshold $y$ 5, then

$y$ 6

If the true CDF vectors are piecewise constant with at most $y$ 7 jumps, this becomes $y$ 8, which is nearly parametric when $y$ 9 (Padilla et al., 14 May 2025).

These results also clarify the role of rearrangement. A rearranged estimator $F_x^{-1}(\alpha)\le F_{x'}^{-1}(\alpha)$ 0 that is explicitly forced to be nondecreasing in the threshold argument $F_x^{-1}(\alpha)\le F_{x'}^{-1}(\alpha)$ 1 can be obtained without sacrificing the CRPS rate, ensuring that the final fitted objects are genuine CDFs (Padilla et al., 14 May 2025).

A major semiparametric extension is the Distributional (Single) Index Model (DIM), which combines a parametric index $F_x^{-1}(\alpha)\le F_{x'}^{-1}(\alpha)$ 2 with nonparametric isotonic distributional regression. The model assumes

$F_x^{-1}(\alpha)\le F_{x'}^{-1}(\alpha)$ 3

where $F_x^{-1}(\alpha)\le F_{x'}^{-1}(\alpha)$ 4 is stochastically ordered in the scalar or vector index. Estimation proceeds in two stages: first estimate $F_x^{-1}(\alpha)\le F_{x'}^{-1}(\alpha)$ 5, then apply IDR to the training pairs $F_x^{-1}(\alpha)\le F_{x'}^{-1}(\alpha)$ 6 (Henzi et al., 2020). Under Lipschitz continuity of $F_x^{-1}(\alpha)\le F_{x'}^{-1}(\alpha)$ 7, density of the design in index space, and a near-parametric uniform rate for a monotone transform of $F_x^{-1}(\alpha)\le F_{x'}^{-1}(\alpha)$ 8, the plug-in IDR estimator is uniformly consistent on interior index regions at rate $F_x^{-1}(\alpha)\le F_{x'}^{-1}(\alpha)$ 9 (Henzi et al., 2020). An important invariance property is that IDR depends only on the order induced by the index: strictly monotone transformations of the index leave the fitted distributions unchanged up to reparametrization (Henzi et al., 2020).

Another extension operates in the space of probability measures rather than scalar responses. In distribution-on-distribution regression via optimal transport maps, both covariates and responses are distributions on a compact interval $\alpha\in(0,1)$ 0. The model assumes

$\alpha\in(0,1)$ 1

for an increasing transport map $\alpha\in(0,1)$ 2, and the Fréchet least-squares estimator of $\alpha\in(0,1)$ 3 reduces after discretization to the weighted isotonic least-squares problem

$\alpha\in(0,1)$ 4

which is solved by PAVA (Ghodrati et al., 2021). This places isotonicity at the level of transport maps rather than conditional CDFs and suggests a broader interpretation of isotonic distributional regression as regression in spaces of distributions under monotone transport structure.

A different multivariate extension is Brenier isotonic regression, where the scalar monotonicity of classical isotonic regression is replaced by cyclic monotonicity for vector-valued outputs. The method uses optimal transport couplings and Brenier potentials to define a multivariate analogue of isotonic regression, with applications to multiclass probability calibration and generalized linear models (Bao et al., 11 Mar 2026). This suggests that isotonic distributional regression can be generalized from stochastic order on scalar distributions to convex-analytic monotonicity in higher-dimensional output spaces.

6. Applications, empirical behavior, and adjacent distributional viewpoints

IDR has been used as a probabilistic postprocessing method for day-ahead electricity prices. In that setting, a point forecast serves as the one-dimensional regressor, and IDR estimates a full predictive distribution under the assumption that larger point forecasts correspond to stochastically larger price outcomes. The method minimizes CRPS under isotonicity, is implemented via an abridged pool-adjacent-violators algorithm together with interpolation, and is used to extract 99 percentiles for operational forecasting (Lipiecki et al., 2024). In the reported experiments, standalone IDR has more variable performance than quantile regression averaging and conformal prediction, but ensemble combinations that include IDR achieve the best overall CRPS and are at par with or better than state-of-the-art distributional deep neural networks over long test periods in German and Spanish electricity markets (Lipiecki et al., 2024). This suggests that IDR contributes complementary distributional information even when it is not the strongest individual forecaster.

Within the DIM framework, IDR has also been applied to intensive-care-unit length-of-stay prediction. A parametric index model is first estimated from clinical covariates, then IDR is applied to the fitted index values to obtain stochastically ordered conditional length-of-stay distributions. In that application, the resulting probabilistic forecasts outperform an empirical CDF benchmark, quantile regression, and Cox regression in mean CRPS across most intensive care units, while retaining good calibration as assessed by PIT histograms and reliability diagrams (Henzi et al., 2020).

There is also a distinct but related sense in which isotonic regression becomes “distributional”: the distributional theory of the estimator itself. For constant underlying sequences and exchangeable noise, the isotonic least-squares estimator can be represented as the slope process of the greatest convex minorant of a random walk, and each coordinate of the estimator has the same distribution as an order statistic of running averages (Soloff et al., 2018). This yields exact finite-sample distribution-free risk formulas such as

$\alpha\in(0,1)$ 5

for exchangeable noise with pairwise correlation $\alpha\in(0,1)$ 6, and $\alpha\in(0,1)$ 7 in the i.i.d. mean-zero, variance-one case (Soloff et al., 2018). That line of work concerns the law of isotonic estimators rather than conditional distribution estimation, but it supplies a complementary distributional perspective on isotonic methods.

Taken together, these strands define isotonic distributional regression as a family of methods in which monotonicity constraints are imposed directly on distributional objects—conditional CDFs, quantile functions, transport maps, or density ratios—and estimated by projection or empirical risk minimization under those constraints. The unifying principle is that stochastic or geometric order replaces parametric specification, while CRPS, Wasserstein, or related distributional losses provide the criterion by which the fitted objects are selected (Henzi et al., 2019).