simulate.metric_graph: Simulate a Whittle-Matérn field on a metric graph

In MetricGraph: Random Fields on Metric Graphs

Simulate a Whittle-Matérn field on a metric graph

Description

Draws unconditional (prior) samples of a Whittle-Matérn field on a metric graph using one of two exact algorithms from Section 6.6 of the paper.

Usage

## S3 method for class 'metric_graph'
simulate(
  object,
  nsim = 1,
  seed = NULL,
  alpha = 1,
  method = c("direct", "kriging", "extended"),
  impl = c("cpp", "R"),
  kappa,
  tau,
  range,
  sigma,
  PtE = NULL,
  type = "manual",
  BC = 1,
  parallel = FALSE,
  n_cores = NULL,
  cluster = NULL,
  ...
)

Arguments

object	A metric_graph object.
nsim	Number of samples. Returns a vector (nsim = 1) or matrix (nsim > 1) of sampled values.
seed	Integer seed for set.seed() before the vertex draw. Pass NULL (default) to use the current RNG state.
alpha	Smoothness parameter (1 or 2).
method	Simulation algorithm: "direct" (Method A, O(m^3)), "kriging" (Method B, O(m)), or "extended" (add simulation locations as graph vertices, single sparse Cholesky draw).
impl	Implementation: "cpp" (default, Rcpp/Eigen) or "R" (pure R reference, for testing/comparison). Ignored for method = "extended".
kappa	Range parameter.
tau	Precision parameter.
range	Practical correlation range (alternative to kappa/tau).
sigma	Marginal standard deviation (alternative to kappa/tau).
PtE	Matrix with columns ⁠(edge_number, normalised_position)⁠. Required when type = "manual" (the default).
type	Location specification: "manual" (use PtE), "mesh" (mesh nodes), or "obs" (observation locations).
BC	Boundary condition for degree-1 vertices: 1 = stationary (default), 0 = Neumann.
parallel	Logical. Use a parallel cluster for the edge loop?
n_cores	Number of parallel workers. Ignored if cluster is supplied. Defaults to detectCores() - 1.
cluster	An already-registered doParallel cluster. If NULL and parallel = TRUE, a cluster is created and stopped automatically.
...	Ignored.

Details

Algorithm A ("direct"): interior values on each edge are drawn from the bridge conditional \mathcal{N}(S_e b_e, \Sigma^*_e) via a dense Cholesky of the m \times m bridge covariance. Cost: O(m^3) per edge. Preferred when m is small (roughly m \lesssim 100).

Algorithm B ("kriging"): an unconditioned Matérn process is simulated on each edge by a linear Markov recursion, then the bridge boundary data are imposed by a kriging correction (Theorem 2 / Proposition 4). Cost: O(m) per edge. Preferred for large m.

Both algorithms produce draws from exactly the same finite-dimensional distribution; they differ only in computational cost.

The two-step procedure is:

1. Draw the vertex (boundary) state from the vertex precision via a sparse Cholesky.

2. For each edge (independently given the vertex state), draw the interior values using the chosen algorithm.

When parallel = TRUE the edge loop is distributed over a doParallel/foreach cluster. Reproducibility requires passing a seed so that per-edge seeds are derived deterministically.

Value

Numeric vector (nsim = 1) or matrix with nsim columns of field values at the requested locations, ordered to match PtE.

See Also

sample_spde()

Examples

# Small square graph, alpha = 1
V <- rbind(c(0,0), c(1,0), c(1,1), c(0,1))
E <- rbind(c(1,2), c(2,3), c(3,4), c(4,1))
g <- metric_graph$new(V = V, E = E)
t_norm <- seq(0.1, 0.9, by = 0.2)
PtE    <- do.call(rbind, lapply(1:4, function(e) cbind(e, t_norm)))
u <- simulate(g, alpha = 1, method = "direct", kappa = 1, tau = 1, PtE = PtE)

# Method B on a larger graph, parallel
u2 <- simulate(g, alpha = 2, method = "kriging", range = 1, sigma = 1,
               PtE = PtE, parallel = TRUE, n_cores = 2)

MetricGraph documentation built on Nov. 11, 2026, 1:09 a.m.