# Easily Computed Marginal Likelihoods for Multivariate Mixture Models Using the THAMES Estimator

Martin Metodiev, Nicholas J Irons, Marie Perrot-Dockès, Pierre Latouche, Adrian E Raftery

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Martin Metodiev, Nicholas J Irons, Marie Perrot-Dockès, Pierre Latouche, Adrian E Raftery. Easily Computed Marginal Likelihoods for Multivariate Mixture Models Using the THAMES Estimator. 2025. hal-05053852v2

HAL Id: hal-05053852

https://hal.science/hal-05053852v2

Preprint submitted on 5 Jul 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Easily Computed Marginal Likelihoods for Multivariate Mixture Models Using the THAMES Estimator

Martin Metodiev* ◦, Nicholas J. Irons† ◦, Marie Perrot-Dockès‡, Pierre Latouche*,§ and Adrian E. Raftery¶

**Abstract.** We present a new version of the truncated harmonic mean estimator (THAMES) for univariate or multivariate mixture models. The estimator computes the marginal likelihood from Markov chain Monte Carlo (MCMC) samples, is consistent, asymptotically normal and of finite variance. In addition, it is invariant to label switching, does not require posterior samples from hidden allocation vectors, and is easily approximated, even for an arbitrarily high number of components. Its computational efficiency is based on an asymptotically optimal ordering of the parameter space, which can in turn be used to provide useful visualisations. We test it in simulation settings where the true marginal likelihood is available analytically. It performs well against state-of-the-art competitors, even in multivariate settings with a high number of components. We demonstrate its utility for inference and model selection on univariate and multivariate data sets.

**MSC2020 subject classifications:** Primary 62F15, 62-04; secondary 62F12.

**Keywords:** Bayesian statistics, bridge sampling, harmonic mean estimator, label switching, marginal likelihood estimation, mixture modelling.

## 1 Introduction

Choosing the right number of components for a mixture model is a topic of much discussion [see 16, 9, for reviews]. The marginal likelihood, also known as the evidence or the normalising constant, is a Bayesian model choice criterion that can be used for this task and that is optimal from a Bayesian point of view, since it minimises the model selection error rate on average over the prior distribution [31].

Unfortunately, the marginal likelihood often does not possess an analytic expression. A variety of algorithms have been proposed for its computation, none of which is clearly seen as a gold standard [see 37, for a review]. This is why the truncated harmonic mean estimator [THAMES, 41] was conceived, with the authors arguing that it is the only marginal likelihood estimator that meets the following three desiderata: it is precise,

<footnote>arXiv: 2504.21812</footnote>

<footnote>*Laboratoire de Mathématiques Blaise Pascal, Université Clermont Auvergne, mar tin.metodiev@doctorant.uca.fr</footnote>

<footnote>†Department of Statistics and Leverhulme Centre for Demographic Science, University of Oxford, nicholas.irons@stats.ox.ac.uk</footnote>

<footnote>‡CNRS, MAP5, Université Paris Cité, marie.perrot-dockees@u-paris.fr</footnote>

<footnote>§Institut universitaire de France, pierre.latouche@uca.fr</footnote>

<footnote>¶Departments of Statistics and Sociology, University of Washington, raftery@uw.edu</footnote>
generic and simple. However, it is not adequate for mixture models. In addition, there
is an almost complete lack in the literature of marginal likelihood estimation techniques
for mixture models that are multivariate and contain a high number of components.
The goal of this article is to present an adaptation of the THAMES that can be used
in this setting.

**Marginal likelihood estimation** Available approaches for marginal likelihood estima-
tion do not satisfy the three desiderata described before. These include bridge sampling
[38] (when using a generic, non-symmetric proposal such as the bridge sampling imple-
mentation of Gronau et al. [25]) and Chib's estimator [11], which have been observed to
be imprecise when used to estimate the marginal likelihood of mixture models [45, 9].
Many precise alternatives are not generic: they either become intractable as the number
of mixture components increases [e.g., the proposals of 6, 34], or require specific Markov
chain Monte Carlo (MCMC) samplers that can obtain samples from hidden allocation
vectors [e.g., the proposals of 47, 54, 50, 28]. A truly generic estimator would not have
to depend on these restrictions. Indeed, there are real-world examples of mixture models
with a substantial number of components. For example, in Bouveyron et al. [7] a dataset
is analysed for which the true number of components is estimated to be 15. There are
also many popular samplers that do not provide samples from the hidden allocation vec-
tors, such as the Hamiltonian Monte Carlo (HMC) algorithm implemented in Stan [58],
the classical Metropolis-Hastings algorithm [42, 29] and the tempered MCMC algorithm
[44].

Finally, a simple estimator has to avoid the label switching problem without simu-
lating from randomly permuted MCMC samples, because many parameter estimators
of interest, such as the arithmetic mean, will not work on these samples. In this paper,
label-switched samples are bypassed entirely by post-process relabelling algorithms such
as Stephens’ algorithm [60] and the equivalence classes representatives (ECR) algorithm
[49]. Estimators that require simulating from randomly permuted MCMC samples [e.g.,
the proposals of 15, 17] are not simple by our definition.

The version of the THAMES that was introduced in Metodiev et al. [41] is not
appropriate in the context of mixture models. This is because it is optimal only in the
case that the posterior distribution of the parameters is asymptotically normal. While
this case is important, it does not occur in mixture models due to the highly irregular
and often non-identifiable structure of the posterior distribution.

**Our proposal** Building on the ideas of Berkhof et al. [6] and Reichl [51], we propose an
adapted version of the THAMES that can be used in multivariate mixture models with
a high number of components, while still adhering to the three desiderata of precision,
generality and simplicity. It is precise, since it is consistent, asymptotically normal, and
optimal in a certain sense, even for the multimodal posterior densities that typically
appear in mixture models. It is also generic, since it can be easily approximated when
the number of components is large, and it is compatible with any MCMC algorithm that
provides a sample from the posterior distribution of the parameters. Most importantly,
it is simple, since it is a symmetric truncated harmonic mean of the unnormalised
posterior density values for a relabelled sample from the posterior parameters.
**Outline of the paper** The paper is organised as follows. In Section 2, we describe previous work on which the THAMES for mixture models is built. Then, we introduce the THAMES for mixture models in Section 3.1 and explain how it can be computed efficiently in Section 3.2, using the principle of overlap graphs to obtain an ordering constraint. In addition, we propose a model choice criterion, the criterion of overlap (CO), and show how to implement an identity from Nobile [46, 47] in Section 3.3 to better deal with the phenomenon of empty components. We compare the THAMES to state-of-the-art estimators in simulation studies, where the true marginal likelihood is known, in Section 4.1, and on real datasets, in Section 4.2. We conclude with a discussion in Section 5.

## 2 Previous work relating to the THAMES

Let $\theta$ denote the model parameter of dimension $R$ and $\mathcal{D}$ the data. The Bayesian framework is considered, where $\pi(\theta)$ is the prior probability density function (PDF) and $L(\theta) = p(\mathcal{D}|\theta)$ is the likelihood.

Let $\theta^{(1)}, \dots, \theta^{(T)}$ be a sample from the posterior distribution $p(\theta|\mathcal{D})$. The truncated harmonic mean estimator [THAMES, 41] of the reciprocal marginal likelihood is defined as

$$ \hat{Z}^{-1} = \frac{1}{T/2} \sum_{\substack{t=T/2+1, \\ \theta^{(t)} \in A}}^{T} \frac{1/V(A)}{\pi(\theta^{(t)})L(\theta^{(t)})}, \qquad (1) $$

where $A$ is the set

$$ A = E_{\hat{\theta}, \hat{\Sigma}, c} = \{\theta | (\theta - \hat{\theta})^{\mathrm{T}} \hat{\Sigma}^{-1} (\theta - \hat{\theta}) < c^2\}, \qquad (2) $$

an ellipsoid centered at an estimate of the posterior mean $\hat{\theta}$ and scaled by an estimate of the posterior covariance matrix $\hat{\Sigma}$. Notice that the posterior sample is split: $\hat{\theta}$ and $\hat{\Sigma}$ are computed using the first half (or portion, more generally), while the THAMES is computed from the second half. Here $V(A)$ denotes the volume of $A$.

Any choice of $A$ within the posterior support leads to an unbiased estimator on the scale of the reciprocal marginal likelihood. Under mild assumptions, an optimal truncation set is given by an $\alpha$-highest-posterior-density (HPD) region

$$ H_{\alpha} = \{\theta | \pi(\theta) L(\theta) > q_{\alpha}\}, $$

where the quantile $q_{\alpha}$ is chosen such that the posterior probability of $H_{\alpha}$ is equal to the tuning-parameter $\alpha$ [for a proof of this statement, see 40, Supplement A]. Unfortunately, the volume of $H_{\alpha}$ is needed for the computation of the THAMES. We know of no simple way to compute this volume in general. Robert et al. [53] suggested using a convex hull which contains $H_{\alpha}$, but approximating the volume of a convex hull is difficult when the dimension of the parameter space is large. This is avoided by the proposal of Reichl [51],
who suggested to instead use a set that is sufficiently close to $H_\alpha$, but whose volume is easy to compute, by taking the intersection of the ellipsoid $E_{\hat{\theta}, \hat{\Sigma}, c}$ with $H_\alpha$:

$$ 
\begin{aligned}
B_{\hat{\theta}, \hat{\Sigma}, c, \alpha} &= E_{\hat{\theta}, \hat{\Sigma}, c} \cap \{\theta | \pi(\theta) L(\theta) > \hat{q}_{\alpha}\} \\
&= \{\theta | (\theta - \hat{\theta})^{\top} \hat{\Sigma}^{-1} (\theta - \hat{\theta}) < c^2, \pi(\theta) L(\theta) > \hat{q}_{\alpha}\},
\end{aligned}
$$

where $\hat{q}_\alpha$ is the empirical $(1-\alpha)$-quantile of the sample of unnormalised log-posterior values.

This proposal requires additional Monte Carlo simulations to compute the volume of $B_{\hat{\theta}, \hat{\Sigma}, c, \alpha}$. This is done by using a technique that was successfully employed by Sims et al. [57] in reciprocal importance sampling: an i.i.d sample $\nu^{(1)}, \dots, \nu^{(N)}$ from the uniform distribution on $E_{\hat{\theta}, \hat{\Sigma}, c}$ is simulated to calculate

$$ V(B_{\hat{\theta}, \hat{\Sigma}, c, \alpha}) \approx V(E_{\hat{\theta}, \hat{\Sigma}, c}) \cdot \frac{1}{N} \sum_{j=1}^{N} 1_{\{\theta | \pi(\theta) L(\theta) > \hat{q}_{\alpha}\}} (\nu^{(j)}). \quad (3) $$

The original version of the THAMES did not require additional Monte Carlo simulations in many cases. However, the estimator of Reichl [51] is arguably a necessary adaptation in the case of mixture models, since it is presumably more robust in the face of irregularities in the posterior distribution. Reichl [51] implemented this estimator and obtained good results in a variety of situations. In particular, in the context of mixture models, it performed well, but only after adding additional MCMC simulations from random permutation sampling via a technique from Frühwirth-Schnatter [14]. This runs against our desideratum of simplicity.

# 3 Methods

## 3.1 Truncated harmonic mean estimation for mixture models

The standard mixture model setting is considered, where the data $D = Y \in \mathcal{M}_{n \times d}(\mathbb{R})$ consist of an i.i.d sample of $n$ different $d$-dimensional observations, $Y = (Y_1, \dots, Y_n)$, with a distribution that is a mixture of distributions from the same family,

$$ Y_i | \xi, \tau \stackrel{\text{i.i.d}}{\sim} \sum_{g=1}^{G} \tau_g f(\cdot | \xi_g), \quad i = 1, \dots, n, \quad (4) $$

where $G > 1$, $\xi = (\xi_1, \dots, \xi_G)$ is a matrix containing the mixture component parameters, $\tau_1, \dots, \tau_G > 0$, $\sum_{g=1}^G \tau_g = 1$ are the mixture proportions represented by the vector $\tau = (\tau_1, \dots, \tau_G)$, and $f(\cdot | \xi_g)$ is a density characterised by $\xi_g$. We are interested in estimating the marginal likelihood $Z(G)$ for different values of the number of mixture components $G$. The full parameter is denoted by $\theta = (\xi_1, \dots, \xi_G, \tau_1, \dots, \tau_{G-1})$.

If the prior distribution of $(\tau_g, \xi_g)$ does not depend on the index $g$, which is commonly the case and an assumption that we make from now on, the posterior distribution is
symmetric. This results in the well-known label-switching problem, where some parts of the MCMC sample randomly switch their labels. Luckily, it is possible to simulate a sample whose labels are not switched by relabelling, e.g., with Stephens' algorithm [60] or the ECR algorithm [49]. In the former, the sample is permuted such that it minimises the posterior expected loss under a class of loss functions, while the latter is based on permuting the sample such that it follows the law of a non-symmetric distribution. The two algorithms give similar results in our experience.

Let $\theta^{(1)*}, \dots, \theta^{(T)*}$ denote the posterior sample after relabelling. A marginal likelihood estimator can be based entirely on this sample if it is symmetric, i.e., if the labels of the parameters do not change the value of the estimator [for a proof, see 40, Supplement A]. By a procedure suggested in Berkhof et al. [6], a symmetric estimator can be obtained by averaging a non-symmetric estimator over the permutations $P_1, \dots, P_G$ of the posterior sample. Hence, we suggest combining the proposal of Reichl [51] with the symmetrisation procedure used by Berkhof et al. [6] to obtain the symmetric (S) THAMES for mixture models,

$$ \hat{Z}_S^{-1}(G) = \frac{1}{G!} \sum_{o=1}^{G!} \frac{1}{T/2} \sum_{\substack{t=T/2+1, \\ P_o(\theta^{(t)*}) \in B_{\hat{\theta}, \hat{\Sigma}, c, \alpha}}} \frac{1/V(B_{\hat{\theta}, \hat{\Sigma}, c, \alpha})}{\pi(\theta^{(t)*})L(\theta^{(t)*})} \quad (5) $$

We compute the volume of $B_{\hat{\theta}, \hat{\Sigma}, c, \alpha}$ by setting $N = T$ in Equation (3). Note that a naive computation of the full sum over the $G!$ permutations is not necessary, since most permutations will evaluate to 0. We elaborate on this in the next section.

The symmetric THAMES (hereafter just THAMES) possesses the same attractive properties as the original proposal of Metodiev et al. [41]: it is unbiased on the scale of the reciprocal marginal likelihood, consistent, and asymptotically normal as the size $T$ of the simulations increases [for a proof, see 40, Supplement A]. Equation (5) also directly implies that this estimator is symmetric, and we show in Metodiev et al. [Supplement A, 40] that this symmetrisation is in some sense optimal over all mixtures of truncated harmonic mean estimators. Only the relabelled posterior sample is required for its computation, with the first half being used to estimate the posterior mean and covariance matrix. It possesses two tuning-parameters ($\alpha$ and $c$) and we illustrate in Metodiev et al. [Supplement B, 40] how these can be chosen in an optimal way. Finally, we show in the next section that the THAMES for mixture models can be computed quickly, even as the number of components $G$ grows large.

## 3.2 Calculating the THAMES for mixture models efficiently

A naive computation of Equation (5) would require $G!$ computation steps, one for each permutation of the $G$ different labels. This rapidly becomes infeasible computationally as $G$ increases. However, in practice, many of these computations can be avoided due to the fact that the support $B_{\hat{\theta}, \hat{\Sigma}, c, \alpha}$ is contained within a very simple geometric object, the ellipsoid $E_{\hat{\theta}, \hat{\Sigma}, c}$. We can hence reduce computation time by selecting a restricted number of components and only summing over those permutations which respect the
<figure><img src="image_3.png" /><figcaption>Figure 1: Example of two overlap graphs; an edge is added whenever there is overlap between 2 components.</figcaption></figure>

order of these components within the ellipsoid $E_{\hat{\theta}, \hat{\Sigma}, c}$. The selection of the components is done via graph-theoretic tools, while the permutations are determined via quadratic discriminant analysis (QDA; Ghojogh and Crowley 2019).

In the next three subsections, we propose a way to calculate the THAMES more efficiently. First, we apply the principle of “overlap graphs” to choose a maximum set of components that do not overlap. Then, we use this set to construct an ordering constraint via QDA. We use this ordering constraint to compute the THAMES efficiently by limiting the number of permutations required in the computation. We emphasise that overlapping components are determined in the parameter space and not in the observation space, since the objective is to reduce the computation cost of Equation (5), which is computed on the parameters from the posterior.

## Visualising the posterior distribution via overlap graphs

Since the radius $c = \sqrt{R+1}$ is taken as suggested by Metodiev et al. [41], $E_{\hat{\theta}, \hat{\Sigma}, c}$ roughly estimates a 50%-HPD-region on the relabelled posterior distribution. For more details, we refer the reader to Metodiev et al. [41]. In such a region, two components should not be overlapping if none of the components is superfluous. It can easily be verified [see 40, Supplement B] if the hyperplane $\{\theta|\xi_{g1} = \xi_{g2}\}$ crosses the ellipsoid $E_{\hat{\theta}, \hat{\Sigma}, c}$. If this is the case, we say that there is overlap between the components $g1$ and $g2$. This could be an indicator of either $g1$ or $g2$ being superfluous since they could be the same component.

Figure 1 shows two “overlap graphs”, where each node is a component and an edge is added between two nodes if there is overlap between the corresponding components. In Model 1, there is no overlap between all components, so there are no edges present. In Model 2, there is overlap between the pairs (1,2), (2,4), (4,3) and (3,1), so there are at most two non-overlapping components present.

Using graph-theoretic techniques on an overlap graph, we can find a maximum independent set, i.e., a maximum set of non-overlapping components. In Model 2 of Figure 1 such a set is, for instance, given by components 1 and 4. This set can be determined via the `buildCoverGreedy` function in the `gor` package [2] and is denoted by $I(G)$. Note that the complexity of `buildCoverGreedy` is quadratic with respect to $G$. As we
are going to see, I(G) will strongly help reducing the computational cost of Equation (5). Interestingly, the size |I(G)| can also be interpreted as an estimate of the number of distinct (i.e., non-overlapping) components. Thus, we define the criterion of overlap (CO) as

$$CO = \#_{\text{components}}^{\text{distinct}} - \#_{\text{components}}^{\text{overlapping}} = |I(G)| - (G - |I(G)|).$$

The CO is related to the goal of determining the number of distinct components (i.e., clusters) instead of the number of true components. The difference between these two different goals was pointed out in Baudry et al. [4]. While determining the number of components G that maximises the marginal likelihood is useful in the latter case, choosing G such that it maximises the CO is useful in the former. The set I(G) that determines the CO is also going to be used to define an ordering constraint to estimate the marginal likelihood: nodes outside of I(G) cannot be ordered together with nodes inside I(G) and should thus not be used directly to define the constraint. The constraint will be constructed via QDA, and will in turn be used to quickly compute the THAMES.

## Applying QDA to the posterior sample

Our aim is to find a classification of a component parameter $\xi_g$ such that the labels of two component parameters differ if, and only if, their classifications differ. This is done using discriminant analysis.

Focusing on the mixture component parameters in $\xi$, discriminant analysis can be performed on the relabelled posterior sample because the labels of the parameters are known by construction. Indeed, $\xi_1$ are the parameters of component 1, $\xi_2$ of component 2, and so on. Thus, a labelled sequence of size $G \cdot T$ from the relabelled posterior sample is given by

$$\xi_1^{(1)*}, \dots, \xi_G^{(1)*}, \dots, \xi_1^{(T)*}, \dots, \xi_G^{(T)*}.$$

In the following, QDA is applied to this sequence. The idea is that QDA should be able to distinguish two samples from two component parameters if, and only if, they are concentrated around two different values.

Let $\hat{\xi}_1, \dots, \hat{\xi}_G$ and $\hat{S}_1, \dots, \hat{S}_G$ denote the estimated component means and covariance matrices on the second half of the posterior sample, respectively. Let MVN$_u(\xi_{g'}; \hat{\xi}_g, \hat{S}_g)$ denote the multivariate normal distribution of dimension $u$ with mean $\hat{\xi}_g$ and covariance matrix $\hat{S}_g$, evaluated at $\xi_{g'}$. Then, any value $\xi_{g'} \in \mathbb{R}^u$ (of which the true component $g'$ is unknown) will be assigned to

$$\hat{g}(\xi_{g'}) = \operatorname*{argmax}_g \hat{w}_g(\xi_{g'}), \quad \hat{w}_g(\xi_{g'}) = \frac{\frac{1}{G} \operatorname{MVN}_u(\xi_{g'}; \hat{\xi}_g, \hat{S}_g)}{\sum_{\hat{g}=1}^{G} \frac{1}{G} \operatorname{MVN}_u(\xi_{g'}; \hat{\xi}_g, \hat{S}_g)}. \quad (6)$$

Note that the cluster proportions are always estimated as $\frac{1}{G}$ by construction. Thus, $\xi_{g'}^{(t)*}$ is assigned the wrong label if $\hat{g}(\xi_{g'}^{(t)*}) \neq g'$.
It should also be noted that, if G is the true number of components, the relabelled posterior distribution should converge to the true value as the size of the data increases. Moreover, the condition needed to apply QDA is asymptotic normality, and posteriors are often asymptotically normal when converging to the truth [30, 22, 56, 43]. Since this condition is often met, we argue that an asymptotically optimal ordering constraint can be defined via QDA.

Defining the ordering constraint

In the literature, ordering constraints have been employed to assure identifiability [see for example 52, 59, 35, 20]. However, the specific choice of the constraint can be of crucial importance [10]. Results from the previous two sections will be applied to construct a constraint that is useful for our purposes, in the sense that the constraint will mostly hold within an approximate 50%-HPD-region.

An approximate 50%-HPD-region was used to define the set I(G) such that a constraint could hold on this region for all of the components of I(G). We now restrict the set of components over which the QDA allocation probabilities are computed. Thus, conditioning the QDA output on I(G) gives the probability $\hat{w}_g(\xi_{g'}|I(G)) = \frac{\hat{w}_g(\xi_{g'})}{\sum_{\bar{g} \in I(G)} \hat{w}_{\bar{g}}(\xi_{g'})}$ for all $g \in I(G), g' \in \{1, \dots, G\}$.

Given normality, the assignment $\hat{g}$ and the uncertainty assessment $\hat{w}$ meet the optimality results provided by QDA [21]. With this in mind, we combine this assignment $\hat{g}(\xi_{g'}|I(G))$ with its probability $\hat{w}_{\hat{g}(\xi_{g'}|I(G))}(\xi_{g'}|I(G))$ to obtain the summary measure $W(\xi_{g'})$:

$$W(\xi_{g'}) = \hat{g}(\xi_{g'}|I(G)) + 1 - \hat{w}_{\hat{g}(\xi_{g'}|I(G))}(\xi_{g'}|I(G)).$$

W distinguishes the components in I(G) by construction: if $\xi_{g'}^{(t)*}$ is assigned to the component $\hat{g}(\xi_{g'}^{(t)*}|I(G))$ by QDA, it is placed within $[\hat{g}(\xi_{g'}^{(t)*}|I(G)), \hat{g}(\xi_{g'}^{(t)*}|I(G)) + 1]$ by W. It is moved further away from $\hat{g}(\xi_{g'}^{(t)*}|I(G))$ if its placement was assigned with low certainty.

W is used to define the subspace

$$\{\xi_1, \dots, \xi_G : W(\xi_1) \le \dots \le W(\xi_G)\}. \qquad (7)$$

This subspace will be used to subdivide the parameter space in G! different subspaces $J_1, \dots, J_G$ created by permuting the ordering constraint, e.g., $J_1 = \{\xi_1, \dots, \xi_G : W(\xi_1) \le W(\xi_2) \le \dots \le W(\xi_G)\}, J_2 = \{\xi_1, \dots, \xi_G : W(\xi_2) \le W(\xi_1) \le \dots \le W(\xi_G)\}$ and so on. This subdivision will be used for the computation of the THAMES in Equation (5).

Easily computing the THAMES

Let $\psi_1, \dots, \psi_G$ denote the G! different ordering algorithms that map $\theta$ to the different permutations $J_1, \dots, J_G$ of the subspace defined in Equation (7), respectively, by permuting the labels of $\theta$. Note that these are different from $P_1, \dots, P_G$, since $\psi_1, \dots, \psi_G$
<figure><img src="image_4.png" /><figcaption>Figure 2: A toy example to illustrate the computation of the THAMES; the points that lie in the ellipsoid can be counted by applying $P_1, P_2$ or $\psi_1, \psi_2$.</figcaption></figure>

always enforce some version of Equation (7) independent of the value of $(\xi_1, \dots, \xi_G)$. For example, $\psi_1$ may order $(\xi_1, \dots, \xi_G)$ to $(\xi_1^{\psi_1}, \dots, \xi_G^{\psi_1})$ such that $W(\xi_1^{\psi_1}) \le W(\xi_2^{\psi_1}) \le \dots \le W(\xi_G^{\psi_1})$ for all values of $(\xi_1, \dots, \xi_G)$, while $\psi_2$ would result into $W(\xi_2^{\psi_2}) \le W(\xi_1^{\psi_2}) \le \dots \le W(\xi_G^{\psi_2})$ for all values of $(\xi_1, \dots, \xi_G)$.

The difference between $\psi_1, \dots, \psi_G!$ and $P_1, \dots, P_G!$ is shown in Figure 8 in the case that $G=2, T=10$ and that there are only two parameters $\theta^{(t)\star} = (\mu_1^{(t)\star}, \mu_2^{(t)\star})$, one for each component, with $W(\mu_g) = \mu_g$. Computing Equation (5) boils down to counting the number of times that $\psi_o(\theta^{(t)\star}), P_o(\theta^{(t)\star})$ falls in the ellipse for any $o=1, \dots, G!$. One can either apply $P_1$ (the identity) and $P_2$ (flipping the x- and y-axis), or $\psi_1$ (ordering the points such that $\mu_1^{(t)\star} \le \mu_2^{(t)\star}$ for all $t=1, \dots, T$) and $\psi_2$ (ordering the points such that $\mu_2^{(t)\star} \le \mu_1^{(t)\star}$ for all $t=1, \dots, T$). The outcome will be the same: 9 points are counted to be in the ellipse for some $\psi_o(\theta^{(t)\star}), P_o(\theta^{(t)\star})$, but computing $\psi_2(\theta^{(t)\star})$ could be avoided entirely since the ellipse does not cross the identity line.

We are going to determine which of $\psi_1, \dots, \psi_G!$ to evaluate as follows: let $\Delta$ denote the adjacency matrix defined by $\Delta_{g_1,g_2} = 1$ if $W(\xi_{g_1}) < W(\xi_{g_2})$ for all $\xi_{g_1}, \xi_{g_2}$ that satisfy the ellipsoidal constraint on $E_{\hat{\theta}, \hat{\Sigma}, c}$. In addition, let $\Omega$ denote the index set over the collection of all topological orderings on the graph defined by $\Delta$. In other words, $o \in \Omega$ if the previously selected inequalities $W(\xi_{g_1}^{\psi_o}) < W(\xi_{g_2}^{\psi_o})$ hold within the support of the THAMES. Thus, the THAMES from Equation (5) is equal to

$$ \hat{Z}_S^{-1}(G) = \frac{1}{G!} \sum_{o \in \Omega} \frac{1}{T/2} \sum_{\substack{t=T/2+1 \\ \psi_o(\theta^{(t)\star}) \in B_{\hat{\theta}, \hat{\Sigma}, c, \alpha}}} \frac{1/V(B_{\hat{\theta}, \hat{\Sigma}, c, \alpha})}{\pi(\theta^{(t)\star})L(\theta^{(t)\star})}. \quad (8) $$

The idea is that we only need to compute those orderings $\psi_o$ whose output does not contradict the constraints that hold in $B_{\hat{\theta}, \hat{\Sigma}, c, \alpha}$. All other orderings would be 0 anyway. For a proof, see Metodiev et al. [Supplement B, 40].

The inequalities $W(\xi_{g_1}) < W(\xi_{g_2})$ can be verified approximately by re-using the Monte Carlo sample from Equation (3). See Metodiev et al. [Supplement B, 40] for details on this computation and the computation of $\Omega$, which can be done in polynomial time with respect to the size of $\Omega$. Note that any measurable function $W$ could have
been used to construct $\Omega$. $W$ was chosen by QDA to target a function for which the inequalities $W(\xi_{g_1}) < W(\xi_{g_2})$ are fulfilled for as many pairs $(g_1, g_2)$ as possible, thus reducing the size of $\Omega$.

The parameters $c$ and $\alpha$ can always be chosen such that the set $\Omega$ is small enough and that the THAMES can be computed in a reasonable amount of time through Equation (17): the smaller the chosen ellipsoid is, the more the overlap decreases. See Metodiev et al. [Supplement B, 40] for more details. In practice we found that this is rarely necessary, since few of the different orderings need to be evaluated. This is likely due to the fact that the components overlap less and less as the size of the data $n$ increases.

## 3.3 The THAMES for high dimensional mixture models in the case of empty components

Aside from the computational challenges already discussed, employing mixture models in high dimensional settings can pose further difficulties. In particular, when the number of components is large, we run the risk of overfitting the data, resulting in empty components assigned to no data points. The posterior distribution of these empty component parameters is often highly diffuse which makes it difficult to fit an ellipsoid around them, a technique that is vital for the precision of the THAMES. In this section, we will show how this problem can be sidestepped.

The component $g$ of a posterior sample entry $\theta^{(t)}$ is identified as possibly being empty if the product of conditional proportions

$$ \prod_{i=1}^{n} (1 - \hat{z}_{i,g}^{(t)}), \text{ where } \hat{z}_{i,g}^{(t)} = \frac{\tau_{g}^{(t)} f(Y_{i}; \xi_{g}^{(t)})}{\sum_{\tilde{g}=1}^{G} \tau_{\tilde{g}}^{(t)} f(Y_{i}; \xi_{\tilde{g}}^{(t)})}, $$

is not close to zero. This may be true for a large value of the total number of components $G$, but is going to be less and less likely as $G$ decreases. Thus, the empty component eventually disappears when simulating from the posterior of a sufficiently small model. In addition, when the probability of having an empty component is large, it is also easily estimated. Suppose that a Dirichlet prior with parameters $(e_1, \dots, e_G) > 0$ is used for the proportions. We use an identity from Nobile [46, 47] to recursively estimate the THAMES as

$$ \hat{Z}_S(G) = \frac{\Gamma\left(\sum_{g=1}^{G} e_g\right) \Gamma\left(n + \sum_{g=1}^{G-1} e_g\right)}{\Gamma\left(n + \sum_{g=1}^{G} e_g\right) \Gamma\left(\sum_{g=1}^{G-1} e_g\right)} \cdot \hat{Z}_S(G-1) \cdot \frac{1}{\hat{p}_0(G)}, \quad (9) $$

with

$$ \hat{p}_0(G) = \frac{1}{G} \sum_{g=1}^{G} \frac{1}{T} \sum_{t=1}^{T} \prod_{i=1}^{n} (1 - \hat{z}_{i,g}^{(t)}). $$

The validity of this approach is proven in Metodiev et al. [Supplement A, 40]. Thus, we can simply check if $\hat{p}_0(G)$ exceeds some upper threshold (e.g., $1/T$), and reduce the
estimation procedure to the estimation of the marginal likelihood of the lower dimensional model if this is true<fnref n="1" />.

# 4 Experiments

We now demonstrate the utility of the THAMES for mixture models on simulated and real datasets. For conciseness, we limit the model descriptions to the definition of the mixture models and the definition of the datasets only. More detailed descriptions, such as the exact prior distributions used, are presented in Metodiev et al. [Supplement C, 40].

## 4.1 Simulations

We perform simulations for two settings in which the marginal likelihood can be either computed analytically or easily approximated with a high degree of accuracy. Formulas for these exact or approximate solutions are given in Metodiev et al. [Supplement D, 40]. In the first, the exact expression of the marginal likelihood is computed in a brute-force manner, which can be done due to the low number of data points. In the second, the marginal likelihood is computed for a high number of extremely distinct components.

### Univariate Gaussian mixtures

**Model** The mixture model is univariate Gaussian with $n = 10$ data points,

$$Y \in \mathbb{R}^n, \quad Y_1, \dots, Y_n | \mu, \sigma, \tau \stackrel{\text{i.i.d}}{\sim} \sum_{g=1}^{G} \tau_g N(\mu_g, \sigma_g^2), \quad (10)$$

where only the parameters $\mu = (\mu_1, \dots, \mu_G)$ are unknown and $\tau = (\tau_1, \dots, \tau_G)$, $\sigma = (\sigma_1, \dots, \sigma_G)$ are known.

**Settings** Simulations are conducted with the fitted and true $G = 2$, but also when **underfitting**, where the true $G = 3$ but the model is fitted for $G = 2$, and **overfitting**, where the converse is true. Component overlap is varied via the parameter $\rho \in [0, 1]$ defined in Metodiev et al. [Supplement C, 40]. It measures how far the posterior means of the component parameters are away from each other. $\rho = 0$ means very little overlap, while $\rho = 1$ implies high overlap.

**Results** The THAMES is compared to the bridge sampling [38] implementation of Gronau et al. [25] and the estimator from Reichl [51] (using the same tuning parameters as the THAMES) for 50 independently generated datasets, with $T = 10,000$, a burn-in

<footnote n="1">¹It should be noted that this technique can only be used if $\xi_1, \dots, \xi_G$ are a priori independent and their prior distribution does not depend on the total number of components. This was true in all multivariate Gaussian mixture models that we used.</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 3: Boxplots of the error log($\hat{Z}$) – log($Z$) of bridge sampling, the estimator from Reichl [51], and the THAMES for 50 different simulations of the data for increasing levels of overlap $\rho = 0, 1/2, 1$ and 3 different scenarios corresponding to the true model, underfitting and overfitting with $G = 2, G = 3$, respectively.</figcaption></figure>

of 2,000 and using ECR for relabelling. Results are shown in Figure 3. Boxplots of the THAMES overlapped with the truth in all cases. However, the other estimators are imprecise and worsen as $\rho$ decreases.

## Multivariate Gaussian mixture models

**Model** Let $Y_1, \dots, Y_n$ be i.i.d values from a mixture of $G$ multivariate normal distributions, namely

$$Y \in \mathbb{R}^{n \times d}, \quad Y_1, \dots, Y_n | \mu, \Sigma, \tau \stackrel{\text{i.i.d}}{\sim} \sum_{g=1}^{G} \tau_g \text{MVN}_d(\mu_g, \Sigma_g). \quad (11)$$

The settings were each chosen such that they correspond to a real dataset.

## Settings

(a) Setting 1 (moderate $G$): the model is $d = 6$-dimensional with $G = 5$ components and $n = 200$ data points.

(b) Setting 2 (large $G$): the model is $d = 5$-dimensional with $G = 15$ components and $n = 345$ data points. The number of free parameters $R$ is very large in this setting ($R = 314$). As is typical in such high dimensional models [see for example 5], it can be decreased by imposing additional constraints on the component covariance matrices. In this case, they are restricted to be diagonal matrices, resulting in much fewer free parameters ($R = 164$).
<figure><img src="image_6.png" /><figcaption>Figure 4: The error log($\hat{Z}$) – log(Z) of the THAMES and bridge sampling for different values of T; on the left: the unconstrained Gaussian mixture model with $G = 5, d = 6$; on the right: the constrained Gaussian mixture model with $G = 15, d = 5$.</figcaption></figure>

**Results** The results are shown in Figure 4. The THAMES and bridge sampling were each computed on an MCMC sample with a burn-in of 2,000, after relabelling via the ECR algorithm. Their evaluations are shown for different posterior sample sizes. The THAMES converges to the true value, while bridge sampling is off by a factor of $G!$. It should be noted that the strategy of simply shifting a non-symmetric marginal likelihood estimator by $G!$ does not work in general. It would work in this specific, artificial setting, only because all components are well separated. It does not work in others, as pointed out by [45].

## 4.2 Real datasets

### The univariate case

The THAMES was evaluated for the datasets used in Richardson and Green [52], namely the “acidity”, “enzyme”, and “galaxy” dataset, available at Ameijeiras-Alonso et al. [1]. Results of other, related marginal likelihood estimators on the same datasets are available in Celeux et al. [9, Section 2.3.2] and were thus used for comparison.

**Settings** All datasets were fitted to a univariate Gaussian mixture model with hierarchical priors, varying $G$ between 2 and 6.

**Results** The THAMES estimates were compared to those provided by Celeux et al. [9] for other estimators under the same circumstances, namely $T = 10,000$ and a burn-in of 2,000. The results are shown in Figure 5. The other estimators are fully symmetric and random permutation bridge sampling (BSF and BSR), importance sampling (ISF and ISR) and reciprocal importance sampling (RIF and RIR). Out of these, it was argued in Celeux et al. [9] that BSF and ISF performed well. Next to them, the median out of
<figure><img src="image_7.png" /><figcaption>Figure 5: Estimates from fully symmetric and random permutation bridge- (BSF, BSR), reciprocal- (RIF, RIR), importance sampling (ISF, ISR) and THAMES for the enzyme, galaxy and acidity datasets.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: Overlap graphs of the galaxies data set.</figcaption></figure>
20 independent THAMES estimates computed on different MCMC samples is shown. The THAMES is close to ISF and BSF in almost all settings, with slight differences occurring when $G = 6$.

One limitation of the ISF and BSF estimators is that they become intractable as $G$ increases, for example if $G = 15$. This is not the case for the THAMES. To illustrate this, it was computed from $G = 2$ to 15 for the galaxy dataset with $T = 100,000$ and a burn-in of 2,000. It peaked at $G = 6$, while the criterion of overlap (CO) was maximised at $G = 3$. The latter result is visualised with overlap graphs, in which the sharp increase of overlapping components at $G > 3$ is apparent (Figure 6). The difference in the CO and the marginal likelihood in this setting is not surprising and can be interpreted as follows: the number of true components is estimated to be 6, while the number of distinguishable components is estimated to be 3. The "correct" result depends on the interests of the practitioner [see 4].

## The multivariate case

Marginal likelihoods are computed for the Swiss banknote dataset, which contains 100 genuine and 100 counterfeit old-Swiss 1000-franc bank notes. They are also computed for the BUPA liver disorders dataset [36], which contains 345 blood samples from different male individuals.

**Settings** The Swiss banknotes dataset consists of $n = 200$ observations and $d = 6$ variables. The BUPA liver disorders dataset consists of $n = 345$ observations and $d = 5$ variables. Multivariate Gaussian mixtures were fitted to both datasets. The covariance matrices of the BUPA liver disorders dataset were restricted to be diagonal matrices to deal with the high dimension of the parameters.

**Results** We set $T = 10,000$ with a burn-in of 2,000 and apply ECR for relabelling. The THAMES was evaluated from $G = 2$ to $G = 5$ for the Swiss banknotes dataset, with the maximum being obtained at $G = 3$. This is consistent with the literature, where the estimated number of components for the Swiss banknotes dataset is often 2 or 3 [12]. For the BUPA liver disorders dataset, the THAMES was evaluated for $G$ between 2 and 15. It was maximised at $G = 4$, with distinct levels of gamma-glutamyl transpeptidase for each component. For more details, we refer the reader to Metodiev et al. [Supplement C, 40].

# 5 Discussion

We have presented an adapted version of the THAMES for mixture models. It is consistent, asymptotically normal, symmetric, and is computed on the relabelled MCMC parameters, without the need of simulations from the hidden allocation vectors. Its simplicity is exemplified by the fact that it only requires two inputs: the unnormalized log-posterior density (i.e., the logarithm of the prior times the likelihood) and the
relabelled MCMC sample from the posterior. A function that takes exactly those two
arguments is implemented in the R-package thamesmix [39].

The THAMES uses the symmetrisation method that was used in Berkhof et al. [6] and is tractable in the number of components, unlike other estimators that use this method. It should be noted that adjustments presented in Berkhof et al. [6] and Lee and Robert [34] already decrease the computation time of symmetrised estimators, but they do not lower it below $O(G!)$, which is the case for the THAMES. For example, we only needed to evaluate less than $10^{-6}$ percent of the permutations over which we sum in the case that $G = 15$ for the galaxy dataset. However, for large $G$, higher overlap implies that the ellipsoid onto which the THAMES is defined will have to be shrunk to ensure that the THAMES can be computed in a reasonable amount of time, thus increasing the variance of the THAMES.

Finally, we want to emphasise that our estimator can be applied to very general high dimensional and multivariate mixture models. This opens up a new avenue of numerous applications. One interesting application may be to further adapt the THAMES for the stochastic block model [62, 48], in which only the complete data likelihood is available.

## Competing interests

No competing interest is declared.

## Data availability

All code is made available at the following link. The acidity, enzyme and galaxy datasets are available at Ameijeiras-Alonso et al. [1], the Swiss banknote dataset is available at Scrucca et al. [55], the BUPA dataset is available at https://doi.org/10.24432/C54G67.

## Acknowledgments

The authors would like to thank Qing Mai from Florida State University for the rec-
ommendation of QDA, as well as Alexander Reisach, Léo Hahn Lecler and Michael
Schaller for fruitful discussions on partial orderings. In addition, they would like to
thank Christian Robert for helpful comments. Raftery's research was supported by NIH
grant R01 HD-070936, the Fondation des Sciences Mathématiques de Paris (FSMP),
Université Paris Cité (UPC), and the Blumstein-Jordan professorship at the University
of Washington. Irons's research was supported by a Shanahan Endowment Fellowship,
an NICHD training grant, T32 HD101442-01, to the Center for Studies in Demography
and Ecology at the University of Washington, the Florence Nightingale Bicentenary
Fellowship in Computational Statistics and Machine Learning from the University of
Oxford Department of Statistics and the Leverhulme Centre for Demographic Science
and the Leverhulme Trust (Grant RC-2018-003). Latouche's research was supported by
the Institut Universitaire de France (IUT).
References

[1] Ameijeiras-Alonso, J., R. M. Crujeiras, and A. Rodríguez-Casal (2021). multimode: An R package for mode assessment. *Journal of Statistical Software* 97(9), 1–32. 13, 16

[2] Asensio, C. (2023). *gor: Algorithms for the Subject Graphs and Network Optimization*. R package version 1.0. 6

[3] Bacharoglou, A. (2010). Approximation of probability distributions by convex mixtures of Gaussian measures. *Proceedings of the American Mathematical Society* 138(7), 2619–2628. 27

[4] Baudry, J.-P., A. E. Raftery, G. Celeux, K. Lo, and R. Gottardo (2010). Combining mixture components for clustering. *Journal of computational and graphical statistics* 19(2), 332–353. 7, 15

[5] Bensmail, H., G. Celeux, A. E. Raftery, and C. P. Robert (1997). Inference in model-based cluster analysis. *statistics and Computing* 7, 1–10. 12, 39

[6] Berkhof, J., I. van Mechelen, and A. Gelman (2003). A Bayesian approach to the selection and testing of mixture models. *Statistica Sinica* 13(2), 423–442. 2, 5, 16

[7] Bouveyron, C., G. Celeux, T. B. Murphy, and A. E. Raftery (2019). *Model-based clustering and classification for data science: with applications in R*, Volume 50. Cambridge University Press. 2

[8] Carpenter, B., A. Gelman, M. D. Hoffman, D. Lee, B. Goodrich, M. Betancourt, M. Brubaker, J. Guo, P. Li, and A. Riddell (2017). Stan: A probabilistic programming language. *Journal of Statistical Software* 76(1), 1–32. 36

[9] Celeux, G., S. Frühwirth-Schnatter, and C. P. Robert (2019). Model selection for mixture models–perspectives and strategies. In *Handbook of mixture analysis*, pp. 117–154. Chapman and Hall/CRC. 1, 2, 13, 39, 40

[10] Celeux, G., M. Hurn, and C. P. Robert (2000). Computational and inferential difficulties with mixture posterior distributions. *Journal of the American Statistical Association* 95(451), 957–970. 8

[11] Chib, S. (1995). Marginal likelihood from the Gibbs output. *Journal of the American Statistical Association* 90(432), 1313–1321. 2

[12] Coraggio, L. and P. Coretto (2023). Selecting the number of clusters, clustering models, and algorithms. a unifying approach based on the quadratic discriminant score. *Journal of Multivariate Analysis* 196, 105181. 15

[13] Fraley, C. and A. E. Raftery (2007). Bayesian regularization for normal mixture estimation and model-based clustering. *Journal of classification* 24(2), 155–181. 39

[14] Frühwirth-Schnatter, S. (2001). Markov chain Monte Carlo estimation of classical and dynamic switching and mixture models. *Journal of the American Statistical Association* 96(453), 194–209. 4
[15] Frühwirth-Schnatter, S. (2004). Estimating marginal likelihoods for mixture and Markov switching models using bridge sampling techniques. *Econometrics Journal* 7, 143–167. 2

[16] Frühwirth-Schnatter, S. (2006). *Finite mixture and Markov switching models*. Springer. 1, 28, 36, 38, 45

[17] Frühwirth-Schnatter, S. (2019). Keeping the balance—bridge sampling for marginal likelihood estimation in finite mixture, mixture of experts and Markov mixture models. *Brazilian Journal of Probability and Statistics* 33(4), 706–733. 2

[18] Gelfand, A. E. and D. K. Dey (1994). Bayesian model choice: asymptotics and exact calculations. *Journal of the Royal Statistical Society: Series B (Methodological)* 56, 501–514. 26

[19] Gelman, A., J. B. Carlin, H. S. Stern, and D. B. Rubin (1995). *Bayesian data analysis*. Chapman and Hall/CRC. 45

[20] Geweke, J. (2007). Interpretation and inference in mixture models: Simple MCMC works. *Computational Statistics & Data Analysis* 51(7), 3529–3550. 8

[21] Ghojogh, B. and M. Crowley (2019). Linear and quadratic discriminant analysis: Tutorial. 8

[22] Ghosal, S. (2000). Asymptotic normality of posterior distributions for exponential families when the number of parameters tends to infinity. *Journal of Multivariate Analysis* 74, 49–68. 8

[23] Goldfarb, D. and A. Idnani (1983). A numerically stable dual method for solving strictly convex quadratic programs. *Mathematical programming* 27(1), 1–33. 35

[24] Goldfarb, D. and A. Idnani (2006). Dual and primal-dual methods for solving strictly convex quadratic programs. In *Numerical Analysis: Proceedings of the Third IIMAS Workshop Held at Cocoyoc, Mexico, January 1981*, pp. 226–239. Springer. 35

[25] Gronau, Q. F., H. Singmann, and E.-J. Wagenmakers (2020). bridgesampling: An R package for estimating normalizing constants. *Journal of Statistical Software* 92(10), 1–29. 2, 11

[26] Gruen, B. (2023). bayesmix: Bayesian Mixture Models with JAGS. R package version 0.7-6. 40

[27] Grün, B. and F. Leisch (2009). Dealing with label switching in mixture models under genuine multimodality. *Journal of Multivariate Analysis* 100(5), 851–861. 27

[28] Hairault, A., C. P. Robert, and J. Rousseau (2022). Evidence estimation in finite and infinite mixture models and applications. 2, 45

[29] Hastings, W. K. (1970). Monte Carlo sampling methods using Markov chains and their applications. *Biometrika* 57(1), 97–109. 2

[30] Heyde, C. C. and I. M. Johnstone (1979). On asymptotic posterior normality for stochastic processes. *Journal of the Royal Statistical Society: Series B (Methodological)* 41, 184–189. 8
[31] Jeffreys, H. (1961). *Theory of Probability* (3rd ed.). Oxford, U.K.: Oxford University Press. 1

[32] Klenke, A. (2013). *Probability theory: a comprehensive course*. Springer Science & Business Media. 25

[33] Knuth, D. E. and J. L. Szwarcfiter (1974). A structured program to generate all topological sorting arrangements. *Information Processing Letters* 2(6), 153-157. 35

[34] Lee, J. E. and C. P. Robert (2016). Importance sampling schemes for evidence approximation in mixture models. *Bayesian Analysis* 11(2), 573–597. 2, 16

[35] Lenk, P. J. and W. S. DeSarbo (2000). Bayesian inference for finite mixtures of generalized linear models with random effects. *Psychometrika* 65(1), 93–119. 8

[36] Liver-Disorders (2016). Liver Disorders. UCI Machine Learning Repository. DOI: https://doi.org/10.24432/C54G67. 15

[37] Llorente, F., L. Martino, D. Delgado, and J. Lopez-Santiago (2023). Marginal likelihood computation for model selection and hypothesis testing: An extensive review. *SIAM Review* 65, 3–58. 1

[38] Meng, X.-L. and W. H. Wong (1996). Simulating ratios of normalizing constants via a simple identity: A theoretical exploration. *Statistica Sinica* 6(4), 831–860. 2, 11

[39] Metodiev, M., N. J. Irons, and M. Perrot-Dockès (2025). *thamesmix: Truncated Harmonic Mean Estimator of the Marginal Likelihood for Mixtures*. R package version 0.1.1. 16

[40] Metodiev, M., N. J. Irons, M. Perrot-Dockès, P. Latouche, and A. E. Raftery (2025). Supplement to "easily computed marginal likelihoods for multivariate mixture models using the thames estimator". 3, 5, 6, 9, 10, 11, 15

[41] Metodiev, M., M. Perrot-Dockès, S. Ouadah, N. J. Irons, P. Latouche, and A. E. Raftery (2024). Easily computed marginal likelihoods from posterior simulation using the THAMES estimator. *Bayesian Analysis (advance publication)*, 1 – 28. 1, 2, 3, 5, 6, 27, 28, 29, 33, 34, 35

[42] Metropolis, N., A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller (1953). Equation of state calculations by fast computing machines. *Journal of Chemical Physics* 21, 1087–1092. 2

[43] Miller, J. W. (2021). Asymptotic normality, concentration, and coverage of generalized posteriors. *The Journal of Machine Learning Research* 22, 7598–7650. 8

[44] Neal, R. M. (1996). Sampling from multimodal distributions using tempered transitions. *Statistics and computing* 6, 353–366. 2

[45] Neal, R. M. (1999). Erroneous results in “marginal likelihood from the Gibbs output”. *Unpublished manuscript*. 2, 13

[46] Nobile, A. (2004). On the posterior distribution of the number of components in a finite mixture. *The Annals of Statistics* 32(5), 2044 – 2073. 3, 10, 32
[47] Nobile, A. (2007). Bayesian finite mixtures: a note on prior specification and posterior computation. *arXiv preprint arXiv:0711.0458*. 2, 3, 10, 32

[48] Nowicki, K. and T. A. B. Snijders (2001). Estimation and prediction for stochastic blockstructures. *Journal of the American statistical association* 96 (455), 1077–1087.
16

[49] Papastamoulis, P. and G. Iliopoulos (2010). An artificial allocations based solution to the label switching problem in Bayesian analysis of mixtures of distributions. *Journal of Computational and Graphical Statistics* 19(2), 313–331. 2, 5

[50] Perrakis, K., I. Ntzoufras, and E. G. Tsionas (2014). On the use of marginal posteriors in marginal likelihood estimation via importance sampling. *Computational Statistics & Data Analysis* 77, 54–69. 2

[51] Reichl, J. (2020). Estimating marginal likelihoods from the posterior draws through a geometric identity. *Monte Carlo Methods and Applications* 26(3), 205–221. 2, 3, 4, 5, 11, 12, 34, 37

[52] Richardson, S. and P. J. Green (1997). On Bayesian analysis of mixtures with an unknown number of components (with discussion). *Journal of the Royal Statistical Society Series B: Statistical Methodology* 59(4), 731–792. 8, 13, 39

[53] Robert, C. P., D. Wraith, P. M. Goggans, and C.-Y. Chan (2009). Computational methods for Bayesian model choice. In *AIP Conference Proceedings*, pp. 251–262. AIP. 3

[54] Rufo, M. J., J. Martín, and C. Pérez (2010). New approaches to compute Bayes factor in finite mixture models. *Computational statistics & data analysis* 54(12), 3324–3335. 2

[55] Scrucca, L., C. Fraley, T. B. Murphy, and A. E. Raftery (2023). *Model-Based Clustering, Classification, and Density Estimation Using mclust in R*. Chapman and Hall/CRC. 16

[56] Shen, X. (2002). Asymptotic normality of semiparametric and nonparametric posterior distributions. *Journal of the American Statistical Association* 97, 222–235. 8

[57] Sims, C. A., D. F. Waggoner, and T. Zha (2008). Methods for inference in large multiple-equation Markov-switching models. *Journal of Econometrics* 146(2), 255–274. 4

[58] Stan Development Team (2022). RStan: the R interface to Stan. R package version 2.21.7. 2

[59] Stephens, M. (1997). *Bayesian methods for mixtures of normal distributions*. Phd thesis, University of Oxford. Available at https://stephenslab.uchicago.edu/assets/papers/DPhilMS.pdf. 8

[60] Stephens, M. (2000). Dealing with label switching in mixture models. *Journal of the Royal Statistical Society: Series B (Statistical Methodology)* 62(4), 795–809. 2, 5
[61] Turlach, B. A. and A. Weingessel (2019). *quadprog: Functions to Solve Quadratic Programming Problems*. R package version 1.5-8. 35

[62] Wang, Y. J. and G. Y. Wong (1987). Stochastic blockmodels for directed graphs. *Journal of the American Statistical Association* 82(397), 8–19. 16

[63] Wilcox, R. R. (2011). *Introduction to robust estimation and hypothesis testing*. Academic press. 34
# Supplementary Material

## Table of contents

<table><tr><td><strong>Supplement A: proofs</strong></td><td style="text-align: right;">23</td></tr><tr><td>A.1 Proving that an optimal truncation set is given by a HPD region</td><td style="text-align: right;">23</td></tr><tr><td>A.2 Proving theoretical properties of symmetric estimators</td><td style="text-align: right;">26</td></tr><tr><td>A.3 Proving that the variance of the THAMES minimises an upper bound</td><td style="text-align: right;">27</td></tr><tr><td>A.4 Proving that we can easily compute the THAMES</td><td style="text-align: right;">29</td></tr><tr><td>A.5 Proving that we can avoid dealing with empty components</td><td style="text-align: right;">31</td></tr><tr><td><strong>Supplement B: details on the algorithm</strong></td><td style="text-align: right;">33</td></tr><tr><td>B.1 Choosing the tuning-parameters</td><td style="text-align: right;">33</td></tr><tr><td>B.2 Choosing the grid for α</td><td style="text-align: right;">34</td></tr><tr><td>B.3 Numerically verifying that ξ<sub>g1</sub> = ξ<sub>g2</sub> lies in E<sub>θ̂,Σ̂,c</sub></td><td style="text-align: right;">35</td></tr><tr><td>B.4 Quickly computing the THAMES</td><td style="text-align: right;">35</td></tr><tr><td><strong>Supplement C: details on the experiments</strong></td><td style="text-align: right;">36</td></tr><tr><td>C.1 Simulations on univariate Gaussian mixtures</td><td style="text-align: right;">36</td></tr><tr><td>C.2 Simulations on multivariate Gaussian mixtures</td><td style="text-align: right;">38</td></tr><tr><td>C.3 Univariate datasets (galaxy, enzyme, acidity)</td><td style="text-align: right;">39</td></tr><tr><td>C.4 Multivariate datasets (Swiss banknotes, BUPA liver)</td><td style="text-align: right;">40</td></tr><tr><td><strong>Supplement D: exact and approximate density computations</strong></td><td style="text-align: right;">43</td></tr><tr><td>D.1 Exact marginal likelihood for a small dataset</td><td style="text-align: right;">43</td></tr><tr><td>D.2 Highly accurate marginal likelihood approximation for an extremely well separated dataset</td><td style="text-align: right;">44</td></tr></table>
## Supplement A: proofs

Please note that the optimality results in Supplements A1 and A3 require additional assumptions on the distribution of the MCMC sample, as pointed out in the main document. However, the results on symmetry, unbiasedness, consistency and asymptotic normality in Supplement A2 make no such assumptions.

### A.1 Proving that an optimal truncation set is given by a HPD region

The following theorem proves that the variance of the estimator

$$ \hat{Z}^{-1} = \frac{1}{T/2} \sum_{\substack{t=T/2+1, \\ \theta^{(t)} \in A}}^{T} \frac{1/V(A)}{\pi(\theta^{(t)})L(\theta^{(t)})}, $$

is minimised by a $\alpha$-HPD region under conditions on which this region is well-defined, if the samples $\theta^{(1)}, \dots, \theta^{(T)}$ independently follow the posterior distribution. The condition on the distribution of $\theta^{(1)}, \dots, \theta^{(T)}$ is rarely satisfied exactly, but should be satisfied approximately, for example if a sufficiently large burn-in period is used and the sample is sufficiently heavily thinned.

**Theorem 5.1.** Assume that $\pi(\theta)L(\theta)$ is continuous on its support and that $\pi(\theta)L(\theta)$ is non-constant, in the sense that the volume of each level set $V(\{\theta|\pi(\theta)L(\theta) = q\}) = 0$ for all $q > 0$. Let

$$ H_{\alpha} = \{\theta | \pi(\theta)L(\theta) > q_{\alpha}\}, \qquad (12) $$

where $q_{\alpha}$ is the constant for which the posterior probability of $\{\theta \in H_{\alpha}\}$ is equal to $\alpha$. If the draws $\theta^{(1)}, \dots, \theta^{(T)} \sim p(\theta|D)$ are independent, there exists a constant $\alpha \in (0, 1]$ such that $A = H_{\alpha}$ minimises the variance of the THAMES over all measurable sets $A$, conditional on $D$.

*Proof.* Let $A$ be an arbitrary set of non-zero, finite volume that is contained in the posterior support. Because of independence and unbiasedness, the variance of $\hat{Z}^{-1}$ is given by

$$ \begin{aligned} \mathrm{Var}[\hat{Z}^{-1}|D] &= \mathrm{Var} \left[ \frac{1}{T/2} \sum_{\substack{t=T/2+1, \\ \theta^{(t)} \in A}}^{T} \frac{1/V(A)}{\pi(\theta^{(t)})L(\theta^{(t)})} \middle| D \right] \\ &= \frac{1}{T/2} \mathrm{Var} \left[ \frac{1_A(\theta^{(1)})/V(A)}{\pi(\theta^{(1)})L(\theta^{(1)})} \middle| D \right] \\ &= \frac{1}{T/2} \mathrm{E} \left[ \left( \frac{1_A(\theta^{(1)})/V(A)}{\pi(\theta^{(1)})L(\theta^{(1)})} \right)^2 \middle| D \right] - \frac{Z^{-2}}{T/2} \\ &= \frac{1}{T/2} \int_A \left( \frac{1/V(A)}{\pi(\theta)L(\theta)} \right)^2 p(\theta|D) \, d\theta - \frac{Z^{-2}}{T/2} \end{aligned} $$
$$ = \frac{Z^{-2}}{T/2} \int_A \frac{1/V(A)^2}{p(\theta|\mathcal{D})} d\theta - \frac{Z^{-2}}{T/2}, $$

where $\mathbb{1}_A(\theta^{(1)})$ denotes the indicator function of A evaluated at $\theta^{(1)}$. We are going to show that there exists a level $\alpha$ such that the THAMES has an equal- or lower variance under $H_\alpha$.

### Proving that there exists an $\alpha$ such that $H_\alpha$ is optimal

Let $H_\alpha = \{\theta | p(\theta|\mathcal{D}) \ge q_\alpha\}$ denote an $\alpha$-HPD region. Since $V(H_\alpha)$ is continuous, we can choose $\alpha$ such that

$$ V(H_\alpha) = V(A). $$

We have just shown that the variance of the THAMES is, up to translation, proportional to

$$ \mathrm{Var}[\hat{Z}^{-1}|\mathcal{D}] \propto \int_A \frac{1/V(A)^2}{p(\theta|\mathcal{D})} d\theta. \qquad (13) $$

Comparing this to the same integral of the THAMES under $H_\alpha$ gives

$$ \begin{aligned} \mathrm{Var}[\hat{Z}^{-1}|\mathcal{D}] &\propto \frac{1}{V(A)^2} \int_A \frac{1}{p(\theta|\mathcal{D})} d\theta \\ &= \frac{1}{V(H_\alpha)^2} \int_{A \cap H_\alpha} \frac{1}{p(\theta|\mathcal{D})} d\theta + \frac{1}{V(H_\alpha)^2} \int_{A \setminus H_\alpha} \frac{1}{p(\theta|\mathcal{D})} d\theta \\ &\ge \frac{1}{V(H_\alpha)^2} \left( \int_{A \cap H_\alpha} \frac{1}{p(\theta|\mathcal{D})} d\theta + \frac{1}{q_\alpha} V(A \setminus H_\alpha) \right). \end{aligned} $$

This is due to the fact that $A \setminus H_\alpha = \{\theta \in A | p(\theta|\mathcal{D}) < q_\alpha\}$. Using this and the fact that $V(A \setminus H_\alpha) = V(H_\alpha \setminus A)$ one can similarly conclude

$$ \begin{aligned} & \frac{1}{V(H_\alpha)^2} \left( \int_{A \cap H_\alpha} \frac{1}{p(\theta|\mathcal{D})} d\theta + \frac{1}{q_\alpha} V(A \setminus H_\alpha) \right) = \\ & \frac{1}{V(H_\alpha)^2} \left( \int_{A \cap H_\alpha} \frac{1}{p(\theta|\mathcal{D})} d\theta + \frac{1}{q_\alpha} V(H_\alpha \setminus A) \right) \ge \\ & \frac{1}{V(H_\alpha)^2} \int_{A \cap H_\alpha} \frac{1}{p(\theta|\mathcal{D})} d\theta + \frac{1}{V(H_\alpha)^2} \int_{H_\alpha \setminus A} \frac{1}{p(\theta|\mathcal{D})} d\theta = \\ & \frac{1}{V(H_\alpha)^2} \int_{H_\alpha} \frac{1}{p(\theta|\mathcal{D})} d\theta. \end{aligned} $$

It follows that for any set A with $V(A) \in (0, \infty)$ that is contained in the posterior support an $\alpha$ exists such that the variance of the THAMES is lower or equal when choosing $H_\alpha$ instead of A. This finishes the proof.

Please note that $H_\alpha$ is well-defined and $V(H_\alpha)$ is strictly increasing and continuous. This follows from the conditions given in the theorem. We could not find a direct proof in the literature (though the result is probably well-known), so we give one here.
### Proving that $H_\alpha$ is well-defined and $V(H_\alpha)$ is strictly increasing and continuous

Let $q \in (0, \infty)$ be arbitrary. Let $P(\cdot|\mathcal{D})$ denote the probability measure of the posterior distribution. It is well-known that any measure is continuous, in the sense that

$$ \lim_{i \to \infty} P(A_i|\mathcal{D}) = P\left(\bigcap_{i=1}^{\infty} A_i \mid \mathcal{D}\right) $$

for any decreasing series of sets $A_1 \supseteq A_2 \supseteq \dots$ [32]. Thus

$$ \begin{aligned} \lim_{\tilde{q} \uparrow q} P(\{\theta | \pi(\theta) L(\theta) > \tilde{q}\} | \mathcal{D}) &= P\left(\bigcap_{\tilde{q} \uparrow q} \{\theta | \pi(\theta) L(\theta) > \tilde{q}\} \mid \mathcal{D}\right) \\ &= P(\{\theta | \exists \tilde{q} > q : \pi(\theta) L(\theta) > \tilde{q}\} | \mathcal{D}) = P(\{\theta | \pi(\theta) L(\theta) > q\} | \mathcal{D}), \end{aligned} $$

for all $q \in \mathbb{R}$ due to the fact that $\{\theta | \pi(\theta) L(\theta) > \tilde{q}\}$ is contained in $\{\theta | \pi(\theta) L(\theta) > q\}$ for any $\tilde{q} > q$. Since $V(\{\theta | \pi(\theta) L(\theta) = \tilde{q}\}) = 0$, it is also the case that

$$ \begin{aligned} \lim_{\tilde{q} \downarrow q} P(\{\theta | \pi(\theta) L(\theta) > \tilde{q}\} | \mathcal{D}) &= 1 - \lim_{\tilde{q} \downarrow q} P(\{\theta | \pi(\theta) L(\theta) \le \tilde{q}\} | \mathcal{D}) \\ &= 1 - \lim_{\tilde{q} \downarrow q} P(\{\theta | \pi(\theta) L(\theta) < \tilde{q}\} | \mathcal{D}), \end{aligned} $$

and we can use the same technique as before to show that

$$ \lim_{\tilde{q} \downarrow q} P(\{\theta | \pi(\theta) L(\theta) > \tilde{q}\} | \mathcal{D}) = P(\{\theta | \pi(\theta) L(\theta) > q\} | \mathcal{D}). $$

Since we have just shown that the function $q \mapsto P(\{\theta | \pi(\theta) L(\theta) > q\} | \mathcal{D})$ is continuous, and since it tends to 0 as $q \to \infty$, to 1 as $q \downarrow 0$, it follows that the equation

$$ P(\{\theta | \pi(\theta) L(\theta) > q\} | \mathcal{D}) = \alpha \qquad (14) $$

has at least one solution. In addition, the function is also strictly decreasing in the interval $(\inf_{\pi(\theta)L(\theta)>0} \pi(\theta)L(\theta), \sup_{\pi(\theta)L(\theta)>0} \pi(\theta)L(\theta))$, since for any $\inf_{\pi(\theta)L(\theta)>0} \pi(\theta)L(\theta) < q_1 < q_2 < \sup_{\pi(\theta)L(\theta)>0} \pi(\theta)L(\theta)$

$$ \begin{aligned} P(\{\theta | \pi(\theta) L(\theta) > q_1\} | \mathcal{D}) - P(\{\theta | \pi(\theta) L(\theta) > q_2\} | \mathcal{D}) \\ = P(\{\theta | \pi(\theta) L(\theta) \in (q_1, q_2]\} | \mathcal{D}) > 0, \end{aligned} $$

because the set $\{\theta | \pi(\theta) L(\theta) \in (q_1, q_2]\}$ has non-zero volume due to the fact that $\pi(\theta)L(\theta)$ is continuous and non-constant. It follows that Equation (14) has only one solution, defined as $q_\alpha$. We can use the exact same technique to show that the volume $V(H_\alpha)$ is strictly increasing and continuous, since the only property of $P(\cdot|\mathcal{D})$ that we used is that it is a measure, and the volume is also a measure. $\square$
## A.2 Proving theoretical properties of symmetric estimators

Let $\theta^{(1)}, \dots, \theta^{(T)}$ denote the sample from the posterior and $\theta^{(1)\star}, \dots, \theta^{(T)\star}$ denote the same sample, but relabelled. Also, let $Q$ be a symmetric estimator, in the sense that

$$Q(P_{j_1}(\theta^{(1)}), \dots, P_{j_T}(\theta^{(T)})) = Q(\theta^{(1)}, \dots, \theta^{(T)}),$$

for any sequence of permutations $P_{j_1}(\theta^{(1)}), \dots, P_{j_T}(\theta^{(T)})$ and for all possible values of the posterior sample. We are going to show that $Q(\theta^{(1)}, \dots, \theta^{(T)})$ and $Q(\theta^{(1)\star}, \dots, \theta^{(T)\star})$ have the same distribution. This implies that any symmetric estimator that is consistent for a sample of the posterior is also going to be consistent for the relabelling of that sample and, on top of that, all theoretical properties of its distribution, such as asymptotic normality, are preserved when using the relabelled posterior sample.

**Theorem 5.2.** $Q(\theta^{(1)}, \dots, \theta^{(T)})$ and $Q(\theta^{(1)\star}, \dots, \theta^{(T)\star})$ have the same distribution.

*Proof.* Since $\theta^{(1)\star}, \dots, \theta^{(T)\star}$ is a relabelling of $\theta^{(1)}, \dots, \theta^{(T)}$, this means that for any realisation of $\theta^{(1)}, \dots, \theta^{(T)}$ there exists a sequence of permutations $P_{j_1}, \dots, P_{j_T}$ such that

$$(\theta^{(1)\star}, \dots, \theta^{(T)\star}) = (P_{j_1}(\theta^{(1)}), \dots, P_{j_T}(\theta^{(T)})).$$

This immediately implies

$$Q(\theta^{(1)\star}, \dots, \theta^{(T)\star}) = Q(\theta^{(1)}, \dots, \theta^{(T)})$$

due to symmetry. Since this is true for every realisation of the posterior, we have equality in distribution. $\square$

## Proving that the THAMES is consistent and asymptotically normal

The THAMES is defined as

$$\hat{Z}_S^{-1}(G) = \frac{1}{G!} \sum_{o=1}^{G!} \frac{1}{T/2} \sum_{\substack{t=T/2+1, \\ P_o(\theta^{(t)\star}) \in B_{\hat{\theta}, \hat{\Sigma}, c, \alpha}}} \frac{1}{\pi(\theta^{(t)\star})L(\theta^{(t)\star})}.$$

This is a reciprocal importance sampling (RIS) estimator [18], where an RIS estimator on the posterior sample $\theta^{(T/2+1)}, \dots, \theta^{(T)}$ is defined as $\frac{1}{T/2} \sum_{t=T/2+1}^{T} \frac{h(\theta^{(t)})}{\pi(\theta^{(t)})L(\theta^{(t)})}$, with $h$ denoting any density function specified independently from this sample whose support is within the posterior support. It is well-known that this type of estimator is unbiased on the scale of the reciprocal marginal likelihood, consistent, asymptotically normal and has finite variance as long as the condition

$$\int \frac{h(\theta)^2}{L(\theta)\pi(\theta)} d\theta < \infty \qquad (15)$$
holds. The symmetric THAMES is in fact a RIS estimator with

$$h(\theta^{(t)\star}) = \frac{1}{G!} \sum_{o=1}^{G!} \frac{\mathbb{1}_{B_{\hat{\theta},\hat{\Sigma},c,\alpha}}(P_o(\theta^{(t)\star}))}{V(B_{\hat{\theta},\hat{\Sigma},c,\alpha})}$$

being a mixture of $G!$ uniform distributions on $B_{\hat{\theta},\hat{\Sigma},c,\alpha}$ and $\mathbb{1}_{B_{\hat{\theta},\hat{\Sigma},c,\alpha}}$ denoting the indicator function on $B_{\hat{\theta},\hat{\Sigma},c,\alpha}$. This follows from changing the order of summation:

$$\frac{1}{T/2} \sum_{t=T/2+1}^{T} \frac{h(\theta^{(t)\star})}{\pi(\theta^{(t)\star})L(\theta^{(t)\star})} = \frac{1}{T/2} \sum_{t=T/2+1}^{T} \frac{\frac{1}{G!} \sum_{o=1}^{G!} \frac{\mathbb{1}_{B_{\hat{\theta},\hat{\Sigma},c,\alpha}}(P_o(\theta^{(t)\star}))}{V(B_{\hat{\theta},\hat{\Sigma},c,\alpha})}}{\pi(\theta^{(t)\star})L(\theta^{(t)\star})} \\ = \frac{1}{G!} \sum_{o=1}^{G!} \frac{1}{T/2} \sum_{t=T/2+1, \quad P_o(\theta^{(t)\star}) \in B_{\hat{\theta},\hat{\Sigma},c,\alpha}}^{T} \frac{1/V(B_{\hat{\theta},\hat{\Sigma},c,\alpha})}{\pi(\theta^{(t)\star})L(\theta^{(t)\star})} = \hat{Z}_{\mathrm{S}}^{-1}(G)$$

Thus, the posterior density is bounded from below whenever $h$ is not 0 and the support of $h$ is finite, so Equation (15) is satisfied due to boundedness. In addition, while the symmetric THAMES is defined on the relabelled, not the full posterior sample, we already showed that this makes no difference for a symmetric estimator. It follows that the THAMES is indeed unbiased on the scale of the reciprocal marginal likelihood, consistent and asymptotically normally distributed.

### A.3 Proving that the variance of the THAMES minimises an upper bound

Consider a convex combination of truncated harmonic mean estimators,

$$\hat{Z}_{\mathrm{Mix}}^{-1} = \sum_{k=1}^{K} \omega_k \frac{1}{T/2} \sum_{\substack{t=T/2+1, \\ \theta^{(t)} \in A_k}}^{T} \frac{1/V(A_k)}{\pi(\theta^{(t)})L(\theta^{(t)})}, \quad (16)$$

where $\omega_1, \dots, \omega_K$ are weights that sum to 1 and $A_1, \dots, A_K$ are truncation sets of the different estimators. Optimality results were obtained in Metodiev et al. [41] in the case that $K = 1$, under the assumption of normality of the posterior. We are going to present results in the case that the posterior is a mixture of normal distributions. This case is important because any density can be approximated by such a mixture, under few assumptions [3].

It is easy to show that the symmetric THAMES corresponds to Equation (16) with the optimal parameters that we suggest so long as the posterior density of the relabelled sample is unimodal. This isn't necessarily the case, since the phenomenon of genuine multimodality, that is, multimodality in the posterior density of the relabelled sample, can be observed [27]. In such cases, where the number of genuine modes is larger than 1, it may be useful to average different versions of the symmetric THAMES, each centered
in a different genuine mode, by applying mclust to the relabelled posterior sample. However, the assumption of only G! symmetric modes may very well be justified by the large sample behaviour of the posterior: given an identifiable (up to permutations) mixture model, the posterior itself can often be observed to converge to a mixture of G! multivariate Gaussian distribution, as pointed out in Frühwirth-Schnatter [16]. This is exactly the condition under which we constructed the THAMES for multimodality.

**Theorem 5.3.** Suppose that the posterior density $p(\theta|\mathcal{D})$ is a Gaussian mixture, $p(\theta|\mathcal{D}) = \sum_{k=1}^K a_k \text{MVN}_R(\theta; m_k, S_k)$, where $\text{MVN}_R(\theta; m_k, S_k)$ denotes the Gaussian law of dimension $R$ with mean $m_k$ and covariance matrix $S_k$. If the draws $\theta^{(1)}, \dots, \theta^{(T)}$ from the posterior distribution are independent, an upper bound on the variance of $\hat{Z}_{\text{Mix}}^{-1}$, conditional on $\mathcal{D}$, is minimised by $A_k = E_{m_k, S_k, c_R}$, $\omega_k = a_k$, where $c_R$ is the minimising constant determined in Metodiev et al. [41], which is asymptotically equal to $\sqrt{R+1}$, i.e., $\lim_{R \to \infty} \frac{c_R}{\sqrt{R+1}} = 1$.

*Proof.* Because of independence, unbiasedness, the Cauchy-Schwartz and the QM-AM inequality, the variance of $\hat{Z}_{\text{Mix}}^{-1}$ is given by

$$
\begin{align*}
\operatorname{Var}[\hat{Z}_{\text{Mix}}^{-1}|\mathcal{D}] &= \operatorname{Var}\left[\sum_{k=1}^{K} \omega_k \frac{1}{T/2} \sum_{\substack{t=T/2+1, \\ \theta^{(t)} \in A_k}}^{T} \frac{1/V(A_k)}{\pi(\theta^{(t)})L(\theta^{(t)})} \middle| \mathcal{D}\right] \stackrel{\text{i.i.d}}{=} \\
&\quad \frac{1}{T/2} \operatorname{Var}\left[\sum_{k=1}^{K} \omega_k \frac{1_{A_k}(\theta^{(1)})/V(A_k)}{\pi(\theta^{(1)})L(\theta^{(1)})} \middle| \mathcal{D}\right] \stackrel{\text{unbiased}}{=} \\
&\quad \frac{1}{T/2} \operatorname{E}\left[\left(\sum_{k=1}^{K} \omega_k \frac{1_{A_k}(\theta^{(1)})/V(A_k)}{\pi(\theta^{(1)})L(\theta^{(1)})}\right)^2 \middle| \mathcal{D}\right] - \frac{Z^{-2}}{T/2} \stackrel{\text{linearity}}{=} \\
&\quad \frac{1}{T/2} \sum_{j,k=1}^{K} \omega_k \omega_j \operatorname{E}\left[\frac{1_{A_k \cap A_j}(\theta^{(1)})/(V(A_k)V(A_j))}{(\pi(\theta^{(1)})L(\theta^{(1)}))^2} \middle| \mathcal{D}\right] - \frac{Z^{-2}}{T/2} \stackrel{\text{Cauchy Schwartz}}{\le} \\
&\quad \frac{1}{T/2} \sum_{j,k=1}^{K} \sqrt{\operatorname{E}\left[\omega_k^2 \frac{1_{A_k}(\theta^{(1)})/V(A_k)^2}{(\pi(\theta^{(1)})L(\theta^{(1)}))^2} \middle| \mathcal{D}\right]} \operatorname{E}\left[\omega_j^2 \frac{1_{A_j}(\theta^{(1)})/V(A_j)^2}{(\pi(\theta^{(1)})L(\theta^{(1)}))^2} \middle| \mathcal{D}\right] \\
&\quad - \frac{Z^{-2}}{T/2} \stackrel{\text{QM-AM}}{\le} \\
&\quad \frac{K^2}{T/2} \sqrt{\frac{1}{K^2} \sum_{j,k=1}^{K} \omega_k^2 \omega_j^2 \operatorname{E}\left[\frac{1_{A_k}(\theta^{(1)})/V(A_k)^2}{(\pi(\theta^{(1)})L(\theta^{(1)}))^2} \middle| \mathcal{D}\right]} \operatorname{E}\left[\frac{1_{A_j}(\theta^{(1)})/V(A_j)^2}{(\pi(\theta^{(1)})L(\theta^{(1)}))^2} \middle| \mathcal{D}\right] \\
&\quad - \frac{Z^{-2}}{T/2} = \\
&\quad \frac{K^2}{T/2} \left( \frac{1}{K} \sum_{k=1}^{K} \omega_k^2 \operatorname{E}\left[\frac{1_{A_k}(\theta^{(1)})/V(A_k)^2}{(\pi(\theta^{(1)})L(\theta^{(1)}))^2} \middle| \mathcal{D}\right] \right) - \frac{Z^{-2}}{T/2} =
\end{align*}
$$
$$ \frac{K}{T/2} \sum_{k=1}^{K} \omega_k^2 E \left[ \frac{1_{A_k}(\theta^{(1)})/V(A_k)^2}{(\pi(\theta^{(1)})L(\theta^{(1)}))^2} \middle| \mathcal{D} \right] - \frac{Z^{-2}}{T/2}. $$

We can use the fact that $p(\theta|\mathcal{D}) \ge a_k \text{MVN}_d(\theta; m_k, S_k)$ for any $k$ to get

$$ \begin{aligned} & \frac{K}{T/2} \sum_{k=1}^{K} \omega_k^2 E \left[ \frac{1_{A_k}(\theta^{(1)})/V(A_k)^2}{(\pi(\theta^{(1)})L(\theta^{(1)}))^2} \middle| \mathcal{D} \right] - \frac{Z^{-2}}{T/2} = \\ & \frac{K}{T/2} \sum_{k=1}^{K} \omega_k^2 \int_{A_k} \frac{p(\theta|\mathcal{D})/V(A_k)^2}{(\pi(\theta)L(\theta))^{2/3}} d\theta - \frac{Z^{-2}}{T/2} = \\ & \frac{K}{T/2} \sum_{k=1}^{K} \omega_k^2 \int_{A_k} \frac{p(\theta|\mathcal{D})/V(A_k)^2}{(Zp(\theta|\mathcal{D}))^2} d\theta - \frac{Z^{-2}}{T/2} = \\ & \frac{KZ^{-2}}{T/2} \sum_{k=1}^{K} \omega_k^2 \int_{A_k} \frac{1/V(A_k)^2}{p(\theta|\mathcal{D})} d\theta - \frac{Z^{-2}}{T/2} \le \\ & \frac{KZ^{-2}}{T/2} \sum_{k=1}^{K} \frac{\omega_k^2}{a_k} \int_{A_k} \frac{1/V(A_k)^2}{\text{MVN}_R(\theta; m_k, S_k)} d\theta - \frac{Z^{-2}}{T/2}. \end{aligned} $$

$\text{MVN}_R(\theta; m_k, S_k)$ is itself a density, and the integral above is exactly the integral found in the proof of Theorem 1, Equation (13). Thus, this integral is proportional to the variance of a THAMES estimator for which the posterior density is exactly normal. This case has been dealt with in Metodiev et al. [41, Theorem 1], and the integral is indeed minimised by $E_{m_k,S_k,c_R}$. The transformation formula implies that the integral is independent of $k$ in this case, since

$$ \begin{aligned} & \frac{K}{T/2} \sum_{k=1}^{K} \frac{\omega_k^2}{a_k} \int_{E_{m_k, S_k, c_R}} \frac{Z^{-2}/V(E_{m_k, S_k, c_R})^2}{\text{MVN}_R(\theta; m_k, S_k)} d\theta - \frac{Z^{-2}}{T/2} = \\ & \frac{K}{T/2} \sum_{k=1}^{K} \frac{\omega_k^2}{a_k} \int_{E_{0_R, I_R, c_R}} \frac{Z^{-2}/V(E_{0_R, I_R, c_R})^2}{\text{MVN}_R(\theta; 0_R, I_R)} d\theta - \frac{Z^{-2}}{T/2} = \\ & \frac{K}{T/2} \int_{E_{0_R, I_R, c_R}} \frac{Z^{-2}/V(E_{0_R, I_R, c_R})^2}{\text{MVN}_R(\theta; 0_R, I_R)} d\theta \left( \sum_{k=1}^{K} \frac{\omega_k^2}{a_k} \right) - \frac{Z^{-2}}{T/2}, \end{aligned} $$

and the quadratic mean is minimised by $\omega_k = a_k$. $\square$

<figure><img src="image_9.png" /></figure>

### A.4 Proving that we can easily compute the THAMES

**Theorem 5.4.** Let $\psi_1, \dots, \psi_G!$ denote the $G!$ different ordering algorithms that map $\theta$ to different permutations of the ordering constraint

$$ W(\xi_1) \le \dots \le W(\xi_G), $$
<figure><img src="image_10.png" /><figcaption>Figure 7: The values of $\mu_g, \sigma_g$ and $W$ for a sample $(\mu_g^{(1)}, \sigma_g^{(1)}), \dots, (\mu_g^{(T)}, \sigma_g^{(T)})$ from the posterior of the galaxy dataset of size $T = 10000$ and different values of $G$; the acceptance regions of QDA are shown via a grid in the background</figcaption></figure>

respectively, by permuting its labels. Let $\Delta$ denote the adjacency matrix defined by $\Delta_{g_1,g_2} = 1$ if $W(\xi_{g_1}) < W(\xi_{g_2})$ for all $\xi_{g_1}, \xi_{g_2}$ that satisfy the ellipsoidal constraint on $E_{\hat{\theta}, \hat{\Sigma}, c}$, and let $\Omega$ denote the index set over the collection of all topological orderings on the graph defined by $\Delta$. The THAMES is equal to

$$ \hat{Z}_S^{-1}(G) = \frac{1}{G!} \sum_{o \in \Omega} \frac{1}{T/2} \sum_{\substack{t=T/2+1 \\ \psi_o(\theta^{(t)*}) \in B_{\hat{\theta}, \hat{\Sigma}, c, \alpha}}}^{T} \frac{1/V(B_{\hat{\theta}, \hat{\Sigma}, c, \alpha})}{\pi(\theta^{(t)*})L(\theta^{(t)*})} \quad (17) $$

The role of the function $W$ can be explained via Figure 7: it shows the different values of $(\mu_1^{(t)}, \sigma_1^{(t)}), \dots, (\mu_G^{(t)}, \sigma_G^{(t)})$ for an MCMC simulation from the posterior of size $T = 10000$ in conjunction with the values of $W$. Ordering by the means would not be enough in the case that $G = 2$, since they overlap, although the variances do not. Ordering by the variances would not be enough when $G = 3$ for the converse reason. However, the values of $W$ differ for most points whenever they are in a different component. The exception is component 4 in the case that $G = 4$: it is superfluous, and was thus not included into the calculation of $W$.
<figure><img src="image_11.png" /><figcaption>Figure 8: A toy example to illustrate the computation of the THAMES; the points that lie in the ellipsoid can be counted by applying $P_1, P_2$ or $\psi_1, \psi_2$.</figcaption></figure>

*Proof.* The key identity for this proof is that for every $t$ and for every function $\iota$ of $\theta$

$$ \sum_{o=1}^{G!} \iota(P_o(\theta^{(t)\star})) = \sum_{o=1}^{G!} \iota(\psi_o(\theta^{(t)\star})). $$

It is true because the functions $\psi_o, P_o$ each assign a different permutation to $\theta$. This can be illustrated via Figure 8: let $T = 10, \theta^{(t)\star} = (\mu_1^{(t)\star}, \mu_2^{(t)\star})$, and $\iota(\theta^{(t)\star}) = 1_{E_{\hat{\theta}, \hat{\Sigma}, c}}(\theta^{(t)\star})$. In this case, we are simply counting the number of times that $P_o(\theta^{(t)\star})$ is in $E_{\hat{\theta}, \hat{\Sigma}, c}$. This can be done by applying $P_1, P_2$ one by one and counting each time (plots A and B), or by ordering every parameter such that $\mu_1^{(t)} > \mu_2^{(t)}$ via $\psi_1$, and then switching the order to $\mu_2^{(t)} > \mu_1^{(t)}$ via $\psi_2$ (plots C and D). In fact, $\psi_2$ can be avoided entirely in Figure 8 since $E_{\hat{\theta}, \hat{\Sigma}, c}$ does not cross the identity line.

Since $\hat{Z}_S^{-1}$ is also a sum over all permutations of $\theta^{(t)}$, we can replace $P_o$ by $\psi_o$ to get

$$ \hat{Z}_S^{-1} = \frac{1}{G!} \sum_{o=1}^{G!} \frac{1}{T/2} \sum_{\substack{t=T/2+1, \\ \psi_o(\theta^{(t)\star}) \in B_{\hat{\theta}, \hat{\Sigma}, c, \alpha}}} \frac{1/V(B_{\hat{\theta}, \hat{\Sigma}, c, \alpha})}{\pi(\theta^{(t)\star})L(\theta^{(t)\star})}. $$

Suppose that $\psi_o$ is such that $W(\xi_{g_1}^{\psi_o}) > W(\xi_{g_2}^{\psi_o})$, but $\Delta_{g_1, g_2} = 1$. This immediately implies that $\psi_o(\theta^{(t)\star})$ will not be included in the sum for any $\theta^{(t)\star}$, since $W(\xi_{g_1}) \le W(\xi_{g_2})$ for all $\xi_{g_1}, \xi_{g_2}$ satisfying the ellipsoidal constraint on $E_{\hat{\theta}, \hat{\Sigma}, c}$, and thus in particular in $B_{\hat{\theta}, \hat{\Sigma}, c, \alpha}$. Thus, the sum only includes orderings which respect the partial ordering implied by $\Delta$. This is how topological orderings are defined. Since $\Omega$ only excludes those parts of the series that would be excluded anyway, this finishes the proof. $\square$

### A.5 Proving that we can avoid dealing with empty components

**Theorem 5.5.** Suppose that a Dirichlet prior with parameters $e_1, \dots, e_G > 0$ is used for the proportions and that $\hat{Z}(G-1)$ is a consistent estimator of $Z(G-1)$. If $\xi_1, \dots, \xi_G$ are
a priori independent and their prior distribution does not depend on the total number of components, then

$$ \hat{Z}(G) = \frac{\Gamma\left(\sum_{g=1}^{G} e_g\right) \Gamma\left(n + \sum_{g=1}^{G-1} e_g\right)}{\Gamma\left(n + \sum_{g=1}^{G} e_g\right) \Gamma\left(\sum_{g=1}^{G-1} e_g\right)} \cdot \hat{Z}(G-1) \cdot \frac{1}{\hat{p}_0(G)}, \quad (18) $$

with

$$ \hat{p}_0(G) = \frac{1}{G} \sum_{g=1}^{G} \frac{1}{T} \sum_{t=1}^{T} \prod_{i=1}^{n} (1 - \hat{z}_{i,g}^{(t)}), \quad \hat{z}_{i,g}^{(t)} = \frac{\tau_{g}^{(t)} f(Y_i; \xi_{g}^{(t)})}{\sum_{\tilde{g}=1}^{G} \tau_{\tilde{g}}^{(t)} f(Y_i; \xi_{\tilde{g}}^{(t)})} $$

is a consistent estimator of $Z(G)$ as long as $\hat{p}_0(G)$ is a consistent estimator of its posterior mean.

Note that $\hat{p}_0(G)$ is symmetric. Thus, it can be evaluated on the relabelled posterior sample with no consequences.

*Proof.* The data has the following latent variable representation:

$$ C_i | \tau \stackrel{\text{i.i.d}}{\sim} \text{Multinom}_G(1, \tau), \\ Y_i | \xi, C_{i,g} = 1 \sim f(\cdot; \xi_g), $$

with the hidden matrix of allocation vectors $C = (C_1, \dots, C_n) \in \mathcal{M}_{n \times G}(\{0, 1\})$.

Nobile [46, 47] showed that

$$ \frac{Z(G)}{Z(G-1)} = \frac{\Gamma\left(\sum_{g=1}^{G} e_g\right) \Gamma\left(n + \sum_{g=1}^{G-1} e_g\right)}{\Gamma\left(n + \sum_{g=1}^{G} e_g\right) \Gamma\left(\sum_{g=1}^{G-1} e_g\right)} \cdot \frac{1}{1 - p(\exists i : C_{i,G} = 1 | Y)}. \quad (19) $$

As pointed out by Nobile [47],

$$ 1 - p(\exists i : C_{i,G} = 1) = p\left(\sum_{i=1}^{n} C_{i,G} = 0 \mid Y\right), $$

the posterior probability that component $G$ is empty. By the law of total probability

$$ p\left(\sum_{i=1}^{n} C_{i,G} = 0 \mid Y\right) = \int p\left(\sum_{i=1}^{n} C_{i,G} = 0 \mid Y, \theta\right) p(\theta | Y) d\theta $$

The distribution of $C_i$ given $Y, \theta$ is i.i.d multinomial with probabilities

$$ z_{i,g} = \frac{\tau_g f(Y_i, \xi_g)}{\sum_{\tilde{g}=1}^{G} \tau_{\tilde{g}} f(Y_i; \xi_{\tilde{g}})}, \quad g = 1, \dots, G. $$
This simplifies the integral to

$$ \int p\left(\sum_{i=1}^{n} C_{i,G} = 0 \mid Y, \theta\right) p(\theta|Y) d\theta = \int \prod_{i=1}^{n} (1 - z_{i,G}) p(\theta|Y) d\theta, $$

which is in turn equal to the posterior expectation of $\prod_{i=1}^{n} (1 - z_{i,G})$. By the law of large numbers,

$$ \frac{1}{T} \sum_{t=1}^{T} \prod_{i=1}^{n} (1 - \hat{z}_{i,G}^{(t)}) \qquad (20) $$

is a valid estimator of this expectation. Due to the symmetry of the posterior, choosing any other component $g = 1, \dots, G$ instead of the last component $G$ would give the exact same result. Averaging over all of these versions produces a symmetric estimator. This estimator is exactly equal to $\hat{p}_0$. Combining the unbiasedness of $\hat{p}_0$ with the consistency of $\hat{Z}(G - 1)$ gives the assertion. □

# Supplement B: details on the algorithm

## B.1 Choosing the tuning-parameters

We choose $c = \sqrt{R+1}$, following the results of Metodiev et al. [41]. The hyperparameter $\alpha$ could be chosen via grid-search over $(0, 1)$ such that it minimises the empirical variance of the THAMES. However, while this approach may work, it is computationally expensive, since we would have to approximate the volume of $B_{\hat{\theta}, \hat{\Sigma}, c, \alpha}$ via a new Monte Carlo sample for each value that $\alpha$ can take on the grid. Instead, we suggest choosing $\alpha$ such that the posterior distribution is well-behaved on $B_{\hat{\theta}, \hat{\Sigma}, c, \alpha}$. We take this to mean that the distribution of the sample from the posterior within $B_{\hat{\theta}, \hat{\Sigma}, c, \alpha}$ is as close as possible to a normal distribution, truncated such that at most 50 percent of the posterior sample is within $B_{\hat{\theta}, \hat{\Sigma}, c, \alpha}$. This is because the tails of the posterior can cause difficulties in reciprocal importance sampling estimators and the tails of the normal distribution are particularly well behaved in the case of the THAMES if they are cut off at 50 percent, as was shown in Metodiev et al. [41].

If the posterior is normal, the negative unnormalised log-posterior values

$$ \eta^{(1)}, \dots, \eta^{(T)} = -\log(\pi(\theta^{(1)})L(\theta^{(1)})), \dots, -\log(\pi(\theta^{(T)})L(\theta^{(T)})) $$

follow, if they are truncated to an HPD region, a truncated, shifted and scaled $\chi^2$-distribution with $d$ degrees of freedom. This is because the quadratic form

$$ (\theta^{(t)} - m)^{\mathrm{T}} S^{-1} (\theta^{(t)} - m) $$

follows a $\chi^2$-distribution if $m$ and $S$ are the true posterior mean and covariance matrix, respectively.
We set $\alpha$ such that it minimises the Kolmogorov distance against this $\chi^2$-distribution. Let $F$ denote the empirical cumulative distribution function (ECDF) of $\eta^{(1)}, \dots, \eta^{(T)}$ conditioned on the $\alpha$-HPD region $H_\alpha$,

$$F(\eta) = \frac{1}{\sum_{t=1}^{T} 1_{H_{\alpha}}(\theta^{(t)})} \sum_{t=1}^{T} 1_{H_{\alpha}}(\theta^{(t)}) \cdot 1_{(-\infty, \eta]}(\eta^{(t)}),$$

with $1$ denoting the indicator function, and let $\chi_{M,s}^2$ denote the cumulative distribution function (CDF) of the shifted and scaled $\chi^2$ distribution truncated on $[0, \max\{\eta^{(1)}, \dots, \eta^{(T)}\}]$ with $d$ degrees of freedom, mean $M$ and variance $v$. The Kolmogorov distance [see for instance 63] on $H_\alpha$ is given by

$$\sup_{\theta \in H_\alpha} \left| \chi_{M_\alpha, \hat{v}_\alpha}^2 (-\log(\pi(\theta)L(\theta))) - F(-\log(\pi(\theta)L(\theta))) \right|, \quad (21)$$

where the mean and variance parameter are computed using the moment-estimators $\hat{M}_\alpha, \hat{v}_\alpha$ on the values of $\eta^{(1)}, \dots, \eta^{(T)}$ for which the corresponding parameters are within $H_\alpha$. We minimise this distance by evaluating Equation (21) on some subset of the grid $\{\theta^{(1)}, \dots, \theta^{(T)}\}$ for all sensible values of $\alpha$, with $\alpha$ being capped at 50 percent for the reasons that we mentioned before (for more details on this technique, see the following subsection).

**Summary** Choosing $\alpha$ such that it minimises Equation (21) has many advantages. An optimal value is easy to compute, since it can be efficiently approximated by evaluating at most $T^2$ different values due to the discreteness of $F$. It also conforms with previous results from Metodiev et al. [41]: if the posterior is exactly normal, the Kolmogorov distance is going to approach 0 if $\alpha = 1$, which implies $B_{\hat{\theta}, \hat{\Sigma}, c, \alpha} = E_{\hat{\theta}, \hat{\Sigma}, c}$. This means that the asymptotically optimal set derived by Metodiev et al. [41] for the case that the posterior is exactly normal is also going to be chosen by our method if this case occurs. Finally, if the posterior is not normal, our method chooses $\alpha$ such that the THAMES is well-behaved, in the sense that the Kolmogorov distance is going to penalise extreme outliers.

Even without this optimisation algorithm, we have observed that choosing $c = \sqrt{R+1}$ and $\alpha = 0.5$ lead to good results, which is in congruence with similar observations from Metodiev et al. [41] and Reichl [51].

## B.2 Choosing the grid for $\alpha$

We want to choose $\alpha$ such that it minimises the Kolmogorov distance

$$\sup_{\theta \in H_\alpha} \left| \chi_{M_\alpha, \hat{s}_\alpha}^2 (-\log(\pi(\theta)L(\theta))) - F(-\log(\pi(\theta)L(\theta))) \right|.$$

We do this by first choosing a grid for the different values of $\alpha$, and then, for each value on the grid, approximating the supremum via another grid search over $\theta$, respectively.
Regarding the choice of the grid for $\alpha$, it is important to note that the above function is itself an estimator of the distance between the true CDF of the normalised log posterior values, approximated by $F$, and the CDF of the $\chi^2$-distribution. For this reason, we have experienced that small values of $\alpha$ caused the function to be very volatile and unreliable, since almost all of the normalised log posterior values are excluded in this case. We have experienced that this can be avoided by only evaluating values of $\alpha$ that are larger than 0.2.

In addition, the function $F$ is constant for most values of $\theta$, excluding the values of the posterior sample itself, $\{\theta^{(1)}, \dots, \theta^{(T)}\}$. In turn, this leaves us with a possibility of $T$ different values of $\alpha$ to check, but we chose to thin the sequence by a step-size of 100, since we made the experience that this greatly reduces computation time with no loss in the efficiency of the THAMES. Regarding the choice of the grid for $\theta$, we chose the full grid, since the computation time was quite small, even for a large grid. If the chosen $\alpha$ is larger than 0.5, it is then truncated to 0.5, to be in accordance with Metodiev et al. [41].

### B.3 Numerically verifying that $\xi_{g_1} = \xi_{g_2}$ lies in $E_{\hat{\theta}, \hat{\Sigma}, c}$

$\xi_{g_1} = \xi_{g_2}$ lies in $E_{\hat{\theta}, \hat{\Sigma}, c}$ if, and only if, the minimum

$$ \min_{\theta: \xi_{g_1} = \xi_{g_2}} (\theta - \hat{\theta})^{\top} \hat{\Sigma}^{-1} (\theta - \hat{\theta}) \quad (22) $$

is smaller or equal to $c^2$. This is a quadratic optimisation problem, which can be solved in polynomial time since $\hat{\Sigma}^{-1}$ is positive definite, for example by using one of the algorithms given in Goldfarb and Idnani [23, 24]. We do this using the quadprog R package [61].

### B.4 Quickly computing the THAMES

The matrix $\Delta$ is constructed as follows: if it has been established that $\xi_{g_1} = \xi_{g_2}$ for $g_1, g_2$ and for some pair of values $(\xi_{g_1}, \xi_{g_2})$ that follows the ellipsoidal constraint given by $E_{\hat{\theta}, \hat{\Sigma}, c}$, then $\Delta_{g_1, g_2} = 0$ since $W(\xi_{g_1}) = W(\xi_{g_2})$ in that case. If not, then we test $W(\xi_{g_1}^{(j)}) < W(\xi_{g_2}^{(j)})$ with $(\xi_{g_1}^{(j)}, \xi_{g_2}^{(j)})$ being from the Monte Carlo sample that is simulated uniformly on $E_{\hat{\theta}, \hat{\Sigma}, c}$. If there exists a $j$ such that the condition is false, then $\Delta_{g_1, g_2}$ is set to 0. It is set to 1 otherwise. The law of large numbers guarantees $\Delta_{g_1, g_2}$ will be evaluated accurately if the size of the Monte Carlo sample $N$ is large enough.

The set $\Omega$ is constructed via an algorithm that returns all topological orderings on the graph defined by $\Delta$. An algorithm that does this in polynomial time with respect to the size of $\Omega$ is given in Knuth and Szwarcfiter [33]. We implemented a recursive version of this algorithm.

A very rough upper bound on the size of $\Omega$ is given by $G! / |I'(G)|!$, where $|I'(G)|$ denotes the number of vertices of the longest path in the graph generated by $\Delta$. This is because there are $|I'(G)| + 1$ possibilities of choosing the position of the first node that is not in $I'(G)$, $|I'(G)| + 2$ possibilities of choosing the position of the second node and so
on, and the order of the nodes within $I'(G)$ is fixed. Whenever this bound passes 50,000, we take this as a cue that the THAMES cannot be evaluated in a reasonable amount of time. In this case, we iteratively shrink the ellipse until the bound passes 50,000. This is done by repeatedly dividing $c$ by 2 and, if the ellipsoid ends up containing no samples due to this shrinkage, setting $\hat{\theta}$ to the mode of the second half of the posterior sample and continuing to shrink $c$. This makes sure that the ellipsoid always contains at least one point, the mode of the second half of the posterior sample. To make sure that $W$ can differentiate this point after $\hat{\theta}$ has been shifted, we train QDA only on those sample points that are within the ellipsoid in this case.

## Supplement C: details on the experiments

Note that the conditions of Theorem 5.5 are not met for the univariate model in Supplement C3, so the theorem was not used for univariate data. In addition, we have noticed that a large amount of the ellipsoid needed to be truncated in the case that its centre was close to the edge of the parameter space, which happened with variance and covariance matrix parameters. To avoid this, these parameters where mapped to the real line and the prior was adjusted via the absolute value of the determinant of the jacobian before evaluating the THAMES. The transformations used are described in Carpenter et al. [8].

### C.1 Simulations on univariate Gaussian mixtures

An analytic solution of the marginal likelihood was computed in the case that the number of parameters and data points is very small [the same idea was used in 16], such that the performance of the THAMES could be compared to other state-of-the-art estimators.

**Model** We consider a univariate Gaussian mixture model,

$$Y \in \mathbb{R}^n, \quad Y_1, \dots, Y_n | \mu, \sigma, \tau \stackrel{\text{i.i.d}}{\sim} \sum_{g=1}^{G} \tau_g N(\mu_g, \sigma_g^2), \qquad (23)$$

where only the parameters $\mu = (\mu_1, \dots, \mu_G)$ are unknown and $\tau = (\tau_1, \dots, \tau_G)$, $\sigma = (\sigma_1, \dots, \sigma_G)$ are known. They are known in the sense that we put a Dirac density on these parameters, where a Dirac density is a discrete density which is only equal to 1 at exactly 1 point. Let $C = (C_1, \dots, C_n) \in \mathcal{M}_{n \times G}(\{0, 1\})$ denote the matrix of hidden allocation vectors. An analytic solution was computed and a Gibbs sampler was constructed using the following theorem. It is proven in Supplement D.

**Theorem 5.6.** If the prior is set to the distribution with the following symmetric densities,

$$\pi(\mu) = \prod_{g=1}^{G} N(\mu_g; 0, 1), \quad \pi(\sigma_g) = \begin{cases} 1 & \sigma_g = 1, \\ 0 & \text{else,} \end{cases} \quad \pi(\tau_g) = \begin{cases} 1 & \tau_g = 1/G, \\ 0 & \text{else,} \end{cases} \qquad (24)$$
then the marginal likelihood is given by

$$ 
\begin{aligned}
p(Y) &= \sum_{C_1, \dots, C_n} p(Y|C_1, \dots, C_n)p(C_1, \dots, C_n) \\
&= \frac{1}{G^n} \sum_{C_1, \dots, C_n} \text{MVN}_n(Y; 0_n, I_n + O),
\end{aligned} 
\quad (25) $$

where O is a matrix such that $O_{i,j} = 1$ if $C_i = C_j$, $O_{i,j} = 0$ else. The conditional distributions of $\mu$ and the latent components $C_1, \dots, C_n$ are

$$ 
\begin{aligned}
\mu_g|C, Y &\stackrel{\text{i.i.d}}{\sim} \mathcal{N}\left(\frac{\sum_{i=1}^n C_{i,g} Y_i}{1 + \sum_{i=1}^n C_{i,g}}, \frac{1}{1 + \sum_{i=1}^n C_{i,g}}\right) \\
C_i|\mu, Y &\stackrel{\text{i.i.d}}{\sim} \text{Multinom}_G\left(1, \left(\frac{\exp\left(-\frac{1}{2}(Y_i - \mu_1)^2\right)}{\sum_{g=1}^G \exp\left(-\frac{1}{2}(Y_i - \mu_g)^2\right)}, \dots, \frac{\exp\left(-\frac{1}{2}(Y_i - \mu_G)^2\right)}{\sum_{g=1}^G \exp\left(-\frac{1}{2}(Y_i - \mu_g)^2\right)}\right)\right).
\end{aligned} 
$$

**Settings** In all of the following scenarios, we set $n = 10$ such that the formula given in Theorem 5.6 could be computed in a reasonable amount of time. In addition, even though our prior assumptions were symmetric, we wanted to test how our estimator performs in the case that the data is simulated from a distribution that is not symmetric, in the sense that its components have distinct means and proportions. For all scenarios, we simulated 50 different datasets from the fixed parameters $(\mu^*, \sigma^*, \tau^*)$. We set $\sigma^* = (1, \dots, 1)$, $\tau^* = (1/3, 2/3)$ if $G = 2$, $\tau^* = (2/6, 1/6, 3/6)$ if $G = 3$. The mean was varied,

$$ \mu^* = \begin{cases} (3, 3)^\mathrm{T} + (2(1 - \rho) + 1)(-1, 1)^\mathrm{T} & G = 2, \\ (3, 3, 3)^\mathrm{T} + (2(1 - \rho) + 1)(-1, 0, 1)^\mathrm{T} & G = 3, \end{cases} \quad \rho = 0, 1/2, 1. $$

This assures linear interpolation between means of $(0, 3, 6)^\mathrm{T}$, $(0, 6)^\mathrm{T}$, whose values are particularly well separated, and means of $(2, 3, 4)^\mathrm{T}$, $(2, 4)^\mathrm{T}$, which are less separate. Note that the true proportion parameters are not included in the class of models that was fitted. This is not a problem since we are not assessing the performance of estimators of the parameters, but of estimators of the marginal likelihood.

**Sampling procedure** A sample from the posterior of size 12,000 was produced via Gibbs sampling initialised in the midpoint of the data. The first 2,000 samples were discarded (burn-in) and the remaining $T = 10,000$ samples were used to compute the bridge sampling estimator, the estimator from Reichl [51] (using the same tuning parameters as the THAMES) as well as the THAMES on the relabelled posterior sample, where all estimators were defined with the prior from Equation (24).

We tested three scenarios: in the first, $G = 2$ and the data is simulated with $G = 2$. In the second, we tested our performance in the case of overfitting by computing the marginal likelihood for $G = 3$ on the data that was simulated with $G = 2$. Analogously, the third scenario corresponds to underfitting in the case that $G = 3$.
## C.2 Simulations on multivariate Gaussian mixtures

In this section, the performance of the THAMES was assessed in the case that the number of components G is large. The analytic expression given by Equation (25) is untractable in this case. However, it can be computed if the mixture components are sufficiently distinct from each other [the same idea was used in 16]. This is because in this case, the likelihood conditional on the allocations is equal to $p(Y|C) \approx 0$ for all but one allocation $C^0$ (and its label-switched versions, which result in the same likelihood). Thus, the marginal likelihood reduces to

$$p(Y) \approx G! \cdot p(Y|C^0) \cdot p(C^0), \qquad (26)$$

where the factor $G!$ comes from the different label-switched versions of $C^0$.

$C^0$ was estimated via the maximum aposteriori (MAP) estimator ($\hat{C}_{i,g}^0 = 1$ if $g = \text{argmax}_{\tilde{g}} \hat{z}_{i,\tilde{g}}$ and $\hat{C}_{i,g}^0 = 0$ else) and a prior distribution was chosen for which Expression (26) has an analytic solution and for which Theorem 5.5 holds. This is the case if the prior on $\tau$ is a Dirichlet distribution and the priors on $\xi_g$ are conjugate and independent. An exact expression of $p(Y|C^0)p(C^0)$ is given in Supplement D.

**Setting** To assure that Equation (23) can be applied, we needed to sufficiently separate the different components. We did this by choosing an absolute difference of 100 for all mean-component entries and a variance vector of ones for each respective component. We then assessed the validity of Equation (23) empirically: if all allocation vectors generated by the Gibbs sampler were identical, we assumed that there was indeed only one allocation vector with posterior probability not numerically equal to 0.

Let MVN$_d(\mu_g, \Sigma_g)$ denote the Gaussian law of dimension $d$ with mean $\mu_g$ and covariance matrix $\Sigma_g$. The mixture model is a multivariate Gaussian,

$$Y \in \mathbb{R}^{n \times d}, \quad Y_1, \dots, Y_n | \mu, \Sigma, \tau \stackrel{\text{i.i.d}}{\sim} \sum_{g=1}^{G} \tau_g \text{MVN}_d(\mu_g, \Sigma_g). \qquad (27)$$

We sampled 1 data set from Equation (27) with identical proportions $\tau_g = \frac{1}{G}$ and covariance matrices $\Sigma_g = I_d$, but different means $\mu_g = 100 \cdot g \cdot 1_d$, where $1_d$ denotes a vector of ones. We chose the priors

$$\begin{align*}
\tau &\sim \text{Dirichlet}(G; e_0, \dots, e_0), \quad \tau \in \mathbb{R}^G, \\
\mu_g | \Sigma_g &\sim \text{MVN}_d(\beta, \Sigma_g / \kappa_0), \quad \mu_1, \dots, \mu_G \in \mathbb{R}^m, \quad g = 1, \dots, G, \\
\Sigma_g &\sim \text{InvWis}(\phi_0, \Lambda), \quad \Sigma_1, \dots, \Sigma_G \in \mathbb{R}^{d \times d}, \quad g = 1, \dots, G, \\
e_0 &= 1, \\
\beta_{r,g} &= \frac{1}{n} \sum_{i=1}^n Y_{r,i}, \\
\kappa_0 &= 0.00001, \\
\phi_0 &= d,
\end{align*}$$
$$ \Lambda = (1 + \phi_0 + d) \frac{1}{G} \sum_{g=1}^{G} \tilde{\Sigma}_g, $$

where $\tilde{\Sigma}_1, \dots, \tilde{\Sigma}_G$ are estimates for the component-wise covariance matrices, computed via mclust.

Out of all of these values, only the value of $\Lambda$ is non-standard. We suggest this choice of $\Lambda$ since it assures that the prior mode of the covariance matrices is positioned within a range of plausible values of the covariance matrices. This needed to be done due to the large range of the data, which precludes more standard choices such as the sample covariance matrix, which was several orders of magnitude away from the true value in our case. We also tried to remedy this by scaling the sample covariance by $G^{-d/2}$, as suggested in Fraley and Raftery [13], but this does in turn result in a value of $\Lambda$ that is very close to the null-matrix, since e.g., for one of the following settings we set $G^{-d/2} = 15^{-5/2}$, which is fairly small. The hyperparameter that we suggest resulted in a posterior distribution that is fairly concentrated around the true value of each covariance matrix. Gibbs sampling was used, with the sampler being initialized in the allocation vector, using mclust. We tried 2 different settings, each chosen such that it exhibits similarity to a real dataset:

**Setting 1 (moderate d and G)** We simulated from a $d = 6$-dimensional model with $G = 5$ components and $n = 200$ data points. Thus, the dimension of the parameter space is $R = d - 1 + Gd + Gd(d+1)/2 = 139$.

**Setting 2 (large d and G)** We simulated from a $d = 5$-dimensional model with $G = 15$ components and $n = 345$ data points. It is of course possible to use the THAMES to determine the marginal likelihood of the full model, such as in the first setting. However, the number of free parameters is very large, and, as is typical in such high dimensional models, it can be decreased by imposing additional constraints on the covariance matrices of the model. Gibbs sampling is also possible in this case, and a variety of alternative parsimonious geometric models with their associated Gibbs sampling algorithm are presented in Bensmail et al. [5]. We choose the following: we impose each covariance matrix to be diagonal with a variance vector $\sigma_{g,1}^2, \dots, \sigma_{g,d}^2$ for each component. We choose independent inverse gamma priors $\text{InvGamma}(\phi_0, \lambda_r)$ for each of the variances, with $\phi_0 = 2, \lambda_r = 2 \frac{1}{G} \sum_{g=1}^{G} \hat{\sigma}_{g,r}^2$ analogous to the multivariate case.

### C.3 Univariate datasets (galaxy, enzyme, acidity)

The performance of the THAMES was evaluated for the datasets used in Richardson and Green [52], namely the “acidity”, “enzyme”, and “galaxy” dataset. This is because results of other, similar marginal likelihood estimators on the same datasets are available in Celeux et al. [9, Section 2.3.2]. We compared the THAMES to these estimators.
<table><caption>Table 1: The log-THAMES and criterion of overlap (CO) for different values of G for the galaxy dataset; the CO is maximized by G = 3, the marginal likelihood by G = 6. Numbers were rounded to the first decimal place.</caption><thead><tr><th></th><th>G=2</th><th>G=3</th><th>G=4</th><th>G=5</th><th>G=6</th><th>G=7</th><th>G=8</th></tr></thead><tbody><tr><th>log-THAMES</th><td>-235.2</td><td><strong>-226.7</strong></td><td>-226.0</td><td>-225.6</td><td><strong>-225.4</strong></td><td>-226.9</td><td>-226.4</td></tr><tr><th>CO</th><td>2</td><td><strong>3</strong></td><td>2</td><td>1</td><td><strong>0</strong></td><td>-1</td><td>-2</td></tr><tr><th></th><th>G=9</th><th>G=10</th><th>G=11</th><th>G=12</th><th>G=13</th><th>G=14</th><th>G=15</th></tr><tr><th>log-THAMES</th><td>-275.0</td><td>-307.6</td><td>-336.1</td><td>-329.7</td><td>-363.5</td><td>-375.2</td><td>-387.6</td></tr><tr><th>CO</th><td>-5</td><td>-4</td><td>-7</td><td>-6</td><td>-9</td><td>-10</td><td>-9</td></tr></tbody></table>

**Settings** All datasets were fitted to a univariate Gaussian mixture model,

$$Y \in \mathbb{R}^n, \quad Y_1, \dots, Y_n | \mu, \sigma, \tau \stackrel{\text{i.i.d}}{\sim} \sum_{g=1}^{G} \tau_g \mathcal{N}(\mu_g, \sigma_g^2),$$

with hierarchical priors

$$\tau \sim \text{Dirichlet}(G; 1, \dots, 1), \quad \tau \in \mathbb{R}^G$$

$$\mu_g \sim \mathcal{N}(x_0, \lambda), \quad g = 1, \dots, G,$$

$$\sigma_g^2 | \zeta \sim \text{InvGamma}(2, \zeta), \quad g = 1, \dots, G,$$

$$\zeta \sim \text{Gamma}(0.2, 10/\lambda^2),$$

where $x_0$ and $\lambda$ are the midpoint and length of the observation interval, respectively.

**Sampling procedure** We used the same sampling procedure as Celeux et al. [9, Section 2.3.2]: Gibbs sampling was conducted to obtain a posterior sample of size 12,000 with a burn-in of 2,000, yielding a sample of size T = 10,000. The values of the galaxy dataset were divided by 100 to be in accordance with Celeux et al. [9, Section 2.3.2]. The bayesmix package [26] was used to implement the Gibbs sampler.

**Additional results** The THAMES was used on a posterior sample of size T = 100,000 with a burn-in of 2,000 to compute the marginal likelihood and criterion of overlap (CO) for different values of G. The results are shown in Table 1. The marginal likelihood is maximised by G = 6, while the CO is maximised by G = 3. As pointed out in the main article, this difference is likely due to different meanings of the CO and marginal likelihood: the number of true components is estimated to be 6, while the number of distinguishable components is estimated to be 3.

## C.4 Multivariate datasets (Swiss banknotes, BUPA liver)

**Settings** The likelihood and priors correspond to those described in Supplement C2.
<table><caption>Table 2: The log-THAMES and criterion of overlap (CO) for different values of G for the Swiss banknote dataset; both model choice criteria were maximised by G = 3. Numbers were rounded to the second decimal place.</caption><thead><tr><th></th><th>G=2</th><th>G=3</th><th>G=4</th><th>G=5</th></tr></thead><tbody><tr><td>log-THAMES</td><td>-926.78</td><td><strong>-909.49</strong></td><td>-922.93</td><td>-941.87</td></tr><tr><td>CO</td><td>2</td><td><strong>3</strong></td><td>2</td><td>0</td></tr></tbody></table>

<table><caption>Table 3: The log-THAMES and criterion of overlap (CO) for different values of G for the BUPA liver disorders dataset; the CO and log-THAMES were both maximised by G = 4. Numbers were rounded to full integers.</caption><thead><tr><th></th><th>G=2</th><th>G=3</th><th><strong>G=4</strong></th><th>G=5</th><th>G=6</th><th>G=7</th><th>G=8</th></tr></thead><tbody><tr><td>log-THAMES</td><td>-6689</td><td>-6655</td><td><strong>-6648</strong></td><td>-6668</td><td>-6695</td><td>-6726</td><td>-6764</td></tr><tr><td>CO</td><td>2</td><td>3</td><td><strong>4</strong></td><td>3</td><td>0</td><td>-1</td><td>-2</td></tr><tr><td></td><td>G=9</td><td>G=10</td><td>G=11</td><td>G=12</td><td>G=13</td><td>G=14</td><td>G=15</td></tr><tr><td>log-THAMES</td><td>-6807</td><td>-6856</td><td>-6910</td><td>-6971</td><td>-7037</td><td>-7109</td><td>-7187</td></tr><tr><td>CO</td><td>-3</td><td>-4</td><td>-5</td><td>-6</td><td>-7</td><td>-8</td><td>-9</td></tr></tbody></table>

**Sampling procedure** Gibbs sampling was used, with the sampler being initialised in the allocation vector, using mclust. We used 12,000 iterations with a burn-in of 2,000, resulting in a posterior sample of size T = 10,000.

**Additional results** The THAMES was used to estimate the log marginal likelihood for G = 2 to G = 5 for the Swiss banknote dataset (Table 2), G = 2 to G = 15 for the BUPA liver disorders dataset (Table 3). G = 3 was found to maximise the CO and marginal likelihood for the Swiss banknote dataset.

The labels “genuine” and “counterfeit” are known for the banknote dataset. Table 4 shows the confusion matrix with respect to these labels, where we clustered the data points for G = 3, since this value was estimated to maximise the marginal likelihood and CO. All but one of the datapoints are assigned the “counterfeit” label if belonging to component 1 or 3, the “genuine” label if belonging to component 2.

The component means estimated for the BUPA liver dataset seem to heavily differ around the variable “gammagt”, denoting the results of a blood test on gamma-glutamyl transpeptidase. This is indicated by the estimates of the posterior mean and standard deviation, computed on the relabelled posterior sample (Table 5).
<table><caption>Table 4: Confusion matrix for the Swiss banknote dataset for G = 3, the value which maximises the marginal likelihood; g denotes the label of component g; components 1 and 3 seem to correspond to counterfeits, while component 2 seems to correspond to genuine banknotes.</caption><thead><tr><th rowspan="2">g</th><th colspan="2">banknote type</th></tr><tr><th>genuine</th><th>counterfeit</th></tr></thead><tbody><tr><td>1</td><td>1</td><td>16</td></tr><tr><td>2</td><td>99</td><td>0</td></tr><tr><td>3</td><td>0</td><td>84</td></tr></tbody></table>

<table><caption>Table 5: Posterior mean estimate +/- the posterior standard deviation estimate of the different component means for the BUPA liver dataset for G = 4, the value which maximised the marginal likelihood; the means correspond to the results of blood tests on the mean corpuscular volume (mcv), alkaline phosphotase (alkphos), alanine aminotransferase (sgpt), aspartate aminotransferase (sgot) and gamma-glutamyl transpeptidase (gammagt). g denotes the label of component g; the component means seem to strongly differ around the level of gammagt. Numbers were rounded to the first decimal place.</caption><thead><tr><th rowspan="2">g</th><th colspan="5">blood tests</th></tr><tr><th>mcv</th><th>alkphos</th><th>sgpt</th><th>sgot</th><th><strong>gammagt</strong></th></tr></thead><tbody><tr><td>1</td><td>90.3 ± 0.5</td><td>75.2 ± 2.2</td><td>29.4 ± 1.1</td><td>24.4 ± 0.6</td><td><strong>37.8 ± 2.4</strong></td></tr><tr><td>2</td><td>90.9 ± 0.8</td><td>75.8 ± 2.4</td><td>52.0 ± 3.2</td><td>36.4 ± 1.9</td><td><strong>91.3 ± 11.1</strong></td></tr><tr><td>3</td><td>89.6 ± 0.3</td><td>64.4 ± 1.2</td><td>21.1 ± 0.6</td><td>19.6 ± 0.4</td><td><strong>17.7 ± 0.9</strong></td></tr><tr><td>4</td><td>95.6 ± 2.1</td><td>65.3 ± 9.2</td><td>122.0 ± 26.1</td><td>66.2 ± 8.8</td><td><strong>133.3 ± 22.9</strong></td></tr></tbody></table>
# Supplement D: exact and approximate density computations

## D.1 Exact marginal likelihood for a small dataset

*Proof of Theorem 5.6.* By the law of total probability

$$p(Y) = \sum_{C_1, \dots, C_n} p(Y|C_1, \dots, C_n)p(C_1, \dots, C_n). \qquad (28)$$

The prior on the clusters is given by

$$p(C_1, \dots, C_n) = \prod_{i=1}^{n} \prod_{g=1}^{G} \left(\frac{1}{G}\right)^{C_{i,g}} = G^{-\sum_{g=1}^{G} \sum_{i=1}^{n} C_{i,g}}, \qquad (29)$$

so we only need to solve for $p(Y|C_1, \dots, C_n)$. Notice that $Y$ can be written as

$$Y = \left( \prod_{g=1}^{G} \mu_{g}^{C_{1,g}}, \dots, \prod_{g=1}^{G} \mu_{g}^{C_{n,g}} \right)^{\top} + \varepsilon$$

where $\varepsilon \sim \text{MVN}(0, I_n)$ is independent from $\mu$. Thus, the marginal distribution of $Y$ given $C_1, \dots, C_n$ is given by

$$p(Y_i|C_1, \dots, C_n) = \text{MVN}_n(Y_i; 0_n, I_n + O). \qquad (30)$$

Plugging (29) and (30) into (28) gives the result.

We are now going to derive the distribution of $\mu$ given $Y$ and $C_1, \dots, C_n$. Let

$$\tilde{\mu} = \left( \prod_{g=1}^{G} \mu_{g}^{C_{1,g}}, \dots, \prod_{g=1}^{G} \mu_{g}^{C_{n,g}} \right)^{\top}.$$

Then, by Bayes' rule

$$\begin{align*}
p(\mu|Y, C_1, \dots, C_n) &\propto p(\mu)p(C_1, \dots, C_n)p(Y|\mu, C_1, \dots, C_n) \\
&\propto p(\mu)p(Y|\mu, C_1, \dots, C_n) \\
&= \left( \prod_{g=1}^{G} \mathcal{N}(\mu_g; 0, 1) \right) \text{MVN}_n(Y; \tilde{\mu}, I_n) \\
&\propto \exp \left( -\frac{1}{2} \sum_{g=1}^{G} \mu_g^2 \right) \exp \left( -\frac{1}{2} (Y - \tilde{\mu})^\top (Y - \tilde{\mu}) \right) \\
&\propto \exp \left( -\frac{1}{2} \sum_{g=1}^{G} \mu_g^2 \right) \exp \left( -\frac{1}{2} (\tilde{\mu}^\top \tilde{\mu} - 2Y^\top \tilde{\mu}) \right)
\end{align*}$$
$$ 
\begin{aligned}
&= \exp \left( -\frac{1}{2} \sum_{g=1}^{G} \mu_g^2 - \frac{1}{2} \sum_{i=1}^{n} (\tilde{\mu}_i^2 - 2\tilde{\mu}_i Y_i) \right) \\
&= \exp \left( -\frac{1}{2} \sum_{g=1}^{G} \mu_g^2 - \frac{1}{2} \sum_{g=1}^{G} \sum_{i:C_{i,g}=1} (\mu_g^2 - 2\mu_g Y_i) \right).
\end{aligned}
$$

We recognise that the density factorises into the densities

$$ 
\begin{aligned}
p(\mu_g|Y, C_1, \dots, C_n) &\propto \exp \left( -\frac{1}{2}\mu_g^2 - \frac{1}{2}\sum_{i:C_{i,g}=1} (\mu_g^2 - 2\mu_g Y_i) \right) \\
&= \exp \left( -\frac{1}{2} \left( \left(1 + \left(\sum_{i=1}^n C_{i,g}\right)\right)\mu_g^2 - 2\mu_g \left(\sum_{i=1}^n C_{i,g}Y_i\right) \right) \right) \\
&\propto \mathcal{N} \left( \mu_g; \frac{\sum_{i=1}^n C_{i,g}Y_i}{1 + \sum_{i=1}^n C_{i,g}}, \frac{1}{1 + \sum_{i=1}^n C_{i,g}} \right).
\end{aligned}
$$

Analogously,

$$ 
\begin{aligned}
p(C_i|Y, \mu, C_{-i}) &\propto p(C_i)p(Y|\mu, C_1, \dots, C_n) \\
&\propto \exp \left( -\frac{1}{2} \sum_{g=1}^{G} \sum_{i:C_{i,g}=1} (\mu_g - Y_i)^2 \right) \\
&\propto \prod_{g=1}^{G} \exp \left( -\frac{1}{2} C_{i,g} (\mu_g - Y_i)^2 \right) \\
&= \prod_{g=1}^{G} \left( \exp \left( -\frac{1}{2} (\mu_g - Y_i)^2 \right) \right)^{C_{i,g}},
\end{aligned}
$$

so the posterior probabilities on $C_i$ given $\mu, Y$ and all other allocation vectors are the normalization of the vectors

$$ \left( \exp \left( -\frac{1}{2} (\mu_1 - Y_i)^2 \right), \dots, \exp \left( -\frac{1}{2} (\mu_G - Y_i)^2 \right) \right). $$

<figure><img src="image_12.png" /></figure>

## D.2 Highly accurate marginal likelihood approximation for an extremely well separated dataset

We know that the marginal likelihood is given by

$$ Z(G) = \sum_{C_1, \dots, C_n} p(Y|C_1, \dots, C_n)p(C_1, \dots, C_n) \approx G! p(Y|C^0)p(C^0) $$
if $p(Y|C^0) \simeq 0$ for all but one allocation matrix $C^0$ (and its permutations). The expression $p(C^0)$ is equal to [see for example 16, 28]

$$ \frac{\Gamma(G)}{\Gamma(G+n)} \prod_{g=1}^{G} \Gamma \left( 1 + \sum_{i=1}^{n} C_{i,g}^{0} \right). $$

First of all, let us note that (denoting $\mu = (\mu_1, \dots, \mu_G)$, $\Sigma = (\Sigma_1, \dots, \Sigma_G)$)

$$ p(Y|C^0, \tau, \mu, \Sigma) = p(Y|C^0, \mu, \Sigma), $$

since the proportions are rendered redundant by the knowledge of the component membership. $p(Y|C^0)$ can thus be further decomposed into

$$ \begin{aligned} p(Y|C^0) &= \int p(Y|C^0, \mu, \Sigma) \pi(\mu, \Sigma) d(\mu, \Sigma) \\ &= \int \prod_{i=1} p(Y_i|C^0, \mu, \Sigma) \pi(\mu, \Sigma) d(\mu, \Sigma) \\ &= \prod_{g=1}^{G} \int p_g(Y|\mu_g, \Sigma_g) \pi_g(\mu_g, \Sigma_g) d(\mu_g, \Sigma_g) = \prod_{g=1}^{G} Z_g(G), \end{aligned} $$

where $\pi_g$ denotes the prior on $(\mu_g, \Sigma_g)$ and

$$ p_g(Y|\mu_g, \Sigma_g) = \prod_{i:C_{i,g}=1} p(Y_i|C^0, \mu_g, \Sigma_g). $$

This is because $Y_1, \dots, Y_n$ are independent given $C^0$, $\mu$, $\Sigma$ and $(\mu_1, \Sigma_1), \dots, (\mu_G, \Sigma_G)$ are a priori independent.

Because $p_g$ denotes a likelihood on those data points $Y_i$ for which $C_{i,g} = 1$ with length $n_g = \sum_{i=1}^n C_{i,g}$, $p_g$ and $\pi_g$ can be interpreted as the prior and likelihood, respectively, of a model with a Gaussian inverse Wishart distribution and the above integrals $Z_g(G)$ can be interpreted as the respective marginal likelihoods of each of these $G$ different models. These marginal likelihoods are tractable due to conjugacy and Bayes rule:

$$ \frac{\int p_g(Y|\mu_g, \Sigma_g) \pi_g(\mu_g, \Sigma_g) d(\mu_g, \Sigma_g) = \frac{p_g(Y|\mu_g, \Sigma_g) \pi_g(\mu_g, \Sigma_g)}{p_g(\mu_g, \Sigma_g|Y)} \left( \prod_{i:C_{i,g}=1} \text{MVN}_d(Y_i; \mu_g, \Sigma_g) \right) \text{InvWis}(\Sigma_g; \phi_0, \Lambda) \text{MVN}_R(\mu_g; \beta, \Sigma_g/\kappa_0) \text{InvWis}(\Sigma_g; d, \Lambda_{n,g}) \text{MVN}_d(\mu_g; \beta_{n,g}, \Sigma_g/\kappa_{n,g})}{\text{InvWis}(\Sigma_g; d, \Lambda_{n,g}) \text{MVN}_d(\mu_g; \beta_{n,g}, \Sigma_g/\kappa_{n,g})} $$

This holds for any $\mu_g, \Sigma_g$, where $p_g(\mu_g, \Sigma_g|Y)$ denotes the posterior distribution with respect to the $g$'th model.

The posterior hyperparameters are for example given in Gelman et al. [19]. Plugging in and factoring out $\mu_g, \Sigma_g$ gives the following (rather long) expression, where $\Gamma_d$ denotes
the multivariate Gamma function (the posterior hyperparameter $\beta_{n,g}$ ends up being irrelevant since it factors out in the calculation but is still given for completeness):

$$
\begin{align*}
p(Y) &\approx G! \cdot p(Y|C^0) \cdot p(C^0) \\
p(C^0) &= \frac{\Gamma(G)}{\Gamma(G+n)} \prod_{g=1}^{G} \Gamma \left( 1 + \sum_{i=1}^{n} C_{i,g}^0 \right) \\
p(Y|C^0) &= \prod_{g=1}^{G} Z_g(G) \\
\log(Z_g(G)) &= (\phi_0/2) \log(\det(\Lambda)) - (\phi_0 d/2) \log(2) - \log(\Gamma_d(\phi_0/2)) \\
&\quad + (\phi_{n,g} d/2) \log(2) + \Gamma_d(\phi_{n,g}/2) + (-\phi_{n,g}/2) \log(\det(\Lambda_{n,g})) \\
&\quad - n_g d \log(\sqrt{2\pi}) - 2 \log(\kappa_{n,g}) \\
\phi_{n,g} &= \phi_0 + n_g \\
\Lambda_{n,g} &= \Lambda + \frac{\kappa_0 n_g}{\kappa_0 + n_g} (\bar{Y}_g - \beta)(\bar{Y}_g - \beta)^\top + \sum_{i:C_{i,g}=1} (Y_i - \bar{Y}_g)(Y_i - \bar{Y}_g)^\top \\
\beta_{n,g} &= \frac{\kappa_0}{\kappa_0 + n_g} \beta + \frac{n_g}{n_g + \kappa_0} \bar{Y}_g \\
\bar{Y}_{g,r} &= \frac{1}{n_g} \sum_{i:C_{i,g}=1} Y_{i,r} \\
\kappa_{n,g} &= \kappa_0 + n_g \\
n_g &= \sum_{i:C_{i,g}=1} 1
\end{align*}
$$