arXiv:1401.0062v4 [] 23 Jun 2016

Bernoulli 22(4), 2016, 2301–2324
DOI: 10.3150/15-BEJ729

The combinatorial structure of beta negative
binomial processes
CREIGHTON HEAUKULANI1 and DANIEL M. ROY2
1

Department of Engineering, University of Cambridge, Trumpington Street, Cambridge, CB2
1PZ, United Kingdom. E-mail: ckh28@cam.ac.uk
2
Department of Statistical Sciences, University of Toronto, 100 St. George Street, Toronto, ON,
M5S 3G3, Canada. E-mail: droy@utstat.toronto.edu
We characterize the combinatorial structure of conditionally-i.i.d. sequences of negative binomial
processes with a common beta process base measure. In Bayesian nonparametric applications,
such processes have served as models for latent multisets of features underlying data. Analogously, random subsets arise from conditionally-i.i.d. sequences of Bernoulli processes with a
common beta process base measure, in which case the combinatorial structure is described by
the Indian buffet process. Our results give a count analogue of the Indian buffet process, which
we call a negative binomial Indian buffet process. As an intermediate step toward this goal,
we provide a construction for the beta negative binomial process that avoids a representation
of the underlying beta process base measure. We describe the key Markov kernels needed to
use a NB-IBP representation in a Markov Chain Monte Carlo algorithm targeting a posterior
distribution.
Keywords: Bayesian nonparametrics; Indian buffet process; latent feature models; multisets

1. Introduction
The focus of this article is on exchangeable sequences of multisets, that is, set-like objects
in which repetition is allowed. Let Ω be a complete, separable metric space equipped with
its Borel σ-algebra A and let Z+ := {0, 1, 2, . . .} denote the non-negative integers. By a
point process on (Ω, A), we mean a random measure X on (Ω, A) such that X(A) is a
Z+ -valued random variable for every A ∈ A. Because (Ω, A) is Borel, we may write
X=

X

δγk

(1.1)

k≤κ

for a random element κ in Z+ := Z+ ∪ {∞} and some – not necessarily distinct – random
elements γ1 , γ2 , . . . in Ω. We will take the point process X to represent the multiset of

This is an electronic reprint of the original article published by the ISI/BS in Bernoulli,
2016, Vol. 22, No. 4, 2301–2324. This reprint differs from the original in pagination and
typographic detail.
1350-7265

c

2016 ISI/BS

2

C. Heaukulani and D.M. Roy

its unique atoms γk with corresponding multiplicities X{γk }. We say X is simple when
X{γk } = 1 for all k ≤ κ, in which case X represents a set.
In statistical applications, latent feature models associate each data point yn in a
dataset with a latent point process Xn from an exchangeable sequence of simple point
processes, which we denote by (Xn )n∈N := (X1 , X2 , . . .). The unique atoms among the
sequence (Xn )n∈N are referred to as features, and a data point is said to possess those
features appearing in its associated point process. We can also view these latent feature
models as generalizations of mixture models that allow data points to belong to multiple,
potentially overlapping clusters [2, 10]. For example, in an object recognition task, a
model for a dataset consisting of street camera images could associate each image with
a subset of object classes – for example, “trees”, “cars”, and “houses”, etc. – appearing
in the images. In a document modeling task, a model for a dataset of news articles could
associate each document with a subset of topics – for example, “politics”, “Europe”, and
“economics”, etc. – discussed in the documents. Recent work in Bayesian nonparametrics
utilizing exchangeable sequences of simple point processes have focused on the Indian
buffet process (IBP) [7, 10], which characterizes the marginal distribution of the sequence
(Xn )n∈N when they are conditionally-i.i.d. Bernoulli processes, given a common beta
process base measure [11, 24].
If the point processes (Xn )n∈N are no longer constrained to be simple, then data
points can contain multiple copies of features. For example, in the object recognition
task, an image could be associated with two cars, two trees, and one house. In the
document modeling task, an article could be associated with 100 words from the politics
topic, 200 words from the Europe topic, and 40 words from the economics topic. In this
article, we describe a count analogue of the IBP called the negative binomial Indian buffet
processes (NB-IBP), which characterizes the marginal distribution of (Xn )n∈N when it
is a conditionally i.i.d. sequence of negative binomial processes [3, 28], given a common
beta process base measure. This characterization allows us to describe new Markov Chain
Monte Carlo algorithms for posterior inference that do not require numerical integrations
over representations of the underlying beta process.

1.1. Results
e0 be a non-atomic, finite measure on Ω, and let Π be a Poisson (point)
Let c > 0, let B
process on Ω × (0, 1] with intensity
e0 (ds).
(ds, dp) 7→ cp−1 (1 − p)c−1 dpB

(1.2)

As this intensity is non-atomic and merely σ-finite, ΠPwill have an infinite number of
∞
atoms almost surely (a.s.), and so we may write Π = j=1 δ(γj ,bj ) for some a.s. unique
random elements b1 , b2 , . . . in (0, 1] and γ1 , γ2 , . . . in Ω. From Π, construct the random
measure
B :=

∞
X
j=1

b j δγj ,

(1.3)

The combinatorial structure of beta negative binomial processes

3

which is a beta process [11]. The construction of B ensures that the random variables
B(A1 ), . . . , B(Ak ) are independent for every finite, disjoint collection A1 , . . . , Ak ∈ A, and
B is said to be completely random or equivalently, have independent increments [14]. 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 γj with probability bj for every j ≥ 1, which is, conditioned on B, called a
Bernoulli process (with base measure B) [24]. In this article, we study the point process
∞
X
(1.4)
ζj δγj ,
X :=
j=1

where the random variables ζ1 , ζ2 , . . . are conditionally independent given B and
ζj |bj ∼ NB(r, bj ),

j ∈ N,

(1.5)

for some parameter r > 0. Here, NB(r, p) denotes the negative binomial distribution with
parameters r > 0, p ∈ (0, 1], whose probability mass function (p.m.f.) is
NB(z; r, p) :=

(r)z z
p (1 − p)r ,
z!

z ∈ Z+ ,

(1.6)

where (a)n := a(a + 1) · · · (a + n − 1) with (a)0 := 1 is the nth rising factorial. Note that,
conditioned on B, the point process X is the (fixed component) of a negative binomial
process [3, 28]. Unconditionally, X is the ordinary component of a beta negative binomial
process, which we formally define in Section 2.
Conditioned on B, construct a sequence of point processes (Xn )n∈N that are i.i.d.
copies of X. In this case, (Xn )n∈N is an exchangeable sequence of beta negative binomial
processes, and our primary goal is to characterize the (unconditional) distribution of the
sequence. This task is non-trivial because the construction of the point process X in
equation (1.4) is not finitary in the sense that no finite subset of the atoms of B determines X with probability one. In the case of conditionally-i.i.d. Bernoulli processes, the
unconditional distributions of the measures remain in the class of Bernoulli processes,
and so a finitary construction is straightforwardly obtained with Poisson (point) processes. Then the distribution of the sequence, which Thibaux and Jordan [24] showed is
characterized by the IBP, may be derived immediately from the conjugacy between the
classes of beta and Bernoulli processes [11, 13, 24]. While conjugacy also holds between
the classes of beta and negative binomial processes [3, 28], the unconditional law of the
point process X is no longer that of a negative binomial process; instead, it is the law of
a beta negative binomial process.
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 [3, 19, 23, 28]. In this work, we

4

C. Heaukulani and D.M. Roy

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 [11]. 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 [21], which has p.m.f. (for parameters r, θ > 0) given by
digamma(z; r, θ) :=

1
(r)z
z −1 ,
ψ(r + θ) − ψ(θ) (r + θ)z

z ≥ 1,

(1.7)

where ψ(a) := Γ′ (a)/Γ(a) denotes the digamma function. In Section 3, we prove the
following:
Theorem 1.1. Let Y be a Poisson process on (Ω, A) with finite intensity
e0 (ds),
ds 7→ c[ψ(c + r) − ψ(c)]B

(1.8)

Pκ
that is, Y = k=1 δγk for a Poisson random variable κ with mean c[ψ(c + r) −
e0 (Ω) and i.i.d. random variables (γk )k∈N , independent from κ, each with distribuψ(c)]B
e
e0 (Ω). Let (ζk )k∈N be an independent collection of i.i.d. digamma(r, c) random
tion B0 /B
variables. Then
d

X=

κ
X

ζk δγk ,

(1.9)

k=1

where X is the beta negative binomial process defined in equation (1.4).
With this construction and conjugacy (the relevant results are reproduced in Section 4),
characterizing the distribution of (Xn )n∈N is straightforward. However, in applications
we are only interested in the combinatorial structure of the sequence (Xn )n∈N , that is,
the pattern of sharing amongst the atoms while ignoring the locations of the atoms
themselves. More precisely, for every n ∈ N, let Hn := Zn+ \ {0n } be the set of all lengthn sequences of non-negative integers, excluding the all-zero sequence. Elements in Hn
are called histories, and can be thought of as representations of non-empty multisets
of [n] := {1, . . . , n}. For every h ∈ Hn , let Mh be the number of elements s ∈ Ω such
that Xj {s} = h(j) for all j ≤ n. By the combinatorial structure of a finite subsequence
X[n] := (X1 , . . . , Xn ), we will mean the collection (Mh )h∈Hn of counts, which together can
be understood as representations of multisets of histories. These counts are combinatorial
in the following sense: Let φ: (Ω, A) → (Ω, A) be a Borel automorphism on (Ω, A), that
is, a measurable permutation of Ω whose inverse is also measurable, and define the
transformed processes Xjφ := Xj ◦ φ−1 , for every j ≤ n, where each atom s is repositioned
to φ(s). The collection (Mh )h∈Hn is invariant to this transformation, and it is in this
sense that they only capture the combinatorial structure. In Section 4, we prove the
following.

The combinatorial structure of beta negative binomial processes

5

Theorem 1.2. The probability mass function of (Mh )h∈Hn is
P{Mh = mh : h ∈ Hn }
P

(cT )
= Q

h∈Hn

h∈Hn

where S(h) :=

mh

mh !
P

exp(−cT [ψ(c + nr) − ψ(c)])

j≤n h(j),

Y

h∈Hn

"

(1.10)
#mh

n
Γ(S(h))Γ(c + nr) Y (r)h(j)
Γ(c + nr + S(h)) j=1 h(j)!

,

e0 (Ω) > 0.
for every h ∈ Hn , and T := B

As one would expect, equation (1.10) is reminiscent of the p.m.f. for the IBP, and
indeed the collection (Mh )h∈Hn is characterized by what we call the negative binomial
Indian buffet process, or NB-IBP. Let beta-NB(r, α, β) denote the beta negative binomial
distribution (with parameters r, α, β > 0), that is, we write Z ∼ beta-NB(r, α, β) if there
exists a beta random variable p ∼ beta(α, β) such that Z|p ∼ NB(r, p). In the NB-IBP, a
sequence of customers enters an Indian buffet restaurant:
• The first customer
– selects Poisson(cγ[ψ(c + r) − ψ(c)]) distinct dishes, taking digamma(r, c) servings
of each dish, independently.
• For n ≥ 1, the (n + 1)st customer
– takes beta-NB(r, Sn,k , c + nr) servings of each previously sampled dish k; where
Sn,k is the total number of servings taken of dish k by the first n customers;
– selects Poisson(cγ[ψ(c + (n + 1)r) − ψ(c + nr)]) new dishes to taste, taking
digamma(r, c + nr) servings of each dish, independently.
The interpretation here is that, for every h ∈ Hn , the count Mh is the number of
dishes k such that, for every j ≤ n, customer j took h(j) servings of dish k. Then the
sum S(h) in equation (1.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 1.2, must be invariant to every permutation
of the customers. We can state this property formally as follows.
Theorem 1.3 (Exchangeability). Let π be a permutation of [n] := {1, . . . , n}, and, for
h ∈ Hn , note that the composition h ◦ π ∈ Hn is given by (h ◦ π)(j) = h(π(j)), for every
j ≤ n. Then
d

(Mh )h∈Hn = (Mh◦π )h∈Hn .

(1.11)

The exchangeability of the combinatorial structure and its p.m.f. in equation (1.10)
allows us to develop Gibbs sampling techniques analogous to those originally developed
for the IBP [7, 17]. In particular, because the NB-IBP avoids a representation of the
beta process underlying the exchangeable sequence (Xn )n∈N , these posterior inference
algorithms do not require numerical integration over representations of the beta process.
We discuss some of these techniques in Section 5.

6

C. Heaukulani and D.M. Roy

2. 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.

2.1. Completely random measures
Let M(Ω, A) denote the space of σ-finite measures on (Ω, A) equipped with the σalgebra generated by the projection maps µ 7→ µ(A) for all A ∈ A. A random measure
ξ on (Ω, A) is a random element in M(Ω, A), and we say that ξ is completely random
or has independent increments when, for every finite collection of disjoint, measurable
sets A1 , . . . , An ∈ A, the random variables ξ(A1 ), . . . , ξ(An ) are independent. Here, we
briefly review completely random measures; for a thorough treatment, the reader should
consult Kallenberg [12], Chapter 12, or the classic text by Kingman [14]. Every completely
random measure ξ can be written as a sum of three independent parts
X
X
pδs
a.s.,
(2.1)
ξ = ξ̄ +
ϑs δs +
s∈A

(s,p)∈η

called the diffuse, fixed, and ordinary components, respectively, where:
1. ξ̄ is a non-random, non-atomic measure;
2. A ⊆ Ω is a non-random countable set whose elements are referred to as the fixed
atoms and whose masses ϑ1 , ϑ2 , . . . are independent random variables in R+ (the nonnegative real numbers);
3. η is a Poisson process on Ω × (0, ∞) whose intensity Eη is σ-finite and has diffuse
projections onto Ω, that is, the measure (Eη)(· × (0, ∞)) on Ω is non-atomic.
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 ξ 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.2. Definitions
By a base measure on (Ω, A), we mean a σ-finite measure B on (Ω, A) such that B{s} ≤ 1
for all s ∈ Ω. For the remainder of the article, fix a base measure B0 . We may write
X
e0 +
B0 = B
b̄s δs
(2.2)
s∈A

e0 ; a countable set A ⊆ Ω; and constants b̄1 , b̄2 , . . . in
for some non-atomic measure B
1
(0, 1]. As discussed in the Introduction, a convenient model for random base measures
1 Note

e0 (in the Introduction) to be merely σ-finite.
that we have relaxed the condition on B

The combinatorial structure of beta negative binomial processes

7

are beta processes, a class of completely random measures introduced by Hjort [11]. For
the remainder of the article, let c: Ω → R+ be a measurable function, which we call a
concentration function (or parameter when it is constant).
Definition 2.1 (Beta process). A random measure B on (Ω, A) is a beta process
with concentration function c and base measure B0 , written B ∼ BPL (c, B0 ), when it is
purely atomic and completely random, with a fixed component
X
ind
(2.3)
ϑs δs ,
ϑs ∼ beta(c(s)b̄s , c(s)(1 − b̄s )),
s∈A

and an ordinary component with intensity measure
e0 (ds).
(ds, dp) 7→ c(s)p−1 (1 − p)c(s)−1 dpB

(2.4)

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.
Definition 2.2 (Negative binomial process). A point process X on (Ω, A) is a
negative binomial process with parameter r > 0 and base measure B0 , written X ∼
NBP(r, B0 ), when it is purely atomic and completely random, with a fixed component
X

ϑs δs ,

ind

ϑs ∼ NB(r, b̄s ),

(2.5)

s∈A

and an ordinary component with intensity measure
e0 (ds).
(ds, dp) 7→ rδ1 (dp)B

(2.6)

The fixed component in this definition was given by Broderick et al. [3] and Zhou et
al. [28] (and by Thibaux [25] for the case r = 1). Here, we have additionally defined an
ordinary component, following intuitions from Roy [20].
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 Lévy–Khinchin formula for Borel spaces, one can show that the Laplace
functional of X is
 Z
 Y
r
1 − b̄s
e0 (ds)
f 7→ E[e−X(f ) ] = exp − (1 − e−f (s) )rB
,
(2.7)
1 − b̄s e−f (s)
s∈A

R
where f ranges over non-negative measurable functions and X(f ) := f (s)X(ds).
Finally, we define beta negative binomial processes via their conditional law.

8

C. Heaukulani and D.M. Roy

Definition 2.3 (Beta negative binomial process). A random measure X on (Ω, A)
is a beta negative binomial process with parameter r > 0, concentration function c, and
base measure B0 , written
X ∼ BNBP(r, c, B0 ),
if there exists a beta process B ∼ BPL (c, B0 ) such that
X|B ∼ NBP(r, B).

(2.8)

This characterization was given by Broderick et al. [3] and can be seen to match
a special case of the model in Zhou et al. [28] (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
 Z 

r 

1−p
−1
c(s)−1
e0 (ds)
E[e−X(f ) ] = exp −
1−
c(s)p
(1
−
p)
dp
B
1 − pe−f (s)
(2.9)
Y Z  1 − p r
beta(p; c(s)b̄s , c(s)(1 − b̄s )) dp,
×
1 − pe−f (s)
s∈A

for f : Ω → R+ measurable, where we note that the factors in the product term take the
form of the Laplace transform of the beta negative binomial distribution.

2.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 Podgórski [15] 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 [16], 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 [10], Titsias [26] 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. [28], who described the

The combinatorial structure of beta negative binomial processes

9

distribution of the number of new dishes sampled by each customer. Recent work by
Zhou, Madrid and Scott [29], 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 [1] and extended
to general spaces by Grégoire [9] and Wolpert and Ickstadt [27]. These measures are
generally Cox processes on (Ω, A) directed by random measures of the form
Z
ν(t, ds)G(dt),
ds 7→
R+

where G is again a gamma process, this time on R+ , and ν is a probability kernel from
Ω to R+ , for example, the Gaussian kernel.

3. 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
λr,θ := ψ(θ + r) − ψ(θ) for some r, θ > 0. Following a representation by Sibuya [21],
we may relate the digamma and beta negative binomial distributions as follows: Let
Z ∼ digamma(r, θ) and define W := Z − 1, the latter of which has p.m.f.
P{W = w} = (θλr,θ )−1

w+r
beta-NB(w; r, 1, θ),
w+1

w ∈ Z+ .

(3.1)

Deriving the Laplace transform of the law of W is straightforward, and because E[e−tW ] =
et E[e−tZ ], one may verify that the Laplace transform of the digamma distribution is given
by

r 
Z 
1−p
−1
−tZ
1−
Ψr,θ (t) := E[e
] = 1 − λr,θ
p−1 (1 − p)θ−1 dp.
(3.2)
1 − pe−t
The form of equation (3.1) suggests the following rejection sampler, which was first
proposed by Devroye [6], Proposition 2, Remark 1: Let r > 0 and let (Un )n∈N be an i.i.d.
sequence of uniformly distributed random numbers. Let
i.i.d.

(Yn )n∈N ∼ beta-NB(r, 1, θ),
and define η := inf{n ∈ N : max{r, 1} · Un <

Yi +r
Yi +1 }.

Then

Yη + 1 ∼ digamma(r, θ),
and
Eη =

max{r, 1}
;
θ[ψ(r + θ) − ψ(θ)]

Eη < max{r, r−1 }.

10

C. Heaukulani and D.M. Roy

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.1 (in the Introduction) to a non-homogeneous process, which also has a fixed
component.
Theorem 3.1. Let r > 0, and let ϑ := (ϑs )s∈A be a collection of independent random
variables with
ϑs ∼ beta-NB(r, c(s)b̄s , c(s)(1 − b̄s )),

s∈A.

(3.3)

Let Y be a Poisson process on (Ω, A), independent from ϑ, with (finite) intensity
e0 (ds).
ds 7→ c(s)[ψ(c(s) + r) − ψ(c(s))]B

(3.4)

Pκ
Write Y = k=1 δγk for some random element κ in Z+ and a.s. unique random elements γ1 , γ2 , . . . in Ω, and put F := σ(κ, γ1 , γ2 , . . .). Let (ζj )j∈N be a collection of random
variables that are independent from ϑ and are conditionally independent given F , and let
ζj |F ∼ digamma(r, c(γj )),

j ∈ N.

(3.5)

Then
X=

X

ϑs δs +

κ
X

ζj δγj ∼ BNBP(r, c, B0 ).

(3.6)

j=1

s∈A

Proof. We have
F

E [e

−X(f )

]=

Y

E[e

−ϑs f (s)

]×

κ
Y

EF [e−ζj f (γj ) ],

(3.7)

j=1

s∈A

for every f : Ω → R+ measurable. For s ∈ Ω, write g(s) = Ψr,c(s) (f (s)) for the Laplace
transform of the digamma distribution evaluated at f (s), where Ψr,θ (t) is given by equation (3.2). We may then write
κ
Y

EF [e−ζj f (γj ) ] =

κ
Y

g(γj ).

(3.8)

j=1

j=1

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

 Z
Y
−X(f )
−ϑs f (s)
e
(3.9)
E[e
]=
E[e
] × exp − (1 − g(s))c(s)λr,c(s) B0 (ds)
Ω

s∈A

=

Y Z 

s∈A

1−p
1 − pe−f (s)

r

beta(p; c(s)b̄s , c(s)(1 − b̄s )) dp



(3.10)

The combinatorial structure of beta negative binomial processes
11


 Z
r 

1−p
−1
c(s)−1
e0 (ds) ,
1−
× exp −
c(s)p
(1
−
p)
dp
B
1 − pe−f (s)
(0,1]×Ω
which is the desired form of the Laplace functional.



A finitary construction for conditionally-i.i.d. sequences of negative binomial processes
with a common beta process base measure now follows from known conjugacy results. In
particular, for every n ∈ N, let X[n] := (X1 , . . . , Xn ). The following theorem characterizes
the conjugacy between the (classes of) beta and negative binomial processes and follows
from repeated application of the results by Kim [13], Theorem 3.3 or Hjort [11], Corollary 4.1. This result, which is tailored to our needs, is similar to those already given by
Broderick et el. [3] and Zhou et al. [28], and generalizes the result given by Thibaux [25]
for the case r = 1.
Theorem 3.2 (Hjort [11], Zhou et al. [28]). Let B ∼ BPL (c, B0 ) and, conditioned on
B, let (Xn )n∈N be a sequence of i.i.d. negative binomial processes with parameter r > 0
and base measure B. Then for every n ∈ N,


c
1
B|X[n] ∼ BPL cn , B0 + Sn ,
(3.11)
cn
cn
P
where Sn := ni=1 Xi and cn (s) := c(s) + Sn {s} + nr, for s ∈ Ω.

Remark 3.1. It follows immediately that, for every n ∈ N, the law of Xn+1 conditioned
on X1 , . . . , Xn is given by


1
c
(3.12)
Xn+1 |X[n] ∼ BNBP r, cn , B0 + Sn .
cn
cn
We may therefore construct this exchangeable sequence of beta negative binomial
processes with Theorem 3.1.

4. Combinatorial structure
We now characterize the combinatorial structure of the exchangeable sequence X[n] in the
e0 ) is non-atomic. In order to make this precise, we
case when c > 0 is constant and B0 (= B
introduce a quotient of the space of sequences of integer-valued measures. Let n ∈ N and
for any pair U := (U1 , . . . , Un ) and V := (V1 , . . . , Vn ) of (finite) sequences of integer-valued
measures, write U ∼ V when there exists a Borel automorphism φ on (Ω, A) satisfying
Uj = Vj ◦ φ−1 for every j ≤ n. It is easy to verify that ∼ is an equivalence relation. Let
[[U ]] denote the equivalence class containing U . The quotient space induced by ∼ is itself
a Borel space, and can be related to the Borel space of sequences of Z+ -valued measures
by coarsening the σ-algebra to that generated by the functionals
Mh (U1 , . . . , Un ) := #{s ∈ Ω : ∀j ≤ n, Uj {s} = h(j)},

h ∈ Hn , j ≤ n,

(4.1)

12

C. Heaukulani and D.M. Roy

where #A denotes the cardinality of A, and Hn := Zn+ \ {0n } is the space of histories
defined in the Introduction. The collection (Mh )h∈Hn of multiplicities (of histories) corresponding to X[n] , also defined in the Introduction, then satisfies Mh = Mh (X[n] ) for
every h ∈ Hn . The collection (Mh )h∈Hn thus identifies a point in the quotient space
induced by ∼. Our aim is to characterize the distribution of (Mh )h∈Hn , for every n ∈ N.
(ℏ)
Let ℏ ∈ Hn , and define Hn+1 := {h ∈ Hn+1 : ∀j ≤ n, h(j) = ℏ(j)} to be the collection
of histories in Hn+1 that agree with ℏ on the first n entries. Then note that
X
Mℏ =
Mh ,
ℏ ∈ Hn ,
(4.2)
(ℏ)

h∈Hn+1

that is, the multiplicities (Mh )h∈Hn+1 at stage n + 1 completely determine the multiplicities (Mℏ )ℏ∈Hn at all earlier stages. It follows that
P{Mh = mh : h ∈ Hn+1 } = P{Mℏ = mℏ : ℏ ∈ Hn }
× P{Mh = mh : h ∈ Hn+1 |Mℏ = mℏ : ℏ ∈ Hn },
where mℏ =

P

(ℏ)

h∈Hn+1

(4.3)

mh for ℏ ∈ Hn . The structure of equation (4.3) suggests an induc-

tive proof for Theorem 1.2.

4.1. The law of Mh for h ∈ H1
Note that H1 is isomorphic to N and that the collection (Mh )h∈H1 counts the number of
atoms of each positive integer mass. It follows from Theorem 1.1 and a transfer argument
[12], Propositions 6.10, 6.11 and 6.13, that there exists:
e0 (Ω) < ∞;
1. a Poisson random variable κ with mean cT λr,c , where T := B
2. an i.i.d. collection of a.s. unique random elements γ1 , γ2 , . . . in Ω;
3. an i.i.d. collection (ζj )j∈N of digamma(r, c) random variables;
all mutually independent, such that
X1 =

κ
X

ζj δγj

a.s.

j=1

It follows that
Mh = #{j ≤ κ: ζj = h(1)}
and κ =

P

h∈H1

a.s., for h ∈ H1 ,

(4.4)

Mh a.s. Therefore,
P{Mh = mh : h ∈ H1 }

 

X
X
mh .
mh P Mh = mh : h ∈ H1 κ =
=P κ=
h∈H1

h∈H1

(4.5)

The combinatorial structure of beta negative binomial processes

13

Because ζ1 , ζ2 , . . . are i.i.d., the collection (Mh )h∈H1 has a multinomial distribution conditioned on its sum κ. Namely, Mh counts the number of times, in κ independent trials,
that the multiplicity h(1) arises from a digamma(r, c) distribution. In particular,


X
P Mh = mh : h ∈ H1 κ =
mh
h∈H1

(4.6)

P
( h∈H1 mh )! Y
m
[digamma(h(1); r, c) h ].
= Q
(m
!)
h
h∈H1
h∈H1

It follows that

P{Mh = mh : h ∈ H1 }
P

(cT λr,c )
= Q

h∈H1

mh

h∈H1 (mh !)

exp(−cT λr,c )

Y

(4.7)
mh

[digamma(h(1); r, c)

].

h∈H1

4.2. The conditional law of Mh for h ∈ Hn+1
Let Sn :=

Pn

j=1 Xj .

Recall that s(ℏ) :=
Sn =

P

j≤n ℏ(j)

Mℏ
X X

for ℏ ∈ Hn . We may write

s(ℏ)δωℏ,j ,

(4.8)

ℏ∈Hn j=1

for some collection ω := (ωℏ,j )ℏ∈Hn ,j∈N of a.s. distinct random elements in Ω. It follows
from Remark 3.1, Theorem 1.1, and a transfer argument that there exists:
1. a Poisson random variable κ with mean cT λr,c+nr ;
2. an i.i.d. collection of a.s. unique random elements γ1 , γ2 , . . . in Ω, a.s. distinct also
from ω;
3. an i.i.d. collection (ζj )j∈N of digamma(r, c + nr) random variables;
4. for each ℏ ∈ Hn , an i.i.d. collection (ϑℏ,j )j∈N of random variables satisfying
ϑℏ,j ∼ beta-NB(r, s(ℏ), c + nr)

for j ∈ N;

all mutually independent and independent of X[n] , such that
Xn+1 =

Mℏ
X X

ℏ∈Hn j=1

ϑℏ,j δωℏ,j +

κ
X

ζj δγj

a.s.

(4.9)

j=1

Conditioned on X[n] , the first and second terms on the right-hand side correspond to the
fixed and ordinary components of Xn+1 , respectively. Let
(0)

Hn+1 := {h ∈ Hn+1 : h(j) = 0, j ≤ n}

(4.10)

14

C. Heaukulani and D.M. Roy

be the set of histories h for which h(n + 1) is the first non-zero element. Then, with
probability one,
Mh = #{j ≤ κ: ζj = h(n + 1)}

(0)

for h ∈ Hn+1 ,

(4.11)

and
Mh = #{j ≤ Mℏ : ϑℏ,j = h(n + 1)}

(ℏ)

for ℏ ∈ Hn and h ∈ Hn+1 .

(4.12)

By the stated independence of the variables above, we have
P{Mh = mh : h ∈ Hn+1 |Mℏ = mℏ : ℏ ∈ Hn }
Y
(0)
(ℏ)
= P{Mh = mh : h ∈ Hn+1 }
P{Mh = mh : h ∈ Hn+1 |Mℏ = mℏ }.

(4.13)

ℏ∈Hn

S
(ℏ)
+
Let Hn+1
:= ℏ∈Hn Hn+1 . For every ℏ ∈ Hn , the random variables ϑℏ,1 , ϑℏ,2 , . . . are
i.i.d., and therefore, conditioned on Mℏ , the collection (Mh )h∈H(ℏ) has a multinomial
n+1

distribution. In particular, the product term in equation (4.13) is given by
Y

(ℏ)

P{Mh = mh : h ∈ Hn+1 |Mℏ = mℏ }

ℏ∈Hn

Q
(mℏ !)
= Q ℏ∈Hn
h∈H+ (mh !)
n+1

Y

[beta-NB(h(n + 1); r, S(h) − h(n + 1), c + nr)

mh

].

h∈H+
n+1

The p.m.f. of the beta negative binomial distribution is given by
beta-NB(z; r, α, β) =

(r)z B(z + α, r + β)
,
z
B(α, β)

z ∈ Z+ ,

(4.14)

for positive parameters r, α, P
and β, where B(α, β) := Γ(α)Γ(β)/Γ(α + β) denotes the beta
function. We have that κ = h∈H(0) Mh a.s., and therefore
n+1

(0)

P{Mh = mh : h ∈ Hn+1 }


X
=P κ=
mh

(4.15)

(0)
h∈Hn+1


(0)
× P Mh = mh : h ∈ Hn+1 κ =

X

(0)

h∈Hn+1


mh .

The combinatorial structure of beta negative binomial processes

15

Because ζ1 , ζ2 , . . . are i.i.d., conditioned on the sum κ, the collection (Mh )h∈H(0) has a
n+1

multinomial distribution, and so

(0)
P Mh = mh : h ∈ Hn+1 κ =
P
( h∈H(0) mh )!
n+1
= Q
(0) (mh !)
h∈H
n+1

It follows that

Y

X

mh

(0)

h∈Hn+1


(4.16)

[digamma(h(n + 1); r, c + nr)

mh

].

(0)

h∈Hn+1

P{Mh = mh : h ∈ Hn+1 |Mℏ = mℏ : ℏ ∈ Hn }
P

(0)

(cT λr,c+nr ) h∈Hn+1
P
=
( h∈H(0) mh )!

mh

exp(−cT λr,c+nr )

n+1

Q
(mℏ !)
× Q ℏ∈Hn
+
(mh !)
h∈H
n+1

(4.17)

Y

[beta-NB(h(n + 1); r, S(h) − h(n + 1), c + nr)

Y

[digamma(h(n + 1); r, c + nr)

mh

]

h∈H+
n+1

P
( h∈H(0) mh )!
n+1
× Q
(0) (mh !)
h∈H
n+1

mh

].

(0)
h∈Hn+1

Proof of Theorem 1.2. The proof is by induction. The p.m.f. P{Mh = mh : h ∈ H1 }
is given by equation (4.7), which agrees with equation (1.10) for the case n = 1. The
conditional p.m.f. P{Mh = mh : h ∈ Hn+1 |Mℏ = mℏ : ℏ ∈ Hn } is given by equation (4.17).
By the inductive hypothesis, the p.m.f. P{Mℏ = mℏ : ℏ ∈ Hn } is given by equation (1.10).
Then by equation (4.3), we have
P{Mh = mh : h ∈ Hn+1 }
(

=

(cT )

P

ℏ∈Hn

Q

mℏ )

h∈H+
n+1

×

Y

h∈H+
n+1

"

(

P

(0)

(cT λr,c+nr ) h∈Hn+1
Q
(mh !) h∈H(0) (mh !)

mh )

exp −cT

j=1

n+1

B(S(h) − h(n + 1), c + nr)

n
Y
(r)h(j)
h(j)!
j=1

× beta-NB(h(n + 1); r, S(h) − h(n + 1), c + nr)
×

Y

(0)
h∈Hn+1

n+1
X

[digamma(h(n + 1); r, c + nr)]

mh

.

#mh

λr,c+(j−1)r

!
(4.18)

16

C. Heaukulani and D.M. Roy

In the first product term on the right-hand side of equation (4.18), note that, for every
+
h ∈ Hn+1
,
n
Y
(r)h(j)
beta-NB(h(n + 1); r, S(h) − h(n + 1), c + nr)
B(S(h) − h(n + 1), c + nr)
h(j)!
j=1

= B(S(h), c + (n + 1)r)

n+1
Y
j=1

(r)h(j)
.
h(j)!

In the second product term, note that
Y
m
[digamma(h(n + 1); r, c + nr)] h
(0)

h∈Hn+1

=

 mh
Y 
(r)h(n+1)
−1
λr,c+nr
B(h(n + 1), c + (n + 1)r)
h(n + 1)!
(0)

h∈Hn+1
−(

P

h∈H

= λr,c+nr

(0)
n+1

mh )

Y

(0)

h∈Hn+1

"

B(h(n + 1), c + (n + 1)r)

n+1
Y
j=1

(r)h(j)
h(j)!

#mh

,

where for the last equality, we have used the fact that h(j) = 0 for every j ≤ n and
P
P
P
(0)
h ∈ Hn+1 . Note that ℏ∈Hn mℏ + h∈H(0) mh = h∈Hn+1 mh . Then equation (4.18) is
n+1

equal to

P

(cT )
Q

h∈Hn+1

mh

h∈Hn+1 (mh !)

×

Y

h∈Hn+1

"

exp −cT

n+1
X

[ψ(c + jr) − ψ(c + (j − 1)r)]

j=1

B(S(h), c + (n + 1)r)

n+1
Y
j=1

(r)h(j)
h(j)!

# mh

!

(4.19)
.

Pn+1
Noting that j=1 [ψ(c + jr) − ψ(c + (j − 1)r)] = ψ(c + (n + 1)r) − ψ(c), we obtain the
expression in equation (1.10) for n + 1, as desired.

By construction, equation (1.10) defines the finite-dimensional
marginal distributions
S
of the stochastic process (Mh )h∈H∞ with index set H∞ := n∈N Hn . The exchangeability
result given by Theorem 1.3 then follows from the exchangeability of the sequence X[n] .

5. 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

The combinatorial structure of beta negative binomial processes

17

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.

5.1. Representations as random arrays/matrices
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 Wi,j associated with the ith row and jth column is the multiplicity/mass
of the atom labeled j in the ith point process Xi .
More carefully, fix n ∈ N and let (Mh )h∈Hn be the combinatorial structure of a sequence
X1 , . . . , Xn of conditionally i.i.d. negative binomial processes, given a shared beta process
e0 of
base measure with concentration
parameter c > 0 and non-atomic base measure B
P
finite mass T . Let κ := h∈Hn Mh be the number of unique atoms among X[n] . Then W
is an n × κ array of non-negative integers such that, for every h ∈ Hn , there are exactly
Mh columns of W equal to h, where h is thought of as a length-n column vector. Note
that W will have no columns when κ = 0.
All that remains is to order the columns of W . Every total order on Hn induces a
unique ordering of the columns of W . Titsias [26] defined a unique ordering in this way,
analogous to the left-ordered form defined by Griffiths and Ghahramani [10] for the IBP.
In particular, for h, h′ ∈ Hn , let  denote the lexicographic order given by: h  h′ if and
only if h = h′ or h(η) < h′ (η), where η is the first coordinate where h and h′ differ. We
say W is left-ordered when its columns are ordered according to . Because there is a
bijection between combinatorial structures (Mh )h∈Hn and their unique representations
by left-ordered arrays, the probability mass function of W is given by equation (1.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 [4] for the IBP. Note that the number of
distinct ways of ordering the κ columns is given by the multinomial coefficient
Q

κ!
h∈Hn

Mh !

,

(5.1)

where the denominator arises from the fact that there are Mh indistinguishable columns
for every history h ∈ Hn . The following result is then immediate:
Theorem 5.1. Let W be a uniform random labeling of (Mh )h∈Hn described above, let
n×k
w ∈ Z+
be an array of non-negative integers with n rows and k ≥ 0 non-zero columns,

18
and for every j ≤ k, let sj :=

C. Heaukulani and D.M. Roy
Pn

i=1 wi,j

be the sum of column j. Then

"
#
k
n
Y
Γ(sj )Γ(c + nr) Y (r)wi,j
(cT )k
. (5.2)
exp(−cT [ψ(c + nr) − ψ(c)])
P{W = w} =
k!
Γ(sj + c + nr) i=1 wi,j !
j=1
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 [10] 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
[7] for the IBP, with the parameter r providing flexibility with respect to the counts in
the array. In particular, the total number of features, κ, is Poisson distributed with mean
cT [ψ(c + nr) − ψ(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 κ, column exchangeable.) Finally, note that
the mean of the digamma(r, c) distribution exists for c > 1 and is given by
r
,
(c − 1)(ψ(r + c) − ψ(c))

(5.3)

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.

The combinatorial structure of beta negative binomial processes

19

Figure 1. Simulated Z+ -valued arrays from the NB-IBP. Dots are positive entries, the magnitudes of which determine the size of the dot. The total mass parameter T is varied along the
top row; the concentration parameter c is varied along the middle row; the negative binomial
parameter r is varied along the bottom row. See the text for a summary of how these parameters
affect the expected number of features in total, features per row, and feature multiplicities.

20

C. Heaukulani and D.M. Roy

5.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 y1 , . . . , yn , where each group
yi is comprised of Di measurements yi = (yi,1 , . . . , yi,Di ). In particular, Wi,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. [22] and Titsias [26] 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 Wi,j models the relative
number of times object class j appears in image i. For every i ≤ n, image yi consists of Di
local patches yj,1 , . . . , yj,Di detected in the image, which are (collections of) continuous
variables representing, for example, color, hue, location in the image, etc. Let κ be the
number of columns of W , that is, the number of features. The P
local patches in image i are
modeled as conditionally i.i.d. draws from a mixture of Si = κj=1 Wi,j Gaussian distributions, where Wi,j of these components are associated with feature j. For k = 1, 2, . . . ,
(j,k)

(j,k)

(j,k)

:= (mi , Σi ) denote the mean and covariance of the Gaussian components
let Θi
associated with feature j for image i. Let zi,d = (j, k) when yi,d is assigned to compo(j,k)
nent k ≤ Wi,j associated with feature j ≤ κ. Conditioned on Θ := (Θi )i≤n,j≤κ,k≤Wi,j
and the assignments z := (zi,d )i≤n,d≤Di , the distribution of the measurements admits a
conditional density
p(y|W, Θ, z) =

Di
n Y
Y

z

z

N (yi,d ; mi i,d , Σi i,d ).

(5.4)

i=1 d=1
(j,k)

To share statistical strength across images, the parameters Θi
Bayesian prior:
(j,k)

Θi

i.i.d.

|Θ(j) ∼ ν(Θ(j) )
Θ

(j) i.i.d.

∼ ν0

are given a hierarchical

for every i and k,

for every j.

(5.5)
(5.6)

A typical choice for ν(·) is the family of Gaussian–inverse-Wishart distributions with
feature-specific parameters Θ(j) drawn i.i.d. from a distribution ν0 . Finally, for every image i ≤ n, conditioned on W , the assignment variables zi,1 , . . . , zi,Di for the local patches
in image n are assumed to form a multivariate Pólya urn scheme, arising from repeated
draws from a Dirichlet-distributed probability vector over {(j, k) : j ≤ κ, k ≤ Wi,j }. The
parameters for the Dirichlet distributions are tied in a similar fashion to Θ. The interpretation here is that local patch d in image i is assigned to one of the Si instances of
the latent objects appearing in the image. The number of object instances to which a

The combinatorial structure of beta negative binomial processes

21

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 [5] and Zhou et al. [28] proposed models for text documents in terms of latent
features representing topics. More carefully, let yi,v be the number of occurrences of word
v in document i. Conditioned on W and a collection of non-negative topic-word weights
Θ := (θj,v )j≤κ,v≤V , the word counts are assumed to be conditionally i.i.d. and
yi,v |W, Θ ∼ Poisson

κ
X
j=1

!

Wi,j θj,v .

(5.7)

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 θj,v is large are likely to appear many times. There are a total of
κ topics that are shared across the documents. The topic-word weights Θ 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 [8], where yi,v represents
the rating a Netflix user i assigns to a film v.

5.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 [3, 19, 23, 28].
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 [7, 10].
The models described in Section 5.2 associate every feature with a latent parameter.
Therefore, conditioned on the number of columns κ, let Θ = (θ1 , . . . , θκ ) be an i.i.d. sequence drawn from some non-atomic distribution νΘ , and assume that the data y admits
a conditional density p(y|W, Θ). We will associate the jth column of W with Θj , and so
the pair (W, Θ) 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
p(W, Θ|y) ∝ p(y|W, Θ) × p(W, Θ),

(5.8)

22

C. Heaukulani and D.M. Roy

where p(W, Θ) is a density for the joint distribution of (W, Θ). 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 Wi,j , conditioned on the remaining
elements of the array (collectively denoted by W−(i,j) ), the data y, and the parameters
Θ. By Bayes’ rule, and the independence of Θ and W given κ, we have
P{Wi,j = z|y, W−(i,j) , Θ}
∝ p(y|{Wi,j = z}, W−(i,j) , Θ) × P{Wi,j = z|W−(i,j) }.

(5.9)

Recall that the array W is row-exchangeable, and so, in the language of the NBIBP, we may associate the ith row with the final customer at the buffet. The count
Wi,j is the number of servings the customer takes of dish j, which has been served
P
(−i)
(−i)
> 0, we have
:= i′ 6=i Wi′ ,j times previously. When Sj
Sj
(−i)

Wi,j |W−(i,j) ∼ beta-NB(r, Sj

, c + (n − 1)r).

(5.10)

Therefore, we can simulate from the unnormalized, unbounded discrete distribution in
equation (5.9) using equation (5.10) 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. [17], the second kernel resamples the number, positions, and
(−i)
values of those singleton columns j ′ such that Sj ′ = 0. Simultaneously, we propose a
corresponding change to the sequence of latent parameters Θ, preserving the relative
ordering with the columns of W . This corresponding change to Θ cancels out the effect
of the κ! term appearing in the p.m.f. of the array W . Let Ji be the number of singleton
columns, that is, let
(−i)

Ji = #{j ≤ κ: Wi,j > 0 and Sj

= 0},

(5.11)

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, Ji may be interpreted as the number of
new dishes sampled by the final customer, in which case, we know that
Ji ∼ Poisson(cT [ψ(c + nr) − ψ(c + (n − 1)r)]).

(5.12)

We therefore propose a new array W ∗ by removing the Ji singleton columns from the
array and insert Ji∗ new singleton columns at positions drawn uniformly at random, where
Ji∗ is sampled from the (marginal) distribution of Ji given in equation (5.12). Like those
columns that were removed, each new column has exactly one non-zero entry in the ith
row: We draw each non-zero entry independently and identically from a 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.

The combinatorial structure of beta negative binomial processes

23

Finally, we form a new sequence of latent parameters Θ∗ by removing those entries
from Θ associated with the Ji columns that were removed from W and inserting Ji∗
new entries, drawn i.i.d. from νΘ , at the same locations corresponding to the Ji∗ newly
∗
introduced columns. Let κ∗ := κ − Ji + Ji∗ , and note that there were Jκ∗ possible ways
i
to insert the new columns. Therefore, the proposal density is
−1
κ∗
q(W , Θ |W, Θ) =
Poisson(Ji∗ ; cT [ψ(c + nr) − ψ(c + (n − 1)r)])
Ji∗
(5.13)
Y
Y
∗
νΘ (θ).
×
digamma(Wi,j
; r, c + (n − 1)r)
∗

∗



j≤κ∗

θ∈Θ∗ \Θ

With manipulations similar to those in the proof of Theorem 1.2, it is straightforward
to show that a Metropolis–Hastings kernel accepts a proposal (W ∗ , Θ∗ ) with probability
min{1, α∗ }, where
α∗ =

p(y|W ∗ , Θ∗ )
.
p(y|W, Θ)

(5.14)

Combined with appropriate Metropolis–Hastings moves that shuffle the columns of W
and resample the latent parameters Θ, we obtain a Markov chain whose stationary distribution is the conditional distribution of W and Θ given the data y.
Another benefit of the characterization of the distribution of W in (5.1) 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 [18]. In the particular case when T is given a gamma prior distribution, say
T ∼ gamma(α, β) for some positive parameters α and β, the conditional distribution
again falls into the class of gamma distributions. In particular, the conditional density is
p(T |W, κ) ∝ T α+κ−1 exp(−cT [ψ(c + nr) − ψ(c)] − βT )
∝ gamma(T ; α + κ, β + cT [ψ(c + nr) − ψ(c)]).

(5.15)
(5.16)

Acknowledgements
We thank Mingyuan Zhou for helpful feedback and for pointing out the relation of our
work to that of Sibuya [21]. 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.

24

C. Heaukulani and D.M. Roy

References
[1] Barndorff-Nielsen, O. and Yeo, G.F. (1969). Negative binomial processes. J. Appl.
Probab. 6 633–647. MR0260001
[2] Broderick, T., Jordan, M.I. and Pitman, J. (2013). Cluster and feature modeling from
combinatorial stochastic processes. Statist. Sci. 28 289–312. MR3135534
[3] Broderick, T., Mackey, L., Paisley, J. and Jordan, M.I. (2014). Combinatorial clustering and the beta-negative binomial process. IEEE Trans. Pattern Anal. Mach. Intell.
37 290–306. Special issue on Bayesian nonparametrics.
[4] Broderick, T., Pitman, J. and Jordan, M.I. (2013). Feature allocations, probability
functions, and paintboxes. Bayesian Anal. 8 801–836. MR3150470
[5] Canny, J. (2004). Gap: A factor model for discrete data. In Proceedings of the 27th Annual
International ACM SIGIR Conference on Research and Development in Information
Retrieval, Sheffield, United Kingdom.
[6] Devroye, L. (1992). Random variate generation for the digamma and trigamma distributions. J. Stat. Comput. Simul. 43 197–216. MR1389440
[7] Ghahramani, Z., Griffiths, T.L. and Sollich, P. (2007). Bayesian nonparametric latent
feature models. In Bayesian Statistics 8. Oxford Sci. Publ. 201–226. Oxford: Oxford
Univ. Press. MR2433194
[8] Gopalan, P., Ruiz, F.J.R., Ranganath, R. and Blei, D.M. (2014). Bayesian nonparametric Poisson factorization for recommendation systems. In Proceedings of the 17th
International Conference on Artificial Intelligence and Statistics, Reykjavik, Iceland.
[9] Grégoire, G. (1984). Negative binomial distributions for point processes. Stochastic Process. Appl. 16 179–188. MR0724064
[10] Griffiths, T.L. and Ghahramani, Z. (2006). Infinite latent feature models and the Indian
buffet process. In Advances in Neural Information Processing Systems 19, Vancouver,
Canada.
[11] Hjort, N.L. (1990). Nonparametric Bayes estimators based on beta processes in models
for life history data. Ann. Statist. 18 1259–1294. MR1062708
[12] Kallenberg, O. (2002). Foundations of Modern Probability, 2nd ed. Probability and Its
Applications (New York). New York: Springer. MR1876169
[13] Kim, Y. (1999). Nonparametric Bayesian estimators for counting processes. Ann. Statist.
27 562–588. MR1714717
[14] Kingman, J.F.C. (1967). Completely random measures. Pacific J. Math. 21 59–78.
MR0210185
[15] Kozubowski, T.J. and Podgórski, K. (2009). Distributional properties of the negative
binomial Lévy process. Probab. Math. Statist. 29 43–71. MR2553000
[16] Lo, A.Y. (1982). Bayesian nonparametric statistical inference for Poisson point processes.
Z. Wahrsch. Verw. Gebiete 59 55–66. MR0643788
[17] Meeds, E., Ghahramani, Z., Neal, R.M. and Roweis, S.T. (2007). Modeling dyadic
data with binary latent factors. In Advances in Neural Information Processing Systems
20, Vancouver, Canada.
[18] Neal, R.M. (2003). Slice sampling. Ann. Statist. 31 705–767. With discussions and a
rejoinder by the author. MR1994729
[19] Paisley, J., Zaas, A., Woods, C.W., Ginsburg, G.S. and Carin, L. (2010). A stickbreaking construction of the beta process. In Proceedings of the 27th International
Conference on Machine Learning, Haifa, Israel.

The combinatorial structure of beta negative binomial processes

25

[20] Roy, D.M. (2014). The continuum-of-urns scheme, generalized beta and Indian buffet
processes, and hierarchies thereof. Preprint. Available at arXiv:1501.00208.
[21] Sibuya, M. (1979). Generalized hypergeometric, digamma and trigamma distributions.
Ann. Inst. Statist. Math. 31 373–390. MR0574816
[22] Sudderth, E.B., Torralba, A., Freeman, W.T. and Willsky, A.S. (2005). Describing
visual scenes using transformed Dirichlet processes. In Advances in Neural Information
Processing Systems 18, Vancouver, Canada.
[23] Teh, Y.W., Görür, D. and Ghahramani, Z. (2007). Stick-breaking construction for the
Indian buffet process. In Proceedings of the 11th International Conference on Artificial
Intelligence and Statistics, San Juan, Puerto Rico.
[24] Thibaux, R. and Jordan, M.I. (2007). Hierarchical beta processes and the Indian buffet
process. In Proceedings of the 11th International Conference on Artificial Intelligence
and Statistics, San Juan, Puerto Rico.
[25] Thibaux, R.J. (2008). Nonparametric Bayesian models for machine learning. Ph.D. thesis,
EECS Department, Univ. California, Berkeley. MR2713095
[26] Titsias, M. (2007). The infinite gamma-Poisson feature model. In Advances in Neural
Information Processing Systems 20.
[27] Wolpert, R.L. and Ickstadt, K. (1998). Poisson/gamma random field models for spatial
statistics. Biometrika 85 251–267. MR1649114
[28] Zhou, M., Hannah, L., Dunson, D. and Carin, L. (2012). Beta-negative binomial process
and Poisson factor analysis. In Proceedings of the 29th International Conference on
Machine Learning, Edinburgh, United Kingdom.
[29] Zhou, M., Madrid, O. and Scott, J.G. (2014). Priors for random count matrices derived
from a family of negative binomial processes. Preprint. Available at arXiv:1404.3331v2.
Received June 2014 and revised March 2015

