spLMPredictJoint | R Documentation |
The function spLMPredictJoint
collects posterior predictive samples
for a set of new locations given a spLM
object from the spBayes
package.
spLMPredictJoint(sp.obj, pred.coords, pred.covars, start = 1,
end = nrow(sp.obj$p.theta.samples), thin = 1, verbose = TRUE, n.report = 100,
noisy = FALSE, method = "eigen")
sp.obj |
An |
pred.coords |
An |
pred.covars |
An |
start |
Specifies the first sample included in the composition sampling. |
end |
Specifies the last sample included in the composition.
The default is to use all posterior samples in |
thin |
A sample thinning factor. The default of 1 considers all
samples between |
verbose |
If |
n.report |
The interval to report sampling progress. |
noisy |
If |
method |
Method used to decompose covariance matrix. Options are "chol", "eigen", and "svd" for the Cholesky, Eigen, and singular value decomposition approaches, respectively. |
This function samples from the joint posterior predictive distribution of a Bayesian spatial linear model. Specifically, it is intended to be similar to the spPredict
function in the spBayes
except that it samples from the joint distribution instead of the marginal distribution. However, it will only work for spLM
objects and should have the same limitations as the spLM
and spPredict
functions. Note that the spRecover
function is called internally to recover the posterior samples form the posterior distribution of the spatial model.
The function returns a np \times B
matrix of posterior predictive samples, where B
is the number of posterior samples. The class is jointPredicitveSample
.
Joshua French
spLM, spPredict, spRecover
# Set parameters
n <- 100
np <- 12
n.samples <- 10
burnin.start <- .5 * n.samples + 1
sigmasq <- 1
tausq <- 0.0
phi <- 1
cov.model <- "exponential"
n.report <- 5
# Generate coordinates
coords <- matrix(runif(2 * n), ncol = 2);
pcoords <- as.matrix(expand.grid(seq(0, 1, len = 12), seq(0, 1, len = 12)))
# Construct design matrices
X <- as.matrix(cbind(1, coords))
Xp <- cbind(1, pcoords)
# Specify priors
starting <- list("phi" = phi, "sigma.sq"= sigmasq, "tau.sq" = tausq)
tuning <- list("phi"=0.1, "sigma.sq"=0.1, "tau.sq"=0.1)
priors.1 <- list("beta.Norm"=list(c(1, 2, 1), diag(100, 3)),
"phi.Unif"=c(0.00001, 10), "sigma.sq.IG"=c(1, 1))
# Generate data
B <- rnorm(3, c(1, 2, 1), sd = 10)
phi <- runif(1, 0, 10)
sigmasq <- 1/rgamma(1, 1, 1)
V <- simple.cov.sp(D = dist1(coords), cov.model, c(sigmasq, 1/phi), error.var = tausq,
smoothness = nu, finescale.var = 0)
y <- X %*% B + rmvnorm(1, rep(0, n), V) + rnorm(n, 0, sqrt(tausq))
# Create spLM object
library(spBayes)
m1 <- spBayes::spLM(y ~ X - 1, coords = coords, starting = starting,
tuning = tuning, priors = priors.1, cov.model = cov.model,
n.samples = n.samples, verbose = FALSE, n.report = n.report)
# Sample from joint posterior predictive distribution
y1 <- spLMPredictJoint(m1, pred.coords = pcoords, pred.covars = Xp,
start = burnin.start, verbose = FALSE, method = "chol")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.