Skip to content

Stochastic path sampling

sample_paths runs Monte Carlo walkers on the MSM until two convergence criteria for each endpoint are simultaneously satisfied:

  1. ESS via Sokal-windowed integrated autocorrelation. For each target b, the indicator 1[walker hit b] is treated as a time series across walker IDs; the integrated autocorrelation time tau_int is estimated with the standard Sokal windowing. The effective sample size is N / (2 * tau_int). Default threshold: ess_min=200.
  2. Gelman-Rubin R-hat across chains. Walkers are split into m >= 4 chains; R-hat is computed on the binary hit indicator. Default threshold: rhat_max=1.01.

Endpoints with empirical hit probability below relevance_threshold (default 1e-3) are excluded from the check, so vanishingly-rare targets do not hold up convergence.

Stochastic walker trajectories traced across the genotype graph from wildtype Stochastic walker trajectories traced across the genotype graph from wildtype

Each walker is one Monte Carlo realization of the evolutionary chain; overlaying many traces shows which routes the ensemble actually takes. Sampling continues until both convergence diagnostics settle:

Sampler convergence: effective sample size rising and R-hat falling toward its threshold Sampler convergence: effective sample size rising and R-hat falling toward its threshold

v1 vs v2

v1 stopped on a Euclidean distance threshold over the empirical path distribution. That criterion has no statistical interpretation: correlated samples push the L2 norm down faster than the marginals actually converge, the distribution is multinomial on a sparse support, and continuing to sample after the marginals stabilize keeps reducing the L2 norm anyway.

v2 replaces it with ESS + R-hat. KL divergence with bootstrap CI is offered as a separate diagnostic in analysis.observables but is not a stopping rule (KL is unstable on sparse empirical multinomials with zeros).

API

from gpvolve import ConvergenceCheck, sample_paths

cc = ConvergenceCheck(
    ess_min=200.0,
    rhat_max=1.01,
    relevance_threshold=1e-3,
    chunk_size=10_000,
    max_walkers=1_000_000,
    n_chains=4,
)
ens = sample_paths(msm, source=0, targets=[31, 63], convergence=cc, seed=0)

The returned PathEnsemble has metadata["convergence"] matching SCHEMA section 6: {ess, rhat, n_walkers, n_chunks, converged}. If the walker budget is exhausted without convergence, ConvergenceError is raised; the last-known ess and rhat are in the exception message.

  • Tiny maps (n <= 64): chunk_size=200, ess_min=50 is usually enough.
  • 2^10 to 2^14 binary maps: defaults are fine.
  • Larger maps with deep metastable traps: increase max_walkers to 5e6 and consider running multiple seeds to verify R-hat across independent runs.