Sample the posterior of a distribution's parameters by MCMC

Fits distribution to data with uniform priors spanning the family's parameter constraints (exactly the C# GetParameterConstraints bounds) and returns the full sampler output: per-chain draws, acceptance rates, the MAP estimate, posterior summaries, and Gelman-Rubin / effective-sample-size diagnostics. Seven of the eight ported samplers are available (Gibbs needs a model-specific conditional proposal and is not exposed here).

Usage

mcmc_sample(
  data,
  distribution = "Normal",
  sampler = c("RWMH", "ARWMH", "DEMCz", "DEMCzs", "HMC", "NUTS", "SNIS"),
  iterations = NULL,
  warmup = NULL,
  chains = NULL,
  thinning = NULL,
  seed = 12345,
  initialize = c("MAP", "Randomize")
)

Arguments

data

numeric vector of observations.

distribution

the distribution family name; see distribution_names().

sampler

one of "RWMH" (the default), "ARWMH", "DEMCz", "DEMCzs", "HMC", "NUTS", or "SNIS".

iterations

iterations per chain (sampler default if NULL).

warmup

warm-up iterations discarded from each chain (sampler default if NULL).

chains

number of chains (sampler default if NULL).

thinning

thinning interval (sampler default if NULL).

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 parameters (parameter names), chains (a list of one draws-by-parameters matrix per chain), acceptance_rates, map (parameter values at the posterior mode), map_fitness, posterior_mean, posterior_sd, posterior_median, posterior_lower_ci, posterior_upper_ci, rhat, and ess (all per-parameter vectors).

Details

Runs the ported serial chain driver with the C# default seed, so a seeded run reproduces the C# sampler stream bit-for-bit (and matches corehydropy exactly). The priors are always the constraint-based uniforms; custom priors are not exposed – use univariate_analysis() and friends for the full Bayesian workflow.

Examples

# \donttest{
x <- dist_random(distribution("Normal", c(100, 15)), 200, seed = 42)
fit <- mcmc_sample(x, "Normal", sampler = "DEMCzs", iterations = 2000)
fit$posterior_mean
#> [1] 98.64253 15.76604
fit$rhat
#> [1] 0.9996262 0.9997092
# }