The combinatorial structure of beta negative binomial processes

Introduction

Let $c>0$, let $\widetilde{B}_{0}$ be a non-atomic, finite measure on $\Omega$, and let $\Pi$ be a Poisson (point) process on $\Omega\times(0,1]$ with intensity

As this intensity is non-atomic and merely $\sigma$-finite, $\Pi$ will have an infinite number of atoms almost surely (a.s.), and so we may write $\Pi=\sum_{j=1}^{\infty}\delta_{(\gamma_{j},b_{j})}$ for some a.s. unique random elements $b_{1},b_{2},\ldots$ in $(0,1]$ and $\gamma_{1},\gamma_{2},\ldots$ in $\Omega$. From $\Pi$, construct the random measure

which is a beta process . The construction of $B$ ensures that the random variables $B(A_{1}),\ldots,B(A_{k})$ are independent for every finite, disjoint collection $A_{1},\ldots,A_{k}\in\mathcal{A}$, and $B$ is said to be completely random or equivalently, have independent increments . We review completely random measures in Section 2.

The conjugacy of the family of beta distributions with various other exponential families carries over to beta processes and randomizations by probability kernels lying in these same exponential families. The beta process is therefore a convenient choice for further randomizations, or in the language of Bayesian nonparametrics, as a prior stochastic process. For example, previous work has focused on the (simple) point process that takes each atom $\gamma_{j}$ with probability $b_{j}$ for every $j\geq 1$, which is, conditioned on $B$, called a Bernoulli process (with base measure $B$) . In this article, we study the point process

where the random variables $\zeta_{1},\zeta_{2},\ldots$ are conditionally independent given $B$ and

for some parameter $r>0$. Here, $\operatorname{NB}(r,p)$ denotes the negative binomial distribution with parameters $r>0$, $p\in(0,1]$, whose probability mass function (p.m.f.) is

where $(a)_{n}:=a(a+1)\cdots(a+n-1)$ with $(a)_{0}:=1$ is the $n$th rising factorial. Note that, conditioned on $B$, the point process $X$ is the (fixed component) of a negative binomial process . Unconditionally, $X$ is the ordinary component of a beta negative binomial process, which we formally define in Section 2.

Existing constructions for beta negative binomial processes truncate the number of atoms in the underlying beta process and typically use slice sampling to remove the error introduced by this approximation asymptotically . In this work, we instead provide a construction for the beta negative binomial process directly, avoiding a representation of the underlying beta process. To this end, note that while the beta process $B$ has a countably infinite number of atoms a.s., it can be shown that $B$ is still an a.s. finite measure . It follows as an easy consequence that the point process $X$ is a.s. finite as well and, therefore, has an a.s. finite number of atoms, which we represent with a Poisson process. The atomic masses are then characterized by the digamma distribution, introduced by Sibuya , which has p.m.f. (for parameters $r,\theta>0$) given by

where $\psi(a):=\Gamma^{\prime}(a)/\Gamma(a)$ denotes the digamma function. In Section 3, we prove the following:

Let $Y$ be a Poisson process on $(\Omega,\mathcal{A})$ with finite intensity

where $X$ is the beta negative binomial process defined in equation (4).

The probability mass function of $(M_{h})_{h\in\mathcal{H}_{n}}$ is

where $S({h}):=\sum_{j\leq n}h(j)$, for every $h\in\mathcal{H}_{n}$, and $T:=\widetilde{B}_{0}(\Omega)>0$.

selects $\operatorname{Poisson}(c\gamma[\psi(c+r)-\psi(c)])$ distinct dishes, taking $\operatorname{digamma}(r,c)$ servings of each dish, independently.

selects $\operatorname{Poisson}(c\gamma[\psi(c+(n+1)r)-\psi(c+nr)])$ new dishes to taste, taking $\operatorname{digamma}(r,c+nr)$ servings of each dish, independently.

The interpretation here is that, for every $h\in\mathcal{H}_{n}$, the count $M_{h}$ is the number of dishes $k$ such that, for every $j\leq n$, customer $j$ took $h(j)$ servings of dish $k$. Then the sum $S({h})$ in equation (10) is the total number of servings taken of dish $k$ by the first $n$ customers. Because the NB-IBP is the combinatorial structure of a conditionally i.i.d. process, its distribution, given in Theorem 2, must be invariant to every permutation of the customers. We can state this property formally as follows.

Let $\pi$ be a permutation of $[n]:=\{1,\ldots,n\}$, and, for $h\in\mathcal{H}_{n}$, note that the composition $h\circ\pi\in\mathcal{H}_{n}$ is given by $(h\circ\pi)(j)=h(\pi(j))$, for every $j\leq n$. Then

Preliminaries

Here, we review completely random measures and formally define the negative binomial and beta negative binomial processes. We provide characterizations via Laplace functionals and conclude the section with a discussion of related work.

Let $\mathcal{M}(\Omega,\mathcal{A})$ denote the space of $\sigma$-finite measures on $(\Omega,\mathcal{A})$ equipped with the $\sigma$-algebra generated by the projection maps $\mu\mapsto\mu(A)$ for all $A\in\mathcal{A}$. A random measure $\xi$ on $(\Omega,\mathcal{A})$ is a random element in $\mathcal{M}(\Omega,\mathcal{A})$, and we say that $\xi$ is completely random or has independent increments when, for every finite collection of disjoint, measurable sets $A_{1},\ldots,A_{n}\in\mathcal{A}$, the random variables $\xi(A_{1}),\ldots,\xi(A_{n})$ are independent. Here, we briefly review completely random measures; for a thorough treatment, the reader should consult Kallenberg , Chapter 12, or the classic text by Kingman . Every completely random measure $\xi$ can be written as a sum of three independent parts

called the diffuse, fixed, and ordinary components, respectively, where:

$\bar{\xi}$ is a non-random, non-atomic measure;

In this article, we will only study purely-atomic completely random measures, which therefore have no diffuse component. It follows that we may characterize the law of $\xi$ by (1) the distributions of the atomic masses in the fixed component, and (2) the intensity of the Poisson process underlying the ordinary component.

2 Definitions

By a base measure on $(\Omega,\mathcal{A})$, we mean a $\sigma$-finite measure $B$ on $(\Omega,\mathcal{A})$ such that $B\{s\}\leq 1$ for all $s\in\Omega$. For the remainder of the article, fix a base measure $B_{0}$. We may write

and an ordinary component with intensity measure

It is straightforward to show that a beta process is itself a base measure with probability one. This definition of the beta process generalizes the version given in the introduction to a non-homogeneous process with a fixed component. Likewise, we generalize our earlier definition of a negative binomial process to include an ordinary component.

and an ordinary component with intensity measure

The fixed component in this definition was given by Broderick et al. and Zhou et al. (and by Thibaux for the case $r=1$). Here, we have additionally defined an ordinary component, following intuitions from Roy .

The law of a random measure is completely characterized by its Laplace functional, and this representation is often simpler to manipulate: From Campbell’s theorem, or a version of the Lévy–Khinchin formula for Borel spaces, one can show that the Laplace functional of $X$ is

Finally, we define beta negative binomial processes via their conditional law.

A random measure $X$ on $(\Omega,\mathcal{A})$ is a beta negative binomial process with parameter $r>0$, concentration function $c$, and base measure $B_{0}$, written

This characterization was given by Broderick et al. and can be seen to match a special case of the model in Zhou et al. (see the discussion of related work in Section 2.3). It is straightforward to show that a beta negative binomial process is also completely random, and that its Laplace functional is given by

3 Related work

The term “negative binomial process” has historically been reserved for processes with negative binomial increments – a class into which the process we study here does not fal – and these processes have been long-studied in probability and statistics. We direct the reader to Kozubowski and Podgórski for references.

One way to construct a process with negative binomial increments is to rely upon the fact that a negative binomial distribution is a gamma mixture of Poisson distributions. In particular, similarly to the construction by Lo , consider a Cox process $X$ directed by a gamma process $G$ with finite non-atomic intensity. So constructed, $X$ has independent increments with negative binomial distributions. Like the beta process (with a finite intensity underlying its ordinary component), the gamma process has, with probability one, a countably infinite number of atoms but a finite total mass, and so the Cox process $X$ is a.s. finite as well. Despite similarities, a comparison of Laplace functionals shows that the law of $X$ is not that of a beta negative binomial process. Using an approach directly analogous to the derivation of the IBP in , Titsias characterizes the combinatorial structure of a sequence of point processes that, conditioned on $G$, are independent and identically distributed to the Cox process $X$. See Section 4 for comments. This was the first count analogue of the IBP; the possibility of a count analogue arising from beta negative binomial processes was first raised by Zhou et al. , who described the distribution of the number of new dishes sampled by each customer. Recent work by Zhou, Madrid and Scott , independent of our own and proceeding along different lines, describes a combinatorial process related to the NB-IBP (following a re-scaling of the beta process intensity).

Finally, we note that another negative binomial process without negative binomial increments was defined on Euclidean space by Barndorff-Nielsen and Yeo and extended to general spaces by Grégoire and Wolpert and Ickstadt . These measures are generally Cox processes on $(\Omega,\mathcal{A})$ directed by random measures of the form

Constructing beta negative binomial processes

Before providing a finitary construction for the beta negative binomial process, we make a few remarks on the digamma distribution. For the remainder of the article, define $\lambda_{r,\theta}:=\psi(\theta+r)-\psi(\theta)$ for some $r,\theta>0$. Following a representation by Sibuya , we may relate the digamma and beta negative binomial distributions as follows: Let $Z\sim\operatorname{digamma}(r,\theta)$ and define $W:=Z-1$, the latter of which has p.m.f.

With digamma random variables, we provide a finitary construction for the beta negative binomial process. The following result generalizes the statement given by Theorem 1 (in the Introduction) to a non-homogeneous process, which also has a fixed component.

Let $r>0$, and let $\vartheta:=(\vartheta_{s})_{s\in\mathscr{A}}$ be a collection of independent random variables with

Let $Y$ be a Poisson process on $(\Omega,\mathcal{A})$, independent from $\vartheta$, with (finite) intensity

Then by the chain rule of conditional expectation, complete randomness, and Campbell’s theorem,

which is the desired form of the Laplace functional. ∎

where $S_{n}:=\sum_{i=1}^{n}X_{i}$ and $c_{n}(s):=c(s)+S_{n}\{s\}+nr$, for $s\in\Omega$.

We may therefore construct this exchangeable sequence of beta negative binomial processes with Theorem 4.

Combinatorial structure

Let $\hslash\in\mathcal{H}_{n}$, and define $\mathcal{H}_{n+1}^{(\hslash)}:=\{h\in\mathcal{H}_{n+1}:\forall j\leq n,h(j)=\hslash(j)\}$ to be the collection of histories in $\mathcal{H}_{n+1}$ that agree with $\hslash$ on the first $n$ entries. Then note that

that is, the multiplicities $(M_{h})_{h\in\mathcal{H}_{n+1}}$ at stage $n+1$ completely determine the multiplicities $(M_{\hslash})_{\hslash\in\mathcal{H}_{n}}$ at all earlier stages. It follows that

where $m_{\hslash}=\sum_{h\in\mathcal{H}_{n+1}^{(\hslash)}}m_{h}$ for $\hslash\in\mathcal{H}_{n}$. The structure of equation (LABEL:eqJointCombStruct) suggests an inductive proof for Theorem 2.

a Poisson random variable $\kappa$ with mean $cT\lambda_{r,c}$, where $T:=\widetilde{B}_{0}(\Omega)<\infty$;

an i.i.d. collection of a.s. unique random elements $\gamma_{1},\gamma_{2},\ldots$ in $\Omega$;

and $\kappa=\sum_{h\in\mathcal{H}_{1}}M_{h}$ a.s. Therefore,

Because $\zeta_{1},\zeta_{2},\ldots$ are i.i.d., the collection $(M_{h})_{h\in\mathcal{H}_{1}}$ has a multinomial distribution conditioned on its sum $\kappa$. Namely, $M_{h}$ counts the number of times, in $\kappa$ independent trials, that the multiplicity $h(1)$ arises from a $\operatorname{digamma}(r,c)$ distribution. In particular,

𝑛1h\in\mathcal{H}_{n+1} Let $S_{n}:=\sum_{j=1}^{n}X_{j}$. Recall that $s(\hslash):=\sum_{j\leq n}\hslash(j)$ for $\hslash\in\mathcal{H}_{n}$. We may write

a Poisson random variable $\kappa$ with mean $cT\lambda_{r,c+nr}$;

an i.i.d. collection of a.s. unique random elements $\gamma_{1},\gamma_{2},\ldots$ in $\Omega$, a.s. distinct also from $\omega$;

all mutually independent and independent of $X_{[n]}$, such that

Conditioned on $X_{[n]}$, the first and second terms on the right-hand side correspond to the fixed and ordinary components of $X_{n+1}$, respectively. Let

be the set of histories $h$ for which $h(n+1)$ is the first non-zero element. Then, with probability one,

By the stated independence of the variables above, we have

Let $\mathcal{H}^{+}_{n+1}:=\bigcup_{\hslash\in\mathcal{H}_{n}}\mathcal{H}_{n+1}^{(\hslash)}$. For every $\hslash\in\mathcal{H}_{n}$, the random variables $\vartheta_{\hslash,1},\vartheta_{\hslash,2},\ldots$ are i.i.d., and therefore, conditioned on $M_{\hslash}$, the collection $(M_{h})_{h\in\mathcal{H}_{n+1}^{(\hslash)}}$ has a multinomial distribution. In particular, the product term in equation (45) is given by

The p.m.f. of the beta negative binomial distribution is given by

for positive parameters $r,\alpha$, and $\beta$, where $\mathcal{B}(\alpha,\beta):=\Gamma(\alpha)\Gamma(\beta)/\Gamma(\alpha+\beta)$ denotes the beta function. We have that $\kappa=\sum_{h\in\mathcal{H}^{(0)}_{n+1}}M_{h}$ a.s., and therefore

Because $\zeta_{1},\zeta_{2},\ldots$ are i.i.d., conditioned on the sum $\kappa$, the collection $(M_{h})_{h\in\mathcal{H}^{(0)}_{n+1}}$ has a multinomial distribution, and so

In the first product term on the right-hand side of equation (50), note that, for every $h\in\mathcal{H}^{+}_{n+1}$,

where for the last equality, we have used the fact that $h(j)=0$ for every $j\leq n$ and $h\in\mathcal{H}^{(0)}_{n+1}$. Note that $\sum_{\hslash\in\mathcal{H}_{n}}m_{\hslash}+\sum_{h\in\mathcal{H}^{(0)}_{n+1}}m_{h}=\sum_{h\in\mathcal{H}_{n+1}}m_{h}$. Then equation (50) is equal to

Noting that $\sum_{j=1}^{n+1}[\psi(c+jr)-\psi(c+(j-1)r)]=\psi(c+(n+1)r)-\psi(c)$, we obtain the expression in equation (10) for $n+1$, as desired.

Applications in Bayesian nonparametrics

In Bayesian latent feature models, we assume that there exists a latent set of features and that each data point possesses some (finite) subset of the features. The features then determine the distribution of the observed data. In a nonparametric setting, exchangeable sequences of simple point processes can serve as models for the latent sets of features. Similarly, exchangeable sequences of point processes, like those that can be constructed from beta negative binomial processes, can serve as models of latent multisets of features. In particular, atoms are features and their (integer-valued) masses indicate multiplicity. In this section, we develop posterior inference procedures for exchangeable sequences of beta negative binomial processes.

A convenient way to represent the combinatorial structure of an exchangeable sequence of point processes is via an array/matrix $W$ of non-negative integers, where the rows correspond to point processes and columns correspond to atoms appearing among the point processes. Informally, given an enumeration of the set of all atoms appearing in $X_{[n]}$, the entry $W_{i,j}$ associated with the $i$th row and $j$th column is the multiplicity/mass of the atom labeled $j$ in the $i$th point process $X_{i}$.

All that remains is to order the columns of $W$. Every total order on $\mathcal{H}_{n}$ induces a unique ordering of the columns of $W$. Titsias defined a unique ordering in this way, analogous to the left-ordered form defined by Griffiths and Ghahramani for the IBP. In particular, for $h,h^{\prime}\in\mathcal{H}_{n}$, let $\preceq$ denote the lexicographic order given by: $h\preceq h^{\prime}$ if and only if $h=h^{\prime}$ or $h(\eta)<h^{\prime}(\eta)$, where $\eta$ is the first coordinate where $h$ and $h^{\prime}$ differ. We say $W$ is left-ordered when its columns are ordered according to $\preceq$. Because there is a bijection between combinatorial structures $(M_{h})_{h\in\mathcal{H}_{n}}$ and their unique representations by left-ordered arrays, the probability mass function of $W$ is given by equation (10).

Other orderings have been introduced in the literature: If we permute the columns of $W$ uniformly at random, then $W$ is the analogue of the uniform random labeling scheme described by Broderick, Pitman and Jordan for the IBP. Note that the number of distinct ways of ordering the $\kappa$ columns is given by the multinomial coefficient

where the denominator arises from the fact that there are $M_{h}$ indistinguishable columns for every history $h\in\mathcal{H}_{n}$. The following result is then immediate:

An array representation makes it easy to visualize some properties of the model. For example, in Figure 1 we display several simulations from the NB-IBP with varying values of the parameters $T,c$, and $r$. The columns are displayed in the order of first appearance, and are otherwise ordered uniformly at random. (A similar ordering was used by Griffiths and Ghahramani to introduce the IBP.) The relationship of the model to the values of $T$ and $c$ are similar to the characteristics described by Ghahramani, Griffiths and Sollich for the IBP, with the parameter $r$ providing flexibility with respect to the counts in the array. In particular, the total number of features, $\kappa$, is Poisson distributed with mean $cT[\psi(c+nr)-\psi(c)]$, which increases with $T$, $c$, and $r$. From the NB-IBP, we know that the expected number of features for the first (and therefore, by exchangeability, every) row is $T$. Because of the ordering we have chosen here, the rows are not exchangeable, despite the sequence $X_{[n]}$ being exchangeable. (In contrast, a uniform random labeling $W$ is row exchangeable and, conditioned on $\kappa$, column exchangeable.) Finally, note that the mean of the $\operatorname{digamma}(r,c)$ distribution exists for $c>1$ and is given by

which increases with $r$ and decreases with $c$. This is the expected multiplicity of each feature for the first row, which again, by exchangeability, must hold for every row. We may therefore summarize the effects of changing each of these parameters (as we hold the others constant) as follows:

Increasing the mass parameter $T$ increases both the expected total number of features and the expected number of features per row, while leaving the expected multiplicities of the features unchanged.

Increasing the concentration parameter $c$ increases the expected total number of features and decreases the expected multiplicites of the features, while leaving the expected number of features per row unchanged.

Increasing the parameter $r$ increases both the expected total number of features and the expected multiplicities of the features, while leaving the expected number of features per row unchanged.

These effects can be seen in the first, second, and third rows of Figure 1, respectively. We note that $r$ has a weak effect on the expected total number of features (seen in the third row of Figure 1), and $c$ has a weak effect on the expected multiplicities of the features (seen in the second row of Figure 1). The model may therefore be effectively tuned with $T$ and $c$ determining the size and density of the array, and $r$ determining the multiplicities. The most appropriate model depends on the application at hand, and in Section 5.3 we discuss how these parameters may be inferred from data.

2 Examples

Latent feature models with associated multiplicities and unbounded numbers of features have found several applications in Bayesian nonparametric statistics, and we now provide some examples. In these applications, the features represent latent objects or factors underlying a dataset comprised of $n$ groups of measurements $y_{1},\ldots,y_{n}$, where each group $y_{i}$ is comprised of $D_{i}$ measurements $y_{i}=(y_{i,1},\ldots,y_{i,D_{i}})$. In particular, $W_{i,j}$ denotes the number of instances of object/factor $j$ in group $i$.

These nonparametric latent feature representations lend themselves naturally to mixture models with an unbounded number of components. For example, consider a variant of the models by Sudderth et al. and Titsias for a dataset of $n$ street camera images where the latent features are interpreted as object classes that may appear in the images, such as “building”, “car”, “road”, etc. The count $W_{i,j}$ models the relative number of times object class $j$ appears in image $i$. For every $i\leq n$, image $y_{i}$ consists of $D_{i}$ local patches $y_{j,1},\ldots,y_{j,D_{i}}$ detected in the image, which are (collections of) continuous variables representing, for example, color, hue, location in the image, etc. Let $\kappa$ be the number of columns of $W$, that is, the number of features. The local patches in image $i$ are modeled as conditionally i.i.d. draws from a mixture of $S_{i}=\sum_{j=1}^{\kappa}W_{i,j}$ Gaussian distributions, where $W_{i,j}$ of these components are associated with feature $j$. For $k=1,2,\ldots,$ let $\Theta_{i}^{(j,k)}:=(m_{i}^{(j,k)},\Sigma_{i}^{(j,k)})$ denote the mean and covariance of the Gaussian components associated with feature $j$ for image $i$. Let $z_{i,d}=(j,k)$ when $y_{i,d}$ is assigned to component $k\leq W_{i,j}$ associated with feature $j\leq\kappa$. Conditioned on $\Theta:=(\Theta_{i}^{(j,k)})_{i\leq n,j\leq\kappa,k\leq W_{i,j}}$ and the assignments $z:=(z_{i,d})_{i\leq n,d\leq D_{i}}$, the distribution of the measurements admits a conditional density

To share statistical strength across images, the parameters $\Theta_{i}^{(j,k)}$ are given a hierarchical Bayesian prior:

A typical choice for $\nu(\cdot)$ is the family of Gaussian–inverse-Wishart distributions with feature-specific parameters $\Theta^{(j)}$ drawn i.i.d. from a distribution $\nu_{0}$. Finally, for every image $i\leq n$, conditioned on $W$, the assignment variables $z_{i,1},\ldots,z_{i,D_{i}}$ for the local patches in image $n$ are assumed to form a multivariate Pólya urn scheme, arising from repeated draws from a Dirichlet-distributed probability vector over $\{(j,k):j\leq\kappa,k\leq W_{i,j}\}$. The parameters for the Dirichlet distributions are tied in a similar fashion to $\Theta$. The interpretation here is that local patch $d$ in image $i$ is assigned to one of the $S_{i}$ instances of the latent objects appearing in the image. The number of object instances to which a patch may be assigned is specific to the image, but components across all images that correspond to the same feature will be similar.

Latent feature representations are also a natural choice for factor analysis models. Canny and Zhou et al. proposed models for text documents in terms of latent features representing topics. More carefully, let $y_{i,v}$ be the number of occurrences of word $v$ in document $i$. Conditioned on $W$ and a collection of non-negative topic-word weights $\Theta:=(\theta_{j,v})_{j\leq\kappa,v\leq V}$, the word counts are assumed to be conditionally i.i.d. and

In other words, the expected number of occurrences of word $v$ in document $i$ is a linear sum of a small number of weighted factors. The features here are interpreted as topics: words $v$ such that $\theta_{j,v}$ is large are likely to appear many times. There are a total of $\kappa$ topics that are shared across the documents. The topic-word weights $\Theta$ are typically chosen to be i.i.d. Gamma random variates, although there may be reason to prefer priors with dependency enforcing further sparsity. This general setup has been applied to other types of data including, for example, recommendations , where $y_{i,v}$ represents the rating a Netflix user $i$ assigns to a film $v$.

3 Conditional distributions

Let $W$ be a uniform random labelling of a NB-IBP as described in Section 5.1. In the applications described above, computing the posterior distribution of $W$ is the first step towards most other inferential goals. Existing inference schemes use stick-breaking representations, that is, they represent (a truncation of) the beta process underlying $W$. This approach has some advantages, including that the entries of $W$ are then conditionally independent negative binomial random variables. On the other hand, the random variables representing the truncated beta process, as well as the truncation level itself, must be marginalized away using auxiliary variable methods or other techniques . Here, we take advantage of the structure of the NB-IBP and do not represent the beta process. The result is a set of Markov (proposal) kernels analogous to those originally derived for the IBP .

The models described in Section 5.2 associate every feature with a latent parameter. Therefore, conditioned on the number of columns $\kappa$, let $\Theta=(\theta_{1},\ldots,\theta_{\kappa})$ be an i.i.d. sequence drawn from some non-atomic distribution $\nu_{\Theta}$, and assume that the data $y$ admits a conditional density $p(y|W,\Theta)$. We will associate the $j$th column of $W$ with $\Theta_{j}$, and so the pair $(W,\Theta)$ can be seen as an alternative representation for an exchangeable sequence $X_{[n]}$ of beta negative binomial processes. By Bayes’ rule, the posterior distributions admits a conditional density

where $p(W,\Theta)$ is a density for the joint distribution of $(W,\Theta)$. We describe two Markov kernels that leave this distribution invariant. Combined, these kernels give a Markov chain Monte Carlo (MCMC) inference procedure for the desired posterior.

The first kernel resamples individual elements $W_{i,j}$, conditioned on the remaining elements of the array (collectively denoted by $W_{-(i,j)}$), the data $y$, and the parameters $\Theta$. By Bayes’ rule, and the independence of $\Theta$ and $W$ given $\kappa$, we have

Recall that the array $W$ is row-exchangeable, and so, in the language of the NB-IBP, we may associate the $i$th row with the final customer at the buffet. The count $W_{i,j}$ is the number of servings the customer takes of dish $j$, which has been served $S_{j}^{(-i)}:=\sum_{i^{\prime}\neq i}W_{i^{\prime},j}$ times previously. When $S_{j}^{(-i)}>0$, we have

Therefore, we can simulate from the unnormalized, unbounded discrete distribution in equation (60) using equation (61) as a Metropolis–Hastings proposal, or we could use inverse transform sampling where the normalization constant is approximated by an importance sampling estimate.

Following Meeds et al. , the second kernel resamples the number, positions, and values of those singleton columns $j^{\prime}$ such that $S_{j^{\prime}}^{(-i)}=0$. Simultaneously, we propose a corresponding change to the sequence of latent parameters $\Theta$, preserving the relative ordering with the columns of $W$. This corresponding change to $\Theta$ cancels out the effect of the $\kappa!$ term appearing in the p.m.f. of the array $W$. Let $J_{i}$ be the number of singleton columns, that is, let

which we note may be equal to zero. Because we are treating the customer associated with row $i$ as the final customer at the buffet, $J_{i}$ may be interpreted as the number of new dishes sampled by the final customer, in which case, we know that

We therefore propose a new array $W^{*}$ by removing the $J_{i}$ singleton columns from the array and insert $J_{i}^{*}$ new singleton columns at positions drawn uniformly at random, where $J_{i}^{*}$ is sampled from the (marginal) distribution of $J_{i}$ given in equation (63). Like those columns that were removed, each new column has exactly one non-zero entry in the $i$th row: We draw each non-zero entry independently and identically from a $\operatorname{digamma}(r,c+(n-1)r)$ distribution, which matches the distribution of the number of servings the last customer takes of each newly sampled dish.

Finally, we form a new sequence of latent parameters $\Theta^{*}$ by removing those entries from $\Theta$ associated with the $J_{i}$ columns that were removed from $W$ and inserting $J_{i}^{*}$ new entries, drawn i.i.d. from $\nu_{\Theta}$, at the same locations corresponding to the $J_{i}^{*}$ newly introduced columns. Let $\kappa^{*}:=\kappa-J_{i}+J_{i}^{*}$, and note that there were ${\kappa^{*}\choose J_{i}^{*}}$ possible ways to insert the new columns. Therefore, the proposal density is

With manipulations similar to those in the proof of Theorem 2, it is straightforward to show that a Metropolis–Hastings kernel accepts a proposal $(W^{*},\Theta^{*})$ with probability $\min\{1,\alpha^{*}\}$, where

Combined with appropriate Metropolis–Hastings moves that shuffle the columns of $W$ and resample the latent parameters $\Theta$, we obtain a Markov chain whose stationary distribution is the conditional distribution of $W$ and $\Theta$ given the data $y$.

Another benefit of the characterization of the distribution of $W$ in (6) is that numerically integrating over the real-valued concentration, mass, and negative binomial parameters $c$, $T$, and $r$, respectively, are straightforward with techniques such as slice sampling . In the particular case when $T$ is given a gamma prior distribution, say $T\sim\operatorname{gamma}(\alpha,\beta)$ for some positive parameters $\alpha$ and $\beta$, the conditional distribution again falls into the class of gamma distributions. In particular, the conditional density is

Acknowledgements

We thank Mingyuan Zhou for helpful feedback and for pointing out the relation of our work to that of Sibuya . We also thank Yarin Gal and anonymous reviewers for feedback on drafts. This research was carried out while C. Heaukulani was supported by the Stephen Thomas studentship at Queens’ College, Cambridge, with funding also from the Cambridge Trusts, and while D.M. Roy was a research fellow of Emmanuel College, Cambridge, with funding also from a Newton International Fellowship through the Royal Society.

References