Skip to contents

Run MCMC inference for a spatial random field model

Source: R/mcmc.R
mcmc.Rd

Performs Markov chain Monte Carlo sampling to estimate the posterior distribution of the spatial field, dependence parameter, and (for hierarchical models) emission parameters. For MDGM models, the DAG structure is also sampled.

Usage

mcmc(
  model,
  y = NULL,
  z_init,
  psi_init,
  theta_init = numeric(0),
  n_iter = 1000L,
  psi_tune = 0.1,
  emission_prior_params = NULL,
  store_z = FALSE,
  seed = NULL,
  nug = NULL
)

Arguments

model

An SrfModel object created by srf_model().

y

Observation data. For hierarchical models, a list of numeric vectors where y[[i]] contains the observations for vertex i. For Bernoulli, values should be 0 or 1; for Poisson, non-negative integers; for Gaussian, any real values. Vertices with no observations should have numeric(0). For standalone models, pass NULL (default).

z_init

Initial color assignment, an integer vector of length n with values in 0, ..., n_colors - 1.

psi_init

Initial value for the dependence parameter (positive).

theta_init

Initial emission parameters. For Bernoulli, a numeric vector of length n_colors (e.g., c(0.2, 0.8)). For Gaussian, a vector of length 2 * n_colors with means then variances (e.g., c(mu_1, mu_2, sigma2_1, sigma2_2)). For Poisson, a vector of length n_colors with rate parameters. Ignored for standalone models.

n_iter

Number of MCMC iterations (default 1000).

psi_tune

Standard deviation of the normal MH proposal for psi (default 0.1).

emission_prior_params

Prior hyperparameters for emission parameters. For Bernoulli: c(a, b) for Beta(a, b) prior. For Gaussian: c(mu_0, sigma2_0, alpha_0, beta_0) where mu_k ~ N(mu_0, sigma2_0) and sigma_k^2 ~ InvGamma(alpha_0, beta_0) independently. For Poisson: c(alpha_0, beta_0) for Gamma(alpha_0, beta_0) prior. Default: c(1, 1) for Bernoulli/Poisson, c(0, 10000, 0.01, 0.01) for Gaussian (non-informative).

store_z

Logical; if TRUE, store the full latent field matrix (n x n_iter). Default is FALSE to conserve memory — on large grids (e.g. 1000x1000 with 10,000 iterations), the z matrix alone requires ~40 GB. When FALSE, per-iteration summary statistics (allocation counts, sufficient statistics, MAP configuration) are still computed.

seed

Optional integer seed for reproducibility.

nug

Optional NaturalUndirectedGraph object. If provided, stored in the result for use by edge_inclusion_probs() and plot().

Value

An SrfResult object containing posterior samples.

See also

srf_model() for model construction, SrfResult for accessing results.

Examples

if (FALSE) { # \dontrun{
# Create a small grid graph
A <- matrix(0, 4, 4)
A[1, 2] <- A[2, 1] <- A[2, 3] <- A[3, 2] <- 1
A[3, 4] <- A[4, 3] <- 1
nug <- nug_from_adj_mat(A, seed = 42L)

# Standalone model
model <- srf_model(nug, spatial = mdgm(dag_type = "spanning_tree"))
result <- mcmc(model, z_init = c(0L, 0L, 1L, 1L),
               psi_init = 0.5, n_iter = 100L)
result$summary()
} # }