| 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) <arXiv:1702.00434v2>). 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, Mark Risser |
| Maintainer: | Daniel Turek <[email protected]> |
| License: | GPL-3 |
| Version: | 0.1.2 |
| Built: | 2025-02-16 05:31:38 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).
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_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_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_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", mu_model = "constant", likelihood = "fullGP", coords, data, constants = list(), monitorAllSampledNodes = TRUE, ... )nsgpModel( tau_model = "constant", sigma_model = "constant", Sigma_model = "constant", 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 |
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, seed = 0, ... )nsgpPredict( model, samples, coords.predict, predict.process = TRUE, constants, 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.
# 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 )# 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 )
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 |
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_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_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 )sgvSetup( coords, coords_pred = NULL, k = 15, seed = NULL, pred.seed = NULL, order_coords = 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. |
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, |