Estimation and uncertainty of reversible Markov models

Benjamin Trendelkamp-Schroer; Hao Wu; Fabian Paul; Frank Noé

Source-linked AI summary

Estimation and uncertainty of reversible Markov models

Benjamin Trendelkamp-Schroer, Hao Wu, Fabian Paul, Frank Noé

arXiv:1507.05990v2Jul 19, 2015physics.chem-phmath-phphysics.comp-phphysics.data-anstat.CO

TL;DR

The paper addresses how to estimate reversible Markov models and quantify uncertainty from finite, correlated molecular-dynamics simulation data. It develops maximum-likelihood estimators, sparse-prior Bayesian inference, and efficient Gibbs samplers, including models with fixed stationary distributions. The resulting priors produce error bars that envelop reference estimates, while the new samplers substantially reduce autocorrelation times relative to a previous Gibbs method.

Problem
Finite molecular-dynamics trajectories provide correlated transition counts and rare events, leaving uncertainty quantification, suitable priors, and efficient reversible-model sampling unresolved.
Method
The paper develops reversible maximum-likelihood estimators, sparse priors for Bayesian inference, and efficient Gibbs or Metropolis-within-Gibbs samplers, including fixed-stationary-distribution models.
Results
Sparse priors produce error bars that envelop reference estimates, while the new sampler's autocorrelation times are more than 5-fold shorter for 233 states and more than 13-fold shorter for 1108 states.
Takeaways & Limitations
Reversible estimation can provide statistically efficient uncertainty-aware Markov models while preserving the observed process's metastable sparsity structure.

Abstract

from arXiv · show

Reversibility is a key concept in Markov models and Master-equation models of molecular kinetics. The analysis and interpretation of the transition matrix encoding the kinetic properties of the model relies heavily on the reversibility property. The estimation of a reversible transition matrix from simulation data is therefore crucial to the successful application of the previously developed theory. In this work we discuss methods for the maximum likelihood estimation of transition matrices from finite simulation data and present a new algorithm for the estimation if reversibility with respect to a given stationary vector is desired. We also develop new methods for the Bayesian posterior inference of reversible transition matrices with and without given stationary vector taking into account the need for a suitable prior distribution preserving the meta- stable features of the observed process during posterior inference. All algorithms here are implemented in the PyEMMA software - http://pyemma.org - as of version 2.0.

I. INTRODUCTION

Reversible Markov models encode physically motivated detailed balance, reducing statistical complexity and ensuring real, interpretable kinetics. The paper develops maximum-likelihood and Bayesian methods to estimate such models from finite, correlated simulation data, including cases with a fixed stationary distribution.

I. INTRODUCTION: Transition probabilities are estimated from trajectory counts, but rare events and correlated transitions make uncertainty quantification and effective counting difficult.Determining statistical inefficiencies for transition-count matrices remains an open problem, although the estimation algorithms can use future effective-count methods.
I. INTRODUCTION: The paper develops simplified maximum-likelihood estimators and Bayesian posterior algorithms for reversible transition matrices, including models with a fixed equilibrium distribution.It also addresses priors and sampling procedures designed for finite simulation data.
A. From microscopic reversibility to discrete-state detailed balance: Finite-data estimates need not satisfy detailed balance exactly, so the paper enforces reversibility during estimation to reflect microscopic reversibility of equilibrium molecular dynamics.The discrete-state relation is π_i p_ij(τ) = π_j p_ji(τ), and it holds independently of lag time τ.
I. INTRODUCTION: Detailed balance reduces the number of independent transition-matrix parameters by roughly half, lowering statistical error and enabling physically meaningful equilibrium-kinetics analyses.Reversible matrices also support analysis tools and matrix routines unavailable or less suitable for nonreversible models.
B. Eigenvalues and eigenvectors: Reversible estimates produce real eigenvalues and can be more accurate and less uncertain than nonreversible estimates in the reported three-state example.For λ3, the reversible estimate is 0.18±0.04 versus 0.07±0.21 for the nonreversible estimate, while nonreversible estimates can have complex eigenvalues.

C. Equilibrium kinetics analyses

The paper formulates likelihood-based Markov-model estimation from transition counts and connects reversible transition matrices to equilibrium-kinetics analyses. Reversibility supports relaxation, metastability, and transition-path analyses, while the paper presents estimation algorithms for known or unknown equilibrium distributions.

Equilibrium kinetics analyses: Reversible transition matrices provide dominant relaxation rates and timescales relevant to kinetic experiments.The largest nontrivial eigenvalues determine relaxation rates, whose inverses are relaxation timescales.
Equilibrium kinetics analyses: Reversibility enables spectral decompositions for molecular observables and supports PCCA+ analysis of metastable sets.PCCA+ is reliably applicable when the transition matrix has real-valued eigenvalues and eigenvectors.
Equilibrium kinetics analyses: In reversible models, transition-path theory yields complementary forward and backward committors and corresponding net-flux expressions.Transition-path statistics can also be computed for nonreversible matrices, but the committor relationship is specific to the reversible case.
Likelihood and estimation: The paper presents likelihood formulations and new maximum-likelihood estimators for reversible Markov models with known or unknown equilibrium distributions.The framework begins from trajectory likelihoods and their count-matrix representation.
Likelihood and counting: The count matrix is sufficient for the Markov-model likelihood, so transition sequences can be reduced to transition counts for estimation.Multiple trajectories contribute by adding their count matrices.

B. Counting

The section explains how transition counts support likelihood estimation and how counting choices affect statistical efficiency and bias. It then develops numerical estimation of reversible transition matrices from generally nonsymmetric count matrices.

Counting: Counting all available lag-τ transitions in a sliding window improves data usage but creates statistically dependent counts when the process is not Markovian at lag 1.Effective counts can correct for dependence through lag- and transition-specific statistical inefficiencies.
Counting: The determination of transition-count statistical inefficiencies remains an open problem, although the estimation algorithms can use future effective-count methods unchanged.A cited approach provides consistent but somewhat too-small uncertainty intervals for practical molecular-dynamics data.
Counting: Symmetrizing forward and backward counts is discouraged because it creates large bias unless trajectories are much longer than the molecule’s slowest timescales.The reversible estimator instead uses the unmodified count matrix while imposing detailed balance during estimation.
Reversible estimation: The reversible MLE enforces detailed balance on the original, generally nonsymmetric forward-only count matrix rather than symmetrizing the data.The optimization has no closed-form solution and is solved numerically, including through a fixed-point iteration.
Numerical algorithms: The quadratic optimizer can converge faster than the direct fixed-point iteration by solving each coordinate’s quadratic subproblem exactly while holding others fixed.Both approaches estimate the reversible matrix from count statistics.
Reversible estimation: The reversible maximum-likelihood estimator is unique because the equivalent log-variable optimization is convex.The fixed-point iteration preserves symmetry and inequality constraints when initialized with a feasible iterate.

3. Reversible estimation for given stationary vector

The paper develops a reversible maximum-likelihood estimator constrained to a specified stationary distribution. Its fixed-point iteration avoids the O(n^6) cost of a standard interior-point solution and converges nearly three orders of magnitude faster than the unknown-equilibrium case in the reported example.

Motivation and formulation: The fixed-stationary-vector algorithm estimates a reversible transition matrix from transition counts while incorporating an externally supplied equilibrium distribution.This targets settings where enhanced sampling informs equilibrium weights more efficiently than unbiased simulations.
Algorithm: O(n^6) interior-point computation becomes infeasible for increasing state-space size, motivating a fixed-point iteration for large systems.The convex formulation has O(n^2) unknowns and O(n^2) equality and inequality constraints.
Estimator properties: The prescribed stationary vector is unique only when the transition matrix is irreducible, so the state space is restricted to the largest connected component of the count graph.This restriction is based on the undirected graph defined by C + C^T.
Algorithm: The fixed-point scheme iterates Lagrange multipliers from positive starting values until convergence, then reconstructs the transition matrix.The inequality constraints are automatically fulfilled when all counts are nonnegative, with special handling required for zero diagonal-count cases.
Convergence: Nearly three orders of magnitude faster convergence is reported for the fixed-stationary-vector estimation than for estimation with an unknown equilibrium distribution.The comparison is described for alanine-dipeptide simulation data and is illustrated alongside convergence results.
Estimator properties: The constrained estimator is unique because it minimizes a strictly convex function over a convex feasible set.The uniqueness argument applies on the support defined by positive symmetrized counts.

IV. BAYESIAN ESTIMATION

Bayesian inference quantifies uncertainty in Markov-model observables, but reliable posterior estimates require priors and count handling that respect sparse metastable dynamics. The section demonstrates that sparse priors avoid artificial pathways and produce credible intervals closer to the true kinetics.

A. Bayes theorem and Monte Carlo sampling: Posterior sampling estimates uncertainty in observables such as relaxation timescales and mean first passage times using Monte Carlo samples of transition matrices.Posterior moments and credible intervals summarize uncertainty from finite observations.
A. Bayes theorem and Monte Carlo sampling: The paper addresses prior choice, correlated transition counts, and sampler efficiency as central obstacles to reliable reversible posterior inference.The proposed approach uses effective uncorrelated counts and new reversible sampling algorithms.
B. Non-reversible sampling: The uniform prior yields a 90% credible interval of [1.9, 2.0] · 10^3 for a true mean first passage time of 2 · 10^5.Artificial nonzero probabilities create pathways around the bottleneck, underestimating the passage time by two orders of magnitude.
B. Non-reversible sampling: The sparse prior b_ij = −1 yields a 90% credible interval of [1.5, 2.7] · 10^5 that covers the true value 2 · 10^5.Its posterior samples preserve the count matrix’s sparsity structure.
B. Non-reversible sampling: The sparse prior’s credible interval reaches the true value with little data, whereas the uniform prior requires simulations orders of magnitude longer than the slowest process.The latter makes inference under the uniform prior impractical in this example.

C. A prior for reversible Markov models

The reversible sampler replaces the Dirichlet prior with a prior defined on symmetric unconditional transition probabilities, enabling detailed-balance-preserving MCMC. With b_ij = −1, unobserved bidirectional transitions are excluded and sampling efficiency is improved through adapted proposals.

C. A prior for reversible Markov models: The new reversible sampler defines the prior on symmetric X matrices, so detailed balance is enforced by sampling symmetric independent variables.Transition matrices are recovered from X after normalization.
C. A prior for reversible Markov models: The sparse prior b_ij = −1 assigns zero posterior probability to transitions with no observed forward or backward counts.This preserves the observed sparsity pattern and prevents unsupported transitions.
C. A prior for reversible Markov models: The algorithm uses Metropolis-within-Gibbs updates with adapted proposal densities, including Gamma approximations for off-diagonal conditional distributions.The proposals match the conditional mode and curvature to obtain acceptance probabilities close to one in typical cases.
C. A prior for reversible Markov models: Transition elements with zero forward and backward counts can be omitted from sampling to account for the sparse prior efficiently.The resulting chain samples only the relevant entries.

E. A prior for reversible Markov models with fixed equilibrium distribution

For a fixed stationary distribution, the reversible sampler updates X while preserving row normalization and symmetry. Its prior again removes unobserved transitions, while diagonal regularization prevents degenerate updates that could break ergodicity.

E. A prior for reversible Markov models with fixed equilibrium distribution: The fixed-stationary-distribution updates preserve π through row normalization and restore symmetry to maintain reversibility.The construction samples independent lower-triangle elements while enforcing the constraints.
E. A prior for reversible Markov models with fixed equilibrium distribution: For unobserved bidirectional transitions, b_kl = −1 forces both transition probabilities to zero in every posterior sample.This concentrates posterior mass around the maximum-likelihood estimate and prevents artificial kinetic pathways.
E. A prior for reversible Markov models with fixed equilibrium distribution: The diagonal prior is selected from the maximum-likelihood estimate so the posterior expectation of p_kk is zero exactly when the estimate is zero.
E. A prior for reversible Markov models with fixed equilibrium distribution: When a diagonal conditional degenerates at zero, choosing b_kk = −1 + ϵ and starting with x_kk > 0 prevents the sampling chain from losing ergodicity.The regularization is used when the corresponding maximum-likelihood diagonal entry is zero.

F. Sampling reversible Markov models with fixed equilibrium distribution

The fixed-equilibrium sampler uses a Metropolis-within-Gibbs procedure with transformed variables and adapted Gamma proposals. These updates target constrained conditional distributions while adding a log-scale move to avoid sticking at extreme values.

F. Sampling reversible Markov models with fixed equilibrium distribution: Algorithm 2 samples each reversible transition pair by transforming x_kl into v and proposing v with an adapted Gamma distribution.The transformation supports bounded transition updates while the proposal approximates the conditional density.
F. Sampling reversible Markov models with fixed equilibrium distribution: A log-normal-style proposal is additionally used to prevent the sampler from becoming stuck at extremely small or large v values.
F. Sampling reversible Markov models with fixed equilibrium distribution: The resulting method is characterized as a Metropolis-within-Gibbs MCMC algorithm with adapted proposal probabilities.
F. Sampling reversible Markov models with fixed equilibrium distribution: The proposal uses a Metropolis acceptance step after mapping the sampled v back to x_kl.The update is applied over index pairs with observed forward or backward counts.

V. RESULTS

The proposed reversible posterior samplers are validated against analytical posteriors and molecular-dynamics reference observables. Sampled distributions agree with analytical results, while alanine-dipeptide posterior means agree with reference values.

Validation: Sampled reversible posteriors agree with analytical probabilities for both unconstrained and fixed-stationary-vector 2 × 2 count matrices.The fixed stationary vector restricts the posterior to the line π1p12 = π2p21.
Alanine dipeptide: The experiments use 10 µs of alanine-dipeptide simulation data discretized on a 20 × 20 dihedral-angle grid.The simulation used OpenMM, Langevin dynamics at 300 K, and a 2 fs time-step.
Alanine dipeptide: Posterior means for alanine-dipeptide implied time-scales agree well with reference values.Tables II and III report posterior means, standard deviations, estimated mean errors, and autocorrelation times.
Alanine dipeptide: Posterior sampling estimates expected hitting times between metastable sets with means in good agreement with reference values.The reported transitions include C5 → Cax, C5 → αL, and C5 → αR.

2. Alanine dipeptide, reversible sampling with fixed equilibrium distribution

Fixing the stationary distribution preserves agreement between posterior means and reference observables while reducing uncertainty. For BPTI, the reversible sampler captures maximum-likelihood relaxation times within its 95% credible intervals and mixes substantially faster than a previous algorithm.

2. Alanine dipeptide, reversible sampling with fixed equilibrium distribution: Fixing the stationary distribution yields posterior means agreeing with reference values while producing smaller standard deviations for all shown observables.The constrained posterior results are reported for implied time-scales and expected hitting times.
3. Bovine pancreatic trypsin inhibitor, reversible sampling: For BPTI, maximum-likelihood relaxation times lie within the reversible sampler’s 95% credible intervals across the evaluated lag times.The analysis used effective count matrices over lag times up to 2 µs.
3. Bovine pancreatic trypsin inhibitor, reversible sampling: The sampler has acceptance probability p = 1.0 for diagonal-element updates, supporting efficient sampling of large reversible matrices.The comparison attributes poor mixing in the previous method to the update that modified the equilibrium distribution.
3. Bovine pancreatic trypsin inhibitor, reversible sampling: The new sampler’s autocorrelation times are more than 5-fold shorter for 233 states and more than 13-fold shorter for 1108 states than the previous sampler.The autocorrelation time increases only mildly with matrix dimension.
VI. CONCLUSION: The paper presents a fixed-equilibrium-distribution maximum-likelihood estimator, sparse priors, and efficient Gibbs samplers for reversible transition matrices.The sparse prior preserves the observed transition-matrix sparsity structure, while the treatment of correlated transition counts remains preliminary.

Appendix A: Details for transition matrix sampling

The transition-matrix sampler updates a selected matrix element using a scale-invariant proposal, renormalizes the resulting matrix, and accepts the update with a posterior-based probability.

Appendix A: Details for transition matrix sampling: A selected transition-matrix element is proposed, while the remaining elements are transformed consistently with the matrix constraints.The proposal-density derivation uses the free variables of the symmetric reversible representation.
Appendix A: Details for transition matrix sampling: The proposal density is scale-invariant, and the proposed matrix is renormalized to retain an element sum of 1.Scale invariance relates proposal densities for matrices differing by a positive multiplicative constant.
Appendix A: Details for transition matrix sampling: The posterior remains invariant when the proposed matrix is accepted with probability min{1, pacc}.The acceptance rule follows from the invariant-distribution argument for the proposed update step.

2. Reversible transition matrix sampling: Efficient proposal densities

The efficient proposal construction uses distribution-specific updates for diagonal and off-diagonal elements, with auxiliary variables sampled from Beta or Gamma distributions.

2. Reversible transition matrix sampling: Efficient proposal densities: The transformed diagonal update can be sampled as s′ ∼ Beta(ckk, ck − ckk) and then back-transformed to obtain x′kk.This construction follows from the derived form of the conditional proposal distribution.
2. Reversible transition matrix sampling: Efficient proposal densities: For off-diagonal elements, the proposal density with bkl = −1 is approximated by a three-parameter family leading to an efficiently sampled Gamma distribution.The parameters α, β, and f0 are matched through second derivatives at the maximum point.
2. Reversible transition matrix sampling: Efficient proposal densities: The maximum point for the proposal approximation is obtained from a quadratic equation, with the negative-root solution excluded because the transformed variable must be non-negative.The quadratic coefficients are denoted a, b, and c.

3. Reversible sampling with fixed stationary distribution: efficient proposal densities

The analysis establishes a maximum-likelihood estimate with at most one positive diagonal entry unsupported by observed self-transitions, then approximates a related density using a Gamma family.

Iteratively transforming a candidate estimate increases its likelihood while preserving maximum-likelihood status, yielding the stated diagonal-entry property.
A maximum-likelihood estimate exists with at most one state having a positive diagonal entry despite zero observed self-transitions.
The density in (A25), for s > 1, is approximated by a Gamma distribution using a three-parameter family.
The resulting approximate density has the same stated desirable properties as the approximation in (A22).