rmnorm_NN_GP: Simulate from a Nearest-Neighbor Gaussian Process (NNGP)

In BayesNSGP: Bayesian Analysis of Non-Stationary Gaussian Process Models

Simulate from a Nearest-Neighbor Gaussian Process (NNGP)

Description

Efficiently simulates a realization from a Gaussian process using the Nearest-Neighbor Gaussian Process (NNGP) approximation as described by Datta (2022).

Usage

rmnorm_NN_GP(n, mu, AD, neighbors.id)

Arguments

n	Integer. Number of samples to generate. Currently only supports n = 1.
mu	Numeric vector (length M). Mean vector of the process.
AD	Numeric matrix (M \times (k+1)). Output from computeAD, where: Columns 1:k contain regression coefficients (A matrix) Column k+1 contains conditional variances (D)
neighbors.id	Integer matrix (M \times k). Neighbor indices for each location, defining the conditioning set N(i). Output from computeNeighbors.

Details

The simulation follows the recursive form derived from the sparse Cholesky factorization of the precision matrix (\Sigma^{-1}). Let z = x - \mu be the zero-mean spatial residuals. The NNGP simulates these residuals sequentially:

z_i = \sum_{j \in N(i)} A_{ij} z_j + \epsilon_i

where:

• N(i) is the neighbor set of location i, such that j < i.

• A_{ij} are the kriging weights computed from the local covariance.

• \epsilon_i \sim \mathcal{N}(0, D_i) are independent innovations.

This sequential approach exploits the sparsity of the Cholesky factor L, where \Sigma^{-1} \approx (I-A)^\top D^{-1} (I-A) and allows for O(Mk^3) simulation, which is significantly faster than the standard O(M^3) Cholesky decomposition for large M.

Value

A numeric vector (length M) representing a simulated realization from the NNGP.

Author(s)

Fabian Ketwaroo

References

Datta, A. (2022). Nearest-neighbor sparse Cholesky matrices in spatial statistics. *Wiley Interdisciplinary Reviews: Computational Statistics*, 14(5), e1574. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1002/wics.1574")}

Examples

M <- 2000      # number of spatial locations
coords <- matrix(runif(2*M), ncol = 2)
Nk <- 15       # number of neighbors
res <- computeNeighbors(coords, k = Nk)
AD <- computeAD(edist = res$edist_sorted,
                nid.dist = res$neighbors_dist,
                neighbors.id = res$neighbor_idx,
                rho = 0.1,
                sigma2 = 0.3,
                k = Nk)
w <- rmnorm_NN_GP(n = 1, mu = rep(0, M),
                 AD = AD[1:M, 1:(Nk+1)],
                 neighbors.id = res$neighbor_idx[1:M, 1:Nk])

BayesNSGP documentation built on Sept. 10, 2026, 5:08 p.m.