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

Description Usage Arguments Examples

View source: R/sFit.R

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
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.