{"id":"0326fd9c-5fae-4e46-a851-c9ec4983ea3f","arxiv_id":"2504.21812","paper_version":3,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":5,"one_line_summary":"A symmetrized, ordering-constrained THAMES estimator accurately computes marginal likelihoods for univariate and multivariate Gaussian mixtures, demonstrated up to 15 components.","lead":"An adapted THAMES estimator computes marginal likelihoods for mixture models with many components, using only relabelled MCMC samples. It offers a generic, label-switching-invariant tool for choosing the number of clusters, with an R package.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Equations (1) and (5) put πL in the numerator, but Supplement A.2's RIS derivation requires dividing by πL; the displayed estimator would not converge to the reciprocal marginal likelihood, so the central definition and its proof are inconsistent.","rationale":"The reader's weakest assumption concerned the unquantified effect of dependence and of estimating the truncation set from the data. While those are legitimate secondary concerns, the most load-bearing issue is internal: the estimator as displayed in Equations (1) and (5) is not the estimator whose theoretical properties are proved in Supplement A.2. The proof relies on the reciprocal importance sampling identity E_p[h(θ)/(π(θ)L(θ))] = 1/Z for any density h, which requires dividing by the unnormalised posterior density, not multiplying by it. The main text's displayed summand 1/V(B) * π(θ)L(θ) has the wrong power of the unnormalised posterior density. If the equations are taken at face value, the estimator is biased and not consistent for the reciprocal marginal likelihood. The simulations and the Supplement's algebra strongly suggest the implementation actually uses the reciprocal form, so the practical method may be sound; however, the manuscript as written contains a mathematical inconsistency in its central definition, and this must be corrected before the theoretical claims can be evaluated. This does not change the reader's CONDITIONAL verdict, because the condition should now explicitly include fixing this inconsistency and confirming that the code matches the corrected equation.","tokens_in":32209,"tokens_out":20265,"duration_ms":212988,"concrete_test":"Inspect the thamesmix R package source code and verify how the unnormalised posterior density enters the estimator: does the implementation sum 1/(V * unnormalized_posterior) or unnormalized_posterior / V? Alternatively, run a controlled simulation with a known marginal likelihood Z: compute the quantity in Equation (5) exactly as printed (πL in the numerator) and the RIS version (1/(πL) in the numerator) on the same posterior draws. The version that converges to 1/Z identifies which expression the proofs and code actually support; the other version will not. Also independently re-derive Equation (5) from the RIS identity in Supplement A.2; the two displayed forms cannot both equal the same estimator.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The paper's central estimator is defined in Equation (5) as a sum over permutations of 1/V(B) times π(θ)L(θ) for relabelled draws inside B. Taken literally, this is the arithmetic mean of the unnormalised posterior density over a truncated region, not a harmonic mean or a reciprocal importance sampling estimator. Its expectation under the posterior is (1/(Z V(B))) ∫_B (πL)^2 dθ, which is not 1/Z. Supplement A.2, however, proves unbiasedness and consistency by identifying the estimator as an RIS estimator of the form (1/(T/2)) Σ h(θ)/(π(θ)L(θ)), and the final algebra in that Supplement gives 1/(V(B) π(θ)L(θ)) in the summand. Thus the displayed definition in the main text and the proof in Supplement A.2 cannot both be correct: either the equations are typographically wrong (missing the reciprocal of πL) or the proof establishes properties of a different estimator. This is more fundamental than the independence or truncation-set issues, because it affects even the idealized estimator with fixed B and known volume.","agreement_with_reader":"disagree"},"referee_report":{"model":"deepseek-v4-flash","summary":"This paper proposes a symmetrized version of the truncated harmonic mean estimator (THAMES) for the marginal likelihood of mixture models. The estimator averages a Reichl-type volume-truncated harmonic mean over all G! label permutations of relabelled posterior draws, and uses an ordering constraint based on quadratic discriminant analysis and overlap graphs to reduce the cost of the permutation sum. The authors prove in a supplement symmetry, unbiasedness, consistency, asymptotic normality, and a variance-optimality property for an idealized RIS formulation, provide an empty-component correction based on Nobile's identity, and evaluate the estimator on simulated and real datasets, including comparisons with bridge sampling and Reichl's estimator. The typeset estimator in Eq. (5) has π(θ)L(θ) in the denominator, matching the RIS derivation in Supplement A.2; the alleged reciprocal inconsistency does not survive a close reading.","tokens_in":32355,"tokens_out":12209,"duration_ms":136755,"significance":"If the theoretical claims were established for the actual algorithm, this would be a valuable addition to the toolkit for model selection in mixture models: it is generic, uses only relabelled posterior draws, and the reported computations appear feasible for G=15. The paper includes useful analytic benchmarks (Theorem 5.6 and Supplement D), a detailed supplement with proofs, and an R package. The central gap is that the implemented estimator differs from the one analyzed: the truncation set, its volume, and the ordering constraint are estimated from the same MCMC sample used in the average, so the stated guarantees are not proven for the actual procedure. The work is therefore promising but requires a substantial revision to close the gap between the ideal estimator and the implementation.","major_comments":[{"comment":"The RIS analysis in Supplement A.2 treats h(θ) = (1/G!)∑_o 1_B(P_o θ)/V(B) as fixed and independent of the sample used in the average. In the implementation, the quadratic discriminant classifier (using component means and covariances estimated on the second half), the ordering set Ω, the threshold q̂_α, and the Monte Carlo volume V̂(B) are all estimated from the same MCMC output. The first/second-half split covers only θ̂ and Σ̂. Hence the estimator actually computed is not the RIS estimator whose unbiasedness and CLT are proved; the data-adaptive construction requires either an independent evaluation sample or a separate theorem quantifying the effect of estimating these inputs.","section":"Section 3.2; Supplement A.2"},{"comment":"Equation (5) is written with the exact volume V(B), but the algorithm replaces it by the Monte Carlo estimate from Eq. (3) with N=T. The proofs in Supplement A.2 use the exact volume. Because E[1/V̂] does not equal 1/V in general, the finite-sample unbiasedness claim fails for the implemented ratio estimator, and consistency of the ratio needs a separate argument, for example uniform convergence of V̂ and independence from the inner average.","section":"Eq. (3) and Eq. (5)"},{"comment":"The variance-optimality theorems (Theorem 5.1 in Supplement A.1 and Theorem 5.3 in Supplement A.3) explicitly assume independent posterior draws, and the main text acknowledges this is approximate. Supplement A.2 states that the symmetry, unbiasedness, consistency, and asymptotic normality results make no such assumptions, but it does not state the Markov-chain conditions used. Since the headline claims are for MCMC samples, the paper should state explicit ergodicity or mixing conditions and clarify that 'unbiased' is an iid-sample or asymptotic statement.","section":"Supplement A.1 and A.2"},{"comment":"The identity replacing the G! permutation sum by the sum over topological orderings Ω is exact only when the adjacency matrix Δ is known exactly. The algorithm estimates Δ by checking W-inequalities on a finite Monte Carlo sample from the ellipsoid, as described in Supplement B.4; a misclassification changes the value of the estimator, not merely its computational cost. No error bound or consistency statement is given for this approximation.","section":"Theorem 5.4 / Eq. (8)"}],"minor_comments":[{"comment":"The typesetting of Equations (1) and (5) should make the reciprocal explicit, for example by writing 1/(V(B)π(θ)L(θ)), since the current line break invites misreading as π(θ)L(θ) in the numerator.","section":"Eq. (1) and Eq. (5)"},{"comment":"The text refers to 'Equation (17)' when discussing computation of the THAMES, but in the main text the relevant display is Equation (8); Equation (17) is in the supplement.","section":"Section 3.2"},{"comment":"The statement that the posterior density is bounded below whenever h is not 0 should explicitly invoke the threshold q̂_α > 0 on B, rather than relying on the boundedness of the support alone.","section":"Supplement A.2"},{"comment":"The captions of Figures 2 and 8 describe counting points in the ellipsoid, but the relevant condition is membership in B = E ∩ {πL > q̂_α}; please clarify the captions.","section":"Figures 2 and 8"},{"comment":"The comparison with bridge sampling would be more informative if the paper specified whether bridge sampling was applied to the relabelled or the raw posterior sample, since non-symmetric estimators applied to relabelled draws are not automatically estimators of the same quantity.","section":"Figure 4"}],"recommendation":"major_revision","confidential_remarks":"The manuscript depends heavily on the authors' earlier THAMES paper [41] and on the supplemental document [40]; the main text is not fully self-contained for the proofs of the central claims. The use of ECR relabelling also sits uneasily with the claim that the method does not require hidden allocation vectors, and this should be clarified in revision. I recommend asking the authors to close the gap between the idealized estimator analyzed in the supplement and the data-adaptive procedure actually implemented."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"The stress-test note does not hold up on close reading. Equations (1) and (5) put πL in the denominator: the expression is (1/V(B)) / (π(θ)L(θ)), not πL in the numerator. I checked the main text and Supplement A.2; they are consistent with each other and with the usual reciprocal importance sampling identity. So the core definition is fine.\n\nWhat is actually new here is a useful combination: symmetrisation over label permutations (borrowed from Berkhof) plus Reichl's ellipsoid-HPD truncation, made computationally tractable via the overlap-graph ordering and QDA-based constraint. The CO criterion is a genuinely new auxiliary tool for choosing the number of distinct components. The univariate simulation, with analytically known marginal likelihoods across 50 datasets and three misspecification scenarios, is solid evidence: THAMES tracks the truth where bridge sampling and the Reichl estimator go off. The multivariate evidence is thinner — one well-separated artificial dataset per setting — but it does show the method scales to G=15, and the real-data examples are sensible.\n\nThe soft spots are real but not fatal. The consistency and variance proofs treat B and its volume as fixed, while the implementation estimates them from the data; the split-sample construction is approximate, and the independence assumption in the optimality theorem is acknowledged but not quantified. The QDA-based ordering relies on approximate normality, which is reasonable for well-identified mixtures but not guaranteed. The 'arbitrarily high number of components' claim in the abstract is hyperbolic; the paper's own discussion notes that high overlap forces a smaller ellipsoid and higher variance. These are addressable with sharper assumptions or more extensive simulations, and they do not undermine the central estimator.\n\nThe paper is honest, clearly written, and ships an R package. It is a methodological contribution that applied Bayesians would find genuinely useful for model selection in multivariate mixtures. I would send it to peer review; a good referee should push on the dependence and estimated-truncation issues, but the work deserves the time.","headline":"The stress-test concern about the estimator's numerator/denominator is a misreading; Equations (1) and (5) are the standard reciprocal form, and the paper's central idea is sound, though the proofs leave the estimated truncation region and MCMC dependence unquantified.","tokens_in":33020,"tokens_out":2372,"would_cite":true,"duration_ms":25782,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["62F15","62-04","62F12"],"pacs":[],"model":"deepseek-v4-flash","headline":"The paper adapts the truncated harmonic mean estimator to multivariate mixture models by symmetrising over component labels and making the permutation sum tractable through a QDA-based ordering constraint.","keywords":["Bayesian statistics","bridge sampling","harmonic mean estimator","label switching","marginal likelihood estimation","mixture modelling","truncated harmonic mean estimator","overlap graphs"],"falsifier":"On a small mixture with an exactly computable marginal likelihood, take one long MCMC run and compute THAMES on the raw chain, on a chain thinned by a factor of 10, and on a set of independent draws; if the three estimates differ by more than their reported Monte Carlo standard errors, the independence and split-sample assumptions are violated in a way that matters.","tokens_in":31882,"feed_emoji":"🧮","tokens_out":10175,"duration_ms":93219,"temperature":0.7,"pith_summary":"The paper claims that a symmetric version of the truncated harmonic mean estimator (THAMES) computes the marginal likelihood of a multivariate mixture model from nothing more than the unnormalised log-posterior and a relabelled posterior parameter sample. The estimator is claimed to be consistent, asymptotically normal and of finite variance, and invariant to label switching, so no simulation from hidden allocation vectors or randomly permuted samples is needed. The naive definition averages over all $G!$ label permutations, which would be intractable, but the paper shows that the sum reduces to a small set of topological orderings defined by a quadratic-discriminant-analysis-based ordering score, so the cost stays manageable even for $G=15$. In simulations with analytically known marginal likelihoods the estimator converges to the truth while bridge sampling can be off by a factor of $G!$, and on real datasets it selects plausible numbers of components. If the claims hold, applied Bayesians obtain a simple and generic tool for model choice in high-dimensional mixture models.","feed_headline":"Symmetrised harmonic mean beats bridge sampling on mixture evidence","feed_subtitle":"It matches exact values for 15-component Gaussian mixtures where bridge sampling is off by a factor of G!","key_machinery":"The load-bearing object is the symmetric THAMES estimator of Equation (5), viewed as a reciprocal importance sampling estimator whose proposal is a mixture of uniform distributions on the symmetrised truncation set. Its tractability comes from the partial order defined by the ordering score $W(\\xi_g) = \\hat{g}(\\xi_g|I(G)) + 1 - \\hat{w}_{\\hat{g}(\\xi_g|I(G))}(\\xi_g|I(G))$, built from quadratic discriminant analysis restricted to a maximum independent set of non-overlapping components; the inequalities $W(\\xi_{g_1}) < W(\\xi_{g_2})$ are verified on the ellipsoid $E_{\\hat{\\theta},\\hat{\\Sigma},c}$, producing an adjacency matrix whose topological orderings are exactly the nonzero terms of the estimator. The graph machinery also yields the criterion of overlap, a separate estimate of the number of distinguishable components. In short, the ordering constraint converts an exponential sum over label permutations into a small sum over graph-theoretic orderings without changing the value of the estimator.","core_discovery":"The central discovery is that the mixture-adapted THAMES, defined as a uniform average over label permutations of a truncated reciprocal-importance-sampling estimator on the set $B_{\\hat{\\theta},\\hat{\\Sigma},c,\\alpha} = E_{\\hat{\\theta},\\hat{\\Sigma},c} \\cap \\{\\theta : \\pi(\\theta)L(\\theta) > \\hat{q}_\\alpha\\}$, is both theoretically sound and computationally tractable. The paper proves unbiasedness on the reciprocal scale, consistency, asymptotic normality and finite variance, and shows that this symmetrisation is optimal among mixtures of truncated harmonic mean estimators. The computational breakthrough is Equation (8): because the ordering constraint $W(\\xi_1) \\le \\cdots \\le W(\\xi_G)$ holds throughout the truncation region for most pairs of components, the sum over $G!$ permutations equals a sum over the topological orderings of a directed graph, and only those orderings need to be evaluated. The score $W$ is derived from quadratic discriminant analysis on the relabelled posterior sample, restricted to a maximum independent set of non-overlapping components found via overlap graphs. An additional identity lets the estimator handle empty components by reducing to a lower-dimensional model, which matters when many components are fitted and some capture no data points.","pith_inferences":["Because the estimator only needs a relabelled posterior and the log posterior, a natural extension is to models beyond Gaussian mixtures, such as t-mixtures or mixtures of experts; the ordering-score construction would need re-testing when the relabelled posterior is not approximately normal.","The split-sample independence assumption could be probed directly: on a long chain, compare THAMES computed from heavily thinned, lightly thinned, and independent draws; the difference in estimates is a measure of how much the proofs rely on the approximation.","The $G!$ discrepancy with bridge sampling in the well-separated simulation suggests a diagnostic for practitioners: when two marginal likelihood estimators disagree by a large factor, test whether one proposal is label-asymmetric, since the paper's symmetrisation is what removes that bias.","For very large numbers of components, the paper's own discussion indicates variance grows as the ellipsoid must be shrunk to keep the ordering set small; a useful extension would be an automated rule for trading off ellipsoid size against variance."],"forward_implications":["Marginal likelihood becomes a routine post-processing step for multivariate Gaussian mixtures with up to at least 15 components, using only relabelled posterior draws and unnormalised log-posterior values.","Because hidden allocation vectors are not required, the estimator works with samplers that do not produce them, widening the class of tractable mixture models.","Model selection for mixtures can be performed with an estimator that is a formal reciprocal importance sampler, so its consistency and asymptotic normality bring the usual machinery of Monte Carlo error assessment.","In well-separated settings the $G!$-factor bias seen in bridge sampling is avoided, so the two estimators can serve as mutual checks in practice.","The criterion of overlap provides a complementary answer to a different question: how many distinguishable clusters exist, alongside the marginal likelihood's answer to how many components the model needs."],"supporting_citations":[{"why":"Supplies the original THAMES definition, the tuning constant c = sqrt(R+1), and the unbiasedness, consistency and asymptotic-normality results that the mixture extension inherits.","marker":"[41]"},{"why":"Introduces the intersection of an ellipsoid with an HPD region and its Monte Carlo volume estimate, which replaces the pure HPD region in the mixture setting.","marker":"[51]"},{"why":"Provides the symmetrisation device of averaging a non-symmetric estimator over label permutations, on which the symmetric THAMES is built.","marker":"[6]"},{"why":"Establishes reciprocal importance sampling as the framework that yields unbiasedness, consistency and asymptotic normality for this estimator.","marker":"[18]"},{"why":"Supplies one of the relabelling algorithms that produce the relabelled posterior sample the estimator requires.","marker":"[60]"},{"why":"Provides the ECR relabelling algorithm used in the simulations and real-data analyses, an alternative to the other relabelling method.","marker":"[49]"},{"why":"Gives the identity linking the marginal likelihood ratio Z(G)/Z(G-1) to the posterior probability that a component is empty, used to handle empty components.","marker":"[46]"},{"why":"Defines bridge sampling, the main baseline against which the THAMES is compared in the simulation studies.","marker":"[38]"},{"why":"Supplies the optimality of quadratic discriminant analysis used to construct the ordering score W.","marker":"[21]"}],"fun_headline_variants":["THAMES estimator computes mixture evidence exactly via graph orderings","Symmetrised harmonic mean gives exact evidence for high-dimensional mixtures","Invariant to label switching, THAMES yields easy marginal likelihood","Mixture model evidence via efficient graph-based permutation sum"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The theoretical guarantees assume the MCMC draws used in the second half are effectively independent of each other and of the first half, which defines the truncation set; the paper states this is approximate, and the proofs do not quantify the effect of dependence or of estimating the truncation set from the data.","fun_headline_variants_meta":{"raw":{"variants":["THAMES estimator computes mixture evidence exactly via graph orderings","Symmetrised harmonic mean gives exact evidence for high-dimensional mixtures","Invariant to label switching, THAMES yields easy marginal likelihood","Mixture model evidence via efficient graph-based permutation sum"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.001183,"raw_usage":{"total_tokens":4888,"prompt_tokens":950,"completion_tokens":3938,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":566,"completion_tokens_details":{"reasoning_tokens":3869}},"tokens_in":566,"tokens_out":3938,"duration_ms":29398,"temperature":1.0,"reasoning_tokens":3869,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-16T04:55:05.022913+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"On a small mixture with an exactly computable marginal likelihood, take one long MCMC run and compute THAMES on the raw chain, on a chain thinned by a factor of 10, and on a set of independent draws; if the three estimates differ by more than their reported Monte Carlo standard errors, the independence and split-sample assumptions are violated in a way that matters.","supporting_citations":[{"cited_title":"Perrot-Dockès, S","cited_arxiv_id":null,"evidence_quote":"Supplies the original THAMES definition, the tuning constant c = sqrt(R+1), and the unbiasedness, consistency and asymptotic-normality results that the mixture extension inherits."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Introduces the intersection of an ellipsoid with an HPD region and its Monte Carlo volume estimate, which replaces the pure HPD region in the mixture setting."},{"cited_title":"van Mechelen, and A","cited_arxiv_id":null,"evidence_quote":"Provides the symmetrisation device of averaging a non-symmetric estimator over label permutations, on which the symmetric THAMES is built."},{"cited_title":"Bayesianmodelchoice:asymptoticsandexact calculations","cited_arxiv_id":null,"evidence_quote":"Establishes reciprocal importance sampling as the framework that yields unbiasedness, consistency and asymptotic normality for this estimator."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Supplies one of the relabelling algorithms that produce the relabelled posterior sample the estimator requires."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Provides the ECR relabelling algorithm used in the simulations and real-data analyses, an alternative to the other relabelling method."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Gives the identity linking the marginal likelihood ratio Z(G)/Z(G-1) to the posterior probability that a component is empty, used to handle empty components."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Defines bridge sampling, the main baseline against which the THAMES is compared in the simulation studies."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Supplies the optimality of quadratic discriminant analysis used to construct the ordering score W."}],"review_version":1}