Distributional regression: CRPS-error bounds for model fitting, model selection and convex aggregation

Authors: clément  dombry
ahmed  zaoui
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-13 10:53:52

Abstract: Distributional regression aims at estimating the conditional distribution of a target variable given explanatory co-variates. It is a crucial tool for forecasting when a precise uncertainty quantification is required. A popular methodology consists in fitting a parametric model via empirical risk minimization where the risk is measured by the Continuous Rank Probability Score (CRPS). For independent and identically distributed observations, we provide a concentration result for the estimation error and an upper bound for its expectation. Furthermore, we consider model selection performed by minimization of the validation error and provide a concentration bound for the regret. A similar result is proved for convex aggregation of models. Finally, we show that our results may be applied to various models such as EMOS, distributional regression networks, distributional nearest neighbours or distributional random forests and we illustrate our findings on two data sets (QSAR aquatic toxicity and Airfoil self-noise).

Main paper: 
 introduction: 
 Motivation and related literature. We consider in this paper the distributional regression problem, where we want to estimate the conditional distribution of a target random variable Y given explanatory co-variates X. We assume (X, Y ) ∈ R d × R and we letEQUATIONbe the conditional cumulative distribution functions (c.d.f.). Estimation is built on a training sample D n = {(X i , Y i ), 1 ≤ i ≤ n} of independent identically distributed (i.i.d.) copies of (X, Y ).Let us emphasize that distributional regression is much more challenging than standard regression where only a point prediction for Y given X = x is provided which typically reduces to an estimation of the conditional expectation. The conditional distribution provides a full account for the variability of Y given X = x and distributional regression is therefore a crucial tool for forecasting when a precise uncertainty quantification is required.Distributional regression is of relevance in various applied fields. To name a few, let us mention statistical post-processing of weather forecast (Matheson and Winkler, 1976; Gneiting et al., 2005) , forecasting of wind gusts (Baran and Lerch, 2015) , of solar irradiance (Schulz et al., 2021) , of ICU length stays during COVID-19 pandemy (Henzi et al., 2021a) , of breast cancer ODX score in oncology (Al Masry et al., 2024) . . .The methodology relies on probabilistic forecast where the forecaster produces a predictive distribution for the quantity of interest (Gneiting and Katzfuss, 2014) . Fitting a distributional regression model typically involves the minimization of a proper scoring rule, the most popular one being the Continuous Ranked Probability Score (CRPS, Matheson and Winkler, 1976; Gneiting and Raftery, 2007) . It compares the actual observation with the predictive distribution in a comprehensive manner. Many different models have been proposed for distributional regression among which: Analog similar to a nearest neighbor method (KNN, Toth, 1989) , Ensemble Model Output Statistics similar to Gaussian heteroscedastic regression (EMOS, Gneiting and Raftery, 2007) , Isotonic Distributional Regression (Henzi et al., 2021b) and more recent machine learning methods such as Distributional Regression Network (DRN, Rasp and Lerch, 2018) or Distributional Random Forest (DRF, Ćevid et al., 2022) .Despite these successful achievements in terms of methods and applications, a sound theory for distributional regression via scoring rule minimization is still missing. Recently, minimax rates of convergence with CRPS-error have been considered for distributional regression (Pic et al., 2023) . The aim of this paper is to provide statistical learning guarantees for model fitting, model selection and convex aggregation based on CRPS minimization.Model selection and convex aggregation are very important techniques in the field of statistics and machine learning and have been used very successfully for regression function estimation, see Tsybakov (2003) ; Bunea et al. (2007) . The methods use two independent samples; the first sample, called training sample, is used to construct the initial estimators which are constituted as a dictionary (a collection of candidates). The second sample, called the validation sample, is used to aggregate them. Model selection enables the selection of the best candidate in the dictionary, while convex aggregation provides the optimal convex combination from these candidates.Contributions. Our main results include a concentration bound for the theoretical risk when a parametric model is fitted via CRPS empirical risk minimization. We also consider model selection and model aggregation via CRPS minimization on a validation set and provide concentration bounds for the regret. Our results are first derived under sub-Gaussianity assumptions and then extended under weaker moment assumptions. We show that they apply to several popular models such as EMOS, DRN, KNN or DRF and provide a short illustration on two different datasets.Structure of the paper. Section 2 first provides some background on distributional regression and proper scoring rules and then presents precisely the main methods and goals. The main results are stated in Section 3: an oracle inequality for the estimation error in model fitting (Theorem 1), and concentration bounds for the regret in model selection (Theorem 2) and model aggregation (Theorem 3). In Section 4, some specific models are introduced and the assumptions for our results to hold are checked. A short illustration on two different data set is provided in Section 5. Finally, an appendix contains all the proofs as well as some additional results.2 Background on distributional regression and main goals 2.1 Probabilistic forecast and its evaluation with scoring rulesWe first consider the simple setting of probabilistic forecast without co-variate where a future observation Y is predicted by a probability distribution F , called predictive distribution. Proper scoring rules are used in order to compare the predictive distribution F and the materializing observation y which are objects of different nature. Let P 0 ⊂ P(R) denote a subset of the set of all probability measures on R, often identified with their c.d.f. A scoring rule2 on P 0 is a function S : P 0 × R → [0, +∞). The quantity S(F, y) is interpreted as the error between the predictive distribution F and the materializing observation y. The mean error when Y has "true" distribution G is denoted byS(F, G) = E Y ∼G [S(F, Y )].The following notion of proper and strictly proper scoring rule is central in the theory. Definition 1. The scoring rule S is said proper on P 0 when S(F, G) ≥ S(G, G), for all F, G ∈ P 0 .(2) It is said strictly proper when equality in Eq. (2) implies F = G.Stated differently, the scoring rule S is strictly proper on P 0 when arg minF ∈P0 S(F, G) = {G}, for all G ∈ P 0 .The interpretation is that, in order to minimize its mean error, the forecaster has to predict the "true" observation distribution G.In this paper, we consider the Continuous Ranked Probability Score (CRPS, Matheson and Winkler 1976) . This scoring rule is defined by the formulaS(F, y) = R 1 {y≤z} -F (z) 2 dz (3)for all F finite absolute moment. Here the subset P 0 ⊂ P(R) is the Wasserstein spaceP 1 (R) = F ∈ P(R) : m 1 (F ) = R |y|F (dy) < ∞of probability measures on R with finite first moment. One can easily check from this definition thatEQUATION)This quantity is nonnegative and vanishes if and only if F = G, ensuring that the CRPS is a strictly proper scoring rule. For discrete predictive distributions of the form F = n i=1 w i δ yi (with δ y the Dirac mass at y), the CRPS can be computed simply (Gneiting and Raftery, 2007) byS(F, y) = n i=1 w i |y i -y| - 1 2 i̸ =j w i w j |y i -y j |.Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?Answer: [Yes] Justification: The main claims presented in the abstract and introduction accurately reflect the contributions and scope of the paper. The abstract outlines the study's focus on distributional regression and its importance for uncertainty quantification in forecasting. It mentions key methodologies such as fitting parametric models using empirical risk minimization with CRPS. The claims about providing concentration results for estimation error, bounds on expectation, model selection, and convex aggregation via validation error minimization, as well as application to various regression models, are all consistent with the paper's contributions and scope as described. Additionally, the examples cited (QSAR aquatic toxicity and Airfoil self-noise datasets) demonstrate practical applications of the discussed methodologies. Therefore, the abstract and introduction effectively summarize the paper's objectives, findings, and relevance.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 be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 
model fitting, model selection and convex aggregation: 
 We briefly present the methods and objectives that we address in this paper.In a regression framework, we observe a sampleD n = {(X i , Y i ), 1 ≤ i ≤ n} of independent copies of (X, Y ) ∈ R d × R.Distributional regression aims at estimating the conditional distribution Y |X = x characterized by its c.d.f. F * x defined in Eq. ( 1). The marginal distribution of X is denoted by P X . The forecaster uses the training sample D n and some algorithm to build a functional estimator Fn : x → Fn,x of the map F * : x → F *x . The accuracy of this estimator is here measured by its theoretical riskR( Fn ) = E S( Fn,X , Y ) .where expectation is taken with respect to the joint law of (X, Y ). This quantity can be seen as the counterpart of the mean squared error in point regression. The excess risk of Fn is defined asR( Fn ) -R(F * ) = E S( Fn,X , Y ) -S(F * X , Y ) = E S( Fn,X , F * X ) -S(F * X , F * X ) ≥ 0.The nonnegativity is ensured by the fact that S is a proper scoring rule according to Definition 1. If S is strictly proper, the excess risk is equal to 0 if and only if Fn,x = F * x almost everywhere (with respect to P X ). For the CRPS, Eq. ( 4) implies that the excess risk can be rewritten asR( Fn ) -R(F * ) = E R Fn,X (u) -F * X (u) 2 du .Our first interest lies in model fitting by empirical risk minimization. Here we consider a parametric family(F θ ) θ∈Θ , Θ ⊂ R K , where F θ : x ∈ R d → F θ,x ∈ P 0 ⊂ P(R).The empirical risk associated with F θ is computed on the training sample D n byRn (F θ ) = 1 n n i=1 S(F θ,Xi , Y i )and is an empirical counterpart of the theoretical risk R(F θ ). Empirical risk minimization consists in finding θn = arg min θ∈Θ Rn (F θ )and proposing the estimator F θn which is thought as almost optimal within the family (F θ ) θ∈Θ . A classical decomposition of the excess risk of the corresponding estimator is given byR(F θn ) -R(F * ) = R(F θn ) -inf θ∈Θ R(F θ ) + inf θ∈Θ R(F θ ) -R(F * ) .where the two terms are called the estimation error and the approximation error respectively. The approximation error is deterministic and depends on the ability of the family (F θ ) θ∈Θ to approximate F * . The estimation error is random as it depends on the training sample D n . Our first goal is the following: Goal 1: provide non asymptotic estimates for the estimation error R(F θn ) -inf θ∈Θ R(F θ ).Our second interest lies in model selection and convex aggregation via validation error minimization.Here we suppose that a validation sampleD ′ N = {(X ′ i , Y ′ i ), 1 ≤ i ≤ N } is available, which is assumed independent of the training sample D n .Model selection. A common situation in machine learning is that we have M algorithms at hand that are trained on D n , resulting in models F 1 n , . . . , F M n . In order to select the best model, we compute the empirical risks on the validation sampleR′ N ( F m n ) = 1 N N i=1 S( F m n,X ′ i , Y ′ i )and select the modelm = arg min 1≤m≤M R′ N ( F m n ).An oracle having access to the theoretical risk would have selectedm * = arg min 1≤m≤M R( F m n ),leading to the definition of the regretR( F m n ) -R( F m * n ) = R( F m n ) -min 1≤m≤M R( F m n ).Goal 2: provide non asymptotic estimates for the regret R(F m n ) -min 1≤m≤M R( F m n ).Convex aggregation. We define the convex aggregation of models F 1 n , . . . , F M n with weights λ 1 , . . . , λ M byF λ n,x = M m=1 λ m F m n,x .Here λ = (λ 1 , . . . , λ M ) is an element of the simplex Λ M = {λ : λ m ≥ 0, 1≤m≤M λ m = 1}. The best weights for convex aggregation are obtained by minimization of the validation error, i.e. λ = arg minλ∈Λ M R′ N ( F λ n ).Goal 3: provide non asymptotic estimates for the regret R(F λ n ) -inf λ∈Λ M R( F λ n ).Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes] Justification: The approach is very general and may apply to parametric and non-parametric methods and does not assume any model assumption. Observations are assumed independent and identically distributed which is a very standard assumption in statistical learning. The results are derived under a sub-Gaussian assumption and a weaker moment condition is also considered.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a 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 on 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 a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not 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 to address problems of privacy and fairness. • 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 acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
main results: 
 We present our main results for model fitting, model selection and convex aggregation in distributional regression. We first focus on concentration bound under sub-Gaussianity assumptions. In a second step, we extend our results under weaker moment assumptions. The definition of sub-Gaussian random variables and distributions is given in Appendix A as well as some useful concentration inequalities. All the proofs are postponed to Appendices B, C and D.We provide concentration results for the estimation error when fitting a parametric model via empirical CRPS-error minimization according to the framework described in Section 2.2.2. Our working assumptions are the following. Assumption 1. (sub-gaussianity)i) the variable Y is β 1 -sub-Gaussian; ii) there exists β 2 > 0 such that m 1 (F θ,X ) is β 2 -sub-Gaussian for all θ ∈ Θ.Assumption 1-i) is classical in a regression setting (Györfi et al., 2002; Biau and Devroye, 2015) while Assumption 1-ii) characterizes the sub-Gaussian behavior of the absolute moment of F θ,X . Importantly, Assumption 1 implies that the variableZ θ = S(F θ,X ) is 2(β 2 1 + β 22 )-sub-Gaussian for all θ ∈ Θ, see Proposition 5 in the appendix.Our second assumption requires compactness of the parameter space and Lipschitz continuity of the model (F θ ) θ∈Θ . We denote by W 1 the Wasserstein distance of order 1 on the space P 1 (R). Assumption 2. (regularity) The parameter space Θ ⊂ R K is compact and there exists a constantL > 0 such that W 1 (F θ1,x , F θ2,x ) ≤ L∥θ 1 -θ 2 ∥ for all θ 1 , θ 2 ∈ Θ, x ∈ R d .This assumption ensures the Lipschitz continuity of the empirical risk θ → R n (F θ ), see Proposition 6 and, by compactness, the existence of θn , the ERM estimator defined in Eq. ( 5). We provide in Section 4 examples of popular models verifying the two assumptions.Our main result provides a concentration bound on the estimation error. We let R > 0 be such that Θ is included in the ball centered at 0 and with radius R. Theorem 1. Under Assumptions 1-2, for all δ ∈ (0, 1), the estimation error satisfies, with probability at least 1 -δ,EQUATIONwith c β = 64(β 2 1 + β 2 2 ) and provided that n is large enough so that n log(2n K /δ) ≥ (48LR) 2 /c β .The inequality (6) is commonly known as an oracle inequality. The convergence rate depends on log(2n K /δ)/n with K the dimension of the parameter space and n the sample size. A key point in the proof is the the combinatorial complexity of Θ in terms of ϵ-net, see Devroye et al. (1996) . A bound in expectation can easily be deduced from Theorem 1 and its proof.Corollary 1. Under Assumptions 1-2, E R(F θn ) -inf θ∈Θ R(F θ ) ≤ 2 c β log(2n K ) n provided n is large enough so that n log(2n K ) ≥ (48LR) 2 /c β .In particular, the ERM approach is weakly consistent when n goes to infinity.We provide concentration results for the regret in model selection and convex aggregation. Recall that the methodology is based on minimization of the test error on a validation sample, see the precise framework in Section 2.2.3. We need the following assumption. Assumption 3. (sub-Gaussianity)i) the variable Y is β 1 -sub-Gaussian;ii) conditionally on D n , there existsβ n = β(D n ) > 0 such that m 1 ( F m n,X) is β n -sub-Gaussian for all m = 1, . . . , M .We will see in Section 4 that for several popular models such as distributional nearest neighbors or distributional random forest, Assumption 3-ii) is satisfied with Vershynin, 2018, Exercise 2.5.8 p.25) .β n = max 1≤i≤n |Y i | which is of order β 1 √ log n (Under this assumption, a control of the regret in model selection is provided by the following theorem.Our results hold conditionally on the training set D n . Theorem 2. Under Assumption 3, for all δ ∈ (0, 1), the regret in model selection satisfiesP R( F m n ) -min 1≤m≤M R( F m n ) ≤ 4 c n log(2M/δ)/N D n ≥ 1 -δ with c n = β 2 1 + β 2 n . Furthermore, E R( F m n ) -min 1≤m≤M R( F m n ) D n ≤ 8 c n log(2M ) N .When selecting the hyperparameter k in nearest neighbor distributional regression or mtry in distributional random forest, one has respectively M = n and M = d if all possible values are considered -see Section 4 for more details. Note that when the response variable Y is bounded, then c n does not depend on n and is constant. We now state a bound for the regret in convex aggregation. Theorem 3. Under Assumption 3, for all δ ∈ (0, 1), the regret in convex aggregation satisfiesP R( F λ n ) -inf λ∈Λ R( F λ n ) ≤ 8 c n log(2N M /δ)/N D n ≥ 1 -δ, with c n = β 2 1 + β 2 n provided N is large enough so that N log(2N M /δ) ≥ 48 2 /c n . Furthermore, E R( F λ n ) -inf λ∈Λ R( F λ n ) D n ≤ 2 c n log(2N M ) N provided N log(2N M ) ≥ 48 2 /c nThe preceding results hold under a strong sub-Gaussianity. We next adapt our results to the following weaker moment condition. Assumption 4. There is p ≥ 2 andD > 0 such that E[|Y | p ] ≤ D and E[|m 1 (F θ,X )| p ] ≤ D for all θ ∈ ΘWe recall that, for p ≥ 1, the L p -norm of a random variable Z is defined by ∥Z∥ L p = E[|Z p |] 1/p . Theorem 4. Under Assumptions 2 and 4, we haveR(F θn ) -inf θ∈Θ R(F θ ) L p ≤ Cn -p/(2(p+K)) .with constant C > 0 depending only on K, p, L, D, R and possibly made explicit from the proof. This implies the bound for the expected estimation errorE R(F θn ) -inf θ∈Θ R(F θ ) ≤ Cn -p/(2(p+K)) .For p large, the rate of convergence tends to the parametric rate n -1/2 obtained (up to a logarithmic factor) in the sub-Gaussian case . We propose additional results for model selection (Theorem 5) and convex aggregation (Theorem 6) which are, for the sake of brevity, postponed to Appendix D.2.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Answer: [Yes] Justification: the different assumptions are clearly stated in the suitable Assumption environment. At the beginning of each Theorem, the relevant assumptions are recalled. Furthermore, the assumptions are discussed and checked on a variety of different examples. All the proofs are provided and carefully explained in the supplementary material.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 crossreferenced. • All assumptions should be clearly stated or referenced in the statement of any theorems.• The proofs 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 to 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 be properly referenced. 
examples and popular models for distributional regression: 
 We present the most popular models for distributional regression for which we want to apply our results. The first two (EMOS and DRN) are parametric, while the last two (distributional k-NN and DRF) are fully non-parametric.EMOS. The EMOS model was designed by Gneiting et al. (2005) for the purpose of statistical post-processing of ensemble weather forecast. In this framework, the predictive distribution takes the form of a discrete distributions m -1 m l=1 δ y l with members y 1 , . . . , y m corresponding to different scenarios obtained from numerical weather predictions. Such forecast typically suffer from bias and underdispersion so that statistical post-processing is needed. The explanatory variable for distributional regression are Ensemble Member Output Statistics such as the ensemble mean ȳ and ensemble variance v 2 y . In its simplest version, EMOS models the predictive distribution as a Gaussian distribution with parameters m = β 0 + β 1 ȳ andσ 2 = β ′ 0 + β ′ 1 v 2 y . This is a parametric model with θ = (β 0 , β 1 , β ′ 0 , β ′ 1 ) in Θ = R 2 × (0, ∞) 2 .Minimum CRPS estimation is used for model fitting as described in 2.2.2. This simple yet successful method has encountered many generalizations (Scheuerer 2013; Scheuerer and Hamill 2015; Baran and Nemoda 2016) and we shall consider the following general setting. Given x ∈ R d , the predictive distribution takes the formF (•; θ) = F ((• -m)/σ) with θ = (α, β, α ′ , β ′ ) ∈ R 1+d × R 1+d and m(x; θ) = α + β ⊤ x σ 2 (x; θ) = softplus(α ′ + β ′ ⊤ x)with softplus(u) = log(1 + e u ). The predictive distribution belongs to a location/scale family and the location and log-scale parameters are linear in the covariate x. The softplus link function ensures positivity of the scale.Distributional Regression Networks (DRN). In order to consider higher dimension co-variates together with non-linear dependence, Rasp and Lerch (2018) introduce distributional regression networks. In a setting similar to EMOS, neural networks are used to model complex response functions m(x; θ) and log σ 2 (x; θ). In the case of a single hidden layer with H units, the equations writeEQUATIONwhere g denotes the activation function acting componentwise, (α, β, α ′ , β ′ ) ∈ R 1+H × R 1+H the parameters for the output layers and (γ, δ) ∈ R H+Hd the biases and weights of the hidden layer.The global parameter θ = α, β, α ′ , β ′ , γ, δ) lies in dimension (d + 3)H + 2. Note that in absence of hidden layers, DRN reduces to EMOS. Extension to a MLP structure with multiple hidden layers is straightforward, see Schulz and Lerch (2022) for a review of DRNs and their applications.We next provide conditions ensuring that our results hold for the EMOS and DRN models.Proposition 1. Let (F θ ) θ∈Θ be a EMOS or DRN model with parameters restricted to a compact subset Θ. If Y is sub-Gaussian, X is bounded and the activation function g is Lipschitz continuous, then Assumptions 1, 2 and 3 are satisfied.Distributional k-Nearest Neighbors (KNN). The predictive distribution is built on a straightforward extension of k-nearest neighbor regression: we setFn,x (y) = 1 k k i=1 1 {Xi∈knn(x)} 1 {Yi≤y}where knn(x) denotes the set of k nearest neighbors of x in the training set (X i ) 1≤i≤n . The main hyperparameter is the number k of neighbors, leading to a model selection problem as considered in 2.2.3. This simple method is known as the Analog method in the framework of statistical postprocessing of weather forecast (Toth, 1989) .Distributional Random Forests (DRF). Random forests for distributional regression have been introduced by Ćevid et al. (2022) and are a powerful nonparametric method. We describe only the main lines of the method. Recall that in standard regression, Breiman's Random Forest (Breiman, 2001) estimates the regression function byμ(x) = 1 B B b=1 T (b) (x) = 1 B B b=1 1 |L b (x)| n i=1 Y i 1 {Xi∈L b (x)} ,where T 1 , . . . , T B denote randomized regression trees built on bootstrap samples of the original data and L b (x) the leaf containing x in the tree T b . Interverting the two sums yieldsμ(x) = n i=1 1 B B b=1 1 {Xi∈L b (x)} |L b (x)| Y i = n i=1 w ni (x)Y i .Using these Random Forest weights, the predictive distribution writesFn,x (y) = n i=1 w ni (x)1 {Yi≤y} .This strategy based on Breiman's regression tree is used by Meinshausen (2006) to construct the quantile regression forest. Different splitting strategies in the tree construction have been considered to better detect changes in distribution rather than changes in mean (Taillardat et al. 2016; Athey et al. 2019) . Note that Ćevid et al. ( 2022) minimizes a scoring rule to construct the splits of the trees. The main hyperparameter of this nonparametric method is the so-called mtry parameter that controls the number of co-variates tested at each split in the trees.For the distributional KNN and DRF models, model fitting with ERM is irrelevant and we only consider model selection and convex aggregation.Proposition 2. Let Fn be a KNN or DRF model fitted on a training set D n . Then Assumption 3-ii) is satisfied withβ n = max 1≤i≤n |Y i |.Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?Answer: [Yes] Justification: All the details of the numerical experiment can be found in the Models and Methodology parts of Section 5. 
numerical results: 
 In this section, we illustrate how model selection and model aggregation work on real data sets. The cross validation methodology, widely used in practice, is justified by the theoretical results obtained in Section 3. The source code for these experiments can be found at https://github. com/ZaouiAmed/Neurips2024_DistributionalRegression.Datasets. We consider two datasets used in the framework of heteroscedastic regression with reject option as detailed in Zaoui et al. (2020) . The first dataset, QSAR aquatic toxicity (Ballabio et al., 2019) is referred to as qsar and comprises 546 observations with 8 numerical features used to predict acute toxicity in Pimephales promelas. The toxicity output ranges from 0.12 to 10.05 with evidence of a low heteroscedasticity. The second dataset, Airfoil Self-Noise (Brooks et al., 2014) is referred to as airfoil and consists of 1503 observations with 5 measured features from aerodynamic and acoustic tests. The output represents the scaled sound pressure level in decibels, ranging from 103 to 140 with evidence of a strong heteroscedasticity..We consider the KNN and DRF models to predict the conditional distribution of the output variable. Our focus lies on hyperparameters selection via minimization of the validation error, where the main hyperparameter is the number k of neighbors and the number mtry of variables considered at each split for KNN and DRF respectively. We utilize the implementation of these methods from the R packages KernelKnn and DRF . For DRF, a preliminary exploration shows that sensible choices for the other parameters are num.trees= 1000, sample.fraction=0.9, min.node.size= 1 and default values for others parameters. Finally, the R package ScoringRule is used for CRPS computation and the optim function based on the Nelder-Mead method is used for parameter optimization in convex aggregation.Methodology. We use the same methodology for the two datasets. We divide the data into three parts: 50% for training, 20% for validation and 30% for testing. In a first stage, the training set is used to train the model (KNN or DRF) for the various hyperparameters (k or mtry) and the validation set is used to select best hyperparameters k or mtry. Furthermore, the validation set also used to choose for model selection (MS), choosing between KNN and RF, and the best convex aggregation (CA). In a second stage, the different models (KNN, RF, MS and CA) are refitted on the union of training and validation sets (70% of data) and evaluated on the test set (30%) by computing the test CRPS-error. This process is repeated 100 times for different random splits of the data into training, validation and testing sets. The distributions of the test error for KNN, DRF, MS and CA is then analyzed in terms of mean, standard error and boxplot.Results. For the sake of brevity, only the results for the qsar dataset are presented, while the results for the airfoil dataset are postponed to Appendix F. The first stage of the procedure where the hyperparameter (k or mtry) is selected by minimization of the validation error is shown on the left plot of Fig. 1 for KNN and on the middle plot for DRF. In both cases, the validation error curve shows a clear minimum allowing to select the hyperparameter k = 8 and mtry = 4 corresponding to CRPS-error of 0.696 and 0.678 on the validation set respectively. In the second stage, KNN, DRF, MS and CA are evaluated on the test set and the right plot of Fig. 1 shows the distribution of the test error over 100 repetitions. We can see from the boxplots that the distribution of the test error is slightly larger for KNN than for DRF, that MS achieves almost the same performance as RF and that CA achieves slightly better performance than DRF. Hence model selection and convex aggregation accomplish their goal. The numerical values of the means together with standard errors are summarized in Table 1 , confirming our analysis from the boxplots.  
conclusion: 
 In this work, we considered distributional regression with error assessed using the CRPS and investigated model fitting via empirical risk minimization. Additionally, we explored classical aggregation procedures: model selection and convex aggregation where two independent samples are used, the first for model construction, the second for model selection or convex aggregation. We derived oracle concentration inequalities for the estimation error and established an upper bound for its expectation within the sub-Gaussian framework and beyond, considering weaker moment assumptions. These new theoretical results are solid mathematical justifications for common practices in the framework of distributional regression. In future work, we will study the minimax convergence rates for our approaches and consider the use of empirical process theory to strengthen our results. 
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: standard error are provided in parenthesis in Table 1 and reveal the statistical significance of the results (e.g. the convex aggregation is significantly better than the other models).Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence 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.) • The assumptions made should be given (e.g., Normally distributed errors).• It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 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. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [NA] Justification: The main purpose of the paper is theoretical and the numerical experiment is very quick to run so that it is not useful to discuss the computer resource.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, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs 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 it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines?Answer: [Yes] Justification: Nothing specific to declare.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 circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration 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: [Yes] Justification: The introduction remind the reader that uncertainty quantification is important for decision making. 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 or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness 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 a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an 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 a 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 as 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 or 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 to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe 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)? Answer: [NA] Justification: [NA] Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a 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 a best faith effort. 12. Licenses for existing assets Question: 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: 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 an equivalent approval/review based on the requirements of your country or institution) were obtained? Answer: [NA] Justification: [NA] 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 NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: The concept of sub-Gaussianity has been extensively studied and is a key tool in deriving concentration inequalities, as documented in the monographs by Vershynin (2018) and Boucheron et al. (2013) . For the convenience of the reader, we recall the definition and basic properties of sub-Gaussian random variables as well as the useful Hoeffding inequality.A.1 Sub-Gaussian random variables Definition 2 (Sub-Gaussian random variable). Let β > 0. The random variable X is said to beβ-sub-Gaussian if E e λ(X-E[X]) ≤ e 1 2 λ 2 β 2 for all λ ∈ R. Proposition 3. Let β, β 1 , β 2 > 0. i) If X is β-sub-Gaussian, then |X| is also β-sub-Gaussian. ii) If X is a bounded random variable such that X ∈ [a, b] almost surely, then X is β-sub- Gaussian with β = (b -a)/2. iii) If X, Y are random variables such that X is β-sub-Gaussian and 0 ≤ Y ≤ X, then Y is also β-sub-Gaussian. iv) If X i is β i -sub-Gaussian for i = 1, 2, the sum Y = X 1 + X 2 is β-sub-Gaussian with β = 2(β 2 1 + β 2 2 ).In this section we gather several technical results which are used to derive the contributions of this work. We start with the general Hoeffding's inequality (Vershynin, 2018 , Theorem 2.6.2 p. 28).Proposition 4 ( Hoeffding (1963) ). Let Z 1 , . . . , Z n be independent random variables such that Z i is β-sub-Gaussian for all i = 1 . . . , n. Then, for all t ≥ 0 we haveP 1 n n i=1 Z i -E[Z i ] ≥ t ≤ 2 exp - nt 2 2β 2 .The following Lemma is often useful to derive upper bound in expectation from concentration inequalities.Lemma 1. Let a ≥ 1 , b > 0 and Z a positive random variable such thatP (Z ≥ t) ≤ exp(a -bt 2 ) for all t ≥ 0.Then we haveE[Z] ≤ 1 + √ π 2 a b ≤ 2 a b .Proof.The assumption implies that P (Z ≥ t) ≤ min 1, exp(a -bt 2 ) for all t ≥ 0. Then,EQUATIONwhere the last inequality uses a ≥ 1. Combining Eq. ( 8) and Eq. ( 9) yields the result.The results from the last two sections are applied to the CRPS and the empirical risk. Proposition 5. Under Assumption 1, for all θ ∈ Θ, the random variable S(F θ,X , Y ) is β-sub- Gaussian with β = 2(β 2 1 + β 2 2 ).Proof. We use following alternative representation of the CRPS, see Gneiting and Raftery (2007) : for fixed F ∈ P 1 (R) and y ∈ R,S(F, y) = E[|Z -y|] - 1 2 E[|Z ′ -Z|]where Z, Z ′ denote independent random variables with distribution F . The triangle inequality then impliesS(F, y) = E[|Z -y|] - 1 2 E[|Z ′ -Z|] ≤ |y| + 1 2 E[|Z + Z| -|Z ′ -Z|] ≤ |y| + E[|Z|] = |y| + m 1 (F ),with m 1 (F ) the absolute moment of F . We deduce that, for all θ ∈ Θ,EQUATIONwhere, by Assumption 1, |Y | and m 1 (F θ,X ) are sub-Gaussian with parameter β 1 and β 2 respectively. Then, Proposition 3 implies that S(F θ,X , Y ) is β-sub-Gaussian with β = 2(β 2 1 + β 2 2 ).Corollary 2. Under Assumption 1, for all θ ∈ Θ, the empirical risk Rn (F θ ) satisfies the concentration inequalityP Rn (F θ ) -R(F θ ) ≥ t ≤ 2 exp - nt 2 2β 2 , for all t ≥ 0 with β = 2(β 2 1 + β 2 2 ).Proof. The random variables Z θ,i = S(F θ,Xi , Y i ), 1 ≤ i ≤ n, are i.i.d. with expectationE[Z θ,i ] = E[S(F θ,Xi , Y i )] = R(F θ )and empirical mean1 n n i=1 Z θ,i = 1 n n i=1 S(F θ,Xi , Y i ) = Rn (F θ ).By Proposition 5, Z θ,i is β-sub-Gaussian with β = 2(β 2 1 + β 2 2 ). We can then apply the general Hoeffding's inequality (Proposition 4) and deduce the result.Before proving Theorem 1, we introduce preliminary results about the Lipschitz regularity of the CRPS and the empirical risk. Lemma 2. For all F 1 , F 2 ∈ P 1 (R) and y ∈ R, we have|S(F 1 , y) -S(F 2 , y)| ≤ 2W 1 (F 1 , F 2 )This Lemma states that the CRPS is 2-Lipschitz in the first variable with respect to the Wasserstein distance W 1 .Proof. By the definition (3) of the CRPS,S(F 1 , y) -S(F 2 , y) = R (1 {y≤z} -F 1 (z)) 2 -(1 {y≤z} -F 2 (z)) 2 dz. Using a 2 -b 2 = (a -b)(a + b), we get |S(F 1 , y) -S(F 2 , y)| ≤ R |F 1 (z) -F 2 (z)||F 1 (z) + F 2 (z) -21 {z≤y} |dz ≤ 2 R |F 1 (z) -F 2 (z)|dz.We recognize the expression of the Wasserstein distance of order 1 for probability measures on R as the L 1 -distance between their cdf, i.e. the formulaR |F 1 (z) -F 2 (z)|dz = W 1 (F 1 , F 2 ),see Bobkov and Ledoux (2019) . The result follows.We can deduce that, under Assumption 2, the theoretical and empirical risks are also Lipschitz continuous. Proposition 6. Under Assumption 2, for all θ 1 , θ 2 ∈ Θ,|R(F θ1 ) -R(F θ2 )| ≤ 2L∥θ 1 -θ 2 ∥ and | Rn (F θ1 ) -Rn (F θ2 )| ≤ 2L∥θ 1 -θ 2 ∥.Proof. Using the definition of the theoretical risk and the Lipschitz continuity of the CRPS, we have|R(F θ1 ) -R(F θ2 )| = |E[S(F θ1,X , Y ) -S(F θ2,X , Y )]| ≤ E[|S(F θ1,X , Y ) -S(F θ2,X , Y )|] ≤ 2E[W 1 (F θ1,X , F θ2,X )].Then, Assumption 2 implies|R(F θ1 ) -R(F θ2 )∥ ≤ 2L∥θ 1 -θ 2 ∥.Similarly for the empirical risk,| Rn (F θ1 ) -Rn (F θ2 )| ≤ 1 n n i=1 S(F θ1,Xi , Y i ) -S(F θ2,Xi , Y i )] ≤ 2 n n i=1 W 1 (F θ1,Xi , F θ2,Xi ) ≤ 2L∥θ 1 -θ 2 ∥.We are now ready for the proof of Theorem 1.Proof of Theorem 1. By compactness of Θ and continuity of θ → R(θ), the theoretical risk reaches a minimum on Θ and we can define θ * = arg min θ∈Θ R(F θ ). Similarly, θn = arg min θ∈Θ Rn (F θ ) is well defined. The estimation error can then be decomposed into three termsEQUATIONBy definition of θn , the second term is non-positive and we deduceEQUATION)Since Θ is compact, it can be included in some centered closed Euclidean ball BR with radius R > 0.Therefore, for all ϵ ≤ R, there exists an ϵ-net Θ ϵ of Θ such that card(Θ ϵ ) ≤ 3R ϵ K , see Devroye et al. (1996) . The term ϵ-net means that, for all θ ∈ Θ, there exists θ ϵ ∈ Θ ϵ such that ∥θ -θ ϵ ∥ ≤ ϵ. Now, for all θ ∈ Θ, we introduce the decompositionRn (F θ ) -R (F θ ) ≤ Rn (F θ ) -Rn (F θϵ ) + Rn (F θϵ ) -R (F θϵ ) + |R (F θϵ ) -R (F θ )| .By the Lipschitz properties stated in Proposition 6, the first and third terms are bounded from above by 2Lϵ. Therefore, we deduceEQUATIONwhich implies, for t ≥ 16Lϵ,P 2 sup θ∈Θ Rn (F θ ) -R (F θ ) ≥ t ≤ P sup θ∈Θϵ Rn (F θ ) -R (F θ ) ≥ t/4 .We then use Corollary 2 stating that, for all θ ∈ Θ,P Rn (F θ ) -R(F θ ) ≥ t ≤ 2 exp - nt 2 2β 2 with β = 2(β 2 1 + β 2 2 ).The union bound and the inequality card(Θ ϵ ) ≤ 3R ϵ K imply, for t ≥ 16Lϵ, P 2 sup θ∈Θ Rn (F θ ) -R (F θ ) ≥ t ≤ θ∈Θϵ P Rn (F θ ) -R (F θ ) ≥ t/4 ≤ 2 3R ϵ K exp - nt 2 32β 2 .Setting ϵ = 3R/n in this inequality and using Eq. ( 12), we deduceEQUATION)The right hand side is equal to δ for t = 32β 2 log(2n K /δ)/n, yielding the result. Note that we need this specific choice of t to satisfy t ≥ 48LR/n, or equivalently n log(2nK /δ) ≥ (48LR) 2 /c β with c β = 64(β 2 1 + β 2 2 ).Proof of Corollary 1. Setting a = log(2n K ) and b = n/c β , Eq. ( 14) can be rewritten asP R(F θn ) -inf θ∈Θ R(F θ ) ≥ t ≤ exp a -bt 2 , t ≥ 48LR/n.For t ≤ a/b, the inequality is trivial because the right hand side is larger than 1. We deduce that the inequality holds for all t ≥ 0 as soon as 48LR/n ≤ a/b which is equivalent to n log(2n K ) ≥ (48LR) 2 /c β , which holds for n large enough. Then we can apply Lemma 1 and deduce the upper bound 2 a/b for the expectation of the estimation error.C Proof of the results of Section 3.2Proof of Theorem 2. A straightforward adaptation of the proof of Proposition 5 shows that Assumption 3 implies that the random variables S( F m n,X , Y ), m = 1, . . . , M , are β-sub-Gaussian conditionally on D n with β = 2(β 2 1 + β 2 n ). Then, similarly as in Corollary 2, Hoeffding inequality impliesEQUATIONfor all t ≥ 0 and m = 1, . . . , M . With a similar reasoning as in Eq. ( 11), the regret in model selection is bounded from above byEQUATION)This together with the union bound and Eq. ( 15) impliesP R( F m n ) -min 1≤m≤M R( F m n ) ≥ t D n ≤ P max 1≤m≤M R′ N ( F m n ) -R( F m n ) ≥ t 2 D n ≤ M m=1 P R′ N ( F m n ) -R( F m n ) ≥ t 2 D n ≤ 2M exp - N t 2 16(β 2 1 + β 2 n ).The right hand side is equal to δ for t = 4 c n log(2M/δ)/N , which proves the first claim. The second claim follows by an application of Lemma 1 with a = log(2M ) and b = N/(16c n ).For the proof of Theorem 3, we need the following Lemma. Lemma 3. Consider cumulative distribution functions (G m ) 1≤m≤M and, for λ ∈ Λ, the convex aggregationG λ = M m=1 λ m G m . Then the following Lipschitz property is satisfied: for all λ 1 , λ 2 ∈ Λ, W 1 (G λ1 , G λ2 ) ≤ max 1≤m≤M m 1 (G m ) ∥λ 1 -λ 2 ∥ 1Proof. In dimension 1, the Wasserstein distance is given byW 1 (G λ1 , G λ2 ) = R |G λ1 (y) -G λ2 (y)|dy.By definition of G λ , we deduceW 1 (G λ1 , G λ2 ) = R | M m=1 (λ 1,m -λ 2,m )(G m (y) -1 {y≥0} )|dy ≤ M m=1 |λ 1,m -λ 2,m | R |G m (y) -1 {y≥0} |dy = M m=1 |λ 1,m -λ 2,m | m 1 (G m ) ≤ max 1≤m≤M m 1 (G m ) ∥λ 1 -λ 2 ∥ 1 .In the first line, we use that M m=1 (λ 1,m -λ 2,m ) = 0 so that we can introduce the indicator function 1 {y≥0} . In the third line, we useR |G m (y) -1 {y≥0} |dy = W 1 (G m , δ 0 ) = m 1 (G m ).Proof of Theorem 3. We work conditionally on D n so that F m n , 1 ≤ m ≤ M, can be seen as deterministic. The empirical risk is computed on the validation set D ′ N with cardinal N . Furthermore, Assumption 3 implies that, conditionally on D n , Y and F λ n satisfy Assumption 1 with Θ replaced by Λ and constant β 2 replaced by β n . Lemma 3 implies that, conditionally on D n , F λ x satisfy Assumption 2 with Θ replaced by Λ M and constant L replaced by√ M max 1≤m≤M m 1 ( F m n,x). As a consequence, we can apply Theorem 1 and deduce the result.Without sub-Gaussian assumptions, we cannot use Hoeffding inequality. Alternatively, we use a consequence of Rosenthal inequality providing a simple control of moments. For a reference, see Petrov (1995, Theorems 2.9 and 2.10). Proposition 7. Let Z, Z 1 , . . . , Z n be i.i.d. random variables such that mp (Z) = E[|Z-E[Z]| p ] < ∞ for some p ≥ 2. Then, E 1 n n i=1 Z i -E[Z i ] p ≤ c(p)m p (Z)n -p/2 ,with c(p) > 0 depending only on p.We deduce a simple moment bound for the empirical risk. Proposition 8. Under Assumption 4, we have, for all θ ∈ ΘE Rn (F θ ) -R(F θ ) p ≤ c ′ (p)Dn -p/2with c ′ (p) = 4 p c(p) depending only on p.Proof. Eq. ( 10) together with the upper bound |a+ b| p ≤ 2 p-1 (|a| p + |b| p ) imply |S(F θ,X , Y )| p ≤ |Y | + m 1 (F θ,X ) p ≤ 2 p-1 |Y | p + m 1 (F θ,X ) p .Taking the expectation and using Assumption 4, we deduce, with the notationZ θ = S(F θ,X , Y ), E[|Z θ | p ] ≤ 2 p-1 |Y | p + m 1 (F θ,X ) p ≤ 2 p D.By Jensen inequality, this implies the upper boundE |Z θ -E[Z θ ]| p ≤ 4 p D.The centered empirical risk can be written asRn (F θ ) -R(F θ ) = 1 n n i=1 Z θ,i -E[Z θ,i ] with Z θ,i = S(F θ,Xi , Y i ), 1 ≤ i ≤ n, i.i.d. random variables.Then the result follows from an application of Proposition 7.Proof of Theorem 4. We use the same notation as in the proof of Theorem 1. The beginning of the proof follows exactly the same lines: according to Eq. 12-13, for all ϵ-net Θ ϵ ⊂ Θ, we haveEQUATION)This impliesE R(F θn ) -inf θ∈Θ R(F θ ) p ≤ 2 p-1 (8Lϵ) p + 2 p E sup θ∈Θϵ Rn (F θ ) -R (F θ ) p .The expectation in the right hand side is bounded from above byE sup θ∈Θϵ Rn (F θ ) -R (F θ ) p ≤ θ∈Θϵ E Rn (F θ ) -R (F θ ) p ≤ (3R) K ϵ -K c ′ (p)Dn -p/2where the last inequality follows from Proposition 8 and the bound card(Θ ϵ ) ≤ (3R/ϵ) K . Hence we deduceE R(F θn ) -inf θ∈Θ R(F θ ) p ≤ 2 p-1 (8Lϵ) p + 2 p (3R) K ϵ -K c ′ (p)Dn -p/2 .Minimizing the right hand side with respect to ϵ according to Lemma 4 below and taking the p-th root, we obtainE R(F θn ) -inf θ∈Θ R(F θ ) p 1/p ≤ Cn -p/(2(p+K))where the constant C = C(K, p, L, D, R) does not depend on n and can be made explicit thanks to Lemma 4. Finally, Jensen's inequality yields the result.Lemma 4. For all a, b > 0 and p, q > 0, inf ϵ>0 aϵ p + bϵ -q = C p,q a q p+q b p p+q , with C p,q = p q q p+q + q p p p+q .Proof. A straightforward analysis of the function ϵ → aϵ p + bϵ -q shows that its derivative vanishesat ϵ = bq ap 1 p+qwhere the minimum is reached.We provide some additional results to Theorem 4 from Section 3.3 regarding model selection and convex aggregation. Our working assumption in this framework is the following.Assumption 5. For p ≥ 2, i) ∥Y ∥ p L p ≤ D; ii) conditionnaly on D n , ∥m 1 ( F m n,X )∥ p L p ≤ D n = D(D n ) for all m = 1, . . . , M .For instance, one can show with a reasoning similar to Proposition 1 that for the EMOS and DRN models, Assumption 5 holds as soon as Y and F have a finite moment of order p and X remains bounded. For the KNN and DRF models, a reasoning similar to Proposition 2 shows that the assumption holds withD n = max 1≤i≤n |Y i | p .Model selection. Here is our result on model selection under moment assumption only.Theorem 5. Under Assumption 5, we haveE R( F m n ) -min 1≤m≤M R( F m n ) D n ≤ 2 c ′ (p) max(D, D n )M 1/p N -1/2 ,with c ′ (p) the constant from Proposition 8.Before proving the Theorem, we establish the following result.Proposition 9. Under Assumption 5, we have, for all 1 ≤ m ≤ M ,E R′ N ( F m n ) -R( F m n ) p D n ≤ c ′ (p) max(D, D n )N -p/2 .Proof. We work conditionally on D n so that F m n can be seen as a deterministic quantity. The empirical risk is computed on the validation set D ′ N with cardinal N . Furthermore, Assumption 5 implies that, conditionally on D n , Y and F m n satisfy 4 with constant D replaced by max(D, D n ). Given these remarks, Proposition 9 follows by an application of Proposition 8 where F θ is replaced by F m n , Rn is replaced by R′ N and the expectation by the conditional expectation given D n .Proof of Theorem 5. Thanks to Eq. ( 16),R( F m n ) -min 1≤m≤M R( F m n ) ≤ 2 max 1≤m≤M R′ N ( F m n ) -R( F m n ) , whence we deduce E R( F m n ) -min 1≤m≤M R( F m n ) p D n ≤ 2 p E max 1≤m≤M R′ N ( F m n ) -R( F m n ) p D n ≤ 2 p M m=1 E R′ N ( F m n ) -R( F m n ) p D n ≤ 2 p M c ′ (p) max(D, D n )N -p/2. The result follows by Jensen's inequality.Here is our main result on convex aggregation under moment assumption only. Theorem 6. Under Assumption 5, we haveE R( F λ n ) -inf λ∈Λ M R( F λ n ) D n ≤ C L K max(D, D n )N -p/2 1/p+K with L = √ M max 1≤m≤M m 1 ( F m n) and constant C > 0 depending only on p and K.Proof of Theorem 6. The proof follows by an application of Theorem 4 exactly in the same way that Theorem 3 was derived from Theorem 1.The following Definition and Lemma about location/scale families are useful for the proof of Section 4. Definition 3. Let F be a probability distribution on R with zero mean and unit variance. The associated location-scale family is defined as {F m,σ , m ∈ R, σ > 0} where F m,σ is the law of m + σZ with Z ∼ F . Proposition 10 (Chhachhi and Teng 2023). For all m 1 , m 2 ∈ R and σ 1 , σ 2 > 0, we haveW 1 (F m1,σ1 , F m2,σ2 ) ≤ |m 1 -m 2 | + m 1 (F ) |σ 1 -σ 2 |.Proof of Proposition 1. Denote by F m(x;θ),σ(x;θ) the output distribution of the DRN with parameter θ ∈ Θ and input x, where m(x; θ) and σ(x; θ) are the location and scale parameters given by Eq. ( 7). Because the activation function g is assumed Lipschitz continuous, when x and θ remain in bounded sets, the functions θ → m(x; θ) and θ → σ(x; θ) are Lipschitz continuous with a Lipschitz constant bounded by some constant C uniformly in x. Then, Proposition 10 implies that the output distributions satisfyW 1 (F m(x;θ1),σ(x;θ1) , F m(x;θ2),σ(x;θ2) ) ≤ |m(x; θ 1 ) -m(x; θ 2 )| + m 1 (F ) |σ(x; θ 1 ) -σ(x; θ 2 )| ≤ C(1 + m 1 (F ))∥θ 1 -θ 2 ∥.for all θ in Θ compact and x in the support of X bounded. It follows that Assumption 2 is satisfied. Furthermore, this estimation of the Wasserstein distance implies that the absolute moment m 1 (F θ,X ) remains bounded for all θ ∈ Θ and x in the support of X. Hence m 1 (F θ,X ) is sub-Gaussian (uniformly in θ) and this proves Assumption 1. For the same reason, the absolute moment m 1 ( Fn ) = m 1 (F θn,X ) remains bounded and this implies Assumption 3 with β(D n ) constant not depending on the training set.Proof of Proposition 2. The two models KNN and DRF can be written in the form Fn,x (y) = n i=1 w ni (x)1 {Yi≤y} for suitable probability weights (w ni (x)) 1≤i≤n . Hence the absolute moment of the predictive distribution satisfiesm 1 ( Fn,x ) = n i=1 w ni (x)|Y i | ≤ max 1≤i≤n |Y i |.We can see that the random variable m 1 ( Fn,X ) remains bounded and is hence β n -sub-Gaussian with β n = max 1≤i≤n |Y i |, which proves Assumption 3.The numerical results for the airfoil data is presented in Figure 2 and Table 2 . The left plot of Figure 2 corresponds to KNN, and the middle plot to DRF. In both instances, the validation curve exhibits a clear minimum, enabling the selection parameter values k = 3 and mtry = 5 associated with CRPS-error of 2.16 and 1.47 on the validation set respectively. Then the models KNN, DRF, MS and CA are tested on the test set, and the right plot of Figure 2 shows the test error evaluated with 100 repetitions. The boxplots reveal a clear difference between the test errors for KNN and DRF in favor of DRF. In this case, MS achieves the same performance as DRF (in fact DRF is systematically selected), and a slight improve is achieved with CA. These results are confirmed by the numerical values in Table 2 . • The answer NA means that the paper does not include experiments.• If the paper includes experiments, a 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 a 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 a 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 it 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 to 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 a 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 submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example (d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe 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., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results. Justification: All the details of the numerical experiment can be found in the Models and Methodology parts of Section 5 (data description, methods and algorithms used, hyperparameters, optimizers . . . ) 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.For the sake of simplicity, we consider only the case of real-valued observation Y and of non-negative scoring rule. A more general definition can be found inGneiting and Raftery (2007). 


Decision: The paper provides concentration inequalities for the empirical risk minimization (ERM) estimator with continuous rank probability score (CRPS) as loss in the context of distributional regression, which aims to estimate the conditional distribution of a target random variable $Y$ given covariates $X$. Concentration bounds on the estimation error are provided for model fitting, model selection, and model averaging. The proof techniques are standard but the paper contains interesting results. The authors argue that the main takeaway of their theoretical results is that "the empirical risk minimization widely used in practice (be it for model fitting or model selection) is mathematically grounded". However, it might improve the paper substantially if the author(s) could use the theoretical results to provide deeper insights for the ERM methods.

Reviews: review: paper_summary: paper_summary: The paper provides theoretical guarantees for distributional regression, which aims to estimate the conditional distribution of a target random variable Y given covariates X.    These theoretical guarantees hold when the regression is learned by minimizing a particular proper scoring rule, the Continuous Rank Probability Score (CRPS), in the presence of i.i.d. data (X,Y).    Concentration bounds on the estimation error are provided for model fitting, model selection, and model averaging.
strength_weakness: strength_weakness: 1.) The experiments do not seem strongly related to the theory, other than that the experiments investigate distributional regression.  I expected the theoretical results to play a stronger role here.   Is there some way in which the theory guides the data analysis?  If not, what is the purpose of including an experiment at all? If nothing else, the experiments would seem to provide an opportunity to perform sanity checks.  For example, can the authors provide empirical evidence of consistency to the theoretical limiting values?  

2.) Generalizing #1, the authors do not provide guidance on what their theoretical results can be used for.
questions: questions: 1.) The authors write: "[...] distributional regression is way stronger and more challenging than standard regression, where only a point prediction for Y given X=x is provided which typically reduces to an estimation of the conditional expectation."  I found this statement a bit puzzling.    For instance, both frequentist and Bayesian GLM's would seem to estimate conditional distributions of a target given explanatory variables, rather than simply providing point estimates.    Could the authors clarify?

Editing notes:

1.) "Way stronger and more challenging" (line 20) is redundant, and the first phrase is too colloquial.  Streamline to "more challenging" or perhaps "much more challenging". 

2.) It is unclear why "we will mostly use the CRPS for discrete predictive distributions" (line 87).  Please clarify.

3.) Can the Wasserstein space `P_1(\R)` (defined between lines 83 and 84) be equivalently described as the set of probability distributions over `\R` with finite first moment? If so, this might be worth stating explicitly.
limitations: limitations: No limitations were provided.
paper_summary: paper_summary: The paper considers the problem of distributional regression, i.e., learning the distribution of Y conditional on X. In particular, the distributional regression problem is formulated as an empirical risk minimization problem, where standard concentration techniques are applied to obtain non-asymptotic bounds on the excessive risks.
strength_weakness: strength_weakness: After formulating the problem as an empirical risk minimization problem, applying concentration techniques appears to be very standard --- the technical contribution of the current paper is relatively weak.
questions: questions: The current paper mainly considers the case of iid data, and correspondingly leverages concentration tools for iid data. I wonder if it could be generalized to adaptively collected data, since there are parallel concentration tools for such data.
limitations: limitations: The paper has discussed its limitation.
paper_summary: paper_summary: This paper considers the problem of conditional distribution prediction. For covariate-response pair $(X, Y)$, the objective is to estimate the conditional distribution of $Y|X=x$ for all $x$. The paper provides concentration bounds for the empirical risk minimization (ERM) estimator with continuous rank probability score (CRPS) as loss. It also provides regret bounds for the model selection and model aggregation with CRPS minimization.
strength_weakness: strength_weakness: 1. The proving technique for the theorems presented in the paper is pretty standard. The work lacks technical contribution.
2. In Theorem 1, the definition of $R$ which denotes the boundary of the parameter space $\Theta$ should be added in the main text. 
3. The theorems have not provided insights for the algorithm design. The paper has not proposed new algorithms, but in line 262 it says "demonstrate the effectiveness of our proposed methods".
4. The entire experimental section does not utilize the theoretical guarantees established earlier. Additionally, there is a lack of experiments to validate the effectiveness and tightness of the theoretical bounds.
questions: questions: As given in the "Weakness" section.
limitations: limitations: As given in the "Weakness" section.
paper_summary: paper_summary: The paper studies distributional regression, a statistical technique that models not just the conditional mean of the response variable given the predictors (as in traditional regression) but the entire conditional distribution. The learning task minimizes the risk function, measured by the Continuous Rank Probability Score (CRPS). The authors aim to estimate the error bounds for the expectation of the predictive distribution. They also investigate model selection and aggregation via CRPS minimization on a validation set.
strength_weakness: strength_weakness: * Typos in the Text: The text should be checked again by the authors to remove the typos (e.g., line 88: "on can show" should be "one can show").

* Math Symbols and Abbreviations: Several math symbols or abbreviations are used in the text without prior definition, reducing the readability of the text. For instance, "EMOS" in the abstract or $\mathcal{G}$ in line 69. Additionally, some symbols are used multiple times with different meanings. For example, $\mathcal{D}$ in line 19 indicates the training set, while $\mathcal{D}$ in line 71 denotes a subset of the set of all probability measures on $\mathcal{R}$. It would be better to use symbols more carefully and ensure consistency throughout the text.

* Comparison with Existing Methods: While the paper focuses on the advantages of the proposed CRPS-error bounds for distributional regression, it would be beneficial to include a comparison with existing techniques, both in terms of theoretical properties and empirical performance. Related works and available baselines have not been properly discussed in the paper. Including such comparisons would provide a clearer context for the contributions of the proposed method and its relative performance.

* Complexity Analysis: A detailed complexity analysis is needed to understand the computational cost of the proposed method, especially in comparison with other techniques.


* Numerical Analysis Clarity: The results presented in the numerical analysis should be clearer, with a more detailed explanation of what the numbers represent and how they support the claims made in the paper. The numerical analysis is not sufficient, and not enough sensitivity analysis has been provided. Only two small datasets have been used, and the simulation does not cover all aspects of the model. It would be better if the authors designed new experiments with different datasets and settings to comprehensively validate the proposed method. This would help demonstrate the robustness and applicability of the approach across various scenarios.
questions: questions: See Weaknesses.
limitations: limitations: The limitations have not been properly discussed in the paper.

Comments: title: Response to the authors' feedback
Further comments by the reviewers?
comment: Thank you for the response. I will remain my initial score.
Thank you for your response.  My initial assessment of strengths and weaknesses persists after the rebuttal.  Therefore, I maintain my original score.
Thank you very much for your detailed response. I will keep my original score.
Thank you for your revisions and the effort you’ve put into improving the manuscript. The gap between theoretical claims and numerical experiments still exists. After careful consideration, I have decided to retain my original score.
Dear all reviewers:

Can you please respond to the rebuttal as soon as possible? Your comments will be greatly appreciated. Many thanks,

AC