A Bayesian hierarchical model for related densities using Polya trees
Abstract
Bayesian hierarchical models are used to share information between related samples and obtain more accurate estimates of samplelevel parameters, common structure, and variation between samples. When the parameter of interest is the distribution or density of a continuous variable, a hierarchical model for distributions is required. A number of such models have been described in the literature using extensions of the Dirichlet process and related processes, typically as a distribution on the parameters of a mixing kernel. We propose a new hierarchical model based on the Polya tree, which allows direct modeling of densities and enjoys some computational advantages over the Dirichlet process. The Polya tree also allows more flexible modeling of the variation between samples, providing more informed shrinkage and permitting posterior inference on the dispersion function, which quantifies the variation among sample densities. We also show how the model can be extended to cluster samples in situations where the observed samples are believed to have been drawn from several latent populations.
1 Introduction
Many statistical applications deal with learning and comparing the distributions of two or more related samples. We may be interested in learning how samples are similar or testing whether they are distinguishably different from each other. Because distributions are complex infinitedimensional objects, classical approaches work with lowdimensional parameterizations of the distribution. Analysis of variance, for example, reduces distributions to a mean and variance, which are sufficient under the assumption of normality. A wide range of other parametric models within both the Bayesian and frequentist inferential frameworks use other parameterizations of the distribution to reduce the dimensionality of the problem. Classical nonparametric approaches use features of the samples such as medians[31], rankbased scores[32], or summaries of the empirical distribution functions, as in the KolmogorovSmirnov[14] and Cramérvon Mises tests[1].
A number of Bayesian nonparametric approaches embrace the infinitedimensional nature of the problem using extensions of Dirichlet processes. The Hierarchical Dirichlet process of Teh et al. [29] builds a hierarchical model using the Dirichlet process, allowing it to share information between samples. The Nested Dirichlet process of Rodriguez, Dunson, and Gelfand [24] takes a different approach, using a Dirichlet process as the base measure of a second Dirichlet process. This induces clustering in the samples, with samples in a cluster being modeled with a single density. Müller, Quintana, and Rosner [23] model each sample density as a mixture of two components: one Dirichlet process mixture of Gaussians representing common structure between samples, and a second Dirichlet process mixture of Gaussians representing the idiosyncratic structure of the given sample. A variety of other dependent Dirichlet Processes[19] have been described in the literature. Teh and Jordon [28] give an overview of hierarchical models based on the Dirichlet process.
While the Dirichlet process has been the basis of most of the work in this area, work has also been done on hierarchical extensions of other priors. For example, Teh [27] defines a Hierarchical PitmanYor process, taking advantage of the more flexible clustering structure of the PitmanYor process over the Dirichlet process. Camerlenghi et. al [3] consider hierarchical models based on Normalized Random Measures [2], which includes both the Dirichlet process and the PitmanYor process as special cases.
An alternative approach for modeling densities within the Bayesian nonparametric framework is to use a model derived from Polya trees [9]. Polya trees are a class of tailfree prior in which an infinite recursive binary partition is placed on the sample space, and probability mass assigned to the elements of the partition by a corresponding infinite sequence of Betadistributed random variables. A special case gives the familiar Dirichlet process [7], but the Polya tree family is considerably more flexible. Under appropriate specification of the prior parameters, the Polya tree assigns probability 1 to the set of absolutely continuous distributions [8]. This property allows it to be used to directly model probability densities without the encumbrance of a mixing kernel.
The Polya tree model allows tractable computation of the marginal likelihood, which has made it popular in Bayesian hypothesis testing of nonparametric density models. Examples include [13, 12]. These approaches are unsatisfactory when estimation and prediction rather than formal hypothesis testing are of primary interest. Statisticians have recognized the benefits of partial shrinkage as far back as Stein’s shrinkage estimator for the mean of a multivariate normal distribution [26]. We may expect that related samples will have similar but not identical distributions. In this situation, neither independence nor complete shrinkage are appropriate models. A partial shrinkage model allows borrowing of information across samples while preserving betweensample variation.
Within the Bayesian inference framework, partial shrinkage is naturally achieved with a hierarchical model. While a simple hierarchy of Polya trees is possible and is described herein, we build on the richer Adaptive Polya tree model of Ma [18]. This model allows us to learn the concentration parameters of the Polya tree, rather than fixing them in the prior. In contrast to Dirichlet processbased models, the richer structure of the Polya tree allows us to model both the means and the betweengroup variation of the densities in a fully nonparametric manner. This not only improves the estimation of the distributions, but also allows us to perform inference on how the variation between groups differs over the sample space.
In Section 2 we describe the model and contrast it with several existing models. We discuss Bayesian posterior inference and computation in Section 3. In Section 4 we give theoretical grounding for the appropriateness of the method. Section 5 describes two ways in which the model allows us to go beyond the capabilities of existing models. We show simulation results in Section 6 and demonstrate application to real data in Section 7. Finally, we conclude with some discussion in Section 8.
2 Model
2.1 Reviewing the Polya tree construction
We begin with a brief sketch of the Polya tree; the reader interested in the mathematical details should refer to [21, 16, 17]. Applications of the Polya tree are discussed in, for example, [10, 30]. The Polya tree consists of an infinite recursive partition of the sample space and a corresponding infinite sequence of Betadistributed random variables which assign mass to the various regions of the partition. Figure 1 illustrates the partitioning sequentially. In this illustration our sample space is the interval on the real line, and our prior mean, shown in the first pane, is the uniform distribution on that interval. At each level of the recursive partition we cut each region in half at the midpoint. Although arbitrary partitions may be used, the dyadic partition described here is convenient and is often used as a default partition in the absence of a compelling reason to use a different one. Other partitions may be more convenient in the presence, for example, of censored data [22]. The second pane shows the result after the first cut and mass allocation, in this case the majority of the mass having been allocated to the righthand side. In the next step we cut each of the two regions from the second pane in half again, and assign the probability mass to each to its children according to a Betadistributed random variable. This results in four regions, shown in the third pane, which are again cut and mass distributed in the fourth pane. The process continues indefinitely.
We denote the Polya tree prior as , where is the centering distribution and is an (infinite dimensional) concentration parameter describing how much is expected to vary from . The parameters of the sequence of Beta distributions from which the mass allocations are drawn are derived very simply from and . For an arbitrary region belonging to the recursive partition , the fraction of the mass allocated to the left child of is given by the random variable
where with the remainder of the mass allocated to the right child . thus determines the mean of the mass allocations (and hence the expectation of the resulting density), while controls the variation of the mass allocations, and hence the dispersion of around . Alternatively, controls the strength of the shrinkage of the posterior mean density from the empirical process towards . With an appropriate choice of , the Polya tree prior almost surely generates an absolutely continuous distribution [15].
2.2 The Hierarchical Polya Tree
It is conceptually straightforward to extend the Polya tree to a hierarchical model. Let be samples arising from related distributions on a complete, separable space . For ease of exposition we again use , though like the Polya tree the model can be applied to more general sample spaces. We model these samples as coming from exchangeable distributions , which are centered at a common underlying measure , itself unknown. Applying Polya tree priors to both and the gives us the hierarchical model
(1) 
Here is the overall prior mean, controls the variation among samples around the common structure , and controls the variation of the common structure from , which determines the smoothness of . This model, which we call the Hierarchical Polya tree, allows nonparametric estimation of the sample distributions and the common structure . The Hierarchical Polya tree model is illustrated in Figure 2. The first row shows the upper level of the hierarchy, which like the basic Polya tree illustrated in Figure 1 is centered at the uniform distribution on . The second row shows the lower level of the hierarchy, the individual group distributions conditioned on . They follow exactly the same Polya tree construction, but each cut is centered on the corresponding cut from , rather than on the uniformed distribution. captures the common structure of the groups, while the model the idiosyncratic structure of each group.
The hierarchy of Polya trees translates directly to the decomposed space as a hierarchical model for Beta random variables. For an arbitrary region , we have
(2)  
The representation of the hierarchy of Polya trees as a hierarchy of Beta random variables allows tractable posterior inference, as described in Section 3.
2.3 The Stochastically Increasing Shrinkage prior on dispersion
The Polya tree’s concentration parameter is traditionally set to increase with depth at a predetermined rate to ensure absolute continuity, with a constant multiplicative term to control the overall level of variation which may be treated as a tuning parameter (as in [13]) or have a prior placed on it. Hanson [11] discusses some of the necessary considerations when placing a prior on this parameter. However, Ma [18] shows that putting a fully nonparametric prior on the concentration parameter allows the Polya tree to learn the true distribution more accurately, particularly when the smoothness of the underlying density varies over the sample space. We can extend the Hierarchical Polya tree model by placing priors on both concentration parameters. In addition to more accurate inference on and the , putting a prior distribution on allows us to learn the variation between samples in a nonparametric way. That is, we can estimate a posterior dispersion function which summarizes the variability between sample densities at any given point in the sample space. Dispersion can be measured in a variety of ways; in Section 5.1 we show how to derive a function which gives the posterior mean variance of the densities, as well as a standardized version using the coefficient of variation to correct for the height of the density. Nonparametric inference on the variation between samples across the sample space is made possible by the flexibility of the Polya tree model. We contrast how several other models treat betweensample variation in Section 2.5.
A simple approach is to put independent priors on the variance of each Beta random variable. However, we expect spatial structure in the dispersion function—depending on how smooth the underlying densities are, the amount of Beta variance is generally smaller for deeper levels of the partition though the decay in the Beta variance can be heterogeneous over the sample space depending on the local smoothness of the densities. While the recursive partitioning allows independent priors to capture some spatial structure, we can do better by introducing dependency between regions in the partition. Ma [18] introduces Markov dependency on the concentration parameter, following the tree topology. The Stochastically Increasing Shrinkage (SIS) prior introduces a state variable supported on a finite set of integers , corresponding to decreasing prior variance and increasing shrinkage for the Beta random variables. For example, we may have with implying with the stochastically ordered and being a point mass at infinity. The number of states and the corresponding distributions can be chosen to balance the flexibility and computational complexity of the model. A simple way to enforce such a stochastic ordering is through partitioning the support of the concentration parameter into disjoint intervals. Given these latent states, the SIS prior adopts a transition probability matrix for , denoted . Ma [18] discusses several prior possibilities for this transition matrix, including an empirical Bayes estimate. We adopt the simplest recommendation there that provides generally adequate flexibility, namely
This uppertriangular transition probability matrix induces stochastically increasing shrinkage, and ensures that the model generates absolutely continuous densities (see Theorem 1). We denote the SIS prior
Figure 3 illustrates the SIS prior in action. As you move down to finer resolutions the shrinkage state tends to increase, eventually reaching complete shrinkage. Because the transitions are stochastic and the state of each node can be learned from the data, the amount of shrinkage or prior variance can vary over the sample space to capture large smooth features in one part of the sample space and smaller scale features in another part. This allows the resulting density to have heterogeneous smoothness across the sample space.
2.4 The Hierarchical Adaptive Polya Tree
Having described the SIS prior, we can adopt this prior for the concentration parameters and , and write the complete model as follows:
We call this model the Hierarchical Adaptive Polya Tree (HAPT). In contrast to existing models, this specification allows fully nonparametric inference on both the densities and the variation between densities. The inclusion of the SIS priors also allows the model to adapt the level of shrinkage or information borrowing in different parts of the sample space to more accurately capture the densities of each group, rather than using fixed uniform shrinkage. Figure 4 shows a graphical representation of the model.
2.5 Comparison to existing models
The most prominant existing nonparametric models for estimation of related distributions are based on the Dirichlet process, including the Hierarchical Dirichlet process (HDP) [29], Nested Dirichlet process (NDP) [24], and the Hierarchichal Dirichlet Process Mixture (HDPM) [23]. The Hierarchical Adaptive Polya tree enjoys several advantages over these methods:

Nonparametric estimation of betweensample variation. The HDP and NDP have a single concentration parameter that controls the dependence between samples. The HDPM has one parameter per sample that controls what proportion of the sample is explained by common structure and how much by idiosyncratic structure.
In contrast, HAPT places a fully nonparametric prior on the variation between samples, which allows it to learn spatially heterogeneous variation between samples. Indeed, the variation between samples at different locations of the sample space may be of primary interest in some scientific applications: learning where common structure is largely preserved between samples and where distributions vary widely may point the way to understanding important underlying phenomena.

Computation. HDP, NDP, and HDPM all rely on MCMC methods to draw from the posterior. The Hierarchical Adaptive Polya tree is not fully conjugate, but the necessary integration can be split into lowdimensional integrals and approximated extremely quickly using adaptive quadrature methods, without concerns about Markov chain convergence. See Section 3.1 for details.

Interpretability. HAPT provides an easily interpretable estimate of the common structure: The posterior estimate of is both the estimate of the mean density across samples, and the expected value of the posterior predictive for a new sample. In contrast, the HDP estimates a discrete instead of continuous distribution, and the NDP does not provide any estimate of common structure. The HDPM provides an estimate of common structure, but it is neither the mean of sample distributions nor a posterior predictive estimate. Interpretation of the common structure in the HDPM model is most straightforward if variation between samples involves contamination of an underlying distribution by an idiosyncratic process for each sample.
Other existing models based on Normalized Random Measures[3], including the Hierarchical PitmanYor model[27], can suffer from drawbacks similar to those of the models based on the Dirichlet process, including not allowing flexible modeling of the dispersion between samples, requiring MCMC samplers, and requiring mixture kernels to model continuous distributions.
3 Bayesian inference and computation
The HAPT model is partially conjugate: the conditional posterior for is a standard Polya tree. Though not conjugate, the remaining joint posterior for can be reliably approximated using adaptive quadrature methods. The computational strategies used are described in Section 3.1
To derive the posterior we use the representation of the Polya trees and in terms of Betadistributed random variables and for each node of the tree. With this notation, The second and third lines in (1) can be written in terms of the and as in Equation 2:
Including the concentration parameters, we can write the posterior for the parameters of a region in the following form conditional on the state parameters :
(3)  
where is the Beta function. The full posterior is the summation of Equation (3) over the possible states of and , with their respective priors factored in. The derivation of this posterior is given in the supplementary material.
3.1 Computation
Posterior computation of the HAPT model requires three distinct computational strategies. We split the model (see Figure 4) into three parts, each of which requires a different approach. Each part of the model is conditioned on all parameters which are further to the right in Figure 4. We describe how to integrate out each of the first two parts, which reduces the problem to evaluating the posterior probabilities of all possible combinations of the state variables .

: The individual sample densities , conditional on all other parameters, are a priori distributed according to a standard Polya tree. The corresponding conditional posterior is therefore also a Polya tree. This allows us to integrate out the when computing the posterior. If individual sample densities are of inferential interest their posteriors can easily be reconstructed after the main posterior computation is completed.

: The remaining continuous parts of the joint model, namely the common structure and the continuous concentration parameters and , conditioned on the discrete state parameters and , are not conjugate and must be integrated numerically. Because all parameter dependence across nodes in the Polya tree topology is restricted to the discrete state parameters, by conditioning on those parameters we are able to compute the posterior of the remaining parameters for each node of the tree independently.
This has two significant implications. First, rather than tackling a very high dimensional integral over the product space of the parameters for all nodes, we have a much more tractable collection of lowdimensional integrals: we need only integrate the threedimensional joint posterior of for each region in the Polya tree. Each of these integrals is tractable using adaptive quadrature techniques. Second, these integrals can be computed in parallel.
An additional observation allows us to further accelerate the adaptive quadrature. We note that the joint posterior distribution for conditional on and , with the other parameters integrated out, can be factored as
This allows us to factor the three dimensional integral:
This factorization effectively reduces the dimensionality of the integral: rather than evaluating the unnormalized posterior at points throughout a 3dimensional space, we need only evaluate it on the union of two 2dimensional spaces.

The posterior distribution of the state parameters and at first appears to be the most intimidating part of the model: It is a distribution over the product space of a large number of discrete parameters, resulting in an enormous number of level combinations. The naive computation of the joint posterior,
is straightforward but needs to be repeated for every possible combination of and for every node in the tree, which is computationally prohibitive. Here the Markov dependency structure comes to our rescue. The shrinkage states constitute a Hidden Markov Model following the tree topology[5], and we can factor the joint distribution and calculate the posterior probabilities using a forwardbackward algorithm in a manner analogous to inference strategies for linear Hidden Markov Models.
During the forwardbackward algorithm we can compute expectations of any function that can be expressed in the form
where is an arbitrary function in . This includes the marginal likelihood, the expected value of the estimated common density or any individual sample density at any given point, expectations of random variables or , and a wide variety of other functions, such as the variance function described in Section 5.1.
This computation is recursive, and we give a brief example of how it is carried out for the marginal likelihood. The previous two computational strategies give us the ability to calculate (up to quadrature approximation) the marginal likelihood, within a given node, of all remaining parameters conditional on and . Combining this with the prior transition parameters specified in we are able to calculate the posterior probabilities for and , and the overall marginal likelihood of the distribution on , by considering the recursivelycalculated marginal likelihoods of the child regions of under each possible state.
Obviously this recursion requires a stopping point. The simplest method is to truncate the tree at a predetermined depth. Hanson[11] offers some guidance on how to choose the depth of a truncated Polya tree based on sample size and other considerations. We recommend using as large a tree as is computationally feasible in order to minimize approximation errors due to truncation. If the data deviate very strongly from the prior distribution—as is common, for example, in highdimensional settings—a more sophisticated approach may be required, such as truncating a branch of the tree when it reaches a depth where the node contains only a few data points.
4 Theoretical results
Theorem 1.
(Absolute Continuity) Let be random measures distributed according to a HAPT model with base measure . If the SIS priors on and each include a complete shrinkage state that absorbs all possible chains in a finite number of steps with probability 1, then with probability 1.
Remark.
A sufficient condition for the complete shrinkage state to absorb all chains in a finite number of steps with probability 1 is that the transition probability from every other state to the complete shrinkage state is bounded away from zero, which is satisfied by our choice of .
Theorem 2.
(Prior support) Let be the probability density functions corresponding to arbitrary absolutely continuous distributions on . Let be corresponding densities from a HAPT model satisfying

and are equivalent measures, that is and ;

, i.e. there are at least two shrinkage states at each level.
Then are in the prior support.
These results follow immediately from [18]. Our final result gives posterior consistency in the case where the groups have equal sample sizes:
Theorem 3.
(Posterior consistency) Let be observed data consisting of independent samples, each of size , from absolutely continuous distributions . Let be a Hierarchical Adaptive Polya tree model on the densities with overall prior mean , and let be the corresponding posterior. If for all then we have posterior consistency under the weak topology. That is,
for any weak neighborhood of the product measure .
Proofs are given in the supplementary material.
5 Methodological applications of HAPT
We present two ways in which the HAPT model can be applied to infer structure that existing models have not been able to capture. The first application is the ability of HAPT to model the “dispersion function” (defined below) on the sample space; we show how to calculate the fitted dispersion function from the parameter. No existing model permits inference on the variation between sample densities in this manner. The second application is clustering samples based on their distributions, while allowing for withincluster variation. While the Nested Dirichlet process clusters distributions, it allows no variation (beyond sampling variation) within clusters. The importance of allowing for variation within clusters was first pointed out by MacEachern [20], who described a dependent Dirichlet process which would incorporate withincluster variation.
5.1 Inferring the betweensample dispersion function
The primary target of inference in problems with multiple samples is often the variation between samples. It is this inference, for example, which lends ANOVA its name, though the model is typically presented in terms of the overall and sample means.
In the HAPT model, variation between samples is captured by , an infinitedimensional parameter which estimates variation at different locations and scales. Rather than trying to provide guidance on how to interpret the multiscale structure in , we show how to recast it into an estimate of the variation between samples at any given point in the sample space, giving us a posterior variance function analogous to the posterior mean function.
Let be the posterior predictive density for a new sample, with corresponding Betadistributed random variables for each region in the recursive partition. Since is random, we wish to estimate the variance function which gives, for any point in the sample space, the variance of conditional on the density of the common structure, . This is analogous to the variance between treatments in traditional ANOVA. We have
where is the density corresponding to the prior mean distribution . We note that
We will need the facts that
and
which together imply
Let . It is clear that contains exactly the information about which is relevant to ; that is, the distribution of is identical to the distribution of . With the above, we can calculate
We note that the expectation in the second term can be factored into a product of expectations, none of which depend on or . The variance of this product is thus zero. We rewrite the remaining line using the identify :
Conditional on and , the terms in both products are mutually independent, both a priori and a posteriori. Thus, we can factor the expectations:
We can now plug in our expressions for the expectations, calculated above:
Since is itself unknown, we take the expectation with respect to it:
We call this expectation the variance function. Finally, we note that conditional on , are independent of for any two distinct regions and . This allows us to rearrange our expectation to obtain
The expectations with respect to and can be estimated during the same quadrature routines used to compute the posterior distributions of , described in Section 3.1. The expectation with respect to can then be calculated during the forwardbackward routine used to calculate the posterior distribution of , as described in the same section.
In addition to the variance we can consider other measures of dispersion. For example, because we naturally expect more absolute variation between samples in areas where the density of all the samples is higher, we can consider a standardized dispersion function measuring the coefficient of variation. The posterior mean coefficient of variation of at any given point is not analytically tractable; we can obtain an estimate by taking the square root of the variance function and dividing by the mean density function. We illustrate the application of these dispersion functions in Section 6.1.
5.2 Dirichlet Process Mixture of HAPT
In many applications we may not believe that the samples collected all share a single common structure. A more appropriate model may be that the samples are drawn from several latent populations, with samples being drawn from the same population having structure in common. In this case we may reconstruct the latent structure by clustering the samples. To learn the clustering of samples without fixing the number of clusters in advance, we add a Dirichlet process component to the model. We can write the model as follows:
where the base measure can be factored as
The Dirichlet process introduces clustering among the samples, so that some set of , belonging to the same cluster, share a corresponding and . Conditional on the clustering structure, the model reduces to a collection of independent HAPT models.
We call this model a Dirichlet Process Mixture of Hierarchical Adaptive Polya Trees, or DPMHAPT. It is comparable to the Nested Dirichlet process in the way it induces clustering among the samples, but is considerably more flexible. While the NDP requires that all samples in a cluster have identical distributions [20], DPMHAPT allows the distributions within cluster to vary according to the HAPT model. In addition, the advantages of the HAPT model discussed earlier, such as flexible modeling of variation in different parts of the sample space, still apply.
Posterior computation for HAPTDPM consists of a combination of standard Dirichlet Process methods and the HAPT posterior calculations described in Section S1: Posterior derivation. As noted above, conditional on the clustering structure, the model consists of a number of independent HAPT models. Although the HAPT model is not fully conjugate, our posterior computation strategy allows us to calculate the marginal likelihood with arbitrary precision. This allows the use of a Dirichlet process algorithm designed for conjugate mixture models. We use the Chinese Restaurant Process representation to sample the clustering structure, using marginal likelihoods calculated from the HAPT model.
This algorithm requires the computation of

marginal likelihoods under the HAPT model for clusters including a single sample, and

marginal posterior predictive likelihoods under the HAPT model of one sample conditional on one or more other samples making up a cluster.
The first item is straightforward. The second is easily achieved by fitting the HAPT model twice, and calculating
6 Simulation results
6.1 Estimation of the dispersion function
In order to obtain data with dispersion that varies across the sample space, we simulate from a mixture of four Beta distributions, such that the variance of the sample densities is much larger on the right half of the space than the left. Details of the simulation setting are described in the supplementary material. One hundred sample densities are plotted in Figure 5(a). The variation in the dispersion of the sample densities can be seen clearly.
Figure 5(b) shows the estimated centroid and sample densities from fitting the HAPT model. The posterior mean variance function of the sample densities is plotted in Figure 5(c); because the variance of the sample densities is naturally tied to the height of the density, the estimated coefficient of variation—obtained by taking the square root of the posterior mean variance and dividing by the posterior mean density value—may be more interpretable, and is plotted in Figure 5(d). While the variance plot is dominated by the spike in variance at the far right of the sample space, the standardized dispersion function on the coefficient of variation scale shows a clearer picture of the variation between samples growing, falling to near zero around 0.8, and then growing again.
6.2 DPMHAPT Simulations
We simulate a simple 1dimensional example to demonstrate the clustering behavior of the DPMHAPT model. The simulation contains 30 samples belonging to three true clusters, with 15, 10, and 5 samples respectively. Each sample is drawn from a mixture of a Uniform(0,1) distribution and a Beta distribution, with the parameters of the Beta varying by cluster: Beta(1,5) for the first cluster, Beta(3,3) for the second cluster, and Beta(5,1) for the third cluster. The weights of the two components are randomized in each sample. The weight of one component is drawn from a Beta(10,10) distribution, which creates weights varying approximately between 0.3 and 0.7, with the actual observed proportions in realized samples varying more widely due to the additional Binomial variation. Each sample contains points. Histograms of the samples from one simulation run are shown in Figure 6.
An MCMC sampler is run using the Chinese Restaurant Process to sample the clustering structure. We summarize the results by looking at how often each of the 30 samples is clustered together with each other sample. These results are plotted in Figure 7(a). We can see in the figure that the DPMHAPT model clearly identified the three clusters. In contrast, the clustering structure from the Nested Dirichlet process applied to the same data is shown in Figure 7(b). This example demonstrates the inability to identify cluster structure when the variation between samples within a cluster is substantially larger than sampling variation alone can account for.
We now consider an example in which the variation between samples differs across the sample space. Samples with heterogeneous variation are drawn from mixtures of betas as described in Section 6.1, with varying parameters for each cluster. We consider 30 samples belonging to three true clusters, as above. Figure 8 (a) shows one hundred draws from each of three clusters used in this simulation, to illustrate the variability between and within clusters and the heterogeneity across the sample space. We can see that the clusters have substantially more withincluster variation on the right half of the space than on the left half.
As in the previous example, we simulate three clusters with 15, 10 and 5 groups respectively. We draw 200 observations from each group, and apply DPMHAPT to cluster the resulting samples. We run the MCMC for 10,000 draws after burnin; estimated coclustering probabilities are plotted in Figure 8 (b). The three true clusters are clearly identified, with a few groups having ambiguous cluster membership.
7 Applications
7.1 DNaseseq profile clustering
DNase sequencing (DNAseseq) is a method used to identify regulatory regions of the genome [25]. DNA is treated with Deoxyribonuclease (DNase) I, an enzyme that cuts the DNA strand. The cut strands are then sequenced and the locations of the cuts are identified and tallied. The vulnerability of the DNA strand to DNase varies by location, resulting in a distribution of cut counts which is nonuniform. The density of this distribution is related to various biological factors of interest; for example, it tends to be high near potential binding sites for transcription factors, since these proteins require access to the DNA strand in much the same way as DNase I, but will be low if a transcription factor already bound at that site blocks access to the DNA strand.
We consider the problem of clustering DNaseseq profiles near potential transcription binding sites, identified by a specific genetic motif. Each sample consists of observed counts in a range of 100 base pairs on either side of one occurrence of the motif. A single motif, consisting of 10–20 base pairs, may appear tens of thousands of time in the genome, with each occurrence presenting one sample for analysis. Many samples, however, have very few cuts observed. For analysis we restrict ourselves to samples which meet a minimum sample size threshold.
Different locations where the transcription factor motif of interest appears may be expected to show different DNase behavior in the region around the motif for a variety of reasons. This makes clustering a more appropriate approach to the problem than treating all the samples as having a single common structure. Identifying clusters of locations which have similar DNaseseq profiles may reveal previously unrecognized factors. We also expect withincluster variation beyond simple sampling variation, which makes the Nested Dirichlet process unsuitable.
Here we present data from the ENCODE project [6] for locations surrounding a motif associated with the RE1silencing transcription factor (REST). REST suppresses neuronal genes in nonneuronal cells [4]. The data includes 48,549 locations where the REST motif appears in the genome. The motif consists of 21 base pairs, and the data includes an additional 100 base pairs on each side, for a total of 221 base pairs. In all, 922,704 cuts were tallied, an average of 19 per location. 468 locations have zero cuts observed. The number of cuts per location is heavily right skewed, with a median of 13 observations, first and third quartiles of 7 and 21 respectively, and a maximum of 2,099 cuts observed in a single sample.
For this analysis we restrict ourselves to locations which have at least 200 observations, a total of 265 samples. These samples include a total of 70,225 observations, an average of 330 observations per sample. The distribution is still quite skewed, with a minimum of 201 observations and a maximum of 2099. The median is 279 and the first and third quartiles are 232 and 366 observations.
We fit the DPMHAPT clustering model to this data, using 200 postburnin draws for inference. The model estimates between 8 and 13 clusters with high probability (see Figure 9(a)). A heatmap of the clustering structure is shown in Figure 9(b). The clusters are quite clearly differentiated and vary in size, with the largest clusters containing about 50 locations.
The estimated mean densities of the three largest clusters are plotted in Figure 10.
One of these clusters includes locations with cuts which are roughly symmetric around the transcription factor binding site, while the other two largest clusters include locations which have cuts heaped up on one side or another of the binding site. Additional clusters show other densities which do not conform to the general patterns of the three largest clusters. The plots show substantial variation around the cluster centroid, much more than can be explained by sampling variation alone. The estimated dispersion functions are quite noisy; this is not surprising given that they are betweensample dispersions estimated with only 40–50 degrees of freedom.
8 Discussion
The HAPT model offers a compelling alternative to existing nonparametric models that share information across multiple samples. The Polya tree’s ability to directly model the density of an absolutely continuous distribution frees us from the necessity of using mixture models to obtain densities, while the computational tractability of the posterior avoids the need to run MCMC chains for posterior inference. The addition of the SIS prior allows us not only to more accurately model the densities of interest, but also to estimate a fully nonparametric dispersion function over the sample space. The model extends easily—both conceptually and computationally—to the setting where we do not believe all our samples have the same common structure, where DPMHAPT allows us to learn both clustering structure and the distributional structure within each cluster.
Although we have presented HAPT in a onedimensional space for the sake of clarity, adoption of the randomized recursive partitioning scheme first introduced in [33] allows the extension of HAPT to model densities in multidimensional spaces. Variables other than simple continuous ones can also be handled naturally—all that is needed is the definition of an appropriate recursive partition. This allows inclusion of categorical and ordinalvalued variables, as well has more exotic possibilities: a continuous variable that lives on the surface of a torus, a partially ordered categorical variable, or a zeroinflated variable with a point mass at zero and a continuous component on the positive halfline.
The Polya tree’s decomposition of the density space into orthogonal Betadistributed random variables, which extends to HAPT, is central to HAPT’s computational efficiency. It also allows the performance of quick online updates in the HAPT model: when a new data point arrives we only need to update the nodes of the Polya tree which contain the new data point, rather than recomputing the entire posterior. In a HAPT truncated at a depth of levels, this means we need to reevaluate the posteriors of only nodes, rather than . HAPT may thus be used in streaming data settings where fast online updates are essential.
The DPMHAPT model may also be easily extended by replacing the Dirichlet process by any process which admits inference by a Chinese Restaurant Processtype algorithm, such as the PitmanYor process. This allows the properties of the clustering process to be adapted if the clustering assumptions implicit in the Dirichlet process are not appropriate.
References
 [1] T. W. Anderson. On the Distribution of the TwoSample Cramervon Mises Criterion. The Annals of Mathematical Statistics, 33(3):1148–1159, September 1962.
 [2] Ernesto Barrios, Antonio Lijoi, Luis E. NietoBarajas, and Igor Prünster. Modeling with Normalized Random Measure Mixture Models. Statistical Science, 28(3):313–334, August 2013.
 [3] Federico Camerlenghi, Antonio Lijoi, and Igor Prünster. Bayesian prediction with multiplesamples information. Journal of Multivariate Analysis, 156(Supplement C):18–28, April 2017.
 [4] Jayhong A Chong, José TapiaRamirez, Sandra Kim, Juan J ToledoAral, Yingcong Zheng, Michael C Boutros, Yelena M Altshuller, Michael A Frohman, Susan D Kraner, and Gail Mandel. REST: A mammalian silencer protein that restricts sodium channel gene expression to neurons. Cell, 80(6):949–957, March 1995.
 [5] M. S. Crouse, R. D. Nowak, and R. G. Baraniuk. Waveletbased statistical signal processing using hidden Markov models. IEEE Transactions on Signal Processing, 46(4):886–902, April 1998.
 [6] ENCODE Project Consortium. An integrated encyclopedia of DNA elements in the human genome. Nature, 489(7414):57–74, September 2012.
 [7] Thomas S. Ferguson. A bayesian analysis of some nonparametric problems. Ann. Statist., 1(2):209–230, 03 1973.
 [8] Thomas S. Ferguson. Prior distributions on spaces of probability measures. Ann. Statist., 2(4):615–629, 07 1974.
 [9] David A. Freedman. On the asymptotic behavior of bayes’ estimates in the discrete case. Ann. Math. Statist., 34(4):1386–1403, 12 1963.
 [10] Timothy Hanson and Wesley O. Johnson. Modeling Regression Error with a Mixture of Polya Trees. Journal of the American Statistical Association, 97(460):1020–1033, 2002.
 [11] Timothy E. Hanson. Inference for Mixtures of Finite Polya Tree Models. Journal of the American Statistical Association, 101(476):1548–1565, 2006.
 [12] Chris C. Holmes, FranÃ§ois Caron, Jim E. Griffin, and David A. Stephens. Twosample bayesian nonparametric hypothesis testing. Bayesian Anal., 10(2):297–320, 06 2015.
 [13] Alessandra Guglielmi James O. Berger. Bayesian and conditional frequentist testing of a parametric model versus nonparametric alternatives. Journal of the American Statistical Association, 96(453):174–184, 2001.
 [14] AN Kolmogorov. Sulla Determinazione Empirica di una Legge di Distribuzione. Giornale dell’Istituto Italiano degli Attuari, 4:83–91, 1933.
 [15] Charles H. Kraft. A class of distribution function processes which have derivatives. Journal of Applied Probability, 1(2):385–388, December 1964.
 [16] Michael Lavine. Some aspects of polya tree distributions for statistical modelling. Ann. Statist., 20(3):1222–1235, 09 1992.
 [17] Michael Lavine. More aspects of polya tree distributions for statistical modelling. Ann. Statist., 22(3):1161–1176, 09 1994.
 [18] Li Ma. Adaptive Shrinkage in Pólya Tree Type Models. Bayesian Analysis, 12(3):779–805, September 2017.
 [19] Steven N MacEachern. Dependent nonparametric processes. In ASA Proceedings of the Section on Bayesian Statistical Science, pages 50–55. Alexandria, Virginia. Virginia: American Statistical Association; 1999, 1999.
 [20] Steven N MacEachern. Discussion of ”The nested Dirichlet process” by A.E. Gelfand, D.B. Dunson and A. Rodriguez. Journal of the American Statistical Association, 103:1149–1151, 2008.
 [21] R. Daniel Mauldin, William D. Sudderth, and S. C. Williams. Polya trees and random distributions. Ann. Statist., 20(3):1203–1221, 09 1992.
 [22] Pietro Muliere and Stephen Walker. A Bayesian Nonparametric Approach to Survival Analysis Using Polya Trees. Scandinavian Journal of Statistics, 24(3):331–340, September 1997.
 [23] Peter Müller, Fernando Quintana, and Gary Rosner. A method for combining inference across related nonparametric bayesian models. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 66(3):735–749, 2004.
 [24] Abel RodrÃguez, David B Dunson, and Alan E Gelfand. The nested dirichlet process. Journal of the American Statistical Association, 103(483):1131–1154, 2008, http://dx.doi.org/10.1198/016214508000000553.
 [25] Lingyun Song and Gregory E. Crawford. DNaseseq: A highresolution technique for mapping active gene regulatory elements across the genome from mammalian cells. Cold Spring Harbor protocols, 2010(2):pdb.prot5384, February 2010.
 [26] Charles Stein. Inadmissibility of the usual estimator for the mean of a multivariate normal distribution. In Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability, Volume 1: Contributions to the Theory of Statistics, pages 197–206, Berkeley, Calif., 1956. University of California Press.
 [27] Yee Whye Teh. A Hierarchical Bayesian Language Model Based on PitmanYor Processes. In Proceedings of the 21st International Conference on Computational Linguistics and the 44th Annual Meeting of the Association for Computational Linguistics, ACL44, pages 985–992, Stroudsburg, PA, USA, 2006. Association for Computational Linguistics.
 [28] Yee Whye Teh and Michael I. Jordan. Hierarchical Bayesian nonparametric models with applications, pages 158–207. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, Apr 2010.
 [29] Yee Whye Teh, Michael I Jordan, Matthew J Beal, and David M Blei. Hierarchical dirichlet processes. Journal of the American Statistical Association, 101(476):1566–1581, 2006, http://dx.doi.org/10.1198/016214506000000302.
 [30] Stephen G. Walker, Paul Damien, PuruShottam W. Laud, and Adrian F. M. Smith. Bayesian Nonparametric Inference for Random Distributions and Related Functions. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 61(3):485–527, January 1999.
 [31] J. Westenberg. Significance test for median and interquartile range in samples from continuous populations of any form. Proceedings Koninklijke Nederlandse Akademie van Wetenshappen, 51:252–261, 1948.
 [32] Frank Wilcoxon. Individual Comparisons by Ranking Methods. Biometrics Bulletin, 1(6):80–83, 1945.
 [33] Wing H. Wong and Li Ma. Optional pólya tree and bayesian inference. Ann. Statist., 38(3):1433–1459, 06 2010.
Supplementary material
S1: Posterior derivation
We use the following notation in this section:

is a set in the recursive partition of implicit in the Polya tree.

and indicate the left and right children of , respectively; these are also sets in the recursive partition of .

is the number of data points across all samples that are contained in ; is the number of datapoints in sample that are contained in . Thus .
We ignore for the time being the priors on and , as they are not important for this derivation and can be reinserted at the end.
The density of the distribution generated from the HAPT model at an observation consists of three factors:

The baseline density

A term indicating how the common structure modifies the density:
Note that under the canonical representation,

A term indicating how the idiosyncratic structure of the particular sample containing modifies the density:
Note that by definition, and .
Altogether this gives us the following likelihood:
Rearranging terms, the likelihood can be written as follows:
which gives the following form for the posterior of for a particular :
Conditional on and the are Beta distributed, and we can integrate them out analytically:
We can then simply multiply by the priors on and to obtain
S2: Proofs of theoretical results
Proof of Theorem 1: Absolute continuity. The result follows directly from repeated application of Theorem 3 in [18]. From one application of that theorem we have that with probability 1; by a second application we have that for each with probability 1. Thus with probability 1.
Proof of Theorem 2: Large L1 Support. The upper level of the hierarchy () follows immediately from Theorem 4 in [18], since its prior coverage is not altered by the addition of the hierarchy. For the second level of the hierarchy, note that by Theorem 4 in [18] we have
Further, by definition of conditional expectation,
is measurable with respect to . Thus
as the integrand is positive with probability 1.
Proof of Theorem 3: Weak Consistancy. By Schwartz’s theorem, it is sufficient to show that jointly lie in the KullbackLeibler support of the prior. The proof proceeds as follows: We consider the product space of all distributions so that we can show joint convergence. We restrict ourselves to a compact set with mass , and define a set , which depends on . We show that has positive prior mass. We then show that by choosing and appropriately, we can make lie within an arbitrarily small KL ball around . This shows that are jointly in the support of the prior, and concludes the proof. We follow closely the proof of Theorem 5 in [18].
Let be the product measure on , and let be the corresponding density. Let