| Title: | Bayesian Analysis of Non-Stationary Gaussian Process Models |
|---|---|
| Description: | Enables off-the-shelf functionality for fully Bayesian, nonstationary Gaussian process modeling. The approach to nonstationary modeling involves a closed-form, convolution-based covariance function with spatially-varying parameters; these parameter processes can be specified either deterministically (using covariates or basis functions) or stochastically (using approximate Gaussian processes). Stationary Gaussian processes are a special case of our methodology, and we furthermore implement approximate Gaussian process inference to account for very large spatial data sets (Finley, et al (2017) <doi:10.48550/arXiv.1702.00434>). Bayesian inference is carried out using Markov chain Monte Carlo methods via the "nimble" package, and posterior prediction for the Gaussian process at unobserved locations is provided as a post-processing step. |
| Authors: | Daniel Turek [aut, cre], Mark Risser [aut] |
| Maintainer: | Daniel Turek <[email protected]> |
| License: | GPL-3 |
| Version: | 0.2.0 |
| Built: | 2026-06-09 07:23:39 UTC |
| Source: | https://github.com/cran/BayesNSGP |
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.
calcQF(u, v, AD, nID)calcQF(u, v, AD, nID)
u |
Vector; left product. |
v |
Vector; right product |
AD |
N x (k+1) matrix; the first k columns are the 'A' matrix, and the
last column is the 'D' vector. Represents the Cholesky of |
nID |
N x k matrix of neighbor indices. |
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.
calculateAD_ns calculates A and D matrices (the Cholesky of the
precision matrix) needed for the NNGP approximation.
calculateAD_ns( dist1_3d, dist2_3d, dist12_3d, Sigma11, Sigma22, Sigma12, log_sigma_vec, log_tau_vec, nID, N, k, nu, d )calculateAD_ns( dist1_3d, dist2_3d, dist12_3d, Sigma11, Sigma22, Sigma12, log_sigma_vec, log_tau_vec, nID, N, k, nu, d )
dist1_3d |
N x (k+1) x (k+1) array of distances in the x-coordinate direction. |
dist2_3d |
N x (k+1) x (k+1) array of distances in the y-coordinate direction. |
dist12_3d |
N x (k+1) x (k+1) array of cross-distances. |
Sigma11 |
N-vector; 1-1 element of the Sigma() process. |
Sigma22 |
N-vector; 2-2 element of the Sigma() process. |
Sigma12 |
N-vector; 1-2 element of the Sigma() process. |
log_sigma_vec |
N-vector; process standard deviation values. |
log_tau_vec |
N-vector; nugget standard deviation values. |
nID |
N x k matrix of neighbor indices. |
N |
Scalar; number of data measurements. |
k |
Scalar; number of nearest neighbors. |
nu |
Scalar; Matern smoothness parameter. |
d |
Scalar; dimension of the spatial domain. |
A N x (k+1) matrix; the first k columns are the 'A' matrix, and the last column is the 'D' vector.
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.
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 )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 )
dist1_3d |
N x (k+1) x (k+1) array of distances in the x-coordinate direction. |
dist2_3d |
N x (k+1) x (k+1) array of distances in the y-coordinate direction. |
dist12_3d |
N x (k+1) x (k+1) array of cross-distances. |
Sigma11 |
N-vector; 1-1 element of the Sigma() process. |
Sigma22 |
N-vector; 2-2 element of the Sigma() process. |
Sigma12 |
N-vector; 1-2 element of the Sigma() process. |
log_sigma_vec |
N-vector; process standard deviation values. |
log_tau_vec |
N-vector; nugget standard deviation values. |
nu |
Scalar; Matern smoothness parameter. |
nID |
N x k matrix of (ordered) neighbor indices. |
cond_on_y |
A matrix indicating whether the conditioning set for each
(ordered) location is on the latent process (y, |
N |
Scalar; number of data measurements. |
k |
Scalar; number of nearest neighbors. |
d |
Scalar; dimension of the spatial domain. |
M |
Scalar; number of prediction sites. |
Returns a sparse matrix representation of the Cholesky of the precision matrix for a fixed set of covariance parameters.
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.
conditionLatentObs(nID, coords_ord, N)conditionLatentObs(nID, coords_ord, N)
nID |
N x k matrix of neighbor indices. |
coords_ord |
N x 2 matrix of locations. |
N |
Scalar; number of locations (observed only!). |
A matrix indicating whether the conditioning set for each location
is on the latent process (y, 1) or the observed values (z, 0).
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.
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 )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 )
Xdists |
N x N matrix of Euclidean distances |
coords |
N x d matrix of coordinate/input locations |
Pcoords |
N x d matrix of coordinate/input locations |
d |
Scalar; dimension of the spatial domain. |
n1 |
Scalar; number of outer products. |
n2 |
Scalar; number of bump functions in each outer product. |
r0 |
Scalar; length-scale of sparse stationary kernel. |
s0 |
Scalar; signal-variance of sparse stationary kernel. |
cstat_opt |
Scalar; determines the compactly supported kernel. See Details. |
normalize |
Logical; should C_sparse have 1's along the diagonal (1 = TRUE) |
bumpLocs |
Array of bump function locations (n2*d x n1) |
rads |
Matrix of bump function radii (n1 x n2; denoted |
ampls |
Matrix of bump function amplitudes (n1 x n2; denoted |
shps |
Matrix of bump function shape parameters (n1 x n2; denoted |
Xdist1_sq |
N x N matrix; contains values of pairwise squared distances in the x-coordinate. |
Xdist2_sq |
N x N matrix; contains values of pairwise squared distances in the y-coordinate. |
Xdist12 |
N x N matrix; contains values of pairwise signed cross-
distances between the x- and y-coordinates. The sign of each element is
important; see |
Sigma11 |
Vector of length N; contains the 1-1 element of the anisotropy process for each station. |
Sigma22 |
Vector of length N; contains the 2-2 element of the anisotropy process for each station. |
Sigma12 |
Vector of length N; contains the 1-2 element of the anisotropy process for each station. |
PSigma11 |
Vector of length N; contains the 1-1 element of the anisotropy process for each station. |
PSigma22 |
Vector of length N; contains the 2-2 element of the anisotropy process for each station. |
PSigma12 |
Vector of length N; contains the 1-2 element of the anisotropy process for each station. |
nu |
Scalar; Matern smoothness parameter. |
log_sigma_vec |
Vector of length N; log of the signal standard deviation. |
Plog_sigma_vec |
Vector of length N; log of the signal standard deviation. |
Returns a sparse matrix (N x 3) of the nonzero elements of the product between the core and sparse kernel.
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.
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 )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 )
dists |
N x N matrix of Euclidean distances |
coords |
N x d matrix of coordinate/input locations |
N |
Scalar; number of data measurements. |
d |
Scalar; dimension of the spatial domain. |
n1 |
Scalar; number of outer products. |
n2 |
Scalar; number of bump functions in each outer product. |
r0 |
Scalar; length-scale of sparse stationary kernel. |
s0 |
Scalar; signal-variance of sparse stationary kernel. |
cstat_opt |
Scalar; determines the compactly supported kernel. See Details. |
normalize |
Logical; should C_sparse have 1's along the diagonal |
bumpLocs |
Array of bump function locations (n2*d x n1) |
rads |
Matrix of bump function radii (n1 x n2; denoted |
ampls |
Matrix of bump function amplitudes (n1 x n2; denoted |
shps |
Matrix of bump function shape parameters (n1 x n2; denoted |
dist1_sq |
N x N matrix; contains values of pairwise squared distances in the x-coordinate. |
dist2_sq |
N x N matrix; contains values of pairwise squared distances in the y-coordinate. |
dist12 |
N x N matrix; contains values of pairwise signed cross-
distances between the x- and y-coordinates. The sign of each element is
important; see |
Sigma11 |
Vector of length N; contains the 1-1 element of the anisotropy process for each station. |
Sigma22 |
Vector of length N; contains the 2-2 element of the anisotropy process for each station. |
Sigma12 |
Vector of length N; contains the 1-2 element of the anisotropy process for each station. |
nu |
Scalar; Matern smoothness parameter. |
log_sigma_vec |
Vector of length N; log of the signal standard deviation. |
lognuggetSD |
Vector of length N; log of the error standard deviation. |
Returns a sparse matrix (N x 3) of the nonzero elements of the product between the core and sparse kernel.
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.
determineNeighbors(coords, k)determineNeighbors(coords, k)
coords |
N x 2 array of N 2-dimensional (x,y) spatial coordinates. |
k |
Scalar; number of neighbors |
An N x k matrix of nearest neighbor indices
coords <- cbind(runif(100), runif(100)) determineNeighbors(coords, 20)coords <- cbind(runif(100), runif(100)) determineNeighbors(coords, 20)
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.
dmnorm_gp2Scale(x, mean, Cov, N, Nnz, log = 1)dmnorm_gp2Scale(x, mean, Cov, N, Nnz, log = 1)
x |
Vector of measurements |
mean |
Vector of mean values |
Cov |
Matrix of size N x N; sparse kernel |
N |
Number of measurements in x |
Nnz |
Number of measurements in x |
log |
Logical; should the density be evaluated on the log scale. |
Returns the Gaussian likelihood using the gp2Scale sparse covariance.
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.
dmnorm_nngp(x, mean, AD, nID, N, k, log)dmnorm_nngp(x, mean, AD, nID, N, k, log)
x |
N-vector of data. |
mean |
N-vector with current values of the mean |
AD |
N x (k+1) matrix; the first k columns are the 'A' matrix, and the last column is the 'D' vector. |
nID |
N x k matrix of neighbor indices. |
N |
Scalar; number of data measurements. |
k |
Scalar; number of nearest neighbors. |
log |
Scalar; should the density be on the log scale (1) or not (0). |
The NNGP approximate density.
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.
dmnorm_sgv(x, mean, U, N, k, log = 1)dmnorm_sgv(x, mean, U, N, k, log = 1)
x |
Vector of measurements |
mean |
Vector of mean valiues |
U |
Matrix of size N x 3; representation of a sparse N x N Cholesky of the precision matrix. The first two columns contain row and column indices, respectively, and the last column is the nonzero elements of the matrix. |
N |
Number of measurements in x |
k |
Number of neighbors for the SGV approximation. |
log |
Logical; should the density be evaluated on the log scale. |
Returns the SGV approximation to the Gaussian likelihood.
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.
inverseEigen(eigen_comp1, eigen_comp2, eigen_comp3, which_Sigma)inverseEigen(eigen_comp1, eigen_comp2, eigen_comp3, which_Sigma)
eigen_comp1 |
N-vector; contains values of the log of the first anisotropy eigenvalue for a set of locations. |
eigen_comp2 |
N-vector; contains values of the log of the second anisotropy eigenvalue for a set of locations. |
eigen_comp3 |
N-vector; contains values of the rescaled logit of the anisotropy rotation for a set of locations. |
which_Sigma |
Scalar; one of |
A vector of anisotropy values (Sigma11, Sigma22, or Sigma12; depends
on which_Sigma) for the corresponding set of locations.
# 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# 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
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.
matern_corr(dist, rho, nu)matern_corr(dist, rho, nu)
dist |
N x N matrix; contains values of pairwise Euclidean distances in the x-y plane. |
rho |
Scalar; "range" parameter used to rescale distances |
nu |
Scalar; Matern smoothness parameter. |
A correlation matrix for a fixed set of stations and fixed parameter values.
# 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)# 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
nimble_sparse_chol(i, j, x, n)nimble_sparse_chol(i, j, x, n)
i |
Vector of row indices. |
j |
Vector of column indices. |
x |
Vector of values in the matrix. |
n |
Length of the vector |
nimble_sparse_chol
nimble_sparse_cholesky(i, j, x)nimble_sparse_cholesky(i, j, x)
i |
Vector of row indices. |
j |
Vector of column indices. |
x |
Vector of values in the matrix. |
nimble_sparse_crossprod
nimble_sparse_crossprod(i, j, x, z, n, subset, transp)nimble_sparse_crossprod(i, j, x, z, n, subset, transp)
i |
Vector of row indices. |
j |
Vector of column indices. |
x |
Vector of values in the matrix. |
z |
Vector to calculate the cross-product with. |
n |
Length of the vector |
subset |
Optional vector of rows to include in the calculation. |
transp |
Optional indicator of using the transpose |
nimble_sparse_solve
nimble_sparse_solve(i, j, x, z)nimble_sparse_solve(i, j, x, z)
i |
Vector of row indices. |
j |
Vector of column indices. |
x |
Vector of values in the matrix. |
z |
Vector to calculate the cross-product with. |
nimble_sparse_crossprod
nimble_sparse_solveMat(i, j, x, z, transp)nimble_sparse_solveMat(i, j, x, z, transp)
i |
Vector of row indices. |
j |
Vector of column indices. |
x |
Vector of values in the matrix. |
z |
Vector to calculate the cross-product with. |
transp |
Optional indicator of using the transpose |
nimble_sparse_tcrossprod
nimble_sparse_tcrossprod(i, j, x, subset)nimble_sparse_tcrossprod(i, j, x, subset)
i |
Vector of row indices. |
j |
Vector of column indices. |
x |
Vector of values in the matrix. |
subset |
Optional vector of rows to include in the calculation. |
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.
nsCorr(dist1_sq, dist2_sq, dist12, Sigma11, Sigma22, Sigma12, nu, d)nsCorr(dist1_sq, dist2_sq, dist12, Sigma11, Sigma22, Sigma12, nu, d)
dist1_sq |
N x N matrix; contains values of pairwise squared distances in the x-coordinate. |
dist2_sq |
N x N matrix; contains values of pairwise squared distances in the y-coordinate. |
dist12 |
N x N matrix; contains values of pairwise signed cross-
distances between the x- and y-coordinates. The sign of each element is
important; see |
Sigma11 |
Vector of length N; contains the 1-1 element of the anisotropy process for each station. |
Sigma22 |
Vector of length N; contains the 2-2 element of the anisotropy process for each station. |
Sigma12 |
Vector of length N; contains the 1-2 element of the anisotropy process for each station. |
nu |
Scalar; Matern smoothness parameter. |
d |
Scalar; dimension of the spatial coordinates. |
A correlation matrix for a fixed set of stations and fixed parameter values.
# 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))# 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))
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.
nsCrosscorr( Xdist1_sq, Xdist2_sq, Xdist12, Sigma11, Sigma22, Sigma12, PSigma11, PSigma22, PSigma12, nu, d )nsCrosscorr( Xdist1_sq, Xdist2_sq, Xdist12, Sigma11, Sigma22, Sigma12, PSigma11, PSigma22, PSigma12, nu, d )
Xdist1_sq |
M x N matrix; contains values of pairwise squared cross-distances in the x-coordinate. |
Xdist2_sq |
M x N matrix; contains values of pairwise squared cross-distances in the y-coordinate. |
Xdist12 |
M x N matrix; contains values of pairwise signed cross/cross-
distances between the x- and y-coordinates. The sign of each element is
important; see |
Sigma11 |
Vector of length N; contains the 1-1 element of the anisotropy process for each observed location. |
Sigma22 |
Vector of length N; contains the 2-2 element of the anisotropy process for each observed location. |
Sigma12 |
Vector of length N; contains the 1-2 element of the anisotropy process for each observed location. |
PSigma11 |
Vector of length N; contains the 1-1 element of the anisotropy process for each prediction location. |
PSigma22 |
Vector of length N; contains the 2-2 element of the anisotropy process for each prediction location. |
PSigma12 |
Vector of length N; contains the 1-2 element of the anisotropy process for each prediction location. |
nu |
Scalar; Matern smoothness parameter. |
d |
Scalar; dimension of the spatial domain. |
A M x N cross-correlation matrix for two fixed sets of stations and fixed parameter values.
# 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))# 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))
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.
nsCrossdist(coords, Pcoords, scale_factor = NULL, isotropic = FALSE)nsCrossdist(coords, Pcoords, scale_factor = NULL, isotropic = FALSE)
coords |
N x 2 matrix; contains x-y coordinates of station (observed) locations. |
Pcoords |
M x 2 matrix; contains x-y coordinates of prediction locations. |
scale_factor |
Scalar; optional argument for re-scaling the distances. |
isotropic |
Logical; indicates whether distances should be calculated
using Euclidean distance ( |
A list of distances matrices, with the following components:
dist1_sq |
M x N matrix; contains values of pairwise squared cross- distances in the x-coordinate. |
dist2_sq |
M x N matrix; contains values of pairwise squared cross- distances in the y-coordinate. |
dist12 |
M x N matrix; contains values of pairwise signed cross- distances between the x- and y-coordinates. |
scale_factor |
Value of the scale factor used to rescale distances. |
# Generate some coordinates coords <- cbind(runif(100),runif(100)) Pcoords <- cbind(runif(200),runif(200)) # Calculate distances Xdist_list <- nsCrossdist(coords, Pcoords)# Generate some coordinates coords <- cbind(runif(100),runif(100)) Pcoords <- cbind(runif(200),runif(200)) # Calculate distances Xdist_list <- nsCrossdist(coords, Pcoords)
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().
nsCrossdist3d( coords, predCoords, P_nID, scale_factor = NULL, isotropic = FALSE )nsCrossdist3d( coords, predCoords, P_nID, scale_factor = NULL, isotropic = FALSE )
coords |
N x d matrix; contains the x-y coordinates of stations. |
predCoords |
M x d matrix |
P_nID |
N x k matrix; contains indices of nearest neighbors. |
scale_factor |
Scalar; optional argument for re-scaling the distances. |
isotropic |
Logical; indicates whether distances should be calculated
separately for each coordinate dimension (FALSE) or simultaneously for all
coordinate dimensions (TRUE). |
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.
# 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)# 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)
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.
nsDist(coords, scale_factor = NULL, isotropic = FALSE)nsDist(coords, scale_factor = NULL, isotropic = FALSE)
coords |
N x 2 matrix; contains the x-y coordinates of stations |
scale_factor |
Scalar; optional argument for re-scaling the distances. |
isotropic |
Logical; indicates whether distances should be calculated
separately for each coordinate dimension (FALSE) or simultaneously for all
coordinate dimensions (TRUE). |
A list of distances matrices, with the following components:
dist1_sq |
N x N matrix; contains values of pairwise squared distances in the x-coordinate. |
dist2_sq |
N x N matrix; contains values of pairwise squared distances in the y-coordinate. |
dist12 |
N x N matrix; contains values of pairwise signed cross- distances between the x- and y-coordinates. |
scale_factor |
Value of the scale factor used to rescale distances. |
# 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)# 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)
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().
nsDist3d(coords, nID, scale_factor = NULL, isotropic = FALSE)nsDist3d(coords, nID, scale_factor = NULL, isotropic = FALSE)
coords |
N x 2 matrix; contains the x-y coordinates of stations. |
nID |
N x k matrix; contains indices of nearest neighbors. |
scale_factor |
Scalar; optional argument for re-scaling the distances. |
isotropic |
Logical; indicates whether distances should be calculated
separately for each coordinate dimension (FALSE) or simultaneously for all
coordinate dimensions (TRUE). |
Arrays with nearest neighbor distances in each coordinate direction.
# Generate some coordinates and neighbors coords <- cbind(runif(100),runif(100)) nID <- determineNeighbors(coords, 10) # Calculate distances nsDist3d(coords, nID)# Generate some coordinates and neighbors coords <- cbind(runif(100),runif(100)) nID <- determineNeighbors(coords, 10) # Calculate distances nsDist3d(coords, nID)
This function sets up and compiles a nimble model for a general nonstationary Gaussian process.
nsgpModel( tau_model = "constant", sigma_model = "constant", Sigma_model = "constant", sparse_model = "none", mu_model = "constant", likelihood = "fullGP", coords, data, constants = list(), monitorAllSampledNodes = TRUE, ... )nsgpModel( tau_model = "constant", sigma_model = "constant", Sigma_model = "constant", sparse_model = "none", mu_model = "constant", likelihood = "fullGP", coords, data, constants = list(), monitorAllSampledNodes = TRUE, ... )
tau_model |
Character; specifies the model to be used for the log(tau)
process. Options are |
sigma_model |
Character; specifies the model to be used for the
log(sigma) process. See |
Sigma_model |
Character; specifies the model to be used for the
Sigma anisotropy process. Options are |
sparse_model |
Character; specifies the model to be used for the sparse
kernel when using the gp2Scale likelihood. Defaults to |
mu_model |
Character; specifies the model to be used for the mu mean
process. Options are |
likelihood |
Character; specifies the likelihood model. Options are
|
coords |
N x d matrix of spatial coordinates. |
data |
N-vector; observed vector of the spatial process of interest |
constants |
A list of constants required to build the model; depends on the specific parameter process models chosen. |
monitorAllSampledNodes |
Logical; indicates whether all sampled nodes
should be stored ( |
... |
Additional arguments can be passed to the function; for example,
as an alternative to the |
A nimbleCode object.
# 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 )# 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 )
nsgpPredict conducts posterior prediction for MCMC samples generated
using nimble and nsgpModel.
nsgpPredict( model, samples, coords.predict, predict.process = TRUE, constants = list(), seed = 0, ... )nsgpPredict( model, samples, coords.predict, predict.process = TRUE, constants = list(), seed = 0, ... )
model |
A NSGP nimble object; the output of |
samples |
A matrix of |
coords.predict |
M x d matrix of prediction coordinates. |
predict.process |
Logical; determines whether the prediction corresponds to
the y(·) process ( |
constants |
An optional list of contants to use for prediction; alternatively, additional arguments can be passed to the function via the ... argument. |
seed |
An optional random seed argument for reproducibility. |
... |
Additional arguments can be passed to the function; for example,
as an alternative to the |
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.
## 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)## 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)
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).
orderCoordinatesMMD(coords, exact = FALSE)orderCoordinatesMMD(coords, exact = FALSE)
coords |
N x 2 array of N 2-dimensional (x,y) spatial coordinates. |
exact |
Logical; |
A list of distances matrices, with the following components:
orderedCoords |
N x 2 matrix; contains the ordered spatial coordinates
as |
orderedIndicesNoNA |
N-vector; contains the ordered indices with any NA values removed. |
coords <- cbind(runif(100), runif(100)) orderCoordinatesMMD(coords)coords <- cbind(runif(100), runif(100)) orderCoordinatesMMD(coords)
R_sparse_chol
R_sparse_chol(i, j, x, n)R_sparse_chol(i, j, x, n)
i |
Vector of row indices. |
j |
Vector of column indices. |
x |
Vector of values in the matrix. |
n |
Length of the vector |
R_sparse_chol
R_sparse_cholesky(i, j, x)R_sparse_cholesky(i, j, x)
i |
Vector of row indices. |
j |
Vector of column indices. |
x |
Vector of values in the matrix. |
nimble_sparse_crossprod
R_sparse_crossprod(i, j, x, z, n, subset = -1, transp = 1)R_sparse_crossprod(i, j, x, z, n, subset = -1, transp = 1)
i |
Vector of row indices. |
j |
Vector of column indices. |
x |
Vector of values in the matrix. |
z |
Vector to calculate the cross-product with. |
n |
Length of the vector |
subset |
Optional vector of rows to include in the calculation. |
transp |
Optional indicator of using the transpose |
nimble_sparse_solve
R_sparse_solve(i, j, x, z)R_sparse_solve(i, j, x, z)
i |
Vector of row indices. |
j |
Vector of column indices. |
x |
Vector of values in the matrix. |
z |
Vector to calculate the cross-product with. |
nimble_sparse_crossprod
R_sparse_solveMat(i, j, x, z, transp = 1)R_sparse_solveMat(i, j, x, z, transp = 1)
i |
Vector of row indices. |
j |
Vector of column indices. |
x |
Vector of values in the matrix. |
z |
Vector to calculate the cross-product with. |
transp |
Optional indicator of using the transpose |
nimble_sparse_tcrossprod
R_sparse_tcrossprod(i, j, x, subset = -1)R_sparse_tcrossprod(i, j, x, subset = -1)
i |
Vector of row indices. |
j |
Vector of column indices. |
x |
Vector of values in the matrix. |
subset |
Optional vector of rows to include in the calculation. |
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.
rmnorm_gp2Scale(n, mean, Cov, N, Nnz)rmnorm_gp2Scale(n, mean, Cov, N, Nnz)
n |
Number of realizations to generate |
mean |
Vector of mean values |
Cov |
Matrix of size N x N; sparse kernel |
N |
Number of measurements in x |
Nnz |
Number of measurements in x |
Not applicable.
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.
rmnorm_nngp(n, mean, AD, nID, N, k)rmnorm_nngp(n, mean, AD, nID, N, k)
n |
N-vector of data. |
mean |
N-vector with current values of the mean |
AD |
N x (k+1) matrix; the first k columns are the 'A' matrix, and the last column is the 'D' vector. |
nID |
N x k matrix of neighbor indices. |
N |
Scalar; number of data measurements. |
k |
Scalar; number of nearest neighbors. |
The NNGP approximate density.
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.
rmnorm_sgv(n, mean, U, N, k)rmnorm_sgv(n, mean, U, N, k)
n |
Vector of measurements |
mean |
Vector of mean valiues |
U |
Matrix of size N x 3; representation of a sparse N x N Cholesky of the precision matrix. The first two columns contain row and column indices, respectively, and the last column is the nonzero elements of the matrix. |
N |
Number of measurements in x |
k |
Number of neighbors for the SGV approximation. |
Not applicable.
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.
sgvSetup( coords, coords_pred = NULL, k = 15, seed = NULL, pred.seed = NULL, order_coords = TRUE, order_coords_pred = TRUE )sgvSetup( coords, coords_pred = NULL, k = 15, seed = NULL, pred.seed = NULL, order_coords = TRUE, order_coords_pred = TRUE )
coords |
Matrix of observed locations. |
coords_pred |
Optional matrix of prediction locations. |
k |
Number of neighbors. |
seed |
Setting the seed for reproducibility of the observed location ordering |
pred.seed |
Setting the seed for reproducibility of the prediction ordering. |
order_coords |
Logical; should the coordinates be ordered. |
order_coords_pred |
Logical; should the coordinates be ordered. |
A list with the following components:
ord |
A vector of ordering position for the observed locations. |
ord_pred |
A vector of ordering position for the prediction
locations (if |
ord_all |
A concatenated vector of |
coords_ord |
A matrix of ordered locations (observed and prediction), included for convenience. |
nID_ord |
A matrix of (ordered) neighbor indices. |
condition_on_y_ord |
A matrix indicating whether the conditioning
set for each (ordered) location is on the latent process (y, |