{"id":"372006fb-7140-401e-831e-3d88417437e2","arxiv_id":"2507.16756","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":4,"one_line_summary":"Bayesian inference for discretely observed CTMCs via a pseudo-likelihood that couples the transition matrix to a spectral generator, with a Gibbs sampler whose cost per iteration is near-invariant in the number of observations.","lead":"A new statistical method estimates how continuous time Markov chains change between regularly spaced snapshot observations, without reconstructing every hidden jump. It runs much faster than existing Bayesian approaches on long data and handles state spaces of tens to hundreds.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The sampler's P-update ignores the spectral penalty in (6), so the implemented algorithm is not shown to target the pseudo-posterior; coherent uncertainty for L is therefore unsupported even within the biorthogonal class.","rationale":"The reader's weakest_assumption focuses on the biorthogonal spectral decomposition excluding defective generators. That is a real limitation, but it is an explicit model restriction acknowledged in the paper. A more damaging issue is that the algorithm as implemented does not appear to sample from the pseudo-posterior that the theory is supposed to justify. Section 3.1 gives the Dirichlet update as an asymptotic approximation and offers an MH correction, but the algorithm and experiments use the approximation directly. The reported simulations use a large penalty coefficient (ν=10^4) and relatively small n, so the approximation's error is not automatically negligible. The theoretical guarantees in Section 4 are for a different object: Theorem 1 is a BvM result for the Dirichlet-multinomial posterior of P, and Theorem 2 is a delta-method for the plug-in spectral estimator. Neither establishes that the pseudo-posterior from (6) concentrates or that its credible intervals have valid frequentist coverage. Because the paper's central contribution is scalable Bayesian inference with coherent uncertainty quantification, the mismatch between the implemented sampler and the claimed target distribution is the most load-bearing concern. I do not think it warrants rejection, since the MH variant is described and could be implemented, but it must be demonstrated before the posterior uncertainty claims can be accepted. This does not change the reader's conditional verdict, hence UNCHANGED.","tokens_in":19793,"tokens_out":4122,"duration_ms":43907,"concrete_test":"Implement the Metropolis–Hastings variant suggested in Section 3.1: propose each row P(p,·) from Dir(α+cp) and accept with probability min(1, L(xn|P_prop,{λk,ϕk,ψk})/L(xn|P_current,{λk,ϕk,ψk})) using the pseudo-likelihood (6), with the same hyperparameters as the §5.1 simulations (ν=10^4, n=10^3, m=4). Report the average acceptance rate and compare posterior means and 95% credible intervals for L to those obtained from the direct-Dirichlet sampler in Algorithm 1. If acceptance rates are materially below 1 or the intervals differ beyond Monte Carlo error, the approximation in (7) is not adequate in this regime and the implemented sampler is not targeting the pseudo-posterior.","verdict_should_be":"UNCHANGED","load_bearing_attack":"Section 3.1 observes that the full conditional of P under the pseudo-likelihood (6) is non-conjugate, but states that as n→∞ or ν→0 it converges to the Dirichlet posterior Dir(α+cp). Algorithm 1 (line 5) samples each row of P directly from this limiting Dirichlet distribution, with the comment 'Or use as proposal in MH if desired.' No MH acceptance step is implemented or reported. Consequently, the deployed sampler does not target the joint pseudo-posterior defined by (6): P is drawn from its Dirichlet-multinomial marginal while the spectral parameters are updated conditionally on P through the spectral penalty. The stationary distribution of this composite update is not specified. The paper's theoretical results do not repair this gap: Theorem 1 is a Bernstein–von Mises theorem for the Dirichlet-multinomial posterior P(P|xn), not for the pseudo-posterior of (6); Theorem 2 is a frequentist delta-method for the plug-in spectral estimator, not a posterior consistency result for the spectral parameters under the pseudo-posterior. Thus the central claim of 'coherent posterior uncertainty' for the generator L is not justified by the stated theory and is potentially violated by the algorithm as implemented. This concern is load-bearing because it affects the method's validity in exactly the finite-sample regime where the paper claims scalable Bayesian inference: the simulations use ν=10^4 and n as small as 10^2, so the spectral penalty in (6) is not negligible, and the approximation underlying the Dirichlet update can be poor.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper proposes a scalable Bayesian approach for inference on discretely observed continuous-time Markov chains. It defines a pseudo-likelihood that combines the observed transition log-likelihood with a Frobenius penalty coupling the transition matrix P to a biorthogonal spectral decomposition of the generator L (Eq. 6). A Gibbs sampler (Algorithm 1) alternates between sampling rows of P from Dirichlet distributions based on transition counts and sampling the spectral parameters from truncated Gaussian full conditionals. The paper claims a Bernstein-von Mises theorem for P and posterior consistency for the spectral parameters of L, and supports the method with simulations and a metastable diffusion application. The central promise is near-cost-invariant inference with coherent posterior uncertainty, bypassing latent path integration.","tokens_in":20092,"tokens_out":5620,"duration_ms":64228,"significance":"If the method worked as claimed, it would be a useful contribution: it avoids latent path sampling, respects embeddability within the assumed spectral class, and appears computationally attractive for state spaces of tens to hundreds of states. The paper also provides code and an R package, and it compares against exact latent-path methods. However, the gap between the theoretical results and the actual sampler is substantial, and the current manuscript does not establish the advertised coherent posterior uncertainty. The biorthogonal assumption and the restriction to real spectra are also underexplored. The core algorithmic idea is promising, but the theoretical framing and implementation claims need significant reworking.","major_comments":[{"comment":"The implemented sampler does not target the pseudo-posterior defined by Eq. (6). The full conditional of P under (6) is non-conjugate; Section 3.1 argues that it converges to Dir(α + c_p) as n → ∞ or ν → 0, and Algorithm 1 line 5 samples P directly from this limiting Dirichlet, with only a parenthetical suggestion to use it as an MH proposal. No Metropolis-Hastings acceptance step is reported, so the stationary distribution of the composite update is not the pseudo-posterior (6). The simulations use ν = 10^4 and n as small as 10^2, so the penalty term is far from negligible in the regime where the approximation is used. This undermines the claim that the sampler yields coherent posterior uncertainty for L in the finite-sample settings presented.","section":"Section 3.1 and Algorithm 1, line 5"},{"comment":"Theorem 2 is a frequentist delta-method for the spectral decomposition of the plug-in empirical generator L̂ = Δ^{-1} log P̂, not a posterior consistency result for the spectral parameters under the pseudo-posterior (6) or under Algorithm 1's sampling distribution. The theorem's statement concerns √n(λ̂_k − λ_k^0), √n(φ̂_k − φ_k^0), and √n(ψ̂_k − ψ_k^0) for the empirical estimator, not the posterior samples produced by the sampler. The abstract's claim of 'posterior consistency for the spectral parameters' is therefore not supported by the displayed result. Additionally, the proof's formula for the left eigenvector, √n(ψ̂_k − ψ_k^0) → ψ_k^0 Q R_k, does not match the theorem statement's √n(ψ̂_k − ψ_k^0) → R_k^T Q^T ψ_k^0; these differ and the dimensions should be checked.","section":"Section 4.2, Theorem 2"},{"comment":"The method as implemented is restricted to generators with real eigenvalues. The prior in Eq. (8) constrains the eigenvalues to 1 = Λ_1 > Λ_2 ≥ ... ≥ Λ_m > 0, and the eigenvectors are modeled as real Gaussian vectors in Eq. (9). However, the biorthogonal representation L = Σ λ_k φ_k ψ_k^T in Section 2.2 allows, in general, complex eigenvalues and eigenvectors for non-reversible CTMCs, and Section 2.2 even orders eigenvalues by Re(λ_k). For a non-reversible generator with a complex conjugate pair, the proposed sampler cannot run. This scope restriction is not stated as a modeling assumption in the abstract or introduction, and the simulation studies only use symmetric generators, so the method's applicability to general CTMCs is overstated.","section":"Section 3.2 and Eq. (8)"},{"comment":"The Bernstein-von Mises theorem in Theorem 1 concerns the posterior P(P | x_n) under the Dirichlet-multinomial likelihood, not the pseudo-posterior of Eq. (6). Even if one accepts the algorithm's replacement of the P-conditional with Dir(α + c_p), Theorem 1 does not account for the spectral penalty in (6); hence it cannot justify the two-stage procedure used in Algorithm 1. What is needed is a statement about the joint pseudo-posterior, or an explicit analysis of the approximate sampler's stationary distribution, neither of which is provided.","section":"Section 4.1 and Theorem 1"}],"minor_comments":[{"comment":"In the m = 4, n = 10^4 row for the 'blandt' method, the reported Frobenius norm is 0.16 with a standard deviation of 0.79; a standard deviation larger than the mean by a factor of five seems implausible and is likely a typographical error.","section":"Table 2"},{"comment":"The covariance matrix of Q is stated as 2 I^{-1}(P_0) without derivation; the factor 2 and the precise form of the Fisher information for a Markov chain with row-wise multinomial counts should be made explicit, since the asymptotic covariance of the empirical transition matrix is not simply the inverse Fisher information of a single transition.","section":"Theorem 2 statement and proof"},{"comment":"The notation uses m for both the state-space size and the number of spectral components; this is acceptable only if it is clearly stated that the spectral truncation uses all m components. A short clarification would help.","section":"Section 2.2"},{"comment":"The phrase 'converges in total variation of the conditional posterior' is imprecise: the convergence in Eq. (7) should be stated as convergence of the conditional posterior distribution of P(p, ·) to the specified Dirichlet distribution, and under which norm (total variation) it holds.","section":"Section 3.1"}],"recommendation":"major_revision","confidential_remarks":"The reader's stress-test concern is well-founded and lands directly on the central claim. The paper needs either to implement a valid Metropolis-Hastings correction so that the sampler targets Eq. (6), or to reframe the method as an approximate pseudo-posterior algorithm with explicit approximation error control in the finite-sample regime used in the simulations. The theoretical claims in Theorem 2 also need to be relabeled as frequentist delta-method results, and the real-spectrum restriction must be stated prominently. I recommend major revision rather than rejection because the core algorithmic idea is promising and the empirical comparisons show competitive accuracy, but the current gap between the claims and the analysis is too large to accept."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"The core idea here is solid and practically valuable: a Gibbs sampler for discretely observed CTMCs that avoids latent path integration by penalizing the transition matrix toward a biorthogonal spectral form. The empirical comparisons are the strength. bigmac matches Bladt–Sørensen accuracy in the tested regimes (m up to 8, n up to 1e5) at a fraction of the cost, and the m=30 metastable diffusion example suggests scalability. The code is public, which is a plus.\n\nWhat is genuinely new is the pseudo-likelihood in (6)—transition log-likelihood plus a Frobenius penalty to the spectral decomposition—and the resulting sampler with state-dependent truncations. The biorthogonal parametrization itself comes from Crommelin and Vanden-Eijnden (2006), and the paper is honest about that. The near-O(1)-in-n cost is the key selling point, and the simulations support it.\n\nBut the stress-test concern is real and load-bearing. Section 3.1 states that the full conditional of P under (6) is not conjugate, and that sampling P from the Dirichlet-multinomial is an approximation valid as n→∞ or ν→0. Algorithm 1 then does exactly that, with no MH correction, while the spectral parameters are updated using the penalty. The stationary distribution of this composite update is not specified. The theory does not rescue the gap: Theorem 1 is a BvM for the Dirichlet-multinomial posterior for P, not for the pseudo-posterior of (6); Theorem 2 is a frequentist delta-method for the plug-in log-P estimator, not posterior consistency for the spectral parameters. The abstract and Section 4.2 overstate these results. In the simulations, ν=10^4 and n as low as 1e2, so the approximation is not obviously benign.\n\nThe biorthogonal assumption is another soft spot—it excludes defective generators, and there are no diagnostics for when it fails. The paper's own truncation-based prior is acknowledged as approximate, which is fine.\n\nWho should read this: practitioners who need fast approximate Bayesian inference for CTMCs, and methodologists interested in approximate MCMC with theory gaps. A serious referee should push the authors to either implement the MH step, provide a bound on the distance between the approximate and exact stationary distributions, or soften the uncertainty quantification claims. The empirical contribution deserves publication after major revision, not desk rejection.","headline":"Genuinely useful computational idea for CTMC inference, but the theoretical claims overstate what is proved and the implemented sampler targets an approximation of the pseudo-posterior; still worth serious refereeing.","tokens_in":20685,"tokens_out":2966,"would_cite":true,"duration_ms":33925,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["62M05","62F15","60J27"],"pacs":[],"model":"deepseek-v4-flash","headline":"A pseudo-likelihood coupling the observed transition matrix to a biorthogonal spectral form of the generator makes Bayesian inference for discretely observed continuous-time Markov chains both scalable and asymptotically justified.","keywords":["continuous-time Markov chain","generator matrix","pseudo-likelihood","Bayesian inference","Bernstein-von Mises theorem","biorthogonal spectral decomposition","embeddability","Gibbs sampling"],"falsifier":"Simulate a three-state CTMC whose generator is defective (a single Jordan block, so no biorthogonal decomposition exists), observe it at regular intervals, and run the proposed Gibbs sampler: if the posterior contracts to a generator that misrepresents the true transition probabilities or the spectral penalty residual remains large, the biorthogonal assumption is doing the work.","tokens_in":19562,"feed_emoji":"⏱️","tokens_out":6112,"duration_ms":59956,"temperature":0.7,"pith_summary":"The paper proposes a Bayesian method for continuous-time Markov chains observed only at discrete time points, where the exact likelihood is intractable. Instead of integrating over unobserved holding times, it couples the transition matrix P with the generator L through a pseudo-likelihood that penalizes differences between P and the spectral form of L. The resulting Gibbs sampler has cost nearly invariant to the number of observations and scales to state spaces with tens to hundreds of states. The authors prove a Bernstein–von Mises theorem for P and posterior consistency for the spectral parameters of L, and demonstrate the method on simulations and a metastable diffusion.","feed_headline":"A shortcut makes Bayesian CTMC inference scale to dozens of states","feed_subtitle":"Pseudo-likelihood ties the transition matrix to the generator's spectral form, skipping expensive latent path sampling.","key_machinery":"The load-bearing object is the biorthogonal spectral decomposition $L = \\sum_{k=1}^m \\lambda_k \\phi_k \\psi_k^\\top$ with $\\psi_j^\\top \\phi_k = \\delta_{jk}$, under which $\\exp(\\Delta L) = \\sum_k \\exp(\\lambda_k \\Delta) \\phi_k \\psi_k^\\top$. The pseudo-likelihood combines the observed-transition log-likelihood with a Frobenius penalty $\\|P - \\sum_k \\exp(\\lambda_k \\Delta) \\phi_k \\psi_k^\\top\\|_F^2$ scaled by $\\nu > 0$. This couples the empirical transition matrix to the generator, and the conditional updates for P, the eigenvalues, and the eigenvectors take conjugate-friendly forms: Dirichlet for the rows of P, and truncated normal or multivariate normal for the spectral parameters. Biorthogonality is encouraged by a shrinkage prior toward the Stiefel manifold and truncations that keep off-diagonal entries of L nonnegative.","core_discovery":"The central claim is that the intractable likelihood of a discretely observed CTMC can be replaced by a pseudo-likelihood that is both tractable and asymptotically justified: the log-likelihood of the observed transitions plus a Frobenius penalty tying P to the biorthogonal spectral representation $L = \\sum_{k=1}^m \\lambda_k \\phi_k \\psi_k^\\top$. Under this pseudo-likelihood, the posterior of P is asymptotically normal at the parametric rate, and the spectral parameters of L converge at $\\sqrt{n}$, with the rate explicitly showing that closely spaced eigenvalues inflate the uncertainty in eigenvectors. The payoff is a Gibbs sampler that avoids latent path integration, obeys embeddability, and returns an embeddable generator with coherent posterior uncertainty.","pith_inferences":["Because the asymptotic theorems fix the regularization strength $\\nu$, the finite-sample bias-variance tradeoff is not addressed; one testable extension is to let $\\nu$ grow slowly with n and examine the posterior coverage of P.","The biorthogonal assumption rules out defective generators; for near-defective systems, the posterior may concentrate on a diagonalizable approximation, and the Frobenius residual could serve as a diagnostic for that mismatch.","The same pseudo-likelihood coupling could be extended to irregularly spaced observations by modeling P as a function of the time lag, rather than fixing a single $\\Delta$.","When eigenvalues are close, the inflated eigenvector uncertainty provides a principled way to flag non-identifiability in the embedding problem."],"forward_implications":["Bayesian inference for discretely observed CTMCs no longer requires sampling latent paths, so computation is nearly invariant to the number of observations.","The posterior distribution of the transition matrix P concentrates at the parametric rate, giving calibrated uncertainty quantification.","The spectral parameters of the generator converge at $\\sqrt{n}$, and the convergence rate reveals that closely spaced eigenvalues make eigenvectors harder to estimate.","The method returns an embeddable generator, avoiding the failure mode where an empirical transition matrix has no valid matrix-logarithm.","The approach scales to medium-to-high dimensional state spaces, demonstrated on a 30-state metastable diffusion example."],"supporting_citations":[{"why":"Defines the embedding problem that motivates the biorthogonal parametrization of the generator.","marker":"Kingman (1962)"},{"why":"Introduces the biorthogonal spectral parametrization and the low-rank caveat that this paper builds on and moves beyond.","marker":"Crommelin & Vanden-Eijnden (2006)"},{"why":"Provides the exact-likelihood latent-path integration baseline whose computational cost the proposed method aims to avoid.","marker":"Bladt & Sørensen (2005)"},{"why":"Supplies the latent-path Gibbs sampling speed-up used in the baseline implementation that is compared against.","marker":"Fearnhead & Sherlock (2006)"},{"why":"Offers the Metropolis–Hastings alternative without latent paths, serving as a second baseline requiring proposal tuning.","marker":"Riva-Palacio et al. (2023)"},{"why":"Establishes the weakly dependent Bernstein–von Mises conditions used to prove Theorem 1.","marker":"Connault (2014)"},{"why":"Provides the concentration inequalities used to verify Conditions (C6) and (C7).","marker":"Paulin (2015)"},{"why":"Supplies the first-order matrix perturbation theory used to propagate uncertainty from P to spectral parameters in Theorem 2.","marker":"Stewart & Sun (1990)"}],"fun_headline_variants":["Pseudo-likelihood unlocks scalable Bayesian CTMC inference","Skip latent paths: efficient Bayesian CTMC inference","Bayesian CTMC inference without latent path sampling","Pseudo-likelihood makes CTMC inference scalable and fast"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The data-generating generator must be diagonalizable with a biorthogonal spectral decomposition, and the observations must be regularly spaced; if the true generator is defective or the observation times are irregular, the stated guarantees do not follow.","fun_headline_variants_meta":{"raw":{"variants":["Pseudo-likelihood unlocks scalable Bayesian CTMC inference","Skip latent paths: efficient Bayesian CTMC inference","Bayesian CTMC inference without latent path sampling","Pseudo-likelihood makes CTMC inference scalable and fast"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000666,"raw_usage":{"total_tokens":3014,"prompt_tokens":892,"completion_tokens":2122,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":508,"completion_tokens_details":{"reasoning_tokens":2058}},"tokens_in":508,"tokens_out":2122,"duration_ms":14972,"temperature":1.0,"reasoning_tokens":2058,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-06T15:02:42.791945+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Simulate a three-state CTMC whose generator is defective (a single Jordan block, so no biorthogonal decomposition exists), observe it at regular intervals, and run the proposed Gibbs sampler: if the posterior contracts to a generator that misrepresents the true transition probabilities or the spectral penalty residual remains large, the biorthogonal assumption is doing the work.","supporting_citations":[{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Defines the embedding problem that motivates the biorthogonal parametrization of the generator."},{"cited_title":"& Vanden-Eijnden, E","cited_arxiv_id":null,"evidence_quote":"Introduces the biorthogonal spectral parametrization and the low-rank caveat that this paper builds on and moves beyond."},{"cited_title":"& S rensen, M","cited_arxiv_id":null,"evidence_quote":"Provides the exact-likelihood latent-path integration baseline whose computational cost the proposed method aims to avoid."},{"cited_title":"& Sherlock, C","cited_arxiv_id":null,"evidence_quote":"Supplies the latent-path Gibbs sampling speed-up used in the baseline implementation that is compared against."},{"cited_title":", Mena, R","cited_arxiv_id":null,"evidence_quote":"Offers the Metropolis–Hastings alternative without latent paths, serving as a second baseline requiring proposal tuning."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Establishes the weakly dependent Bernstein–von Mises conditions used to prove Theorem 1."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Provides the concentration inequalities used to verify Conditions (C6) and (C7)."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Supplies the first-order matrix perturbation theory used to propagate uncertainty from P to spectral parameters in Theorem 2."}],"review_version":1}