Bayesian MCMC fit

Fit a model with a Bayesian MCMC analysis and return a fit object carrying the raw chains, the posterior summary (mean/median/sd/credible interval, R-hat, effective sample size), and the usual Bayesian goodness-of-fit scalars. Wraps the shared C++ BayesianAnalysis ported from USACE-RMC RMC.BestFit – the same estimator univariate_analysis() and estimation_diagnostics() build on.

Usage

fit_bayesian(
  model,
  distribution = NULL,
  sampler = "DEMCz",
  chains = 4L,
  iterations = 3000L,
  warmup = NULL,
  output_length = 10000L,
  thinning_interval = -1L,
  seed = 12345L,
  point_estimator = NULL,
  credible_level = 0.9,
  ...
)

Arguments

model: a model_univariate() (or any model_*()) object, or a plain numeric vector of observations together with distribution. A model can bring censored observations (see analysis_data()), nonstationary trends (see trend()), and parameter bounds or priors (see model_parameter()).
distribution: distribution family name, required only when model is a numeric vector.
sampler: MCMC sampler: "DEMCz" (default), "DEMCzs", "ARWMH", or "NUTS". BayesianAnalysis can only construct these four; RWMH, HMC and SNIS are real MCMC samplers but need mcmc_sample() instead.
chains: number of parallel Markov chains.
iterations: number of post-warmup MCMC iterations, per chain.
warmup: number of warmup (burn-in) iterations. Defaults to max(50, iterations %/% 2) when omitted: BayesianAnalysis's own class default (1500) would otherwise silently trip its sampler's warmup <= iterations / 2 guard for any iterations below about 3000.
output_length: number of posterior draws retained (thinned down from the raw chains) for the summary and any downstream uncertainty quantification.
thinning_interval: MCMC thinning interval; -1 (default) keeps the sampler's own default.
seed: PRNG seed for the sampler (fixed for reproducibility – a seeded call returns identical draws in R and Python).
point_estimator: which posterior summary $parameters reports: "PosteriorMean" or "PosteriorMode" (the MAP). NULL (default) leaves BayesianAnalysis's own class default, which is "PosteriorMean" – so by default $parameters is bit-identical to $summary$mean. Pass "PosteriorMode" to report the MAP point instead.
credible_level: width of the posterior credible interval $summary's lower/upper columns report, between 0 and 1. Defaults to 0.9, BayesianAnalysis's own class default. confint() reuses these bounds when its level matches and otherwise re-runs the identical seeded chain, so setting it here is how a bare confint(f) avoids that rebuild.
...: sampler-specific tuning knobs. Passing a knob the chosen sampler does not use is an error: "DEMCz" accepts jump, jump_threshold, noise; "DEMCzs" additionally accepts snooker_threshold; "ARWMH" accepts scale, beta; "NUTS" accepts max_tree_depth.

Value

An object of class corehydro_fit with method == "BayesianAnalysis". $draws is a 3-D array [iteration, chain, parameter], the axis order posterior::as_draws_array() expects, so it can be handed to the posterior or coda packages with no reshaping. $posterior is the thinned draw matrix the analyses consume, $posterior_rows by length(coef(f)), with the parameter names on its columns. $summary is a data frame (one row per parameter, named to match coef()) with columns mean, median, sd, lower, upper, rhat, ess. $map and $posterior_mean are the two posterior point estimates, named to match coef(). $mean_log_likelihood has one entry per iteration, the chain-averaged log-likelihood trace. $acceptance_rates has one entry per chain. $warmup records the warmup actually used and $credible_level the width $summary's lower/upper were computed at. $dic, $waic, $looic are the usual Bayesian goodness-of-fit scalars, $looic_se the standard error of $looic, and $waic_pd/$loo_pd their effective parameter counts. See fit_diagnostics() for leverage/influence diagnostics off a Bayesian fit.

Accepted ranges

BayesianAnalysis rejects a configuration outside these bounds, so fit_bayesian() checks them up front and names the offending value: chains between 4 and 20, iterations between 100 and 1,000,000, warmup between 50 and 100,000, output_length between 100 and 1,000,000, seed not negative, and credible_level strictly between 0 and 1.

Examples

peaks <- c(12500, 15300, 8900, 22100, 18700, 14200, 9800, 28500, 17400, 11600)
f <- fit_bayesian(model_univariate("Normal", peaks),
  sampler = "DEMCz", iterations = 200, output_length = 500, seed = 12345
)
f$summary
#>                  mean    median       sd    lower    upper      rhat      ess
#> Mean (µ)    15530.223 15631.282 2207.107 11923.93 19177.86 0.9998316 284.0929
#> Std Dev (σ)  6641.725  6347.736 1674.511  4278.37  9714.25 1.0267234 326.1511

Usage

Arguments

Value

Accepted ranges

See also

Examples