A sequential algorithm for fast fitting of Dirichlet process mixture models

1.6 MB·English

by  David Nott

A sequential algorithm for fast fitting of Dirichlet process mixture
models
BY DAVID J. NOTT1, XIAOLE ZHANG2, CHRISTOPHER YAU2 & AJAY JASRA1
1DepartmentofStatistics&AppliedProbability,NationalUniversityofSingapore,Singapore,117546,SG.
E-Mail:[email protected], [email protected]
2DepartmentofMathematics,ImperialCollegeLondon,London,SW72AZ,UK.
E-Mail:[email protected], [email protected]
3 Abstract
1
0 In this article we propose an improvement on the sequential updating and greedy search (SUGS) algorithm
2
[20] for fast fitting of Dirichlet process mixture models. The SUGS algorithm provides a means for very fast

n
a approximateBayesianinferenceformixturedatawhichisparticularlyofusewhendatasetsaresolargethatmany
J
standard Markov chain Monte Carlo (MCMC) algorithms cannot be applied efficiently, or take a prohibitively
4
1
long time to converge. In particular, these ideas are used to initially interrogate the data, and to refine models

]
such that one can potentially apply exact data analysis later on. SUGS relies upon sequentially allocating data
O
C to clusters and proceeding with an update of the posterior on the subsequent allocations and parameters which
.
t assumesthisallocationiscorrect. Ourmodificationsoftensthisapproach,byprovidingaprobabilitydistribution
a
t
s over allocations, with a similar computational cost; this approach has an interpretation as a variational Bayes
[
procedureandhencewetermitvariationalSUGS(VSUGS).ItisshowninsimulatedexamplesthatVSUGScan
1
v out-perform, in terms of density estimation and classification, the original SUGS algorithm in many scenarios.
7
9 Inaddition,wepresentadataanalysisforflowcytometrydata,andSNPdataviaathree-classdirichletprocess
8
2 mixture model illustrating the apparent improvement over SUGS.
.
1 Key-words: Approximate Bayesian Inference; Mixture Modelling; Variational Bayes; Density Estimation.
0
3
1
: 1 Introduction
v
i
X
r The demands of fitting models to large data-sets have exploded over the last decade. Increasingly complex data
a
sets are available, which has placed demands on statisticians to develop realistic models to represent these data.
Inevitably, for many classes of models, this places a further emphasis on being able to fit such models accurately
and in a reasonable time-frame.
Inthisarticle,weconsiderfastBayesianstatisticalinferenceforDirichletprocessmixture(DPM)models[1,13].
This particular class of models have proven to be popular in the literature as a tool for both clustering and density
estimation and there are a wide variety of elegant MCMC and sequential Monte Carlo algorithms; see e.g. [15, 19].
Such algorithms provide exact inference from DPM models, but can be very computationally demanding when
trying to analyze extremely large data-sets and even more, exact statistical inference from mixtures is notoriously
difficult; see [11]. As mentioned above, this issue often leads to researchers resorting to approximate inference to
1
browse or interrogate the data, so as to refine model specifications for an exact analysis; one particular important
and interesting method in this direction is the SUGS algorithm.
The SUGS algorithm is a procedure for fast approximate fitting of DPM models. It relies on an approximation
of the exact posterior distribution on the allocation of data to components and the component specific parameters,
by sequentially adding data-points to the model. These data are allocated to a given mixture component and this
allocation is frozen and taken as the truth when updating the posterior for new data points. The method can be
sensitive to the ordering of the data, but [20] provide procedures to select this ordering in a systematic way. In this
paper,wereinterprettheSUGSalgorithmwithinavariationalBayesframework. Thisallowsonetoderivedifferent
approximations of the posterior distribution.
Inparticular,ourinterpretationdoesnotmeanthatoneneedstoallocatedatatoaclusterand,instead,provides
aprobabilitydistributionontheseallocations;thisisdoneataminorincreaseincomputationalcost. Theadvantage
ofthisgeneralization,whichwecallVSUGS,isapparentwhenonefitstheapproximationtodatawhosecomponents
are close in some sense. In this scenario, we have consistently found (and as illustrated in Section 6) that VSUGS
outperformsSUGSinavarietysenses;thisisimportantas,whenthemixturecomponentsarehighlyseparated,such
initial browsing or interrogation of the data is less important. Moreover, our variational approximation provides
a lower-bound on the log-marginal likelihood; we empirically find that this can be used as a technique for model
selection (e.g. selecting the order in which the data arrive) which did not seemingly work well in [20].
This article is structured as follows. We begin with a motivating example in Section 2. In Section 3 we give
a basic summary of DPM models. In Section 4 we discuss SUGS and in Section 5 our generalization VSUGS. In
Section 6 we give numerical examples; both a simulation study and real data analyses associated to flow cytometry
dataandSNPdataviaathree-classdirichletprocessmixturemodel. InSection7weconcludethearticle,discussing
avenues for future work.
2 Motivating Application
A motivating application for us is the problem of genotyping single nucleotide polymorphisms (SNPs) from SNP
genotyping microarray data [9]. Figure 6 illustrates an example dataset. The statistical problem is to characterise
the three genotype classes AA, AB and BB and classify each data point into one of these three classes. This
is a straightforward three-way classification problem that can be approached using hierarchical mixture modelling
where,forexample,theclass-conditionaldensitiesaremodelledusingmultivariateNormalorStudentt-distributions.
Although, these models work well in practice, it is clear that the class-conditional densities are not Normal (or
Student). We can obtain increased accuracy through the use of semi-parametric models for the class-conditional
densities using Dirichlet Process Mixtures. However, the size of the data sets presents a massive challenge for
this type of modelling approach. For a single experiment (individual), modern genotyping microarray produces
2
300,000-5,000,000two-dimensionalmeasurements. Eachstudymayconsistofhundredstothousandsofindividuals.
The data sizes here prohibit the use of Monte Carlo inference and motivate approximate approaches that are able
to scale to the size of problems encountered.
3 Dirichlet process mixture model
Consider a Dirichlet process mixture model of the form, for i∈N:
Y |θ˜ i.∼i.d.P(·|θ˜) θ˜ ∼P P ∼DP(αP ) (1)
i i i i 0
where Yi ∈ Y ⊆ Rdy, θ˜i ∈ Θ ⊆ Rdθ are observation specific parameters, P(·|θ˜i) is a conditional probability which
admits a density p(·|θ˜) w.r.t. a single dominating σ−finite measure for each θ˜ (which is often Lebesgue), P is
i i
an unknown mixing distribution, and P ∼ DP(αP ) indicates that the prior for P is a Dirichlet process [8] with
0
precision parameter α ∈ R and base measure P . P is known and we also consider α to be fixed. We remark
+ 0 0
that, in connection to subsequent methodology to be presented, [20] consider a way of handling unknown α which
can also be used in all the extensions we consider but for simplicity we do not consider this below.
A well known property of the Dirichlet process is that a distribution drawn from it will put all its mass on a
countable set of points. Following the notation of [20] we will write θ ={θ }∞ for the set of distinct values in the
j j=1
sequence θ˜={θ˜ }∞ where points in θ are labelled according to their order of appearance in θ˜. Next, let δ =j if
j j=1 i
θ˜ =θ and write δ ={δ }∞ . Using the P´olya urn characterization of the Dirichlet process [3] we can rewrite (1)
i j j j=1
in the form
i.i.d.
Y |δ ,θ ∼ P(·|θ ) p(δ,θ)=p(δ)p(θ)
i i δi δi
where p(θ) = (cid:81)∞ p (θ ), p (·) is the probability density associated to P and p(δ) is as follows, (writing δ =
j=1 0 j 0 0 1:i
(δ ,...,δ ), with the convention δ is the null vector), for i∈{2,3,...}
1 i 1:0

p(δ =j|δ ) =  α+n(jii−)1 j ∈{1,...,ni} (2)
i 1:i−1
 α+αi−1 j =ni+1
where p(δ = 1) = 1, n(i) = Card({δ : δ = j,1 ≤ m ≤ i−1}), and n = n (δ ) is the maximum value in
1 j m m i i 1:i−1
δ (i.e. the number of components “seen” in the data up to to time i−1). This representation of the model
1:i−1
where the unknown measure P is integrated out is important for many Monte Carlo sampling schemes for fitting
DP mixtures [5, 7, 14, 15].
LaterwewillworkwithatruncatedDirichletprocessmixturemodel,whichisoftenconvenientforcomputations
[10]. SuchtruncationsarebasedonthestickbreakingrepresentationoftheDirichletprocess[18]. Supposewelimit
the number of distinct values appearing in the sequence {θ˜ }∞ to an upper truncation limit T > 1. Then
j j=1
3
generalizing the P´olya urn representation we can consider the truncated Dirichlet process with p(δ,θ) = p(δ)p(θ)
wherenowp(θ)=(cid:81)T p (θ )andp(δ)isdefinedrecursivelybyp(δ =1)=1andsimilarlyto(2)(see,forexample,
j=1 0 j 1
[10, Section 2.2.2]) for i∈{2,3,...}

p(δ =j|δ ) =  n(jαi)++i−α/1T j ∈{1,...,ni} (3)
i 1:i−1
 α(α1−+ni−i/1T) j =ni+1
We note that truncations have also been used in the context of variational approximations for Dirichlet process
mixtures [4] although these authors consider the truncation point as a variational parameter without truncating
the original model. The algorithm of [4], although related to ours, is not a sequential algorithm however - here
we are interested in very fast sequential algorithms related to the SUGS method of [20]. [20] compare their
approach with a variety of other fast computational methodologies for Dirichlet process mixtures, and show that
theiralgorithmiscompetitivewithotherfastapproximationmethodologies. Inthisworkwefocusonlyoncomparing
our new approach with the original SUGS algorithm, and refer the reader to [20] for information about the relative
performance of alternative approximations to SUGS.
4 The SUGS algorithm
SUGSisarecursivealgorithmthattakesattimei−1anestimateδˆ ofδ andanapproximationofp(θ|y )
1:i−1 1:i−1 1:i−1
and produces an estimate δˆ of δ and an approximation of p(θ|y ). To start the recursion we use δˆ = 1 and
1:i 1:i 1:i 1
p(θ|y )=p(θ |y ,δˆ )(cid:81) p (θ ).Hereandinwhatfollowsatermδˆ intheconditioningmeansδ =δˆ andsimilarly
1 1 1 1 j>1 0 j i i i
for δˆ . Suppose at time i−1 we have an estimate δˆ of δ . Consider the posterior distribution p(δ ,θ|y )
1:i 1:i−1 1:i−1 i 1:i
and approximate this by p(δ ,θ|y ,δˆ ). That is to say, we initially consider our approximation of p(δ ,θ|y ),
i 1:i 1:i−1 i 1:i
p (δ ,θ|y ), as
(cid:98)i i 1:i
p (δ ,θ|y ) ∝ p(δ ,θ|y ,δˆ )p(y |θ )
(cid:98)i i 1:i i 1:i−1 1:i−1 i δi
 
= p(δ |δˆ )(cid:89)p(θ |δˆ ,y )p(y |θ )
i 1:i−1 j 1:i−1 1:i−1 i δi
j∈N 
Using this approximation and integrating out θ,
(cid:90)
p (δ |y ) ∝ p(δ |δˆ ) p(θ |δˆ ,y )p(y |θ )dθ .
(cid:98)i i 1:i i 1:i−1 δi 1:i−1 1:i−1 i δi δi
The SUGS algorithm sets:
(cid:104) (cid:105)
δˆ :=argmax p(δ |y ) (4)
i δi∈{1,...,ni(δˆ1:i−1)+1} (cid:98) i 1:i
4
and then one replaces
p(δ |δˆ )p(θ |δˆ ,y )p(y |θ )
i 1:i−1 δi 1:i−1 1:i−1 i δi
by
I (δ )p(θ |δˆ ,y )p(y |θ )
{δˆi} i δi 1:i−1 1:i−1 i δi
to form the approximation:
 
p (δ ,θ|y )=I (δ )(cid:89) p(θ |δˆ ,y )p(θ |δˆ ,y ). (5)
(cid:98)i i 1:i {δˆi} i  j 1:i−1 1:i−1  δˆi 1:i 1:i
j(cid:54)=δˆi
Inthisapproximationthecomponentsθ areindependentfordifferentj,theposteriorforθ forj (cid:54)=δˆ isunchanged
j j i
from time i−1, and the posterior for θ is updated by assuming that δ =δˆ so that y represents an observation
δˆi i i i
from this mixture component. If the mixture components are from the exponential family and conjugate priors
are used, the densities such as p(θ |δˆ ,y ) can be calculated in closed form and sufficient statistics are updated
j 1:i 1:i
recursively, leading to a very efficient update.
[20] consider a number of further innovations in their algorithm. First, since the fitting algorithm is sequential
and there is a dependence of the fit on the ordering of the data, they suggest running their algorithm for different
random orderings and then choosing the best according to a pseudo-likelihood criterion. Secondly, for model
comparison they suggest the approximation
p(y )=p(y |δˆ )
(cid:98) 1:i 1:i 1:i
andshowthatthiscrudeapproximationcanbeusefulfortaskssuchascomparisonofparametricandnonparametric
models. Thirdly, they suggest a way of dealing with unknown α in the Dirichlet process prior.
5 An improvement of the SUGS algorithm
HerewesuggestasimpleimprovementoftheSUGSalgorithmwhichwecallVSUGS(variationalSUGS).Webegin
with a brief introduction to variational Bayes (VB) methods.
5.1 Variational Bayes
Suppose we have a parameter ξ ∈ Ξ ⊆ Rdξ and data y, p(ξ) is the prior density, p(y|ξ) the likelihood and p(ξ|y)
denotestheposteriordensityw.r.t.Lebesguemeasure(whichweuseforpresentationalpurposesonly). InVB[12,2]
we split ξ into blocks ξ =ξ =(ξ ,...,ξ ), ξ ∈Ξ , with Ξ ×···×Ξ =Ξ and seek to find a good approximation
1:k 1 k j j 1 k
5
to p(ξ |y) of the form
1:k
k
(cid:89)
q(ξ )= q(ξ )
(cid:98) 1:k j
j=1
where each q(ξ ) is a probability density w.r.t. the appropriate dimensional Lebesgue measure. Given known
j
probability densities for q(ξ ), j (cid:54)=i, the optimal choice for q(ξ ) for minimizing the Kullback-Leibler divergence
j i
(cid:90) (cid:18) q(ξ ) (cid:19)
KL(q||p):= log (cid:98) 1:k q(ξ )dξ (6)
(cid:98) p(ξ |y) (cid:98) 1:k 1:k
1:k
is
q(ξ ) ∝ exp{E [logp(Ξ)p(y|Ξ)]} (7)
(cid:98) i −q(cid:98)(ξi)
where E [·] denotes expectation w.r.t. (cid:81) q(ξ ). Hence there is a gradient descent algorithm for minimizing
−q(cid:98)(ξi) j(cid:54)=i(cid:98) j
(7)basedonchoosinginitialvalues forthe factors inq(η)andtheniterativelyupdatingeachtermaccordingto(7).
Minimizing (6) is equivalent to maximizing
(cid:90) (cid:18)p(ξ )p(y|ξ )(cid:19)
L(q) := log 1:k 1:k q(ξ )dξ (8)
q(ξ ) (cid:98) 1:k 1:k
(cid:98) 1:k
(cid:82)
and (8) is a lower bound on the log marginal likelihood logp(y) where p(y) = p(ξ )p(y|ξ )dξ . p(y) is a
1:k 1:k 1:k
key quantity in Bayesian model selection and the lower bound is tight, L(q) = logp(y), when q(ξ ) = p(ξ |y).
(cid:98) 1:k 1:k
Generally L(q) is often used as an approximation to logp(y) for model selection in the VB framework.
5.2 The VSUGS algorithm
Aswehaveseenin(5)theSUGSalgorithmrecursivelyapproximatesp(δ ,θ|y ). Consideringthisinavariational
1:i 1:i
framework, suppose we have an approximation to the posterior p(δ ,θ|y ) of the form
1:i−1 1:i−1
(cid:26)i−1 (cid:27)(cid:26) (cid:27)
(cid:89) (cid:89)
q (δ ) q (θ ) .
i−1 j i−1 j
j=1 j∈N
In the above expression, we omit certain conditionings (as will become apparent below) to reduce the subsequent
notational burdens. We will suggest a way to update this approximation of p(δ ,θ|y ), using variational ideas.
1:i 1:i
The approximation will be of the form {(cid:81)i q (δ )}{(cid:81) q (θ )}. We start the recursion with q (δ = 1) = 1,
j=1 i j j∈N i j (cid:98)1 1
q (θ )=p(θ |y ,δ =1), q (θ )=p (θ ), j ∈{2,3,...}.
(cid:98)1 1 1 1 1 (cid:98)1 j 0 j
The idea is to make a particular fixed choice at time i for q (δ ) and not to revisit that choice at future times.
(cid:98)i i
That is, the solution to the (partial) variational optimization at time i−1 is used to initialize the optimization at
timei. TheoriginalSUGSalgorithmchoosesq (δ )=I (δ )whereδˆ isdefinedin(4)butitispossibletomakea
(cid:98)i i {δˆi} i i
6
better choice than this without sacrificing the attractive computational properties of the original SUGS algorithm.
In particular, at time i, we set q (δ )=q (δ ), j ∈{1,...,i−1} and choose
(cid:98)i j (cid:98)i−1 j
(cid:90)
q (δ =j) = q q (θ )p(y |θ )dθ (9)
(cid:98)i i ij (cid:98)i−1 δi i δi δi
for j ∈{1,...,i∧T} where T is a truncation point for the number of mixture components and

 (cid:80)(ki=−11)q(cid:98)i−1(δk=j)+α/T j ∈{1,...,(i−1)∧T}
q = α+i−1
ij
 α(1−[(i−1)∧T]/T) j =(i−1)∧T +1.
α+i−1
q (δ =j)willbechosenasanapproximationtop(δ |y ). Toprovidesomeintuitionforthisselectionofq (δ =j),
(cid:98)i i i 1:i (cid:98)i i
we remark that
p(δ |y ) ∝ p(δ |y )p(y |δ ,y )
i 1:i i 1:i−1 i i 1:i−1
(cid:90)
= p(δ |y ) p(y |θ )p(θ |y )dθ (10)
i 1:i−1 i δi δi 1:i−1 δi
Next, note that
p(δ |y )=E[p(δ |y ,δ )|y ]=E[p(δ |δ )|y ]
i 1:i−1 i 1:i−1 1:i−1 1:i−1 i 1:i−1 1:i−1
where the expectation is w.r.t. p(δ |y ). This suggests approximating p(δ |y ) by taking the expectation
1:i−1 1:i−1 i 1:i−1
in the above expression with respect to the variational posterior q (δ ). Although this approximation is still
(cid:98)i−1 1:i−1
not easy to work with, if we we condition on n =(i−1)∧T in (3) and then take the expectation with respect to
i
q (δ ), we get q . So q is an approximation to p(δ |y ) in (10) and the term
(cid:98)i−1 1:i−1 ij ij i 1:i−1
(cid:90)
q (θ )p(y |δ ,θ )dθ
(cid:98)i−1 δi i i δi δi
in (9) simply approximates the integral in (10) by replacing p(θ |y ) with the variational posterior q (θ ).
δi 1:i−1 (cid:98)i−1 δi
As noted earlier, the original SUGS algorithm can be placed in our framework by using q (δ )=I (δ ), a “hard”
(cid:98)i i {δˆi} i
rather than “soft” allocation to clusters which tends to result in greater under-estimation of uncertainty than in
our approach.
Next, using q (θ ) as the prior for θ at time i for processing the data point y , the optimal choice for q (θ )
(cid:98)i−1 j j i (cid:98)i j
is, via (7)
q (θ ) ∝ q (θ )exp(cid:32)E (cid:32)(cid:88)Ti I (δ )log(cid:0)p(y |Θ )(cid:1)(cid:33)(cid:33)
(cid:98)i j (cid:98)i−1 j −q(cid:98)i(θj) {k} i i k
k=1
∝ q(cid:98)i−1(θj)p(yi|θj)q(cid:98)i(δi=j) (11)
7
where E denotes expectation w.r.t. {(cid:81)i q (δ )}{(cid:81) q (θ )} and T = i ∧ T. Again, with the choice
−q(cid:98)i(θj) h=1(cid:98)i h h(cid:54)=j (cid:98)i h i
q (δ ) = I (δ ), exponential family mixture components and conjugate priors this reduces to the SUGS update.
(cid:98)i i {δˆi} i
If the mixture components are normal then the generalized update above can be done in closed form – we will
give details of this below. Note the attractive form of the above update. The likelihood contribution from the ith
observation is split among different mixture components j according to the weight q (δ =j) rather than assuming
(cid:98)i i
the most likely allocation as in the original SUGS algorithm. We can also use the variational lower bound (8) to
approximate the marginal likelihood more accurately than with log(p(y |δ )) in the SUGS algorithm. Details of
1:i 1:i
this are given in the next section for the case of normal mixture components.
5.3 VSUGS for DP mixtures of normals
[20] consider the case of DP mixtures of normals in detail. In this case, θ = (µ ,ζ ) where µ is the mean for the
i i i i
ith component and ζ is the precision and we have Y |δ = j,θ i.∼i.d. N(µ ,ζ−1). They consider a normal inverse-
i i i δi j j
gamma prior for θ , p (θ ) = p (µ |m,νζ−1)p (ζ |a,b) where p (µ |m,νζ−1) is a normal density with associated
j 0 j 0 j j 0 j 0 j j
distribution with mean ρ and variance νζ−1 and ρ and ν are known hyperparameters, and p (ζ |a,b) is a gamma
j 0 j
density with known parameters a and b.
OurVSUGSalgorithmresultsinq (θ )beingnormalinverse-gammaalso,q (θ )=q (µ |ρ(i),ν(i)ζ−1)q (ζ |a(i),b(i))
(cid:98)i j (cid:98)i j (cid:98)i j j j j (cid:98)i j j j
where q (µ |ρ(i),ν(i)ζ−1) is the normal density with mean ρ(i) and variance ν(i)ζ−1, and q (ζ |a(i),b(i)) is a gamma
(cid:98)i j j j j j j j (cid:98)i j j j
density with parameters a(i) and b(i). The parameters ρ(i), ν(i), a(i) and b(i) are updated recursively by (c.f. [20,
j j j j j j
p. 204])
(cid:110) (cid:111)−1
ν(i) = (ν(i−1))−1+q (δ =j)
j j (cid:98)i i
(cid:26)(cid:110) (cid:111)−1 (cid:27)
ρ(i) = ν(i) ν(i−1) ρ(i−1)+q (δ =j)y
j j j j (cid:98)i i i
q (δ =j)
a(i) = a(i−1)+ (cid:98)i i
j j 2
1 ρ(i)2 ρ(i−1)2
b(i) = b(i−1)+ q (δ =j)y2+ j − j
j j 2(cid:98)i i i ν(i) ν(i−1) 
j j
To calculate the terms q (δ =j) in the VSUGS algorithm, we also need to evaluate the integral
(cid:98)i i
(cid:90)
q (θ )p(y |θ )dθ .
(cid:98)i−1 j i j j
In the normal case with the priors we have chosen this integral evaluates to a t density, t (y ;ρ(i−1),b(i−1)/
2a(i−1) i j j
j
a(i−1)(ν(i−1)+1)) where t (y;m,s2) denotes a t density for y with d degrees of freedom, location parameter m and
j j d
scale parameter s.
An approximate variational lower bound on logp(y ) can also be computed recursively. We can think of the
1:i
8
posterior at stage i−1 as the prior to be updated by the likelihood contribution for the ith observation:
p(δ ,θ|y )∝p(δ ,θ|y )p(y |θ ).
i 1:i i 1:i−1 i δi
Approximating p(δ |y ) by q as we did previously and approximating p(θ|y ) by (cid:81)Ti q (θ ) and calcu-
i 1:i−1 ij 1:i−1 j=1(cid:98)i−1 j
lating the lower bound (8) using these priors for the likelihood contribution p(y |θ ) gives
i δi
T
L(q) = (cid:88)(cid:110)(a(i)−a(i−1))ψ(a(i))−log(cid:0)Γ(a(i))(cid:1)+log(cid:0)Γ(a(i−1))(cid:1)+a(i−1)(cid:16)log(cid:0)b(i)(cid:1)−log(cid:0)b(i−1)(cid:1)(cid:17)
j j j j j j j j
j=1
(b(i−1)−b(i)) (ρ(i)−ρ(i−1))2a(i) 1(cid:32) ν(i) (cid:18) ν(i) (cid:19)(cid:33)(cid:41)
+a(i) j j + j j j + j −1−log j
j b(i) 2ν(i−1) b(i) 2 ν(i−1) ν(i−1)
j j j j j
+(cid:88)T q (δ =j)(cid:40)1ψ(a(i))− 1log(cid:0)b(i)(cid:1)− 1log(cid:0)2π(cid:1)− 1(cid:40)ν(i)+(y −ρ(i))2b(ji)(cid:41)(cid:41)
(cid:98)i i 2 j 2 j 2 2 j j j a(i)
j=1 j
(cid:88)Ti (cid:0) (cid:1) (cid:88)Ti (cid:0) (cid:1)
− q (δ =j)log q (δ =j) + q (δ =j)log q
(cid:98)i i (cid:98)i i (cid:98)i i ij
j=1 j=1
where ψ(·) denotes the digamma function.
6 Numerical examples
6.1 Density Estimation
Setup. We generated 100 data sets for each setting and all results are averaged over that for each data set.
2 3 3
i.i.d.
Y ∼ N(−dµ,0.25)+ N(0,0.5)+ N(dµ,2),i=1,...,N,
i 5 10 10
for 0≤dµ≤5 with N =500. Examples are shown in Figure 2.
WecomparedtherelativeerrorsofSUGStoVSUGSunderdifferentvaluesofα,i.e. α∈{0.1, 0.5, 1, 5, 10...45, 50}.
Throughout T =200 for VSUGS. We used 50 different (but random) orderings of the data and chose the ordering
with the maximal variational lower-bound for VSUGS in Section 5.3 and the best ordering for SUGS as in [20].
We also used a standard Collapsed Gibbs Sampling method [17] for posterior inference on some of the datasets for
comparison. To assess the performance in density estimation we compute the value of
N
e=(cid:88)(fˆ(y )−f(y ))2 (12)
j j
j=1
where fˆ(y ), f(y ), are the estimated (predictive) and true density for the data, evaluated at data-point y and
j j j
examine the relative errors.
9
Results. Figure 4 shows example predictive density estimates from SUGS and VSUGS. Figure 3(a) shows that
for large values of α and closely spaced clusters dµ < 1, VSUGS provides more accurate density estimates than
SUGS.However, fordµ>1andα<20, i.e.well-separatedclusters, thedensityestimatesfromSUGSarerelatively
more accurate.
The computation time for VSUGS is constant for given truncation level T as we use a fixed maximum number
of mixture components. In contrast, the computation time required for SUGS is variable and depends both on the
data set and the order in which the data is processed. Figure 3(b) considers the computation burden for the two
methods. In particular, for large values of α and more mixture components, SUGS can be computationally quite
demanding due to the excessive numbers of mixture components that are realised. Whilst in practice, one might
estimate α, this value is not known and hence SUGS could both be significantly less accurate and computationally
more expensive in many situations.
We compared the SUGS and VSUGS predictive densities with those obtained from Collapsed Gibbs Sampling,
we considered the case dµ=0.2,α=0.1 and show results in Table 1 for different data sizes N and the truncation
parameter T. Using Collapsed Gibbs Sampling as a “gold standard”, we find that VSUGS consistently provides
betterpredictivedensityestimates. ExamplecomputationaltimesforN =500were4secondsforSUGS,12seconds
for VSUGS (T =150) and 193 seconds for Collapsed Gibbs Sampling.
6.2 Density estimation for flow cytometry data
We analyzed the flow cytometry data example, which has been studied thoroughly by [16]. Flow cytometers detect
fluorescent reporter markers that typically correspond to specific cell surface or intracellular proteins on individual
cells, and can assay millions of such cells in a fluid stream in minutes. These data points are associated with one
(or more) componentsof a Gaussian mixture model ([6]) and arefrom human peripheral blood cells, with 6 marker
measurementseach: ForwardScatter(measureofcellsize),SideScatter(Measureofcellgranularity),CD4(marker
for helper T cells), IFNg+IL-2 (effector crytokines), CD8 (a marker for cytotoxic T cells), CD3 (marker for all T
cells); that is, the observations are 6 dimensional (the priors are modified to Normal-inverse Wishart, which leads
to a similar derivation of the VSUGS algorithm as in Section 5.3, in this multivariate scenario). Our objective is to
compare the performance of VSUGS to SUGS and Collapsed Gibbs Sampling for clustering and density estimation
in this multivariate, large data setting.
Data. The size of the whole data is 50,000 with 6 dimensions and [16] state the components of these data
are centered closely. In the following simulations, we adopted a Gamma prior for α, i.e. α ∼ G(1,1) for the three
approaches. Whenconsideringαasunknownweusetheapproachtohandlinguncertaintyinαdescribedin[20]for
all algorithms. The Collapsed Gibbs sampler was run for a 300 iteration burn-in followed by 1000 iterations. This
low number is adopted due to the size and complexity of the data; these type of data scenarios are exactly those
which motivate the development of SUGS and VSUGS algorithms. For the VSUGS approximation, the truncation
10