[remislp]

Bayesian inference over a mechanism's rate constants with the exact HJC missed-events likelihood, following Epstein, Calderhead, Girolami & Sivilotti (2016) Biophys J 111:333–348. It is a new pure-Python module, HJCFIT.likelihood.mcmc, driven by the same fit specification as hjcfit fit. There are no C++ changes.

[mcmc]
n = 20000
burnin = 5000
chains = 4

hjcfit sample my-fit.toml -o post

What's in it

HJCFIT.likelihood.mcmc is duck-typed on a mechanism like HJCFitter, and does not import scalcs.

• LogPosterior is ln prior + ln likelihood over one or more Records. It wraps HJCFitter but never clips a rate to its limits. A point outside the prior gives −∞ without computing the likelihood. A likelihood that can't be computed also gives −∞ and is counted.
• Priors: UniformPrior (the paper's) and LogUniformPrior. UniformPrior.from_mechanism() reads each free rate's limits, and scalcs' default limits are exactly the paper's prior.
• mwg_sample is the paper's pilot: Metropolis-within-Gibbs on log rates, with the Jacobian applied and step sizes tuned during burn-in only.
• adaptive_sample is adaptive Metropolis (Haario et al. 2001). The mixture can be the paper's weighted sum or the Roberts & Rosenthal (2009) random choice, and an optional initial covariance can shape the first proposals.
• Chain holds samples, proposals, acceptance and scales, saves to .npz, and has mode() and kept().
• Diagnostics: effective_sample_size (Geyer's initial monotone sequence, as the paper), significant_lags, autocorrelation and split-R̂ (potential_scale_reduction).
• gaussian_approximation gives N(mode, (−H)⁻¹) from a numerical Hessian. It uses relative steps and one Richardson extrapolation, and reports an error estimate with every Hessian (relative_error). These are the paper's red curves and the ML standard errors.

Fitting layer

• HJCFitter(solver=...) passes the root-finding options (SOLVER_OPTIONS) to every record's Log10Likelihood.

Fit specification and command

• [likelihood] table: those same options, for fit and sample alike.
• [mcmc] table:
  • start = "fit" (the default): search first, then start the chains at the ML estimate, with the Hessian as the first proposal covariance. If the Hessian can't be used, the run says so and falls back to the paper's small step.
  • start = "guess": the paper's pilot from the initial guess.
  • Also sampler, n, burnin, prior, chains, seed, and a few others.
• runner.sample(): chains after the first start from draws of the Gaussian approximation.
  • Parallel chains run as spawned processes, each rebuilding its posterior from the spec (the C++ likelihood can't be pickled) and each limited to its share of the OpenMP threads.
  • It writes <prefix>.json (summary, spec, provenance) and one .npz per chain.
• hjcfit sample SPEC [-o PREFIX] [--chains N] [--processes P] exits with status 2 when split-R̂ > 1.1.

Why the solver options matter

On the paper's three AChR records, its settings (nmax=2, tolerances 1e-12) give an ln L 0.38 higher than HJCFIT's defaults at the same rates. With those settings, a scalcs mechanism through HJCFIT reproduces the value from the paper's own MATLAB code (BICME) exactly: 237158.0646.

Validation

• Unit tests: 288 pass locally (124 new).
  • The samplers are tested on targets whose answers are known: a correlated Gaussian, a log-normal and a box.
  • Removing the Jacobian from the log-space walk makes its tests fail, with the means shifted by exactly −σ² as theory predicts.
  • The statistical tests also pass for ten other seeds.
  • The Hessian is tested on functions with known second derivatives, including rates eight decades apart.
• Reproduction of the paper, kept outside this repo since it uses unpublished data. Figs 3–10 were reproduced. On identical data, our chains agree with the paper's stored chains rate by rate (|z| ≤ 2.3 for the HJC comparison). Two headline results:
  • Real AChR data: the MCMC and Hessian standard deviations agree to 0.96–1.04 for eight rates. α1b and β1b come out non-Gaussian (SD ratios 1.32 and 2.68), as in the paper's Fig. 7.
  • Synthetic data: all ten true rates lie inside their 95% intervals. The adaptive sampler gives about 4× the effective samples per iteration of the pilot (the paper: 3×).
• Performance, measured:
  • 99% of each evaluation is the C++ likelihood: 57, 36 and 29 ms per AChR evaluation on 1, 2 and 4 OpenMP threads.
  • On a 4-core laptop, four parallel chains ran 1.3–1.6× faster than four in turn.

Two things a reviewer should know

• Nothing here was built locally (there is no compiler on the machine). The Python was tested against the 1.1.0 wheel, whose compiled part this branch does not change. This PR's CI is the first real build. mcmc.py is added to PYTHON_MODULES in likelihood/python/CMakeLists.txt.
• A trap in scalcs: mec.theta_unsqueeze() does not rebuild mec.Q; only set_eff() does. HJCFitter and LogPosterior are safe because they call set_eff for every record. Any new code that reads mec.Q directly after setting rates must call it too.

🤖 Generated with Claude Code