Priors for Random Count Matrices Derived from a Family of Negative Binomial Processes

Introduction

The need to model a random count matrix arises in many settings, from linguistics to marketing to ecology. For example, in text analysis, we often observe a document-term matrix, whose rows record how many times word $k$ appeared in a given document. In a biodiversity study, we may observe a site-species matrix, where each row records the number of times species $k$ was observed at a given site. Similar applications arise in a wide variety of fields; for examples, see Cameron and Trivedi (1998), Chib et al. (1998), Canny (2004), Buntine and Jakulin (2006), Winkelmann (2008), Titsias (2008), and Zhou et al. (2012).

Nonparametric Bayesian analysis provides a natural setting in which to study random matrices, especially those with no natural upper bound on the number of rows or columns. Yet while there is a wide selection of nonparametric Bayesian models for random count vectors and random binary matrices, prior distributions over random count matrices are relatively underdeveloped. Moreover, a major conceptual problem in modeling a random count matrix arises when new rows are added sequentially. For example, as new documents are collected and processed in text analysis, each new document (represented by a new row of the matrix) may contain previously unseen words (features). This requires that new columns be added to the existing count matrix. But it is not obvious how to define the predictive distribution of this new row of a random count matrix, if the row contains previously unseen features. This is especially important in natural language processing, where a common application is to build a naive Bayes model for classifying new documents. Without having a predictive distribution that accounts for new features, one must often use a predetermined vocabulary and simply ignore the previously unseen terms appearing in a new document.

We directly address these issues by investigating a family of nonparametric Bayesian priors for random count matrices constructed from stochastic processes: the gamma-Poisson process, the gamma-negative binomial process (GNBP), and the beta-negative binomial process (BNBP). We show that all these processes lead to random count matrices with independent and identically distributed (i.i.d.) columns, which can be constructed by drawing all the columns at once, or by adding one row at a time. In addition, we show the gamma-Poisson process, and for special cases of the GNBP and BNBP with common row-wise parameters, the generated random count matrices are exchangeable in both rows and columns.

Our derivation exactly marginalizes out the underlying stochastic processes to arrive at a probability mass function (PMF) for a column-i.i.d. random count matrix. In contrast to existing techniques that take the infinite limit of a finite-dimensional model, this novel procedure allows for the construction and analysis of much more flexible nonparametric priors for random matrices, and highlights certain model properties that are not evident from the finite-model limit. The argument relies upon a novel combinatorial analysis for calculating the number of ways to map a column-i.i.d. random count matrix to a structured random count matrix whose columns are ordered in a certain manner. This is a key step in deriving the predictive distribution of a new random count vector under a random count matrix.

As an application of our proposed framework, we construct a naive-Bayes text classification model. The approach does not require a predefined list of terms (features), and naturally accounts for documents with previously unseen terms. This also implies that random count matrices of different categories can be updated, analyzed, and tested completely in parallel. Moreover, the algorithm requires neither feature selection nor parameter tuning. Following Crammer et al. (2012), the algorithm may also be conveniently extended to an online learning setting. Empirical results suggest that both the proposed GNBP and BNBP models lead to substantially better out-of-sample classification performance, versus both the gamma-Poisson model and the multinomial model with Laplace smoothing. They also clearly outperform the text classification algorithms that first learn lower-dimensional feature vectors for documents and then train a multi-class classifier, and have comparable performance to the state-of-the-art discriminatively trained text classification algorithms, whose features need to be carefully constructed and parameters carefully selected.

2 Connections with existing work

Our paper is in the spirit of existing work on nonparametric Bayesian priors for random count vectors and random binary matrices. To model a random count vector, one may use the Chinese restaurant process, or any one of many other stochastic processes characterized by exchangeable partition probability functions (EPPFs) or sample-size dependent EPPFs; see, for example, Blackwell and MacQueen (1973), Pitman (2006), Lijoi and Prünster (2010), and Zhou and Walker (2014). Likewise, to model a random binary matrix, one may use the Indian buffet process (Griffiths and Ghahramani, 2005, Teh and Gorur, 2009). These well-studied nonparametric Bayesian priors, however, are not directly useful for describing random count matrices. To address this gap, we investigate a family of nonparametric Bayesian priors for random count matrices, each based on a previously proposed stochastic process that has not been thoroughly studied: the gamma-Poisson process (Lo, 1982, Titsias, 2008), the gamma-negative binomial process, or GNBP (Zhou and Carin, 2015); and the beta-negative binomial process, or BNBP (Zhou et al., 2012, Broderick et al., 2015).

All three models can be derived as the marginal distribution of a suitably defined stochastic process with respect to a traditional sampling model for integer-valued counts. This parallels the construction of the models for count vectors or binary matrices mentioned previously. For example, the Chinese restaurant process describes a random count vector as the marginal of the Dirichlet process (Ferguson, 1973) under multinomial sampling. Likewise, the Indian buffet process describes a random binary matrix as the marginal of the beta process (Hjort, 1990) under Bernoulli sampling (Thibaux and Jordan, 2007). Similarly, we present the negative binomial process as the marginal of the gamma process under Poisson sampling, the GNBP as the marginal of the gamma process under negative binomial sampling, and the BNBP as the marginal of the beta process under negative binomial sampling.

The remainder of the paper is organized as follows. After some preliminary definitions and notation, we introduce in Section 2 three distinct nonparametric Bayesian priors for random count matrices. In Section 3, we construct nonparametric Bayesian naive Bayes classifiers to classifier a count vector to one of several existing count matrices and demonstrate their use in document categorization. The details for deriving the random count matrix priors from their underlying hierarchical stochastic processes are provided in the Supplementary Material.

3 Notation and preliminaries

A beta process (Hjort, 1990) $B\sim\mbox{BP}(c,B_{0})$ on the product space $\times\Omega$, is also defined by two parameters: a finite and continuous base measure $B_{0}$ over a complete separable metric space $\Omega$, and a concentration parameter $c>0$. The Lévy measure of the beta process in this paper is defined as

As $\int_{\times\Omega}\nu(dpd\omega)=\infty$ and $\int_{\times\Omega}\min\{p,1\}\nu(dpd\omega)<\infty$, a draw from $B\sim\mbox{BP}(c,B_{0})$ can be represented as $B=\sum_{k=1}^{\infty}p_{k}\delta_{\omega_{k}},~{}\omega_{k}\sim g_{0},$ where $\gamma_{0}=B_{0}(\Omega)$ is the mass parameter and $g_{0}(d\omega)=B_{0}(d\omega)/\gamma_{0}$ is the base distribution.

Our convention is that a prior for a random count matrix is named by the stochastic process used to generate each of its rows. In this paper, we study three hierarchical stochastic processes, all in the family of negative binomial processes. Each such stochastic process is defined by the prior for an almost-surely discrete random measure, together with a sampling model for generating counts. We denote the distribution of such a matrix as ${{\bf N}}\sim\mbox{ProcessM}(\boldsymbol{\theta})$, where “Process” is the name of the underlying hierarchical stochastic process, “M” stands for matrix, and $\boldsymbol{\theta}$ encodes the parameters of the process.

For example, to construct a gamma-Poisson or negative binomial process random count matrix, ${{\bf N}}_{J}\sim\mbox{NBPM}(\gamma_{0},c)$, we draw a random measure $G\sim\Gamma\mbox{P}(G_{0},1/c)$ from a gamma process. Then for each row of the matrix, we independently draw $X_{j}\mid G\sim\mbox{PP}(G)$: a Poisson process such that $X_{j}(A)\sim\mbox{Pois}[G(A)]$ for all $A\subset\Omega$. As $G=\sum_{k=1}^{\infty}r_{k}\delta_{\omega_{k}}$ is atomic, we have $X_{j}=\sum_{k=1}^{\infty}n_{jk}\delta_{\omega_{k}},~{}n_{jk}\sim\mbox{Pois}(r_{k})$. Although $\{X_{j}\}_{1,J}$ contains countably many atoms, we will show in later sections that only a finite number of them have nonzero counts. The count matrix ${{\bf N}}_{J}$ is constructed by organizing all the nonzero column count vectors, $\{{\boldsymbol{n}}_{:k}\}_{k:n_{\boldsymbol{\cdot}k}>0}$, in an arbitrary order into a random count matrix. Thus the statistical features we care about, such as words or species, are identified with the atoms of the underlying random measure.

The notation $u\sim\mbox{Log}(p)$ denotes a random variable having a logarithmic distribution (Quenouille, 1949) with PMF

A related distribution, called the sum-logarithmic, is defined as follows. Let $u_{t}\sim\mbox{Log}(p)$, and let $n=\sum_{t=1}^{l}u_{t}$. The marginal distribution of $n$ is a sum-logarithmic distribution (Zhou and Carin, 2015), expressed as $n\sim\mbox{SumLog}(l,p)$, with PMF

where $|s(n,l)|$ are unsigned Stirling numbers of the first kind. These are related to gamma functions by

The joint distribution of $n\sim\mbox{SumLog}(l,p)$ and $l\sim\mbox{Pois}[-r\ln(1-p)]$ is described as the Poisson-logarithmic bivariate distribution in Zhou and Carin (2015), with PMF

The marginalization of $l$ from this compound Poisson representation leads to the negative binomial distribution $n\sim\mbox{NB}(r,p)$, with PMF

We describe in the Supplementary Material several other useful distributions, including the logarithmic mixed sum-logarithmic (LogLog), the negative binomial mixed sum-logarithmic, the gamma-negative binomial (GNB), the beta-negative binomial (BNB), the digamma distribution, and the logbeta distributions.

Nonparametric Priors for Random Count Matrices

Let ${{\bf N}}_{J}\sim\mbox{NBPM}(\gamma_{0},c)$ denote a gamma-Poisson or negative binomial process (NBP) random count matrix, parameterized by a mass parameter $\gamma_{0}$ and a concentration parameter $c$. This prior arises from marginalizing out the gamma process $G\sim\Gamma\mbox{P}(G_{0},1/c)$ from $J$ conditionally independent Poisson process draws $X_{j}\mid G\sim\mbox{PP}(G)$, with the rows of ${{\bf N}}_{J}$ corresponding to the $X_{j}$’s and the columns of ${{\bf N}}_{J}$ corresponding to the atoms with at least one nonzero count.

As $\{X_{j}\}_{1,J}$ are i.i.d. given $G$, they are exchangeable according to de Fennetti’s theorem. With a draw from the gamma-Poisson process expressed as $X_{j}=\sum_{k=1}^{\infty}n_{jk}\delta_{\omega_{k}},~{}n_{jk}\sim\mbox{Pois}(r_{k})$, where $r_{k}=G(\omega_{k})$ is the weight of the atom $\omega_{k}$ of the gamma process $G\sim\Gamma\mbox{P}(G_{0},1/c)$, we may write the likelihood of $\{X_{j}\}_{1,J}$, given $G$, as

Fix an arbitrary labeling of the indices of the atoms in $\mathcal{D}_{J}$ from $1$ to $K_{J}$. We now appeal to the definition of a gamma process and rewrite the conditional likelihood of $\{X_{j}\}_{1,J}$ as

where $G(\Omega\backslash\mathcal{D}_{J}):=\sum_{k:n_{k}=0}r_{k}$ is the total mass of the rest of the (absolutely continuous) space. The idea is to first marginalize out $G$ from (4) to obtain the marginal distribution $p(\{X_{j}\}_{1,J}\mid\gamma_{0},c)$, whose derivation using the Palm formula is provided in the Supplementary Material, and then use combinatorial argument to find the marginal distribution of the random count matrix ${{\bf N}}_{J}$ organized from $\{X_{j}\}_{1,J}$.

1.2 Marginal distribution and combinatorial analysis

One of our main results is that the PMF of ${{\bf N}}_{J}\sim\mbox{NBPM}(\gamma_{0},c)$, with $J$ rows and a random $K_{J}$ number of columns, is

where the unordered column vectors $\{{\boldsymbol{n}}_{:k}\}_{1,K_{J}}$ of the count matrix ${{\bf N}}_{J}$ represent a draw from the underlying stochastic process, and the normalization constant of $1/K_{J}!$ arises from the fact that the mapping from a realization of $\{X_{j}\}_{1,J}$ to ${{\bf N}}_{J}$ is one-to-many, with $K_{J}!$ distinct column orderings.

By construction, the rows of a NBP random count matrix are exchangeable. Moreover, one may verify by direct calculation that a NBP random count matrix with PMF (5) can be generated column by column as i.i.d. count vectors:

It is clear from (2.1.2) that the columns of ${{\bf N}}_{J}$ are independent multivariate count vectors, which all follow the same logarithmic-multinomial (mixture) distribution. Thus the NBP random count matrix ${{\bf N}}_{J}$ is row-column exchangeable (see, e.g. Hoover, 1982, Aldous, 1985, Orbanz and Roy, 2014, for a general treatment of row-column exchangeable matrices).

where $\boldsymbol{\theta}:=\{\gamma_{0},c\}$ and $p({{\bf N}}^{+}_{j+1}\!\mid\!{{\bf N}}_{j},\boldsymbol{\theta}):=f({{\bf N}}_{j+1}\!\mid\!\boldsymbol{\theta})/f({{\bf N}}_{j}\!\mid\!\boldsymbol{\theta})$ is the prediction rule to add the new part brought by row $(j+1)$ into the matrix ${{\bf N}}_{j}$. Direct calculations using (2.1.2) yield the following form for this prediction rule, expressed in terms of familiar PMFs:

The normalizing constant $(K_{J}!\ K^{+}_{J+1}!)/K_{J+1}!$ in (2.1.2) plays a key role in our combinatorial analysis, and will appear again in both the gamma- and beta- negative binomial processes. It emerges directly from the calculations, and can also be interpreted in the following way. After drawing $K_{J+1}^{+}$ new columns, we must insert them into the original $K_{J}$ columns while keeping the relative orders of both the original and new columns unchanged. This is a one-to-many mapping, with the number of such order-preserving insertions given by the binomial coefficient. For example, if the original ${{\bf N}}_{J}$ has two columns and the new row $J+1$ introduces two more columns, then we construct ${{\bf N}}_{J+1}$ by rearranging the two old columns 1 and 2 and the two new columns iii and iv in one of $\binom{4}{2}=6$ possible ways: (1 2 iii iv), (1 iii 2 iv), (iii 1 2 iv), (1 iii iv 2), (iii 1 iv 2), and (iii iv 1 2), where (1 2 iii iv) represents the construction appending the new columns to the right of the original matrix.

It is instructive to compare (2.1.2), which generates a NBP random matrix by drawing all its columns at once, with (2.1.2), which generates an identically distributed random matrix one row at a time. The matrix generated with (2.1.2) has i.i.d. columns. The matrix generated with (2.1.2) adds $K^{+}_{J+1}$ new columns when it adds the $(J+1)$th row, and if the newly added columns are inserted into random locations among original columns with their relative orders preserved, then we arrive at an identically distributed column-i.i.d. random count matrix. If the newly added columns are inserted in a particular way, then the distribution of the generated random matrix would be different up to a multinomial coefficient. For example, if we generate row vectors ${\boldsymbol{n}}_{j}$ from $j=1$ to $j=J$ and each time we append the new columns to the right of the original matrix, then this ordered matrix $\widetilde{{{\bf N}}}_{J}$ will appear with probability

Shown in the first row of Figure 1 are three NBP random count matrices simulated in this manner. We note that the gamma-Poisson process is related to the model of Lo (1982), as well as the model of Titsias (2008), which can be considered as a special case of the NBP with the concentration parameter $c$ fixed at one.

1.3 Inference for parameters

Although the marginal likelihood alone is not amenable to posterior analysis, the NBP parameters can be conveniently inferred using both the conditional and marginal likelihoods. To complete the model, we let $\gamma_{0}\sim\mbox{Gamma}(e_{0},{1}/{f_{0}})$ and $c\sim\mbox{Gamma}(c_{0},1/d_{0})$. With (4), (5) and $G(\Omega):=G(\Omega\backslash\mathcal{D}_{J})+\sum_{k=1}^{K_{J}}r_{k}$, we sample the parameters in closed form as

Similar strategies will be used to infer the parameters of the other two stochastic processes. Having closed-form update equations for parameter inference via Gibbs sampling is a unique feature shared by all the nonparametric Bayesian priors proposed in this paper.

2 The gamma-negative binomial process

Let ${{\bf N}}_{J}\sim\mbox{GNBPM}(\gamma_{0},c,p_{1},\ldots,p_{J})$ denote a gamma-negative binomial process (GNBP) random count matrix, parameterized by a mass parameter $\gamma_{0}$, a concentration parameter $c$, and $J$ row-specific probability parameters $\{p_{j}\}_{1,J}$. This random count matrix is the direct outcome of marginalizing out the gamma process $G\sim\Gamma\mbox{P}(G_{0},1/c)$, with data augmentation, from $J$ conditionally independent negative binomial process draws $X_{j}\mid G\sim\mbox{NBP}(G,p_{j})$, which are defined such that $X_{j}(A)\sim\mbox{NB}\left(G(A),p_{j}\right)$ for each $A\subset\Omega$.

As directly marginalizing out the gamma process under negative binomial sampling is difficult, our construction is based on the compound-Poisson representation of the negative binomial, described in Section 1.3. Specifically, consider the joint distribution of ${{\bf N}}_{J}$ and a latent count matrix ${{\bf L}}_{J}$, whose dimension and locations of nonzero counts are the same as those of ${{\bf N}}_{J}$. These two matrices parallel the scalar $n$ and $l$ given in the joint PMF of the Poisson-logarithmic distribution (3). This joint distribution is defined as

where $\boldsymbol{\theta}:=\{\gamma_{0},c,p_{1},\ldots,p_{J}\}$, $q_{j}:=-\ln(1-p_{j})$ and $q_{\boldsymbol{\cdot}}:=\sum_{j=1}^{J}q_{j}$. The detailed derivation is in the Supplementary Material.

Similar to the analysis in Section 2.1 for the NBP, we show in the Supplementary Material that the GNBP random count matrix can be constructed by either drawing its i.i.d. columns at once or adding one row at a time, and it has closed-form Gibbs sampling update equations for model parameters. Different from the NBP random count matrix that is row-column exchangeable, the GNBP random count matrix no longer maintains row exchangeability if its row-wise probability parameters $p_{j}$ are set differently for different rows.

Shown in the second row of Figure 1 are three sequentially constructed GNBP random count matrices, with the new columns introduced by each row appended to the right of the matrix. Similar to the combinatorial arguments that lead to (8), this particularly structured matrix and its auxiliary matrix appear with probability $\binom{K_{J}}{K^{+}_{1},\ldots,K^{+}_{J}}f({{\bf N}}_{J},{{\bf L}}_{J}\!\mid\!\boldsymbol{\theta})$.

3 The beta-negative binomial process

Let ${{\bf N}}_{J}\sim\mbox{BNBPM}(\gamma_{0},c,r_{1},\ldots,r_{J})$ denote a beta-negative binomial process (BNBP) random count matrix, parameterized by a mass parameter $\gamma_{0}$, a concentration parameter $c$, and $J$ row-specific dispersion parameters $\{r_{j}\}_{1,J}$, whose PMF is defined as

where $\boldsymbol{\theta}:=\{\gamma_{0},c,r_{1},\ldots,r_{J}\}$. The PMF is the direct outcome of marginalizing out the beta process $B\sim\mbox{BP}(c,B_{0})$ from $J$ conditionally independent negative binomial process draws $X_{j}\!\mid\!B\sim\mbox{NBP}(r_{j},B)$, which are defined such that $X_{j}(A)=\sum_{k:\omega_{k}\in A}n_{jk},~{}n_{jk}\sim\mbox{NB}(r_{j},p_{k})$ for each $A\subset\Omega$, where $p_{k}=B(\omega_{k})$ is the weight of atom $k$. The detailed derivation is provided in the Supplementary Material.

Similar to the analysis in Section 2.1 for the NBP, we show in the Supplementary Material that the BNBP random count matrix can be constructed by either drawing its i.i.d. columns at once or adding one row at a time using an “ice cream” buffet process, and it has closed-form Gibbs sampling update equations for all model parameters except for the concentration parameter $c$. The BNBP random count matrix no longer maintains row exchangeability if its row-wise dispersion parameters $r_{j}$ are set differently for different rows.

Shown in the last row of Figure 1 are three sequentially constructed BNBP random count matrices, with the new columns introduced by each row appended to the right of the matrix. Similar to the combinatorial arguments that lead to (8), this particularly structured matrix appears with probability $\binom{K_{J}}{K^{+}_{1},\ldots,K^{+}_{J}}f({{\bf N}}_{J}\!\mid\!\boldsymbol{\theta})$.

4 The predictive distribution of a new row count vector

It is critical to note that the prediction rule $p({{\bf N}}^{+}_{J+1}\!\mid\!{{\bf N}}_{J},\boldsymbol{\theta})$ of the NBP shown in (2.1.2) is for sequentially constructing a column-i.i.d. random count matrix, but it is not the predictive distribution for a new row count vector. The $1\times K_{J}$ submatrix of ${{\bf N}}^{+}_{J+1}$ orders its column in the same way as ${{\bf N}}_{J}$ does, and the $(J+1)\times K^{+}_{J+1}$ submatrix of ${{\bf N}}^{+}_{J+1}$ also maintains a certain order of its columns; however, the indexing of these $K^{+}_{J+1}$ columns are in fact arbitrarily chosen from $K^{+}_{J+1}!$ possible permutations. Therefore, the predictive distribution of a row vector ${\boldsymbol{n}}_{J+1}$ that brings $K^{+}_{J+1}$ new columns shall be

The normalizing constant $1/{K^{+}_{J+1}!}$ in (12) arises because a realization of ${{\bf N}}^{+}_{J+1}$ to ${\boldsymbol{n}}_{J+1}$ is one-to-many, with $K^{+}_{J+1}!$ distinct orderings of these new columns brought by the $(J+1)$th row. Our experimental results show that omitting this normalizing term may significantly deteriorate the out-of-sample prediction performance.

An equivalent representation in (13) shows that one may first consider the distribution of a matrix constructed by appending the new columns brought by ${\boldsymbol{n}}_{J+1}$ to the right of ${{\bf N}}_{J}$, which is $\frac{K_{J+1}!}{K_{J}!K^{+}_{J+1}!}f({{\bf N}}_{J+1}\!\mid\!\boldsymbol{\theta})$, and then apply the Bayes’ rule to derive the conditional distribution of this particularly ordered ${\boldsymbol{n}}_{J+1}$ given ${{\bf N}}_{J}$. The normalizing constant ${K_{J}!}/{K_{J+1}!}$ in (13) can be interpreted in the following way. We need to insert the $K^{+}_{J+1}$ new columns one by one into the original matrix. The first, second, $\ldots$, and last new columns can choose from $K_{J}+1$, $K_{J}+2$, $\ldots$, and $K_{J}+K^{+}_{J+1}$ possible locations, respectively, thus there are $\prod_{i=1}^{K^{+}_{J+1}}(K_{J}+i)!={K_{J+1}!}/{K_{J}!}$ ways to insert the $K^{+}_{J+1}$ new columns into the original ordered $K_{J}$ columns, which is again a one-to-many mapping. The same combinatorial analysis applies to both the GNBP and BNBP. For the GNBP, to compute the predictive likelihood of ${\boldsymbol{n}}_{J+1}$, one will need to take extra care as the computation involves ${{\bf L}}_{J}$, an auxiliary random count matrix that is not directly observable. In Section 3, we will discuss in detail how to compute the predictive likelihood via Monte Carlo integration.

5 Comparison

In the Supplementary Material, we provide further details on the construction of random count matrices from the negative binomial process, as well as those derived from both the gamma-negative binomial process (GNBP) and beta-negative binomial process (BNBP). While the PMFs for all three proposed nonparametric priors are complicated, their relationship and differences become evident once we show that they all govern random count matrices with a Poisson-distributed number of i.i.d. columns. Table 1 shows the differences among the three priors’ row-wise sequential construction, and the following list shows the variance-mean relationship for each prior for the counts at existing columns. Together, these provide additional insights on how the priors differ from each other.

The NBP can be used to generate a row-column exchangeable random count matrix with a potentially unbounded number of columns. However, as shown in (2.1.2), to model the total count of a column $n_{\boldsymbol{\cdot}k}$, the NBP uses the logarithmic distribution, which has only one free parameter, always has the mode at one, and monotonically decreases. In addition, each column sum $n_{\boldsymbol{\cdot}k}$ is assigned to the $J$ rows with a multinomial distribution that has a uniform probability vector $(1/J,\ldots,1/J)$. Furthermore, as shown in Table 1, for out-of-sample prediction, it models counts at existing columns using $\mbox{NB}\left[n_{(J+1)k};n_{\boldsymbol{\cdot}k},{1}/{(J+c+1)}\right]$, whose variance-mean relationship (14) may be restrictive in modeling highly overdispersed counts. Finally, the expected number of new columns brought by a row, equal to $\gamma_{0}\ln[1+{1}/{(J+c)}]$, monotonically decreases. These constraints limit the potential use of the NBP model.

The variance-mean relationships expressed by (14)-(16) show that the GNBP and BNBP can model much more overdispersed counts than the NBP. This fact is borne out by the simulated random count matrices in Figure 1, which provide some intuition for the practical differences among the models. The parameters for the three priors have been chosen so that each random matrix has the same expected total count. Yet the counts in the NBP random count matrices have small dynamic ranges, whereas the counts in both the GNBP and BNBP matrices can contain values that are significantly above the average.

6 Parameter inference

An appealing feature of all three negative binomial process random count matrix priors is that their parameters can be inferred with closed-form Gibbs sampling update equations, by exploiting both the conditional and marginal distributions, together with the data augmentation and marginalization techniques unique to the negative binomial distribution. Parameter inference for the NBP is provided in Section 2.1.3. The details of parameter inference for both the GNBP and BNBP are provided in the Supplementary Material.

Negative Binomial Process Naive Bayes Classifiers

Given a random count matrix, finding the predictive distribution of a row count vector, which may bring additional columns, involves interesting and challenging combinatory arguments that have been throughly addressed in this paper. With these combinatorial structures carefully analyzed, we are ready to construct a NBP, a GNBP, and a BNBP naive Bayes classifiers. We do so as follows. First, for each category, the training row count vectors are summarized as a random count matrix ${{\bf N}}_{J}$, each column of which must contain at least one nonzero count (i.e. columns with all zeros are excluded). Second, Gibbs sampling is used to infer the parameters $\boldsymbol{\theta}$ that generate ${{\bf N}}_{J}$. To represent the posterior of $\boldsymbol{\theta}$, $S$ MCMC samples $\{\boldsymbol{\theta}^{[s]}\}_{1,S}$ are collected. For the GNBP, a posterior MCMC sample ${{\bf L}}_{J}^{[s]}$ for the auxiliary random matrix is also collected when $\boldsymbol{\theta}^{[s]}$ is collected. Finally, to test a row count vector ${\boldsymbol{n}}_{J+1}$, its predictive likelihood given ${{\bf N}}_{J}$ is calculated via Monte Carlo integration using

for the GNBP. Although a larger $S$ shall lead to a more accurate calculation of the predictive likelihood, the computational complexity for testing is a linear function of $S$. It is therefore of practical importance to find out how the value of $S$ impacts the performance of the proposed nonparametric Bayesian naive classifiers. Below we consider experiments on document categorization, for which we will show that $S=1$ performs essentially just as well as selecting a much larger $S$ in terms of the categorization accuracy.

2 Experiment settings

We consider the example of categorizing the 18,774 documents of the 20 newsgroups datasethttp://qwone.com/$\sim$jason/20Newsgroups/, where each bag-of-words document is represented as a word count vector under a vocabulary of size $V=$ 61,188. We also consider the TDT2 corpushttp://www.cad.zju.edu.cn/home/dengcai/Data/TextData.html ( NIST Topic Detection and Tracking corpus): with the documents appearing in two or more categories removed, this subset of TDT2 consists of 9,394 documents from the largest 30 categories, with a vocabulary of size $V=$ 36,771; this dataset was used to compare document clustering algorithms in Cai et al. (2005). We train all three negative binomial processes using 10%, 20%, $\ldots$, or 80% of the documents in each newsgroup of the 20 newsgroups dataset, and in each category of the TDT2 corpus. We then test on the remaining documents. We report our results based on five random training/testing partitions.

To make comparison to other commonly used text categorization algorithms, we also consider a default setting for the 20 newsgroups dataset: using the first 11,269 documents for training and the other 7,505 documents collected at later times for testing. For this setting, we reports our results based on five independent runs with random initializations. This allows us to compare our performance to many other papers that have proposed text classification algorithms and benchmarked their methods using this same split of the 20 newsgroups dataset.

We collect $S$ MCMC samples of model parameters and auxiliary variables to compute the predictive likelihood for a new row count vector. In this paper, we run $S$ independent Markov chains and collect the 2500th sample of each chain. Note that one may also consider collecting $S$ samples at a certain interval from a single Markov chain after the burn-in period. We consider non-informative hyper-parameters as $a_{0}=b_{0}=\ldots=f_{0}=0.001$. For the BNBP, we set $c_{0}=d_{0}=1$. The document-term training count matrix of the $i$th newsgroup is modeled as ${{\bf N}}^{(i)}_{J^{(i)}}\sim\mbox{NBPM}(\gamma^{(i)}_{0},c^{(i)})$, {{\bf N}}^{(i)}_{J^{(i)}}\sim\mbox{GNBPM}\big{(}\gamma^{(i)}_{0},c^{(i)},p^{(i)}_{1},\ldots,p^{(i)}_{J^{(i)}}\big{)}, and {{\bf N}}^{(i)}_{J^{(i)}}\sim\mbox{BNBPM}\big{(}\gamma^{(i)}_{0},c^{(i)},r^{(i)}_{1},\ldots,r^{(i)}_{J^{(i)}}\big{)} under the three priors respectively.

Note that we are facing typical “small $n$ and large $p$” problems as the number of rows of a document-term count matrix is typically much smaller than the number of columns. For example, the first newsgroup of the 20 newsgroups dataset contains 798 documents with 12,665 unique words, which is summarized as a $798\times 12665$ count matrix; and the 30th category of the TDT2 subset contains 52 documents with 2904 unique words, which is summarized as a $52\times 2904$ count matrix. As the number of unique terms in a category might be significantly smaller than the vocabulary size of the whole corpus, our approach for both training and testing could be much faster than the approach that considers all the terms in the vocabulary of the corpus. In addition, our approach provides a principled, model-based way to handle terms that appear in a testing document but not in the training documents. By contrast, many traditional approaches have to discard these terms not present in training.

3 Training and posterior predictive checking

It is clear that the NBP is restrictive, in that the generated random matrix looks the least similar to the observed count matrix. This is unsurprising, as the NBP has a limited ability to model highly overdispersed counts, does not model row-heterogeneity, and can barely adjust the number of new columns brought by a row. On the other hand, both the generated GNBP and BNBP random count matrices resemble the original count matrix much more closely. This is expected, since both priors use heavy-tailed count distributions to model highly overdispersed counts, and have row-wise probability or dispersion parameters to model row-heterogeneity and to control the number of new columns brought by each row. Note that the observed matrix has 2904 columns, but each of the generated random count matrices has a different (random) number of columns. This is because there are one-to-one correspondences between their row indices, but not their column indices.

4 Out-of-sample prediction and categorization for count vectors

For out-of-sample prediction on a new row vector, we first compute that vector’s likelihood under different categories’ training count matrices. We then use these likelihoods in a naive-Bayes classifier to categorize the new vector. For example, for testing row count vector ${\boldsymbol{n}}_{j^{\prime}}$ under category $i$, we will first match the column indices (features) of this row count vector to those of the training count matrix ${{\bf N}}^{(i)}_{J^{(i)}}$; each feature that belongs to one of the $K^{(i)}_{J^{(i)}}$ features of ${{\bf N}}^{(i)}_{J^{(i)}}$ but not present in ${\boldsymbol{n}}_{j^{\prime}}$ will be assigned a zero count; and the ${K^{+}_{j^{\prime}}}^{(i)}$ features that are present in vector $j^{\prime}$ but not in ${{\bf N}}^{(i)}_{J^{(i)}}$ will be treated as new features brought by vector $j^{\prime}$ to to ${{\bf N}}^{(i)}_{J^{(i)}}$. For the the GNBP, we first find an estimate of $p^{(i)}_{j^{\prime}}$ as $p^{(i)}_{j^{\prime}}={(a_{0}+n^{(i)}_{j^{\prime}\boldsymbol{\cdot}})}/{[a_{0}+b_{0}+n^{(i)}_{j^{\prime}\boldsymbol{\cdot}}+G^{(i)}(\Omega)]}$. For the BNBP, we first find an expectation-maximization estimate of $r_{j^{\prime}}$ by running the updates

iteratively for 20 iterations, where for a testing row vector with all zeros, we let $l^{(i)}_{j^{\prime}\boldsymbol{\cdot}}=1$. Given the column sums of ${{\bf N}}^{(i)}$ and the inferred model parameters (together with auxiliary variables for the GBNB), the predictive likelihoods of a new row count vector are calculated using (17) for both the NBP and BNBP and with (18) for the GNBP.

Note that when the predictive distributions are used to calculate the likelihoods, the models are not constrained under a predetermined vocabulary. But if we are given a vocabulary of size $V$ that includes all the important terms, exploiting that information might further improve the performance. Thus to test document $j^{\prime}$, we also consider using

as the likelihood for the GNBP, and using

as the likelihood for the BNBP. Note that for this testing procedure we also compute $p({\boldsymbol{n}}_{j^{\prime}}\mid{{\bf N}}_{J}^{(i)})$ using Monte Carlo integration based on $S$ posterior MCMC samples. In contrast to its truly nonparametric Bayesian counterpart with an infinite vocabulary, this testing procedure is expected to have higher computational complexity, but may produce better out-of-sample prediction if the predetermined finite vocabulary fits the testing documents well. Below we show the results produced by both testing procedures.

For comparison, we consider the multinomial naive Bayes classifier with Laplace smoothing (McCallum and Nigam, 1998, Manning et al., 2008), where a test document $j^{\prime}$ has the likelihood under newsgroup $i$ as

The results of some other commonly used text classification algorithms will also be included as benchmarks. Note that all these classifiers require the same predefined finite vocabulary for both training and testing. Thus any new terms in a testing document that are not listed in that vocabulary must be discarded.

5 Example results

We first consider choosing $S=10$ in (17) and (18) to compute the predictive likelihood $p({\boldsymbol{n}}_{j^{\prime}}\!\mid\!{{\bf N}}_{J}^{(i)})$ for test document $j^{\prime}$. Assuming a uniform prior for all the $C$ categories, we assign document $j^{\prime}$ to category $i$ with probability

and categorize document $j^{\prime}$ to the category under which its word count vector ${\boldsymbol{n}}_{j^{\prime}}$ has the highest probability. As shown in Figures 4 and 4, the NBP has the worst categorization accuracy. Both the BNBP and GNBP clearly outperform the NBP and the multinomial naive-Bayes classifier with Laplace smoothing, especially when the number of training documents is small. Both for fitting the training count matrix and making out-of-sample prediction, the NBP is the most restrictive, as it has only two free parameters $\gamma_{0}$ and $c$. In addition to these two parameters, the GNBP (BNBP) has a probability (dispersion) parameter for each row count vector. Moreover, as both the GNB and BNB distributions are mixed negative-binomial distributions, they have heavier tails that may help model the burstiness of words in documents (Church and Gale, 1995, Madsen et al., 2005, Clinchant and Gaussier, 2008).

For the 20 newsgroups dataset, with the 7,505 documents collected at later times used for testing, our NBP, BNBP, and GNBP with an infinite vocabulary and $S=10$ achieve categorization accuracies of 61.9%, 78.7%, and 80.9%, respectively. With a finite vocabulary they achieve accuracies of 61.7%, 79.1%, and 80.9%, respectively. Despite the simplicity of the model, this performance meets or exceeds that of other competing methods, which we briefly describe. The multinomial naive Bayes classifier with Laplace smoothing achieves an accuracy of 78.1%. Lan et al. (2009) consider a range of reweighted term-frequency features in a $k$-nearest neighbors ($k$NN) classifier. Under an optimal choice of $k$ and set of features, they achievs an accuracy of 69.1%. The same authors report that a support vector machine (SVM) classifier achieves an accuracy of 80.8%. Larochelle et al. (2012) use restricted Boltzmann machine for classification, with an optimized training strategy and cross-validated model parameters. They report an accuracy of 76.2% using binary features for the 5000 most frequent words. The accuracy increases to 79.1% when using binary features for the 25247 most frequent words, but the algorithm is too computationally intensive to include more word features.

We also note that text categorization performance significantly deteriorates if one trains a multi-class classifier on the lower-dimensional features extracted using unsupervised feature learning algorithms, such as latent Dirichlet allocation (LDA) (Blei et al., 2003) or the deep Boltzmann machine (Srivastava et al., 2013). As shown in Srivastava et al. (2013), even with tuned parameters, neither LDA nor deep Boltzmann machines combined with a multinomial logistic regression classifier can achieve an accuracy above 70% on this data set. It is also shown in Zhu et al. (2012) that LDA plus an SVM classifier fails to achieve an accuracy above 65%. The performance of LDA could be improved by using a supervised training strategy (Blei and Mcauliffe, 2008). However, as shown in Zhu et al. (2012), the maximum-entropy discrimination LDA (MedLDA), a state-of-the-art supervised LDA algorithm, still does not achieve an accuracy above 80%, despite the fact that the number of topics and model parameters are carefully tuned through cross validation and complex inference and heavy computations are employed to learn the latent features. Both the BNBP and GNBP naive classifiers, while being tuning-free and fast and simple to train using the raw counts, compare favorably to the state-of-the-art text classification algorithms that often rely on heavy computation and carefully selected features and parameters.

Note that for the proposed naive Bayes classifiers, a larger $S$ usually leads to a more accurate computation of the predictive likelihood via Monte Carlo integration, but may not necessarily lead to a clear gain in accuracy for document categorization. This is confirmed by examining the experimental results with $S$ set as small as one (i.e. a single MCMC sample) on both the 20 newsgroups and TDT2 datasets, which are found to be very similar to the results with $S=10$ that are shown in Figures 4 and 4. This is not surprising since it is not the absolute magnitude of the category-specific predictive likelihoods, only their relative rankings, that determine the categorization accuracy.

To further elaborate on this point, we consider the CNAE-9 datasethttps://archive.ics.uci.edu/ml/datasets/CNAE-9 of Ciarelli and Oliveira (2009), which contains 1080 documents of free text business descriptions of Brazilian companies divided into nine categories, with a vocabulary size of $V=856$; and we randomly select 20% of documents from each category as training, and calculate each test document’s predictive probabilities under the nine categories, using the GNBP naive Bayes classifier with $S=1000$ samples, each of which is the 2500th MCMC sample of an independent Markov chain. As shown in Figure 5 (a), in most cases, there is a little ambiguity on which category a test document should be assigned to. Hence letting $S=1000$ or $S=1$ make little practical difference in terms of categorization accuracy. In Figure 5 (b), from the left to right, we show the boxplot of 1000 accuracies produced by 1000 independent runs of the same testing procedure, each of which is calculated with $S=1$ MCMC sample; the boxplot of 250 accuracies with $S=4$; the boxplot of 100 accuracies with $S=10$; and the boxplot of 20 accuracies with $S=50$. It is clear from Figure 5 (b) that the larger the $S$ is, the less the categorization accuracy varies, which is expected as the error of Monte Carlo integration decreases with $\sqrt{N}$. However, there is no substantial improvement for the mean of the accuracies as $S$ increases. Even with $S=1$, the worst categorization accuracy is not too far from its mean. Therefore, in practice one may simply choose a small $S$ to compute the predictive likelihoods for the purpose of document categorization.

As opposed to the conventional multinomial naive-Bayes classifier that estimates the probability of each word in the vocabulary by normalizing the word counts, the proposed negative binomial processes provide new methods that directly analyze the raw counts and take into account the total length of a document. Moreover, there is no need to predetermine the vocabulary, as new features not present in the training data have been taken care of by the nonparametric Bayesian predictive distributions of the negative binomial processes that are discussed in Section 2.4.

Conclusions

This paper fills a gap in the nonparametric Bayesian literature, deriving a family of probability mass functions for random count matrices by exploiting the gamma-Poisson, gamma-negative binomial, and beta-negative binomial processes. The resulting random count matrices have a random number of i.i.d. columns, and their parameters can be inferred with closed-form update equations. Any random count matrix in this family can be constructed by generating all its i.i.d. columns at once, or by adding one row at a time. Our results also allow us to define the predictive distribution of an infinite-dimensional random count vector under any of the proposed priors, leading to three nonparametric Bayesian naive Bayes classifiers for count vectors. The proposed classifiers, which directly operate on the raw counts and require no parameter tuning, alleviate the need to predetermine a shared finite vocabulary, and can account for features not present in the training data. Example results on document categorization show that the proposed gamma-negative binomial process and beta-negative binomial process clearly outperform both the negative binomial process and the multinomial naive Bayes classifier with Laplace smoothing, and have comparable performance to other state-of-the-art discriminatively-trained text classification algorithms. We are currently extending the techniques developed here to construct nonparametric Bayesian priors for a random count matrix, which has an unbounded number of columns and each row of which sums to a fixed integer; this extension can be used to construct nonparametric Bayesian discrete latent variable models, whose feature usages are represented with infinite random count matrices that are not directly observable.

Acknowledgements

The authors thank the editor, associate editor, and two anonymous referees, whose invaluable comments and suggestions have helped us to improve the paper substantially. JGS acknowledges the support of CAREER grant DMS-1255187 from the U.S. National Science Foundation.

References

Appendix A The Negative Binomial Process: Details

To generate a random count matrix, we construct a gamma-Poisson process as

Zhou and Carin (2015) derives the marginal distribution of $X=\sum_{j=1}^{J}X_{j}$ and calls it as the negative binomial process (NBP), a draw from which is represented as an exchangeable random count vector. We do not consider that simplification in this paper and consequently our definition of the NBP, a draw from which is represented as a row-column exchangeable random count matrix, differs from the one in Zhou and Carin (2015).

The conditional likelihood in (4) can be re-written as

Appendix B Gamma-Negative Binomial Process: Details

Given the gamma process $G\sim\Gamma{\mbox{P}}(G_{0},1/c)$, we define $X\mid G\sim\mbox{NBP}(G,p)$ as a negative binomial process such that $X(A)\sim\mbox{NB}(G(A),p)$ for each $A\subset\Omega$. Replacing the Poisson processes in (A.1) with the negative binomial processes defined in this way yields a gamma-negative binomial process (GNBP):

With a draw from the gamma process $G\sim\Gamma{\mbox{P}}(G_{0},1/c)$ expressed as $G=\sum_{k=1}^{\infty}r_{k}\delta_{\omega_{k}}$, a draw from $X_{j}\mid G\sim\mbox{NBP}(G,p_{j})$ can be expressed as $X_{j}=\sum_{k=1}^{\infty}n_{jk}\delta_{\omega_{k}},~{}n_{jk}\sim\mbox{NB}(r_{k},p_{j}).$ The GNBP employs row-specific probability parameters $p_{j}$ to model row heterogeneity, and hence $X_{j}$ are conditionally independent but not identically distributed if $p_{j}$ at different rows are set differently. Note that the GNBP is previously proposed in Zhou and Carin (2015), which focuses on finding the conditional posterior of $G$, without considering the marginalization of $G$.

The GNBP hierarchical construction is conceptually simple, but to obtain a random count matrix, we have to marginalize out the gamma process $G\sim\Gamma{\mbox{P}}(G_{0},1/c)$. As it is difficult to directly marginalize $G$ out of the conditional likelihood of the observed $J$ rows as

where ${\boldsymbol{p}}:=(p_{1},\ldots,p_{J})$, we first augment each $n_{jk}\sim\mbox{NB}(r_{k},p_{j})$ under its compound Poisson representation as $n_{jk}\sim\mbox{SumLog}(l_{jk},p_{j}),~{}l_{jk}\sim\mbox{Pois}(r_{k}q_{j}).$

Define $X\sim\mbox{SumLogP}(L,p)$ as a sum-logarithmic process such that $X(A)\sim\mbox{SumLog}(L(A),p)$ for each $A\subset\Omega$. With $X_{j}\sim\mbox{NBP}(G,p_{j})$ augmented as $X_{j}\sim\mbox{SumLogP}(L_{j},p_{j}),~{}L_{j}\sim\mbox{PP}(q_{j}G)$, we may express the joint likelihood of $X_{j}$ and $L_{j}$ as

With $l_{\boldsymbol{\cdot}k}:=\sum_{j=1}^{J}l_{jk}$, similar to the analysis in Section A, we can reexpress the likelihood as

Although not obvious, one may verify that (10) defines the PMF of a compound random count matrix, which can be generated via

Let $\sigma(1),\ldots,\sigma(J)$ denote a random permutation of the column indices. If $p_{j}$ are set differently for different rows, then $\mbox{Mult}(l_{\boldsymbol{\cdot}k},{q_{\sigma(1)}}/{q_{\boldsymbol{\cdot}}},\ldots,{q_{\sigma(J)}}/{q_{\boldsymbol{\cdot}}})\buildrel d\over{\neq}\mbox{Mult}(l_{\boldsymbol{\cdot}k},{q_{1}}/{q_{\boldsymbol{\cdot}}},\ldots,{q_{J}}/{q_{\boldsymbol{\cdot}}})$ and hence the introduced random count matrix no longer maintains row exchangeability.

Comparing (B.3) with (2.1.2), one may identify several key differences between the GNBP and NBP random count matrices. First, one may increase $p_{j}$ to encourage the $j$th row to have larger counts than the others. Second, both $n_{jk}$ and the column sum $n_{\boldsymbol{\cdot}k}$ are generated from compound distributions. In fact, if we let $p_{j}\equiv 1-e^{-1}$, then the matrix $\{l_{jk}\}_{jk}$ in (B.3) is exactly a NBP random count matrix, and the GNBP builds its random matrix using $n_{jk}\sim\mbox{SumLog}(l_{jk},p_{j})$.

The sequential construction of a GNBP random count matrix can be intuitively explained as drawing dishes, drawing tables at each dish, and then drawing customers at each table. Similar to the definition of ${{\bf N}}^{+}_{J+1}$, we let ${{\bf L}}^{+}_{J+1}$ represent the new row and columns added to ${{\bf L}}_{J}$. Using (10), following the analysis in Section 2.1, one may show with direct calculation that

Thus to add a new row, we first draw $\mbox{NB}[l_{\boldsymbol{\cdot}k},{q_{J+1}}/{(c+q_{\boldsymbol{\cdot}}+q_{J+1})}]$ tables at existing columns (dishes); we then draw $K^{+}_{J+1}\sim\mbox{Pois}\{\gamma_{0}[\ln(c+q_{\boldsymbol{\cdot}}+q_{J+1})-\ln(c+q_{\boldsymbol{\cdot}})]\}$ new dishes, each of which is associated with $\mbox{Log}[{q_{J+1}}/{(c+q_{\boldsymbol{\cdot}}+q_{J+1})}]$ tables; we further draw $\mbox{Log}(p_{J+1})$ customers at each table and aggregate the counts across the tables of the same dish as $n_{(J+1)k}=\sum_{t=1}^{l_{(J+1)k}}n_{(J+1)kt};$ and in the final step, we insert the $K^{+}_{J+1}$ new columns into the $K_{J}$ original columns without reordering, which again is a one to $K_{J+1}!/\left(K_{J}!\ K^{+}_{J+1}!\right)$ mapping. We emphasize that the number of tables (customers) for a new dish, which follows a logarithmic (sum-logarithmic) distribution, must be at least one; the implication is that there are infinite many dishes that have not yet been ordered by any of the tables seated by existing customers. The sequential construction provides a convenient way to construct a GNBP random count matrix one row at a time.

With the latent counts $l_{(J+1)k}$ marginalized out, one may show that the predictive distribution for ${{\bf N}}^{+}_{J+1}$, given ${{\bf N}}_{J}$ and ${{\bf L}}_{J}$, can be expressed in terms of the Poisson, LogLog and GNB distributions as

B.2 Inference for parameters

Both the GNB and LogLog distributions have complicated PMFs involving Stirling numbers of the first kind and it seems difficult to infer their parameters. Fortunately, using the likelihoods (B.1) and (10) and the data augmentation techniques developed for the negative binomial distribution (Zhou and Carin, 2015), we are able to derive closed-form conditional posteriors for the GNBP. To complete the model, we let $\gamma_{0}\sim\mbox{Gamma}(e_{0},{1}/{f_{0}})$, $p_{j}\sim\mbox{Beta}(a_{0},b_{0})$ and $c\sim\mbox{Gamma}(c_{0},1/d_{0})$. We sample the model parameters as

Appendix C Beta-Negative Binomial Process: Details

The GNBP generalizes the NBP by replacing the Poisson process in (A.1) using a negative binomial process and shares the negative binomial dispersion parameters across rows. Exploiting an alternative strategy that shares the negative binomial probability parameters across rows, we construct a BNBP as

where $p_{k}=B(\omega_{k})$ is the weight of the atom $\omega_{k}$ of the beta process $B\sim\mbox{BP}(c,B_{0})$, and $X_{j}\mid B\sim\mbox{NBP}(r_{j},B)$ is a negative binomial process such that $X_{j}(A)=\sum_{k:\omega_{k}\in A}n_{jk},~{}n_{jk}\sim\mbox{NB}(r_{j},p_{k})$ for each $A\subset\Omega$.

With $\boldsymbol{r}:=(r_{1},\ldots,r_{J})$, similar to the analysis in Appendix B, the likelihood of the BNBP can be expressed as

where $p_{*}$ denotes the sum over all the atoms in the absolutely continuous space $\Omega\backslash\mathcal{D}_{J}$ as

and $r_{\boldsymbol{\cdot}}:=\sum_{j=1}^{J}r_{j}$. Using the Lévy-Khintchine theorem and (1), the Laplace transform of $p_{*}$ can be expressed as

where $\psi(x)={\Gamma^{\prime}(x)}/{\Gamma(x)}$ is the digamma function; we define such a random variable as the logbeta random variable

where the PMFs of both the Dirichlet-multinomial (DirMult) and digamma distributions are shown in the Appendix. Note that if $r_{j}$ are set differently for different rows, then $\mbox{DirMult}(n_{\boldsymbol{\cdot}k},r_{\sigma(1)},\ldots,r_{\sigma(J)})\buildrel d\over{\neq}\mbox{DirMult}(n_{\boldsymbol{\cdot}k},r_{1},\ldots,r_{J})$ and hence the corresponding random count matrix no longer maintains row exchangeability.

The sequential construction of a BNBP random count matrix can be intuitively understood as an “ice cream” buffet process (ICBP). Using (11), similar to the analysis in Section 2.1, we have

A related marked BNBP of Zhou et al. (2012), Zhou and Carin (2012) attaches an independent negative binomial dispersion parameter $r_{k}$ for each atom of the beta process, and infers its values under a finite approximation of the beta process; another related BNBP of Broderick et al. (2015) uses a single dispersion parameter $r$ and sets its value empirically. None of these papers, however, marginalize out the beta process to define a prior on column-i.i.d. random count matrices, a challenge tackled in this paper.

Independently of our work, Heaukulani and Roy (2013) also describe the marginalization of the beta process from the negative binomial process, where the obtained BNBP is called the negative binomial Indian buffet process. Although the idea of marginalizing out the beta process is shared by both papers, the techniques and combinatorial arguments used are quite different. Their paper focuses on a special case of the BNBP where a single dispersion parameter $r$ is used for all the $X_{j}$’s. Our model allows row-specific dispersion parameters $r_{j}$, develops an efficient inference scheme for all model parameters, derives the predictive distribution of a new row count vector under a BNBP random count matrix, and also situates the BNBP in the larger family of count-matrix priors derived from negative-binomial processes.

C.2 Inference for parameters

For all the atoms in the absolutely continuous part of the space, $\Omega\backslash\mathcal{D}_{J}$, we have that

Thus the Laplace transform of $(p_{*}|-)$ can be expressed as

and hence we have $(p_{*}|-)\sim\mbox{logBeta}(\gamma_{0},c+r_{\boldsymbol{\cdot}}).$ With its Laplace transform, we sample $(p_{*}|-)$ using the method proposed in Ridout (2009). To complete the model, we let $\gamma_{0}\sim\mbox{Gamma}(e_{0},{1}/{f_{0}})$, $r_{j}\sim\mbox{Gamma}(a_{0},b_{0})$ and $c\sim\mbox{Gamma}(c_{0},1/d_{0})$. Using both the conditional likelihood (C.1) and the marginal likelihood (11), and the data augmentation techniques developed in Zhou and Carin (2015), we sample the model parameters as

as the proposal distribution in an independence chain Metropolis-Hastings sampling step. One may also sample $c$ using a griddy-Gibbs sampler (Ritter and Tanner, 1992).

Appendix D Some useful distributions

Direct calculation shows that the logarithmic mixed sum-logarithmic (LogLog) distribution, expressed as n\sim\mbox{SumLog}(l,p),~{}l\sim\mbox{Log}\Big{(}\frac{-\ln(1-p)}{c-\ln(1-p)}\Big{)}, has PMF

for $n\in\{1,2,\ldots\}$; and the negative binomial mixed sum-logarithmic distribution, expressed as $n\sim\mbox{SumLog}(l,p),~{}l\sim\mbox{NB}\left(e,\frac{-\ln(1-p)}{c-\ln(1-p)}\right),$ has PMF

for $n\in\{0,1,\ldots\}$. The iterative calculation of ${|s(n,l)|}/{n!}$ under the logarithmic scale is described in Appendix E. Using (2), one may show that the negative binomial mixed sum-logarithmic distribution shown above is equivalent to a gamma mixed negative binomial (GNB) distribution, generated by $n\sim\mbox{NB}(r,p),~{}r\sim\mbox{Gamma}(e,1/c)$. Note that $n\sim\mbox{LogLog}(c,p)$ is the limit of $n\sim\mbox{GNB}(e,c,p)$ as $e\rightarrow 0$, conditioning on $n>0$, thus it can be considered as a truncated GNB distribution.

The Dirichlet-multinomial (DirMult) distribution (Mosimann, 1962, Madsen et al., 2005) is a Dirichlet mixed multinomial distribution, with PMF

and the digamma distribution (Sibuya, 1979) has PMF

where $n=1,2,\ldots$. Since the beta-negative binomial (BNB) distribution has PMF

one may show that conditioning on $n>0$, $n\sim\mbox{BNB}(r,e,c)$ becomes $n\sim\mbox{Digam}(r,c)$ as $e\rightarrow 0$. Thus the digamma distribution can be considered as a truncated BNB distribution.

Since the Laplace transform of the logbeta random variable $p_{*}\sim\mbox{logBeta}(\gamma_{0},c)$ can be reexpressed as

we can generate $p_{*}\sim\mbox{logBeta}(\gamma_{0},c)$ as an infinite sum of independent compound Poisson random variables as

Appendix E Calculating Stirling Numbers of the First Kind

The unsigned Stirling numbers of the first kind $|s(n,l)|$ appear in the predictive distribution for the GNBP. It is numerically unstable to recursively calculate $|s(n,l)|$ based on $|s(n,l)|=(n-1)|s(n-1,l)|+|s(n-1,l-1)|$, as $|s(n,l)|$ would rapidly reach the maximum value allowed by a finite precision machine as $n$ increases. Denoting

we iteratively calculate $g(n,l)$ with $g(n,1)=\ln(n-1)-\ln({n})+\ln g(n-1,1)$, $g(n,n)=g(n-1,n-1)-\ln n$, and

for $2\leq l\leq n-1$. This approach is found to be numerically stable.