mcmc_posterior

mcmc_posterior(
    log_likelihood,
    priors,
    sampler='RWMH',
    proposal=None,
    gradient=None,
    iterations=None,
    warmup=None,
    chains=None,
    thinning=None,
    output_length=None,
    seed=12345,
    initialize='MAP',
)

Sample your own posterior by MCMC.

Runs any of the ported MCMC samplers against a log-likelihood you write and priors you choose. This is the constructor the upstream C# library itself exposes, MCMCSampler(priorDistributions, logLikelihoodFunction): :func:corehydropy.mcmc_sample can only fit a built-in distribution family under uniform priors spanning its parameter constraints, while this function takes any model you can write down.

log_likelihood is called with one list of floats, as long as priors, and must return a single number. Following the upstream contract, it should return the log of the data likelihood PLUS the log prior density; the priors list is used for the feasible parameter bounds and for chain initialization, and is never added to your value behind your back. A parameter outside its prior’s support is rejected before your function sees it, so a flat (uniform) prior needs no term of its own.

Parameters

Name	Type	Description	Default
log_likelihood	callable	Takes a list of floats, returns one float.	required
priors	Distribution or sequence of Distribution	One prior per parameter, in the order log_likelihood reads them.	required
sampler	('RWMH', 'ARWMH', 'DEMCz', 'DEMCzs', 'HMC', 'NUTS', 'SNIS', 'Gibbs')	The MCMC sampler. "Gibbs" requires proposal.	"RWMH"
proposal	callable	The conditional proposal function "Gibbs" samples with, and the only sampler that takes one. It is called as proposal(parameters, rng) and must return a sequence as long as priors (a bare number when there is one parameter): the next state of the chain, which Gibbs accepts unconditionally. Ordinarily that is a draw from the full conditional of the model. rng is a :class:corehydropy.Rng handle on the generator this chain is running on – draw from it with rng.uniform() and rng.integers(), not with :mod:random or :mod:numpy.random, or the seeded run stops being reproducible.	None
gradient	callable	An analytic gradient of log_likelihood for "HMC" and "NUTS", the only samplers that take one. It is called as gradient(parameters) and must return a sequence as long as priors (a bare number when there is one parameter). Left unset, both samplers use the ported bound-aware finite-difference gradient, which costs two extra log_likelihood calls per parameter per leapfrog step; an analytic gradient is usually a large saving and is always more accurate.	None
iterations	int	Iterations per chain (sampler default if omitted). "SNIS" needs at least 10000: its ported validation requires iterations to be at least the output length, whose default is 10000, so lowering output_length lowers the floor with it.	None
warmup	int	Warm-up iterations discarded from each chain (sampler default if omitted). When iterations is given and warmup is not, half of iterations is used, matching :func:corehydropy.mcmc_sample. "SNIS" is the exception: its ported validation rejects any warm-up at all, so none is derived for it.	None
chains	int	Number of chains (sampler default if omitted). "DEMCz" and "DEMCzs" require at least three; "SNIS" draws independently and runs one.	None
thinning	int	Thinning interval (sampler default if omitted).	None
output_length	int	Total number of retained draws across all chains (sampler default, 10,000, if omitted). The ported sampler collects ceil(output_length / chains) draws per chain after the iteration loop, so this is the second setting – with thinning – that multiplies how many times your function is called. The ported floor is 100 and it is refused below that.	None
seed	int	PRNG seed; 12345 is the C# default.	12345
initialize	('MAP', 'Randomize')	Chain initialization: from the posterior-mode estimate (the C# default) or randomized draws from the priors.	"MAP"

Returns

Name	Type	Description
	dict	The same fields :func:corehydropy.mcmc_sample returns: parameters (parameter names, p1, p2, … since your model names nothing), chains (a list of one (n_draws, n_params) array per chain), acceptance_rates, map, map_fitness, posterior_mean, posterior_sd, posterior_median, posterior_lower_ci, posterior_upper_ci, rhat, and ess.

Notes

Reproducing a run in R. A seeded :func:corehydropy.mcmc_sample run is bit-identical between Python and R, because every arithmetic operation happens in the shared C++ core. That guarantee is WEAKER here, and it is worth stating plainly. The draws still come from the core’s seeded Mersenne Twister, but the log-density is your own Python code, and Python and R do not guarantee identical rounding for the same formula. MCMC amplifies a single differing bit far harder than an optimizer does: one flipped accept-or-reject changes every state after it, so the two chains diverge outright rather than drift apart slowly.

A seeded run reproduces across the two languages if and only if your function returns bit-identical values. Arithmetic (+ - * /) is IEEE-deterministic and does reproduce; log, exp, gamma and friends come from each platform’s own math library and are not guaranteed to. Note also that R’s sum() accumulates in extended precision while Python’s does not, so a plain loop is the portable spelling on both sides.

Performance. Every evaluation calls back into Python, and there are far more of them than iterations suggests: the count is (iterations + output_length / chains) * thinning * chains, and the ported defaults are 4 chains, a thinning interval of 20 and an output length of 10,000. So iterations=10000 is a million crossings, not ten thousand.

That is affordable. All the figures here come from one machine, one session and one problem – a 50-point Gaussian log-density at iterations=10000, so exactly the million evaluations above – and corehydror’s help page reports the same experiment, so the two can be read side by side. The run took 7.1 seconds, against 1.4 seconds for the :func:corehydropy.mcmc_sample call on the same data: the callback path is about five times the built-in one, and almost none of that is the boundary. The same function called a million times directly from Python takes 6.4 seconds, so your own code is nearly the whole total. The crossing itself is about 0.5 microseconds, measured by replacing the log-density with lambda p: 0.0, which brings the run to 0.6 seconds.

Two settings dominate the count and neither is obvious. thinning multiplies it, so thinning=1 turned the same run into 0.33 seconds. And initialize="MAP", the C# default, runs the DifferentialEvolution optimizer over your function before the first chain iteration; initialize="Randomize" skips the optimizer.

"Gibbs" is the sampler whose defaults surprise: the ported constructor sets one chain, no thinning, and 100,000 iterations on top of a 10,000-draw output block, so an iterations you do not set is 110,000 iterations of BOTH your log-likelihood and your proposal. Set iterations.

Interrupting a long run. Ctrl-C returns control with a KeyboardInterrupt, but not instantly. The ported samplers have no cancellation hook, so the chain runs to the end of its loop – rejecting every remaining point without calling your function again – before the interrupt surfaces. The wait is therefore set by how much of the run is left, not by the interrupt: measured on a 200,000-iteration chain interrupted one second in, 6.4 seconds from the signal to the exception reaching the caller. corehydror behaves the same way and reports the same measurement at 100,000 iterations.

See Also

corehydropy.mcmc_sample : a built-in family under constraint-based priors, which is faster and bit-identical across languages.

Examples

>>> import corehydropy as ch
>>> data = [4.9, 5.1, 5.0, 5.2, 4.8]
>>> # A plain loop and ``+ - * /`` alone, the portable spelling described above.
>>> def ll(p):
...     acc = 0.0
...     for x in data:
...         acc += (x - p[0]) * (x - p[0])
...     return -0.5 * acc
>>> fit = ch.mcmc_posterior(ll, [ch.Distribution("Uniform", [0.0, 10.0])],
...                         iterations=200, seed=12345)