Distributional Random Forests: Heterogeneity Adjustment and Multivariate Distributional Regression

Random Forest (Breiman, 2001) is a successful and widely used regression and classification algorithm. Part of its appeal and reason for its versatility is its (implicit) construction of a kernel-type weighting function on training data, which can also be used for targets other than the original mean estimation. We propose a novel forest construction for multivariate responses based on their joint conditional distribution, independent of the estimation target and the data model. It uses a new splitting criterion based on the MMD distributional metric, which is suitable for detecting heterogeneity in multivariate distributions. The induced weights define an estimate of the full conditional distribution, which in turn can be used for arbitrary and potentially complicated targets of interest. The method is very versatile and convenient to use, as we illustrate on a wide range of examples. The code is available as Python and R packages drf.


Introduction
In practice, one often encounters heterogeneous data, whose distribution is not constant, but depends on certain covariates. For example, data can be collected from several different sources, its distribution might differ across certain subpopulations or it could even change with time, etc. Inferring valid conclusions about a certain target of interest from such data can be very challenging as many different aspects of the distribution could potentially change.
As an example, in medical studies, the effectiveness of a certain treatment might not be constant throughout the population but depend on certain patient characteristics such as age, race, gender, or medical history. Another issue could be that different patient groups were not equally likely to receive the same treatment in the observed data.
Obviously, pooling all available data together can result in invalid conclusions. On the other hand, if for a given test point of interest one only considers similar training data points, i.e. a small homogeneous subpopulation, one may end up with too few samples for accurate statistical estimation. In this paper, we propose a method based on the Random Forest algorithm (Breiman, 2001) which in a data-adaptive way determines for any given test point which training data points are relevant for it. This in turn can be used for drawing valid conclusions or for accurately estimating any quantity of interest.
Let Y = (Y 1 , Y 2 , . . . , Y d ) ∈ R d be a multivariate random variable representing the data of interest, but whose joint distribution is heterogeneous and depends on some subset of a potentially large number of covariates X = (X 1 , X 2 , . . . , X p ) ∈ R p . Throughout the paper, vector quantities are denoted in bold. We aim to estimate a certain target object τ (x) that depends on the conditional distribution P(Y | X = x) = P(Y | X 1 = x 1 , . . . , X p = x p ), where x = (x 1 , . . . , x p ) is an arbitrary point in R p . The estimation target τ (x) can range from simple quantities, such as the conditional expectations E[f (Y) | X] (Breiman, 2001) or quantiles Q α [f (Y) | X] (Meinshausen, 2006) for some function f : R d → R, to some more complicated aspects of the conditional distribution P(Y | X = x), such as conditional copulas or conditional independence measures. Given the observed data {(x i , y i )} n i=1 , the most straightforward way of estimating τ (x) nonparametrically would be to consider only the data points in some neighborhood N x around x, e.g. by considering the k nearest neighbors according to some metric. However, such methods typically suffer from the curse of dimensionality even when p is only moderately large: for a reasonably small neighborhood, such that the distribution P(Y | X ∈ N x ) is close to the distribution P(Y | X = x), the number of training data points contained in N x will be very small, thus making the accurate estimation of the target τ (x) difficult. The same phenomenon occurs with other methods which locally weight the training observations such as kernel methods (Silverman, 1986), local MLE (Fan et al., 1998) or weighted regression (Cleveland, 1979) even for the relatively simple problem of estimating the conditional mean E[Y | X = x] for fairly small p. For that reason, more importance should be given to the training data points (x i , y i ) for which the response distribution P(Y | X = x i ) at point x i is similar to the target distribution P(Y | X = x), even if x i is not necessarily close to x in every component.
In this paper, we propose the Distributional Random Forest (DRF) algorithm which estimates the multivariate conditional distribution P(Y | X = x) in a locally adaptive fashion. This is done by repeatedly dividing the data points in the spirit of the Random Forest algorithm (Breiman, 2001): at each step, we split the data points into two groups based on some feature X j in such a way that the distribution of Y for which X j ≤ l, for some level l, differs the most compared to the distribution of Y when X j > l, according to some distributional metric. One can use any multivariate two-sample test statistic, provided it can detect a wide variety of distributional changes. As the default choice, we propose a criterion based on the Maximal Mean Discrepancy (MMD) statistic (Gretton et al., 2007a) with many interesting properties. This splitting procedure partitions the data such that the distribution of the multivariate response Y in the resulting leaf nodes is as homogeneous as possible, thus defining neighborhoods of relevant training data points for every x. Repeating this many times with randomization induces a weighting function w x (x i ) as in Jeon (2002, 2006), described in detail in Section 2, which quantifies the relevance of each training data point x i for a given test point x. The conditional distribution is then estimated by an empirical distribution determined by these weights (Meinshausen, 2006). This construction is data-adaptive as it assigns more weight to the training points x i that are closer to the test point x in the components which are more relevant for the distribution of Y.
Our forest construction does not depend on the estimation target τ (x), but it rather estimates the conditional distribution P(Y | X = x) directly and the induced forest weights can be used to estimate τ (x) in a second step. This approach has several advantages. First, only one DRF fit is required to obtain estimates of many different targets, which has a big computational advantage. Furthermore, since those estimates are obtained from the same forest fit, they are mutually compatible. For example, if the conditional correlation matrix {Cor(Y i , Y j | X = x)} d i,j=1 is estimated componentwise using some other method, the resulting matrix might not be positive semidefinite, and as another example, the CDF estimatesP(Y ≤ y | X = x) might not be monotone in y, see Figure 6. Finally, it could be extremely difficult to tailor forest construction to some complex targets τ (x). The induced weighting function can thus be used not only for obtaining simple distributional aspects such as, for example, the conditional quantiles, conditional correlations, or joint conditional probability statements, but also to obtain more complex objectives, such as conditional independence tests (Zhang et al., 2011), heterogeneous regression (see also Section 4.4 for more details) (Künzel et al., 2019;Wager and Athey, 2018) or semiparametric estimation by fitting a parametric model for Y, having nonparametrically adjusted for X (Bickel et al., 1993). Representation of the conditional distribution via the weighting function has a great potential for applications in causality such as causal effect estimation or as a way of implementing do-calculus (Pearl, 2009) for finite samples, as we discuss in Section 4.4.
Therefore, DRF is used in two steps: in the first step, we obtain the weighting function w x (·) describing the conditional distribution P(Y | X = x) in a target-and model-free way, which is then used as an input for the second step. Even if the method used in the second step does not directly support weighting of the training data points, one can easily resample the data set with the sampling probabilities equal to {w x (x i )} n i=1 . This two-step approach is visualized in the following diagram: 1) get w x (·) with DRF objective 2) compute from w x (·) induced estimator

Related work and our contribution
Several adaptations of the Random Forest algorithm have been proposed for targets beyond the original one of the univariate conditional mean E[Y | X = x]: for survival analysis (Hothorn et al., 2006), conditional quantiles (Meinshausen, 2006), density estimation (Pospisil and Lee, 2018), CDF estimation (Hothorn and Zeileis, 2021) or heterogeneous treatment effects (Wager and Athey, 2018). Almost all such methods use the weights induced by the forest, as described in Section 2, rather than averaging the estimates obtained per tree. This view of Random Forests as a powerful adaptive nearest neighbor method is well known and dates back to Jeon (2002, 2006). It was first used for targets beyond the conditional mean in Meinshausen (2006), where the original forest construction with univariate Y was used (Breiman, 2001). However, the univariate response setting considered there severely restricts the number of interesting targets τ (x) and DRF can thus be viewed as an important generalization of this approach to the multivariate setting.
In order to be able to perform certain tasks or to achieve a better accuracy, many forest-based methods adapt the forest construction by using a custom splitting criterion tailored to their specific target, instead of relying on the standard CART criterion. In Zeileis et al. (2008) and Hothorn and Zeileis (2021), a parametric model for the response Y | X = x ∼ f (θ(x), ·) is assumed and recursive splitting is performed based on a permutation test which uses the user-provided score functions. Similarly,  estimate certain univariate targets for which there exist corresponding score functions defining the local estimating equations. The data is split so that the estimates of the target in resulting child nodes differ the most. This is different, though, to the target-free splitting criterion of DRF, which splits so that the distribution of Y in child nodes is as different as possible.
Since the splitting step is extensively used in the algorithm, its complexity is crucial for the overall computational efficiency of the method, and one often needs to resort to approximating the splitting criterion (Pospisil and Lee, 2018; to obtain good computational run time. We propose a splitting criterion based on a fast random approximation of the MMD statistic (Gretton et al., 2012a;Zhao and Meng, 2015), which is commonly used in practice for two-sample testing as it is able to detect any change in the multivariate distribution of Y with good power (Gretton et al., 2007a). DRF with the MMD splitting criterion also has interesting theoretical properties as shown in Section 3 below.
The multivariate response case has not received much attention in the Random Forest literature. Most of the existing forest-based methods focus on either a univariate response Y or on a certain univariate target τ (x). One interesting line of work considers density estimation (Pospisil and Lee, 2018) and uses aggregation of the CART criteria for different response transformations. Another approach (Kocev et al., 2007;Segal and Xiao, 2011;Ishwaran and Kogalur, 2022) is based on aggregating standard univariate CART splitting criteria for Y 1 , . . . , Y d and targets only the conditional mean of the responses, a task which could also be solved by separate regression fits for each Y i . In order to capture any change in the distribution of the multivariate response Y, one needs to not only consider the marginal distributions for each component Y i , but also to determine whether their dependence structure changes, see e.g. Figure 8.
There is an increasing number of methods that nonparametrically estimate the joint multivariate conditional distribution P(Y | X = x) in the statistics and machine learning literature. In addition to a few simple classical methods such as k-nearest neighbors and kernel regression, there exist methods based on normalizing flows such as Inverse Autoregressive Flow (Kingma et al., 2016) or Masked Autoregressive Flow (Papamakarios et al., 2017) and also conditional variants of several popular generative models such as Conditional Generative Adversarial Networks (Mirza and Osindero, 2014) or Conditional Variational Autoencoder (Sohn et al., 2015). The focus of these methods is more on the settings with large response dimension d and small covariate dimension p, such as image or text generation. Another interesting and related line of research focuses on estimating the conditional mean embedding (CME), as described e.g., in Song et al. (2009); Muandet et al. (2017); Song et al. (2013); Park and Muandet (2020), rather than estimating the conditional distribution directly. CMEs generalize the concept of embedding (marginal) probability distributions into a Reproducing Kernel Hilbert Space (RKHS) to the conditional case. Interestingly, DRF with the MMD-based splitting criterion can also be viewed as a method for estimating the CME, as discussed in Section 3 below. This viewpoint provides a natural connection between the Random Forest and kernel embedding literature. A comparison of DRF with the methods for distributional estimation listed above can be found in Section 4.1.
Our contribution, resulting in the proposal of the Distributional Random Forest (DRF), can be summarized as follows: First, we introduce the idea of forest construction based on sequential multivariate two-sample test statistics. It does not depend on a particular estimation target and is completely nonparametric, which makes its implementation and usage very simple and universal. Not only does it not require additional user input such as the log-likelihoods or score functions, but it can be used even for complicated targets for which there is no obvious forest construction. Furthermore, it has a computational advantage as only a single forest fit is needed for producing estimates of many different targets that are additionally compatible with each other. Second, we propose an MMD-based splitting criterion with good statistical and computational properties, for which we also derive interesting theoretical results in Section 3. It underpins our implementation, which we provide as R and Python packages drf. Finally, we show on a broad range of examples in Section 4 how many different statistical estimation problems, some of which not being easily tractable by existing forest-based methods, can be cast to our framework, thus illustrating the usefulness and versatility of DRF.

Method
In this section we describe the details of the Distributional Random Forest (DRF) algorithm. We closely follow the implementations of the grf  and ranger (Wright and Ziegler, 2017) R-packages. A detailed description of the method and its implementation and the corresponding pseudocode can be found in the Appendix A.

Forest Building
The trees are grown recursively in a model-free and target-free way as follows: For every parent node P , we determine how to best split it into two child nodes of the form C L = {X j ≤ l} and C R = {X j > l}, where the variable X j is one of the randomly chosen splitting candidates and l denotes its level based on which we perform the splitting. The split is chosen such that we maximize a certain (multivariate) two-sample test statistic which measures the difference of the empirical distributions of the data Y in the two resulting child nodes C L and C R . Therefore, in each step we select the candidate predictor X j which seems to affect the distribution of Y the most, as measured by the metric D(·, ·). Intuitively, in this way we ensure that the distribution of the data points in every leaf of the resulting tree is as homogeneous as possible, which helps mitigate the bias caused by pooling the heterogeneous data together. A related idea can be found in GRF , where one attempts to split the data so that the resulting estimatesτ L andτ R , obtained respectively from data points in C L and C R , differ the most: where we write n P = |{i | x i ∈ P }| and n L , n R are defined analogously. One could construct the forest using any metric D(·, ·) for empirical distributions. However, in order to have a good accuracy of the overall method, the corresponding twosample test using D(·, ·) needs to have a good power for detecting any kind of change in distribution, which is a difficult task in general, especially for multivariate data (Bai and Saranadasa, 1996;Székely and Rizzo, 2004). Another very important aspect of the choice of distributional metric D(·, ·) is the computational efficiency; one needs to be able to sequentially compute the values of D ({y i | x i ∈ C L } , {y i | x i ∈ C R }) for every possible split very fast for the overall algorithm to be computationally feasible, even for moderately large data sets. Below, we propose a splitting criterion based on the MMD two-sample test statistic (Gretton et al., 2007a) which has both good statistical and computational properties.
In contrast to other forest-based methods, we do not use any information about our estimation target τ in order to find the best split of the data, which comes with a certain trade-off. On one hand, it is sensible that tailoring the splitting criterion to the target should improve the estimation accuracy; for example, some predictors might affect the conditional distribution of Y, but not necessarily the estimation target τ and splitting on such predictors unnecessarily reduces the number of training points used for estimating τ . On the other hand, our approach has multiple benefits: it is easier to use as it does not require any user input such as the likelihood or score functions and it can also be used for very complicated targets for which one could not easily adapt the splitting criterion. Furthermore, only one DRF fit is necessary for producing estimates of many different targets, which has both computational advantage and the practical advantage that the resulting estimates are mutually compatible (see e.g. Figure 5).
Interestingly, sometimes it could even be beneficial to split based on a predictor which does not affect the target of estimation, but which affects the conditional distribution. This is illustrated by the following toy example. Suppose that for a bivariate response (Y 1 , Y 2 ) we are interested in estimating the slope of the linear regression of Y 2 on Y 1 conditionally on p = 30 predictors X, i.e. our target is τ ( . This is one of the main use cases for GRF and its variant which estimates this target is called Causal Forest (Wager and Athey, 2018;. Let us assume that the data has the following distribution: i.e. X 1 affects only the mean of the responses, while the other p − 1 predictors have no effect. In Figure 1 we illustrate the distribution of the data when n = 300, p = 30, σ = 0.2, together with the DRF and GRF splitting criteria. The true value of the target is τ (x) = 0, but when σ is not too big, the slope estimatesτ on pooled data will be closer to 1. Therefore, the difference ofτ L andτ R between the induced slope estimates for a candidate split, which is used for splitting criterion (2) of GRF, might not be large enough for us to decide to split on X 1 , or the resulting split might be too unbalanced. This results in worse forest estimates for this toy example, see Figure 1.  (2) splitting criteria (right) at the root node. The curves of different colors correspond to different predictors, with X 1 denoted in black. Top right: Comparison of the estimates of DRF and Causal Forest  which respectively use those splitting criteria. Test points were randomly generated from the same distribution as the training data. Black dashed line indicates the correct value of the target quantity.

Weighting Function
Having constructed our forest, just as the standard Random Forest (Breiman, 2001) can be viewed as the weighted nearest neighbor method (Lin and Jeon, 2002), we can use the induced weighting function to estimate the conditional distribution at any given test point x and thus any other quantity of interest τ (x). This approach is commonly used in various forest-based methods for obtaining predictions, see e.g., Hothorn and Zeileis (2021); Pospisil and Lee (2018); .
Suppose that we have built N trees T 1 , . . . , T N . Let L k (x) be the set of the training data points which end up in the same leaf as x in the tree T k . The weighting function w x (x i ) is defined as the average of the corresponding weighting functions per tree (Lin and Jeon, 2006): The weights are positive and add up to 1: In the case of equally sized leaf nodes, the assigned weight to a training point x i is proportional to the number of trees where the test point x and x i end up in the same leaf node. This shows that forest-based methods can in general be viewed as adaptive nearest neighbor methods. The sets L k (x) of DRF will contain data points (x i , y i ) such that P(Y | X = x i ) is close to P(Y | X = x), thus removing bias due to heterogeneity of Y caused by X. On the other hand, since the trees are constructed randomly and are thus fairly independent (Breiman, 2001), the leaf sets L k (x) will be different enough so that the induced weights w x (x i ) are not concentrated on a small set of data points, which would lead to high estimation variance. Such good bias-variance tradeoff properties of forest-based methods are also implied by their asymptotic properties (Biau, 2012;Wager, 2014), even though this is a still active area of research and not much can be shown rigorously.
One can estimate the conditional distribution P(Y | X = x) from the weighting function by using the corresponding empirical distribution: where δ y i is the point mass at y i . Two-step approach using weights. The weighting function w x (x i ) can directly be used for any target τ (x) in a second step and not just for estimating the conditional distribution. For example, the estimated conditional joint CDF is given bŷ It is important to point out that using the induced weighting function for locally weighted estimation is different than the approach of averaging the noisy estimates obtained per tree (Wager and Athey, 2018), originally used in standard Random Forests (Breiman, 2001). Even though the two approaches are equivalent for conditional mean estimation, the former approach is often much more efficient for more complex targets , since the number of data points in a single leaf is very small, leading to large variance of the estimates.
For the univariate response, the idea of using the induced weights for estimating targets different than the original target of conditional mean considered in Breiman (2001) dates back to Quantile Regression Forests (QRF) (Meinshausen, 2006), where a lot of emphasis is put on the quantile estimation, as the number of interesting targets is quite limited in the univariate setting. In the multivariate case, on the other hand, many interesting quantities such as, for example, conditional quantiles, conditional correlations or various conditional probability statements can easily be directly estimated from the weights. By using the weights as an input for some other method, we can accomplish some more complicated objectives, such as conditional independence testing, causal effect estimation, Figure 2: Top: the characteristics of the important training sites, for a fixed test site whose position is indicated by a black star and whose characteristics are indicated in the title. The total weight assigned corresponds to the symbol size. Bottom: estimated joint conditional distribution of two pollutants NO 2 and PM2.5, where the weights correspond to the transparency of the data points. Green area corresponds to 'Good' air quality category (AQI ≤ 50). semiparametric learning, time series prediction or tail-index estimation in extreme value analysis. As an example, suppose that our data Y come from a certain parametric model, where the parameter θ is not constant, but depends on X instead, i.e. Y | X = x ∼ f (θ(x), ·), see also Zeileis et al. (2008). One can then estimate the parameter θ(x) by using weighted maximum likelihood estimation: Another example is heterogeneous regression, where we are interested in the regression fit of an outcome Y ∈ R on certain predicting variables W ∈ R s conditionally on some event {X = x}. This can be achieved by weighted regression of Y on W, where the weights w x (x i ) assigned to each data point (w i , y i ) are obtained from DRF with the multivariate response (Y, W) ∈ R s+1 and predictors X ∈ R p , for an illustration see Section 4.4. The weighting function of DRF is illustrated on the air quality data in Figure 2. Five years (2015 − 2019) of air pollution measurements were obtained from the US Environmental Protection Agency (EPA) website. Six main air pollutants (nitrogen dioxide (NO 2 ), carbon monoxide (CO), sulphur dioxide (SO 2 ), ozone (O 3 ) and coarse and fine particulate matter (PM10 and PM2.5)) that form the air quality index (AQI) were measured at many different measuring sites in the US for which we know the longitude, latitude, elevation, location setting (rural, urban, suburban) and how the land is used within a 1/4 mile radius. Suppose we would want to know the distribution of the pollutant measurements at some new, unobserved, measurement site. We train DRF with the measurements (intraday maximum) of the two pollutants PM2.5 and NO 2 as the responses, and the site longitude, latitude, elevation, land use and location settings as the predictors and choose two decommissioned measurement sites as test points. For each test point we obtain the weights to all training measurements. We further combine the weights for all measurements corresponding to the same site. The top row illustrates for a given test site, whose characteristics are indicated in the plot title, how much weight in total is assigned to the measurements from a specific training site. We see that the important sites share many characteristics with the test site and that DRF determines the relevance of each characteristic in a data-adaptive way. The bottom row shows the corresponding estimates of the joint conditional distribution of the pollutants (we choose 2 of them for visualization purposes), where the transparency of each training point reflects the assigned weight. One can clearly see how the estimated pollution levels are larger for the suburban site than for the rural site. The forest weights can be used, for example, for estimating the joint density (whose contours can be seen in the plot) or for estimating the probability that the AQI is below a certain value by summing the weights in the corresponding region of space.

Distributional Metric
In order to determine the best split of a parent node P , i.e. such that the distributions of the responses Y in the resulting child nodes C L and C R differ the most, one needs a good distributional metric D(·, ·) (see Equation (1)) which can detect change in distribution of the response Y when additionally conditioning on an event {X j > l}. Testing equality of distributions from the corresponding samples is an old problem in statistics, called twosample testing problem. For univariate data, many good tests exist such as Wilcoxon rank test (Wilcoxon, 1946), Welch's t-test (Welch, 1947, Wasserstein two-sample testing (Ramdas et al., 2017), Kolmogorov-Smirnov test (Massey Jr, 1951) and many others, but obtaining an efficient test for multivariate distributions has proven to be quite challenging due to the curse of dimensionality (Friedman and Rafsky, 1979;Baringhaus and Franz, 2004).
Additional requirement for the choice of distributional metric D(·, ·) used for data splitting is that it needs to be computationally very efficient as splitting is used extensively in the algorithm. If we construct N trees from n data points and in each node we consider mtry candidate variables for splitting, the complexity of the standard Random Forest algorithm (Breiman, 2001) in the univariate case is O(N × mtry × n log n) provided our splits are balanced. It uses the CART splitting criterion, given by: where y L = 1 n L x i ∈C L y i and y R is defined analogously. This criterion has an advantage that not only it can be computed in O(n P ) complexity, but this can be done for all possible splits {X j ≤ l} as cutoff level l varies, since updating the splitting criterion when moving a single training data point from one child node to the other requires only O(1) computational steps (most easily seen by rewriting the CART criterion as in (13)).
If the time complexity of evaluating the DRF splitting criterion (1) for a single splitting candidate X j and all cutoffs l of interest (usually taken to range over all possible values) is at least n c for some c > 1, say O(f (n)) for some function f : R → R, then by solving the recursive relation we obtain that the overall complexity of the method is given by O(N × mtry × f (n)) (Akra and Bazzi, 1998), which can be unfeasible even for moderately large n if f grows too fast.
The problem of sequential two-sample testing is also central to the field of change-point detection (Wolfe and Schechtman, 1984;Brodsky and Darkhovsky, 2013), with the slight difference that in the change-point problems the distribution is assumed to change abruptly at certain points in time, whereas for our forest construction we only are interested in finding the best split of the form {X j ≤ l} and the conditional distribution P(Y | {X ∈ P } ∩ {X j ≤ l}) usually changes gradually with l. The testing power and the computational feasibility of the method play a big role in change-point detection as well. However, the state-of-the-art change-point detection algorithms (Li et al., 2019;Matteson and James, 2014) are often too slow for our purpose as sequential testing is done O(N × mtry × n) times for forest construction, much more frequently than in change-point problems.

MMD splitting criterion
Even though DRF could in theory be constructed with any distributional metric D(·, ·), as a default choice we propose splitting criterion based on the Maximum Mean Discrepancy (MMD) statistic (Gretton et al., 2007a). Let (H, ·, · H ) be the RKHS of real-valued functions on R d induced by some positive-definite kernel k, and let ϕ : R d → H be the corresponding feature map satisfying that k(u, v) = ϕ(u), ϕ(v) H .
The MMD statistic D MMD(k) (U, V ) for kernel k and two samples U = {u 1 , . . . , u |U | } and V = {v 1 , . . . , v |V | } is given by: MMD compares the similarities, described by the kernel k, within each sample with the similarities across samples and is commonly used in practice for two-sample testing. It is based on the idea that one can assign to each distribution P its embedding µ(P) into the RKHS H, which is the unique element of H given by The MMD two-sample statistic (8) can then equivalently be written as the squared distance between the embeddings of the empirical distributions with respect to the RKHS norm · H : recalling that δ y is the point mass at y. As the sample sizes |U | and |V | grow, the MMD statistic (10) converges to its population version, which is the squared RKHS distance between the corresponding embeddings of the data-generating distributions of U and V . Since the embedding map µ is injective for a characteristic kernel k, we see that MMD is able to detect any difference in the distribution. Even though the power of the MMD two sample test also deteriorates as the data dimensionality grows, since the testing problem becomes intrinsically harder , it still has good empirical power compared to other multivariate two-sample tests for a wide range of k (Gretton et al., 2012a).
Fast random splitting criterion approximation. The O((|U | + |V |) 2 ) complexity for computing D MMD(k) (U, V ) from (8) is nevertheless too large for many applications. For that reason, several fast approximations of MMD have been suggested in the literature (Gretton et al., 2012a,b;Zaremba et al., 2013;Chwialkowski et al., 2015;Jitkrittum et al., 2016). As already mentioned, the complexity of the distributional metric D(·, ·) used for DRF is crucial for the overall method to be computationally efficient, since the splitting step is used extensively in the forest construction. We therefore propose splitting based on an MMD statistic computed with an approximate kernelk, which is also a fast random approximation of the original MMD statistic (Zhao and Meng, 2015).
Bochner's theorem (see e.g. Wendland (2004, Theorem 6.6)) gives us that any bounded shift-invariant kernel can be written as i.e. as a Fourier transform of some measure ν. Therefore, by randomly sampling the frequency vectors ω 1 , . . . , ω B from normalized ν, we can approximate our kernel k by another kernelk (up to a scaling factor) as follows: where we definek(u, v) = ϕ(u), ϕ(v) C B as the kernel function with the feature map given by which is a random vector consisting of the Fourier features ϕ ω (u) = e iω T u ∈ C (Rahimi and Recht, 2008). Such kernel approximations are frequently used in practice for computational efficiency (Rahimi and Recht, 2009;Le et al., 2013). As a default choice of k we take the Gaussian kernel with bandwidth σ, since in this case we have a convenient expression for the measure ν and we sample ω 1 , . . . , The bandwidth σ is chosen as the median pairwise distance between all training responses {y i } n i=1 , commonly referred to as the 'median heuristic' (Gretton et al., 2012c).
From the representation of MMD via the distribution embeddings (10), we can obtain that MMD two-sample test statistic D MMD(k) using the approximate kernelk is given by Interestingly, D MMD(k) is not only an MMD statistic on its own, but can also be viewed as a random approximation of the original MMD statistic D MMD(k) (8) using kernel k; by using the kernel representation (11), it can be written as Finally, our DRF splitting criterion D(·, ·) (1) is then taken to be the (scaled) MMD statistic n L n R n 2 with the approximate random kernel k used instead of k, which can thus be conveniently written as: where we recall that n P = |{i | x i ∈ P }| and n L , n R are defined analogously. The additional scaling factor n L n R n 2 P in (12) occurs naturally and compensates the increased variance of the test statistic for unbalanced splits; it also appears in the GRF (2) and CART (see representation (13)) splitting criteria. The main advantage of the splitting criterion based on D MMD(k) is that by using the representation (1) it can be easily computed for every possible splitting level l in O(Bn P ) complexity, whereas the MMD statistic D MMD(k) using kernel k would require O(n 2 P ) computational steps, which makes the overall complexity of the algorithm O (B × N × mtry × n log n) instead of much slower O N × mtry × n 2 .
We do not use the same approximate random kernelk for different splits; for every parent node P we resample the frequency vectors {ω b } B b=1 defining the corresponding feature mapφ. Using differentk at each node might help to better detect different distributional changes. Furthermore, having different random kernels for each node agrees well with the randomness of the Random Forests and helps making the trees more independent. Since the MMD statistic D MMD(k) used for our splitting criterion is not only an approximation of D MMD(k) , but is itself an MMD statistic, it inherits good power for detecting any difference in distribution of Y in the child nodes for moderately large data dimensionality d, even when B is reasonably small. One could even consider changing the number of random Fourier features B at different levels of the tree, as n P varies, but for simplicity we take it to be fixed.
Relationship to CART. There is some similarity of our MMD-based splitting criterion (12) with the standard variance reduction CART splitting criterion (7) when d = 1, which can be rewritten as: The derivation can be found in Appendix B. From this representation, we see that the CART splitting criterion (7) is also equivalent to the GRF splitting criterion (2) when our target is the univariate conditional mean τ (x) = E[Y | X = x] which is estimated for C L and C R by the sample meansτ L = y L andτ R = y R . Therefore, as it compares the means of the univariate response Y in the child nodes, the CART criterion can only detect changes in the response mean well, which is sufficient for prediction of Y from X, but might not be suitable for more complex targets. Similarly, for multivariate applications, aggregating the marginal CART criteria (Kocev et al., 2007;Segal and Xiao, 2011) across different components Y i of the response can only detect changes in the means of their marginal distributions. However, it is possible in the multivariate case that the pairwise correlations or the variances of the responses change, while the marginal means stay (almost) constant. For an illustration on simulated data, see Figure 7. Additionally, aggregating the splitting criteria over d components of the response Y can reduce the signal size if only the distribution of a few components change. Our MMD-based splitting criterion (12) is able to avoid such difficulties as it implicitly inspects all aspects of the multivariate response distribution.
If one takes a trivial kernel k id (y i , y j ) = y i y j with the identity feature map ϕ id (y) = y, the distributional embedding (9) is given by µ(P) = E Y ∼P [Y ] and thus the corresponding splitting criterion based on D MMD(k id ) (10) is exactly equal to the CART splitting criterion (7), which can be seen from its equivalent representation (13). Interestingly, Theorem 1 in Section 3 shows that the MMD splitting criterion with general kernel k can also be viewed as the abstract version of the CART criterion in the RKHS H corresponding to k (Fan et al., 2010), with the response variable being the feature map ϕ(Y) ∈ H. Therefore, DRF with the MMD splitting criterion can also be viewed as a forest-based method for estimation of the conditional embedding, which further justifies the proposed method. In Section 3 below, we use this relationship to derive interesting theoretical properties of DRF with the MMD splitting criterion.

Theoretical Results
In this section we first use the properties of the kernel mean embedding in order to relate DRF with the MMD splitting criterion to an abstract version of the standard Random Forest with the CART splitting criterion (Breiman, 2001), where the response is taking values in the corresponding RKHS. This representation reveals that DRF with the MMD splitting criterion can be viewed as a Random Forest estimator of the conditional mean embedding (CME) (Park and Muandet, 2020), similarly as the standard Random Forest estimates the conditional mean. This relationship is further exploited to adapt the existing theoretical results from the Random Forest literature to show that our estimate (5) of the conditional distribution of the response is consistent with respect to the MMD metric for probability measures and with a good rate. Finally, we show that this implies consistency of the induced DRF estimates for a range interesting targets τ (x), such as conditional CDFs or quantiles. The proofs of all results can be found in the Appendix B.

Casting DRF as a Random Forest in an RKHS
Recalling the notation from above, let (H, ·, · H ) be the Reproducing kernel Hilbert space induced by the positive definite kernel k : R d × R d → R and let ϕ : R d → H be its corresponding feature map. The kernel embedding function µ : M b (R d ) → H maps any bounded signed Borel measure P on R d to an element µ(P) ∈ H defined by see also (9). Boundedness of k ensures that µ is indeed defined on all of M b (R d ), while continuity of k ensures that H is separable (Hsing and Eubank, 2015). By considering the kernel embedding µ(·) and using its linearity, the embedding of the distributional estimate µ(P(Y | X = x)) of DRF (5) can be written as the average of the embeddings of the empirical distributions of Y in the leaves containing x over all trees: This is analogous to the prediction of the response for the standard univariate Random Forest, but where we average the embeddings µ(δ y i ) = ϕ(y i ) ∈ H instead of the response values y i themselves.
Furthermore, one can relate the MMD splitting criterion to the original CART criterion (7), which measures the mean squared prediction error for splitting a certain parent node P into children C L and C R . On one hand, from Equation (13) we see that the CART criterion also measures the squared distance between the response averages 1 n L x i ∈C L y i and 1 n R x i ∈C R y i in the child nodes, but on the other hand, Equation (10) shows that the MMD splitting criterion measures the RKHS distance between the embeddings of the empirical response distributions in C L and C R . This is summarized in the following theorem, which not only shows that the MMD splitting criterion can be viewed as the abstract CART criterion in the RKHS H (Fan et al., 2010), but also that DRF with the MMD splitting criterion can be viewed asymptotically as a greedy minimization of the average squared MMD distance between our estimateP(Y | X = x) and the truth P(Y | X = x): Theorem 1 For any split of a parent node P into child nodes C L and C R , letP split (x) = j∈{L,R} 1(x ∈ C j ) 1 n j x i ∈C j δ y i denote the resulting estimate of the distribution P(Y | X = x) when x ∈ P . Then the MMD splitting criterion can be viewed as the version of the CART criterion (7) on H: Moreover, for any node P and any fixed distributional estimatorP(Y | X = x), we have: In conclusion, from the above results we see that by applying the kernel embedding (14), we can shift the perspective to the RKHS H and view DRF as the analogue of the original Random Forest for estimation of the CME µ(P(Y | X = x)) = E[ϕ(Y) | X = x] in an abstract Hilbert space H. Like some traditional CME estimators (Song et al., 2009(Song et al., , 2013Muandet et al., 2017;Park and Muandet, 2020) it is also of the form It can be shown that, since DRF produces nonnegative weights that sum to one in (16), there is one-to-one correspondence between the resulting estimate in H and the empirical probability distribution in M b (R d ). Thus DRF can be seen as a CME estimator through (16), or directly as an estimator for the conditional distribution through (5). By contrast, other CME estimators of the form (16) have weights that are unconstrained and can be negative. Finding an appropriate distribution on R d for a given mean-embedding (sometimes referred to as "distributional inverse image problem", see e.g. McCalman et al. (2013); Muandet et al. (2017)) is not straightforward in general. For certain tasks, such as sampling from the estimated conditional distributionP(Y | X = x) or obtaining the plug-in estimates of some target functionals τ (P(Y | X = x)), this is crucial.

Convergence of Conditional Distribution Estimates
As we have seen, DRF can be viewed as the abstract version of the standard Random Forest when the response takes value in an RKHS. In principle, one could thus derive properties of DRF by adapting any existing theoretical result from the literature to the RKHS case. However, a lot of care is needed for making the results rigorous in this abstract setup, as many useful properties of R need not hold for infinite-dimensional H. This section is inspired by the results from Wager and Athey (2018).
We suppose that the forest construction satisfies the following properties, which significantly facilitate the theoretical considerations of the method and ensure that our forest estimator is well behaved, as stated in Wager and Athey (2018): (P1) (Data sampling) The bootstrap sampling with replacement, usually used in forest-based methods, is replaced by a subsampling step, where for each tree we choose a random subset of size s n out of n training data points. We consider s n going to infinity with n, with the rate specified below.
(P2) (Honesty) The data used for constructing each tree is split into two partOn the pitfalls of Gaussian scoring for causal discoverys; the first is used for determining the splits and the second for populating the leaves and thus for estimating the response.
(P3) (α-regularity) Each split leaves at least a fraction 0 < α ≤ 0.2 of the available training sample on each side. Moreover, the trees are grown until every leaf contains between κ and 2κ − 1 observations, for some fixed tuning parameter κ ∈ N.
(P4) (Symmetry) The (randomized) output of a tree does not depend on the ordering of the training samples.
(P5) (Random-split) At every split point, the probability that the split occurs along the feature X j is bounded below by π/p, for some π > 0 and for all j = 1, . . . , p.
The validity of the above properties are easily ensured by the forest construction used. From Equation (15), the prediction of DRF for a given test point x can be viewed as an element of H. If we denote the i-th training observation by Z i = (x i , µ(δ y i )) ∈ R p × H, then by (15) we estimate the embedding of the true conditional distribution µ(P(Y | X = x)) by the average of the corresponding estimates per tree: where Z k is a random subset of {Z i } n i=1 of size s n chosen for constructing the j-th tree T j and ε j is a random variable capturing all randomness in growing T j , such as the choice of the splitting candidates. T (x; ε, Z) denotes the output of a single tree: i.e. the average of the terms µ(δ Y i ) over all data points Z i contained in the leaf L(x) of the tree constructed from ε and Z.
Since one can take the number of trees N to be arbitrarily large, we consider an "idealized" version of our estimator, as done in Wager and Athey (2017), which we denote asμ n (x): where the sum is taken over all n sn possible subsets of {Z i } n i=1 .We have that µ(P(Y | X = x)) →μ n (x) as N → ∞, while keeping the other variables constant, and thus we assume for simplicity that those two quantities are the same.
Our main result shows that, under similar assumptions as in Wager and Athey (2017), the embedding of our conditional distribution estimatorμ n (x) = µ(P(Y | X = x)) consistently estimates µ(x) := µ(P(Y | X = x)) with respect to the RKHS norm with a certain rate: Theorem 2 Suppose that our forest construction satisfies properties (P1)-(P5). Assume additionally that k is a bounded and continuous kernel and that we have a random design with X 1 , . . . , X n independent and identically distributed on [0, 1] p with a density bounded away from 0 and infinity. If the subsample size s n is of order n β for some 0 < β < 1, the mapping Remark 3 The rate in (18) is analogous to the one from Wager and Athey (2018), who used it further to derive the asymptotic normality of the Random Forest estimator in R.
Unfortunately, this alone is not enough to establish asymptotic normality of (μ n (x)−µ(x))/σ n as an element of H. To do so, one needs to prove a functional central limit theorem with a Gaussian limiting process in the Hilbert space H. This then allows to deduce asymptotic normality of smooth real-valued functionals. We will provide the detailed derivations in future work.

Convergence of the Induced Estimates
The above result shows that DRF estimateP(Y | X = x) converges fast to the truth P(Y | X = x) in the MMD distance, i.e. the RKHS distance between the corresponding embeddings. Even though this is interesting on its own, ultimately we want to relate this result to estimation of certain distributional targets τ (x) = τ (P(Y | X = x)).
For any f ∈ H, we have that the DRF estimate of the target where we recall the weighting function w x (·) induced by the forest (4). Therefore, the consistency result (18) in Theorem 2 directly implies that i.e. that the DRF consistently estimates the targets of the form τ ( (18) we also obtain the rate of convergence when s n n β : for γ as in Theorem 2. When k is continuous, it is well known that all elements of H are continuous, see e.g. Hsing and Eubank (2015). Under certain assumptions on the kernel and its input space, holding for several popular kernels, (e.g. the Gaussian kernel) (Sriperumbudur, 2016), we can generalize the convergence result (19) to any bounded and continuous function f : R d → R, as the convergence of measuresP(Y | X = x) → P(Y | X = x) in the MMD metric will also imply their weak convergence, i.e. k metrizes weak convergence (Sriperumbudur, 2016;Simon-Gabriel and Schölkopf, 2018;Simon-Gabriel et al., 2020): Corollary 4 Assume that one of the following two sets of conditions holds: (a) The kernel k is bounded, (jointly) continuous and has Moreover, y → k(y 0 , y) is vanishing at infinity, for all y 0 ∈ R d .
(b) The kernel k is bounded, shift-invariant, (jointly) continuous and ν in the Bochner representation in (11) is supported on all of R d . Moreover, Y takes its values almost surely in a closed and bounded subset of R d .
Then, under the conditions of Theorem 2, we have for any bounded and continuous function f : Recalling the Portmanteau Lemma on separable metric spaces, see e.g. Dudley (2002, Chapter 11), this has several other interesting consequences, such as the consistency of CDF and quantile estimates; Let F Y | X = x (·) be the conditional CDF of Y and for any index be the corresponding DRF estimates via weighting function (6). Then we have the following result: Corollary 5 Under the conditions of Corollary 4, for any 1 ≤ i ≤ d, we havê

Applications and Numerical Experiments
The goal of this section is to demonstrate the versatility and applicability of DRF for many practical problems. We show that DRF can be used not only as an estimator of the multivariate conditional distribution, but also as a two-step method to easily obtain out-of-the box estimators for various, and potentially complex, targets τ (x).
Our main focus lies on the more complicated targets which cannot be that straightforwardly approached by conventional methods. However, we also illustrate the usage of DRF for certain applications for which there already exist several well-established methods. Whenever possible in such cases, we compare the performance of DRF with the specialized, task-specific methods to show that, despite its generality, there is at most a very small loss of precision. However, we should point out that for many targets such as, that can not be written in a form of a conditional mean or a conditional quantile, for example, conditional correlation, direct comparison of the accuracy is not possible for real data, since no suitable loss function exists and the ground truth is unknown. Finally, we show that, in addition to directly estimating certain targets, DRF can also be a very useful tool for many different applications, such as causality and fairness.
Detailed descriptions of all competing methods, data sets and the corresponding analyses can be found in Appendix C, and some additional simulations can be found in the Appendix D.

Estimation of Conditional Multivariate Distributions
In order to provide good estimates for any target τ (x) = τ (P(Y | X = x)), our method needs to estimate the conditional multivariate distribution P(Y | X = x) well. Therefore, we first investigate here the accuracy of the DRF estimate (5) of the full conditional distribution and compare its performance with the performance of several existing methods. In addition to a few simple methods such as the k-nearest neighbors or the kernel regression, which locally weight the training points, we also consider the CME estimator of Park and Muandet (2020) and several advanced machine learning methods such as the Conditional Generative Adversarial Network (CGAN) (Mirza and Osindero, 2014;Aggarwal et al., 2019), Conditional Variational Autoencoder (CVAE) (Sohn et al., 2015) and Masked Autoregressive Flow (Papamakarios et al., 2017). It is worth mentioning that the focus in the machine learning literature has been more on applications where d is very large (e.g. pixels of an image) and p is very small (such as image labels). Even though some methods do not provide the estimated conditional distribution in a form as simple as DRF, one is still able to sample from the estimated distribution and thus perform any subsequent analysis and make fair comparisons between the methods. For the CME estimator we simply set the negative weights to zero and renormalize, such that the weights are nonnegative and sum to one.
We first illustrate the estimated distributions by the above methods on a toy example where n = 1000, p = 10, d = 2 and (21) In the above example X 1 affects the mean of Y 1 , whereas X 2 affects the both mean and variance of Y 2 , and X 3 , . . . , X p have no impact. The results can be seen in Figure 3. We see j u r a s l u m  Table 1: NLPD loss computed on out-of-sample observations for the estimated conditional distributions obtained by several different methods (corresponding to rows) for many real data sets (corresponding to columns). The best method is indicated in bold.
that, unlike some other methods, DRF is able to balance the importance of the predictors X 1 and X 2 and thus to estimate the distributions of Y 1 and Y 2 well.
One can do a more extensive comparison on a collection of real data sets. We use the benchmark data sets from the multi-target regression literature (Tsoumakas et al., 2011) together with some additional ones created from the data sets described throughout this paper. The performance of DRF is compared with the performance of other existing methods for nonparametric estimation of multivariate distributions by using the Negative Log Predictive Density (NLPD) loss, which evaluates the logarithm of the induced multivariate density estimate (Quinonero-Candela et al., 2005). As the number of test points grows to infinity, NLPD loss becomes equivalent to the average KL divergence between the estimated and the true conditional distribution and is thus able to capture how well one estimates the whole distribution, instead of only its mean.
In addition to the methods mentioned above, we also include some methods that are intended only for mean prediction, by assuming that the distribution of the response around its mean is homogeneous, i.e. that the conditional distribution P (Y − E[Y | X] | X = x) does not depend on x. This is fitted by regressing each component of Y separately on X and using the pooled residuals. We consider the standard nonparametric regression methods such as Random Forest (Breiman, 2001), XGBoost (Chen and Guestrin, 2016), and Deep Neural Networks (Goodfellow et al., 2016).
The results are shown in Table 1. We see that DRF performs well for a wide range of sample size and problem dimensionality, especially in problems where p is large and d is moderately big. It does so without the need for any tuning or involved numerical optimization.

Estimation of Statistical Functionals
Figure 4: Scatter plot of predictions of the 0.1, 0.5 and 0.9 quantiles against X 1 for randomly generated 500 test data points X test ∼ U (−1, 1) p . The true values of the quantiles are displayed by black dashed lines. The columns corresponds to different methods DRF (red), GRF (green), QRF (blue), TRF (purple). The rows correspond to different simulation scenarios. The first two are taken from .
We first investigate the performance for the classical problem of univariate quantile estimation on simulated data. We consider the following three data generating mechanisms with p = 40, n = 2000 and X i i.i.d. (1)   In Figure 4 we can see the corresponding estimates of the conditional quantiles for DRF, Quantile Regression Forest (QRF) (Meinshausen, 2006), which uses the same forest construction with CART splitting criterion as the original Random Forest (Breiman, 2001) but estimates the quantiles from the induced weighting function, Generalized Random Forests (GRF)  with a splitting criterion specifically designed for quantile estimation and Transformation Forests (TRF) (Hothorn and Zeileis, 2021). We see that DRF is performing very well even compared to methods that are specifically tailored to quantile estimation.
The multivariate setting is however more interesting, as one can use DRF to compute much more interesting statistical functionals τ (x). We illustrate this in Figure 5 for the air quality data set, described in Section 2.2. The left plot shows one value of the estimated multivariate CDF, specifically the estimated probability of the event that the air quality index (AQI) is at most 50 at a given test site. This corresponds to the "Good" category and means that the amount of every air pollutant is below a certain threshold determined by the EPA. Such probability estimates can be easily obtained by summing the weights of the training points belonging to the event of interest. For both plots in Figures 2 and 5, we train the single DRF with the same set of predictor variables and take the three pollutants O 3 , SO 2 and PM2.5 as the responses. In this way we still have training data from many different sites.
In order to investigate the accuracy of the conditional CDF obtained by DRF, we compare the estimated probabilities with estimates of the standard univariate classification forest (Breiman, 2001) with the response 1(AQI ≤ 50). In the left plot of Figure 6, we can see that the DRF estimates of the P(AQI ≤ 50 | X = x) (also visualized in Figure 5) are quite similar to the estimates of the classification forest predicting the outcome 1(AQI ≤ 50). Furthermore, the cross-entropy loss evaluated on the held-out measurements equals 0.4671 and 0.4663 respectively, showing almost no loss of precision. In general, estimating the simple functionals from the weights provided by DRF comes usually at a small to no loss compared to the classical methods specifically designed for this task.
In addition to the classical functionals τ (x) in the form of an expectation for some function f : R d → R, which can also be computed by solving the corresponding one-dimensional problems, additional interesting statistical  Figure 5) and by the classification forest. Right: Example how the CDF estimated by using the classification forest (blue) need not be monotone, whereas the DRF estimates (red) are well-behaved.
functionals with intrinsically multivariate nature that are not that simple to estimate directly are accessible by DRF, such as, for example, the conditional correlations Cor As an illustration, the estimated correlation of the sulfur dioxide (SO 2 ) and fine particulate matter (PM2.5) is shown in the right plot of Figure 5. The plot reveals also that the local correlation in many big cities is slightly larger than in its surroundings, which can be explained by the fact that the industrial production directly affects the levels of both pollutants.
A big advantage of the target-free forest construction of DRF is that all subsequent targets are computed from same the weighting function w x obtained from a single forest fit. First, this is computationally more efficient, since we do not need for every target of interest to fit the method specifically tailored to it. For example, estimating the CDF with classification forests requires fitting one forest for each function value. Secondly and even more importantly, since all statistical functionals are plug-in estimates computed from the same weighting function, the obtained estimates are mathematically well-behaved and mutually compatible. For example, if we estimate Cor , and Var(Y j | X = x), one can not in general guarantee the estimate to be in the range [−1, 1], but this is possible with DRF. Alternatively, the correlation or covariance matrices that are estimated entrywise are guaranteed to be positive semi-definite if one uses DRF. As an additional illustration, Figure 6 shows that the estimated (univariate) CDF using the classification forest need not be monotone due to random errors in each predicted value, which can not happen with the DRF estimates.

Conditional Copulas and Conditional Independence Testing
One can use the weighting function not only to estimate certain functionals, but also to obtain more complex objects, such as, for example, the conditional copulas. The well-known Sklar's theorem (Sklar, 1959)  for 1 ≤ i ≤ d, as follows: Copulas capture the dependence of the components Y i by the joint distribution of the corresponding quantile levels of the marginal distributions: Decomposing the full multivariate distribution to marginal distributions and the copula is a very useful technique used in many fields such as risk analysis or finance (Cherubini et al., 2004). Using DRF enables us to estimate copulas conditionally, either by fitting certain parametric model or nonparametrically, directly from the weights.
To illustrate this, consider an example where the 5-dimensional Y is generated from the equicorrelated Gaussian copula ∼ U (0, 1) p , where p = 30 and n = 5000. All Y i have a N (0, 1) distribution marginally, but their conditional correlation for i = j is given by Figure 7 shows that DRF estimates the full conditional distribution at different test points x quite accurately and thus we can obtain a good nonparametric estimate of the conditional copula as follows. First, for each component Y i , we compute the corresponding marginal CDF estimateF Y i | X = x (·) from the weights. Second, we map Figure 7 in order to visualize the copula.
If we want to instead estimate the copula parametrically, we need to find the choice of parameters for a given model family which best matches the estimated conditional distribution, e.g. by weighted maximum likelihood estimation (MLE). For the above example, the correlation parameter of the Gaussian copula can be estimated by computing the weighted correlation with weights {w x (x i )} n i=1 . The left plot in Figure 8 shows the resulting estimates of the conditional correlation Cor (Y 1 , Y 2 | X = x) obtained from DRF MMD , which uses the MMD splitting criterion (12) described in Section 2.3.1, and DRF CART , which aggregates the marginal CART criteria (Kocev et al., 2007;Segal and Xiao, 2011). We see that DRF MMD is able to detect the distributional heterogeneity and provide good estimates of the conditional correlation. On the other hand, DRF CART cannot detect the change in distribution of Y caused by X 1 that well. The distributional heterogeneity can not only occur in marginal distribution of the responses (a case extensively studied in the literature), but also in their interdependence structure described by the conditional copula C x , as one can see from decomposition (22). Since DRF MMD relies on a distributional metric for its splitting criterion, it is capable of detecting any change in distribution (Gretton et al., 2007a), whereas aggregating marginal CART criteria for Y 1 , . . . , Y d in DRF CART only captures the changes in the marginal means.
This is further illustrated for a related application of conditional independence testing, where we compute some dependence measure from the obtained weights. For example, we can test the independence Y 1 ⊥ ⊥ Y 2 conditionally on the event X = x by using the Hilbert Schmidt Independence Criterion (HSIC) (Gretton et al., 2007b), which measures the difference between the joint distribution and the product of the marginal distributions. The right plot of Figure 8 shows that the DRF MMD estimates are quite close to the population value of the HSIC, unlike the ones obtained by DRF CART .

Heterogeneous Regression and Causal Effect Estimation
In this and the following section, we illustrate that, in addition to direct estimation of certain targets, DRF can also be a useful tool for complex statistical problems and applications, such as causality.
Suppose we would like to investigate the relationship between some (univariate) quantity of interest Y and certain predictors W from heterogeneous data, where the change in distribution of (W, Y ) can be explained by some other covariates X. Very often in causality applications, W is a (multivariate) treatment variable, Y is the outcome, which is commonly, but not necessarily, binary, and X is a set of observed confounding variables for which we need to adjust if we are interested in the causal effect of W on Y . This is illustrated by the following causal graph: The problem of nonparametric confounding adjustment is hard; not only can the marginal distributions of Y and W be affected by X, thus inducing spurious associations due to confounding, but the way how W affects Y can itself depend on X, i.e. the treatment effect might be heterogeneous. The total causal effect can be computed by using the adjustment formula (Pearl, 2009): In general, implementing do-calculus for finite samples and potentially non-discrete data might not be straightforward and comes with certain difficulties. In this case, the standard approach would be to estimate the conditional mean E[Y | W = w, X = x] nonparametrically by regressing Y on (X, W) with some method of choice and to average out the estimates over different x sampled from the observed distribution of X. Using DRF for this approach is not necessary, but has an advantage that one can easily estimate the full interventional distribution P(Y | do(W = w)) and not only the interventional mean E[Y | do(W = w)].
Another way of computing the causal effect, which allows to add more structure to the problem, is explained in the following: We use DRF to first fit the forest with the multivariate response (W, Y ) and the predictors X. In this way, one can for any point of interest x obtain the joint distribution of (W, Y ) conditionally on the event X = x and then the weights {w x (x i )} n i=1 can be used as an input for some regression method for regressing Y on W in the second step. This conditional regression fit might be of an independent interest, but it can also be used for estimating the causal effect E[Y | do(W = w)] from (23), by averaging the estimates E[Y | W = w, X = x] over x, where x is sampled from the empirical observation of X. In this way one can efficiently exploit and incorporate any prior knowledge of the relationship between W and Y , such as, for example, monotonicity, smoothness or that it satisfies a certain parametric regression model, without imposing any assumptions on the effect of X on (W, Y ). Furthermore, one might be able to better extrapolate to the regions of space where P(W = w, X = x) is small, compared to the standard approach which computes E[Y | W = w, X = x] directly, by regressing Y on (W, X). Extrapolation is crucial for causal applications, since for computing E[Y | do(W = w)] we are interested in what would happen with Y when our treatment variable W is set to be w, regardless of the value achieved by X. However, it can easily happen that for this specific combination of X and W there are very few observed data points, thus making the estimation of the causal effect hard (Pearl, 2009).
As an illustration, we consider the following synthetic data example, with continuous outcome Y , continuous univariate treatment W , n = 5000 and p = 20: A visualization of the data can be seen on the left side of Figure 9; treatment W affects Y nonlinearly, X 2 is a confounding variable that affects the marginal distributions of Y and W and X 1 makes the treatment effect heterogeneous. The middle part of Figure 9 shows the conditional regression fits, i.e. the estimates of E[Y | W = w, X = x] as w varies and x is fixed. In general, the conditional regression fit is related to the concept of the conditional average treatment effect (CATE) as it quantifies the effect of W on Y for the subpopulation for which X = x. We see that combination of DRF with response (Y, W ) and predictors X with the smoothing splines regression of Y on W (blue curve) is more accurate than the estimates obtained by standard Random Forest (Breiman, 2001) with response Y and predictors (W, X) (red curve). Furthermore, we see that the former approach can extrapolate better to regions with small number of data points, which enables us to better estimate the causal effect E[Y | do(W = w)] from (23), by averaging the corresponding estimates of E[Y | W = w, X = x] over observed x, as shown in the right plot of Figure 9. There exist many successful methods in the literature for estimating the causal effects and the (conditional) average treatment effects for a wide range of settings (Abadie and Imbens, 2006;Chernozhukov et al., 2018;Wager and Athey, 2018;Künzel et al., 2019). However, some methods are not designed for the most general case and make certain modeling assumptions or are designed specifically for the (very common) case where the treatment variable is univariate or even binary. Due to its versatility, DRF can easily be used when the underlying assumptions of conventional methods are violated, when some additional structure is given in the problem or for the general, nonparametric, settings (Imbens, 2004;Ernest and Bühlmann, 2015;Kennedy et al., 2017). Appendix D contains additional comparisons with some existing methods for causal effect estimation.

Births data
We further illustrate the applicability of DRF for causality-related problems on the natality data obtained from the Centers for Disease Control and Prevention (CDC) website, where we have information about all recorded births in the USA in 2018. We investigate the relationship between the pregnancy length and the birthweight, an important indicator of baby's health. Not only is this relationship complex, but it also depends on many different factors, such as parents' race, baby's gender, birth multiplicity (single, twins, triplets...) etc.
In the left two plots of Figure 10 one can see the estimated joint distribution of birthweight and pregnancy length conditionally on many different covariates, as indicated in the plot. The black curves denote the subsequent regression fit, based on smoothing splines. In addition to the estimate of the mean, indicated by the solid curve, we also include the estimates of the conditional 0.1-and 0.9-quantiles, indicated by dashed curves, which is very useful in practice for determining whether a baby is large or small for its gestational age. Notice how DRF assigns less importance to the mother's race when the point of interest is a twin (middle plot), as in this case more weight is given to twin births, regardless of the race of the parents.
Suppose now we would like to understand how a twin birth T causally affects the birthweight B, but ignoring the obvious indirect effect due to shorter pregnancy length L. For example, sharing of resources between the babies might have some effect on their birthweight. We additionally need to be careful to adjust for other confounding variables X, such as, for example, the parents' race, which can affect B, T and L. We assume that this is represented by the following causal graph: In order to answer the above question, we investigate the causal quantity P(B | do(T = t, L = l)). Even though one cannot make such do-intervention in practice, this quantity describes the total causal effect if the birth multiplicity and the length of the pregnancy could be manipulated and thus for a fixed pregnancy length l, we can see the difference in birthweight due to T . We compute this quantity as above, by using DRF with subsequent regression fits, which has the advantage of better extrapolating to regions with small probability, such as long twin pregnancies (see the middle plot of Figure 10). In the right plot of Figure 10 we show the mean and quantiles of the estimated interventional distribution and we see that, as one might expect, a twin birth causes smaller birthweight on average, with the difference increasing with the length of the pregnancy.

Fairness
Being able to compute different causal quantities with DRF could prove useful in a range of applications, including fairness (Kusner et al., 2017). We investigate the data on approximately 1 million full-time employees from the 2018 American Community Survey by the US Census Bureau from which we have extracted the salary information and all covariates that might be relevant for salaries. In the bottom left plot of Figure 11 one can see the distribution of hourly salary of men and women (on the logarithmic scale). The overall salary was scaled with working hours to account for working part-time and for the fact that certain jobs have different working hours. We can see that men are paid more in general, especially for the very high salaries. The difference between the median hourly salaries, a commonly used statistic in practice, amounts 17% for this data set.
We would like to answer whether the observed gender pay gap in the data is indeed unfair, i.e. only due to the gender, or whether it can at least in part be explained by some other factors, such as age, job type, number of children, geography, race, attained education level and many others. Hypothetically, it could be, for example, that women have a preference for jobs that are paid less, thus causing the gender pay gap.
In order to answer this question, we assume that the data is obtained from the following causal graph, where G denotes the gender, W the hourly wage and all other factors are denoted by X: G X W i.e. G is a source node and W is a sink node in the graph. In order to determine the direct effect of the gender on wage that is not mediated by other factors, we would like to compute the distribution of the nested counterfactual W (male, X(female)), which is interpreted as the women's wage had they been treated in same way as men by their employers for determining the salary, but without changing their propensities for other characteristics, such as the choice of occupation (Chernozhukov et al., 2013). Therefore, it can be obtained from the observed distribution as follows: Put in the language of the fairness literature, it quantifies the unfairness when all variables X are assumed to be resolving (Kilbertus et al., 2017), meaning that any difference in salaries directly due to factors X is not viewed as gender discrimination. For example, one does not consider unfair if people with low education level get lower salaries, even if the gender distribution in this group is not balanced. There are several ways how one can compute the distribution of W (male, X(female)) from (25) with DRF. The most straightforward option is to take W as the response and (G, X) as predictors in order to compute the conditional distribution P (W | G = male, X = x). However, with this approach it could happen that for predicting P (W | G = male, X = x) we also assign weight to training data points for which G = female. This happens if in some trees we did not split on variable G, which is likely, for example, if P(G = male | X = x) is low. Using salaries of both genders to estimate the distribution of men's salaries might be an issue if our goal is to objectively compare how women and men are paid.
Another approach is to take (W, G) as a multivariate response and X as the predictors for DRF and thus obtain joint distribution of (W, G) conditionally on the event X = x. In this way we can also quantify the gender discrimination of a single individual with characteristics x by comparing his/her salary to the corresponding quantile of the salary distribution of people of the opposite gender with the same characteristics x (Plečko and Meinshausen, 2020). This is interesting because the distribution of salaries, and thus also the gender discrimination, can be quite different depending on other factors such as the industry sector or job type, as illustrated for a few choices of x in the top row of Figure 11.
Finally, by averaging the DRF estimates of P (W | X = x, G = male), conveniently represented via the weights, over different x sampled from the distribution P(X | G = female), we can compute the distribution of the nested counterfactual W (male, X(female)) (Chernozhukov et al., 2013). In the middle panel in the bottom row of Figure 11 a noticeable difference in the means, also called natural direct effect in the causality literature (Pearl, 2009), is still visible between the observed distribution of women's salaries and the hypothetical distribution of their salaries had they been treated as men, despite adjusting for indirect effects of the gender via covariates X. By further matching the quantiles of the counterfactual distribution P (W (male, X(female))) with the corresponding quantiles of the observed distribution of men's salaries in the bottom right panel of Figure 11, we can also see that the adjusted gender pay gap even increases for larger salaries. Median hourly wage for women is still 11% lower than the median wage for the hypothetical population of men with exactly the same characteristics X as women, indicating that only a minor proportion of the actually observed hourly wage difference of 17% can be explained by other demographic factors.

Conclusion
We have shown that DRF is a flexible, general and powerful tool, which exploits the well-known properties of the Random Forest as an adaptive nearest neighbor method via the induced weighting function. Not only does it estimate multivariate conditional distributions well, but it constructs the forest in a model-and target-free way and is thus an easy to use out-of-the-box algorithm for many, potentially complex, learning problems in a wide range of applications, including also causality and fairness, with competitive performance even for problems with existing tailored methods.

A. Implementation Details
Here we present the implementation of the Distributional Random Forests (DRF) in detail. The code is available in the R-package drf and the Python package drf. The implementation is based on the implementations of the R-packages grf  and ranger (Wright and Ziegler, 2017). The largest difference is in the splitting criterion itself and the provided user interface. Algorithm 1 gives the pseudocode for the forest construction and computation of the weighting function w x (·).
• Every tree is constructed based on a random subset of size s (taken to be 50% of the size of the training set by default) of the training data set, similar to Wager and Athey (2018). This differs from the original Random Forest algorithm (Breiman, 2001), where the bootstrap subsampling is done by drawing from the original sample with replacement.
• The principle of honesty (Biau, 2012;Denil et al., 2014;Wager and Athey, 2018) is used for building the trees (line 4), whereby for each tree one first performs the splitting based on one random set of data points S build , and then populates the leaves with a disjoint random set S populate of data points for determining the weighting function w x (·). This prevents overfitting, since we do not assign weight to the data points which we used to built the tree.
• We borrow the method for selecting the number of candidate splitting variables from the grf package . This number is randomly generated as min(max(Poisson(mtry), 1), p), where mtry is a tuning parameter. This differs from the original Random Forests algorithm, where the number of splitting candidates is fixed to be mtry.
• The number of trees built is N = 2000 by default.
• The factor variables in both the responses and the predictors are encoded by using the one-hot encoding, where we add an additional indicator variable for each level l of some factor variable X k . This implies that in the building step, if we split on this indicator variable, we divide the current set of data points in the sets where X k = l and X k = l. This works well if the number of levels is not too big, since otherwise one makes very uneven splits and the dimensionality of the problem increases significantly.
Handling of categorical problems is a general challenge for the forest based methods and is an area of active research (Johannemann et al., 2019). We will leave improving on this approach for the future development.
• We try to enforce splits where each child has at least a fixed percentage (chosen to be 10% as the default value) of the current number of data points. In this way we achieve balanced splits and reduce the computational time. However, we cannot enforce this if we are trying to split on the variable X i with only a few unique values, e.g. indicator variable for a level of some factor variable.
• All components of the response Y are scaled for the building step (but not when we populate the leaves). This ensures that each component of the response contributes return w 40: end procedure equally to the kernel values, and consequently to the MMD two-sample test statistic. Plain usage of the MMD two-sample test would scale the components of Y at each node. However, this approach favors always splitting on the same variables, even though their effect will diminish significantly after having split several times.
• By default, in step 20 of the Algorithm 1, we use the MMD-based splitting criterion given by is used as the default choice, with the bandwidth σ chosen as the median pairwise distance between all training responses {y i } n i=1 , commonly referred to as the 'median heuristic' (Gretton et al., 2012c). However the algorithm can be used with any choice of kernel, or in fact with any two-sample test.
• The number B of random Fourier features is fixed and taken to be 20 by default. The performance of the trees empirically shows stability for large range of B. Smaller values of B help making the trees more independent which could improve the performance. One could even use an adaptive strategy of choosing B, possibly increasing B as the depth of the tree increases, but we decided to keep B fixed for simplicity.
• We compute variable importance similarly as for the original Random Forest algorithm (Breiman, 2001;Wright and Ziegler, 2017), by sequentially permuting each variable and investigating the decrease in the performance. However, since we target the full conditional distribution of the multivariate response, as our performance measure we use for every test point (x, y) the MMD distance between the estimated joint distributionP(Y | X = x), described by the DRF weights, and the point mass δ y .

B. Derivations and Proofs
In this section we present proofs and further details for the results in Sections 2.3.1 and 3, in the order in which they appear. Sections B.1 and B.3 derive the splitting criterion presented in Equation (12) in the main text, with Section B.3 showing that the CART criterion can be written in an analogous way. Section B.4 provides background and proves to the statements in Section 3.

B.1 Expressing MMD test statistic as an integral in the feature space
The biased MMD two-sample statistic is given as Assume that the kernel k is bounded and shift-invariant, then by Bochner's theorem there exist a measure ν such that k can be written as k( where ϕ ω (y) = e iω T y ∈ C are the corresponding Fourier features, which is what we wanted to show.

B.2 Approximate kernel and its MMD
When the kernel k is bounded and shift invariant, we have seen that it can be written as k(x, y) = R d e iω T (x−y) dν(ω). This integral can be approximated by sampling from ν: Let ω 1 , . . . , ω B ∼ ν be a random sample from the measure ν. Then we can write where ϕ(u) = (ϕ ω 1 (u), . . . , ϕ ω B (u)) T is a random complex vector consisting of the Fourier features ϕ ω (u) = e iω T u ∈ C. The kernelk is analogous to the kernel k, but where the measure ν is replaced by the empirical measureν = 1 Analogously as in the section B.1, we can now write the MMD for the kernelk as: which can also additionally be interpreted as the approximation of D MMD . Therefore, our splitting criterion is obtained as the MMD of the random approximate kernelk: The scaling factor n L n R n 2 P occurs naturally and penalizes the increased variance of the sample MMD statistic when n L or n R are small: it appears when we rewrite the CART criterion in the related form, see section B.3.
This representation of the MMD is the key why we use the approximate kernelk instead of k. This splitting criterion can be computed in O(Bn P ) complexity, by updating the sums (1) computations, whereas this is not possible for D MMD .

B.3 CART criterion rewritten
Standard CART criterion used in Random Forests (Breiman, 2001) is the following: we repeatedly choose to split the parent node P of size n P in two children C L and C R , of sizes n L and n R respectively, such that the expression Similarly we obtain which gives us that the CART criterion (26) can be written as since n L + n R = n P . Since the first term depends only on the parent node and not on the chosen split, we conclude that minimizing the CART criterion (26) is equivalent to maximizing the following expression This equivalent criterion can be interpreted as comparing the difference in the means of the resulting child nodes, i.e. we will choose the split such that the means in the child nodes are as heterogeneous as possible. The scaling factor n L n R n 2 P appears naturally, penalizing uneven splits due to the increased variance of Y L or Y R .

B.4 Proofs for Section 3
Preliminaries. We first set notation and define basic probabilistic concepts on the separable Hilbert space (H, , ·, H ). We thereby mostly refer to Hsing and Eubank (2015) and Pisier (2016). The initial results derived here parallel some of the results derived in Park and Muandet (2020), but where derived independently. Let (Ω, A, P) be the underlying probability space. Let (H, ·, · H ) be the Hilbert space induced by the kernel k and µ : M b (R d ) → H be the embedding function of, µ(P ) ∈ H for all bounded signed Borel measures P on R d . Throughout we assume that k is bounded and continuous in its two arguments. Boundedness of k ensures that µ is indeed defined on all of M b (R d ), while continuity of k : R d × R d → R ensures H is separable. Thus measurability issues can be avoided, in particular, a map ξ : (Ω, A) → (H, B(H)) is measurable iff ξ, f H is measurable for all f ∈ H. Moreover, a quick check reveals that µ(P ) is linear on where the integral is meant in a Bochner sense. Separability and E[ ξ H ] < ∞ mean this integral is well-defined and moreover for any continuous linear function F : Define moreover for q ≥ 1, and ξ, ξ 1 , ξ 2 ∈ L 2 (Ω, A, H), 1. Here and later F (ξ) is meant to mean F (ξ(ω)) for all ω ∈ Ω.
It is well-known, that (L q , · L q (H) ) is a Banach space, with We can then also define conditional expectation. For a sub σ− algebra F ⊂ A, ξ ∈ L 1 (Ω, A, H), E[ξ | F] is the (a.s.) unique element such that See e.g. Umegaki and Bharucha-Reid (1970) again we refer to (Umegaki and Bharucha-Reid, 1970). We note that, as with conditional expectation on R, E[ξ | F] is only defined uniquely a.s. As such all (in)equalitie statements hold only a.s. However, we will often not explicitly write this going forward. We The following Proposition shows that this notion is well-defined and some further properties of Hilbert space-valued conditional expectation, in addition to (C1) -(C3): Proposition 6 Let (H 1 , ·, · 1 ), (H 2 , ·, · 2 ) be two separable Hilbert spaces, X, X 1 , X 2 ∈ L 1 (Ω, A, H 1 ) and ξ 1 , ξ 2 , ξ ∈ L 1 (Ω, A, H 2 ). 2 (C4) There exists a measurable function h : (H 1 , B(H Proof We will use the following fact in the proof: Under the assumption of separability, all relevant notions of measurability are the same, see e.g. Hsing and Eubank (2015, Chapter 2, 7). In particular, ξ ∈ L q (Ω, F, H 2 ), for any q ≥ 1, means that ξ : (Ω, F) → (H 2 , B(H 2 )) is measurable, which in turns means there exists a sequence of simple functions with g k ∈ H 2 and A k ∈ F for all k and such that f n → ξ a.s. on H 2 and even f n −ξ L q (H) → 0, see e.g. Pisier (2016, Proposition 1.2).
2. We again note that all equalities technically only hold a.s.
We first show (C5) for simple functions and then extend this to L 2 (Ω, σ(X), H 1 ), using the fact at the beginning of the proof. Let thus f n be of the form (28). Then for all F ∈ σ(X), from the properties of real-valued conditional expectation. As additionally E[f n | X], ξ 2 H is clearly σ(X) measurable, (C1) and (C2) are met for this candidate. Since conditional expectation is (a.s.) uniquely defined by (C1) and (C2), (C5) holds true for the special case of simple functions. For general ξ 1 ∈ L 2 (Ω, A, H 1 ), let f n have f n − ξ 1 L 2 (H) → 0. The goal is to show that We can bound both terms by the same quantity, using Cauchy-Schwarz: as the random variable 1 F ξ 2 H is σ(X) measurable by assumption and The result thus follows from the Hölder inequality, Finally (C6) and (C7) is easily proven with the technique of "scalarization" (Pisier, 2016, Chapter 1). We do the full argument for (C6), (C7) can be shown analogously. That is, we show that for any F : H 2 → R linear and continuous, for almost all ω and some representation of E[ξ | X 1 , X 2 ]. This then immediately implies the result. Now, using the property of real-valued conditional expectations , we obtain (31).
(C8) If ξ and X 2 are conditionally independent given X 1 , then E[ξ | X 1 , Proof We prove (C8) again using the scalarization trick: For any F : H 2 → R continuous and linear, it holds that is a bounded and measurable function. Thus by assumption, s. an application of Fatou's Lemma for (real-valued) conditional expectation (see e.g., Dudley (2002, Problem 10 Thus, we have shown that for all F : H 2 → R continuous and linear, F (E[ξ | X 1 , X 2 ]) = F (E[ξ | X 1 ]), proving the claim.
Finally, combining (C5), (C7) and (C8), we obtain for ξ 1 , ξ 2 conditionally independent given X Let for x ∈ X , P x be the conditional distribution of Y given x on R d and similarly with P X (i.e. the regular conditional probability measure). Note that P x ∈ H, while P X is a random element mapping into H.
As in Wager and Athey (2018), we define for f , g two functions, with lim inf s→∞ g(s) > 0, for some C > 0. If C = 1, then we write f (s) g(s). For a sequence of random variables X n : Ω → R, and a n ∈ (0, +∞), n ∈ N, we write as usual i.e. if X n is bounded in probability. We write X n = o p (a n ), if a −1 n X n converges in probability to zero. Similarly, for (S, d) a separable metric space, X n : (Ω, A) → (S, B(S)), n ∈ N and X : (Ω, A) → (S, B(S)) measurable, we write X n p → X, if d(X n , X) = o p (1). Finally let X ∈ L 2 (Ω, A, H 1 ), ξ ∈ L 2 (Ω, A, H 2 ) and assume that A ⊂ Ω depends on X, A = A(X). Thus for X fixed to a certain value, A is a fixed set. If P(A | X) > 0 almost everywhere, we define It then holds by construction that Let again µ(x) := µ(P x ) be the embedding of the true conditional distribution into H. We first state 3 preliminary results: Proof We first note that P x exists and is a probability measure on R d . Since k is bounded µ(P x ) ∈ H exists and is uniquely defined by the relation On the other hand, by the above E[µ(δ Y ) | X] ∈ L 2 (Ω, A, H) exists and for all f ∈ H, since F : H → R, F (g) = g, f H defines a continuous linear function. In particular, for all f ∈ H, Since µ(δ Y ) = k(Y, ·), Lemma 8 in fact corresponds to Lemma 3.2 in Park and Muandet (2020). For a more compact notation in the following Lemma, let N = {1, . . . , n} and let for A ⊂ N and k ≤ |A|, C k (A) be the set of all subsets of size k drawn from A without replacement, with C 0 := ∅.
Write X n = (X 1 , . . . , X n ) and let T : and it holds that and Proof Composition (33) was proven in a range of different ways for real-valued T , see e.g. Hoeffding (1961), Dynkin and Mandelbaum (1983), Efron and Stein (1981), or Vitale (1992). We consider and slightly extend the elegant proof of Dynkin and Mandelbaum (1983) to also prove (34). See also Vitale (1992). Let We note that T does not depend on the exact indices in A thanks to the assumed symmetry of T . Since E[T (X n ) | X 1 , . . . , X n ] = T (X n ) is part of T n (X n ), this leads to a telescoping sum, already proving (33). Adapting the approach of Dynkin and Mandelbaum (1983), consider now L 2 (Ω, σ(X n ), H) and let Q i be the projection operator into L 2 (Ω, σ(X 1 , . . . , X i−1 , X i+1 , . . . , X n ), H). That is by (C3). Now as in Vitale (1992), In particular, Moreover, it holds that Expanding the identity as in Vitale (1992) and using Q i (1 − Q i ) = 0, proves (33). Furthermore the above implies that for any subset A l ⊂ N that intersects with A = {i 1 , . . . , i }, we must have E[T (X i 1 , . . . , X i ) | π N A l (X n )] = 0. Indeed assume A l ∩ A = {i 2 , . . . , i }, then since the elements of X n are independent, by (C6), i.e. all elements outside the intersection are irrelevant. Moreover, and clearly this projection can only be 0. The same argument can be made for any other intersection set A l ∩ A , even if it is the empty set.
Thus combining the above with (C5), for A = A l , it holds that In other words, the covariance between any two elements in the decomposition of T (X n ) − E[T (X)] in (33) is uncorrelated. Finally then In order to proceed, we first prove Theorem 1, which is somewhat separate from the remainder of this section: Theorem 1 For any split of a parent node P into child nodes C L and C R , let n j x i ∈C j δ y i denote the resulting estimate of the distribution P(Y | X = x) when x ∈ P . Then the MMD splitting criterion can be viewed as the version of the CART criterion (7) on H: Moreover, for any node P and any fixed distributional estimatorP(Y | X = x), we have: Proof The first part of the Theorem is shown analogously to the proof in Section B.3 of the appendix, but where we replace the standard dot product in R with the inner product ·, · H associated with (H, k) and use the induced RKHS norm · H . Also, since k is bounded, the embedding µ(D) into RKHS H exists for any distribution D, so everything is well-defined.
For the second statement of the Theorem 1, note that n P ∼ Binomial(π, n), where π := P(X ∈ P ) > 0. LetP x =P(Y | X = x) be a fixed conditional distribution estimator and recall that P x = P(Y | X = x). We now write Then it holds that It follows with Lemma 8 and (C5) that, Combining the three equations, this means using that E[g(X) | X ∈ P ] = E[g(X)1{X ∈ P }]/P(X ∈ P ). We now show that the difference between 1 2 H and the expectation on the left of Equation (35) is O p (n −1/2 ), using standard CLT arguments. Define K = sup z,z | k(z, z ) | < ∞, as we have assumed that k is bounded. For any two distributions D 1 , D 2 we now obtain where Z 1 , Z 1 ∼ D 1 and Z 2 , Z 2 ∼ D 2 are independent random variables. Thus, implying that both first and second moments of the random variable . . , n are its i.i.d. realizations, it follows directly from the CLT that: By multiplying the above equation with n/n P = (1/π + o p (1)) = O p (1), it also holds that Thus, Now both terms in (37) are O p (1): For the first term this follows from (36). For the second which in turn is true by another application of the CLT on random variables 1{X i ∈ P } and the fact that n/n P = O p (1): Combining the fact that the expression in (37) is O p (1) with (35) gives the result.
Letμ n (x) be defined as in (17): where the sum is taken over all n sn possible subsamples Z i 1 , . . . , Z is n of Z 1 , . . . Z n and s n → ∞ with n. Moreover, a single tree is given as Given these preliminary results, the proofs we use are for the most part analogous to the ones in Wager and Athey (2018). Thus the proofs are given mostly for completeness and sometimes omitted altogether. We introduce the following additional notation, similar to Section 3: Let Z s = (Z 1 , . . . , Z s ) collect s i.i.d. copies of Z and define for j = 1, . . . , s n , The index s will take the role of s n or n, depending on the situation. We note that, due to i.i.d. sampling, it doesn't matter for variance or expectation what kind of subset Z i 1 , . . . , Z is we are considering. In particular, we might just take Z s each time.
Before going on to the main proofs, we repeat here the assumed properties of the trees for completeness: (P1) (Data sampling) The bootstrap sampling with replacement, usually used in forest-based methods, is replaced by a subsampling step, where for each tree we choose a random subset of size s n out of n training data points. We consider s n going to infinity with n, with the rate specified below.
(P2) (Honesty) An observation Z = (X, µ(δ Y )) is either used to place the splits in a tree or to estimate the response, but never both.
(P3) (α-regularity) Each split leaves at least a fraction α ≤ 0.2 of the available training sample on each side. Moreover, the trees are grown until every leaf contains between κ and 2κ − 1 observations, for some fixed tuning parameter κ ∈ N.
(P4) (Symmetry) The (randomized) output of a tree does not depend on the ordering of the training samples.
(P5) (Random-split) At every split point, the probability that the split occurs along the feature X j is bounded below by π/p, for some π > 0 and for all j = 1, . . . , p.
Note that assumption (P2) differentiates from assumption (P2) in the main text. We refer to it as (P2') (Double Sampling) The data used for constructing each tree is split into two parts; the first is used for determining the splits and the second for populating the leaves and thus for estimating the response.
(P2) will allow us to assume that all s n observations of the tree are used to estimate the response, as in (39) and the trees are built with some auxiliary data. This is done for simplicity of exposition, the results can be extended to hold in case of (P2') as well. In fact, since the (random) division into the two data sets can be seen as part of ε k , the adaptation simply involves changing s n to s n /2. Now, we may directly apply Lemma (9) to the U-statisticsμ n (x): Lemma 10 Letμ n (x) be as in (38) and assume T satisfies (P4) and Var(T ) < ∞.
Proof Using the composition in (33) on T (Z i 1 , . . . , Moreover, it holds by symmetry and i.i.d. sampling, that for A 1 = A 2 ⊂ N , |A 1 | = |A 2 | = , T (π A 1 (Z n )) = T (π A 2 (Z n )). Thus we obtain, where (s n ) j = s n (s n − 1) · (s n − (j − 1)) = s n !/(s n − j)!. In particular n s n −1 n − 1 Consequently, with covariances between terms equal to 0, as in Lemma 9. Thus Var(μ n (x)) = s 2 n n 2 nVar( where the last step followed from (34). This proves the first inequality. On the other hand, we have from Lemma 9 that leading to the second inequality.
Lemma 11 [ Lemma 2 from Wager and Athey (2018)] Let T is a tree satisfying (P3), (P5) trained on Z s = (ξ 1 , X 1 ), . . . (ξ s , X s ) and let L(x, Z s ) be the leaf containing x. Suppose that X 1 , . . . , X s are i.i.d. on [0, 1] p independently with a density f bounded away from 0 and infinity. Then, Lemma 12 Let T be a tree satisfying (P1), (P2) and L(x, Z s ) be the leaf containing x. Then and Var(T (Z s )) ≤ sup Proof We want to prove Let for the following N x = s i=1 1{X i ∈ L(x, Z s )}. Then due to i.i.d. sampling: The last expression can be broken into Now, by honesty, given the knowledge that X 1 ∈ L(x, Z s ), ξ 1 is independent of N x , thus: Now it holds by i.i.d. sampling that, We focus on the second term. For the first, the bound follows by analogous arguments. Similar as before, By the same argument as above Thus combining (46) and (47), where in the last step we used independence of (ξ 1 , 1{X 1 ∈ L(x, Z s )}), (ξ 2 , 1{X 2 ∈ L(x, Z s )}) conditionally on L(x, Z s ) and (C9). Finally, proving the claim.
Corollary 13 Under the conditions of Lemma 11, assume is Lipschitz and that the trees T in the forest satisfy (P2) and (P3). Then and If moreover, is Lipschitz, then: Proof By (42), it holds as in Wager and Athey (2018) π p away from 0 and infinity. If the subsample size s n is of order n β for some 0 < β < 1, the mapping we obtain the consistency w.r.t. the RKHS norm: Proof We first note that sup (43) implies Var(T ) < ∞. Thus, from Markov's inequality and Lemma 10, for γ ≤ (1 − β)/2. In particular, it goes to zero for any ε > 0, if γ < (1 − β)/2. Since, the result follows as soon as the second expression goes to zero. Now from Theorem 13, with This goes to zero provided that, 1/2βC α π p > γ.
To prove Corollaries 4 and 5, we first need another auxiliary result: Lemma 14 Let (S, d) be a separable metric space and X n : (Ω, A) → (S, B(S)), n ∈ N and X : (Ω, A) → (S, B(S)) be measurable. Then X n p → X if and only if for every subsequence n(k) there exists a further subsequence n(k(l)) such that X n(k(l)) → X a.s.
Proof If X n p → X, then so does any subsequence X n(k) . By well-known results, see e.g. van der Vaart (1998, Chapter 2), this implies that there exists a further subsequence X n(k(l)) such that a.s. convergence holds.
As is well known, there exists a metric ρ on P(S) such that ρ(X n , X) → 0 iff X n p → X, see e.g., Dudley (2002, Chapter 11). Now assume that for any subsequence we can find a further subsequence such that (51) holds, but the overall sequence does not converge in probability. Then we can build a subsequence such that for some ε > 0, ρ(X n(k) , X) ≥ ε for all elements of that subsequence. Thus any further subsequence will also not convergence in probability and consequently cannot converge a.s. This proves the claim.
We note that the set A with P (A) = 1 on which (51) holds is allowed to depend on the subsequence. Corollary 4 and 5 are finally proven jointly in the following Corollary. The proof is motivated by the tools used in Berti et al. (2006).
Corollary 15 Assume that one of the following two sets of conditions holds: (a) The kernel k is bounded, (jointly) continuous and has Moreover, y → k(y 0 , y) is vanishing at infinity, for all y 0 ∈ R d .
(b) The kernel k is bounded, shift-invariant, (jointly) continuous and ν in the Bochner representation in (11) is supported on all of R d . Moreover, Y takes its values almost surely in a closed and bounded subset of R d .
Then, under the conditions of Theorem 2, we have for any bounded and continuous function f : R d → R that DRF consistently estimates the target τ ( for all points of continuity t ∈ R d and t ∈ R of F Y | X = x (·) and F −1 Proof As shown in (Sriperumbudur, 2016, Theorem 3.2), (a) implies that k metrizes weak convergence. Similarly, from Theorem 9 in Sriperumbudur et al. (2010), it follows that k is characteristic on the compact subspace of R d in which Y takes its value almost surely. Thus, ignoring the Null set, Theorem 23 in Sriperumbudur et al. (2010) implies that k metrizes the weak convergence in this case as well. Thus in both cases μ n(k(l)) (x) − µ(x) H → 0 implies weak convergence ofμ n(k(l)) (x) to µ(x). From Theorem 14, for any subsequence, we can choose a further subsequence, such that μ n(k(l)) (x) − µ(x) H → 0 , a.s. and since it is assumed that · H metrizes weak convergence,μ n(k(l)) (x) converges weakly to µ(x) on a set A, depending on the subsequence, with P(A) = 1. Let C b (R d ) denote the space of all bounded continuous functions on R d . By the Portmanteau theorem (see e.g. Dudley (2002)), this implies that on A where we omitted the dependence on the subsequence. But, since the subsequence n(k) was arbitrary, this immediately implies for the overall sequence. On the other hand, (II) implies that on A, for the given subsequence n(k(l)), is continuous. Using for each ω ∈ A the arguments in van der Vaart (1998, Chapter 21), this implies that on A for the given subsequence. Again, as the subsequence n(k) was arbitrary, this implies the result.

C. Simulation Details
In this section we describe in detail all our simulations shown in the main paper, together with the data used in the analysis. The data sets are available in the R-package drf as well.

C.1 Air quality data
Data. This data is obtained from the website of the Environmental Protection Agency website (https://aqs.epa.gov/aqsweb/airdata/download_files.html). We have daily measurements for 5 years of data (2015-2019) for 6 'criteria' pollutants that form the Air Quality Index (AQI): • O 3 -ground ozone (8 hours' average, expressed in pieces per million (ppm)) • SO 2 -sulfur dioxide (1 hour average, expressed in pieces per billion (ppb)) • CO -carbon monoxide (8 hours' average, expressed in pieces per million (ppm)) • NO 2 -nitrogen dioxide (1 hour average, expressed in pieces per billion (ppb)) • PM2.5 -fine particulate matter smaller than 2.5 micrometers (24 hours' average, expressed in νg/m 3 ) • PM10 -large particulate matter, smaller than 10 micrometers (24 hours' average, expressed in µg/m 3 ) For the above quantities, we have the maximal and mean value within the same day. In our analysis we have used only the maximal intraday values. The pollutants are measured at different measurement sites. For each site we have information about • site address (street, city, county, state, zip code) • site coordinates (longitude and latitude) • site elevation • location setting (rural, urban, suburban) • how the land is used within a 1/4 mile radius (agricultural, forest, desert, industrial, commercial, residential, blighted area, military reservation, mobile) • date when the measurement site was put in operation • date when the measurement site was decommissioned (NA if the site is still operational) We have information about 19 739 sites, much more than the number of 2 419 sites from which we have measurements in years 2015-2019, since many sites were only operating in the past and are decommissioned. In total there is 5 305 859 pollutant measurements. Many pollutants are measured at the same site, but it is important to note that not every site measures every pollutant, so there is a lot of 'missing' measurements. It can also occur that there are several measuring devices for the same pollutant at the same site, in which case we just average the measurements across the devices and do not report those measurements separately.
Analysis. Since we have a lot of missing data, we use only the data points (identified by the measurement date and the measurement site) for which we have measurements of all the pollutants chosen as the responses. For that reason we also do not train DRF with all 6 pollutants as the responses, but only those that we are interested in, since only 64 sites measure all pollutants. For computational feasibility, we only use 50 000 of the available measurements for the training step. We also omit the states Alaska and Hawaii and the US territories for plotting purposes.
To obtain the results displayed in Figure 2, we train the DRF with the measurements (intraday maximum) of the two pollutants PM2.5 and NO 2 as the responses, and the site longitude, latitude, elevation, land use and location settings as the predictors. We manually choose two decommissioned measurement sites (for which we have no measurements in years 2015-2019) as the test points. For each test point we obtain the weights to all training measurements. We further combine the weights for all measurements corresponding to the same site, which is represented by the symbol size in the top row. The bottom row shows the estimated distribution of the response, where the transparency (alpha) each training point corresponds to the assigned weight. We also add some estimated contours.
For all plots in Figures 2 and 5, we train the single DRF with the same set of predictor variables and take the three pollutants O 3 , SO 2 and PM2.5 as the responses. In this way we still have training data from many different sites (see the above discussion on missing data) and moreover, those are the 3 pollutants that most likely cross the threshold for the "Good" AQI category set by the EPA. Carbon monoxide (CO), for example, almost never crosses this threshold.
In left plot of Figure 6, we compare the estimated CDF value with the standard classification forest which has the indicator 1(O 3 < 0.055ppm, SO 2 < 36ppb, PM2.5 < 12.1µg/m 3 ) as the univariate response. In the right plot, we obtain the estimated CDF by fitting for each threshold a separate classification forest with an indicator 1(O 3 ≤ threshold).
We pick a test point such that the classification performs bad, just to illustrate that its estimated CDF need not be monotone, which cannot happen with DRF. In most of the cases, the estimated CDFs are very similar, as can also be seen from the left plot in Figure 6.

C.2 Benchmark Analysis
In this part we compare the resulting distributional estimates of DRF with several benchmark methods on a number of data sets. Because our target of estimation is now the whole conditional distribution one needs to use a distributional loss, and there appears to be no well-established choice in the literature. Furthermore, for any test point x i we only have one observation y i from P(Y | X = x), which makes performance evaluation of our estimatorP(Y | X = x) hard. We thus use the following performance measure: • (NLPD loss) For a fixed conditional distribution estimatorP(Y | X = x i ), we sample a set of m = 500 observations from which we estimate the conditional density via a Gaussian kernel estimator, using the L 2 loss with scale components and the median heuristic for the choice of the bandwidth parameter. We then evaluate the negative loglikelihood of the test observation y i implied by the kernel estimate of the distribution and average over the test set (consisting of multiple pairs (x i , y i )). To reduce the dependence of these results on single large values of the log-likelihood, we use an 0.05−trimmed mean to average the losses over the training set.
This loss definition provides a fair way to compare the ability to estimate the conditional distribution since most of the candidate methods only allow for sampling from the estimated conditional distributionP(Y | X = x i ).

C.2.1 Competing methods
We compare DRF that uses the MMD splitting criterion with many existing methods that can be used for estimation of the conditional distribution.
• Nearest Neighbor (k-NN): The standard k-nearest neighbors algorithm with the Euclidean metric. An estimated conditional distributionP(Y|X = x) at a test point x is defined by a uniform distribution over the k nearest observations in the training set. k is chosen to be the square root of the training set size.
• Gaussian kernel (kernel): The estimate of the conditional distributionP(Y|X = x) at a test point x is obtained by assigning to each training observation (x i , y i ) the weight proportional to the Gaussian kernel k(x, x i ), analogously to usual kernel estimation methods. Median heuristic is used for bandwidth selection.
• Homogeneous distribution models: This method makes the homogeneity assumption that the residuals have constant distribution and only the conditional mean changes. The estimate of the conditional distributionP(Y|X = x) is obtained by first fitting a regression method of choice, computing the residuals, assigning the same weight to every residual and then adding those residuals to the predicted mean. We chose three different methods for the mean estimation: 1. Random Forests (RF), a classical univariate regression forest is fitted independently for each response component; 2. Extreme Gradient Boosting (XGBoost), a tree gradient boosting model (as described in Chen and Guestrin (2016)) is fitted independently for each response; 3. Deep Neural Network (DNN), a single deep neural network is fitted to predict the conditional mean of each response.
• Conditional Generative Adversarial Neural Network (CGAN): The estimated conditional distributionP(Y|X = x) is obtained through sampling from the discriminator with conditional feature x. The implementation of the CGAN is taken from Aggarwal et al. (2019). The architecture of the neural networks was taken to be the best one for the considered data sets among a set of candidates.
• Conditional Variational Auto-Encoder (CVAE): The estimated conditional dis-tributionP(Y|X = x) is obtained through sampling from the decoder of the CVAE with conditional feature x. The implementation of the CVAE follows the one in Sohn et al. (2015). The architecture of the neural networks was taken to be the best one for the considered data sets among a set of candidates.
• Masked Autoregressive Flow (MAF): The estimated conditional distribution P(Y|X = x) is obtained through sampling from the normalizing flow model with conditional feature X. The implementation of the model follows the one presented in Papamakarios et al. (2017). The number of layers is chosen to be the best value from a set {5, 10} for the considered data set.
• Conditional Mean Embedding (CME): The CME is calculated as in (16) with the weights given as in e.g., Song et al. (2009);Muandet et al. (2017); Park and Muandet (2020). We choose both kernels to be Gaussian kernels with σ = 0.01, as in Park and Muandet (2020) and set λ = 0.01. The estimated conditional distribution P(Y|X = x) is obtained through sampling from the obtained weights, renormalized such that they lie in [0, 1] and sum up to one.

C.2.2 Benchmark data sets
Many benchmark data sets used come from the multiple target regression literature, where only the conditional means of the multivariate response is considered. We have used the data sets: jura, slump, wq, enb, atp1d, atp7d, scpf, sf1 and sf2 collected in the Mulan (Tsoumakas et al., 2011) library. Description about the dimensionality of the data sets, together with the descriptions of the outcomes and the regressors can be found in Tsoumakas et al. (2011) with links to the relevant papers introducing these data sets. In each data set categorical variables have been represented by the one-hot dummy encoding, the observations with missing data were removed together with constant regressors.
We additionally added 5 data sets obtained from the data sets used in the main paper: • copula: Simulated Gaussian copula example where the response Y is bivariate and whose marginal distribution is N (0, 1), but the correlation between Y 1 and Y 2 depends on X 1 .
• birth1: This data set is created from the CDC natality data and contains many covariates as predictors and the pregnancy length and birthweights as the responses.
• birth2: This data set is similar as the above one, but we take pregnancy length as the predictor and add 3 more measures of baby's health as the response: APGAR score measured 5 minutes after birth and indicators whether there were any abnormal conditions and congenital anomalies.
• wage: This data set is created from the 2018 American Community Survey. We take the logarithmic hourly wage and gender as the response, as it was done in the fairness example in the main paper.
• air: This data set is obtained from the EPA air quality data. All six pollutants were taken as the response and we add both the information about the measuring site (location, which setting it is in, etc.), as well as the temporal information when the measurement has taken place (month, day of the week).

C.3 Births data
Data. This data set is obtained from the CDC Vital Statistics Data Online Portal (https: //www.cdc.gov/nchs/data_access/vitalstatsonline.htm) and contains the information about the ≈ 3.8 million births in 2018. However, as we do not need this many data points, we subsample 300 000 of them. Even though the original data contains a lot of variables, we have taken only the following variables from the source data: • mother's age, height, weight before the pregnancy and BMI before pregnancy • mother's race (black, white, asian, NHOPI, AIAN or mixed), marital status (married or unmarried) and the level of education (in total 8 levels) • father's age, race and education level • month and year of birth • plurality of the birth (how many babies were born at once) • whether and when the prenatal care started • length of the pregnancy • delivery method (vaginal or C-section) • birth order -the total number of babies born by the same mother (including the current one) • birth interval -number of months passed since last birth (NA if this is the first child) • number of cigarettes smoked per day on average during the pregnancy • birthweight (in grams) and gender of the baby • APGAR score (taken after 5min and 10min) • indicators whether baby had any abnormal condition or some congenital anomalies Analysis. After removing the data points with any missing entries and taking only the data points where the race of both parents is either black, white or Asian (for nicer plotting), we are left with 183 881 data points. We use randomly chosen 100 000 data points for training the DRF. We take the birthweight and the pregnancy length as the bivariate response and for the predictors we take: mother's age, race, education, marital status, height, BMI; father's age, race and education level; birth plurality, birth order, delivery method, baby's gender, number of cigarettes and indicator whether prenatal care took place. For arbitrary test points from the data we can get the estimated weights by the fitted DRF, thus estimating the joint distribution of birthweight and pregnancy length conditional on all other variables mentioned above. Two such distributions are shown in Figure 10. In addition we use the weights to fit a parametric model for the mean and 0.1 and 0.9 quantiles. This is done as follows: • We slightly upweight the data points where the pregnancy length is significantly above or below the usual range. This is to avoid the bulk of the data points to dominate the fit obtained for very long or short pregnancies.
• We apply the transformation f (·) = log(log(·)) on both the pregnancy length and the birthweight since then the scatterplots look much nicer.
• We estimate the mean with smoothing splines with a small manually chosen number of degrees of freedom.
• The fitted mean is subtracted from the response (birthweight). The residuals seem well behaved with maybe slight, seemingly linear trend in standard deviation.
• We fit the 0.1 and 0.9 quantiles as the best linear functions that minimize the sum of quantile losses, by using the quantreg package (Koenker et al., 2012).
• The data is transformed back on the original scale by using the function f −1 (·) = exp(exp(·)).
For the right plot in Figure 10, we have the following causal graph, as mentioned in the main paper: We want to determine the direct effect (indicated in bold) of the twin pregnancy T on the birthweight B that is due to sharing of resources by the babies (space, food etc.) and is not due to the fact that twin pregnancy causes shorter pregnancy length L, which in turn causes the smaller birthweight. Another big issue is that we have confounding factors Z which can directly affect B, L and T . For example, the number of twin pregnancies significantly depends on the parents' race, but so do the pregnancy length and the birthweight, e.g. black people have more twins, shorter pregnancies and smaller babies. We take all other variables as the potential confounders Z and adjust for all of them (mother's age, race, education, marital status, height, BMI; father's age, race and education level; birth plurality, birth order, baby's gender, number of cigarettes and indicator whether prenatal care took place). In order to do it, we fit the same DRF as before, where Z and T are the predictors and (B, L) is the bivariate response for which we can fit the parametric model described above. We compute then the interventional distribution P(B | do(T = t, L = l)) for all values of t and l, by using the do-calculus to adjust the confounding Z via the backdoor criterion (Pearl, 2009), where we also use the obtained parametric regression fit. In this way we can generalize the fit well, which is important when doing the do-calculus, since we are interested in some hypothetical combinations of covariates which might not occur frequently in the observed data, such as very long twin pregnancies.

C.4 Wage data
Data. The PUMS (Public Use Microdata Area) data from the 2018 1-Year American Community Survey is obtained from the US Census Bureau API (https://www.census.gov/ content/dam/Census/data/developers/api-user-guide/api-guide.pdf). The survey is sent to ≈ 3.5 million people annually and aims to give more up to date data than the official census that is carried out every decade. The 2018 data set has 3 214 539 anonymized data points for the 51 states and District of Columbia. Even though the original survey contains many questions, we have retrieved only the subset of variables that might be relevant for the salaries: • person's gender, age, race (AIAN, black, white, asian, mix, NHOPI, other), indicator of hispanic origin, state of residence, US citizenship indicator (5 ordered levels), indicator whether the person is foreign-born • person's marital status, number of own children in the same household and the number of family members in the same household • person's education level (24 ordered levels) and level of English knowledge (5 ordered levels) • person's employment status (employed, not at work, not in workforce, unemployed) • for employed people we have annual salary earnings, number of weeks worked in a year and average number of hours worked per week • for employed people we have employer type (government, non-profit company, for-profit company, self-employed), occupation (530 levels), industry where the person works (271 levels) and the geographical unit where the person works (59 levels) • statistical weight determined by the US Census Bureau which aims to correct sampling bias For our purposes, since we want to analyze the unfairness of the gender pay gap, we consider only employed people that are at least 17 years of age, have worked full-time (at least 48 weeks in a year) and have worked at least 16 hours a week on average. We also omit the self-employed persons, since they often report zero annual salary and the pay gap there, if exists, cannot be called unfair as the salary is not determined by any employer. Since there are no missing data which would need to be omitted, we finally end up with 1 071 866 data points.
Analysis. We scale the salary with the amount of time spent working (determined from the number of weeks worked and average hours worked per week) to compute the logarithm of the hourly wages. The scaling with the time spent working is necessary, since full-time employed men spend on average 11% more time working than women. The logarithmic transformation is used since the salaries are very skewed (positively) and logarithmic wages show nice behavior.
We also reduce the large number of levels of some of the categorical variables: for the occupation we use the group of 530 jobs into 20 categories provided in the SOC system (https://www.bls.gov/soc/); for the industry information we group the 271 possibilities in 23 categories as is done in the NIACS classification (https://www.bls.gov/bls/naics.htm); for the work place we group the 59 US states and foreign territories into 9 economic regions (including the "abroad" category), as determined by the Bureau of the Economic Analysis (https://apps.bea.gov/regional/docs/regions.cfm).
We want to investigate how the logarithmic hourly wage W is affected by the gender G, depending on the other factors Z: age, race, hispanic origin, citizenship, being foreignborn, marital status, family size, number of children, education level, knowledge of English, occupation, industry type and place of work. To do this, we train DRF with bivariate response (W, G) and predictors Z on a subsample of 300 000 data points. With it we can answer the following: For fixed values of the covariates Z = z, what are the distributions of salaries of men and women. In addition, we can determine the "propensities", i.e. the proportion of men and women corresponding to Z = z. This information is displayed in the top row of Figure 11 for a combination of covariates corresponding to some person in the left-out data. It illustrates how the distribution of salaries and their relationship can vary with different covariates Z.
We do not only want to determine how different covariates Z affect the salary distribution, but we want to quantify the overall fairness of the pay, after appropriate adjustments. In Figure 11, we can see that the observed salaries of men and women differ noticeably, and this difference in the logarithmic wages means that an average woman has 17% smaller salary than an average men. However, the question is how much of this difference is "fair". For example, the effect of the gender on the salary can be mediated through some variables such as, for example, the occupation, workplace or the level of education and we are only interested in the direct effect. This is illustrated in the following causal graph: G Z W If we assume that people have the freedom to choose such variables themselves, the pay gap which arises from such different choices for men and women is fair and those variables are resolving variables (Kilbertus et al., 2017). Another way that the pay gap can be explained is that some of the variables are not statistically independent of the gender in the population of full-time employed people (e.g. the race or the age), but they themselves have an effect on the salary.
In order to address those issues, we compute the distribution of the nested counterfactual W (male, Z(female)), corresponding to the wages of a person that has characteristics Z as a woman, but which was treated as a man for obtaining the salary. Such distribution can be computed from the DRF, as described in the main paper: we randomly draw a female person and for its characteristics z we obtain the conditional distribution of wages of men with those characteristics P(W | G = male, Z = z) via the weights. Those distributions are averaged over random draw of 1 000 women (that were not used in the training step of the DRF). In case that the difference in salary is fair, the distribution of the counterfactual salary W (male, Z(female)) should be exactly the same as the observed distribution of women's wages. However, we can see that this is not the case and that the median salaries of the two distributions differ by 11%. Even though this is smaller than the 17% we obtain by comparing only the observational distributions, it still shows that women are paid less compared to men.
The univariate response (case d = 1) is by far the most studied case in the regression literature. However, at the level of the whole conditional distribution and compared to the multivariate case, the range of practically interesting targets τ (x) is quite reduced, e.g. conditional mean of some functional E(f (Y ) | X) or conditional quantiles Q α (Y | X).
In Figure 4, we have compared the performance of DRF (which uses the MMD splitting criterion) with 3 different tree-based univariate methods that can estimate the conditional quantiles in the univariate case: • QRF: the quantile regression forest introduced in Meinshausen (2006), which is equivalent to DRF CART in the univariate case and uses the standard forest construction (Breiman, 2001) to get the weights.
• GRF: the quantile forest proposed in  based on the generalized random forest algorithm.
Additionally to the visual inspection of the performance given in Figure 4, we present here a formal performance comparison for the three simulation scenarios also described in the main paper. The first two scenarios correspond exactly to the examples given in  for the quantile version of the GRF, which serve to illustrate its advantage compared to the conventional quantile regression forest (QRF) (Meinshausen, 2006). Scenario 3, in addition, aims at assessing the ability to detect a change of distribution that does not relate to a change in the first two moments.  Table 2: Average quantile losses for scenarios 1 (left), 2 (middle), 3 (right) over the repeated out-of-sample validations.
The performance of each method is evaluated as follows: We consider the quantile (pinball) loss for the resulting quantile estimates provided by each candidate method for the different percentiles α ∈ {0.1, 0.3, 0.5, 0.7, 0.9}. The losses are presented and computed based on repeated (10 times) out-of-sample validation (with a 70 − 30% ratio between the training and testing sets sizes). The results are presented respectively for each scenario in Table 2. We additionally include the estimates obtained by k-nearest neighbor algorithm for several different values of k.
Furthermore, Table 3 shows non-inferiority of DRF compared to the standard Random Forest for the classical task of estimating the conditional mean. We observe that DRF has a good relative performance that makes it on par with existing algorithms, some of which specially designed for the problem of estimating conditional quantiles. Furthermore, it seems that the MMD splitting criterion improves the CART criterion for distributional regression in a general heterogeneous case (see e.g. scenarios 2 and 3), since the CART criterion is suitable only for detecting the change in the conditional mean, unlike MMD. Dependence of the estimated quantiles on X 1 for each method (except the k-nearest neighbors) is displayed in the main paper in Figure 4. In addition, the estimates of 2-Wasserstein distance to the true conditional distribution, quantifying the difference in the estimated CDFs, are shown in Figure 12. Figure 12: Scatter plot of discrete estimates of the 2-Wasserstein distance between the estimated and true conditional distribution against X 1 for a grid of test points of the form (x 1 , 0, . . . , 0). The 2-Wasserstein distance is estimated over a grid of 100 quantiles with levels equally spaced on [0, 1]. Different colors corresponds to different methods: DRF (red), GRF (green), QRF (blue), TRF (purple).

D.2 Heterogeneous regression and causal effects
We explore here the performance of DRF on the synthetic data for the setup of heterogeneous regression, where we want to obtain the regression fit of Y on the explanatory (or treatment) variables W, but where this fit might change depending on values of X. This can be done by DRF by using X as predictors and (W, Y ) as the response and then using some standard regression method for regressing Y on W in the second step, having already obtained the weights that describe the conditional distribution P((W, Y ) | X = x).
The most important such setup is when the data come from the following causal graph: In the case of such a causal graph, X are confounding variables, which we need to adjust for to understand the causal effect of W on Y . Not only can the marginal distributions of Y and W be affected by X, but also the regression fit (e.g. the regression coefficients).

D.2.1 CATE and ATE
One special case of this setup that is intensively studied in the causal literature is when W is a (univariate) binary treatment variable. In this case we are interested in the distribution of the potential outcomes Y (W = 0) and Y (W = 1) and especially in their difference. It is commonly measured by using the Conditional Average Treatment Effect Competing methods. We will compare the performance of the DRF with the following methods, specially designed for estimation of the CATE (or ATE) • Double Machine Learning (DML) of Chernozhukov et al. (2018), which assumes the model Y = m(X) + W θ + with constant treatment effect and can thus only be used for estimating ATE and not CATE.  (2018)) In order to make the comparison fair, we use the local centering approach for DRF as well.
Data. We will use the following data models for our simulations, where the first three are taken directly from Athey et al. (2019): 1. In this model X 3 is a confounder affecting both W and Y : X ∼ U (0, 1) p , W | X ∼ Bernoulli 1 4 (1 + β 2,4 (X 3 )) , Y | X, W ∼ 2 X 3 − 1 2 + N (0, 1), where β a , b(x) is the density of the beta random variable with parameters a and b.
Results. For every model we generate n data points (X 1 , . . . , X p , W, Y ) i=1,...,n . We run all methods and compute the root mean squared error of the obtained CATE estimate on a randomly generated test set X test containing 1000 data points. CATE corresponds to the coefficient of W in the data generating mechanism of Y . We repeat the same procedure 100 times and report the average result. For methods other than the DML, we estimate ATE by averaging the CATE estimates over the randomly generated test set. The results can be seen in Table 4 and Figure 13. Even though DRF is performing less well in general compared to the methods that are specially designed for the task of estimating CATE, we can still see that its estimates are fairly good.

D.2.2 Continuous W , linear treatment effect
When the treatment variable W is continuous, many methods designed for binary treatments, such as the X-learner (Künzel et al., 2019) cannot be used. However, many important real-word examples fall within this framework. As an example, we might be interested in how the amount of medicine W affects some biological parameter of interest Y (conditionally on X). When W affects Y linearly conditionally on X, one can still use the Causal Forest (CF) Wager and Athey, 2018) method, which makes the splits based on the slope of the conditional linear fit Y ∼ W . Due to its generality and versatility, DRF can trivially be used in such setting as well.
To illustrate this, we consider the Model 3, as described in the previous section, which is also taken from , but where we change the distribution of the binary    treatment variable so that it is continuous and it has a normal distribution with the same mean and variance, which depend on X. In this model W affects Y linearly, which is a  (right), which uses the local centering approach. The data is generated from Model 3 described in the previous section, with continuous treatment variable W with the same mean and variance conditionally on X.
crucial assumption for the CF approach to work. We take n = 10000 and p = 10. The concept of CATE does not exist in this form in such setup and therefore we consider how the forest obtained by each method estimates both the intercept and the slope of the fit Y ∼ W , conditionally on X. The results can be seen in Figure 14. We see that the estimate of the slope for DRF is slightly worse than for the CF, whose forest construction is specially designed for estimating the conditional slope. However, DRF estimates the intercept significantly better than the CF, especially in combination the local centering approach, which uses the centered data Y − Y (X) and W − W (X) instead. In this example the slope depends only on X 2 , whereas the intercept depends on both X 1 and X 2 . Since CF targets only the slope for forest construction, it will split mostly on X 2 and not on X 1 , which leads to poor estimate of the intercept term. On the other hand, DRF splits both on X 1 and X 2 , depending on the size of their effect on the joint distribution of (W, Y ). For many applications, especially in causality (see also the example in next section), it is essential to know the whole conditional distribution P(Y | W, X) so the DRF approach might be more beneficial than the CF.

D.2.3 Nonlinear treatment effect
There are very few methods that can estimate the treatment effect when the treatment variable W is continuous and affects Y nonlinearly, as it is commonly the case in real world settings. Here we demonstrate that DRF can be easily used in such setup as well, as opposed to the CF, which assumes a linear, though heterogeneous, treatment effect of W on Y . We further show how the obtained regression fits Y ∼ W conditional on X can be used to estimate the causal effect E[Y | do(W = w)], as it is done in the main paper in the birth data example: We compare the performance of DRF with the straightforward and commonly used approach, where we first regress Y on (W, X) and use this regression fit which estimates E[Y | W = w, X = x] together with the above formula to estimate the causal effect.
Therefore, W affects Y highly nonlinearly through a sine function. X 3 is a confounding variable that affects the marginal distributions of Y and W . X 2 regulates the error level for Y , whereas X 1 makes the treatment effect heterogeneous. This is illustrated in the following plot: Figure 15: Visualization how X 1 , X 2 , X 3 affect the conditional nonlinear regression fit Y ∼ W . X 1 changes the effect size, X 2 changes the noise level, whereas X 3 is a confounding variable which affects the means of W and Y .

Results
In Figure 16 we can see the estimated joint distribution of (Y, W ) conditionally on X = x, where the values of X 1 and X 3 vary, while the rest are fixed (even though X 2 also affects the conditional distribution, the effect is much weaker than for X 1 of X 3 , see Figure  16). We see that the estimated distribution matches the true regression line, denoted in red, very well. The estimated distribution induced by the DRF weights enables us to fit some specialised regression method for regressing Y on W for every fixed value of X. The blue line indicates the fit obtained by using smoothing splines. Compared to the green line, which shows the predicted values for regression Y ∼ (W, X), it is nicer looking and also is able to extrapolate much better to values of W which have low probability conditionally on X. This extrapolation is crucial for causal applications, since for computing E[Y | do(W = w)] we are interested in what would happen with Y when our treatment variable W is fixed to be w, regardless of which values are achieved by X. However, it can easily happen that for this specific combination of X and W there are very few observed data points, which makes the estimation hard (Pearl, 2009). In this example, W tends to be small for small values of X 3 and vice-versa and thus is hard to say what would happen with Y when X 3 is large and W is set to a small value by an outside intervention.
In Figure 17, we indeed see that the estimates of the causal effect E[Y | do(W = w)] by DRF are much better. One can still see that the error increases for the border values of W , which have small probability for some values of X, since the estimation there is much harder, but this error is much less pronounced for DRF than for the standard regression approach.