Bayesian inference

Fitting a mechanism 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

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

class HJCFIT.likelihood.mcmc.LogPosterior(mec, records, prior=None, solver=None)

The log posterior density of a mechanism’s free rate constants.

Calling it with a vector of free rates returns \(\ln p(\theta) + \ln L(\theta)\). The prior is evaluated first, and the likelihood is not computed at all where the prior is zero.

Parameters:

• mec – A mechanism carrying its constraints. It is modified: every evaluation puts its rates on the mechanism, as fitting does.

• records – A sequence of Record, fitted simultaneously, each at its own concentration.

• prior – Anything with logpdf(rates) and k. Defaults to UniformPrior.from_mechanism().

• solver (dict) – Root-finding options for the likelihood; see SOLVER_OPTIONS.

__call__(rates)

\(\ln p(\theta) + \ln L(\theta)\).

property k

Number of free parameters.

log_likelihood(rates)

\(\ln L(\theta)\) summed over the records; -inf if it cannot be computed. Rates are used as given, never clipped.

log_prior(rates)

\(\ln p(\theta)\).

nevals

Likelihood evaluations requested, failures included.

nfailures

Likelihood evaluations that could not be computed.

Two things differ from fitting on purpose:

• Nothing is clipped. 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. LogPosterior.nfailures counts these, so a mechanism that the likelihood struggles with is visible.

Priors

class HJCFIT.likelihood.mcmc.UniformPrior(lower, upper, names=None)

Independent uniform distributions on the free rate constants.

This is the prior of Epstein et al. (2016). It is flat between the limits a maximum-likelihood fit would also impose, so the mode of the posterior is the maximum-likelihood estimate. The paper uses U(1e-2, 1e6) s^-1 for every rate and U(1e-2, 1e10) M^-1 s^-1 for association rates, which are the limits scalcs gives its rates by default.

Parameters:

• lower – Lower bound of each free rate.

• upper – Upper bound of each free rate.

• names – Names of the free rates, in the same order. Optional.

contains(rates)

Whether rates lie inside the bounds, boundaries included.

classmethod from_mechanism(mec)

The prior bounded by each free rate’s own limits.

Raises:

ValueError – if a free rate has no limits.

property k

Number of parameters.

logpdf(rates)

Log density at rates: -inf outside the bounds.

sample(rng, size=None)

Draws from the prior.

Parameters:

• rng – A numpy.random.Generator.

• size – Number of draws; None for a single vector.

class HJCFIT.likelihood.mcmc.LogUniformPrior(lower, upper, names=None)

Independent log-uniform distributions on the free rate constants.

Uniform in the logarithm of each rate between the bounds: every decade is equally likely a priori. This is the scale-free alternative to UniformPrior, which puts almost all of its mass in the top decade. The bounds must be positive.

The default prior is 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^-1, and U(1e-2, 1e10) M^-1 s^-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 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:

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 \(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).

HJCFIT.likelihood.mcmc.mwg_sample(log_density, x0, n, burnin, rng=None, *, initial_scale=1.0, tune_every=50, acceptance=(0.1, 0.5), tune_step=0.1, log_space=True, callback=None, callback_every=1000)

Metropolis-within-Gibbs: a random walk on one parameter at a time.

Each iteration visits every parameter in turn and proposes \(y_i' = y_i + s_i z\), \(z \sim N(0, 1)\), on the logarithm of the rate (a multiplicative step on the rate itself), accepting it by the Metropolis rule. This is the pilot sampler of Epstein et al. (2016): multiplicative steps suit rates spanning several decades, and it locates the posterior mode reliably from a poor starting point.

During burn-in, every tune_every iterations, each \(s_i\) is multiplied by 1 - tune_step if that parameter’s acceptance over the window is below acceptance[0], or by 1 + tune_step if above acceptance[1].

Parameters:

• log_density – Callable returning the log posterior at a vector of rates, such as a LogPosterior.

• x0 – Starting rates; the density there must be positive.

• n (int) – Iterations, each a sweep over all parameters.

• burnin (int) – Iterations during which the step sizes are tuned.

• rng – A numpy.random.Generator, or a seed.

• initial_scale – Starting \(s_i\), a scalar or one per parameter. On log rates, 1 is a step of about a factor of e.

• log_space (bool) – Walk the logarithms of the rates (default) or the rates themselves.

• callback – Called as callback(iteration, chain) every callback_every iterations with the chain so far.

Returns:

A Chain with per-parameter accepted and scales.

HJCFIT.likelihood.mcmc.adaptive_sample(log_density, x0, n, burnin, rng=None, *, adapt_start=100, beta=0.05, initial_step=0.1, initial_covariance=None, mixture='sum', tune_every=50, acceptance=(0.1, 0.5), tune_step=0.1, log_space=False, callback=None, callback_every=1000)

Adaptive Metropolis: a block random walk that learns its covariance.

For the first adapt_start iterations the proposal is \(x + (\epsilon/\sqrt{k})\,L_0 z\), where \(\epsilon\) is initial_step and \(L_0 L_0^T\) is initial_covariance, or the identity if none is given. From then on \(\hat\Sigma\), the covariance of every sample so far, takes over. It is scaled by \(2.38^2/k\), which is optimal for a Gaussian target (Roberts & Rosenthal 2001). It is also mixed with the small step, so that the chain cannot lock onto a degenerate \(\hat\Sigma\):

• mixture='sum', as in Epstein et al. (2016): \(x' = x + (1-\beta)\,L z_1 + \beta\,(\epsilon/\sqrt{k})\,L_0 z_2\), where \(L L^T = (2.38^2/k)\,s\,\hat\Sigma\).

• mixture='choice', as in Roberts & Rosenthal (2009): the small step with probability \(\beta\), otherwise the learned one.

Both are symmetric. During burn-in, a global multiplier \(s\) on \(\hat\Sigma\) is tuned from the acceptance rate every tune_every iterations, as in mwg_sample().

By default the walk is on the rates themselves, as in that paper; start it at the mode mwg_sample() found. With log_space=True it walks their logarithms instead.

Parameters:

• log_density – Callable returning the log posterior at a vector of rates, such as a LogPosterior.

• x0 – Starting rates; the density there must be positive.

• n (int) – Iterations.

• burnin (int) – Iterations during which s is tuned. The covariance keeps adapting afterwards. The weight of each new sample in it falls as \(1/n\), which is what keeps the chain ergodic (Haario et al. 2001).

• rng – A numpy.random.Generator, or a seed.

• adapt_start (int) – Iterations before the learned covariance is used.

• beta (float) – Weight of the small step.

• initial_covariance – \(L_0 L_0^T\). The inverse Hessian at a maximum-likelihood fit is a good choice. It shapes the small step throughout.

• callback – Called as callback(iteration, chain) every callback_every iterations with the chain so far.

Returns:

A Chain; proposal_covariance is the learned covariance at the end, scaled as it was being used.

class HJCFIT.likelihood.mcmc.Chain(samples: ndarray, log_posterior: ndarray, proposals: ndarray, accepted: ndarray, scales: ndarray, burnin: int, names: tuple = None, sampler: str = '', settings: dict = <factory>, seconds: float = 0.0, nevals: int = 0, nfailures: int = 0, proposal_covariance: ndarray = None)

What a sampler produced.

Parameters:

• samples – (n, k) parameter values, one row per iteration, as rates even when the sampler walked their logarithms.

• log_posterior – (n,) log density at each row of samples.

• proposals – (n, k) the points proposed, accepted or not, as rates. Plotting them against samples shows how well a proposal fits the posterior (Epstein et al. 2016, Fig. 3).

• accepted – Whether each proposal was accepted: (n,) for a block sampler, (n, k) for one parameter at a time.

• scales – The step-size multipliers in use at each iteration, (n,) or (n, k) like accepted.

• burnin – Iterations during which the sampler tuned itself. They are not draws from the posterior; kept() drops them.

• names – The parameters’ names, if the target had them.

• sampler – Which sampler produced the chain.

• settings – The sampler’s settings, for provenance.

• seconds – Wall-clock time of the run.

• nevals – Calls to the log density.

• nfailures – Of those, how many the target counted as failures.

• proposal_covariance – The proposal covariance at the end of an adaptive run, in the space it walked.

acceptance_rate(after_burnin=True)

Fraction of proposals accepted; per parameter for one-at-a-time sampling.

kept()

The samples after burn-in: the draws from the posterior.

classmethod load(path)

Read a chain written by save().

mode()

The sample with the highest log posterior, and that value. This is where Epstein et al. (2016) start the adaptive sampler.

property n

Iterations, burn-in included.

save(path)

Write the chain to a NumPy .npz file.

How much a chain is worth

Samples from a Markov chain are correlated, so \(N\) of them carry less information than \(N\) independent draws. The effective sample size is how many independent draws they are worth.

HJCFIT.likelihood.mcmc.effective_sample_size(x, max_lag=None)

Effective sample size by Geyer’s (1992) initial monotone sequence.

Autocorrelations are summed in adjacent pairs, \(\Gamma_m = \rho_{2m} + \rho_{2m+1}\). The pairs are made non-increasing and kept while positive, and then \(\text{ESS} = N / (-1 + 2\sum_m \Gamma_m)\). The denominator is floored at 1, so the estimate never exceeds N. This is the estimator of Epstein et al. (2016).

Parameters:

• x – A series, or an (n, k) array of them.

• max_lag – Autocorrelations beyond this lag are not used. None takes each series’ own significant_lags(). That paper truncates every parameter at the significant lags of \(\alpha_2\); pass that number to do the same.

Returns:

A float, or one per column.

HJCFIT.likelihood.mcmc.significant_lags(x, max_lag=200, nsd=2.0)

The first lag whose autocorrelation falls below \(\text{nsd}/\sqrt{N}\), the approximate bound for a series with no autocorrelation; max_lag if none does.

Epstein et al. (2016, Table 2) report this as the number of significant lags, and truncate the effective sample size there.

HJCFIT.likelihood.mcmc.autocorrelation(x, max_lag)

Sample autocorrelation of a series at lags 0 .. max_lag.

The biased estimator: autocovariances divided by the length of the series, then by the variance. It is the usual estimator for MCMC output, and the one MATLAB’s autocorr computes.

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. 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.

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

HJCFIT.likelihood.mcmc.gaussian_approximation(log_density, mode, rel_step=0.003, steps=None)

The GaussianApproximation to a posterior at its mode.

Parameters:

• log_density – Callable returning the log posterior at a vector of rates, such as a LogPosterior.

• mode – The mode: a fit’s maximum-likelihood estimate, or the best sample of a chain (Chain.mode()). The Hessian is evaluated there, so a point that is not a maximum is refused.

• rel_step – See hessian().

• steps – See hessian().

Raises:

ValueError – if the negative Hessian is not positive definite, so the point is not a maximum; or see hessian().

class HJCFIT.likelihood.mcmc.GaussianApproximation(mode: ndarray, covariance: ndarray, hessian: ndarray, hessian_error: ndarray, log_density: float, names: tuple = None)

A normal distribution matched to the posterior at its mode.

The mean is the mode and the covariance is \((-H)^{-1}\), the inverse of the negative Hessian of the log posterior there. With a flat prior, the mode is the maximum-likelihood estimate and this is the asymptotic distribution of that estimate: the basis of maximum-likelihood standard errors and correlations (Colquhoun, Hatton & Hawkes 2003). Epstein et al. (2016, Figs 4 and 7) overlay it on the marginals the sampler finds. Where the two differ, the standard errors misstate the uncertainty.

Parameters:

• mode – Where the approximation was made.

• covariance – \((-H)^{-1}\).

• hessian – \(H\).

• hessian_error – Estimated error in each element of hessian.

• log_density – The log density at the mode.

• names – Parameter names, if the density had them.

property correlation

The correlation matrix.

marginal_pdf(parameter, x)

The normal density of one parameter, by index or name, at x.

property relative_error

The largest estimated error in the Hessian, relative to the scale of its row and column: \(\max_{ij} \epsilon_{ij} / \sqrt{|H_{ii} H_{jj}|}\). Small (1e-4 or less) when the Hessian is well determined. Near 1, a direction is barely curved at all, and the noise of the likelihood decides it. A single-concentration record often does this, because it does not identify every rate.

sample(rng, size=None)

Draws from the approximation.

property sd

Standard deviations: square roots of the covariance diagonal.

HJCFIT.likelihood.mcmc.hessian(f, x, rel_step=0.003, steps=None)

Hessian of a scalar function by central differences with one Richardson extrapolation.

Second differences with steps \(h\) and \(h/2\) each carry an \(O(h^2)\) error. \((4 D_{h/2} - D_h) / 3\) cancels it, leaving \(O(h^4)\), and the difference between that and \(D_{h/2}\) estimates the error left.

Steps are relative, \(h_i = \text{rel\_step} \times |x_i|\), because rate constants span eight decades: one absolute step cannot suit a rate of 1 s^-1 and one of 10⁸ M^-1 s^-1 at once.

The step is a trade-off. Too large, and the extrapolation has not yet removed the curvature error. Too small, and noise in f is divided by \(h^2\). The HJC likelihood carries noise from its root finding: about 5e-8 in ln L with the default tolerances. On the three AChR records of Epstein et al. (2016), steps from 3e-2 down to 1e-3 give the same standard deviations, with the smallest error estimate near the default. Tightening the likelihood’s solver tolerances lowers the noise floor.

Parameters:

• f – Scalar function of a vector.

• x – Where to differentiate.

• rel_step (float) – Step relative to each coordinate.

• steps – Absolute steps, one per coordinate, instead.

Returns:

(H, error), both k x k.

Raises:

ValueError – if f is not finite at any point used. A difference across a prior boundary or a failed likelihood would otherwise give a meaningless number without complaint.

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 A fit from a file), 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-\(\hat R\), and writes post.json (all of that, the specification and the provenance) and one post_chain<i>.npz per chain (Chain.load() reads them). It exits with status 2 when \(\hat R\) exceeds 1.1: the chains plainly disagree, so run them longer.

HJCFIT.likelihood.mcmc.potential_scale_reduction(chains)

Split-\(\hat R\) (Gelman et al. 2013) for one parameter.

Each chain is split in half, and the variance between the halves’ means is compared with the variance within them. Near 1 when every chain samples the same distribution; values above about 1.01 mean the chains have not yet agreed. Splitting makes a single chain testable too, against a drift between its first and second halves.

Parameters:

chains – A sequence of 1-d series, one per chain, burn-in already removed. They are truncated to the shortest.

Returns:

\(\hat R\).