Sparse Density Trees and Lists: An Interpretable Alternative to High-Dimensional Histograms

1. Introduction

A histogram is a piecewise constant density estimation model. There are good reasons that the histogram is among the first techniques taught to any student dealing with data (Chakrabarti et al. 2006); (i) histograms are easy to visualize, (ii) they are accurate as long as there are enough data in each bin, and (iii) they have a logical structure that most people find interpretable. A downside of the conventional histogram is that all of these properties fail in more than two or three dimensions, particularly for binned binary or categorical data. One cannot easily visualize a conventional higher-dimensional bar plot or histogram. For binary data, this would require us to visualize bins on a high-dimensional hypercube. Worse, there may not be enough data in each bin, so the estimates would cease to be accurate. In terms of interpretability, for a higher-dimensional histogram, a large set of logical conditions characterizing each bin ceases to be an interpretable representation of the data and can easily obscure important relationships between variables. Considering all of the marginals is often useless for binary variables because there are only two bins (zero and one). Not only do histograms become uninterpretable in high dimensions, other high-dimensional density estimation methods are also uninterpretable; flexible nonparametric approaches, such as kernel density estimation (KDE), simply produce a formula, and the estimated density landscape cannot be easily visualized without projecting it to one or two dimensions, in which case we would lose substantial information. The question we ask is how to construct a piecewise constant density estimation model (like a histogram) for categorical data that have the three properties mentioned; (i) it can be visualized, (ii) it is accurate, and (iii) it is interpretable.

In this paper, we present three methods for constructing tree- and list-based density estimation models. These types of models are alternatives to bar plots or variable bin-width histograms (e.g., see Scott 1979, Wand 1997). In our models, a leaf is analogous to a histogram bin (i.e., a probability mass function bin) and is defined by conditions on a subset of variables (e.g., if $\mathbf{x}$ is the feature vector, the conditions can be “the second component of $\mathbf{x}$ is zero” and “the first component of $\mathbf{x}$ is one”), and the density is estimated to be constant with each leaf (that is, regardless of what the other components of $\mathbf{x}$ are, the density is constant). Our approaches use only a subset of the variables to define the bins, making them more interpretable.

Our density estimation models can be useful in multiple domains to detect new patterns or errors in the data. For example, Figure 1 shows the sparse density tree for the COCO-Stuff (Caesar et al. 2018) labels. The data set contains 118,000 training samples over 91 stuff categories. Although the labels are sparse, our method finds interesting combinations, such as mirror and blanket or railing, skyscraper, and ground, that are shown in Figure 2.

Let us give a second example, this time illustrating how each bin is modeled to be of constant density. We are interested in understanding the set of passengers on the Titanic. Specifically, our goal is to create a density model to understand how common or unusual the particular details of a passenger might be given three categorical features: passenger class (first, second, third, crew), whether someone is an adult (adult, child), and gender (male, female). A leaf (histogram bin) in our model might be if passenger class is third class and a person is child, then the probability of belonging to this leaf is $p(\text{state})$ is 0.038. That is, the total density in the bin where these conditions hold is 0.038, so 3.8% of passengers are children traveling in third class. The density is constant within the leaf, so for an additional variable gender not in the tree (with values male and female), each of these values would be equally probable in the leaf, with probability $0.038/2=0.019$. We described just one bin, whereas a full tree is in Figure 3.

Each of our methods aims to globally optimize a Bayesian posterior possessing a sparsity prior, which acts as a regularization term. (Unlike past work on density trees (Ram and Gray 2011), our three methods are not constructed using greedy tree induction; they are optimized instead.) In this way, Bayesian priors control the shape of the tree or list. This helps with both generalization and interpretability. For the first of our three methods, the prior parameter controls the number of leaves; its value should be set to be the user’s desired number of leaves. For the second method, the prior controls the desired number of branches for nodes of the tree. For the third method, which creates lists (one-sided trees), the prior controls the desired number of leaves and also, the length of the leaf descriptions.

The structure of these sparse tree models aims to fix the three issues with conventional histograms: (i) visualization (we need only write down the conditions we used in the density tree (or density rule list) to visualize the model); (ii) accuracy (the prior encourages the density estimation model to be smaller, which means the bins are larger and generalize better); and (iii) interpretability (the prior encourages sparsity and encourages the model to obey a user-defined notion of interpretability).

In what follows, we provide related works, and then, we describe our three methods in Section 3. In Section 4, we discuss how we optimized the posteriors for three methods. Section 5 provides experiments and examples, Section 6 provides a run-time analysis and a study using simulated data sets. Section 7 provides a consistency proof, and we conclude in Section 8.

2. Related Work

2.1. Nonparametric Density Estimation

Density estimation is a classic topic in statistics and unsupervised machine learning. Without using domain-specific generative assumptions, the most useful techniques have been nonparametric, mainly variants of KDE (Akaike 1954, Rosenblatt et al. 1956, Parzen 1962, Cacoullos 1966, Nadaraya 1970, Rejtö and Révész 1973, Silverman 1986, Devroye 1991, Wasserman 2006, Mahapatruni and Gray 2011, Cattaneo et al. 2020, Varet et al. 2023). KDE is highly tunable, is not domain dependent, and can generalize well, but it does not have the interpretable logical structure of histograms. Similar alternatives include mixtures of Gaussians (Ormoneit and Tresp 1996, 1998; Zhuang et al. 1996; Li and Barron 2000; Chen et al. 2006; Seidl et al. 2009), forest density estimation (Liu et al. 2011), RODEO (Liu et al. 2007), and other nonparametric Bayesian methods (Müller and Quintana 2004), which have been proposed for general purpose (not interpretable per se) density estimation. Jebara (2008) provides a Bayesian treatment of latent directed graph structure for non independent and identically distributed (i.i.d.) data but does not focus on sparsity. Pólya trees have been generated probabilistically for real-valued features (Wong and Ma 2010) and could be used as priors for our method. Friedman et al. (1984) uses a projection pursuit method to perform density estimation. Another task related to density estimation is level set estimation, where the goal is to determine whether the density at a leaf is higher than a prespecified value γ; Willett and Nowak (2007) address the problem using tree representations, and Holmström et al. (2015) use a discretized kernel to construct level set trees. Some works (e.g., Sasaki and Hyvärinen 2018, Liu et al. 2021) use neural networks to perform density estimation that do not aim to be interpretable. In Luo et al. (2019), a smoothing spline is used to perform density estimation. In Rehn et al. (2018), a nonparametric density estimator called “FRONT” segments a data stream through a periodically updated linear transformation.

The most closely related paper to ours is on density estimation trees (DETs) (Ram and Gray 2011) and their extensions. DETs are constructed in a top-down greedy way. This gives them a disadvantage in optimization, often leading to lower-quality trees. They also do not have a generative interpretation, and their parameters do not have a physical meaning in terms of the shape of the trees (unlike the methods defined in this work). DET was used by Wu et al. (2018), leveraging ideas from Lu et al. (2013) with random forests to perform density estimation. DET has also been applied to high-energy physics (Anderlini 2015). Techniques to avoid overfitting in tree-based density estimation models have been discussed by Anderlini (2016). Other top-down greedy approaches (e.g., Yang and Wong 2014a, b; Li et al. 2016) use discrepancy, negative log likelihood, or MISE (Ooi 2012) as splitting criteria. A distinction between our work and existing work is that we place priors directly on the shape of the trees that we desire using a Bayesian approach. We do not rely on greedy splitting and pruning; the splits are optimized instead.

2.2. Bayesian Tree Models

Bayesian tree models are commonly used for tasks other than density estimation (i.e., classification and regression). Some examples include Bayesian CART (Wu et al. 2007), Bayesian additive regression trees (BARTs) (Chipman et al. 2010), and Bayesian rule lists (Letham et al. 2015, Yang et al. 2017). Bayesian CART and BART use priors that specify the probability that a node is terminal and a uniform probability distribution over the choices for a split. Our priors function differently. Often, we have priors over global properties of the trees, such as the number of total leaves (our method I and method III). Also, we have prior parameters governing the number of branches at a node (our method II), which is different from Bayesian CART or BART where there are only two branches at every node. Our rule list density approach (method III) has a prior on the number of conditions used in each rule, which is similar to Bayesian rule lists but not similar to Bayesian CART or BART, which have only one condition defining each split.

3. Methods

We use a Bayesian approach to achieve sparsity by introducing priors on the shape of the trees. In particular, in method I, we define a prior on the number of leaves in the tree, then calculate the likelihood of the data having been generated by a particular tree, and multiply the prior and the likelihood to create a posterior to be optimized over all trees. In method II, we instead choose a prior over the number of branches for each split in the tree, preferring a small number of branches. In method III, we switch to rule lists, where the prior prefers models with a smaller number of rules and a smaller number of conjunctions per rule.

Before the introduction of the three methods, we first present notation. We will focus on the problem of estimating the unknown distribution f with tree-structured approximations given a set of n data points $X=\{x_i,\dots ,x_n\}$ drawn i.i.d. from f on $\mathcal{X}\subset {\mathbb{R}}^p$.

For the tree-structured estimations, we express the path to a leaf as the set of conditions on each of p features along the path. For instance, for a particular leaf (leaf t in Figure 4), we might see conditions on the first feature (e.g., $x_{.1}\in \{4,5,6\}$) and the second feature ($x_{.2}\in \{100,101\}$). Thus, the leaf is defined by the set of all feature values that obey these conditions: that is, the leaf could be

This implies that there is no restriction on $x_{.3},x_{.4},\dots ,x_{.p}$ for observations within the leaf. Notationally, a condition on the jth feature is denoted $x_{.j}\in {\sigma }_j(l)$, where ${\sigma }_j(l)$ is the set of allowed values for feature j along the path to leaf l. If there are no conditions on feature j along the path to l, then ${\sigma }_j(l)$ includes all possible values for feature j. Thus, leaf l includes feature values x obeying:

For categorical data, the volume in a leaf l, which is required for computing the density, is calculated as

We give example of this volume computation next.

The data are categorical. Possible values for $x_{.1}$ are $\{1,2,3,4,5,6,7\}$. Possible values for $x_{.2}$ are $\{100,101,102,103\}$. Possible values for $x_{.3}$ are $\{10,11,12,13,14,15\}$, and those for $x_{.4}$ are $\{8,9,10\}$. Consider the tree in Figure 4.

We compute the volume for leaf l. Here, ${\sigma }_1(l)=\{4,5\}$ because l requires both $x_{.1}\in \{4,5,6\}$ and $x_{.1}\in \{4,5\}$. ${\sigma }_2(l)=\{100,101\},$${\sigma }_3(l)=\{10,15\},$ and ${\sigma }_4(l)=\{8,9,10\}$ because there is no restriction on $x_{.4}$. So, ${\mathbb{V}}_l={\prod }_j|{\sigma }_j(l)|=2·2·2·3=24.$

Our notation handles only categorical data (for ease of exposition) but can be extended to handle ordinal and continuous data. For ordinal data, the definition is the same as for categorical, but σ_j can (optionally) include only contiguous values (e.g., $\{3,4,5\}$ but not $\{3,4,6\}$). For continuous variables, σ_j is the “volume” of the continuous variables: for example, for node condition $x_{.j}\in (0,0.5),{\sigma }_j=0.5-0$.

In the next subsections, we present the three methods: the leaf-sparse modeling approach, the branch-sparse modeling approach, and the density list approach. We define the prior and likelihood for each of the three modeling methods, combining them to get posteriors.

3.1. Method I: Leaf-Sparse Density Trees

Let us define the prior, likelihood, and posterior for method I.

3.1.1. Prior.

For this modeling method, the main prior on tree T is on the number of leaves K_T. This prior we choose to be Poisson with a particular scaling (which will make sense later on), where the Poisson is centered at a user-defined parameter λ, which is the user’s desired number of leaves in the tree prior to seeing data. Notation $N_{K_T}$ is the number of trees with K_T leaves. It can be calculated directly given K_T. The prior over trees is

Thus, λ allows the user to control the number of leaves in the tree. The number of possible trees is finite; thus, the distribution can be trivially normalized. Among trees with K_T leaves, tree T is chosen uniformly, with probability $1/N_{K_T}$. This means that the probability to choose a particular tree T is Poisson:

We place a uniform prior over the probabilities for a data point to land in each of the leaves. To do this, we start from a Dirichlet distribution with equal parameters ${\alpha }_1=\cdots ={\alpha }_{K_T}=\alpha \in {\mathbb{Z}}^{+}$, where $\alpha >1$. The hyperparameter α is a pseudocount that is typically chosen to be a small number (one or two) to avoid a zero value for the estimated densities. We denote the vector with K_T equal entries $[\alpha ,\dots ,\alpha ]$ as ${\mathit{\alpha }}_{K_T}$. We draw multinomial parameters $\mathit{\theta }=[{\theta }_1,\dots ,{\theta }_{K_T}]$ from Dir$({\mathit{\alpha }}_{K_T})$, which govern the prior on the popularity of each leaf. Thus, the prior on feature values given the tree structure is

Here, $\text{Beta}({\mathit{\alpha }}_{K_T})=\frac{{\prod }_{l=1}^{K_T}\mathrm{\Gamma }({\alpha }_{K_T,l})}{\mathrm{\Gamma }({\sum }_{l=1}^{K_T}{\alpha }_{K_T,l})}=\frac{{(\mathrm{\Gamma }(\alpha ))}^{K_T}}{\mathrm{\Gamma }(K_T\alpha )}$ is the multinomial beta function, which is also the normalizing constant for the Dirichlet distribution.

Thus, the first part of our generative modeling method is as follows given hyperparameters λ and α. The number of leaves in T is $K_T\propto \text{scaled}\hspace{0.17em}\text{Poisson}(\lambda )\hspace{0.17em}(\text{i}.\text{e}.,N_{K_T}\times \text{Poisson}(K_T,\lambda )$), the tree shape is T $\propto \text{Uniform}\hspace{0.17em}\text{over}\hspace{0.17em}\text{trees}\hspace{0.17em}\text{with}\hspace{0.17em}{K}_T \text{leaves}$, and $\text{Prior}\hspace{0.17em}\text{distribution}\hspace{0.17em}\text{over}\hspace{0.17em}\text{leaves}: \mathit{\theta }\propto \text{Dir}({\mathit{\alpha }}_{K_T}).$ As usual, the prior is regularization that can be overwhelmed with a large amount of data. However, it takes a substantial amount of data and an enormous number of feature-value pairs to reduce the prior’s influence on the posterior.

3.1.2. Likelihood.

Here, we calculate the likelihood of the data to have arisen from a particular tree. Let n_l denote the number of points captured by the lth leaf, and denote ${\mathbb{V}}_l$ as the volume of that leaf, defined in (1). The probability to land at any specific value within leaf l is $\frac{{\theta }_l}{{\mathbb{V}}_l}$. The likelihood for the full data set is thus

3.1.3. Posterior.

The posterior, which is (by definition) the likelihood times the prior, can be written as follows, where we have substituted the distributions from the prior into the formulas:

where $P(T|\lambda )$ is simply Poisson($K_T,\lambda$) as discussed earlier. For numerical stability, we maximize the log posterior, which is equivalent to maximizing the posterior.

As compared with full Bayesian approaches, by maximizing the posterior, we leverage relatively faster computation time and optimize for a single model, which can be important for interpretability. However, in turn, we miss out on uncertainty information, which we could get by modeling the posterior density.

For the purposes of prediction, we are required to estimate the density that is being assigned to leaf l. This is calibrated to the data simply as $\widehat{f}=n_l/(n{\mathbb{V}}_l)$, where n is the total number of training data points and n_l is the number of training data points that reside in leaf l. The formula implicitly states that the density in the leaf is uniformly distributed over the features whose values are undetermined within the leaf (features for which σ_j contains all possible values for feature j). That is, the probability mass function is the same for any point within the leaf.

3.2. Method II: Branch-Sparse Density Trees

In method I, we regularized the number of leaves of the density tree. In method II, we instead regularize the number of branches at each node to simplify the model by avoiding models with too many branches. In other words, method I aims to have a small number of leaves. Method II aims to have a small number of branches.

In method I, a Dirichlet distribution is drawn only over the leaves. In method II, a Dirichlet distribution is drawn at every internal node to determine branching. Similar to the previous method, we choose the tree that optimizes the posterior.

The prior is composed of two pieces: the part that creates the tree structure and the part that determines how data propagate through it.

3.2.1. Tree Structure Prior.

For tree T, we let $B_T=\{b_i|i\in I\}$ be a multiset, where each element is the count of branches from a node of the tree. For instance, if in tree T with three nodes, the three nodes have three branches, two branches, and two branches, respectively, then $B_T=\{3,2,2\}$. We let $N_{B_T}$ denote the number of trees with the same multiset B_T. Note that B_T is unordered, so {3, 2, 2} is the same multiset as {2, 3, 2} or {2, 2, 3}.

Let I denote the set of internal nodes of tree T, and let L denote the set of leaves. As before, we let ${\mathbb{V}}_l$ denote the volume of leaf l.

In the generative model for the branch-sparse density tree, a Poisson distribution with parameter λ is used at each internal node in a top-down fashion to determine the number of branches. Iteratively, for node i, the number of branches, b_i, obeys $b_i\sim \text{Poisson}(\lambda ),$ where the parameter λ is the desired number of branches (before seeing data). Hence, at any node i, with probability $\exp (-\lambda )\frac{{\lambda }^{b_i}}{b_i!}$, there are b_i branches from node i. This implies that with probability $\exp (-\lambda )$, the node is a leaf. In summary,

Among trees with multiset B, tree T is chosen uniformly, with probability $\frac{1}{N_B}.$ This means that the probability to choose a particular tree shape is

3.2.2. Tree Propagation Prior.

After the tree structure is determined, we need a generative process for how the data propagate through each internal node. We denote θ_l as the probability to land in leaf l. We denote ${\tilde{\theta }}_{ij}$ as the probability to traverse to node j from internal node i. Notation $\mathit{\theta }$ is the vector of leaf probabilities (the θ_l’s), $\tilde{\mathit{\theta }}$ is the set of all ${\tilde{\theta }}_{ij}$’s, and ${\tilde{\mathit{\theta }}}_i$ is the set of all internal node transition probabilities from node i (the ${\tilde{\theta }}_{ij}$’s).

We compute $P(\widetilde{{\mathit{\theta }}_i}|\alpha ,T)$ for all internal nodes i of tree T. As before, α is a pseudocount to avoid zero-valued estimated densities. At each internal node, we draw a sample from a Dirichlet distribution with parameter $[\alpha ,\dots ,\alpha ]$ (of size equal to the number of branches b_i of i) to determine the proportion of data, ${\tilde{\theta }}_{i,j}$, that should go along the branch leading to each child node j from the internal parent node i. Thus, $\widetilde{{\mathit{\theta }}_i}\sim \text{Dir}(\mathit{\alpha })$ for each internal node i: that is,

where ${\text{Beta}}_k(\mathit{\alpha })$ is the normalizing constant for the Dirichlet distribution with parameter $\alpha$ and k categories and C_i are the indices of the children of i. Thus,

Thus, the prior is $P(T|\lambda )·P(\tilde{\mathit{\theta }}|\alpha ,T)$, where $P(T|\lambda )$ is in (2) and $P(\tilde{\mathit{\theta }}|\alpha ,T)$ is in (3).

In summary, the prior of method II is as follows given hyperparameters λ and α:

• multiset of branches $B_T\propto N_{B_T}\left[{\prod }_{i\in I}{e}^{-\lambda }\frac{{\lambda }^{b_i}}{b_i!}\right]\left[{\prod }_{l\in L}{e}^{-\lambda }\right],$

• tree shape $T\sim \text{Uniform}\hspace{0.17em}\text{over}\hspace{0.17em}\text{trees}\hspace{0.17em}\text{with}\hspace{0.17em}\text{branches}\hspace{0.17em}{B}_T,$ and

• prior distribution over each branch ${\tilde{\mathit{\theta }}}_i\sim \text{Dir}(\alpha ).$

The likelihood is the same as that for method I.

3.2.3. Posterior.

The posterior is proportional to the prior times the likelihood terms. Here, we are integrating over the ${\tilde{\mathit{\theta }}}_i$ terms for each of the internal nodes i:

where $c_1,\dots ,c_{b_i}\in C_i$ in the second to last expression. We used the equation ${\sum }_{i\in I}{b}_i=|L|+|I|-1$ for a tree in the last line.

3.2.4. Possible Extension.

We can include an upper layer of the hierarchical Bayesian model to control (regularize) the number of features d (of a total of p dimensions) that are used in the model. This would introduce an extra multiplicative factor within the posterior of $(\begin{array}{c}p\\ d\end{array}){\gamma }^d{(1-\gamma )}^{p-d}$, where γ is a parameter between zero and one and a smaller value favors a simpler model. γ corresponds to the probability that a feature is chosen to be included in the model. For example, the value 0.5 corresponds to the case where the user prefers to have half of the features to be chosen. The posterior would become

3.3. Method III: Sparse Density Rule List

Rather than producing a general tree, an alternative approach is to produce a rule list, which is a one-sided tree. Rule lists are easier to optimize than trees. Each tree can be expressed as a rule list by creating a rule for each leaf, where the conditions defining the leaf also define the rule. By using lists, we implicitly hypothesize that the full space of trees may not be necessary and that simpler rule lists may suffice.

An example of a sparse density rule list is as follows: if x obeys a₁then density$(x)=f_1$ else if x obeys a₂then density$(x)=f_2$$\dots$ else if x obeys a_m then density$(x)=f_m$ else density$(x)=f_0$.

Here, as with the trees, the density is the probability mass function, which is constant for the entire portion of the feature space that falls into the leaf.

The antecedents a₁,…,a_m are Boolean assertions that are either true or false for each data point x_i. They are chosen from a large premined collection of possible antecedents called A. We define A to be the set of all possible antecedents of size at most H, where the user chooses H. The size of A is $|A|={\sum }_{j=0}^H{A}_j,$ where A_j is the number of antecedents of size j:

where feature i consists of q_i categories.

For example, say the features consist of two, three, four, and five categories, respectively. If H = 2, then the total number of elements of A is $|A|=$ 1 (no feature is chosen) + 2 (because there are two categories for the first feature) + 3 (because there are three categories for the second feature) + 4 (because there are four categories from the third feature) + 4 (because there are five for the fourth feature) + six (possible combinations of feature 1 and feature 2) + $\dots$ + 20 (possible combinations of feature 3 and feature 4).

3.3.1. Generative Process.

We now sketch the generative process for the tree from the observations $X=\{x_i\}$ and antecedents $A=\{a_j\}$. Prior parameter λ is the user’s preference of the length of the density list (in the absence of data), and η is the user’s preference for the number of conjunctions in each subrule a_j.

Define $a_{<j}$ as the antecedents before j in the rule list if there are any. For example $a_{<3}=\{a_1,a_2\}$. Similarly, let c_j be the cardinalities of the antecedents before j in the rule list. Let d denote the rule list. Following the exposition of Letham et al. (2015), we use a prior over rule lists to encourage sparsity. The generative process is described in Algorithm 1. It depends on the input prior parameters $\lambda ,\eta$, and $\mathit{\alpha }$, which is a user-chosen vector of size m + 1 (as before, usually all elements in $\mathit{\alpha }$ are the same and indicate pseudocounts).

Algorithm 1

(Density Rule Lists Generation Procedure)

Input: Prior parameters λ and η, pseudocount $\mathit{\alpha }$, observations $X=\{x_i\}$, antecedents $A=\{a_j\}$

Output: Density rule list

• 1: Sample a decision list length $m\sim P(m|A,\lambda )$

• 2: for decision list rule $j=1,\dots ,m$ do

• 3:  Sample the cardinality of antecedent a_j in d as $c_j\sim P(c_j|c_{<j},A,\eta )$.

• 4:  Sample a_j of cardinality c_j from $P(a_j|a_{<j},c_j,A)$.

• 5: end for

• 6: for observation $i=1,\dots ,n$ do

• 7:  Find the antecedent a_j in d that is the first that applies to x_i.

• 8:  If no antecedents in d applies, set j = 0.

• 9: end for

• 10: Sample parameter $\mathit{\theta }\sim$ Dirichlet ($\mathit{\alpha }$) for the probability to be in each of the leaves

• 11: for each leaf l in the rule list do

• 12:  Compute volume V_l according to (1)

• 13:  Compute density $f_l=\frac{{\theta }_l}{{\mathbb{V}}_l}$

• 14: end for

For instance, in step 1 of Algorithm 1’s generation process, we might find out that the rule list is of size m = 4. Then, in steps 2–5, we sample the cardinality of each rule, so we may find that the four rules are of cardinality $2,1,2,3$. Then, in steps 6–9, we sample to determine which rules are actually used; for instance, the first rule might be $a_1:\{x_{14}=\text{blue}\hspace{0.17em}\hspace{0.17em}\text{\&}\hspace{0.17em}\hspace{0.17em}{x}_{21}=\text{large}\}$, the second rule could be $a_2:\{x_{52}=\text{feathered}\}$, the third rule might be $a_3:\{x_{13}=\mathrm{10,110}\hspace{0.17em}\hspace{0.17em} \&\hspace{0.17em}\hspace{0.17em}{x}_3=\text{angled}\}$, and the fourth rule might be $a_4:\{x_{55}=000\hspace{0.17em}\hspace{0.17em}\&\hspace{0.17em}\hspace{0.17em}{x}_8=\text{hypothetical}\&\hspace{0.17em}\hspace{0.17em}{x}_{21}=\text{small}\}$. The rest of the feature space (that does not fall into any of these rules) would go to the default rule, again having constant density within the rule.

3.3.2. Prior.

The distribution of m in step 1 is the Poisson distribution truncated at the total number of preselected antecedents:

When $|A|$ is huge, we can use the approximation $P(m|A,\lambda )\approx {\lambda }^m/m!$ as the denominator would be close to one.

For step 2, we let $R_j(c_1,\dots ,c_j,A)$ be the set of antecedent cardinalities that are available after drawing antecedent j, and we let $P(c_j|c_{<j},A,\eta )$ be a Poisson truncated to remove values for which no rules are available with that cardinality:

We use a uniform distribution over antecedents in A of size c_j excluding those in a_j:

The sparse prior for the antecedent lists is thus

The prior distribution over the leaves $\mathit{\theta }=[{\theta }_1,\dots ,{\theta }_m,{\theta }_0]$ is drawn from Dir(${\mathit{\alpha }}_{m+1}$):

It is straightforward to sample an ordered antecedent list d from the prior by following the generative process that we just specified, generating rules from the top down.

3.3.3. Likelihood.

Similar to method I, the probability to land at any specific value within leaf l is $\frac{{\theta }_l}{{\mathbb{V}}_l}$. Hence, the likelihood for the full data set is thus $P(X|\mathit{\theta },d)={\prod }_{l=0}^m{\left(\frac{{\theta }_l}{{\mathbb{V}}_l}\right)}^{n_l}.$

3.3.4. Posterior.

The posterior can be written as

where the last equality uses the standard Dirichlet-multinomial distribution derivation.

4. Numerical Methods to Optimize the Objective Function

In the previous section, we presented the posterior functions for three generative modeling methods. Because the search space of our problems is large, we use simulated annealing, a metaheuristic optimization algorithm that allows us to approximate global solutions. More specifically, in this section we describe a simulated annealing scheme that we implemented in order to find the optimal tree that maximizes the posterior for method I and method II as well as discuss method III’s optimization details.

4.1. Simulated Annealing for Tree-Based Methods

A successful simulated annealing scheme requires us to create a useful definition of a neighborhood. We define our neighborhood such that each move explores a neighboring tree where we are able to extend or shrink the tree.

At each iteration, we need to determine which neighboring tree to move to. To decide which neighbor to move to, we fix a parameter $ϵ>0$ beforehand, where ϵ is small (approximately 0.01 in our experiments). The parameter ϵ is the probability that we will perform a structural change to jump out of a possible local minimum. All other actions are taken with equal probability. Thus, at each time, we generate a number from the uniform distribution on (0, 1). Then, we do the following.

1. (Shrink at leaf) If the number is smaller than $\frac{1-ϵ}{4}$, we select uniformly at random a “parent” node that has leaves as its children and remove its children. This is always possible unless the tree is the root node itself, in which case we cannot remove it and this step is skipped.

2. (Expand) If the random number is between $\frac{1-ϵ}{4}$ and $\frac{1-ϵ}{2}$, we pick a leaf randomly and a feature randomly. If it is possible to split on that feature, then we create children for that leaf. (If the feature has been used up by the leaf’s ancestors, we cannot split, and we then skip this round.)

3. (Regroup) If the random number is between $\frac{1-ϵ}{2}$ and $\frac{3(1-ϵ)}{4}$, we pick a node randomly, delete its descendants, and split the node, creating two child nodes where the splitting is done on subsets of the node’s feature values. Sometimes this is not possible: for example, if we pick a node where all the features have been used up by the node’s ancestors or if the node has only one category. In that case, we skip this step.

4. (Merge sibling nodes) If the random number is between $\frac{3(1-ϵ)}{4}$ and $(1-ϵ)$, we choose two nodes that share a common parent, delete all their descendants, and merge the two nodes (e.g., black, white, red, and green become black or white, red, and green).

5. (Structural change) If the random number is more than $1-ϵ$, we perform a structural change operation, where we remove all the children of a randomly chosen node of the tree.

See Algorithm 2 for the pseudocode of the simulated annealing procedure. The last three actions avoid problems with local minima. The algorithms can be warm started using solutions from other algorithms (e.g., DET trees). We found it useful to occasionally reset to the best tree encountered so far or the trivial root node tree.

Algorithm 2

(Simulated Annealing for Tree-Based Methods)

Input: Prior parameters, θ, maximum number of iterations, N, ϵ, cooling schedule Cool(iteration) for the simulated annealing

Output: Optimal density tree

• 1: Initialize the initial tree T to be a single node, compute the posterior using the objective function and the prior parameters. Set iteration number to be zero.

• 2: while iteration number $<N$ do

• 3:  Draw a random number r, from Uni(0, 1).

• 4:  if $r<\frac{1-ϵ}{4}$ then

• 5:   Perform shrink at leaf operation.

• 6:  else if $r<\frac{2(1-ϵ)}{4}$ then

• 7:   Perform expand operation.

• 8:  else if $r<\frac{3(1-ϵ)}{4}$ then

• 9:   Perform regroup operation.

• 10:  else if $r<1-ϵ$ then

• 11:   Perform merge sibling nodes operation.

• 12:  else

• 13:   Perform structural change.

• 14:  end if

• 15:  Compute the objective value for modified tree.

• 16:  if a better objective value is obtained than current best then

• 17:   Update T to be the current tree.

• 18:  end if

• 19:  if the current tree is worse than the current best tree T then

• 20:   With probability defined by the cooling schedule Cool(iteration), update T to be the current tree. This will always be a small probability.

• 21:  end if

• 22: end while

• 23: Return T

4.2. Sparse Density Rule List Optimization

To search for optimal density rule lists that fit the data, we use local moves (adding rules, removing rules, and swapping rules) and use the Gelman–Rubin convergence diagnostic applied to the log posterior function.

A technical challenge that we need to address in our problem is the computation of the volume of a leaf. Volume computation is not needed in the construction of a decision list classifier like that of Letham et al. (2015), but it is needed in the computation of density list. There are multiple ways to compute the volume of a leaf of a density rule list.

4.2.1. Approach 1: Analytical Computation.

Use the inclusion-exclusion principle to directly compute the volume of each leaf. Consider computing the volume of the ith leaf in a density rule list. Let $V_{a_i}$ denote the volume induced by the rule a_i: that is, the number of points in the domain that satisfy a_i. To belong to that leaf, a data point has to satisfy a_i and not earlier rules $a_{<i}$. Hence, the volume of the ith leaf is equal to the volume obeying a_i alone minus the volume that has been used by earlier rules. Using notation $a_k^c$ to denote the complement of condition a_k, we have the following:

where the last expression is because of the inclusion-exclusion principle, and it only involves the volume resulting from conjunctions. The volume resulting from conjunctions can be easily computed from data. Without loss of generality, suppose we want to compute the volume of $V_{a_1\wedge \cdots \wedge a_k}$. For each feature that appears, we examine if there is any contradiction between the rules; for example, if feature 1 is present in both a₁ and a₂ where rule a₁ specifies feature 1 to be zero, whereas a₂ specifies feature 1 to be one, then we have found a contradiction, and the volume of the intersection of a₁ and a₂ should be zero. If there is no contradiction, then the volume is equal to the product of the number of distinct categories of all the features that are not used. If all features are used, then the volume is one. By using the inclusion-exclusion principle, we reduce the problem to just computing a volume of conjunctions as in (4). Note that for this approach, we still need to iterate over all conjunctions for each volume computation, which can be computationally expensive.

4.2.2. Approach 2: Uniform Sampling.

Create uniform data over the whole domain, and count the number of points that satisfy the antecedents. This approach would be expensive when the domain is huge but easy to implement for smaller problems.

4.2.3. Approach 3: MCMC.

Use an Markov chain Monte Carlo (MCMC) sampling approach to sample the whole domain space. This approach is again not practical when the domain size is huge as the number of samples required will increase exponentially because of the curse of dimensionality.

We use the analytical computation approach 1 in our implementation.

Some works (Angelino et al. 2018, Rudin and Ertekin 2018) have achieved provable optimality on minimization of objectives over a set of premined rules, although for supervised classification (not density estimation). They have also noted that randomized methods, such as those of Yang et al. (2017) (or those considered in the present work), tend to produce models that are close to these optimal solutions, leading us to believe that our search methods may actually achieve close-to-optimal solutions fairly often, particularly for problems where there may be a large “Rashomon” set of good models (Fisher et al. 2019, Semenova et al. 2022). None of these earlier works consider density estimation, but nonetheless, we have reason to hypothesize that simulated would also yield near-optimal models for density estimation.

5. Experiments

Our experimental setup is as follows. We considered five methods: the leaf-sparse density tree, the branch-sparse density tree, the sparse density rule list, regular histograms, and DETs (Wu et al. 2018). Note that the implementation of DET is meant for continuous variables, but we use it anyway for comparison. To our knowledge, these methods can serve to represent the full set of logic-based, high-dimensional density estimation methods. To ascertain uncertainty, we split the data in five folds and assessed test log likelihood (i.e., out-of-sample performance) and sparsity of the trees for each method on every fold. A model with fewer bins and higher test likelihood is a better model.

Let us discuss how the baselines were implemented. For the standard high-dimensional histogram baseline, we treated each possible set of feature values (e.g., $x_{.1}=1,x_{.2}=0,\dots ,x_{.10}=1$) as a separate bin. (We call a set of feature values a configuration; it is a point in our feature space.) Histograms have the disadvantage that they create a large number of bins and thus, may not generalize well to the test set; they are also not interpretable because they cannot be visualized in a tree or list.

DET was designed for continuous data, which meant that the computation of volume needed to be adapted for discrete data; it is the number of configurations in the bin rather than the lengths of each bin multiplied together. Thus, we used our own computations for the density in each leaf following volume computations in (1).

The DET method has two parameters: the minimum allowable support in a leaf and the maximum allowable size of a leaf. We originally planned to use a minimum of zero and a maximum of the size of the full data set, but the algorithm often produced trivial models when we did this (i.e., models with one leaf). Therefore, we tried values $\{1,3,5\}$ for the minimum size of a leaf and $\left\{5,10,50,100,\lfloor \frac{n}{10}\rfloor ,\lfloor \frac{n}{5}\rfloor \right\}$ for the maximum size of a leaf, where n is the number of training data points. We use nested crossvalidation over five folds. For each fold, we optimize parameters for validation log likelihood and report out-of-distribution log likelihood. DET has the disadvantage of being a greedy algorithm, and the available implementation of DET is not designed for categorical data (see Appendix B); thus, DET may not produce trees that are as useful or sparse as those from method I, II, or III.

For the data sets in the following two sections for the leaf-sparse density tree model (method I), the mean of the Poisson prior was chosen from the set $\{5,8,10\}$ using nested crossvalidation. For the branch-sparse density tree model (method II), the parameter to control the number of branches was chosen from the set $\{2,3\}$. The parameter α was set to be two for the experiment. This corresponds to a pseudocount of two data points in each bin (to prevent bins with zero data points). For the sparse density rule list (method III), the parameters $\lambda ,\eta$, and α were chosen among $[3,5,7],1$, and 1, respectively. We provide a summary of parameters and their suggested values in Appendix D.

5.1. Titanic Data Set

As discussed earlier, a sparse density tree or list would help us understand the distribution of people on board the Titanic. The Titanic data set has an observation for each of the 2,201 people aboard the Titanic. There are three features: gender, whether someone is an adult, and the class of the passenger (first class, second class, third class, or crew member) (Figures 5 and 6).

Figure 6 shows the results, both for out-of-sample log likelihood (on the y axis) and sparsity (on the x axis), for each method for each of the five test folds. The histogram method had the most leaves (by design) and thus, tended to overfit. Our methods performed well; arguably, the sparse density rule list method performed slightly better in the likelihood-sparsity trade-off. In general, the results are consistent across folds; the histogram produces too many bins, the sparse density rule list and density tree methods performs well, and DET has worse log likelihood.

Figure 3 shows one of the density models generated by the leaf-sparse density tree method. The reason for the split is clear; there were fewer children than adults, the distributions of the males and females were different (mainly because of the fact that the crew was mostly male), and the number of crew members was very different than the numbers of first, second, and third class passengers.

Figure 5 shows density rule lists for the Titanic data set by choosing parameter λ, which indicates preferred list length, from set {2, 4, 7}; we change the density rule lists as well as its length. Lower values of the parameter lead to shorter density rule lists, whereas larger preferred lengths correspond to longer density rule lists.

5.2. Crime Data Set

The house break-in data used for this experiment are from the Cambridge Police Department, Cambridge, Massachusetts. The motivation is to understand the common types of modus operandi characterizing house break-ins, which is important in crime analysis. The data consist of 3,739 separate house break-ins occurring in Cambridge between 1997 and 2012 inclusive. The six categorical features for the crime data set are as follows: (1) location of entry (“window,” “door,” “wall,” and “basement”); (2) means of entry (“forceful” (cut, broke, cut screen, etc.), “open area,” “picked lock,” “unlocked,” and “other”); (3) whether the resident was inside; (4) whether the premise was judged to be ransacked by the reporting officer; (5) whether the entry happened on a “weekday” or “weekend”; and (6) type of premise. The first category is “residence” (including apartment, residence/unknown, dormitory, single-family house, two-family house, garage (personal), porch, apartment hallway, residence unknown, apartment basement, condominium). The second category is nonmedical, nonreligious “workplace” (commercial unknown, accounting firm, research, school). The third group “medical” consists of halfway houses, nursing homes, medical buildings, and assisted living. The fourth group “parking” consists of parking lots and parking garages, and the fifth group “social” consists of Young Women’s Christian Associations, Young Men’s Christian Associations, and social clubs. The last groups are “storage,” “construction site,” “street,” and “church,” respectively.

The experiments in Figure 6(b) show that our approaches dominate DET for the crime data set and are sparser than DET trees. The standard histogram’s results were not reported because they involve too many bins (1,440) to fit on the figure, and they are thus not competitive.

One of the trees obtained from the leaf-sparse density tree method is in Figure 7. It states that most burglaries happen at residences; the nonresidential density has values less than 2 × 10⁻⁴. Given that a crime scene is a residence, most crimes happened on weekends. For residential crimes, burglary is more likely to happen when the owner is not inside (density 0.0046 if weekday and 0.2 if weekend, the premise is judged to be not ransacked, and there is forceful entry through the door or window). When the premise is judged to be ransacked, the crime is more likely to happen with the door as the location of entry (density 0.0042) compared with wall, window, and basement (density 4.72 × 10⁻⁴). On weekends, for residential and not-ransacked premises, doors and windows are more common locations of entry. If the entry is not forceful, unlocked windows and doors are the most common means of entry (density is 0.0361). If the means of entry is picked lock, the density is 0.0011, and if the area is open, the density is 0.0068.

Important aspects of the modus operandi are within the leaves of the tree: for instance, that the owner of a residence was not inside, that the house was not ransacked, and that the entry was forceful through the door or window. If this approach would have been performed using a regular histogram, it would require 1,440 different markers (discrete states), whereas the crime tree groups the crimes into just 20 bins.

These types of results can be useful for crime analysts to assess whether a particular modus operandi is unusual. For instance, according to the tree, it is clearly more unusual for the owner to be inside during the break-in, as shown by the smaller density values in the leaves when the owner is inside. Also, according to the densities in the leaves, it is more common for the means of entry to be forceful and for the location of entry to be windows and doors. A density list for these data is presented in Figure 8. The preferred length of the list was chosen from the set {3, 5, 7}.

6. Empirical Performance Analysis

The experiments are designed to provide insight into how the methods operate.

6.1. Run-Time Analysis

We studied the performance of the density estimation methods on data sets of different complexity. We chose 17 data sets, including financial data sets (Bank-Full, Telco Customer Churn, home equity line of credit (HELOC)), recidivism risk score data (COMPAS) (Larson et al. 2016), UCI repository data sets (e.g., Car, Mushroom, U.S. Census data) (Markelle Kelly, Rachel Longjohn, Kolby Nottingham, The UCI Machine Learning Repository, https://archive.ics.uci.edu), and detection labels of COCO-Stuff + thing data. The complexity of the data set, in this case, is defined based on the number of samples and/or the number of feature-value pairs (the sum of categories for each feature). For the considered data sets, the number of samples ranged from 625 to around 2.5 million. The number of feature-value pairs ranged from 9 to around 400, and there were from 3 to 68 features. See Table C.1 in Appendix C for data set statistics and the preprocessing steps taken.

For data sets with fewer than 100,000 samples and fewer than 200 feature-value pairs, the tree-based methods (methods I and II) performed in under a minute, and the rule list method (method III) performed in under seven minutes. The most complicated data set in terms of both the number of samples and feature-value pairs is the U.S. Census data set (∼2.09 million training samples, ∼400 feature-value pairs), which took around 1.75 hours for method I, 2.5 hours for method II, and 8 and $\frac{2}{3}$ hours for method III (note that majority of the time for method III went into data loading and preprocessing). In our implementation, we represent data through bit vectors; thus, an increase in the number of samples causes a relatively smaller increase in the run time compared with the run-time increase caused by a larger number of feature categories. We also found that the leaf-sparse density tree method (method I) performed fastest on average for all data sets considered. In Figure 9, we provide a visualization of the time taken to estimate the density of each data set given its complexity for all three methods. Table C.2 in Appendix C shows more details on the timing. All time measurements are averaged over five repeats.

6.2. Recovering a Sparse Tree That Generates Data

This experiment is a test to see whether we can recover a tree that actually generates a data set. Specifically, we generated a data set that arises from a tree with six leaves, involving three features. The data consist of 1,000 data points, where 100 points are tied at value (1, 2, 1), 100 points are at (1, 2, 2), 100 points are at (2, 1, 1), 400 points are at (2, 1, 2), and 300 points are at (2, 2, 2).

We trained the methods on half of the data set and tested on the other half. Figure 10(c) shows the scatterplot of out-of-sample performance and sparsity. This is a case where the DET failed badly to recover the true model. It produced a model that was too sparse, with only three leaves. The leaf-sparse density tree method recovered the full tree. Figure 10(a) is the density tree that generated the data set, and we present the tree that we obtained in Figure 10(b), where it is recovered automatically.

7. Consistency

A consistent method has estimates that converge to the real densities as the size of the training set grows. Consistency of conventional histograms is well studied (for example, Abou-Jaoude 1976, Devroye and Györfi 1983). More generally, consistency for general rectangular partitions has been studied by Zhao et al. (1991) and Lugosi and Nobel (1996).

Typical consistency proofs (e.g., Devroye et al. 1996, Ram and Gray 2011) require the leaf diameters to become asymptotically smaller as the size of the data grows. In our case, if the ground truth density is a tree, we do not want our methods to asymptotically produce smaller and smaller bin sizes; we would rather they reproduce the ground truth tree. This means we require a new type of consistency proof.

Definition 1.

Trees have a single root, and there are conditions on each branch. A density value, f_l is associated with each leaf l of the tree.

Definition 2.

Two trees, T₁ and T₂, are equivalent with respect to density f if they assign the same density values to every data point on the domain: that is, $f_{T_1}(x)=f_{T_2}(x)$ for all x. We denote the class of trees that are equivalent to T as ${[T]}_f$.

Theorem 1.

Let Θ be the set of all trees. Consider these conditions.

1. $T_n\in {\text{arg}\hspace{0.17em}\text{max}}_{T}Obj(T)$. The objective function can be decomposed into $\text{Obj}(T)=\ln q_n(T|X)+\ln g_n(T|X)$, where ${\text{arg}\hspace{0.17em}\text{max}}_T[\ln q_n(T|X)+\ln g_n(T|X)]\equiv {\text{arg}\hspace{0.17em}\text{max}}_T\ln g_n(T|X)$ as $n\to \infty$.

2. $\ln g_n(T|X)$ converges in probability, for any tree T, to the empirical log likelihood that is obtained by the maximum likelihood principle, ${\widehat{l}}_n(T|X)=\frac{1}{n}{\sum }_{i=1}^n\ln {\widehat{f}}_n(x_i|T)$.

3. ${\sup }_{T\in \mathrm{\Theta }}|{\widehat{l}}_n(T|X)-l(T)|\stackrel{P}{\to }0$,

   where $l(T)={\mathbb{E}}_x(\ln (f(x|T)))$.

4. $T_{\text{MLE}}^{*}\in {\text{arg}\hspace{0.17em}\text{max}}_{T}l(T)$ is unique up to equivalence among elements of ${[T_{\text{MLE}}^{*}]}_f$.

If these conditions hold, then the trees T_n that we learned, $T_n\in {\text{arg}\hspace{0.17em}\text{max}}_T\mathit{\text{Obj}}(T)$, obey $T_n\in {[T_{\text{MLE}}^{*}]}_f$ for n > M for some M.

The proof of the result is as follows.

Proof of Theorem 1.

The first condition and the second condition are true any time we use a Bayesian method. They are also true any time we use regularized empirical likelihood, where the regularization term’s effect fades with the number of observations. Note that the third condition is automatically true by the law of large numbers. The last condition is not automatically true and requires regularity conditions for identifiability. The result states that our learned trees are equivalent to maximum likelihood trees when there are enough data.

From definition of T_n, T_n is an optimal value of the log-objective function, and hence, it is also an optimal solution to g_n as n is sufficiently large because of the first condition. We have that

$$\ln \text{Obj}(T_n|X)-\ln \text{Obj}(T_{\text{MLE}}^{*}|X)\ge 0$$

by definition of T_n as the maximizer of Obj. Because Obj becomes close to g_n, we have that

$$\ln g_n(T_n|X)-\ln g_n(T_{\text{MLE}}^{*})\ge 0$$(5)

as n is sufficiently large.

From condition (2), we know that $\ln g_n(T|X)-{\widehat{l}}_n(T|X)\stackrel{P}{\to }0$, and from condition (3), we have ${\widehat{l}}_n(T|X)-l(T)\stackrel{P}{\to }0.$ Adding this up using the fact that convergence in probability is preserved under addition, we know that $\ln g_n(T|X)\stackrel{P}{\to }l(T)$.

Hence, by taking the limit of (5) as n grows, we have that

$$\underset{n\to \infty }{\lim } l(T_n|X)\ge l(T_{\text{MLE}}^{*}).$$

Because $T_{\text{MLE}}^{*}$ is optimal for l(T) by definition and by condition (4), we conclude that T_n stays in $[T_{\text{MLE}}^{*}]$ when n is sufficiently large. □

8. Conclusion

In this work, we have presented a Bayesian approach to density estimation using sparse piecewise constant estimators for categorical (or binary) data. Our methods have nice properties, including that their prior encourages sparsity, which permits interpretability.

For tree-based methods, the prior is the user’s desired number of leaves or branches in each node in the tree, whereas for the density rule list, the prior regularizes the length of the list. We designed a simulated annealing scheme, which alongside the inclusion-exclusion principle and efficient data representation via bit vectors, allows us to find an optimal solution relatively fast.

Our methods outperform existing baselines. They do not have the pitfalls of other nonparametric density estimation methods like density estimation trees, which are top-down greedy. Further, they are consistent, without needing to asymptotically produce infinitesimally small leaves.

The interpretability of density trees and rule lists allows for easier visualization of the estimated density values, aiding in the understanding of the data distribution, detection of outliers and errors, and model selection, and it could assist with decision making. The approaches presented here have given us insight into real data sets (including the house break-in data set from the Cambridge Police and COCO-Stuff labels) that we could not have obtained reliably in any other way.