Package {BayesNSGP}

Calculate sparse kernel, core kernel, and determine nonzero entries

Description

Cy_sm calculates the normalized sparse kernel for a fixed set of bump function hyperparameters and returns the nonzero entries. Note that the matrix is calculated and returned in dense format.

Usage

Cy_sm(
  dists,
  coords,
  N,
  d,
  n1,
  n2,
  r0,
  s0,
  cstat_opt,
  normalize,
  bumpLocs,
  rads,
  ampls,
  shps,
  dist1_sq,
  dist2_sq,
  dist12,
  Sigma11,
  Sigma22,
  Sigma12,
  nu,
  log_sigma_vec,
  lognuggetSD
)

Arguments

Value

Returns a sparse matrix (N x 3) of the nonzero elements of the product between the core and sparse kernel.

Posterior predictive sampling for NNGP

Description

Generates a posterior predictive draw at a new location x_0 using the Nearest-Neighbor Gaussian Process (NNGP) approximation. This function automatically handles coordinate sorting, neighbor selection, and distance calculations internally following Algorithm 2 of Finley et al. (2019).

Usage

NNGP.pred(x0, coords_sorted, rho, sigma2, w.x0.all, k)

Arguments

Details

The function merges the prediction point x_0 into the training coordinates and identifies the k nearest neighbors from the set of points that precede it in the x-axis ordering.

The predictive distribution is Gaussian:

w(x_0) \mid \mathbf{w}_{N_0} \sim \mathcal{N}(m, v)

where m is the **Kriging mean** and v is the **Kriging variance** defined as:

• m = c^\top C_{N_0}^{-1} \mathbf{w}_{N_0}

• v = C(x_0, x_0) - c^\top C_{N_0}^{-1} c

Here:

• N_0 is the set of indices of the k nearest neighbors of x_0 among the observed locations.

• \mathbf{w}_{N_0} is the vector of realizations of the spatial process at those neighbor locations (passed as w.s.x0).

• c is the k \times 1 covariance vector between x_0 and its neighbors in N_0.

• C_{N_0} is the k \times k covariance matrix among the neighbors in N_0.

The exponential covariance function is used: C(d) = \sigma^2 \exp(-d / \rho). A small nugget (10^{-6}) is added to the diagonal of C_{N_0} for numerical stability.

The internal subsetting of w.x0.all ensures that \mathbf{w}_{N_0} correctly represents the process realizations at the identified neighbor locations.

Value

A numeric scalar representing a posterior predictive draw at location x_0.

Author(s)

Fabian Ketwaroo

References

Finley, A. O., Datta, A., Cook, B. D., Morton, D. C., Andersen, H. E., & Banerjee, S. (2019). Efficient algorithms for Bayesian nearest neighbor Gaussian processes. *Journal of Computational and Graphical Statistics*, 28(2), 401–414. doi:10.1080/10618600.2018.1537924

Examples

# Setup training data
M <- 50
coords <- matrix(runif(M * 2), ncol = 2)
w_all <- rnorm(M) # Mock spatial realizations

# Predict at a new location using k=5 neighbors
new_loc <- c(0.5, 0.5)
pred_draw <- NNGP.pred(x0 = new_loc, 
                       coords_sorted = coords, 
                       rho = 0.2, 
                       sigma2 = 1.0, 
                       w.x0.all = w_all,
                       k = 5)

Helper to extract NNGP local neighborhood nodes for custom MCMC sampling

Description

Pre-calculates the local graph structure, target indices, and specific coordinate strings within the AD matrix for a single target node. This forms the control list payload needed by sampler_RW_NN_GP.

Usage

RWNNGP_setup(node_id, AD, neighbors.id, Rneighbors.id, N.neighbors, k = k)

Arguments

Details

This utility is designed to run in an R loop when configuring an MCMC specification prior to model compilation. By converting matrix indexing operations into explicit character node paths (e.g., "AD[5,1]"), it shifts the burden of matrix searching from the runtime execution loop of the C++ compiled sampler into a one-time R setup cost.

It maps out how a change in node_id will ripple through its forward neighbors and back through the reverse neighbors whose conditional distributions depend directly on the target node's value.

Value

A named list containing specific structural elements for the target node:

• update_id: Integer vector combining the node_id and its reverse neighbors.

• Fneighbors.id: Integer vector containing non-zero forward neighbor indices.

• AFnodes: Character vector of parsed string addresses pointing to forward coefficients.

• ARnodes: Character vector of parsed string addresses pointing to reverse coefficients.

• A.neighbors: Flattened character vector of matrix addresses representing all relevant neighborhood coefficients needed to re-evaluate residuals during the loop.

Author(s)

Fabian Ketwaroo

R_sparse_chol

Description

R_sparse_chol

Usage

R_sparse_chol(i, j, x, n)

Arguments

R_sparse_chol

Description

R_sparse_chol

Usage

R_sparse_cholesky(i, j, x)

Arguments

nimble_sparse_crossprod

Description

nimble_sparse_crossprod

Usage

R_sparse_crossprod(i, j, x, z, n, subset = -1, transp = 1)

Arguments

nimble_sparse_solve

Description

nimble_sparse_solve

Usage

R_sparse_solve(i, j, x, z)

Arguments

nimble_sparse_crossprod

Description

nimble_sparse_crossprod

Usage

R_sparse_solveMat(i, j, x, z, transp = 1)

Arguments

nimble_sparse_tcrossprod

Description

nimble_sparse_tcrossprod

Usage

R_sparse_tcrossprod(i, j, x, subset = -1)

Arguments

Calculate the Gaussian quadratic form for the NNGP approximation

Description

calcQF calculates the quadratic form in the multivariate Gaussian based on the NNGP approximation, for a specific parameter combination. The quadratic form is t(u)C^{-1}v.

Usage

calcQF(u, v, AD, nID)

Arguments

Value

A list with two components: (1) an N x 2 array containing the same spatial coordinates, ordered by MMD, and (2) the same thing, but with any NA values removed.

Calculate A and D matrices for the NNGP approximation

Description

calculateAD_ns calculates A and D matrices (the Cholesky of the precision matrix) needed for the NNGP approximation.

Usage

calculateAD_ns(
  dist1_3d,
  dist2_3d,
  dist12_3d,
  Sigma11,
  Sigma22,
  Sigma12,
  log_sigma_vec,
  log_tau_vec,
  nID,
  N,
  k,
  nu,
  d
)

Arguments

Value

A N x (k+1) matrix; the first k columns are the 'A' matrix, and the last column is the 'D' vector.

Calculate the (sparse) matrix U

Description

calculateU_ns calculates the (sparse) matrix U (i.e., the Cholesky of the inverse covariance matrix) using a nonstationary covariance function. The output only contains non-zero values and is stored as three vectors: (1) the row indices, (2) the column indices, and (3) the non-zero values. NOTE: this code assumes the all inputs correspond to the ORDERED locations.

Usage

calculateU_ns(
  dist1_3d,
  dist2_3d,
  dist12_3d,
  Sigma11,
  Sigma22,
  Sigma12,
  log_sigma_vec,
  log_tau_vec,
  nu,
  nID,
  cond_on_y,
  N,
  k,
  d,
  M = 0
)

Arguments

Value

Returns a sparse matrix representation of the Cholesky of the precision matrix for a fixed set of covariance parameters.

Compute NNGP regression coefficients and conditional variances

Description

Computes the Nearest-Neighbor Gaussian Process (NNGP) factorization components for each location: the regression coefficients (A) and conditional variances (D) following the algorithms described by Finley et al. (2019)

Usage

computeAD(edist, nid.dist, neighbors.id, rho, sigma2, k)

Arguments

Details

This function computes the parameters of the NNGP factorization by expressing the joint distribution as a product of conditional Gaussian densities:

p(\mathbf{x}) = \prod_{i=1}^M p(x_i \mid x_{N(i)})

where N(i) denotes the set of neighbors of location i.

For each location i, the conditional distribution is:

x_i \mid x_{N(i)} \sim \mathcal{N}\left( A_i x_{N(i)}, \; D_i \right)

where:

• A_i is a vector of regression coefficients

• D_i is the conditional variance

These quantities are obtained from the covariance structure:

• A_i = C_{i,N(i)} C_{N(i),N(i)}^{-1}

• D_i = C_{i,i} - C_{i,N(i)} C_{N(i),N(i)}^{-1} C_{N(i),i}

where:

• C_{N(i),N(i)} is the covariance matrix among neighbors

• C_{i,N(i)} is the covariance between location i and its neighbors

The covariance structure is defined using an exponential covariance function.

These parameters are used to efficiently evaluate the NNGP likelihood and simulate from the process.

Value

A numeric matrix (M \times (k+1)) where:

• First k columns: regression coefficients (A matrix)

• Last column: conditional variances (D diagonal entries)

Author(s)

Fabian Ketwaroo

References

Finley, A. O., Datta, A., Cook, B. D., Morton, D. C., Andersen, H. E., & Banerjee, S. (2019). Efficient algorithms for Bayesian nearest neighbor Gaussian processes. *Journal of Computational and Graphical Statistics*, 28(2), 401–414. doi:10.1080/10618600.2018.1537924

Examples

coords <- matrix(runif(40), ncol = 2)
Nk = 5
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)

Construct local covariance matrix for NNGP conditional distributions

Description

Constructs the local covariance matrix associated with a focal location and its neighbor set under the Nearest-Neighbor Gaussian Process (NNGP).

Usage

computeC(fdist, ndist, rho, sigma2)

Arguments

Details

This function constructs the local covariance matrix used in the NNGP approximation for conditional Gaussian distributions. The covariance structure is defined using an exponential covariance function.

The structure of the matrix is:

• Top-left block: covariance among neighbors

• Last row/column: covariance between focal location and neighbors

• Bottom-right: marginal variance of the focal location

A small nugget term (1e-6) is added to the diagonal for numerical stability.

Value

A numeric matrix of dimension (k+1) \times (k+1) representing the local covariance matrix for the focal location and its neighbors.

Author(s)

Fabian Ketwaroo

Compute neighbor sets for a Nearest-Neighbor Gaussian Process (NNGP)

Description

Computes ordered neighbor sets based on Vecchia's approximation for a set of spatial coordinates. Locations are first ordered by their x-coordinate, and neighbors are selected from previously ordered locations.

Usage

computeNeighbors(coords, k, proj)

Arguments

Details

Neighbor selection follows a Vecchia-type ordering. For each location i, neighbors are chosen from the set of previously ordered locations 1, \dots, i-1.

Distances are computed using sf::st_distance and depend on the CRS.

The first location has no neighbors. For the first k+1 locations, all previous locations are used as neighbors. For subsequent locations, the k nearest neighbors among previously ordered points are selected.

Value

A list with components:

coords_sorted

Matrix of coordinates sorted by x-coordinate

edist_sorted

Matrix (M \times M) of pairwise distances between sorted coordinates

neighbors

Binary matrix (M \times M) indicating neighbor relationships

neighbor_idx

Matrix (M \times k) of neighbor indices for each location

neighbors_dist

Matrix (M \times k) of distances to neighbors

Author(s)

Fabian Ketwaroo

See Also

st_distance

Examples

coords <- matrix(runif(40), ncol = 2)
res <- computeNeighbors(coords, k = 5, proj = 4326)

Compute quadratic form for NNGP log-likelihood

Description

Evaluates the quadratic form appearing in the Gaussian log-likelihood under the Nearest-Neighbor Gaussian Process (NNGP).

Usage

computeQF(u, v, AD, neighbors.id)

Arguments

Details

This function computes the quadratic form associated with the NNGP precision matrix without explicitly constructing the full covariance or precision matrix.

The expression is obtained by decomposing the joint Gaussian density into a product of conditional densities:

p(\mathbf{x}) = \prod_{i=1}^M p(x_i \mid x_{N(i)})

where N(i) denotes the set of neighbors of location i.

Each term corresponds to a squared, standardized conditional residual:

r_i = u_i - A_i u_{N(i)}

and the quadratic form is:

\sum_i \frac{r_i^2}{D_i}

This provides an efficient way to evaluate the Gaussian log-likelihood under the NNGP approximation.

Value

A scalar value representing the quadratic form:

\sum_{i=1}^M \frac{\left(u_i - A_i u_{N(i)}\right) \left(v_i - A_i v_{N(i)}\right)}{D_i}

Author(s)

Fabian Ketwaroo

See Also

computeAD

Assign conditioning sets for the SGV approximation

Description

conditionLatentObs assigns q_y(i) vs q_z(i) following Section 5.1 in Katzfuss and Guinness (2018). This function only needs to be run once per SGV analysis.

Usage

conditionLatentObs(nID, coords_ord, N)

Arguments

Value

A matrix indicating whether the conditioning set for each location is on the latent process (y, 1) or the observed values (z, 0).

Calculate sparse kernel, core kernel, and determine nonzero entries

Description

Cy_sm calculates the normalized sparse kernel for a fixed set of bump function hyperparameters and returns the nonzero entries. Note that the matrix is calculated and returned in dense format.

Usage

crossCy_sm(
  Xdists,
  coords,
  Pcoords,
  d,
  n1,
  n2,
  r0,
  s0,
  cstat_opt,
  normalize,
  bumpLocs,
  rads,
  ampls,
  shps,
  Xdist1_sq,
  Xdist2_sq,
  Xdist12,
  Sigma11,
  Sigma22,
  Sigma12,
  PSigma11,
  PSigma22,
  PSigma12,
  nu,
  log_sigma_vec,
  Plog_sigma_vec
)

Arguments

Value

Returns a sparse matrix (N x 3) of the nonzero elements of the product between the core and sparse kernel.

Determine the k-nearest neighbors for each spatial coordinate.

Description

determineNeighbors returns an N x k matrix of the nearest neighbors for spatial locations coords, with the ith row giving indices of the k nearest neighbors to the ith location, which are selected from among the 1,...(i-1) other spatial locations. The first row is -1's, since the first location has no neighbors. The i=2 through i=(k+1) rows each necessarily contain 1:i.

Usage

determineNeighbors(coords, k)

Arguments

Value

An N x k matrix of nearest neighbor indices

Examples

coords <- cbind(runif(100), runif(100))
determineNeighbors(coords, 20)

Plot neighbor sets for a Nearest-Neighbor Gaussian Process (NNGP)

Description

Visualizes the neighbor set for a given location based on an ordered Nearest-Neighbor Gaussian Process (NNGP) constructed using Vecchia's approximation.

Usage

displayNeighbors(ind, coords.sorted, neighbor_idx, proj)

Arguments

Details

The function assumes that coordinates and neighbor indices are based on an ordered representation of the data (e.g., sorted by x-coordinate or another ordering used in Vecchia's approximation).

The neighbor set for location i consists of previously ordered locations, as defined by the NNGP construction.

Value

A ggplot2 object showing:

• All spatial locations (black markers)

• The focal location (red point)

• Its neighbors (blue points)

• A vertical reference line at the focal location's x-coordinate

Author(s)

Fabian Ketwaroo

See Also

computeNeighbors

Examples

coords <- matrix(runif(40), ncol = 2)
res <- computeNeighbors(coords, k = 5, proj = 4326)
p <- displayNeighbors(
  ind = 10,
  coords.sorted = res$coords_sorted,
  neighbor_idx = res$neighbor_idx,
  proj = 4326
)
print(p)

Compute Nearest-Neighbor Gaussian Process (NNGP) log-density

Description

Evaluates the log-density of a multivariate normal distribution under the Nearest-Neighbor Gaussian Process (NNGP) approximation as defined by Data et al 2016.

Usage

dmnorm_NN_GP(x, mu, AD, neighbors.id, log = 1)

Arguments

Details

The NNGP approximation factorizes the joint density as:

p(\mathbf{x}) \approx \prod_{i=1}^M p(x_i \mid x_{N(i)})

leading to the log-density:

-\frac{1}{2} \left( N \log(2\pi) + \sum_i \log D_i + Q \right)

where:

• D_i are conditional variances

• Q is the quadratic form computed by computeQF

The quadratic form corresponds to:

Q = \sum_{i=1}^M \frac{(x_i - \mu_i - A_i(x_{N(i)} - \mu_{N(i)}))^2}{D_i}

This formulation avoids constructing the sparse covariance matrix, enabling scalable likelihood evaluation for large spatial datasets.

Value

A scalar value corresponding to the (log-)density.

Author(s)

Fabian Ketwaroo

References

Datta, A., Banerjee, S., Finley, A.O. and Gelfand, A.E., 2016. Hierarchical nearest-neighbor Gaussian process models for large geostatistical datasets. *Journal of the American Statistical Association*, 111(514), 800-812. doi:10.1080/01621459.2015.1044091

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])
loglike <- dmnorm_NN_GP(w, mu = rep(0, M),
                       AD = AD[1:M, 1:(Nk+1)],
                       neighbors.id = res$neighbor_idx[1:M, 1:Nk],
                       log = TRUE)

Function for the evaluating the Gaussian likelihood with gp2Scale sparse covariance.

Description

dmnorm_gp2Scale (and rmnorm_gp2Scale) calculate the usual Gaussian likelihood for a fixed set of parameters (but with sparse matrices). Finally, the distributions must be registered within nimble.

Usage

dmnorm_gp2Scale(x, mean, Cov, N, Nnz, log = 1)

Arguments

Value

Returns the Gaussian likelihood using the gp2Scale sparse covariance.

Function for the evaluating the NNGP approximate density.

Description

dmnorm_nngp (and rmnorm_nngp) calculate the approximate NNGP likelihood for a fixed set of parameters (i.e., A and D matrices). Finally, the distributions must be registered within nimble.

Usage

dmnorm_nngp(x, mean, AD, nID, N, k, log)

Arguments

Value

The NNGP approximate density.

Function for the evaluating the SGV approximate density.

Description

dmnorm_sgv (and rmnorm_sgv) calculate the approximate SGV likelihood for a fixed set of parameters (i.e., the U matrix). Finally, the distributions must be registered within nimble.

Usage

dmnorm_sgv(x, mean, U, N, k, log = 1)

Arguments

Value

Returns the SGV approximation to the Gaussian likelihood.

Exponential covariance function for Gaussian processes

Description

Computes the exponential covariance matrix.

Usage

expcov(edists, rho, sigma2)

Arguments

Details

The covariance function is defined as:

C_{i,j} = \sigma^2 \exp(-d_{i,j} / \rho)

where d_{ij} is the distance between locations i and j.

Larger values of \rho imply slower decay of correlation with distance (i.e., stronger spatial autocorrelation). Larger values of \sigma^2 increase the marginal variance of the process.

A small nugget term (1e-6) is added to the diagonal for numerical stability.

Value

A numeric matrix (M \times M) representing the covariance between spatial locations.

Author(s)

Fabian Ketwaroo

Examples

coords <- matrix(runif(10), ncol = 2)
d <- as.matrix(dist(coords))
expcov(edists = d, rho = 0.5, sigma2 = 1)

Find reverse neighbors for a single NNGP node

Description

Identifies all "children" (reverse neighbors) for a specific target location j in a Nearest-Neighbor Gaussian Process (NNGP). It searches the neighbor matrix to find all locations k that condition upon location j.

Usage

get_single_reverse_neighbors(target_node, neighbor_idx)

Arguments

Details

This function optimizes the reverse lookup for a single node by avoiding the full O(M) list construction loop. It uses vectorization to scan rows of the neighbor matrix where the row index k is greater than the target index j, enforcing the Directed Acyclic Graph (DAG) structure of the NNGP.

Value

An integer vector containing the indices of all locations k that have target_node as a neighbor, maintaining the NNGP ordering constraint (j < k).

Author(s)

Fabian Ketwaroo

Examples

nn_matrix <- matrix(c(0, 0,
                      1, 0,
                      1, 2,
                      1, 3), nrow = 4, byrow = TRUE)

# Find which nodes have node 1 as a neighbor
get_single_reverse_neighbors(target_node = 1, neighbor_idx = nn_matrix)

Calculate covariance elements based on eigendecomposition components

Description

inverseEigen calculates the inverse eigendecomposition – in other words, the covariance elements based on the eigenvalues and vectors. For a 2x2 anisotropy (covariance) matrix, we parameterize the three unique values in terms of the two log eigenvalues and a rotation parameter on the rescaled logit. The function is coded as a nimbleFunction (see the nimble package) but can also be used as a regular R function.

Usage

inverseEigen(eigen_comp1, eigen_comp2, eigen_comp3, which_Sigma)

Arguments

Value

A vector of anisotropy values (Sigma11, Sigma22, or Sigma12; depends on which_Sigma) for the corresponding set of locations.

Examples

# Generate some eigendecomposition elements (all three are real-valued)
eigen_comp1 <- rnorm(10)
eigen_comp2 <- rnorm(10)
eigen_comp3 <- rnorm(10)
inverseEigen( eigen_comp1, eigen_comp2, eigen_comp3, 2) # Return the Sigma22 values

Calculate a stationary Matern correlation matrix

Description

matern_corr calculates a stationary Matern correlation matrix for a fixed set of locations, based on a range and smoothness parameter. This function is primarily used for the "npGP" and "approxGP" models. The function is coded as a nimbleFunction (see the nimble package) but can also be used as a regular R function.

Usage

matern_corr(dist, rho, nu)

Arguments

Value

A correlation matrix for a fixed set of stations and fixed parameter values.

Examples

# Generate some coordinates
coords <- cbind(runif(100),runif(100))
nu <- 2
# Calculate distances -- can use nsDist to calculate Euclidean distances
dist_list <- nsDist(coords, isotropic = TRUE)
# Calculate the correlation matrix
corMat <- matern_corr(sqrt(dist_list$dist1_sq), 1, nu)

nimble_sparse_chol

Description

nimble_sparse_chol

Usage

nimble_sparse_chol(i, j, x, n)

Arguments

nimble_sparse_chol

Description

nimble_sparse_chol

Usage

nimble_sparse_cholesky(i, j, x)

Arguments

nimble_sparse_crossprod

Description

nimble_sparse_crossprod

Usage

nimble_sparse_crossprod(i, j, x, z, n, subset, transp)

Arguments

nimble_sparse_solve

Description

nimble_sparse_solve

Usage

nimble_sparse_solve(i, j, x, z)

Arguments

nimble_sparse_crossprod

Description

nimble_sparse_crossprod

Usage

nimble_sparse_solveMat(i, j, x, z, transp)

Arguments

nimble_sparse_tcrossprod

Description

nimble_sparse_tcrossprod

Usage

nimble_sparse_tcrossprod(i, j, x, subset)

Arguments

Calculate a nonstationary Matern correlation matrix

Description

nsCorr calculates a nonstationary correlation matrix for a fixed set of locations, based on vectors of the unique anisotropy parameters for each station. Since the correlation function uses a spatially-varying Mahalanobis distance, this function requires coordinate- specific distance matrices (see below). The function is coded as a nimbleFunction (see the nimble package) but can also be used as a regular R function.

Usage

nsCorr(dist1_sq, dist2_sq, dist12, Sigma11, Sigma22, Sigma12, nu, d)

Arguments

Value

A correlation matrix for a fixed set of stations and fixed parameter values.

Examples

# Generate some coordinates and parameters
coords <- cbind(runif(100),runif(100))
Sigma11 <- rep(1, 100) # Identity anisotropy process
Sigma22 <- rep(1, 100)
Sigma12 <- rep(0, 100)
nu <- 2
# Calculate distances
dist_list <- nsDist(coords)
# Calculate the correlation matrix
corMat <- nsCorr(dist_list$dist1_sq, dist_list$dist2_sq, dist_list$dist12,
                 Sigma11, Sigma22, Sigma12, nu, ncol(coords))

Calculate a nonstationary Matern cross-correlation matrix

Description

nsCrosscorr calculates a nonstationary cross-correlation matrix between two fixed sets of locations (a prediction set with M locations, and the observed set with N locations), based on vectors of the unique anisotropy parameters for each station. Since the correlation function uses a spatially-varying Mahalanobis distance, this function requires coordinate- specific distance matrices (see below). The function is coded as a nimbleFunction (see the nimble package) but can also be used as a regular R function.

Usage

nsCrosscorr(
  Xdist1_sq,
  Xdist2_sq,
  Xdist12,
  Sigma11,
  Sigma22,
  Sigma12,
  PSigma11,
  PSigma22,
  PSigma12,
  nu,
  d
)

Arguments

Value

A M x N cross-correlation matrix for two fixed sets of stations and fixed parameter values.

Examples

# Generate some coordinates and parameters
coords <- cbind(runif(100),runif(100))
Sigma11 <- rep(1, 100) # Identity anisotropy process
Sigma22 <- rep(1, 100)
Sigma12 <- rep(0, 100)
Pcoords <- cbind(runif(200),runif(200))
PSigma11 <- rep(1, 200) # Identity anisotropy process
PSigma22 <- rep(1, 200)
PSigma12 <- rep(0, 200)
nu <- 2
# Calculate distances
Xdist_list <- nsCrossdist(coords, Pcoords)
# Calculate the correlation matrix
XcorMat <- nsCrosscorr(Xdist_list$dist1_sq, Xdist_list$dist2_sq, Xdist_list$dist12,
   Sigma11, Sigma22, Sigma12, PSigma11, PSigma22, PSigma12, nu, ncol(coords))

Calculate coordinate-specific cross-distance matrices

Description

nsCrossdist calculates coordinate-specific cross distances in x, y, and x-y for use in the nonstationary cross-correlation calculation. This function is useful for calculating posterior predictions.

Usage

nsCrossdist(coords, Pcoords, scale_factor = NULL, isotropic = FALSE)

Arguments

Value

A list of distances matrices, with the following components:

Examples

# Generate some coordinates
coords <- cbind(runif(100),runif(100))
Pcoords <- cbind(runif(200),runif(200))
# Calculate distances
Xdist_list <- nsCrossdist(coords, Pcoords)

Calculate coordinate-specific cross-distance matrices, only for nearest neighbors and store in an array

Description

nsCrossdist3d generates and returns new 3-dimensional arrays containing the former dist1_sq, dist2_s1, and dist12 matrices, but only as needed for the k nearest-neighbors of each location. these 3D matrices (dist1_3d, dist2_3d, and dist12_3d) are used in the new implementation of calculateAD_ns().

Usage

nsCrossdist3d(
  coords,
  predCoords,
  P_nID,
  scale_factor = NULL,
  isotropic = FALSE
)

Arguments

Value

Arrays with nearest neighbor distances in each coordinate direction. When the spatial dimension d > 2, dist1_3d contains squared Euclidean distances, and dist2_3d and dist12_3d are empty.

Examples

# Generate some coordinates and neighbors
coords <- cbind(runif(100),runif(100))
predCoords <- cbind(runif(200),runif(200))
P_nID <- FNN::get.knnx(coords, predCoords, k = 10)$nn.index # Prediction NN
# Calculate distances
Pdist <- nsCrossdist3d(coords, predCoords, P_nID)

Calculate coordinate-specific distance matrices

Description

nsDist calculates x, y, and x-y distances for use in the nonstationary correlation calculation. The sign of the cross-distance is important. The function contains an optional argument for re-scaling the distances such that the coordinates lie in a square.

Usage

nsDist(coords, scale_factor = NULL, isotropic = FALSE)

Arguments

Value

A list of distances matrices, with the following components:

Examples

# Generate some coordinates
coords <- cbind(runif(100),runif(100))
# Calculate distances
dist_list <- nsDist(coords)
# Use nsDist to calculate Euclidean distances
dist_Euclidean <- sqrt(nsDist(coords, isotropic = TRUE)$dist1_sq)

Calculate coordinate-specific distance matrices, only for nearest neighbors and store in an array

Description

nsDist3d generates and returns new 3-dimensional arrays containing the former dist1_sq, dist2_sq, and dist12 matrices, but only as needed for the k nearest-neighbors of each location. these 3D matrices (dist1_3d, dist2_3d, and dist12_3d) are used in the new implementation of calculateAD_ns().

Usage

nsDist3d(coords, nID, scale_factor = NULL, isotropic = FALSE)

Arguments

Value

Arrays with nearest neighbor distances in each coordinate direction.

Examples

# Generate some coordinates and neighbors
coords <- cbind(runif(100),runif(100))
nID <- determineNeighbors(coords, 10)
# Calculate distances
nsDist3d(coords, nID)

NIMBLE code for a generic nonstationary GP model

Description

This function sets up and compiles a nimble model for a general nonstationary Gaussian process.

Usage

nsgpModel(
  tau_model = "constant",
  sigma_model = "constant",
  Sigma_model = "constant",
  sparse_model = "none",
  mu_model = "constant",
  likelihood = "fullGP",
  coords,
  data,
  constants = list(),
  monitorAllSampledNodes = TRUE,
  ...
)

Arguments

Value

A nimbleCode object.

Examples

# Generate some data: stationary/isotropic
N <- 100
coords <- matrix(runif(2*N), ncol = 2)
alpha_vec <- rep(log(sqrt(1)), N) # Log process SD
delta_vec <- rep(log(sqrt(0.05)), N) # Log nugget SD
Sigma11_vec <- rep(0.4, N) # Kernel matrix element 1,1
Sigma22_vec <- rep(0.4, N) # Kernel matrix element 2,2
Sigma12_vec <- rep(0, N) # Kernel matrix element 1,2
mu_vec <- rep(0, N) # Mean
nu <- 0.5 # Smoothness
dist_list <- nsDist(coords)
Cor_mat <- nsCorr( dist1_sq = dist_list$dist1_sq, dist2_sq = dist_list$dist2_sq,
                   dist12 = dist_list$dist12, Sigma11 = Sigma11_vec,
                   Sigma22 = Sigma22_vec, Sigma12 = Sigma12_vec, nu = nu )
Cov_mat <- diag(exp(alpha_vec)) %*% Cor_mat %*% diag(exp(alpha_vec))
D_mat <- diag(exp(delta_vec)^2)
set.seed(110)
data <- as.numeric(mu_vec + t(chol(Cov_mat + D_mat)) %*% rnorm(N))
# Set up constants
constants <- list( nu = 0.5, Sigma_HP1 = 2 )
# Defaults: tau_model = "constant", sigma_model = "constant", mu_model = "constant",
# and Sigma_model = "constant"
Rmodel <- nsgpModel(likelihood = "fullGP", constants = constants, coords = coords, data = data )

Posterior prediction for the NSGP

Description

nsgpPredict conducts posterior prediction for MCMC samples generated using nimble and nsgpModel.

Usage

nsgpPredict(
  model,
  samples,
  coords.predict,
  predict.process = TRUE,
  constants = list(),
  seed = 0,
  ...
)

Arguments

Value

The output of the function is a list with two elements: obs, a matrix of J posterior predictive samples for the N observed locations (only for likelihood = "SGV", which produces predictions for the observed locations by default; this element is NULL otherwise); and pred, a corresponding matrix of posterior predictive samples for the prediction locations. Ordering and neighbor selection for the prediction coordinates in the SGV likelihood are conducted internally, as with nsgpModel.

Examples

## Not run: 
# Generate some data: stationary/isotropic
N <- 100
coords <- matrix(runif(2*N), ncol = 2)
alpha_vec <- rep(log(sqrt(1)), N) # Log process SD
delta_vec <- rep(log(sqrt(0.05)), N) # Log nugget SD
Sigma11_vec <- rep(0.4, N) # Kernel matrix element 1,1
Sigma22_vec <- rep(0.4, N) # Kernel matrix element 2,2
Sigma12_vec <- rep(0, N) # Kernel matrix element 1,2
mu_vec <- rep(0, N) # Mean
nu <- 0.5 # Smoothness
dist_list <- nsDist(coords)
Cor_mat <- nsCorr( dist1_sq = dist_list$dist1_sq, dist2_sq = dist_list$dist2_sq,
                   dist12 = dist_list$dist12, Sigma11 = Sigma11_vec,
                   Sigma22 = Sigma22_vec, Sigma12 = Sigma12_vec, nu = nu )
Cov_mat <- diag(exp(alpha_vec)) %*% Cor_mat %*% diag(exp(alpha_vec))
D_mat <- diag(exp(delta_vec)^2)
set.seed(110)
data <- as.numeric(mu_vec + t(chol(Cov_mat + D_mat)) %*% rnorm(N))
# Set up constants
constants <- list( nu = 0.5, Sigma_HP1 = 2 )
# Defaults: tau_model = "constant", sigma_model = "constant", mu_model = "constant",
# and Sigma_model = "constant"
Rmodel <- nsgpModel(likelihood = "fullGP", constants = constants, coords = coords, data = data )
conf <- configureMCMC(Rmodel)
Rmcmc <- buildMCMC(conf)
Cmodel <- compileNimble(Rmodel)
Cmcmc <- compileNimble(Rmcmc, project = Rmodel)
samples <- runMCMC(Cmcmc, niter = 200, nburnin = 100)
# Prediction
predCoords <- as.matrix(expand.grid(seq(0,1,l=10),seq(0,1,l=10)))
postpred <- nsgpPredict( model = Rmodel, samples = samples, coords.predict = predCoords )

## End(Not run)

Order coordinates according to a maximum-minimum distance criterion.

Description

orderCoordinatesMMD orders an array of (x,y) spatial coordinates according to the "maximum minimum distance" (MMD), as described in Guinness, 2018. (Points are selected to maximize their minimum distance to already- selected points).

Usage

orderCoordinatesMMD(coords, exact = FALSE)

Arguments

Value

A list of distances matrices, with the following components:

Examples

coords <- cbind(runif(100), runif(100))
orderCoordinatesMMD(coords)

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

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. 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])

Function for the evaluating the SGV approximate density.

Description

dmnorm_gp2Scale (and rmnorm_gp2Scale) calculate the usual Gaussian likelihood for a fixed set of parameters (but with sparse matrices). Finally, the distributions must be registered within nimble.

Usage

rmnorm_gp2Scale(n, mean, Cov, N, Nnz)

Arguments

Value

Not applicable.

Function for the evaluating the NNGP approximate density.

Description

dmnorm_nngp (and rmnorm_nngp) calculate the approximate NNGP likelihood for a fixed set of parameters (i.e., A and D matrices). Finally, the distributions must be registered within nimble.

Usage

rmnorm_nngp(n, mean, AD, nID, N, k)

Arguments

Value

The NNGP approximate density.

Function for the evaluating the SGV approximate density.

Description

dmnorm_sgv (and rmnorm_sgv) calculate the approximate SGV likelihood for a fixed set of parameters (i.e., the U matrix). Finally, the distributions must be registered within nimble.

Usage

rmnorm_sgv(n, mean, U, N, k)

Arguments

Value

Not applicable.

NNGP-Specific Random Walk Metropolis-Hastings Sampler

Description

A custom NIMBLE sampler for updating individual spatial random effects under a Nearest-Neighbor Gaussian Process (NNGP) approximation. This sampler uses a factorized likelihood approach to perform highly efficient local Metropolis-Hastings sampling with a normal proposal distribution (Metropolis, 1953), implementing the adaptation routine given in Shaby and Wells (2011).

Usage

sampler_RW_NN_GP(model, mvSaved, target, control)

Arguments

Details

Instead of calculating the full NNGP log-likelihood, which is O(M) (where M is the total number of spatial locations), this sampler exploits the Directed Acyclic Graph (DAG) structure of the NNGP to reduce the complexity of a single-node update to O(k^2), where k is the number of neighbors.

During the MCMC configuration stage, the sampler automatically queries the local graph architecture using get_single_reverse_neighbors and RWNNGP_setup to resolve dependencies internally. When a single node w_i is updated, only its own conditional density and the conditional densities of its "children" (the reverse neighbors that depend on it) are affected.

The local log-Metropolis-Hastings ratio (\log MHR) is calculated as:

\log MHR = -\frac{(r_i^{*2} - r_i^2)}{2D_i} + \sum_{j \in \mathcal{R}(i)} -\frac{(r_j^{*2} - r_j^2)}{2D_j}

where:

• r_i and r_i^* are the current and proposed conditional residuals for the target node.

• \mathcal{R}(i) is the set of reverse neighbors (indices j such that i \in N(j)).

• r_j and r_j^* are the current and proposed residuals for those reverse neighbors.

• D represents the conditional variances.

The residuals are updated efficiently using the difference \delta = w_i^* - w_i, such that r_j^* = r_j - A_{ji}\delta, where A_{ji} is the NNGP regression coefficient.

Author(s)

Fabian Ketwaroo

References

Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., and Teller, E. (1953). Equation of state calculations by fast computing machines. *The Journal of Chemical Physics*, 21(6), 1087-1092.

Shaby, B. A. and Wells, M. T. (2011). Exploring an adaptive Metropolis-Hastings algorithm. *Department of Statistical Science, Cornell University Tech Report*.

One-time setup wrapper function for the SGV approximation

Description

sgvSetup is a wrapper function that sets up the SGV approximation. Three objects are required: (1) ordering the locations, (2) identify nearest neighbors, and (3) determine the conditioning set. This function only needs to be run once per SGV analysis.

Usage

sgvSetup(
  coords,
  coords_pred = NULL,
  k = 15,
  seed = NULL,
  pred.seed = NULL,
  order_coords = TRUE,
  order_coords_pred = TRUE
)

Arguments

Value

A list with the following components: