CRPS, Quantile regression, Decision trees, Random forest, Distributional regression, Uncertainty quantification

1 Introduction↩︎

Unlike standard regression methods which target the conditional mean of the response to explanatory variables, distributional regression aims at estimating the entire conditional distribution of the outcome. Distributional regression is becoming more and more popular since understanding variability and uncertainty is crucial in fields like econometrics, climatology, and biostatistics. [1] laid foundational work with their generalized additive models for location, scale, and shape, in which a parametric distribution is assumed for the response variable but the parameters of this distribution can vary according to explanatory variables using linear, nonlinear or smooth functions. Another line of work relies on density estimation techniques to learn the conditional distributions with methods ranging from mixture models [2] to Bayesian ones [3]. More recently, the deep learning area gave birth to generative models achieving remarkable results, especially in the field of imaging [4]. These generative models aim at learning to sample from the training data distribution via a transformation of a simple distribution. One can mention normalizing flows [5] or diffusion models [6]. These methods have been widely applied to images and text, but recent works [7] adapt these methods to classical regression. We refer to a recent review of distributional regression approaches in [8].

While distributional regression provides a comprehensive view of the response distribution, there are cases where practitioners are interested in specific aspects of this distribution, rather than the entire probability landscape. For instance, understanding the behavior of a response variable at its extremes or at specific percentiles can be crucial for applications such as risk assessment or resource allocation. This is where quantile regression [9], [10] becomes a powerful complementary tool. Quantile regression has become a standard technique and is usually based on the minimization of the so-called pinball loss. The pinball loss \(\ell_{\tau}\) for the quantile of level \(\tau \in (0,1)\) is defined by \[\ell_{\tau}(\xi) = (\tau- \mathbb{1}_{\xi <0})\xi,\quad \xi \in \mathbb{R}\]and is such that \[q_{\tau}(x) \in {\arg \min}_{q \in \mathbb{R}} \mathbb{E}_P[ \ell_{\tau} (Y -q) \mid X=x],\] where \(q_{\tau}(x)\) is the quantile of level \(\tau\) of the distribution of \(Y\) conditional on \(X=x\). While the pinball loss is not differentiable at zero, optimization can still be performed using subgradient methods or by employing a smoothed approximation, such as the Huber loss for the \(\ell^1\) loss. This approach makes it possible to use of standard regression models to estimate specific conditional quantiles of interest. Notable examples include the use of neural networks [11], kernel methods [12], Gaussian processes [13] or splines [14]. We refer to the survey paper [15] for a more exhaustive list of quantile regression methods.

Simultaneous multiple quantile regression provides a natural bridge between full distributional regression and single quantile regression. Estimating a set of quantiles rather than a single one can be important in various applications such as decision making [16], personalized health [17] or finance [18] where professionals often need to understand the potential losses or gains at various quantiles to assess the risk of different investments or portfolios. A straightforward approach to multiple quantile regression involves fitting separate models for each quantile of interest. However, this method can be suboptimal in practice, as it fails to leverage the shared information between the models, potentially leading to poorer performance. Additionally, it may result in the undesirable phenomenon of quantile crossing, where estimated quantiles are inconsistent, namely when higher level quantiles are lower than lower level quantiles. To address these issues, regularized approaches [19] or constrained methods [20] have been proposed to couple different quantile models. To improve computational efficiency, a promising strategy is to optimize the sum of pinball losses across multiple quantiles in a single model [21], rather than fitting separate models for each quantile. The sum of pinball losses is known as the weighted interval score (WIS) in the probabilistic forecasting community [22]. Given \(M\) quantile levels \(0<\tau_1<\dots <\tau_M<1/2\), the WIS loss is defined by: \[\ell_{WIS(\tau_{1:M})}(\{q_{\tau_m}\}_{m=1}^{2M+1}, y)=\frac{2}{2M+1} \sum_{m=1}^{2M+1} \ell_{\tau_m}(y-q_{\tau_m}),\] where \(\tau_{M+1}=1/2\) and where \(\tau_{M+2}=1-\tau_{M}\), ..., \(\tau_{2M+1}=1-\tau_{1}\). The WIS is primarily a loss function on quantiles, but can actually be interpreted as a loss function over a cumulative distribution function \(F\) by rewriting it as \[\ell_{WIS(\tau_{1:M})}(F, y)=\frac{2}{2M+1} \sum_{m=1}^{2M+1} \ell_{\tau_m}(y-F^{-1}(\tau_m)).\] By increasing the number of quantiles considered, one intuitively expects the WIS loss to be sensitive to any change in the cumulative distribution function \(F\). This intuition can be formalized, as shown in [23], and it can be shown that the WIS loss converges asymptotically to the continuous ranked probability score (CRPS) as \(M\to \infty\): \[\mathrm{CRPS}(F,y) := \int_{-\infty}^{+\infty} \big( F(s) - \mathbb{1}_{s\geq y}\big)^2 ds.\] This convergence bridges the gap between quantile-based evaluation and full distributional assessment.

Both the CRPS and WIS are proper scoring rules, i.e. loss functions that are minimized, in expectation, when the forecasted distribution (for CRPS) or the forecasted quantiles (for WIS) match the corresponding true conditional distribution or quantiles [24]. While these proper scoring rules encourages well-calibrated probabilistic forecasts, their adoption in practice remains limited due to several challenges. Firstly, optimizing the risk associated to the CRPS loss can be computationally challenging since explicit solutions are only available for limited parametric cases [25] and gradients of the CRPS are only tractable under some strong assumption on the data distribution [26]. Secondly, most existing methods for minimizing the risk associated to the WIS loss require either adding constraints during training or implementing post-processing adjustments to enforce non-crossing quantiles. These additional procedures not only complicate the optimization process but might also negatively impact overall model performance. Thirdly, and perhaps most critically, many of the recent deep learning-based methods for distributional learning, while powerful, suffer from a lack of interpretability. This is a significant drawback in high-stakes applications, such as the medical field, where understanding model decisions is crucial for trust and accountability. There is, therefore, a growing need for interpretable models that can optimize CRPS or WIS losses while maintaining strong predictive performance and computational efficiency.

Among interpretable models, one of the most popular is regression trees. This can be attributed to their intuitive representation of decision-making processes, which appeals to both experts and non-experts. Their versatility in handling diverse data types renders them applicable from finance to healthcare, image recognition, and natural language processing. The construction of decision trees entails an iterative process of selecting features and thresholds to recursively partition the input space, guided by an entropy measure, aka impurity score, quantifying a level of error in each node . At each node, the algorithm assesses potential splitting thresholds for each feature, choosing the one that minimizes the overall impurity. Commonly used impurity measures are in fact directly associated with a convex loss function \(\ell\) and defined as the value of the risk evaluated at its minimizer. Mathematically, for a loss function \(\ell\), the associated impurity measure or generalized entropy \(\mathbb{H}_{\ell}( P)\) is equal to \[\label{eq:entropy95loss}\mathbb{H}_{\ell}( P) := \inf_{a} \mathbb{E}_{P}[\ell(a,Y)].\tag{1}\] This connection is well-known in the literature [27],[28], [29]. The generalized entropy associated with a specific loss function is also sometimes referred to as the Bayes risk [29], [30]. In classification, the choice of entropies remains relatively limited. In binary classification, the usual choices for entropy include the Gini entropy, the Shannon entropy, or the 0-1 entropy [31] corresponding respectively to the squared loss, the log loss and the 0-1 loss. For regression, common practice involves opting for a splitting criterion based on the squared loss or the absolute error.

The main contribution of this work lies in the development of an efficient algorithm for learning distributional regression trees that directly optimize the CRPS loss without making any assumptions on the underlying data distribution. To our knowledge, this is the first approach that enables efficient, assumption-free CRPS optimization at scale for trees. In addition, we propose efficient training algorithms for the WIS loss based on an improvement over the approach of [32]. Our algorithms update multiple conditional quantile estimates concurrently as increasing thresholds on a feature are considered. This makes it easy to define common splits for all quantiles, which produce non-crossing quantile estimates directly. For both CRPS- and pinball-based trees, instead of relying on classical empirical information gain estimates we propose to use better estimates of the population information gain that comes at no extra computational cost, based either on Leave-One-Out (LOO) or, for CRPS-based trees, on a Mallows \(C_p\)-style correction. This LOO or Mallows criteria penalize complexity and provide a data-driven stopping rule, automatically halting tree growth when splits creating leaf nodes with too few observations do not improve the estimated population risk. Together, these contributions provide practical, scalable, and interpretable methods for calibrated, non-parametric probabilistic forecasting. The proposed methods are particularly relevant for conformal prediction, where a natural choice of prediction sets is based on quantiles of predefined levels [33] or on estimated conditional distributions [34].

1.0.0.1 Related work.

To our knowledge, two approaches to learn a conditional quantile model using regression trees have been proposed in the literature:

• The first approach is the so-called Quantile Regression Forest (QRF) method [35]. This algorithm, using the squared loss, benefits from well-established and optimized packages for training regression trees. Additionally, these trees can be used to estimate quantiles of any level with the resulting quantiles guaranteed to be non-crossing (cf. [33]). However, a key drawback is that the splitting criterion used to construct the tree is well adapted to capture the fluctuations of the conditional mean and therefore not necessarily those of the conditional quantiles of interest, which potentially decreases statistical efficiency.

• [32] consider the problem of learning a single conditional quantile model using a single quantile regression trees trained with the quantile loss function, also known as the pinball loss. The entropy associated to the pinball loss is \[\begin{align} \mathbb{H}_{\ell_{\tau}}( P) &= \min_{q \in \mathbb{R}} \mathbb{E}_P[ \ell_{\tau} (Y -q)]\\ &= (1-\tau) \mathbb{E}_P[(q-Y)_+] + \tau \mathbb{E}_P[(Y-q)_+], \end{align}\] where for any real number \(x\), \((x)_+=\max(x,0)\). In the empirical setting where we have observations \(\boldsymbol{y}=\{ y_i \}_{i=1}^n\) sampled from \(P\), the empirical entropy is \[\begin{align} \label{eq:empi95pinball} H_{\ell_{\tau}}(\boldsymbol{y}) :=\min_{q \in \mathbb{R}} \frac{1}{n}\sum_{i=1}^n \ell_{\tau} (y_i -q)=\frac{1}{n}\sum_{i=1}^n \ell_{\tau} (y_i -\hat{q}_{\tau}), \end{align}\tag{2}\] where \(\hat{q}_{\tau} = y_{(\lceil\tau n\rceil)}\) is the empirical quantile of level \(\tau\) of the sequence \(\boldsymbol{y}=(y_i)_{i=1}^n\). The tree is then trained using the standard CART algorithm [36] with the generalized entropy from Eq.@eq:eq:empi95pinball : at each node, the data is split by selecting the feature and the threshold that minimize the entropy within the resulting subsets. This process iteratively creates nodes that represent the most homogeneous possible groups by reducing overall entropy. The procedure to find the data split maximizing the information gain defined with the entropy from Eq.@eq:eq:empi95pinball at a given node of the tree with \(n\) samples and \(d\) features has a time complexity of order \(\mathcal{O}(dn\log n)\) (cf. [32]).

In this paper, we leverage similar ideas as those of [32] to propose two efficient algorithms to learn distributional regression trees for the WIS and CRPS losses.

1.0.0.2 Outline.

In this work, we first introduce two novel algorithms to learn distributional regression trees (Sec. 2): the first minimizes the empirical entropy associated to the WIS loss, which we call Pinball Multiple Quantiles Regression Trees (PMQRT); the second algorithm learns distributional regression trees (CRPS-RT) for the CRPS loss. To enhance the algorithm, we derive unbiased leave-one-out estimates of the population risk for both PMQRT and CRPS, and further provide a Mallows \(C_p\)-style correction for CRPS; notably, these corrections incur essentially no additional computational cost (Sec. 2.4). Detailed algorithmic descriptions, along with illustrative examples, are provided in the supplementary material. In Section 3, we evaluate the empirical performance of random forests of PMQRTs (PMQRF) and CRPS-RTs (CRPS-RF) through a series of experiments. We demonstrate the scalability of the proposed algorithms and that by using a CRPS-RF, we can develop a non-parametric method for learning conditional distribution functions that rivals existing approaches.

1.0.0.3 Software Package and Algorithms.

This paper builds upon several fundamental data structures, namely Weight-Balanced Binary Trees, Fenwick Trees, and Heaps. To assist readers who may be less familiar with these tools, we provide in Section 5 concise definitions and concrete examples illustrating how each data structure is used within our algorithms.

All methods presented in this work are available through a publicly released software package called DisTreebution¹, which includes complete documentation and usage examples. Section 7 of the Appendix provides a detailed description of the key features and configurable options of our Python package.

2 Learning distributional regression trees↩︎

This section introduces algorithms designed to construct regression trees with common splits for modeling multiple quantiles or the entire conditional distribution. Unlike traditional approaches that treat different quantiles independently, our methods seek to identify shared split points to improve interpretability and computational efficiency. We propose two tree-learning algorithms tailored to distinct loss functions: one for the WIS and another for the CRPS. The following subsections outline these methods and address the challenges inherent in optimizing these distinct objectives.

2.1 Notations↩︎

To evaluate potential splits based on a single input feature, the vector of output values must first be sorted by increasing order of that input feature. For each possible split index \(s\), we then consider the subset of output values corresponding to the first \(s\) elements of the sorted sequence and compute an empirical quantile - that is, an order statistic - of these labels. This motivates the introduction of the following notations: Consider an arbitrary sequence \(\boldsymbol{a} = (a_1,\dots,a_n)\) of real values of size \(n\) (here typically label values sorted according to a given input feature). For any \(s\in[n]\), we denote by \(\boldsymbol{a}^{(s)}:=(a_1,\dots,a_s)\) the prefix sequence of \(\boldsymbol{a}\). For any \(j\in [n]\), \(a_{(j)}\) corresponds to the \(j\)-th largest value of the sequence \(\boldsymbol{a}\), while for any \(i \in [s]\), \(a^{(s)}_{(i)}\) refers to the \(i\)-th largest value of the prefix sequence \(\boldsymbol{a}^{(s)}\).

2.2 Pinball multi-quantile regression tree↩︎

2.2.0.1 Extracting multiple quantiles from a single tree.

When estimating multiple quantiles, a natural option is to use QRF, since a single forest can be queried at any level. However, because QRF relies on trees learned with the square loss, it is primarily tailored to conditional mean estimation and may yield suboptimal quantile estimates. Another option is to train a separate tree for each quantile level \(\tau\), following [32], but this does not guarantee non-crossing quantiles and produces a different partition of the feature space for each \(\tau\), which may reduce statistical efficiency and is undesirable for applications such as group conformal prediction (cf. Section 3.4).

To bypass these limitations, we present an approach that takes the best of the ones above-mentioned. We propose learning a tree that produces all conditional quantiles estimates in each leaf and that most importantly relies on a common split for all quantiles considering the sum of pinball losses \(\ell_{\tau_1,\ldots,\tau_m}\) defined for any \({\boldsymbol{\tau}},\mathbf{q} \in \mathbb{R}^M\) as \[\ell_{\tau_1,\ldots,\tau_m}(\mathbf{q},y):=\sum_{m=1}^M \ell_{\tau_m}(y-q_m).\] The multi-dimensional optimization problem over \(\mathbf{q}\in\mathbb{R}^M\) that defines the empirical entropy associated with \(\ell_{\tau_1,\ldots,\tau_m}\) is separable, which leads to \[\label{eq:multi95entropy}H_{\ell_{\tau_1,\ldots,\tau_m}}(\mathbf{y})=\sum_{m=1}^M H_{\ell_{\tau_m}}(\mathbf{y}).\tag{3}\] Relying on the CART algorithm, training a regression tree considering the loss \(\ell_{\tau_1,\ldots,\tau_m}\) requires to compute the entropies \(H_{\ell_{\tau_1,\ldots,\tau_m}}(\boldsymbol{y}^{(s)})\) associated with any prefix sequence \(\boldsymbol{y}^{(s)}:= (y_1,\dots, y_s)\) from a given sequence of values \(\boldsymbol{y}=(y_1,\dots,y_n)\). We end up with a problem similar to [32] but with multiple quantiles at the same time.

2.2.0.2 Efficient pinball loss node splitting with max heaps.

The algorithm proposed in [32] to learn a quantile regression tree that minimizes the empirical entropy associated with a single pinball loss \(\ell_{\tau}\) uses two heaps (cf. [37]). At a given node in the tree, with data \(\{(\boldsymbol{x}_i, y_i)\}_{i=1}^n\), the goal is to find the threshold that provides the best information gain using the entropy from Eq.@eq:eq:empi95pinball for each feature \(j\). Assuming without loss of generality that \(x_{i,j} \leq \dots \leq x_{n,j}\) for some feature \(j\), identifying the optimal threshold reduces to computing the entropies for all prefix and suffix sequences of \(\boldsymbol{y}\). Since the suffix sequences of \(\boldsymbol{y}\) correspond to prefix sequences of \((y_{(n)}, \dots, y_{(1)})\), it is sufficient to explain how [32] computes efficiently the entropies of the sequences \(\boldsymbol{y}^{(s)}\) for all \(s \in [n]\). As shown in [32], the entropy \(H_{\ell_{\tau}}(\boldsymbol{a})\) for any nondecreasing sequence \(\boldsymbol{a}\) of length \(s\) can be computed in \(\mathcal{O}(1)\) time using \(a_{\lceil\tau(s-1)\rceil}\), \(a_{\lceil\tau s\rceil}\), \(\sum_{i=1}^{\lceil\tau(s-1)\rceil} a_{i}\), and \(\sum_{i=1}^s a_i\). The authors show that these quantities can be efficiently tracked by initializing two empty heaps - one max-heap \(\mathcal{H}_{low}\) and one min-heap \(\mathcal{H}_{high}\) - and by sequentially adding elements \(y_s\) from \(s=1\) to \(n\) to them. At each step \(s\), the update of the two heaps ensures that \(\mathcal{H}_{low}\) contains all elements \((y_{(1)}, \dots, y_{(\lceil\tau(s-1)\rceil)})\) and \(\mathcal{H}_{high}\) contains all elements \((y_{(\lceil\tau s\rceil)}, \dots, y_{(s)})\). Querying the maximum value in \(\mathcal{H}_{low}\) and the minimum value in \(\mathcal{H}_{high}\) can be done in constant time while the heaps update require \(\mathcal{O}(\log(s))\) time.

From Eq. 3 , the empirical entropy of a sum of pinball losses \(\ell_{\tau_1, \dots, \tau_M}\) is the sum of the empirical entropies for each individual pinball loss. As a consequence, a straightforward approach to estimate quantiles at multiple levels using a single tree is to apply the method from [32] independently for each quantile level. However, this brute-force approach incurs a space complexity of \(\mathcal{O}(Mn)\), making it inefficient as the number of quantile levels \(M\) increases.

To address this limitation, we propose a new algorithm whose space complexity remains of the order of the size of the dataset, regardless of the number of estimated quantiles.

2.2.0.3 Efficient pinball multi-quantile loss node splitting with a set of min-max heaps.

Using the decomposition Eq. 3 , one can deduce from [32] that the empirical entropy \(H_{\ell_{\tau_1,\ldots,\tau_m}}(\boldsymbol{y}^{(s)})\) of any prefix sequence \(\boldsymbol{y}^{(s)}\) can be obtained by keeping track of \(i)\) the empirical quantiles \(\widehat q_{\tau_1},\dots,\widehat q_{\tau_M}\) of level \(\tau_1,\dots,\tau_M\) of \(\boldsymbol{y}^{(s)}\) and \(ii)\) of the cumulative sums \(\left( \sum_{i\leq \ceil{\tau_m s}}y^{(s)}_{(i)}\right)_{m=1}^M\) and \(\sum_{i\leq s}y^{(s)}_{(i)}\). We propose an efficient algorithm to keep track of these quantities relying on \(M+1\) min-max heaps: \(\mathcal{H}_{1},\dots,\mathcal{H}_{M}, \mathcal{H}_{M+1}\). Let us recall that a min-max heap is a variant of a binary heap allowing for constant-time access to both the minimum and maximum elements. On even levels, the heap adheres to min-heap properties, while odd levels follow max-heap properties. This ensures that the root always holds the minimum, and its children represent the maximum elements. Insertions and deletions have logarithmic time complexity, providing a balance between efficient extremal element access and effective heap modification. As for standard regression trees, we increase the index \(s\) from 1 to \(n\) to cover all the prefix sequences \(\boldsymbol{y}^{(s)}\) of interest to determine the position of the split at the current node. The \(M+1\) min-max heaps are updated at each iteration \(s\) to ensure that \(\mathcal{H}_{1}\) contains all the elements \(\{y^{(s)}_{(1)},\dots y^{(s)}_{(\ceil{\tau_1 s})} \}\), \(\mathcal{H}_{2}\) contains all the elements \(\{y^{(s)}_{(\ceil{\tau_1 s}+1)},\dots y^{(s)}_{(\ceil{\tau_2 s})}\}\), \(\dots\) , \(\mathcal{H}_{M}\) contains all the elements \(\{y^{(s)}_{(\ceil{\tau_{M-1} s}+1)},\dots y^{(s)}_{(\ceil{\tau_M s})}\}\), and \(\mathcal{H}_{M+1}\) contains all the elements \(\{y^{(s)}_{(\ceil{\tau_M s}+1)},\dots y^{(s)}_{(s)}\}\). Querying the quantile of level \(\tau_m\) simply involves accessing the maximum of the heap \(\mathcal{H}_m\), or the minimum of the heap \(\mathcal{H}_j\) for \(j:=\arg \min \{i >m\mid \mathcal{H}_i\neq \emptyset\}\) in case \(\mathcal{H}_m\) is empty. Our min-max heap structures are also designed to store the total sums of elements in each heap, facilitating access to cumulative sums of the form \(\sum_{i\leq \lceil\tau_m s\rceil}y^{(s)}_{(i)}\) in \(\mathcal{O}(m)\) time by querying the total sums of heaps \(\mathcal{H}_1\) to \(\mathcal{H}_m\). Overall, the worst time complexity for a node split with \(n\) data points and \(d\) features is \(\mathcal{O}(Md n\log n)\) as illustrated in Table 1. The detailed algorithm can be found in Section 2.1 of the supplement.

2.3 CRPS distributional regression tree↩︎

In the previous section, we presented a computationally efficient method for learning a pinball multi-quantile regression tree. Such a tree has the advantage of leveraging information from different parts of the conditional distribution to build the partition, ensuring non-crossing quantile estimates, and offering a flexible tool for modern uncertainty quantification, ideal for applications such as group conformal methods, see Section 3.4. However, a Pinball multi-quantile regression tree focuses only on specific quantiles, limiting its ability to represent the full conditional distribution. As a result, it may overlook key aspects of the data distribution, affecting predictive accuracy. Furthermore, certain methods, such as the distributional conformal prediction method proposed in [34], require querying quantiles at arbitrary levels, and constraining the model to a fixed set of predefined quantiles may hinder its effectiveness.

To bypass these limitations, we propose to learn trees that minimize the empirical entropy associated with the Continuous Ranked Probability Score (CRPS), namely \[H_{\mathrm{CRPS}}(\boldsymbol{y}):=\min_F \frac{1}{n} \sum_{i=1}^n \mathrm{CRPS}(F,y_i) = \frac{1}{n} \sum_{i=1}^n \mathrm{CRPS}(\widehat F_n,y_i) , \label{eq:CRPS95empirical}\tag{4}\] where \(\widehat F_n\) is the empirical CDF of the observed data \(\boldsymbol{y}\). The CRPS is a strictly proper scoring rule, as detailed in [24], incentivizing forecasters to provide accurate and well-calibrated uncertainty estimates, making it a valuable tool in probabilistic forecasting evaluations.

2.3.0.1 CRPS-Based tree training: Theoretical foundations and design.

To construct a distributional regression tree based on the generalized entropy induced by the CRPS loss, one must identify the split at a given leaf that results in the largest impurity gain. Given a leaf with data \(\left\{(\boldsymbol{x}_i, y_i)\right\}_{i=1}^n\), this process involves computing the entropies for every feature \(j\) using every prefix sequence \(\boldsymbol{y}^{(s)} = \left\{y_1, \dots, y_s\right\}\) and every suffix sequence \((y_{s+1}, \dots, y_n)\) of the data sequence \(\boldsymbol{y}\), assuming for simplicity that \(x_{1,j} \leq \dots \leq x_{n,j}\). Relying on the entropy expression from Eq. 4 , this requires constructing the empirical CDFs for all prefix and suffix sequences of \(\boldsymbol{y}\) for each feature \(j\). This process is computationally intensive: specifically, the construction of empirical CDFs for the prefix sequences of a feature \(j\) requires sequentially inserting a new element \(y_{s+1}\) into the current sorted version of the list \(\boldsymbol{y}^{(s)}\), with a time complexity of \(\mathcal{O}(s)\). Consequently, the basic construction of a distributional regression tree based on the CRPS loss results in a node split with a time complexity of \(\mathcal{O}(d n \log(n) \sum_{i=1}^n s) = \mathcal{O}(d n^3 \log n)\), which is not feasible for practical use.

To efficiently train a tree using the CRPS loss, reducing computations is key. Due to the sequential structure of the computation of entropies associated with prefix (resp. suffix) sequences to perform a node split, one can investigate the possibility to rely on iterative schemes to reduce the computational burden of the algorithm. Theorem 1 presents an iterative formula allowing us to compute the entropies associated to all prefix sequences for a given feature \(j\). Note that a similar formula is directly obtained from Theorem 1 for suffix sequences by flipping the sign of each entry of the feature vectors \(\left\{\boldsymbol{x}_i\right\}_{i=1}^n\). In what follows, we denote by \(h^{(s)}_{\text{mirror}}\) the counterpart of the entropies \(h^{(s)}\) from Theorem 1 computed under this sign-flipped transformation of the input features.

Theorem 1 shows that training a distributional regression tree using the CRPS loss is achievable if the insertion ranks \(\big(r_s\big)_{s \in [n]}\) and cumulative sums \(S^{(s)}_{r_s}\), \(S^{(s)}\), and \(W^{(s)}\) can be computed efficiently. To this end, we introduce the CRPS-Optimal Split algorithm, which leverages tailored data structures to compute the entropy of all candidate splits for any feature \(j\) in \(\mathcal{O}(n \log n)\) time. The procedure is presented in Algorithm 2 and described in detail in the remainder of this section. To efficiently compute the ranks \(r_s\), the algorithm uses a Weight-balanced tree. This is a type of self-balancing binary search tree that keeps its shape roughly balanced as elements are inserted or removed. Because of this balance, operations like insertion, deletion, and searching can all be done in \(\mathcal{O}(\log n)\) time.

In our case, we use the Weight-balanced tree to determine the rank of each element as it is inserted — that is, how many elements in the current set are smaller than the new one. This can be viewed as inserting each element into a dynamically maintained sorted list while simultaneously determining its sorted position (rank), without requiring a full sort of the list. This rank computation is also performed in \(\mathcal{O}(\log n)\) time thanks to the structure of the tree.

After computing the ranks \((r_s)_s\), we use a Fenwick tree (also known as a binary indexed tree) to compute the cumulative sums \(S^{(s)}_{r_s}\). A Fenwick tree is a clever data structure that allows us to efficiently update values and calculate prefix sums — that is, the sum of values up to a certain index. Unlike a regular array, which would need up to \(n\) steps to compute a prefix sum, a Fenwick tree can compute it in just \(\mathcal{O}(\log n)\) time. It achieves this by storing partial sums in a compact binary representation, enabling efficient updates and queries. An illustration of how the Fenwick Tree works is shown in Figure 1.

2.3.0.2 CRPS-Based tree training: Detailed procedure.

For a given feature \(j\), the process begins by sorting the list of feature entries \(\left\{x_{i,j}\right\}_{i=1}^n\), which requires \(\mathcal{O}(n \log n)\) time. Once sorted, the insertion ranks \(\big(r_s\big)_{s \in [n]}\) can be pre-computed using a weight-balanced tree [38]. Weight-balanced trees offer logarithmic performance similar to that of red-black and AVL trees but also provide the advantage of indexing elements in logarithmic time. Lemma 1 shows that these ranks can be pre-computed in \(\mathcal{O}(n \log n)\) time.

Once the ranks \(\big(r_s\big)_{s \in [n]}\) are computed, we need to obtain the entropies \(H_{\mathrm{CRPS}}(\boldsymbol{y}^{(s)})\) for any prefix sequence \(\boldsymbol{y}^{(s)}\) using the iterative scheme presented in Theorem 1. For this, it is sufficient to query efficiently the sums \(S^{(s)}_{r_s}\), \(S^{(s)}\) and \(W^{(s)}\). Tracking \(S^{(s)} = \sum_{i=1}^s y_i\) across iterations is straightforward, as it simply involves adding the new element \(y_s\) to \(S^{(s-1)}\) at each step \(s\). Moreover, we observe that \[\label{eq:key95update95CRPS}W^{(s+1)} = W^{(s)} +r_{s+1}y_{s+1} +\sum_{i=r_{s+1}}^s y^{(s)}_{(i)}=W^{(s)} +(1+r_{s+1})y_{s+1} +S^{(s)}-S^{(s+1)}_{r_{s+1}},\tag{5}\] where we used that \(\sum_{i=r_{s+1}}^s y^{(s)}_{(i)}=\sum_{i=r_{s+1}}^{s+1} y^{(s+1)}_{(i)}-y_{r_{s+1}}^{(s+1)}=S^{(s+1)}_{r_{s+1}}-y_{s+1}\). Thus, maintaining \(W^{(s)}\) comes at no additional cost if partial sums of the form \(S^{(s)}_{r_s}\) can be queried efficiently. This highlights that efficiently tracking these sums is the key to designing an algorithm that constructs regression trees for the CRPS loss. Denoting by \(\sigma\) a permutation sorting the list \(\boldsymbol{y}\) in non decreasing order and having in mind that \(y_s=y^{(s)}_{(r_s)}\), the crucial observation is: \[\label{eq:motiv95fenwick}S^{(s)}_{r_s}=\underbrace{\sum_{i=1}^{r_s} y^{(s)}_{(i)}}_{\text{r_s smallest elements of \mathbf{y}^{(s)}}} = \underbrace{\sum_{i=1}^s \mathrm 1_{\sigma^{-1}(i)\leq \sigma^{-1}(s)} y_i}_{\text{Elements of \mathbf{y}^{(s)} smaller than y_s}}=\underbrace{\sum_{j=1}^{\sigma^{-1}(s)} \mathrm 1_{\sigma(j)\leq s} y_{\sigma(j)}}_{\substack{\text{Elements of \boldsymbol{y} smaller}\\\text{than y_s with index \leq s}}}.\tag{6}\] Building on Eq.@eq:eq:motiv95fenwick , we rely on a customized Fenwick tree to efficiently track the sums \(S^{(s)}_{r_s}\). A Fenwick tree, as shown in Figure 1, is a binary-indexed data structure that efficiently supports prefix sum and insertion in logarithmic time. In Definition 1, we introduce the tailored Fenwick tree \(\mathcal{F}\), which serves as the cornerstone of our algorithm, enabling efficient node splits.

By construction at any iteration \(s\), \(\mathcal{F}\) is a Fenwick tree for the list \((\mathrm 1_{\sigma(i)\leq s} y_{\sigma(i)} )_{i\in [n]}\). In Lemma 2, we demonstrate that \(\mathcal{F}\) enables tracking the cumulative sums of the form \(S^{(s)}_{r_s}\) and \(S^{(s)}\) which, as a by product, allows for monitoring \(W^{(s)}\) throughout the iterative process. The proof of Lemma 2 is provided in Section 3.2 of the supplement.

In Section 2.2 of the supplement, we present in details our algorithm to perform a node split to grow a CRPS-RT. Overall, the computational time required to perform a node split scales as \(\mathcal{O}(dn\log n)\) where \(n\) is the number of data points in the node and \(d\) is the number of features. The time complexity of the two algorithms presented in this paper are given in Table 1.

2.4 Estimating the population entropies↩︎

To train PQRTs, PMQRTs or CRPS-RTs, we proposed in the previous sections to rely on the standard CART algorithm [40] which is based on empirical entropies. Unfortunately, these are biased estimates of the population entropies (i.e., associated with the population risk), because the same data is used both to estimate the empirical quantiles or empirical cumulative distribution functions and to evaluate the corresponding in-sample empirical risk. As a result, the empirical criterion computed at each step of the tree-growing procedure is typically an optimistically biased estimate of the target prediction risk, and therefore induces a biased estimate of the associated information gain. A naive way to reduce this bias would be to split the sample into two disjoint parts, one used for estimation and the other for evaluation. Such a strategy is, however, statistically inefficient, since only a fraction of the available observations is used for each task at every candidate split.

To address this issue, we consider two complementary bias-correction strategies: leave-one-out (LOO) and a Mallows \(C_p\)-style optimism correction. Using the LOO or Mallows estimate of the information gain provides a data-driven stopping criterion: if all corrected information gains are non-positive, the current node is not split. In contrast, this situation can only arise in trivial settings when using the standard information gain estimate, because its optimistic bias ensures that the gain values are always non-negative.

Despite their conceptual appeal, the use of LOO and Mallows corrections in regression trees could, in principle, introduce a computational overhead. Indeed, applying these corrections at each candidate split could be costly, since the information gain would need to be evaluated for all possible splits when determining the best node split. In this section, we show that a LOO (resp. Mallows) estimate of the information gains can be obtained at no additional computational cost for the WIS and CRPS losses (resp. CRPS loss). Before presenting the formula for computing these estimates, we make a brief and subtle clarification regarding the distinction between the LOO and Mallows targets.

2.4.0.1 LOO vs Mallows \(C_p\)-style estimate of the information gains.

The targets of the LOO and Mallows estimates differ slightly. Consider a node containing \(n\) observations \((y_1,\dots,y_n)\) and let \(\hat{f}_{k:j}^{\mathrm{ERM}}\) denote the empirical risk minimizer for the loss \(\ell\) fitted on the subsample \((y_k,\dots,y_j)\). In the PMQRT setting, \(\hat{f}_{k:j}^{\mathrm{ERM}}\) is the vector of empirical quantiles associated with the prescribed quantile levels, whereas in the CRPS-RT setting it is the empirical cumulative distribution function of the subsample. We write \[\widehat R_{k:j}(\hat{f}_{k:j}^{\mathrm{ERM}}):=H_\ell(\mathbf{y}_{k:j})\] for the corresponding empirical risk (or impurity), and \(R(\hat{f}_m^{\mathrm{ERM}})\) for the population risk of the same learning procedure trained on \(m\) observations. For a candidate split producing children of sizes \(k\) and \(n-k\), a Mallows \(C_p\)-style correction directly targets the expected information gain associated with the full-sample models, namely \[\mathbb{E}\!\left[R(\hat{f}_n^{\mathrm{ERM}})-\bigl(R(\hat{f}_k^{\mathrm{ERM}})+R(\hat{f}_{n-k}^{\mathrm{ERM}})\bigr)\right],\] by correcting the empirical information gain with the appropriate optimism terms. By contrast, LOO evaluates each node with one observation removed and therefore targets the analogous quantity for models trained on \(n-1\), \(k-1\), and \(n-k-1\) observations. Hence, LOO is not exactly unbiased for the target full-sample information gain, although the discrepancy is typically small when node sizes are large. A more detailed discussion is provided in the supplementary material. Our numerical experiments described in Section 3.1 have shown that building trees using either LOO estimates of the information gains or a Mallows correction leads to very consistent results.

2.4.0.2 LOO estimate of the information gains.

For both the CRPS loss and the sum of pinball losses, we derive a LOO correction that yields an unbiased estimate of the corresponding empirical entropy while remaining computationally cost-free. In particular, the naive implementation of LOO would appear prohibitively expensive, since it requires recomputing the impurity criterion after removing each observation in turn. Nevertheless, Propositions 1 and 2 show that, in both cases, the resulting LOO criterion admits a simple closed-form expression and can be evaluated with no additional computational overhead.

Proposition 2 shows that the LOO estimate of empirical entropy for the pinball loss, \(H_{\ell_{\tau}}^{\mathrm{LOO}}(\boldsymbol{y})\), can be derived from \(H_{\ell_{\tau}}(\boldsymbol{y})\) by adding a debiasing term that depends on \(y_{(r^*-1)}\), \(y_{(r^*)}\), and \(y_{(r^*+1)}\), where \(r^* := \lceil n\tau \rceil\). Since these values can be retrieved in constant time using the min-max heap structures described in Section 2.2, a node split, whether using LOO or not, is still performed in \(\mathcal{O}(dMn \log n)\) time when considering a sum of \(M\) pinball losses. The proofs of Propositions 1 and 2 are provided in the supplement.

2.4.0.3 Mallows-style correction for the information gains.

In the specific case of the CRPS loss, we additionally derive a Mallows \(C_p\)-style optimism correction. As presented in Proposition 3 (see the Supplement for the proof), this correction can be computed at essentially no additional computational cost and provides an unbiased estimate of the CRPS information gain at each candidate split. By contrast, for the sum of pinball losses, no computationally cheap Mallows-type correction is currently available to the best of our knowledge. While a Mallows-style correction for the pinball loss might exist in principle, deriving one without introducing substantial additional computational overhead appears substantially more challenging, and we therefore leave this question for future work.

3 Experiments↩︎

In this section, we investigate the empirical performance of the different algorithms introduced in this paper when they are used within random forests based on the principle of bagging. In a first experiment, we assess the usefulness of the LOO and Mallows estimates of the population entropies. In a second part, we study the empirical computational efficiency of Algorithm 2. In a third part, we compare our methods with several alternative approaches, in particular QRF [35]. Finally, we present an application of the distributional regression trees introduced in this work to group conformal prediction.

When using trees, there are two natural aggregation schemes, depending on whether the primary focus is on conditional quantiles or on conditional CDFs. Given a new input \(x\), a first approach follows the aggregation strategy of QRF, where a conditional CDF \(\hat{F}(\cdot \mid X=x)\) is obtained by averaging the empirical CDFs \(\hat{F}_k(\cdot \mid X = x)\) produced by each tree \(k\) in the forest. More precisely, for each tree \(k\), one considers the empirical CDF of the observations contained in the leaf reached by dropping \(x\) from the root of the tree, and defines \[\hat{F}(y \mid X = x) = \frac{1}{K} \sum_{k=1}^K \hat{F}_k(y \mid X = x).\] We refer to this first approach as Distributional Bagging (DB).

A second approach consists in averaging the conditional quantiles produced by each tree. Specifically, for each tree \(k\), let \(\hat{q}_{k,\beta}(Y \mid x)\) denote the empirical quantile of level \(\beta\) computed from the observations in the leaf reached by dropping \(x\) from the root of tree \(k\). One then defines \[\hat{q}_{\beta}(Y \mid x) = \frac{1}{K} \sum_{k=1}^K \hat{q}_{k,\beta}(Y \mid x).\] We refer to this second approach as Quantile Bagging (QB).

Although DB is naturally formulated as an aggregation scheme for conditional CDFs and QB as an aggregation scheme for conditional quantiles, both approaches can in fact be used to recover either object. In particular, DB directly provides an estimate of the conditional CDF \(\hat{F}(\cdot \mid X=x)\), from which conditional quantiles can be obtained by querying \(\hat{F}\). Conversely, QB directly provides conditional quantile estimates as a function of the level \(\beta\), and a conditional CDF can be recovered by considering the corresponding stepwise inverse representation induced by the collection of conditional quantiles.

Both aggregation schemes were implemented and compared in our experiments. In all experiments involving PMQRF or CRPS-RF, we used the QB aggregation scheme. This is motivated by preliminary experiments that show that QB performs comparably or slightly better (cf. Section 6).

3.1 CDF estimation improvement with an unbiased CRPS risk estimate↩︎

To demonstrate the importance of using LOO for calculating information gain, we consider a simple example where the covariate \(X\) is uniformly sampled from the interval \((0, 10)\), and \[\begin{align} \label{eq:simu95gamma} Y \mid X &\sim \text{Gamma}(\text{shape} = \sqrt{X}, \text{scale} = \min\{\max\{X, 1\}, 6\}) \end{align}\tag{7}\] Figure 3 shows the estimated conditional CDFs using a forest of 100 CRPS-RTs, comparing the results with and without LOO used for entropy calculation. On this example, applying LOO produces a significantly better estimate of the conditional CDFs. We perform a similar experiment for PMQRF, demonstrating a significant improvement in quantile estimates achieved through the use of LOO. The results of these simulations are presented in Section 4 of the supplement.

We extend the analysis by quantifying the statistical implications of using LOO and Mallows-style corrections to compute the information gain. We consider CRPS-RF models trained with 100 trees on \(600\) data points sampled from Eq. 7 , exploring various values of the stopping-criterion hyperparameters \(max\_depth\) (maximum depth of trees) and \(min\_samples\_split\) (the minimum number of samples required to split an internal node). For each configuration, we train three versions of CRPS-RF: one using LOO, one using the Mallows \(C_p\)-style correction, and one without any correction.

Once trained, we query the estimated quantiles \(\hat{q}_{0.05}\) and \(\hat{q}_{0.95}\) (corresponding to the 5% and 95% quantiles, respectively) on a test set of 1,000 data points, and compute the marginal coverage induced by the predictive intervals \([\hat{q}_{0.05}, \hat{q}_{0.95}]\). The results are reported in Table 2.

Empirically, when the stopping rules are sufficiently restrictive (i.e., small \(max\_depth\) and/or large \(min\_samples\_split\)), the marginal coverages of CRPS-RF models trained with LOO, Mallows, or without correction are comparable. Conversely, when the stopping rules are less restrictive (i.e., large \(max\_depth\) and small \(min\_samples\_split\)), models trained without correction tend to exhibit marginal coverages that deviate more strongly from the target level of 90%, whereas models trained with either LOO or Mallows maintain much more stable coverage.

In the absence of correction, tree growth generally continues until a predefined stopping criterion is satisfied, since each node split typically produces a positive information gain, except in trivial cases. In contrast, both LOO and Mallows corrections yield improved estimates of the population entropies. Once the corrected information gains become non-positive for all possible splits of a node, the splitting process naturally terminates, even if the explicit stopping criteria have not been reached. Consequently, incorporating either LOO or Mallows-type corrections effectively implements a data-driven, automatic stopping rule.

3pt

3.2 Computational time↩︎

In this section, we conduct experiments to assess the scalability of our method to learn regression trees minimizing the CRPS loss. First, we consider random vectors \((y_1,\dots, y_n)\) with i.i.d.entries sampled from a white noise Gaussian distribution and we compute the empirical entropy \(H_{\mathrm{CRPS}}(\boldsymbol{y}^{(s)})\) associated to all prefix sequences \(\boldsymbol{y} ^{(s)}= (y_1,\dots,y_s)\) for all \(s \in [n]\), namely \[H_{\mathrm{CRPS}}(\boldsymbol{y}^{(s)}) := \frac{1}{s} \sum_{i=1}^s\mathrm{CRPS}(\widehat F^{(s)},y_i),\] where \(\widehat F^{(s)}\) is the empirical CDF of sequence \(\boldsymbol{y} ^{(s)}\). We compare with a method computing all CRPS independently using available Python packages. Figure 4 demonstrates the computational time required for a node split during the training of a regression tree that minimizes the CRPS loss. When naively computing all CRPS values, the time complexity scales empirically as \(O(n^{2.55})\). In contrast, the method proposed in this paper has a computational time complexity that scales as \(O(n \log n)\), as theoretically demonstrated.

3.3 Comparing CRPS-RF and PMQRF with QRF↩︎

In this section, we describe the numerical experiments conducted to evaluate the performance of our models. We selected six datasets from the UCI Machine Learning Repository to cover a variety of real-world applications. The total number of samples \(n\) and features \(d\) for each dataset are summarized in Table 3, with URLs pointing to the corresponding publicly available resources given in the Supplement.

For each dataset, we randomly sampled 1,000 data points without replacement to generate the training sets, and 3,000 other samples used as the test set to compute the prediction error. We consider forests of PMQRTs or CRPS-RTs with 100 trees. Trees in the forests are trained by sampling uniformly at random without replacement \(60\%\) of the training data. This procedure was repeated for 300 different training sets, each generated by sampling the data points randomly. We computed the mean and standard deviation of the CRPS on the test set across these 300 experiments, which are reported in Table 4. For each method, the CRPS was computed on the test set by querying quantiles with level \(\mathcal{Q}_{0.05}:=\{0.05\times k, \quad k\in[19]\}\) and building the corresponding CDF. In the following, we present the different models considered in our experiments:

• We train CRPS-RFs using the algorithm presented in Section 2.3.

• The PMQRFs are trained using the list of quantile levels from \(\mathcal{Q}_ {0.05}\).

• To benchmark our models, we compared their performance with that of the Quantile Regression Forest (QRF) method from [35].

An important aspect of the tree-based methods such as QRF, PMQRF, and CRPS-RF is that they rely on trees learning a single partition of the feature space for all quantile levels. This common partitioning has certain advantages, such as facilitating group conformal prediction (see Section 3.4 for more details) and ensuring that the estimated quantiles are non-crossing. However, for some data distributions, the use of a single partition for all quantiles might lead to suboptimal results. To assess the cost of using a common partition for all quantile levels, we compared the aforementioned methods with:

• We train 95 PQRTs where for each quantile level \(\tau\) in \(\mathcal{Q}_{0.05}\), \(5\) PQRTs are trained specifically using the pinball loss \(\ell_{\tau}\). At inference time, we query an estimate \(\hat{q}_{\tau}\) by aggregating the forest of \(5\) trees trained for the quantile level \(\tau\). We then rearrange the entire set of esimated quantiles following the method proposed in [41].

• An intermediate approach interpolating between the (PQRTs-5) and (PMQRF) methods, in which half of the trees (\(50\) trees) are trained using the sum of pinball losses \(\sum_{k=1}^{9}\ell_{0.05\times k}\) and the other half (\(50\) trees) with \(\sum_{k=10}^{19}\ell_{0.05\times k}\).

• A quantile regression model from the Scikit-learn library, where an independent linear model is trained for each quantile level. This allows for more flexible partitioning of the feature space, though it comes at the cost of losing the benefits of a shared partition.

Table 4 presents the comparative results obtained across all datasets. Overall, the findings reveal several consistent patterns regarding the relative performance and robustness of the evaluated methods.

3.3.0.1 Limitations of parametric models.

The SKLearn model achieves the best performance only on the Turbine dataset, where it significantly outperforms all other methods. However, it ranks as the worst model on every other dataset, often by a large margin. This behavior highlights the strong dependence of SKLearn on its underlying linearity assumption: when the true conditional quantiles can be well approximated by a linear model, SKLearn performs remarkably well, but its accuracy deteriorates substantially when this assumption is violated.

3.3.0.2 CRPS-RF and PMQRF as robust, assumption-free methods.

Across datasets, the CRPS-RF consistently provides the most accurate or at least competitive results among all tree-based approaches. Similarly, the PMQRF demonstrates strong and stable performance. These two methods share an important property: they are assumption-free and trained using interval-based scores (namely, the WIS or the CRPS; see [22]). This formulation allows them to model a wide variety of conditional distributions without relying on restrictive parametric forms, which explains their robustness across heterogeneous datasets.

3.3.0.3 Behavior of PQRTs-5 and related variants.

The PQRTs-5 method underperforms compared to both PMQRF and CRPS-RF, reflecting its limited robustness when the loss function focuses on a restricted set of quantiles. This behavior suggests that extreme quantiles are more difficult to estimate accurately in isolation, while jointly learning several quantiles tends to stabilize training and yield more reliable predictions. The intermediate variant—where half of the trees in the forest are trained on lower quantiles (\(\sum_{k=1}^{9}\ell_{0.05\times k}\)) and the other half on upper quantiles (\(\sum_{k=10}^{19}\ell_{0.05\times k}\))—does not lead to better results than PMQRF. Nevertheless, we believe that this strategy could be beneficial for datasets where the behaviors of the lower and upper tails differ substantially.

3.3.0.4 Assumption-based versus assumption-free models.

It is worth emphasizing that both SKLearn and QRF can be viewed as assumption-based models. SKLearn assumes linearity in the quantile relationship, while QRF relies on standard regression trees that implicitly assume a Gaussian conditional distribution through the use of variance as the impurity measure. This variance corresponds to the generalized entropy (negative log-likelihood) of a Gaussian model, which effectively optimizes a logarithmic score. In contrast, PMQRF and CRPS-RF are trained without such distributional assumptions, relying instead on interval-based scoring rules. These scores are known to be more robust than logarithmic scores [22], particularly when the true data distribution deviates from Gaussianity. Two main conclusions emerge from these experiments:

• Assumption-based models (e.g., SKLearn, QRF) can perform extremely well when their underlying assumptions hold—such as linearity or Gaussianity—but their performance deteriorates markedly when those assumptions are violated.

• PMQRF and CRPS-RF may not always be optimal for specific data distributions, yet they tend to yield competitive and stable results across diverse datasets. Their robustness stems from both their assumption-free nature and their reliance on interval-based scoring rules.

3.4 Adaptive group conformal prediction with CRPS random trees↩︎

In this section, we show how CRPS-RTs can be used to achieve group conformal prediction. Inspired by the simulated dataset from [42], we consider data \(\{(\boldsymbol{x}_i,y_i)\}_{i=1}^n\) where \[\label{eq:data}\forall i \in [n],\quad y_i = f(\boldsymbol{x}_i^{\top}\beta) + \epsilon_i \sqrt{1+\big(\boldsymbol{x}_i^{\top} \beta\big)^2},\tag{8}\] with \(\beta = [1\,, \,1] \in \mathbb{R}^2\), \(\boldsymbol{x}_i \sim \mathrm{Unif}([0,1]^2)\) are independent vectors, and \((\epsilon_i)_{i\in [n]}\) are independent standard normal random variables.

We train a single CRPS-RT with a training set of size \(n=2000\). We define four regions corresponding to the partition of the feature space obtained from the CRPS-RT at depth 2. We then evaluated different conformalized predictive inference methods, based on CRPS-RF or QRF with 50 trees, with both marginal and group-based conformalization considering the four previously defined groups. Our goal is to assess how well these methods achieve the desired coverage while maintaining narrow prediction intervals across the different regions.

We use the split conformal framework to conformalize the different methods. To do so, we consider a calibration set with \(1000\) data points and a test set of \(5000\) points. We conformalize the CRPS-RF method by considering the family of nested sets presented in [34], namely of the form \[[\hat{q}_t(\boldsymbol{x}),\hat{q}_{1-t}(\boldsymbol{x})],\] where \(\hat{q}_t(\boldsymbol{x})\) denotes the empirical quantile of level \(t\) obtained from our CRPS-RT by dropping the feature vector \(\boldsymbol{x} \in [0,1]^2\) from the root in the CRPS-RT and by taking the empirical quantile of level \(t\) from the data points in the leaf where we end up. We compare our method with a similar split conformal method obtained by training QRF and conformalizing using the CQR method from [33]. In both case, we target a marginal coverage of \(90\%=1-\alpha\).

Figure 5 illustrates the impact of different conformalization strategies on coverage and mean prediction interval width across the four regions identified by the CRPS-RT partitioning. Notably, the QRF-CQR method without group conformalization results in highly unbalanced coverage, failing to maintain the desired coverage level uniformly across the regions. This suggests that standard marginal conformalization struggles to adapt to the heterogeneity present in the learned feature space partition. In contrast, the CRPS-RF method with a distributional nested set and marginal calibration naturally provides better balance in coverage across the four regions. This improved calibration highlights the advantage of leveraging CRPS-based probabilistic forests for predictive inference.

Among all the methods considered, the group-conformalized CRPS-RF approach emerges as the most effective. It not only achieves the target coverage across all four regions but also maintains the narrowest marginal mean prediction interval width. This demonstrates that incorporating group-based conformalization within CRPS-RF leads to both improved calibration and sharper predictive intervals, making it the most reliable approach in this setting.

4 Conclusion↩︎

We have shown that it is possible to construct tree-based algorithms that minimize either the sum of \(M\) pinball losses or the CRPS, with computational complexity for one node split of, respectively, \(O(Mdn\log n)\) and \(O(dn\log n)\). In this setting, we also established that improved estimates of the population entropies can be obtained at no additional computational cost, by using LOO or Mallows-type corrections, leading to better generalization performance for the resulting trees and random forests. Empirically, PMQRFs and CRPS-RFs were shown to perform similarly to, or better than, existing approaches such as QRF.

These trees are particularly useful in the context of conformal prediction, since informative conformal intervals can be constructed from at least pairs of quantiles [33], or even from conditional CDFs [34]. Moreover, shallow trees produced by our algorithms can be used to define an interpretable partition of the feature space, thereby yielding natural groups on which group conformal prediction methods can be applied.

Disclosure statement: The authors report there are no competing interests to declare.

5 Background on Data Structures and Computational Tools↩︎

In this section, we provide a brief, self-contained introduction to the key data structures and computational tools used in our distributional regression tree algorithms.

We also provide concrete examples illustrating how our algorithms for learning PMQRT and CRPS-RT update their respective data structures — min–max heaps and Fenwick trees — throughout the iterations. We present the updates for selected steps to highlight the key mechanisms, while the full, step-by-step versions of these examples are available in the documentation of our package (section Algorithms unraveled) for interested readers.

5.1 Heaps↩︎

5.1.0.1 Presentation of the data structures.

A heap is a specialized tree-based data structure that allows efficient access to extreme values in a dataset. In a max-heap, each parent node is greater than or equal to its children, ensuring that the maximum element is always at the root. Conversely, in a min-heap, each parent node is less than or equal to its children, so the minimum element is at the root. Both max-heaps and min-heaps are useful in statistical algorithms when we need to track the largest or smallest values efficiently, such as maintaining quantile thresholds or top-k observations. A min-max heap generalizes these ideas to provide fast access to both the minimum and maximum elements simultaneously. In a min-max heap, nodes at even levels (with the root at level 0) are smaller than all their descendants, while nodes at odd levels are larger than all their descendants. This guarantees that the minimum element is at the root and the maximum element is at one of its children. A concrete example of min-max heap is presented in Figure 6.

5.1.0.2 Practical role in this work.

In our work, min–max heaps are used to easily maintain and update order statistics within each node during split evaluation, allowing us to efficiently perform a node splits for the generalized entropy associated to the WIS loss. Figure 7 gives a concrete example on the way the min-max heaps are updated at a given iteration to perform a node split in a PMQRT. We consider as loss function the sum of two pinball losses with quantile levels respectively \(0.3\) and \(0.7\). We assume that we need to identify the best split for a given feature \(j\) for which \(x_{1,j}\leq \dots\leq x_{7,j}\) with corresponding observations \(\boldsymbol{y}= [0,1,2,3,-1,-2,-3]\). In Figure 7, we illustrate the evolution of the three min-max heaps \(\mathcal{H}_1\), \(\mathcal{H}_2\) and \(\mathcal{H}_3\) from iteration \(5\) to \(6\) to compute the entropies associated to the prefix sequence \(\boldsymbol{y}^{(6)}\). Next to each heap, the blue value represents the sum of all elements within the corresponding heap.

5.2 Prefix and Suffix Sequences↩︎

5.2.0.1 Presentation of the data structures.

A prefix sequence is the subsequence starting from the first element up to a given position. Conversely, a suffix sequence starts from a given position and continues to the last element. Figure 8 illustrates the prefix and suffix sequences on a concrete example.

5.2.0.2 Practical role in this work.

Prefix and suffix sequences are convenient notations in our algorithms since, when evaluating candidate splits, we must compute the entropy for both the prefix sequence (samples that fall to the left of the split) and the suffix sequence (samples that fall to the right). More precisely, in a given leaf containing data points \(\{(x_i, y_i)\}_{i=1}^n\), for which we want to compute the information gains associated with all possible splits on feature \(j\), we assume without loss of generality that \[x_{1,j} \leq x_{2,j} \leq \dots \leq x_{n,j}.\] We then need to evaluate the entropies of the prefix sequence \((y_1, \dots, y_s)\) and the suffix sequence \((y_{s+1}, \dots, y_n)\) for all \(s \in [n]\) in order to compute the information gain corresponding to feature \(j\) and split position \(s\) in the current leaf.

5.3 Weight-Balanced Binary Trees↩︎

5.3.0.1 Presentation of the data structures.

A weight-balanced binary tree is a type of binary search tree in which the sizes of left and right subtrees are maintained roughly equal to ensure balanced height. This guarantees that insertion, deletion, and query operations can be performed in logarithmic time.

5.3.0.2 Practical role in this work.

In our work, weight-balanced trees are used in the construction of CRPS-RTs. When computing the information gains for a feature \(j\) within a given leaf, we iteratively evaluate the gains for all possible splits \(s \in [0, n]\) by computing the entropies of the corresponding prefix and suffix sequences, \((y_1, \dots, y_s)\) and \((y_{s+1}, \dots, y_n)\) (where we consider without loss of generality that \(x_{1,j} \leq x_{2,j} \leq \dots \leq x_{n,j}\)).

To perform these iterations efficiently, we need to determine the insertion rank of the element \(y_{s+1}\) within the sorted version of the list \(\mathbf{y}^{(s)} = (y_1, \dots, y_s)\). This rank is crucial for efficiently updating the Fenwick tree, which we use to compute the information gains associated with all possible splits for a given feature \(j\).

A naive approach would involve searching for the insertion rank of \(y_{s+1}\) in the sorted version of \(\mathbf{y}^{(s)}\) at each iteration, resulting in prohibitive computational cost. To mitigate this, given a feature \(j\), we first iteratively insert the elements \(y_1, y_2, \dots, y_n\) into a Weight-Balanced Binary Tree. This allows us to obtain, in \(\mathcal{O}(n \log n)\) time, the rank of insertion of \(y_{s+1}\) in the sorted version of \(\mathbf{y}^{(s)}\) for any \(s \in [n]\).

These precomputed insertion ranks, denoted by \(\{r_s\}_{s=1}^n\), are then used to update the Fenwick tree efficiently during the computation of information gains. In the next subsection, we provide a concrete example illustrating how the Fenwick tree is updated using these precomputed ranks.

Figure 9 shows how a Weight-Balanced Binary Tree can be used to determine the rank of a new element in a sorted list.

5.4 Fenwick Trees↩︎

5.4.0.1 Presentation of the data structures.

A Fenwick tree, also called a binary indexed tree, is a data structure designed to efficiently compute and update prefix sums over a sequence of numbers. Unlike a simple prefix sum array, a Fenwick tree allows efficient dynamic updates: if an element in the sequence changes, all relevant prefix sums can be updated in logarithmic time, without recomputing sums for the entire array.

Formally, given a sequence \((x_1, x_2, \dots, x_n)\), the Fenwick tree stores partial sums that allow two operations in \(\mathcal{O}(\log n)\) time:

1. Query: compute the sum of the first \(i\) elements, \(\sum_{k=1}^i x_k\).

2. Update: add a value to a particular element \(x_i\) and update all affected prefix sums.

We refer to Figure 1 for a visualization of the structure of a Fenwick tree.

5.4.0.2 Practical role in this work.

We provide a concrete example to illustrate how our algorithm to learn CRPS-RT update the Fenwick tree along the iterations. Figure 10 illustrates the updates of the Fenwick tree to perform a node split in a CRPS-RT. Similarly to Figure 7, we assume that we need to identify the best split for a given feature \(j\) for which \(x_{1,j}\leq \dots\leq x_{6,j}\) with corresponding observations \(\boldsymbol{y}= [2,1,3,-1,-3,-2]\).

6 Study of different aggregation strategies↩︎

We evaluate the performance of the two aggregation methods described in the introduction of Section 3, namely distributional bagging (DB) and quantile bagging (QB). For each dataset introduced in Section 3, we randomly subsample 1,000 points to train a CRPS-RF consisting of 100 trees. The maximum depth was fixed at 10, the \(min\_samples\_split\) parameter was set to 30, and LOO was used to compute the information gains. We then query the \(5\%\) and \(95\%\) quantiles on a test set of size 1,000 and compute both the marginal width of the predictive intervals \([\hat{q}_{0.05},\hat{q}_{0.95}]\) and the corresponding marginal coverage. This procedure is repeated over 100 different subsamples of the original dataset, allowing us to report the mean and standard deviation of the marginal width and coverage across these scenarios. Table 5 summarizes the results, showing that the QB aggregation scheme achieves a marginal coverage close to the nominal \(90\%\), while the DB scheme tends to be more conservative, resulting in wider predictive intervals and higher coverage than expected. For this reason, the QB aggregation scheme was preferred in the experiments presented in our paper.

7 Package description↩︎

Table ¿tbl:tab:config-options? summarizes the main configuration options available in the DisTreebution package. The parameters are organized into three conceptual layers reflecting the structure of the modeling framework. The first layer defines the model type, allowing users to choose between CRPS regression trees (CRPS-RTs) and pinball-loss-based quantile regression trees (PMQRTs). The second layer specifies the tree training settings, which control how individual trees are constructed and aggregated within the ensemble. These include standard hyperparameters such as tree depth and minimum samples per split, as well as more advanced options like the use of a leave-one-out estimate of the information gain for data-driven split selection, and the strategy for aggregating leaf predictions. The final layer relates to the conformalization settings, which govern the calibration of predictive intervals. Users can choose between a distributional or conformalized quantile regression (CQR) approach, optionally applying conformalization at the group level. Together, these parameters offer a flexible yet principled way to configure DisTreebution for a range of uncertainty quantification tasks.

7.0.0.1 Quantile-level configuration for PMQRTs.

When using PMQRTs, the dictionary \(treeID2quantiles\_train\) specifies the set of quantile levels on which each tree is trained. More concretely, setting \(\left\{i:[\beta,1-\beta]\, \text{for} \, i \in [nTrees]\right\}\) trains all trees in the forest using the combined pinball loss \(\ell_{\beta}+\ell_{1-\beta}\). Alternatively, defining \(\left\{i:[\beta], \text{for} , i \in [nTrees/2]\right\}\cup\left\{i:[1-\beta]\, \text{for} \, i \in \{nTrees/2,\dots,nTrees\}\right\}\) results in half of the trees being optimized with \(\ell_{\beta}\) and the other half with \(\ell_{1-\beta}\). By default, each tree can be queried only for the quantile levels it has been trained on, although this behavior can be manually adjusted by the user. We refer readers to the package documentation for additional configuration options. This design provides flexibility, allowing PMQRFs to distribute modeling effort across different quantile levels.

8 List of notations↩︎

9 Algorithms pseudo-codes for PMQRTs and CRPS-RTs↩︎

9.1 Computations of the entropies for all splits for multiple quantiles↩︎

We consider \(M\) quantile levels: \(\tau_1<\dots <\tau_M\). In this section, we provide the detailed algorithm to train regression trees to minimize the empirical entropy associated with a sum of the pinball losses with quantile levels \(\tau_1,\dots,\tau_M\). By convention, querying the minimum or the maximum of an empty heap returns a \(\text{None}\) value. Moreover, each min-max heap data structure keeps track of the total sum of the current elements within the heap.

Algorithm 11 describes the procedure to compute the information gain induced by all the possible splits at a given node.

9.2 Computations of the entropies for all splits with the CRPS↩︎

Theorem 1 ensures that one can compute the information gain at a given node for the possible split of the data points in that node by:

• computing the insertion ranks \((r_s)_{s\in[n]}\) before running the iterative scheme. This can be done in time complexity \(\mathcal{O}(n\log n)\) using a weight-balanced binary tree.

• updating along the iterative process a Fenwick tree \(\mathcal{F}\) that will allow us at iteration \(s\) to query \(S^{(s)}_{r_s}, S^{(s)}\) and \(W^{(s)}\) in \(\mathcal{O}(\log n)\) time.

Algorithm 16 describes the procedure to compute the information gain induced by all the possible splits at a given node.

10 Proofs↩︎

10.1 Proof of Theorem 1↩︎

The CRPS is defined by \[\mathrm{CRPS}(F,y) = \int_{-\infty}^{+\infty} \big( F(s) - \mathbb{1}_{s\geq y}\big)^2 ds.\] The associated entropy is \[\min_F \mathbb{E}[ \mathrm{CRPS}(F,Y)].\] Considering its empirical counterpart, the minimizer is the empirical CDF. \[\hat{F}:s\mapsto \frac{1}{n} \sum_{i=1}^n \mathbb{1}_{s\geq y_i}.\] We deduce that the empirical entropy is given by \[\begin{align} &\frac{1}{n} \sum_{i=1}^n\mathrm{CRPS}(\widehat F,y_i)\\ & = \frac{1}{n} \sum_{i=1}^n \int_{-\infty}^{+\infty} \big( \hat{F}(s) - \mathbb{1}_{s\geq y_i}\big)^2 ds\\ &= \frac{1}{n} \sum_{i=1}^n \int_{-\infty}^{+\infty} \big( \frac{1}{n} \sum_{j=1}^n \mathbb{1}_{s\geq y_j} - \mathbb{1}_{s\geq y_i}\big)^2 ds\\ &= \frac{1}{n^3} \sum_{i,j,k} \int_{-\infty}^{+\infty} \big( \mathbb{1}_{s\geq y_j} - \mathbb{1}_{s\geq y_i}\big)\big( \mathbb{1}_{s\geq y_k} - \mathbb{1}_{s\geq y_i}\big) ds\\ &= \frac{2}{n^3} \sum_{i<j<k} \int_{-\infty}^{+\infty} \big( \mathbb{1}_{s\geq y_{(j)}} - \mathbb{1}_{s\geq y_{(i)}}\big)\big( \mathbb{1}_{s\geq y_{(k)}} - \mathbb{1}_{s\geq y_{(i)}}\big) ds\\ &\qquad + \frac{2}{n^3} \sum_{j<i<k} \int_{-\infty}^{+\infty} \big( \mathbb{1}_{s\geq y_{(j)}} - \mathbb{1}_{s\geq y_{(i)}}\big)\big( \mathbb{1}_{s\geq y_{(k)}} - \mathbb{1}_{s\geq y_{(i)}}\big) ds\\ &\qquad + \frac{2}{n^3} \sum_{j<k<i} \int_{-\infty}^{+\infty} \big( \mathbb{1}_{s\geq y_{(j)}} - \mathbb{1}_{s\geq y_{(i)}}\big)\big( \mathbb{1}_{s\geq y_{(k)}} - \mathbb{1}_{s\geq y_{(i)}}\big) ds\\ &\qquad + \frac{2}{n^3} \sum_{j<i} \int_{-\infty}^{+\infty} \big( \mathbb{1}_{s\geq y_{(j)}} - \mathbb{1}_{s\geq y_{(i)}}\big)^2 ds\\ &= \frac{2}{n^3} \sum_{i<j<k} (y_{(j)}-y_{(i)}) \quad + 0\quad + \frac{2}{n^3} \sum_{j<k<i} (y_{(i)}-y_{(k)})+ \frac{2}{n^3} \sum_{j<i} (y_{(i)}-y_{(j)}) \\ &= \frac{2}{n^3} \sum_{i=1}^{n-2} \sum_{j=i+1}^{n-1} (n-j)(y_{(j)}-y_{(i)}) - \frac{2}{n^3} \sum_{j=1}^{n-2} \sum_{k=j+1}^{n-1} \sum_{i=k+1}^{n}(y_{(k)}-y_{(i)})+\frac{2}{n^3} \sum_{i<j} (y_{(j)}-y_{(i)})\\ &= \frac{2}{n^3} \sum_{i=1}^{n-1} \sum_{j=i+1}^{n} (n-j+1)(y_{(j)}-y_{(i)}) - \frac{2}{n^3} \sum_{j=1}^{n-2} \sum_{k=j+1}^{n-1} \sum_{i=k+1}^{n}(y_{(k)}-y_{(i)})\\ &= - \frac{2}{n^3} \sum_{i=1}^{n-1} [n(n-i) -\frac{n(n-1)}{2}+\frac{i(i-1)}{2}] y_{(i)} + \frac{2}{n^3}\sum_{j=2}^{n}(n-j+1)(j-1)y_{(j)}\\ &\qquad -\frac{2}{n^3} \sum_{i=3}^{n} \sum_{k=2}^{i-1} \sum_{j=1}^{k-1} (y_{(k)}-y_{(i)})\\ &= - \frac{2}{n^3} \sum_{i=1}^{n-1} \frac{1}{2}(n-i)(n-i+1) y_{(i)} + \frac{2}{n^3}\sum_{j=2}^{n}(n-j+1)(j-1)y_{(j)}-\frac{2}{n^3} \sum_{i=3}^{n} \sum_{k=2}^{i-1} (k-1) (y_{(k)}-y_{(i)})\\ &= - \frac{1}{n^3} \sum_{i=1}^{n-1} (n-i)(n-i+1) y_{(i)} + \frac{2}{n^3}\sum_{j=2}^{n}(n-j+1)(j-1)y_{(j)}\\ &\qquad +\frac{2}{n^3} \sum_{i=3}^{n} \frac{(i-2)(i-1)}{2} y_{(i)} - \frac{2}{n^3} \sum_{k=2}^{n-1} (n-k)(k-1)y_{(k)} \\ &= - \frac{1}{n^3} \sum_{i=1}^{n-1} (n-i)(n-i+1) y_{(i)} +\frac{1}{n^3} \sum_{i=3}^{n} (i-2)(i-1) y_{(i)} +\frac{2}{n^3} \sum_{k=2}^{n} (k-1)y_{(k)}\\ &= - \frac{1}{n^3} \sum_{i=1}^{n} (n-i)(n-i+1) y_{(i)} +\frac{1}{n^3} \sum_{i=3}^{n} (i-2)(i-1) y_{(i)} +\frac{2}{n^3} \sum_{k=1}^{n} (k-1)y_{(k)}\\ &= \frac{1}{n^3} \sum_{i=1}^{n} \left[(i-2)(i-1)+2(i-1)-(n-i)(n-i+1) \right] y_{(i)} \\ &= \frac{1}{n^3} \sum_{i=1}^{n} \left[i(i-1)-(n-i)(n-i+1) \right] y_{(i)} \\ &= \frac{1}{n^3} \sum_{i=1}^{n} (i-1)i (y_{(i)}-y_{(n-i+1)}) . \end{align}\] We denote by \(\boldsymbol{y}^{(n)}\) the vector of size \(n\) corresponding to the \(n\) first entries of the vector of observations \(\boldsymbol{y}\) (where we consider without loss of generality that elements in \(\boldsymbol{y}\) have been ordered considering the permutation sorting in increasing order \(\big( x_{i,k}\big)_{i: \boldsymbol{x}_i \in R_j}\) for some given node \(k\) and some feature \(j\)). From the previous computations we have \[\frac{1}{n} \sum_{i=1}^n\mathrm{CRPS}(\widehat F,y_i) = \frac{1}{n^3} \left(h^{(n)}_{\uparrow} -h^{(n)}_{\downarrow} \right),\] where \[\begin{align} h^{(n)}_{\uparrow}&:=\sum_{i=1}^{n} (i-1)iy^{(n)}_{(i)}\quad \text{and} \quad h^{(n)}_{\downarrow}:=\sum_{i=1}^{n} (i-1)i y^{(n)}_{(n-i+1)}. \end{align}\] In the following, we derive iterative formula for both \(h^{(n)}_{\uparrow}\) and \(h^{(n)}_{\downarrow}\).

10.1.0.1 Iterative formula for \(h^{(n)}_{\uparrow}\)

.

\[\begin{align} h^{(n+1)}_{\uparrow}-h^{(n)} _{\uparrow}&=\sum_{i=1}^{n+1} (i-1)iy^{(n+1)}_{(i)} -\sum_{i=1}^{n} (i-1)iy^{(n)}_{(i)} \\ &=\sum_{i=1}^{n+1} (i-1)i y^{(n+1)}_{(i)} -\sum_{i=1}^{r_{n+1}-1} (i-1)iy^{(n+1)}_{(i)} -\sum_{i=r_{n+1}}^{n} (i-1)i y^{(n+1)}_{(i+1)} \\ &=\sum_{i=1}^{n+1} (i-1)iy^{(n+1)}_{(i)} -\sum_{i=1}^{r_{n+1}-1} (i-1)i y^{(n+1)}_{(i)} -\sum_{i=r_{n+1}+1}^{n+1} (i-2)(i-1) y^{(n+1)}_{(i)} \\ &=(r_{n+1}-1) r_{n+1} y^{(n+1)}_{(r_{n+1})} + 2\sum_{i=r_{n+1}+1}^{n+1} (i-1) y^{(n+1)}_{(i)}\\ &=[(r_{n+1}-1) r_{n+1} - 2(r_{n+1}-1)] y^{(n+1)}_{(r_{n+1})} + 2 W^{(n+1)}\\ &-2\sum_{i=1}^{r_{n+1}-1} i y^{(n+1)}_{(i)}-2S^{(n+1)}+2\sum_{i=1}^{r_{n+1}-1} y^{(n+1)}_{(i)}\\ &=[(r_{n+1}-1) r_{n+1} - 2(r_{n+1}-1)] y^{(n+1)}_{(r_{n+1})} + 2 W^{(n+1)}\\ &-2\sum_{i=1}^{r_{n+1}-1} i y^{(n+1)}_{(i)}-2S^{(n+1)}+2S^{(n+1)}_{r_{n+1}-1}. \end{align}\]

10.1.0.2 Iterative formula for \(h^{(n)}_{\downarrow}\)

.

\[\begin{align} h^{(n+1)}_{\downarrow}-h^{(n)}_{\downarrow} &=\sum_{i=1}^{n+1} (i-1)i y^{(n+1)}_{(n-i+2)} -\sum_{i=1}^{n} (i-1)i y^{(n)}_{(n-i+1)} \\ &=\sum_{i=1}^{n+1} (n-i+1)(n-i+2) y^{(n+1)}_{(i)} -\sum_{i=1}^{n} (n-i)(n-i+1) y^{(n)}_{(i)} \\ &=\sum_{i=1}^{n+1} (n-i+1)(n-i+2)y^{(n+1)}_{(i)} -\sum_{i=1}^{r_{n+1}-1} (n-i)(n-i+1) y^{(n+1)}_{(i)}\\ &\quad -\sum_{i=r_{n+1}}^{n} (n-i)(n-i+1) y^{(n+1)}_{(i+1)} \\ &=\sum_{i=1}^{n+1} (n-i+1)(n-i+2) y^{(n+1)}_{(i)} -\sum_{i=1}^{r_{n+1}-1} (n-i)(n-i+1) y^{(n+1)}_{(i)}\\ &\quad -\sum_{i=r_{n+1}+1}^{n+1} (n-i+1)(n-i+2)y^{(n+1)}_{(i)} \\ &= 2\sum_{i=1}^{r_{n+1}-1} (n-i+1)y^{(n+1)}_{(i)} + (n-r_{n+1}+1) (n-r_{n+1}+2) y^{(n+1)}_{(r_{n+1})} \\ &= 2(n+1) S^{(n+1)}_{r_{n+1}-1}-2 \sum_{i=1}^{r_{n+1}-1} i y^{(n+1)}_{(i)} + (n-r_{n+1}+1) (n-r_{n+1}+2) y^{(n+1)}_{(r_{n+1})} . \end{align}\]

10.1.0.3 Conclusion.

Using the previous computations, we obtain that \[\begin{align} h^{(n+1)}_{\uparrow}-h^{(n+1)}_{\downarrow} &=h^{(n)}_{\uparrow}-h^{(n)}_{\downarrow} + n\big[2r_{n+1} -3-n\big] y^{(n+1)}_{(r_{n+1})} + 2 W^{(n+1)}-2S^{(n+1)}-2nS^{(n+1)}_{r_{n+1}-1}. \end{align}\] Noticing that \(S^{(n+1)}_{r_{n+1}-1}=S^{(n+1)}_{r_{n+1}}-y_{n+1}\) concludes the proof.

10.2 Proof of Lemma 2↩︎

10.2.0.1 Querying \(S^{(s)}\) in \(\mathcal{O}(\log n)\) time.

Let us recall that by construction, \(\mathcal{F}\) is at iteration \(s\) a Fenwick tree for the sequence \(\big( \mathrm 1_{\sigma(i)\leq s} y_{\sigma(i)}\big)_{i\in[n]}\). One can thus easily see that getting the sum of all elements in \(\mathcal{F}\) gives \(S^{(s)}=\sum_{i=1}^s y^{(s)}_{(i)}=\sum_{i=1}^s y_i\). Indeed, up to ordering of the elements, the list \(\big( \mathrm 1_{\sigma(i)\leq s} y_{\sigma(i)}\big)_{i\in[n]}\) is equal to the one obtained by applying the permutation \(\sigma^{-1}\) namely to \(\big( \mathrm 1_{i\leq s} y_{i}\big)_{i\in[n]}\) whose sum is equal to \(S^{(s)}\). Hence, querying the sum of all elements of \(\mathcal{F}\) at iteration \(s\) gives \(S^{(s)}\) and this operation requires \(\mathcal{O}(\log n)\) time since \(\mathcal{F}\) is a Fenwick tree.

10.2.0.2 Querying \(S^{(s)}_{r_s-1}\) in \(\mathcal{O}(\log n)\) time.

Let us now prove that querying the sum of all elements in \(\mathcal{F}\) from position \(1\) to \(\sigma^{-1}(s)\) at iteration \(s\) is equal to \(S^{(s)}_{r_s}\). This will prove that \(S^{(s)}_{r_s}\) can be obtained in \(\mathcal{O}(\log n)\) time since \(\mathcal{F}\) is a Fenwick tree. One can then get \(S^{(s)}_{r_s-1}\) using that \(y^{(s)}_{(r_s)}=y_s\) and thus that \[S^{(s)}_{r_s-1}=S^{(s)}_{r_s} -y_s.\] Indeed, \(S^{(s)}_{r_s-1}=\sum_{i=1}^{r_s-1}y^{(s)}_{(i)}\) where \(r_s\) is the position of insertion of \(y_s\) in the non-decreasing sorted version of the list \(\mathbf{y}^{(s-1)}\).

Since \(y^{(s)}_{(r_s)}=y_s\), the list \((y^{(s)}_{(1)}, \dots,y^{(s)}_{(r_s)})\) exactly corresponds to all the elements in \((y_{1}, \dots, y_{s})\) smaller than \(y_s\) i.e. to the list of all \(y_k\) for \(k\leq s\) such that \(\sigma^{-1}(k)\leq \sigma^{-1}(s)\). We deduce that \[\label{eq:lemma2} S^{(s)}_{r_s}=\sum_{i=1}^{r_s} y^{(s)}_{(i)} = \sum_{i=1}^s \mathrm 1_{\sigma^{-1}(i)\leq \sigma^{-1}(s)} y_i=\sum_{j\in \{ \sigma^{-1}(k) \mid k\in[s]\}} \mathrm 1_{j\leq \sigma^{-1}(s)} y_{\sigma(j)}=\sum_{j=1}^{\sigma^{-1}(s)} \mathrm 1_{\sigma(j)\leq s} y_{\sigma(j)},\tag{9}\] where in the penultimate equality we made the change of variable \(i=\sigma(j)\) and where in the last equality we used that \[\{ \sigma^{-1}(k) \mid k\in[s], \sigma^{-1}(k)\leq \sigma^{-1}(s)\}= \{ j\mid j\in[\sigma^{-1}(s)], \sigma(j)\leq s\}.\] Since \(\mathcal{F}\) is at iteration \(s\) a Fenwick tree for the sequence \(\big( \mathrm 1_{\sigma(i)\leq s} y_{\sigma(i)}\big)_{i\in[n]}\), we deduce from Eq.@eq:eq:lemma2 that querying the sum of elements from position 1 to \(\sigma^{-1}(s)\) in \(\mathcal{F}\) gives \(S^{(s)}_{r_s}\).

As a remark, one can notice that a direct consequence of this result is that querying the sum of all elements from position 1 to \(\sigma^{-1}(s+1)\) at iteration \(s\) from \(\mathcal{F}\) gives \(S^{(s)}_{r_{s+1}-1}\).

10.2.0.3 Getting \(W^{(s+1)}\) in constant time given \(W^{(s)}\), \(S^{(s)}\) and \(S^{(s+1)}_{r_{s+1}-1}\).

The total rank-weighted sum \(W^{(s)}\) can be efficiently tracked along the iterative process using the formula \[\label{eq:key95update95CRPS95apdx}W^{(s+1)} = W^{(s)} +r_{s+1}y_{s+1} +\sum_{i=r_{s+1}}^s y^{(s)}_{(i)}=W^{(s)} +r_{s+1}y_{s+1} +S^{(s)}-S^{(s)}_{r_{s+1}-1},\tag{10}\] where \[S^{(s)}_{r_{s+1}-1} = S^{(s+1)}_{r_{s+1}-1},\] since \(r_{s+1}\) is the position of insertion of \(y_{s+1}\) in the list \(( y^{(s)}_{(1)}, \dots, y^{(s)}_{(s)})\).

11 Unbiased estimate of information gains↩︎

11.1 Reminder: Risk estimation for the quadratic loss↩︎

Consider a dataset \(\mathcal{D}_n={(x_i,y_i)}_{i=1}^n\) with squared loss \(\ell(\hat{y},y)=(y-\hat{y})^2\). Let \(\hat{f}_{1:n}\) denote the estimator trained on the full dataset and \(\hat{f}^{-i}_{1:n}\) the estimator trained on the sample with observation \(i\) removed.

11.1.0.1 Mallows \(C_p\) and optimism correction for the quadratic loss.

A natural estimator of the prediction risk of \(\hat{f}_{1:n}\) is the empirical risk \[\widehat R_{1:n}(\hat{f}_{1:n})=\frac{1}{n}\sum_{i=1}^n (y_i-\hat{f}_{1:n}(x_i))^2 .\] However, this quantity is optimistically biased because the same data are used both to fit the model and evaluate its performance. In expectation, \[\mathbb{E}[\widehat R_{1:n}(\hat{f}_{1:n})] \le \mathbb{E}[R(\hat{f}_{1:n})],\]

where \(R(\hat{f}_{1:n})=\mathbb{E}_{(X,Y)}[(Y-\hat{f}_{1:n}(X))^2]\) denotes the population risk.

A general way to correct this bias is to introduce the optimism term \[\mathrm{pen}_{n} := \mathbb{E}\left[ R(\hat{f}_{1:n})-\widehat R_{1:n}(\hat{f}_{1:n}) \right].\] This quantity measures the expected gap between the true prediction risk and the training error. If an unbiased estimate \(\widehat{\mathrm{pen}}_{1:n}\) of \(\mathrm{pen}_n\) is available, one can construct an unbiased estimator of the prediction risk: \[\widehat R_{1:n}(\hat{f}_{1:n}) + \widehat{\mathrm{pen}}_{1:n} .\] Mallows \(C_p\) provides such an estimator in the case of squared loss under the classical homoscedastic regression model. For linear estimators with effective degrees of freedom \(d\) and noise variance \(\sigma^2\), the optimism admits the simple expression \[\mathrm{pen}_n = \frac{2\sigma^2}{n} d .\] This leads to the well-known \(C_p\) criterion \[C_p = \widehat R_{1:n}(\hat{f}_{1:n}) + \frac{2\sigma^2}{n} d .\]

11.1.0.2 Leave-One-Out (LOO) estimate.

The leave-one-out estimate of the prediction risk is defined as \[\widehat R^{\text{LOO}}_{1:n} := \frac{1}{n}\sum_{i=1}^{n} (y_i - \hat{f}^{-i}_{1:n}(x_i))^2 .\] Each term evaluates the prediction error of a model trained without observation \(i\) on that left-out observation. \(\widehat R^{\text{LOO}}\) is an unbiased estimate of the expected population risk of the learning algorithm trained on \(n-1\) observations, \[\mathbb{E}\big[ R(\hat{f}_{1:n-1}) \big],\] where \(R(\hat{f}_{1:n-1})=\mathbb{E}_{(X,Y)}[(Y-\hat{f}_{1:n-1}(X))^2]\) denotes the population risk of the estimator \(\hat{f}_{1:n-1}\) trained with \(n-1\) observations.

Importantly, this means that LOO does not provide an unbiased estimate of the risk of \(\hat{f}_{1:n}\) (the estimator trained on the full dataset). The two risks are typically very close when \(n\) is large, but they correspond to different training sample sizes.

11.2 Unbiased estimate of information gain in regression trees.↩︎

Several authors have recognized that the standard CART algorithm can produce biased estimates of the information gain because splits are evaluated on the same data used to estimate node parameters. In particular, S. Athey et al.² already noted that one should use strategies, such as Mallows-style corrections, to adjust for the additional sampling variance in small leaves. Implementing such a correction in practice requires efficiently estimating the corresponding optimism term for each possible split and feature.

In the case of the quadratic loss, the standard Mallows \(C_p\) correction terms does not induce additional computational complexity since the estimate of \(\mathrm{pen}_{n}\) given by \(\widehat{\mathrm{pen}}_{1:n}:=\frac{2\widehat \sigma^2}{n}d\) where \(\widehat \sigma^2\) is the standard unbiased estimate of the noise variance. Similarly, the LOO estimate of the information gain for the quadratic loss does not cause any computational overhead.

The choice between Mallows-style correction and LOO should be guided by their respective theoretical properties and computational costs. Suppose we are evaluating a node split for a leaf containing \(n\) observations. We denote by \(\hat{f}_{k:j}^{\mathrm{ERM}}\) the empirical risk minimizer for the loss \(\ell\) (typically the CRPS of the WIS) computing from the data points \(k\) to \(j\) (with \(1\leq k< j\leq n\)). With the notations introduced in the main paper, we have \(\widehat R_{k:j}(\hat{f}^{\mathrm{ERM}}_{k:j})=H_{\ell}(\mathbf{y}_{k:j})\). Rather than relying on a standard empirical estimate of information gain, we can enhance the decision procedure using one of the following approaches:

• Mallows-style correction provides an unbiased estimate of the expected information gain for the model trained on the same number of data points as the one the would be used at inference time. Namely, for a candidate split leading to the children nodes with data \((y_1,\dots,y_k)\) and \((y_{k+1},\dots,y_n)\), the standard information gain \[\begin{align} \widehat { \mathrm{IG}}_k&:=\widehat R_{1:n}(\hat{f}_{1:n}^{\mathrm{ERM}})-\left(\widehat R_{1:k}(\hat{f}_{1:k}^{\mathrm{ERM}}) + \widehat R_{k:n}(\hat{f}_{k:n}^{\mathrm{ERM}})\right)\\ &= H_{\ell}(\mathbf{y}_{1:n}) - (H_{\ell}(\mathbf{y}_{1:k}) + H_{\ell}(\mathbf{y}_{k:n})) \end{align}\] can be corrected by adding the appropriate optimism terms, namely one adds \(\widehat{\mathrm{pen}}_{1:n}\) to \(H_{\ell}(\mathbf{y}_{1:n})\), \(\widehat{\mathrm{pen}}_{1:k}\) to \(H_{\ell}(\mathbf{y}_{1:k})\) and \(\widehat{\mathrm{pen}}_{k:n}\) to \(H_{\ell}(\mathbf{y}_{k:n})\). It holds \[\mathbb{E}\left[\widehat{\mathrm{IG}}_k+\widehat{\mathrm{pen}}_{1:n}-\left(\widehat{\mathrm{pen}}_{1:k}+\widehat{\mathrm{pen}}_{k:n}\right) \right]=\mathbb{E}\left[R(\hat{f}_{n}^{\mathrm{ERM}}) - \left(R(\hat{f}_{k}^{\mathrm{ERM}}) + R(\hat{f}_{n-k}^{\mathrm{ERM}})\right)\right].\label{eq:target95IG}\tag{11}\]

• LOO estimates of the information gains yield unbiased estimates of the risk for models trained one less data point in each child leaf. Stated otherwise, the expected value of the LOO estimate of the information gain is given by: \[\mathbb{E}\left[ R(\hat{f}_{n-1}^{\mathrm{ERM}}) - \left( R(\hat{f}_{k-1}^{\mathrm{ERM}}) + R(\hat{f}_{n-k-1}^{\mathrm{ERM}})\right)\right]\] Consequently, LOO provides a biased estimate for the quantity we are targeting, namely the one given in Eq.@eq:eq:target95IG .

  However, this bias tends to penalize splits that create very small leaves, because removing a single observation from a small leaf has a larger effect on the LOO error. As a result, LOO-based selection often favors splits that produce larger, more balanced leaves, leading to more robust trees and potentially better generalization performance. Our numerical experiments, presented in the main paper, indicate that the results obtained using either LOO or Mallows corrections are very consistent.

In the rest of this section, we show that LOO estimates of the information gains are available for both the sum of pinball losses and the CRPS loss without additional computational cost. We also show that a Mallows-style correction can be derived for the CRPS loss at essentially no extra cost. Deriving a similar Mallows-type correction for the sum of pinball losses appears to be more challenging. We therefore do not provide such a result here, and to the best of our knowledge we are not aware of any existing work that derives such a correction.

11.2.1 LOO for PMQRT↩︎

We propose using a leave-one-out method to obtain an unbiased estimate of the information gain used to build the PQRT or PMQRT. In the following, we show that such LOO estimate of the entropy associated with a pinball loss \(\ell_{\tau}\) with data \(\boldsymbol{y}=(y_1,\dots,y_n)\) can be computed efficiently.

Let us denote by \(\mathcal{H}_n\) the entropy of the sequence \(\boldsymbol{y}\) and by \(\mathcal{H}_{n,LOO}\) the leave-one-out estimate of the entropy of \(\boldsymbol{y}\) for the pinball loss \(\ell_{\tau}\), namely \[\mathcal{H}_n = \frac{1}{n}\sum_{i=1}^n \ell_{\tau}(y_i, \hat{q}_{\tau}),\quad \mathcal{H}_{n,LOO} = \frac{1}{n}\sum_{i=1}^n \ell_{\tau}(y_i, \hat{q}_{\tau}^{-i}),\] where \(\hat{q}_{\tau}\) (resp. \(\hat{q}_{\tau}^{-i}\)) is the empirical quantile of order \(\tau\) of the sequence \(\boldsymbol{y}\) (resp. \(\boldsymbol{y}^{-i}:=(y_1,\dots, y_{i-1},y_{i+1},\dots, y_n)\)).

Let us consider a permutation \(\sigma\) sorting the sequence \(\boldsymbol{y}\) by increasing order. We further denote \(r_i:=\sigma^{-1}(i)\) which corresponds to the position in the sorted list of the element \(y_i\). Defining \(r^-:=\ceil{\tau (n-1)}\) and \(r^*:=\ceil{n\tau}\), we can distinguish two possible cases:

• If \(r^-=r^*\):

  • If \(r_i>r^-\): \(\ell_{\tau}(y_i, \hat{q}^{-i}_{\tau})=\ell_{\tau}(y_i, \hat{q}_{\tau}) = \tau (y_i - y_{(r^*)}).\)

  • Else: \(\ell_{\tau}(y_i, \hat{q}^{-i}_{\tau})=(1-\tau)(y_{(r^*+1)}-y_i).\)

  Hence: \[\begin{align} \mathcal{H}_{n,LOO} &= \mathcal{H}_n - \frac{(1-\tau)}{n} \sum_{i \mid r_i\leq r^*} (y_{(r^*)} -y_i) + \frac{(1-\tau)}{n} \sum_{i \mid r_i\leq r^*} (y_{(r^*+1)} -y_i) \\ &= \mathcal{H}_n + \frac{(1-\tau) r^*}{n} (y_{(r^*+1)} - y_{(r^*)} ). \end{align}\]

• Else (i.e. if \(r^-+1=r^*\)):

  • If \(r_i>r^-\): \(\ell_{\tau}(y_i, \hat{q}^{-i}_{\tau})=\ell_{\tau}(y_i, y_{(r^*-1)}) = \tau (y_i - y_{(r^*-1)}).\)

  • Else: \(\ell_{\tau}(y_i, \hat{q}^{-i}_{\tau})=(1-\tau)(y_{(r^*)}-y_i).\)

  Hence: \[\begin{align} \mathcal{H}_{n,LOO} &= \mathcal{H}_n - \frac{\tau}{n} \sum_{i \mid r_i\leq r^*} (y_i-y_{(r^*)} ) + \frac{\tau}{n} \sum_{i \mid r_i\leq r^*} (y_i-y_{(r^*-1)} )\\ &= \mathcal{H}_n + \frac{\tau (n-r^*+1)}{n} (y_{(r^*)} - y_{(r^*-1)} ). \end{align}\]

All quantities \(y_{(r^*-1)}, y_{(r^*)}\) and \(y_{(r^*+1)}\) can be queried in \(\mathcal{O}(\log(n))\) time using the min-max heap structures used to train PMQRT (or PQRT). Overall, a node split using or not these LOO estimates of entropies is performed in \(\mathcal{O}(d n \log(n))\) time where \(d\) is the number of features and \(n\) is the number of data points in the leaf. Therefore, no additional cost is induced by using LOO.

Figure 19 illustrates on a numerical example the difference on the estimated conditional quantiles when learning PMQRTs using or not without.

11.2.2 LOO for CRPS-RT↩︎

It is noteworthy that using leave-one-out to obtain an unbiased estimate of the information gain at each node for determining the best split is very simple. In fact, the computations presented below show that the leave-one-out estimate of the information gain can be efficiently derived by appropriately scaling the standard empirical entropies.

Denoting by \(\hat{F}^{-i}\) the empirical CDF obtained from the sequence \((y_1,\dots,y_{i-1},y_{i+1},\dots,y_n)\), we have \[\begin{align} &\frac{1}{n} \sum_{i=1}^n\mathrm{CRPS}(\widehat F^{-i},y_i)\\ & = \frac{1}{n} \sum_{i=1}^n \int_{-\infty}^{+\infty} \big( \hat{F}^{-i}(s) - \mathbb{1}_{s\geq y_i}\big)^2 ds\\ &= \frac{1}{n} \sum_{i=1}^n \int_{-\infty}^{+\infty} \big( \frac{1}{n-1} \sum_{j\neq i} \mathbb{1}_{s\geq y_j} - \mathbb{1}_{s\geq y_i}\big)^2 ds\\ &= \frac{1}{n(n-1)^2} \sum_{i=1}^n \sum_{j\neq i} \sum_{k\neq i}\int_{-\infty}^{+\infty} \big( \mathbb{1}_{s\geq y_j} - \mathbb{1}_{s\geq y_i}\big)\big( \mathbb{1}_{s\geq y_k} - \mathbb{1}_{s\geq y_i}\big) ds\\ &= \frac{1}{n(n-1)^2} \sum_{i=1}^n \sum_{j=1}^n \sum_{k=1}^n\int_{-\infty}^{+\infty} \big( \mathbb{1}_{s\geq y_j} - \mathbb{1}_{s\geq y_i}\big)\big( \mathbb{1}_{s\geq y_k} - \mathbb{1}_{s\geq y_i}\big) ds\\ &= \frac{n^2}{(n-1)^2} \left(\frac{1}{n} \sum_{i=1}^n\mathrm{CRPS}(\widehat F,y_i) \right). \end{align}\]

11.2.3 Optimism correction for CRPS-RT↩︎

We derive here the analogue of Mallows’ \(C_p\) for the CRPS risk when the predictive distribution is the empirical distribution function \[\widehat F(s)=\frac{1}{n}\sum_{i=1}^n \mathbb{1}_{\{Y_i\le s\}},\] constructed from an i.i.d.sample \(\mathcal{D}_n=(Y_1,\dots,Y_n)\) with common distribution function \(F\).

For a fixed predictive CDF \(\widehat F\), the empirical CRPS criterion is \[\begin{align} \widehat { R}_{1:n}(\widehat F)=H_{\mathrm{CRPS}}(\mathbf{y}) &= \frac{1}{n} \sum_{i=1}^n \int_{-\infty}^{+\infty} \left(\widehat F(s) - \mathbb{1}_{s\geq y_i}\right)^2 ds \\ &= \int_{-\infty}^{+\infty} \left( \widehat F(s)^2 -2 \widehat F(s) \frac{1}{n}\sum_{i=1}^n\mathbb{1}_{s\geq y_i} + \frac{1}{n}\sum_{i=1}^n\mathbb{1}_{s\geq y_i}\right) ds \\ &= \int_{-\infty}^{+\infty} \left( \widehat F(s)-\widehat F(s)^2 \right) ds . \end{align}\] The corresponding out-of-sample (prediction) risk is \[\begin{align} R(\widehat F) &= \mathbb{E}_{Y,\mathcal{D}_n}\left[\int_{-\infty}^{+\infty} \left(\widehat F(s) - \mathbb{1}_{s\geq Y}\right)^2 ds\right] \\ &= \int_{-\infty}^{+\infty} \mathbb{E}_{\mathcal{D}_n}\left[ \widehat F(s)^2 - 2 \widehat F(s) \mathbb{E}_{Y}[\mathbb{1}_{s\geq Y} ] + \mathbb{E}_{Y}[\mathbb{1}_{s\geq Y}^2 ] \right]ds\\ &= \int_{-\infty}^{+\infty} \mathbb{E}_{\mathcal{D}_n}\left[ \widehat F(s)^2 - 2 \widehat F(s) F(s) + F(s) \right]ds\\ &= \int_{-\infty}^{+\infty}\left( \mathbb{E}_{\mathcal{D}_n}\left[ \widehat F(s)^2 \right] - 2 F(s)^2 + F(s) \right) ds. \end{align}\] where \(Y\) is an independent future draw from \(F\). We used that \(\mathbb{E}[\mathbb{1}_{\{Y\le s\}}]=F(s)\) and \(\mathbb{E}[\mathbb{1}_{\{Y\le s\}}^2]=F(s)\). Therefore the expected optimism is \[\begin{align} \mathrm{pen}_n&=\mathbb{E}_{\mathcal{D}_n}\left[ R(\widehat F) - \widehat { R}_{1:n}(\widehat F) \right] \\ &= \int_{-\infty}^{+\infty}\left(2\mathbb{E}_{\mathcal{D}_n}\left[\widehat F(s)^2 \right]-2F(s)^2 \right)ds = \int_{-\infty}^{+\infty} 2 \mathrm{Var}_{\mathcal{D}_n}(\widehat F(s))ds\\ &= \int_{-\infty}^{+\infty} 2 \mathrm{Var}_{\mathcal{D}_n}(\frac{1}{n} \sum_{i=1}^n \mathbb{1}_{s\geq y_i}) ds\\ &= \int_{-\infty}^{+\infty} \frac{2}{n} \mathrm{Var}_{\mathcal{D}_n}( \mathbb{1}_{s\geq y_1}) ds\\ &= \int_{-\infty}^{+\infty} \frac{2}{n} F(s)(1-F(s))ds. \end{align}\]

This suggests the plug-in estimator of \(\mathrm{pen}_n\) given by: \[\widehat{\mathrm{pen}}_{1:n}: =\int_{-\infty}^{+\infty} \frac{2}{n-1} \widehat F(s)(1-\widehat F(s))ds =\frac{2}{n-1} \widehat R _{1:n}(\widehat F)=\frac{2}{n-1} H_{\mathrm{CRPS}}(\mathbf{y}).\]

Normalization by \(n-1\) ensures that the plug-in estimator is exactly unbiased. Indeed, for every \(s\in\mathbb{R}\), \[\mathbb{E}_{\mathcal{D}_n}[\widehat F(s)]=F(s) \qquad\text{and}\qquad \mathrm{Var}_{\mathcal{D}_n}(\widehat F(s))=\frac{1}{n}F(s)\bigl(1-F(s)\bigr),\] so that \[\begin{align} \mathbb{E}_{\mathcal{D}_n}\!\left[\widehat F(s)\bigl(1-\widehat F(s)\bigr)\right] &= \mathbb{E}_{\mathcal{D}_n}[\widehat F(s)]-\mathbb{E}_{\mathcal{D}_n}[\widehat F(s)^2] \\ &= F(s)-\Bigl(\mathrm{Var}_{\mathcal{D}_n}(\widehat F(s))+\mathbb{E}_{\mathcal{D}_n}[\widehat F(s)]^2\Bigr) \\ &= F(s)-\frac{1}{n}F(s)(1-F(s))-F(s)^2 \\ &= \left(1-\frac{1}{n}\right)F(s)\bigl(1-F(s)\bigr). \end{align}\] Integrating over \(s\) yields \[\mathbb{E}_{\mathcal{D}_n}\!\left[H_{\mathrm{CRPS}}(\mathbf{y})\right] = \left(1-\frac{1}{n}\right)\int_{-\infty}^{+\infty}F(s)\bigl(1-F(s)\bigr)\,ds = \frac{n-1}{n}\int_{-\infty}^{+\infty}F(s)\bigl(1-F(s)\bigr)\,ds.\] Therefore, \[\begin{align} \mathbb{E}_{\mathcal{D}_n}\!\left[\frac{2}{n-1}H_{\mathrm{CRPS}}(\mathbf{y})\right] &= \frac{2}{n-1}\cdot \frac{n-1}{n}\int_{-\infty}^{+\infty}F(s)\bigl(1-F(s)\bigr)\,ds \\ &= \frac{2}{n}\int_{-\infty}^{+\infty}F(s)\bigl(1-F(s)\bigr)\,ds = \mathrm{pen}_n, \end{align}\] which proves that the normalization by \(n-1\) yields an exactly unbiased estimator of the optimism term.

Hence the counterpart of Mallows \(C_p\) for the CRPS loss is given by: \[\begin{align} &H_{\mathrm{CRPS}}(\mathbf{y}) + \frac{2}{n-1} H_{\mathrm{CRPS}}(\mathbf{y}) = \frac{n+1}{n-1}H_{\mathrm{CRPS}}(\mathbf{y}). \end{align}\] In other words, CRPS-RT admits a simple closed-form optimism correction, so that model selection can be based on a direct multiplicative rescaling of the in-sample CRPS criterion.

In the regression tree setting, this correction is applied locally at each candidate split: the distribution function \(F\) should therefore be understood as the conditional cumulative distribution function of the response \(Y\) given that the covariate vector \(X\) belongs to the node (or child node) defined by the sequence of previous splits together with the split currently under consideration, so that the resulting Mallows-style correction provides an unbiased estimate of the corresponding CRPS information gain.

12 Datasets↩︎

Table 7 lists the datasets used in our experiments, along with their corresponding URLs.

References↩︎

---

1. github.com/quentin-duchemin/DisTreebution/↩︎

2. Footnote 6: Wager, S., & Athey, S. (2018). Estimation and inference of heterogeneous treatment effects using random forests. Journal of the American Statistical Association, 113(523), 1228-1242.↩︎