# Theory of Posterior Concentration for
Generalized Bayesian Additive Regression Trees

Enakshi Saha
Department of Biostatistics
Harvard T.H. Chan School of Public Health
Boston, MA
<EMAIL_ADDRESS>

###### Abstract

Bayesian Additive Regression Trees (BART) are a powerful semiparametric
ensemble learning technique for modeling nonlinear regression functions.
Although initially BART was proposed for predicting only continuous and binary
response variables, over the years multiple extensions have emerged that are
suitable for estimating a wider class of response variables (e.g. categorical
and count data) in a multitude of application areas. In this paper we describe
a Generalized framework for Bayesian trees and their additive ensembles where
the response variable comes from an exponential family distribution and hence
encompasses a majority of these variants of BART. We derive sufficient
conditions on the response distribution, under which the posterior
concentrates at a minimax rate, up to a logarithmic factor. In this regard our
results provide theoretical justification for the empirical success of BART
and its variants.

_K_ eywords Bayesian additive regression trees, BART, Posterior concentration,
Minimax rate, Exponential family, Generalized linear models

##  1 Introduction

Additive ensemble of Bayesian trees [1, 2], more popularly known as Bayesian
additive regression trees (BART) [3] is a flexible semiparametric tool that
has been extremely successful in a multitude of high dimensional
classification and regression tasks. Aided by efficient software
implementations, (BART R package of [4], bartMachine R package of [5],
parallel BART of [6] and XBART of [7]), BART has thrived in a wide range of
application areas, including causal inference [8, 9, 10], interaction
detection [11], survival analysis [12], time series analysis [13, 14] and
variable selection [5, 15, 16, 17], to name a few. Even though BART was
initially proposed for predicting univariate continuous and binary response
variables, due to its flexibility and impressive performance, multiple
extensions have emerged over the subsequent years, that are suitable for both
univariate and multivariate prediction problems where the response variable is
of a wider variety (e.g. categorical and count data [18], heteroscedastic
responses [19]) and / or the target regression surface is of a constrained
nature (e.g. monotone BART [20], varying coefficient BART [14], BART with
targeted smoothing [21] etc.).

Despite a long history of empirical success, theoretical studies on Bayesian
trees and forests is a relatively new area of research. Recently emerging
results along this line are geared towards providing a theoretical perspective
on why these models have been so successful in a wide range of classification
and regression problems. Among the initial developments, [22] and [23]
demonstrated that the posterior concentration rate of BART equals to the
minimax rate up to a logarithmic factor for various tree priors. Built on
these findings, [24] derived a semiparametric Bernstein von-Mises theorem for
the BART estimator. Extensions of BART, adapted to various special function
types have also been studied from a theoretical perspective: [25] studied a
version of BART suitable for smooth function estimation; [26] conducted a
multiscale analysis of BART and [27] derived posterior concentration results
for anisotropic functions. In this paper we study the posterior concentration
rates of a generalized version of BART, thereby supplementing this newly
emerging area of research.

We formulate a Generalized BART (G-BART) model that extends the existing
theoretical developments in several directions. Firstly while existing results
focus on Gaussian response variables, we allow the response to come from an
exponential family distribution. Hence G-BART can be regarded as
semiparametric extensions of the widely popular ‘Generalized Linear Models’
(GLM) [28]. Many prominent Bayesian CART and BART models used in practice [2,
3, 18], including the traditional BART model [3], can be viewed as a special
case of this generalized extension. Therefore theoretical properties of these
conventional adaptations of BART can be studied as direct corollaries of
analogous properties for the G-BART model.

Secondly, existing results [23, 22, 25] build upon the assumption that the
underlying regression function is Hölder continuous. However given the
efficacy of BART models in a multitude of prediction problems with varying
degrees of complexity, the assumption of Hölder continuity seems too
restrictive. In this paper we demonstrate that similar posterior optimality
results can be obtained for non-smooth functions as well, such as step
functions and monotone functions, thus extending the theoretical findings on
BART beyond the assumption of Hölder continuity.

Finally, the BART model [3] approximate the regression functions through step
functions and assume that these step heights come from a Gaussian
distribution. All subsequent theoretical and empirical developments have
adopted this specification. In the G-BART setup we assume that the
distribution of these step heights belong to a broader family of distributions
that include both the Gaussian distribution and also some thicker tailed
distributions like Laplace. We demonstrate that the BART model maintains a
near-minimax posterior concentration rate, if the step heights come from any
of the distributions belonging to this broader family, thus providing a wide
range of distributional choices without sacrificing fast posterior
concentration. The theory also shows how important modelling choices such as
link functions can impact performance of the posterior and hence can serve as
a guide for empirical implementations as well.

This paper is organized as follows. In Section 2 we describe the generalized
BART model with the associated priors. Section 3 discusses the notion of
posterior concentration, followed by the main theoretical results on G-BART in
Section 4. Broader implications of these results are described in Section 5.
Finally, Section 6 concludes with a discussion. Proofs of the main theoretical
results are provided in the supplementary material.

###  1.1 Our contributions

To summarize our previous discussion, we now briefly highlight our key
contributions.

##### Response distribution:

We assume that the response variable comes from an exponential family
distribution and derive sufficient conditions on the response density under
which the posterior concentration rate of the BART model adapted to this
particular response type would be almost equal to the minimax rate. This
extends the existing theoretical results on BART for Gaussian regression.

##### Step Size distribution:

Instead of assigning a Gaussian distribution on the step-heights associated to
the BART model, we impose sufficient conditions on the cumulative distribution
function that guarantee a near-optimal posterior concentration rate. The
resulting family of distributions encompasses the Gaussian distribution along
with several thicker tailed distributions like Laplace, thus widening modeling
choices for empirical applications.

##### Types of functions:

The objective of BART models is to estimate unknown functions that
characterize the relationship between the response and the covariates. All
existing results on BART assume this function to be Hölder continuous. We
extend these results to the situations where the underlying function to be
estimated is either a monotone function or a step function supported on an
axes-paralleled partition. The results on step functions are particularly
important because posterior concentration rates for more general class of
functions can be built upon these, aided by the “simple function approximation
theorem” [29].

##### Empirical implications:

As we will see in Section 5, specific model choices such as the choice of link
functions can influence the posterior concentration rate of the G-BART model.
The results discussed in this paper can provide useful insights into selecting
link functions that provide faster concentration rates of the posterior,
possibly leading to better empirical performance.

###  1.2 Notations:

For any two real numbers $a$ and $b$, $a\vee b$ will denote the maximum of $a$
and $b$. The notations $\gtrsim$ and $\lesssim$ will stand for “greater than
or equal to up to a constant” and “less than or equal to up to a constant”,
respectively. The symbol $P_{f}$ will abbreviate $\int fdP$ and
$\mathbb{P}_{f}^{(n)}=\prod_{i=1}^{n}\mathbb{P}_{f}^{i}$ will denote the
$n$-fold product measure of the $n$ independent observations, where the $i$-th
observation comes from the distribution $P_{f}^{i}$. Let
$h(f,g)=\left(\int(\sqrt{f}-\sqrt{g})^{2}d\mu\right)^{1/2}$ and $K(f,g)=\int
f\log(f/g)d\mu$ denote the Hellinger distance and the Kullback-Leibler
divergence, respectively between any two non-negative densities $f$ and $g$
with respect to a measure $\mu$. We define another discrepancy measure
$V(f,g)=\int f\left(\log(f/g)\right)^{2}d\mu$. Finally, for any set of real
vectors $\boldsymbol{X}_{1},\dots,\boldsymbol{X}_{n}\in\mathbb{R}^{q}$ of size
$n$, define the average discrepancy measures
$H_{n}(f,g)=\frac{1}{n}\sum_{i=1}^{n}H\left(f(\boldsymbol{X}_{i}),g(\boldsymbol{X}_{i})\right)$,
$K_{n}(f,g)=\frac{1}{n}\sum_{i=1}^{n}K\left(f(\boldsymbol{X}_{i}),g(\boldsymbol{X}_{i})\right)$
and
$V_{n}(f,g)=\frac{1}{n}\sum_{i=1}^{n}V\left(f(\boldsymbol{X}_{i}),g(\boldsymbol{X}_{i})\right)$,
where $f(\theta)$ and $g(\theta)$ denote the densities $f$ and $g$ with
respect to parameter $\theta$. Also, for any $L_{p}$ norm
$\left\|\cdot\right\|_{p}$, define the average norm
$\left\|f-g\right\|_{p,n}=\frac{1}{n}\sum_{i=1}^{n}\left\|f-g\right\|_{p}$.

##  2 The Generalized BART Prior

The BART method of [3] is a prominent example of Bayesian ensemble learning,
where individual shallow trees are entwined together into a forest, that is
capable of estimating a wide variety of nonlinear functions with exceptional
accuracy, while simultaneously accounting for different orders of interactions
among the covariates. Building upon BART, we describe a generalized model,
where the response variable is assumed to come from an exponential family
distribution. For continuous Gaussian response variables, this generalized
BART model reduces to the original BART prior of [3].

The data setup under consideration consists of
$\boldsymbol{Y}_{i}=(y_{i1},\dots,y_{ip})^{\prime}\in\mathbb{R}^{p}$, a set of
$p$-dimensional outputs, and
$\boldsymbol{X}_{i}=(x_{i1},\dots,x_{iq})^{\prime}\in[0,1]^{q}$, a set of $q$
dimensional inputs for $1\leq i\leq n$. We assume $\boldsymbol{Y}$ follows
some distribution in the exponential family with density of the following
form:

$P_{f_{0}}(\boldsymbol{Y}\,|\>\boldsymbol{X})=h(\boldsymbol{Y})g\left[f_{0}(\boldsymbol{X})\right]\exp\left[\eta\left(f_{0}(\boldsymbol{X})\right)^{T}T(\boldsymbol{Y})\right],$
(1)

where $h:\mathbb{R}^{p}\rightarrow\mathbb{R}$,
$g:\mathbb{R}\rightarrow\mathbb{R}$,
$\eta:\mathbb{R}^{p}\rightarrow\mathbb{R}^{J}$,
$T:\mathbb{R}^{p}\rightarrow\mathbb{R}^{J}$ for some integer $J$ and
$f_{0}:\mathbb{R}^{q}\rightarrow\mathbb{R}^{D}$, for some integer $D$, are all
real valued functions. Among these functions, $h$, $g$, $\eta$ and $T$ are
usually known depending on the nature of the response $\boldsymbol{Y}$. The
function $f_{0}$, connecting the input $\boldsymbol{X}$ with the output
$\boldsymbol{Y}$, is the only unknown function and estimating this function is
the primary objective of the G-BART estimator.

We assume that $f_{0}$ is an unconstrained function, i.e. the range of $f_{0}$
is the entire space $\mathbb{R}^{D}$ for some integer $D$. A suitable link
function $\Psi(\cdot)$ is used to transform $f_{0}$ to the natural parameter
of the distribution of $\boldsymbol{Y}$, which is often constrained. For
example, for the binary classification problem, $\boldsymbol{Y}\sim
Bernoulli\left(p(\boldsymbol{X})\right)$. Here the natural parameter
$p(\boldsymbol{X})\in(0,1)$ is restricted and hence we can use
$\Psi(z)=\frac{1}{1+\exp(-z)}$, the logistic function (or a probit function,
as in [3]) to map the unconstrained function $f_{0}(\boldsymbol{X})$ to the
natural parameter $p(\boldsymbol{X})$. There are usually several different
choices for the link function. As we will see in Section 5, the BART estimator
might have different posterior concentration rates depending on which link
function is used to transform the function $f_{0}$ to the natural parameter of
the response distribution.

The univariate regression and the two-class classification problem considered
in the original BART paper [3] and many of its important extensions, such as
the multi-class classification and the log-linear BART [18] for categorical
and count responses can be formulated as special cases of (1). The specific
forms of the functions $h,g,\eta$ and $T$ for continuous regression and multi-
class classification are given in Table 1.

Table 1:  Univariate Regression (column 2) and Multi-class Classification (column 3), as special cases of the Generalized BART model. $\Phi$ denotes the $Softmax$ function and $\mathcal{M}(\cdot)$ denotes the $Multinomial(1;\cdot)$ distribution. $\left(\\{\mathbb{I}\\{Y=i\\}\\}_{i=1}^{p}\right)^{\prime}$ denotes the row vector where the $i$-th coordinate equals to one if $\boldsymbol{Y}$ belongs to class $i$ and zero otherwise.  Response ($\boldsymbol{Y}$) | Continuous | Categorical
---|---|---
Dist.($\boldsymbol{Y}$) | $\mathcal{N}\left(f_{0}(\boldsymbol{X}),\sigma^{2}\right)$ | $\mathcal{M}\left(\Phi(f_{0}(\boldsymbol{X}))\right)$
$h(\boldsymbol{Y})$ | $1/\sqrt{2\pi\sigma}$ | 1
$g\left(f_{0}(\boldsymbol{X})\right)$ | $\exp\left(-f_{0}(\boldsymbol{X})^{2}/\sigma^{2}\right)$ | 1
$\eta\left(f_{0}(\boldsymbol{X})\right)$ | $\left(f_{0}(\boldsymbol{X}),1\right)$ | $f_{0}(\boldsymbol{X})$
$T(\boldsymbol{Y})$ | $\left(2Y/\sigma^{2},-Y^{2}/\sigma^{2}\right)$ | $\left(\\{\mathbb{I}\\{Y=i\\}\\}_{i=1}^{p}\right)^{\prime}$
$f_{0}(\boldsymbol{X})$ | $\mathbb{R}^{q}\rightarrow\mathbb{R}$ | $\mathbb{R}^{q}\rightarrow\mathbb{R}^{p-1}$

Next a regression tree is used to reconstruct the unknown function
$f_{0}:\mathbb{R}^{q}\rightarrow\mathbb{R}^{D}$ via a mapping
$f_{\mathcal{T},{\boldsymbol{\beta}}}:[0,1]^{q}\rightarrow\mathbb{R}^{D}$ so
that $f_{\mathcal{T},{\boldsymbol{\beta}}}(\boldsymbol{X})\approx
f_{0}(\boldsymbol{X})$ for
$\boldsymbol{X}\notin\\{\boldsymbol{X}_{i}\\}_{i=1}^{n}$. Each such mapping is
essentially a step function of the form

$f_{\mathcal{T},{\boldsymbol{\beta}}}(\boldsymbol{X})=\sum_{k=1}^{K}\beta_{k}\mathbb{I}(\boldsymbol{X}\in\Omega_{k})$
(2)

supported on a tree-shaped partition $\mathcal{T}=\\{\Omega_{k}\\}_{k=1}^{K}$
and specified by a vector of step heights
$\boldsymbol{\beta}=(\beta_{1},\dots,\beta_{K})^{\prime}$. The vector
$\beta_{k}\in\mathbb{R}^{p}$ represents the value of the expected response
inside the $k$-th cell of the partition $\Omega_{k}$.

Bayesian additive trees consist of an ensemble of multiple shallow trees, each
of which is intended to be a weak learner, geared towards addressing a
slightly different aspect of the prediction problem. These trees are then
woven into an additive forest mapping of the form

$f_{\mathcal{E},{\boldsymbol{B}}}(\boldsymbol{x})=\sum_{t=1}^{T}f_{\mathcal{T}_{t},{\boldsymbol{\beta}}_{t}}(\boldsymbol{x}),$
(3)

where each $f_{\mathcal{T}_{t},{\boldsymbol{\beta}}_{t}}(\boldsymbol{x})$ is
of the form (2), $\mathcal{E}=\\{\mathcal{T}_{1},\dots,\mathcal{T}_{T}\\}$ is
an ensemble of $T$ trees and
$\boldsymbol{B}=\\{\boldsymbol{\beta}_{1},\dots,\boldsymbol{\beta}_{T}\\}^{\prime}$
is a collection of jump sizes corresponding to the $T$ trees.

Since each individual member of the approximating space is a step function of
the form (3), supported on a Bayesian additive forest, the prior distribution
should include three components: (i) a prior $\pi(T)$ on the number of trees
$T$ in the ensemble, (ii) a prior on individual tree partitions
$\pi(\mathcal{T})$ and their collaboration within the ensemble and (iii) given
a single tree partition $\mathcal{T}$, a prior
$\pi(\boldsymbol{\beta}\,|\>\mathcal{T})$ has to be imposed on the individual
step heights $\boldsymbol{\beta}$.

In this paper we follow the recommendation by [3] and assume the number of
trees $T$ to be fixed at a large value (e.g. $T=200$ for regression and $T=50$
for classification). This is equivalent to assigning a degenerate prior
distribution on $T$, where all probability mass is concentrated on a single
positive integer. Alternatively, one can also assign a prior with higher
dispersion, as in [23] and [25] and replicate the steps of the proofs provided
in the appendix with minor modifications.

Given the total number of trees in the ensemble, individual trees are assumed
to be independent and identically distributed with some distribution
$\pi(\mathcal{T})$. This reduces the prior on the ensemble to be of the form

$\pi(\mathcal{E},\boldsymbol{B})=\prod_{t=1}^{T}\pi(\mathcal{T}_{t})\pi(\boldsymbol{\beta}_{t}\,|\>\mathcal{T}_{t}),$
(4)

where $\pi(\mathcal{T}_{t})$ is the prior probability of a partition
$\mathcal{T}_{t}$, while $\pi(\boldsymbol{\beta}_{t}\,|\>\mathcal{T}_{t})$ is
the prior distribution over the jump sizes. The specific forms of the priors
$\pi(\mathcal{T})$ and $\pi(\boldsymbol{\beta}\,|\>\mathcal{T})$ are described
below.

###  2.1 Prior on partitions

We consider two distinct prior distributions on the partitions
$\pi(\mathcal{T})$ proposed by [1] and [2] respectively. The posterior
concentration results discussed in Section 4 are applicable to both these
priors. [1] specifies the prior over trees implicitly as a tree generating
stochastic process, described as follows:

  1. 1.

Start with a single leaf (a root node) encompassing the entire covariate
space.

  2. 2.

Split a terminal node, say $\Omega$, with a probability

$p_{split}(\Omega)\propto\alpha^{-d(\Omega)}\text{ for some $0<\alpha<1/2$.}$
(5)

where $d(\Omega)$ is the depth of the node $\Omega$ in the tree architecture.
This choice, motivated by [22], is slightly different from the original prior
of [1]111The reason behind this modification is that the original BART prior
of [3] does not decay at a fast enough rate. However since we examine only
sufficient (but not necessary) conditions for optimal posterior concentration,
our results do not guarantee that the original prior is inherently worse than
the modified prior In fact, empirical results indicate otherwise. The original
BART prior will be examined in future work.

  3. 3.

If the node $\Omega$ splits, assign a splitting rule and create left and right
children nodes. The splitting rule consists of picking a split variable $j$
uniformly from available directions $\\{1,\dots,p\\}$ and picking a split
point $c$ uniformly from available data values $x_{1j},\dots,x_{nj}$.

A description of the prior proposed by [2] is given in Section A.1 in the
supplementary material.

###  2.2 Prior on step heights

We impose a broad class of priors on the step heights that incorporate the
corresponding component of the classical BART model as a special case. Given a
tree partition $\mathcal{T}_{t}$ with $K_{t}$ steps, [3] considers identically
distributed independent Gaussian jumps with mean $0$ and variance
$\sigma^{2}$. In the G-BART set-up we assume that the $j$-th step height of
the $t$-th tree, $\beta_{tj}\stackrel{{\scriptstyle i.i.d}}{{\sim}}F_{\beta}$,
where $F_{\beta}$ is any general distribution with the following property: for
some constants $C_{1},C_{2},C_{3}$ such that $C_{1}>0$, $0<C_{2}\leq 2$ and
$C_{3}>0$,

$F_{\beta}(\left\|\beta\right\|_{\infty}\leq
t)\gtrsim\left(e^{-C_{1}t^{C_{2}}}t\right)^{p}\quad\text{for $0<t\leq 1$}$ (6)

and

$F_{\beta}(\left\|\beta\right\|_{\infty}\geq t)\lesssim
e^{-C_{3}t}\quad\text{for $t\geq 1$}$ (7)

where $\left\|\cdot\right\|_{\infty}$ represents the $L_{\infty}$ norm and
$F_{\beta}(\left\|\beta\right\|_{\infty}\geq t)$ denotes the tail probability
of the distribution on the step heights $\beta\in\mathbb{R}^{p}$. Both the
multivariate Gaussian and the multivariate Laplace distribution come from this
family of distributions and so do any sub-Gaussian distributions. A proof of
these statements is provided in the appendix. We will see in Section 4.1 and
Section 4.3 that these conditions are sufficient to guarantee that the G-BART
estimator has a near-optimal posterior concentration rate.

However we should note that the conditions (6)-(7), although sufficient, are
not necessary conditions and distributional assumptions on the step sizes that
do not satisfy these conditions, might still guarantee a near-optimal
posterior concentration rate. For such an example, please refer to the
‘classification with Dirichlet step-heights’ in Section 5.

##  3 Posterior Concentration

Posterior concentration statements are a prominent artifact in Bayesian
nonparametrics, where the primary motivation is to examine the quality of a
Bayesian procedure, by studying the learning rate of its posterior, i.e. the
rate at which the posterior distribution, centralizes around the truth as the
sample size $n\rightarrow\infty$. In empirical settings, posterior
concentration results have often influenced the proposal and fine-tuning of
priors. Oftentimes seemingly unremarkable priors give rise to capricious
outcomes, specially in the infinite-dimensional parameter spaces, such as the
one considered here ([30], [31]) and designing well-behaved priors turn out to
be of utmost importance, thus further reinstating the importance of posterior
concentration statements.

The Bayesian approach proceeds by imposing a prior measure $\Pi(\cdot)$ on
$\mathcal{F}$, the set of all estimators of $f_{0}$. For the G-BART models
this corresponds to the set of all step functions supported on an additive
ensemble of Bayesian trees. Given observed data
$\boldsymbol{Y}^{(n)}=(Y_{1},\dots,Y_{n})^{\prime}$, the inference about
$f_{0}$ is solely dependent on the posterior distribution

$\Pi(A\mid\boldsymbol{Y}^{(n)})=\frac{\int_{A}\prod_{i=1}^{n}\Pi_{f}(Y_{i}\mid\boldsymbol{X}_{i})\mathrm{d}\,\Pi(f)}{\int\prod_{i=1}^{n}\Pi_{f}(Y_{i}\mid\boldsymbol{X}_{i})\mathrm{d}\,\Pi(f)}\quad\forall
A\in\mathcal{B}$

where $\mathcal{B}$ is a $\sigma$-field on $\mathcal{F}$ and where
$\Pi_{f}(Y_{i}\,|\>\boldsymbol{X}_{i})$ is the conditional likelihood function
for the output $Y_{i}$, given the covariates $\boldsymbol{X}_{i}$, under the
parameterization $f$.

Ideally under a suitable prior, the posterior should put most of its
probability mass around a small neighborhood of the true function and as the
sample size increases, the diameter of this neighborhood should go to zero at
a fast pace. Formally speaking, for a given sample size $n$, if we examine an
$\varepsilon_{n}$-neighborhood of the true function
$\mathcal{A}_{\varepsilon_{n}}$, for some $\varepsilon_{n}\rightarrow 0$ and
$n\varepsilon_{n}^{2}\rightarrow\infty$, we should expect

$\Pi(\mathcal{A}_{\varepsilon_{n}}^{c}\,|\>\boldsymbol{Y}^{(n)})\rightarrow
0\quad\text{in $\mathbb{P}_{f_{0}}^{(n)}$-probability as
$n\rightarrow\infty$,}$ (8)

where $\mathcal{A}_{\varepsilon_{n}}^{c}$ denotes the complement of the
neighborhood $\mathcal{A}_{\varepsilon_{n}}$.

In the context of G-BART, given observed data
$\boldsymbol{Y}^{(n)}=(\boldsymbol{Y}_{1},\dots,\boldsymbol{Y}_{n})^{\prime}$,
we are interested in evaluating whether the posterior concentrates around the
true likelihood $\mathbb{P}_{f_{0}}^{(n)}=\prod_{i=1}^{n}P_{f_{0}}^{i}$ at a
near-minimax rate, where
$P_{f_{0}}^{i}=P_{f_{0}}(\boldsymbol{Y}_{i}\,|\>\boldsymbol{X}_{i})$ is of the
form (1), for $i=1,\dots,n$. Following the suggestions of [32], we look at the
smallest $H_{n}$-neighborhoods around $\mathbb{P}_{f_{0}}^{(n)}$ that contain
the bulk of the posterior probability. Specifically, for a diameter
$\varepsilon>0$ define

$\mathcal{A}_{\varepsilon}=\\{f\in\mathcal{F}:H_{n}(P_{f},P_{f_{0}})\leq\varepsilon\\}$
(9)

Theorem 4 of [32] demonstrates that the statement (8) can be proved by
verifying three sufficient conditions. The first condition, henceforth
referred to as the “entropy condition” specifies that

$\displaystyle\sup_{\varepsilon>\varepsilon_{n}}\log
N\left(\tfrac{\varepsilon}{36};\mathbb{F}_{n}\cap\mathcal{A}_{\varepsilon};H_{n}\right)\lesssim
n\,\varepsilon_{n}^{2},$ (C1)

where $N(\varepsilon;\Omega;d)$ denotes the $\varepsilon$-covering number of a
set $\Omega$ for a semimetric $d$, i.e. the minimal number of $d$-balls of
radius $\varepsilon$ needed to cover the set $\Omega$ and
$\\{\mathbb{F}_{n}\\}_{n\geq 1}$ denotes an increasing sequence of
approximating sieves. The sequence of sieves used in this paper is described
in the appendix.

The second condition requires that the prior puts enough mass around the true
likelihood $\mathbb{P}_{f_{0}}^{(n)}$, meaning that for a given sample size
$n\in\mathbb{N}\setminus\\{0\\}$ and for some $d>2$,

$\displaystyle{\Pi(f\in\mathcal{F}:K_{n}(f,f_{0})\vee
V_{n}(f,f_{0})\leq\varepsilon_{n}^{2})}\gtrsim
e^{-d\,n\,\varepsilon_{n}^{2}},$ (C2)

where $K_{n}$ and $V_{n}$ are the Kullback-Leibler divergence and the
variation, averaged over the observed data points.

The final condition, referred to as the “prior decay rate condition”
stipulates that the sequence of sieves $\mathbb{F}_{n}\uparrow\mathcal{F}$
captures the entire parameter space with increasing accuracy, in the sense
that the complementary space $\mathcal{F}\backslash\mathbb{F}_{n}$ has
negligible prior probability mass for large values of $n$.

$\displaystyle\Pi(\mathcal{F}\backslash\mathbb{F}_{n})=o(e^{-(d+2)\,n\,\varepsilon_{n}^{2}})$
(C3)

The results of type (8) quantify not only the typical distance between a point
estimator (posterior mean/median) and the truth, but also the typical spread
of the posterior around the truth and hence are stronger than ‘posterior
consistency’ statements. These results are usually the first step towards
further uncertainty quantification statements such as semiparametric
Bernstein-von Mises theorem [33].

##  4 Main Results

In this section we describe our main theoretical findings, which describe the
posterior concentration rates of the generalized Bayesian trees and their
additive ensembles (G-BART), when the true function $f_{0}$ connecting the
response $\boldsymbol{Y}$ with the covariates $\boldsymbol{X}$, is either (a)
a step function (Theorem 4.1), or (b) a monotone function (Theorem 4.3), or
(c) a $\nu$-Hölder continuous function with $0<\nu\leq 1$ (Theorem 4.4). We
make two important assumptions: the first assumption (subsequently referred to
as Assumption 1), given below restricts the distribution of the response
variable $\\{\boldsymbol{Y}_{1},\dots,\boldsymbol{Y}_{n}\\}\in\mathbb{R}^{p}$
to a specific class of exponential family distributions while the second
assumption (subsequently referred to as Assumption 2) concerns the spread of
the covariates
$\\{\boldsymbol{X}_{1},\dots,\boldsymbol{X}_{n}\\}\in\mathbb{R}^{q}$.

##### Assumption 1:

Let $\boldsymbol{Y}_{1},\dots,\boldsymbol{Y}_{n}\sim P_{f}$, where $P_{f}$
denotes a probability density function of the form (1), such that, $\eta(z)=z$
and there exists strictly increasing positive sequences
$\\{C^{n}_{g}\\}_{n\geq 1}$ and $\\{C^{n}_{\beta}\\}_{n\geq 1}$, such that

$\left|\frac{\nabla g(\boldsymbol{\beta})}{g(\boldsymbol{\beta})}\right|\leq
C^{n}_{g}\textbf{1}_{p},\quad\forall\boldsymbol{\beta}\in
B_{n}=\left\\{\boldsymbol{\beta}:\left\|\boldsymbol{\beta}\right\|_{\infty}\leq
C^{n}_{\beta}\right\\},$ (10)

where $\textbf{1}_{p}=(1,\dots,1)\in\mathbb{R}^{p}$ denotes a $p$-dimensional
vector of ones and $\nabla g$ denotes the vector of partial derivatives. We
assume $\\{C^{n}_{g}\\}\vee\\{C^{n}_{\beta}\\}\lesssim n^{M}$ for some $M>0$.
The significance is that the function $g(\cdot)$ should not change too
rapidly, and the higher the sample size the larger the rate of change is
allowed. The above assumption is satisfied by most distributions commonly used
in the regression and classification settings, as will be demonstrated in
Section 5.

##### Assumption 2:

For a k-d tree partition,
$\smash{\widehat{\mathcal{T}}}=\\{\smash{\widehat{\Omega_{k}}}\\}$, with
$K=2^{ps}$-many leaves, the dataset
$\\{\boldsymbol{X}_{1},\dots,\boldsymbol{X}_{n}\\}$ satisfies the following
condition: for any nonnegative integer $s$, there exists some large enough
constant $M>0$ such that

$\max_{1\leq k\leq
K}\text{diam}(\smash{\widehat{\Omega_{k}}})<M\sum_{k=1}^{K}\mu(\Omega_{k})\text{diam}(\smash{\widehat{\Omega_{k}}}),$
(11)

where
$\mu(\Omega_{k})=\frac{1}{n}\sum_{i=1}^{n}\mathbb{I}\\{\boldsymbol{X}_{i}\in\Omega_{k}\\}$
denotes the proportion of observations in the cell $\Omega_{k}$ and
$\text{diam}(\smash{\widehat{\Omega_{k}}})=\max_{\boldsymbol{x},\boldsymbol{y}\in\Omega_{k}}\left\|\boldsymbol{x}-\boldsymbol{y}\right\|_{2}$
denotes the spread of the cell $\Omega_{k}$ with respect to the $L_{2}$-norm.

###  4.1 Results on Step-Functions

Let us suppose $f_{0}$ is a step function supported on an axes-paralleled
partition $\\{\Omega_{k}\\}_{k=1}^{K_{0}}$. For any such step function
$f_{0}$, we define the complexity of $f_{0}$, as the smallest $K$ such that
there exists a partition $\\{\Omega_{k}\\}_{k=1}^{K}$ with $K$ cells, for
which the step function
$f(x)=\sum_{k=1}^{K}\beta_{k}\mathbb{I}\\{x\in\Omega_{k}\\}$ can approximate
$f_{0}$ without any error, for some step heights
$(\beta_{1},\dots,\beta_{K})\in\mathbb{R}^{K}$. This complexity number,
denoted by $K_{f_{0}}$, depends on the true number of step $K_{0}$, the
diameter of the intervals $\\{\Omega_{k}\\}_{k=1}^{K_{0}}$, and the number of
covariates $q$. The actual minimax rate for approximating such piecewise-
constant functions $f_{0}$ with $K_{0}>2$ pieces, is
$n^{-1/2}\sqrt{K_{0}\log{(n/K_{0})}}$ [34]. The following theorem shows that
the posterior concentration rate of G-BART is almost equal to the minimax
rate, except that $K_{0}$ gets replaced by $K_{f_{0}}$. The discrepancy is an
unavoidable consequence of the fact that the true number of steps $K_{0}$ is
unknown. Had this information been available, the G-BART estimator would have
attained the exact minimax rate.

######  Theorem 4.1.

If we assume that the distribution of the step-sizes satisfies (6) and (7),
then under Assumptions 1 and 2 with $q\lesssim\sqrt{\log n}$, the generalized
BART estimator satisfies the following property:

If $f_{0}$ is a step-function, supported on an axes-paralleled partition, with
complexity $K_{f_{0}}\lesssim\sqrt{n}$ and
$\left\|f_{0}\right\|_{\infty}\lesssim\sqrt{\log n}$, then with
$\varepsilon_{n}=n^{-1/2}\sqrt{K_{f_{0}}\log^{2\gamma}{\left(n/K_{f_{0}}\right)}}$
and $\gamma>1/2$,

$\Pi\left(f\in\mathcal{F}:H_{n}(\mathbb{P}_{f},\mathbb{P}_{f_{0}})>\varepsilon_{n}\mid\boldsymbol{Y}^{(n)}\right)\to
0,$

in $\mathbb{P}_{f_{0}}^{(n)}$-probability, as $n,q\to\infty$.

###### Proof.

Proof is given in the appendix. ∎

###  4.2 Results on Monotone Functions

An important implication of Theorem 4.1 is that posterior concentration
results on step functions can potentially build the foundation for similar
results on broader class of functions, aided by the “simple function
approximation theorem” [29], which states that for any measurable function $f$
on $\mathcal{E}\subseteq\mathbb{R}^{q}$, there exists a sequence of step
functions $\\{f_{k}\\}$ which converges point-wise to $f$ almost everywhere
[29]. As a corollary to this theorem, we can derive the following result on
the set of all monotone functions. A function
$f_{0}:\mathbb{R}^{q}\rightarrow\mathbb{R}$ is defined as monotone increasing
(or decreasing) if $f_{0}(\boldsymbol{x}_{1})\geq f_{0}(\boldsymbol{x}_{2})$
(or, $f_{0}(\boldsymbol{x}_{1})\leq f_{0}(\boldsymbol{x}_{2})$) for all
$\boldsymbol{x}_{1},\boldsymbol{x}_{2}$ such that every coordinate of
$\boldsymbol{x}_{1}$ is greater than or equal to the corresponding coordinate
of $\boldsymbol{x}_{2}$.

######  Lemma 4.2.

Any monotone bounded function $f_{0}$ can be approximated with arbitrary
precision $\varepsilon$, by a step function supported on a k-d tree partition
with number of leaves $K_{f_{0}}(\varepsilon)\geq\lceil 1/\varepsilon\rceil$.
We define $K_{f_{0}}(\varepsilon)$ to be the complexity of the monotone
function $f_{0}$ with respect to $\varepsilon>0$.

The complexity $K_{f_{0}}(\varepsilon)$ also depends on the dimension of the
domain $q$ as well as on the magnitude of the true function
$\left\|f_{0}\right\|_{\infty}$. This paves the way for deriving the posterior
concentration rate of G-BART when the true function $f_{0}(\cdot)$ connecting
the covariates $\boldsymbol{X}$ with a univariate response $\boldsymbol{Y}$ is
a monotone function. The minimax rate of estimation for such densities is
$n^{-1/(2+q)}$ [35]. The following theorem states that the posterior
concentration rate of G-BART equals to this optimum rate up to a logarithmic
function, provided that the magnitude of the true function $f_{0}$ is not “too
large”.

######  Theorem 4.3.

If we assume that the distribution of the step-sizes satisfies (6) and (7),
then under Assumptions 1 and 2 with $q\lesssim\sqrt{\log n}$, the generalized
BART estimator satisfies the following property:

If the true function $f_{0}:\mathbb{R}^{q}\rightarrow\mathbb{R}$ is monotonic
on every coordinate, with $\left\|f_{0}\right\|_{\infty}\lesssim\sqrt{\log
n}$, then with $\varepsilon_{n}=n^{-1/(2+q)}\sqrt{\log n}$,

$\Pi\left(f\in\mathcal{F}:H_{n}(\mathbb{P}_{f},\mathbb{P}_{f_{0}})>\varepsilon_{n}\mid\boldsymbol{Y}^{(n)}\right)\to
0,$

in $\mathbb{P}_{f_{0}}^{(n)}$-probability, as $n,q\to\infty$.

###### Proof.

The first step of the proof involves finding an approximating step-function
$\smash{\widehat{f_{0}}}$ by Lemma 4.2, such that
$\left\|f_{0}-\smash{\widehat{f}}_{0}\right\|_{2,n}<\varepsilon_{n}/2$. The
rest of the proof follows by retracing the steps as in the proof of Theorem
4.4 given in the supplementary material. ∎

The above result demonstrates that the Generalized BART model adapts to
monotonic patterns in the true function $f_{0}$, without any additional prior
assumptions.

###  4.3 Results on Hölder Continuous Functions

This section describes the posterior concentration results on G-BART when the
true function $f_{0}$ connecting $\boldsymbol{X}$ with $\boldsymbol{Y}$ is a
$\nu$-Hölder continuous function with $0<\nu\leq 1$. [23] and [22] proved that
the posterior concentration rates of the BART model (under the priors of [2]
and [3] respectively) are equal to $n^{-\alpha/(2\alpha+q)}$, the minimax rate
of estimation for such functions [36], except for a logarithmic factor. These
results can be derived as direct corollaries of the following theorem for
G-BART, when $\boldsymbol{Y}$ is a univariate continuous response and the
step-sizes are assumed to follow a Gaussian distribution.

######  Theorem 4.4.

If we assume that the distribution of the step-sizes satisfies (6) and (7),
then under Assumptions 1 and 2 with $q\lesssim\sqrt{\log n}$, the generalized
BART estimator satisfies the following property:

If $f_{0}$ is a $\nu$-Hölder continuous function with $0<\nu\leq 1$, where
$\|f_{0}\|_{\infty}\lesssim\sqrt{\log n}$, then with
$\varepsilon_{n}=n^{-\alpha/(2\alpha+q)}\sqrt{\log n}$,

$\Pi\left(f\in\mathcal{F}:H_{n}(\mathbb{P}_{f},\mathbb{P}_{f_{0}})>\varepsilon_{n}\mid\boldsymbol{Y}^{(n)}\right)\to
0,$

in $\mathbb{P}_{f_{0}}^{(n)}$-probability, as $n,q\to\infty$.

###### Proof.

Proof is given in the appendix. ∎

##### Remark:

Interestingly, the posterior concentration rates derived in Theorems 4.1-4.4,
do not depend on the number of trees $T$ in the generalized BART ensemble. In
other words the concentration rate is equally valid for a single tree (i.e.
$T=1$), as well as for tree ensembles (i.e. $T>1$), when the true regression
function $f_{0}$ is $\nu$-Hölder continuous with $0<\nu\leq 1$. However as has
been seen in multiple empirical applications [3], Bayesian forests consisting
of multiple trees provide superior out-of-sample predictive performance,
compared to a single tree, the reason being that multiple weak tree learners,
when woven together into a forest, can accommodate a wider class of
partitions, as opposed to a single tree.

This phenomenon can be reinforced by theoretical results, such as Theorem 6.1
of [23]. When the true function $f_{0}$ is of the form
$f_{0}=\sum_{t=1}^{T_{0}}f_{0}^{t}$, where $f_{0}^{t}$ is a $\nu_{t}$-Hölder
continuous function, with $0\leq\nu^{t}\leq 1$, a forest with multiple trees
have a posterior concentration rate equal to
$\varepsilon_{n}^{2}=\sum_{t=1}^{T_{0}}n^{-2\nu_{t}/(2\nu_{t}+p)}\log n$,
provided $T_{0}\lesssim n$, whereas single regression trees fail to recognize
the additive nature of the true function and attain a slower concentration
rate. A similar result is presented in Theorem 4 of [25], under a kernel-
smoothed version of the BART prior.

Although the BART prior considered by [23] is fundamentally different from the
classical BART prior [3] considered here, their result on additive functions
can be replicated in the present set up as well, provided we allow the number
of trees $T$ in the BART ensemble to be stochastic. In particular, we might
assume that $\pi(T)\propto e^{-C_{T}T}$, for $T\in\mathbb{N}\setminus\\{0\\}$,
with $C_{T}>\log 2$, thus enabling the number of trees in the forest to adapt
to unknown $T_{0}$, as $n,p\rightarrow\infty$.

##  5 Implications

The primary significance of Theorems 4.1, 4.3 and 4.4 is that these results
provide a frequentist theoretical justification for superior empirical
performance of generalized Bayesian trees and forests, claiming that the
posterior concentrates around the truth at a near-optimal learning rate. As
demonstrated below, we can show that the original BART model [3], along with
some of its commonly used variants (such as BART for multi-class
classification and regression on count data) have near-optimal posterior
concentration rates, as direct corollaries of Theorems 4.1 \- 4.4. Another
important consequence of these results is that (see Section A.5 of the
supplementary material), they show that the posterior distribution on the
number of leaves in a generalized Bayesian tree does not exceed the optimal
number of splits by more than a constant multiple and hence are resilient to
overfitting.

Below we demonstrate the breadth of applicability of Theorems 4.1, 4.3 and 4.4
in proving analogous theoretical results for a wide range of commonly used
BART models.

##### Continuous Regression:

For a (multivariate) continuous regression, assume that the response
$\boldsymbol{Y}\,|\>\boldsymbol{X}\sim\mathcal{N}_{p}(\boldsymbol{\mu}(\boldsymbol{X}),\Sigma)$,
for some positive definite $\Sigma$. The function
$g(f_{0}(\boldsymbol{X}))=g(\boldsymbol{\mu})=e^{-\boldsymbol{\mu}^{T}\Sigma^{-1}\boldsymbol{\mu}/2}$
satisfies (10) with $B_{n}=[-n,n]^{p}$ and $C_{g}^{n}=n\lambda(\Sigma)$, where
$\lambda(\Sigma)$ denotes the maximum eigenvalue of $\Sigma$. Hence from
Theorems 4.1, 4.3 and 4.4, we can conclude that for continuous regression, the
G-BART estimator has a near-minimax posterior concentration rate, provided
that the true function $f_{0}$ connecting the input $\boldsymbol{X}$ with the
output $\boldsymbol{Y}$ is either a step function, a monotone function or a
$\nu$-Hölder continuous function with $0<\nu\leq 1$.

##### Classification with Gaussian Step Heights:

For a $p$-class classification the response $\boldsymbol{Y}$ can be written as
a $p$ dimensional binary vector that has $1$ at the $l$-th coordinate if
$\boldsymbol{Y}$ belongs to category $l\in\\{1,\dots,p\\}$ and $0$ elsewhere.
We can assume
$\boldsymbol{Y}\,|\>\boldsymbol{X}\sim\text{Multinomial}(1;\boldsymbol{\pi}(\boldsymbol{X}))$
for some $\boldsymbol{\pi}:\mathbb{R}^{q}\in(0,1)^{p}$ such that
$\boldsymbol{\pi}^{\prime}\textbf{1}_{p}=1$. The unrestricted function
$f_{0}(\boldsymbol{X})$ can be transformed to the natural parameter
$\pi(\boldsymbol{X})$ by a logistic (softmax) or an inverse-probit link
function [3] denoted by $\Psi(\cdot)$, so that
$\pi(\boldsymbol{X})=\Psi(f_{0}(\boldsymbol{X}))$. In either case, the
function $g(f_{0}(\boldsymbol{X}))=1$ trivially satisfies condition (10).
Hence from Theorem 4.1 and Theorem 4.4, we can conclude that the BART model
for multi-class classification has a near-minimax posterior concentration
rate.

##### Classification with Dirichlet Step-Heights

For the same multi-class classification problem with $p$ classes described
above, an alternative prior specification is recommended by [2]. The
parameters $\boldsymbol{\pi}(\boldsymbol{X})$ can be approximated by
multivariate step functions of the form
$f_{\mathcal{T},{P}}(\boldsymbol{x})=\sum_{k=1}^{K}P_{k}\mathbb{I}(\boldsymbol{x}\in\Omega_{k})$
on a tree-partition $\\{\Omega_{k}\\}_{k=1}^{K}$. [2] assumes that
$P_{k}=\left(P_{k1},\dots,P_{kp}\right)\stackrel{{\scriptstyle
i.i.d}}{{\sim}}\text{Dirichlet}(\alpha_{1},\dots,\alpha_{p})$, where
$\alpha_{l}>0,\quad\forall l\in\\{1,\dots,p\\}$. For example, in a binary
classification ($p=2$) problem, we can assign prior
$P_{k}\stackrel{{\scriptstyle i.i.d}}{{\sim}}\text{Beta}(2,2)$ on the step-
heights. The prior $\text{Beta}(2,2)$ violates condition (6). But we can show
that this estimator has a near-optimal posterior concentration rate, even if
we cannot conclude this from the results discussed in Section 4. A proof is
given in the supplementary material. This demonstrates that the assumptions we
make in Section 4 are merely sufficient but not necessary conditions for
proving that the generalized Bayesian tree estimator has a near-minimax
posterior concentration rate.

##### Count Regression:

For count response variable, $\boldsymbol{Y}\sim
Poisson\left[\lambda(\boldsymbol{X})\right]$ with $\lambda(\boldsymbol{X})>0$.
There are several choices for the link function $\Psi(\cdot)$ to map the
unconstrained function $f_{0}(\boldsymbol{X})$ to the constrained parameter
$\lambda(\boldsymbol{X})$. The posterior concentration rate of the Generalized
Bayesian tree estimator might differ depending on which link function is used.
For example, if we use $\Psi(z)=\log\left(1+\exp(z)\right)$, the softplus link
function, then
$g(f_{0}(\boldsymbol{X}))=1/(1+\exp\left(f_{0}(\boldsymbol{X})\right)$,
trivially satisfies condition (10) and we can conclude that the generalized
tree estimator has a near-minimax concentration rate from Theorems 4.1, 4.3
and 4.4.

In contrast, if we use $\Psi(z)=\exp(z)$ as the link function, then
$g(f_{0}(\boldsymbol{X}))=\exp\left(-\exp(f_{0}(\boldsymbol{X}))\right)$ does
not satisfy the condition (10), when the true function $f_{0}$ is a
$\nu$-Hölder continuous function. Therefore we cannot apply Theorem 4.4
anymore to imply that the generalized tree estimator has a near-optimal rate
of posterior concentration. When $f_{0}$ is a step function with complexity
$K_{f_{0}}$, the condition (10) is satisfied with $B_{n}=[-K_{f_{0}}\log
n,K_{f_{0}}\log n]$ and $C_{g}^{n}=n^{K_{f_{0}}}$. The posterior concentration
rate becomes
$\varepsilon_{n}=n^{-\frac{1-\alpha}{2}}\sqrt{K_{f_{0}}\log^{2\eta}(n/K_{f_{0}})}$
under the assumption $K_{f_{0}}\lesssim n^{\alpha}$ for some $0<\alpha<1$.
This is slower than the near-optimal concentration rate
$n^{-\frac{1}{2}}\sqrt{K_{f_{0}}\log^{2\eta}(n/K_{f_{0}})}$, if we use
$\Psi(z)=\log\left(1+\exp(z)\right)$, the softplus link function, instead.
This demonstrates the need for choosing suitable link functions in empirical
applications.

##  6 Discussion

In this paper we have examined a general framework for Bayesian Additive
Regression Tree Models that encapsulates various conventional BART models
adapted to a wide range of regression and classification tasks. We
demonstrated that these models have a near-minimax posterior concentration
rate for a wide range of functions, thus corroborating the empirical success
of BART and its variants, from a theoretical perspective. These results also
build the foundation for uncertainty quantification statements for a wide
variety of BART models, opening up an interesting avenue for future research.
Among empirical implications, we have established the need for careful
modeling choices such as selecting appropriate link functions. The theoretical
results also substantiate the scope of a wider variety of distributions on
approximating step-heights, that can prove advantageous for applications where
the response distribution has a thicker tail. These theoretical findings also
provide strong motivation for exploring novel application areas for flexible
BART-like models.

## References

  * [1] 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.
  * [2] David GT Denison, Bani K Mallick, and Adrian FM Smith.  A Bayesian CART algorithm.  Biometrika, 85(2):363–377, 1998.
  * [3] Hugh A Chipman, Edward I George, Robert E McCulloch, et al.  Bart: Bayesian additive regression trees.  The Annals of Applied Statistics, 4(1):266–298, 2010.
  * [4] Rodney Sparapani, Charles Spanbauer, and Robert McCulloch.  The BART R package, 2019.
  * [5] Justin Bleich, Adam Kapelner, Edward I George, and Shane T Jensen.  Variable selection for bart: an application to gene regulation.  The Annals of Applied Statistics, pages 1750–1781, 2014.
  * [6] 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.
  * [7] Jingyu He, Saar Yalov, and P Richard Hahn.  XBART: Accelerated Bayesian additive regression trees.  In The 22nd International Conference on Artificial Intelligence and Statistics, pages 1130–1138, 2019.
  * [8] Jennifer L Hill.  Bayesian nonparametric modeling for causal inference.  Journal of Computational and Graphical Statistics, 20(1):217–240, 2011.
  * [9] Jennifer Hill and Yu-Sung Su.  Assessing lack of common support in causal inference using Bayesian nonparametrics: Implications for evaluating the effect of breastfeeding on children’s cognitive outcomes.  The Annals of Applied Statistics, pages 1386–1420, 2013.
  * [10] P Richard Hahn, Jared S Murray, and Carlos Carvalho.  Bayesian regression tree models for causal inference: regularization, confounding, and heterogeneous effects.  arXiv preprint arXiv:1706.09523, 2017.
  * [11] Junliang Du and Antonio R Linero.  Interaction detection with Bayesian decision tree ensembles.  In The 22nd International Conference on Artificial Intelligence and Statistics, pages 108–117. PMLR, 2019.
  * [12] Rodney A Sparapani, Brent R Logan, Robert E McCulloch, and Purushottam W Laud.  Nonparametric survival analysis using Bayesian additive regression trees (BART).  Statistics in medicine, 35(16):2741–2753, 2016.
  * [13] Matthew A Taddy, Robert B Gramacy, and Nicholas G Polson.  Dynamic trees for learning and design.  Journal of the American Statistical Association, 106(493):109–123, 2011.
  * [14] Sameer K Deshpande, Ray Bai, Cecilia Balocchi, and Jennifer E Starling.  VC-BART: Bayesian trees for varying coefficients.  arXiv preprint arXiv:2003.06416, 2020.
  * [15] Antonio R Linero.  Bayesian regression trees for high-dimensional prediction and variable selection.  Journal of the American Statistical Association, 113(522):626–636, 2018.
  * [16] Yi Liu, Veronika Rockova, and Yuexi Wang.  Abc variable selection with Bayesian forests.  arXiv preprint arXiv:1806.02304, 2018.
  * [17] Yi Liu and Veronika Rockova.  Variable selection via thompson sampling.  arXiv preprint arXiv:2007.00187, 2020.
  * [18] Jared S Murray.  Log-linear bayesian additive regression trees for multinomial logistic and count regression models.  Journal of the American Statistical Association, 116(534):756–769, 2021.
  * [19] MT Pratola, HA Chipman, EI George, and RE McCulloch.  Heteroscedastic BART via multiplicative regression trees.  Journal of Computational and Graphical Statistics, pages 1–13, 2019\.
  * [20] Hugh A Chipman, Edward I George, Robert E McCulloch, and Thomas S Shively.  High-dimensional nonparametric monotone function estimation using bart.  arXiv preprint arXiv:1612.01619, 2016.
  * [21] Jennifer E Starling, Jared S Murray, Carlos M Carvalho, Radek K Bukowski, James G Scott, et al.  Bart with targeted smoothing: An analysis of patient-specific stillbirth risk.  Annals of Applied Statistics, 14(1):28–50, 2020.
  * [22] Veronika Rockova and Enakshi Saha.  On theory for BART.  In The 22nd International Conference on Artificial Intelligence and Statistics, pages 2839–2848, 2019.
  * [23] Veronika Rockova, Stéphanie van der Pas, et al.  Posterior concentration for Bayesian regression trees and forests.  Annals of Statistics, 48(4):2108–2131, 2020.
  * [24] Veronika Rockova.  On semi-parametric inference for BART.  In International Conference on Machine Learning, pages 8137–8146. PMLR, 2020.
  * [25] Antonio Ricardo Linero and Yun Yang.  Bayesian regression tree ensembles that adapt to smoothness and sparsity.  arXiv preprint arXiv:1707.09461, 2017.
  * [26] Ismael Castillo and Veronika Ro$v$cková.  Uncertainty quantification for bayesian cart.  The Annals of Statistics, 49(6):3482–3509, 2021.
  * [27] Seonghyun Jeong and Veronika Rockova.  The art of BART: On flexibility of Bayesian forests.  arXiv preprint arXiv:2008.06620, 2020.
  * [28] John Ashworth Nelder and Robert WM Wedderburn.  Generalized linear models.  Journal of the Royal Statistical Society: Series A (General), 135(3):370–384, 1972.
  * [29] Elias M Stein and Rami Shakarchi.  Real analysis: measure theory, integration, and Hilbert spaces.  Princeton University Press, 2009.
  * [30] Dennis D Cox.  An analysis of bayesian inference for nonparametric regression.  The Annals of Statistics, pages 903–923, 1993.
  * [31] Persi Diaconis and David Freedman.  On the consistency of bayes estimates.  The Annals of Statistics, pages 1–26, 1986.
  * [32] Subhashis Ghosal, Aad Van Der Vaart, et al.  Convergence rates of posterior distributions for noniid observations.  The Annals of Statistics, 35(1):192–223, 2007.
  * [33] Ismaël Castillo, Richard Nickl, et al.  On the bernstein–von mises phenomenon for nonparametric bayes procedures.  The Annals of Statistics, 42(5):1941–1969, 2014.
  * [34] Chao Gao, Fang Han, and Cun-Hui Zhang.  Minimax risk bounds for piecewise constant models.  arXiv preprint arXiv:1705.06386, 2017.
  * [35] Gérard Biau and Luc Devroye.  On the risk of estimates for block decreasing densities.  Journal of multivariate analysis, 86(1):143–165, 2003.
  * [36] Charles J Stone.  Optimal global rates of convergence for nonparametric regression.  The annals of statistics, pages 1040–1053, 1982.

##  Appendix A Appendix

###  A.1 Bayesian CART Prior by Denison et al. (1998)

We describe the Bayesian CART prior proposed by [2]. The prior on individual
Bayesian trees is assigned conditional on the number of terminal nodes/ leaves
$K$ and all prior probability is concentrated on the set of all valid tree
partitions, as defined below (Definition 3.1 of [23]):

######  Definition A.1.

Denote by $\mbox{\boldmath$\Omega$}=\\{\Omega\\}_{k=1}^{K}$, a partition of
$[0,1]^{p}$, We say that $\Omega$ is valid if

$\mu(\Omega_{k})\geq\frac{C}{n}\quad\forall k=1,\dots,K$ (12)

for some $C\in\mathbb{N}\setminus\\{0\\}$.

Valid partitions have non-empty cells, where the allocation does not need to
be balanced. Now the prior on tree partitions is specified as follows:

  1. 1.

The number of leaves in a tree $K$ follows a Poisson distribution with
parameter $\lambda>0$

$P(K)=\frac{\lambda^{K}}{(e^{\lambda}-1)K!},\qquad k=1,2,\dots$ (13)

  2. 2.

Given the number of leaves $K$, a tree is chosen uniformly at random from the
set of all available valid tree-partitions with $K$ leaves. Number of valid
tree partitions is given by

$\Delta(V_{K})=\frac{q^{K-1}n!}{(n-K+1)!}$ (14)

This is a slightly modified version of the original prior proposed by [2].
This modified version was used by [23] to derive posterior concentration rates
for the BART estimator under this prior.

  3. 3.

At each node, the splitting rule consists of picking a split variable $j$
uniformly at random from the available directions $\\{1,\dots,q\\}$ and
picking a split point $c$, also uniformly at random from the available data
values $x_{1j},\dots,x_{nj}$.

###  A.2 Preliminary Results with Proof

######  Lemma A.2.

The multivariate Gaussian $\mathcal{N}_{p}({\bf 0},\mathbb{I}_{p})$ and the
multivariate Laplace $\mathcal{L}_{p}({\bf 0},\mathbb{I}_{p})$ distribution
belong to the general family of distributions with CDF $F_{\beta}$ that has
the following property: For some $C_{1}>0$, $0<C_{2}\leq 2$ and $C_{3}>0$ and
any $t>0$,

$\displaystyle F_{\beta}(\left\|\beta\right\|_{\infty}\leq
t)\gtrsim\left(e^{-C_{1}t^{C_{2}}}t\right)^{p}\quad\text{for $t>0$}$ (15)
$\displaystyle F_{\beta}(\left\|\beta\right\|_{\infty}\geq t)\lesssim
e^{-C_{3}t}\quad\text{for $t\geq 1$}$ (16)

###### Proof.

If $F_{\beta}=\mathcal{N}_{p}({\bf 0},\mathbb{I}_{p})$, then for any $t>0$,

$F_{\beta}(\left\|\beta\right\|_{\infty}\leq
t)\gtrsim\left(e^{-t^{2}/2}\int_{-t}^{t}d\beta\right)^{p}\gtrsim
e^{-pt^{2}/2}t^{p}$

For $t\geq 1$

$F_{\beta}(\left\|\beta\right\|_{\infty}\geq
t)\lesssim\left(e^{-t^{2}/4}2\int_{t}^{\infty}e^{-z^{2}/4}dz\right)^{p}\lesssim
e^{-C_{3}t}$

If $F_{\beta}=\mathcal{L}_{p}({\bf 0},\mathbb{I}_{p})$, then for any $t>0$,

$F_{\beta}(\left\|\beta\right\|_{\infty}\leq
t)\gtrsim\left(e^{-t}\int_{-t}^{t}d\beta\right)^{p}\gtrsim e^{-pt}t^{p}$

Also, for any $t\geq 0$,

$F_{\beta}(\left\|\beta\right\|_{\infty}\geq t)=\frac{e^{-pt}}{2}<e^{-pt}$

∎

######  Lemma A.3.

Let $f$ and $f_{0}$ denote step functions of the form
$f(\boldsymbol{X})=\sum_{k=1}^{K}\beta_{k}\mathbb{I}(\boldsymbol{X}\in\Omega_{k})$
and
$f_{0}(\boldsymbol{X})=\sum_{k=1}^{K}\beta^{0}_{k}\mathbb{I}(\boldsymbol{X}\in\Omega_{k})$
respectively, on a tree-shaped partition $\\{\Omega_{k}\\}_{k=1}^{K}$. Let
$P_{f}$ and $P_{f_{0}}$ denote two probability densities belonging to an
Exponential family distribution of the form

$P_{f}(\boldsymbol{Y}\,|\>\boldsymbol{X})=h(\boldsymbol{Y})g\left[f(\boldsymbol{X})\right]\exp\left[\eta\left(f(\boldsymbol{X})\right)^{T}T(\boldsymbol{Y})\right],$
(17)

with parameters $f$ and $f_{0}$ respectively. If
$\left|\frac{\nabla^{T}g(\boldsymbol{\beta})}{g(\boldsymbol{\beta})}\right|\leq
C^{n}_{g}\textbf{1}_{p}$, for some positive sequence $\\{C^{n}_{g}\\}_{n\geq
1}$, then

$\displaystyle K_{n}(P_{f},P_{f_{0}})\vee V_{n}(P_{f},P_{f_{0}})\lesssim
C^{n}_{g}\sum_{k=1}^{K}\left\|\beta_{k}-\beta^{0}_{k}\right\|_{1}$ (18)
$\displaystyle H_{n}(P_{f},P_{f_{0}})\lesssim
C^{n}_{g}\sum_{k=1}^{K}\left\|\beta_{k}-\beta^{0}_{k}\right\|_{1}$ (19)

###### Proof.

Denoting $f_{i}=f(\boldsymbol{X}_{i})$ and $f_{i0}=f_{0}(\boldsymbol{X}_{i})$,
we can write

$\displaystyle K_{n}(P_{f},P_{f_{0}})$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}g(f_{i})\int
h(\boldsymbol{Y})\exp\left(f_{i}T(\boldsymbol{Y})\right)[\log\frac{g(f_{i})}{g(f_{i0})}$
$\displaystyle+\exp\left[(f_{i}-f_{i0})^{T}T(\boldsymbol{Y})\right]]d\boldsymbol{Y}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\left[\log\frac{g(f_{i})}{g(f_{i0})}+(f_{i}-f_{i0})^{T}\mathbb{E}\left[T(\boldsymbol{Y})\right]\right]$
$\displaystyle=\sum_{k=1}^{K}\mu(\Omega_{k})\left[\log\frac{g(\beta_{k})}{g(\beta_{k}^{0})}-\frac{\nabla^{T}g(\beta_{k})}{g(\beta_{k})}\left(\beta_{k}-\beta^{0}_{k}\right)\right]$

By triangle inequality and Taylor series approximation of $\log\nabla
g(\beta_{k})$ about $\beta^{0}_{k}$, we get

$\displaystyle K_{n}(P_{f},P_{f_{0}})$
$\displaystyle\lesssim\sup\left|\frac{\nabla^{T}g(\cdot)}{g(\cdot)}\right|\sum_{k=1}^{K}\left\|\beta_{k}-\beta_{k}^{0}\right\|_{1}$
$\displaystyle=C^{n}_{g}\sum_{k=1}^{K}\left\|\beta_{k}-\beta_{k}^{0}\right\|_{1},$

Similar technique works for $V_{n}(P_{f},P_{f_{0}})$

Also, Since Hellinger metric is bounded from above by Kullback-Leibler
divergence, $H_{n}$ satisfies,

$H_{n}(P_{f},P_{f_{0}})\lesssim
C^{n}_{g}\sum_{k=1}^{K}\left\|\beta_{k}-\beta^{0}_{k}\right\|_{1}$

∎

######  Lemma A.4.

Any bounded monotone function $f_{0}$ can be approximated with arbitrary
precision $\varepsilon_{n}$, by a step function supported on a k-d tree
partition with $\smash{\widehat{K}}\geq\lceil 1/\epsilon_{n}\rceil$ leaves.

###### Proof.

Without loss of generality, assume $0\leq f_{0}(\cdot)\leq 1$. Partition
interval $[0,1]$ by $0=y_{0}<y_{1}<\dots<y_{k}<\dots<y_{K-1}<y_{K}=1$, with
$K=\lceil 1/\epsilon_{n}\rceil$. Then
$\,|\>y_{k}-y_{k-1}\,|\><\varepsilon_{n}$ and we can approximate
$f_{0}(\boldsymbol{X})$ by the step function:

$f(\boldsymbol{X})=\sum_{k=1}^{K}y_{k}\mathbb{I}\\{\boldsymbol{X}\in\Omega_{k}\\}$

, where $\Omega_{k}=f^{-1}[y_{k-1},y_{k}]$.

If $f$ is monotone, $\Omega_{k}=\prod_{j=1}^{q}\\{x_{j}\in\mathrm{I}_{j}\\},$
where $\mathrm{I}_{j}$ is an interval and $x_{j}$ denotes the $j$-th
coordinate of $\boldsymbol{X}\in\mathbb{R}^{q}$.

Since any step function supported on an axis-paralleled partition has an
equivalent step function supported on a k-d tree, we can approximate the axis
paralleled partition $\\{\Omega_{k}\\}_{k=1}^{K}$ by a recursive binary tree
partition $\\{\smash{\widehat{\Omega}}_{k}\\}_{k=1}^{\smash{\widehat{K}}}$
with number of leaves $\smash{\widehat{K}}\geq K$. ∎

###  A.3 Proof of Main Results

In this section we prove Theorem 4.1 and Theorem 4.3. Most steps in the proofs
are identical and hence for simplicity we describe the common steps of the
proofs together and mark the steps that are different by the corresponding
theorem number. We need to prove three conditions: entropy condition (C1),
prior concentration condition (C2) and prior decay rate condition (C3). The
steps of the proofs for each of these conditions are described below.

####  A.3.1 Entropy Condition

Define

$\mathcal{F}_{n}=\\{f_{\mathcal{T},{\boldsymbol{\beta}}}(\boldsymbol{X})\text{
of the form \eqref{eq:tree_mapping} with }K=k_{n}\text{ and
}\left\|\beta\right\|_{\infty}\leq C^{n}_{\beta}\\},$

where $k_{n}\propto n\varepsilon_{n}^{2}/\log n$ and $C^{n}_{\beta}$ is
defined in Assumption 1.

Since $\left\|\boldsymbol{z}\right\|_{1}\leq
Kp\left\|\boldsymbol{z}\right\|_{\infty}$ for any
$\boldsymbol{z}\in\mathbb{R}^{Kp}$, by the bound (19) and by definition of
$\mathcal{F}_{n}$, we can write

$\displaystyle N\left(\frac{\varepsilon_{n}}{36},\mathcal{F}_{n},H_{n}\right)$
$\displaystyle\lesssim\sum_{K=1}^{k_{n}}N\left(\frac{\varepsilon_{n}}{36C^{n}_{g}Kp},\\{\beta:\left\|\beta\right\|_{\infty}\leq
C^{n}_{\beta}\\},\left\|\cdot\right\|_{\infty}\right)$
$\displaystyle\lesssim\sum_{K=1}^{k_{n}}\left(\frac{36C_{\beta}^{n}C_{g}^{n}Kq}{\varepsilon_{n}}\right)^{Kq}$

Therefore the LHS of (C1) can be bounded from above by

$\displaystyle(k_{n}+1)p\left[\log 36+\log(C_{\beta}^{n}C_{g}^{n})+\log
k_{n}+\log p-\log\varepsilon_{n}\right]$

Since $C_{\beta}^{n}C_{g}^{n}\lesssim n^{M}$ for some $M>0$, ignoring smaller
terms, proving condition (C1) reduces to proving

$\displaystyle(k_{n}+1)p\log n\lesssim n\varepsilon_{n}^{2}$ (20)

##### Theorem 4.1:

When $f_{0}$ is a step function with complexity $K_{f_{0}}$ we can prove (20)
by replacing
$\varepsilon_{n}=n^{-1/2}\sqrt{K_{f_{0}}\log^{2\eta}(n/K_{f_{0}})}$ and
$k_{n}\propto\frac{n\varepsilon_{n}^{2}}{p\log(n/K_{f_{0}})}=K_{f_{0}}\log^{2\theta-1}(n/K_{f_{0}})$
for some $\theta>1/2$.

##### Theorem 4.3:

When $f_{0l}$ is a $\nu$-Hölder continuous function with $0<\nu\leq 1$ for all
$l=1,\dots,p$, replacing $\varepsilon_{n}=n^{-\nu/(2\nu+q)}\sqrt{\log n}$ and
$k_{n}\propto\frac{n\varepsilon_{n}^{2}}{\log n}=n^{q/(2\nu+q)}$ proves (20).

####  A.3.2 Prior Concentration Condition

Let
$\widetilde{f}_{0}=\left(f_{\mathcal{T},{\boldsymbol{B}_{1}^{0}}}(\boldsymbol{x}),\dots,f_{\mathcal{T},{\boldsymbol{B}_{q-1}^{0}}}(\boldsymbol{x})\right)$
denote the projection of $f_{0}$ onto a balanced k-d tree partition with
$a_{n}$ leaves, where $a_{n}$ is chosen so that
$\left\|f_{0}-\widetilde{f}_{0}\right\|_{2,n}<\varepsilon_{n}/2$.

##### Theorem 4.1:

If $f_{0}$ is a step function, $a_{n}=K_{f_{0}}$

##### Theorem 4.3:

If $f_{0}$ is a $\nu$-Hölder continuous function, $a_{n}$ is chosen by the
following lemma, which is analogous to Lemma 3.2 of [23].

######  Lemma A.5.

Denote $f=\\{f_{l}\\}_{l=1}^{p}$ and assume $f_{l}\in\mathcal{H}^{\nu_{l}}$
where $\nu_{l}\leq 1$ for all $l=1,\dots,p$ and $\mathcal{X}$ is regular. Then
there exists tree structured step functions
$\hat{f}=\\{f_{\mathcal{T},\boldsymbol{B}_{l}}\\}_{l=1}^{p}\in\mathcal{F}_{K}$
for some given tree partition $\mathcal{T}$ with $K\in\mathbb{N}$ leaves such
that for some constant $C>0$,

$\left\|\hat{f}-f\right\|_{2,n}\leq
Cd\sum_{l=1}^{p}\left(\frac{1}{K^{\nu_{l}/q}}\left\|f_{l}\right\|_{\mathcal{H}^{\nu_{l}}}\right)\leq
C\frac{q}{K^{\nu/q}}\sum_{l=1}^{p}\left(\left\|f_{l}\right\|_{\mathcal{H}^{\nu_{l}}}\right),$

where $\nu=\min_{l=1}^{p}\nu_{l}$.

As a corollary, replacing
$C_{0}=C\left(\sum_{l=1}^{p}\left\|f_{l}\right\|_{\mathcal{H}^{\nu}}\right)$,
$a_{n}$ satisfies

$\left(\frac{2C_{0}q}{\varepsilon_{n}}\right)^{q/\nu}\leq
a_{n}\leq\left(\frac{2C_{0}q}{\varepsilon_{n}}\right)^{q/\nu}+1$ (21)

Using (18) and by triangle inequality, we can bound the LHS of (C2) from below
by

$\displaystyle C\pi(a_{n})\Pi\left(\beta\in
B_{n}^{a_{n}}:\left\|\beta-\beta^{0}\right\|_{1}\leq\frac{\epsilon_{n}^{2}}{2C_{g}^{n}}\right)$

For the prior by [3], $C=1$ and $\pi(a_{n})\gtrsim e^{-a_{n}\log a_{n}}$ (by
Corollary 5.2 of [22]).

For the prior by [2],
$C=\frac{1}{\,|\>F_{a_{n}}\,|\>}>(a_{n}dn)^{-a_{n}}>e^{-a_{n}\log a_{n}}$ (by
Lemma 3.1 of [23]) and $\pi(a_{n})\gtrsim e^{-a_{n}\log a_{n}}$ (by proof of
Theorem 4.1 of [23]).

Thus for both priors $C\pi(a_{n})\gtrsim e^{-2a_{n}\log a_{n}}$.

Next we bound $\Pi\left(\beta\in
B_{n}^{a_{n}}:\left\|\boldsymbol{\beta}-\boldsymbol{\beta}^{0}\right\|_{1}\leq\frac{\epsilon_{n}^{2}}{2C_{g}^{n}}\right)$,
up to a constant, from below by

$\displaystyle\Pi\left(\boldsymbol{\beta}:\left\|\boldsymbol{\beta}\right\|_{\infty}\leq
C_{\beta}^{n},\quad\left\|\boldsymbol{\beta}-\boldsymbol{\beta}^{0}\right\|_{\infty}\leq\frac{\epsilon_{n}^{2}}{2a_{n}qC_{g}^{n}}\right)$

Since $C^{n}_{g}$ and $C^{n}_{\beta}$ both are increasing with $n$, for
sufficiently large $n$, the above expression is bounded below by

$\displaystyle\Pi\left(\boldsymbol{\beta}:\left\|\boldsymbol{\beta}-\boldsymbol{\beta}^{0}\right\|_{\infty}\leq\frac{\varepsilon_{n}^{2}}{2a_{n}pC_{g}^{n}}\right)$
$\displaystyle\gtrsim$ $\displaystyle
e^{-C_{1}a_{n}p\left(\left\|\beta_{0}\right\|_{\infty}+\frac{\varepsilon_{n}^{2}}{2a_{n}pC_{g}^{n}}\right)^{C}_{2}}\left({\left\|\beta_{0}\right\|_{\infty}+\frac{\varepsilon_{n}^{2}}{2a_{n}pC_{g}^{n}}}\right)^{a_{n}p}$

Since $\varepsilon_{n}^{2}\rightarrow 0$ and both $a_{n}$ and $C_{g}^{n}$ are
both increasing with $n$, assuming
$\left\|f_{0}\right\|_{\infty}\lesssim\sqrt{\log n}$, the above bound reduces
to

$\displaystyle
e^{-C_{1}a_{n}p\log^{C_{2}/2}n}\left\|\beta_{0}\right\|_{\infty}^{a_{n}p/2}\gtrsim\log\left[-C_{1}a_{n}p\log^{C_{2}/2}n\right]$

We can prove $e^{-a_{n}\log n}\gtrsim e^{-n\varepsilon_{n}^{2}}$ for Theorem
4.1 and Theorem 4.3 separately by replacing appropriate values of
$\varepsilon_{n}$. Since $C_{2}\leq 2$, this would complete the proof.

####  A.3.3 Prior Decay Rate Condition

##### Theorem 4.1:

When $f_{0}$ is a step-function with complexity $K_{f_{0}}$,

$\displaystyle\Pi(\mathcal{F}\setminus\mathcal{F}_{n})$
$\displaystyle\leq\Pi(\mathcal{F}\setminus\bigcup_{K=1}^{k_{n}}F_{K})+\Pi(\bigcup_{K\leq
k_{n}}\\{f\in F_{K}:\left\|\beta\right\|_{\infty}>C_{\beta}^{n}\\})$
$\displaystyle\leq\Pi(\bigcup_{K>k_{n}}F_{K})+e^{-K_{f_{0}}\log n/2}$
$\displaystyle=\Pi(\bigcup_{K>k_{n}}F_{K})+o(e^{-n\varepsilon_{n}^{2}})$

The last line is due to the fact $C_{\beta}^{n}\gtrsim K_{f_{0}}\log n$ when
$f_{0}$ is a step-function with complexity $K_{f_{0}}$.

##### Theorem 4.3:

When $f_{0}$ is a $\nu$-Hölder continuous function, the LHS of condition (C3)
can be bounded from above by

$\displaystyle\Pi(\mathcal{F}\setminus\mathcal{F}_{n})$
$\displaystyle\leq\Pi(\mathcal{F}\setminus\bigcup_{K=1}^{k_{n}}F_{K})+\Pi(\bigcup_{K\leq
k_{n}}\\{f\in F_{K}:\left\|\beta\right\|_{\infty}>C_{\beta}^{n}\\})$
$\displaystyle\leq\Pi(\bigcup_{K>k_{n}}F_{K})+\sum_{K=1}^{k_{n}}\Pi(\\{\beta:\left\|\beta\right\|_{\infty}>C_{\beta}^{n}\\})$
$\displaystyle\leq\Pi(\bigcup_{K>k_{n}}F_{K})+\sum_{K=1}^{k_{n}}e^{-C^{n}_{\beta}},\quad\text{by
condition \eqref{eq:beta_tail2}}$
$\displaystyle\leq\Pi(\bigcup_{K>k_{n}}F_{K})+k_{n}e^{-C^{n}_{\beta}}$
$\displaystyle=\Pi(\bigcup_{K>k_{n}}F_{K})+o(e^{-n\varepsilon_{n}^{2}})$

The last line is due to the fact $C_{\beta}^{n}\gtrsim n$, when $f_{0}$ is a
$\nu$-Hölder continuous functions.

Therefore it is enough to show that

$\Pi(\bigcup_{K>k_{n}}F_{K})\lesssim e^{-n\varepsilon_{n}^{2}}$

This condition is satisfied for both priors under consideration. This follows
from section 8.3 of [23] for the prior by [2] and from Corollary 5.2 of [22]
for the prior by [3].

###  A.4 Classification with Dirichlet Step Heights

For a multi-class classification problem with $p$ classes, where the response
variable $\boldsymbol{Y}$ is a categorical random variable with $p$
categories, $\boldsymbol{Y}$ can be written as a $p$ dimensional binary vector
that has $1$ at the $l$-th coordinate if $\boldsymbol{Y}$ belongs to category
$l\in\\{1,\dots,p\\}$ and $0$ elsewhere. G-BART assumes

$\boldsymbol{Y}\,|\>\boldsymbol{X}\sim\text{Multinomial}(1,\boldsymbol{f_{0}}(\boldsymbol{X})),$
(22)

where
$\boldsymbol{f_{0}}=\left(f_{01},\dots,f_{0p}\right)^{\prime}:\mathbb{R}^{q}\rightarrow(0,1)^{p}$
is a constrained function with
$\boldsymbol{f_{0}}(\boldsymbol{X})^{\prime}\textbf{1}_{p}=1$ for any
$\boldsymbol{X}\in\mathbb{R}^{q}$. Each $f_{0l}(\cdot)$ can be approximated by
a step function of the form

$f_{\mathcal{T},{P}}(\boldsymbol{x})=\sum_{k=1}^{K}P_{k}\mathbb{I}(\boldsymbol{x}\in\Omega_{k})$
(23)

on a tree-partition $\\{\Omega_{k}\\}_{k=1}^{K}$. [2] assumes

$P_{k}=\left(P_{k1},\dots,P_{kp}\right)\stackrel{{\scriptstyle
i.i.d}}{{\sim}}\text{Dirichlet}(\alpha_{1},\dots,\alpha_{p}),$ (24)

where $\alpha_{l}>0,\quad\forall l\in\\{1,\dots,p\\}$.

######  Theorem A.6.

If we assume that the distribution of the step-sizes satisfies (24), then
under Assumptions 1 & 2 described in section 4 of the manuscript, the Bayesian
Tree estimator satisfies the following property,:

(i) If $f_{0}$ is $\nu$-Hölder continuous with $0<\nu\leq 1$ where
$\|f_{0}\|_{\infty}\lesssim\log^{1/2}n$, then with
$\varepsilon_{n}=n^{-\alpha/(2\alpha+p)}\log^{1/2}n$, and

(ii) If $f_{0}$ is step-function with complexity $K_{f_{0}}\lesssim\sqrt{n}$,
then with
$\varepsilon_{n}=n^{-1/2}\sqrt{K_{f_{0}}p\log^{2\nu}{\left(n/K_{f_{0}}p\right)}}n$,

$\Pi\left(f\in\mathcal{F}:H_{n}(\mathbb{P}_{f},\mathbb{P}_{f_{0}})>M_{n}\,\varepsilon_{n}\mid\boldsymbol{Y}^{(n)}\right)\to
0,$

for any $M_{n}\to\infty$ in $\mathbb{P}_{f_{0}}^{(n)}$-probability, as
$n,p\to\infty$.

The above statement is true for both tree priors considered in this paper: the
prior by [2] and a modified version of the prior by [1] with
$p_{split}(\Omega_{t})=\alpha^{d(\Omega_{t})}$ for some $1/n\leq\alpha<1/2$.

###### Proof.

We need to prove three conditions: entropy condition (C1), prior concentration
condition (C2) and prior decay rate condition (C3). Among these (C1) and (C3)
can be proved by the same technique as in section A.3. Therefore we will only
prove Condition (C3) here. We need to show, for some $c>0$

$\Pi\left(f\in\mathcal{F}:max\\{K_{n}(f,f_{0}),V_{n}(f,f_{0})\\}\leq\varepsilon_{n}^{2}\right)\gtrsim
e^{-cn\varepsilon_{n}^{2}}$ (25)

Let
$\widetilde{f}_{0}=\left(f_{\mathcal{T},{P_{1}^{0}}}(\boldsymbol{x}),\dots,f_{\mathcal{T},{P_{q}^{0}}}(\boldsymbol{x})\right)$
denote the projection of $f_{0}$ onto a balanced k-d tree partition
$\mathcal{T}$ with $a_{n}$ leaves, where $a_{n}$ is chosen so that
$\left\|f_{0}-\widetilde{f}_{0}\right\|_{2,n}<\varepsilon_{n}/2$. If $f_{0}$
is a step function, $a_{n}=K_{f_{0}}$. If $f_{0}$ is a $\nu$-Hölder continuous
function, $a_{n}$ is chosen by Lemma 3.2 of [23], where replacing
$C_{0}=C\left(\sum_{l=1}^{p}\left\|f_{l}\right\|_{\mathcal{H}^{\nu}}\right)$
we get

$\left(\frac{2C_{0}q}{\varepsilon_{n}}\right)^{q/\nu}\leq
a_{n}\leq\left(\frac{2C_{0}q}{\varepsilon_{n}}\right)^{q/\nu}+1$ (26)

$f_{\mathcal{T},{P^{0}_{l}}}(\boldsymbol{x})$ is of the form (23) for some
tree topology $\mathcal{T}$ with $a_{n}$ leaves and
$P^{0}_{l}=\\{P_{kl}^{0}\\}_{k=1}^{a_{n}}$ for $l=1,\dots,p$. We assume there
exists some $\delta_{0}>0$ such that $\min{f_{0l}}>\delta_{0}$ for all
$l=1,\dots,q$. This implies $P_{lk}^{0}>\delta_{0}$ for all $l=1,\dots q$ and
all $k=1,\dots,K$. Therefore by (18), we can bound the LHS of (25) from above
by

$\displaystyle
C\pi(a_{n})\Pi\left(P\in[0,1]^{a_{n}p}:\left\|P-P^{0}\right\|_{1}\leq\delta_{0}\varepsilon_{n}^{2}/2\right)$

For the prior by [1], $C=1$ and for the prior by [2],
$C=\frac{1}{\,|\>F_{a_{n}}\,|\>}>(a_{n}dn)^{-a_{n}}>e^{-a_{n}\log a_{n}}$ (by
Lemma 3.1 of [23]). By Corollary 5.2 of [22] for the prior by [1] and by proof
of Theorem 4.1 of [23] for the prior by [2], we can show $\pi(a_{n})\geq
e^{-a_{n}\log a_{n}}$. Thus for both priors,

$C\pi(a_{n})>e^{-2a_{n}\log a_{n}}$ (27)

Since $P_{k}\sim\text{Dirichlet}(\alpha_{1},\dots,\alpha_{p})$ for all
$k=1,\dots,K$ and $P_{lk}^{0}>\delta_{0}$, for all $l=1,\dots,p$ and all
$k=1,\dots,K$, we can bound
$\Pi\left(P\in[0,1]^{a_{n}q}:\left\|P-P^{0}\right\|_{1}\leq\delta_{0}\epsilon_{n}^{2}/2\right)$
from above by

$\displaystyle\Pi\left(P\in[0,1]^{a_{n}p}:\left\|P-P^{0}\right\|_{\infty}\leq\frac{\delta_{0}\epsilon_{n}^{2}}{2a_{n}p}\right)\gtrsim
C_{\alpha}(\frac{\delta_{0}\epsilon_{n}^{2}}{a_{n}p})^{a_{n}p},$ (28)

where $C_{\alpha}$ is a constant that depends on the Dirichlet parameters
$\alpha=(\alpha_{1},\dots,\alpha_{q})$. Combining (27) and (28) completes the
proof. ∎

###  A.5 Parsimony of G-BART

As a by-product of the theoretical results discussed in section 4, we can the
following statements which support the empirical observation that generalized
Bayesian trees are resilient to overfitting.

  * (i)

Under the assumptions of Theorem 4.1 we have $\Pi\left(K\gtrsim
K_{f_{0}}\mid\boldsymbol{Y}^{(n)}\right)\to 0$ in
$\mathbb{P}_{f_{0}}^{(n)}$-probability, as $n,q\to\infty$.

  * (ii)

Under the assumptions of Theorem 4.3 we have $\Pi\left(K\gtrsim
n^{q/(2+q)}\mid\boldsymbol{Y}^{(n)}\right)\to 0$ in
$\mathbb{P}_{f_{0}}^{(n)}$-probability, as $n,q\to\infty$.

  * (iii)

Under the assumptions of Theorem 4.4 we have $\Pi\left(K\gtrsim
n^{q/(2\nu+q)}\mid\boldsymbol{Y}^{(n)}\right)\to 0$ in
$\mathbb{P}_{f_{0}}^{(n)}$-probability, as $n,q\to\infty$.

###### Proof.

The proofs of (i), (ii) and (iii) follow from Lemma 1 of [32], in conjunction
with the proofs of Theorems 4.1, 4.3 and 4.4 respectively. ∎