Bayesian inference: sampling the posterior of a mechanism's rate constants (Epstein et al. 2016)

documentation/source/api/python.rst

@@ -15,6 +15,7 @@ Python API Reference
    python/optimization.rst
    python/fitting.rst
    python/fitspec.rst
+   python/mcmc.rst
 
 
 Extras

documentation/source/api/python/fitspec.rst

@@ -18,6 +18,7 @@ the result.
     hjcfit template -o my-fit.toml     a specification to edit
     hjcfit check my-fit.toml           say what it would do; run nothing
     hjcfit fit my-fit.toml -o out.json run it, and keep the result
+    hjcfit sample my-fit.toml -o post  sample the posterior of the rates
 
 ``examples/fit_template.ipynb`` is the same thing as a notebook, over
 ``examples/CH82.toml``.
@@ -185,6 +186,11 @@ The specification
    :members: from_dict, validate
 
 .. autoclass:: SearchSpec
+
+.. autoclass:: LikelihoodSpec
+   :members: solver
+
+.. autoclass:: MCMCSpec
    :members: from_dict, validate
 
 .. autoexception:: SpecError
@@ -215,6 +221,14 @@ Carrying it out
 
 .. autofunction:: write_result
 
+.. autofunction:: sample
+
+.. autoclass:: SampleOutcome
+
+.. autofunction:: samples_as_dict
+
+.. autofunction:: write_samples
+
 .. autoexception:: RunnerError
 
 The command
@@ -230,5 +244,5 @@ environment on Windows.
 
 Exit status is 0 on success, 1 on anything the command can explain -- one line
 beginning with ``hjcfit:`` -- and 2 from ``fit`` when the search did not
-converge. The estimates are still printed in that case; it is the status that
+converge, or from ``sample`` when the chains disagree (split-R-hat above 1.1). The estimates are still printed in that case; it is the status that
 says so, because that is what a script reads.

documentation/source/api/python/fitting.rst

@@ -102,6 +102,12 @@ The fitter
 
 .. autodata:: FAILURE_COST
 
+.. autodata:: SOLVER_OPTIONS
+
+``HJCFitter(..., solver=dict(nmax=2, xtol=1e-12))`` passes them to every
+record's likelihood. Options left out keep the defaults of
+:py:class:`~HJCFIT.likelihood.Log10Likelihood`.
+
 The result
 """"""""""
 

documentation/source/api/python/mcmc.rst

@@ -0,0 +1,195 @@
+.. _python_mcmc_api:
+
+Bayesian inference
+------------------
+
+.. currentmodule:: HJCFIT.likelihood.mcmc
+
+:ref:`python_fitting_api` finds the rate constants that maximise the HJC
+likelihood. This module treats them as random variables instead, and describes
+their posterior distribution: the likelihood of the records multiplied by a
+prior. The approach follows Epstein, Calderhead, Girolami & Sivilotti (2016)
+*Biophys J* 111:333-348, who sampled that posterior by Markov chain Monte
+Carlo using the same exact missed-events likelihood, and it uses their two
+samplers.
+
+Like the fitter, it needs a mechanism from SCALCS (``pip install
+hjcfit[fitting]``), and it is duck-typed on one, so it does not import scalcs
+itself.
+
+The posterior
+"""""""""""""
+
+.. code-block:: python
+
+    import HJCFIT
+    from HJCFIT.likelihood.fitting import Record
+    from HJCFIT.likelihood.mcmc import LogPosterior
+    from scalcs.samples import samples
+
+    bursts = HJCFIT.read_idealized_bursts("CH82", tau=1e-4, tcrit=4e-3)
+    record = Record(conc=100e-9, tres=1e-4, tcrit=4e-3,
+                    groups=tuple(tuple(b) for b in bursts))
+
+    mechanism = samples.CH82()
+    posterior = LogPosterior(mechanism, [record])
+
+    posterior(mechanism.theta())      # ln prior + ln likelihood
+
+.. autoclass:: LogPosterior
+   :members: log_prior, log_likelihood, __call__, k, nevals, nfailures
+
+Two things differ from fitting on purpose:
+
+* **Nothing is clipped.** :py:class:`~HJCFIT.likelihood.fitting.HJCFitter`
+  resets a rate that leaves its limits, which suits a search. For a
+  posterior, a point outside the prior has zero density. Moving it to the
+  boundary instead would pile probability there.
+* **A likelihood that cannot be computed gives** ``-inf``, **not a penalty.**
+  A sampler rejects the move. :py:attr:`LogPosterior.nfailures` counts these,
+  so a mechanism that the likelihood struggles with is visible.
+
+Priors
+""""""
+
+.. autoclass:: UniformPrior
+   :members: from_mechanism, logpdf, contains, sample, k
+
+.. autoclass:: LogUniformPrior
+
+The default prior is :py:meth:`UniformPrior.from_mechanism`: flat between the
+limits each free rate already carries. Those are the limits a
+maximum-likelihood fit resets against, so the mode of the posterior is the
+maximum-likelihood estimate. They are also the prior of Epstein et al.
+(2016): scalcs' default limits are U(1e-2, 1e6) s\ :sup:`-1`, and
+U(1e-2, 1e10) M\ :sup:`-1` s\ :sup:`-1` for association rates.
+
+Solver settings matter
+""""""""""""""""""""""
+
+Pass ``solver=`` to reproduce a published value. On that paper's three AChR
+records, its root-finding settings (``nmax=2``, tolerances 1e-12) give a
+natural log-likelihood 0.38 higher than the defaults at the same rates. See
+:py:data:`~HJCFIT.likelihood.fitting.SOLVER_OPTIONS`.
+
+Sampling
+""""""""
+
+Two stages, as in that paper. A Metropolis-within-Gibbs pilot walks the
+logarithm of one rate at a time and finds the mode. An adaptive Metropolis
+chain, started at that mode, learns the posterior covariance and samples
+along correlations that the pilot crosses only slowly:
+
+.. code-block:: python
+
+    from HJCFIT.likelihood.mcmc import (
+        adaptive_sample, effective_sample_size, mwg_sample)
+
+    pilot = mwg_sample(posterior, mechanism.theta(), n=10000, burnin=5000,
+                       rng=1)
+    start, _ = pilot.mode()
+    chain = adaptive_sample(posterior, start, n=100000, burnin=50000, rng=2)
+
+    draws = chain.kept()                 # (50000, k) rates, burn-in dropped
+    effective_sample_size(draws)         # one value per rate
+
+Step sizes are tuned only during burn-in, so the samples that are kept come
+from a chain with a fixed kernel. The adaptive covariance keeps learning
+afterwards, but with a weight that falls as :math:`1/n`, which is what keeps
+the chain ergodic (Haario et al. 2001).
+
+A good first check of any chain is to plot ``chain.log_posterior`` against
+the iteration number, and ``chain.proposals`` against ``chain.samples`` for a
+correlated pair. Proposals falling well outside the cloud of accepted samples
+mean the proposal does not fit the posterior (Epstein et al. 2016, Fig. 3).
+
+.. autofunction:: mwg_sample
+
+.. autofunction:: adaptive_sample
+
+.. autoclass:: Chain
+   :members: n, kept, acceptance_rate, mode, save, load
+
+How much a chain is worth
+"""""""""""""""""""""""""
+
+Samples from a Markov chain are correlated, so :math:`N` of them carry less
+information than :math:`N` independent draws. The effective sample size is how
+many independent draws they are worth.
+
+.. autofunction:: effective_sample_size
+
+.. autofunction:: significant_lags
+
+.. autofunction:: autocorrelation
+
+The Gaussian approximation
+""""""""""""""""""""""""""
+
+Maximum-likelihood standard errors come from the curvature of the log
+likelihood at its maximum. The estimate is taken to be normally distributed
+about the maximum, with covariance equal to the inverse of the negative
+Hessian. :py:func:`gaussian_approximation` computes that distribution at any
+mode: a fit's estimate, or the best sample of a chain. Overlaid on the
+marginals a sampler finds, it shows where standard errors would misstate the
+uncertainty. For the AChR data in Epstein et al. (2016, Fig. 7), that happens
+for the monoliganded opening and shutting rates of site B.
+
+.. code-block:: python
+
+    from HJCFIT.likelihood.mcmc import gaussian_approximation
+
+    mode, _ = chain.mode()
+    approx = gaussian_approximation(posterior, mode)
+    approx.sd, approx.correlation      # standard errors and correlations
+    approx.relative_error              # how well the Hessian is determined
+
+.. autofunction:: gaussian_approximation
+
+.. autoclass:: GaussianApproximation
+   :members: sd, correlation, relative_error, marginal_pdf, sample
+
+.. autofunction:: hessian
+
+A Hessian needs a record that determines every rate. On the CH82 sample
+record, a single concentration fitted with eight free rates, one direction is
+so flat that the likelihood's own numerical noise decides its curvature.
+``relative_error`` is then large, and at small steps the point may not even
+appear to be a maximum. That is a statement about identifiability, not about
+the differentiation: fitting several concentrations at once removes it.
+
+From a specification
+""""""""""""""""""""
+
+Everything above is available from a fit specification (see
+:ref:`python_fitspec_api`), with no Python written. Add an ``[mcmc]`` section
+and, to reproduce a published value, a ``[likelihood]`` section with the
+root-finding settings it was computed with::
+
+    [likelihood]
+    nmax = 2
+    xtol = 1e-12
+    rtol = 1e-12
+
+    [mcmc]
+    n = 20000        # iterations per chain
+    burnin = 5000    # of which tuning, not kept
+    chains = 4       # in parallel processes
+
+and run it::
+
+    hjcfit sample my-fit.toml -o post
+
+By default this fits first and starts every chain at the maximum-likelihood
+estimate. The inverse Hessian there shapes the first proposals, at the optimal
+scale. ``start = "guess"`` runs the paper's own pilot sampler from the initial
+guess instead. Chains after the first start from draws of the Gaussian
+approximation, so that chains which agree have had a chance not to. The
+command prints, for each rate, the posterior mean, standard deviation, 95%
+interval, effective sample size and split-:math:`\hat R`, and writes
+``post.json`` (all of that, the specification and the provenance) and one
+``post_chain<i>.npz`` per chain (:py:meth:`Chain.load` reads them). It exits
+with status 2 when :math:`\hat R` exceeds 1.1: the chains plainly disagree, so
+run them longer.
+
+.. autofunction:: potential_scale_reduction

likelihood/python/CMakeLists.txt

@@ -70,7 +70,8 @@ if(pythonBindings)
   # Named rather than globbed, because CMake will not re-run on a new file
   # appearing and a glob would make the wheel depend on when cmake last ran.
   set(PYTHON_MODULES __init__.py random.py _methods.py optimization.py
-                     simplex_hjc.py fitting.py fitspec.py runner.py cli.py)
+                     simplex_hjc.py fitting.py fitspec.py runner.py cli.py
+                     mcmc.py)
 
   install(FILES  ${PYTHON_MODULES}
                  ${CMAKE_CURRENT_BINARY_DIR}/likelihood.py

likelihood/python/cli.py

@@ -22,6 +22,7 @@
         hjcfit template -o ch82.toml     a specification to edit
         hjcfit check ch82.toml           say what it would do; run nothing
         hjcfit fit ch82.toml -o out.json run it, and keep the result
+        hjcfit sample ch82.toml -o post  sample the posterior ([mcmc])
 
     Also reachable as ``python -m HJCFIT.likelihood.cli`` when the console
     script is not on the path, which is the usual state of affairs inside a
@@ -80,6 +81,20 @@ def _parser():
     fit.add_argument('-q', '--quiet', action='store_true',
                      help='print the result and nothing before it')
 
+    sample = sub.add_parser(
+        'sample', help='sample the posterior of the rates, as [mcmc] says')
+    _add_spec_argument(sample)
+    sample.add_argument('-o', '--output', metavar='PREFIX', default=None,
+                        help='write PREFIX.json (summary, provenance) and '
+                             'PREFIX_chain<i>.npz (the chains)')
+    sample.add_argument('--chains', type=int, default=None,
+                        help='number of chains, overriding [mcmc] chains')
+    sample.add_argument('--processes', type=int, default=None,
+                        help='worker processes (default: one per chain, up '
+                             'to the number of CPUs; 1 runs them in turn)')
+    sample.add_argument('-q', '--quiet', action='store_true',
+                        help='print the result and nothing before it')
+
     show = sub.add_parser(
         'show', help='read a specification and write it back, normalised')
     _add_spec_argument(show)
@@ -162,7 +177,8 @@ def _check(args):
     # only thing here that proves the records, the mechanism and the dead time
     # are mutually possible. A guess at which the likelihood cannot be
     # computed is a guess no search can start from.
-    fitter = HJCFitter(mec, records, log_params=spec.search.log_params)
+    fitter = HJCFitter(mec, records, log_params=spec.search.log_params,
+                       solver=spec.likelihood.solver() or None)
     print('')
     print('log10L at the guess: {0:.4f}'.format(fitter.log10_likelihood()))
     print('nothing was fitted; run: hjcfit fit {0}'.format(args.spec))
@@ -191,6 +207,36 @@ def _fit(args):
     return 0 if outcome.result.success else 2
 
 
+def _sample(args):
+    from dataclasses import replace
+
+    from .runner import sample, write_samples
+
+    spec = _load(args.spec)
+    if args.chains is not None:
+        if args.chains < 1:
+            raise SpecError('--chains: must be at least 1, not {0}'
+                            .format(args.chains))
+        spec = replace(spec, mcmc=replace(spec.mcmc, chains=args.chains))
+    if not args.quiet:
+        print(spec)
+        if not spec.mcmc.as_dict():
+            print('  sampling: {0} (the defaults; no [mcmc] section)'.format(
+                spec.mcmc))
+        print('')
+    outcome = sample(spec, processes=args.processes, verbose=not args.quiet)
+    if not args.quiet:
+        print('')
+    print(outcome)
+    if args.output:
+        paths = write_samples(args.output, outcome)
+        print('')
+        print('wrote {0}'.format(", ".join(paths)))
+    # As for fit: numbers are printed either way, and the status says whether
+    # to trust them. Chains that plainly disagree have not converged.
+    return 0 if not outcome.max_rhat > 1.1 else 2
+
+
 def _show(args):
     sys.stdout.write(_load(args.spec).to_toml())
     return 0
@@ -202,7 +248,8 @@ def main(argv=None):
         :param argv: Arguments, without the program name. None takes
             ``sys.argv``.
         :returns: 0 on success, 1 on anything this module can explain, and 2
-            from ``fit`` when the search did not converge.
+            from ``fit`` when the search did not converge or from ``sample``
+            when the chains disagree (split-R-hat above 1.1).
     """
     parser = _parser()
     args = parser.parse_args(argv)
@@ -215,7 +262,7 @@ def main(argv=None):
         return 0
 
     handlers = {'template': _template, 'check': _check, 'fit': _fit,
-                'show': _show}
+                'sample': _sample, 'show': _show}
     try:
         return handlers[args.command](args)
     except SpecError as error: