sFit: Fit a spatially mean-zero spatial Gaussian process model

In jmhewitt/bisque: Approximate Bayesian Inference via Sparse Grid Quadrature Evaluation (BISQuE) for Hierarchical Models

Description

Uses a Gibbs sampler to estimate the parameters of a Matern covariance function used to model observations from a Gaussian process with mean 0.

Usage

sFit(
  x,
  coords,
  nSamples,
  thin = 1,
  rw.initsd = 0.1,
  inits = list(),
  C = 1,
  alpha = 0.44,
  priors = list(sigmasq = list(a = 2, b = 1), rho = list(L = 0, U = 1), nu = list(L = 0,
    U = 1))
)

Arguments

x	Observation of a spatial Gaussian random field, passed as a vector
coords	Spatial coordinates of the observation
nSamples	(thinned) number of MCMC samples to generate
thin	thinning to be used within the returned MCMC samples
rw.initsd	initial standard devaition for random walk proposals. this parameter will be adaptively tuned during sampling
inits	list of initial parameters for the MCMC chain
C	scale factor used during tuning of the random walk proposal s.d.
alpha	target acceptance rate for which the random walk proposals should optimize
priors	parameters to specify the prior distributions for the model

Examples

library(fields)

simulate.field = function(n = 100, range = .3, smoothness = .5, phi = 1){
  # Simulates a mean-zero spatial field on the unit square
  #
  # Parameters:
  #  n - number of spatial locations
  #  range, smoothness, phi - parameters for Matern covariance function
  
  coords = matrix(runif(2*n), ncol=2)
  
  Sigma = Matern(d = as.matrix(dist(coords)), 
                 range = range, smoothness = smoothness, phi = phi)
  
  list(coords = coords,
       params = list(n=n, range=range, smoothness=smoothness, phi=phi),
       x = t(chol(Sigma)) %*%  rnorm(n))
}

# simulate data
x = simulate.field()

# configure gibbs sampler  
it = 100

# run sampler using default posteriors
post.samples = sFit(x = x$x, coords = x$coords, nSamples = it)

# build kriging grid
cseq = seq(0, 1, length.out = 10)
coords.krig = expand.grid(x = cseq, y = cseq)

# sample from posterior predictive distribution
burn = 75
samples.krig = sKrig(x$x, post.samples, coords.krig = coords.krig, burn = burn)

jmhewitt/bisque documentation built on Feb. 9, 2020, 2:36 a.m.