# Density Regression with
Bayesian Additive Regression Trees

Vittorio Orlandi
Dept. of Statistical Science
Duke University
Durham, NC 27708
<EMAIL_ADDRESS>
& Jared Murray
Dept. of Information,
Risk, and Operations Management
University of Texas, Austin
Austin, TX 78712
<EMAIL_ADDRESS>Antonio Linero
Dept. of Statistics and Data Sciences
University of Texas, Austin
Austin, TX 78712
<EMAIL_ADDRESS>&Alexander Volfovsky
Dept. of Statistical Science
Duke University
Durham, NC 27708
<EMAIL_ADDRESS>

###### Abstract

Flexibly modeling how an entire density changes with covariates is an
important but challenging generalization of mean and quantile regression.
While existing methods for density regression primarily consist of covariate-
dependent discrete mixture models, we consider a continuous latent variable
model in general covariate spaces, which we call DR-BART. The prior mapping
the latent variable to the observed data is constructed via a novel
application of Bayesian Additive Regression Trees (BART). We prove that the
posterior induced by our model concentrates quickly around true generative
functions that are sufficiently smooth. We also analyze the performance of DR-
BART on a set of challenging simulated examples, where it outperforms various
other methods for Bayesian density regression. Lastly, we apply DR-BART to two
real datasets from educational testing and economics, to study student growth
and predict returns to education. Our proposed sampler is efficient and allows
one to take advantage of BART’s flexibility in many applied settings where the
entire distribution of the response is of primary interest. Furthermore, our
scheme for splitting on latent variables within BART facilitates its future
application to other classes of models that can be described via latent
variables, such as those involving hierarchical or time series data.

_Keywords_ Bayesian Nonparametrics $\cdot$ Conditional Density Estimation
$\cdot$ Posterior Concentration $\cdot$ Latent Variables $\cdot$
Heteroscedasticity

##  1 Introduction

Data analysis frequently concerns itself with associating the change in a
function of some response variable $y$ with a set of covariates $\mathbf{x}$.
Arguably the most common tool for this is mean regression, which focuses on
the expectation $\mathrm{E}[y\mid\mathbf{x}]$ and foregoes inference about
other parts of the conditional density $p(y\mid\mathbf{x})$, a much more
general quantity. This inflexibility has been recognized as problematic in
many modern applications (see, e.g. Wittman, 2009; Karabatsos and Walker,
2011; Gneiting and Katzfuss, 2014). It is then natural to ask if other
functionals of $p(y\mid\mathbf{x})$ are more appropriate, with two immediate
candidates: quantile regression and density regression.

Quantile regression models specific quantiles of the conditional response
distribution. This can help address settings where some quantiles carry more
probative value. For example, Black et al. (2007) note that the impact of a
job training program on the upper quantiles of the distribution is considered
to be more important by policy makers than that on the lower quantiles. While
this approach is more general than mean regression, one problem with its
application is that computing estimates of functionals of the quantiles is not
always straightforward. Another problem is that estimated quantiles oftentimes
do not obey the monotonicity constraint inherently satisfied by the true
distribution. While there are approaches for joint modeling of quantiles (e.g.
Kadane and Tokdar (2012); Sangnier et al. (2016)) or post-hoc reordering of
estimates (e.g. Chernozhukov et al. (2010)), the fundamental limitation with
the approach is that individual quantiles are being targeted as proxies for
features of the distribution as a whole.

By modeling the entire probability distribution of the response, density
regression methods perform a substantially harder task than mean or even
quantile regression. In doing so, however, they are able to compute coherent
point estimates and perform uncertainty quantification for arbitrary
functionals of the distribution that may be of interest. Studies of income
inequality, for example, typically take into account the impact of variables
on the entire income distribution (e.g. Gerfin, 1994; Daly and Valletta,
2006). To date, there is a large literature on density regression. We
contribute to this body of work by providing a general, yet reasonably
structured, formulation of the problem that enjoys important theoretical
guarantees and yields an efficient sampler with strong empirical performance.

In this work, we consider the following model for the conditional density:

$p(y\mid\mathbf{x},f,\sigma)=\int_{0}^{1}\frac{1}{\sigma(\mathbf{x},u)}\phi\left({\frac{y-f(\mathbf{x},u)}{\sigma(\mathbf{x},u)}}\right)du,$
(1)

which is parameterized by the latent variable $u$. This model generalizes the
covariate dependent mixture model discussed in detail below. At the same time,
it posits a fairly specific form by which $f$ and $\sigma$ interact to
generate a conditional density. Various forms of this model have been
considered in Pati et al. (2011), Kundu and Dunson (2014), and Zhou et al.
(2017), where Gaussian process priors were placed on $f(x,u)$ and inverse
gamma priors on either $\sigma(x,u)=\sigma$ or $\sigma^{2}(x,u)=\sigma^{2}$.
Instead, we propose placing BART-type priors on both $f$ and $\sigma$. While
BART was designed for mean regression, we also consider a modification
allowing it to model variance function to grant the model added flexibility in
finite sample settings, and generalize both to accommodate latent variables.
Asymptotically, we show that the posterior specified by our model concentrates
around a true underlying density, provided its log is $\alpha$-Hölder for
$0<\alpha\leq 2$. The rate of concentration is removed from the minimax rate
by a factor of $\alpha/(\alpha+1)$ because BART generates piecewise constant
functions; however, a near–minimax rate can be attained by using the SBART
model of Linero and Yang (2018) and restricting to a slightly smaller function
class. In addition to yielding a more efficient sampler, we also show that
this model outperforms its counterpart employing Gaussian process priors, as
well as a variety of other methods for density regression in several empirical
evaluations.

We also apply our method to two real world datasets. In the first, we compute
quantile growth targets for students in elementary and middle school
mathematics classes using data provided in Betebenner et al. (2011). Mean
regression approaches fail to identify interesting aspects of the data — such
as that conditional distributions of test scores become more skewed over time
— and quantile regression approaches can suffer from quantile crossing and a
limited description of the uncertainty in their estimates. Our model addresses
both these problems simultaneously. In the second, we study returns to
education from US census microdata originally compiled by Angrist et al.
(2006). The returns are a nonlinear functional of quantiles of the wage
distribution and therefore well suited for analysis by density regression
methods, which fully capture uncertainty about their estimates.

The paper proceeds as follows. The remainder of the introduction discusses
past work on density regression models. Section 2 gives a brief overview of
BART and its application to modeling mean and variance functions. In Section
3, we motivate our use of BART for modeling components of the conditional
density and state our full model for density regression, which we refer to as
DR-BART. Section 4 outlines theoretical results concerning our model, namely
the rate at which the posterior contracts about a true density. Section 5
compares DR-BART to other models for density regression on a variety of
simulated datasets and Section 6 applies DR-BART to two real world datasets
from education and economics. Section 7 concludes.

###  1.1 Related Work

Our proposal generalizes a common approach to density regression that uses
covariate dependent mixture models. In these models, the conditional density
is given by:

$p(y\mid\mathbf{x},\theta,\pi)=\sum_{h=1}^{k}\pi_{h}(\mathbf{x})\mathcal{K}(y;\theta_{h}(\mathbf{x})),$
(2)

allowing for $k=\infty$, where $\mathcal{K}(\cdot,\cdot)$ is a positive
definite kernel function and $\pi_{h}(\mathbf{x})$ are covariate-dependent
mixture weights. We restrict attention to normal kernels:

$p(y\mid\mathbf{x},\theta\equiv(\mu,\sigma),\pi)=\sum_{h=1}^{k}\pi_{h}(\mathbf{x})\phi_{\sigma_{h}(\mathbf{x})}(y-\mu_{h}(\mathbf{x})),$
(3)

where $\phi_{\sigma}(z)=(1/\sigma)\phi(z/\sigma)$ and $\phi$ denotes the
standard normal pdf – but the extension to other kernels is straightforward.
In practice $\pi$ or $\theta$ may not vary with $\mathbf{x}$, or may vary in a
limited way. (e.g., by taking
$\mu_{h}(\mathbf{x})=\mathbf{x}^{\prime}\beta_{h}$
$\sigma_{h}(\mathbf{x})=\sigma_{h}$). Our proposed model in (1) recovers that
in (3) when $f$ and $\sigma$ are step functions with the same points of
discontinuity (which can depend on the covariates).

Models of the form in (3) appear in the machine learning literature as
“mixtures of experts” (Jacobs et al., 1991; Jordan and Jacobs, 1994) where the
initial focus was on using these models for flexible mean regression or
classification. Geweke and Keane (2007) and Villani et al. (2009) study models
of this form for semiparametric density regression, using finite $k$ and
multinomial probit and logit regression models for $\pi(\mathbf{x})$. While it
is possible to get consistency properties for large classes of conditional
densities (see e.g. Norets et al. (2010); Pati et al. (2013); Norets and
Pelenis (2014) and the monograph Norets and Pati (2014)), practical experience
in finite samples suggests that there can be value in allowing the kernel
variance to depend on $\mathbf{x}$, as this can reduce the number of clusters
required for an accurate approximation. Villani et al. (2009) provide
simulated examples and discussion in the case of fixed $k$, and we will
revisit this point in the context of the models introduced here (Section 5).

Another approach proposes to leverage the joint model for $(y,\mathbf{x})$ as
a convenient device for inducing a particular conditional model as
$p(y\mid\mathbf{x},\theta)=p(y,\mathbf{x}\mid\theta)/p(\mathbf{x}\mid\theta)$
(as in e.g. Muller et al. (1996); Park and Dunson (2010), among others). There
is some cost to the joint modeling approach – which itself has commanded a
large literature (see West and Escobar (1993); Muller et al. (1996) and
variations in Shahbaba and Neal (2009); Taddy and Kottas (2010); Molitor et
al. (2010); Wade et al. (2011); Dunson and Bhattacharya (2011); Hannah et al.
(2011); Wade et al. (2014)) – in terms of computation and accuracy as the
dimension of the covariate vector grows (see Hannah et al. (2011); Wade et al.
(2014)). The posterior can also depend on the distribution of the covariates,
even when care is taken to separate the parameter spaces in the prior, as the
auxiliary joint model assumes a common clustering for the response and the
covariates (see Griffin’s discussion of Dunson and Bhattacharya (2011); also,
Walker and Karabatsos (2013) and Wade et al. (2014)). Thus, other
nonparametric Bayesian models focus explicitly on the conditional
distributions of interest; these date to (at least) MacEachern’s seminal work
on dependent Dirichlet processes (DDPs) (MacEachern, 1999, 2000). A DDP is a
prior for a collection of distributions such that at each covariate value the
process is marginally a DP. Models in this class include De Iorio et al.
(2004); Griffin and Steel (2006); Dunson and Peddada (2008); De Iorio et al.
(2009); Wang and Dunson (2011), and numerous other specializations to
spatiotemporal or hierarchical models. Barrientos et al. (2012) characterize
the DDP in terms of copulas and provide results about its support and about
kernel mixtures using the DDP.

##  2 Heteroscedastic Regression with BART priors

###  2.1 Bayesian Additive Regression Trees (BART)

Bayesian Additive Regression Trees (BART) were introduced by Chipman et al.
(2010) (henceforth CGM) as a nonparametric prior over a regression function
$f$ designed to capture complex, nonlinear relationships and interactions.
Specifically, for observed data pairs
$\mathcal{D}=\\{(y_{i},\mathbf{x}_{i});1\leq i\leq n\\}$, CGM propose the
regression model:

$y_{i}=f(\mathbf{x}_{i})+\epsilon_{i},\quad\epsilon_{i}\overset{iid}{\sim}N(0,\sigma^{2}).$
(4)

The BART prior represents $f$ as the sum of many piecewise constant regression
trees. Each tree $T_{h},\;1\leq h\leq m$ consists of a set of interior
decision nodes (where decisions are generally of the form $x_{j}<c$ for some
value $c$) and a set of $b_{h}$ terminal nodes. The terminal nodes have
associated parameters $M_{h}=(\mu_{h1},\mu_{h2},\dots\mu_{hb_{h}})^{\prime}$.
For each tree there is a partition of the covariate space
$\\{\mathcal{A}_{h1},\dots,\mathcal{A}_{hb_{h}}\\}$ with each element of the
partition corresponding to a terminal node. A tree and its associated
parameters define step functions:

$g(\mathbf{x},T_{h},M_{h})=\mu_{hb}\text{ if
}\mathbf{x}\in\mathcal{A}_{hb}\text{ (for $1\leq b\leq b_{h}$)}.$ (5)

These functions are additively combined to obtain $f$:

$f(\mathbf{x})=\sum_{h=1}^{m}g(\mathbf{x},T_{h},M_{h}).$ (6)

This model has been shown empirically and theoretically (see, e.g. Chipman et
al. (2010); Jeong and Rockova (2020)) to be accurate, highly flexible, and
robust to the presence of irrelevant covariates. The default prior parameters
work very well in practice and the model admits an efficient sampler. More
details on BART can be found in the Supplement.

###  2.2 A BART prior for variance functions

As mentioned in the introduction, to allow for flexible density regression in
finite samples, it may be useful to allow the variance of the process to
depend on $\mathbf{x}$ as well. To this end, we adopt the log-linear BART
prior of Murray (2021) for the log-variance. Specifically, consider the
heteroskedastic regression problem
$y_{i}=f(\mathbf{x})+\sigma^{2}(\mathbf{x})\epsilon_{i}$ where
$\epsilon_{i}\overset{iid}{\sim}N(0,\sigma_{0}^{2})$. Murray places a log-
linear BART prior on $\sigma^{2}(\cdot)$:

$\log[\sigma^{2}(\mathbf{x})]=\sum_{h=1}^{m_{v}}g(\mathbf{x},T^{(v)}_{h},M^{(v)}_{h}),$

where $\\{(T^{(v)}_{h}$, $M^{(v)}_{h})\\}$ are trees and parameters for the
variance function. The prior for the exponentiated leaf parameters
$\exp(\mu_{hb}^{(v)})$ is conjugate, symmetric on the log scale, and can be
calibrated to match the expected prior range of the log-variance process. More
details on BART’s extension to log-linear models can be found in the
Supplement.

##  3 Density Regression with BART priors

Even with a flexible variance function, the normality assumption of
heteroscedastic BART may be too restrictive; for example, at any covariate
value it yields symmetric predictive distributions. In this section, we extend
the heteroscedastic BART model to general density regression problems by
introducing a continuous latent variable $U$, which is treated as an omitted
variable independent of $\mathbf{x}$. Before introducing the model in full
generality we will motivate the use of continuous latent variables in the
density estimation setting.

###  3.1 Continuous latent variable priors for a single density

Consider the following generalized location model for estimating a single
density:

$p(y)=\int_{0}^{1}\phi_{\sigma}(y-f(u))\,du.$ (7)

An equivalent representation in terms of a latent variable is

$Y=f(U)+\epsilon,\quad U\sim U(0,1),\quad\epsilon\sim N(0,\sigma^{2})$ (8)

where (7) is obtained on marginalizing over $U$. In the limit as
$\sigma\rightarrow 0$, $Y\,{\buildrel d\over{=}}\,f(U)$ where $U\sim U(0,1)$.
This class of models can be quite broad, depending on the prior for $f$; if
$f$ is the quantile function of a distribution $P$, then $Y\sim P$. We do not
restrict $f$ to be monotone; while it would be possible to do so, this
substantially increases the computational burden and since subsequent
inference is on the induced density for $y$ (or $y\mid\mathbf{x}$ below) or
its functionals rather than $f$ itself, the monotonicity constraint is not
necessary.

Discrete location mixtures arise as a special case of this model when $f$ is a
step function. Suppose $f(u)=\mu_{h}$ for $u\in[\nu_{h},\nu_{h+1})$, where
$\nu$ is an increasing sequence on $[0,1)$ such that $\nu_{1}=0$ and
$\sum_{h=0}^{\infty}(\nu_{h+1}-\nu_{h})=1$. Then we have:

$p(y)=\int_{0}^{1}\phi_{\sigma}(y-\mu_{h}){\mathbf{1}}{\left(u\in[\nu_{h},\nu_{h+1})\right)}d\,u=\sum_{h=1}^{\infty}(\nu_{h+1}-\nu_{h})\phi_{\sigma}(y-\mu_{h}).$

This representation is intimately related to the augmented model used for
slice sampling infinite mixture models (Walker, 2007; Kalli et al., 2011),
where a prior on $f$ is induced via the prior on mixture component weights
$\pi_{h}\equiv\nu_{h+1}-\nu_{h}$.

Priors on mixture weights are only one way to induce the prior on $f$. Kundu
and Dunson (2014) proposed placing a Gaussian process prior directly on $f$,
suggesting models centered on a prior guess of the quantile function and using
a squared exponential covariance function. Theoretically, this is a flexible
choice (see Pati et al., 2011; Kundu and Dunson, 2014); however, it introduces
computational difficulties, requiring a discretization of the space that may
reduce the quality of the subsequent inference. Furthermore, it is not
immediately clear how to introduce multiple or categorical covariates into
this framework.

The continuous latent variable model is appealing, however. In many contexts
it is more intuitive to think of distributional features as arising from some
omitted or latent continuous variables, and not from heterogeneity due to
multiple independent subpopulations. Adapting BART to this setting yields
priors for $f$ which incorporate covariates flexibly, are approximately
smooth, and do not require discretization of the latent variable _a priori_.

###  3.2 Density Regression with BART (DR-BART)

Figure 1: Representing a step function (left) as a binary tree (right). Figure
2: Modifying the tree in Fig. 1 to incorporate covariates

To motivate the use of tree-based priors, recall that the normal location
mixture model corresponds to a step function for $f$, which can also be
represented as a binary tree (Fig. 1). In the tree-based representation,
incorporating covariates is simply a matter of introducing additional
covariate-based splitting rules; Figure 2 modifies the tree to split on
$\mathbf{x}$. Marginalizing over the latent variable in this example gives the
two conditional densities:

$\displaystyle p(y\mid\mathbf{x}\in A)=$
$\displaystyle\pi_{1}\phi_{\sigma}(y-\theta_{1})+\pi_{2}\phi_{\sigma}(y-\theta_{2})+(\pi_{3}+\pi_{4})\phi_{\sigma}(y-\theta_{3})$
(9) $\displaystyle p(y\mid\mathbf{x}\not\in A)=$
$\displaystyle\pi_{1}\phi_{\sigma}(y-\theta_{1})+\pi_{2}\phi_{\sigma}(y-\theta_{2})+\pi_{3}\phi_{\sigma}(y-\theta_{4})+\pi_{4}\phi_{\sigma}(y-\theta_{5})$
(10)

The resulting model has some interesting properties. There are shared
components and shared weights (e.g. $\theta_{1},\theta_{2}$ and
$\pi_{1},\pi_{2}$) but also components and weights that are unique to each
conditional distribution (e.g. $\pi_{3}+\pi_{4}$ versus $\pi_{3},\pi_{4}$).
This allows borrowing of information across the covariate space, with the
degree of borrowing controlled by the tree structure. Trees with multiple
interior splits on $\mathbf{x}$ allow the model to capture multiscale
structure, as the degree of borrowing varies across the covariate space. In
the extreme case, a tree that splits on $\mathbf{x}$ before splitting on $U$
will yield two independent density estimates, while a tree that doesn’t split
on $U$ yields a standard Bayesian CART model.

In practice, a single tree model will probably be inadequate for many
applications. With many covariates, nearly smooth mean functions, or highly
skewed/multimodal distributions the tree will have to grow quite large.
Additionally, the presence of the latent variable $U$ will tend to yield poor
mixing and further complicate the design of MCMC algorithms for the single
tree model (see Chipman et al. (1998); Denison et al. (1998) for discussion of
the complications arising in MCMC for single tree models _without_ latent
variables).

A BART prior for $f$ is a natural alternative to a single tree model and
additional flexibility can be obtained by modeling the variance function as
well. For density regression we modify both BART priors from Section 2 to
include a latent variable $U$. The most general model for density regression
with BART (DR-BART) is:

$\displaystyle U\sim U[0,1],$ $\displaystyle\quad\epsilon\sim
N(0,\sigma_{0}^{2})$ (11) $\displaystyle Y=f(\mathbf{x},U)$
$\displaystyle+\exp[v(\mathbf{x},U)/2]\epsilon.$ (12)

Letting $\sigma(\mathbf{x},u)=\exp[v(\mathbf{x},u)/2]\sigma_{0}$, the density
function at $\mathbf{x}$ is

$p(y\mid\mathbf{x})=\int_{0}^{1}\phi_{\sigma(\mathbf{x},u)}(y-f(\mathbf{x},u))\,du,$
(13)

a location-scale mixture of normals. In addition to the full DR-BART model,
two reduced models are potentially of interest:

  1. 1.

Location-only mixture (DR-BART-L): Constant bandwidth,
$\sigma(\mathbf{x},u)\equiv\sigma_{0}$

  2. 2.

Location mixture with heteroscedasticity (DR-BART-LH): Covariate-dependent
bandwidth parameter,
$\sigma(\mathbf{x},u)\equiv\sigma(\mathbf{x})=\exp[v(\mathbf{x})/2]\sigma_{0}$

Priors for the DR-BART parameters are specified as follows:

  * •

$f$, the location function: The CGM BART prior with $k=2$ and $m=250$ trees.
We also require that the leaves of each tree contain at least 5 observations.
In the context of DR-BART this condition is related to priors on mixture
models for a single density that require all the components to be occupied
(see e.g. Diebolt and Robert (1994)).

  * •

$v$, the bandwidth function: The BART variance prior from Section 2, with 100
trees and $a_{0}$ (a hyperparameter for $\sigma(\mathbf{x},u)$ defined in the
Supplement) calibrated to a reasonable range, as described in the Supplement.

  * •

Depending on the model, $\sigma_{0}^{2}$ is either fixed (DR-BART, DR-BART-LH)
or given a further inverse gamma prior (DR-BART-L): $\sigma_{0}^{2}\sim
IG(\nu_{0}/2,\nu_{0}\xi_{0}/2)$, with $\xi_{0}$ a prior guess at an
appropriate bandwidth. In DR-BART-L, posterior inferences are insensitive to
reasonable specifications of $\xi_{0}$, but allowing $\nu_{0}\rightarrow 0$
puts too much mass near 0. Since $\sigma_{0}^{2}$ is not identified separately
from the flexible variance function $\nu$ in the full DR-BART/DR-BART-LH
models, the prior should be very informative and in practice fixing it at a
sensible guess (some fraction of the sample variance or of the variance of the
OLS residuals) seems to work well. It could also be elicited more formally.

Results appear to be more or less insensitive to the numbers of trees in $f$
and $v$ provided they are large enough, and the values chosen here reflect
experience with mean regression BART and the belief that variance functions
are less complex than location functions.

Before describing posterior sampling we will describe the properties of $f$
and $v$ and provide some intuition for their roles in the model.

####  3.2.1 The location function $f$

With a BART prior for $f$, $f(\mathbf{x},u)$ decomposes as

$\displaystyle f(\mathbf{x},u)$
$\displaystyle=\sum_{h=1}^{m}g((\mathbf{x},u),T_{h},M_{h})\equiv
f_{x}(\mathbf{x})+f_{xu}(\mathbf{x},u)+f_{u}(u).$ (14)

The three functions in (14) are defined (from left to right) as the sum of the
trees splitting only on $\mathbf{x}$, on both $\mathbf{x}$ and $u$, and only
on $u$. These terms can capture covariate effects that are pure location
shifts, covariate-dependent distributional features, and distributional
features common across covariate space (respectively). Using posterior samples
of the trees to try to infer which variables influence the responses and in
what manner is somewhat fraught, however. For example, it is possible that a
tree in $f_{xu}(\mathbf{x},u)$ split trivially on $U$ in that the leaf
parameters on either side of the split are nearly identical or cancelled by
the contribution of other trees. This is essentially the same difficulty
reported by CGM in doing variables selection in BART mean regression by
counting the trees splitting on a particular variable. But (14) does give some
insight into how the model _can_ capture the complex effects that $\mathbf{x}$
may have on the distribution of $Y$, even with a constant bandwidth.

Since $f$ and $v$ are step functions, this model is equivalent to a discrete
mixture of normal distributions. But the prior is much different than the
usual priors in covariate-dependent mixture models. First, the number of
distinct mixture components with positive probability varies across covariate
space like in the single tree model. Second, unlike in the single tree model,
the components are correlated _a priori_ ; given a fixed set of $m$ trees $T$
we have:

$\displaystyle\mathrm{Cov}(f(\mathbf{x},u),f(\mathbf{x}^{\prime},u^{\prime}))$
$\displaystyle=\sigma^{2}_{\mu}\mathcal{N}[(\mathbf{x},u),(\mathbf{x}^{\prime},u^{\prime})]$
(15)
$\displaystyle\mathrm{Cor}(f(\mathbf{x},u),f(\mathbf{x}^{\prime},u^{\prime}))$
$\displaystyle=\mathcal{N}[(\mathbf{x},u),(\mathbf{x}^{\prime},u^{\prime})]/m$
(16)

where $\mathcal{N}[(\mathbf{x},u),(\mathbf{x}^{\prime},u^{\prime})]$ is the
number of trees where $(\mathbf{x},u)$ and $(\mathbf{x}^{\prime},u^{\prime})$
are in the same leaf. It follows from (16) that if
$|u-u^{\prime}|>|u-u^{\prime\prime}|$, then

$\displaystyle\mathrm{Cov}(f(\mathbf{x},u),f(\mathbf{x},u^{\prime}))$
$\displaystyle\geq\mathrm{Cov}(f(\mathbf{x},u),f(\mathbf{x},u^{\prime\prime}))$
(17) $\displaystyle\mathrm{Cor}(f(\mathbf{x},u),f(\mathbf{x},u^{\prime}))$
$\displaystyle\geq\mathrm{Cor}(f(\mathbf{x},u),f(\mathbf{x},u^{\prime\prime})).$
(18)

Given the potentially strong correlation, it is misleading to think of the
steps in $f$ as “mixture components” in the usual sense.

####  3.2.2 The bandwidth function $v$

A covariate-adaptive bandwidth parameter will often be important in these
models. Since $U$ and $\epsilon$ are independent, in the DR-BART-L model
$\mathrm{Var}(Y\mid\mathbf{x})\geq\mathrm{Var}(\sigma_{0}\epsilon)=\sigma_{0}^{2}$.
On the other hand, in the DR-BART-LH model
$\mathrm{Var}(Y\mid\mathbf{x})\geq\mathrm{Var}(\exp[v(\mathbf{x})/2]\epsilon)=\sigma(\mathbf{x})^{2}$.
In a model without a covariate adaptive bandwidth, $\sigma^{2}_{0}$ must be at
least as small as the most concentrated predictive density to avoid
oversmoothing, yielding much rougher density estimates elsewhere. The
efficiency of the MCMC sampler suffers as well, as smaller bandwidths imply
more concentrated distributions for $U_{i}$. A covariate dependent bandwidth
might be preferable for this reason, even if a single bandwidth seems like a
reasonable simplification.

The full DR-BART model has an additional degree of flexibility due to scale
mixing. It can allow $v(\mathbf{x},u)$ to grow large in some areas of of
$U$-space, effectively “turning off” portions of $f(\mathbf{x},u)$ or
capturing relatively flat areas of the density. The behavior of the different
models is easiest to understand with an example, presented in Section 5.

####  3.2.3 Posterior Sampling

Generating samples from the posterior with MCMC is straightforward.
Conditional on values for the latent variables $u=(u_{1},u_{2},\dots,u_{n})$,
DR-BART reduces to the heteroscedastic BART model so that sampling for the
other parameters proceeds as described in the Supplement. The latent $u_{i}$
have full conditionals

$p(u_{i}\mid-)\propto\phi_{\sigma(\mathbf{x},u)}(y-f(\mathbf{x},u)){\mathbf{1}}{\left(u\in\mathcal{B}(T,T^{(v)})\right)},$
(19)

where $\mathcal{B}(T,T^{(v)})$ is the set of possible $u$-values; that is,
those that do not yield trees with leaves having fewer than 5 observations.
Since $f$ and $\sigma$ are step functions, (19) is piecewise constant, so
$u_{i}$ can be updated with a Gibbs step: If
$u^{*}_{1}<u^{*}_{2}<\dots<u^{*}_{k}$ are the points of discontinuity of (19)
and $u^{*}_{0}=0,\ u^{*}_{k-1}=1$, the Gibbs step first samples an interval
from

$\Pr(u_{i}\in(u^{*}_{h},u^{*}_{h+1})\mid-)\propto\phi_{\sigma(\mathbf{x},\tilde{u}_{h})}(y-f(\mathbf{x},\tilde{u}_{h})){\mathbf{1}}{\left(u\in\mathcal{B}(T,T^{(v)})\right)},$
(20)

where $\tilde{u}_{h}=(u^{*}_{h+1}-u^{*}_{h})/2$ (or any other point in the
interval) and then $u_{i}$ is sampled uniformly from the selected interval.
Note that ${\mathbf{1}}{\left(u\in\mathcal{B}(T,T^{(v)})\right)}$ necessarily
equals either 0 or 1 on the entire interval $(u^{*}_{h},u^{*}_{h+1})$ for each
$h$ by construction.

While conceptually simple, the Gibbs sampling update can require a large
number of likelihood evaluations ($k$ is often well into the hundreds). On the
other hand, a Metropolis step is difficult to tune because (19) is in general
multimodal. An efficient alternative that doesn’t require tuning is slice
sampling which introduces latent variables $\omega_{i}$ so that

$p(\omega_{i},u_{i}\mid-)\propto{\mathbf{1}}{\left(\omega_{i}<\phi_{\sigma(\mathbf{x},u)}(y-f(\mathbf{x},u)\right)}{\mathbf{1}}{\left(u\in\mathcal{B}(T,T^{(v)})\right)}.$
(21)

Sampling proceeds using the techniques developed in Neal (2003). The slice
sampler is much more efficient, which tends to make up for any loss in
theoretical efficiency or mixing.

##  4 Theory

Here, we present some properties of the DR-BART model, showing that our
proposed prior generates trees that are almost surely finite and upper
bounding the rate at which DR-BART estimates the conditional density
$p(y\mid\mathbf{x})$. We focus on the special case where the predictors
$\mathbf{x}_{i}\in[0,1]^{p}$ are continuous. This is a common assumption when
studying theoretical properties of density regression models (e.g. Chung and
Dunson, 2012; Pati et al., 2013; Li et al., 2020). The assumption is even more
innocuous here, as BART is invariant to monotone transformations of the
covariates. Proofs can be found in the Supplement.

######  Theorem 1.

Assume that the prior over a binary tree is as in CGM, but with a continuous
uniform prior on splitting locations and no restriction to nonempty leaves. A
tree sampled from this prior has finite depth with probability one.

Thus, introducing a latent $u$ into the prior does not affect the finite depth
of the trees. We now focus on providing upper bounds for the posterior
concentration rate $\epsilon_{n}$ of the posterior. A rate
$\epsilon_{n}\downarrow 0$ is said to be a _rate of convergence_ of the
posterior with respect to a divergence measure $h$ if there exists a positive
constant $M$ such that $\Pi\\{h(p_{0},p_{f,\sigma})\geq
M\epsilon_{n}\mid\mathcal{D}_{n}\\}\to 0$ in $F_{0}$-probability, where
$\mathcal{D}_{n}=\\{(\mathbf{X}_{i},Y_{i}):i=1\ldots,n\\}$ and
$(\mathbf{X}_{i},Y_{i})\overset{iid}{\sim}F_{0}(d\mathbf{x},dy)=p_{0}(y\mid\mathbf{x})\
dy\ F_{\mathbf{X}}(d\mathbf{x})$.

The conditional density of $y$ induced by our model is given by the
convolution

$\displaystyle
p_{f,\sigma}(y\mid\mathbf{x})=\int\phi_{\sigma(\mathbf{x},u)}(y-f(\mathbf{x},u))\
du,$

and the limit $\sigma(\mathbf{x},u)\to 0$ is associated with a random variable
with quantile function $f(\mathbf{x},u)$ if $f(\mathbf{x},u)$ is monotonically
increasing in $u$. This suggests that a reasonable strategy for establishing
that $p_{f,\sigma}$ is close to $p_{0}$ is to show that $f(\mathbf{x},u)$ is
close to the true conditional quantile function $f_{0}(\mathbf{x},u)$.

We characterize the concentration of DR-BART with respect to the integrated
Hellinger distance, defined:
$h(p,q)=\big{(}\int(\sqrt{p(y|\mathbf{x})}-\sqrt{q(y|\mathbf{x})})^{2}\ dy\
F_{\mathbf{X}}(d\mathbf{x})\big{)}^{1/2}$. Two other divergence measures which
will be useful for us are the _generalized Kullback-Leibler divergences_
$\operatorname{KL}(p_{1}\|p_{2})=\int p_{1}\log(p_{1}/p_{2})\ dy\
F_{\mathbf{X}}(d\mathbf{x})$ and $V(p_{1}\|p_{2})=\int
p_{1}\log^{2}(p_{1}/p_{2})\ dy\ F_{\mathbf{X}}(d\mathbf{x})$.

To study the posterior concentration of DR-BART, we make use of results from
(i) Jeong and Rockova (2020) involving the concentration of BART in a
regression setting and (ii) Pati et al. (2011); Zhou et al. (2017) who
leverage similar results about Gaussian processes to show convergence in a
latent variable model similar to the one considered here. Our proof extends
both of these works in fundamental ways. We extend Jeong and Rockova (2020) by
introducing latent variables into the tree structure. While having immediate
implications for DR-BART this also lays the groundwork for concentration
results for latent variable BART models that we will consider in the future.
Compared to Pati et al. (2011), we not only introduce covariates into the
setup, but also place a non-trivial DR-BART prior over the bandwidth, compared
to their choice of a parametric, covariate-independent prior.

Next, we outline conditions required for our proof. Throughout, we will write
$a\lesssim b$ to mean that there exists a positive constant $C$, possibly
depending on $p_{0}$ and on hyperparameters but otherwise independent of $n$,
$p$, or any other variables, such that $a\leq Cb$.

##### Condition F (on $p_{0}$)

We assume that $\log p_{0}(y\mid\mathbf{x})$ is $\alpha$-Hölder as a map from
$[0,1]^{p+1}$ to $\mathbb{R}$ for some $0<\alpha\leq 2$. Additionally, we
assume $\log p_{0}$ is $d_{0}$-sparse in the sense that it depends on
$(y,\mathbf{x})$ only through the coordinates in
$S_{0}\subseteq\\{1,\ldots,p+1\\}$ where $|S_{0}|=d_{0}$, $d_{0}=o(\log(n))$,
and $d_{0}\log p=o(n)$. Lastly, we assume that
$||f_{0}||_{\infty}\lesssim\sqrt{\log n}$.

######  Remark 1.

The assumption that $\log p_{0}$ is $\alpha$-Hölder implies that $p_{0}$ is
bounded and bounded away from $0$. Condition F is used both to ensure that
$p_{0}$ can be well-approximated with convolutions and to ensure that
$\|f-f_{0}\|_{\infty}$ is small with sufficiently large prior probability.

##### Condition P (on $\Pi$)

Let $S\subseteq\\{1,\ldots,p+1\\}$ denote the coordinates of $\mathbf{x}$
which the trees $f(\mathbf{x},u)$ and $v(\mathbf{x},u)$ split on.

  * (P1)

The support set $S$ of $(f,\sigma)$ has prior
$\pi(S)=\binom{p+1}{D}^{-1}\pi_{D}(D)$ where $D\equiv|S|$ and $\pi_{D}(d)$ is
an _exponentially decaying prior_ satisfying
$a_{1}(p+1)^{-a_{3}}\pi_{D}(d-1)\leq\pi_{D}(d)\leq
a_{2}(p+1)^{-a_{4}}\pi_{D}(d-1)$ for some positive constants
$a_{1},a_{2},a_{3},a_{4}$ and $d=1,\ldots,p+1$.

  * (P2)

Given $S$, each tree $T_{h},h=1,\ldots,m$ and $T_{h}^{(v)},h=1,\ldots,m_{v}$
is assigned the branching process prior with splitting proportion
$q(d)=\nu^{d}$ for some $\nu\in(0,1/2)$.

  * (P3)

The leaf node parameters $\mu_{hl}$ of $T_{h}$ are assigned independent
$N(0,\sigma^{2}_{\mu})$ priors.

  * (P4)

The log-variance function is given by
$v(\mathbf{x},u)=\xi+\sum_{h=1}^{m_{v}}g(\mathbf{x},u;T_{h}^{(v)},M_{h}^{(v)})$
where $e^{-\xi}\sim\operatorname{Gam}(a_{\sigma},b_{\sigma})$ and
$\mu_{hl}^{(v)}\overset{iid}{\sim}\pi_{v}$ where $\pi_{v}$ is a strictly
positive density supported on an interval $[-V,V]$.

  * (P5)

Splits in the tree ensemble can occur only at a number $b_{n}$ of candidate
split points $\mathcal{Z}_{n}\subseteq[0,1]^{p+1}$, which are selected from
uniformly. Additionally, $\log b_{n}\lesssim\log n$.

  * (P6)

For each $n$ there exists a decision tree $\widehat{T}$ and leaf node values
$\widehat{M}$ built from the candidate split-points in $\mathcal{Z}_{n}$ such
that the regression tree
$f^{\star}(\mathbf{x},u)=g(\mathbf{x},u;\widehat{T},\widehat{M})$ satisfies
$\|f_{0}-f^{\star}\|_{\infty}\lesssim(\log n/n)^{\beta/(2\beta+d_{0})}$ where
$\beta=\min\\{\alpha,1\\}$.

######  Remark 2.

Though P4 constrains the support of $\pi_{v}$, $\xi$ is unbounded, allowing
the variance function to have arbitrary scale even as the trees grant
arbitrary flexibility in its shape.

######  Remark 3.

The only assumption which is seemingly beyond our direct control is P6, which
asserts that $f_{0}$ can be uniformly approximated with a single decision tree
using the candidate split points $\mathcal{Z}_{n}$. Jeong and Rockova (2020)
give several valid configurations of split points for which P6 would hold; for
example, when $\mathcal{Z}_{n}$ is a regular grid of size $b_{n}\asymp n^{cp}$
for $c$ a sufficiently large constant if Condition F holds.

Under these assumptions, we have the following theorem:

######  Theorem 2.

Assume that Condition F and Condition P hold. Then, there exists a positive
constant $M$ such that $\Pi\\{h(p_{0},p_{f,\sigma})\geq
M\epsilon_{n}\mid\mathcal{D}_{n}\\}\to 0$ in $F_{0}$-probability, where
$\epsilon_{n}=(n/\log
n)^{-\frac{\alpha}{\alpha+1}\times\frac{\beta}{2\beta+d_{0}}}+\sqrt{d_{0}\log(p+1)/n}$,
where $\beta=\min\\{\alpha,1\\}$.

######  Remark 4.

The gap between what is attainable by DR-BART-LH is larger than might be
expected, as the rate is removed by a factor of $\alpha/(\alpha+1)$ from the
minimax rate. The reason this occurs is that Condition F implies
$f_{0}(\mathbf{x},u)$ has Hölder smoothness of $\alpha+1>1$ in $u$, whereas
BART is not known to be able to adapt to smoothness levels higher than 1. If
we modify Condition F to state that $f_{0}(\mathbf{x},u)$ is $\alpha+1$-smooth
as a function of $(u,\mathbf{x})$ (as opposed to just in $u$) then it is
possible to show that replacing the BART model with the SBART model of Linero
and Yang (2018) gives the rate
$n^{\alpha/(2\alpha+d_{0})}\log(n)^{-\alpha(d_{0}+1)/(2\alpha+d_{0})}+\sqrt{d_{0}\log(p+1)/n}$
adaptively over $\alpha$ and $S_{0}$ for $\alpha\in(0,2]$; extending these
results to higher $\alpha$ using results of Plummer et al. (2021) is deferred
to future work.

As argued by Li et al. (2020), the posterior rate of convergence is (bounded
by) a sequence $\epsilon_{n}$ with $n\epsilon_{n}^{2}\to\infty$ if we can find
positive constants $C_{1},\ldots,C_{4}$ such that, for every sufficiently
large $n$, there exists a set $\mathcal{G}_{n}$ of conditional densities
satisfying the following:

  * (G1)

$\Pi\\{p_{f,\sigma}\in\operatorname{KL}_{p_{0}}(C_{1}\epsilon_{n})\\}\geq\exp\\{-C_{2}n\epsilon_{n}^{2}\\}$,
with
$\operatorname{KL}_{p_{0}}(\epsilon)=\\{p:\operatorname{KL}(p_{0}\|p),V(p_{0}\|p)\leq\epsilon^{2}\\}$

  * (G2)

$\Pi(\mathcal{G}_{n}^{c})\leq C_{3}\exp\\{-(C_{2}+4)n\epsilon_{n}^{2}\\}$.

  * (G3)

$\log N(\mathcal{G}_{n},\bar{\epsilon}_{n},h)\leq C_{4}n\epsilon_{n}^{2}$
where $\bar{\epsilon}_{n}$ is a constant multiple of $\epsilon_{n}$ and
$N(\mathcal{G}_{n},\epsilon,d)$ denotes the $\epsilon$-covering number of
$\mathcal{G}_{n}$ with respect to $d$ (see, e.g. Ghosal et al., 2000).

The following two lemmas proved in the course of establishing Theorems 2 and 5
of Jeong and Rockova (2020) play a key role in establishing our results. The
first ensures the prior on $f$ places sufficient mass around $f_{0}$ and is
essential in establishing G1.

######  Lemma 1.

Suppose that Condition F and Condition P hold and let $\delta_{n}=(\log
n/n)^{\beta/(2\beta+d_{0})}$. Then for sufficiently large $n$ we have

$\displaystyle-\log\Pi(\|f-f_{0}\|_{\infty}\leq\delta_{n}\mid S=S_{0})\lesssim
n\delta_{n}^{2}.$

The second lemma ensures that the support of the BART prior is “small” in a
suitable sense, allowing us to verify G2 and G3.

######  Lemma 2.

Let $\mathcal{F}$ denote the collection of decision tree ensembles with $m$
trees which (i) split on no more than $d$ variables, (ii) have at most $K$
leaf nodes per tree, (iii) have at most $b_{n}$ candidate split points, and
(iv) satisfy $\sup_{hl}|\mu_{hl}|\leq U$. Then

$\displaystyle\log N(\mathcal{F},\epsilon,\|\cdot\|_{\infty})\lesssim
d\log(p+1)+K\log\left(\frac{d^{m}b_{n}^{m}KU}{\epsilon}\right).$

Finally, Lemma 3 connects $h(p,q)$, $\operatorname{KL}(p\|q)$, and $V(p\|q)$
to the supremum norm.

######  Lemma 3.

Suppose that Condition F holds. Then there exists a constant
$C_{\operatorname{KL}}$ independent of $(n,p)$ such that, for sufficiently
small $\epsilon$, we have

$\displaystyle\operatorname{KL}_{p_{0}}(C_{\operatorname{KL}}\epsilon)\supseteq\\{\sigma(\mathbf{x},u)\equiv\sigma\textnormal{
is
constant},\sigma\in(\epsilon^{1/\alpha},2\epsilon^{1/\alpha}),\|f-f_{0}\|_{\infty}\leq\epsilon^{1+1/\alpha})\\}.$

Additionally, for any bounded measurable functions
$f_{1},f_{2},\log\sigma_{1},\log\sigma_{2}:[0,1]^{p+1}\to\mathbb{R}$:

$\displaystyle
h(p_{f_{1},\sigma_{1}},p_{f_{2},\sigma_{2}})\lesssim\sqrt{\|\log\sigma_{1}-\log\sigma_{2}\|_{\infty}}+\frac{\|f_{1}-f_{2}\|_{\infty}}{\inf_{\mathbf{x},u}\sigma_{1}(\mathbf{x},u)\wedge\sigma_{2}(\mathbf{x},u)}.$

Straightforward application of Lemmas 1 – 3 suffices to verify G1 – G3. The
first part of Lemma 3 allows us to decompose the probability of the KL ball in
G1 into pieces that can be bounded by Lemma 1. To verify G2, we consider a
sieve defined by individual sieves for $f$, $v$, and $\xi$; the second part of
Lemma 3 allows us to bound the Hellinger distance between a conditional
density and an element in this sieve and Lemma 2 then ensures that the entropy
is appropriately bounded. Given the sieve, G3 follows given Condition P.

##  5 Simulations

Here, we evaluate how well DR-BART and other methods estimate conditional
densities and appropriately express uncertainty about these estimates. We do
so via variants of a challenging univariate example. Additional simulation
details can be found in the Supplement.

###  5.1 Simulation 1: Contrasting DR-BART Models

Here, we introduce our basic simulation setup and gain insight into the
differences between DR-BART-L, DR-BART-LH, and the full DR-BART model, before
comparing to other methods. We consider a challenging example with a single
regressor:

$\displaystyle
X_{i}\overset{iid}{\sim}U(0,1),\,Y_{i}=f_{0}(X_{i})+\epsilon_{i}(X_{i})$ (22)

where

$f_{0}(\mathbf{x})=5\exp[15(x-0.5)]/(1+\exp[15(x-0.5)])-4x$

and $\epsilon_{i}(X_{i})$ is given by

$\displaystyle p(\epsilon\mid
X=x)=\lambda(x)N(\epsilon;2x-0.6,0.3^{2})+(1-\lambda(x))p_{G}(\exp(\epsilon),0.5+x^{2},1.0)\exp(\epsilon),$

where $\lambda(x)=\exp[-10(x-0.8)^{2}]$ and the second component of
$p(\epsilon\mid x)$ is a log-Gamma distribution with scale $1$ and shape
$0.5+x^{2}$. Figure 3 shows selected quantile processes and conditional
densities. The log-Gamma component is skewed and heteroscedastic, and the
normal component is much more concentrated. The conditional distributions are
nearly all unimodal, but for $x$ values around 0.4 – 0.6 the density is quite
peaked around the mode with a heavy left tail. At $x=0.8$, $\epsilon$ has
exactly a $N(1,0.3^{2})$ distribution.

Figure 3: (Left) Conditional quantiles of the example in Section 5. Bands are
equal-tailed 98%, 90%, and 50% intervals as a function of $x$. The solid line
is the median process. (Right) The conditional pdf for selected $\mathbf{x}$
values.

Representative results for a single dataset are shown in Figure 4, which shows
the estimated predictive density at $x=0.1,0.5,0.8.$ All three models do well
when $x=0.8$; BART captures the nonlinear mean function particularly well
since the noise is low and symmetric, and the variance is well-estimated since
a peak of approximately this width is present in many of the conditional
densities. However, DR- BART-L predictably struggles at $x=0.1$ when the peak
disappears and the density becomes much more diffuse and skewed. Some of this
severe multi-modality can be mitigated by adjusting the prior to include more
/ larger trees. But the fundamental problem is with the single bandwidth
parameter: it must be low to capture the peak, which means that the spread in
$p(y|x=0.1)$ has to be captured by $f(0.1,u)$ varying greatly in $u$, yielding
rougher densities.

Figure 4: Representative true (dotted) and estimated (solid) densities from
Simulation 1.

DR-BART-LH and DR-BART perform much better, with the full model doing slightly
better, particularly when $x=0.1$, since it is a location-scale mixture and
can capture the long left tail by splitting on $u$ in the variance function to
add high-variance “components”. The differences are fairly modest though, and
the increase in computation time is $\scriptstyle\sim$25%. Given its
comparable performance, DR-BART-LH is a viable alternative to the full model,
especially when most of the conditional densities are not severely multimodal
or skewed.

###  5.2 Competing Methods

We now compare the full DR-BART model to a variety of other methods for
conditional density estimation. First, we compare to the Probit Stick-Breaking
Process Mixture (PSBPM) in Chung and Dunson (2012). Chung et al. model the
conditional density $p(y|\mathbf{x})$ via an infinite mixture of normal linear
regressions according to mixing distributions $P_{x}$:

$p(y|\mathbf{x})=\int\phi_{\tau}(y-\mathbf{x}^{\prime}\beta)\,dP_{\mathbf{x}}$

where the prior over the $P_{\mathbf{x}}$ is defined by a covariate-dependent
stick-breaking prior.

Second, we compare to the Soft BART Density Sampler (SBART-DS) of Li et al.
(2020) which models conditional densities by modulating a base model
$h(y|\mathbf{x},\theta)$ via a link function $\Phi(\mu):$
$p(y|\mathbf{x},\theta)\propto{h(y|\mathbf{x},\theta)\Phi\\{r(y,\mathbf{x})\\}}$.
There are many possible choices for $h$, $r$, and $\Phi$; Li et al. choose to
center their conditional densities on a normal linear regression:
$h(y|\mathbf{x},\theta)=\phi_{\sigma_{\theta}}(y-\mathbf{x}^{\prime}\beta_{\theta})$
and $\Phi$ to be a probit link. They take $r$ to be the weighted sum of Soft
Additive Regression Trees (Linero and Yang, 2018) with random Fourier
expansions approximating a Gaussian process in the leaf nodes, as in Starling
et al. (2020).

Third, we compare to a Dirichlet Process Mixture Model (DPMM) of Jara et al.
(2011) that models $(\mathbf{x},y)$ as jointly normal and computes the implied
conditional of $y$ on $\mathbf{x}$.

Lastly, we compare to a generalization of Zhou et al. (2017) that incorporates
covariates into the transfer function $\mu$, which has a Gaussian process
prior (DR-GP). This approach is similar to ours in its use of a latent $u$ to
perform conditional density estimation. It differs in its use of 1. a Gaussian
process to map $(u,\mathbf{x})$ to the observed response and 2. a
homoscedastic, inverse Gamma prior on $\sigma$, instead of our BART priors on
each.

###  5.3 Simulation 2: Univariate Example

Here, we consider the same simulation as above. We simulate $n=800$ data
points from model (22) and fit all methods. Fits and credible intervals for a
representative simulation are shown in Figure 5. We see that DR-BART estimates
the mean well with appropriate uncertainty quantification across all
conditional densities, though it struggles to fully capture the peak of the
normal density at $x=0.8$ given how concentrated it is. Most the other methods
are able to fit the densities at $x=0.1$ and $x=0.8$ reasonably well, but have
difficulty with capturing both the peak and strong skew present in $p(y\mid
x=0.5)$. Results aggregated across 100 simulations are shown in Table 1. As
suggested by Figure 5, DR-BART outperforms all other methods when
$x\in\\{0.1,0.5\\}$, but falls short of always capturing the peak at $x=0.8$.
DPMM estimates $p(y\mid x=0.8)$ particularly well, which is to be expected
given that it is built upon a normal specification.

We next assess how well the 95% credible bands of each method cover the true
$p(y\mid x)$. Table 1 also displays the proportion of simulations where the
credible bands fully contain the true density within the true 95% HDR
interval. As Figure 5 suggests, DR-BART is able to do so consistently for
$x=0.1$ and $x=0.5$, but not always for $x=0.8$. For this and the next two
simulations, the Supplement also contains information on credible band width
and predictive coverage, which was consistently close to nominal for all
methods but PSBPM.

Figure 5: Representative true (solid) and estimated (dashed) densities from
Simulation 2.

DR-BART DPMM SBART-DS DR-GP PSBPM Error Coverage Error Coverage Error Coverage
Error Coverage Error Coverage $x=0.1$ 1 0.93 1.67 0.75 1.65 0.48 2.74 0.02
9.80 0.00 $x=0.5$ 1 0.61 1.01 0.62 2.20 0.00 1.46 0.00 2.14 0.00 $x=0.8$ 1
0.19 0.28 0.85 0.79 0.00 1.03 0.00 2.27 0.00

Table 1: Error and Coverage for Simulation 2. Error is Wasserstein distance
between true and estimated densities (normalized at each $x$ for
interpretability), averaged over 100 simulations. Coverage is the proportion
of simulations in which 95% credible bands cover the truth within the true 95%
HDR region.

###  5.4 Simulation 3: Irrelevant Covariates

Here, we generate data as above, but supply an additional 14 irrelevant,
uniformly distributed variables, with pairwise correlations of 0.3, to each
method. For each method, we evaluate the predictive density with the
irrelevant variables fixed to $0.5$. Figure 6 and Table 2 illustrate the
results. While the quality of DR-BART’s fit decreases, it compensates by
increasing the width of its credible bands to capture the truth. DPMM and DR-
GP, which in the previous simulation fit better at $x=0.1,0.8$ do not,
however, adjust their estimates of uncertainty to counter the substantial
degradation in their fits. Perhaps most notable is that DR-BART still has good
pointwise coverage, whereas the other methods are overconfident in their
overfitting to the irrelevant covariates. We also note that the SBART-DS
credible bands are not more erratic than in Simulation 1, reflecting BART’s
ability to perform variable selection.

Figure 6: Representative true (solid) and estimated (dashed) densities from
Simulation 3.

DR-BART DPMM SBART-DS DR-GP PSBPM Error Coverage Error Coverage Error Coverage
Error Coverage Error Coverage $x=0.1$ 1 0.78 1.10 $0.00$ 0.78 0.46 2.63 0.00
5.58 0.00 $x=0.5$ 1 0.94 2.96 0.00 1.97 0.00 2.23 0.00 1.97 0.00 $x=0.8$ 1
1.00 1.60 0.00 0.56 0.02 2.52 0.00 1.70 0.00

Table 2: Error and Coverage for Simulation 3. Error is Wasserstein distance
between true and estimated densities (normalized at each $x$ for
interpretability), averaged over 50 simulations. Coverage is the proportion of
simulations in which 95% credible bands cover the truth within the true 95%
HDR region.

###  5.5 Simulation 4: Insufficient Data

Next, we explore the ability of DR-BART to appropriately express uncertainty
about conditional distributions in regions of $x$-space where there is little
data. To do so, we simulate $y$ given $x$ as in the previous experiments;
however, instead of simulating $x$ uniformly, we simulate it from a mixture of
uniforms: $0.475\times U_{[0,0.4]}+0.05\times U_{[0.4,0.6]}+0.475\times
U_{[0.6,1]}$, where $U_{[a,b]}$ denotes the density of a $U(a,b)$ random
variable. With only 5% of the mass lying between $x=0.4$ and $x=0.6$,
estimation of $p(y\mid x=0.5)$ is much harder. With this experiment, we expect
the performance of all methods to degrade; we are interested in whether their
uncertainty grows appropriately to account for the lack of data.

Figure 7: Representative true (solid) and estimated (dashed) densities from
Simulation 4.

Figure 7 shows a representative simulation. Most of the methods perform better
in estimating $p(y\mid x=0.8)$, as there is now substantially more data in
that region of the space. At $x=0.5$, only DR-BART and DPMM increase their
uncertainty in their estimation; the other methods have relatively poor fits
with insufficient uncertainty about them. Table 3 summarizes this information
across 50 simulations. DR-BART performs better in predicting $p(y\mid x=0.5)$,
but the relative performance for $x\in\\{0.1,0.8\\}$ is comparable to that in
Section 5.3. Table 3 also supports the observation that DR-BART and DPMM are
unique in increasing the uncertainty in their estimates in response to the
lack of data.

DR-BART DPMM SBART-DS DR-GP PSBPM Error Coverage Error Coverage Error Coverage
Error Coverage Error Coverage $x=0.1$ $1$ $0.96$ $1.55$ $0.94$ $1.59$ $0.36$
$2.87$ $0.02$ 10.35 0.00 $x=0.5$ $1$ $0.54$ $1.34$ $0.56$ $1.91$ $0.02$ $1.34$
$0.00$ 1.26 0.00 $x=0.8$ $1$ $0.16$ $0.29$ $0.80$ $0.79$ $0.02$ $1.01$ $0.00$
2.69 0.00

Table 3: Error and Coverage for Simulation 4. Error is Wasserstein distance
between true and estimated densities (normalized at each $x$ for
interpretability), averaged over 50 simulations. Coverage is the proportion of
simulations in which 95% credible bands cover the truth within the true 95%
HDR region.

###  5.6 Simulation 5: Comparison of BART Models

Here, we further compare DR-BART to SBART-DS. Because both are built off of
(S)BART models, we expect them to share many desirable properties, such as
high flexibility and the ability to identify interaction effects and perform
variable selection. SBART-DS, however, is centered on a base model, which may
degrade the fit when the truth is far from the base model. Furthermore, there
are cases in which prior information suggesting a reasonable base model may be
lacking. To explore reliance on the base model, we run the original univariate
simulation, but change the mean to be $f_{0}(x)=a(x-0.5)^{2}$. For small $a$,
the mean is approximately linear for $x\in[0,1]$, as specified by SBART-DS.
But as $a$ increases, it becomes increasingly nonlinear and we expect
performance to degrade. Using $n=1500$, we run 50 simulations for each
$a\in\\{0,1,5,15,25,35\\}$, summarized in Table 4. For small $a$, both methods
perform well. But the increasing nonlinearity of the mean as $a$ increases
hinders SBART-DS’ performance deteriorates, while DR-BART is unaffected.

$a=0$ $a=1$ $a=5$ $a=15$ $a=25$ $a=35$ DR-B S-DS DR-B S-DS DR-B S-DS DR-B S-DS
DR-B S-DS DR-B S-DS $x=0.1$ 0.39 0.56 0.20 0.55 0.37 0.50 0.39 0.62 0.82 0.97
0.37 1.04 $x=0.5$ 0.76 1.09 0.74 0.85 0.89 0.91 0.39 1.50 0.71 1.57 0.80 2.20
$x=0.8$ 1.11 0.94 1.32 0.96 1.05 1.31 1.05 2.53 1.26 2.90 1.04 3.25

Table 4: Wasserstein distance between true, estimated densities in Simulation
5, averaged over 50 simulations. As $a$ increases and the true mean becomes
increasingly nonlinear, SBART-DS (S-DS) is hindered by its reliance on a
linear base model, while DR-BART (DR-B) is unaffected .

##  6 Applications

We consider two applications of DR-BART. In both cases, we fit DR-BART-LH; the
full DR-BART model gave similar results since both applications involve fairly
well-behaved densities. For all models, $\sigma_{0}^{2}$ is set to be half the
standard deviation of the OLS residuals and results seem to be insensitive to
reasonable choices of this value. We also set $a_{0}=\log(\sqrt{4})^{-2}$. For
each model, 10,000 MCMC samples are collected after a burn-in of 10,000
iterations. The estimands of interest in both sections are functions of
quantiles; we will use $Q(s\mid\mathbf{x})$ to denote the quantile function of
$P(y\mid\mathbf{x})$.

###  6.1 Calculating Student Growth

We consider the estimation of student growth from a series of test scores,
using anonymized mathematics test scores provided in the SGP R package
(Betebenner et al., 2011). Castellano and Ho (2013) gives an overview of
current mean and quantile regression approaches to this problem, as well as
some consideration to the substantive problem of measuring student growth.
Statistically, the problem reduces to estimating a series of conditional
densities $p(y_{t}\mid y_{t-1},y_{t-2},\dots,y_{0})$. These distributions are
heteroscedastic and tend to become skewed as they approach the extremes. The
current state of the art is the quantile regression methodology implemented in
the SGP package and detailed in Betebenner (2009). The quantile regression
models are specified as

$Q_{y_{t}}(s\mid y_{t-1},\dots y_{0})=\sum_{r=0}^{t-1}\eta_{r}^{(s)}(y_{r}).$
(23)

where each function $\eta_{r}^{(s)}(\cdot)$ is given a B-spline basis
expansion. This model is fit at each of the 99 percentiles to approximate the
conditional distribution. This procedure does not yield valid estimates for
percentiles of the conditional distribution in general, due to quantile
crossing (though this can be corrected post-hoc (Chernozhukov et al., 2010)).
Assessing uncertainty in this framework is challenging as well, and analysis
generally relies solely on point estimates. Bayesian density regression
addresses both problems simultaneously.

We fit independent DR-BART-LH models to scores from grades 4-7, each
conditional on all the previous scores as well as the grade 3 scores. For
simplicity we took a subsample of 3,000 students with complete data from
grades 3-7. An interesting feature of this data is the interactions between
previous test scores on the predictive distributions. One mechanism for this
is regression to the mean: large jumps in test scores are unlikely to be
sustained. Figure 8 gives an example, displaying posterior predictive
densities for grade 5 scores when fixing grade 3 and 4 scores at all
combinations of their marginal quartiles. In the absence of interactions the
shift between the pairs of solid and dashed densities would be equal. Figure 9
further suggests the presence of an interaction effect. When $y_{4}=455$,
$y_{3}$ has a smaller effect on the quantiles of $p(y_{5}\mid y_{4},y_{3})$
than when $y_{4}=538$. Note that the additive model in (23) cannot capture
such an interaction and requires the curves to be equal.

Figure 8: Four posterior predictive distributions for grade 5 scores obtained
by setting grade 3 and 4 scores at each combination of their marginal
quartiles. The effect of the grade 3 score on the predictive distribution
clearly depends on the value of the grade 4 score, an important interaction
that DR-BART is well-equipped to capture. Figure 9: Quantile effects
$Q_{y_{5}}(s\mid y_{4},y_{3}=523)-Q_{y_{5}}(s\mid y_{4},y_{3}=425)$ for
$y_{4}=455,538$ along with 90% credible intervals. Under Betebenner (2009)’s
additive model these two curves are forced to be equal. This reveals a likely
interaction between grade 3 and grade 4 scores: the posterior mean of the 538
curve is nearly outside the 455 credible interval.

A key objective in student growth modeling is computing growth quantile
targets. A growth quantile is a quantile of the predictive distribution for
the current test score given the score history. Intuitively, it provides a
measure of how well a student performed on the test relative to academic peers
(the hypothetical population of students with identical test history). A
growth quantile target is the level of consecutive quantile growth that would
be required to achieve pre-established achievement targets. For these data,
there are four achievement levels used to define achievement targets:
Unsatisfactory, Partially Proficient, Proficient, or Advanced. Growth quantile
targets answer questions like “What level of sustained growth is necessary for
a student with grade 3 test score $y_{3}$ to be Proficient by grade 7?”. This
is different than simulating score trajectories and computing the probability
that a student reaches a target. Growth quantile targets are intended to
promote a “what will it take” attitude over a fatalistic “where will s/he be”
attitude (Betebenner, 2011).

The current methodology for computing growth quantile targets uses a series of
point estimates, ignoring uncertainty in the quantile curve estimates. But
this uncertainty is generally not negligible, particularly for students with
extreme test scores. To illustrate, we computed posterior samples of growth
quantile curves for a hypothetical third grade student who scored at the cusp
of Unsatisfactory/Partially Proficient on the third grade math test. Figure 10
plots this student’s growth quantile curves and 90% credible intervals for
grades 4-7. It is easy to read off growth quantile targets from these charts:
For example, the bottom right panel of Figure 10 shows that if the student
sustains $60^{th}$ percentile growth in grades 5,6 and 7, there is a 95%
probability of remaining partially proficient in grade 7. Having a 50% chance
of reaching proficiency in grade 7 would require sustained growth at the
$80^{th}$ percentile, which is probably unattainable without some new
intervention. Similar curves can be computed starting with any grade and any
score history.

Figure 10: Projected scores as a function of quantile growth for a student at
the cusp of Unsatisfactory/Partial Proficiency in grade three. The solid line
is the posterior mean and dashed lines are pointwise 90% credible intervals.
Shaded regions correspond to achievement thresholds: Unsatisfactory (U),
Partially Proficient (PP), Proficient (P).

In particular, we analyze a subset of public use microdata from the U.S.
Census, originally compiled by Angrist et al. (2006) to estimate returns to
education across the income distribution. For the 1980, 1990 and 2000 samples
they extracted all U.S. born white and black men aged 40-49 with positive
annual earnings and positive hours worked in the year prior to the Census.
Records with imputed values were excluded, and wages were adjusted to 1989
dollars. Angrist et al. (2006) and its supplementary materials contain details
of how the data were obtained and cleaned. The response is log monthly wages,
and the covariates include years of education, experience (defined as $age-
education-12$) and race (white or black). The study objective was to estimate
returns to education as a function of quantile index. Taking $Q$ as the
quantile function of wages, the return to education at quantile index $s$ is

$100\times\frac{Q(s\mid\mathbf{x}_{2})-Q(s\mid\mathbf{x}_{1})}{Q(s\mid\mathbf{x}_{1})},$
(24)

the predicted percentage change in monthly wage from modifying education in
$\mathbf{x}_{1}$ to yield $\mathbf{x}_{2}$. While Angrist et al. (2006) fit a
linear quantile regression to estimate the effect of each additional year of
education, we fit DR-BART-LH to each census sample separately, subsampling
5000 units from each, and study the difference between 12 and 16 years of
education. In general the effect of a single additional year of schooling is
probably heterogenous. For example, the difference in earnings between 11 and
12 years of schooling should be larger than the difference between 10 and 11,
because 12 years of education in the U.S. typically indicates that the
respondent completed high school. BART readily accommodates non-smooth
features in the regression function, while usual econometric analyses assume
that regression functions are linear or quadratic in education (or
experience).

Figure 11 shows the return to education for a 45 year old white male,
comparing 12 to 16 years of schooling. In 1990, returns were highest at low
and high quantiles, whereas in 2000 the returns are actually increasing as a
function of the quantile index. Other covariate vectors show similar patterns,
although the exact estimates vary (and are somewhat unstable for black men due
to small sample size). Angrist et al. (2006) found a similar pattern in the
returns at the population level and verified it using additional data from the
Current Population Survey, so this seems to be a robust finding.

Figure 11: Returns to education (16 years versus 12 years) for a 45 year old
white man.

##  7 Conclusion

In this paper, we introduced a new nonparametric model, DR-BART, extending
Bayesian Additive Regression Trees to the novel and challenging setting of
density regression. This model has the appeal of being flexible yet easy to
specify and understand, with few prior parameters to set. This distinguishes
DR-BART from other nonparametric Bayesian methods for density regression,
which often include collections of infinite dimensional regression parameters.
Inference via MCMC is fast, with acceptable mixing in our examples obtained in
a matter of minutes for datasets with thousands of observations. We showed DR-
BART to empirically outperform a variety of other density regression methods
in its ability to point estimate conditional densities and to express
uncertainty about these estimates. Lastly, we extended previous work on
posterior concentration results for BART and for density estimation via a
latent variable model to prove concentration rates for DR-BART.

Having introduced the latent variable as a modeling device a natural question
is whether it might represent some real structure in the scientific problem,
like measurement error in $\mathbf{x}$, some combination of omitted variables,
or a latent construct like ability or motivation. Any of these seem plausible
in the applications presented here. In the educational testing application we
treated the latent variables as independent. This is a useful “saturated”
model for the joint distribution of test scores, but models for dependence
across these variables or lower-dimensional representations based on shared
latent variables may be more scientifically meaningful or efficient. Further
modeling of latent variables within BART is a promising area for future work,
even beyond the setting of density regression.

## References

  * Wittman [2009] D. Wittman.  What Lies Beneath: Using $p(z)$ To Reduce Systematic Photometric Redshift Errors.  _The Astrophysical Journal_ , 700(2):L174–L177, jul 2009.
  * Karabatsos and Walker [2011] George Karabatsos and Stephen G. Walker.  Bayesian unimodal density regression for causal inference.  Technical report, Society for Research on Educational Effectiveness, May 2011.
  * Gneiting and Katzfuss [2014] Tilmann Gneiting and Matthias Katzfuss.  Probabilistic forecasting.  _Annual Review of Statistics and Its Application_ , 1:125–151, 2014.
  * Black et al. [2007] Dan A. Black, Jose Galdo, and Jeffrey A. Smith.  Evaluating the worker profiling and reemployment services system using a regression discontinuity approach.  _The American Economic Review_ , 97(2):104–107, 2007.  ISSN 00028282.
  * Kadane and Tokdar [2012] Joseph B. Kadane and Surya T. Tokdar.  Simultaneous Linear Quantile Regression: A Semiparametric Bayesian Approach.  _Bayesian Analysis_ , 7(1):51 – 72, 2012.  doi:10.1214/12-BA702.
  * Sangnier et al. [2016] Maxime Sangnier, Olivier Fercoq, and d’Alché Buc Florence.  Joint Quantile Regression in Vector-Valued RKHSs.  In _Advances in Neural Information Processing Systems_ , 2016.
  * Chernozhukov et al. [2010] Victor Chernozhukov, Iván Fernández-Val, and Alfred Galichon.  Quantile and probability curves without crossing.  _Econometrica_ , 78(3):1093–1125, 2010.
  * Gerfin [1994] Michael Gerfin.  Income distribution, income inequality and life cycle effects - a nonparametric analysis for Switzerland.  _Swiss Journal of Economics and Statistics_ , 130:509–522, 1994.
  * Daly and Valletta [2006] Mary C. Daly and Robert G. Valletta.  Inequality and poverty in United States: The effects of rising dispersion of men’s earnings and changing family behaviour.  _Economica_ , 73:75–98, 2006.
  * Pati et al. [2011] Debdeep Pati, Anirban Bhattacharya, and David Dunson.  Posterior Convergence Rates in Non-linear Latent Variable Models, 2011.
  * Kundu and Dunson [2014] S. Kundu and D. B. Dunson.  Latent factor models for density estimation.  _Biometrika_ , 101(3):641–654, June 2014.  ISSN 0006-3444.  doi:10.1093/biomet/asu019.
  * Zhou et al. [2017] Shuang Zhou, Debdeep Pati, Anirban Bhattacharya, and David Dunson.  Adaptive posterior convergence rates in non-linear latent variable models, 2017.
  * Linero and Yang [2018] A. R. Linero and Y. Yang.  Bayesian regression tree ensembles that adapt to smoothness and sparsity.  _Journal of the Royal Statistical Society: Series B_ , 80(5), 2018.
  * Betebenner et al. [2011] Damian W Betebenner, Adam VanIwaarden, Ben Domingue, Yi Shang, Jonathan Weeks, John Stewart, Jinnie Choi, Xin Wei, Hi Shin Shim, Xiaoyuan Tan, et al.  SGP: An R package for the calculation and visualization of student growth percentiles & percentile growth trajectories.  2011\.
  * Angrist et al. [2006] Joshua Angrist, Victor Chernozhukov, and Iván Fernández-Val.  Quantile regression under misspecification, with an application to the US wage structure.  _Econometrica_ , 74(2):539–563, 2006.
  * Jacobs et al. [1991] Robert A Jacobs, Michael I Jordan, Steven J Nowlan, and Geoffrey E Hinton.  Adaptive mixtures of local experts.  _Neural Computation_ , 3(1):79–87, 1991.
  * Jordan and Jacobs [1994] Michael I Jordan and Robert A Jacobs.  Hierarchical mixtures of experts and the em algorithm.  _Neural computation_ , 6(2):181–214, 1994.
  * Geweke and Keane [2007] John Geweke and Michael Keane.  Smoothly mixing regressions.  _Journal of Econometrics_ , 138(1):252–290, 2007\.
  * Villani et al. [2009] Mattias Villani, Robert Kohn, and Paolo Giordani.  Regression density estimation using smooth adaptive Gaussian mixtures.  _Journal of Econometrics_ , 153(2):155–173, 2009\.
  * Norets et al. [2010] Andriy Norets et al.  Approximation of conditional densities by smooth mixtures of regressions.  _The Annals of Statistics_ , 38(3):1733–1766, 2010.
  * Pati et al. [2013] Debdeep Pati, David B Dunson, and Surya T Tokdar.  Posterior consistency in conditional distribution estimation.  _Journal of Multivariate Analysis_ , 116:456–472, 2013\.
  * Norets and Pelenis [2014] Andriy Norets and Justinas Pelenis.  Posterior consistency in conditional density estimation by covariate dependent mixtures.  _Econometric Theory_ , 30(03):606–646, 2014.
  * Norets and Pati [2014] Andriy Norets and Debdeep Pati.  Adaptive Bayesian estimation of conditional densities, 2014.
  * Muller et al. [1996] P. Muller, A. Erkanli, and M. West.  Bayesian curve fitting using multivariate normal mixtures.  _Biometrika_ , 83(1):67–79, March 1996.  ISSN 0006-3444.  doi:10.1093/biomet/83.1.67.
  * Park and Dunson [2010] Ju-Hyun Park and David B Dunson.  Bayesian generalized product partition model.  _Statistica Sinica_ , 20(20):1203–1226, 2010\.
  * West and Escobar [1993] Mike West and Michael D Escobar.  _Hierarchical priors and mixture models, with applications in regression and density estimation_.  Duke University, 1993.
  * Shahbaba and Neal [2009] Babak Shahbaba and Radford Neal.  Nonlinear Models Using Dirichlet Process Mixtures.  _The Journal of Machine Learning Research_ , 10:1829–1850, December 2009.  ISSN 1532-4435.
  * Taddy and Kottas [2010] Matthew A. Taddy and Athanasios Kottas.  A Bayesian Nonparametric Approach to Inference for Quantile Regression.  _Journal of Business & Economic Statistics_, 28(3):357–369, July 2010.  ISSN 0735-0015.  doi:10.1198/jbes.2009.07331.
  * Molitor et al. [2010] John Molitor, Michail Papathomas, Michael Jerrett, and Sylvia Richardson.  Bayesian profile regression with an application to the National survey of children’s health.  _Biostatistics_ , 11(3):484–98, July 2010.  ISSN 1468-4357.  doi:10.1093/biostatistics/kxq013.
  * Wade et al. [2011] Sara Wade, Silvia Mongelluzzo, and Sonia Petrone.  An Enriched Conjugate Prior for Bayesian Nonparametric Inference.  _Bayesian Analysis_ , 6(3):359–385, 2011.
  * Dunson and Bhattacharya [2011] David B Dunson and Abhishek Bhattacharya.  Nonparametric Bayes Regression and Classification Through Mixtures of Product Kernels.  In Jose M. Bernardo, M. J. Bayarri, James O. Berger, A. P. Dawid, David Heckerman, Adrian F. M. Smith, and Mike West, editors, _Proceedings of 9th Valencia International Conference on Bayesian Statistics_. Oxford University Press, 2011.
  * Hannah et al. [2011] Lauren A. Hannah, David M. Blei, and Warren B. Powell.  Dirichlet Process Mixtures of Generalized Linear Models.  _The Journal of Machine Learning Research_ , pages 1923–1953, July 2011.  ISSN 1532-4435.
  * Wade et al. [2014] Sara Wade, DB Dunson, Sonia Petrone, and Lorenzo Trippa.  Improving Prediction from Dirichlet Process Mixtures via Enrichment.  _Journal of Machine Learning Research_ , 15:1041–1071, 2014\.
  * Walker and Karabatsos [2013] Stephen G. Walker and George Karabatsos.  Revisiting Bayesian curve fitting using multivariate normal mixtures.  In Paul Damien, Petros Dellaportas, Nicholas G. Polson, and David A. Stephens, editors, _Bayesian Theory and Applications_. Oxford University Press, 01 edition, 3 2013.  ISBN 9780199695607.
  * MacEachern [1999] Steven N MacEachern.  Dependent nonparametric processes.  In _ASA proceedings of the section on Bayesian statistical science_ , pages 50–55, 1999.
  * MacEachern [2000] Steven N MacEachern.  Dependent Dirichlet processes.  _Unpublished manuscript, Department of Statistics, The Ohio State University_ , 2000.
  * De Iorio et al. [2004] Maria De Iorio, Peter Müller, Gary L Rosner, and Steven N MacEachern.  An ANOVA model for dependent random measures.  _Journal of the American Statistical Association_ , 99(465):205–215, 2004.
  * Griffin and Steel [2006] Jim E Griffin and MF J Steel.  Order-based dependent Dirichlet processes.  _Journal of the American Statistical Association_ , 101(473):179–194, 2006.
  * Dunson and Peddada [2008] David B Dunson and Shyamal D Peddada.  Bayesian nonparametric inference on stochastic ordering.  _Biometrika_ , 95(4):859–874, 2008.
  * De Iorio et al. [2009] Maria De Iorio, Wesley O Johnson, Peter Müller, and Gary L Rosner.  Bayesian nonparametric nonproportional hazards survival modeling.  _Biometrics_ , 65(3):762–771, 2009.
  * Wang and Dunson [2011] Lianming Wang and David B Dunson.  Bayesian isotonic density regression.  _Biometrika_ , 98(3):537–551, 2011.
  * Barrientos et al. [2012] Andrés F. Barrientos, Alejandro Jara, and Fernando A. Quintana.  On the support of MacEachern’s dependent Dirichlet processes and extensions.  _Bayesian Analysis_ , 7(2):277–310, 06 2012\.  doi:10.1214/12-BA709.
  * Chipman et al. [2010] Hugh A. Chipman, Edward I. George, and Robert E. McCulloch.  BART: Bayesian additive regression trees.  _The Annals of Applied Statistics_ , 4(1):266–298, March 2010.  ISSN 1941-7330.
  * Jeong and Rockova [2020] Seonghyun Jeong and Veronika Rockova.  The art of BART: On flexibility of Bayesian forests, 2020.
  * Murray [2021] Jared S Murray.  Log-linear bayesian additive regression trees for multinomial logistic and count regression models.  _Journal of the American Statistical Association_ , pages 1–14, 2021\.
  * Walker [2007] Stephen G. Walker.  Sampling the Dirichlet Mixture Model with Slices.  _Communications in Statistics - Simulation and Computation_ , 36(1):45–54, January 2007.  ISSN 0361-0918.  doi:10.1080/03610910601096262.
  * Kalli et al. [2011] Maria Kalli, JE Griffin, and SG Walker.  Slice sampling mixture models.  _Statistics and Computing_ , 21(1):93–105, September 2011.  ISSN 0960-3174.  doi:10.1007/s11222-009-9150-y.
  * Chipman et al. [1998] Hugh A Chipman, Edward I George, and Robert E McCulloch.  Bayesian CART model search.  _Journal of the American Statistical Association_ , 93(443):935–948, 1998.
  * Denison et al. [1998] David GT Denison, Bani K Mallick, and Adrian FM Smith.  A Bayesian CART algorithm.  _Biometrika_ , 85(2):363–377, 1998.
  * Diebolt and Robert [1994] Jean Diebolt and Christian P Robert.  Estimation of finite mixture distributions through Bayesian sampling.  _Journal of the Royal Statistical Society. Series B_ , pages 363–375, 1994.
  * Neal [2003] Radford M Neal.  Slice sampling.  _Annals of Statistics_ , pages 705–741, 2003.
  * Chung and Dunson [2012] Yeonseung Chung and David Dunson.  Nonparametric bayes conditional distribution modeling with variable selection.  _Journal of the American Statistical Association_ , 2012.
  * Li et al. [2020] Yinpu Li, Antonio R. Linero, and Jared S. Murray.  Adaptive conditional distribution estimation with Bayesian decision tree ensembles, 2020.
  * Plummer et al. [2021] Sean Plummer, Shuang Zhou, Anirban Bhattacharya, David Dunson, and Debdeep Pati.  Statistical guarantees for transformation based models with applications to implicit variational inference.  In _The 24th Intl. Conference on Artificial Intelligence and Statistics, AISTATS 2021_ , 2021.
  * Ghosal et al. [2000] Shen Ghosal, Jayanta K. Ghosh, and A.W. van der Vaart.  Convergence Rates of Posterior Distributions.  _Annals of Statistics_ , 28(2):500––531, 2000\.
  * Starling et al. [2020] Jennifer Starling, Jared Murray, Carlos Carvalho, Radek Bukowski, and James Scott.  BART with targeted smoothing: An analysis of patient-specific stillbirth risk.  _Annals of Applied Statistics_ , 14(1), 2020.
  * Jara et al. [2011] Alejandro Jara, Timothy Hanson, Fernando A. Quintana, Peter Múller, and Gary L. Rosner.  DPpackage: Bayesian Semi- and Nonparametric Modeling in R.  _Journal of Statistical Software_ , 40(5), 2011.
  * Castellano and Ho [2013] Katherine Elizabeth Castellano and Andrew Dean Ho.  Contrasting OLS and quantile regression approaches to student “growth” percentiles.  _Journal of Educational and Behavioral Statistics_ , 38(2):190–215, 2013.
  * Betebenner [2009] D Betebenner.  Norm-and criterion-referenced student growth.  _Educational Measurement: Issues and Practice_ , 28(4):42–51, 2009.
  * Betebenner [2011] Damian W Betebenner.  A technical overview of the student growth percentile methodology: Student growth percentiles and percentile growth projections/trajectories. the national center for the improvement of educational assessment, 2011.
  * Pratola et al. [2014] Matthew T. Pratola, Hugh A. Chipman, James R. Gattiker, David M. Higdon, Robert McCulloch, and William N. Rust.  Parallel Bayesian additive regression trees.  _Journal of Computational and Graphical Statistics_ , 23(3):830–852, 2014.  doi:10.1080/10618600.2013.841584.
  * Bleich and Kapelner [2014] Justin Bleich and Adam Kapelner.  Bayesian Additive Regression Trees With Parametric Models of Heteroskedasticity.  _arXiv preprint arXiv:1402.5397_ , February 2014.
  * Ghosal and van der Vaart [2007] Shen Ghosal and Aad van der Vaart.  Convergence Rates of Posterior Distributions for Noniid Observations.  _Annals of Statistics_ , 35:192––223, 2007.
  * Stan Development Team [2020] Stan Development Team.  RStan: the R interface to Stan, 2020.  R package version 2.21.2.

##  8 Supplement

###  8.1 BART

In this section, we review the basic details of Bayesian Additive Regression
Trees (BART). We refer the reader to Chipman et al. [2010] for a full
exposition.

####  8.1.1 Model and Priors

Chipman et al. consider the regression model

$y_{i}=f(\mathbf{x}_{i})+\epsilon_{i},\epsilon_{i}\stackrel{{\scriptstyle
iid}}{{\sim}}N(0,\sigma^{2})$

and place a BART prior on $f$. Such a prior represents the function as the sum
of $m$ many piecewise constant regression trees. Each tree $T_{h},\;1\leq
h\leq m$ consists of a set of interior decision nodes (where decisions are
generally of the form $x_{j}<c$ for some value $c$) and a set of $b_{h}$
terminal nodes. The terminal nodes have associated parameters
$M_{h}=(\mu_{h1},\mu_{h2},\dots\mu_{hb_{h}})^{\prime}$. For each tree there is
a partition of the covariate space
$\\{\mathcal{A}_{h1},\dots,\mathcal{A}_{hb_{h}}\\}$ with each element of the
partition corresponding to a terminal node. A tree and its associated
parameters define step functions:

$g(\mathbf{x},T_{h},M_{h})=\mu_{hb}\text{ if
}\mathbf{x}\in\mathcal{A}_{hb}\text{ (for $1\leq b\leq b_{h}$)}.$ (25)

These functions are additively combined to obtain $f$:

$f(\mathbf{x})=\sum_{h=1}^{m}g(\mathbf{x},T_{h},M_{h}).$ (26)

In the spirit of boosting, each term in the sum is constrained by a strong
prior to be a “weak learner"; that is, the prior on $(T_{h},M_{h})$ strongly
favors small trees and leaf parameters that are near zero. Each tree
independently follows the prior described by Chipman et al. [1998], where the
probability that a node at depth $d$ splits (is not terminal) is given by

$\alpha(1+d)^{-\beta},\;\;\alpha\in(0,1),\;\beta\in[0,\infty).$ (27)

A variable to split on, and a cutpoint to split at, are then selected
uniformly at random from the available splitting rules. We follow CGM
throughout by taking $\alpha=0.95$ and $\beta=2$. Traditionally the prior on
cutpoints for the $j^{th}$ variable is a discrete uniform distribution over a
uniformly spaced grid or some collection of quantiles of $\\{x_{ij}:1\leq
i\leq n\\}$. Most implementations of tree-based models also require that the
terminal nodes not be empty, or contain at least 5 observations.

To set the priors on $M$, CGM suggest scaling the data to lie in $\pm\,0.5$
and assign the leaf parameters independent priors:

$\mu_{hb}\sim N(0,\sigma^{2}_{\mu})\;\;\textrm{ where
}\sigma_{\mu}=1/(2k\sqrt{m}).$ (28)

CGM recommend $1\leq k\leq 3$, with $k=2$ as a reasonable default choice. This
prior shrinks $g_{l}(x)$ strongly toward zero, while ensuring that the induced
prior for $f(x)$ is centered at zero and puts approximately 95% of the prior
mass within $\pm 0.5$. Larger values of $k$ imply increasing degrees of
shrinkage. Performance is fairly insensitive to the number of trees $m$, as
long as it is large enough. In practice, $m=200$ is a common choice.

####  8.1.2 Posterior Sampling

Chipman et al. [2010] provide full details of the MCMC algorithm used to fit
BART. A key ingredient of the sampler is performing blocked updates for
$(T_{h}M_{h}\mid\\{(T_{l},M_{l}\\}_{l\neq h},-)$. This “backfitting” step
utilizes the fact that the full conditional for $(T_{h},M_{h})$ depends on
$\\{(T_{l},M_{l}\\}_{l\neq h}$ and $y$ only through the residuals

$R_{hi}=\left(y_{i}-\sum_{l\neq h}^{m}g(\mathbf{x}_{i},T_{l},M_{l})\right)\sim
N(g(\mathbf{x}_{i},T_{h},M_{h}),\sigma^{2}),$ (29)

These residuals follow the single tree model studied by Chipman et al. [1998]
with parameters $(T_{h},M_{h},\sigma^{2})$, so the Chipman et al. [1998]
Metroplis-Hastings update can be embedded within the BART Gibbs sampler to
sample $(T_{h},M_{h})$ jointly from their full conditional. Since $M_{h}$ has
a normal prior, the integrated likelihood $p(R_{h}\mid T_{h},\sigma^{2})$ is
available in closed form. Therefore $(T_{h},M_{h})$ can be updated in a block
by sampling $T_{h}$ marginally over $M_{h}$ with a Metropolis step and then
sampling $M_{h}$ from its full conditional given the new $T_{h}$. Complete
details including proposal distributions are given in Chipman et al. [1998]
(the results in this paper use a smaller set of propsal distributions,
outlined in Pratola et al. [2014]). This blocked update enhances mixing and
obviates the need for transdimensional MCMC algorithms due to the changing
dimensionality of $M_{h}$.

###  8.2 Details on Heteroscedastic Regression with BART priors

####  8.2.1 Model and Priors

Here, we give an overview of the Murray [2021] extension of the BART prior to
model variance functions. Begin with the loglinear model

$v(x)=\sum_{h=1}^{m_{v}}g(\mathbf{x},T^{(v)}_{h},M^{(v)}_{h})$ (30)

where $\\{(T^{(v)}_{h}$, $M^{(v)}_{h})\\}$ are trees and parameters for the
variance functions. To collect posterior samples, a variant of the backfitting
MCMC algorithm is possible. Suppose we are updating
$(T^{(v)}_{h},M^{(v)}_{h})$. The first step is constructing the scaled
residuals

$R^{(v)}_{hi}=\frac{(Y_{i}-f(x_{i}))}{\exp[\sum_{l\neq
h}^{m}g(\mathbf{x},T^{(v)}_{l},M^{(v)}_{l})]}.$ (31)

Now $R^{(v)}_{hi}\sim
N(0,\sigma^{2}\exp[g(\mathbf{x},T^{(v)}_{h},M^{(v)}_{h})])$ so we have

$p(R^{(v)}_{h}\mid
T^{(v)}_{h},M^{(v)}_{h},\sigma^{2})=\prod^{b^{(v)}_{h}}_{l=1}\prod_{i:\mathbf{x}_{i}\in\mathcal{A}_{l}}\frac{1}{\sigma_{0}\exp[\mu^{(v)}_{hl}/2]}\phi\left(\frac{R^{(v)}_{ih}}{\sigma_{0}^{2}\exp[\mu^{(v)}_{hl}]}\right)$
(32)

As a normal prior for $\mu^{(v)}_{hb}$ is no longer conditionally conjugate,
Murray introduces another prior which is symmetric about $0$ on the log scale
and also admits closed form representations for the integrated likelihood and
full conditional distribution. This prior is specified as an equal probability
mixture of gamma and inverse gamma distributions with the same parameters. Let
$\tau^{(v)}_{hb}:=\exp(\mu_{hb}^{(v)})$. Then, the prior is:

$p\left(\tau^{(v)}_{hb}\right)=\frac{1}{2}\frac{b^{a}}{\Gamma(a)}\left(\tau^{(v)}_{hb}\right)^{-a-1}\exp\left(-b/\tau^{(v)}_{hb}\right)+\frac{1}{2}\frac{b^{a}}{\Gamma(a)}\left(\tau^{(v)}_{hb}\right)^{a-1}\exp\left(-b\tau^{(v)}_{hb}\right)$
(33)

Throughout this paper we take $a=b=a_{0}m_{v}$. To motivate this choice, note
that if $\tau^{(v)}_{hb}$ has density (33) with $a=b=a_{0}m_{v}$ then
$\log(\tau^{(v)}_{hb})\,{\buildrel d\over{=}}\,(-1)^{Z_{h}}\log(W_{h})$ where
$Z_{h}\overset{iid}{\sim}Bern(0.5)$ and
$W_{h}\overset{iid}{\sim}Ga(a_{0}m_{v},a_{0}m_{v})$. The first two moments of
$\log(\tau^{(v)}_{hb}))$ are 0 and $(a_{0}m_{v})^{-1}$ and for large
$a_{0}m_{v}$ its distribution is nearly normal since the Gamma distributions
are so concentrated near 1. Thus, marginally, $v(x)$ will be approximately
distributed $N(0,1/a_{0})$ and can be scaled to match the _a priori_ expected
range of the log-variance process. Even for relatively small $a_{0}m_{v}$ this
approximation is excellent; see Figure 12. The marginal prior for $\log(v(x))$
is nearly normal for $a_{0}m_{v}=5$ and is indistinguishable from the normal
density when $a_{0}m_{v}$ is greater than 25. In practice $m_{v}$ will
generally be at least 100. To set $a_{0}$, note that
$a_{0}=[\log(\sqrt{d})]^{-2}$ implies that
$\Pr(\exp[v(\mathbf{x})]\in(\sigma^{2}_{0}/d,d\sigma_{0}^{2}))\approx 0.95$.
For example, taking $a_{0}=1.5$ makes $d$ about $5$. In practice, values of
$d$ around 2-4 tend to work well; much larger and the risk of overfitting and
degenerate bandwidths increases.

Figure 12: Marginal prior on $v(x)$ for $a_{0}=1$ and a range of $m_{v}$. The
blue dashed line is the standard normal density.

####  8.2.2 Posterior Sampling

Let
$n_{hb}=\sum_{i=1}^{n}{\mathbf{1}}{\left(\mathbf{x}_{i}\in\mathcal{A}^{(v)}_{hb}\right)}$,
the number of observations in leaf $b$ of $T^{(v)}_{h}$, and let
$r^{2}_{hb}=\sum_{i=1}^{n}[R^{(v)}_{hi}]^{2}{\mathbf{1}}{\left(\mathbf{x}_{i}\in\mathcal{A}^{(v)}_{hb}\right)}$
be the sum of squared residuals in each leaf. Under (33) the full conditional
distribution for $\tau^{(v)}_{hb}$ is another mixture distribution:

$\begin{split}p(\tau^{(v)}_{hb}\mid
R^{(v)}_{h},\sigma^{2})\propto\frac{\Gamma(n_{hb}/2+a)}{(b+r^{2}_{hb}/2)^{n_{hb}/2+a}}p_{IG}(\tau^{(v)}_{hb},n_{hb}/2+a,b+r_{hb}^{2}/2)\\\
+\frac{2K_{(n_{hb}/2-a)}(\sqrt{2br_{hb}^{2}})}{(2b/r_{hb}^{2})^{(a-n_{hb}/2)/2}}p_{GIG}(a-n_{hb}/2,b,r_{hb}^{2}/2)\end{split}$
(34)

where $K_{p}(x)$ is the modified Bessel function of the second kind, $p_{IG}$
is the pdf of an inverse gamma random variable with rate $b$, and $p_{GIG}$ is
the pdf of a generalized inverse Gaussian random variable:

$p_{GIG}(\tau;\lambda,\psi,\chi)=\frac{(\psi/\chi)^{(\lambda/2)}}{2K_{\lambda}(\sqrt{\psi\chi})}\tau^{\lambda-1}\exp\left[{-\frac{1}{2}}\left(\psi\tau+\chi/\tau\right)\right]$
(35)

The integrated likelihood is

$P(R^{(v)}_{h}\mid
T^{(v)}_{h},\sigma^{2})=\prod^{b^{(v)}_{h}}_{l=1}\frac{1}{2(2\pi)^{n_{hl}/2}}\frac{b^{a}}{\Gamma(a)}\left(\frac{\Gamma(n_{hl}/2+a)}{(b+r^{2}_{hl}/2)^{n_{hl}/2+a}}+\frac{2K_{(n_{hl}/2-a)}(\sqrt{2br^{2}_{hl}})}{(2b/r_{hl}^{2})^{(a-n_{hl}/2)/2}}\right).$
(36)

Block updating $(T^{(v)}_{h},M^{(v)}_{h})$ proceeds by substituting the
integrated likelihood (36) into the Metropolis step for $T^{(v)}_{h}$ and
drawing $M^{(v)}_{h}$ by sampling each $\tau_{hb_{h}}^{(v)}$ from (36). Given
$v$, updates for the rest of the parameters are straightforward although the
integrated likelihood and full conditionals also have to be adjusted in the
update for $(T_{h},M_{h})$ to account for the heteroscedastic residuals. See
Bleich and Kapelner [2014] for details.

###  8.3 Proofs of Theorems

#### Proof of Theorem 1

Let $N_{d}$ be the number of nodes at depth $d$ of a tree sampled from this
prior. There are $2^{d-1}$ pairs of possible nodes at depth $d$. Label these
pairs by $1\leq i\leq 2^{d-1}$ and define $U_{i}^{(d)}$ as 1 if the $i^{th}$
pair of nodes exists and 0 otherwise, so that
$N_{d}=\sum_{i=1}^{2^{d-1}}2U_{i}^{(d)}$. Now $U_{i}^{(d)}=1$ if and only if
all its parents are splitting nodes, so
$\Pr(U_{i}^{(d)}=1)=\mathrm{E}(U_{i}^{(d)})=\prod_{s=0}^{d-1}\alpha(1+s)^{-\beta}=\alpha^{d}(d!)^{-\beta}$.
Then
$E(N_{d})=\sum_{i=1}^{2^{d-1}}2\alpha^{d}(d!)^{-\beta}=(2\alpha)^{d}(d!)^{-\beta}$.
Since $N_{d}$ is nonnegative, $\Pr(N_{d}\geq 1)\leq\mathrm{E}(N_{d})$. Since
$\sum_{d=0}^{\infty}\mathrm{E}(N_{d})=\sum_{d=0}^{\infty}(2\alpha)^{d}(d!)^{-\beta}<\infty$,
$\Pr(N_{d}\geq 1\ i.o.)=0$ so the tree is almost surely finite.

#### Proof of Theorem 2

As discussed in Section 4, we must find constants $C_{1}$, …, $C_{4}$ and a
sieve $\mathcal{G}_{n}$ that allow us to verify conditions G1, G2, and G3. We
first prove the conditions using Lemmas 1 – 3; in the next subsection we prove
a series of Lemmas allowing us to establish Lemma 3.

###### Proof of G1.

Let $\delta_{n}=(\log n/n)^{\beta/(2\beta+d_{0})}$ and let
$\widetilde{\epsilon}_{n}=(\log
n/n)^{\frac{\alpha}{\alpha+1}\times\frac{\beta}{2\beta+d_{0}}}$.

Applying Lemma 3 with $\epsilon=\widetilde{\epsilon}_{n}$, for sufficiently
large $n$ we have

$\displaystyle\begin{split}-\log\Pi\\{\operatorname{KL}_{p_{0}}(C_{\operatorname{KL}}\widetilde{\epsilon}_{n})\\}&\leq-\log\Pi(S=S_{0})-\log\Pi(\|f-f_{0}\|_{\infty}\leq\delta_{n}\mid
S=S_{0})\\\ &-\log\Pi\\{\sigma(\mathbf{x},u)\equiv\sigma\text{ is
constant},\sigma\in[\widetilde{\epsilon}_{n}^{1/\alpha},2\widetilde{\epsilon}_{n}^{1/\alpha}]\mid
S=S_{0}\\}.\end{split}$ (37)

As argued by Jeong and Rockova [2020], P1 implies $-\log\Pi(S=S_{0})\lesssim
d_{0}\log(p+1)\lesssim n\epsilon_{n}^{2}$. Next, by Lemma 1 we have
$-\log\Pi(\|f-f_{0}\|_{\infty}\leq\delta_{n}\mid S=S_{0})\lesssim
n\delta_{n}^{2}\lesssim n\epsilon_{n}^{2}$. Finally, note that with
$\sigma_{n}=\widetilde{\epsilon}_{n}^{1/\alpha}$, by P4 and standard
properties of the inverse gamma distribution we can bound the final term of
(37) by

$\displaystyle-\log\Pi\\{e^{\xi}\in[1.2\sigma_{n},1.8\sigma_{n}],\sup_{hl}|\mu_{hl}^{(v)}|\leq
0.1/m_{v},\text{all trees have depth 0}\mid S=S_{0}\\}$
$\displaystyle\quad\lesssim-\log\Pi\\{e^{\xi}\in[1.2\sigma_{n},1.8\sigma_{n}]\\}\lesssim\sigma_{n}^{-1}\lesssim
n\epsilon_{n}^{2}.$

Combining these facts, we get
$-\log\Pi\\{\operatorname{KL}_{p_{0}}(C_{\operatorname{KL}}\widetilde{\epsilon}_{n})\\}\lesssim
n\epsilon_{n}^{2}$, and because $\widetilde{\epsilon}_{n}\leq\epsilon_{n}$
this implies that G1 holds for some choice of $C_{1}$ and $C_{2}$. ∎

###### Proof of G2 and G3.

Let $\kappa$ be a large constant to be determined later and define
$\mathcal{F}_{f}$ as in Lemma 2 with the choices $d=\kappa
n\epsilon_{n}^{2}/\log(p+1)$, $K=\lfloor\kappa n\epsilon_{n}^{2}/\log
n\rfloor$, $U^{2}=\kappa n\epsilon_{n}^{2}$, and
$\epsilon=\epsilon_{n}/(n\epsilon_{n}^{2})$. Let $\tilde{v}$ denote the tree-
based contribution to the variance function, without the intercept:
$\tilde{v}:=v(\mathbf{x},u)-\xi$. Then, we similarly define
$\mathcal{F}_{\tilde{v}}$ but with the choice $\epsilon=\epsilon_{n}^{2}$ and
$U=V$. Lastly, define $\ell=1/(\kappa n\epsilon_{n}^{2})$ and $u=e^{\kappa
n\epsilon_{n}^{2}}$. We take the set $\mathcal{G}$ to be given by

$\displaystyle\mathcal{G}=\\{p_{f,\sigma}:f\in\mathcal{F}_{f},\tilde{v}\in\mathcal{F}_{\tilde{v}},\xi\in[\log\ell,\log
u]\\}.$

Let $\widehat{F}_{f}$ and $\widehat{F}_{\tilde{v}}$ denote
$\epsilon_{n}/(n\epsilon_{n}^{2})$ and $\epsilon_{n}^{2}$ nets of
$\mathcal{F}_{f}$ and $\mathcal{F}_{\tilde{v}}$ respectively and let
$\widehat{R}$ denote an $\epsilon_{n}^{2}$-net of $[\log\ell,\log u]$. Given a
$p_{f_{1},\sigma_{1}}\in\mathcal{G}$ we can find $f_{2}\in\widehat{F}_{f}$,
$\tilde{v}_{2}\in\widehat{F}_{\tilde{v}}$, and $\xi_{2}\in\widehat{R}$ such
that, by Lemma 3:

$\displaystyle h(p_{f_{1},\sigma_{1}},p_{f_{2},\sigma_{2}})\leq
C_{h}\left\\{\sqrt{2\epsilon_{n}^{2}}+e^{m_{v}V}\frac{\epsilon_{n}(\kappa
n\epsilon_{n}^{2})}{n\epsilon_{n}^{2}}\right\\}\leq C_{\kappa}\epsilon_{n},$

for a global constant $C_{h}$ and a constant $C_{\kappa}$ depending only on
the constant $\kappa$ and the prior. Hence

$\displaystyle\log N(\mathcal{G},C_{\kappa}\epsilon_{n},h)\leq\log
N(\mathcal{F}_{f},\epsilon_{n}/(n\epsilon_{n}^{2}),\|\cdot\|_{\infty})+\log
N(\mathcal{F}_{\tilde{v}},\epsilon_{n}^{2},\|\cdot\|_{\infty})+\log\frac{\log
u+\log\ell^{-1}}{\epsilon_{n}^{2}}.$

By Lemma 2, each term is easily verified to be bounded by a constant multiple
of $\kappa n\epsilon_{n}^{2}$ so that $\log
N(\mathcal{G},C_{\kappa}\epsilon_{n},h)\lesssim\kappa n\epsilon_{n}^{2}$.

Next we compute the complementary probability bound. First, by the union
bound:

$\displaystyle\begin{split}\Pi(p\notin\mathcal{G})&\leq\Pi(\xi<\log\ell)+\Pi(\xi>\log
u)+\Pi(D>d)\\\ &\quad+\Pi(f\notin\mathcal{F}_{f}\mid D\leq
d)+\Pi(v\notin\mathcal{F}_{\tilde{v}}\mid D\leq d).\end{split}$ (38)

Using standard properties of the inverse gamma distribution we have

$\displaystyle\log\Pi(\xi<\log\ell)$ $\displaystyle\lesssim-\ell^{-1}=-\kappa
n\epsilon_{n}^{2},$ $\displaystyle\log\Pi(\xi>\log u)$
$\displaystyle\lesssim-\log u=-\kappa n\epsilon_{n}^{2}.$

By P1 we also have $\log\Pi(D>d)\lesssim-d\log(p+1)=-\kappa
n\epsilon_{n}^{2}$. Next, if $D\leq d$ but $f\notin\mathcal{F}_{f}$ then
either (i) $K_{m}>K$ for some $m$ where $K_{m}$ is the number of leaf nodes in
tree $m$ or (ii) $K_{m}\leq K$ for all $m$, but $\sup_{h\ell}|\mu_{h\ell}|>U$.
Hence, by the tail properties of Gaussian random variables,

$\displaystyle\Pi(f\notin\mathcal{F}_{f}\mid D\leq d)$
$\displaystyle\leq\Pi(K_{m}>K\text{ for some
$m$})+\Pi(\sup_{h\ell}|\mu_{h\ell}|>U\mid K_{m}\leq K\text{ for all $m$})$
$\displaystyle\leq m\Pi(K_{1}>K)+Kme^{-U^{2}/(2\sigma^{2}_{\mu})}.$

As noted in the proof of Theorem 2 of Jeong and Rockova [2020] we have
$\log\Pi(K_{1}>K)\lesssim-K\log K\lesssim-\kappa n\epsilon_{n}^{2}$ while
$-U^{2}/(2\sigma^{2}_{\mu})\lesssim-\kappa n\epsilon_{n}^{2}$ by our choice of
$U$. Hence $\log\Pi(f\notin\mathcal{F}_{f}\mid D\leq d)\lesssim-\kappa
n\epsilon_{n}^{2}$ as well, provided that $n$ is sufficiently large. Finally,
the only way for $\tilde{v}\notin\mathcal{F}_{\tilde{v}}$ to occur when $D\leq
d$ is for at least one tree to have more than $K$ leaf nodes, and we have
already seen that the log of this probability is bounded by a constant
multiple of $-\kappa n\epsilon_{n}^{2}$. Putting all of these facts together
and bounding each term of (38) by the maximum, we have

$\displaystyle\Pi(p\notin\mathcal{G})\leq\exp\\{-C_{\mathcal{G}}\kappa
n\epsilon_{n}^{2}\\}$

for some small constant $C_{\mathcal{G}}$ depending only on the prior. By
taking $\kappa$ larger than $(C_{2}+4)/C_{\mathcal{G}}$ we have
$\Pi(p\notin\mathcal{G})\leq\exp\\{-(C_{2}+4)n\epsilon_{n}^{2}\\}$ and $\log
N(\mathcal{G},\bar{\epsilon}_{n},h)\leq C_{5}n\bar{\epsilon}_{n}^{2}$ for some
constant $C_{5}$ and $\bar{\epsilon}_{n}$ a constant multiple of
$\epsilon_{n}$, as desired. ∎

#### Proof of Supporting Lemmas

As before, we use $a\lesssim b$ to mean there is a positive constant $C$ which
can be computed from the prior and $p_{0}$, independent of $n$ and $p$, such
that $a\leq Cb$. Unless otherwise stated, the constant is assumed to be
universal in the sense that if we write $f(x;\xi)\lesssim g(x)$ then we have
$f(x;\xi)\leq Cg(x)$ for all $(x,\xi)$ (unless $\xi$ is part of the prior or a
function of $p_{0}$).

Throughout this section, $\phi_{\sigma}*p_{0}(y\mid\mathbf{x})$ denotes the
convolution $\int\phi_{\sigma}(y-z)\,p_{0}(z\mid\mathbf{x})\ dz$ where
$\phi_{\sigma}(z)=(2\pi\sigma^{2})^{-1/2}e^{-z^{2}/(2\sigma^{2})}$. When
$\sigma(\mathbf{x},u)$ is a constant we note that
$\phi_{\sigma}*p_{0}(y\mid\mathbf{x})=p_{f_{0},\sigma}$.

######  Lemma 4.

Suppose $p_{0}$ satisfies Condition F. Then for $\sigma<1$, we have

$\displaystyle
p_{0}(y\mid\mathbf{x})\lesssim\phi_{\sigma}*p_{0}(y\mid\mathbf{x})$

###### Proof.

By Condition F, $p_{0}(y\mid\mathbf{x})$ is uniformly bounded and bounded away
from 0 on $[0,1]$. First, if $y\notin[0,1]$ then we have
$\phi_{\sigma}*p_{0}(y\mid\mathbf{x})>0=p_{0}(y\mid\mathbf{x})$ so that
$\phi_{\sigma}*p_{0}(y\mid\mathbf{x})>Cp_{0}(y\mid\mathbf{x})$ for any
positive constant $C$; hence we can assume without loss of generality that
$y\in[0,1]$. Let $M=\inf_{y,\mathbf{x}}p_{0}(y\mid\mathbf{x})$ and write

$\displaystyle\phi_{\sigma}*p_{0}(y\mid\mathbf{x})$
$\displaystyle=\int_{0}^{1}\phi_{\sigma}(y-z)\,p_{0}(z\mid\mathbf{x})\ dz$
$\displaystyle\geq M\int_{0}^{1}\phi_{\sigma}(y-z)\ dz$
$\displaystyle=M\int_{y-1}^{y}\phi_{\sigma}(u)\ du.$ $\displaystyle\geq
M\int_{0}^{1}\phi_{\sigma}(u)\ du\geq
M/3\geq\frac{M}{3\sup_{y,\mathbf{x}}p_{0}(y\mid\mathbf{x})}p_{0}(y\mid\mathbf{x}),$

where the second inequality follows from the fact
$U(y;\sigma)=\Phi\\{y/\sigma\\}-\Phi\\{(y-1)/\sigma\\}$ is minimized at $y=1$
for $y\in[0,1]$ and all $\sigma$, and the third inequality follows from the
fact that $U(1;\sigma)>1/3$ for all $\sigma<1$. Hence
$\phi_{\sigma}*p_{0}(y\mid\mathbf{x})\gtrsim p_{0}(y\mid\mathbf{x})$ for
$\sigma<1$. ∎

######  Lemma 5.

If $p_{0}$ satisfies Condition F and $\sigma(\mathbf{x},u)$ is constant with
$\sigma<1$ then

$\displaystyle\log\|p_{0}/p_{f,\sigma}\|_{\infty}\leq
C_{\inf}+\frac{\|f-f_{0}\|^{2}_{\infty}}{\sigma^{2}},$

for some $C_{\inf}$.

###### Proof.

The proof is essentially the same as the proof of Lemma 6.2 of Zhou et al.
[2017]. We write

$\displaystyle p_{f,\sigma}(y\mid\mathbf{x})$
$\displaystyle=\frac{1}{\sqrt{2\pi}\sigma}\int_{0}^{1}\exp\left\\{-\frac{(y-f(\mathbf{x},u))^{2}}{2\sigma^{2}}\right\\}$
$\displaystyle\geq\sqrt{2}\exp\\{-\|f-f_{0}\|_{\infty}^{2}/\sigma^{2}\\}\int_{0}^{1}\frac{1}{\sqrt{2\pi\sigma^{2}/2}}\exp\left\\{-\frac{(y-f_{0}(\mathbf{x},u))^{2}}{\sigma^{2}}\right\\}$
$\displaystyle=\sqrt{2}\exp\\{-\|f-f_{0}\|_{\infty}^{2}/\sigma^{2}\\}\phi_{\sigma/\sqrt{2}}*p_{0}(y\mid\mathbf{x}).$

Therefore, using the fact that
$\phi_{\sigma/\sqrt{2}}*p_{0}(y\mid\mathbf{x})\gtrsim p_{0}(y\mid\mathbf{x})$
we have

$\displaystyle\log\frac{p_{0}(y\mid\mathbf{x})}{p_{f,\sigma}(y\mid\mathbf{x})}\leq
C_{\inf}+\frac{\|f-f_{0}\|_{\infty}^{2}}{\sigma^{2}}$

for some constant $C_{\inf}$. ∎

######  Lemma 6.

Let $p_{0}$ satisfy Condition F. Then
$\|p_{0}-\phi_{\sigma}*p_{0}\|_{\infty}\lesssim\sigma^{\alpha}$ for
sufficiently small $\sigma$.

###### Proof.

Suppose $\alpha\leq 1$. Let
$L=\sup_{(y,\mathbf{x})\neq(y^{\prime},\mathbf{x}^{\prime})}\frac{|p_{0}(y\mid\mathbf{x})-p_{0}(y^{\prime}\mid\mathbf{x}^{\prime})|}{\|(y-y^{\prime},\mathbf{x}-\mathbf{x}^{\prime})\|^{\alpha}}$,
which is finite by Condition F. Then:

$\displaystyle|p_{0}(y\mid\mathbf{x})-\phi_{\sigma}*p_{0}(y\mid\mathbf{x})|$
$\displaystyle=|p_{0}(y\mid\mathbf{x})-\int\phi_{\sigma}(u)p_{0}(y-u\mid\mathbf{x})\,du|$
$\displaystyle=|\int
p_{0}(y\mid\mathbf{x})\phi_{\sigma}(u)\,du-\phi_{\sigma}(u)p_{0}(y-u\mid\mathbf{x})\,du|$
$\displaystyle\leq\int|p_{0}(y\mid\mathbf{x})\phi_{\sigma}(u)-p_{0}(y-u\mid\mathbf{x})\phi_{\sigma}(u)|\,du$
$\displaystyle\leq L\int|u|^{\alpha}\,\phi_{\sigma}(u)\
du\lesssim\sigma^{\alpha}.$

For $1<\alpha\leq 2$, let
$\dot{p}_{0}(y\mid\mathbf{x})=\frac{d}{dy}p_{0}(y\mid\mathbf{x})$ and let
$L=\sup_{(y,\mathbf{x})\neq(y^{\prime},\mathbf{x}^{\prime})}\frac{|\dot{p}_{0}(y\mid\mathbf{x})-\dot{p}_{0}(y^{\prime}\mid\mathbf{x}^{\prime})|}{\|(y,\mathbf{x})-(y^{\prime},\mathbf{x}^{\prime})\|^{\alpha-1}}$,
which is finite by Condition F.
Then:

$\displaystyle
p_{0}(y-u\mid\mathbf{x})=p_{0}(y\mid\mathbf{x})-u\dot{p}_{0}(y^{\prime}\mid\mathbf{x})+u\dot{p}_{0}(y)-u\dot{p}_{0}(y\mid\mathbf{x})$

for some $y^{\prime}$ between $y$ and $y-u$, depending on $(y,\mathbf{x},u)$.
Integrating with respect to $\phi_{\sigma}(u)$, the term
$u\dot{p}_{0}(y\mid\mathbf{x})$ drops; then taking the absolute value we have:

$\displaystyle\bigg{|}\int(p_{0}(y-u\mid\mathbf{x})-p_{0}(y\mid\mathbf{x}))\phi_{\sigma}(u)\,du\bigg{|}$
$\displaystyle=\bigg{|}\int
u(\dot{p}_{0}(y\mid\mathbf{x})-\dot{p}(y^{\prime}\mid\mathbf{x}))\phi_{\sigma}(u)\,du\bigg{|}$

The left hand side becomes
$|p_{0}(y\mid\mathbf{x})-\phi_{\sigma}*p_{0}(y\mid\mathbf{x})|$ and the right
hand side is upper bounded by:

$\displaystyle\int|u(\dot{p}_{0}(y)-\dot{p}_{0}(y^{\prime}))\phi_{\sigma}(u)|\,du\leq
L\int|u|^{\alpha}\phi_{\sigma}(u)\,du\lesssim\sigma^{\alpha}$

yielding the stated result. ∎

######  Lemma 7.

Suppose Condition F holds. Then for sufficiently small $\sigma$ we have

$\displaystyle\int
p_{0}(y\mid\mathbf{x})\log\frac{p_{0}(y\mid\mathbf{x})}{\phi_{\sigma}*p_{0}(y\mid\mathbf{x})}\
dy\ F_{\mathbf{X}}(d\mathbf{x})\lesssim\sigma^{2\alpha}.$

###### Proof.

For fixed $\mathbf{x}$ we can bound the Kullback-Leibler divergence in terms
of the chi-squared divergence as

$\displaystyle\int
p_{0}(y\mid\mathbf{x})\log\frac{p_{0}(y\mid\mathbf{x})}{\phi_{\sigma}*p_{0}(y\mid\mathbf{x})}\
dy\
F_{\mathbf{X}}(d\mathbf{x})\leq\int\frac{(p_{0}(y\mid\mathbf{x})-\phi_{\sigma}*p_{0}(y\mid\mathbf{x}))^{2}}{\phi_{\sigma}*p_{0}(y\mid\mathbf{x})}\
dy\ F_{\mathbf{X}}(d\mathbf{x}).$

By Lemma 4 and Condition F, we have that the denominator is bounded away from
0 for small enough $\sigma$. Lemma 6 further tells us that the numerator is
$O(\sigma^{2\alpha})$, giving the result. ∎

#### Proof of Lemma 3

We first establish the second statement. Let
$h^{2}_{\mathbf{x}}(p,q)=\int\\{p(y\mid\mathbf{x})^{1/2}-q(y\mid\mathbf{x})^{1/2}\\}^{2}\
dy$ denote the usual Hellinger distance between $p(\cdot\mid\mathbf{x})$ and
$q(\cdot\mid\mathbf{x})$. We begin by noting that, by Cauchy-Schwarz,

$\displaystyle
p_{f_{1},\sigma_{1}}(y\mid\mathbf{x})\,p_{f_{2},\sigma_{2}}(y\mid\mathbf{x})\geq\left\\{\int_{0}^{1}\sqrt{\phi_{\sigma_{1}}(y-f_{1}(\mathbf{x},u))\phi_{\sigma_{2}}(y-f_{2}(\mathbf{x},u))}\right\\}^{2}.$

By Fubini’s Theorem and the usual formula for the affinity between normal
distributions we have

$\displaystyle\begin{split}h^{2}_{\mathbf{x}}(p_{f_{1},\sigma_{1}},p_{f_{2},\sigma_{2}})&\leq
1-\int_{0}^{1}\int\sqrt{\phi_{\sigma_{1}}(y-f_{1}(\mathbf{x},u))\phi_{\sigma_{2}}(y-f_{2}(\mathbf{x},u))}\
dy\ du\\\
&=\int_{0}^{1}1-\sqrt{\frac{2\sigma_{1}(\mathbf{x},u)\sigma_{2}(\mathbf{x},u)}{\sigma_{1}(\mathbf{x},u)^{2}+\sigma_{2}(\mathbf{x},u)^{2}}}\exp\left\\{-\frac{(f_{1}(\mathbf{x},u)-f_{2}(\mathbf{x},u))^{2}}{4(\sigma_{1}(\mathbf{x},u)^{2}+\sigma_{2}(\mathbf{x},u)^{2})}\right\\}\
du\end{split}$ (39)

By the triangle inequality we have
$h_{\mathbf{x}}^{2}(p_{f_{1},\sigma_{1}},p_{f_{2},\sigma_{2}})\leq
2\\{h_{\mathbf{x}}^{2}(p_{f_{1},\sigma_{1}},p_{f_{1},\sigma_{2}})+h^{2}_{\mathbf{x}}(p_{f_{1},\sigma_{2}}p_{f_{2},\sigma_{2}})\\}$.
Now, note from the inequality $1-e^{-|x|}\leq|x|$ that for any
$\sigma_{1},\sigma_{2}$ we have

$\displaystyle
1-\sqrt{\frac{2\sigma_{1}\sigma_{2}}{\sigma_{1}^{2}+\sigma_{2}^{2}}}\leq
1-\sqrt{\frac{\sigma_{1}^{2}\wedge\sigma_{2}^{2}}{\sigma_{1}^{2}\vee\sigma_{2}^{2}}}=1-\exp\\{-|\log\sigma_{1}-\log\sigma_{2}|\\}\leq|\log\sigma_{1}-\log\sigma_{2}|.$

Hence using (39) with $f_{1}$ in place of $f_{2}$ gives

$\displaystyle
h_{\mathbf{x}}^{2}(p_{f_{1},\sigma_{1}},p_{f_{1},\sigma_{2}})\leq\int_{0}^{1}|\log\sigma_{1}(\mathbf{x},u)-\log\sigma_{2}(\mathbf{x},u)|\
du\leq\|\log\sigma_{1}-\log\sigma_{2}\|_{\infty}.$

Next, we use (39) with $\sigma_{2}$ in place of $\sigma_{1}$ to get

$\displaystyle h_{\mathbf{x}}^{2}(p_{f_{1},\sigma_{2}},p_{f_{2},\sigma_{2}})$
$\displaystyle\leq\int_{0}^{1}1-\exp\left\\{-\frac{(f_{1}(\mathbf{x},u)-f_{2}(\mathbf{x},u))^{2}}{4\sigma_{2}(\mathbf{x},u)^{2}}\right\\}\
du$
$\displaystyle\leq\int_{0}^{1}\frac{(f_{1}(\mathbf{x},u)-f_{2}(\mathbf{x},u))^{2}}{4\sigma_{2}(\mathbf{x},u)^{2}}\
du$
$\displaystyle\leq\frac{\|f_{1}-f_{2}\|^{2}_{\infty}}{4\inf_{u}\sigma_{2}(\mathbf{x},u)^{2}}.$

Putting these together we get

$\displaystyle
h_{\mathbf{x}}^{2}(p_{f_{1},\sigma_{1}},p_{f_{2},\sigma_{2}})\leq
2\left\\{\|\log\sigma_{1}-\log\sigma_{2}\|_{\infty}+\frac{\|f_{1}-f_{2}\|^{2}_{\infty}}{4\inf_{u}\sigma_{1}(\mathbf{x},u)^{2}\wedge\sigma_{2}(\mathbf{x},u)^{2}}\right\\}.$

Integrating with respect to $F_{\mathbf{X}}(d\mathbf{x})$ and taking the
square root gives the result.

To prove the first bound, fix $\mathbf{x}$ and suppose $p_{f,\sigma}$ is such
that $\sigma(\mathbf{x},u)$ is constant,
$\sigma\in(\epsilon^{1/\alpha},2\epsilon^{1/\alpha})$, and
$\|f-f_{0}\|_{\infty}\leq\epsilon^{1+1/\alpha}$. Note that
$\phi_{\sigma}*p_{0}=p_{f_{0},\sigma}$. For any $p_{f,\sigma}$, applying the
triangle inequality we have

$\displaystyle h_{\mathbf{x}}^{2}(p_{0},p_{f,\sigma})\leq
2\\{h_{\mathbf{x}}^{2}(p_{0},p_{f_{0},\sigma})+h^{2}_{\mathbf{x}}(p_{f_{0},\sigma},p_{f,\sigma})\\}.$

Using the fact that $h_{\mathbf{x}}^{2}(p,q)\leq\int p\log(p/q)\ dy$, the
proof of Lemma 7 gives
$h_{\mathbf{x}}^{2}(p_{0},p_{f_{0},\sigma})\lesssim\sigma^{2\alpha}\lesssim\epsilon^{2}$
for sufficiently small $\epsilon$. Next, we have

$\displaystyle h^{2}_{\mathbf{x}}(p_{f_{0},\sigma},p_{f,\sigma})\leq
1-\exp\left\\{-\frac{\|f-f_{0}\|^{2}_{\infty}}{8\sigma^{2}}\right\\}\lesssim\frac{\|f-f_{0}\|^{2}_{\infty}}{\sigma^{2}}\lesssim\frac{\epsilon^{2+2/\alpha}}{\epsilon^{2/\alpha}}\lesssim\epsilon^{2}.$

Combining these facts and integrating with respect to
$F_{\mathbf{X}}(d\mathbf{x})$ gives
$h^{2}(p_{0},p_{f\sigma})\lesssim\epsilon^{2}$. An application of Lemma 8 of
Ghosal and van der Vaart [2007] and Lemma 5 gives
$K(p_{0}\|p_{f,\sigma})\lesssim\epsilon^{2}$ and
$V(p_{0}\|p_{f,\sigma})\lesssim\epsilon^{2}$ as well; let
$C_{\operatorname{KL}}$ denote the constant which makes
$\operatorname{KL}(p_{0}\|p_{f,\sigma})\vee V(p_{0}\|p_{f,\sigma})\leq
C^{2}_{\operatorname{KL}}\epsilon^{2}$. Hence
$p_{f,\sigma}\in\operatorname{KL}(C_{\operatorname{KL}}\epsilon)$, proving the
result.

###  8.4 Simulation Details

####  8.4.1 Implementation

Code for conducting these simulations can be found at:
https://github.com/vittorioorlandi/DR-BART.

Implementation details are as below:

  * •

DR-BART DR-BART was implemented in R using the Rcpp package.

  * •

DR-GP We implemented DR-GP via the rstan package on CRAN for interfacing with
Stan [Stan Development Team, 2020]. We followed Kundu and Dunson [2014], used
a squared exponential kernel for the Gaussian Process, and introduced
covariates appropriately.

  * •

SBART-DS The code for SBART-DS was graciously provided by Li et al.

  * •

PSBPM This method was run via the implementation on the author’s GitHub at:
https://github.com/david-dunson/probit-stick-breaking. Prior parameters were
as suggested in Chung and Dunson [2012].

  * •

DPMM The code for DPMM was taken from version 1.17 of the archived CRAN
package DPpackage, which can be found here:
https://cran.r-project.org/src/contrib/Archive/DPpackage/. Prior parameters
were left at their defaults.

All methods were run for 25,000 iterations of burn-in, after which 25,000
posterior samples were saved.

####  8.4.2 Predicted Coverage and Credible Band Width

Here, we provide additional information on 1. the predictive coverage and 2.
the credible band width of each method across various simulation settings. To
compute the predictive coverage, an additional $n=1000$ test points were
generated for each run of a simulation. The predicted coverage for a run is
the proportion of the test points that were contained in the 95% HDR
intervals; the reported values are averages across all runs of a simulation.
The reported widths are average widths – across runs of a simulation – of 95%
posterior credible bands within the 95% HDR region of the true density. That
is, we evaluate credible band width in a high density region of the data. The
predictive coverage below shows that all methods expect for PSBPM consistently
have good predictive coverage. The information on credible bands is useful in
conjunction with the coverage results in the main text, as it helps show which
of the methods that undercover tend to do so because they are unreasonably
confident in their estimates (e.g. DR-GP and PSBPM) versus those that simply
do not capture the shape of the density well enough (e.g. SBART-DS).

$x$ | DR-BART | DPMM | SBART-DS | DR-GP | PSBPM
---|---|---|---|---|---
0.1 | 0.95 | 0.95 | 0.93 | 0.96 | 0.63
0.5 | 0.95 | 0.95 | 0.96 | 0.96 | 0.75
0.8 | 0.99 | 0.97 | 0.97 | 0.99 | 0.00
Table 5: Mean predictive coverage of 95% predictive intervals, across runs of Simulation 1. $x$ | DR-BART | DPMM | SBART-DS | DR-GP | PSBPM
---|---|---|---|---|---
0.1 | 0.12 | 0.13 | 0.07 | 0.08 | 0.07
0.5 | 0.26 | 0.20 | 0.17 | 0.10 | 0.12
0.8 | 0.67 | 0.33 | 0.27 | 0.20 | 0.00
Table 6: Mean 95% credible band width within the true 95% HDR region, averaged across runs of Simulation 1. $x$ | DR-BART | DPMM | SBART-DS | DR-GP | PSBPM
---|---|---|---|---|---
0.1 | 0.31 | 0.1 | 0.07 | 0.05 | 0.08
0.5 | 0.53 | 0.1 | 0.17 | 0.07 | 0.12
0.8 | 1.55 | 0.3 | 0.27 | 0.20 | 0.00
Table 7: Mean 95% credible band width within the true 95% HDR region, averaged across runs of Simulation 2. $x$ | DR-BART | DPMM | SBART-DS | DR-GP | PSBPM
---|---|---|---|---|---
0.1 | 0.92 | 0.94 | 0.93 | 0.90 | 0.64
0.5 | 0.94 | 0.97 | 0.95 | 0.97 | 0.73
0.8 | 1.00 | 0.99 | 0.97 | 1.00 | 0.00
Table 8: Mean predictive coverage of 95% predictive intervals, across runs of Simulation 2. $x$ | DR-BART | DPMM | SBART-DS | DR-GP | PSBPM
---|---|---|---|---|---
0.1 | 0.95 | 0.95 | 0.94 | 0.96 | 0.61
0.5 | 0.95 | 0.97 | 0.96 | 0.97 | 0.78
0.8 | 1.00 | 0.96 | 0.97 | 0.99 | 0.00
Table 9: Mean predictive coverage of 95% predictive intervals, across runs of Simulation 3. $x$ | DR-BART | DPMM | SBART-DS | DR-GP | PSBPM
---|---|---|---|---|---
0.1 | 0.11 | 0.13 | 0.07 | 0.08 | 0.07
0.5 | 0.33 | 0.34 | 0.24 | 0.16 | 0.17
0.8 | 0.63 | 0.33 | 0.25 | 0.20 | 0.00
Table 10: Mean 95% credible band width within the true 95% HDR region,
averaged across runs of Simulation 3.

###  8.5 Applications

Code for performing these analyses can be found at:
https://github.com/vittorioorlandi/DR-BART. The data for the application on
returns to education can be found here. The data for the student growth
application can be found in the SGPdata R package.