higher posterior variance, systematic over-estimation of the number of
clusters, and worse computational efficiency. The CRP-Oracle has better
performance, but does clearly worse than p-CRP, and there is still a tendency
for over-estimation. Figure 21 suggests that CRP will have larger probability
on larger cluster numbers especially when the sample size increases, while
pCRP tends to have larger probability on the true cluster number as the sample
size increases. For example, in Sim 1, the probability of selecting three
clusters increases from 0.55 to 0.68 in pCRP when $N$ increases from 300 to
2000 and the probability for all the other cluster number decreases. However,
the probability of finding four clusters stabilizes around 0.37 and 0.38 in
CRP-Oracle when $N$ increases from 300 to 2000. CRP has increased probability
of selecting larger number of clusters (say 5, 6, 7, 8 clusters) when $N$
increases from 300 to 2000. In fact, the proposed pCRP has the largest
concentration probability on the true number of clusters among all the three
methods including CRP-Oracle, and this observation is consistent between
$N=300$ and $N=2000$.

(a) Sim 1, $N$=300

(b) Sim 2, $N$=300

(c) Sim 1, $N$=2000

(d) Sim 2, $N$=2000

Figure 21: Estimated posterior of the number of clusters in observed data for
CRP-Oracle (red x), CRP (blue circle) and pCRP (green star).

Table 5 provides numerical summaries of this simulation. We can see all three
methods lead to similar NMI, but pCRP consistently gives the highest value.
Furthermore, pCRP leads to the lowest value of VI in most tests. The
parsimonious effect of pCRP discussed above is further confirmed by the
average and maximum number of clusters; see the columns $K$ and $K_{\max}$ in
the table.

The posterior density plots in Figure 20 show that there is one small
unnecessary cluster in CRP-Oracle and two small unnecessary clusters in CRP,
while all three methods capture the general shape of the true density thus
provide good fitting performance. The over-clustering effect of CRP is much
reduced by pCRP as seen in Figure 20(c).

$N=300$
---
Method | NMI (SE $\times 10^{-3}$) | VI (SE $\times 10^{-3}$) |  $K$ (SE $\times 10^{-2}$) | $K_{\max}$
Ground truth (Sim 1) | 1.0 | 0.0 | 3 | -
CRP-Oracle (Sim 1) | 0.800 (1.1) | 0.669 (4.5) | 4.2 (2.3) | 8
CRP (Sim 1) | 0.773 (1.2) | 0.795 (5.4) | 5.3 (3.3) | 12
pCRP (Sim 1) | 0.827 (0.7) | 0.580 (4.4) | 3.6 (1.7) | 7
Ground truth (Sim 2) | 1.0 | 0.0 | 2 | -
CRP-Oracle (Sim 2) | 0.211 (1.0) | 1.803 (6.4) | 3.5 (2.7) | 8
CRP (Sim 2) | 0.189 (0.9) | 2.164 (7.8) | 6.2 (4.2) | 13
pCRP (Sim 2) | 0.228 (1.0) | 1.518 (2.1) | 2.4 (1.3) | 6
$N=2000$
Method | NMI (SE $\times 10^{-4}$) | VI (SE $\times 10^{-3}$) |  $K$ (SE $\times 10^{-2}$) | $K_{\max}$
Ground truth (Sim 1) | 1.0 | 0.0 | 3 | -
CRP-Oracle (Sim 1) | 0.812 (5.3) | 0.610 (2.6) | 4.0 (2.3) | 10
CRP (Sim 1) | 0.782 (8.5) | 0.732 (4.0) | 5.8 (3.6) | 12
pCRP (Sim 1) | 0.823 (6.6) | 0.869 (7.3) | 3.5 (1.6) | 7
Ground truth (Sim 2) | 1.0 | 0.0 | 2 | -
CRP-Oracle (Sim 2) | 0.258 (7.0) | 1.537 (5.3) | 3.5 (2.6) | 8
CRP (Sim 2) | 0.238 (7.1) | 1.755 (6.8) | 6.8 (4.7) | 15
pCRP (Sim 2) | 0.258 (4.2) | 1.368 (0.7) | 2.2 (1.0) | 5
Table 5: Comparison of CRP and pCRP on Sim 1 and Sim 2. $K$ is the average
number of found clusters. $K_{\max}$ is the maximum number of clusters during
sampling. SE is the standard error of mean. Ground truth is calculated using
the true assignments.

(a) True clustering when $N$=3000

(b) CRP-Oracle when $N$=3000

(c) CRP when $N=3000$

(d) pCRP when $N$=3000

Figure 22: Results of clustering 3000 randomly sampled digits from 1 to 4 in
spectral space. Observations in the same color represent the same digit. CRP-
Oracle and CRP seem to over-fit the noise (the red cluster). $N=1000$
---
Method | NMI (SE $\times 10^{-4}$) | VI (SE $\times 10^{-3}$) |  $K$ (SE $\times 10^{-2}$) | $K_{\max}$
Ground truth | 1.0 | 0 | 4 | -
CRP-Oracle | 0.651 (3.3) | 1.382 (1.4) | 4.37 (1.3) | 7
CRP | 0.651 (3.3) | 1.386 (1.4) | 4.58 (1.6) | 8
pCRP | 0.651 (3.3) | 1.382 (1.4) | 4.08 (0.6) | 6
$N=3000$
Method | NMI (SE $\times 10^{-4}$) | VI (SE $\times 10^{-3}$) |  $K$ (SE $\times 10^{-2}$) | $K_{\max}$
Ground truth | 1.0 | 0.0 | 4 | -
CRP-Oracle | 0.651 (2.0) | 1.400 (1.1) | 5.17 (1.2) | 8
CRP | 0.651 (2.0) | 1.402 (1.1) | 5.44 (1.6) | 9
pCRP | 0.652 (1.9) | 1.389 (1.1) | 4.57 (1.2) | 7
Table 6: Comparison of CRP-Oracle, CRP and pCRP on a 4 digits subset of MNIST.
$K$ is the average number of found clusters. $K_{\max}$ is the maximum number
of clusters during sampling.

#####  11.10.8 Digits 1-4

In this experiment, we cluster 1000 and 3000 digits of the classes 1 to 4 in
MNIST data set (LeCun et al., 2010), where the four clusters are approximate
equally distributed. From cross validation on a different set of 1000 samples,
we obtain the power value $r=1.05$. The concentration parameter $\alpha$ in
CRP-Oracle is calculated as 0.58 ($N=1000$) and 0.5 $(N=3000)$.

Figure 22 shows the clustering result of all the three methods for $N=3000$.
Both CRP and CRP-Oracle seem to over-fit the data by introducing a small
cluster (in red), while pCRP gives a cleaner clustering result with four
clusters. This comparison is further confirmed by Table 6, where the average
posterior cluster number in CRP apparently increases when $N$ grows to 3000.
In contrast, pCRP is closer to the true situation by reducing the over-
clustering effect, even compared to CRP-Oracle; see the columns of $K$ and
$K_{\max}$. All methods lead to similar NMI but pCRP gives lower VI.

(a) The distribution of Old Faithful Geyser after standardization

(b) Clustering results using CRP-Oracle, CRP, pCRP and manual clustering.

Figure 23: Clustering result for Old Faithful Geyser

#####  11.10.9 Old Faithful Geyser

The Old Faithful Geyser data ($N=272$) are widely used to illustrate the
performance of clustering algorithms. We use a test sample of 100 in CV
leading to the power value $r=1.11$. We compare all methods on the other 172
data points. A manual clustering that consists of two Gaussian components is
viewed as the ground truth. The concentration parameter is 0.39 in CRP-Oracle.
Figure 23(b) shows the size of each component obtained from all methods and
the manual clustering. We can see that there are two mixture components in
CRP-Oracle and pCRP, and four mixture components in the CRP method. In this
case where the sample size is relatively small, we again see that pCRP
successfully suppresses small components and generate results closer to the
ground truth than CRP.

####  11.11 Some issues in finite and infinite Gaussian mixture model

#####  11.11.1 Non-identifiability due to overfitting finite mixture model or
infinite mixture model

In finite Gaussian mixture model, when the number of components is unknown,
the analyst can intentionally or unintentionally propose an overfitting model,
i.e., one with more components than the true component the data come from. And
in infinite Gaussian mixture model, the model does not assume the upper bound
of the number of components. Thus, the problem of non-identifiability in
estimation of overfitted mixture model and infinite mixture model is well
known. For example, (Frühwirth-Schnatter, 2006) observed that identifiability
will be violated as either one of the component weights is 0 or two of the
component parameters are equal.

More precisely, as in Section 10.1, assume we have $N$ observations
$\mathcal{X}=\\{\bm{x}_{1},\bm{x}_{2},\ldots,\bm{x}_{N}\\}$ sampled i.i.d.,
from a finite mixture distribution with density

$f_{0}(\bm{x}|\bm{\pi},\boldsymbol{\gamma})=\sum_{k=1}^{K_{0}}\pi_{k}^{0}\Phi(\bm{x}|\boldsymbol{\gamma}_{k}^{0}),$
(144)

with $\boldsymbol{\gamma}_{k}\in\Gamma$, $\Gamma$ is the metric space of the
parameter for some kernel function, $\Phi$ is the kernel function of each
component and $K_{0}$ is the true component number both in finite and finite
Gaussian mixture model. We wish to make Bayesian inference for the model
parameters $\boldsymbol{\theta}=\\{\bm{\pi}^{0},\boldsymbol{\gamma}^{0}\\}$.
In such cases the model is non-identifiable since all values of the parameter
in the form

$\boldsymbol{\theta}=\\{\pi_{1}^{0},\ldots,\pi_{K_{0}}^{0},0;\boldsymbol{\gamma}_{1}^{0},\ldots,\boldsymbol{\gamma}_{K_{0}}^{0},\boldsymbol{\gamma}\\},$
(145)

for all $\boldsymbol{\gamma}\in\Gamma$, and all values of the parameter in the
form
$\boldsymbol{\theta}=\\{\pi_{1}^{0},\ldots,\pi_{j},\ldots,\pi_{K_{0}}^{0},\pi_{K_{0}+1};$
$\boldsymbol{\gamma}_{1}^{0},\ldots,\boldsymbol{\gamma}_{K_{0}}^{0},\boldsymbol{\gamma}_{j}^{0}\\}$
with $\pi_{j}+\pi_{K_{0}+1}=\pi_{j}^{0}$ satisfy
$f_{0}=f_{\boldsymbol{\theta}}$. As stated in (Rousseau and Mengersen, 2011),
this non-identifiability is much more tough to deal with than the non-
identifiability corresponding to permutation of the labels in the mixture
representation. Interesting readers can refer to (Rousseau and Mengersen,
2011) for more details and the references there in.

#####  11.11.2 Read on

If you are interested in a more through understanding of the hyperprior in
mixture model, (Rasmussen, 1999) gives some ideas how we can put hyperprior on
NIW such that release the burden to select hype-parameter $\boldsymbol{\beta}$
for NIW.

####  11.12 Pruning methods for Dirichlet process mixture model

Further to power Chinese restaurant process, we introduce some pruning methods
for Dirichlet proces mixture model.

As a recap, some problem of DP mixture models have been brought about when
apply them to practical problems. For example, they always produce more
components than that the real data should have. The small mixture components
are mainly caused by noise. Some approaches have been proposed for solving
this problem. In (McCullagh et al., 2008), an upper bound of the number of
components is fixed in advance to limit the number in modeling. In (Escobar
and West, 1995), components with little data points are simply discarded, and
these data points are reassigned to other existing components. However, this
two methods based on simple upper bound or thresholds can not be directly used
for real world data, because when you choose a larger bound, DP mixture models
can still results in small clusters, and choosing the best thresholds is
usually difficult. In this thread, we focus on how to shrink small clusters
during sampling.

#####  11.12.1 Constrained sampling (cSampling)

During Gibbs sampling, remove small number of clusters every $s$ (e.g., 20)
iterations, reassign these data to big clusters by the probability of each
cluster. Further extension on cSampling is that when we reassign these small
clusters, we can use the assignment method similar to K-means. However, the
basic idea is the same.

Important points:

  * •

Different to the method in (Escobar and West, 1995), our cSampling method does
not need to choose best threshold to discard unuseful clusters. A small
threshold is good. It aims to remove very small clusters that can be easily
recognized as “noise”.

  * •

When we remove small cluster during sampling, it will have influence on later
sampling iterations. That is where the name “constrained” come from.

#####  11.12.2 loss based sampling (lSampling)

In (Kulis and Jordan, 2011), the authors introduce an algorithm called DP-
means. We here briefly review DP-means. The authors considered asymptotic
behavior of DP mixture models, obtaining a hard clustering algorithm that
behaves similarly to K-means with the exception that a new cluster is formed
whenever a point is farther than $distance$ away from existing cluster
centroid. The $distance$ is very hard to decide, the authors used cross-
validation to decide it. However, this violates the setting of unsupervised
learning.

Inspired by DP-means, when we doing sampling during DP mixture models, we can
shrink out small cluster by some metric, for example: marginal of data and
component assignment $p(\mathcal{X},\bm{z}|\alpha,\boldsymbol{\beta})$, where
$\alpha$ is the concentration parameter on Dirichlet Process,
$\boldsymbol{\beta}$ is the prior parameter on kernel (see Section 12.1). And
also we can use inertia, or so called the within-cluster sum of squares
criterion (see Section 12.3). Again, during Gibbs sampling, remove small
number of clusters every $s$ (e.g., 20) iterations by applying to the chosen
metric, i.e., if removing the small clusters will get smaller loss, we remove,
otherwise, we keep them as they are. In our proposal, we recommend to use the
following loss function

$\sum_{k=1}^{K}\sqrt{\sum_{j:j\in
C_{k}}^{N_{k}}||\bm{x}_{j}-\overline{\bm{x}}_{k}||^{2}},$ (146)

where $C_{k}$ is the data samples in the $k^{th}$ cluster. The reason we use a
square root over each cluster is that it can overcome identifiablity issue.
One can imagine that if two clusters have same center value, the square root
operation will force the two clusters into one cluster. This idea on the loss
function comes from (Petralia et al., 2012), in which case they put a
“repulsive” prior on the mixture components, thus overcomes identifiability
issue in some sense.

#####  11.12.3 Posterior inference

Although the proposed cSampling and lSampling is generic, we focus on its
application in Gaussian mixture models for concreteness. Again, we develop a
collapsed Gibbs sampling algorithm (Neal, 2000) for posterior computation.

Let $\mathcal{X}$ be the observations, assumed to follow a mixture of
multivariate Gaussian distributions. We use a conjugate normal-inverse-Wishart
(NIW) prior $p(\boldsymbol{\mu},\boldsymbol{\Sigma}|\boldsymbol{\beta})$ for
the mean vector $\boldsymbol{\mu}$ and covariance matrix $\boldsymbol{\Sigma}$
in each multivariate Gaussian component, where $\boldsymbol{\beta}$ consists
of all the hyperparameters in NIW. A key quantity in a collapsed Gibbs sampler
is the probability of each customer $i$ sitting with table $k$:
$p(z_{i}=k|\mathbf{z}_{-i},\mathcal{X},\alpha,\boldsymbol{\beta})$, where
$\mathbf{z}_{-i}$ are the seating assignments of all the other customers and
$\alpha$ is the concentration parameter in CRP. This probability is calculated
as follows:

$\displaystyle
p(z_{i}=k|\mathbf{z}_{-i},\mathcal{X},\alpha,\boldsymbol{\beta})$
$\displaystyle\varpropto
p(z_{i}=k|\mathbf{z}_{-i},\alpha,\cancel{\boldsymbol{\beta}})p(\mathcal{X}|z_{i}=k,\mathbf{z}_{-i},\cancel{\alpha},\boldsymbol{\beta})$
(147)
$\displaystyle=p(z_{i}=k|\mathbf{z}_{-i},\alpha)p(\mathbf{x}_{i}|\mathcal{X}_{-i},z_{i}=k,\mathbf{z}_{-i},\boldsymbol{\beta})p(\mathcal{X}_{-i}|\cancel{z_{i}=k},\mathbf{z}_{-i},\boldsymbol{\beta})$
$\displaystyle\varpropto
p(z_{i}=k|\mathbf{z}_{-i},\alpha)p(\mathbf{x}_{i}|\mathcal{X}_{-i},z_{i}=k,\mathbf{z}_{-i},\boldsymbol{\beta})$
$\displaystyle\varpropto
p(z_{i}=k|\mathbf{z}_{-i},\alpha)p(\mathbf{x}_{i}|\mathcal{X}_{k,-i},\boldsymbol{\beta}),$

where $\mathcal{X}_{k,-i}$ are the observations in table $k$ excluding the
$i^{th}$ observation. Algorithm 6 and 7 give the pseudo code of the collapsed
Gibbs sampler to implement cSampling and lSampling in Gaussian mixture models.

input : Choose an initial $\bm{z}$, set constrained step $s$,
threshold=$thres$;

for  _$t\leftarrow 1$ to $T$ iterations_ do

for  _$i\leftarrow 1$ to $N$_ do

Remove $\mathbf{x}_{i}$’s statistics from component $z_{i}$ ;

for  _$k\leftarrow 1$ to $K$_ do

Calculate
$p(z_{i}=k|\mathbf{z}_{-i},\alpha)=\frac{\textit{N}_{k,-i}}{N+\alpha-1}$;

Calculate $p(\mathbf{x}_{i}|\mathcal{X}_{k,-i},\boldsymbol{\beta})$;

Calculate
$p(z_{i}=k|\mathbf{z}_{-i},\mathcal{X},\alpha,\boldsymbol{\beta})\propto
p(z_{i}=k|\mathbf{z}_{-i},\alpha)p(\mathbf{x}_{i}|\mathcal{X}_{k,-i},\boldsymbol{\beta})$;

end for

Calculate
$p(z_{i}=k^{\star}|\mathbf{z}_{-i},\alpha)=\frac{\alpha}{N+\alpha-1}$;

Calculate $p(\mathbf{x}_{i}|\boldsymbol{\beta})$;

Calculate
$p(z_{i}=k^{\star}|\mathbf{z}_{-i},\mathcal{X},\alpha,\boldsymbol{\beta})\propto
p(z_{i}=k^{\star}|\mathbf{z}_{-i},\alpha)p(\mathbf{x}_{i}|\boldsymbol{\beta})$;

Sample $k_{new}$ from
$p(z_{i}|\mathbf{z}_{-i},\mathcal{X},\alpha,\boldsymbol{\beta})$ after
normalizing;

Add $\mathbf{x}_{i}$’s statistics to the component $z_{i}=k_{new}$ ;

If any component is empty, remove it and decrease $K$.

end for

if  _t == s_ then // constrain step

Get cluster index $k$ where cluster number $n_{k}>thres$, put these indexes
into set $\mathcal{C}$;

for  _$i\leftarrow 1$ to $N$_ do

Remove $\mathbf{x}_{i}$’s statistics from component $z_{i}$ ;

for  _$k$ in $\mathcal{C}$_ do

Calculate
$p(z_{i}=k|\mathbf{z}_{-i},\alpha)=\frac{\textit{N}_{k,-i}}{N+\alpha-1}$;

Calculate $p(\mathbf{x}_{i}|\mathcal{X}_{k,-i},\boldsymbol{\beta})$;

Calculate
$p(z_{i}=k|\mathbf{z}_{-i},\mathcal{X},\alpha,\boldsymbol{\beta})\propto
p(z_{i}=k|\mathbf{z}_{-i},\alpha)p(\mathbf{x}_{i}|\mathcal{X}_{k,-i},\boldsymbol{\beta})$;

end for

Sample $k_{new}$ from
$p(z_{i}|\mathbf{z}_{-i},\mathcal{X},\alpha,\boldsymbol{\beta})$ after
normalizing;

Add $\mathbf{x}_{i}$’s statistics to the component $z_{i}=k_{new}$ ;

If any component is empty, remove it and decrease $K$.

end for

end if

end for

Algorithm 6 Collapsed Gibbs sampler for constrained sampling: reassign these
data to big clusters by the probability of each cluster.

input : Choose an initial $\bm{z}$, set loss-based step $s$;

for  _$t\leftarrow 1$ to $T$ iterations_ do

for  _$i\leftarrow 1$ to $N$_ do

Remove $\mathbf{x}_{i}$’s statistics from component $z_{i}$ ;

for  _$k\leftarrow 1$ to $K$_ do

Calculate
$p(z_{i}=k|\mathbf{z}_{-i},\alpha)=\frac{\textit{N}_{k,-i}}{N+\alpha-1}$;

Calculate $p(\mathbf{x}_{i}|\mathcal{X}_{k,-i},\boldsymbol{\beta})$;

Calculate
$p(z_{i}=k|\mathbf{z}_{-i},\mathcal{X},\alpha,\boldsymbol{\beta})\propto
p(z_{i}=k|\mathbf{z}_{-i},\alpha)p(\mathbf{x}_{i}|\mathcal{X}_{k,-i},\boldsymbol{\beta})$;

end for

Calculate
$p(z_{i}=k^{\star}|\mathbf{z}_{-i},\alpha)=\frac{\alpha}{N+\alpha-1}$;

Calculate $p(\mathbf{x}_{i}|\boldsymbol{\beta})$;

Calculate
$p(z_{i}=k^{\star}|\mathbf{z}_{-i},\mathcal{X},\alpha,\boldsymbol{\beta})\propto
p(z_{i}=k^{\star}|\mathbf{z}_{-i},\alpha)p(\mathbf{x}_{i}|\boldsymbol{\beta})$;

Sample $k_{new}$ from
$p(z_{i}|\mathbf{z}_{-i},\mathcal{X},\alpha,\boldsymbol{\beta})$ after
normalizing;

Add $\mathbf{x}_{i}$’s statistics to the component $z_{i}=k_{new}$ ;

If any component is empty, remove it and decrease $K$.

end for

if  _t == s_ then // loss-based step

while  _$K\geq 2$_ do

Get cluster index $k$ where cluster number $n_{k}$ is minimal, put these
indexes which are not equal to $k$ into set $\mathcal{C}$;

for  _$i\leftarrow 1$ to $N$_ do

Remove $\mathbf{x}_{i}$’s statistics from component $z_{i}$ ;

for  _$k$ in $\mathcal{C}$_ do

Calculate
$p(z_{i}=k|\mathbf{z}_{-i},\alpha)=\frac{\textit{N}_{k,-i}}{N+\alpha-1}$;

Calculate $p(\mathbf{x}_{i}|\mathcal{X}_{k,-i},\boldsymbol{\beta})$;

Calculate
$p(z_{i}=k|\mathbf{z}_{-i},\mathcal{X},\alpha,\boldsymbol{\beta})\propto
p(z_{i}=k|\mathbf{z}_{-i},\alpha)p(\mathbf{x}_{i}|\mathcal{X}_{k,-i},\boldsymbol{\beta})$;

end for

Sample $k_{new}$ from
$p(z_{i}|\mathbf{z}_{-i},\mathcal{X},\alpha,\boldsymbol{\beta})$ after
normalizing;

Add $\mathbf{x}_{i}$’s statistics to the component $z_{i}=k_{new}$ ;

If any component is empty, remove it and decrease $K$.

end for

Calculate current loss $l$ ;

end while

Roll back to the status with minimal loss;

end if

end for

Algorithm 7 Collapsed Gibbs sampler for loss-based sampling: remove small
clusters if it results in smaller loss.

#####  11.12.4 Examples

We conduct experiments to demonstrate the main advantages of the proposed
pruning sampling methods using both synthetic and real data.

In all experiments, we run the Gibbs sampler 20,000 iterations with a burn-in
of 10,000. The sampler is thinned by keeping every 5th draw. We use the same
concentration parameter $\alpha=1$ for both CRP and pruning sampling methods
in all scenarios. In addition, we equip CRP with an unfair advantage to match
the magnitude of its prior mean $\alpha\log(N)$ to the true number of
clusters, termed CRP-Oracle. In order to measure overall clustering
performance, we use normalized mutual information (NMI) (McDaid et al., 2013)
and variation of information (VI) (Meilă, 2003), which measures the similarity
between the true and estimated cluster assignments. Higher NMI and lower VI
indicate better performance. If applicable, metrics using the true clustering
are calculated to provide an upper bound for all methods, coded as ‘Ground
Truth’. The metrics are discussed in Section 12. Feel free to skip this
section for a first reading.

#####  11.12.5 Simulation experiments

The parameters of the simulations are as follows:

  * •

Sim 1: $K_{0}=3$, $\bm{\pi}$={0.35, 0.4, 0.25}, $\bm{\mu}$={0, 2, 5} and
$\bm{\Sigma}$={0.5, 0.5, 1};

  * •

Sim 2: $K_{0}=2$, $\bm{\pi}$={0.65, 0.35}, $\bm{\mu}$={1, 1} and
$\bm{\Sigma}$={10, 1};

In practice, we think that a trivial threshold for cSampling is $4\%$, which
means that we consider cluster with samples smaller than $4\%$ of total
samples can be regarded as noise. For lSampling, we shrink by the proposed
loss function every 20 steps. Table 7 shows the posterior summary for these
tests. We notice that cSampling gives better results than CRP-Oracle but a
little bit worse than pCRP. lSampling gives worse results than CRP-Oracle.

$N=300$
---
Method | NMI (SE $10^{-3}$) | VI (SE $10^{-3}$) |  $K$ (SE $10^{-2}$) | $K_{\max}$ | $K_{mode}$
Ground truth (Sim 1) | 1.0 | 0.0 | 3 | - | -
CRP-Oracle (Sim 1) | 0.800 (1.1) | 0.669 (4.5) | 4.2 (2.3) | 8 | 4
CRP (Sim 1) | 0.773 (1.2) | 0.795 (5.4) | 5.3 (3.3) | 12 | 5
pCRP (Sim 1) | 0.827 (0.74) | 0.580 (4.4) | 3.6 (1.7) | 7 | 3
cSampling (Sim 1) | 0.829 (0.75) | 0.695 (6.9) | 3.3 (1.3) | 7 | 3
lSampling (Sim 1) | 0.791 (1.6) | 0.682 (5.2) | 4.3 (3.0) | 11 | 3
Ground truth (Sim 2) | 1.0 | 0.0 | 2 | - | -
CRP-Oracle (Sim 2) | 0.211 (1.0) | 1.803 (6.4) | 3.5 (2.7) | 8 | 3
CRP (Sim 2) | 0.189 (0.9) | 2.164 (7.8) | 6.2 (4.2) | 13 | 6
pCRP (Sim 2) | 0.228 (1.0) | 1.518 (2.1) | 2.4 (1.3) | 6 | 2
cSampling (Sim 2) | 0.231 (1.0) | 1.526 (2.2) | 2.5 (2.0) | 7 | 2
lSampling (Sim 2) | 0.218 (1.1) | 1.707 (5.7) | 4.1 (4.3) | 13 | 2
$N=2000$
Method | NMI (SE $10^{-4}$) | VI (SE $10^{-3}$) |  $K$ (SE $10^{-2}$) | $K_{\max}$ | $K_{mode}$
Ground truth (Sim 1) | 1.0 | 0.0 | 3 | - | -
CRP-Oracle (Sim 1) | 0.812 (5.3) | 0.610 (2.6) | 4.0 (2.3) | 10 | 4
CRP (Sim 1) | 0.782 (8.5) | 0.732 (4.0) | 5.8 (3.6) | 12 | 5
pCRP (Sim 1) | 0.823 (6.6) | 0.869 (7.3) | 3.5 (1.6) | 7 | 3
cSampling (Sim 1) | 0.825 (2.6) | 0.552 (3.1) | 3.3 (1.4) | 7 | 3
lSampling (Sim 1) | 0.815 (5.5) | 0.580 (1.7) | 4.4 (3.0) | 10 | 3
Ground truth (Sim 2) | 1.0 | 0.0 | 2 | - | -
CRP-Oracle (Sim 2) | 0.258 (7.0) | 1.537 (5.3) | 3.5 (2.6) | 8 | 3
CRP (Sim 2) | 0.238 (7.1) | 1.755 (6.8) | 6.8 (4.7) | 15 | 6
pCRP (Sim 2) | 0.258 (4.2) | 1.368 (0.74) | 2.2 (1.0) | 5 | 2
cSampling (Sim 2) | 0.286 (4.0) | 1.351 (0.81) | 2.5 (2.0) | 7 | 2
lSampling (Sim 2) | 0.278 (4.4) | 1.396 (1.6) | 4.1 (4.2) | 11 | 2
Table 7: Posterior summary for pruning methods. $K_{\max}$ is the maximum
number of clusters during sampling. $K_{mode}$ is the most probability of
cluster number during sampling.

#####  11.12.6 Digits 1-4

In this experiment, we cluster 1000 digits of the classes 1 to 4 in MNIST data
set (LeCun et al., 2010), where the four clusters are approximate equally
distributed. The concentration parameter $\alpha$ in CRP-Oracle is calculated
as 0.5 ($N=3000$).

Both CRP and CRP-Oracle seem to over-fit the data by introducing a small
cluster, while cSampling gives a cleaner clustering result with four clusters.
This comparison is further confirmed by Table 8. In contrast, cSampling and
lSampling are closer to the true situation by reducing the over-clustering
effect, even compared to CRP-Oracle; see the columns of $K$ and $K_{\max}$.
All methods lead to similar NMI. We also observe similar result as in
simulation test that lSampling gives worse results than cSampling.

(a) True clustering when $N$=3000

(b) cSampling when $N$=3000

(c) lSampling $N$=3000

Figure 24: Results of clustering 3000 randomly sampled digits from 1 to 4 in
spectral space. Observations in the same color represent the same digit.
$N=3000$
---
Method | NMI (SE $\times 10^{-4}$) | VI (SE $\times 10^{-3}$) |  $K$ (SE $\times 10^{-2}$) | $K_{\max}$
Ground truth | 1.0 | 0.0 | 4 | -
CRP-Oracle | 0.651 (2.0) | 1.400 (1.1) | 5.17 (1.2) | 8
CRP | 0.651 (2.0) | 1.402 (1.1) | 5.44 (1.6) | 9
pCRP | 0.652 (1.9) | 1.389 (1.1) | 4.57 (1.2) | 7
cSampling | 0.659 (1.7) | 1.353 (0.7) | 4.25 (1.2) | 7
lSampling | 0.653 (3.6) | 1.379 (1.2) | 4.81 (1.9) | 9
Table 8: Comparison of CRP-Oracle, CRP and pCRP on a 4 digits subset of MNIST.
$K$ is the average number of found clusters. $K_{\max}$ is the maximum number
of clusters during sampling.

###  12 Some metrics

In order to evaluate the Gibbs sampling procedure and to ensure that mixing is
taking place, it is useful to have some metrics to calculate over the sampling
iterations. We consider two kinds of metrics, one is label-related, in which
case we use the true label of clustering to evaluate the process; the other
one is non-label-related, in which case we do not use the true label to
evaluate the process.

####  12.1 Marginal of data and component assignment

#####  12.1.1 In Bayesian finite Gaussian mixture model

Marginal of the data and component assignments
$p(\mathcal{X},\bm{z}|\boldsymbol{\alpha},\boldsymbol{\beta})$ is useful for
evaluating the Gibbs sampling process since it captures both changes in the
likelihood of the data under the current assignments through
$p(\mathcal{X}|\bm{z},\boldsymbol{\beta})$, as well as the probability of the
current component assignments $p(\bm{z}|\boldsymbol{\alpha})$. This marginal
of data and component assignments can be calculated as follows

$\displaystyle p(\mathcal{X},\bm{z}|\boldsymbol{\alpha},\boldsymbol{\beta})$
$\displaystyle=p(\mathcal{X}|\bm{z},\boldsymbol{\beta})p(\bm{z}|\boldsymbol{\alpha})$
(148)
$\displaystyle=\left(\prod_{k=1}^{K}p(\mathcal{X}_{k}|\boldsymbol{\beta})\right)p(\bm{z}|\boldsymbol{\alpha}),$

where $\mathcal{X}_{k}$ is the set of data observations assigned to
component/cluster $k$. The terms in the product in Equation (148) can each be
calculated using Equation (38) and Equation (77).

#####  12.1.2 In Bayesian infinite Gaussian mixture model

Similar to the finite case, the marginal of the data and component assignments
$p(\mathcal{X},\bm{z}|\alpha,\boldsymbol{\beta})$ can be used as evaluation of
the sampling process

$\displaystyle p(\mathcal{X},\bm{z}|\alpha,\boldsymbol{\beta})$
$\displaystyle=p(\mathcal{X}|\bm{z},\boldsymbol{\beta})p(\bm{z}|\alpha)$ (149)
$\displaystyle=\left(\prod_{k=1}^{K}p(\mathcal{X}_{k}|\boldsymbol{\beta})\right)p(\bm{z}|\alpha),$

The terms in the product in Equation (149) can each be calculated using
Equation (38) and Equation (119). The only difference between Equation (149)
and Equation (148) is that in Equation (148) the marginal probability of
assignments depends on a vector $\boldsymbol{\alpha}$, while in Equation
(149), it depends on a scalar $\alpha$. Also in infinite case, the cluster
number $K$ can increase or decrease at each iteration.

####  12.2 Mixture likelihood

The marginal likelihood can be used as a metric to evaluate the sampling
iterations since it captures both the likelihood of data under current
assignment through $p(\mathcal{X}|\bm{z},\boldsymbol{\beta})$ , and the
probability of the current component assignment
$p(\bm{z}|\boldsymbol{\alpha})$. In this sense, an alternative metric can be
utilized from Equation (60) after we sample out the distribution parameters
$\boldsymbol{\mu}_{k}$’s and $\boldsymbol{\Sigma}_{k}$’s for multivariate
Gaussian distributions:

$\displaystyle p(\bm{\pi},\boldsymbol{\gamma}|\mathcal{X},\bm{z})$
$\displaystyle=\prod_{i=1}^{N}\pi_{z_{i}}\mathcal{N}(\bm{x}_{i}|\boldsymbol{\gamma}_{z_{i}})$
$\displaystyle=\prod_{k=1}^{K}\pi_{k}^{N_{k}}\left[\prod_{i:i\in
C_{k}}\mathcal{N}(\bm{x}_{i}|\boldsymbol{\gamma}_{k})\right],$

where
$\boldsymbol{\gamma}_{k}=\\{\boldsymbol{\mu}_{k},\boldsymbol{\Sigma}_{k}\\}$,
and $C_{k}$ is the data samples in the $k^{th}$ cluster. By evaluating with
this metric, the mixture model tends to select maximum likelihood estimates.

####  12.3 Inertia

Inertia 121212The name is following from Scikit-learn document, see also
http://scikit-learn.org/stable/modules/clustering.html, or within-cluster sum-
of-squares is mostly used in K-means, in which it aims to choose centroids
that minimize the inertia function

$\sum_{k=1}^{K}\sum_{j:j\in
C_{k}}^{N_{k}}||\bm{x}_{j}-\overline{\bm{x}}_{k}||^{2}.$ (150)

where $C_{k}$ is the data samples in the $k^{th}$ cluster and
$\overline{\bm{x}}_{k}$ is the mean vector for cluster $k$. This is exactly
the same metric used in K-means. Inertia is not a normalized metric, so we
just know that lower values are better and zero is optimal. But in very high-
dimensional spaces, Euclidean distances tend to become inflated (this is an
instance of the so-called “curse of dimensionality”).

####  12.4 Squared inertia

We propose the following loss function

$\sum_{k=1}^{K}\sqrt{\sum_{j:j\in
C_{k}}^{N_{k}}||\bm{x}_{j}-\overline{\bm{x}}_{k}||^{2}}.$ (151)

where $C_{k}$ is the data samples in the $k^{th}$ cluster and
$\overline{\bm{x}}_{k}$ is the mean vector for cluster $k$. The square root
has an important impact in favoring a smaller nunber of clusters; for example,
inducing a price to be paid for introducing two clusters with the same mean.
One can imagine that if two clusters have same center value, the square root
operation will force the two clusters into one cluster. For example, if we
have the loss value for two clusters 100 and 30 respectively. If we use this
square root operation, we will get $\sqrt{100}+\sqrt{30}>\sqrt{100+30}$, thus
favoring small cluster number. However, if we do not use square root
operation, we will get $(100)+(30)=(100+30)$.

####  12.5 Label-related metrics

With the increasing popularity of algorithms for clustering, given a set of
true cluster assignments, and the set of clusters found by an algorithm, these
sets of cluster assignment can be compared to see how similar or different the
sets are. We call this as label-related metrics. A normalized measure is
desirable in many contexts, for example assigning a value of 0 where the two
sets are totally dissimilar, and 1 where they are identical (McDaid et al.,
2013). We first introduce two un-normalized measures, and a normalized measure
is described, all of which come from information theory, a field has deep
links to statistics and machine learning. A Python implementation is available
online. 131313https://github.com/junlulocky/infopy

Figure 25: Venn diagram illustrating the relation between information
entropies, mutual information and variation of information. The area
containing both circles is the joint entropy $H(X,Y)$. The circle on the left
(green and grey) is the individual entropy $H(X)$, with the green being the
conditional entropy $H(X|Y)$. The circle on the right (purple and grey) is
$H(Y)$, with the purple being $H(Y|X)$. The grey is the mutual information
$I(X;Y)$. Figure is due to wikipedia.

#####  12.5.1 Mutual information and variation of information

Formally, the mutual information of two discrete random variables $X$ and $Y$
can be defined as:

$MI(X;Y)=\sum_{y\in Y}\sum_{x\in
X}p(x,y)\log\left(\frac{p(x,y)}{p(x)p(y)}\right),$ (152)

where $p(x,y)$ is the joint probability distribution function of $X$ and $Y$,
and $p(x)$ and $p(y)$ are the marginal probability distribution functions of
$X$ and $Y$ respectively 141414See also https://nlp.stanford.edu/IR-
book/html/htmledition/evaluation-of-clustering-1.html.

Then the variation of information between the two discrete random variables is
defined as

$VI(X;Y)=-\sum_{y\in Y}\sum_{x\in
X}p(x,y)\left[\log\frac{p(x,y)}{p(x)}+\log\frac{p(x,y)}{p(y)}\right].$ (153)

Noted that unlike the mutual information, however, the variation of
information is a true metric, in that it obeys the triangle inequality.

The relationship between mutual information and variation of information can
be shown in Figure 25, where the definition of the entropy of a random
variable $X$ is $H(X)=\sum_{x\in X}p(x)\log p(x)$. From which we have the
relationship between mutual information and variation of information:
$VI(X;Y)=H(X)+H(Y)-2MI(X;Y)$.

In clustering, each clustering algorithm $C=\\{C_{1},C_{2},\ldots,C_{K}\\}$
defines the probability distribution $P_{C}$

$P_{C}(k)=\frac{n_{k}}{N},$ (154)

where $n_{k}$ is the number of points in the $k^{th}$ cluster $C_{k}$ and $N$
is the total number of points in the data set. Different clustering algorithms
can determine different number of clusters.

For any clustering distributions $P_{C_{1}}=(p_{1},p_{2},\ldots,p_{n})$ and
$P_{C_{2}}=(q_{1},q_{2},\ldots,q_{m})$, define the probability distribution
$P_{C_{1}}$ and $P_{C_{2}}$ and joint probability distribution $R(i,j)$

$\displaystyle P_{C_{1}}(i)$ $\displaystyle=\frac{n_{i}}{N},$ (155)
$\displaystyle P_{C_{2}}(j)$ $\displaystyle=\frac{m_{j}}{N},$ $\displaystyle
R(i,j)$ $\displaystyle=\frac{N_{i,j}}{N}\triangleq\frac{|n_{i}\cap
m_{j}|}{N}.$

where $|n_{i}\cap m_{j}|$ is the number of observations that is both in
cluster $i$ of $C_{1}$ and cluster $j$ of $C_{2}$. Thus, in clustering
algorithms, we set $p(x,y)=R(i,j)$, $p(x)=P_{C_{1}}(i)$ and
$p(y)=P_{C_{2}}(j)$. Then we arrive at the definition of mutual information
(MI) and variation of information (VI) (Meilă, 2003) in clustering

$MI(C_{1};C_{2})=\sum_{i,j}\frac{N_{i,j}}{N}\log\frac{\frac{N_{i,j}}{N}}{\frac{n_{i}}{N}\cdot\frac{m_{j}}{N}}=\sum_{i,j}\frac{N_{i,j}}{N}\log\frac{N\cdot
N_{i,j}}{n_{i}\cdot m_{j}},$ (156)

and

$VI(C_{1};C_{2})=-\sum_{i,j}\frac{N_{i,j}}{N}[\log\frac{\frac{N_{i,j}}{N}}{\frac{n_{i}}{N}}+\log\frac{\frac{N_{i,j}}{N}}{\frac{m_{j}}{N}}]=-\sum_{i,j}\frac{N_{i,j}}{N}[\log\frac{N_{i,j}}{n_{i}}+\log\frac{N_{i,j}}{m_{j}}].$
(157)

#####  12.5.2 Normalized mutual information

A normalized measure is desirable in many contexts, for example assigning a
value of 0 where the two sets are totally dissimilar, and 1 where they are
identical. From Figure 25, we find that the mutual information
$MI(X;Y)<[H(X)+H(Y)]/2$. Thus we normalized the mutual information to get the
normalized mutual information

$NMI(X;Y)=\frac{MI(X;Y)}{[H(X)+H(Y)]/2},$ (158)

and in clustering, we have

$NMI(C_{1};C_{2})=\frac{MI(C_{1};C_{2})}{[H(C_{1})+H(C_{2})]/2}.$ (159)

In most situations, we need to compare the clustering algorithm to a true
clustering situation, in which case we just set $C_{1}$ to be the true
clustering label. And thus, we expect the higher the mutual information (or
normalized mutual information) the better; and the lower variation of
information the better.

###  A Deriving the Dirichlet distribution

####  A.1 Derivation

Let $X_{1},X_{2},\ldots,X_{K}$ be i.i.d., random variables drawn from the
Gamma distribution such that $X_{k}\sim\mathrm{Ga}(\alpha_{k},1)$ for
$k\in\\{1,2,\ldots,K\\}$. The joint p.d.f., of $X_{1},X_{2},\ldots,X_{K}$ is
given by

$f_{X_{1},X_{2},\ldots,X_{K}}(x_{1},x_{2},\ldots,x_{K})=\left\\{\begin{aligned}
&\prod_{k=1}^{K}\frac{1}{\Gamma(\alpha_{k})}x_{k}^{\alpha_{k}-1}\exp(-x_{k}),&\mathrm{\,\,if\,\,}x_{k}\geq
0.\\\ &0,&\mathrm{\,\,if\,\,}\text{otherwise}.\end{aligned}\right.$

Define variables $Y_{k}$’s as follows

$\displaystyle Y_{k}$
$\displaystyle=\frac{X_{k}}{\sum_{k=1}^{K}X_{k}},\qquad\forall\,\,k\in\\{1,2,\ldots,K-1\\}$
(160) $\displaystyle Y_{K}$
$\displaystyle=\frac{X_{K}}{\sum_{k=1}^{K}X_{k}}=1-\sum_{k=1}^{K-1}Y_{k},$

and

$Z_{K}=\sum_{k=1}^{K}X_{k}.$ (161)

Let $\bm{X}=[X_{1},X_{2},\ldots,X_{K}]$,
$\bm{Y}=[Y_{1},Y_{2},\ldots,Y_{K-1},Z_{K}]$,
$\bm{x}=[x_{1},x_{2},\ldots,x_{K}]$, and
$\bm{y}=[y_{1},y_{2},\ldots,y_{K-1},z_{K}]$. By multidimensional
transformation of variables, we have

$f_{\bm{Y}}(\bm{y})=f_{\bm{X}}(g^{-1}(\bm{y}))\left|\det\left[J_{g^{-1}}(\bm{y})\right]\right|,$

where

$\begin{bmatrix}x_{1}\\\ x_{2}\\\ \vdots\\\ x_{K-1}\\\
x_{K}\end{bmatrix}=g^{-1}(\bm{y})=g^{-1}\left(\begin{bmatrix}y_{1}\\\ y_{2}\\\
\vdots\\\ y_{K-1}\\\ z_{K}\end{bmatrix}\right)=\begin{bmatrix}y_{1}\cdot
z_{K}\\\ y_{2}\cdot z_{K}\\\ \vdots\\\ y_{K-1}\cdot z_{K}\\\ y_{K}\cdot
z_{K}\end{bmatrix},$

and the Jacobian matrix is given by

$\displaystyle J_{g^{-1}}(\bm{y})$
$\displaystyle=\begin{bmatrix}\frac{\partial}{\partial
y_{1}}g_{1}^{-1}(\bm{y})&\cdots&\frac{\partial}{\partial
y_{K-1}}g_{1}^{-1}(\bm{y})&\frac{\partial}{\partial
z_{K}}g_{1}^{-1}(\bm{y})\\\ \vdots&\ddots&\cdots&\vdots\\\
\frac{\partial}{\partial
y_{1}}g_{K}^{-1}(\bm{y})&\cdots&\frac{\partial}{\partial
y_{K-1}}g_{K}^{-1}(\bm{y})&\frac{\partial}{\partial
z_{K}}g_{K}^{-1}(\bm{y})\\\ \end{bmatrix}$
$\displaystyle=\begin{bmatrix}z_{K}&0&\cdots&0&y_{1}\\\
0&z_{K}&\cdots&0&y_{2}\\\ \vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&\cdots&z_{K}&y_{K-1}\\\
-z_{K}&-z_{k}&\cdots&-z_{K}&(1-\sum_{k=1}^{K-1}y_{k})\end{bmatrix}=z_{K}^{K-1}.$

This implies the joint p.d.f, of $\bm{Y}$ is

$f_{\bm{Y}}(\bm{y})=f_{\bm{X}}(g^{-1}(\bm{y}))z_{K}^{K-1}=\frac{y_{1}^{\alpha_{1}-1}y_{2}^{\alpha_{2}-1}\ldots
y_{K-1}^{\alpha_{K-1}-1}y_{K}^{\alpha_{K}-1}}{\prod_{k=1}^{K}\Gamma(\alpha_{k})}\exp(-z_{K})z_{K}^{\alpha_{1}+\alpha_{2}+\ldots+\alpha_{K}-1}.$

We realize that the righthand size of above equation is proportional to a
p.d.f. of Gamma distribution and

$\int\exp(-z_{K})z_{K}^{\alpha_{1}+\alpha_{2}+\ldots+\alpha_{K}-1}dz_{K}=\Gamma(\alpha_{1}+\alpha_{2}+\ldots+\alpha_{K}).$

Let $\alpha_{+}=\alpha_{1}+\alpha_{2}+\ldots+\alpha_{K}$, this implies

$f(y_{1},y_{2},\ldots,y_{K-1})=\frac{\Gamma(\alpha_{+})}{\prod_{k=1}^{K}\Gamma(\alpha_{k})}\prod_{k=1}^{K}y_{k}^{\alpha_{k}-1}.$

We notice that $Y_{k}$’s are defined that $0<Y_{k}<1$ for all
$k\in\\{1,2,\ldots,K\\}$, and $\sum_{k=1}^{K}Y_{k}=1$. This implies the above
equation is the p.d.f., of the Dirichlet distribution. The construction shown
above can be utilized to generate random variables from the Dirichlet
distribution.

####  A.2 Properties of Dirichlet distribution

Suppose
$\bm{Y}=[Y_{1},Y_{2},\ldots,Y_{K}]\sim\mathrm{Dirichlet}(\boldsymbol{\alpha})$
with $\boldsymbol{\alpha}=[\alpha_{1},\alpha_{2},\ldots,\alpha_{K}]$, we here
show the moments and properties of the Dirichlet distribution.

###### Mean of Dirichlet distribution

Write out the expectation:

$\displaystyle\mathrm{E}[Y_{1}]$ $\displaystyle=\int\cdots\int
y_{1}\cdot\mathrm{Dirichlet}(\bm{y}|\boldsymbol{\alpha})dy_{1}dy_{2}\cdots
dy_{K}$ $\displaystyle=\int\cdots\int
y_{1}\frac{\Gamma(\alpha_{+})}{\prod_{k=1}^{K}\Gamma(\alpha_{k})}\prod_{k=1}^{K}y_{k}^{\alpha_{k}-1}dy_{1}dy_{2}\cdots
dy_{K}$
$\displaystyle=\frac{\Gamma(\alpha_{+})}{\prod_{k=1}^{K}\Gamma(\alpha_{k})}\int\cdots\int
y_{1}^{\alpha_{1}+1-1}\prod_{k=2}^{K}y_{k}^{\alpha_{k}-1}dy_{1}dy_{2}\cdots
dy_{K}$
$\displaystyle=\frac{\Gamma(\alpha_{+})}{\prod_{k=1}^{K}\Gamma(\alpha_{k})}\frac{\Gamma(\alpha_{1}+1)\prod_{k=2}^{K}\Gamma(\alpha_{k})}{\Gamma(\alpha_{+}+1)}$
$\displaystyle=\frac{\Gamma(\alpha_{+})}{\Gamma(\alpha_{1})}\frac{\Gamma(\alpha_{1}+1)}{\Gamma(\alpha_{+}+1)}$
$\displaystyle=\frac{\alpha_{1}}{\alpha_{+}},$

where the last equality comes from the fact that $\Gamma(x+1)=x\Gamma(x)$.

###### Variance of Dirichlet distribution

Write out the variance
$\mathrm{Var}[Y_{i}]=\mathrm{E}[Y_{i}^{2}]-\mathrm{E}[Y_{i}]^{2}$. Similarly
from the proof of the mean, we have

$\mathrm{E}[Y_{i}^{2}]=\frac{\Gamma(\alpha_{+})}{\Gamma(\alpha_{+}+2)}\frac{\Gamma(\alpha_{i}+2)}{\Gamma(\alpha_{i})}=\frac{(\alpha_{i}+1)\alpha_{i}}{(\alpha_{+}+1)\alpha_{+}}.$

This implies

$\mathrm{Var}[Y_{i}]=\mathrm{E}[Y_{i}^{2}]-\mathrm{E}[Y_{i}]^{2}=\frac{(\alpha_{i}+1)\alpha_{i}}{(\alpha_{+}+1)\alpha_{+}}-(\frac{\alpha_{i}}{\alpha_{+}})^{2}=\frac{\alpha_{i}(\alpha_{+}-\alpha_{i})}{\alpha_{+}^{2}(\alpha_{+}+1)}.$

###### Covariance of Dirichlet distribution

Write out the covariance
$\mathrm{Cov}[Y_{i}Y_{j}]=\mathrm{E}[Y_{i}Y_{j}]-\mathrm{E}[Y_{i}]\mathrm{E}[Y_{j}]$.
Again, similarly from the proof of the mean, for $i\neq j$, we have

$\mathrm{E}[Y_{i}Y_{j}]=\frac{\Gamma(\alpha_{+})}{\Gamma(\alpha_{+}+2)}\frac{\Gamma(\alpha_{i}+1)}{\Gamma(\alpha_{i})}\frac{\Gamma(\alpha_{j}+1)}{\Gamma(\alpha_{j})}=\frac{\alpha_{i}\alpha_{j}}{\alpha_{+}(\alpha_{+}+1)}.$

This implies

$\mathrm{Cov}[Y_{i}Y_{j}]=\mathrm{E}[Y_{i}Y_{j}]-\mathrm{E}[Y_{i}]\mathrm{E}[Y_{j}]=\frac{\alpha_{i}\alpha_{j}}{\alpha_{+}(\alpha_{+}+1)}-\frac{\alpha_{i}\alpha_{j}}{\alpha_{+}^{2}}=\frac{-\alpha_{i}\alpha_{j}}{\alpha_{+}^{2}(\alpha_{+}+1)}.$

###### Marginal distribution of $Y_{i}$

By definition in Equation (160) and Equation (161), we have
$Z_{K}-X_{i}\sim\mathrm{Ga}(\alpha_{+}-\alpha_{i},1)$. This implies

$Y_{i}=\frac{X_{i}}{Z_{K}}=\frac{X_{i}}{X_{i}+(Z_{K}-X_{i})}\sim\mathrm{Beta}(\alpha_{i},\alpha_{+}-\alpha_{i}).$

which is from the fact about the p.d.f., of two independent Gamma random
variables. 151515Suppose $X\sim\mathrm{Ga}(a,\lambda)$ and
$Y\sim\mathrm{Ga}(b,\lambda)$, then $\frac{X}{X+Y}\sim\mathrm{Beta}(a,b)$.

###### Aggregation property

Suppose
$[Y_{1},Y_{2},\ldots,Y_{K}]\sim\mathrm{Dirichlet}([\alpha_{1},\alpha_{2},\ldots,\alpha_{K}])$,
Then, Let $M=Y_{i}+Y_{j}$, it follows that

$\displaystyle\,\,\,\,\,\,\,\,\,[Y_{1},\ldots
Y_{i-1},Y_{i+1},\ldots,Y_{j-1},Y_{j+1},\ldots,Y_{K},M]$
$\displaystyle\sim\mathrm{Dirichlet}([\alpha_{1},\ldots,\alpha_{i-1},\alpha_{i+1},\ldots,\alpha_{j-1},\alpha_{j+1},\ldots,\alpha_{K},\alpha_{i}+\alpha_{j}]).$

Proof  We realize that $M\sim\mathrm{Ga}(\alpha_{i}+\alpha_{j},1)$. Again by
the multidimensional transformation of variables as shown in the beginning of
this section, we conclude the result.
The results can be extended to a more general case. If
$\\{A_{1},A_{2},\ldots,A_{r}\\}$ is a partition of $\\{1,2,\ldots,K\\}$, then

$\left[\sum_{i\in A_{1}}Y_{i},\sum_{i\in A_{2}}Y_{i},\ldots,\sum_{i\in
A_{r}}Y_{i}\right]\sim\mathrm{Dirichlet}\left(\left[\sum_{i\in
A_{1}}\alpha_{i},\sum_{i\in A_{2}}\alpha_{i},\ldots,\sum_{i\in
A_{r}}\alpha_{i}\right]\right).$

###### Condition distribution

Let $Y_{0}=\sum_{k=3}^{K}Y_{i}$ and
$\alpha_{0}=\alpha_{+}-\alpha_{1}-\alpha_{2}$, then
$[Y_{1},Y_{2},Y_{0}]\sim\mathrm{Dirichlet}([\alpha_{1},\alpha_{2},\alpha_{0}])$.
Therefore

$f_{Y_{1},Y_{2}}(y_{1},y_{2})=\frac{\Gamma(\alpha_{1}+\alpha_{2}+\alpha_{0})}{\Gamma(\alpha_{1})\Gamma(\alpha_{2})\Gamma(\alpha_{0})}y_{1}^{\alpha_{1}-1}y_{2}^{\alpha_{2}-1}(1-y_{1}-y_{2})^{\alpha_{0}-1}.$

Similarly, we have

$f_{Y_{2}}(y_{2})=\frac{\Gamma(\alpha_{1}+\alpha_{2}+\alpha_{0})}{\Gamma(\alpha_{2})\Gamma(\alpha_{1}+\alpha_{0})}y_{2}^{\alpha_{2}-1}(1-y_{2})^{\alpha_{1}+\alpha_{0}-1}=\mathrm{Beta}(y_{1}|\alpha_{2},\alpha_{1}+\alpha_{0}),$

which is a p.d.f., of a Beta distribution. Therefore, the conditional p.d.f.,
of $Y_{1}|Y_{2}=y_{2}$ is given by

$f_{Y_{1}|Y_{2}=y_{2}}(y_{1}|y_{2})=\frac{f_{Y_{1},Y_{2}}(y_{1},y_{2})}{f_{Y_{2}}(y_{2})}=\frac{\Gamma(\alpha_{1}+\alpha_{0})}{\Gamma(\alpha_{1})\Gamma(\alpha_{0})}\left(\frac{y_{1}}{1-y_{2}}\right)^{\alpha_{1}-1}\left(1-\frac{y_{1}}{1-y_{2}}\right)^{\alpha_{0}-1}\frac{1}{1-y_{2}},$

which implies

$\frac{1}{1-y_{2}}Y_{1}|Y_{2}=y_{2}\sim\mathrm{Beta}(\alpha_{1},\alpha_{0}).$

Apply this procedure, we will have

$\bm{Y}_{-i}|Y_{i}\sim(1-y_{i})\mathrm{Dirichlet}(\alpha_{-i}),$

where $Y_{-i}$ is all the $K-1$ variables except $Y_{i}$, and similarly for
$\alpha_{-i}$.

###  B Cholesky decomposition

######  Theorem 24 (Cholesky Decomposition)

Every positive definite matrix $\bm{A}\in{}^{n\times n}$ can be factored as

$\bm{A}=\bm{R}^{\top}\bm{R},$

where $\bm{R}$ is an upper triangular matrix with positive diagonal elements.
This decomposition is known as Cholesky decomposition of $\bm{A}$. $\bm{R}$ is
known as the Cholesky factor or Cholesky triangle of $\bm{A}$.

####  B.1 Existence of the Cholesky decomposition

Before showing the existence of Cholesky decomposition, we need the following
definitions and lemmas.

######  Definition 25 (Positive Definite and Positive Semidefinite)

A matrix $\bm{A}\in{}^{n\times n}$ is positive definite if
$\bm{x}^{\top}\bm{A}\bm{x}>0$ for all nonzero $\bm{x}\in{}^{n}$. And a matrix
$\bm{A}\in{}^{n\times n}$ is positive semidefinite if
$\bm{x}^{\top}\bm{A}\bm{x}\geq 0$ for all $\bm{x}\in{}^{n}$.

######  Lemma 26 (Positive Diagonals of Positive Definite Matrices)

The diagonal elements of a positive definite matrix $\bm{A}$ are all positive.

Proof [of Lemma 26] From the definition of positive definite matrix, we have
$\bm{x}^{\top}\bm{A}\bm{x}>0$ for all nonzero $\bm{x}$. In particular, let
$\bm{x}=\bm{e}_{i}$ where $\bm{e}_{i}$ is the $i$-th unit vector with the
$i$-th entry equal to 1 and other entries equal to 0. Then,

$\bm{e}_{i}^{\top}\bm{A}\bm{e}_{i}=\bm{A}_{ii}>0,\qquad\forall
i\in\\{1,2,\cdots,n\\}.$

This completes the proof.

######  Lemma 27 (Schur Complement of Positive Definite Matrices)

For any positive definite matrix $\bm{A}\in{}^{n\times n}$, its Schur
complement of $\bm{A}_{11}$ is
$\bm{S}_{n-1}=\bm{A}_{2:n,2:n}-\frac{1}{\bm{A}_{11}}\bm{A}_{2:n,1}\bm{A}_{2:n,1}^{\top}$
and it is also positive definite.

Note that the subscript $n-1$ of $\bm{S}_{n-1}$ means it is of size
$(n-1)\times(n-1)$ and it is a Schur complement of a $n\times n$ positive
definite matrix. We will use this notation in the following section.

Proof [of Lemma 27] For any nonzero vector $\bm{v}\in{}^{n-1}$, we can
construct a vector $\bm{x}\in{}^{n}$

$\bm{x}=\begin{bmatrix}-\frac{1}{\bm{A}_{11}}\bm{A}_{2:n,1}^{\top}\bm{v}\\\
\bm{v}\end{bmatrix},$

which is nonzero. Then

$\displaystyle\bm{x}^{\top}\bm{A}\bm{x}$
$\displaystyle=[-\frac{1}{\bm{A}_{11}}\bm{v}^{\top}\bm{A}_{2:n,1}\qquad\bm{v}^{\top}]\begin{bmatrix}\bm{A}_{11}&\bm{A}_{2:n,1}^{\top}\\\
\bm{A}_{2:n,1}&\bm{A}_{2:n,2:n}\end{bmatrix}\begin{bmatrix}-\frac{1}{\bm{A}_{11}}\bm{A}_{2:n,1}^{\top}\bm{v}\\\
\bm{v}\end{bmatrix}$
$\displaystyle=[-\frac{1}{\bm{A}_{11}}\bm{v}^{\top}\bm{A}_{2:n,1}\qquad\bm{v}^{\top}]\begin{bmatrix}0\\\
\bm{S}_{n-1}\bm{v}\end{bmatrix}$
$\displaystyle=\bm{v}^{\top}\bm{S}_{n-1}\bm{v}.$

Since $\bm{A}$ is positive definite, we have
$\bm{x}^{\top}\bm{A}\bm{x}=\bm{v}^{\top}\bm{S}_{n-1}\bm{v}>0$ for all nonzero
$\bm{v}$. Thus, $\bm{S}$ is positive definite.

A word on the Schur complement: It can be easily proved that this Schur
complement
$\bm{S}_{n-1}=\bm{A}_{2:n,2:n}-\frac{1}{\bm{A}_{11}}\bm{A}_{2:n,1}\bm{A}_{2:n,1}^{\top}$
is also nonsingular if $\bm{A}$ is nonsngular and $\bm{A}_{11}\neq 0$.
Similarly, the Schur complement of $\bm{A}_{nn}$ in $\bm{A}$ is
$\bar{\bm{S}}_{n-1}=\bm{A}_{1:n-1,1:n-1}-\frac{1}{\bm{A}_{nn}}\bm{A}_{1:n-1,n}\bm{A}_{1:n-1,n}^{\top}$
which is also positive definite if $\bm{A}$ is positive definite. This
property can help prove the leading principle minors of positive definite
matrices are all positive. See Appendix C for more details.

We then prove the existence of Cholesky decomposition using these lemmas.

Proof [of Theorem 24: Existence of Cholesky Decomposition] For any positive
definite matrix $\bm{A}$, we can write out (since $\bm{A}_{11}$ is positive)

$\displaystyle\bm{A}$
$\displaystyle=\begin{bmatrix}\bm{A}_{11}&\bm{A}_{2:n,1}^{\top}\\\
\bm{A}_{2:n,1}&\bm{A}_{2:n,2:n}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\sqrt{\bm{A}_{11}}&\boldsymbol{0}\\\
\frac{1}{\sqrt{\bm{A}_{11}}}\bm{A}_{2:n,1}&\bm{I}\end{bmatrix}\begin{bmatrix}\sqrt{\bm{A}_{11}}&\frac{1}{\sqrt{\bm{A}_{11}}}\bm{A}_{2:n,1}^{\top}\\\
\boldsymbol{0}&\bm{A}_{2:n,2:n}-\frac{1}{\bm{A}_{11}}\bm{A}_{2:n,1}\bm{A}_{2:n,1}^{\top}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\sqrt{\bm{A}_{11}}&\boldsymbol{0}\\\
\frac{1}{\sqrt{\bm{A}_{11}}}\bm{A}_{2:n,1}&\bm{I}\end{bmatrix}\begin{bmatrix}1&\boldsymbol{0}\\\
\boldsymbol{0}&\bm{A}_{2:n,2:n}-\frac{1}{\bm{A}_{11}}\bm{A}_{2:n,1}\bm{A}_{2:n,1}^{\top}\end{bmatrix}\begin{bmatrix}\sqrt{\bm{A}_{11}}&\frac{1}{\sqrt{\bm{A}_{11}}}\bm{A}_{2:n,1}^{\top}\\\
\boldsymbol{0}&\bm{I}\end{bmatrix}$
$\displaystyle=\bm{R}_{1}^{\top}\begin{bmatrix}1&\boldsymbol{0}\\\
\boldsymbol{0}&\bm{S}_{n-1}\end{bmatrix}\bm{R}_{1}.$

where

$\bm{R}_{1}=\begin{bmatrix}\sqrt{\bm{A}_{11}}&\frac{1}{\sqrt{\bm{A}_{11}}}\bm{A}_{2:n,1}^{\top}\\\
\boldsymbol{0}&\bm{I}\end{bmatrix}.$

Since we proved the Schur complement $\bm{S}_{n-1}$ is positive definite. We
can factor it in the same way

$\bm{S}_{n-1}=\hat{\bm{R}}_{2}^{\top}\begin{bmatrix}1&\boldsymbol{0}\\\
\boldsymbol{0}&\bm{S}_{n-2}\end{bmatrix}\hat{\bm{R}}_{2}.$

We then have

$\displaystyle\bm{A}$
$\displaystyle=\bm{R}_{1}^{\top}\begin{bmatrix}1&\boldsymbol{0}\\\
\boldsymbol{0}&\hat{\bm{R}}_{2}^{\top}\begin{bmatrix}1&\boldsymbol{0}\\\
\boldsymbol{0}&\bm{S}_{n-2}\end{bmatrix}\hat{\bm{R}}_{2}.\end{bmatrix}\bm{R}_{1}$
$\displaystyle=\bm{R}_{1}^{\top}\begin{bmatrix}1&\boldsymbol{0}\\\
\boldsymbol{0}&\hat{\bm{R}}_{2}^{\top}\end{bmatrix}\begin{bmatrix}1&\boldsymbol{0}\\\
\boldsymbol{0}&\begin{bmatrix}1&\boldsymbol{0}\\\
\boldsymbol{0}&\bm{S}_{n-2}\end{bmatrix}\end{bmatrix}\begin{bmatrix}1&\boldsymbol{0}\\\
\boldsymbol{0}&\hat{\bm{R}}_{2}\end{bmatrix}\bm{R}_{1}$
$\displaystyle=\bm{R}_{1}^{\top}\bm{R}_{2}^{\top}\begin{bmatrix}1&\boldsymbol{0}\\\
\boldsymbol{0}&\begin{bmatrix}1&\boldsymbol{0}\\\
\boldsymbol{0}&\bm{S}_{n-2}\end{bmatrix}\end{bmatrix}\bm{R}_{2}\bm{R}_{1}.$

The same formula can be recursively applied. This process gradually continues
down to the bottom-right corner giving us the decomposition

$\displaystyle\bm{A}$
$\displaystyle=\bm{R}_{1}^{\top}\bm{R}_{2}^{\top}\cdots\bm{R}_{n}^{\top}\bm{R}_{n}\cdots\bm{R}_{2}\bm{R}_{1}$
$\displaystyle=\bm{R}^{\top}\bm{R},$

where $\bm{R}_{1},\bm{R}_{2},\cdots\bm{R}_{n},$ are upper triangular matrices
with positive diagonal elements and
$\bm{R}=\bm{R}_{1}\bm{R}_{2}\cdots\bm{R}_{n}$ is also an upper triangular
matrix with positive diagonal elements.
The process in the proof can also be used to compute the Cholesky
decomposition. In next section, we use another point of view to do the
computation.

######  Corollary 28 ($\bm{R}^{\top}\bm{R}$ is PD)

For any upper triangular matrix with positive diagonal elements, then
$\bm{A}=\bm{R}^{\top}\bm{R}$ is positive definite.

Proof [of Corollary 28] If an upper triangular matrix $\bm{R}$ has positive
diagonal, it is full column rank, and the null space of $\bm{R}$ is 0. As a
result, $\bm{R}\bm{x}\neq\boldsymbol{0}$ for any nonzero vector $\bm{x}$. Thus
$\bm{x}^{\top}\bm{A}\bm{x}=||\bm{R}\bm{x}||^{2}>0$ for any nonzero vector
$\bm{x}$.
This corollary can be extended to any $\bm{R}$ with independent columns.

####  B.2 Computing the Cholesky decomposition

To compute Cholesky decomposition, we write out the equality
$\bm{A}=\bm{R}^{\top}\bm{R}$:

$\displaystyle\bm{A}=\left[\begin{matrix}\bm{A}_{11}&\bm{A}_{1,2:n}\\\
\bm{A}_{2:n,1}&\bm{A}_{2:n,2:n}\end{matrix}\right]$
$\displaystyle=\left[\begin{matrix}\bm{R}_{11}&0\\\
\bm{R}_{1,2:n}^{\top}&\bm{R}_{2:n,2:n}^{\top}\end{matrix}\right]\left[\begin{matrix}\bm{R}_{11}&\bm{R}_{1,2:n}\\\
0&\bm{R}_{2:n,2:n}\end{matrix}\right]$
$\displaystyle=\left[\begin{matrix}\bm{R}_{11}^{2}&\bm{R}_{11}\bm{R}_{1,2:n}\\\
\bm{R}_{11}\bm{R}_{1,2:n}^{\top}&\bm{R}_{1,2:n}^{\top}\bm{R}_{1,2:n}+\bm{R}_{2:n,2:n}^{\top}\bm{R}_{2:n,2:n}\end{matrix}\right],$

which allows to determine the first row of $\bm{R}$

$\bm{R}_{11}=\sqrt{\bm{A}_{11}},\qquad\bm{R}_{1,2:n}=\frac{1}{\bm{R}_{11}}\bm{A}_{1,2:n}.$

Let $\bm{A}_{2}=\bm{R}_{2:n,2:n}^{\top}\bm{R}_{2:n,2:n}$. The equality
$\bm{A}_{2:n,2:n}=\bm{R}_{1,2:n}^{\top}\bm{R}_{1,2:n}+\bm{R}_{2:n,2:n}^{\top}\bm{R}_{2:n,2:n}$
gives out

$\displaystyle\bm{A}_{2}=\bm{R}_{2:n,2:n}^{\top}\bm{R}_{2:n,2:n}$
$\displaystyle=\bm{A}_{2:n,2:n}-\bm{R}_{1,2:n}^{\top}\bm{R}_{1,2:n}$
$\displaystyle=\bm{A}_{2:n,2:n}-\frac{1}{\bm{A}_{11}}\bm{A}_{1,2:n}^{\top}\bm{A}_{1,2:n}$
$\displaystyle=\bm{A}_{2:n,2:n}-\frac{1}{\bm{A}_{11}}\bm{A}_{2:n,1}\bm{A}_{1,2:n}\qquad(\bm{A}\mbox{
is positive definite}).$

$\bm{A}_{2}$ is the Schur complement of $\bm{A}_{11}$ in $\bm{A}$ of size
$(n-1)\times(n-1)$. And to get $\bm{R}_{2:n,2:n}$ we must compute the Cholesky
decomposition of matrix $\bm{A}_{2}$ of $(n-1)\times(n-1)$. Again, this is a
recursive algorithm and formulated in Algorithm 8.

1:positive definite matrix $\bm{A}$ with size $n\times n$;

2:Calculate first of $\bm{R}$ by
$\bm{R}_{11}=\sqrt{\bm{A}_{11}},\bm{R}_{1,2:n}=\frac{1}{\bm{R}_{11}}\bm{A}_{1,2:n}$,
($n$ flops);

3:Compute the Cholesky decomposition of the $(n-1)\times(n-1)$ matrix
$\bm{A}_{2}=\bm{R}_{2:n,2:n}^{\top}\bm{R}_{2:n,2:n}=\bm{A}_{2:n,2:n}-\frac{1}{\bm{A}_{11}}\bm{A}_{2:n,1}\bm{A}_{1,2:n},\qquad\text{($n^{2}-n$
flops).}$

Algorithm 8 Cholesky Decomposition

Further, this process can be used to determine if a matrix is positive
definite or not. If we try to factor a non positive definite matrix, at some
point, we will encounter a nonpositive element in entry (1,1) of
$\bm{A},\bm{A}_{2},\bm{A}_{3},\cdots$.

###### Theorem 29

Algorithm 8 requires $\sim(1/3)n^{3}$ flops to compute a Cholesky
decomposition of an $n\times n$ positive definite matrix.

Proof [of Theorem 29] Step 1 takes 1 square root and $(n-1)$ division which
takes $n$ flops totally.

For step 2, Note that
$\frac{1}{\bm{A}_{11}}\bm{A}_{2:n,1}\bm{A}_{1,2:n}=(\frac{1}{\sqrt{\bm{A}_{11}}}\bm{A}_{2:n,1})(\frac{1}{\sqrt{\bm{A}_{11}}}\bm{A}_{1,2:n})=\bm{R}_{1,2:n}^{\top}\bm{R}_{1,2:n}$.
If we calculate the complexity directly from the equation in step 2, we will
get the same complexity as LU decomposition. But since
$\bm{R}_{1,2:n}^{\top}\bm{R}_{1,2:n}$ is symmetric, the complexity of
$\bm{R}_{1,2:n}^{\top}\bm{R}_{1,2:n}$ reduces from $(n-1)\times(n-1)$
multiplications to $1+2+\cdots+(n-1)=\frac{n^{2}-n}{2}$ multiplications. The
cost of matrix division reduces from $(n-1)\times(n-1)$ to
$1+2+\cdots+(n-1)=\frac{n^{2}-n}{2}$ as well. So it costs $n^{2}-n$ flops for
step 2.

Simple calculation will show the total complexity is
$\frac{2n^{3}+3n^{2}+n}{6}$ flops which is $(1/3)n^{3}$ flops if we keep only
the leading term.

###  C Leading Principle Minors of PD Matrices

In Lemma 27, we proved for any positive definite matrix $\bm{A}\in{}^{n\times
n}$, it’s Schur complement of $\bm{A}_{11}$ is
$\bm{S}_{n-1}=\bm{A}_{2:n,2:n}-\frac{1}{\bm{A}_{11}}\bm{A}_{2:n,1}\bm{A}_{2:n,1}^{\top}$
and it is also positive definite. This is also true for its Schur complement
of $\bm{A}_{nn}$, i.e.,
$\bm{S}_{n-1}^{\prime}=\bm{A}_{1:n-1,1:n-1}-\frac{1}{\bm{A}_{nn}}\bm{A}_{1:n-1,n}\bm{A}_{1:n-1,n}^{\top}$
is also positive definite.

We then claim all the leading principle minors of a positive definite matrix
$\bm{A}\in{}^{n\times n}$ are positive.

Proof  We will prove by induction. Since all the components on the diagonal of
positive definite matrices are all positive (see Lemma 26). The case for $n=1$
is trivial that $\det(\bm{A}_{11})>0$.

Suppose all the leading principle minors for $k\times k$ matrices are all
positive. If we could prove this is also true for $(k+1)\times(k+1)$ matrices,
then we complete the proof.

For a $(k+1)\times(k+1)$ matrix $\bm{M}=\begin{bmatrix}\bm{A}&\bm{b}\\\
\bm{b}^{\top}&d\end{bmatrix}$, where $\bm{A}$ is a $k\times k$ submatrix. Then
its Schur complement of $d$,
$\bm{S}_{k}=\bm{A}-\frac{1}{d}\bm{b}\bm{b}^{\top}$ is also positive definite
and its determinant is positive from the assumption. And
$\det(\bm{M})=\det(d)\det(\bm{A}-\frac{1}{d}\bm{b}\bm{b}^{\top})$= 161616By
the fact that if matrix $\bm{M}$ has a block formulation:
$\bm{M}=\begin{bmatrix}\bm{A}&\bm{B}\\\ \bm{C}&\bm{D}\end{bmatrix}$, then
$\det(\bm{M})=\det(\bm{D})\det(\bm{A}-\bm{B}\bm{D}^{-1}\bm{C})$.
$d\cdot\det(\bm{A}-\frac{1}{d}\bm{b}\bm{b}^{\top})>0$, which completes the
proof.

###  D Convexity results

We prove that $x\mapsto\Gamma(Kx)/(x^{K-1}[\Gamma(x)]^{K})$ is strictly log-
convex, $x\mapsto\Gamma(Kx)/[\Gamma(x)]^{K}$ is strictly log-concave and the
function $x\mapsto\Gamma(x+t)/\Gamma(x)$ is also strictly log-concave where
$\Gamma(x)$ is the Gamma function.

###### Theorem 30

Define

$F(x)=\frac{\Gamma(Kx)}{x^{K-1}[\Gamma(x)]^{K}},\quad
G(x)=\frac{\Gamma(Kx)}{[\Gamma(x)]^{K}}.$ (162)

For $x>0$ and $K$ is an arbitrary positive integer, the function F is strictly
log-convex and the function G is strictly log-concave.

Proof  Follow from (Abramowitz et al., 1966) we get
$\Gamma(Kx)=(2\pi)^{\frac{1}{2}(1-K)}K^{Kx-\frac{1}{2}}\prod_{i=0}^{K-1}\Gamma(x+\frac{i}{K})$.
Then

$\log F(x)=\frac{1}{2}(1-K)\log(2\pi)+(Kx-\frac{1}{2})\log
K+\sum_{i=0}^{K-1}\log\Gamma(x+\frac{i}{K})-K\log\Gamma(x)-(K-1)\log x$ (163)

and

$[\log F(x)]^{\prime}=K\log
K+\sum_{i=0}^{K-1}\Psi(x+\frac{i}{K})-K\Psi(x)-(K-1)x^{-1},$ (164)

where $\Psi(x)$ is the Digamma function, and

$\Psi^{\prime}(x)=\sum_{h=0}^{\infty}\frac{1}{(x+h)^{2}}.$ (165)

Thus

$[\log
F(x)]^{\prime\prime}=\left[\sum_{i=0}^{K-1}\Psi^{\prime}(x+\frac{i}{K})\right]-K\Psi^{\prime}(x)+\frac{K-1}{x^{2}}>0,\quad(x>0).$
(166)

The last inequality comes from (165) (also, we can find the derivative of
Digamma function is monotone decreasing). Easily, we can get

$\log G(x)=\frac{1}{2}(1-K)\log(2\pi)+(Kx-\frac{1}{2})\log
K+\sum_{i=0}^{K-1}\log\Gamma(x+\frac{i}{K})-K\log\Gamma(x)$ (167)

and

$[\log G(x)]^{\prime}=K\log K+\sum_{i=0}^{K-1}\Psi(x+\frac{i}{K})-K\Psi(x),$
(168)

Thus

$[\log
G(x)]^{\prime\prime}=\left[\sum_{i=0}^{K-1}\Psi^{\prime}(x+\frac{i}{K})\right]-K\Psi^{\prime}(x)<0,\quad(x>0).$
(169)

This concludes the theorem.
This theorem is a general case of Theorem 1 in (Merkle, 1997).

###### Theorem 31

Define

$H(x)=\frac{\Gamma(x+t)}{\Gamma(x)}.$ (170)

For $x>0$ and $t$ is a constant that $x+t>0$ (or for simplicity we can let
$t>0$), the function H is strictly log-concave.

Proof  We can easily get

$\log H(x)=\log\Gamma(x+t)-\log\Gamma(x)$ (171)

and

$[\log H(x)]^{\prime}=\Psi(x+t)-\Psi(x),$ (172)

where $\Psi(x)$ is the Digamma function, and

$\Psi^{\prime}(x)=\sum_{h=0}^{\infty}\frac{1}{(x+h)^{2}}.$ (173)

Thus

$[\log
H(x)]^{\prime\prime}=\Psi^{\prime}(x+t)-\Psi^{\prime}(x)<0,\quad(x>0,t>0).$
(174)

The last inequality comes from (173) which is monotone decreasing and
concludes the theorem.

### References

  * Abramowitz et al. (1966) Milton Abramowitz, Irene A Stegun, et al.  Handbook of mathematical functions.  _Applied mathematics series_ , 55(62):39, 1966\.
  * Anderson (1991) John R Anderson.  The adaptive nature of human categorization.  _Psychological Review_ , 98(3):409, 1991.
  * Anderson (2003) Theodore Wilbur Anderson.  An introduction to multivariate statistical analysis.  Technical report, Wiley New York, 2003.
  * Andrieu et al. (2003) Christophe Andrieu, Nando De Freitas, Arnaud Doucet, and Michael I Jordan.  An introduction to mcmc for machine learning.  _Machine learning_ , 50(1):5–43, 2003.
  * Antoniak (1974) Charles E Antoniak.  Mixtures of dirichlet processes with applications to bayesian nonparametric problems.  _The annals of statistics_ , pages 1152–1174, 1974.
  * Bernardo et al. (2003) JM Bernardo, MJ Bayarri, JO Berger, AP Dawid, D Heckerman, AFM Smith, and M West.  Bayesian clustering with variable and transformation selections.  In _Bayesian Statistics 7: Proceedings of the Seventh Valencia International Meeting_ , page 249. Oxford University Press, USA, 2003.
  * Bishop (2006) Christopher M Bishop.  Pattern recognition.  _Machine learning_ , 128(9), 2006.
  * Blei and Frazier (2011) David M Blei and Peter I Frazier.  Distance dependent Chinese restaurant processes.  _Journal of Machine Learning Research_ , 12:2461–2488, August 2011.
  * Bonawitz (2008) Keith Allen Bonawitz.  _Composable Probabilistic Inference with Blaise_.  PhD thesis, Massachusetts Institute of Technology, 2008.
  * Box and Draper (1987) George EP Box and Norman R Draper.  _Empirical model-building and response surfaces._ John Wiley & Sons, 1987.
  * Chipman et al. (2001) Hugh Chipman, Edward I George, Robert E McCulloch, Merlise Clyde, Dean P Foster, and Robert A Stine.  The practical implementation of bayesian model selection.  _Lecture Notes-Monograph Series_ , pages 65–134, 2001.
  * Das (2014) Rajarshi Das.  Collapsed gibbs sampler for dirichlet process gaussian mixture models (dpgmm).  _Talk_ , 2014.
  * De Blasi et al. (2015) Pierpaolo De Blasi, Stefano Favaro, Antonio Lijoi, Ramsés H Mena, Igor Prünster, and Matteo Ruggiero.  Are Gibbs-type priors the most natural generalization of the Dirichlet process?  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 37(2):212–229, 2015.
  * Duchi et al. (2011) John Duchi, Elad Hazan, and Yoram Singer.  Adaptive subgradient methods for online learning and stochastic optimization.  _Journal of Machine Learning Research_ , 12(Jul):2121–2159, 2011.
  * Escobar and West (1995) Michael D Escobar and Mike West.  Bayesian density estimation and inference using mixtures.  _Journal of the american statistical association_ , 90(430):577–588, 1995.
  * Fahrmeir et al. (2007) Ludwig Fahrmeir, Thomas Kneib, Stefan Lang, and Brian Marx.  _Regression_.  Springer, 2007.
  * Ferguson (1973) Thomas S Ferguson.  A Bayesian analysis of some nonparametric problems.  _The Annals of Statistics_ , 1(2):209–230, 1973\.
  * Ferguson (1974) Thomas S Ferguson.  Prior distributions on spaces of probability measures.  _The Annals of Statistics_ , pages 615–629, 1974.
  * Fox et al. (2007) Emily B Fox, Erik B Sudderth, and Alan S Willsky.  Hierarchical dirichlet processes for tracking maneuvering targets.  In _Information Fusion, 2007 10th International Conference on_ , pages 1–8. IEEE, 2007.
  * Fraley and Raftery (2007) Chris Fraley and Adrian E Raftery.  Bayesian regularization for normal mixture estimation and model-based clustering.  _Journal of classification_ , 24(2):155–181, 2007\.
  * Franzén (2006) Jessica Franzén.  nbayesian inference for a mixture model using the gibbs sam’pler.  _MResearch Report_ , 1, 2006.
  * Frigyik et al. (2010) Bela A Frigyik, Amol Kapila, and Maya R Gupta.  Introduction to the dirichlet distribution and related processes. department of electrical engineering, university of washignton.  Technical report, UWEETR-2010-0006, 2010.
  * Frühwirth-Schnatter (2006) Sylvia Frühwirth-Schnatter.  _Finite mixture and Markov switching models_.  Springer Science & Business Media, 2006.
  * Gelman et al. (2013) Andrew Gelman, John B Carlin, Hal S Stern, David B Dunson, Aki Vehtari, and Donald B Rubin.  _Bayesian data analysis_.  CRC press, 2013.
  * Gelman et al. (2014) Andrew Gelman, John B Carlin, Hal S Stern, and Donald B Rubin.  _Bayesian data analysis_ , volume 2.  Chapman & Hall/CRC Boca Raton, FL, USA, 2014.
  * Geman and Geman (1984) Stuart Geman and Donald Geman.  Stochastic relaxation, gibbs distributions, and the bayesian restoration of images.  _IEEE Transactions on pattern analysis and machine intelligence_ , (6):721–741, 1984.
  * Gentle (2007) James E Gentle.  _Matrix algebra: theory, computations, and applications in statistics_.  Springer Science & Business Media, 2007.
  * Gershman and Blei (2012) Samuel J Gershman and David M Blei.  A tutorial on bayesian nonparametric models.  _Journal of Mathematical Psychology_ , 56(1):1–12, 2012.
  * Geyer (2011) Charles Geyer.  Introduction to markov chain monte carlo.  _Handbook of markov chain monte carlo_ , pages 3–48, 2011.
  * Gilks and Wild (1992) Walter R Gilks and Pascal Wild.  Adaptive rejection sampling for gibbs sampling.  _Applied Statistics_ , pages 337–348, 1992.
  * Gnedin and Pitman (2005) Alexander Gnedin and Jim Pitman.  Exchangeable Gibbs partitions and Stirling triangles.  _Zap. Nauchn. Sem. St Peterburg. Otdel. Mat. Inst. Steklov._ , 325:83–102, 2005.
  * Görür and Edward Rasmussen (2010) Dilan Görür and Carl Edward Rasmussen.  Dirichlet process gaussian mixture models: Choice of the base distribution.  _Journal of Computer Science and Technology_ , 25(4):653–664, 2010.
  * Graves (2011) Alex Graves.  Practical variational inference for neural networks.  _Advances in neural information processing systems_ , 24, 2011.
  * Griffiths and Ghahramani (2005) Thomas L Griffiths and Zoubin Ghahramani.  Infinite latent feature models and the Indian buffet process.  In _Advances in Neural Information Processing Systems_ , volume 18, pages 475–482, 2005.
  * Hastings (1970) W Keith Hastings.  Monte carlo sampling methods using markov chains and their applications.  _Biometrika_ , 57(1):97–109, 1970.
  * Hoff (2009) Peter D Hoff.  _A first course in Bayesian statistical methods_.  Springer Science & Business Media, 2009.
  * Hoffman et al. (2013) Matthew D Hoffman, David M Blei, Chong Wang, and John Paisley.  Stochastic variational inference.  _Journal of Machine Learning Research_ , 14(5), 2013.
  * Ishwaran and James (2001) Hemant Ishwaran and Lancelot F James.  Gibbs sampling methods for stick-breaking priors.  _Journal of the American Statistical Association_ , 96(453):161–173, 2001.
  * Ishwaran and James (2002) Hemant Ishwaran and Lancelot F James.  Approximate dirichlet process computing in finite normal mixtures: smoothing and prior information.  _Journal of computational and graphical statistics_ , 11(3):508–532, 2002.
  * Ishwaran and James (2003) Hemant Ishwaran and Lancelot F James.  Generalized weighted chinese restaurant processes for species sampling mixture models.  _Statistica Sinica_ , pages 1211–1235, 2003.
  * Ji et al. (2010) Yangfeng Ji, Tong Lin, and Hongbin Zha.  Cdp mixture models for data clustering.  In _Pattern Recognition (ICPR), 2010 20th International Conference on_ , pages 637–640. IEEE, 2010.
  * Jordan and Bishop (2004) Michael I Jordan and Chris Bishop.  An introduction to graphical models, 2004.
  * Jordan et al. (1999) Michael I Jordan, Zoubin Ghahramani, Tommi S Jaakkola, and Lawrence K Saul.  An introduction to variational methods for graphical models.  _Machine learning_ , 37(2):183–233, 1999.
  * Kamper (2013) Herman Kamper.  Gibbs sampling for fitting finite and infinite gaussian mixture models, 2013.
  * Kruschke (2014) John Kruschke.  Doing bayesian data analysis: A tutorial with r, jags, and stan.  2014\.
  * Kulis and Jordan (2011) Brian Kulis and Michael I Jordan.  Revisiting k-means: New algorithms via bayesian nonparametrics.  _arXiv preprint arXiv:1111.0352_ , 2011.
  * LeCun et al. (2010) Yann LeCun, Corinna Cortes, and Christopher JC Burges.  MNIST handwritten digit database.  _AT &T Labs [Online]. Available: http://yann. lecun. com/exdb/mnist_, 2010.
  * Lian (2010) Heng Lian.  Sparse bayesian hierarchical modeling of high-dimensional clustering problems.  _Journal of Multivariate Analysis_ , 101(7):1728–1737, 2010.
  * Lijoi and Prünster (2010) Antonio Lijoi and Igor Prünster.  Models beyond the Dirichlet process.  In Nils Lid Hjort, Chris Holmes, Peter Müller, and Stephen G Walker, editors, _Bayesian Nonparametrics_ , volume 28, pages 80–136. Cambridge Univ. Press, Cambridge, 2010.
  * Lijoi et al. (2007) Antonio Lijoi, Ramsés H Mena, and Igor Prünster.  Bayesian nonparametric estimation of the probability of discovering new species.  _Biometrika_ , 94(4):769–786, 2007.
  * Lo (2005) Albert Y Lo.  Weighted chinese restaurant processes.  _Cosmos_ , 1(01):107–111, 2005.
  * Lu (2017a) Jun Lu.  Hyperprior on symmetric dirichlet distribution.  _arXiv preprint arXiv:1708.08177_ , 2017a.
  * Lu (2017b) Jun Lu.  Revisit log-convexity of a ratio of gamma functions.  _EPFL_ , 2017b.
  * Lu (2017c) Jun Lu.  Robust model-based clustering for big and complex data.  2017c.
  * Lu (2021a) Jun Lu.  Numerical matrix decomposition and its modern applications: A rigorous first course.  _arXiv preprint arXiv:2107.02579_ , 2021a.
  * Lu (2021b) Jun Lu.  A rigorous introduction for linear models.  _arXiv preprint arXiv:2105.04240_ , 2021b.
  * Lu (2021c) Jun Lu.  A rigorous introduction for linear models.  _arXiv preprint arXiv:2105.04240_ , 2021c.
  * Lu et al. (2018) Jun Lu, Meng Li, and David Dunson.  Reducing over-clustering via the powered chinese restaurant process.  _arXiv preprint arXiv:1802.05392_ , 2018.
  * Ma et al. (2014) Zhanyu Ma, Pravin Kumar Rana, Jalil Taghia, Markus Flierl, and Arne Leijon.  Bayesian estimation of dirichlet mixture model with variational inference.  _Pattern Recognition_ , 47(9):3143–3157, 2014\.
  * Mandt and Blei (2014) Stephan Mandt and David Blei.  Smoothed gradients for stochastic variational inference.  _arXiv preprint arXiv:1406.3650_ , 2014.
  * McCullagh et al. (2008) Peter McCullagh, Jie Yang, et al.  How many clusters?  _Bayesian Analysis_ , 3(1):101–120, 2008.
  * McDaid et al. (2013) Aaron F McDaid, Derek Greene, and Neil Hurley.  Normalized mutual information to evaluate overlapping community finding algorithms.  _arXiv preprint arXiv:1110.2515v2_ , 2013.
  * Meilă (2003) Marina Meilă.  Comparing clusterings by the variation of information.  In Bernhard Schölkopf and Manfred K. Warmuth, editors, _Learning Theory and Kernel Machines_ , pages 173–187. Springer Berlin Heidelberg, 2003.
  * Merkle (1997) Milan Merkle.  On log-convexity of a ratio of gamma functions.  _Publikacije Elektrotehničkog fakulteta. Serija Matematika_ , pages 114–119, 1997.
  * Metropolis et al. (1953) Nicholas Metropolis, Arianna W Rosenbluth, Marshall N Rosenbluth, Augusta H Teller, and Edward Teller.  Equation of state calculations by fast computing machines.  _The journal of chemical physics_ , 21(6):1087–1092, 1953.
  * Miller and Harrison (2013) Jeffrey W Miller and Matthew T Harrison.  A simple example of Dirichlet process mixture inconsistency for the number of components.  In _Advances in Neural Information Processing Systems_ , pages 199–206, 2013.
  * Murphy (2007) Kevin P Murphy.  Conjugate bayesian analysis of the gaussian distribution.  _def_ , 1(2$\sigma$2):16, 2007.
  * Murphy (2012) Kevin P Murphy.  _Machine learning: a probabilistic perspective_.  MIT press, 2012.
  * Neal (2000) Radford M Neal.  Markov chain sampling methods for Dirichlet process mixture models.  _Journal of Computational and Graphical Statistics_ , 9(2):249–265, 2000.
  * Ollivier (2015) Yann Ollivier.  Laplace’s rule of succession in information geometry.  In _International Conference on Geometric Science of Information_ , pages 311–319. Springer, 2015.
  * Perman et al. (1992) Mihael Perman, Jim Pitman, and Marc Yor.  Size-biased sampling of Poisson point processes and excursions.  _Probability Theory and Related Fields_ , 92(1):21–39, 1992.
  * Petralia et al. (2012) Francesca Petralia, Vinayak Rao, and David B Dunson.  Repulsive mixtures.  In _Advances in Neural Information Processing Systems_ , pages 1889–1897, 2012.
  * Ranganath et al. (2014) Rajesh Ranganath, Sean Gerrish, and David Blei.  Black box variational inference.  In _Artificial intelligence and statistics_ , pages 814–822. PMLR, 2014.
  * Rasmussen (1999) Carl Edward Rasmussen.  The infinite gaussian mixture model.  In _Advances in Neural Information Processing Systems_ , volume 12, pages 554–560, 1999.
  * Rasmussen (2003) Carl Edward Rasmussen.  Gaussian processes in machine learning.  In _Summer school on machine learning_ , pages 63–71. Springer, 2003\.
  * Rousseau and Mengersen (2011) Judith Rousseau and Kerrie Mengersen.  Asymptotic behaviour of the posterior distribution in overfitted mixture models.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 73(5):689–710, 2011.
  * Ruder (2016) Sebastian Ruder.  An overview of gradient descent optimization algorithms.  _arXiv preprint arXiv:1609.04747_ , 2016.
  * Seeger (2004) Matthias Seeger.  Low rank updates for the cholesky decomposition.  Technical report, 2004.
  * Shen et al. (2013) Weining Shen, Surya T Tokdar, and Subhashis Ghosal.  Adaptive Bayesian multivariate density estimation with Dirichlet mixtures.  _Biometrika_ , 100(3):623–640, 2013.
  * Socher et al. (2011) Richard Socher, Andrew L Maas, and Christopher D Manning.  Spectral Chinese restaurant processes: Nonparametric clustering based on similarities.  In _Fourteenth International Conference on Artificial Intelligence and Statistics (AISTATS)_ , pages 698–706, 2011.
  * Stephens (1997) Matthew Stephens.  Bayesian methods for mixtures of normal distributions, 1997.
  * Teh (2007) Yee Whye Teh.  Exponential families: Gaussian, gaussian-gamma, gaussian-wishart, multinomial, 2007.
  * Teh (2011) Yee Whye Teh.  Dirichlet process.  In _Encyclopedia of Machine Learning_ , pages 280–287. Springer, 2011\.
  * Tierney (1998) Luke Tierney.  A note on metropolis-hastings kernels for general state spaces.  _Annals of applied probability_ , pages 1–9, 1998.
  * Turchin (1971) Valentin F Turchin.  On the computation of multidimensional integrals by the monte-carlo method.  _Theory of Probability & Its Applications_, 16(4):720–724, 1971.
  * van Havre et al. (2015) Zoé van Havre, Nicole White, Judith Rousseau, and Kerrie Mengersen.  Overfitting Bayesian mixture models with an unknown number of components.  _PloS one_ , 10(7):e0131739, 2015.
  * West (1992) Mike West.  _Hyperparameter estimation in Dirichlet process mixture models_.  Duke University ISDS Discussion Paper# 92-A03, 1992.
  * West and Escobar (1993) Mike West and Michael D Escobar.  _Hierarchical priors and mixture models, with application in regression and density estimation_.  Institute of Statistics and Decision Sciences, Duke University, 1993.
  * Wild and Gilks (1993) Pascal Wild and WR Gilks.  Algorithm as 287: Adaptive rejection sampling from log-concave density functions.  _Journal of the Royal Statistical Society. Series C (Applied Statistics)_ , 42(4):701–709, 1993.
  * Xu et al. (2016) Yanxun Xu, Peter Müller, and Donatello Telesca.  Bayesian inference for latent biologic structure with determinantal point processes (dpp).  _Biometrics_ , 2016.
  * Yang and Dunson (2014) Yun Yang and David B Dunson.  Minimax optimal bayesian aggregation.  _arXiv preprint arXiv:1403.1345_ , 2014.
  * Zeiler (2012) Matthew D Zeiler.  Adadelta: an adaptive learning rate method.  _arXiv preprint arXiv:1212.5701_ , 2012.