Risk Bounds For Distributional Regression
Authors: Carlos Misael Madrid Padilla
OSCAR HERNAN MADRID PADILLA
Sabyasachi Chatterjee
Venue: NeurIPS 2025 poster ()
Keywords: Nonparametric distribution estimation
isotonic
trend filtering
dense ReLU networks.
Primary area: theory

Abstract: This work examines risk bounds for nonparametric distributional regression estimators. For convex-constrained distributional regression, general upper bounds are established for the continuous ranked probability score (CRPS) and the worst-case mean squared error (MSE) across the domain. These theoretical results are applied to isotonic and trend filtering distributional regression, yielding convergence rates consistent with those for mean estimation. Furthermore, a general upper bound is derived for distributional regression under non-convex constraints, with a specific application to neural network-based estimators. Comprehensive experiments on both simulated and real data validate the theoretical contributions, demonstrating their practical effectiveness.

Risk Bounds For Distributional Regression
Carlos Misael Madrid Padilla
Oscar Hernan Madrid Padilla
Department Statistics and Data Science Washington University in St Louis
Department of Statisties University of California, Los Angeles
Louis, MO 63130
Los Angele, CA 90095 oscar www.mmddddddaaaaommmmmmmmm ucla. edu
carlosmiiseel@@ulsst. edu
Sabyasachi Chatterjee
Department of Statistics
University of Illinois at Urbana-Champaign
Champaign, IL 61820 @c17060illinoix
Abstract
This work examines risk bounds for nonparametric distributional regression esti- mators. For convex-constraine distributional regression, general upper bounds are established for the continuous ranked probability score (CRPS) and the worst-case mean squared error (MSE) across the domain. These theoretical results are applied to isotonic and trend filtering distributional regression, yielding convergence rates consistent with those for mean estimation. Furthermore, a general upper bound is derived for distributional regression under non-convex constraints, with specific application to neural network based estimators. Comprehensive experiments on both simulated and real data validate the theoretical contributions, demonstrating their practical effectiveness.
Introduction
While regression methods are widely popular across statistics and machine learning, it is well recognized that the conditional mean alone often fails to capture the full relationship between response variable and set of covariates. As noted by Shaked and Shanthikumar [57]. "the ultimate goal of regression analysis is to obtain information about the conditional distribution of response given a set of explanatory variables." A common framework is quantile regression, which estimates
conditional quantiles to provide a more detailed view of the response distribution (35). A more direct approach is distributional regression, which estimates the conditional distribution of the response given the covariates.
Distributional regression has found applications in diverse areas, including electricity spot price analysis 33 understanding income determinants 35 modeling weather data 63 and improving
precipitation forecasts 28 55
The estimation of distributions of random variables under structural constraints is fundamental problem in many statistical and machine learning tasks, including nonparametric regression, density estimation, and probabilistic forecasting [26] 62 23). Consider the sequence distributional model, which we observe independent random variables y1....: Yn E R, each drawn from an unknown distribution Fi..... ......, F", respectively. The objective of this paper is to estimate the vector F* (I) ( for IR. Here, (1) Ply, < () represents the cumulative distribution function (CDF) at a specified for each observation yı.
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

In this paper, we will explore different structural constraints to estimate (t).. These constraints not only ensure interpretable and robust estimators but also prevent overfitting, making them essential in domains such as signal processing, medical diagnostics, and probabilistic weather forecasting. For example, in survival analysis, monotonicity reflects the cumulative nature of survival probabilities, while in genomics, smoothness helps capture gradual trends in gene expression data B 62]. To
rigorously evaluate the quality of the estimators, we employ the continuous ranked probability Score (CRPS), widely used metric for assessing the accuracy of probabilistic forecasts [22]. By quantifying the distance between the estimated and true CDFs, CRPS provides an interpretable and robust framework for comparing estimators under various structural constraints.
1.1 Summary of Results
We now provide a brief summary of the contributions in this paper.
Unified Framework for Estimation: We study unified framework for estimating (* by minimizing a quadratic loss over a convex set K, C R. The convex set K, C R enforces structural constraints on the parameter (t). For the resulting estimator, we provide rigorous theoretical guarantees, including non-asymptotic bounds on both the mean squared error (MSE) and the CRPS.
Applications to Monotonicity and Bounded Variation: We demonstrate the applicability of our
framework to convex constraints arising in monotonicity (9 2 [2] and bounded total variation
These examples illustrate the flexibility and practical utility of the proposed approach for structured parameter estimation.
42
62 23
General Theory for Distributional Regression: We establish a general theory for distributional regression with constraints encoded by arbitrary sets. The main result provides uniform bound the empirical fe error of F* across I, with the upper bound explicitly dependent on the approximation error and the complexity of the sets Kt.
Convergence Rates for Neural Networks: Exploiting the general results for arbitrary sets, we derive convergence rates for distributional regression using dense neural networks. This extends the framework of Kohler and Langer [37 to the context of distributional regression.
1.2 Other Related Work
Distributional regression involves modeling the cumulative distribution function (CDF) of a random variable, whereas quantile regression focuses on estimating the inverse CDF. Koenker et al. [36] provides comprehensive review of the distinctions between these approaches. For instance, the outcome variable is income and the covariate is a binary variable representing educational attainment, distributional regression would model the probability that income falls below certain threshold for each educational group. In contrast, quantile regression would estimate income differences between individuals ranked the same quantile within the two groups. For further discussion on these differences, see also peracchi [50). The approach we adopt in this paper is based on modeling the mean of the random variables of the form 1 (w: which represent the indicator of the events {yr S t} for 1,...,,, n. This idea first introduced in (00 and has since been explored in various contexts, including Firpo and Sunao [19]. Rothe [52] Rothe and Wied (53) and Chernozhukov al. [15] More recently, distributional regression has been studied in diverse settings, such as isotonic regression (28), random forests [10], and neural networks (88 31. In
particular, Henzi [27] develop consistent estimators of conditional CDFs under increasing concave and convex stochastic orders, providing a flexible framework that accommodates heterogeneous variance structures and distributional crossings situations where traditional stochastic dominance assumptions may fail. These works highlight the versatility and applicability of distributional regression across a range of methodologies and problem domains. Finally, other nonparametric approaches to distributional regression include Dunson et al. [7 and Hall et al. 24
2 Notation
Througout, for vector E R", we denote by ||e|| and [| the la and & norms, respectively. Thus, ||a|| = V of and ||o|]] E |vy]- Furthermore, given two functions G, H R + R, define CRPS(G, [C() dt. dt. Also, for for n 0 and E R", we write

---

(2) == fuER" 1 u| S n). For metric space (X, d), let K be subset of A' and r > 0 be a positive number. Let (x, d) be the ball of radius with center X. We say that subset C' C A'' is an r-external covering of K if K C Urec B, (2r, d). Then external covering number of K, written as (r, K d), is defined as the minimum cardinality of any r-external covering of K. Furthermore, for function f X' R. we define its &c norm as |of == sup=ex [f(x)]- Also, for two sequences on and bn we write an X b, if on S cb,, for a positive constant c, and if an & b, and bn X an then we write an x bn. The indicator function of set A is denoted as (t),, which takes value if E A and 0 otherwise. Finally, the indicator an event A is 1A which takes value 1 if A holds and otherwise.
2.1 Outline
The paper is organized as follows: Section 3 establishes a unified framework for distributional regression under structural constraints. In particular, Section 3.T introduces the general methodol- ogy. Section 3.2 derives statistical risk guarantees for the convex case. Sections 3.2.1 and 3.2.2 provide concrete examples of convex estimators, focusing on isotonic regression and trend filtering, respectively. Section 3.3 extends the framework to the non-convex setting. Specifically, Section 3.3.T developp the general theory for non-convex estimators, and Section 3.3.2 explores an estimator present simulation studies and real-data experiments, respectively, comparing the proposed methods to state-of-the-an competitors. Finally, Section concludes with discussion of future research directions, and the Appendix provides additional theoretical results and experimental findings, including two additional real-data applications.
based on deep neural networks. Section [
3 Theory
3.1
General Result for Constrained Estimators
We begin this section by addressing general problems in distributional regression under structural constraints. Using the notation introduced in Section [ our goal to estimate the vector of dis- tribution functions evaluations (t). To achieve this, we adopt empirical risk minimization framework based on the continuous ranked probability score (CRPS), widely used tool for evaluating distributional forecasts
4 22
Specifically, we consider the estimator
F
arg min (A) F()EK for all
CRPS(F, (p.ss))
(I)
where F(t) = (P) (O)) E IR" and K C IR" is a set encoding the structural constraint.
Related CRPS-based formulations have been proposed in the context of isotonic distributional regression (28].
This formulation can be viewed as form of M-estimation over function-valued parameters, where the loss is defined via proper scoring rule. Specifically, the estimator F minimizes the empirical CRPS
loss, computed relative to the observed point masses 1 (....)- Because CRPS is strictly proper [22] the population version of this loss-where expectations minimized at the true distribution functions F.- Hence, the estimator in ( can be interpreted as an justified. The constraint Ftt
taken over the distribution of each y:- y:-is is
empirical risk minimizer, and is statistically and K for all imposes additional structure on the estimator, promoting regularity and interpretability. We now show that the solution the empirical CRPS minimization problem can be obtained via
simple projection estimator. For each E R, define w(t) (1 (t) (mm E IR", and consider the projection
F(t) == arg min (llw(t) -1
(2)
6EK
Lemma 1. For K C R", the function-valued estimator defined by 2 with F((() := (F())), solves solves the empirical CRPS minimization problem in (T).
Thus, for each / E R, the solution F()) in problem 1 is obtained by projecting the empirical vector wi()) onto the convex set K, making the estimation process both computationally efficient and conceptually transparent. In the more general case where the structural constraint may vary with t,
3

---

i.., K2 C R", we extend the estimator to
F(t)
== arg min [||u(() - 0EK,
This projection-based formulation allows to define flexible estimators tailored to varying structural assumptions.
3.2 Risk bounds for convex case
In this subsection, we focus the special case where each constraint set Kr C R is convex and satisfies (t) E Kr- As described in (3), the corresponding estimator F(I) is obtained by projecting the empirical vector vo(() onto K,, which amounts to minimizing the empirical E--loss under convex constraints.
Our next result provides an upper bound on the expected value of the CRPS error. This a consequence of modified version of Theorem A.I in Guntuboyina et al. 23] see Theorem 5 in the Appendix. While our result is conceptually related to Theorem A. in Guntuboyina et al. [23] our setting and proof strategy differ in several important ways; Appendix B for detailed discussion comparing both results.
Theorem 1. Suppose that (0, ...,, 0), (1, ...,,,.................... E Kr, and Ply, E .2) = 1, for all i,
(4)
and for some fixed compact set Ω C R. Then, there exists constant C > O such that
Cn²
(R,s(F,FF))
<
(5)
n2
n
for every >> satisfying
S n/²
sup E tER
sup 110 [10-F*(IISS
g F (t))
(6)
BEK,
where ~ NO, In), universal constant L>> 0.
Thus, Theorem I shows that bounding the expected CRPS error can be reduced analyzing the local Gaussian complexity, as given on the left-hand side of Equation (6).
As direct consequence of Theorem 1 we can derive the same upper bound as in for a rear- rangement (see e.g. Lorentz [39], Bennett and Sharpley [6]) of a truncated version of F, which is non-decreasing by construction. This result is formalized in the following theorem.
Corollary 1. Suppose that (O,..., 0), (1, ...., 1) E K2 for all t E R. and that (-) is continuous for all i. Given i E (1,....n), let F+((1) (1) max (1)) Let 7(1) .... S U() be the order statistics of y. Define ai F (wt) for 1...., 1 and sort the vector a,,
if < y()
ai 2 ... 2 aijn-1s and let F. be defined as Fi(1) =
as Masi-1onne (omm)))))) if y() $ < Uin)
if u(n) St.
where up = H(m) and = V(m) Ek- - voo) for I,..., 1. Then, with the
Cn2
notation from Theorem I we have that
E
CRPS(F.,F)) F7)
The function (t) defined above non-decreasing and is constructed by modifying the original estimator ((( through rearrangement of values. Specifically, this is achieved by applying a change of variable to (t), followed by the Hardy-Littlewood decreasing rearrangement. This classical construction 25 preserves the level set measures and minimizes the L distance among all decreasing equimeasurable functions. An example of F.(1) and (t(() is given in Figure 2 in the Appendix. We now present the final result from this section.

---

Theorem 2. Suppose that K₁ C. K for all t. For any 17 > 0 it holds that
P
sup Σ
(F()) 2n²
1 1p J.
C
Cylogn
for positive constant C.
Theorem 2 provides a high-probability concentration bound for the MSE in estimating F* holding uniformly over all R. This sets it apart from standard sub-Gaussian bounds, which typically yield control at a fixed evaluation point f, see for example [4]. To ensure that the bound in Theorem 2 small, it suffices to upper bound the local entropy of K K, where K an upper set that contains the sets Kr.
3.2.1 Isotonic Regression
In this subsection, we present the first application of our general theory from Section focusing on distributional isotonic regression-a topic that has garnered significant attention in the literature
The most relevant works to our results are Mösching and Dümbgen AA which examined distributional isotonic regression under smoothness constraints, and Henzi et al. (28). which proposed an interpolation method equivalent the formulation in with K, enforcing monotonicity constraint.
16
18
29
32
Setting
(8) we now consider the case of isotonic distributional regression assuming that (* E K,, which equivalent to P(Y1 S t) I < P(Y, I t). With the constraint sets as in (8). the resulting estimator
K2 =
in (3) can be found with the pool adjacent violators algorithm from Robertson (51.
Corollary 2. Consider the estimators {Fe [FiheR defined in
3
with K, defined as 8 (f) EK
and 44) holds for some fixed compact set Q, then
CRPS(P,FF) F.) S Cn 2/3
E
(9)
for some positive constant C > O. Moreover, (9) holds replacing F with the corresponding F as defined in Corollary 7
provided that each function F. continuous. Finally,
(P:(0) (P()))(())) ==:00=== = Op 2/33
log
sup
(10)
EER
n
where 70p holds without requiring nor continuity of the Fi's.
The result in Corollary Sacallaryessae establishes that distributional isotonic regression achieves an estimation rate of n 2/3 for both the expected average CRPS and the worst-case MSE, as shown in (10). This result improves upon Theorem 3 in Henzi al. 28 in the univariate case, which only demonstrated convergence in probability for isotonic distributional regression. However, we emphasize that Henzi al. 28 study the more general setting of multivariate covariates, while our analysis restricted to the univariate case (d 1).
To further contextualize theoretical guarantees, compare them with Theorem 3.3 in Mösching and Dümbgen [4]] which establishes uniform consistency for estimating the conditional distribution function. Considering the fixed design case in their formulation and assuming without loss of generality that in their notation X1 < ... < Xx, E R. the goal is to estimate an unknown family of distributions (F2 (FZ))ERR where for each fixed E H. the map x Ft(t) is assumed be non- decreasing and a- a-Hölder continuous with constant C > O that is the same across t. Additionally, the design assumed be asymptotically dense- i.e., the covariate values X, < < X,, sufficiently cover the domain. Translating their setup into our notation, their target Fx (t) == P(Y; < X;) corresponds to F (4) := Ply: < I) in our sequence model. Their assumption that x + (2()) is non- decreasing for each aligns with our isotonic regression framework, where we impose monotonicity of the sequence (1)
< F. (f). However, in contrast to their setting, we do not require any smoothness assumptions such Hölder continuity the sequence {F) [t) nor do require
5

---

the covariates to be dense. Despite that, our method achieves faster convergence rate of order 2/3 for both the average CRPS risk and the worst-case MSE, compared to their rate n 222//(2a+)) (up logarithmic factors) for the worst-case MSE, with a (0, 1]. established under their stronger regularity assumptions. A summary of these comparisons is provided in Table 2 Appendix
We also show in Appendix that faster rates achievable under additional structural assumptions on F (t). Specifically, if F [t) has few strict e.c.. if it is piecewise constant with a small number of jumps then the estimator can attain nearly parametric risk rates up to logarithmic factors.
3.2.2 Trend Filtering
In this subsection, we apply the theory from Section 3 to distributional regression under a total variation constraint. Total variation- methods were independently introduced by Rudin et al. [54], Mammen and Van De Geer [42], and Tibshirani et al. [6]. These methods have been extensively
studied in various contexts within the statistics literature, including univariate settings [62 38 23
46), grid graphs 30 and general graphs [65 48
Before establishing our proposed total variation estimators, we introduce some additional notation. For a vector 0 E IR",, define (( (0) = 0, (0) = (02 01..... On - 0n- and ()), for
r 2 2, is recursively defined as ( = D (D( where ()) E RR-r. With this notation, for > 1, the rth order total variation of vector is given as
(V) n'''n'"" (
(11)
The concept of the rth total variation can be understood as follows. Consider the evaluations of an r times differentiable function f [0, 1] - I on the grid (1/n, 2/n,..., n//.. In this case, a Riemann approximation of the integral Sp... (t) corresponds precisely to TV(* (0). where
f(r) denotes the rth derivative of f. Therefore, for natural instances of 0, it is reasonable to expect that TV(*) (0) O()). The above discussion motivates us to define the sets
K2 = E R." TV v.].
(12)
for some V; > O, and consider the corresponding estimator in This estimator is referred to as trend filtering following standard terminology in the literature (see, for example, [65] 62] 61). where methods based on bounded total variation of order are collectively known as trend filtering estimators. The intuition here is that if (1) E K, then the probabilities F (L),..., Fx (1) change smoothly over i in the sense that F (1) has bounded rth total variation. We also highlight that the set K, defined in Equation 12 is convex; see Appendix [G. a formal proof. The resulting set in (12) allows us to define the trend filtering distributional regression estimator subject of our next corollary which follows from the results Section 3.1
Refined risk bounds under additional sparsity assumptions are presented in Appendix A where we show that trend filtering estimators can achieve near near-parame ear-paramettic rates when the signal both smooth and piecewise sparse. These results extend recent adaptive risk bounds in trend filtering, see, for example, Guntuboyina al. 233
Corollary 3. Consider the estimator in 3 with K, as in (72) for an integer r satisfying r 2 1. If (I) E Kg 44 holds for some fixed compact set Q and sup, V S V, then
a
(T)
-
log
E
S
(13)
n
for positive constant C. Moreover, the upper bound in (13) also holds for the corresponding sorted estimators F. as defined in Corollary 1
if in addition each function Př is continuous. Finally,
V
log
sup
n (Fi(t) (F))) = -0
+
(14)
ER
n2 2r++
n
where (14g holds without requiring nor continuity of the F. 's.

---

Corollary 3 establishes that the constrained version of trend filtering for distributional regression achieves the rate / ignoring logarithmic factors, for both the CRPS and the worst-case MSE. This result aligns with the convergence rate trend filtering one-dimensional regression, where the same rate is attained when the regression function has rth-order total variation summary of these comparisons is provided in Table 2 in Appendix E Additionally, per Corollary in Appendix B the penalized version of trend filtering for distributional regression achieves the same rate in terms of the worst-case MSE, further reinforeing its consistency with classical trend filtering results.
42 62 23
3.3
Risk Bounds for the General Case
3.3.1
General Result
This subsection aims to present our main result on constrained distributional regression in scenarios where the constraint sets K, are arbitrary, not necessarily convex, and potentially misspecified for (t).
Theorem Let F() be the estimator defined in (3) for all I but with K. not necessarily convex and with F (1) not necessarily in Kr. Suppose that sup sup [IF()) S B for some constant
SERF(t)EK:
B 1, and Kr K for all and some set K. Let G(t) be defined as G(l) arg min [F(() F(jEK,
F' (t)]loo.. Then, for v) > 1, with K(n) = (K- - K] n (0),, we have that
PP
sup EER
[I(()) - F+]I + (t supvnl
G(
tER
2//2
Cylogn
Cyn
Vlog ((e, (()), || I1)de +
sup ||G(t) (t)ll)
20-2
1
rER
for some positive constant C > 0, and where J Ban(2nB//')
Tog
The intuition behind Theorem iss that captures the estimation error, which depends on the local covering complexity of the sets (v)) The second term, suptea y [ ( corresponds to the approximation error, measuring how well the true F* (t) can be approximated within the model class K,. While such a decomposition into estimation and approximation error standard in nonparametric theory, our result provides a uniform guarantee over all E R. in contrast to classical bounds that control error only at a fixed evaluation point, in 47. This distinction particularly relevant for distributional regression problems, where the goal is to control the entire CDF path. The proof of Theorem 3 provided Appendix G.8 relies on peeling argument and extends techniques originally developed for the convex case.
3.3.2 Dense ReLU Networks
We now turn to the application of Theorem 3 to the problem of distributional regression using dense neural networks. The results in this section add to the literature statistical theory for rectified linear unit (ReLU) networks as in Bauer and Kohler [4]. Schmidt-Hieber (56), Kohler and Langer 37], Padilla al. [49, Ma and Safikhani (40) Zhang et al. (67). Padilla al. [47].
Before presenting our main result, we first introduce some notation. Suppose we are given i.i.d. data { y.) [0, 1]doo x R and let
(x,,, == P(y: S fl., x) fi
(15)
for functions f (0, 1] + R for all f. We set the conditional edf (t)
To define the constraint set Kr, we follow Kohler and Langer [37] and assume all hidden layers have the same width. Let F(L, v) denote the set of neural networks with depth L, width and ReLU activation, restricting the functions in F(L, v) to satisfy || < 1. The structure of these networks is described in Appendix C Then, for all E R, define
K₁ = K = 10 R" 0, f(x.), i= 1.... n, for some f F(L,
(16)
7

---

We next summarize the main assumptions underlying the analysis, stated formally in Appendix C Specifically, for each threshold we assume that the true conditional distribution function G (- t) belongs to a hierarchical composition model class HHl, PP) as defined in Kohler and Langer [37 Here, I denotes the number of hierarchical composition levels, and P C (0, oo) x N specifies allowable pairs (p, M) corresponding to the smoothness p and input dimension M of each component function. These assumptions, ensure that the neural network class F(L, v/) can approximate G* (., t) at the
rate @n, = max(p.M n M2) which determines the convergence behavior established the result below.
We now present our main result concerning distributional regression with ReLU neural networks.
Corollary 4. Let F(t) be the estimator from (3) with the set Kr as in 16) for all t E IR with F (t) not necessarily in K₁- Suppose that Assumption [
described in Appendix C holds. Let
6n maX M)P n Under the choices of L and specified in equations (51 or (2) in Appendix G.9 the estimator satisfies
1 (R(0) (Z()) = =Op (logn (1ogm + n
sup CER
Φn log' n)
(17)
n
Corollary 4 demonstrates that dense ReLU neural network estimators for distributional regression uniformly achieve the rate o... up to logarithmic factors, in terms of the worst-case MSE for estimating the true parameters (F* (F*)t) (t)), provided these parameters belong to a hierarchical composition class. Importantly, while this rate matches that of Kohler and Langer (37 for mean regression under sub-Gaussian error assumptions, thresholds / E R, rather than for a fixed target. A summary of this comparison is provided in Table 2 in Appendix E This extension essential for distributional learning tasks where uniform control required, such as CRPS-based risk bounds.
result strengthens the guarantee by holding uniformly over all
4 Simulated data analysis
We evaluate the performance of the proposed methods against state-of-the-art approaches across diverse simulation settings that reflect various practical challenges and structural assumptions. Specif- ically, six distinct scenarios are considered to evaluate different aspects of the distributional regression problem. The implementation of all experiments is available /https://// : / thub. com/ coom/cmadrddp. UnifDR We refer to our proposed approach as UnifDR which adapts different methods based on the scenario. In the first two scenarios, UnifDR applies the isotonic regression method from Section in the next two, it uses the trend filtering approach from Section 3.2.2, and in the final two scenarios, employs the Dense ReLU Networks method described in Section 3.3.2
Scenario 1 (S1). We generate data y ~ (44. 1) where 14; = 1 i/n for i = 1,..., n.
Scenario 2 (S2). We consider y Unif(a;, b,), where = (n- (0--1) in/n and by = + 1.
Scenario 3 (S3). The true CDFs are modeled F (f) exp( (op(-t/mus), uz + 0.5 (n.5sin(2m1/11
Scenario (S4). Consider (t) = Gamma(shape 0.7, scale = 14), where 14 = 6 1 ((15n/1) +
2 /////4 + 8 [n/22<153/44 + 4 1ii33n/4)-
Scenario 5 (S5). Let x ~ Unif( (0., 1 The true CDFs are given by (t) 2(( h(x;) 1)/0.))
where h(x;)) 3r. (1) + log(1 + x (2)) + x²{²) + 5x(()) + (a (o) with (1(1) denoting the j th
coordinate of z,. The function & represents the standard normal cumulative distribution function.
Scenario 6 (86). Let x; [ni 1j10) and y/
x² (x:)) Here, h(xx))
log -0. 0.5 Σ sin(Tx.))) (3) 0.5 E - 6) + 0.5 cos(2 (o))) +2). with a!!) denoting the
th coordinate z;-
To implement our proposed approach UnifDR we proceed as follows. The isotonic method introduced in Section 3.2.1 I is implemented in using the pool adjacent violators algorithm (PAVA) from Robertson [51]. For this approach there are no direct competitors in the distributional regression problem. The Trend Filtering estimator in Section 3.2.2 is implemented using the trendfilter function from the glmgen package in R, and compare it with additive smoothing splines (AddSS)
8

---

CRPS results for Scenarios S3 and 54 in A2
CRPS results for Scenario 55 in A2
0.04
Competitors
0.05
S3-UnifR
0.03
54-UnifDR
53 AddSS
0.04
DRT MARS
54 AddSS
CRP5
0.0.0
0.00
0.02
0.01
0.01
0.00
n=000
n=800
n=1600
n=400
n=800
n=1600
Sample Size
Sample Size
Figure 1: Box plots of CRPS results in Ag. The left plot corresponds to S3 and S4, while the right plot displays the results for S5. Table Evaluation metries for UnifDR (Trend Filtering approach) and its competitor AddSS on the 2015 Chicago crime dataset.
METHOD CRPS (MEAN + STD)
MSD (MEAN + STD)
UNIFDR
0.097610.0017
0.2509-0.000 0.2850=0.0391
ADDSS 0.1223=0.0078
via the smooth. spline function in R. For the Dense ReLU Networks method in Section 3.3.2 use fully connected feedforward architecture with input layer, two hidden layers (64 units each), and an output layer. The network is implemented in Python and trained using the Adam optimizer with learning rate of 0.001. In this case, compare the proposed UnifDR method with five benchmark methods. First, we consider Classification and Regression Trees (CART) (8) implemented in R via the rpart package, with the complexity parameter used for tuning. Second, we evaluate Multivariate Adaptive Regression Splines (MARS) (21]. available in the earth package, where the penalty parameter serves as the tuning parameter. Third, we assess Random Forests (RF) 77). implemented in R. via the randomFores orest package, using 500 trees and tuning the minimum terminal node size. Additionally, we consider two recent methods. Distributional Random Forests (DRF) [10 is implemented via the drf package in Python. DRF employs tree-based ensemble with tuning parameters including the splitting rule and the number of trees. Lastly, Engression (EnG) [58] is implemented using the engression package in Python. EnG utilizes hierarchical structured neural networks, and we adopt the same training hyperparameters as our deep learning approach to ensure optimization consistency. EnG also requires sampling procedure, with the number of samples set to 1000 for accurate distribution estimation.
Performance Evaluation: For each scenario, datasets with sample sizes E (400, 800, 1600) are generated, with each experiment repeated 100 times using Monte Carlo simulations. Evaluations conducted at 100 evenly spaced points t from three fixed intervals: A = [-1, 0.4], A2 = 2- 2]. and Ag = [0.8, 10). Each dataset is randomly split into 75% training and 25% test sets. Competing models undergo f-fold cross-validation on the training data for hyperparameter tuning, with performance assessed on the test set. For the isotonic regression method, test set predictions are obtained via naive nearest neighbor interpolation. The accuracy of the estimated CDFs (t()) relative to the true CDFs F (t) is evaluated using the following performance metrics, averaged over 100 Monte Carlo repetitions. CRPS: CRPS evaluates the overall fit of (t) to F across all evaluation points Since the evaluations in experiments are performed over finite set of 100
in IR, see Section 2
values in A, CRPS is approximated via a Riemann sum: CRPS
Σie
100 EtEA (f((t)
Test
(t)) where [est]] = 4 the size of the test set. Maximum Squared Difference (MSD): MSD
captures the worst-case discrepancy between Ft(t) and F. (t), and is approximated as: MSD
max
Teas
EiETesi
100-Faa (z)
The results below focus on the CRPS metric for Scenarios S3, S4, and 85, where the CDFs are evaluated the points in Ag. Additional results, including MSD performance and CRPS evaluations at A₁ and Ag across all scenarios. as well as the A₁ -bssed CRPS for Scenarios S1. S2, and S6, are provided in Appendix Figure I presentss the performance of UnifDR S3 and S4, where we compare the Trend Filtering approach with AddSS. The same figure also includes results for S5, where UnifDR utilizes the Dense ReLU Networks method against five state-of-the-art competitors: CART, MARS, RF, DRE, and EnG. In all scenarios, UnifDR consistently outperforms competing

---

methods, with its performance superiority becoming more pronounced as the sample size increases. This dominance further confirmed by the extended evaluations in Appendix [F] reinforcing the robustness of UnifDR across various conditions.
5 Real data application
In this section, we evaluate the performance of the proposed UnifDR method using both the Trend Filtering and Dense ReLU network procedures on real-world datasets.
5.1 Chicago crime data
We analyze the 2015 Chicago crime dataset, available at lttps:://ata. goveoopn--oooopp--ao-ooeoo gov/ which records reported crimes in Chicago throughout the year. Following Tansey et al. [60], the spatial domain is discretized into a 100 x 100 grid, where each grid cell aggregates crime counts within its spatial boundary. The response variable is defined as the log-transformed total crime counts per grid cell. Grid cells with zero observed crimes are excluded, yielding a final dataset of 3,844 grid cells. The UniIDR' Trend Filtering procedure does not covariates. Instead, the spatial grid is treated as an ordered sequence, assuming a smooth spatial trend. The grid cells are ordered lexicographically. The dataset is randomly partitioned into 100 train (75%)-test (25%) splits, and evaluation is conducted at evenly spaced points A in the interval [-1,6]. Performance is assessed using the Continuous Ranked Probability Score (CRPS) and Maximum Squared Difference (MSD) metrics, comparing estimated CDFs (t() against empirical indicators (t). where & A. AddSS is used as benchmark for the Trend Filtering approach on the same dataset. Table [ presents CRPS and MSD metries, demonstrating UnifDR superior performance. Furthermore, Figure 4 in Appendix provides a visualization of (1) at L = 3 for both competing methods. The same data set ondlyzed using UnifDR with the Dense ReLU network framework, with further details in Appendix
F111
5.2 Other real data examples
Beyond the Chicago crime dataset, UnifDR is further evaluated on California housing prices and daily Ozone measurements. Detailed descriptions, pre-processing steps, and results are provided in the Appendix F.
6 Conclusion
This paper introduced a unified framework for nonparametric distributional regression under convex and non-convex structural constraints. We established theoretical risk bounds for the estimation of cumulative distribution functions (CDFs) in various settings, including isotonic regression, trend filtering, and deep neural networks. Our analysis leveraged continuous ranked probability scores (CRPS) and worst-case mean squared error (MSE) to quantify estimation accuracy, demonstrating that structured constraints such as monotonicity, bounded total variation, and hierarchical function composition lead to improved estimation accuracy. The resulting upper bounds are general and apply broad class of convex and non-convex estimators. While we do not derive lower bounds at this level of generality, the rates obtained for specific cases such as isotonic and trend filtering estimators-mmn match known minimax results for mean estimation, providing evidence for the tightness and near-optimality of our bounds in these canonical examples (see Table 2 In particular, we establish explicit convergence guarantees for isotonic, strengthening prior findings. Moreover, the trend filtering estimator achieves rates consistent with classical one-dimensional regression results. For deep neural networks, we show that our estimator achieve comparable rates under hierarchical composition constraints, aligning with existing results in structured regression. Experiments on simulated and real datasets further validated the theoretical guarantees, with our proposed methods consistently outperforming alternative approaches. In particular, UnifDR, the distributional regression framework we study, demonstrated superior performance across all considered settings.
An important avenue for future research is extending our theory to dependent data settings. Many real-world applications involve time-series data, spatial data, or sttorr--structured data, where dependencies among observations must be accounted for.
10

---

References
[1] Dennis Amelunxen, Martin Lotz, Michael B McCoy, and Joel Tropp. Living on the edge: A geometric theory of phase transitions in convex optimization. 2013.
[2] Richard E Barlow and Hugh D Brunk. The isotonic regression problem and its dual. Journal of the American Statistical Association, 37(337):140-1444 1972.
[3] Richard E. Barlow, David Bartholomew, John Bremmer, and H Dale Brunk. Statistical inference under order restrictions: The theory and application of isotonic regression. Wiley, 1972.
[4] Benedikt Bauer and Michael Kohler. On deep learning remedy for the curse of dimension- ality in nonparametric regression. 2019.
[5] Pierre C Bellec, Sharp oracle inequalities for least squares estimators in shape restricted regression. The Annals of Statisfics, 46(2):745-780, 2018.
[6] Colin Bennett and Robert Sharpley. Interpolation of operators. Academic press, 1988.
[7] Leo Breiman. Random forests. Machine learning, 45:5-33 32, 2001.
[8] Leo Breiman, Jerome Friedman, and Charles Stone. e richard olshen, 1984.
[9] Hugh Brunk. Estimation of isotonic regression. University of Missouri-Columne 1969.
[10] Domagoj Cevid, Loris Michel, Jeffrey Näf, Peter Bühlmann, and Nicolai Meinshausen. Distri- butional random forests: Heterogeneity adjustment and multivariate distributional regression. Journal of Machine Learning Research, 23(333):1-79, 2022.
[11] Sabyasachi Chatterjee and Subhajit Goswami. New risk bounds for 2d total variation denoising. IEEE Transactions on Information Theory, 67(6):4060-44 2021.
[12] Sabyasachi Chatterjee, Adityanand Guntuboyina, and Bodhisattva Sen. On risk bounds in isotonic and other shape restricted regression problems. 2015.
[13] Sourav Chatterjee. A new perspective least squares under convex constraint. 2014.
[14] Sourav Chatterjee. Matrix estimation by universal singular value thresholding. 2015.
[15] Victor Chernozhukov, Iván Fernández-Vän Val, and Blaise Melly, Inference on counterfactual distributions. Econometrica, 81(6):2205-222) 2013.
[16] Ori Davidov and George Iliopoulos. Estimating a distribution function subject to a stochastic order restriction: comparative study. Journal of Nonparametric Statistics, 24(4):923-933, 2012.
[17] David B Dunson, Natesh Pillai, and Ju-Hyun Park. Bayesian density regression. Journal of the Royal Statistical Society Series B: Statistical Methodology, 69(2):163- 183, 2007.
[18] Hammou El Barmi and Hari Mukerjee. Inferences under stochastic ordering constraint: the k-sample case. Journal of the American Statistical Association, 100(469);252-261, 2005.
[19] Sergio Firpo and Stefanie Sunao. Decomposition methods. In Handbook of Labor, Human Resources and Population Economics, pages 1-32. Springer, 2023.
[20] Silverio Foresi and Franco Peracchi. The conditional distribution of excess returns: An empirical analysis. Journal of the American Statistical Association, 90(430):451-466, 1995.
[21] Jerome H Friedman. Multivariate adaptive regression splines. The annals of statistics, 19(1): 1-67, 1991.
[22] Tilmann Gneiting and Adrian E. Raftery. Strictly proper scoring rules, prediction, and estimation. Journal of the American Statistical Association, 102(477):35-377 2007.

---

[23] Adityanand Guntuboyina, Donovan Lieu, Sabyasachi Chatterjee, and Bodhisattva Sen. Adaptive risk bounds in univariate total variation denoising and trend filtering. The Annals of Statistics, 48:205-229, 2020.
[24] Peter Hall, Rodney CL Wolff, and Qiwei Yao. Methods for estimating conditional distribution function. Journal of the American Statistical association, 94(445):154 1999.
[25] Godfrey Harold Hardy and John Edensor Littewood. Some properties of fractional integrals. i. Mathematische Zeitschrift, 27(1):6656666, 1928.
[26] Trevor Hastie, Robert Tibshirani, and Jerome Friedman, The Elements Staristical Learning: Data Mining, Inference, and Prediction. Springer, New York, 2nd edition, 2009. ISBN 778-0387848570.
[27] Alexander Henzi. Consistent estimation of distribution functions under increasing concave and convex stochastic ordering. Journal of Business & Economic Statistics, (((4):1203-124,, 2023.
[28] Alexander Henzi, Johanna Ziegel, and Tilmann Gneiting. Isotonic distributional regression. Journal of the Royal Statistical Society Series B: Statistical Methodology, 33(5):963-992) 2021.
[29] Robert V Hogg. On models and hypotheses with restricted alternatives. Journal of the American Statistical Association, 60(312):1 11T3::113-116 1965.
[30] Aan-Christian Hütter and Philippe Rigollet. Optimal rates for total variation denoising. In Conference on Learning Theory, pages 1115-1146, PMLR, 2016.
[31] Ehsan Imani and Martha White. Improving regression performance with distributional losses. In International conference on machine learning, pages 2157-2166. PMLR, 2018.
[32] Javier Rojo Jiménez and Hammou El Barmi. Estimation of distribution functions under second order stochastic dominance. Statistica Sinica, pages 903-926, 2003.
[33] Nadja Klein. Distributional regression for data analysis. Annual Review of Statistics and Its Application, 11, 2024.
[34] Thomas Kneib, Alexander Silbersdorff, and Benjamin Säfken. Rage against the mean-a review of distributional regression approaches. Econometrics and Statistics, 26:99-123, 2023.
[35] Roger Koenker and Gilbert Bassett Jr. Regression quantiles. Econometrica: journal of the Econometric Society, pages 33-50, 1978.
[36] Roger Koenker, Samantha Leorato, and Franco Peracchi. Distributional vs. quantile regression. 2013.
[37] Michael Kohler and Sophie Langer. On the rate of convergence of fully connected deep neural network regression estimates. The Annals of Statistics, 49(4):2231-2244 2021.
(38) Kevin Lin, James L. Sharpnack, Alessandro Rinaldo, and Ryan J Tibshirani. A sharp error analysis for the fused lasso, with application approximate changepoint screening. Advances in neural information processing systems, 30, 2017.
[39] Georg Gunther Lorentz. An inequality for rearrangements. The American Mathematical Monthly, 60(3):176-179, 1953.
[40] Mingliang Ma and Abolfazl Safikhani. Theoretical analysis of deep neural networks for temporally dependent observations. Advances in Neural Information Processing Systems, 35: 37324-37334, 2022.
[41] Oscar Herman Madrid Padilla and Sabyasachi Chatterjee. Risk bounds for quantile trend filtering. Biometrika, 109(3)7511-66, 2022.
[42] Enno Mammen and Sara Van De Geer. Locally adaptive regression splines. The Annals of Statistics, 25(1):387-413, 1997.
12

---

[43] James E Matheson and Robert L Winkler. Scoring rules for continuous probability distributions. Management science, 228100:0077-106, 1976.
[44] Alexandre Mösching and Lutz Dümbgen. Monotone least squares and isotonic quantiles. 2020.
[45] Kelley Pace and Ronald Barry. Sparse spatial autoregressions. Statistics Probability Letters, 33(3):291-297, 1997.
[46] Carlos Misael Madrid Padilla, Oscar Hernan Madrid Padilla, and Daren Wang. Temporal-spatial model via trend filtering. arXiv preprint arXiv:2308 161722223...77,,,... 2023.
[47] Carlos Misael Madrid Padilla, Oscar Hernan Madrid Padilla, Yik Lun Kei, Zhi Zhang, and Yanzhen Chen. Confidence interval construction and conditional variance estimation with dense relu networks. arXiv preprint arXiv:2412.20335,, 2024.
[48] Oscar Hernan Madrid Padilla, James Sharpnack, James G Scott, and Ryan Tibshirani. The dis fused lasso: Linear-time denoising over general graphs. Journal of Machine Learning Research, 18(176):1-36, 2018.
[49] Oscar Hernan Madrid Padilla, Wesley Tansey, and Yanzhen Chen. Quantile regression with relu networks: Estimators and minimax rates. Journal of Machine Learning Research, 23(247): 1-42, 2022.
[50] Franco Peracchi. On estimating conditional quantiles and distribution functions. Computational statistics & data analysis, 98(4):433-444) 2002.
[51] Tim Robertson. Order restricted statistical inference/(et al.). 1988.
[52] Christoph Rothe. Partial distributional policy effects. Econometrica, 80(5):229--2301 2012.
[53] Christoph Rothe and Dominik Wied. Misspecification testing in a class of conditional distribu- tional models. Journal of the American Statistical Association, 108(501):314-32) 2013.
[54] Leonid Rudin, Stanley Osher, and Emad Fatemi. Nonlinear total variation based noise removal algorithms. Physica D: nonlinear phenomena, 60(1-4):259-25 1992.
[55] Lisa Schlosser, Torsten Hothorn, Reto Stauffer, and Achim Zeileis. Distributional regression forests for probabilistic precipitation forecasting in complex terrain. 2019.
[56] Johannes Hiiee.. Nonparametric regression using deep neural networks with relu activation function. 2020.
[57] Moshe Shaked and J George Shanthikumar Stochastic orders. Springer, 2007.
[58] Xinwei Shen and Nicolai Meinshausen. Engression: extrapolation through the lens of distribu- tional regression. Journal of the Royal Statistical Society Series B: Statistical Methodology, page gkac 108, 2024.
[59] Michel Talagrand. The generic chaining: upper and lower bounds of stochastic processes. Springer Science & Business Media, 2005.
[60] Wesley Tansey, Jesse Thomason, and James Scott. Maximum-variance total variation denoising for interpretable spatial smoothing. Proceedings of the AAAI Conference on Artificial Intelligence, volume 32, 2018.
[61] Robert Tibshirani, Michael Saunders, Saharon Rosset, Ji Zhu, and Keith Knight. Sparsity and smoothness via the fused lasso. Journal of the Royal Statistical Society: Series (Statistical Methodology), 67(1): 111)8991108.... 2005.
[62] Ryan J Tibshirani. Adaptive piecewise polynomial estimation via trend filtering. 2014.
[63] Nikolaus Umlauf and Thomas Kneib. A primer bayesian distributional regression. Statistical Modelling, 18:2-4) 2018.
[64] Martin J Wainwright. High-dimensional statistics: non-asymptotic viewpoint, volume 48. Cambridge university press, 2019.
13

---

[65] Yu- Yu-Xiang Wang, James Sharpnack, Alexander J Smola, and Ryan J Tibshirani. Trend filtering on graphs. Journal of Machine Learning Research, 17(105):1-41, 2016.
[66] Steven Siwei Ye and Oscar Hernan Madrid Padilla. Non-parametric quantile regression via k-nn fused lasso. Journal of Machine Learning Reseurch, 22(111):1-333 2021.
[67] Zhi Zhang, Carlos Misael Madrid Padilla, Xiaokai Luo, Oscar Hernan Madrid Padilla, and Daren Wang. Dense relu neural networks for cemporal-spatial model. arXiv preprint 0rXiv:2411.09961 2024.
14

---

NeurIPS Paper Checklist
Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
See Section S and 3.3
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results be expected to generalize other settings.
is fine include aspirational goals as motivation as long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: See Section 0
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification. asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low images are taken in low lighting. Or speech-to-text system might be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach address problems of privacy and faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question:
For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
15

---

Justification: See Section and 3.3 for the results with assumptions, and Appendix G for the corresponding proofs.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proof's can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: See section 4 Section 5 and Appendix
Guidelines:
The answer NA means that the paper does include experiments. If the paper includes experiments, No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm. (b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome deseribe the particular way they provide for reproducibility, In the case of closed-source models, it may be that access to the model is limited in some way (e.g., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
16

---

reproducible, they Answer: [Yes]
Justification: See the supplementary material. It includes all code and implementation details for the experiments. We also include the Real data URL to access our examples in Section 5 and Appendix
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (kttps://nips.co cc/
public/guiden htwps///ooeSuuuusiiinnneeeee
for more details.
While encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https
7 7/nips. //nips.cc/pub. httppss.ccppbblc/guuueeee wttpeee/s/oooSSeeesseeiiiiii
ssioiiolioiiioiocyicc for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (c.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Answer: [Yes]
Justification: See section 4 Section 5 and Appendix
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: All our experiments include box plots or tables with standard errors reported. See Figures 20 and tables 7
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results accompanied by error bars, confi- dence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
17

---

The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: See Section f and Section 5
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental rums as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code of ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics httpp/////.aaopps.cccocc//////ia
Answer: [Yes]
Justification:
We follow the Code of Ethics stablished here https://neurips.ca
Mublic/Ethicsuuider
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., if there special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
This paper presents work whose goal is to advance the field of Machine Learning. There many potential societal consequences of work, none which we feel must be specifically highlighted here
Justification:
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
18

---

Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), faireess considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper deseribe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
No risk.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [NA]
Justification: not use existing assets.
Guidelines:
The answer NA means that the paper does use existing assets. The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website). the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms use in package should be provided. For popular datasets, paperswithcode com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
19

---

For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification:
The paper does not release new assets.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moddl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: paper does not involve any crowdsourcing or research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Bthies, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: paper does not involve research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, not include any information that would break anonymity (if applicable), such the institution conducting the review.
20

---

16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology. scientific rigorousness, originality of the research, declaration is not required.
Answer: [NA]
Justification: LLMs were not used as an original or non-standard component of the core methods in this research.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, or non-standarr components.
Please refer to our LLM policy ftttps://neuriiee /t//Cooffeeeccss/22222 for what should or should be described.
21

---

A
Fast Rates
This appendix develops refined risk bounds for distributional regression estimators under additional structural assumptions on the true signal. While our general theory establishes minimax-optinn rates under convex constraints such as monotonicity or bounded total variation, certain low-complexity signal classes allow for significantly faster convergence.
We provide two canonical examples illustrating this phenomenon. The first focuses isotonic regression, where we show that if the true distribution functions have only a small number of strict increases, the estimator achieves nearly parametric rate. The second concerns trend filtering, where we demonstrate that sparsity in the higher-order differences of the signal-combined with a minimum segment length condition leads to similarly fast convergence.
These results extend analogous adaptive risk bounds from classical point estimation 12 to the distributional regression setting.
A.1 Isotonic Regression
We begin by revisiting the isotonic case. In the main text, C general risk bound of order 2-3 under minimal monotonicity assumptions. Here, we refine this analysis by showing that much faster rates, nearly parametric, are achievable when the true distribution function F* exhibits low complexity, in the sense of having few strict increases.
established
To formalize this, consider the isotonic constraint set
K (0 IR" :1 S S..
(18)
and assume (t) E K for all R. Define
k(t) = | (1..... 1] FT F. (t)
This quantity counts the number of strict increases in (t), and hence controls the complexity of the signal.
Corollary Let F. be the estimator defined in with K2 = K as in 18 Assume F (* E and that 4 holds for a compact set Qt Then
CRPS(F,,FF)
1+ ++t(t)
en
E
Csup YER
log
1 k1))(t))
n 1=1
for some constant C > O. Moreover, this bound also applies to the sorted estimator F from Corollary D provided each F. is continuous.
This result improves on the general n 2/3 bound when the number of strict increases k(t) is small. For example, if F* (t) is piecewise constant with at most jumps, then k(t) S and the bound becomes nearly achieving the parametric rate when = 0(1). The proof of Corollary is
1+8 log
provided in Appendix G.10
A.2 Trend Filtering
We now turn to the trend filtering setting. In the main text, Corollary Sacollaryyestallliin esttablise general risk bound of order n (rr((++))) under a total variation constraint of order r, matching the minimax-optimal rate for function estimation under such constraints 42 62
Here, we show that significantly faster rates are achievable in the distributional regression setting when the signal exhibits additional structure. In particular, when the rth-order differences of the vector F* (t) are sparse and well-separated, the estimator can attain nearly parametric accuracy up to logarithmic factors.
Define
K, { (00
(19)
where VI" (F*)(F(()).
22

---

Corollary 6. Consider the estimator in (3) with Ke as in (19 for an integer satisfying 2 1. Suppose s = = (1)llo and - { (D()) F (1), + 0 for all t. Let ju < J1 <... < js+1
r and Js the elements of S. With this notation define Tho Mi++11 0. Then for j E S define ni to be if (Di (t)) (Di--)) (())) (t))
be such that jo 1, jst1
otherwise set πΙ = 1. Suppose that (t) satisfies the following minimum length assumption
en
min [E(s] (14))) #9/t11
(141 - ji) 2
8 8+1
for some constant e. satisfying 0 S cs 1. Then, for sup, V? S VV'
E Σ CRPS(P. - C
(s+ log
en
(20)
for positive constant C. Moreover, the upper bound in (20g also holds for the corresponding sorted estimators F, as defined in Corollary 1
provided that each function F is continuous and V S V
for all t.
This result shows that, under a sparse difference structure and a minimum segment length condition, trend filtering estimators for distributional regression can achieve nearly parametric accuracy (up to logarithmic factors). The bound in 20) mirrors the adaptive rates derived for standard trend filtering in point estimation problems; see, for example, Guntuboyina et al. (3]. Our extension demonstrates that the same refined rate behavior persists in the more general context of distributional regression, where the goal is to estimate full conditional distributions rather than scalar means.
The proof of Corollary 6 is provided in Appendix G.11
B General Result for Penalized Estimators
A natural alternative to shape-constraineed estimators is the use of penalized estimators, where the penalty term promotes a desired behavior the signal being estimated. Motivated by this approach, we present general result for distributional regression using penalized estimators in this subsection. Specifically, consider estimators of the form
F(t)
== arg min BER
(w() - 0112 + ).pen,
(21)
where de > 0 is a tuning parameter and pen R" - BR is a penalty function. We now present our main result for the penalized estimator defined in (2I).
Theorem 4. Suppose that pen, (.) is convex for all t and it is a semi-norm. In addition, assume that sup pen, (F)( V. Let K == {0 E R" pe),(0) I 6V}. Then for any nj 0, it holds that tER
IP
sup IER
(7()))(()= 2n²
S
C
n/4
v Vlog N(, (,, 1I) + Cvlogn n
1p
1
(22)
for some constant C > O, proivide that we set X = /2/4pen(F*(())
Theorem demonstrates that achieving uniform upper bound on the MSE can be accomplished by controlling the covering number of sets of the form K n (0), as outlined on the right-hand side of
(22).
We now
statistical guarantee for the penalized version of trend filtering in distributional regression. This result follows directly from Theorem and involves a calculation analogous to the proof of (14.
Corollary 7. Consider the estimator in 21 with pen, (0) = (r))
(8) and x chosen as in Theorem
4 Then
(F((0) = =0;
V
log n
sup EER
2r +
(23)
1=1
where = sup pen,(F (I). tER
23

---

C Dense ReLU Networks: assumption and definitions
In this appendix, we provide additional details for Section 3.3.2 Before outlining our assumptions on the functions (t), introduce notation related dense ReLU networks. To that end, describe a dense neural network with architecture (L, k) employing the ReLU activation function
given as p(s) max (0, 8} for any E IR. Such network is represented real-valued function f: I2d R satisfying the following properties:
k
((x)
Σ
9
(4)
(x)
(24)
for weights (L)
Ci, 9 E R and for f! 's recursively defined by
h
(z) =
Go
(6-1) (1--) (x) + 0
(0-1)
(25)
for some (a-1)
E R. (,,..., DJ. and (x) P (E x000 1.k.1
1--1)
(0)
1,0 .....
Cig
+ C,,0
with (0)
od E R.
Next we provide some notation necessary for define the class of signals where the (*)''s belong.
Definition 1 ((p, C')smmothneess) Let + s for some E N = Z+ U (00 and 0 < I 1. We say that a function g + R is (p. -msooothh for every a = (ai,..., ad) E N with E zt where Σj=1 Oj = the partial derivative 0rgg/ (Ou, du exists and
ggg
09
(u)
(v)
s Clllu e||
0u
Ou
Ou
Ou
for all E Rdd.
Let us now define the generalized hierarchical interaction models H(I,P)
Definition 2 (Space of Hierarchical Composition Models, (37). For | = and smoothness constraint
P S (0, oo) x N. the space of hierarchical composition models is defined
H(1, P) = (h: Rd + R h(a) = m (a(a())...... ((m(M))) where
m : RM + IR is (p, -mmsooothh for some (p, M) E P and म (1,..., M} + (1,...,d)}
For l 1. we
H1I, P) = {h: Rdd RR h(xx = m (f) far(a), where
m RMR.M. + R is (p, -mssmooth for some (p, M) P and fe E 1(1-
With the notation above, are ready to state assumption the true signals in the spirit of B7 Assumption 1. Suppose that for all the function G* is in the class H(l, P) in Definition In addition, assume that each function g² in the definition of G (., t) can have different smoothness for qg² E N, 8g* E (0, 1}, and of potentially different input dimension Mgs, so that 1pgt,Mgi E P. Let Mmax be the largest input dimension and Pmax the largest smoothness of any of the functions g² for all Suppose that for each g¹ all the partial derivatives of order less than equal to 9g* are uniformly bounded by constant Smoooh and each function g is Lipschitz continuous with Lipschitz constant CLip 2 Also, assume that max max(pmaxs M Mmax)
2
Pyt 4g* + 8gt, for
O(1)
D Rearrangement of estimates
Figure shows an example the original P, and rearrangement F, as Corollary i which is guaranteed to be non-decreasing.
24

---

1.00
1.00
0.75
0.75
0.50
0.50
0.25
0.25
0.00
0.00
Figure 2: The plot in the left shows display of an example of function P. and the right panel shows the corresponding rearrangement F; described in Corollary
E
Summary of Results and Theoretical Comparisons with Prior Work
For a summary of all theoretical results presented in the paper and comparisons with prior work such as [28. 44 42 62 23 37 10 58). see Table 2 This table synthesizes the convergence rates and underlying assumptions for our estimators (isotonic, trend filtering, and dense ReLU networks) alongside those of key existing methods discussed throughout the manuscript.
In addition, Appendix E. provides a complementary comparison that focuses specifically on methods designed for conditional distribution estimation. In particular, we contrast our Dense ReLU approach from Section [3..22] with the main competing estimators evaluated in the simulation section, Section 4 including those from 100 and [58] highlighting differences in assumptions, estimation strategies, and theoretical guarantees.
Finally, this appendix also includes dedicated subsection (Appendix E.2 that clarifies the concep- tual relation and technical distinctions between our Theorem [ and Theorem A. of Guntuboyina et al. 223] This comparison emphasizes how our framework extends the convex-projection setting to distributional regression with Bernoulli-type observations and CRPS loss, relying on different concentration and normalization arguments.
E.1 Error Bounds for Conditional Distribution Estimation
The result in Corollary establishes that the dense ReLU network-based estimator achieves convergence rate of order
log 72 n
o., log" n,
where o. max(p.Mee n (pp(2p+ for the worst-case squared error over thresholds t E IP, as shown in (17). This result provides high-probability bound on the squared & error between the estimated and true conditional CDF vectors
F(I) (f. ((..........))
and
(*) where (T) = (ye S1 wx) G( l), uniformly across E Here ((.,)))) C [0, 1] x R are i.i.d. data.
To contextualize this result, we compare it with the result obtained in [10 and 58
E.1.1 Comparison with Distributional Random Forests (DRF) 10
Work [10] studies conditional distribution estimation using Distributional Random Forests (DRF). Their framework assumes access to i.i.d. samples {(z,, }:) [0, 1,2 RRd, where each [0, 1|. denotes a covariate vector and yı Rd a multivariate response. The goal is to estimate the conditional distribution P(Y X = 2), where Y E IRd and X E [0, 1]", using a weighted empirical
measure of the form
iF x) = Σ wo( ,ng
25

---

where 6g. is the Dirac measure at ys. assigning unit mass to the observation y: E R Moreover, w, (2,) denotes the DRF-induced weight assigned training point z,. Specifically, after constructing forest of N randomized trees, each test point x E [0, defines, for every tree Tı. a leaf (*(x) consisting of training points 2i that fall in the same terminal node as x. The DRF weighting function is then defined as
Th (xy)
1 [mjEcr(m) N k=1
Σ
[xx(x)
This leads the following plug-in in estimator of the conditional CDF at threshold t R
PDRF
(t) = wx (xz) ((us)) (to X
which computes weighted empirical CDF based on the forest-derived weights. To align with our setup, consider the univariate case d = 1 with = dg. where each observation (I y1) E [0. 1]
R. In this setting, their goal becomes estimating the conditional CDF Fy
P(Y S t
= x) for any threshold t E R. The DRF estimator in this case simplifies to PDDF y|X== (t)
Σ==1 wx (1²) 1 {リ, S tt, providing weighted empirical estimate of the conditional distribution function at point Viewed from our perspective, where the focus is the sequence model F (1) (F)(),.... (() with (L) == IP (y: / x4), the DRF method yields estimates
ADRE (t) Σ j=1
Wx. (xz) 1 {uy S.
where the weight vector Wz,() (-) is centered around the design point 2;. These estimates EDRF (t) are directly comparable to our estimator (t). In terms of assumptions, their theoretical guarantce relies on several conditions about how the weights wz, (1y) are constructed. One key requirement is that each tree in the forest is built using a random subsample of the data, rather than sampling with replacement. The size of this subsample 8n is allowed to grow with n, specifically with s, n" for some 3 < They also assume that the covariate distribution on [0, 1] has density bounded away from zero and infinity. This contrasts with our setting, where no such condition on the distribution of the covariates is required. Another assumption is that for each threshold the conditional CDF function x + Fy)x==(() is Lipschitz. Translated to our notation, this corresponds to requiring that for every E R. the function t), where F (t) G( t). is Lipschitz in In terms of our assumption that G(.,t) E H(l, P). stated Assumption this equivalent to assuming that G.,,t) H(1, {11, do))), that is, G(,,t) t) belongs to Lipschitz class, which is strict
[
subset of the more general hierarchical composition model we allow.
Under these conditions, the theoretical guarantees in 10] for the DRF method are limited to con- vergence in probability. Specifically, Corollary in their paper establishes that the DRF estimator j.DRE (t) converges in probability to the true conditional CDF P7 (t) at fixed thresholds & I. However, this result does not provide explicit rates of convergence in mean squared error or uniform guarantees over
In contrast, our result provides a bigh-probability bound on the & error over the full vector F*(t). with convergence rate
(1(1)) (R(1--700 (P(0))7(0))= =Op loog + on log n n),
sup ER
where = matches the approximation rate for the Lipschitz class assumed in their setup. More generally, our framework also provides rates under weaker smoothness conditions, allowing each threshold function G(.,t) to belong the hierarchical composition model class (11 P), as defined in Assumption]
summary of this comparison is provided in Table 2
E.1.2 Comparison with Engression (EnG) 58
The work on (58 introduces method for estimating conditional distributions using energy score minimization. Their framework assumes access to i.i.d. data W.)R C. [0, 1]d x E. where
26

---

each E [0, 1] is covariate vector and E R is real-valued response. The main objective is to estimate the full conditional distribution P(Y 2).
To solve this problem, Engression models the conditional distribution Y X = x: as a transformation of independent noise variable Unif(0, 1]. i.e., it posits the model
Y = EX,,)
for some unknown function g²: [0, 1J [0, R.
The goal to learn an approximation ge g", parameterized by a deep neural network, such that conditional distribution of ex(X X matches that of Y X. The model is trained by minimizing an empirical energy score loss, which quantifies the discrepancy between the observed responses and the model-generated samples. For each data point (æ,, 34), the method samples independent noise
variables Ei,1...,, Ei,m ~ Unif(0, 1]. The resulting empirical loss function is
iid.
min
y g0 Esg)l 2m((m 1)
1
Σ j71
(00( go (1 E
m
Once the network go is trained, the conditional distribution at test point z can be approximated by generating . independent samples E (1....,M) M).
iid. Unif[0, and evaluating go Em) for each
To align with our setting, we interpret their covariate dimension d as our notation dg. so that each design point E [0, 1]de In our sequence-based setting, this leads to the following estimator of the conditional cumulative distribution function:
F.nnG
(t)
1 M
Σ 1gg0(ss, S st.. m=1
which provides approximation to the true conditional CDF (t) P(y S t x;). This estimator is directly comparable to our method's output (t). and both aim recover the conditional distribution function evaluated at x; for any threshold & E R.
Specifically, under the assumptions that the true function g" is Lipschitz continuous in both arguments and that the marginal distribution of covariates has a density bounded away from zero and infinity on its support, [58] shows that the estimator BBB (t) converges in probability to F (t) for any fixed threshold IR.
Translating these conditions to our notation, where (t) = G*,, (1r,,) t). the assumption that g" is Lipschitz in its first argument implies that the map + G (x, t) Lipschitz for every fixed In Assumption posits that G(,,t) E H(l,P), rich class of
contrast, our structural assumption hierarchical compositions of functions. The Lipschitz class assumed in Engression corresponds to the special case H(1, {[1, do))). Therefore, our framework subsumes the setting considered in Engression and allows for more expressive function classes beyond global Lipschitz continuity. Additionally, we do not require the covariate distribution to be bounded away from zero. Importantly, estimator comes with high-probability guarantees the & error over the full CDF vector Ft(t) =
Specifically, we establish that
(A((t) - (((00)=(())) = =Op (logn + On log n m))
sup AER
where the term &n depends on the complexity of the class H(1, P). In particular, when G(: () is Lipschitz, we obtain 6 = n 2/(2
A summary of this comparison is provided in Table 2
E.2 Relation between Theorem 1 and Theorem A.1 in Guntuboyina et al. 23
Our Theorem I conceptually related to, but technically distinct from, the high-probability bound established in Theorem A.l of Guntuboyina et al. 23 (restated Theorem Their setting considers
27

---

the estimation of a mean vector E I" based on noisy Gaussian observations Y = + Z. where Z ~ (0, I.). The estimator f is the Euclidean projection of Y onto a convex set K R", and the goal is to bound the squared error 7 = in expectation, which is then strengthened to
high-pne brobbbbbbbbiity bound.
To achieve this, they define the random process
8²2
H(t) =
sup PEK-U-uHES
Z (v- - 11)
2
which is strictly concave in t and attains its unique maximum at t" = -177 -11l2. The corresponding deterministic function
fu(() = [t)
is also concave and attains its maximum at tu.. which characterizes the typical estimation scale. The key analytical step in Guntuboyina et al. [23] is to show that t* concentrates sharply around ty. via a peeling argument: within a local window [t, - Ivtt tu + x /ttll. stochastic fluctuations of (t) - (t) are controlled by Gaussian concentration, and outside this window, the deterministic function (t) exhibits quadratic decay. This interplay between curvature and Gaussian concentration yields strong probabilistic control of the projection error.
Our framework preserves the geometric structure of convex projection but differs several essential respects. First, the observations in our setting are
w(t) = (1 S =1))
which are Bernoulli random variables that play the role of Y. The associated noise vector e(l)
w(() - F (L) has independent, mean-zero, heteroskedastic components with variances P7 (4)(1 (Q). This noise is sub-Gaussian but not Gaussian, and thus lacks the rotational invariance and linearity properties exploited in Guntuboyina et al. [23]
Second, while their goal is control the squared Euclidean error - - our analysis focuses on bounding the expected CRPS loss between the true conditional distribution functions F and the estimators F,, which serve as the natural analogues of p. and ş in our context. The non-Gaussiynity of the noise, the use of the CRPS loss, and the absence of rotational symmetry make the direct application of their argument infeasible.
In Theorem 1 our proof reduces the problem to controlling the same central quantity, namely the squared error ||F Fl,,,, which corresponds to ||ž - μ|l2 in their notation. The key scale parameter in our analysis, denoted balances stochastic variability and approximation bias, serving a role analogous to ty in their framework. We then show that || F* concentrates sharply around n using a peeling-sty seylln argument similar in spirit to theirs.
However, because of the heteroskedastic and Gau-Gansii nature of our noise, the technical realiza- tion of this step is entirely different. Our analysis relies on
1. a lab-Gaussian concentration inequality for Lipschitz and separately convex functions, used to control the local empirical process supsce (1)))) E(t) (0 (t)); and
2. a star-shaped self- noffaiormaliztion argument that extends this local control uniformly over the full constraint set K,.
In summary, both approaches employ pecling argument to control the projection error, but they differ fundamentally in probabilistic structure and analytical tools. Theorem A.I of Guntuboyina et al.
exploits Gaussian concavity and symmetry, whereas our proof adapts non-Euclidean empirical process methods to accommodate aub-Gaussian, heteroskedasttic noise within the CRPS- based
2
framework of distributional regression. The novelty of our result lies in extending the probabiiiyy convex-projection analysis to this more general and statistically heterogeneous setting.
Therefore, Theorem I can be viewed as CRPS- analogue of Theorem A. 1 in Guntuboyina et al. (23. adapted to the non-Gaussian, heteroscedastic, and distribbutionallregrsssoo setting considered in this paper.
28

---

Table 2: Comparison of convergence rates for isotonic, trend filtering, and dense ReLU network methods in distributional regression. The "Old" rows refer to existing methods in the literature, while "Our" rows describe the corresponding results introduced in this work.
ISOTONIC
METHOD
ASSUMPTIONS
CONVERGENCE RATE
K, MONOTONE CONE (t( E Kg
28
o(1)
PRIOR WORK
K, MONOTONE CONE
F (L) E Kr. DENSE DESIGN
Tt 222/(20t 1) WORSE-CASE MSE)
HOLDER CONTINUITY
GENERAL RATE
K₁ MONOTONE CONE F (L) E Kr
n n-2/3 (CRPS, WORSE-CASE MSE)
THIS
WORK
FEW STRICT INCREASES WITH THE NUMBER OF STRICT INCREASES k(L) &
FAST RATE
1+MO log T+R(0)
en
(CRPS)
TREND FILTERING
METHOD
ASSUMPTIONS
CONVERGENCE RATE
F = - 10ER" TV) W f WITH f THE REGRESSION
(2/(2++)) 13//(3r11) (MSE)
442
62, 23
PRIOR WORK
FUNCTION Ys = f (x.) +
SPARSE D
2
+11 log
(MSE)
MIN SEGMENT LENGTH
K₁ k== f0e (O : (V) (T()) E Kr
22/(2r 11)
THIS WORK
GENERAL RATE
20/(0+11
(CRPS, WORSE-CASE MSE)
SPARSE D (*()) MIN SEGMENT LENGTH
FAST RATE
atd log ++1
(CRPS)
DENSE RELU NETWORKS
METHOD
ASSUMPTIONS
CONVERGENCE RATE
REGRESSION yi = f(xi) +
du (MSE)
H
(t) = Go (x,, = P (y/ S f x4) (' E (1,
PRIOR WORK
100
o(1)
FOREST-WEIGHTED EMPIRICAL CDF
(t) G (x,,t) (M, St x))
58
(' H1,, do))
o(1)
NEURAL NET-BASED SAMPLING
THIS WORK
(t) G (x,, (v. st)x)) (.tt)) E H(IP))
logn on log" n
GENERAL RATE
(CRPS, WORSE-CALE MSE)
Additional Numerical Results
This appendix provides an extensive evaluation the proposed method, UnifDR, including additional results and analyses omitted from Sections 4 and These supplementary results further demonstrate the effectiveness and robustness of our methods across diverse settings. The appendix presents:
Additional real data applications to illustrate the practical utility of the proposed methods (Appendix F.T.
Comprehensive evaluations alternative evaluation sets (A₁ and Ag) supplementing the results for A2 Section 4 This section also includes missing evaluation results for A2 related to the CRPS metric (Appendix F2
29

---

Performance results based on the Maximum Squared Difference (MSD) metric across all scenarios, which were included in Section (Appendix] F3)
F.1 Additional Real Data Applications
This appendix presents two additional real-world data applications to further demonstrate the effective- ness of the proposed methods. These examples span different domains, illustrating the versatility and robustness of our approach. Each case study includes a description of the dataset, the experimental setup, and a comparative performance analysis.
For each dataset, the data is randomly divided into training subset (75%) and a testing (25%) subset. Model performance is evaluated using the empirical cumulative distribution function (CIF),, (t)(), computed over 100 evenly spaced points a predefined set A. The performance of each competitor is assessed using the Continuous Ranked Probability Score (CRPS) and the Maximum Squared Difference (MSD) metrics, defined as follows:
CRPS
Σ
(P((1)
(0000-w)())))) mi(1)
Test
100
YETest
ΛEA
and
MSD
max
Σ TeTe
[i(t) w;(()
TEA
Test
The proposed UnifDR method is implemented in two variants: Trend Filtering, which captures smooth variations through total variation regularization, and Dense ReLU Networks, which leverages a deep neural network to incorporate covariate information.
F111 Chicago Crime Data with Dense ReLU Networks Approach
We revisit the 2015 Chi Caic crime dataset, previously analyzed in the main text using the Trend Filtering approach. We remember that the dataset contains reported crimes in Chicago throughout 2015. As before, the spatial domain is discretized into 100 100 grid, where each grid cell represents an aggregated crime count. The response variable remains the log-transformed total crime counts per grid cell, and grid cells with zero observed crimes are excluded, yielding a final dataset of 3,844 grid cells.
Unlike the Trend Filtering approach which assumes a smooth index trend, the Dense ReI.U network approach leverages covariate information for modeling crime intensity. The following covariates are included. Latitude and Longitude Bins, encapsulating spatial crime patterns. Day of the Week, represented using dummy variables for each weekday (Monday through Sunday). Beat, a categorical identifier for Chicago's ' policing districts. Arrest Indicator, a binary variable denoting whether an arrest was made (1) or not (0).
The dataset is randomly split into 100 train (75%) test (25%) partitions, with evaluation conducted evenly spaced points A ranging from to 6. The Dense ReLUU Networks approach employs fully connected feedforward architecture with five hidden layers of 64 neurons each, using ReL.U activations. The model is trained using the Adam optimizer with learning rate of 0.001 over 1,000 epochs, minimizing the Binary Cross-Entropy (BCE) loss function for improved CDF estimation.
Performance is assessed using the Continuous Ranked Probability Score (CRPS) and Maximum Squared Difference (MSD) metrics, comparing estimated CDFs F((() against empirical indicators w))()), where /. E A. The methods CART, MARS, RF, DRF, and EnG serve as competitors for the Dense ReLU Networks approach. Table summarizes the results demonstrating the superior performance of the Dense ReLU network relative to classical nonparametric regression methods. Additionally, Figure 3 visualizes Ft(t) for t = 3 across test grid cells for all competitors.
F.1.2 California Housing Prices
We evaluate the effectiveness of the proposed UnifDR method by analyzing the 1990 Califor- nia housing dataset, which contains demographic and economic information from various neigh- borhoods across California. Originally introduced in Pace and Barry 45), the dataset com- prises 20,640 observations and is publicly available via the Carnegie Mellon StatLib reposi-
30

---

Predictions of the CDF of Chicago Crime Data (Dense ReLU Networks) UnifDR
CART
MARS
Latitude
60
Random Forest
DRF
Engression
80
Latitude
Re
20
Longitude
Longitude
Longitude
1.0
Prediction
Figure [1) for = 3 and all Test, for all competitors.
Predictions of the CDF of Chicago Crime Data (Trend Filtering) UnifDR
Add55
100
Latitude
60
40
20
01
20
40
60
B0
100
0
20
40
60
80
100
Longitude
Longitude
0.891
0.892
0.893 Prediction
0.894
Figure 4: F.(() for / = and all E Test, for all competitors for the example in Section
tory at httt:://iibb.ttatt cmu wwt...m.atttttttttttsss// as well as the lattps://ww..cc..ffff
Itorgo/Rogressionn cal www.hhh..inn.................... html portal.
Following the approach of Ye and Padilla [66], the geographic area is discretized into 200 200 spatial grid based on latitude and longitude coordinates. The response variable derived by applying logarithmic transformation the median house values within each grid cell to enhance numerical stability and interpretability. Grid cells lacking valid data are excluded from further analysis, resulting final sample size of 3,165 grid cells.
In the Trend Filtering approach no covariates are included. Instead, the spatial grid is treated as an or- dered sequence, allowing total variation regularization to capture smooth spatial variations in housing prices. In contrast, the Dense ReLU Networks method integrates spatial features such as latitude and longitude with socioeconomic attributes like median_ income and average_occuppapeeeeeeeeeeeeeeee occupancy, which is computed as population divided by households. This approach enables the model to capture complex relationships influencing housing prices. The Dense ReLU neural network consists of three hidden layers, each containing 30 neurons followed by a ReL. activation function.
Evaluations performed over 100 evenly spaced points within the range A 15, 15).
31

---

Table 3: Evaluation metrics for UnitDR (Dense ReLU networks) and its competitors on the Chicago crime dataset.
METHOD
CRPS (MEAN + STD)
MSD (MEAN + STD)
UNIFDR
0.08111 0.0018
0.2133 = 0.0031
CART
0.0951 0.0017
0.2622 = 0.0071
DRF
0.0906 0.0018
0.2477 0.0042
ENG
0.1014 0.0028
0.2652 = 0.0045
MARS
0.0974 0.0019
0.2732 + 0.0047
RF
0.0934 + 0.0020
0.2581 + 0.0031
Table 4: Evaluation metrics for UnifDR (Trend Filtering) and its competitor AddSS the California housing dataset.
METHOD
CRPS (MEAN + STD)
MSD (MEAN + STD)
UNIFDR
0.0343++ 0.0013
0.2505 = 0.0097
ADDSS
0.0357 + 0.0028
0.2652 = 0.0407
Figure presents the estimated cumulative distribution functions (t))) at t = 12, comparing the performance of UnifDR (Trend Filtering) against its competitor, AddSS. Similarly, Figure fillustrasss FigureBillustrates the estimated cumulative distribution functions at the same evaluation point, this time comparing UniTDR (Dense ReLU networks) against CART, MARS, RF, DRE, and EnG. The corresponding evaluation metrics, summarized in Tables and demonstrate that UnifDR consistently outperforms all competitors in terms of both CRPS and MSD.
Predictions of the CDF of California Housing Prices (Trend Filtering) UnifDR
AddSS
200
150
Latitude
100
0-
50
100
150
200
0
50
100
150
200
Longitude
Longitude
0.958
0.960
0,962
0.964
0.966
0.968
Prediction
Figure 5: Estimated distribution function (t) at = 12 for all grid cells in the test set, comparing UnifDR (Trend Filtering) and the AddSS competitor.
F.1.3 Ozone Data Analysis
We further evaluate the effectiveness of the proposed UnifDR method by analyzing ozone concentra- tion data collected from the Environmental Protection Agency (EPA) Regional dataset. This dataset consists of daily ozone measurements collected across various monitoring stations in diffeerntt regions of the United States for the year 2024. The dataset includes measurements of ozone concentration along with associated variables such as Air Quality Index (AQI), wind speed, temperature, latitude, and longitude. The available variables include: State Code, County Code. Site Number, Latitude, Longitude, Date (Year, Month, Day), Ozone concentration (in parts per million). AQI, Wind Speed (in miles per hour), Temperature (in Fahrenheit), Observation Percentage (percentage of valid ob- servations for given day). First Maximum Value (highest ozone level recorded in a day). First Maximum Hour (time at which the maximum value was recorded), and Observation Count (number of valid ozone measurements per day). The data publicly accessible via the EPA's AirData portal
32

---

Predictions of the CDF of California Housing Prices (Dense ReLU Networks) UnifDR
CART
MARS
150
Latitude
50-
Random Forest
DRF
Engression
290
150
Latiudd
50
150
200
150
200
200
Longitude
Longitude
Longitude
0.0
0.2
0.4
0.6
0.8
1.0
Prediction
Figure 6: Ft(t) for t = 12 and all i Test, for all competitors.
Table Evaluation metrics for UnifDR (Dense ReLU networks) and its competitors on the California housing price dataset.
METHOD
CRPS (MEAN + STD)
MSD (MEAN 1 STD)
UNIEDR
0.02091 0.00055
0.1450 0.0047
CART
0.0258 0.00053
0.1984 0.0067
DRF
0.0229 = 0.00051
0.1755 = 0.0045
ENG
0.0244 = 0.00059
0.1924 = 0.0057
MARS
0.0233 + 0.00057
0.1798 + 0.0099
RF
0.0221 0.00046
0.1650 0.0047
at httpp:://agg epa. spanning multiple years can be downloaded.
html where historical records
Following the approach of previous studies, the geographic region of interest is discretized based on latitude and longitude coordinates, specifically within the range 30 N 50 and -153 W to W70% In the Trend Filtering setup, each monitoring site within these bounds is uniquely identified using a lexicographic ordering its latitude and longitude values. The response variable y; for each site i is defined as the mean of the daily recorded ozone levels, ensuring robustness against short-term fluctuations and missing data. This approach results total of 1,189 unique monitoring sites, which serve as the basis for spatial trend estimation.
In contrast, the Dense ReLU network method method allows for flexible spatial fit by incorporating spatial location data as covariates, capturing complex relationships that influence ozone concentration. Specifically, the model uses latitude, longitude, mean AQI, mean percentage of valid observations, mean 1st maximum ozone value, mean 1st maximum hour, and mean observation count as input features. The Dense ReLU network consists of two hidden layers, each containing 100 neurons followed by ReLU activation function.
Evaluations are performed over 100 evenly spaced points within the range A [0, 1]. Figures and g
present the estimated cumulative distribution functions (t( at 0.03, comparing the performance of UnifDR using Trend Filtering and Dense ReLU network against their competitors. Tables [ and 7 summarize the evaluation metrics, demonstrating that UnifDR achieves superior performance in terms of both CRPS and MSD.
33

---

Predictions for Ozone Concentration CDF in the USA (Trend Filtering) UnifDR
AddSs
0.18
0.19
0.20
0.22
0.23
0.24
0.25
Ozone Concentration
Figure 7: Estimated distribution function (t()) for all monitoring sites in the test set, comparing UnifDR (Trend Filtering) and the AddSS competitor.
Table 6: Evaluation metrics for UnifDR (Trend Filtering) and its competitor AddSS on the ozone concentration dataset.
METHOD
CRPS (MEAN + STD)
MSD (MEAN + STD)
UNITDR
0.0027.1 0.0002 0.0035 + 0.0002
0.1581 + 0.0124 0.1987 0.0127
ADDSS
Predictions for Ozone Concentration CDF in the USA (Dense ReLU Networks) UnifDR
CART
MARS
Random Forest
DRF
Engression
0.0
0.2
0.4
0.6
0.8
1.0
Ozone Concentration
Figure 8: Ft(t) for 0.03 and all i Test, for all competitors.
Table Evaluation metrics for UnifDR (Dense ReLU networks) and competitors on the Ozone dataset.
METHOD
CRPS (MEAN + STD)
MSD (MEAN STD)
UNIEDR
0.0016 0.00021
0.0959 0.0136
CART
0.0026 1 0.00028
0.1529 + 0.0209
DRF
0.0031 = 0.00016
0.1593 = 0.0118
ENG
0.0022 0.00099
0.1097 0.0603
MARS
0.0022 0.00039
0.1288 0.0376
RF
0.0021 = O. 0.00023
0.1012 = 0.0157
34

---

CRPS results for Scenarios S1 and S2 in A
CRPS results for Scenarios S1 and 52 in Ag
0.0200
Scenarios
0.015
51
51
0.0175
CRPS
0.010
CRPS
0.0150
0.0125
0.005
0.0100
0.0075
0.0050
n-400
n-800
n-
n-1600
n-400
800 Sample Size
n-1600
Sample Size
CRPS results for Scenarios S1 and S2 in A3
Scenarios
0.006
52
+
CRPS
0.004
0.002
4=000
n=800
n= n=1600
Sample Size
Figure 9: Box plots for CRPS results in Scenarios S1 and S2. and A2 (right), while the bottom row displays results for Ag-
The top row shows results for A₁ (left)
F.2 CRPS results on evaluation sets A,, and Aa, and missing results for S1, S1 and S6 in Ag
Section presented results for the evaluation A2 using the Continuous Ranked Probability Score (CRPS). However the results for scenarios S1, S2 and S6 were omitted due to space constraints. Moreover, it is worth noting that the evaluation set A2 represents a balanced range of values centered around zero. To provide a more comprehensive analysis, this appendix includes the omitted CRPS results for Ag, and the CRPS results for the alternative evaluation sets A₃ and Ag, which emphasize distinct distributional regions.
The new evaluation sets are defined follows:
A1: : 100 points evenly spaced between -1 and 0.4, focusing on the lower and middle ranges of the distribution.
Ay: 100 points evenly spaced between 0.8 and 10, capturing the upper tail of the distribution.
This extended analysis provides deeper insights into the robustness of the proposed methods in varying distributional regimes, including regions with lower densities and heavier tails. As described in Section 4 for each evaluation (A1, Ag and Ag). the CRPS is computed and averaged over 100 Monte Carlo repetitions. Figures 9 through t and evaluation regions. Specifically, Figure correspond to Scenarios S1 and A₁, A2, and Ag- Figure T0 presents Scenarios S3 and $4 in A and Ag. Figure T1 focuses on Scenario S5 in A₁ and Ag. Figures
14 summarize the CRPS results across different scenarios
12
13
and 14 consider Scenario S6 in A1, Ag, and Ag. respectively.
Across all scenarios, UnitDR is implemented using different estimation methods: the isotonic estimator (Section Ib b for S1 and S2, the trend filtering estimator (Section 3.2.2) for S3 and S4, and the Dense ReLU Networks method (Section 3.3.2p for S5 and S6. Our proposed framework UnifDR outperforms competing approaches across all scenarios, evaluation regions, and sample sizes, regardless of the specific estimation method employed. These findings highlight the versatility and robustness of UnifDR adapting to diverse structural patterns within the data. Overall, performance trends remain consistent across A1, A2, and Ag. with minor variations due to differences in the underlying data distributions. This extended analysis further reinforces the effectiveness of the proposed methods under varying evaluation conditions.
35

---

CRPS results Scenarios and S4
CRPS results for Scenarios S3 and Ag
Competitors
Competttors
0.06
53/UnIDRR
53-JnfbRR
53-Ade55
AddS
54-UnFDR
54-Adu55
CRPS
0.04
54 AddSS
CRPS
0.04
0.02
0.02
0.00
0.00
n=000
800
n=1600
n=400
n=800
n=1600
Sample Size
Sample Size
Figure 10: Box plots for simulation results of S3-S4 for the CRPS metric. The row shows results for A1 (left) and A2 (right).
CRPS results for Scenario 55 in A
CRPS results for Scenario 55 in Aa
Computioon
Competitors
0.04
0.04
DeF
0.03
BrG MARS
0.03
EnG MARS
CRPS
CRPS
0.02
0.02
0.01
0.01
0.00
0.00
n-400
n- 800
n-1600
40-00
n- 800
n-1600
Sample Size
Sample Size
Figure 11: Box plots for simulation results of S5 for the CRPS metric. The row shows results for A (left) and A2 (right).
CRPS results Scenario A
CRPS results Scenario in Ag
0.06
Compet tors
CART UniT
0.0150
Compotitr's
CART UnimA Der
0.0125
End
0.04
CRPS
0.0075
0.02
0.0050
0.0025
0.00
0.0000
n=400
n=800
n=1600
n=400
n=800
1600
Sample
Sample Size
CRPS results for Scenario S6 Aq
CRPS results for Scenario $6 in A
0.004
0.003
0.0008
CRPS
0.002
CRPS
0.0006
0.0004
0.001
0.0002
0.000
0.0000
n=1600
n=400
n=800 Sample Size
n=1600
Sample Size
Figure 12: Box plots for simulation results of S6 for the CRPS metric for the set A1. The top row shows results for the all the competitors (left) and, competitors with median below 0.01 (right). The bottom displays results for competitors with median below 0.0025 (left), and best two competitors (right).
36

---

CRPS results for Scenario S6 in Az
CRPS results for Scenario S6 in Ag
Compatiaan
0.025
Competroar
0.06
Uniinn
0.020
CRPS
0.04
CRPS
0.015
0.010
0.02
0.005
0.00
0.000
n-400
n-
n-1600
n-400
n-800
n- 1600
Sample Size
Sample Size
CRPS results for Scenario S6 in A/
CRPS results for Scenario 56 in A2
0.0125
0.0025
0.0100
0.0020
MARS
0.0077
0.0015
0.0050
0.0010
0.0025
0.0005
0.0000
400
n=800
n=1600
n=400
n=800
n=1600
Sample Size
Sample Size
Figure 13: Box plots for CRPS results for S6 Ag. The top row shows results for the all the competitors (left) and, competitors with median below 0.02 (right). The bottom displays results for competitors with median below 0.01 (left), and best two competitors (right).
CRPS results for Scenario S6 in A3
CRPS results for Scenario S6 A3
Comoetitors
0.005
Compettrrs
0.025
CART
UnioOl
UnFOR
DR
DRF
0.020
EnG
CRPS
0.015
CRPS
0.003
0.010
0.002
0.005
0.001
0.000
0.000
n=400
n=800
n=1600
n=000
n=800
n=1600
Sample Size
Sample Size
Figure 14: Box plots for simulation results of S6 for the CRPS metric using evaluation Ag. The left plot corresponds to all competitors performance, while the right plot corresponds to best two competitors.
F.3 Additional Results for Maximum Squared Difference (MSD) Metric
This appendix extends the results presented in Section and Appendix F.2 by providing evaluations of the Maximum Squared Difference (MSD) metric across all scenarios for the evaluation sets A A2, and A3. The MSD) metric measures the worst-case discrepancy between the estimated and true cumulative distribution functions (CDFs), offering a stringent assessment of model accuracy and robustness.
Figure 15 displays box plots of the MSD results for Scenarios and S2 across evaluation sets. As outlined in Section 4 UnifDR employs isotonic regression for these scenarios, which lacks direct competitors. The results indicate a decreasing trend MSD values as the sample size increases, demonstrating the consistency and improved accuracy of isotonic regression with larger datasets. Moreover, variations in MSD across A,, A2, and Ag reflect natural differences in the underlying data distributions.
Figure 16 presents MSD results Scenarios S3 and $4, where UnifDR employs trend filtering and compared against the additive smoothing splines (AddSS) method. Across all evaluation sets and
37

---

MSD results for Scenarios S1 and 52 As
MSD results for Scenarios S1 and S2 in Az
0.10
Scenarios
51
51
0.08
MSD
0.020
0.06
0.04
0.015
0.02
n-400
n-
n-1600
n-400
n--800
n-1600
Sample Size
Sample Size
MSD results for Scenarios S1 and $2 in Aa
0.10
Scenarfas
52
0.08
MS0
2.066
0.04
0.02
n=4
n= B00
1600
Sample Size
Figure 15: Box plots for MSD results in Scenarios S1 and S2. The top row shows results for A₁ (left) and A2 (right), while the bottom row displays results for Ag-
sample sizes, UnitDR consistently outperforms AddSS, demonstrating its ability to adapt to complex structural variations. Notably, trend filtering exhibits particularly strong performance in regions with sparse or heavy-tailed data distributions.
Figures 17 to 20 summarize MSD results for Scenarios S5 and S6, where UnifDR is implemented Dense ReLU Networks (Section 3.3.2). The results further confirm the superiority of UnifDR over all competing methods, maintaining its advantage across different sample sizes and evaluation sets. As observed with the CRPS metric, UnifDR effectively captures diverse structural patterns and remains robust in challenging scenarios involving data sparsity and heavy tails.
The inclusion of MSD results provides comprehensive performance evaluation of UnifDR, reinfore- ing its effectiveness and reliability across various experimental conditions.
38

---

results for Scenarios
results Scenarios
S4
0.5
Compet tors
Camoetitor
53 UniDR
0.125
53-UnFDR
0.4
53-1dd55
53-Ado55
0.100
M50
0.3
Add5S
0.075
54-Ade5s
0.2
0.050
0.1
0.025
n=400
n=800
1-1600
800
n-1600
Sample Size
Sample
MSD results Scenarios and in Ag
0.150
Competitors
0.125
$3-1n RDR
53-to055
0.100
$4 Aco55
W50
0.075
0.025
Sample
Figure 16: Box plots for MSD in S3 and S4. The top row shows results for A₁ (left) and A2 (right), while the bottom row displays results for Ag.
MSD results for Scenario $5 in Ag
MSD results for Scenario $5 in Av
Competitors
0.08
CART
LNDDR
CART
DEF
0.125
DAF
0.06
0.100
0.04
0.075
RF
0.02
0.025
0.00
0.000
n=400
n=800
1=100
n=400
n=800
n=1600
Sample
Sample
MSD results for Scenario A3
0.150
Competitors
CART
0.125
UnitDA
Fri
M5D
MARS
0.007
0.050
0.025
0.000
n=400
n=800
1=00
Sample Size
Figure 17: Box plots for MSD in S5. The top row shows results for A₁ (left) and A2 (right), while the bottom row displays results for Ag.
39

---

MSD results for Scenario S6 in A
MSD results for Scenario S6 in A1
0.150
Competitors
0.125
0.006
UnifDR
0.100
En6
0.004
MS
0.075
MSD
XI
0.050
0.025
0.000
0.000
n=400
n=800
n=1600
n=400
n=800
n=1600
Sample
Sample Size
Figure 18: Box plots for MSD in S6 using evaluation set A₁ The left plot corresponds to all competitors performance, while the right plot corresponds best two competitors.
MSD results for Scenario S6 in
MSD results for Scenario S6 in
0.150
Competitors
CART
0.020
Compettors
UnfOR
0.125
UnitDR
DRF
0.015
0.100
EnG
MARS
0.075
0.010
0.050
0.005
0.000
0.000
n=400
n=800
n= 11600
n=400
n=800
n= 1600
Sample Size
Sample Size
Figure 19: Box plots for MSD in S6 using evaluation set A2. The left plot corresponds all competitors performance, while the right plot corresponds to best two competitors.
MSD results for Scenario 56 Ag
MSD results Scenario in
Competiors UnroR
0.15
CART
UnifDR
0.025
DRF
0.10
0.020
MSD
MARS
0.015
0.05
0.010
0.005
0.00
0.000
4=000
n= 800
n=1600
n=400
n=800
n=1600
Sample Size
Sample Size
Figure 20: Box plots for MSD in S6 using evaluation set Ag. The left plot corresponds all competitors performance, while the right plot corresponds to best two competitors.
40

---

G Proofs
G.1 Proof of Lemma
Proof. Let V(1) S 912) ... S V()) the order statistics of y. Notice that
L(F)
Σ CRPS(F,, 1146 1)
(B((1) 11 S 14)
"
Σ Σ (B() 18ys S 13
where A = (-00, M(1)): A₂ = (()) 1(2)):- ..... An ()) m(n)); An-1 [V(() 00). However, for every , 1] and t, AA we have that
min F(()EK
(P,(( 114, s 0))2²
min F(D)E
1 1/y st))²
Hence, letting t, be an element of Ay. we obtain that minimizing L(F) with the constraints F(f) K for all is equivalent to solving the independent problems
Σ (F((t)) - 1iv S t2 1=1
min FOLJER
and the claim follows.
G.2 Proof of Theorem
Theorem 5. [Theorem A.I in [23]] There exists a universal positive constant C such that for
every t
Cmax(17; max Ff(t) - (+))
(((t)-F(())
iI....,
E
<
for every n satisfying
E
sup AEK, : (fllsn
(0 - F*))
$
(26)
where e(t) w(t) (t).
In the followig we first present the proof of Theorem S
Proof. Let & E R. Define x² = max (f) - Fý (t)). We consider the following two cases separately based on the value of a:
s=1...,
1. σ
2. 0 +
Case =0
By definition, o² = max:=1.....
- - 0. This implies that for all = 1, 1,..., n,
F (t) - Fi (t)) = O. Since F (t) Ft) = 0. it follows that either F (t) 0 or (t)
for each i. Now observe that for each i E [l,...,,1} either (t) 0 or (t) = 1. Given that
((w;(t)) F." (t), it follows that (t) (t). Therefore, by the definition of F in Equation 3
we have F F*. this case, holds that
(()) - F (1) O,
41

---

and the result is obtained trivially.
Case 2: + O.
Denote by (c) == (0- (t) E p : 110 F (1)11 S 7), where pp.(1) = {0 - F (L)
0 E Kr}. Notice that the function L : R given by
|= (0 F (t))
sup
(
is 1--iischitz. Moreover observe that L is separately convex. In fact, for any k E {1,...}} we have that
(f1, Ek t + (1 t)(k .et 6n)
(2)
= (le1,..., tell........111 16 (1) .1Ex+11..,,en)
+(11 - ()E1, (1 - --))E-11 1) (2) ,1 - ()ER+t) - 1)mm)
-leR,1 + -1)Ex,2
Therefore,
L (k., tc (1) + (1--)cp (2) (.11...,, Gn)
[tec,1 + (1 - -)) T (0- - (t)
sup ()
[((k,]] (0--F(())) (1)) (1 +)((kx)) ((0k,)))(0001)))))) (0-FF(())) F*()))
sup #e0r-()())
sup ()
(*) + (1 8) ( - F
where the inequality is followed by triangle incquality. Using the fact that the supremum of the sum is bounded by the sum of the supremums, it follows that
L ((e......E-II (2) + +11-))E (2) (..11,....,) 6n)
<t sup (((.1)(( (0- (0-F + (1 -t) sup (E0r(((1)
((Ex,2) ((-F:2)
=th [(1,..,EHI1, ( , en) + t) (E1...,,E--I, ext1. (.411.)),044111)))))
Thus, is separately convex. Notice that (1) by Theorem 3.4 in Wainwright (64, for any 8 >>
(1 s - (1)) satisfies |c(t) S Hence,
sup (
F (t))
SE E
sup (
le(t) (-0)) onp
(27)
with probability at least e -6602/16 Next, we have that
[e(t) (0 (0-F S max 110 - (t) ,I)((n)
(28)
17
for any 0 E Ke, where Σ(n) supoeo (e)))) |e(t) (0 F* (t))]- This conclusion is derived based on the subsequent line of reasoning. If Kr and [|0=F(()]]l S n,. then F- (1) E O (n) and inconsequence (e(t) (0- - (t)) - L(r). by definition of L(n). I E K, and |0 F (t)l2 >
then
[0-F-00] 1) (n) because OF.(1) is star-shaped, given that K, is convex. Also 0- (+)
(1F'((t)
0 [0]](00] 17 17. Hence,
((t) (
[10-F* ( - )
0 F*(1) (0-F(())
< L(n).
which implies,
le(t) (( < 10- F (t)ll, I(m),
42

---

for any 0 E Kn. Then we observe that by the basic inequality
||=(t) - F(t)1 [((w(t) F (t)
This implies that
1 ス 11((1) F ( (F()) FT))
Given that F(t) E K2 it follows from inequality 28 that
1 -1((()) - F S max
[F((t) F (t)ll2 1n,(m),
"
and by inequality (27).
(1((()) F" (1)12 < m
[F(()- FF (F)))
" 11 E sup
|e(1) (() (000--)))))) + =o1o)
(e)
for any 8 > 0. with probability at least --006" Thus, for any 5
[|((t) - Ft) S max |2G(1, 5)
7
with probability at least onjó. Next, by inequality 26).
/18 where G(n, 8) = E (supoc@prcoso)) ((tt) (0-()()--)))))
V 2G ())(n.8 6) 5)
172 + 2amo
V12 2078 ) max {n 206, vn(n 200 } S 205.
max
17
"
In consequence for any 8 > 0,
[1 - F" (t)ll2 Sn + 206.
(29)
with probability at least --02 /266/16 /16
Finally, observe that
B(1I((1) - F* (011))
/ (11(1) - (1) (
8) ds P (I
(0+220)
F 5) ds + s. P (11(()-F'(())))
F* (t) > 8) ds
=11 + 12.
To analyze the term / we observe that P (10(1)- F S 1, and therefore
1 5
(30)
For the term Jz. we perform a change of variables & = (1 + 2 to obtain
P (( (O11 > (1 + 200)² 40 + 200)
1 R
and by inequality (2)..
1a < 40 / e
2,2
(1 + 200)d6
(31)
Moreover,
4g
2,2 e
(n 200)dd Ano
222
15 dõ + 80² f (:0
- 8d8
6-
222
To
-2n
+ 4
<2n crfc((o) + 440-0"
43

---

where erfe(2) = 2 v J e dő. From inequality (30) 31 and the fact that 17 > 1, we conclude
that
-() (1)1)) sCi( + o²),
for an absolute positive constant C1. Finally observe that,
n7² + o² max (n², o²))
Taking C 2C1 the result is achieved.
Now we are ready to start the proof of Theorem
Proof. Notice that for all i, it holds that Et(t) (*) for all inf(a e f2} and (t = (t) = for all t sup(a a 9t. Hence.
1
CRPS F,, (F)
(F((t) - (t)
=
72
71
1
(F((t) PT
71
(R(()) (()) dt
E
n1
C C'max11,19
<
dl
C max
[ dt
where the inequality follows from Theorem
5 by noticing that (6 and Lemma iimply (26).
G.3 Proof of Corollary
Proof.
Throughout we use the notation from Definitions 5 6
First, notice that (t( = (t) = 0 for all < 9(1) and for all i. Similarly, F((t) = Fi(t) = 1 for all
t 2 V(n) and for all i. Therefore,
/ (A+())-F(()) = (aL (F((0-
(F)
V(n)
(32)
Next, define
2.(1) = i7 ()((()) m(1) + vy)) for l E [0, 1),
otherwise.
Clearly, G,(0) = F,t (1)) and Gi(1) (v(r). Moreover, recalling that for E (u<))> Mo))), we
can write
n-1
(1) Σ
aig (o)
then for / E [0, 1), it holds that
G,(t) = Σ
aijj 1 (1sg (1:
(33)
k=1
where
U(1) 9(1)
41 =
V(n)) 4(1)
for E n).
44

---

Furthermore, let
G]((1 { P7 [F ()( M(n) + vay) if E [0, 1),
0
otherwise.
Now, observe that
J ())))
[. - ((1)(())))))) - y(1)) + 8(1) F (1)(m) v(1)) + ((a)))
(((n)
(F (s) P. (s))
(34)
M(x) M()
by making the change of variable Furthermore, by Lemma
(1 (U(()) 3(1) + 901-
((G))(t) ! 00
(t)
/ (t) L/¹ [G)(t) V(m) 9(1) Ju1
[ [FT ds)
((n)
(35)
and
((G.)(()) / 0 (D(G.)())2ar = l0 (2,(0))aa
L (.(0))= V(m) V(1) Jan
((n) |8 ds)
(36)
Also, by Lemma
D(G)(() D(G;)(I)dl -- 1o D(GS1)-DDGTT00
a,(() -)(0)) = -;(()tt=- =- G.(1) GG
(37)
which implies
[ D(Ga()) D(G))(1)d
1
1960
(T(()) (1)dt
(38)
U(m) 1(1)
Combining 35
36 and 37). we obtain that
| (D(G.)() - (C))
1
B(n)
(f) (t)
(39)
U(m) M(1) Jym
However, since G; is decreasing and continuous in [0, 1). then by Lemma holds that D(G;) = (t) for all t 2 0. Thus, from (39..
L (D(G))0) (t) < M(n) 4(1) M03
U(n)
- (t) (t
(40)
Now, by Lemma and 33 we obtain that
((G.)()) Σ
ai.ji 1 Im (mm))
(41)
45

---

where
U1++11 - B(ia) 1= 1,..., 1,
mi
M(n) 9(1)
and with mo = 0.
With 41 in hand, we let Ht = D(G))(1 (t- ((1))(U(()) - B(1))) for E (1(1)) Thus, can write
H1(t) Σ
aish 1 1/.00-))(t
1=1
for all t E (1(1), M()) where to and t1 M(()) Eh- (M()+11) ((())) for all / = 1,.... 1. Also, (t) H((1 ())()))) v(1)) + 3(1)). Hence, from 40) we obtain that
(iin) ( (H+(1) FI (1)) =
[ (D(G (1)
V(n) y(1) V(a)
s 1 (9
- ())
V(m) 1(1)
and as a result,
M (F((t) (t) = / (Hi(t) - FT &
((n)
M(n - (t)) ((
(42)
since (1) for all E (()):Mo)))
Combining
(32 and we obtain,
(F((t) (t) <
1 (f:(0) (t
R
The claim then follows.
G.4 Proof of Theorem
Proof. First, by the basic inequality we have that
1 Σ (P()) (F))) S -1 1 (P F((O)
2
1=1
for all F E {RF* + (1 (1- (1- E (0, 1].. Hence, for E1. . En are independent Rademacher
random variables, we have that
IP
sup IER 1=1
Σ (((()) - F1(1)) > 2n?
sup
sup
(F
FT) > n?
(EX BEK; F10=F= (0llSa i=1
F7( Σ (F)(() (m: (0))-F(())) > 1?
sup
sup
(ER BEK ||0 (MIl
(F) - 1 - FT(())
sup
sup
T
VER DE K 10-F-((1 2=1
E
sup
sup
(F7(1) - 1unset)0
1/2
EER K K: 10|
sup
sup
1/2
HER 0 K K 10155
2
E
E
sup
sup
Σ i=1
&d (11
172
HER OEK K: 50
"
max
sup
112
IE (V)
BEK-K ||0||
46

---

where the second inequality follows from Markov's inequality, the fourth by simmetrization. Next, notice that for fixed and y, the random variables (C: 1 are subGaussian(1) Hence, by Lemma 3
n Σ
s c (' ( Vlog (K K) (K)) B,),,, || []]de
I
max LE (Mi dEK |0 ==1
sup
+ +nnylogn
for some constant
0. The claim then follows.
G.5 Proof of Corollary 2
Proof. Following the proof of Theorem 2.2 in Chatterjee 1 it holds, for ~ [00, In), that
13
we obtain that for any positive integer
E
sup
( - (()))
E
sup
gT ((0-F((t) F*()))
BEKs ( St
AEK (FFK::00=F(F)]]
(0)lsn
(43)
112
S
2 2 2 nm
1/4
21-1
for positive constant Cj. Next, let L the constant
(6).
We now choose large enough such that
Cin²
n72
S
21-1 2L
Furthermore, for this choice of I, we can choose 7/ as 11 x 1/6 such that
n²²
2m 1/4 <
2L
Thus, for a choice of y satisfying x 1/6 we obtain
((--F*(t))
sup E sER
sup OEK: 1|0 (1) (8) Snn
L'
and (9) follows from Theorem follows from Corollary [
Furthermore, the corresponding conclusion for (F())ER
Finally, notice that for some positive constant C2
J 1/1
n/4
log Nee (K [a, ba n [a, b)) n B,) 11)de
$ 2
log /eg I])
7/4
2
/og.N(e/2, (,), 11)de
30
7/4
20 2C2v/n(1
2
de
2 2C2 - 2/2
where the third inequality follows from Lemma 4.20 in Chatterjee
[13 Therefore, the claim in 10
follows form Theorem 2 by taking satisfying 77 (b a) 1/3 Vlog n.
47

---

G.6 Proof of Theorem 4
Proof. First we observe that by the basic inequality, for all E IR,
[[(() - (1]1 aa()) F* (t)) + (F(een(FF(tt) (F(t) n), where a(t) wxt) (t) for all t E R and i {1, (i),..,n., and where the inequality holds for all
(44)
F(I) A(I) = {T(()) + (1- (1-s) SE 11] CR.
Therefore,
(F() [[((t) - (t)11 ((t)T(F(() - (t))
pen, (F())
+ pen, (Fn, (t))
de
for all F(t) E (t) and t R. Hence, by the properties of pen,() we have that pen, (F - (t))
(F(tt) + pen, (t)) a(t) (F(t) (+)) 2pen, (F (F*(tt) (t))
(45)
for F(I) E A() and I R. Now suppose that there exists F(f) E A(t) such that
[[((() (1)l s 1/
and pen, (F(t)) 2 (F Then
pen,(F(t) - (B))
(F(t(F(t)) (F (t)) (F (t))
Hence, we let
8: 4 penu (F* (pen)(F*()))
pen, (t)) (0
Then set
F(() = sF()) + (1 (1- (1- F)) E A(L)
As result,
| (t) F(O)12 = (1)(())F(()) F() sn²
Also,
(F Ft) (sF((sF(t) + (1 (1-8) s)F - (())
(sF(() sF)(t)
spen, (F(t) - F Ft)))
(F"
Therefore,
4pen, (F (t))
pen, (F(t) (F(Ft)F((tt)
a(t) (F(t)) FF)) 2pen, ( (F* (1))
d
where the inequality follows from 45 This implies that
2pen(F*(t)) - a(t) (F(t) Fr
(46)
X
Hence, we let
n
X 4pen(F* (F)
Then from (46 we obtain that
n² < a(t) T (F((t) - F
48

---

It follows that the events
Ω
U ER
sup P10EEN():1P00)-00000
(F()))()) 2 (F
and
sup Eelt F())A((): (f)))) sup pon(F"00--FF))) a(t) (FF() (F(t)-F((t) 2 1/²
92
2
satisfy that 2 C22 Hence, (P()) < (92))
Next, we observe that if [() > 7/. then there exists F(1) E A(1) such that |IF()=F((()]
7.. This implies, by (44), that
7 S (F())FFF + du (p) pen pen(((()))
2
and so from our choice of d
n¹
S (()))(F()) (()(t)=F((t))
(47)
4
Thus, 47 holds for some F(t) E A(t) with [F(t) (t) - provided that [|((t) - F* >
Therefore,
IP sup [F((t) - (+)l > >1)
x P
- - > 77 n + P(1))
sup
SER
TEB
P
sup
sup pen (F00))
(F(t) FF(
HER
P((22
S 2P
sup
sup (Fppe))))))
(F(t)-F*(0)
IEk (()EEM(:11P((0-
E
sup
sup IER F(O): (FFF(O)is (O)
(F(t)
172
V
E
sup
sup
a(t)
EER
#EK:
I|
(48)
Hence, for E1, ......:55 En
independent Rademacher variables independent of y,, have that
sup []((1) - (000-F((M)) F > n)
<
16
sup
sup
Σ & (l)s0))
1ER
1/2
HER 0EK
16
(1(( v))
E
E
sup
sup
112
AEE OEK
|0llsn i=1
16
Σ E.l (use00 ())
n
maxx
sup 8EK
17
where the first inequality follows from
48
and symmetrization. However, by Lemma
J n/4
max tE (V
sup B [10 S
E
< c
/logN(e Kn K (,()), dl)
+ Cnylogn
for some constant C > 0. The claim then follows.
49

---

G.7 Proof of Corollary 3
Proof. We start by noting that the set K = f0r Rº: TV) S V} is convex. In- deed, let ,20)0(3) 022) E K, and define 0la) = a0l1) + (1 a)) for E [0, Then (p())))))) n- | = n | + - a By convexity of the 4-norm, (g(a)) z aTV( + (1 (1- 2 V,. Hence, (la) E K, proving that K is convex.
Let g ~ (O, [x) and set
K = {
R" R TV [v)
(49)
Hence,
(0- (1))
E
sup BEK. 110-F (0)lsn
$ E
sup
T F*)))
BEK 1
(0 1/2r
(50)
S
C,
"
mI logn
where C, is constant the second inequality follows from Lemma B.1 in Guntuboyina et al.
23
Next, notice that
C.n VvVV 1/2r
n?
2L
holds if
(2LC,) (n (22))) S
Also,
17²
C,nylogn logn S
2L
provided that 2LC,vllogn $ Hence, taking
max (
2LC, Vognn
we obtain
The claim for F(t) follows from Corollary
To show 14 we observe that by the proof of Theorem B. in Guntuboyina et al. 23 for 0 < 1) < n, have
1/4 La
B,,0), |l I1)ss 2 Nog (K K) (,)
Vlog (e/2, n B, || 1l)
1/2r
S Ōr
vTV
+
Mylg
for positive constant č, Hence, (14) follows with the argument above.
G.8 Proof Theorem
Proof. Notice that for any e IR we have that
[[F(() - H()) < 2nB
for all F(t), (t) e Ke, This implies that suPtER [I) G(t)ll 2nB. To control the probability of large deviations, we restrict attention to the event
[(()) G() > r2
sup [F(1) G(t) s 2n1) tER
sup tER
n
50

---

which confines the deviations to the range (1²2 Bn We then partition this interval into dyadic subintervals (2 1") for 1, I,....J,, where J the smallest integer such that 2"17² 2 2nB. Applying a union bound over these dyadic shells yields
sup [I(()) - G(DII 11 tE
IP
P
sup EER
[()) G(I] > n/, sup
[F(()) G(I) - 2nB)
FER
[[((t) C(t)]? > 2' 1/
[|((t) - G(t)]l $
sup YER
sup
tER
n
Σ sup 2
Σ (F:(1) - ((nss) -f())) > 20-11?
YE
[(1)) - G(t)]2 & 2'12))
sup tE
<
sup
sup
Σ
(G))))
==1
IER
2=1
2 212
where the first inequality follows from the union bound and the second from the basic inequality. However,
Σ (G.(() - (mmss))(((()) - R((t)) > 21-
sup LER F(DEK: G(0))
sup
1=1
E
sup
sup
[G.(() - ((nss))(((() -i((t)
2 2
1ER P(DEK)
1=1
21 212
1
IF.
sup
sup
(F
F.(1))
(ER F(t)EK -F)) G(00))250 i=1
1
Σ
(G,(tt) - ()) -
2 277 1ER F(t)EK C(00) i=1
sup
sup
20 21
E
sup
sup
(F) 1 ((ss))))) (G((()-F((())
IER F(()KK (())KK;]P(())(UU G010-50 i=1
4v7i
sup |[G(I) - F
2//2n 21
tel
for some constant C > where the first inequality follows from Markov's inequality, and the last inequality holds by Cauchy-Schwar inequality. Furthermore,
E
sup 1EB F()) ((00)5500
sup
1=1 (1) (,s
EE
sup
sup
(t) ((nS)))((())
tER (TtJEK,KK,-1 F()))222
1=1
2/22
< C
/logN(e, -K-K) (0), IIdde + C22/21//oogn
51

---

where the last inequality follows as in the proof of Theorem 2 Therefore,
2//
IP
sup SE
[IF(t) - G(t)ll n) s
25 21
C
log N(e (K (K) ll IDd=
0
4vn
+
C21/2 1) / logn
sup tE
- |G((t)
21
C J
22/3///4
log N(,, (K K) n (,)) 11 I1d= +
2
jml
4C 4Cvlog T
4
sup FF EER
Σ j=1
21/2
21/2
11
"
j=1
C
J
2//21//4
log (e, (K K) n (,) (1-0),1--)dd-------------- 11)ds +
1/
21-2
j=1
4Cvlog
2-1/2
4yn
2-1/2
sup [[G(t) F*) 1 2-1/2
2-1 1/2
tek
and the claim follows noticing that
sup [F() G(\II > n) IER
P sup IEB
[[(()) (1) > + (1IP()) -(0)))))-(())))))
IP
52

---

G.9 Proof of Corollary 4
We begin by restating the corollary be proved.
Corollary 8. Let F(t) be the estimator from not necessarily in K2. Suppose that Assumption described in Appendix holds. Let
(3) with the set K, as in (16 for all t E with (t)
22 n P
On
maxx (p.M)EP
There exists positive constants e₁ and 0g such that if
I. = [ei log n] and
02 max
nE 212 37y
(51)
().
or
L = ci max
4 n 3 log n
and v == [022
(52)
(p.M)EP
then
(Ri(t) (F(0)=F(())*= Or ( (logn
sup
on log' m)
(53)
ER
Then the proof of such result is provided.
Proof. Throughout the proof, condition the covariates 's and omit this dependence from the notation for simplicity. We proceed by using Theorem 3 First, for a vector v E IR,, we let ||v|| = /oll/vn. Then, by Theorem 3 in Kohler and Langer (3), it holds that
sup (F*(() G(t)llo S Civon. IER
(54)
for positive constant C1.
Furthermore, as in the proof of Theorem Zhang et al. 67 see also Lemma 19 in Kohler and Langer 37. we have that
log N(E, v). (1x) & L3v²
for E (0, 1). Therefore, for e (0,2v/). we have that log ((L,F || II) S L' log( Vvn).
Therefore, for 17 Vn. with
log (2m/1²)) log 2
J =
holds that
25
Cylogn "
+ CVñ
IIdds +
sup [[((t) (t)l eR
2
"
n//2
Vog N(e, K, |1 -1]de Cvlogn CVn
sup [G((t F (t)
"
ER
2//4
V Nlo I1)dd + CCIVlog Cynon
2//2
Vlggn
Vnenn
W
log N(E/2, 1ide
"
"
21/2
(2222og((V)) log( Viide +
Vlog v
y nio
2
==1
Lv Luvlog((L) Lv
7/4
Vlog n
/ n
Vlog(n) + de Id- +
2
11
1
53

---

Moreover,
LV Luvlog(Lyy Lv)
2//9
Vlogn n
Vnon
Vlog(n) + 4e +
1/
20-2
0
1
==1
Ia L LwV) 172
2/13
Vlogn
Vnon
(Vlog(n) 2E
1--11/dde
11
LN Logloo(TV) LT))
Lu Loglog(Ly
21/2n/4
log n
V Vnon
2c
11/2
+
21/2
n
2
17
Luy /lglog(( Lv) logn n
Lu Lu Lv)
1
21/2n/4
n/4
Vlogn
y vnon
2c 1/2de
7
nf
20-2
7
7
Log/oogLLu)oog
Lu Loglog(LL
/1I
V log n
Vnon
24
112
231/4-2
11
Luv/lg((Ly) log
Luvlog(Lyy (L))
Vlg
1/3/2
"1
(55)
Moreover, by our choice of and v,
Lv = (logn) vnon-
(56)
Therefore, from 55 and (56)
C
2//4
v No, K(n). || Idd + Cvlogn
CVn
sup [G(() - F (t)lloo
2--2
2
" tER
1=1
n log²
vo log(n)
Vlogn
Vnon n
7
1/3/2
(57)
Hence, the claim follows by taking
n = Vlogn n + v Vno log? n.
54

---

G.10 Proof of Corollary 5
Proof. By Theorem 1
it is enough bound
g0 - (0))
sup E KER
sup |EK=:110-P<0015::
where where g (V. Let, E R and g NO In) be a standard Gaussian vector in R", and let Kr C. R be a constraint set. Now we analyze the quantity
(g. 0 F (1))
E
sup P
(01ls
Since K₁ is convex cone it is closed under positive scalar multiplication; that is, K. K₁ for any ? > Using this, have
F((t)) B(..1)
(K, - F (t) =
In consequence,
E
sup A
gg.0 (*))
=q.
sup DE (E
(g. 03
(58)
Sn
(=()) )oB(..1)
Now, we consider the tangent cone convex set Kr at point F*(1) E K₃. This is defined as,
(t)
cl { ( (0- ('(t) (000))))) :0 kn,
Tk.
11
Now, observe that K, = {0 E Rn 01 S ... S 0n) is closed convex cone, and set K; F
contained in the tangent cone K₁ at (() that is,
(1)
(t)
K,
C TK.
11
71
Therefore,
E
sup F*())
(g
SE
sup 9ETR, F=(()) 11
(g.0)
(59)
AE Ki
(BB(0.1)
Furthermore, sup WeTk.
(+ IB (g, 0) = SEP
(F(()) nss (g, 0), where 5 denotes
the unit Euclidean sphere. Then, by Equation 59
sup r
(g 0)
S E
sup OETH.
(g.0)
(60)
6E
K.
B(U.))
(() nga-1
Now, we define the Gaussian width of subset
(AC) == E
sup
(0)
LHEKnSn-1
Therefore, from Equations (58) and 60
(t)
E
sup
(g. 0 - (0)F((())
Sn w TKe
(61)
PEK. | e- F- (ollsn
Let 8(K) denote the statistical dimension of closed convex cone K. RR defined
8(K) = EgD(g; (C),
where D(y; = Σ
0,(,n; K)
Ny.
55

---

and g .V (0, Here, (yy; denotes the Euclidean projection of E R onto K, formally defined as
((u KC) arg min ||0 yl1². BEK
Using Proposition 10.1 in [T. the Gaussian width is controlled by the statistical dimension:
w (TK. (F(1) $ Vo (T(TK Tk.
(F(())
17
We now analyze the statistical dimension of this tangent cone. Let
k(1) = ( E (1....,
(T) < F (4)
be the number of strict increases in F* (t). Let 0 in < il < ... < (x(() < (x(0)+1 n denote the jump points such that (t) is constant over each segment {iz-1 + 1..... ie}. Then, the tangent cone admits the decomposition
+(1)+1
(t)
Tk.
M
where ny == 1g
t=1
with each M decomposition follows from Remark 2.1 in Chatterjee et al. [1]. Therefore, the statistical dimension satisfies
(0
Re
01 S .... 0nz) denoting the isotonic cone in T This
(TK., TK.
(F())=
k(0+1
Σ 8(M")). E=1
From Example 2.2 in 12 have that 8(Mm)) I log(em). Moreover, applying Jensen's inequality, we conclude
(F(())
en
8 TK.
= (1 + k(t)) log
1 k(t))
1
Hence, the Gaussian width is bounded as
F. (F()) )) & (1
w
Tk.
(1+ k(t))
log
eni I+ k(t)
Therefor, taking sup
(1 + k(t)) log
em (I+(()) from Theorem and Equation 61 we conclude
IER
that
E
CRPS(F,,FF)) F())
Csup IER
1 1+ l(())
log
en
n
+ k<tt)
56

---

G.11 Proof of Corollary 6
Proof. By Theorem 1
it is enough bound
g0 O F0(1)
sup E EER
sup P
where where g (N, L.). Let, R and (0, In) be a standard Gaussian vector in R2 and let Kr C R be a constraint set. Now we analyze the quantity
E
(g. - FM (1) (
sup
To that end, remember that
K2 0 R" |
V
n=-1
and using Definition 4 we have that supae FF (g. F (*)) is the same as the Gaussian complexity of E K; (t) |00 S n). this
R K; F* (t) ||||| n)),
see Equation (63 for such notation. Then, following the same line of arguments as in the proof of Theorem 2.2 see Equation (129) in their Supplementary Material, follows that
R e K; (t) |||| S n)) sC
(s + 1)
s s+11
and the desired result is obtained.
57

---

G.12 Auxiliary Lemmas
Definition Let K C E and 0. subset of K is packing of K if P and the ser (Bi(z))eet DEEP is pairwise disjoint. Then, the t-packing number of K, denoted as (t. K, -I II). is defined as the cardinality of any maximum t-packing.
Lemma 2. [Variant of Dudley's inequality. Let S C R be a finite set and c5)) E B be vector of mean zero independent SubGaussian random variables, for m. Suppose that O E S and E S implies that ||||| I Da/22for some Dr, > O. Then there exists constant C >> 0 such that
max max TT (6)) Sc CC(D.Vlogm
Du/A p
E
Vlog M (,,5,11 dō
j= wES
where M S1 -II1-1) is the packing number of S with respect to the metric || Il--
Proof. For 2 1, let Sn be maximal 2 2,,2-'- seeparaed subset of S, i.e.,
min ||0 u|l > D., B.88
By construction, [S|] M(D S, II). Clearly, because of the maximality,
max min ||U ull S D,,²" vE8 wES:
Furthermore, St for large enough l. Hence, we let
N min (1 2 1 S = S}.
Also, for 1 2 1, let be the function that assigns v E S to the point in S closest to v. By definition,
(T) e s D2
for all a E S and E N. We also write So = [0} and so Tw( = 0 for all v E Next, we observe that
v Te()) (m(()) मパ (i)
for all U E S. It follows that
N
max max U T (1) j= WES
Σ (m) - (-1(r)) l.....m WES
max max
max (m(v) ( T VES
max
and so
N
E max max T Tc0) j=I.... WES
max max (T((U) (1)) Tcbi)
E
jl..... wES
1=1
However, notice that for all 0,
2²
P ((((0)-=-1(0))) (1) 2 1)
< 2exp (1) ग. [w)
and
[m(()) ---((=) [| w|l + || w|l 3D,2-
which implies, by the subGaussian maximal inequality, that
3CDn
E
max max (I)(w) 3 (x))
s 2 Mog(2m|Si][Ss-)
j= vES
3CD, Vlog(2m|Si(3)
21
3CD,
/log(2mM(D22 II))
2
58

---

for some constant C > O. Therefore,
N
Dn 2
E
max max v (Te6)) m vES
3C'
logm + (Dg II))
1=1
3C log
Da
3C
Σ
N
Dn
log
(2M ID)
2
2
1=1
1=1
Dn
3CD, Vlgm + 3C
Σ
V /log (2,
2
D/2
3CDn, / logm + 6C
Vlog (2M (2,
Dn./2it1
D/2
= 3CD, log + 6C
log (2M 1, 11)drr
D./2N11
D=/2
3CD, / logm + 6C
Vlog (rM(r, S, 1))drr
Dn/4
< 3CD, v llgm + 60
log (2M(r, S, 1))dr
Dz/4
+6C
Vlog (2M
2 8, 11))dr
Dn./4
3CD, 3CD, Vlogm 12C
Vlog (2M (2M S, 11,11-111 (,1-,1111))) 1D)
Dn/4
3CDnVlog + 24C
Vlgg (M Sr l))dr
and the claim follows.
Lemma 3. With the notation and conditions of Lemma then
8 E arbitrary (not necessarily finite),
max (up SC (DVVlgm j= REg
(or Da/1
v log M (6, 8, -1168
Proof. Let 5 S be a countable set such that
max sup u T(() max sup v (()). j= vES
j= weS
Let Ši the set consisting of the first elements of S. Without loss of generality let's assume that 0 E for all I. Then by Lemma
2
we have that
E max sup u oTe() j= ve
S C Dr DnVlogm So Da/4
V log M (s, St, |l 11)
hlog
(62)
Da/4
< C Da DaVlogm-
Vlog M (8, S. Mod
for all Hence, the claim follows by letting 1 in 62 and applying the Monotone Convergence Theorem.
L
Definition 4. For & TR let <(t) w(() F*(1) RR
Then for set V C RR define
oies (eE(())
R&(V) = E
sup
VEV
1=1
59

---

Furthermore, define the Gaussian complexity V V as
R(V)
E sup
(63)
where g (0, l).
Lemma 4, There exists a universal constant such that for any set V C R it holds that sup R₂(V) S IR(V),
tER
where L is universal constant.
Proof. Fix t E R. Then for v E V, let Y, v c(t) and X, = _ g. Next, observe that
sup u.eeV
y Yul
sup u.eeV
(Ya Yy) supY,, + sup Y
eEV
Hence, for any vo V,
E
sup Y - Yol
=
supY, u WEV
+
sup Yo BEV
sup u WEV
+ (-Yoo)
(64)
2
supY, WEV
R (V).
Now, we observe that e1( is sub-Gaussian() Hence, by Theorem 2.1.5 from Talagrand 55 we have that
E
sup
Ya
-Yol
LE(sup T 9) LR(V).
wEV for universal constant that does not depend on Hence, the claim follows.
u..EV
Definition For measurable function f + R we define its distribution with respect to the
Lebesgue measure the function P [0, oc) - RR given as
(y(A) p({= If(x)) A})
where y is the Lebesgue measure.
Definition For measurable function D(f) (0, + R given as
+ R we define its decreasing rearrangement
D()()) inf 0 1y(2) 5
where μ is the Lebesgue measure.
Lemma 5. Suppose that I can be written
n-
f(t) = Σ (t)
for bi .... bn-1 and measurable sets Eq C R that are pairwise disjoint. Then the decreasing rearrangement of f is given by
n-1
D()(1
Σ bel mm) (1)
where
mi Σ p(Ex), 1 I,...,
and with mo 0.
60

---

Proof. This is Example 1.6 in Bennett and Sharpley [6]
Lemma 6, For any integrable functions f and g the following hold:
1.
110)) I. (P()(t))
(65)
R
2. [G.H. Hardy and J.E. Littlewood].
J (5())()))
{o- (() (s)ds.
(66)
3. Suppose that is decreasing and continuous in [0, a) for some 0. and f() otherwise. Then [(1) = (()) for all E [0, oo).
Proof. The claim in 665 follows from Proposition 1.8 in Bennett and Sharpley (6. The inequality in (66) is the well-known G.H Hardy and J.E Littlewood inequality, see for instance Theorem 2.2 in Bennett and Sharpley [6].
We now prove the final claim. Let 2 Suppose that E [0, Then
(()(t)
inf ( 2 sup x : > inf A 2 0 sup x 2 0: f(z) > X} st}.
(67)
Hence, for 0 & X < S(),, the continuity of f in a), implies that there exists l' E (1, a) such that f(t') & f(t). Thus, supf# x 2 0 f(x) > X} > t. On the other hand, if (0) 2 À 2 f(t). then,
also by the continuity of in (0,a).
sup(# (r 2 0 f(x) > \ = i= 20: f(x) A}.
Hence, from (67) we obtain f(t (t) for E (0,a).
Suppose now that & E co). Then
(y S a ast.
Hence, (t))(t) 0 and the claim follows.
61


Reviews:
[{"review_id": "Jizf8WoX9e", "reviewer": "Reviewer_c7JP", "summary": "This paper studies the problem of estimating multiple cumulative distribution functions (CDFs) using samples from those distributions under structural constraints. Specifically, given $n$ points sampled from $n$ distributions, the authors propose a distributional regression method which minimizes the averaged CRPS scores between the estimating functions and the one-sample empirical CDFs under functional constraints. They show that the solution corresponds to the $\\ell_2$ projection of the empirical CDF values to the constraint set. Upper bounds of the CRPS and mean squared error (MSE) between the estimator and the true CDFs are proved under various settings of the constraints (convex and non-convex, containing the true CDFs and not containing the true CDFs). They also apply their method to the specific case of isotonic regression, trend filtering, and dense ReLU network by identifying the corresponding constraints to establish risk upper bounds. Finally, they compare the performance of their distributional regression method with other existing methods on simulated and real-world datasets.", "questions": "1. It would be better to explicitly summarize the existing bounds to compare to the proved bounds, preferably in a table. How does the general results in Theorem 1, 2, and 3 compare to existing results of CDF estimation? For Theorem 1, it seems that the results follow from Theorem A.1 in [23]. Then, what is the technical challenges of this paper considering the results in [23]? What improvements have been made comparing to the results in [23]? I would appreciate more discussion on the comparison to [23] in the paper. \n\n2. Are there lower bounds for the considered problem? Without lower bounds, it is unclear how tight the upper bounds are. \n\n3. I do not find explicit statement that $K_t$ is convex in Section 3. Please emphasize on this point in Section 3. \nI am curious about the convexity of $K_t$ in eq. (12). Is it because TV${}^{(r)}$ satisfies triangle inequality? If so, please include a proof or a reference. \nBesides, I would recommend the authors to explain the benefits of convexity in the proofs and how the proof for convex constraints extends to nonconvex constraints. \n\n4. For Theorem 2, $N(\\epsilon, (K-K)\\cap B_\\epsilon(0))=1$. It seems from the proof that this quantity should be $N(\\epsilon, (K-K)\\cap B_\\eta(0))$. Similar issue exists for Theorem 3.\n\n5. In the conclusion, the authors claim that \"structured constraints\" leads to \"improved estimation accuracy\", but it seems that those results under the listed structured constraints are just direct corollaries of the general results. Then, why are there improved estimation accuracy?\n\n6. For the experiments, please explain more on how to apply the proposed method to test data. For training, we can just estimate CDFs of those training variables. Then, how do we use the trained functions in test data?\n\n7. There are some typos and unclear parts. (1) In line 42, since $F^*(t)\\in\\mathbb{R}^n$, we should have $K_t\\subset\\mathbb{R}^n$. \n\n   (2) In line 132, it should be \"for some fixed\". \n\n   (3) Do we need the assumption that $F^*(t)\\in K_t$ for all $t$ for Theorem 1, Corollary 1, and Theorem 2? Is so, please make this assumption explicit in the statement whenever it is needed. \n\n   (4) Does $C$ in line 132 depends on $\\Omega$ and other variables? If so, please make this explicit. \n\n   (5) In line 153 and 158, what is $K-K$?\n\n   (6) Please explain the notation $O_{\\mathbb{P}}$. \n\n   (7) In line 245, a wrong spelling \"unctions\". \n\n   (8) Please explain the notation $\\mathcal{P}$ in line 253.\n\n   (9) Please add a reference to the \"hierarchical composition class\" in line 257.", "limitations": "Yes.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 2, "strengths_and_weaknesses": "I do not check the details of the proofs, but the paper seems to be technically sound. Both theoretical and empirical results of the paper are complete. \n\nThe paper is clearly written and well-organized, but it would be better to include some technical overview of the proofs and more details of the experimental setting. See my questions below for details.\n\nThe results of the paper are significant for CDF estimation and its downstream application. It allows various function classes through the flexibility in the choice of the constraints and thus considerable potential for application. It would be better to make the comparisons with existing theoretical results explicit by summarizing their bounds. Since lower bounds are missing, it is unclear how tight the results of this paper are. See my questions below for details.\n\nThough the idea of distributional regression based on one-sample empirical CDFs are not new, the authors provide a general and practical recipe for this problem with various application settings. The theoretical work of this paper is based on existing results and standard techniques of covering numbers and uniform convergence. \n\nFor other weaknesses, please refer to my questions below.", "was_revised": true, "final_justification": "I appreciate it that the authors could revise the paper in response to my suggestions and questions. As more comparisons with existing work on CDF estimation and explanations of the technical challenges compared to the existing theoretical results are needed to revise the current submission, I will maintain my current score.", "rebuttal": "We thank the reviewer for careful reading with valuable feedback. Below we would like to address the comments and questions raised by the reviewer point-by-point.\n\n\n---\n\n**Q1 & Weakness 2: Comparison and Novelty Relative to Prior Work**\n\nThank you for all your valuable comments. We expanded and clarified the theoretical comparisons in the revision. Table~2 summarizes the convergence rates and assumptions (partially shown here due to space constraints).\n\n**Bounds**\n\nDue to space constraints in this response, we briefly comment on the new comparisons we have added.\n\nOur isotonic result and Theorem 3.3 of Mösching and Dümbgen [44] are compared. Translating their setup into our notation, their quantity $F_{X_i}(t) := \\mathbb{P}(Y_i \\leq t \\mid X_i)$ corresponds to our $F_{i}^{\\ast}(t) := \\mathbb{P}(y_i \\leq t)$. The monotonicity assumption on $x \\mapsto F_x(t)$ aligns with our use of an isotonic constraint on the sequence $F_1^{\\ast }(t) \\leq \\cdots \\leq F_n^{\\ast }(t)$. However, our result does not assume Hölder continuity for the map $x \\mapsto F_x(t)$, or dense covariates in its domain. Despite these weaker requirements, we achieve a faster rate of $n^{-2/3}$ (up to logarithmic factors) for both average CRPS and worst-case MSE, compared to their $n^{-2\\alpha/(2\\alpha+1)}$ rate for the worst-case MSE, where $\\alpha \\in (0,1]$.\n\nWe compared our Dense ReLU estimator with DRF [10] and EnG [57], two methods explicitly designed for conditional distribution estimation.\n\nDRF estimates the conditional CDF using forest-based weights to form a weighted empirical distribution. EnG instead models $Y = g^{\\ast }(X, \\varepsilon)$, with $\\varepsilon \\sim \\text{Unif}[0,1]$, and trains a neural network via energy score minimization; the estimated CDF is obtained by sampling from the learned model. Translating their notation into our notation, both methods assume that $G^{\\ast }(x, t) := \\mathbb{P}(Y \\leq t \\mid X = x)$ is Lipschitz in $x$ for each fixed $t$, and that the covariate density is bounded away from zero and infinity.\n\nOur method instead assumes $G^{\\ast }(\\cdot, t) \\in \\mathcal{H}(l, \\mathcal{P})$, a broader class that includes the Lipschitz case ($\\mathcal{H}(1, \\\\{(1, d_0)\\\\})$) and mitigates the curse of dimensionality, which typically affects Lipschitz-based methods in high-dimensional settings. Moreover, we do not require assumptions on the covariate density. While DRF and EnG show convergence in probability at fixed $t$, our method yields high-probability bounds on the worst-case MSE uniformly over $t$, with explicit rate $\\frac{\\log n}{n} + \\phi_n \\log^4 n$.\n\n**Theorem 1-3 in CDF literature:**\n\nThank you for the question. Theorems 1–3 adapt classical techniques from regression, where projection-based estimators have been studied in both convex (e.g., [5, 14]) and non-convex (e.g., [46]) settings. These prior results focus on controlling mean squared error at a fixed point and were not developed for CDF estimation, which requires uniform control over all $t \\in \\mathbb{R}$. To our knowledge, our extensions provide the first such results for CDF estimation.\n\n**Clarification regarding Theorem A.1 in [23]:**\n\nOur result in Theorem 1 is conceptually related to Theorem A.1 of [23], but technically distinct. In their setting, the observations are Gaussian $Y = \\mu + Z$, with $Z \\sim \\mathcal{N}(0, \\alpha^2 I_n)$. The estimator $\\hat{\\mu}$ is the Euclidean projection of $Y$ onto a convex set $K$. Their goal is to bound the squared error $\\vert\\vert \\hat{\\mu} - \\mu\\vert\\vert_2^2 = t^{\\ast }$ in expectation, which is reduced to be controlled with high probability. For this, the \"KEY\" step is to show that the squared error $t^{\\ast }$ concentrates around a deterministic benchmark $t_\\mu$, which characterizes the scale of estimation error. This is proved using a peeling argument that combines the strict concavity of an associated process with Gaussian concentration.\n\nOur setting also relies on convex projection, but differs in several important ways. We observe Bernoulli vectors $w(t) = 1_{\\\\{ y_i \\leq t \\\\}}$ (their $Y$), with sub-Gaussian, heteroscedastic noise $\\epsilon(t) = w(t) - F^{\\ast }(t)$ (their $Z$), and our target is the CRPS rather than squared error. Despite these differences, our analysis also hinges on addressing the same \"KEY\" step: showing that the squared $\\ell_2$ error $\\vert\\vert \\widehat{F}(t) - F^{\\ast }(t)\\vert\\vert_2^{2}$ (their $t^{\\ast }$) concentrates around a complexity-driven benchmark $\\eta$ (their $t_\\mu$).\n\nHowever, our proof takes a completely different route. We use (i) a sub-Gaussian concentration inequality for Lipschitz and separately convex functions to control a local empirical process, and (ii) a star-shaped self-normalization argument to extend this control to the full constraint set $K_t$. Because of the differences in observations, loss function, and techniques, our result is not a direct application of theirs.\n\n**Q2: Lower bounds**\n\nThank you for the question. Theorems 1–3 are general results that apply to a wide range of constraint sets, both convex and non-convex, and are not tailored to any specific estimation problem. For this reason, we do not provide lower bounds at the general level.\n\nHowever, in specific examples like isotonic regression, or total variation, the rates we obtain using these general theorems match known minimax rates for mean estimation with Gaussian noise. In this sense, the bounds we obtain are optimal for these concrete examples.\n\n**Q3: Convexity discussion**\n\nThank you for pointing this out.\n\n**Statement of convexity in Section 3:**\n\nWe have now added an explicit statement in Section 3.1 clarifying that the sets $K_t \\subset \\mathbb{R}^n$ are assumed to be convex throughout the section. This assumption was previously implicit, and we agree it is important to make it clear.\n\n**Convexity of set in (12):**\n\nWe added a short formal proof of the convexity of the set $K_t$ in Equation (12) in Appendix F.7, which contains the proof of Corollary 3 in the revised manuscript. As you pointed out, the key argument is that $\\mathrm{TV}^{(r)}(\\theta) = n^{r-1}\\\\| D^{(r)} \\theta\\\\|_1$ is convex. We copy it here for your convenience:\n\nWe start by noting that the set $K := \\\\{ \\theta \\in \\mathbb{R}^n : \\mathrm{TV}^{(r)}(\\theta) \\leq V \\\\}$ is convex. Indeed,\nlet $\\theta^{(1)}, \\theta^{(2)} \\in K$, and define $\\theta^{(\\alpha)} := \\alpha \\theta^{(1)} + (1 - \\alpha)\\theta^{(2)}$ for $\\alpha \\in [0,1]$. Then\n\n$\\mathrm{TV}^{(r)}(\\theta^{(\\alpha)}) = n^{r-1} \\\\| D^{(r)} \\theta^{(\\alpha)} \\\\|_1 = n^{r-1} \\\\| \\alpha D^{(r)} \\theta^{(1)} + (1 - \\alpha) D^{(r)} \\theta^{(2)} \\\\|_1$.\n\nBy convexity of the $\\ell_1$-norm,\n\n$\\mathrm{TV}^{(r)}(\\theta^{(\\alpha)}) \\leq \\alpha \\mathrm{TV}^{(r)}(\\theta^{(1)}) + (1 - \\alpha)\\mathrm{TV}^{(r)}(\\theta^{(2)}) \\leq V$.\n\nHence, $\\theta^{(\\alpha)} \\in K_t$, proving that $K$ is convex.\n\n\n**Benefits of convexity and extensions to nonconvex settings:**\n\nThank you for this helpful comment. Convexity plays a central role in both defining the estimator and analyzing its error. In the convex case (Theorems 1 and 2), the convexity guarantees that interpolation between the estimator and the truth remains within the constraint set. This allows us to apply a basic inequality and focus the analysis on a neighborhood around $F^*(t)$, leading to clean bounds based on sub-Gaussian complexity in Theorem 1 and on entropy in Theorem 2.\n\nIn the nonconvex case (Theorem 3), these tools are not available. The proof instead uses a nearby point $G(t) \\in K_t$ and separates the total error into two parts: the estimation error around $G(t)$, and the approximation error between $G(t)$ and $F^*(t)$. The estimation is then controlled using peeling. The final bound includes an extra approximation term that does not appear in the convex case.\n\nWe have added a paragraph in the revision to clarify these points.\n\n**Q6: Applying the Method to Test Data**\n\nThank you for the question. We have clarified this point in the revision. For the Dense ReLU estimator in Section 4.2, the network is trained on the training set, and predictions on test data are obtained by evaluating the fitted model at the test covariates to estimate the conditional CDF at each threshold. These are then compared against the indicators of the responses.\n\n**Q5: Structured Constraints and Accuracy**\n\nThank you for the helpful comment. We have clarified the conclusion, which now reads: “In different settings, when structured assumptions are met, our results suggest that estimation accuracy can be achieved. This is the case, for instance, in isotonic regression if we assume that we are estimating a monotonic signal, and similarly for total variation.”\n\n**Q4 and Q7: Typos**\n\nThank you for these detailed comments. We have corrected the typos, clarified the notation, and made the relevant assumptions explicit in the revised version.\n\n**Q4 and Q7: Table 2**\n\nBelow, we present Table 2 partially, showing only the summary of results for the Isotonic case due to space constraints.\n\n### **Isotonic**\n\n| Work         | Method        | Assumptions                                                   | Convergence Rate                          |\n|--------------|---------------|----------------------------------------------------------------|-------------------------------------------|\n| Prior Work   | [27]          | $K_t$ monotone cone, $F^*(t)\\in K_t$                          | $o(1)$                                    |\n| Prior Work   | [44]          | $K_t$ monotone cone, $F^*(t)\\in K_t$, dense design, Hölder    | $n^{-2\\alpha/(2\\alpha+1)}, \\alpha\\le 1$ (Worst MSE)    |\n| This Work    | General rate  | $K_t$ monotone cone, $F^*(t)\\in K_t$                          | $n^{-2/3}$ (CRPS, Worst MSE)              |\n| This Work    | Fast rate     | Few number of strict increases, $k(t)\\ll n$                                         |$\\frac{1 + k(t)}{n} \\log \\left( \\frac{en}{1 + k(t)} \\right)$ (CRPS)             |"}, {"review_id": "SjPdjhl9AQ", "reviewer": "Reviewer_iGJs", "summary": "This paper presents a unified theoretical and algorithmic framework for nonparametric distributional regression under both convex and non-convex constraints. It analyzes risk bounds based on Continuous Ranked Probability Score (CRPS) and worst-case MSE, with applications to isotonic regression, trend filtering, and deep neural networks. The authors prove convergence rates and support their claims through simulation and real-data experiments.", "questions": "1. While the paper proposes a unified framework for risk bounds in distributional regression, it is not always clear which results represent new theoretical advances versus careful syntheses or extensions of existing work. For instance, the projection-based estimator (Equation 3) and its analysis under convex constraints closely follow the structure of M-estimation with shape constraints, and the CRPS risk bounds leverage tools like Gaussian complexity and entropy bounds that are standard in nonparametric theory. Similarly, the neural network analysis builds directly on Kohler and Langer (2021) with minimal modification to the underlying techniques. The authors should clarify: Which lemmas or theorems are original? Which contributions extend prior work (and how)?\n2. Theorem 3's current form of the bound may be of limited practical use. In particular, the deviation bound involves a threshold eta and achieves small probability only if eta is of order sqrt{n}, which would imply an unacceptably large estimation error. This weakens the utility of the bound, especially compared to existing uniform concentration results in empirical process theory.\n3. Corollary 4 presents convergence rates for neural network-based distributional regression. However, the result depends on Assumption 1 and specific architecture-related conditions that are only described in Appendix C and F.9. Unfortunately, these assumptions are not summarized or even mentioned in the main text, making it difficult to evaluate the generality or strength of the result.", "limitations": "The authors did not discuss the limitations or potential negative societal impacts of their work. While the technical contributions are clearly presented, the paper omits critical reflections on the following:\n1. The methodology assumes prior knowledge of structural constraints (e.g., monotonicity, bounded variation, or function class for neural nets), which may not hold in practice. The authors should acknowledge this and discuss how misspecification could affect estimation accuracy.\n2. Although projection-based estimators are efficient in low dimensions, scalability—particularly for trend filtering and neural networks—can be a concern in large-scale or high-dimensional settings. This should be addressed.\n3. The method is primarily validated on simulations and a few benchmark datasets. Real-world data often involves missing values, measurement errors, or distribution shifts. A discussion of robustness in such settings is missing.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strengths:\n1. The paper presents non-asymptotic risk bounds for distributional regression under convex and non-convex structural constraints. This includes uniform bounds over the domain.\n2. The framework can add a range of constraints including monotonicity and bounded variation, and extends to neural networks.\n3. The bounds are derived using projection-based estimators and Gaussian complexity tools, and they generalize prior results (e.g., isotonic distributional regression).\n\nWeaknesses\n(a) While the paper proposes a general projection-based estimator, this method essentially reduces to projection onto known constraint sets (e.g., isotonic cones or TV balls), which are well-studied in nonparametric regression. The core estimator is structurally straightforward and is similar to existing techniques (e.g., M-estimators under constraints).\n(b) The theoretical analysis heavily relies on known tools—e.g., Gaussian complexity and covering numbers for the convex case, and peeling arguments for the non-convex case. \n(c) Theorem 3 is quite problematic. It seems that eta should be in the magnitude of sqrt{n} to ensure that the probability to be small. This implies that the error should be quite significant. Can you improve it?", "was_revised": true, "final_justification": "The referee appreciates the authors' efforts. The original concerns are still valid. It would take efforts to address them. The paper is solid and has its own merits. Nevertheless, the referee retains the same positive score for this manuscript.", "rebuttal": "We thank the reviewer for careful reading with valuable feedback. Below we would like to address the comments and questions raised by the reviewer point-by-point.\n\n---\n\n**Q1 and Weakness (a) and (b): Simplicity of projection-based estimator and use of standard theoretical tools**\n\nThank you for the thoughtful comments. We agree that projection-based estimators onto known constraint sets, such as isotonic cones or total variation balls, are well-studied in nonparametric regression under squared loss. Our contribution lies not in proposing a new estimator form, but in adapting this framework to the setting of CDF estimation.\n\nIn particular, Theorems 1–3 build on classical techniques from regression, where projection estimators have been studied in both convex (e.g., [5, 14]) and non-convex (e.g., [36]) settings. These prior results focus on pointwise mean squared error control and were not developed for CDF estimation, which requires uniform control over all $t \\in \\mathbb{R}$. To our knowledge, our results provide the first such uniform risk bounds for projection-based estimators in the distributional regression setting.\n\nWhile the analysis uses standard tools, such as Gaussian complexity bounds in the convex case and peeling arguments in the non-convex case, the novelty lies in the adaptation of these techniques to the CDF estimation setting and in the derivation of uniform-in-$t$ risk bounds.\n\nMoreover, the general theorems we establish allow us to analyze concrete examples, such as isotonic regression and total variation denoising, where we show that the resulting rates match known minimax rates for squared loss under Gaussian noise. In this sense, the bounds we obtain are optimal for these structured cases, but now in the more general context of distributional estimation. For a summary comparing our convergence rates with those in the literature, see Table~2 below, which we have incorporated into the revised manuscript.\n\n---\n\n**Weakness (c) and Q2: Concerns about Theorem 3 and scaling of $\\eta$**\n\nThank you for highlighting this aspect. We clarify that the appearance of the $\\sqrt{n}$ factor in Theorem 3 is not due to looseness in the proof. Rather, it reflects the fact that Theorem 3 controls the $\\ell_2$ error over $n$ coordinates, i.e., $\\vert\\vert\\widehat{F}(t)-F^*(t)\\vert\\vert_2$. In applications, the relevant quantity is the average per-sample error, which is normalized by $\\sqrt{n}$.\n\nMoreover, the $\\sqrt{n}$ term appears inside the event in the probability bound as a multiplier of the approximation error term $\\sup_t \\vert\\vert F^{\\ast }(t) - G(t)\\vert\\vert_\\infty$. This quantity measures how well the true distribution vector $F^{\\ast }(t)$ can be approximated by elements in the constraint set $K_t$. A similar approximation term appears in related results in the literature, such as Kohler and Langer [36]. After normalizing by $\\sqrt{n}$, the final bound depends on both the approximation error and $\\eta/\\sqrt{n}$. To ensure the bound holds with high probability, it suffices to choose $\\eta$ of order $\\sqrt{\\log n} + \\sqrt{n} \\cdot \\sup_t \\vert\\vert F^*(t) - G(t)\\vert\\vert_\\infty$.\n\nThis behavior is confirmed in Corollary 4 for our Dense ReLU network example. There, we assume that $G^{\\ast }(\\cdot, t) \\in \\mathcal{H}(l, \\mathcal{P})$ where $F_i^{\\ast }(t) = G^{\\ast }(x_i, t)$, see Assumption 1 in Appendix C, and define the constraint set\n\n$$\nK_t = \\\\{ \\theta \\in \\mathbb{R}^n : \\theta_i = f(x_i), \\text{ for some } f \\in \\mathcal{F}(L, \\nu) \\\\},\n$$\n\nwhere $\\mathcal{F}(L, \\nu)$ is a class of ReLU networks with appropriate number of Layers and Neurons. In this case, the approximation error scales as $\\phi_n = \\max_{(p, M) \\in \\mathcal{P}} n^{-p/(2p + M)}$, and the final rate becomes\n\n$$\n\\sup_{t \\in \\mathbb{R}} \\frac{1}{n} \\sum_{i=1}^n\\left(\\widehat{F}_i(t)-F_i^*(t)\\right)^2 = O_P\\left(\\frac{\\log n}{n} + \\phi_n \\log^4 n \\right).\n$$\n\nThus, the resulting rate is meaningful and matches known convergence rates for mean estimation under Gaussian noise using Dense ReLU networks. This shows that our bound is both general and sharp in concrete, practically relevant settings.\n\n---\n\n**Q3: Assumptions behind Corollary 4 and missing summaries in main text**\n\nThank you for pointing this out. In the revision, we have moved the key elements of Assumption 1 and the architecture-related conditions from the Appendix to the main text to improve clarity and allow readers to better assess the scope and generality of Corollary 4.\n\n---\n\n**Limitations: Assumptions, scalability, and robustness**\n\nThank you for the helpful suggestion. In the revision, we have added a dedicated discussion of limitations at the end of the main text. This includes comments on structural constraint assumptions, scalability, and robustness to real-world data issues.\n\n\n---\n***Table 2***\nBelow, we present a summary of our convergence rates compared to those in the literature.\n\n### **Isotonic**\n\n| Work         | Method        | Assumptions                                                   | Convergence Rate                          |\n|--------------|---------------|----------------------------------------------------------------|-------------------------------------------|\n| Prior Work   | [27]          | $K_t$ monotone cone, $F^*(t)\\in K_t$                          | $o(1)$                                    |\n| Prior Work   | [44]          | $K_t$ monotone cone, $F^*(t)\\in K_t$, dense design, Hölder    | $n^{-2\\alpha/(2\\alpha+1)}, \\alpha\\le 1$ (Worst MSE)    |\n| This Work    | General rate  | $K_t$ monotone cone, $F^*(t)\\in K_t$                          | $n^{-2/3}$ (CRPS, Worst MSE)              |\n| This Work    | Fast rate     | $K_t$ monotone cone, $F^*(t)\\in K_t$, with few strict increases, $k(t)\\ll n$                             | $\\frac{1 + k(t)}{n} \\log \\left( \\frac{en}{1 + k(t)} \\right)$ (CRPS) |\n\n### **Trend Filtering**\n\n| Work         | Method        | Assumptions                                                                 | Convergence Rate                          |\n|--------------|---------------|------------------------------------------------------------------------------|-------------------------------------------|\n| Prior Work   | [41,61,23]    | $F_r := \\\\{ \\theta \\in \\mathbb{R}^n : \\mathrm{TV}^{(r)}(\\theta) \\leq V \\\\}$, $f \\in F_r$ with $f$ the regression function | $\\frac{V^{2/(2r+1)}}{n^{2r/(2r+1)}}$ (MSE) |\n| Prior Work   | [23]          | Sparse $D^{(r)}f$, min segment length                                       | $\\\\{\\max{\\frac{V^*}{n^{r-1}}, 1} \\\\}\\cdot \\frac{s+1}{n} \\log\\left(\\frac{en}{s+1}\\right)$ (MSE) |\n| This Work    | General rate  | $K_t := \\\\{ \\theta \\in \\mathbb{R}^n : \\mathrm{TV}^{(r)}(\\theta) \\leq V_t \\\\}$, $F^*(t) \\in K_t$ | $\\frac{V^{2/(2r+1)}}{n^{2r/(2r+1)}} + \\frac{\\log n}{n}$ (CRPS, Worst MSE) |\n| This Work    | Fast rate     | Sparse $D^{(r)}F^*(t)$, min segment length                                 | $\\\\{\\max{\\frac{V^*}{n^{r-1}}, 1}\\\\} \\cdot \\frac{s+1}{n} \\log\\left(\\frac{en}{s+1}\\right)$ (CRPS) |\n\n### **Dense ReLU Networks**\n\n| Work         | Method        | Assumptions                                                                                       | Convergence Rate                          |\n|--------------|---------------|----------------------------------------------------------------------------------------------------|-------------------------------------------|\n| Prior Work   | [36]          | Regression $y_i = f(x_i) + \\epsilon_i$, $f \\in \\mathcal{H}(l, \\mathcal{P})$                        | $\\phi_n$ (MSE)                            |\n| Prior Work   | [10]          | $F_i^{\\ast }(t) = G^{\\ast }(x_i, t) = \\mathbb{P}(y_i \\leq t \\mid x_i)$, $G^{\\ast }(\\cdot, t) \\in \\mathcal{H}(1, \\\\{(1,d_0)\\\\})$ Forest-Weighted empirical CDF | $o(1)$                                    |\n| Prior Work   | [57]          | $F_i^{\\ast }(t) = G^{\\ast }(x_i, t) = \\mathbb{P}(y_i \\leq t \\mid x_i)$, $G^{\\ast }(\\cdot, t) \\in \\mathcal{H}(1, \\\\{(1,d_0)\\\\})$ Neural Net-based sampling                                             | $o(1)$                                    |\n| This Work    | General rate  | $F_i^{\\ast}(t) = G^{\\ast}(x_i, t)$, $G^{\\ast}(\\cdot, t) \\in \\mathcal{H}(l, \\mathcal{P})$                          | $\\frac{\\log n}{n} + \\phi_n \\log^4 n$ (CRPS, Worst MSE) |\n\n\n---\n***References:***\n\n[5] Pierre C Bellec. Sharp oracle inequalities for least squares estimators in shape restricted regression. The Annals of Statistics, 46(2):745–780, 2018.\n\n[10] Domagoj Cevid, Loris Michel, Jeffrey Näf, Peter Bühlmann, and Nicolai Meinshausen. Distributional random forests: Heterogeneity adjustment and multivariate distributional regression. Journal of Machine Learning Research, 23(333):1–79, 2022.\n\n[14] Sourav Chatterjee. Matrix estimation by universal singular value thresholding. 2015.\n\n[23] Adityanand Guntuboyina, Donovan Lieu, Sabyasachi Chatterjee, and Bodhisattva Sen. Adaptive risk bounds in univariate total variation denoising and trend filtering. The Annals of Statistics, 48:205–229, 2020.\n\n[27] Alexander Henzi, Johanna F Ziegel, and Tilmann Gneiting. Isotonic distributional regression. Journal of the Royal Statistical Society Series B: Statistical Methodology, 83(5):963–993, 2021.\n\n[36] Michael Kohler and Sophie Langer. On the rate of convergence of fully connected deep neural network regression estimates. The Annals of Statistics, 49(4):2231–2249, 2021.\n\n[41] Enno Mammen and Sara Van De Geer. Locally adaptive regression splines. The Annals of Statistics, 25(1):387–413, 1997.\n\n[57] Xinwei Shen and Nicolai Meinshausen. Engression: extrapolation through the lens of distributional regression. Journal of the Royal Statistical Society Series B: Statistical Methodology, page qkae108, 2024.\n\n[61] Ryan J Tibshirani. Adaptive piecewise polynomial estimation via trend filtering. 2014"}, {"review_id": "A5bF9riPlV", "reviewer": "Reviewer_LnJT", "summary": "This paper investigates the risk bounds of nonparametric distribution regression estimators. By establishing general upper bounds for the continuous ranked probability score (CRPS) and the worst-case mean squared error (MSE) across the domain, it proposes methods with theoretical guarantees. Comprehensive experiments on both simulated and real data validate the theoretical contributions and demonstrate their practical effectiveness.", "questions": "1. What are the limitations of the proposed theoretically guaranteed constraints in real-world applications?  \n2. Is the proposed method effective when applied to large-scale datasets and models?  \n3. What challenges or open problems remain in distribution regression currently?\n4. Is the proposed upper bound sufficiently tight?", "limitations": "Yes.", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n1. This paper has clear theoretical guarantees.\n2. This paper  is well-organized and easy to understand.\n3. The experimental results support the theory.", "was_revised": true, "final_justification": "My main concern about this work was the limitations of the theory in the real world, and the author's rebuttal effectively addressed my doubts. Therefore, I ultimately gave a positive score.", "rebuttal": "We thank the reviewer for the evaluation and thoughtful questions regarding our work. We appreciate the opportunity to clarify and expand on several important aspects of the paper. Our detailed responses follow below.\n\n---\n\n**Q1: Limitations of theoretically guaranteed constraints in real-world applications**\n\nThank you for bringing this up. In the revision, we have added a dedicated discussion on the limitations of our theoretically motivated constraints in real-world applications. In particular, we note that structural assumptions such as monotonicity or bounded variation may not hold exactly in practice, and we discuss how misspecification may affect performance. We also comment on the scalability of certain constraint classes and potential robustness issues in the presence of noise or distributional shifts.\n\n---\n\n**Q2: Scalability of the proposed method to large-scale datasets and models**\n\nThank you for pointing this out. We confirm that the proposed method is scalable. While our estimator computes a separate projection $\\widehat{F}(t)$ for each threshold $t$, these projections can be evaluated in parallel. Moreover, the structural constraints we consider, such as isotonicity, trend filtering, and neural network-based function classes, are individually scalable and well supported by existing optimization tools. In particular, our neural network implementation leverages architectures commonly used in large-scale nonparametric regression, ensuring practical feasibility in high-dimensional settings.\n\n---\n\n**Q3: Challenges and open problems in distribution regression**\n\nThank you for the question. An open problem in distributional regression is how to extend the proposed framework to the case of multivariate responses. This remains a challenging setting due to the lack of a natural ordering in multiple dimensions and the need to estimate joint conditional distributions. While this extension is beyond the scope of the current paper, we agree it is an important direction for future research. We have added a note about this in the discussion of limitations in the revised version.\n\n---\n\n**Q4: Tightness of the proposed upper bound**\n\nThank you for the question. Theorems 1--3 are general results that apply to a wide range of constraint sets, both convex and non-convex, and are not tailored to any specific estimation problem. For this reason, we do not provide lower bounds at the general level.\n\nHowever, in specific examples such as isotonic regression and total variation constraints, the convergence rates we obtain using these general theorems match known minimax rates for mean estimation under Gaussian noise. In this sense, the bounds are optimal for these concrete examples.\n\nTo better illustrate this, we have added a dedicated discussion in the revision and created Table~2, which summarizes the convergence rates and key assumptions across our proposed estimators (including isotonic, trend filtering, and Dense ReLU networks) alongside relevant prior work.\n\nDue to space constraints in this response, we briefly comment on the new discussion we incorporated.\n\nFor instance, our isotonic result and Theorem 3.3 of Mösching and Dümbgen [44] are compared. Translating their setup into our notation, their quantity $F_{X_i}(t) := \\mathbb{P}(Y_i \\leq t \\mid X_i)$ corresponds to our $F_i^{\\ast }(t) := \\mathbb{P}(y_i \\leq t)$. The monotonicity assumption on $x \\mapsto F_x(t)$ aligns with our use of an isotonic constraint on the sequence $F_1^{\\ast }(t) \\leq \\cdots \\leq F_n^{\\ast }(t)$. However, our result does not assume Hölder continuity for the map $x \\mapsto F_x(t)$, or dense covariates in its domain. Despite these weaker requirements, we achieve a faster rate of $n^{-2/3}$ (up to logarithmic factors) for both average CRPS and worst-case MSE, compared to their $n^{-2\\alpha/(2\\alpha+1)}$ rate for the worst-case MSE, where $\\alpha \\in (0,1]$.\n\nWe compared our Dense ReLU estimator with DRF [10] and EnG [57], two methods explicitly designed for conditional distribution estimation.\n\nDRF estimates the conditional CDF using forest-based weights to form a weighted empirical distribution. EnG instead models $Y = g^{\\ast }(X, \\varepsilon)$, with $\\varepsilon \\sim \\text{Unif}[0,1]$, and trains a neural network via energy score minimization; the estimated CDF is obtained by sampling from the learned model. Translating their notation into our notation, both methods assume that $G^{\\ast }(x, t) := \\mathbb{P}(Y \\leq t \\mid X = x)$ is Lipschitz in $x$ for each fixed $t$, and that the covariate density is bounded away from zero and infinity.\n\nOur method instead assumes $G^{\\ast }(\\cdot, t) \\in \\mathcal{H}(l, \\mathcal{P})$, a broader class that includes the Lipschitz case ($\\mathcal{H}(1, \\\\{(1, d_0)\\\\})$) and mitigates the curse of dimensionality, which typically affects Lipschitz-based methods in high-dimensional settings. Moreover, we do not require assumptions on the covariate density. While DRF and EnG show convergence in probability at fixed $t$, our method yields high-probability bounds on the worst-case MSE uniformly over $t$, with explicit rate $\\frac{\\log n}{n} + \\phi_n \\log^4 n$.\n\n\n---\n***Table 2***\nBelow, we present a summary of our convergence rates compared to those in the literature.\n\n### **Isotonic**\n\n| Work         | Method        | Assumptions                                                   | Convergence Rate                          |\n|--------------|---------------|----------------------------------------------------------------|-------------------------------------------|\n| Prior Work   | [27]          | $K_t$ monotone cone, $F^*(t)\\in K_t$                          | $o(1)$                                    |\n| Prior Work   | [44]          | $K_t$ monotone cone, $F^*(t)\\in K_t$, dense design, Hölder    | $n^{-2\\alpha/(2\\alpha+1)},\\alpha\\le1$ (Worst MSE)    |\n| This Work    | General rate  | $K_t$ monotone cone, $F^*(t)\\in K_t$                          | $n^{-2/3}$ (CRPS, Worst MSE)              |\n| This Work    | Fast rate     | $K_t$ monotone cone, $F^*(t)\\in K_t$, with few strict increases, $k(t)\\ll n$                             | $\\frac{1 + k(t)}{n} \\log \\left( \\frac{en}{1 + k(t)} \\right)$ (CRPS) |\n\n### **Trend Filtering**\n\n| Work         | Method        | Assumptions                                                                 | Convergence Rate                          |\n|--------------|---------------|------------------------------------------------------------------------------|-------------------------------------------|\n| Prior Work   | [41,61,23]    | $F_r := \\\\{ \\theta \\in \\mathbb{R}^n : \\mathrm{TV}^{(r)}(\\theta) \\leq V \\\\}$, $f \\in F_r$ with $f$ the regression function | $\\frac{V^{2/(2r+1)}}{n^{2r/(2r+1)}}$ (MSE) |\n| Prior Work   | [23]          | Sparse $D^{(r)}f$, min segment length                                       | $\\\\{\\max{\\frac{V^*}{n^{r-1}}, 1} \\\\}\\cdot \\frac{s+1}{n} \\log\\left(\\frac{en}{s+1}\\right)$ (MSE) |\n| This Work    | General rate  | $K_t := \\\\{ \\theta \\in \\mathbb{R}^n : \\mathrm{TV}^{(r)}(\\theta) \\leq V_t \\\\}$, $F^*(t) \\in K_t$ | $\\frac{V^{2/(2r+1)}}{n^{2r/(2r+1)}} + \\frac{\\log n}{n}$ (CRPS, Worst MSE) |\n| This Work    | Fast rate     | Sparse $D^{(r)}F^*(t)$, min segment length                                 | $\\\\{\\max{\\frac{V^*}{n^{r-1}}, 1}\\\\} \\cdot \\frac{s+1}{n} \\log\\left(\\frac{en}{s+1}\\right)$ (CRPS) |\n\n### **Dense ReLU Networks**\n\n| Work         | Method        | Assumptions                                                                                       | Convergence Rate                          |\n|--------------|---------------|----------------------------------------------------------------------------------------------------|-------------------------------------------|\n| Prior Work   | [36]          | Regression $y_i = f(x_i) + \\epsilon_i$, $f \\in \\mathcal{H}(l, \\mathcal{P})$                        | $\\phi_n$ (MSE)                            |\n| Prior Work   | [10]          | $F_i^{\\ast }(t) = G^{\\ast }(x_i, t) = \\mathbb{P}(y_i \\leq t \\mid x_i)$, $G^{\\ast }(\\cdot, t) \\in \\mathcal{H}(1, \\\\{(1,d_0)\\\\})$ Forest-Weighted empirical CDF | $o(1)$                                    |\n| Prior Work   | [57]          | $F_i^{\\ast }(t) = G^{\\ast }(x_i, t) = \\mathbb{P}(y_i \\leq t \\mid x_i)$, $G^{\\ast }(\\cdot, t) \\in \\mathcal{H}(1, \\\\{(1,d_0)\\\\})$ Neural Net-based sampling                                             | $o(1)$                                    |\n| This Work    | General rate  | $F_i^{\\ast}(t) = G^{\\ast}(x_i, t)$, $G^{\\ast}(\\cdot, t) \\in \\mathcal{H}(l, \\mathcal{P})$                          | $\\frac{\\log n}{n} + \\phi_n \\log^4 n$ (CRPS, Worst MSE) |\n\n\n---\n***References:***\n\n[10] Domagoj Cevid, Loris Michel, Jeffrey Näf, Peter Bühlmann, and Nicolai Meinshausen. Distributional random forests: Heterogeneity adjustment and multivariate distributional regression. Journal of Machine Learning Research, 23(333):1–79, 2022.\n\n[23] Adityanand Guntuboyina, Donovan Lieu, Sabyasachi Chatterjee, and Bodhisattva Sen. Adaptive risk bounds in univariate total variation denoising and trend filtering. The Annals of Statistics, 48:205–229, 2020.\n\n[27] Alexander Henzi, Johanna F Ziegel, and Tilmann Gneiting. Isotonic distributional regression. Journal of the Royal Statistical Society Series B: Statistical Methodology, 83(5):963–993, 2021.\n\n[36] Michael Kohler and Sophie Langer. On the rate of convergence of fully connected deep neural network regression estimates. The Annals of Statistics, 49(4):2231–2249, 2021.\n\n[41] Enno Mammen and Sara Van De Geer. Locally adaptive regression splines. The Annals of Statistics, 25(1):387–413, 1997.\n\n[57] Xinwei Shen and Nicolai Meinshausen. Engression: extrapolation through the lens of distributional regression. Journal of the Royal Statistical Society Series B: Statistical Methodology, page qkae108, 2024.\n\n[61] Ryan J Tibshirani. Adaptive piecewise polynomial estimation via trend filtering. 2014"}, {"review_id": "hpfB2IRnEx", "reviewer": "Reviewer_6oF7", "summary": "This paper studies risk bounds for distributional regression and proposes a unified framework that accommodates both convex and non-convex constraints. Theoretical convergence guarantees are established, and the framework is shown to be broadly applicable. Extensive numerical experiments are conducted to validate its performance.", "questions": "See the questions in weakness. I am happy to increase the score if the authors can solve my question.", "limitations": "See the questions in weakness.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 2, "strengths_and_weaknesses": "**Strengths**: The paper is well-written and addresses a fundamental and interesting problem. It extends existing frameworks to the context of distributional regression, with broad applicability. The theoretical bounds derived appear sound, and the authors support their claims with extensive numerical results.\n\n**Weakness**: while the paper is well-written, I have some concerns regarding its novelty.\n\nThe authors should more clearly highlight the contributions of this work in comparison to existing literature. It would be helpful to include a structured comparison—such as a table or bullet-point list—detailing differences in problem formulations, theoretical results, and methodological innovations.\n\nI am particularly concerned about the originality of some key components. For instance, the projection technique used in Lemma 1 appears to be standard in this line of research. Moreover, Theorem 1 plays a central role in the paper but seems closely related to Theorem A.1 in [23]. Notably, a large portion of the proof of Theorem 1—approximately 80%—appears to be directly adapted from the proof of Theorem A.1 in [23]. This raises questions about the specific challenges addressed and the novel contributions of Theorem 1.", "was_revised": true, "final_justification": "Based on all the discussions, I choose to keep my positive score", "rebuttal": "Thank you very much for the comments and appreciation of our work. We provide clarifications and a detailed breakdown of the theoretical and methodological advances in our responses below.\n\n\n---\n\n**Q1: Contributions and structured comparison with prior work**\n\nThank you for the helpful suggestion. We have revised the manuscript to more clearly organize and expand our theoretical comparisons, including both existing content and several new additions. We have also created a consolidated summary (Rate, assumption, method) in Table 2, partially presented below.\n\n**Regarding the bounds**\n\nWe now explicitly compare the convergence rates obtained for our isotonic, trend filtering, and Dense ReLU network examples to those of corresponding methods in the CDF estimation literature, with a focus on both assumptions and statistical accuracy.\n\nAs an example, we reproduce below the updated discussion following Corollary 2 concerning our isotonic estimator:\n\n> The result in Corollary 2 establishes that distributional isotonic regression achieves an estimation rate of $n^{-2/3}$ for both the expected average CRPS and the worst-case MSE, as shown in (10). This result improves upon Theorem 3 in Hänzi et al. [28] in the univariate case, which only demonstrated convergence in probability for isotonic distributional regression. However, we emphasize that Hänzi et al. [28] study the more general setting of multivariate covariates, while our analysis is restricted to the univariate case ($d = 1$).\n\nDue to space constraints in this response, we briefly comment on the new comparisons we have added.\n\nOur isotonic result and Theorem 3.3 of Mösching and Dümbgen [44] are compared. Translating their setup into our notation, their quantity $F_{X_i}(t) := \\mathbb{P}(Y_i \\leq t \\mid X_i)$ corresponds to our $F_i^{\\ast }(t) := \\mathbb{P}(y_i \\leq t)$. The monotonicity assumption on $x \\mapsto F_x(t)$ aligns with our use of an isotonic constraint on the sequence $F_1^{\\ast }(t) \\leq \\cdots \\leq F_n^{\\ast }(t)$. However, our result does not assume Hölder continuity for the map $x \\mapsto F_x(t)$, or dense covariates in its domain. Despite these weaker requirements, we achieve a faster rate of $n^{-2/3}$ (up to logarithmic factors) for both average CRPS and worst-case MSE, compared to their $n^{-2\\alpha/(2\\alpha+1)}$ rate for the worst-case MSE, where $\\alpha \\in (0,1]$.\n\nWe compared our Dense ReLU estimator with DRF [10] and EnG [58], two methods explicitly designed for conditional distribution estimation.\n\nDRF estimates the conditional CDF using forest-based weights to form a weighted empirical distribution. EnG instead models $Y = g^{\\ast }(X, \\varepsilon)$, with $\\varepsilon \\sim \\text{Unif}[0,1]$, and trains a neural network via energy score minimization; the estimated CDF is obtained by sampling from the learned model. Translating their notation into our notation, both methods assume that $G^{\\ast }(x, t) := \\mathbb{P}(Y \\leq t \\mid X = x)$ is Lipschitz in $x$ for each fixed $t$, and that the covariate density is bounded away from zero and infinity.\n\nOur method instead assumes $G^{\\ast }(\\cdot, t) \\in \\mathcal{H}(l, \\mathcal{P})$, a broader class that includes the Lipschitz case ($\\mathcal{H}(1, \\\\{(1, d_0)\\\\})$) and mitigates the curse of dimensionality, which typically affects Lipschitz-based methods in high-dimensional settings. Moreover, we do not require assumptions on the covariate density. While DRF and EnG show convergence in probability at fixed $t$, our method yields high-probability bounds on the worst-case MSE uniformly over $t$, with explicit rate $\\frac{\\log n}{n} + \\phi_n \\log^4 n$.\n\n**Regarding the tightness of the bounds**\n\nIn the revision, we have incorporated a discussion regarding the tightness of the upper bounds. Specifically, the following has been added:\n\nTheorems 1–3 are general results that apply to a wide range of constraint sets, both convex and non-convex, and are not tailored to any specific estimation problem. For this reason, we do not provide lower bounds at the general level.\n\nHowever, in specific examples like isotonic regression, or total variation, the rates we obtain using these general theorems match known minimax rates for mean estimation with Gaussian noise. In this sense, the bounds we obtain are optimal for these concrete examples. \n\n**Regarding the novelty of the general Theorem 1, 2, and 3 in CDF literature**\n\nThe following has been added:\n\nTheorems 1–3 adapt classical techniques from regression, where projection-based estimators have been studied in both convex (e.g., [5, 14]) and non-convex (e.g., [46]) settings. These prior results focus on controlling mean squared error at a fixed point and were not developed for CDF estimation, which requires uniform control over all $t \\in \\mathbb{R}$. To our knowledge, our extensions provide the first such results for CDF estimation.\n\n---\n\n**Q2: Originality of the projection technique in Lemma 1**\n\nThank you for the question. We agree that Lemma 1 is technically straightforward. However, its purpose is not to introduce a novel projection technique per se, but rather to establish a conceptual and formal bridge between the original CRPS-based formulation (1) and the projection estimator (2) that underlies our methodological framework. This connection enables us to reinterpret the constrained distributional regression problem as a sequence of constrained projection problems across $t$, which is key to the development of our general theory. In particular, it provides the foundation for extending the estimator to time-varying constraint sets $K_t$, and for analyzing its performance under different structural assumptions, including those defined by neural networks, trend filtering, and monotonicity.\n\n---\n\n**Q3: Novelty of Theorem 1 relative to Theorem A.1 in [23]**\n\nThank you very much for pointing this out. Our result in Theorem 1 is conceptually related to Theorem A.1 of [23], but technically distinct. In their setting, the observations are Gaussian $Y = \\mu + Z$, with $Z \\sim \\mathcal{N}(0, \\alpha^2 I_n)$. The estimator $\\hat{\\mu}$ is the Euclidean projection of $Y$ onto a convex set $K$. Their goal is to bound the squared error $\\vert\\vert\\hat{\\mu} - \\mu\\vert\\vert_2^2 = t^{\\ast }$ in expectation, which is reduced to be controlled with high probability. For this, the \"KEY\" step is to show that the squared error $\\vert\\vert\\hat{\\mu} - \\mu\\vert\\vert_2^2 = t^{\\ast }$ concentrates around a deterministic benchmark $t_\\mu$, which characterizes the scale of estimation error. This is proved using a peeling argument that combines the strict concavity of an associated process with Gaussian concentration.\n\nOur setting also relies on convex projection, but differs in several important ways. We observe Bernoulli vectors $w(t) = 1_{\\{ y_i \\leq t \\}}$ (their $Y$), with sub-Gaussian, heteroscedastic noise $\\epsilon(t) = w(t) - F^{\\ast }(t)$ (their $Z$), and our target is the CRPS rather than squared error. Despite these differences, our analysis also hinges on addressing the same \"KEY\" step: showing that the squared $\\ell_2$ error $\\vert\\vert\\widehat{F}(t) - F^{\\ast }(t)\\vert\\vert_2^2$ (their $t^{\\ast }$) concentrates around a complexity-driven benchmark $\\eta$ (their $t_\\mu$).\n\nHowever, our proof takes a completely different route. We use (i) a sub-Gaussian concentration inequality for Lipschitz and separately convex functions to control a local empirical process, and (ii) a star-shaped self-normalization argument to extend this control to the full constraint set $K_t$. Because of the differences in observations, loss function, and techniques, our result is not a direct application of theirs.\n\n---\n***Table 2***\nBelow, we present Table 2 partially, showing only the summary of results for the Isotonic and Dense ReLU Network examples due to space constraints.\n\n### **Isotonic**\n\n| Work         | Method        | Assumptions                                                   | Convergence Rate                          |\n|--------------|---------------|----------------------------------------------------------------|-------------------------------------------|\n| Prior Work   | [27]          | $K_t$ monotone cone, $F^*(t)\\in K_t$                          | $o(1)$                                    |\n| Prior Work   | [44]          | $K_t$ monotone cone, $F^*(t)\\in K_t$, dense design, Hölder    | $n^{-2\\alpha/(2\\alpha+1)}, \\alpha\\le 1$ (Worst MSE)    |\n| This Work    | General rate  | $K_t$ monotone cone, $F^*(t)\\in K_t$                          | $n^{-2/3}$ (CRPS, Worst MSE)              |\n| This Work    | Fast rate     | $K_t$ monotone cone, $F^*(t)\\in K_t$, with few strict increases, $k(t)\\ll n$                             | $\\frac{1 + k(t)}{n} \\log \\left( \\frac{en}{1 + k(t)} \\right)$ (CRPS) |\n\n### **Dense ReLU Networks**\n\n| Work         | Method        | Assumptions                                                                                       | Convergence Rate                          |\n|--------------|---------------|----------------------------------------------------------------------------------------------------|-------------------------------------------|\n| Prior Work   | [36]          | Regression $y_i = f(x_i) + \\epsilon_i$, $f \\in \\mathcal{H}(l, \\mathcal{P})$                        | $\\phi_n$ (MSE)                            |\n| Prior Work   | [10]          | $F_i^{\\ast }(t) = G^{\\ast }(x_i, t) = \\mathbb{P}(y_i \\leq t \\mid x_i)$, $G^{\\ast }(\\cdot, t) \\in \\mathcal{H}(1, \\\\{(1,d_0)\\\\})$, Forest-Weighted empirical CDF | $o(1)$                                    |\n| Prior Work   | [57]          | As above but Neural Net-based sampling approach                                            | $o(1)$                                    |\n| This Work    | General rate  | $F_i^{\\ast}(t) = G^{\\ast}(x_i, t)$, $G^{\\ast}(\\cdot, t) \\in \\mathcal{H}(l, \\mathcal{P})$                          | $\\frac{\\log n}{n} + \\phi_n \\log^4 n$ (CRPS, Worst MSE) |"}, {"review_id": "XX67Vb6XMK", "reviewer": "Reviewer_WHfB", "summary": "The paper introduces a general procedure to learn a set of cdfs under constraints, and derive general estimation guarantees for this procedure. They then show several special cases of this approach and their theory, thereby strengthening existing results. Finally they close the paper by showing that their new approach outperforms state of the art methods in real and simulated data.", "questions": "- I dont fully understand in the first part whether (3) is new in this full generality of whether this has been done before?\n- There seems to be something wrong in the Probability statement of Theorem 3 (also in the proof in the appendix). I guess the > \\eta needs to shifted to the right?\n- In the proof of Theorem 3, how do you get from <= 2nB to <= 2^j \\eta^2?\n- It would be great to have a bit more details in Section 4.2. Now the K_t depend on a neural network that I guess needs to be optimized? How is this handled in practice? I am not that interested in the structure of the neural net, but more how the optimization in (3) works together with the optimization of the neural net.\n- Maybe a stupid question, but why is cdf learning with a regularization on TV distance called trend filtering?", "limitations": "Yes", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 4, "originality": 4, "strengths_and_weaknesses": "Strengths\n---------------\n- It appears to me the theoretical results are very strong\n- The paper is dense but well written. The proofs could benefit from a bit more explanation, but overall it is reasonable.\n- Though it is somewhat hidden in the appendix, the simulation and application section appears thorough as well\n\nWeaknesses\n----------------\nIf I had any critique of the paper than that it is almost too dense and appears to follow the trend in statistics to just throw as many results into a paper as possible. In particular, I am not sure what Appendix B is for. The results are great, but it appears to stand completely on it's own just making the paper longer. Similarly, while the main paper has a clear thread and the results have their place, I would wish for a bit more explanation in Section 4.2. Now K_t depends on a neural net as far as I understand; how exactly is this handled? (see also my questions below)", "was_revised": true, "final_justification": "I think this is a great paper with solid theory, generalizing previous work and providing general estimation guarantees. The authors addressed my questions adequately and also replied to my proof questions in details, which was interesting to read. I also see no unresolved issue and I am thus deciding to leave my score.", "rebuttal": "We thank the reviewer for the endorsement of our work and for the helpful suggestions to improve clarity and presentation. Below, we address the general concerns and specific questions in turn.\n\n---\n\n**Q1: Novelty of Equation (3)**\n\nThank you for the question. To the best of our knowledge, Equation (3), which formulates CDF estimation in distributional regression as a family of projection estimators indexed by $t$ and constrained to general sets $K_t \\subset \\mathbb{R}^n$, is novel in this level of generality. While projection-based methods have been used in other contexts, such as mean estimation under monotonicity or trend filtering constraints, their application to distributional regression problems has been limited. Our framework unifies these structural constraints within a general formulation for CDF estimation, allowing for a broad class of constraint sets including isotonic, total variation, and neural network-based sets, to be incorporated into a single projection-based estimator. This unified view enables new theoretical analysis and practical extensions for distributional regression.\n\n---\n\n**Q2: Probability statement in Theorem 3**\n\nThank you for pointing this out. We confirm that the inequality inside the probability statement in Theorem 3 is correct as written, though we agree it may initially appear unintuitive. The bound controls the $\\ell_2$ error $\\vert\\vert\\widehat{F}(t) - F^{\\ast}(t)\\vert\\vert_2$, and the probability is taken over the event that this exceeds a threshold $\\eta$ and an additive term involving $\\sqrt{n} \\cdot \\sup_t \\vert\\vert F^{\\ast}(t) - G(t)\\vert\\vert_\\infty$. In applications, the relevant quantity is the average per-sample error, which is normalized by $\\sqrt{n}$. Therefore, after normalizing by $\\sqrt{n}$, the final bound depends on both the approximation error, $\\sup_t \\vert\\vert F^{\\ast}(t) - G(t)\\vert\\vert_\\infty$, and $\\eta/\\sqrt{n}$. In consequence, it suffices to choose $\\eta$ of order $\\sqrt{\\log n} + \\sqrt{n} \\cdot \\sup_t \\vert\\vert F^{\\ast}(t) - G(t)\\vert\\vert_\\infty$.\n\nThis behavior is confirmed in Corollary 4 for our Dense ReLU network example. There, we assume that $G^{\\ast}(\\cdot, t) \\in \\mathcal{H}(l, \\mathcal{P})$ where $F_i^{\\ast}(t) = G^{\\ast}(x_i, t) $, see Assumption 1 in Appendix C, and define the constraint set\n\n\n$$\nK_t = \\\\{ \\theta \\in \\mathbb{R}^n : \\theta_i = f(x_i), \\text{ for some } f \\in \\mathcal{F}(L, \\nu) \\\\}\n$$\n\n\nwhere $\\mathcal{F}(L, \\nu)$ is a class of ReLU networks with appropriate choice for the number of layers and neurons. In this case, the approximation error scales as $\\phi_n = \\max_{(p, M) \\in \\mathcal{P}} n^{-p/(2p + M)}$, and the final rate becomes\n\n$$\n\\sup_{t \\in \\mathbb{R}} \\frac{1}{n} \\sum_{i=1}^n\\left(\\widehat{F}_i(t)-F_i^*(t)\\right)^2 = O_P\\left(\\frac{\\log n}{n} + \\phi_n \\log^4 n \\right).\n$$\n\nThus, the resulting rate is meaningful and matches known convergence rates for mean estimation under Gaussian noise using Dense ReLU networks. This shows that our bound is both general and sharp in concrete, practically relevant settings.\n\n---\n\n**Q3: Step from $\\leq 2nB$ to $\\leq 2^j \\eta^2$ in the proof of Theorem 3**\n\nThank you for this question. The step involving the transition from the upper bound $\\sup_t \\vert\\vert\\widehat{F}(t) - G(t)\\vert\\vert^2 \\leq 2nB$ to the summation over intervals of the form $(2^{j-1}\\eta^2, 2^j \\eta^2]$ corresponds to a peeling argument.\n\nWe briefly clarify the logic. First, we restrict attention to the event\n\n$$\n \\\\{\\sup_{t \\in \\mathbb{R}} \\vert\\vert\\widehat{F}(t) - G(t)\\vert\\vert^2 > \\eta^2 \\\\}\n\\cap\n \\\\{\\sup_{t \\in \\mathbb{R}} \\vert\\vert\\widehat{F}(t) - G(t)\\vert\\vert^2 \\leq 2nB\\\\},\n$$\n\nwhich allows us to consider only deviations above $\\eta^2$ and below the known upper bound $2nB$.\n\nWe then partition the interval $[\\eta^2, 2nB]$ into dyadic subintervals $(2^{j-1} \\eta^2, 2^j \\eta^2]$ for $j = 1, \\ldots, J$, where $J$ is the smallest integer such that $2^J \\eta^2 \\geq 2nB$. This allows us to write:\n\n$$\n\\mathbb{P}\\left( \\sup_{t \\in \\mathbb{R}} \\vert\\vert\\widehat{F}(t) - G(t)\\vert\\vert^2 > \\eta^2, \n\\sup_{t \\in \\mathbb{R}} \\vert\\vert\\widehat{F}(t) - G(t)\\vert\\vert^2 \\leq 2nB \\right) \n\\leq \n\\sum_{j=1}^J \\mathbb{P}\\left( \\sup_{t \\in \\mathbb{R}} \\vert\\vert\\widehat{F}(t) - G(t)\\vert\\vert^2 \n\\in (2^{j-1}\\eta^2, 2^j \\eta^2] \\right).\n$$\n\nEach term in the sum corresponds to the event that the supremum falls within a specific dyadic shell.\n\nWe have added a brief clarification in the proof to make this peeling step more transparent to the reader.\n\n---\n\n**Q4: Practical implementation and optimization of neural network-based $K_t$**\n\nThank you for the question. In practice, we fix the network architecture in advance by selecting the number of layers and neurons ($L$ and $\\nu$). For each threshold value $t$, we then solve the projection problem in Equation (3) by training a neural network that maps covariates $x_i$ to estimates $F_i(t)$, subject to the specified architecture.\n\nThis training procedure is repeated independently for each $t$, which makes it highly parallelizable and scalable in practice. While the structure of the network remains fixed across $t$, the parameters are learned separately for each projection, effectively allowing us to approximate the conditional CDF at each threshold using a specialized network.\n\n---\n\n**Q5: Terminology — Why “trend filtering” is used for CDF learning with total variation regularization**\n\nThank you for the question. The reason we refer to this estimator as “trend filtering” is because the terminology is standard in the literature, see [41, 61, 23], when the constraint set $K_t$ is defined via a bound on total variation of order $r$. We will add a clarifying sentence in the revised manuscript to make this connection explicit.\n\n---\n\n***References***\n\n[23] Adityanand Guntuboyina, Donovan Lieu, Sabyasachi Chatterjee, and Bodhisattva Sen. Adaptive risk bounds in univariate total variation denoising and trend filtering. The Annals of Statistics, 48:205–229, 2020.\n\n[41] Enno Mammen and Sara Van De Geer. Locally adaptive regression splines. The Annals of Statistics, 25(1):387–413, 1997.\n\n[61] Ryan J Tibshirani. Adaptive piecewise polynomial estimation via trend filtering. 2014"}]

Author rebuttal:

Meta-review: The paper develops risk bounds for distributional regression (estimating conditional CDFs) under structural constraints. It gives general upper bounds on CRPS and worst-case MSE under convex constraints (with applications to isotonic and trend filtering) and also provides a general upper bound for certain non-convex constraints, illustrated with neural network estimators. Experiments support the theory.

Reviewers highlighted that it is
* clearly written and well-organised;
* presents significant results addressing an important problem; and
* is theoretically sound.

Reviewers asked for clarifications on the novelty relative to prior work and on the scaling and tightness of the bounds. The authors’ detailed responses addressed these points satisfactorily. I concur with the reviewers and recommend acceptance.
Decision: Accept (poster)