tulpa_nuts_spde: Sample an SPDE GLM via NUTS, optionally jointly over Matern...

In tulpa: Templated Unified Library for Posterior Approximation in Bayesian Hierarchical Models

Sample an SPDE GLM via NUTS, optionally jointly over Matern hypers

Description

Bayesian counterpart to fit_spde() / laplace_spde_at(). Routes the SPDE-augmented GLM through tulpa's full NUTS backend via the generic LikelihoodSpec interface. Two modes:

• Legacy / fixed-hyper (joint = FALSE, the default): conditions on the supplied ⁠(range, sigma)⁠. Samples the latent block ⁠(beta, w_mesh, log_phi)⁠ jointly. Useful as the inner step of an outer CCD / nested-Laplace integration over hypers.

• Joint hypers (joint = TRUE): samples ⁠(beta, z_mesh, log_kappa, log_tau, log_phi)⁠ jointly via the non-centered transform ⁠w = L(theta)^{-T} z⁠. A PC prior (Fuglstad et al. 2019 JASA) on ⁠(range, sigma)⁠ is taken from prior_range, prior_sigma. Returns additional w_draws, range_draws, sigma_draws, kappa_draws, tau_draws columns/vectors on top of the raw parameter draws.

Usage

tulpa_nuts_spde(
  y,
  X,
  spatial,
  family = c("gaussian", "poisson", "binomial", "gamma", "neg_binomial_2", "beta"),
  n_trials = NULL,
  joint = FALSE,
  range = NULL,
  sigma = NULL,
  prior_range = NULL,
  prior_sigma = NULL,
  log_kappa_init = NULL,
  log_tau_init = NULL,
  beta_prior = .tulpa_default_beta_prior("spde_nuts"),
  log_phi_prior_sd = 3,
  log_phi_init = 0,
  control = list()
)

Arguments

y	Response vector. Family-specific: gaussian: any real poisson, neg_binomial_2: non-negative integers binomial: non-negative integers in ⁠[0, n_trials]⁠ gamma: strictly positive reals beta: strictly in ⁠(0, 1)⁠
X	Fixed-effects design matrix.
spatial	A tulpa_spatial object from spatial_spde() / spatial_spde_custom() – supplies the FEM matrices (C0, G1) and projection (A) plus the smoothness nu.
family	One of "gaussian", "poisson", "binomial", "gamma", "neg_binomial_2", "beta".
n_trials	Integer vector for family = "binomial" (else ignored).
joint	Logical. FALSE (default) conditions on fixed ⁠(range, sigma)⁠. TRUE activates joint sampling of ⁠(log_kappa, log_tau)⁠ with the PC prior from prior_range, prior_sigma.
range, sigma	Fixed-hyper mode only. Matern range and marginal SD on the field. Default to the SPDE prior medians (spatial$prior_range[1], spatial$prior_sigma[1]).
prior_range, prior_sigma	Joint mode only. PC prior anchors as c(value, alpha) pairs: prior_range = c(r0, a_r) encodes P(range < r0) = a_r prior_sigma = c(s0, a_s) encodes P(sigma > s0) = a_s Default to the spec's anchors (spatial$prior_range, spatial$prior_sigma) when those are length-2 PC anchors.
log_kappa_init, log_tau_init	Joint mode only. Initial values for the hyper slots. Default to the value implied by the PC anchor's ⁠(r0, s0)⁠ pair via ⁠kappa = sqrt(8 nu) / r0⁠, ⁠tau = 1 / (sqrt(4 pi) * kappa * s0)⁠.
beta_prior	Fixed-effect prior as list(mean, sd): a mean-zero (mean = 0) Gaussian on each coefficient with SD sd (default the engine default, prior_normal(0, 2.5)).
log_phi_prior_sd	Prior SD on log(phi). Role of phi is family-specific: gaussian: phi is the residual SD (sampled jointly) gamma: phi is the Gamma shape (sampled jointly) neg_binomial_2: phi is the NB size r (sampled jointly) beta: phi is the Beta precision (sampled jointly) poisson, binomial: log_phi is held tight and ignored downstream.
log_phi_init	Starting value for log(phi).
control	A named list of numerical / sampler knobs (statistical arguments stay in the signature): n_iter (default 2000), n_warmup (default 1000), max_treedepth (default 10), adapt_delta (default 0.8), seed (NULL draws from the session RNG), verbose (default FALSE), noncenter (fixed-hyper only, default TRUE: sample the mesh field in a non-centered ⁠v = L^{-T} z⁠ parameterisation – the same target density, a priori isotropised; ignored when joint = TRUE), and mass_matrix (NUTS metric: "auto" (default) picks a dense metric at ⁠<= 200⁠ parameters and diagonal above, capturing the fixed-effect / field cross-curvature that corrects the intercept marginal; "diag", "dense", "block_diag" force the choice).

Value

A list with draws (matrix ⁠n_samples x n_params⁠), means, phi_summary (where applicable), accept_prob, divergent, treedepth, epsilon, joint_hypers, plus the supplied spatial spec. In joint = TRUE mode additionally: w_draws (transformed z -> w per draw), range_draws, sigma_draws, kappa_draws, tau_draws, and range_summary, sigma_summary (mean/median/5%/95% quantiles).

See Also

fit_spde() for the Laplace counterpart and the nested-Laplace path over (range, sigma).

Examples

# Gaussian response with a smooth spatial trend, sampled at fixed
# (range, sigma) on a coarse mesh.
set.seed(1)
n <- 60L
coords <- cbind(runif(n), runif(n))
x <- runif(n, -1, 1)
y <- 0.5 - 0.8 * x + 0.4 * sin(3 * coords[, 1]) + rnorm(n, 0, 0.3)
spec <- spatial_spde(coords, max_edge = c(0.3, 0.6), cutoff = 0.1,
                     prior_range = c(0.5, 0.5), prior_sigma = c(1, 0.5))
fit <- tulpa_nuts_spde(y, cbind(1, x), spatial = spec, family = "gaussian",
                       range = 0.4, sigma = 0.5, log_phi_init = log(0.3),
                       control = list(n_iter = 300L, n_warmup = 150L,
                                      seed = 1L))
colMeans(fit$draws[, c("beta[1]", "beta[2]")])
fit$phi_summary

tulpa documentation built on Oct. 19, 2026, 5:09 p.m.