| Title: | GPU-Portable Vendor-Agnostic Metropolis-Hastings Sampler |
|---|---|
| Description: | A generic random-walk Metropolis-Hastings sampler for Markov chain Monte Carlo. The user declares a model by writing its log-likelihood and log-prior as one-sided formulas; the package compiles them to a stack-machine bytecode that a single portable kernel interprets, so the same model runs on the CPU and on the GPU from one source. The sampler advances many independent chains in one batched pass. A distributional equivalence harness based on the two-sample Kolmogorov-Smirnov test and the split R-hat statistic of Gelman and Rubin (1992) <doi:10.1214/ss/1177011136> checks the GPU output against the CPU output. Dispatching the kernel to a GPU helps when many chains are run or when the log-density is expensive to evaluate over a large data set, and does not help for small models run with few chains. |
| Authors: | Pedro Carvalho Brom [aut, cre, cph] (ORCID: <https://orcid.org/0000-0002-1288-7695>) |
| Maintainer: | Pedro Carvalho Brom <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.7.0 |
| Built: | 2026-11-10 11:25:21 UTC |
| Source: | https://github.com/pcbrom/gpumetropolis |
Estimates the effective sample size of a set of chains: the number of independent draws that would carry the same information as the autocorrelated Markov chain Monte Carlo output. Each chain is reduced by Geyer's initial positive sequence estimator and the per-chain values are summed.
ess(x, warmup = NULL)ess(x, warmup = NULL)
x |
A |
warmup |
Number of initial iterations discarded before the estimate.
When |
A single numeric value, the summed effective sample size.
Geyer, C. J. (1992). Practical Markov chain Monte Carlo. Statistical Science 7(4), 473-483. doi:10.1214/ss/1177011137.
set.seed(1) x <- rnorm(800, mean = 0, sd = 1) fit <- metropolis_gaussian_mean(x, sigma = 1, n_iter = 2000) ess(fit)set.seed(1) x <- rnorm(800, mean = 0, sd = 1) fit <- metropolis_gaussian_mean(x, sigma = 1, n_iter = 2000) ess(fit)
Returns the analytic posterior of the mean mu for observations following
Normal(mu, sigma^2) with known sigma and a flat prior on mu. The
posterior is Normal(mean(data), sigma^2 / length(data)). It provides the
ground truth used to check that metropolis_gaussian_mean() recovers known
parameters.
gaussian_mean_posterior(data, sigma)gaussian_mean_posterior(data, sigma)
data |
Numeric vector of observations. |
sigma |
Positive finite scalar, the known standard deviation. |
A list with the posterior mean and standard deviation sd.
set.seed(1) gaussian_mean_posterior(rnorm(500, mean = 3, sd = 2), sigma = 2)set.seed(1) gaussian_mean_posterior(rnorm(500, mean = 3, sd = 2), sigma = 2)
Advances many independent random-walk Metropolis chains over a model
declared with gpum_model(). The log-density kernel runs on the chosen
backend; the data is uploaded once and stays resident across the run.
gpu_metropolis( model, data = NULL, init = NULL, proposal_sd = 0.1, n_iter = 2000L, n_chains = 4L, warmup = NULL, adapt = TRUE, method = c("rwm", "pt", "de", "mala"), temperatures = NULL, swap_every = NULL, gamma = NULL, de_noise = 0.001, de_every = NULL, de_sync = FALSE, seed = 1L, backend = c("auto", "cpu", "cuda", "vulkan") )gpu_metropolis( model, data = NULL, init = NULL, proposal_sd = 0.1, n_iter = 2000L, n_chains = 4L, warmup = NULL, adapt = TRUE, method = c("rwm", "pt", "de", "mala"), temperatures = NULL, swap_every = NULL, gamma = NULL, de_noise = 0.001, de_every = NULL, de_sync = FALSE, seed = 1L, backend = c("auto", "cpu", "cuda", "vulkan") )
model |
A |
data |
A named list or data frame with one entry per data column named in the model. Ignored for a model with no data term. |
init |
Optional numeric matrix of starting values, |
proposal_sd |
Initial scale of the Gaussian random-walk proposal.
Accepts a scalar (recycled to all parameters and chains), a length-
|
n_iter |
Iterations the sampler runs per chain, counting warmup. Default 2000. |
n_chains |
Number of chains. Used only when |
warmup |
Warmup iterations to discard before returning, following the
convention of Stan and nimble. Must lie in |
adapt |
Whether to adapt the per-chain proposal scale during
warmup. When |
method |
Sampling method: |
temperatures |
Numeric vector of length |
swap_every |
Integer, the iteration count between swap proposals
in |
gamma |
Scale of the difference vector in |
de_noise |
Non-negative factor for the per-dimension jitter added
to the |
de_every |
Integer, the iteration count between refreshes of the
frozen population snapshot used for the difference vectors in
|
de_sync |
Logical. When |
seed |
Integer seed. Each chain advances its own counter-based stream from a triple32 hash; the seed is itself hashed, so runs with consecutive integer seeds get independent streams. |
backend |
Compute backend: |
An object of class gpum_fit: a list with draws (an
n_iter - warmup by n_chains by n_params array of post-warmup
samples), accept_rate, n_iter (kept count), n_iter_total (raw
count), warmup and the rest of the run metadata.
set.seed(1) y <- rnorm(2000, mean = 3, sd = 2) m <- gpum_model(~ -((y - mu)^2) / 8, params = "mu", data = "y") fit <- gpu_metropolis(m, data = list(y = y), proposal_sd = 0.05, n_iter = 1000, n_chains = 4) rhat(fit$draws[, , 1])set.seed(1) y <- rnorm(2000, mean = 3, sd = 2) m <- gpum_model(~ -((y - mu)^2) / 8, params = "mu", data = "y") fit <- gpu_metropolis(m, data = list(y = y), proposal_sd = 0.05, n_iter = 1000, n_chains = 4) rhat(fit$draws[, , 1])
Computes the Bayes factor B_10 = p(y | model1) / p(y | model0) from the
log marginal likelihoods of the two models, each by gpum_evidence(). The
Bayes factor is the formal Bayesian weight of evidence for one model over
another; both models must carry a proper prior.
gpum_bayes_factor(model1, model0, data, ...)gpum_bayes_factor(model1, model0, data, ...)
model1, model0
|
Two |
data |
The data, shared by both models. |
... |
Passed to |
The Bayes factor is sensitive to the priors, sharply so for a point null
under a diffuse prior, the Jeffreys-Lindley effect: a more diffuse prior on
the larger model penalises it regardless of the data. Report the priors with
the factor, and prefer gpum_waic() or gpum_loo() when the goal is
predictive comparison rather than a weight of evidence.
An object of class gpum_bayes_factor with bf10, log_bf10, the
two log evidences and a verbal interpretation on the Jeffreys scale.
gpum_evidence(), gpum_waic(), gpum_loo()
Infers the dependence structure of two variables by fitting a bivariate
copula to their pseudo-observations (average ranks rescaled to the open
unit square), so the fit is invariant to the marginal distributions and
isolates the dependence. The copula parameter is sampled with
gpu_metropolis(); the copula log-density is a normalised density on the
unit square, so WAIC and PSIS-LOO compare families directly.
gpum_copula( data, family = "auto", criterion = c("waic", "loo"), method = "mala", warmup = "auto", n_iter = 8000L, n_chains = 4L, seed = 1L, ... )gpum_copula( data, family = "auto", criterion = c("waic", "loo"), method = "mala", warmup = "auto", n_iter = 8000L, n_chains = 4L, seed = 1L, ... )
data |
A two-column data frame or matrix; the columns are the two variables whose dependence is modelled. |
family |
One of |
criterion |
Selection criterion for |
method |
Sampler passed to |
warmup |
Warmup passed to |
n_iter, n_chains, seed, ...
|
Passed to |
With family = "auto" every candidate family is fit and the workflow
returns the one preferred by predictive comparison together with the full
ranking, since the choice of family is a hypothesis about tail behaviour.
Gumbel carries upper-tail dependence, Clayton lower-tail, Frank symmetric
and tail-independent, Gaussian elliptical and tail-independent; a
selected Gumbel is itself the statement that co-movement concentrates in
the upper tail.
The three Archimedean families (Gumbel, Clayton, Frank) are parametrised for positive association in this release; the Gaussian family carries the sign of the dependence through its correlation. For data with negative rank correlation, the Gaussian family applies directly and the others do not.
An object of class gpum_copula: a list with the selected
family, the gpum_fit, the posterior summaries of the dependence
parameter and of Kendall's tau and the tail-dependence coefficients,
the pseudo-observations, and (for "auto") the comparison table.
gpu_metropolis(), kendall_tau(), tail_dependence()
set.seed(1) u <- runif(300) v <- pnorm(0.8 * qnorm(u) + 0.6 * rnorm(300)) fit <- gpum_copula(data.frame(x = qexp(u), y = qgamma(v, 2)), family = "auto", n_iter = 4000, n_chains = 4) fitset.seed(1) u <- runif(300) v <- pnorm(0.8 * qnorm(u) + 0.6 * rnorm(300)) fit <- gpum_copula(data.frame(x = qexp(u), y = qgamma(v, 2)), family = "auto", n_iter = 4000, n_chains = 4) fit
Computes the observed Fisher information of the model log-likelihood at the
posterior mean by a central-difference Hessian, inverts it to the
Cramer-Rao lower bound on the covariance of an unbiased estimator, and
returns it beside the empirical posterior covariance for comparison. This
is an optional reference, not the headline of the fit: the bound is a
frequentist asymptotic object and is meaningful only where the model is
regular. The function applies four guards and reports
applicable = FALSE with a reason when any fails.
gpum_crlb(fit, data = NULL, at = NULL, rhat_max = 1.05)gpum_crlb(fit, data = NULL, at = NULL, rhat_max = 1.05)
fit |
A |
data |
The same data passed to |
at |
Optional numeric vector, the point at which to evaluate the information. Defaults to the posterior mean (the cold chain for a parallel-tempering fit, the pooled chains otherwise). |
rhat_max |
Convergence guard: when the largest split R-hat over the parameters exceeds this, the comparison is flagged not applicable, since the asymptotic-normal reference assumes a unimodal converged posterior. Default 1.05. |
Under the Bernstein-von Mises theorem, for a regular model with a flat or weak prior the posterior covariance approaches the inverse Fisher information; an agreement is then a check that the sampler recovered the information-bound geometry. The identification is asymptotic and assumes an interior true parameter, identifiability, a finite non-singular information and a smooth correctly specified likelihood. It does not transfer to multimodal posteriors or to parameters on a boundary, and an informative prior legitimately makes the posterior tighter than the bound.
The comparison rests on estimator theory: under the same regularity the posterior mean is a consistent and asymptotically efficient estimator, that is it concentrates on the true parameter and attains the Cramer-Rao bound as the data grow, so the bound is the right yardstick for the posterior spread. Consistency is the weaker and more robust property and can hold where asymptotic efficiency does not; this function reports the bound only where the stronger normal-approximation conditions also hold.
An object of class gpum_crlb: a list with information (the
observed Fisher information matrix), crlb (its inverse, the bound on the
covariance), crlb_sd, crlb_cor, the empirical posterior_cov,
posterior_sd and posterior_cor, the evaluation point at, a logical
applicable and a character note. When not applicable the matrix
entries are NA and note states why.
The reference is not restricted to i.i.d. data. For independent,
non-identically distributed rows (regression-type models, every row
with its own covariates) the observed information is the sum of the
per-row contributions and the asymptotic normality that backs the
comparison holds under the Lindeberg-Feller conditions and local
asymptotic normality (Le Cam; van der Vaart 1998, chapter 7). For
dependent data declared through the conditional factorisation of
gpum_ts_model(), the per-row terms are conditional log-densities,
their score is a martingale difference sequence, and the summed
conditional curvature is again the right information: the martingale
central limit theorem (Hall and Heyde 1980) replaces Lindeberg-Feller,
with the Bernstein-von Mises theorem for ergodic Markov processes
(Borwanker, Kallianpur and Prakasa Rao 1971) backing the posterior
side. The stationarity and ergodicity of the process are model
assumptions the guards below cannot check.
gpu_metropolis(), gpum_diagnose(), gpum_ts_model()
Overlays the kernel density of the observed data with the kernel density of one or more generated samples on a single set of axes, the visual check that a fitted model reproduces the data it was fit to. For competing models the better fit, especially in the relevant tail, is read off directly.
gpum_density_compare( observed, generated, labels = NULL, main = "Observed against generated densities", xlab = "value", colors = NULL, ... )gpum_density_compare( observed, generated, labels = NULL, main = "Observed against generated densities", xlab = "value", colors = NULL, ... )
observed |
Numeric vector of the observed data. |
generated |
A numeric vector, or a named list of numeric vectors, of
generated or posterior-predictive samples (see |
labels |
Optional character labels for the generated samples; taken
from the names of |
main, xlab
|
Plot title and x-axis label. |
colors |
Optional colours for the generated densities. |
... |
Passed to |
Invisibly, a list with the observed and generated density objects.
Produces a per-parameter table (mean, standard deviation, 2.5%, 50%
and 97.5% quantiles, split R-hat, effective sample size and Monte
Carlo standard error) and a convergence verdict from the asymptotic
canonical thresholds (R-hat below 1.01 in every parameter and ESS at
or above 400). When plot = TRUE, opens a multi-panel plot per
parameter showing the trace, the pooled density, the running mean
per chain and the pooled autocorrelation; when the fit carries
adaptation book-keeping, an extra panel shows the per-chain
acceptance rate by warmup batch with the asymptotic optimum drawn as
a horizontal reference.
gpum_diagnose(fit, plot = TRUE, return_data = FALSE, crlb = NULL)gpum_diagnose(fit, plot = TRUE, return_data = FALSE, crlb = NULL)
fit |
A |
plot |
Whether to open the diagnostic plot. Default |
return_data |
Whether to return the structured stats invisibly.
Default |
crlb |
Optional |
When return_data = TRUE, a list with summary (one row
per parameter), verdict (the convergence diagnosis and a
suggested next step), adaptation (passed through from
fit$adaptation) and adaptation_hint (a character string when
the last warmup batch closed below 80% of the asymptotic target
acceptance, NULL otherwise). Otherwise NULL invisibly.
gpu_metropolis(), rhat(), ess()
Estimates the log marginal likelihood (the model evidence) by
thermodynamic integration along a path of power posteriors
p_beta(theta) proportional to p(y | theta)^beta p(theta) from the prior
(beta = 0) to the posterior (beta = 1). The log evidence is the integral
over beta of the expected log-likelihood under p_beta (Gelman and Meng
1998; Friel and Pettitt 2008), evaluated on a power-law grid that is dense
near the prior, where the integrand is most curved. Each rung is sampled by
the package's own sampler, so no kernel change is needed.
gpum_evidence( model, data, n_rungs = 24L, power = 5, n_iter = 4000L, n_chains = 8L, warmup = NULL, method = "de", seed = 1L, backend = "cpu" )gpum_evidence( model, data, n_rungs = 24L, power = 5, n_iter = 4000L, n_chains = 8L, warmup = NULL, method = "de", seed = 1L, backend = "cpu" )
model |
A |
data |
The data for the model. |
n_rungs |
Number of intervals in the |
power |
Exponent of the power-law grid |
n_iter, n_chains, warmup, method, seed, backend
|
Passed to
|
The evidence requires a proper prior: the marginal likelihood is an integral against the prior, undefined for a flat or absent prior. The function stops when the model has no prior term. The estimate is sensitive to the prior, a property of the marginal likelihood itself, not of the sampler.
An object of class gpum_evidence with log_evidence, the betas
grid and the expected_loglik per rung.
gpum_bayes_factor(), gpum_waic()
Detects the support and the modality of a numeric column, fits every
eligible parametric family with the sampler, and ranks them by predictive
comparison. The winner's fitted cumulative distribution is the
probability-integral transform used to feed gpum_copula(), so the
explicit marginals and the copula together model the full joint
distribution (Sklar's theorem).
gpum_fit_catalog( data, catalog = "auto", ranking = c("waic", "loo"), n_iter = 8000L, n_chains = 4L, seed = 1L, backend = "cpu", ... )gpum_fit_catalog( data, catalog = "auto", ranking = c("waic", "loo"), n_iter = 8000L, n_chains = 4L, seed = 1L, backend = "cpu", ... )
data |
A numeric vector, or a one-column data frame or matrix. |
catalog |
Either |
ranking |
Selection criterion, |
n_iter, n_chains, seed, backend, ...
|
Passed to |
Every family is fit in an unconstrained parametrisation, so the sampler
never proposes an invalid parameter; the natural parameters are recovered
for reporting. A column detected multimodal brings the two-component
Gaussian mixture into the eligible set, which is fit with parallel
tempering (method = "pt") and relabeled so the lower-location component
is first, since a random walk cannot cross the modes.
An object of class gpum_catalog: a list with the best family
name and its fit, the comparison ranking table, the support and
modality verdict, and the fitted CDF of the winner.
gpu_metropolis(), gpum_copula(), marginal_cdf()
set.seed(1) y <- rgamma(400, shape = 3, rate = 2) result <- gpum_fit_catalog(y, n_iter = 4000, n_chains = 4, backend = "cpu") resultset.seed(1) y <- rgamma(400, shape = 3, rate = 2) result <- gpum_fit_catalog(y, n_iter = 4000, n_chains = 4, backend = "cpu") result
Reports the posterior probability that a parameter lies in an interval, below it and above it, directly from the draws. This is an estimation decision: it needs no marginal likelihood and inherits the robustness of the posterior the sampler returned.
gpum_hypothesis(fit, parameter, lower = -Inf, upper = Inf)gpum_hypothesis(fit, parameter, lower = -Inf, upper = Inf)
fit |
A |
parameter |
Name of the parameter. |
lower, upper
|
Bounds of the hypothesis interval. Defaults span the
real line, so |
An object of class gpum_hypothesis with prob (mass in the
interval), prob_below, prob_above and the inputs.
Estimates the expected log predictive density of one-step-ahead
forecasts by leave-future-out cross-validation with an expanding
window: for each held-out time t in the evaluation range, the model
is refit on rows 1:(t - 1) and scored on the conditional log
predictive density of row t, integrated over the posterior by
averaging the pointwise likelihood across the draws. This is the
evaluation scheme that respects temporal ordering; gpum_waic() and
gpum_loo() assume exchangeable pointwise contributions and are not
valid for dependent rows (Burkner, Gabry and Vehtari 2020). The refits
are exact rather than importance-sampling approximations: the sampler
is fast enough that the approximation of Burkner, Gabry and Vehtari
(2020) buys nothing here.
gpum_lfo(ts_model, L = NULL, n_iter = 4000L, n_chains = 4L, seed = 1L, ...)gpum_lfo(ts_model, L = NULL, n_iter = 4000L, n_chains = 4L, seed = 1L, ...)
ts_model |
A |
L |
Minimum training rows before the first evaluation: rows
|
n_iter, n_chains, seed, ...
|
Passed to |
A list of class gpum_lfo with elpd_lfo (the summed
expected log predictive density), se (its standard error over the
evaluated points), and pointwise (one row per evaluated time).
gpum_ts_model(), gpum_waic(), gpum_loo()
Builds the conjugate specification of the Gaussian linear model
y = X beta + eps, eps ~ N(0, sigma^2 I), under the
Normal-inverse-Gamma prior
beta | sigma^2 ~ N(b0, sigma^2 B0^{-1}),
sigma^2 ~ InvGamma(c0 / 2, d0 / 2),
with B0 a prior precision matrix. Under this prior the joint posterior
of (beta, sigma^2) is again Normal-inverse-Gamma in closed form, so
gpu_metropolis() samples it exactly and independently, with no warmup,
no proposal and no rejection; the effective sample size equals the number
of draws. The marginal likelihood is also available in closed form and is
returned by the fit as log_marginal.
gpum_lm(formula, data, b0 = 0, B0 = 1e-06, c0 = 0.002, d0 = 0.002)gpum_lm(formula, data, b0 = 0, B0 = 1e-06, c0 = 0.002, d0 = 0.002)
formula |
A standard R model formula, e.g. |
data |
A data frame holding the variables of |
b0 |
Prior mean of |
B0 |
Prior precision of |
c0 |
Prior shape times two of |
d0 |
Prior rate times two of |
The prior must be proper (B0 positive definite, c0 > 0, d0 > 0) for
the marginal likelihood to exist; the sampling itself only needs the
posterior precision X'X + B0 to be positive definite. The defaults are
weakly informative and proper.
An object of class gpum_conjugate holding the design matrix,
the response and the prior, ready for gpu_metropolis().
gpu_metropolis(), gpum_model()
m <- gpum_lm(mpg ~ wt, data = mtcars) fit <- gpu_metropolis(m, n_iter = 4000, n_chains = 4, seed = 1) fit$log_marginalm <- gpum_lm(mpg ~ wt, data = mtcars) fit <- gpu_metropolis(m, n_iter = 4000, n_chains = 4, seed = 1) fit$log_marginal
Computes PSIS-LOO cross-validation by delegating to the loo package, which
implements the Pareto-smoothed importance-sampling estimator and its
diagnostics (Vehtari, Gelman and Gabry 2017). Like WAIC it estimates
out-of-sample predictive accuracy without a marginal likelihood, and the
Pareto k diagnostic flags observations where the estimate is unreliable.
As with gpum_waic(), the leave-one-out question is only meaningful for
independent rows; for the dependent rows of a gpum_ts_model() use
gpum_lfo().
gpum_loo(fit, data, max_draws = 4000L)gpum_loo(fit, data, max_draws = 4000L)
fit |
A |
data |
The same data passed to |
max_draws |
Cap on the number of posterior draws used. Default 4000. |
The loo object returned by loo::loo().
Compiles a log-likelihood and an optional log-prior, declared as one-sided
formulas, into the bytecode the sampler runs. The log-likelihood is a
per-observation expression: the sampler sums it over the data. The formulas
may use +, -, *, /, ^, unary -, and exp, log, sqrt; any
other symbol or function is rejected with a clear error.
gpum_model(loglik, params, data = character(0), prior = NULL)gpum_model(loglik, params, data = character(0), prior = NULL)
loglik |
One-sided formula, the per-observation log-likelihood, up to an additive constant. It may reference the parameter names and the data column names. |
params |
Character vector of parameter names. |
data |
Character vector of data column names. Empty for a model with no data term. |
prior |
One-sided formula, the joint log-prior over the parameters, up
to an additive constant. It may reference only the parameter names. |
An object of class gpum_model.
# Gaussian mean with known sd = 2 and a flat prior on mu. m <- gpum_model( loglik = ~ -((y - mu)^2) / 8, params = "mu", data = "y" )# Gaussian mean with known sd = 2 and a flat prior on mu. m <- gpum_model( loglik = ~ -((y - mu)^2) / 8, params = "mu", data = "y" )
A matrix of panels over the parameters: the marginal posterior density with
its highest density interval on the diagonal, and the bivariate posterior
with its credible-region contour off the diagonal. When a gpum_crlb object
is supplied and applicable, the Cramer-Rao ellipse is overlaid on the
off-diagonal panels, so the recovered convergence region can be read against
the information-bound reference.
gpum_pairs(fit, crlb = NULL, level = 0.95, max_points = 3000L, n = 60L, ...)gpum_pairs(fit, crlb = NULL, level = 0.95, max_points = 3000L, n = 60L, ...)
fit |
A |
crlb |
Optional |
level |
Mass of the highest density interval and credible region. Default 0.95. |
max_points |
Cap on plotted draws per panel. Default 3000. |
n |
Grid resolution of the bivariate density. Default 60. |
... |
Unused. |
The fit, invisibly.
gpum_region(), gpum_surface(), gpum_crlb()
Draws n parameter vectors from the posterior (the cold chain under
parallel tempering) and, for each, calls generate() to simulate from the
model, returning the pooled posterior-predictive sample. The package does
not yet generate from an arbitrary likelihood on its own; the user supplies
the family's one-line simulator, which they already know from the model they
declared. An automatic generator from the model itself is the goal of the
synthesis release on the roadmap.
gpum_ppc(fit, generate, n = 4000L)gpum_ppc(fit, generate, n = 4000L)
fit |
A |
generate |
A function of one named parameter vector (one posterior draw) returning one or more simulated observations. It may close over the data to simulate a covariate-dependent response. |
n |
Number of posterior draws to simulate from. Default 4000. |
A numeric vector, the pooled posterior-predictive sample.
Returns the highest posterior density region of two parameters at level
level, as contour polygons of the bivariate posterior density. The region
is the set the parameters concentrate in, the convergence region of the
pair. The result composes onto any plot with lines().
gpum_region(fit, params, level = 0.95, n = 80L)gpum_region(fit, params, level = 0.95, n = 80L)
fit |
A |
params |
Length-two character vector of parameter names. |
level |
Probability mass of the region. Default 0.95. |
n |
Grid resolution of the density estimate. Default 80. |
An object of class gpum_region with the contours (a list of
x, y polygons), the params and the level. Plot it directly or add
it to an existing plot with lines().
Applies the region-of-practical-equivalence rule (Kruschke 2018): the parameter is practically equivalent to the null when its highest density interval lies inside the ROPE, practically different when the interval lies entirely outside it, and undecided otherwise. The decision is from the draws and needs no marginal likelihood.
gpum_rope(fit, parameter, rope, null = 0, ci = 0.95)gpum_rope(fit, parameter, rope, null = 0, ci = 0.95)
fit |
A |
parameter |
Name of the parameter. |
rope |
Either the half-width of a symmetric ROPE around |
null |
Centre of a symmetric ROPE. Ignored when |
ci |
Mass of the highest density interval. Default 0.95. |
An object of class gpum_rope with the hdi, the rope interval,
the share pct_in_rope of the posterior inside it, and a decision.
Plots the bivariate posterior density of two parameters as level curves in two dimensions and as a three-dimensional surface with the credible-region contours projected on the floor. The level curves are the highest posterior density regions at the requested masses.
gpum_surface( fit, params, type = c("both", "contour", "persp"), levels = c(0.5, 0.8, 0.95), n = 80L, theta = 30, phi = 25, ... )gpum_surface( fit, params, type = c("both", "contour", "persp"), levels = c(0.5, 0.8, 0.95), n = 80L, theta = 30, phi = 25, ... )
fit |
A |
params |
Length-two character vector of parameter names. |
type |
One of |
levels |
Probability masses of the level curves. Default
|
n |
Grid resolution. Default 80. |
theta, phi
|
Viewing angles of the 3D surface. |
... |
Passed to |
The fit, invisibly.
Builds a gpum_model() for dependent data with a Markov dependence of
order p. The user writes the conditional log-density of one
observation given its p predecessors, referring to the current value
by the series name and to the lagged values as <name>_lag1 through
<name>_lagp; gpum_ts_model() assembles the lagged design so the
engine's row sum computes the conditional log-likelihood
sum_{t = p + 1}^{n} log f(y_t | y_{t-1}, ..., y_{t-p}).
gpum_ts_model( loglik, params, series, order = 1L, covariates = NULL, prior = NULL )gpum_ts_model( loglik, params, series, order = 1L, covariates = NULL, prior = NULL )
loglik |
One-sided formula, the conditional log-density of the
current observation. Refer to the series by its name in |
params |
Character vector of parameter names. |
series |
Named list with one entry, the observed series
( |
order |
Markov order |
covariates |
Optional named list of exogenous covariate vectors,
each the same length as the series; row |
prior |
Optional one-sided formula for the log-prior (plus, when wanted, the stationary log-density of the initial state), evaluated once per parameter vector. |
The likelihood conditions on the first p observations, the standard
conditional-likelihood treatment (Billingsley 1961; Hamilton 1994,
section 5.2); for n moderately large the information in the first p
points is negligible. A stationary initial density, when available in
closed form, can be added through prior, which is evaluated once per
parameter vector.
Inference downstream is licensed by the dependent-data asymptotics: the
conditional score is a martingale difference sequence, so the observed
information of the conditional factorisation is the right curvature
(martingale CLT, Hall and Heyde 1980) and the posterior satisfies a
Bernstein-von Mises theorem for ergodic Markov processes (Borwanker,
Kallianpur and Prakasa Rao 1971). Model comparison by gpum_waic() or
gpum_loo() assumes exchangeable pointwise terms and does not transfer
to dependent rows; use gpum_lfo() instead.
A list of class gpum_ts with the compiled model (a
gpum_model), the assembled data list ready for
gpu_metropolis(), the order, and the effective sample rows
n_rows = length(y) - order.
gpum_model(), gpu_metropolis(), gpum_lfo()
set.seed(1) y <- as.numeric(arima.sim(list(ar = 0.6), n = 300)) ts_m <- gpum_ts_model( ~ -0.5 * ((y - mu - phi * (y_lag1 - mu)) / exp(ls))^2 - ls, params = c("mu", "phi", "ls"), series = list(y = y), order = 1, prior = ~ -0.5 * (mu / 10)^2 - 0.5 * (phi / 2)^2 - 0.5 * (ls / 2)^2 ) fit <- gpu_metropolis(ts_m$model, data = ts_m$data, n_iter = 2000, n_chains = 4, seed = 1)set.seed(1) y <- as.numeric(arima.sim(list(ar = 0.6), n = 300)) ts_m <- gpum_ts_model( ~ -0.5 * ((y - mu - phi * (y_lag1 - mu)) / exp(ls))^2 - ls, params = c("mu", "phi", "ls"), series = list(y = y), order = 1, prior = ~ -0.5 * (mu / 10)^2 - 0.5 * (phi / 2)^2 - 0.5 * (ls / 2)^2 ) fit <- gpu_metropolis(ts_m$model, data = ts_m$data, n_iter = 2000, n_chains = 4, seed = 1)
Computes WAIC from the pointwise log-likelihood of the fit (Watanabe 2010; Vehtari, Gelman and Gabry 2017). WAIC estimates out-of-sample predictive accuracy without a marginal likelihood, so it is robust where Bayes factors are prior-sensitive; lower WAIC is preferred when comparing models on the same data.
gpum_waic(fit, data, max_draws = 4000L)gpum_waic(fit, data, max_draws = 4000L)
fit |
A |
data |
The same data passed to |
max_draws |
Cap on the number of posterior draws used, to bound the pointwise matrix. Default 4000. |
The comparison is valid only when the log-likelihood is normalised. The
model DSL evaluates the log-density up to an additive constant, which is
enough for sampling a fixed model but not for comparing models with
different likelihoods, since the dropped constants differ. To compare
models with WAIC or LOO, write the full normalised per-observation
log-likelihood in the formula, including the normalising constant such as
-0.5 * log(2 * pi * sigma^2) for a Gaussian.
The pointwise decomposition also assumes the rows are exchangeable
contributions, which holds for independent data (identically distributed
or not) but fails for the dependent rows of a gpum_ts_model()
conditional factorisation: leaving out a past observation that later
rows condition on is not a predictive question. For dependent data use
gpum_lfo(), which respects the temporal ordering.
An object of class gpum_waic with waic, elpd_waic (the
expected log pointwise predictive density), p_waic (the effective number
of parameters), the standard error se and the pointwise terms.
gpum_loo(), gpum_bayes_factor()
The shortest interval that contains a probability mass ci of the
samples, found by scanning every window of the sorted draws.
hdi(x, ci = 0.95)hdi(x, ci = 0.95)
x |
Numeric samples. |
ci |
Probability mass in the interval. Default 0.95. |
A named numeric vector c(lower, upper).
Returns the posterior summary of Kendall's tau, the rank-correlation
functional of the copula parameter, on its natural [-1, 1] scale. Tau
is comparable across families, unlike the family-specific dependence
parameter.
kendall_tau(object)kendall_tau(object)
object |
A |
A named numeric vector: mean, q2.5, median, q97.5.
gpum_copula(), tail_dependence()
Pools the post-warmup draws of two samplers and applies the two-sample Kolmogorov-Smirnov test. It is the equivalence gate for the package: a CPU and a GPU sampler are treated as equivalent when the test does not reject the hypothesis that their draws come from the same distribution.
ks_equivalence(x, y, warmup = NULL, alpha = 0.05, thin = TRUE)ks_equivalence(x, y, warmup = NULL, alpha = 0.05, thin = TRUE)
x, y
|
Each a |
warmup |
Number of initial iterations discarded from |
alpha |
Significance level of the test. Default 0.05. |
thin |
Controls thinning before the test. |
The Kolmogorov-Smirnov test assumes independent draws, while Markov chain
Monte Carlo output is autocorrelated; feeding it raw draws makes the test
reject far too often. By default the pooled draws are therefore thinned down
to the effective sample size returned by ess() before the test is applied.
Not rejecting is evidence consistent with equivalence, not a proof that the
two distributions are identical. Report the statistic and the p-value, and
read equivalent as "no detected difference at level alpha".
A list with the test statistic, the p_value, the chosen alpha,
a logical equivalent (TRUE when p_value > alpha) and the post-thinning
pooled sample sizes n_x and n_y.
set.seed(1) x <- rnorm(800, mean = 1, sd = 1) a <- metropolis_gaussian_mean(x, sigma = 1, n_iter = 2000, seed = 1) b <- metropolis_gaussian_mean(x, sigma = 1, n_iter = 2000, seed = 2) ks_equivalence(a, b)set.seed(1) x <- rnorm(800, mean = 1, sd = 1) a <- metropolis_gaussian_mean(x, sigma = 1, n_iter = 2000, seed = 1) b <- metropolis_gaussian_mean(x, sigma = 1, n_iter = 2000, seed = 2) ks_equivalence(a, b)
Returns the probability-integral transform of the winning family at the
posterior-mean parameters, either as a function or evaluated at the data.
This is the transform that carries a column into the unit interval for
gpum_copula(), the explicit-marginal counterpart of the rank transform.
marginal_cdf(object, at = NULL)marginal_cdf(object, at = NULL)
object |
A |
at |
Optional numeric vector at which to evaluate the CDF; when omitted the function itself is returned. |
A function, or the CDF values at at.
gpum_fit_catalog(), gpum_copula()
Samples the posterior of the mean mu of observations assumed to follow
Normal(mu, sigma^2) with known sigma and a flat prior on mu. Many
independent chains are advanced together: at every iteration the candidate
state of each chain is evaluated by a single batched call to the log-density
kernel.
metropolis_gaussian_mean( data, sigma, n_iter = 2000L, n_chains = 4L, init = NULL, proposal_sd = NULL, seed = 1L )metropolis_gaussian_mean( data, sigma, n_iter = 2000L, n_chains = 4L, init = NULL, proposal_sd = NULL, seed = 1L )
data |
Numeric vector of observations. Must be finite and non-empty. |
sigma |
Positive finite scalar, the known standard deviation of the observations. |
n_iter |
Number of iterations recorded per chain. Default 2000. |
n_chains |
Number of independent chains. Used only when |
init |
Optional numeric vector of starting values, one per chain. When
supplied, its length sets the number of chains. When |
proposal_sd |
Positive finite scalar, the standard deviation of the
Gaussian random-walk proposal. When |
seed |
Single integer seed. Each chain runs an independent PCG64 stream
derived from |
This is the CPU reference sampler. Its log-density kernel runs on the CPU.
Later package versions dispatch the same kernel to a GPU and are validated
for distributional equivalence against this function with ks_equivalence()
and rhat(). The sequential dependence within a chain (x_{t+1} depends on
x_t) is intrinsic to Markov chain Monte Carlo and is not parallelised; the
parallel axes are the chains and the sum over the data.
An object of class gpumetropolis_fit: a list with draws (an
n_iter by n_chains numeric matrix), accept_rate (per-chain
acceptance rate) and the run metadata.
rhat(), ess(), ks_equivalence(), gaussian_mean_posterior()
set.seed(1) x <- rnorm(500, mean = 3, sd = 2) fit <- metropolis_gaussian_mean(x, sigma = 2, n_iter = 1000) rhat(fit)set.seed(1) x <- rnorm(500, mean = 3, sd = 2) fit <- metropolis_gaussian_mean(x, sigma = 2, n_iter = 1000) rhat(fit)
Draws a three-panel diagnostic for the winning family: the fitted density over the data histogram (observed against generated), the QQ plot of the empirical against the fitted quantiles, and the WAIC ranking of the candidate families. Visual companion of the printed text diagnostic.
## S3 method for class 'gpum_catalog' plot(x, ...)## S3 method for class 'gpum_catalog' plot(x, ...)
x |
A |
... |
Ignored. |
x, invisibly.
Draws a four-panel diagnostic: the pseudo-observations with the fitted copula density contours over them, the posterior of the dependence parameter, the posterior of Kendall's tau with the empirical value marked, and the observed-against-generated check on the unit square. This is the visual companion of the text diagnostic printed by the object.
## S3 method for class 'gpum_copula' plot(x, n_grid = 40L, ...)## S3 method for class 'gpum_copula' plot(x, n_grid = 40L, ...)
x |
A |
n_grid |
Grid resolution for the density contours. Default 40. |
... |
Ignored. |
x, invisibly. Called for the plot.
gpum_copula(), gpum_density_compare()
Computes the split potential scale reduction factor for a set of chains. Each chain is split in half and the halves are treated as separate chains, which makes the statistic sensitive to non-stationarity within a chain. A value near 1 is consistent with the chains having reached the same distribution; values above roughly 1.01 to 1.1 indicate that they have not.
rhat(x, warmup = NULL)rhat(x, warmup = NULL)
x |
A |
warmup |
Number of initial iterations discarded before the split. When
|
A single numeric value, the split R-hat statistic.
Gelman, A. and Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences. Statistical Science 7(4), 457-472. doi:10.1214/ss/1177011136.
set.seed(1) x <- rnorm(500, mean = 0, sd = 1) fit <- metropolis_gaussian_mean(x, sigma = 1, n_iter = 1000) rhat(fit)set.seed(1) x <- rnorm(500, mean = 0, sd = 1) fit <- metropolis_gaussian_mean(x, sigma = 1, n_iter = 1000) rhat(fit)
Returns the posterior summaries of the lower and upper tail-dependence coefficients, the probabilities of joint extremes in each tail. Gumbel has upper-tail dependence and no lower, Clayton the reverse, Frank and Gaussian neither.
tail_dependence(object)tail_dependence(object)
object |
A |
A list with lower and upper, each a named numeric vector.