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

Usage

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

Arguments

log_likelihood: a function taking a numeric parameter vector and returning a single number.
priors: a list of distribution() objects, one per parameter, in the order log_likelihood reads them. A single distribution is accepted for a one-parameter model.
sampler: one of "RWMH" (the default), "ARWMH", "DEMCz", "DEMCzs", "HMC", "NUTS", "SNIS", or "Gibbs". "Gibbs" requires proposal.
proposal: 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 numeric vector as long as priors: 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 handle on the generator this chain is running on – draw from it with rng_uniform() and rng_integers(), not with stats::runif(), or the seeded run stops being reproducible.
gradient: 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 numeric vector as long as priors. Left NULL, 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.
iterations: iterations per chain (sampler default if NULL). "SNIS" needs at least 10000 unless output_length is lowered too: its ported validation requires iterations to be at least the output length, whose default is 10000.
warmup: warm-up iterations discarded from each chain (sampler default if NULL). When iterations is given and warmup is not, half of iterations is used, matching mcmc_sample(). "SNIS" is the exception: its ported validation rejects any warm-up at all, so none is derived for it.
chains: number of chains (sampler default if NULL). "DEMCz" and "DEMCzs" require at least three; "SNIS" draws independently and runs one.
thinning: thinning interval (sampler default if NULL).
output_length: total number of retained draws across all chains (sampler default, 10,000, if NULL). The ported sampler collects ceiling(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.
seed: PRNG seed; 12345 is the C# default.
initialize: chain initialization: "MAP" (from the posterior-mode estimate, the C# default) or "Randomize" (draws from the priors).

Value

A list with the same fields mcmc_sample() returns: parameters (parameter names, p1, p2, ... since your model names nothing), chains (a list of one draws-by-parameters matrix per chain), acceptance_rates, map, map_fitness, posterior_mean, posterior_sd, posterior_median, posterior_lower_ci, posterior_upper_ci, rhat, and ess.

Details

log_likelihood is called with one numeric vector, 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.

Reproducing a run in Python

A seeded mcmc_sample() run is bit-identical between R and Python, 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 R code, and R and Python 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 an explicit loop or Reduce() is the portable spelling.

Performance

Every evaluation calls back into R, 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 corehydropy's help page reports the same experiment, so the two can be read side by side. The run took 5.0 seconds, against 1.3 seconds for the mcmc_sample() call on the same data: the callback path is about four times the built-in one, and most of that is not the boundary. The same closure called a million times directly from R takes 2.7 seconds, so your own code is over half the total. The crossing itself is about 1.5 microseconds, measured by replacing the log-density with function(p) 0, which brings the run to 1.7 seconds.

Two settings dominate the count and neither is obvious. thinning multiplies it, so thinning = 1 turned the same run into 0.23 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 an interrupt condition, 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 100,000-iteration chain interrupted one second in, 3.2 seconds from Ctrl-C to the prompt. corehydropy behaves the same way and reports the same measurement at 200,000 iterations.

Examples

# \donttest{
set.seed(1)
x <- rnorm(50, mean = 5)
# A plain loop and `+ - * /` alone, the portable spelling described above.
ll <- function(p) {
  acc <- 0
  for (xi in x) acc <- acc + (xi - p[1]) * (xi - p[1])
  -0.5 * acc
}
fit <- mcmc_posterior(ll, list(distribution("Uniform", c(0, 10))),
                      iterations = 500, seed = 12345)
fit$posterior_mean
#> [1] 5.100834
# }