arXiv:2011.03121v3 [stat.ME] 7 Dec 2021

Hidden Markov Pólya trees for high-dimensional
distributions
Naoki Awaya and Li Ma
Department of Statistical Science
Duke University, Durham, NC 27708, USA
December 9, 2021
Abstract
The Pólya tree (PT) process is a general-purpose Bayesian nonparametricmodel that has found wide application in a range of inference problems. It
has a simple analytic form and the posterior computation boils down to betabinomial conjugate updates along a partition tree over the sample space. Recent
development in PT models shows that performance of these models can be substantially improved by (i) allowing the partition tree to adapt to the structureof the underlying distributions and (ii) incorporating latent state variables thatcharacterize local features of the underlying distributions. However, important
limitations of the PT remain, including (i) the sensitivity in the posterior inference with respect to the choice of the partition tree, and (ii) the lack ofscalability with respect to dimensionality of the sample space. We consider amodeling strategy for PT models that incorporates a flexible prior on the partition tree along with latent states with Markov dependency. We introduce a
hybrid algorithm combining sequential Monte Carlo (SMC) and recursive message passing for posterior sampling that can scale up to 100 dimensions. While
our description of the algorithm assumes a single computer environment, it has
the potential to be implemented on distributed systems to further enhance thescalability. Moreover, we investigate the large sample properties of the tree
structures and latent states under the posterior model. We carry out extensive
numerical experiments in density estimation and two-group comparison, whichshow that flexible partitioning can substantially improve the performance of
PT models in both inference tasks. We demonstrate an application to a mass
cytometry data set with 19 dimensions and over 200,000 observations.

11Introduction

The Pólya tree (PT) (Freedman, 1963; Ferguson, 1974; Lavine, 1992) is a stochastic
process that generates random probability measures and is introduced as a prior for
Bayesian nonparametric inference. While the PT generalizes the Dirichlet process(DP) (Ferguson, 1973) as it yields the DP under certain hyperparameters (Ferguson,
1974), the statistical properties and practical applications of the PT are very different.
While the DP is most frequently used as a mixing distribution that induces clustering
structures, the PT is often adopted for directly modeling probability densities.
The PT is defined generatively on a recursive partition—or a partition tree—over
the sample space through coarse-to-fine sequential probability assignment at each tree
split. In a classical (univariate) PT, the tree is dyadic and the conditional probability
assigned to the two children nodes at each tree split arises from independent beta priors, which leads to analytic simplicity and ease in computing the posterior. Obtaining
the posterior is straightforward from beta-binomial conjugacy and incurs a computational budget that scales only linearly with the sample size, making the PT one of the
few nonparametric models applicable to data with massive sample size. Moreover,the posterior computation is embassingly parallelizable across the tree nodes.
The PT has been applied in various contexts beyond the original application ofdensity estimation. A far-from-exhaustive list includes survival analysis (Muliere andWalker, 1997; Neath, 2003), imputing missing values (Paddock, 2002), goodness-offit tests (Berger and Guglielmi, 2001), two-group comparison (Ma and Wong, 2011;
Holmes et al., 2015; Chen and Hanson, 2014; Soriano and Ma, 2017), density regression (Jara and Hanson, 2011), ANOVA (Ma and Soriano, 2018), testing independence(Filippi et al., 2017), and hierarchical modeling (Christensen and Ma, 2020). The PT
has also been utilized in semi-parametric analyses such as in (generalized) linear
models (Walker et al., 1999; Hanson and Johnson, 2002; Walker and Mallick, 1997).
2

Early developments of the PT are based on an a priori fixed partition tree on thesample space. The resulting inference can be sensitive to the choice of the partition
points defining the tree. In particular, the resulting process, both a priori and aposteriori, can be jumpy at these points. In hypothesis testing and model choice,this sensitivity is also reflected in the sometimes substantial change in the marginal
likelihood/Bayes factor when the partition points are slightly varied. To remedy theissue, Paddock et al. (2003) modified the PT model so that observations are generated
from the PT model with slightly different partition points. Hanson and Johnson(2002) and Hanson (2006) proposed a mixture of PTs by defining partition points
along quantiles of a parametric model endowed with a hyperprior to allow model
averaging on the partition points. This strategy does not allow individual partition
points to adapt to local features of the distribution but only the whole set of points
to the global structure of the distribution, and is most effective when the underlying
density is close to the specified parametric model. Nieto-Barajas and Müller (2012)
took a different approach by modeling the probability assignments within each level ofthe tree in a correlated manner to smooth out the random measure over the boundaries
of partitioning. While these approaches alleviate the sensitivity to partition points inlow-dimensional settings, they are not easily applicable (though in principle possible)
to problems with even just a handful of dimensions. Moreover, Bayesian inferencewith these models generally require MCMC, whose effectiveness can (in fact often
does) still suffer from the sensitivity with respect to the partition points.
Another related issue regarding the partitioning scheme of the PT is that in multivariate problems, traditionally the partition tree is constructed by dividing all dimensions of the sample space at each split. For example, for a d-dimensional sample
space, each time a tree node is divided, it is split into 2d children nodes, and probability assignment over these 2d child nodes is modeled by independent Dirichlet priors.

3

Wong and Ma (2010) noted that such a “symmetric” partition scheme is undesirable
as the dimensionality increases, in which case due to the exponential growth of thepartition blocks, the vast majority of the blocks are barely, if at all, populated by data.
As such they propose to incorporate adaptivity into the partitioning strategy withrespect to the structure of the underlying distribution through adopting a Bayesian
CART-like prior (Chipman et al., 1998) on the space of dyadic partition trees.
However, in order to maintain the analytic simplicity of the posterior and achieving MCMC-free exact Bayesian inference with a linear computational budget, theBayesian CART-like prior has to be restricted to only divide at the middle point (orotherwise a pre-determined fixed point) on one of the dimensions at each tree split.
Not only does this hamper the model’s ability to adapt to distributional structures,but it makes the model suffer from the same sensitivity with respect to the partition
points. Moreover, even with this restriction, the inference algorithm (based on recursive message passing) is only computationally practical for up to about 10 dimensions
on continuous sample spaces.
In a different vein, recent developments have demonstrated that aside from enhancing the partitioning strategy, the PT can also be substantially improved by adopting
more flexible priors (as opposed to independent betas) on the probability assignment
at each tree split (Jara and Hanson, 2011; Nieto-Barajas and Müller, 2012; Ma, 2017).
One strategy for enriching the PT in this regard is by introducing latent state variables at each tree split and adopt priors on the probability assignment given thesestates. When the latent states are discrete and modeled with Markov dependency,
analytical simplicity is preserved and exact Bayesian inference can proceed throughrecursive message passing that maintains the linear computational budget (Ma, 2017).
Given these developments, we are motivated by the following questions: Is it possible to incorporate into the PT a very flexible partition tree prior, such as the general4Bayesian CART (i.e., without the restriction to partition at middle points), that will(i) enhance its adaptivity to distributional structures in multivariate settings; (ii) resolve its sensitivity to the choice of partition points; and (iii) allow a tractable form of
the joint posterior and a posterior inference algorithm that is scalable to moderately
high-dimensional problems (e.g., up to 100 dimensions)? Moreover, should such a
strategy exist, can the resulting model and inference algorithm be made compatible
with incorporating (possibly Markov dependent) latent states on the tree nodes?
The goals of making the partition tree prior more flexible while enhancing thecomputational scalability appear at odds with each other. Large tree spaces are wellknown to be very hard to compute over. In moderate to high dimensional settingsexact inference involving flexible tree structures is beyond reach and even the mostadvanced MCMC approaches tailor-made for trees encounter substantial difficulty due
to the pervasive multi-modality of distributions in such spaces. Recent advances insequential Monte Carlo (SMC) for regression tree models (Lakshminarayanan et al.,2013; Lu et al., 2013), however, suggest that efficient inference is possible in moderately high-dimensional settings (up to about 100 dimensions). Moreover, once thepartition tree is sampled, the conditional posterior for the rest of the model can becomputed analytically through recursive message passing. We will therefore exploit a
hybrid strategy that uses a new SMC sampler to efficiently sample from the marginal
posterior of the partition tree structure, along with recursive message passing to compute the exact conditional posterior of the latent state variables given the tree.
Beyond the methodological development, we will also investigate the theoreticalproperties of the posterior on the partition tree and the latent states. Previous theoretical literature on the PT and related models have mostly focused on establishing
the posterior consistency and the contraction rate of the random measure induced under these models (Castillo, 2017a; Castillo and Randrianarisoa, 2021). In multivariate

5settings, however, the partition tree itself is highly informative about the underlying
distribution. Moreover, in applications involving model choice and hypothesis testing,it is often the latent states, not the random measures, that are of direct interest. As
such, we focus on studying the asymptotic behavior of the marginal posterior on thepartition tree and latent states, establishing consistency results on their convergence
toward the trees and states that most closely characterize the underlying truth.
The rest of the paper is organized as follows. In Section 2 we describe a flexible
prior on the partition tree structure that relaxes the restriction of “dividing in the
middle” on partition points and present a general form of PT models that adopt thisprior along with latent states associated with the tree nodes with a Markov dependency structure. In Section 3, we present our hybrid computational strategy that canwork effectively up to 100 dimensions consisting of an SMC algorithm for sampling on
the marginal posterior of the partition tree and a recursive message passing algorithm
for obtaining the exact conditional posterior of the latent states and the random measure given the sampled trees. In Section 4 we investigate the asymptotic propertiesof the tree structures and latent states identified under the posterior model. In Section 5, we carry out extensive numerical experiments to examine the performance ofour method in the context of two important applications of PTs—density estimationand the two-group comparison, followed by an application to a data set from mass
cytometry in Section 6. In Section 7 we conclude with a brief discussion. All proofs
are provided in Supplementary Materials C.

2

MethodWe first review the PT model (Ferguson, 1973; Lavine, 1992) on a dyadic recursive
partition in Section 2.1. The model, while defined on a general multivariate sample
space, differs from a traditional multivariate PT which adopts a multi-way symmetric6

recursive partitioning. Then we introduce a new class of PT models that incorporates
both the flexible partition prior and latent states with Markov dependency.

2.1

Pólya trees defined on recursive dyadic partitions

Without loss of generality, we consider a continuous sample space represented as a ddimensional rectangle Ω = (0, 1]d . For unbounded sample spaces such as Rd , one cantransform each margin to [0,1] by applying, say, a cumulative distribution function
transform or by standardizing the data based on its observed range of values. We
use µ to denote the Lebesgue measure on Ω. A (dyadic) recursive partitioning T on
Ω is a sequence of partitions of Ω such that the partition blocks at each level of thepartitioning are obtained by dividing each block in the previous level into two children
Skk
blocks. Formally, we can write T = ∞k=0 A , where A is a partition of Ω in the kthlevel. More specifically, A0 = {Ω}, and A ∈ Ak (k = 0, 1, 2, . . . ) is divided into Al
and Ar , which satisfy Al , Ar ∈ Ak+1 , Al ∪ Ar = A, and Al ∩ Ar = ∅. (Throughout the
paper, a subscript “l” or “r” on a node indicates the left or right child node.) For
example, when d = 1 and if the tree is recursively divided at the middle point of eachnode, then nodes in level k are of the form (l/2k , (l+1)/2k ] for some l ∈ {0, . . . , 2k −1}.
Another common strategy is to define the tree based on the quantiles of a probability)] for l ∈ 0, . . . , 2k −1.
measure F so that A ∈ Ak is of the form A = (F −1 ( 2lk ), F −1 ( l+12kGiven a partition tree T , we can define a random measure Q by putting a prioron the conditional probability θ(A) = Q(Al |A) = 1 − Q(Ar |A) at each A ∈ T .
Under the PT prior, the parameters θ(A) follow independent beta distributions
Beta(αl (A), αr (A)), where αl (A) and αr (A) are hyperparameters. The corresponding posterior, given an i.i.d. sample x1 , . . . , xn from Q, is again a PT with a simple
conjugate update on the conditional proabilities:
θ(A) | x1 , . . . , xn ∼ Beta(αl (A) + n(Al ), αr (A) + n(Ar )),7where n(A) represents the number of observations in a set A ⊂ Ω. Though the treeneeds to be infinitely deep to ensure full support of the PT, for practical purposes, one
typically sets a sufficiently large maximum depth (or resolution) of T and computethe posteriors of θ(A)’s defined on this finite tree structure (Hanson and Johnson,
2002). We shall refer to a node in the deepest level as a “leaf” or “terminal node”.
On a leaf, the conditional distribution can be set to a baseline F (·|A), such as theuniform distribution µ(·|A). In Section 3 when we present inference algorithms, weshall adopt this practical strategy and assume T is finite and use N (T ) and L(T ) todenote the collection of the non-terminal nodes and the leaf nodes, respectively.

2.2

Incorporating flexible partition points

We incorporate a Bayesian-CART like prior on T by randomizing both the dimension
in which to divide a node and the location to divide. Our prior relaxes the “alwaysdivide-in-the-middle” restriction imposed in Wong and Ma (2010). This prior on thepartition tree T differs from that in the mixture of PTs of Hanson (2006), which does
not randomize over the dimension to divide, but generates the boundaries of the tree
nodes jointly using quantiles of a parametric family.
Our prior can be described iteratively as a generative process that recursively
divides the sample space. Specifically, suppose we have a node A in the rectanglar
form, A = (a1 , b1 ]×· · ·×(ad , bd ]. We divide A into two rectangular children by cutting
along a randomly chosen dimension at a random location. The dimension to divide
D(A) ∈ {1, 2, . . . , d}, and the (relative) location to divide L(A) ∈ (0, 1) are given
independent priors of the following forms:
D(A) ∼ Mult(λ1 (A), . . . , λd (A)) and L(A) ∼

NX
L −1

βl (A)δl/NL (·),

(1)

l=1

where δx (·) represents the unit point mass at x, and NL − 1 is the total number ofgrid points along (0, 1). Both {λi (A)}i=1,...,d and {βl (A)}l=1,...,NL −1 sum to 1. In theabove, we have adopted a uniform grid over (0, 1) for notational simplicity, but it8

does not have to be as such. With D(A) = j and L(A) = l/NL , the two children
nodes Al and Ar are
Al = (a1 , b1 ] × · · · × (aj , aj + l/NL · (bj − aj )] × · · · × (ad , bd ],Ar = (a1 , b1 ] × · · · × (aj + l/NL · (bj − aj ), bj ] × · · · × (ad , bd ].
In principle one could adopt a continuous prior on the partition location L(A). A
discretized prior is helpful, however, because it will substantially simplify posterior
computation. In practice, as long as the grid is dense enough, the discrete priorwill be practically just as flexible. Indeed we have verified in extensive numericalexperiments that when NL is large enough (more than 30 to 50) over a uniform grid,posterior inference no longer improves in any noticeable way.
For the prior on D(A), we set λj (A) = 1/d for all nodes A as a default choice.
When L(A) is given a weak prior widely spread over (0, 1), the resulting inference
can be sensitive to the “tail” behavior of the distribution in the node, resulting inhigh posteriors of L(A) near the extreme values 0 and 1. A detailed discussion on
this phenomenon will be provided in Section 5.1.1. This issue can be effectively
addressed by making the prior of L(A) depend on the sample size n(A) so that it
encourages more balanced divisions at large sample sizes. More specifically, we adopt
the following prior with an exponentially decaying tailP (L(A) = l/NL ) = βl ∝ exp [−ηn(A)f (|l/NL − 0.5|)] , l = 1, . . . , NL − 1,(2)

where η ≥ 0 is a hyperparameter and f : [0, ∞) → [0, ∞) is an increasing function
with f (0) = 0. In the following, we shall use a function f (x) = x, and so our prioron L(A) is a (discretized) Laplace distribution. We provide theoretical justification
for adopting this prior with exponential tails in Section 4.
Another generalization of the prior on L(A) is to incorporate a spike-and-slab
set-up with a spike at the middle point 1/2. In particular, one can adopt a dependent spike prior among the nodes such that once a node A is divided exactly at the9

middle point, so are its descendants. This generalization will substantially reduce
the amount of computation in regions of the sample space where the data are either
sparse or lack interesting structure, e.g., close to the uniform distribution. We implement the spike-and-slab in our software but defer the details of this generalization to
Supplementary Materials A to avoid distracting the reader from the main ideas.
Given the tree prior, our PT model now consists of the two components—tree
generation and conditional probability assignment. Figures 1(a) and 1(b) present agraphical model representation for each.

(a) Partition tree generation(A0 : All potential nodes)(b) PT without latent states(c) PT with latent states

Figure 1: Graphical representation of PT models given the tree T . The hyperparameters are hidden for simplicity.

2.3

Hidden Markov Pólya tree models

2.3.1

General framework

Next we extend the above model to accommodate two recent developments in thePT literature: (i) incorporating latent state variables along the tree structure and (ii)joint modeling of multiple groups of observations. Incorporating latent variables allow
more flexibly characterizing distributional features through adding prior dependency.
As in recent literature, we consider incorporating discrete state variables that follow
a Markov process along the tree structure. Because the description in this section
always pertains to the model given the randomly generated partition tree T , forbrevity we shall not keep stating “given T ”.
We generalize our notation to allow observing one or more groups of i.i.d. observations. Let G be the number of groups of i.i.d. observations. For the gth group10(g = 1, 2, . . . , G), let Qg be the sampling measure for that group. Let Q denote the collection of all G sampling measures. That is, Q = {Qg }G
g=1 . Let xg = (xg,1 , . . . , xg,ng )denote the observations in the gth group, which are i.i.d. given Qg , where ng , thesample size for the group, is allowed to differ across the groups. We use x = {xg }G
g=1
to denote the collection of all observations from all groups.
Next we specify a prior on Q in terms of a joint prior on the conditional probability
on each A ∈ T , θg (A) = Qg (Al | A) = 1 − Qg (Ar | A). We use latent variable
modeling to incorporate prior dependency among the tree nodes. Specifically, let{V (A) : A ∈ T } denote a collection of latent state variables, one for each A, andwithout loss of generality, assume that V (A) takes discrete values from {1, . . . , I}.
(In practice, the number of states I can differ among A.) Joint priors of θg (A) for all
g and A are then defined conditionally on these latent states.
Existing literature has exploited these latent states to characterize both the withingroup structure of each distribution Qg and the between-group relationship among
the Qg . An example of within-sample structures is the smoothness of each underlying
distribution, which is explored in the context of density estimation (Ma, 2017). An
example of between-group structures is the difference between two (or more) distributions (Soriano and Ma, 2017).
Dependent modeling of the latent states over the partition tree is desirable asa priori one would expect interesting structures (both within-group and betweengroup) to exhibit themselves in a correlated manner over the sample space. For
example, functions tend to have similar smoothness over adjacent locations, and twogroup difference tend to be clustered in space. A computationally efficient strategy
for modeling such dependency over the tree is by a hidden Markov process along
the tree (Crouse and Baraniuk, 1997), which starts from the root node, A = Ω,and sequentially generates the latent states in a coarse-to-fine fashion according to11

(possibly node-specific) transition matrices ξ(A) whose (i, i0 )th element is
ξi,i0 (A) = P (V (A) = i0 | V (Ap ) = i) for i, i0 ∈ {1, . . . , I},
where Ap denotes A’s parent. (We shall use superscript “p” to indicate the parentof a node in T .) For A = Ω, since Ω has no parent, we can simply let ξi,i0 (Ω) be
constant over i, representing the initial state probabilities on Ω.
Given the V (A)’s, {θg (A)}g=1,...,G can then be modeled as conditionally independent a priori. The specific choices of these conditional priors are problem-dependent.
We give two examples below. Figure 1(c) presents a graphical model representation
for the latent state modeling on G probability distributions by PTs given T , whichalong with our generalized prior on the partition tree T presented in Figure 1(a) forms
the most general version of the model we consider in this work.
Example 1: Density estimation with adaptive smoothness
An example of within-group structures that the latent state V (A) can characterize is
the smoothness of the density functions for the random measures. For example, Ma(2017) proposed the adaptive Pólya tree (APT) model which incorporates latent statesto allow different levels of local smoothness in the underlying distribution. This isachieved by modeling the θ(A)’s as Beta(m(A)ν(A), (1 − m(A))ν(A)), where m(A) isthe prior mean and ν(A) the precision parameter which characterizes the smoothnessof the random measure with larger ν(A) corresponding to more smoothness, andthen modeling the precision parameters conditional on the latent state V (A) with ahyperprior ν(A) | V (A) = i ∼ Fi with the (conditional) hyperprior Fi given V (A) = i.
The hyperpriors F1 , F2 , . . . , FI are chosen to be stochastically increasing F1 ≺ F2 ≺· · · ≺ FI . That is, Fi ([x, ∞)) ≤ Fi+1 ([x, ∞)) for all x > 0. The rate at which thelatent states transition into a higher latent state along each subbranch of the partition
tree determines the local smoothness of the underlying density.

12

Example 2: Two-group comparison
In two-group comparison, we are interested in testing and identifying differences between two measures Q = {Qg }g=1,2 based on an i.i.d. sample from each. The “global”
testing problem can be formulated as testing the following null and alternative hypotheses: H0 : Q1 = Q2 vs H1 : Q1 6= Q2 . Noting that two-group differences may
exist in parts of the sample space and not others, the coupling OPT (Ma and Wong,2011) and the multi-resolution scanning (MRS) model (Soriano and Ma, 2017) are
PT-based models that allow the measures to differ on some nodes A ∈ T and notothers. This more “local” perspective on two-group comparison enables these models
to not only test for H0 vs H1 , but to identify regions on which the two measures differ.
To achieve this, these models incorporate state variables that characterize whether
the conditional probabilities on each A are equal:
V (A) = 1 ⇔ Q1 (Al | A) 6= Q2 (Al | A), V (A) = 2 ⇔ Q1 (Al | A) = Q2 (Al | A). (3)When V (A) = 1, the two corresponding conditional probabilities are given independent beta priors, whereas if V (A) = 2, they are tied and given a single beta prior.
Markov dependency among the states on different nodes are incorporated to induce
the desired spatial correlation of cross-group differences. Additional latent states
can be further incorporated to reflect more complex relationships between the distributions. In fact, the MRS also incorporates an additional state V (A) = 3, which
introduces the same coupled prior as V (A) = 2, but works as an absorbing state thatonce V (A) = 3, all descendants of A will remain in that state, corresponding to thecase that the conditional distributions Q1 (·|A), and Q2 (·|A) are completely equal.

3

Bayesian inference

In sum, the models we consider all share a common structure consisting of the following components: (i) the partition tree T defined by the dimension and location

13

variables D’s and L’s, which follow the priors given in Eq. (1); (ii) the latent statevariables V (A) given T which follow a Markov prior; (iii) the conditional probabilities along the given tree T , {θg (A)}G
g=1 , whose joint prior are specified independently
across the nodes on T given the latent states; and finally (iv) given the random measures Qg defined by T and θg (A)’s, we observe an i.i.d. sample xg from each Qg ,independently across g. Formally, we have the following full hierarchical model:
T | λ, η ∼ p(T | λ, η){V (A) : A ∈ T } | ξ, T ∼ Markov(ξ)
ind(θ1 (A), . . . , θG (A)) | V (A), T ∼ p(θ1 (A), . . . , θG (A) | V (A)) for A ∈ T
iidxg = (xg,1 , xg,2 , . . . , xg,ng ) | Qg ∼ Qg for g = 1, 2, . . . , G.
The key to Bayesian inference is the ability to either compute or sample from thejoint posterior (T, V, θ) given all data x = (x1 , . . . , xG ), where V and θ represent the
totality of all latent states and conditional probabilities given T respectively. Whilein some problems such as density estimation one may mainly be interested in justthe marginal posterior of the Qg ’s, in others such as two-group comparison where one
wants to characterize the between-group relationships among the distributions, thelatent states (along with T ), which characterizes such relationhips, are often of prime
interest. In multivariate and even high-dimensional problems, the tree structure T is
also of great interest as it sheds light on the underlying structures in the distributions.
To this end, we shall take advantage of recent developments in sequential Monte
Carlo (SMC) sampling for tree-based models (Lakshminarayanan et al., 2013; Lu
et al., 2013) and advances in message passing algorithms for PT models with Markov
dependency (Ma, 2017).

We introduce a hybrid algorithm that combines thesetwo computational strategies to effectively sample from the joint posterior in highdimensional spaces. Overall, the hybrid algorithm consists of two stages:
1. Sampling from the marginal posterior of the partition tree14

We design an SMC sampler—that is, a particle filter—to sample a collectionof tree structures T 1 , . . . , T M by growing each tree from coarse to fine scales.
It uses one-step look-ahead message passing to construct proposal distributions
for D(A) and L(A), one node at a time.
2. Computing the conditional posterior given the sampled trees
Given each tree sampled by the SMC, we analytically compute the exact conditional posteriors of V (A)’s and θ(A)’s using recursive message passing. (Thedetailed algorithm will be given in Section 3.2.)

3.1

SMC to sample from tree posterior

In the SMC stage to sample the trees, each particle stores a realized form of a finite
tree structure, and one node of each tree is divided at each step of the SMC sampling.
Suppose Tt is the finite tree obtained after dividing the sample space t times in aparticle, and for this tree we define the target distribution
πt (Tt ) = P (Tt | x) ∝ P (Tt )P (x | Tt ).
Here P (Tt ) is the joint prior of the variables D(A)’s and L(A)’s for the non-leaf nodes
of Tt , and P (x | Tt ) is the marginal likelihood given the tree Tt under the hierarchicalmodel, in which V and θ are integrated out. To sample from this target distribution,m
we sequentially construct a set of M particles {Ttm , Wtm }M
m=1 , where Tt is a realized

tree and Wtm is the associated importance weight for the mth particle. Examples ofgenerated trees are given in Figure 2, where the sample space has been divided threetimes, and in the next step, new partition boundaries will be added in the gray nodes.
Following Lakshminarayanan et al. (2013), we adopt in each step of the SMC
a breadth-first tree-growth strategy by dividing the “oldest active” leaf node—that
is, the one generated in the earliest step and is yet to be terminated in division.
Further division of a node is terminated once the number of observations in that15

node is below a pre-set threshold (e.g., 5 in our software implementation) to avoid
overfitting.

Otherwise a node is bisected along a boundary whose dimension andlocation are randomly drawn from a proposal distribution. For each particle, a finite
tree Tt is formed by a sequence of decisions {Js }ts=1 , where Js = (Ds , Ls ) correspondto all of the variables D(A) and L(A) at the sth step of the SMC.

Figure 2: An example of realized finite trees in the particle system obtained after the
step t = 3. The numbers in the squares indicate in which step the boundaries are
drawn. Among the current leaf nodes, the nodes colored gray are the oldest nodes
generated in the earliest step, so they are split in the next step.
As we will see in Proposition 3.1, the target distribution πt (Tt ) has a decomposition
πt (Tt ) = Ct πt−1 (Tt−1 )πt (Jt | Tt−1 )wt (Tt−1 ), where Ct is a constant independent of Tt ,πt (Jt | Tt−1 ) a conditional distribution on Jt given Tt−1 , and wt (Tt−1 ) a function of
Tt−1 . We will choose πt (Jt | Tt−1 ) as the proposal for Jt under which the corresponding
importance weight will simply be wt (Tt−1 ), independent of Jt .
More specifically, suppose at the current step t, we are to divide At ∈ Tt−1 , into
At,l and At,r with decision Jt . Let Mi (At | Jt ) be the marginal likelihood on the node
At under the decision Jt evaluated based on the observations in At . That is,#
Z "YG
Mi (At | Jt ) =θg (At )ng (At,l ) (1 − θg (At ))ng (At,r ) dP (θ1 (At ), . . . , θG (At )|V (At ) = i),g=1(4)
where ng (A) is the number of observations of the gth group included in A. To avoid
cumbersome notation, we suppress in our notation the dependency of Mi (At | Jt ) on
the observations x. For example, if the {θg (At )}G
g=1 follow independent beta priors
16written as Beta(αli (At ), αri (At )) given V (At ) = i, then the marginal likelihood has theQ
B(αil (A)+ng (At,l ),αir (A)+ng (At,r ))following expression Mi (At | Jt ) = G, where B(·, ·) is
g=1
B(αi (At ),αi (At ))lr

the beta function. Based on the values of Mi (At | Jt ), we can analytically compute theproposal and the importance weight using a general recursive algorithm, as describedin the following proposition.
Proposition 3.1. For every possible decision Jt and states i = 1, . . . , I, let ϕi (At )be a function defined recursively: PIξ1,i (Ω)Mi (Ω|Jt )j=1 ξ1,j (Ω)Mj (Ω|Jt )ϕi (At ) =PIϕj (Apt )ξj,i (At )Mi (At |Jt ) PI j=1PIp
k=1

j=1 ϕj (At )ξj,k (At )Mi (At |Jt )if At = Ω(5)Otherwise,

where Apt is At ’s parent node. Also, let h(Jt | At ) be a function of Jt defined as
( I
)
I
X
X
µ(At,l )−n(At,l ) µ(At,r )−n(At,r ), (6)h(Jt | At ) =ϕj (Apt )ξj,i (At ) Mi (At | Jt )−n(At )
µ(A
)
t
i=1
j=1
where n(A) denotes the total number of observations included in a node A. Then thetarget distribution πt (Tt ) can be expressed in terms of πt−1 (Tt−1 ) asπt (Tt ) = Ct πt−1 (Tt−1 )πt (Jt | Tt−1 )wt (Tt−1 ),where Ct is a constant andX
P (Jt )h(Jt | At ), wt (Tt−1 ) =P (jt )h(jt | At )πt (Jt | Tt−1 ) = P
jt P (jt )h(jt | At )jt
The summation over jt is taken over all possible decisions.
Corollary 3.1. Let h(Jt | At ) be the function defined in Proposition 3.1. Then the
proposal distribution πt (Dt | Tt−1 ) is given by
πt (Dt | Tt−1 ) = πt (Dt | Tt−1 )πt (Lt | Dt , Tt−1 ), where171. πt (Dt | Tt−1 ) is Mult(λ̃1 (At ), . . . , λ̃d (At )) withλ̃j (At ) ∝

NX
L −1

πt ((j, l/NL ) | Tt−1 ) ∝ λj (At )

l=1

NX
L −1

βl (At )h((j, l/NL ) | At ).

l=1

2. Given D(At ) = j, the conditional posterior of L(At ) isπt (Lt = l/NL | Dt = j, Tt−1 ) =

NX
L −1

β̃l (At )δl/NL (·),

l=1

for j = 1, 2, . . . , I and l = 1, . . . , NL − 1 withβ̃l (At ) ∝ β(At )h(j, l/NL | Tt−1 ).
We also have an analytical expression of the incremental weight:
wt (Tt−1 ) =

d NX
L −1
X
j=1

λj (At )βl (At )h((j, l/NL ) | At ).

l=1

Remark: The recursive function ϕi (At ) can be computed based on ϕi (Apt ) with thefixed computational cost. Hence, the optimal proposal πt (Jt | Tt−1 ) and the incremental weight wt (Tt−1 ), which are functions of h(Jt | At ), can be obtained at eachstep with constant computational cost with complexity O(I 2 NL d n(At )). As such, ourinference algorithm scales linearly in both the dimensionality and the sample size.
m M
m
}m=1
, Wt−1
We summarize the algorithm in updating the particle system from {Tt−1

to {Ttm , Wtm }M
m=1 below. All operations are repeated for m = 1, . . . , M .

1. Choosing the current nodem
From Tt−1, choose the the node generated in the earliest step among the current

leaf nodes that have not been terminated, which is denoted by At .
2. Obtaining the information of the parent node
Locate At ’s parent node, Apt , and fetch the values of ϕi (Apt ) for i = 1, . . . , I.
3. Computing the necessary quantities
For all possible Jt = (Dt , Lt ), compute Mi (A | Jt ) (i = 1, . . . , I) and h(Jt | At ).

18

4. Dividing the current nodeCompute the parameters λ̃j (At ) for j = 1, . . . , d and sample
Dtm ∼ Mult(λ̃1 (At ), . . . , λ̃d (At )).
Given Dtm , compute the parameters β̃l (At ) for l = 1, . . . , NL − 1 and sample
Lm
t ∼

NX
L −1

β̃l (At )δl/NL (·).

l=1
m
Divide the current node At with Jtm = (Dtm , Lmt ) to update the tree Tt .

5. Updating the importance weightmCompute the incremental weight wt (Tt−1) and update the importance weights

W m wt (T m )Wtm = PM t−1 m0 t−1 m0 .
m0 =1 Wt−1 wt (Tt−1 )PM
If the effective sample size 1/ m=1 (Wtm )2 is less than some prespecified threshold (e.g., M/10), resample the particles.
6. Computing the information on the current node for its descendants
Given Jtm , compute ϕi (At ) for i = 1, . . . , I.
In the algorithm, we stop dividing At if either the depth of At is equal to a pre-set
maximum resolution K (e.g., 15) or the number of observations in At is less thana pre-set threshold (e.g., 5). The SMC algorithm terminates when all the nodes of
all the particles have been stopped. The maximum resolution K controls the levelof local details that the model allows to infer, and larger values of K require morecomputational time. In a wide range of applications we have found that setting K to
beyond 15 to 20 leads to minimal changes in the resulting inference.
A common technique in SMC is to resample the particles according to the importance weights {Wtm }M
m=1 when the effective sample size of the particles drops below alevel. In sampling trees, however, the importance weights are affected by the choiceof nodes to divide in multiple steps, and so the standard resampling scheme can be19too “short-sighted” and often results in sacrificing promising trees prematurely. To
address this issue we follow the strategy proposed in Lu et al. (2013) by resamplingm κ
the particles according to weights amt ∝ (Wt ) for some κ ∈ (0, 1], and compute the

new importance weights proportional to Wtm /am
t . We generally recommend using amoderate choice of κ such as 0.5, which we have found satisfactory in a variety ofnumerical experiments, and will be our default choice in all of our later examples.

3.2

Posterior computation given sampled tree structures

The second stage of our inference strategy is to compute the posterior distributionsof the latent states V (A) and the conditional probabilities θg (A) given each sampled
tree. We shall focus on deriving a generic recipe for computing the posterior of thelatent states given the tree that works for all models under consideration. Given boththe tree and the latent states, the posterior of θg (A) boils down to the corresponding posterior of standard PT models on a dyadic tree, which is problem-specific asprovided in the literature on each such model.
Now suppose the SMC algorithm has produced a collection of finite trees {T m }M
m=1malong with the importance weights {W m }M
m=1 . Given each tree T , it is possible

to analytically calculate the exact posterior of {V (A)}A∈T m with recursive message
passing (a form of dynamic programming), which we describe below.
For A ∈ N (T m ), let φA (i) be the marginal likelihood on A given that V (A) = i.
(As before we suppress the dependency of the marginal likelihood on the data tosimplify notation.) Specifically,Z
GYYφA (i) = q(x | A)P (dq | V (A) = i), where q(x | A) =qg (z | A).

(7)

g=1 z∈xg (A)In Eq. (7), taking the integration with respect to P (dq | V (A) = i) is equivalent tointegrating out θg (A) as well as the θg (A0 ) and V (A0 ) terms for all descendants A0 of
A. Another useful quantity is the marginal likelihood on a node A given the state of20

its parent node V (Ap ) = i, which we denote as ΦA (i) and is given byQ x∈x(A) µ(x | A)
if A is a leaf node,ΦA (i) =P Ii0 =1 ξi,i0 (A)φA (i0 ) if A is a non-leaf node.

(8)

Note that the ΦA (i) and φA (i) terms are related byφA (i) = Mi (A | J(A))ΦAl (i)ΦAr (i),(9)

where Mi is the marginal likelihood defined in (4) given under the decision J(A) =(D(A), L(A)) to divide A into Al and Ar . By iteratively computing Eqs. (8) and (9)in a bottom-up fashion (i.e., starting from the leaves all the way to the root), we cancompute the pair {(φA (i), ΦA (i)) : A, i} for all nodes in the tree, and these pairs are
the “messages” passed along the tree from leaf to root.
Given the values of {(φA (i), ΦA (i)) : A, i}, we can now obtain the posteriorMarkov transition probability matrices of the latent states given the tree T , ξ̃(A) =(ξ̃i,i0 (A))I×I where the (i, i0 )th element ξ̃i,i0 (A) = P (V (Ω) = i0 | V (Ap ) = i, x, T ) andthe posterior marginal probabilities of the latent states given the tree T ,γ̃(A) = (γ̃i (A))Ii=1 = (P (V (Ω) = i | x, T ))Ii=1 .
Specifically, by the Bayes’ theorem, ξ̃(A) = D1−1 (A)ξ(A)D2 (A), where D1 (A) andD2 (A) are diagonal matrices with D1 (A)i,i = ΦA (i) and D2 (A)i,i = φA (i). Aftercomputing these transition matrices, we can compute γ̃(A) (the feedback “message”)
in the top-down manner (i.e., starting from the root and down to the leaves) asγ̃(Ω) = ξ̃1,· (Ω) and γ̃(A) = γ̃(Ap )ξ̃(A) for A 6= Ω.
We note that ΦΩ (1) computed in the recursive algorithm is the overall marginal
likelihood given the tree T , P (x | T ), which can be used to find the maximum
a posteriori (MAP) tree among the sampled trees, i.e., the sampled tree T m thatmaximizes P (T m | x) ∝ P (T m )P (x | T m ). We can use this “representative” tree,
along with the conditional posterior of the latent states given this tree, to visualize

21and summarize the posterior inference in an interpretable way.
Now that we have completely described our inference algorithm, we provide two
specific examples to demonstrate how one may use the output of the algorithm—
namely the sampled trees along with the conditional posterior given the trees—tocarry out inference. The first example is density estimation which involves learning
the within-sample distributional features while the second is two-group comparison
whose focus is on learning the between-sample structures. The inference strategies
for these quintessential examples are generalizable to a variety of other tasks.
Example 1: Density estimation
We consider the problem of estimating an unknown density based on a set of i.i.d.
observations from that density. This corresponds to G = 1 and so we can drop thesubscript g to simplify the notation. We shall use the posterior mean density, also
called the predictive density—E[q(·) | x]—as an estimate for the density q = dQ/dµ.
It can be computed by integrating out the random trees and the latent variables
based on the SMC sample. Specifically, for a given tree T , we define for all A ∈ T
and i ∈ {1, 2, . . . , I}, the following quantityeA (i) := E[Q(A)I{V (A)=i} | x].
All of these eA (i) terms can be computed together by a single top-down (i.e., rootto-leaf) recursion on the tree as given in the following proposition.
Proposition 3.2. For the root node, eΩ (i) = γ̃1,i (Ω). For a non-root node A, eA (i)can be computed recursively aseA (i0 ) =p

ϑ(A ) =

IX

ξ̃i,i0 (A)E[ϑ(Ap ) | V (Ap ) = i, T m , x] eAp (i), where

i=1



θ(Ap )

if A is the left child of Ap ,
1 − θ(Ap )if A is the right child of Ap .
22

Now with a recipe for obtaining all ei (A)’s for a given tree, we can obtain theconditional predictive measure given the tree T ,I
X µ(B ∩ A) X
eA (i) for any Borel set B.
E[Q(B) | x, T ] =µ(A) i=1
A∈L(T )

Now, given an SMC sample of M trees and weights, the posterior predictive density
can be computed by
E[q(x) | x] ≈

M
X

Wm

m=1

E[Q(B m (x)) | x, T m ],µ(B m (x))

where B m (x) ∈ L(T m ) the leaf node to which x belongs.
Example 2: Two-group comparison
If we are interested in comparing two groups of observations using generalizations to
the PT models described in Section 2.3.1, we shall compute the posterior probabilityof the two hypotheses H0 and H1 . For example, when V (A) is defined as in Eq (3),the posterior probability of the “global” null hypothesis H0 : Q1 = Q2 is given byX
P (H0 | x) =P (V (A) 6= 1 for all A ∈ N (T ) | T, x)P (T | x)
T ∈T

≈

M
X

W m P (V (A) 6= 1 for all A ∈ N (T m ) | T m , x),m=1

where the sum over T in the first row is over all finite trees with maximum resolution
K and the quantity P (V (A) 6= 1 for all A ∈ N (T m ) | T m , x) again is available
analytically by message passing (details given in Supplementary Materials B).
In addition to testing the existence of any difference between groups, it is usuallyof interest to detect where and how the underlying distributions differ. To this end,
we can compute the “posterior marginal alternative probability” (PMAP) on eachnode A, along any sampled tree T m :P (θ1 (A) 6= θ2 (A) | T m , x) = P (V (A) = 1 | T m , x) = γ̃i (A).

23

Reporting the PMAPs along a representative tree such as the MAP among the sampled trees can be a particularly useful visualizing tool to help understand the nature of
the underlying difference. One can also report on each A the estimated magnitude ofthe difference using a notion of “effect size” based on the log-odds ratio (Soriano andihihθ2 (A)
θ1 (A)−
log. In particular, one can report theMa, 2017), eff(A) = log 1−θ
1−θ2 (A)
1 (A)posterior expected effect size E[eff(A) | x, T ], which can be computed using a standard
Monte Carlo (not MCMC) sample from the exact posterior given the representativetree. We will demonstrate this using a mass cytometry data set in Section 6

4Theoretical Properties

Next we investigate the theoretical properties of the proposed model. Previous theoretical analysis on the PT had mostly focused on establishing the marginal posteriorconsistency and contraction of the random measures Qg with respect to an unknown
fixed truth (Walker and Hjort, 2001; Castillo, 2017b). We shall take a different perspective and instead investigate the asymptotic behavior of the marginal posteriorof the partition tree T and the latent states as these are often critical quantities ofpractical importance in applications. We note that once given the tree and the latentstates, our model reduces to standard PTs on a dyadic tree and thus the posteriorconsistency of the random measures Qg ’s will follow from previous results once weestablish the posterior consistency of the tree and the latent states. As the sample
size increases, the two key theoretical questions of interest here are:(1) What tree structures does the marginal posterior of T concentrate around?
(2) How does the posterior of the latent states given the tree behave?
These two questions have broad relevance in inference using PT models, andpreviously several authors have investigated the second question in the two-groupcomparison context for their variants of the PT model (Holmes et al., 2015; Soriano
24and Ma, 2017). In addressing the second question more generally, we aim to provide
results that encompass these previous analyses as special cases. According to ourlimited knowledge, we are not aware of previous studies on the first question.
We will address each of the two questions in turn. Throughout this section, we
consider finite PTs with maximum depth of the trees set to some value K. We use
T K to denote this collection of trees. Also, the asymptotic results are derived under
the prior for L(A) provided in Section 2.2 which can depend on the (finite) sample
size. The case of an uniform priors on L(A) independent of the sample size is included
as a special case where the hyperparameter η = 0. Finally, we consider models thatsatisfy Assumption 1 and Assumption 2 described below. The models discussed
in Section 2.3.1 all meet this requirement.
Assumption 1. For each group g ∈ {1, ..., G}, let ng be the sample size and Pg thetrue probability measure from which the observations are generated. We assume that(i) There exists ζg ∈ (0, 1] such that ζg = limn→∞ nng for g ∈ {1, . . . , G}, wheren = n1 + · · · + nG is the total number of observations across all groups.
(ii) The sampling distribution Pg satisfies Pg  µ, and the density pg = dPg /dµ is
positive almost everywhere.
Additionally, given the tree T and the latent states, the parameters {θg (A)}G
g=1 aregiven one of the following priors (the model can adopt a mix of these priors fordifferent combinations of A and V (A) values):
Prior A : θg (A) independently follow a beta prior.
Prior B : θ1 (A) = · · · = θG (A) and follow a beta prior.
Prior C : θ1 (A) = · · · = θG (A) ≡ c(A), some constant in (0, 1).
Establishing the theoretical properties also requires a condition on the latentstates. In particular, under some states, the support of the prior of the parameters {θg (A)}G
g=1 needs to include the true conditional probabilities. To describe this25

requirement, given a tree T ∈ T K , let Si (A | T ) be the support of the prior on(θ1 (A), . . . , θG (A)) under the state V (A) = i. Then, let τ (A | T ) denote the collection of “feasible states” on A. (A state is “feasible” if the true conditional probabilities
are in the support of the corresponding prior given the state.) That is,τ (A | T ) := {i ∈ {1, . . . , I} : (P1 (Al | A), . . . , PG (Al | A)) ∈ Si (A | T )}.
The next assumption states that the prior for the latent states must give positive
probability for all the latent states to all simultaneously be feasible.
Assumption 2. For every T ∈ T K , P (V (A) ∈ τ (A | T ) for all A ∈ T ) > 0.
With these assumptions, we next derive asymptotic properties for the marginal posteriors for the tree and the state variables. In the following, we use the notation xn
instead of x for the data to indicate the total sample size.
In order to describe the posterior convergence of the partition trees, we introduce
a notion for “tree-based approximation for probability measures”. Let T be a finite
tree and H a probability measure. Then the “tree-based approximation of H underP. for any B ∈ B(Ω).
T ”, denoted by H|T , is defined as H|T (B) = A∈L(T ) H(A) µ(B∩A)µ(A)
The following theorem then characterizes the trees the posterior concentrates on as
the sample size grows.
Theorem 4.1. Let TMK be the collection of trees under which the tree-based approximation of the measures Pg minimizes the Kullback-Leibler divergence from the Pg ’splus a penalty term on unbalanced splits. That is,
TMK = arg min
T ∈T K

where
Bg (T ) =GX

ζg {KL(Pg ||Pg |T ) + ηBg (T )} ,

g=1

X


Pg (A)fA∈N (T )µ(Al )
− 0.5
µ(A)


.

Then the marginal posterior of T concentrates on TMK . That is, as n → ∞,p

P (T ∈ TMK | xn ) →
− 1.
26

(10)

For the state variables, it is desirable that their posterior distribution concentrates
on a collection of feasible states. Moreover, when multiple configurations of the states
are feasible, it is desirable that the posterior concentrates around such configurations
that provide the most parsimonious representation of the true distributions. For
example, if the true conditional distribution on a node is uniform, a model that
introduces a possible non-uniform structure on this node is feasible but redundant.
White and Ghosal (2011) and Li and Ghosal (2014) showed that, in quite general
settings of multi-resolution inference, the posterior probability of such redundant
models tends to concentrate its mass on 0. By adapting their techniques, we showthat the same property holds in the case of our model.
To formally describe the results, we need to define the complexity of the model
specified by the latent states. Given the state V (A) = i, the complexity of the{θg (A)}G
g=1 , in other words, the number of free parameters of the prior distribution
under the ith state is denoted by Ci (A). For example, for two-group comparison,
1 if θ1 (A) = θ2 (A)Ci (A) =2 if θ1 (A) 6= θ2 (A).
Next we introduce the complexity of a combination of states on the tree T . Givena tree T , let V denote a combination of the state variables {V (A)}A∈N (T ) and let
v = {v(A)}A∈N (T ) (v(A) ∈ {1, . . . , I}) be one of the possible realizations of V. Then
we define the model complexity under v as follows:X
C(v) =
Cv(A) (A).

(11)

A∈N (T )The next theorem shows that the posterior distribution of the states given thetree will concentrate on those that are feasible and most parsimonious.
Theorem 4.2. For T ∈ T K , let VT = {v : v(A) ∈ τ (A | T ) for all A ∈ N (T )}.
p

Then P ({V ∈ VT } ∩ {C(V) = minv∈VT C(v)} | T, xn ) →− 1.
Remark: Consistency results for several existing models are special cases of this27

theorem. For example, we derive the consistency of the MRS model for two-groupcomparison as a corollary in Supplementary Materials C.

5

Experiments

In this section, we carry out simulation studies to examine the performance of ourmodel and inference algorithm. In particular, we are interested in (i) understanding
how the model with the flexible tree prior compares to those with a “divide-in-themiddle” restriction, and (ii) verifying the linear scalability of our inference algorithmwith respect to increasing dimensionality. We again consider the two quintessential
examples—(i) density estimation and (ii) the two-group comparison—for inferring
within-group and between-group structures respectively.
We shall consider both low-dimensional settings where the underlying structure
is easy to interpret and software for existing PT models are available, and highdimensional settings for which existing implementation of PT models is not applicable and we use our SMC algorithm to carry out inference for both our model and theearlier models with fixed partitioning points (which are special cases of our model).
Throughout the experiments, the parameters NL and M are fixed to 32 and 1000
respectively. We note that larger NL values can also be adopted at a linear computational cost but did not lead to noticeable change in performance in our examples.
Details such as the settings of hyper-parameters and simulated data sets are provided
in Supplementary Materials F unless explicitly described in this section.

5.1

Density estimation

We first consider 2D examples to observe what kind of tree structures are obtained
under the flexible model and how prior specification in Eq. (2) influences the performance. After that, we move to higher dimensional cases to examine the scalability of
our new SMC method and the effect of incorporating the flexible partition. For this28

task we compare our model with the APT model (Ma, 2017) which also incorporates
a prior on the dimension to divide but restricts partitioning at middle points. Itsposterior computation is implemented by the apt function in the R package PTT.
5.1.1

Two-dimensional cases

Simulated data are generated from the three scenarios with the densities visualized
in the first row of Figure 3. (Details on the simulation settings are provided Supplementary Materials F.1.2.) Also presented in Figure 3 are examples of the posteriormean densities E[q | x] as well as the partition blocks under the MAP tree. Note thatthe posterior mean is computed by integrating out the unknown tree, and the MAP
tree is presented to visualize key distributional features. The results for the firstscenario confirms that our more flexible model is much more effective in capturing
the discontinuous boundaries of the true density. For the second scenario, our modeltends to draw the boundaries that surround the true clusters. In the trees given
under the different values of η , however, we can see that fewer nodes were divided
inside the clusters when η = 0.01. In contrast, when η = 0.1, the representative treedraws outlines of the clusters and divides regions inside of the clusters at the sametime. A similar phenomenon is observed in the third scenario—under our model withflexible partitioning points, partition lines are formed around the region with highdensity, when η = 0.1 for the boundaries were also drawn within the high probability region. The quantitative comparison based on the KL divergence is provided in
Supplementary Materials G, which is consistent with the explanation above.
5.1.2

Higher-dimensional cases

We generate d-dimensional i.i.d. observations from a density with independent pairsQof margins, i.e., f (x1 , x2 , . . . , xd ) = d/2
j=1 fj (x2j−1 , x2j ) where29Figure 3: The posterior means of the densities and the representative trees obtained
under n = 1000. Each column corresponds to a simulation scenario. The first row
shows that true densities, the second row corresponds to the APT model (with fixed
partition), and the third and fourth rows correspond to our model with flexible partitioning with parameters η = 0.01 and 0.1 respectively.

30

fj (x2j−1 , x2j ) =pj Beta(x2j−1 | 0.25, 1) × Beta(x2j | 0.25, 1)+(1 − pj )Beta(x2j−1 | 50/j, 50/j) × Beta(x2j | 50/j, 50/j),with pj = 0.25 + 0.7/j. We consider two different situations: (i) the dimensiond = 6, and the sample size n changes from 5,000 to 50,000; and (ii) the sample sizen = 10, 000 and the dimensionality changes from 10 to 100. For our method, themaximum depth K is set to 15. For the algorithm of the original APT (implemented
by the apt function in the PTT package), in the first case with d = 6, the maximum
resolution is fixed to 9 because setting higher values leads to insufficient memory.
Also, because the apt function does not scale if the dimension is beyond d ≈ 10,
in the second case with large d, we used the proposed SMC algorithm to carry out
inference for the original APT model, which corresponds to setting NL = 2, namely,
fixing all the partition points to the middle. In addition to the original APT, we show
the comparison with the classical PT method in which each node is split into 2d andthe maximum depth is 15, the Dirichlet process Gaussian mixture model (Escobar
and West, 1995; Müller et al., 1996) implemented by the PYdensity function in the
R package BNPmix (Corradin et al., 2021), and for (i) the Gaussian kernel density
estimation using the kde function in the R package ks (Duong et al., 2007).

Figure 4: The wall time under five different data sets. The flexible model with
η = 0.01 is used. The black dashed lines indicate the average times.
Figure 4 presents the computational time for five different data sets. To obtain the31result, we used a singe-core environment using Intel Xeon Gold 6154 (3.00 GHz) CPU.
The computational time is linear in both the sample size and the dimensionality.

Figure 5: Predictive performance of 6 methods. Each point corresponds to the averageof the predictive score in Eq. (12) based on 50 data sets. Each interval is formed byadding and subtracting the standard deviation. In the right plot, the predictive scoresof the DPM model for the over 60 dimensional cases are below the displayed range.
Because in the high-dimensional settings we cannot obtain the KL divergence
between the estimated density and the true density through numerical approximation,
we compare the models based on predictive scores, which is in essence an empiricalestimate of the KL divergence. Specifically, for each training set x, we can generatea new test set denoted by x∗ = (x∗1 , . . . , x∗n ) from the same true model and computen1X
log p̂(x∗i | x),
(12)
n i=1
where p̂ is the posterior mean of the density E[q | x]. We repeat this computation for
50 test/training set pairs and take the average. The results, given in Figure 5, showthat our model substantially outperforms the competitors by this criteria both whend = 6 with varying sample size and when n is fixed with varying dimensionality. We
also investigate the performance under small sample sizes, and the results are similar.
(See Supplementary Materials G for details.)

32

5.2

Two-group comparison

Next we consider the two-group comparison problem, evaluate the performance of ourmodel, and compare it to the original MRS with the “divide-in-the-middle” restriction. We use three scenarios (“Local location shift”, “Local dispersion difference”,
and “Correlation”) to generate 50-dimensional data sets. (Details of the scenarios
are provided in Supplementary Materials F.2.2.) The first two scenarios involve
two-group difference that lies in only parts of the sample space, or “local” differences.
which will help demonstrate the usefulness of inferring the partition tree in identifying
the nature of the differences. The sample size is n1 = n2 = 2, 000 in all scenarios.
The original algorithm for inference under the MRS model by message passing,which is implemented by the mrs function in the R package MRS, is not scalable beyond
about 10 dimensions even with fixed partition locations. Hence we compute theposterior for both our model and the original MRS in all scenarios with our SMC andmessage passing hybrid algorithm. For the two models, the maximum resolution K
is fixed to 15. We compare the performance using receiver operating characteristic(ROC) curves computed based on 200 simulated data sets under each scenario.

Figure 6: The ROC curves for the 50-dimensional examples.
Figure 6 presents the ROC curves. For the location shift and dispersion differences,the model with flexible partitioning results in substantially higher sensitivity. For the33

Figure 7: The node with the highest PMAP P (V (A) = 1 | x) under the three scenarios for the 50-dimensional example, estimated by the MRS with flexible partitioningand η = 0.1. The solid lines mark the boundaries of the nodes and the partition linethat divides them into the two children nodes. The red triangle points and the blue
circle points are the observations of the two groups included in the node. Gray points
indicate the observations outside the node.
correlation scenario, the model with fixed partitioning locations performed slightly
better. This is not surprising since in this scenario the difference exists smoothly overentire ranges of the dimensions without natural “optimal” division points, and so theperformance gap is the cost for searching over more possible partition locations, noneof which improves the model fit than the middle point. It is worth noting again thatwhile the model with fixed partitioning performs well here, it is only with our newcomputational algorithm that it can be fit to data of such dimensionality.
To demonstrate the posterior model can help understand the nature of the differences, we present under each scenario the node with the highest PMAP, or P (V (A) =1 | x) = P (θ1 (A) 6= θ2 (A) | x), in Figure 7. In the location shift and dispersion difference scenarios the boundaries are away from the middle point to characterize thedifference, which partly explains the sensitivity gain in adopting the flexible tree prior.

6

Application to a mass cytometry data set

Finally, we apply our model for two-group comparison to a mass cytometry data
set collected by Kleinsteuber et al. (2016). The data set records 19 different mea-

34

surements including physical measurements and biomarkers on single cells in blood
samples from a group of HIV patients as well as in reference samples from healthy
donors. For demonstration, we compare the sample from an individual patient sample(Patient #1) and to that from a healthy donor to identify differences in immune cell
profiles from these samples. The sample sizes are 29, 226 for the health donor andand 228, 498 for the patient, with each observation corresponding to a cell. We set
η = 0.1 and the maximum depth K to 25.
Given the large sample sizes, the posterior probability for the global alternative
P (Q1 6= Q2 | x) is almost 1 and so is of less interest. Our focus is instead on
identifying the cell subsets on which the samples differ and on quantifying such differences. To this end, we identify a representative tree and report the “effect size”(i.e., the posterior expected log-odds ratio between the two samples) on each node ina representative tree—the MAP among the sampled trees.
The estimated eff(A)’s on the MAP tree up to level 9 is visualized in Figure 8, andthe full tree is provided in Supplementary Materials G. We note that the nodes
on which there is significant evidence for two-group differences, as well as those withlarge estimated effect sizes tend to be nested or clustered in subbranches of the tree,
which is consistent with our intuition that there is spatial correlation in the two-groupdifferences, and justifies the hidden Markov structure embedded in the MRS model.
Figure 14 in Supplementary Materials G presents the 20 nodes with the largest
values of estimated eff(A). In this figure, many of the nodes are in very deep levelsof the tree. We adopted a spike-and-slab with higher spike probability in very deeptree levels to further speed up the computation (details given in Supplementary
Materials A) and that explains why many of the very deep, small nodes plottedhave partition lines in the middle under the MAP tree.

35Figure 8: The MAP tree for the mass cytometry data set visualized up to the 9th
level. The size and the color indicate the estimated eff(A), and the number above anode indicates dimension in which it is split. Only the nodes with the sample size
larger than 50 are drawn.

7

Concluding Remarks

We have proposed a general framework for the PT model that incorporates a flexible
prior on the partition tree and can accommodate latent state variables with Markov
dependency along the partition tree. We have proposed a sampling algorithm thatcombines SMC and recursive message passing that can scale up to moderately highdimensional (∼ 100-dim) problems. Our numerical experiments confirm that oursampling algorithm scales linearly in the sample size and the flexible partitioning
tree prior can result in substantial gain in performance in some settings. Though
we have mainly used two inference tasks—namely density estimation and two-groupcomparison—to demonstrate our model and algorithm, our approach can be readily
applied to other PT models with a hidden Markov structure.
One notable limitation of our model—and in fact all CART-like models—is that
we only consider trees in which the node boundaries are all parallel to the axes. This
could lead to inefficiency in inference. For example, when there is a strong correlation
between several variables, drawing boundaries slanted according to the correlation
structure would be more effective in characterizing the underlying distribution. Such36

trees will need to be represented by more than just the D(A) and L(A) used in ourmodel, and how to efficiently compute the posterior distribution is an open problem.
Finally, our proposed algorithm is currently designed to be run in a single computer environment, so though the computational cost is linear to the sample size n,
dealing with huge data with huge n (e.g., > 1012 ) is not yet realistic. It is of future interest to develop parallelized versions of the algorithm for distributed systems, which
could explore either the parallel structure over nodes or parallel SMC algorithms.

Software
An R package for our method is available at https://github.com/MaStatLab/SMCMP.

Acknowledgment
LM’s research is partly supported by NSF grants DMS-2013930 and DMS-1749789.
NA is partly supported by a fellowship from the Nakajima Foundation.

References
Berger, J. O. and Guglielmi, A. (2001). Bayesian and conditional frequentist testingof a parametric model versus nonparametric alternatives. Journal of the American
Statistical Association, 96(453):174–184.
Castillo, I. (2017a). Pólya tree posterior distributions on densities. In Annales del’Institut Henri Poincaré, Probabilités et Statistiques, volume 53, pages 2074–2102.
Institut Henri Poincaré.
Castillo, I. (2017b). Pólya tree posterior distributions on densities. In Annales del’Institut Henri Poincaré, Probabilités et Statistiques, volume 53, pages 2074–2102.
Institut Henri Poincaré.

37

Castillo, I. and Randrianarisoa, T. (2021). Optional Pólya trees: posterior rates anduncertainty quantification. arXiv preprint arXiv:2110.05265.
Chen, Y. and Hanson, T. E. (2014). Bayesian nonparametric k-sample tests for
censored and uncensored data. Computational Statistics & Data Analysis, 71:335–
346.
Chipman, H. A., George, E. I., and McCulloch, R. E. (1998). Bayesian cart model
search. Journal of the American Statistical Association, 93(443):935–948.
Christensen, J. and Ma, L. (2020). A Bayesian hierarchical model for related densities
by using Pólya trees. Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 82(1):127–153.
Corradin, R., Canale, A., and Nipoti, B. (2021). Bnpmix: an r package for bayesiannonparametric modelling via pitman-yor mixtures. Journal of Statistical Software.
Crouse, M. S. and Baraniuk, R. G. (1997). Contextual hidden markov models forwavelet-domain signal processing. In Conference Record of the Thirty-First Asilomar Conference on Signals, Systems and Computers, volume 1, pages 95–100.
Duong, T. et al. (2007). ks: Kernel density estimation and kernel discriminant analysis
for multivariate data in r. Journal of Statistical Software, 21(7):1–16.
Escobar, M. D. and West, M. (1995). Bayesian density estimation and inference using
mixtures. Journal of the american statistical association, 90(430):577–588.
Ferguson, T. S. (1973). A Bayesian analysis of some nonparametric problems. Annalsof statistics, pages 209–230.
Ferguson, T. S. (1974). Prior distributions on spaces of probability measures. Annalsof statistics, 2(4):615–629.
Filippi, S., Holmes, C. C., et al. (2017). A Bayesian nonparametric approach to testingfor dependence between random variables. Bayesian Analysis, 12(4):919–938.
Freedman, D. A. (1963). On the asymptotic behavior of bayes’ estimates in thediscrete case. Annals of Mathematical Statistics, pages 1386–1403.

38

Hanson, T. and Johnson, W. O. (2002). Modeling regression error with a mixture ofPolya trees. Journal of the American Statistical Association, 97(460):1020–1033.
Hanson, T. E. (2006). Inference for mixtures of finite Polya tree models. Journal of
the American Statistical Association, 101(476):1548–1565.
Holmes, C. C., Caron, F., Griffin, J. E., Stephens, D. A., et al. (2015). Two-sample
Bayesian nonparametric hypothesis testing. Bayesian Analysis, 10(2):297–320.
Jara, A. and Hanson, T. E. (2011). A class of mixtures of dependent tail-free processes.
Biometrika, 98(3):553–566.
Kleinsteuber, K., Corleis, B., Rashidi, N., Nchinda, N., Lisanti, A., Cho, J. L., Medoff,
B. D., Kwon, D., and Walker, B. D. (2016). Standardization and quality control
for high-dimensional mass cytometry studies of human samples. Cytometry PartA, 89(10):903–913.
Lakshminarayanan, B., Roy, D., and Teh, Y. W. (2013). Top-down particle filteringfor Bayesian decision trees. In International Conference on Machine Learning,
pages 280–288.
Lavine, M. (1992). Some aspects of Polya tree distributions for statistical modelling.
Annals of statistics, 20(3):1222–1235.
Li, M. and Ghosal, S. (2014). Bayesian multiscale smoothing of gaussian noised
images. Bayesian analysis, 9(3):733–758.
Lu, L., Jiang, H., and Wong, W. H. (2013). Multivariate density estimation by
Bayesian sequential partitioning. Journal of the American Statistical Association,
108(504):1402–1410.
Ma, L. (2017). Adaptive shrinkage in Pólya tree type models. Bayesian Analysis,
12(3):779–805.
Ma, L. and Soriano, J. (2018). Analysis of distributional variation through graphical multi-scale beta-binomial models. Journal of Computational and Graphical
Statistics, 27(3):529–541.

39

Ma, L. and Wong, W. H. (2011). Coupling optional Pólya trees and the two sample
problem. Journal of the American Statistical Association, 106(496):1553–1565.
Muliere, P. and Walker, S. (1997). A Bayesian non-parametric approach to survival
analysis using Polya trees. Scandinavian Journal of Statistics, 24(3):331–340.
Müller, P., Erkanli, A., and West, M. (1996). Bayesian curve fitting using multivariate
normal mixtures. Biometrika, 83(1):67–79.
Neath, A. A. (2003). Polya tree distributions for statistical modeling of censoreddata. Advances in Decision Sciences, 7(3):175–186.
Nieto-Barajas, L. E. and Müller, P. (2012). Rubbery Polya tree. Scandinavian Journalof Statistics, 39(1):166–184.
Paddock, S. M. (2002). Bayesian nonparametric multiple imputation of partiallyobserved data with ignorable nonresponse. Biometrika, 89(3):529–538.
Paddock, S. M., Ruggeri, F., Lavine, M., and West, M. (2003). Randomized Polyatree models for nonparametric Bayesian inference. Statistica Sinica, pages 443–460.
Schwarz, G. (1978). Estimating the dimension of a model. Annals of statistics,
6(2):461–464.
Soriano, J. and Ma, L. (2017). Probabilistic multi-resolution scanning for two-sample
differences. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 79(2):547–572.
Walker, S. and Hjort, N. L. (2001). On Bayesian consistency. Journal of the Royal
Statistical Society: Series B (Statistical Methodology), 63(4):811–821.
Walker, S. G., Damien, P., Laud, P. W., and Smith, A. F. (1999). Bayesian nonparametric inference for random distributions and related functions. Journal of theRoyal Statistical Society: Series B (Statistical Methodology), 61(3):485–527.
Walker, S. G. and Mallick, B. K. (1997). Hierarchical generalized linear models andfrailty models with Bayesian nonparametric mixing. Journal of the Royal StatisticalSociety: Series B (Statistical Methodology), 59(4):845–860.
White, J. T. and Ghosal, S. (2011). Bayesian smoothing of photon-limited images40

with applications in astronomy. Journal of the Royal Statistical Society: Series B(Statistical Methodology), 73(4):579–599.
Wilks, S. S. (1938). The large-sample distribution of the likelihood ratio for testing
composite hypotheses. Annals of Mathematical Statistics, 9(1):60–62.
Wong, W. H. and Ma, L. (2010). Optional Pólya tree and Bayesian inference. Annalsof Statistics, 38(3):1433–1459.

41

A

Spike-and-slab type prior for L(A)A.1

Introducing an auxiliary variable

The location variable L(A) follows a spike-and-slab type prior which is expressed with
an auxiliary variable R(A) as

L(A) | R(A) ∼ 1{R(A)=1} δ1/2 (·) + 1{R(A)=0}

NX
L −1

β̂l (A)δl/NL (·),

l=1,l/NL 6=1/2

where 1 is the indicator function and the sum of the parameters βl (A) is 1. Under thisprior, L(A) follows the prior degenerated at 1/2 if R(A) = 1 and otherwise follows the
distribution on grid points other than the middle point. R(A) follows an asymmetrichidden Markov process

P (R(A) = 1 | R(Ap ) = 0) = r(A),P (R(A) = 1 | R(Ap ) = 1) = 1,

where r(A) ∈ [0, 1]. R(A) = 1 is the absorbing state, so once A is divided at the
middle point, L(A0 ) = 1/2 for every A’s descendant node A0 . In the estimation we
especially set the parameters as follows:
βl (A)
r(A) = βNL /2 , β̂l = PNL −1

l=1,l/NL 6=1/2 βl (A),where {βl (A)}l=1,...,NL is given in (2). Under this setting the prior of L(A) satisfies
P (L(A) = l | R(Ap ) = 0) = βl (A), l = 1, . . . , NL − 1.

42

Hence, L(A) follows the same prior as defined in (1) unless A’s parent node is divided
at the middle point, so the spike-and-slab prior can be seen as a natural extension.

A.2

SMC algorithm

In the SMC algorithm, we sample values of R(A) in addition to D(A) and L(A).
If R(A) = 1, which is equivalent to L(A) = 1/2, is sampled, we conclude there isno interesting structure on the node A so fix L to 1/2 for all the subsequent nodes.
Hence, we need to generalize the SMC algorithm discussed in Section 3.1 to sample
from the joint posterior distribution of the finite trees and the auxiliary variables R.
To describe this joint posterior, let Tt denote the finite tree structure, which isdetermined by the sequence of decisions J1:t dividing the nodes A1:t , and let R1:t bea sequence of the re-fixing variables for A1:t . Then the target distribution we want tosample from in the SMC is defined as

πt (Tt , R1:t ) = P (Tt , R1:t | x) ∝ P (Tt , R1:t )P (x | Tt , R1:t ).

The prior P (Tt , R1:t ) = P (J1:t , R1:t ) have a Markov chain structure on the tree, andits transition probability is decomposed as

P (Jt , Rt | Jtp , Rtp ) = P (Rt | Rtp )P (Jt ).
where Jtp = J(Apt ) and Rtp = R(Apt ) (Apt is the parent node of At ). On the other hand,because R1:t are conditionally independent of the observations given Tt , the likelihood
only depends on Tt as follows:

P (x | Tt , R1:t ) = P (x | Tt ).

43

The likelihood has the same form as in the original case without the auxiliary variables
R’s. Thus we can obtain the following proposition as a generalization of Proposition
3.1.
Proposition A.1. Let h(Jt | At ) be a function of Jt defined ash(Jt | At ) =

( I
I
X
X
i=1

)
ϕj (Apt )ξj,i (At )Mi (At | Jt )j=1µ(x(At,l ) | At,l )µ(x(At,r ) | At,r ).
µ(x(At ) | At )Then the target distribution πt (Tt , R1:t ) is expressed with πt (Tt−1 , R1:t−1 ) asπt (Tt , Rt ) = Cπt (Tt−1 , R1:t )πt (Rt | Tt−1 , Rtp )πt (Jt | Tt−1 , Rt )wt (Tt−1 , Rtp ),where C is a constant and
P (Rt | Rtp )P
P (jt | Rt )h(jt | At )njtPo,πt (Rt | Tt−1 , Rtp ) = P
p
P
(R
=
i
|
R
)
P
(j
|
R
=
i)h(j
|
A
)
t
t
t
t
t
t
jt
i=0,1
P (Jt | Rt )h(Jt | At ),πt (Jt | Tt−1 , Rt ) = P
jt P (jt | Rt )h(jt | At ))(X
X
P (jt | Rt = i)h(jt | At ) .
wt (Tt−1 , Rtp ) =P (Rt = i | Rtp )
i=0,1

jt

The summation with jt is taken over all possible decisions.
Its proof is essentially the same as Proposition 3.1 and Corollary 3.1 so it is omitted
in this material. The conditional posteriors πt (Rt | Rtp , Tt−1 ) and πt (Jt | Rt , Tt−1 ) are
analytically obtained as follows. First, if Rtp = 0, πt (Rt | Rtp , Tt−1 ) is Bernoulli(r̃(At )),44

wherer̃(At ) = r(At )JX

λj (At )h((j, 1/2) | At )j=1"
× r(At )JX

λj (At )h((j, 1/2) | At ) + (1 − r(At ))

j=1

J NX
L −1
X
j=1

#−1
λj (At )βl (At )h((j, l/NL ) | At )l=1

If Rtp = 1, then Rt is fixed to 1. Second, if Rt = 0, the posterior of Dt and Lt is
the same distribution given in Section 3.1. On the other hand, if Rt = 1, Lt is fixed
to 1/2, and the conditional posterior πt (Jt | Lt , Rt , Tt−1 ) is Mult(λ̂1 (At ), . . . , λ̂d (At )),where

λ̂j (At ) ∝ λj (At )h(j, 1/2 | At ).
After sampling (Rt , Jt ), the incremental weight wt (Tt−1 , Rtp ) (Rtp = 0, 1) is computedaswt (Tt−1 , 0)
= r(At )JX

λj (At )h((j, 1/2) | At ) + (1 − r(At ))j=1

wt (Tt−1 , 1) =

J NX
L −1
X
j=1

J
X

λj (At )βl (At )h((j, l/NL ) | At ),l=1

λj (At )h((j, 1/2) | At ),j=1

with which we update the importance weight Wt as Wt ∝ Wt−1 wt (Tt−1 , Rtp ).

mm M
The procedure to update the particle system {Tt−1, Wt−1}m=1 to obtain {Ttm , Wtm }M
m=1

is described in the following algorithm. The operations involving the index m is repeated for m = 1, . . . , M .

45,

1. Choosing the current nodem
, choose the oldest note from the leaf nodes, which is denoted by At .
From Tt−1

2. Obtaining the information of the parent node
Pick up At ’s parent node, which is denoted by Apt , and load the values of ϕi (Apt )for i = 1, . . . , I and Rtm,p = R(Apt ).
3. Computing the necessary quantities
If Rtm,p = 0, compute Mi (A | j, l/NL ) (i = 1, . . . , I) and h(j, l/NL | At ) for
j = 1, . . . , d and l = 1, . . . , NL − 1.
If Rtm,p = 1, compute Mi (A | j, 1/2) (i = 1, . . . , I) and h(j, 1/2 | At ) for
j = 1, . . . , d.
4. Deciding whether to fix the partition or not
If Rtm,p = 0, compute the parameter r̃(At ) and draw Rtm ∼ Bernoulli(r̃(At )).
If Rtm,p = 1, set Rtm to be 1.
5. Dividing the current node
Sample Jtm = (Dtm , Lmt ) as follows:• If Rtm = 0, compute the parameters λ̃j (At ) for j = 1, . . . , d and sample
Dtm ∼ Mult(λ̃1 (At ), . . . , λ̃d (At )).
Given Dtm , compute the parameters β̃l (At ) for l = 1, . . . , NL −1 and sample

Lm
t ∼

NX
L −1
l=1

46

β̃l (At )δl/NL (·).

• If Rtm = 1, compute the parameters λ̂j (At ) for j = 1, . . . , d and sample
Dtm ∼ Mult(λ̂1 (At ), . . . , λ̂d (At )),and set Lm
t = 1/2.
m
Divide the current node At with Jtm = (Dtm , Lmt ) to update the tree Tt .

6. Storing the information of the state’s posterior
Given Jtm , compute ϕi (At ) for i = 1, . . . , I and store them to the memory.
7. Updating the importance weightm, Rtm,p ) and update the importanceCompute the incremental weight wt (Tt−1weights as
W m wt (T m , Rm,p )Wtm = PM t−1 m0 t−1 m0 t m0 ,p .
)m0 =1 Wt−1 wt (Tt−1 , Rt
If the effective sample size 1/

PMm 2
m=1 (Wt ) is less than some prespecified thresh-old (M/10, say), resample the particles.

B

Additional algorithm for the MRS model

To describe the algorithm proposed in Soriano and Ma (2017), we keep using the same
notations in Section 3.2. Given the tree structure T , we compute functions ψ̃(A) for47

A ∈ N (T ) in the bottom-up (from the leaf nodes to the root node) manner as follows:ψ̃(A) =


ξ̃2,2 (A) + ξ̃2,3 (A)

if A ∈ L(T ),ξ̃2,2 (A)ψ̃(Al )ψ̃(Ar ) + ξ̃2,3 (A)






ξ̃1,2 (A)ψ̃(Al )ψ̃(Ar ) + ξ̃1,3 (A)if A ∈ N (T ) \ {Ω},if A = Ω.

Recall that only the first row of ξ(Ω) is meaningful as the initial distribution. Then
we obtain ψ̃(Ω) = P (H0 | T, x).

C

ProofsC.1

Posterior computation

Lemma C.1. For the finite tree Tt , let As be a node whose children nodes are leaf
nodes. Then we haveπt (Vs = i | Tt ) =πt (Tt , Vs = i)
= ϕi (As ),πt (Tt )where ϕi (As ) is defined in (5).
(Proof) Suppose that As belongs to the kth layer of Tt . Then there is a sub-sequence{ρ(l)}kl=1 such that Aρ(l) belongs to the lth layer andΩ = Aρ(1) ⊃ Aρ(2) ⊃ · · · ⊃ Aρ(k) = As .

By the definition of πt (Tt , V1:t ), for a sequence {vl }kl=1 such that vl ∈ {1, . . . , I}, we48obtain the expression of the conditional posterior of {Vρ(l) }kl=1 as

πt ({Vρ(l) }kl=1 = {vl }kl=1 | Tt ) ∝ P ({Vρ(l) }kl=1 = {vl }kl=1 )kYMvl (Aρ(l) | Jρ(l) )

l=1

=

k
Y

ξvl−1 ,vl (Aρ(l) )Mvl (Aρ(l) | Jρ(l) ),l=1

where v0 = 1. We show that for every k = 1, . . . , K

πt (Vρ(k) = vk | Tt ) ∝

I
X
v1 =1

I
X

( k
Y

vk−1 =1

l=1

···

)
ξvl−1 ,vl (Aρ(l) )Mvl (Aρ(l) | Jρ(l) )∝ ϕvk (As )(13)holds by induction. First, if k = 1, which is equivalent to s = 1, ρ(1) = 1, andAs = Ω, the posterior of V (Ω) is written as

πt (V (Ω) = v1 | Tt ) ∝ ξ1,v1 (Ω)Mv1 (A1 | J1 ) ∝ ϕv1 (Ω).

Second, assume that (13) holds for k = k̄. Then, if k = k̄ + 1, we haveπt (Vk̄ = vk̄ , Vk̄+1 = vk̄+1 | Tt )
∝

∝

I
X

···

I
X

k̄+1
πt ({Vρ(l) }k̄+1l=1 = {vl }l=1 | Tt )v1 =1

vk̄−1 =1

I
X

I
X

( k̄
Y

vk̄−1 =1

l=1

v1 =1

···

)
ξvl−1 ,vl (Aρ(l) )Mvl (Aρ(l) | Jρ(l) ) ξvk̄ ,vk̄+1 (As )Mvk̄+1 (As | Js )∝ ϕvk̄ (Aρ(k̄) )ξvk̄ ,vk̄+1 (As )Mvk̄+1 (As | Js ),49

from which we obtain

πt (Vk̄+1 = vk̄+1 | Tt ) =IX

πt (Vk̄ = vk̄ , Vk̄+1 = vk̄+1 | Tt )vk̄ =1∝I
X

ϕvk̄ (Aρ(k̄) )ξvk̄ ,vk̄+1 (As )Mvk̄+1 (As | Js )vk̄ =1∝ ϕvk̄+1 (As ).

Proof of Proposition 3.1
Let the finite tree Tt consist of a sequence of decisions J1:t = {Js }ts=1 , which sequentially divides nodes A1:t = {As }ts=1 . To derive the proposition for the marginal
posterior of Tt , we first consider the joint posterior of Tt and a sequence of the statevariables V1:t which are defined for the nodes A1:t . From the structure of the model,
the joint posterior is written as

πt (Tt , V1:t ) = P (J1:t , V1:t | x)1
= P (J1:t )P (V1:t )ZttY

MVs (As | Js )s=1,As ∈N (Tt )tYµ(x(As ) | As ),s=1,As ∈L(Tt )(14)
where Zt is the normalizing constant, and As,l and As,r are the children nodes of As .

For Tt−1 and Tt , since At is divided into At,l and At,r , we have

N (Tt ) = N (Tt−1 ) ∪ {At },L(Tt ) = L(Tt−1 ) \ {At } ∪ {At,l , At,r }.

50

With the expression of the joint posterior in (14), we obtain

πt (Tt , V1:t ) =µ(x(As,l ) | At,l )µ(x(As,r ) | As,r )Zt
πt−1 (Tt−1 , V1:t−1 )P (Jt )P (Vt | V1:t−1 )MVt (At | Jt ).
Zt−1µ(x(As ) | As )(15)

Let Apt denote the parent node of At and Vtp = V (Apt ). Then, since the state variables
follow the hidden Markov process, P (Vt | V1:t−1 ) = ξVtp ,Vt (At ). Integrating out V1:t−1 \Vtp in (15) gives
πt (Tt , Vtp , Vt ) =µ(x(As,l ) | At,l )µ(x(As,r ) | As,r )Zt.
πt−1 (Tt−1 , Vtp )P (Jt )ξVtp ,Vt (At )MVt (At | Jt )Zt−1µ(x(As ) | As )Because At is a leaf node of Tt−1 , by Lemma C.1, we haveπt−1 (Tt−1 , Vtp = j) = πt−1 (Tt−1 )πt−1 (Vtp = j | Tt−1 )= πt−1 (Tt−1 )ϕj (Apt ).

Hence, we obtain the expression of the marginal distribution of Tt as

πt (Tt ) =

I X
I
X

πt (Tt , Vtp = j, Vt = i)

i=1 j=1

( I
)
I
X
X
µ(x(At,l ) | At,l )µ(x(At,r ) | At,r )Zt
πt−1 (Tt−1 )P (Jt ),=ϕj (Apt )ξj,i (At ) Mi (At | Jt )
Zt−1
µ(x(A
)
|
A
)
t
t
i=1
j=1
which completes the proof.

51

Proof of Corollary 3.1
For πt (Dt | Tt−1 ), by Proposition 3.1, we obtainπt (Dt = j | Tt−1 ) =
∝

NX
L −1
l=1
NX
L −1

πt ((j, l/NL ) | Tt−1 )P (Dt = j, Lt = l/NL )h((j, l/NL ) | At )l=1

∝ λj (At )

NX
L −1

βl (At )h((j, l/NL ) | At ).

l=1

On the other hand, the conditional probability of Lt is obtained as follows:

πt (Lt = l/NL | Dt = j, Tt−1 ) ∝ P (Dt = j, Lt = l/NL )h((j, l/NL ) | At )∝ β(At )h(j, l/NL | Tt−1 ).

The expression of wt (Tt−1 ) immediately follows Proposition 3.1.
Proof of Proposition 3.2
In this discussion, we suppress x and T m in the expectation for simplicity. First,
when A = Ω, by the definition eΩ (i) = γ̃1,i (Ω). Next, if A is not the root node, we
can decompose eA (i0 ) as0eA (i ) =I
X

E[Q(A)I[V (A) = i0 ]I[V (Ap ) = i]].

i=1

52

For the summand, because θ(Ap ) and V (A) are conditionally independent given
V (Ap ), we obtain

E[Q(A)I[V (A) = i0 ]I[V (Ap ) = i]]= E[E[ϑ(Ap )I[V (A) = i0 ] | V (Ap )]I[V (Ap ) = i]Q(Ap )]= E[E[ϑ(Ap )I[V (A) = i0 ] | V (Ap ) = i]I[V (Ap ) = i]Q(Ap )]= E[E[ϑ(Ap ) | V (Ap ) = i]P (V (A) = i0 | V (Ap ) = i)I[V (Ap ) = i]Q(Ap )]= ξ̃i,i0 E[ϑ(Ap ) | V (Ap ) = i]eAp (i).

Therefore, we obtain0eA (i ) =IX

ξ̃i,i0 (A)E[ϑ(Ap ) | V (Ap ) = i]eAp (i).

i=1

C.2

Asymptotic properties of the tree posteriors

Lemma C.2. Let {Xni }n=1,2,... (i = 1, . . . , L) be sequences of random variables thatsatisfy the following conditions:
1. Xni > 0 for every i and n.
2. As n → ∞, the following convergence occurs
log Xni p→− ci ,nwhere c1 ≥ c2 ≥ · · · ≥ cL .
Then we have
log

PL

i
i=1 Xnn

53

p→− c1 .

(Proof) We only discuss the case where there exists l such that cl > cl+1 . Let  > 0.
Then, by the first condition,
log

P

!


i
X
log Xn1ni=1< c1 −  ≤ P< c1 −  → 0nn

PL

as n → ∞. On the other hand,

P

log

!
i
X
n
i=1
< c1 + 
n

PL

!
 \
l
L
\
log L maxi=1,...,L Xni< c1 +  ∩{Xni > Xnj }ni=1 j=l+1
!


l
L
\
\
log L maxi=1,...,l Xni{Xni > Xnj } → 1.
< c1 +  ∩ni=1 j=l+1


≥P

≥P

Lemma C.3. For T ∈ T K and A ∈ N (T ), if i ∈ τ (A | T ), thenG

log Mi (A | jA ) p X
ζg Pg (A) [Pg (Al | A) log Pg (Al | A) + Pg (Ar | A) log Pg (Ar | A)] ,
→
−
n
g=1
where jA is the splitting rule that divides A into Al and Ar .
(Proof) By the result of Schwarz (1978), since the parameter θ(A) follow the beta
distribution, which belongs to a continuous exponential family, the log of the marginal
likelihood is written asrilog Mi (A | T ) = ˆlA (i, T ) − log n(A) + Op (1),2
"G#YˆlA (i, T ) = log
θ̂g (A)ng (Al ) (1 − θ̂g (A))ng (Ar )
g=1

=

G h
X

i
ng (Al ) log θ̂g (A) + ng (Ar ) log(1 − θ̂g (A)) ,

g=1

54

(16)

where the definition of θ̂g (A) (the MLE) and ri (the number of parameters) depend
on which type of priors in Assumption 1 is introduced by the state i:θ̂g (A) =ng (Al )ng (A)
n(Al )

n(A)






 Pg (Al )Pg (A)(Prior A),(Prior B),(Prior C),G (Prior A),
ri = 1 (Prior B),






0 (Prior C).

(17)

Notice that, for Prior C, the constant c(A) and the true measures Pg must satisfyc(A) = Pg (Al | A)

for every g = 1, . . . , G because if this does not hold, the state i is not included in theset of feasible states τ (A | T ).
p

Since i ∈ τ (A | T ), the law of large numbers gives θ̂g (A) →− Pg (Al | A). Hence, we
obtain the limit of ˆlA (i, T )/n asGˆlA (i, T ) X
ng (Ω) ng (A) ng (Al )ng (Ar )=log θ̂g (A) +log(1 − θ̂g (A))nn(Ω) ng (Ω) ng (A)ng (A)
g=1
p

→
−

G
X

ζg Pg (A) [Pg (Al | A) log Pg (Al | A) + Pg (Ar | A) log Pg (Ar | A)] .

g=1

Lemma C.4. For T ∈ T K , A ∈ N (T ), i ∈ τ (A | T ) and j ∈ {1, . . . , I}, we havelog Mi (A | T ) − log Mj (A | T ) p→− ci,j ,nwhere ci,j = 0 if j ∈ τ (A | T ) and ci,j > 0 if j ∈ {1, . . . , I} \ τ (A | T ).

55

(Proof) If j ∈ τ (A | T ), obtaining the result
log Mi (A | T ) − log Mj (A | T ) p→− 0n
is straightforward from the proof of Proposition C.3. Hence we consider the case ofj ∈ {1, . . . , I}\τ (A | T ). Under the state j, for every g, the estimator θ̂g (A) is definedp

as in (17), and there exists Cg ∈ (0, 1) such that θ̂g (A) →− Cg . By the definition of
τ (A | T ), there exists g ∗ such that Cg∗ 6= Pg∗ (Al | A). As in the proof of Proposition
C.3, for the difference of the marginal likelihoods, we obtainG

log Mi (A | T ) − log Mj (A | T ) p X→−
ζg Pg (A)Λg ,ng=1
Λg = Pg (Al | A) log

Pg (Al | A)
Pg (Ar | A)
+ Pg (Ar | A) log.
Cg
1 − Cg

Because Λg is the KL divergence of the two discrete distributions, Λg ≥ 0 for all gand Λg∗ > 0. By Assumption 1, this result implies thatGX

ζg Pg (A)Λg > 0.

g=1

Proof of Theorem 4.1 (non-informative prior) and Theorem 4.2
This section provides the proofs of Theorem 4.1 under the assumption that the priorof the location variables is independent of the sample size, that is, η = 0 in Eq (2).
The general case is discussed in the next section.
In this proof, we modify the notation for the marginal likelihood defined in Eq. (4)and use Mi (A | T ) to represent the likelihood on A of the tree T under the ith state56

to reflect its dependency on the tree structure.
Let T ∈ T K and V denote a set of a combination of the states for all of the nonleaf nodes of T . Notice that an element of V does not need to satisfy P (V = v) > 0,
where V is the totality of the state variables. In the following proof, for v ∈ V, v(A)denotes a state on a node A. Let l(v, T ) denote the log of the joint likelihood functionl(v, T ) = log P (xn | T, v) =X

lA (v(A), T ) +A∈N (T )X

log µ(xn (A) | A),(18)

A∈L(T )

where lA (v(A), T ) = log MA (v(A) | T ). By Schwarz (1978), this likelihood lA has thefollowing expressionrv(A)log n(A) + Op (1),lA (v(A), T ) = ˆlA (v(A), T ) −2
where ˆlA and ri are defined in (16). Let v̄ ∈ V be a collection of states such that, forall A ∈ N (T ), θg (A) is fixed to µ(Al )/µ(A). For v̄, we havelA (v̄, T ) =X n(Al ) log

A∈N (T )µ(Al )µ(A)+ n(Ar ) logµ(Ar )
µ(A)


+

X

log µ(xn (A) | A)

A∈L(T )= log µ(x) = 0.

Hence l(v, T ) is rewritten asl(v, T ) = l(v, T ) − l(v̄, T ) =oX nˆlA (v(A), T ) − ˆlA (v̄(A), T ) − C(v) log n + Op (1).
2A∈N (T )For the part inside of the braces, when v is replaced with vT ∈ VT , where VT = {v :57

v(A) ∈ τ (A | T ) for all A ∈ N (T )}, the definition of ˆlA gives"#Gˆl(vT (A), T ) − ˆl(v̄(A), T ) X
ng (A) ng (Al )θ̂g (A)ng (Ar )1 − θ̂g (A)
=
log
+
log
n
n
n
(A)
µ(A
|
A)
n
(A)
µ(Ar | A)
g
l
g
g=1


G
X
Pg (Al | A)
Pg (Ar | A)p→−
ζg Pg (A) Pg (Al | A) log+ Pg (Ar | A) log.
µ(Al | A)µ(Ar | A)
g=1
For all A ∈ L(T ), there exists an unique sequence of nodesΩ = BA,0 ⊃ BA,1 ⊃ · · · ⊃ BA,K = A,

(19)

where BA,k ∈ T (k = 0, . . . , K) is a node in the kth level. With this sequence, we
obtain the limit of the scaled log-likelihood asl(vT , T )
n


G
X X
Pg (Al | A)
Pg (Ar | A)p
ζg Pg (A) Pg (Al | A) log→−+ Pg (Ar | A) logµ(Al | A)µ(Ar | A)
g=1
A∈N (T )

=

=

=

G
X

ζg


X 
Pg (Ar | A)
Pg (Al | A)
+ Pg (Ar ) log
Pg (Al ) logµ(Al | A)µ(Ar | A)

g=1

A∈N (T )

G
X

X

ζg

g=1

A∈L(T )

G
X

X

g=1

ζg

A∈L(T )Pg (BA,1 | BA,0 )Pg (BA,K | BA,K−1 )Pg (A) log+ · · · + logµ(BA,1 | BA,0 )µ(BA,K | BA,K−1 )G
Pg (A) X
ζg KL(Pg |T ||µ).
Pg (A) log
=
µ(A)
g=1

Because Pg |T admits the density function

pg |T (x) =X

1A (x)A∈L(T )58

Pg (A)
, x∈Ω
µ(A)



(20)

the KL divergence in (20) is rewritten as

KL(Pg |T ||µ) =X

Pg (A) log

A∈L(T )

Z
=

pg

X

Pg (A)µ(A)1A logA∈L(T )Pg (A)dµµ(A)Zpg log pg |T dµ=Z=pg (x)
pg (x) log
dµ(x) −µ(x)Z
pg log

pg
dµ
pg |T

= KL(Pg ||µ) − KL(Pg ||Pg |T ).

Because KL(Pg ||µ) is independent of T , we obtain another expression of TMK in (10)asGXKζg KL(Pg |T ||µ)TM = arg max
T ∈T K
g=1

By Lemma C.4 and (18), for v ∈ V \ VT , we can show that
X lA (v, T ) − lA (vT , T )l(v, T ) − l(vT , T )= p-lim
> 0,
n
n
n→∞
n→∞

p-lim

(21)

A∈N (T )and p-limn→∞ l(v, T )/n exists. Hence, for TM ∈ TMK , v0 ∈ VTM , T ∈ T K \ TMK andv ∈ V, we haveGGXl(v0 , TM ) − l(v, T ) X
p-lim≥
ζg KL(Pg |TM ||µ) −ζg KL(Pg |T ||µ) > 0.
n
n→∞
g=1
g=1

Hence, for such TM and v0 , we obtainP
X exp(l(v, T ))P (T )P (v)P (T ) v∈V exp(l(v, T ))P (v)P (xn | T )
p
P
=
≤
→
− 0.
0
0
P (xn | TM )P (TM ) v∈V exp(l(v, TM ))P (v) v∈V exp(l(v , TM ))P (TM )P (v )59pThis result implies p(T ∈ TMK | xn ) →− 1, which completes the proof of Theorem 4.1.
To prove Theorem 4.2, we fix T ∈ T K and define a set ST as

ST =

0



v ∈ VT | v ∈ arg min C(v ) .
v0 ∈VTp

Then we want to show P (V ∈ ST | T, xn ) →− 1. The result (21) impliespp(V ∈ VT | T, xn ) →− 1,

so we only need to compare the elements of VT . Let v ∈ VT \ ST and v0 ∈ ST . For
the difference of the log likelihoods, we havel(v0 , T ) − l(v, T ) =i0X hˆlA (v0 (A), T ) − ˆlA (v(A), T ) + C(v) − C(v ) log n + Op (1),2A∈N (T )where ˆlA and C is defined in (17) and (11), respectively. If v(A) and v0 (A) introduce
the same type of the prior (e.g., Prior A and Prior A), because the corresponding
estimators θ̂g (A) have the same form,ˆlA (v0 (A), T ) − ˆlA (v(A), T ) = 0.

On the other hand, if v(A) and v0 (A) introduce different types of the prior (e.g., Prior
A and Prior B), because they are the maximized log-likelihood under the two nested
hypotheses,−2[ˆlA (v0 (A), T ) − ˆlA (v(A), T )]60

weakly converges to the χ2 distribution (Wilks, 1938). Hence, we obtainl(v0 , T ) − l(v, T ) p C(v) − C(v0 )
→
−
> 0,
log n
2
p

which implies P (V ∈ ST | T, xn ) →− 1.
Proof of Theorem 4.1 (informative prior)This section provides the proofs of Theorem 4.1 for the case of the informative prior(η > 0). In the proof, we often use the results provided in the previous section on the
non-informative prior and use the same notations.
We first discuss the asymptotic behavior of the prior of trees, which is in this casedependent on the sample size n. We only need to focus on the prior of the location
variables L since this is only the component that depends on the sample size as given
in Eq (2). By the fact thatGn(A) p X
ζg Pg (A)
→
−
n
g=1

(22)

and Lemma C.2, we obtain the limitGX
log P (L(A) = µ(Al )/µ(A)) pµ(Al )→− −ηζg Pg (A)f
− 0.5 .
nµ(A)g=1
Since given T ∈ T K , the location variables are all independent and the prior of thedimension variables is independent of the sample size n, the log-tree prior l(T ) =log P (T ) has a limit as follows:GX
l(T ) pζg Bg (T ),
→
− −η
n
g=1

61

which is the penalty term introduced in Theorem 4.1.
We next define ψ(T ) as

ψ(T ) =GX

ζg {KL(Pg |T ||µ) − ηBg (T )} .

g=1

By the the proof of Theorem 4.1 for the non-informative case, this is a limit underv ∈ VT (this VT , and V, a collection of possible combination of the state variables on
the tree, are defined in the proof for the non-informative case), that is,l(T ) + l(v, T ) p→− ψ(T )n
as n → ∞. On the other hand, for v ∈ V \ VT , the limit exists andp-liml(T ) + l(v, T )< ψ(T ).
nAs discussed in the proof of Theorem 4.1, we have the equivalencearg max
T ∈T KGX

ζg KL(Pg |T ||µ) = arg min
T ∈T K

g=1

GX

ζg KL(Pg ||Pg |T ),

g=1

which implies thatKTM
= arg max ψ(T ).
T ∈T K

Hence, for TM ∈ TMK , v0 ∈ VTM , T ∈ T K \ TMK , and v ∈ V, we have(l(TM ) + l(v0 , TM )) − (l(T ) + l(v, T ))≥ ψ(TM ) − ψ(T ) > 0.
n
n→∞

p-lim

62

Hence, for such T, TM and v0 , it follows thatP
P (T ) v∈V P (v) exp(l(v, T ))P (T )P (xn | T )P (T, xn )P==
P (TM , xn )P (TM )P (xn | TM )P (TM ) v∈V P (v) exp(l(v, TM ))X exp(l(T ) + l(v, T ))P (v)
p
≤
→
− 0.
0
0
exp(l(TM ) + l(v , TM ))P (v )v∈VpTherefore, P (T ∈ TMK | xn ) →− 1.

DConsistency for the MRS model

To describe the consistency, for a possible node A, we define a variable Z(A) as
follows:

Z(A) =


1

if V (A) = 1,0

if V (A) ∈ {2, 3}.

Hence, θ1 (A) = θ2 (A) if Z(A) = 0 and θ1 (A) 6= θ2 (A) with probability one if Z(A) =1. Then we can obtain the following consistency result.
Corollary D.1. Let Z = {Z(A)}A∈N (T ) and z = {z(A)}A∈N (T ) be a collection ofZ(A) on T ∈ T K and one of its realizations, respectively. If P (Z = z) > 0 for any
possible z, then p
P Z(A) = 1{P1 (Al |A)6=P2 (Al |A)} for all A ∈ N (T ) | T, xn →− 1,

where 1 is the indicator function, and p− 1,
P Z(A) = 1{P1 (Al |A)6=P2 (Al |A)} for all A ∈ N (T ) | xn →63

where T is random.
(Proof) In this case, VT in Theorem 4.2 is written as

VT = {v | v(A) = 1 if P1 (Al | A) 6= P2 (Al | A)} .

We additionally define ṼT as

ṼT = {v | v(A) = 2 if P1 (Al | A) = P2 (Al | A)} .

Then, under the condition that v ∈ VT , the complexity C(v) is minimized if and only
if v ∈ ṼT . Hence, by Theorem 4.2 we obtain

P Z(A) = 1{P1 (Al |A)6=P2 (Al |A)} for all A ∈ N (T ) | T, xnp= P V ∈ VT ∩ ṼT | T, xn →− 1.

We can show the second result by using Theorem 4.1 as follows:

P V ∈ VT ∩ ṼT | xn =

X





P V ∈ VT ∩ ṼT | T, xn P (T | xn )T ∈T K≥ arg min P V ∈ VT ∩ ṼT | TM , xn P (T ∈ TMK | xn )K
TM ∈TM

p

→
− 1.

64

E

Numerical evaluation of the density estimators
under different complexity of data generatingprocesses

In this section, we observe the behavior of our proposed HMPT algorithm in density
estimation in a case where observed points concentrate around certain points or theboundary of the sample space, which occurs in real data analysis in many cases. We
consider the following scenario: the variables Xj (j = 1, . . . , d) independently follow
the beta distribution
Beta(2a , 2a ),where the parameter a takes values from R. When a = 0, the data is uniformly
distributed in the sample space (0, 1]d , and when a is positive, the distribution concentrates on the central point (0.5, 0.5, ..., 0.5). When a is negative, the distribution
concentrates on the vertices of the hyper-cube (0, 1]d . For the cases of a < 0, many
values generated from the beta distribution tend to be indistinguishable from 0 and 1,
and we found it often made the behavior of our HMPT algorithm, which is designed
for continuous distributions without mass points, unstable. As such, for the negativea’s, we use an approximative distribution in which the conditional distributions on(0, 10−5 ] and (0.9 − 10−5 , 1.0] are replaced with the uniform.
In the estimation, we used the APT model as described in Section 5.1. The number of particles is 1,000, and the maximum resolution is set to 10 or 15. The data is
generated under a = (−3, −2, −1, 0, 1, 2, 3), n = 102 , 103 , 104 , and d = 1, 2, 3, 4, andfor each combination, 50 data sets under different random seeds are generated. The
performance is assessed based on L1 distance between the estimated density (posterior mean) and the true distribution. Since the L1 distance is difficult to compute65

analytically, we use the Monte Carlo approximation with 10000 values generated fromthe true distributions.
Figure 9 provides the average L1 distance under the different settings. Whenthe distribution is uniform (a = 0), we can see that the L1 distance is very small,and this can be understood as an advantage of using the APT model to learn theshrinkage level adaptively. On the other hand, especially when a < 0 and d is large,
the distance is larger for a 6= 0 though the performance is slightly improved if the
maximum resolution K is larger. One possible reason to explain this result is that
the APT model, with which we learn the smoothness of unknown densities, is notdesigned to capture such extremely spiky distributions. Thus the performance would
be improved by modifying the prior distributions but this discussion is out of thisresearch’s scope because our main contribution is proposing the algorithm that works
for many types of PT-based models.
Figure 10 compares the number of nodes included in the MAP trees. From thisresult, we can see that the number is larger for the larger maximum resolution K = 15
as naturally expected but tends to be smaller when the parameter a is positive. The
latter phenomenon is explained by the rule that the SMC algorithm stops splitting
nodes when the number of included observations is below the threshold such as 5.
When a is positive, since the distribution concentrates around the central point,
nodes with only a few observations tend to be generated in early levels and thus
are no longer split, resulting in a smaller number of nodes on the tree. Hence thisnumerical result clarifies our proposed algorithm’s tendency of “ignoring” nodes with
only a few observations “zooming up” nodes including many observations.

66

Figure 9: The average L1 distance from the true distribution with error bars thatindicate the standard deviation obtained under different sample size, dimensionality,and values of the parameter a.

67

Figure 10: The average number of nodes in the MAP trees with error bars thatindicate the standard deviation obtained under different sample size, dimensionality,and values of the parameter a.

68

F

Details of the experiments

F.1

Density estimation

F.1.1

Hyper-parameter settings

For the HMPT model and the original APT model, we use the common settings
following Ma (2017). The transition matrix for the latent states on each A thatcharacterize different smoothness levels of the density is given byeβ(i−i0 ) if i ≤ i0 ,
ξi,i0 (A) ∝0
if i > i0 ,where β = 0.1 and the number of states I = 5. For the i(< I)th state, given V (A) = i,the precision ν(A) follows the priorlog10 ν(A) ∼ Unif(a(i), a(i + 1)],
where a(i) = L + (i − 1)(U − L)/(I − 1) with L = −1 and U = 4. For the Ith state,
FI = 1∞ , so ν(A) is fixed to ∞. This is equivalent to stopping the partition andputting the conditional distribution µ(· | A) on A. In the computation, this uniform
distribution is approximated by 5 evenly spaced grid points.

F.1.2

Simulation scenarios in the two-dimensional cases

The data sets are simulated from the following distributions.
1. “Blocks”:

1
1
1[0.1,0.45]×[0.35,0.9] (x1 , x2 ) + 1[0.2,0.8]×[0.45,0.5] (x1 , x2 )3
3
1
+ 1[0.7,0.9]×[0.05,0.6] (x1 , x2 )369

2. “Clusters”:
3
1
Beta(x1 | 1, 1) × Beta(x2 | 1, 1) + Beta(x1 | 15, 45) × Beta(x2 | 15, 45)
10
10
3
3
+ Beta(x1 | 45, 15) × Beta(x2 | 22.5, 37.5) + Beta(x1 | 37.5, 22.5) × Beta(x2 | 45, 15)
10
10
3. “Smooth”:

F.1.3

Beta(x1 | 10, 20)Beta(x2 | 10, 20)

Implementation details for the higher dimensional cases

For the algorithm of the original APT (implemented by the apt function in the PTT
package), in the first case with d = 6, the maximum resolution is fixed to 9 because
setting higher values lead to insufficient memory. Also, because the apt function does
not scale if the dimension is beyond d ≈ 10, in the second case with large d, we used
the proposed SMC algorithm to carry out inference for the original APT model as
well, which corresponds to setting NL = 2.
In the classical PT method, the Dirichlet prior is set to Dir(0.1k 2 , ..., 0.1k 2 ) (k:the depth of the node), and the maximum depth is set to 15. For the Dirichlet process
Gaussian mixture model, we used the PYdensity function in the R package BNPmix(Corradin et al., 2021), in which the strength and discounting parameters are set to
10 and 0, respectively, the covariance matrices are all diagonal (model = "DLS"), andthe size of the burn-in period and the sampling is both set to 1000. For the Gaussian
kernel density estimation, it is implemented by the kde function in the R package ks.
We set the bandwidth matrix by using Hpi.diag function which selects the optimal
diagonal matrix.

70

F.2

Two-group comparison

F.2.1

Hyper-parameter settings

For the transition matrix ξ(A), we use the form proposed in Soriano and Ma (2017)
for incorporating multiple testing control(1 − ρ)(1 − γ)ρ  (1 − ρ)γ−kξ(A) = (1 − ρ)(1 − γ2−k ) ρ  (1 − ρ)γ2
,


0
0
1
where γ ∈ (0, 1), and ρ ∈ (0, 1), and k is the depth of A, and we set (γ, ρ) = (0.3, 0.3)
following recommendations in that paper.

F.2.2

Simulation scenarios

1. “Local location shift”: For j = 1, . . . , 25,
3X
11
N (µl , Σ),(X1,2(j−1)+1 , X1,2j ) ∼ N (µ1 , Σ) +
3
3
l=2
3
X
1
1
(X2,2(j−1)+1 , X2,2j ) ∼ N (µ1 + δj , Σ) +N (µl , Σ),33l=2

where δj = −0.5 for j = 1, . . . , 5 and 0 for j = 6, . . . , 25.
2. “Local dispersion difference”: For j = 1, . . . , 25,
3
X
1
1(X1,2(j−1)+1 , X1,2j ) ∼ N (µ1 , Σ) +N (µl , Σ),
3
3
l=2
3
X
1
1
(X2,2(j−1)+1 , X2,2j ) ∼ N (µ1 , Σ + ∆j ) +N (µl , Σ),33l=2

where ∆j = −0.4 for j = 1, . . . , 5 and 0 for j = 6, . . . , 25.

71

3. “Correlation”: For j = 1, . . . , 25,


 



 0   1 0  (X1,2(j−1)+1 , X1,2j ) ∼ N   , 00 1
  

 0   1 δj (X2,2(j−1)+1 , X2,2j ) ∼ N   , 
 ,
0
δj 1
where δj = 0.75 for j = 1, . . . , 5 and δj = 0 for j = 6, . . . , 25 .
In the “local location shift” and “local dispersion difference”, the parameters are 0.5 0 
µ1 = (−2.5, 1.0), µ2 = (1.0, −2.0), µ3 = (2.0, 2.5), Σ = 
.
0 0.7

G

Additional figures

Figure 11: The average KL divergences between the estimated density and the true
density.

72

Figure 12: The comparison of the predictive performance. Each point corresponds
to the average of the predictive score in Eq. (12) based on 50 data sets. In the rightplot, the predictive scores of the DPM model for the over 30 dimensional cases are
below the displayed range.

Figure 13: The MAP tree for the mass cytometry data set. The size and the color
indicate the estimated eff(A). Only the nodes with the sample size larger than 50 are
drawn. Since there are a huge number of nodes, the information on the dimension isomitted in this figure.

73

Figure 14: The solid lines delineate the nodes with the highest values of eff(A) andtheir two children. The red triangle points and the blue circle points are the observations from the two samples in the node. The observations outside the node are ingray. In this figure, the nodes with ng (A) ≥ 10 (g = 1, 2) are chosen.

74