Recursive partitioning and multi-scale modeling on conditional densities

We introduce a nonparametric prior on the conditional distribution of a (univariate or multivariate) response given a set of predictors. The prior is constructed in the form of a two-stage generative procedure, which in the first stage recursively partitions the predictor space, and then in the second stage generates the conditional distribution by a multi-scale nonparametric density model on each predictor partition block generated in the first stage. This design allows adaptive smoothing on both the predictor space and the response space, and it results in the full posterior conjugacy of the model, allowing exact Bayesian inference to be completed analytically through a forward-backward recursive algorithm without the need of MCMC, and thus enjoying high computational efficiency (scaling linearly with the sample size). We show that this prior enjoys desirable theoretical properties such as full $L_1$ support and posterior consistency. We illustrate how to apply the model to a variety of inference problems such as conditional density estimation as well as hypothesis testing and model selection in a manner similar to applying a parametric conjugate prior, while attaining full nonparametricity. Also provided is a comparison to two other state-of-the-art Bayesian nonparametric models for conditional densities in both model fit and computational time. A real data example from flow cytometry containing 455,472 observations is given to illustrate the substantial computational efficiency of our method and its application to multivariate problems.


Introduction
In recent years there has been growing interest in nonparametrically modeling probability densities based on multi-scale partitioning of the sample space. A prime example in the Bayesian nonparametric literature is the Pólya tree (PT) [12,22,31] and its extensions [17,18,45,21,27]. In particular, Wong and Ma [45] introduced randomization into the partitioning component (involving both random selection of partition directions as well as optional stopping) of the PT framework, which enhances the model's ability to approximate the shape and smoothness of the underlying density. A PT model with these features is called an optional Pólya tree (OPT).
A further desirable feature of the PT and its relatives such as the OPT and the more recently introduced adaptive Pólya tree (APT) [27] is the computational ease for carrying out inference. In turns out that the extra component of randomized partitioning such as that employed in the OPT does not impair the conjugacy enjoyed by the PT. For example, after observing i.i.d. data, the corresponding posterior of an OPT is still an OPT, that is, the same generative procedure for random probability distributions with its parameters updated to their posterior values. Moreover, the corresponding posterior parameter values can be computed exactly through a sequence of recursive computations, which is in essence a forward-backward algorithm [25]. This, together with the constructive nature of these models, allows one to draw samples from the exact posterior directly without resorting to Markov Chain Monte Carlo (MCMC) procedures, and to compute various summary statistics of the posterior analytically. Furthermore, the marginal posterior of the random partitioning adapts to the underlying structure of the data-the sample space will with high posterior probability be more finely divided in places where the underlying distribution has richer structure, i.e. less uniform topological shape.
Motivated by the computational efficiency and statistical properties of the OPT, which is tied to its use of recursive random partitioning, we aim to further exploit the random recursive partitioning idea in the context of multi-scale density modeling, and build such a model for conditional densities for a response (vector) Y given a predictor (vector) X. The objective is to construct a flexible nonparametric model for conditional distributions that maintain all of the desirable statistical and computational properties of PT and OPT.
A variety of inference tasks involve the estimation, prediction, and testing regarding conditional distributions, and nonparametric inference on conditional densities has been studied from both frequentist and Bayesian perspectives. Many frequentist works are based on kernel estimation methods [10,16,11], and they achieve proper smoothing through bandwidth selection, which often involves resampling procedures such as cross-validation [2,19,11] and the bootstrap [16]. An alternative frequentist strategy introduced more recently is to employ the so-called block-wise shrinkage [8,9]. In Bayesian nonparametrics, inference on conditional distributions is often referred to as covariate-dependent distribution modeling, and existing methods fall into two categories. The first category is methods that construct priors for the joint distribution of the response and the predictors, and then use the induced conditional distribution for inference. Some examples are [32,37,33,41], which propose using mixtures of multivariate normals as the model for joint distributions, along with different priors for the mixing distribution. The other category is methods that construct conditional distributions directly without specifying the marginal distribution of the predictors. Many of these methods are based on extending the stick breaking construction for the Dirichlet Process (DP) [39]. Some notable examples, among others, are proposed in [29,20,13,15,7,4,36,1]. Some recent works in this category do not utilize stick breaking. In [43], the authors propose to use the logistic Gaussian process [23,42] together with subspace projection to construct smooth conditional distributions. In [21], the authors incorporate covariate dependency into tail-free processes by generating the conditional tail probabilities from covariate-dependent logistic Gaussian processes, and propose a mixture of such processes as a way for modeling conditional distributions. The authors of [24] introduce dependent normalized complete random measures. In [44] the authors introduce the covariate-dependent multivariate Beta process, and use it to generate the conditional tail probabilities of Pólya trees. More recently, in [40] the authors use the tensor product of B-splines to construct a prior for conditional densities, and incorporate a variable selection feature. While many of these nonparametric models on conditional distributions enjoy desirable theoretical properties, inference using these priors generally relies on intense MCMC sampling, and can take substantial computing time even when both the response and the covariate are one-dimensional.
We introduce a new prior, called the conditional optional Pólya tree, for the conditional density of Y given X, in the form of a two-stage generative procedure consisting of first randomly partitioning the predictor space Ω X , and then for each predictor partition block, generates the response distribution on each block using an OPT, which implicitly employs a further random partitioning of the response space Ω Y . We show that this new prior is a fully nonparametric model and yet achieves extremely high computational efficiency even for multivariate responses and covariates. It enjoys all of the desirable theoretical properties of the PT and the OPT priors-namely large support, posterior consistency, and posterior conjugacy, and its posterior parameters can also be computed exactly through forward-backward recursion. Under this two-stage design, the posterior distribution on the partitions reflect the structure of the conditional distribution at two levelsfirst, the predictor space will be partitioned finely in parts where the conditional distribution changes most abruptly, shedding light on how the conditional distribution depends on the predictors; second, the response space will be divided adaptively for different locations of the predictor space, to capture the local structure of the conditional density through adaptive smoothing.
The rest of the paper is organized as follows. In Section 2 we introduce our two-stage prior and show that it is fully nonparametric-with full (integrated) L 1 support-for conditional densities. In addition, we make a connection to Bayesian CART and show that our method can be considered a nonparametric version of the latter. In Section 3 we show the full conjugacy of the model, derive the exact form of the posterior through forward-backward recursion, and thereby provide a recipe for carrying out Bayesian inference using the prior. We also prove the posterior consistency of such inference. In Section 4 we discuss practical computational issues in implementing the inference. In Section 5 we provide four simulation examples to illustrate the work of our method. The first two are for estimating conditional densities, and the last two concern model selection and hypothesis testing. In Section 6 we apply the proposed method to estimating conditional densities in a flow cytometry data set involving a large number (455,472) of observations, and demonstrate the computational efficiency of the method and its application when both the response and the predictor are multivariate. Section 7 concludes with some discussions. All proofs are given in the Appendix.

Conditional optional Pólya trees
In this section we introduce our proposed prior constructively in terms of a two-stage generative procedure that produces random conditional densities. First we introduce some notions and notations that will be used throughout. Let each observation be a predictor-response pair (X, Y ), where X denotes the predictor (or covariate) vector and Y the response (vector) with Ω X being the predictor space and Ω Y the response space. In this work we consider sample spaces that are either finite spaces, compact Euclidean rectangles, or a product of the two, and Ω X and Ω Y do not have to be of the same type. (See for instance Example 3.) Let µ X and µ Y be the "natural" measures on Ω X and Ω Y . (That is, the counting measure for finite spaces, the Lebesgue measure for Euclidean rectangles, and the corresponding product measure if the space is a product of the two.) Let µ = µ X × µ Y be the "natural" product measure on the joint sample space Ω X × Ω Y .
A partition rule R on a sample space Ω specifies a collection of possible ways to divide any subset A of Ω into a number of smaller sets. For example, for Ω = [0, 1] k , the unit rectangle in R k , the coordinate-wise dyadic mid-split rule allows each rectangular subset A of Ω whose sides are parallel to the k coordinates to be divided into two halves at the middle of the range of each coordinate. For simplicity, in this work we only consider partition rules that allow a finite number of ways for dividing each set. Such partition rules are said to be finite. (Interested readers can refer to [28,Sec. 2] for a more detailed treatment of partition rules and to Examples 1 and 2 in [45] for examples of the coordinate-wise dyadic mid-split rule over Euclidean rectangles and 2 k contingency tables.) We are now ready to introduce our prior for conditional densities as a two-stage constructive procedure. It is important to note that the following describes the generation of conditional densities under our prior and not the operational steps for inference under the prior, which will be addressed Section 3 and Section 4.
Stage I. Predictor partition: We randomly partition Ω X according to a given partition rule R X on Ω X in the following recursive manner. Starting from A = Ω X , draw a Bernoulli variable Let K j (A) be the number of child sets that arise from this partition, and let A j 1 , A j 2 , . . . , A j K(A) denote these children. We then repeat the same partition procedure, starting from the drawing of a stopping variable, on each of these children.
The following lemma, first proved in [45], states that as long as the stopping probabilities are (uniformly) away from 0, this random recursive partitioning procedure will eventually terminate almost everywhere and produce a well-defined partition of Ω X . Lemma 1. If there exists a δ > 0 such that the stopping probability ρ(A) > δ for all A ⊂ Ω X that could arise after a finite number of levels of recursive partition, then with probability 1 the recursive partition procedure on Ω X will stop µ X a.e.
Stage II. Generating conditional densities: Next we move onto the second stage of the procedure to generate the conditional density of the response Y on each of the predictor partition blocks generated in Stage I. Specifically, for each stopped subset A on Ω X produced in Stage I, we let the conditional distribution of Y given X = x be the same across all x ∈ A, and generate this (conditional) distribution on Ω Y , denoted as q 0,A Y , from a "local" prior. When the response space Ω Y is finite, q 0,A Y is simply a multinomial distribution, and so a simple choice of such a local prior is the Dirichlet prior: represents the pseudo-count hyperparameters of the Dirichlet. In this case, we note that the two-stage prior essentially reduces to a version of the Bayesian CART proposed by Chipman et al in [3] for the classification problem. When Ω Y is infinite (or finite but with a large number of elements), one may restrict q 0,A Y to be from a parametric family. For example, when Ω Y = R, one may require q 0,A Y to be normal with some mean µ A and variance σ 2 A , and let µ A |σ 2 A ∼ N(µ 0 , σ 2 ) and σ 2 A ∼ inverse-Gamma(ν/2, νκ/2). In this case the two-stage prior again reduces to a Bayesian CART, this time for the regression problem [3].
The focus of our current work, however, is on the case when no parametric assumptions are placed on the conditional density. To this end, one can draw q 0,A Y from a nonparametric prior. A desirable choice for the local prior, which will result in analytic simplicity and computational efficiency as we will later show, is a Pólya tree type model [27], and in particular an optional Pólya tree (OPT) distribution [45]: independently across As given the partition, where R A Y denotes a partition rule on Ω Y and ρ A Y , λ A Y , and α A Y are respectively the stopping, selection, and pseudo-count hyperparameters [45]. In general we allow the partition rule for these "local" OPTs to depend on A as indicated in the superscript, but adopting a common partition rule on Ω Y -that is to let R A Y ≡ R Y for all A-will suffice for most problems. In the rest of the paper, unless stated otherwise we assume that a common rule R Y is adopted.
This completes the description of our two-stage procedure. We now formally define the resulting prior.
Definition 1. A conditional distribution that arises from the above two-stage procedure is said to have a conditional optional Pólya tree (cond-OPT) distribution. The hyperparameters are the predictor partition rule R X , the response partition rule R Y , the stopping probability ρ(A), the partition selection probabilities λ(A), and the local parameters (ρ A Y , λ A Y , α A Y ) for all A ⊂ Ω X that could arise during the predictor partition under R X .
Remark I: To ensure that this definition is meaningful, one must check that the two-stage procedure will in fact generate a well-defined conditional distribution with probability 1. To see this, first note that because the collection of all potential sets A on Ω X that can arise during Stage I is countable, by Theorem 1 in [45], with probability 1, the two-stage procedure will generate an absolutely continuous conditional distribution of Y given X = x for x in the stopped part of Ω X , provided that ρ A Y is uniformly away from 0. The twostage generation procedure for the conditional density of Y can then be completed by letting Y given X be uniform on Ω Y for the µ X -null subset of Ω X on which the recursive partition in Stage I never stops.
Remark II: While the cond-OPT prior involves many hyperparameters, one can appeal to very simple symmetry and self-similarity principles for choosing their values. Specifically, such considerations lead to the simple choice: following the default choices in [45]. We note that when useful prior knowledge about the structure of the underlying distribution is not available or when one is unwilling to assume particular structure over the distribution, it is desirable to specify the prior parameters in a symmetric and self-similar way. The common stopping probability ρ should not be too close to 0 or 1, but taking a moderate value between 0.1 and 0.9. A sensitivity analysis for such choices demonstrating the robustness of such choices in the context of OPTs is provided in [28]. As for the partition rules, the coordinate-wise dyadic mid-split rule can serve as a simple default choice for both R X and R Y . We will adopt such a specification in all of our numerical examples.
Remark III: One constraint in the cond-OPT is that given the random partition generated in Stage I, the generation of the conditional distribution across different predictor blocks is independent, i.e., in a similar manner as that for Bayesian CART. As we shall see, this constraint is key to the tremendous computational 4 efficiency of the model. It is important to note however that due to the randomized partitioning incurred in Stage I, the marginal prior for the conditional distributions on nearby values of X are in fact dependent, thereby achieving smoothing over Ω X to some extent. More flexible smoothing could be achieved through modeling the "local" priors jointly, but that would incur the need for MCMC sampling and the most desirable feature of PT type models would be lost.
We have emphasized that the cond-OPT prior imposes no parametric assumptions on the conditional distribution. One may wonder whether this prior is truly "nonparametric" in the sense that it can generate all possible conditional densities. Our next theorem confirms this-under mild conditions on the parameters, which the default specification satisfies, the cond-OPT will place positive probability in arbitrarily small L 1 neighborhoods of any conditional density. (A definition of an L 1 neighborhood for conditional densities is also implied in the statement of the theorem.) Theorem 2 (Large support). Suppose q(·|·) is a conditional density function that arises from a cond-OPT prior whose parameters ρ(A) and λ(A) for all A that could arise during the recursive partitioning on Ω X are uniformly away from 0 and 1, and the local OPTs all have full L 1 support on the densities on Ω Y . Moreover, suppose that the underlying partition rules R X and R Y both satisfy the following "fine partition criterion": ∀ > 0, there exists a partition of the corresponding sample space such that the diameter of each partition block is less than . Then for any conditional density function f (·|·) : Ω Y × Ω X → [0, ∞), and any τ > 0, Furthermore, let f X (x) be any density function on Ω X w.r.t. µ X . Then we have ∀τ > 0, Remark: Sufficient conditions for OPTs to have full L 1 support on densities is given in Theorem 2 of [45].

Bayesian inference with cond-OPT
Next we investigate how Bayesian inference on conditional densities can be carried out using this prior. First, we note that Chipman et al [3] and Denison et al [6] each proposed MCMC algorithms that enable posterior inference for Bayesian CART. These sampling and stochastic search algorithms can be applied directly here as the local OPT priors can be marginalized out and so the marginal likelihood under each partition tree that arises in Stage I of the cond-OPT is available in closed form [45,28]. However, as noted in [3] and other works, due to the multi-modal nature of tree structured models, the mixing behavior of the MCMC algorithms is often undesirable. This problem is exacerbated in higher dimensional settings. Chipman et al [3] suggested using MCMC as a tool for searching for good models rather than a reliable way of sampling from the actual posterior.
The main result of this section is that under simple partition rules such as the coordinate-wise dyadic mid-split rule, Bayesian inference under a cond-OPT prior can be carried out in an exact manner in the sense that the corresponding posterior distribution can be computed in closed form and directly sampled from, without resorting to MCMC algorithms. Not only is the computation feasible for multivariate sample spaces of moderate dimensions, but it is in fact highly efficient, scaling linearly with the number of observations. First let us investigate what the posterior of a cond-OPT prior is. Suppose we have observed (x, y) = {(x 1 , y 1 ), (x 2 , y 2 ), . . . , (x n , y n )} where given the x i 's, the y i 's are independent with some density q(y|x). We assume that q(·|·) has a cond-OPT prior denoted by π. Further, for any A ⊂ Ω X we let Then conditional on the event that A arises during the recursive partition procedure on Ω X , we can write q(A) recursively in terms of S(A), J(A), and q A Y as follows the likelihood from the data with x ∈ A if the partitioning stops on A. Equivalently, we can write Integrating out the randomness over both sides of Eq. (3.1), we get is defined to be the marginal likelihood from data with x ∈ A given that A arises during the recursive partitioning on Ω X , whereas is the marginal likelihood from the data with x ∈ A if the recursive partitioning procedure stops on A and the integration is taken over the local OPT( We note that Eqs. 2) provides a recursive recipe for calculating Φ(A) for all A. It is recursive in the sense that Φ(A) is computed based on the value of Φ(·) on A's children. (Of course, to complete the calculation the recursion must eventually terminate everywhere on Ω X . We shall describe the terminal conditions in the next section.) This recursive algorithm is a special case of the forward-backward algorithm [27].
The next theorem establishes the posterior conjugacy of cond-OPT.
Theorem 3 (Conjugacy). After observing {(x 1 , y 1 ), (x 2 , y 2 ), . . . , (x n , y n )} where given the x i 's, the y i 's are independent with density q(y|x), which has a cond-OPT prior, the posterior of q(·|·) is again a cond-OPT (with the same partition rules on Ω X and Ω Y as the prior). Moreover, for each A ⊂ Ω X that could arise during the recursive partitioning, the posterior parameters are given as follows.
1. Stopping probability: 2. Selection probabilities: 3. The local parameters:ρ A Y ,λ A Y , andα A Y are the corresponding posterior parameters for the local OPT after updating using the observed values for the response y(A), This theorem shows that a posteriori our knowledge about the underlying conditional distribution of Y given X can again be represented by the same two-stage procedure that randomly partitions the predictor space and then generates the response distribution accordingly on each of the predictor blocks, except that now the parameters that characterize this two-stage procedure have been updated to reflect the information contained in the data. Moreover, the theorem also provides a recipe for computing these posterior parameters based on Φ(A) and M (A). Given this exact posterior, Bayesian inference can then proceed-samples can be drawn from the posterior cond-OPT directly through vanilla Monte Carlo (as opposed to MCMC) and summary statistics calculated.
In the next section, we provide more details on how to implement such inference in practice. Before that, we present our last theoretical result about the cond-OPT prior-its posterior consistency, which assures the statistician that the posterior cond-OPT distribution will "converge" in some sense to the truth as the amount of data increases. To this end, we first need a notion of neighborhoods for conditional densities under which such convergence holds. We adopt the notion discussed in [35] and [34], by which a (weak) neighborhood of a conditional density function is defined in terms of a (weak) neighborhood of the corresponding joint density. More specifically, for a conditional density function f 0 (·|·) : Ω Y × Ω X → [0, ∞), weak neighborhoods with respect to a marginal density f 0 X (·) on Ω X are collections of conditional densities of the form where the g i 's are bounded continuous functions on Ω X × Ω Y .
Theorem 4 (Weak consistency). Let (x 1 , y 1 ), (x 2 , y 2 ), . . . be independent identically distributed vectors from a probability distribution on Suppose the conditional density f (·|·) is generated from a cond-OPT prior for which the conditions in Theorem 2 all hold.

Practical implementation
Next we address some practical issues in computing the posterior and implementing the inference. For simplicity, from now on we shall refer to a set A ⊂ Ω X that can arise during the (Stage I) recursive partitioning procedure as a "node" (i.e., as a node in the partition tree). A prerequisite for applying Theorem 3 is the availability of the Φ(A) terms, which can be determined recursively through Eq. (3.2). Of course, to carry out the computation of Φ(A) one must specify terminal conditions on Eq. (3.2), or in other words, on what kind of A's the recursion should terminate. We call such nodes terminal nodes.
There are two kinds of nodes for which the value of Φ(A) is available directly according to theory, and thus recursion can terminate on them. They are (i) nodes that cannot be further divided under the partition rule R X , and (ii) nodes that contain no more than one data point. For a node A that cannot be further divided, we must have ρ(A) = 1 and so Φ(A) = M (A). For a node A with no data point, it has no contribution to the likelihood and so Φ(A) = 1. For a node A with exactly one data point, Φ(A) is the predictive density of the local OPT on A evaluated at that data point, which is exactly the density of the prior mean of the local OPT and is directly known when the default symmetric and self-similar prior specification for the local OPTs is adopted as recommended in [45].
Note that with these two types of "theoretical" terminal nodes, in principle the recursion will eventually terminate if one divides the predictor space deep enough. In practice, however, it is unnecessary to take the recursion all the way down to these theoretical terminal nodes. Instead, one can adopt early termination by imposing a technical limit-such as a minimum size (or maximum depth) of the nodes either in terms of the natural measure µ X (A) or the number of observations therein n(A)-to end the recursion. Nodes that are 7 smaller than the chosen size threshold are forced to be terminal, which is equivalent to setting ρ(A) = 1 and thus Φ(A) = M (A) for these nodes. We call these nodes "technical" terminal nodes. With these theoretical and technical terminal nodes, one can then compute Φ(A) through the recursion formula (3.2), and compute the posterior according to Theorem 3. Putting all the pieces together, we can summarize the procedure to carry out Bayesian inference with the cond-OPT prior as a four-step recipe: For the last step, direct simulation from the posterior is straight-forward, but we have not discussed what summary statistics to compute and how to do that. This is problem-specific and will be illustrated in our numeric examples in Section 5.

Examples
In this section we provide four examples to illustrate inference using the cond-OPT prior. The first two illustrate the estimation of conditional densities, the latter two are for model selection and hypothesis testing.
In these examples, the partition rules used on both Ω X and Ω Y are always the coordinate-wise dyadic midsplit rule. We adopt the same prior specification across all the examples: the prior stopping probability on each non-terminal node is always set to 0.5, the prior partition selection probability is always evenly spread over the possible ways to partition each set, and the probability assignment pseudo-counts for the local OPTs are all set to 0.5. For continuous sample spaces, nodes at 12 levels down the partition tree, i.e., with µ X (A) = µ(Ω X )/2 12 , are set to be the technical terminal nodes.
Example 1 (Estimating conditional density with abrupt changes over predictor values). In this example we simulate (X, Y ) pairs according to the following distributions. We generate data sets of three different sample sizes, n = 100, n = 500, and n = 2, 500, and place the cond-OPT prior on the distribution of Y given X. Following the four-step recipe given in the previous section, we can compute the posterior cond-OPT and sample from it. A representative summary of the posterior partitioning mechanism is the so-called hierarchical maximum a posteriori (hMAP) [45] partition tree, which can be computed from the posterior analytically [45] and is plotted in Figure 1 for the different sample sizes. (Chipman et al [3] and Wong and Ma [45] both discussed reasons why the commonly adopted MAP is not a good summary for tree-structured posteriors due to their multi-level nature. See [45,Sec. 4.2] for further details and reasons why the hMAP is often preferred to the MAP.) In Figure 1, within each "leaf" node we plot the corresponding posterior mean of the local OPT. Also plotted for each node is the posterior stopping probability. Even with only 100 data points, the posterior suggests that Ω X should be divided into three pieces-[0,0.25], [0.25,0.5], and [0.5,1]-within which the conditional distribution of Y |X is homogeneous across X. Note that the posterior stopping probabilities on those three intervals are large, in contrast to the near 0 values on the larger sets. Reliably estimating the actual conditional density function on these sets nonparametrically appears to require more than 100 data points. In this example, a sample size of 500 already does a decent job. We compare both the model fit and the computing speed of our cond-OPT prior to two existing Bayesian nonparametric models for conditional densities-namely the linear dependent Dirichlet process mixture of normals (LDDP) [5] and the linear dependent Dirichlet process mixture of Bernstein polynomials (LDBP) [1], both available in the DPpackage in R. In this example and the next, for LDDP and LDBP, we draw 1,000 posterior samples from the MCMC with a 2,000 burn-in period and a thinning interval of 3, and used prior specification given in the examples of the DPpackage. For details, please see the documentation for these two functions in the DPpackage manual on CRAN.
To evaluate model fit, we generate an additional testing data set from the true distribution of (X, Y ), and calculate the log-p score (i.e., the log predictive likelihood of the testing set) for the three methods. Table 1 presents the log-p score for the three methods from a typical simulated data set and the corresponding computing time on the same laptop computer with an Intel Core-i7 CPU using a single core without parallelization. A surprising phenomenon is that the performance of LDBP, in terms of the log-p score for the testing sample, is not always monotone increasing in the sample size-that is, a larger training sample does not always lead to better fit on the testing set. In the particular simulation reported in Table 1, the preformance of LDBP is actually monotone decreasing with sample size. The cause for this is likely to be that under those models the conditional density is assumed to be smoothly varying over the predictors, and so as the true conditional density involves abrupt changes, the misspecified models can be consistently wrong even with large sample sizes. The previous example favors our method because (1) there are a small number of clear boundaries of change for the underlying conditional distribution, and to a lesser extent (2) those boundaries-namely 0.25 and 0.5-lie on the potential partition points of the partition rule. In the next example, we examine the case in which the conditional distribution changes smoothly across a continuous X without any boundary of abrupt change.
Example 2 (Estimating conditional densities that vary smoothly with predictor values). In this example we generate (X, Y ) from a bivariate normal distribution. We generate a data set of size n = 2, 000, and apply the cond-OPT prior on the distribution of Y given X as we did in the previous example. Again we compute the posterior cond-OPT following our four-step recipe. The hMAP tree and the posterior mean estimate of the conditional density given the random partition is presented in Figure 2. Because the underlying predictor space Ω X is unbounded, for simplicity in the above we used the empirically observed range of X as Ω X , which happens to be Ω X = [0.24, 0.92] for our simulated example. (Other ways to handle this situation include transforming X to have a compact support such as through a CDF or rank transform. One interesting observation is that the "leaf" nodes in Figure 2 have very large (close to 1) posterior stopping probability. This may seem surprising as the underlying conditional distribution is not the same for any neighboring values of X. The large posterior stopping probabilities indicate that on those sets, where the sample size is not large, the gain in achieving better estimate of the common features of the conditional distribution for nearby X values outweighs the loss in ignoring the difference among them.
Again, to compare the model fit and computational efficiency with LDDP and LDBP, we repeat a set of simulations with different sample sizes n = 100, 500, and 2500, and again use the log-p score on a testing sample of size 100 to evaluate the performance. The results are summarized in Table 2, and they mostly confirm our intuition-the smooth priors overall outperform our model, especially for small sample sizes. The performance difference vanishes as the sample size increases.
In other words, three predictors X 5 , X 20 and X 30 impact the response in an interactive manner. Our interest is in recovering this underlying interactive structure (i.e. the "model"). To illustrate, we simulate 500 data points from this scenario and place a cond-OPT prior on Y |X, and consider predictor partitions up to four levels deep. This is achieved by setting ρ(A) = 1 for A that arises after four steps of partitioning, and it allows us to search for models involving up to four-way interactions. We again carry out the four-step recipe to get the posterior and calculate the hMAP. The hMAP tree structure along with the predictive conditional density for Y |X within each stopped set given the random partition is presented in Figure 3. The posterior concentrates on partitions involving X 5 , X 20 and X 30 out of the 30 variables. While the predictive conditional density for Y |X is very rough given the limited number of data points in the stopped sets, the posterior recovers the exact interactive structure of the predictors with little uncertainty.
In addition, we sample from the posterior and use the proportion of times each predictor appears in the sampled models to estimate the posterior marginal inclusion probabilities. Our estimates based on 1,000 draws from the posterior are presented in Figure 4(a). Note that the sample size 500 is so large that the posterior marginal inclusion probabilities for the three relevant predictors are all close to 1 while those for The hMAP tree structure on Ω X and the posterior mean estimate of Y |X given the random partition in each of the stopped sets for Example 3. The bold arrows indicate the "true model"-predictor combinations that correspond to "non-null" Y |X distributions. For each node, ρ indicates the posterior stopping probability for each node, λ represents the posterior selection probability for the most probable direction if the partition does not stop on the node, and n represents the number of data points in each node.  the other predictors are close to 0. We carry out the same simulation with a reduced sample size of 200, and plot the estimated posterior marginal inclusion probabilities in Figure 4(b). We see that with a sample size of 200, one can already use the posterior to reliably recover the relevant predictors.
Example 4 (Test of independence). In this example, we illustrate an application of the cond-OPT prior for hypothesis testing. In particular, we use it to test the independence between X and Y . To begin, note that ρ(A|x, y) in Theorem 3 gives the posterior probability for the conditional distribution of Y to be constant over all values of X in A, or in other words, for Y to be independent of X on A. Hence, one can consider ρ(Ω X |x, y) as a score for the statistical significance of dependence between the observed variables. A permutation null distribution of this statistic can be constructed by randomly pairing the observed x and y values, and based on this, permutation p-values can be computed for testing the null hypothesis of independence.
To illustrate, we simulate X = (X 1 , X 2 , . . . , X 10 ) for a sample of size 400 under the same Markov Chain model as in the previous example, and simulate a response variable Y as follows.
In particular, Y is dependent on X but there is no mean or median shift in the conditional distribution of Y over different values of X. Figure 5 gives the histogram of ρ(Ω X |x, y) for 1,000 permuted samples where the vertical dashed line indicates the ρ(Ω X |x, y) for the original simulated data, which equals 0.0384. For this particular simulation, 7 out of the 1,000 permuted samples produced a more extreme test statistic. Remark I: Note that by symmetry one can place a cond-OPT prior on the conditional distribution of X given Y as well and that will produce a corresponding posterior stopping probability ρ(Ω Y |y, x). One can thus alternatively use min{ρ(Ω X |x, y), ρ(Ω Y |y, x)} as the test statistic for independence.
Remark II: Testing using the posterior stopping probability ρ(Ω X |x, y) is equivalent to using a Bayes factor (BF). To see this, note that the BF for testing independence under the cond-OPT can be written as with A = Ω X where the numerator is the marginal conditional likelihood of Y given X if the conditional distribution of Y is not constant over X (i.e. Ω X is divided) and the denominator is that if the conditional distribution of Y is the same for all X (i.e. Ω X is undivided). By Eq. (3.2) and Theorem 3, which is in a one-to-one correspondence to ρ(Ω X |x, y) given the prior parameters.
6. Application to real data: multivariate conditional density estimation in flow cytometry In flow cytometry experiments for immunological studies, a number (typically 4 to 10) of biomarkers are measured on large numbers of blood cells. Estimated densities and conditional densities of such data can be used for tasks such as automatic classification of the cells [30]. We apply cond-OPT to estimate the conditional density of markers "CD4" and "CD8" given two other markers "FSC-H" and "FSC-W" in a flow cytometry data set. So in this case both Ω X and Ω Y are two-dimensional. This particular data set contains n = 455, 472 cells. Flow cytometry experiments often involve large numbers of cells, and thus practical 13 methods must scale well in computing time and memory usage with respect to the number of observations. This poses great challenge to existing nonparametric models that require intense MCMC computation. The values of the four markers are measured in the range of [0,1]. We use maximum level of partitioning to 10 on both the predictor space Ω X and the response space Ω Y but otherwise the same prior specification as before. Figure 6 presents the posterior mean of the conditional density of CD4 and CD8 given FSC-H and FSC-W under the cond-OPT model given the random partition on the predictor space being the one induced under the hMAP tree, which splits the space into 50 pieces. A vast majority, in fact 44 out of the 50 predictor blocks are in fact not technical terminal regions, and so the model indeed smooths the conditional density over the predictor space. Because the number of predictor blocks is relatively large, we present the estimates for only 16 blocks in Figure 6. The entire computation of the full posterior, the hMAP partition, as well as the conditional posterior expectation of the conditional density given the hMAP tree, took about 360 seconds to complete on a single 3.6GHz Intel Core-i7 3820 desktop core without parallelization and required about 8.2 Gbs of RAM. (Reducing the maximum level of partitions from 10 to 8 will reduce computing time to about 116 seconds and RAM to about 0.6 Gbs.)

Discussion
In this work we have introduced a Bayesian nonparametric prior on the space of conditional densities. This prior, which we call the conditional optional Pólya tree, is constructed based on a two-stage procedure that first divides the predictor space Ω X and then generates the conditional distribution of the response through local OPT processes. We have established several important theoretical properties of this prior, namely large support, conjugacy and posterior consistency, and have provided a practical recipe for Bayesian inference using this prior.
The construction of this prior does not depend on the marginal distribution of X. One particular implication is that one can transform X before applying the prior on Y |X without invalidating the posterior inference. (Note that transforming X is equivalent to choosing a different partition rule on Ω X .) In certain situations it is desirable to perform such a transformation on X. For example, if the data points are very unevenly spread over Ω X , then some parts of the space may contain a very small number of data points. There the posterior is mostly dominated by the prior specification and does not provide much information about the underlying conditional distribution. One way to mitigate this problem is to transform X so that the data are more evenly distributed over Ω X . When Ω X is one-dimensional, for example, this can be achieved by a rank transformation on X. Another situation in which a transformation of X may be useful is when the dimensionality of X is very high. In this case a dimensionality reduction transformation can be applied on X before carrying out the inference. Of course, in doing so one often loses the ability to interpret the posterior conditional distribution of Y directly in terms of the original predictors. An alternative approach when X is high-dimensional is through variable selection that imposes certain sparsity assumptions, i.e., only a small number of predictors are affecting the conditional density. Exact calculation of full posterior and the marginal inclusion probabilities as we have carried out in Example 3 is impractical when the number of predictors is large (> 25 ∼ 30). One strategy to overcome this difficulty is through sequential importance sampling as the one proposed in [26].
A general limitation of CART type randomized partitioning methods require a natural ordering of the space to be partitioned on. General partitioning strategies can be designed for unordered spaces, but then the computational efficiency of the proposed model would be lost.
Finally, we note that while we have used recursive partitioning in conjunction with the OPT to build a model for conditional density, one can build such models by replacing the OPT with other multi-scale density models in the family of Pólya tree type models, such as the more recently introduced adaptive Pólya tree (APT) [27].

Software
The proposed model has been implemented in the R package PTT (for Pólya tree type models) as the function cond.opt. A variant of the model that replaces the OPT with an APT is also implemented in the package as 14 function cond.apt. This package is currently available for download at https://github.com/MaStatLab/ PTT and will be submitted to CRAN.

Acknowledgment
The flow cytometry data set was provided by EQAPOL (HHSN272201000045C), an NIH/NIAID/DAIDSsponsored, international resource that supports the development, implementation, and oversight of quality assurance programs (Sanchez PMC4138253).

Appendix: Proofs
Proof of Lemma 1. The proof of this lemma is very similar to that of Theorem 1 in [45]. Let T k 1 be the part of Ω X that has not been stopped after k levels of recursive partitioning. The random partition of Ω X after k levels of recursive partitioning can be thought of as being generated in two steps. First suppose there is no stopping on any set and let J * (k) be the collection of partition selection variables J generated in the first k levels of recursive partitioning. Let A k (J * (k) ) be the collection of sets A that arise in the first k levels of non-stopping recursive partitioning, which is determined by J * (k) . Then we generate the stopping variables S(A) for each A ∈ A k (J * (k) ) successively for k = 1, 2, . . ., and once a set is stopped, let all its descendants be stopped as well. Now for each A ∈ A k (J * (k) ), let I k (A) be the indicator for A's stopping status after k levels of recursive partitioning, with I k (A) = 1 if A is not stopped and = 0 otherwise.
Proof of Theorem 2. We prove only the second result as the first follows by choosing f X (x) ≡ 1/µ X (Ω X ). Also, we consider only the case when Ω X and Ω Y are both compact Euclidean rectangles, because the cases when at least one of the two spaces is finite follow as simpler special cases. For x ∈ Ω X and y ∈ Ω Y , let f (x, y) := f X (x)f (y|x) denote the joint density. First we assume that the joint density f (x, y) is uniformly continuous. In this case it is bounded on Ω X × Ω Y . We let M := sup f (x, y) and |f (x 1 , y 1 ) − f (x 2 , y 2 )|.
By uniform continuity, we have δ( ) ↓ 0 as ↓ 0. In addition, we define Note that in particular the continuity of f (x, y) implies the continuity of f X (x). Let σ > 0 be any positive constant. Choose a positive constant (σ) such that δ X ( (σ)) = δ( (σ))µ Y (Ω Y ) < max(σ/2, σ 3 /2). Because all the parameters in the cond-OPT are uniformly bounded away from 0 and 1, there is positive probability that Ω X will be partitioned into Ω X = ∪ K i=1 B i where the diameter of each B i is less than (σ), and the partition stops on each of the B i 's. (The existence of such a partition follows from the fine partition criterion.) Let A i = B i ∩ {X : f X (x) ≥ σ}, P (X ∈ A i ) = Ai f X (x)µ X (dx), and f i (y) := Ai f (x, y)µ X (dx)/µ X (A i ) if µ X (A i ) > 0, and 0 otherwise. Let I ⊂ {1, 2, . . . , K} be the set of indices i such that µ X (A i ) > 0. Then Let us consider each of the four terms on the right hand side in turn. First, Note that for each i ∈ I, f i (y)µ X (A i )/P (X ∈ A i ) is a density function in y. Therefore by the large support property of the OPT prior (Theorem 2 in [45]), with positive probability, for all i ∈ I. Also, for any x ∈ A i , by the choice of (σ), ≤ δ X ( (σ)) σ(σ − δ X ( (σ)) ≤ σ 3 /2 σ 2 /2 = σ. Thus Finally, again by the choice of (σ), |f i (y) − f (x, y)| ≤ δ( (σ)) < σ, and so Therefore for any τ > 0, by choosing a small enough σ, we can have |q(y|x) − f (y|x)|f X (x)µ(dx × dy) < τ with positive probability. This completes the proof of the theorem for continuous f (x, y). Now we can approximate any density function f (x, y) arbitrarily close in L 1 distance by a continuous onef (x, y). The theorem still holds because |q(y|x) − f (y|x)|f X (x)µ(dx × dy) ≤ q(y|x)|f X (x) −f X (x)|µ(dx × dy) + |q(y|x) −f (y|x)|f X (x)µ(dx × dy) + |f (x, y) − f (x, y)|µ(dx × dy).
≤ |q(y|x) −f (y|x)|f X (x)µ(dx × dy) wheref X (x) andf (y|x) denote the corresponding marginal and conditional density functions forf (x, y). Thus the conjugacy of the prior and the posterior updates for ρ, λ j and OPT(R A Y ; ρ A Y , λ A Y , α A Y ) follows from Bayes' Theorem and the posterior conjugacy of the standard optional Pólya tree prior (Theorem 3 in [45]).

Proof of Theorem 3. Given that a set
Proof of Theorem 4. By Theorem 2.1 in [34], which follows directly from Schwartz's theorem (see [38] and [14,Theorem 4.4.2]), we just need to prove that the prior places positive probability mass in arbitrarily small Kullback-Leibler (K-L) neighborhoods of f (·|·) w.r.t f X . Here a K-L neighborhood w.r.t f X is defined to be the collection of conditional densities for some > 0.
To prove this, we just need to show that any conditional density that satisfies the conditions given in the theorem can be approximated arbitrarily well in K-L distance by a piecewise constant conditional density of the sort that arises from the cond-OPT procedure. We first assume that f (·|·) is continuous. Following the proof of Theorem 2, let δ( ) denote the modulus of continuity of f (·|·). Let Ω X = ∪ K i=1 A i be a reachable partition of Ω X such that the diameter of each partition block A i is less than . Next, for each A i , let Ω Y = ∪ N j=1 B ij be a partition on Ω Y allowed under OPT(R Y ; ρ Ai Y , λ Ai Y , α Ai Y ) such that the diameter of each B ij is also less than . Let g ij = sup x∈Ai,y∈Bij