sprates {nimblewomble}R Documentation

Posterior samples for rates of change

Description

Posterior samples for rates of change

Usage

sprates(
  coords = NULL,
  grid = NULL,
  model = NULL,
  kernel = c("matern1", "matern2", "gaussian")
)

Arguments

coords

coordinates

grid

grid for sampling the rates of change

model

posterior samples of Z(s), \phi, \sigma^2

kernel

choice of kernel; must be one of "matern1", "matern2", "gaussian"

Value

A list containing MCMC samples for gradients and curvatures and the associated estimates and 95

Author(s)

Aritra Halder <aritra.halder@drexel.edu>,
Sudipto Banerjee <sudipto@ucla.edu>

Examples


require(nimble)
require(nimblewomble)

set.seed(1)
# Generated Simulated Data
N = 1e2
tau = 1
coords = matrix(runif(2 * N, -10, 10), ncol = 2); colnames(coords) = c("x", "y")
y = rnorm(N, mean = 20 * sin(sqrt(coords[, 1]^2  + coords[, 2]^2)), sd = tau)

# Create equally spaced grid of points
xsplit = ysplit = seq(-10, 10, by = 1)[-c(1, 21)]
grid = as.matrix(expand.grid(xsplit, ysplit), ncol = 2)
colnames(grid) = c("x", "y")

####################################
# Process for True Rates of Change #
####################################
# Gradient along x
true_sx = round(20 * cos(sqrt(grid[,1]^2 + grid[,2]^2)) *
                grid[,1]/sqrt(grid[,1]^2 + grid[,2]^2), 3)
# Gradient along y
true_sy = round(20 * cos(sqrt(grid[,1]^2 + grid[,2]^2)) *
                grid[,2]/sqrt(grid[,1]^2 + grid[,2]^2), 3)
# Curvature along x
true_sxx = round(20 * cos(sqrt(grid[,1]^2 + grid[,2]^2))/
                  sqrt(grid[,1]^2 + grid[,2]^2) -
                 20 * cos(sqrt(grid[,1]^2 + grid[,2]^2)) *
                 grid[,1]^2/(grid[,1]^2 + grid[,2]^2)^(3/2) -
                 20 * sin(sqrt(grid[,1]^2 + grid[,2]^2)) *
                 grid[,1]^2/(grid[,1]^2 + grid[,2]^2), 3)
# Mixed Curvature
true_sxy = round(-20 * (cos(sqrt(grid[,1]^2 + grid[,2]^2)) -
                 sin(sqrt(grid[,1]^2 + grid[,2]^2))) * grid[,1]
                  * grid[,2]/(grid[,1]^2 + grid[,2]^2), 3)
# Curvature along y
true_syy = round(20 * cos(sqrt(grid[,1]^2 + grid[,2]^2))/
                  sqrt(grid[,1]^2 + grid[,2]^2) -
                 20 * cos(sqrt(grid[,1]^2 + grid[,2]^2)) *
                 grid[,2]^2/(grid[,1]^2 + grid[,2]^2)^(3/2) -
                 20 * sin(sqrt(grid[,1]^2 + grid[,2]^2)) *
                 grid[,2]^2/(grid[,1]^2 + grid[,2]^2), 3)
# Create the plots
p1 = sp_ggplot(data_frame = data.frame(coords, z = y))
p2 = sp_ggplot(data_frame = data.frame(grid[-which(is.nan(true_sx)),],
               z = true_sx[-which(is.nan(true_sx))]))
p3 = sp_ggplot(data_frame = data.frame(grid[-which(is.nan(true_sy)),],
               z = true_sy[-which(is.nan(true_sy))]))
p4 = sp_ggplot(data_frame = data.frame(grid[-which(is.nan(true_sxx)),],
               z = true_sxx[-which(is.nan(true_sxx))]))
p5 = sp_ggplot(data_frame = data.frame(grid[-which(is.nan(true_sxy)),],
               z = true_sxy[-which(is.nan(true_sxy))]))
p6 = sp_ggplot(data_frame = data.frame(grid[-which(is.nan(true_syy)),],
               z = true_syy[-which(is.nan(true_syy))]))

##########################
# Fit a Gaussian Process #
##########################
# Posterior samples for theta
mc_sp = gp_fit(coords = coords, y = y, kernel = "matern2")
# Posterior samples for Z(s) and beta
model = zbeta_samples(y = y, coords = coords,
                      model = mc_sp$mcmc,
                      kernel = "matern2")

###################
# Rates of Change #
###################
gradients = sprates(grid = grid,
                    coords = coords,
                    model = model,
                    kernel = "matern2")
p8 = sp_ggplot(data_frame = data.frame(grid,
                z = gradients$estimate.sx[,"50%"],
               sig = gradients$estimate.sx$sig))
p9 = sp_ggplot(data_frame = data.frame(grid,
               z = gradients$estimate.sy[,"50%"],
               sig = gradients$estimate.sy$sig))
p10 = sp_ggplot(data_frame = data.frame(grid,
                 z = gradients$estimate.sxx[,"50%"],
                sig = gradients$estimate.sxx$sig))
p11 = sp_ggplot(data_frame = data.frame(grid,
                 z = gradients$estimate.sxy[,"50%"],
                 sig = gradients$estimate.sxy$sig))
p12 = sp_ggplot(data_frame = data.frame(grid,
                z = gradients$estimate.syy[,"50%"],
                sig = gradients$estimate.syy$sig))


[Package nimblewomble version 0.1.0 Index]