nimblewomble (version 0.1.0)

sprates: 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")
)

Value

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

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"

Author

Aritra Halder <[e-mail]>,
Sudipto Banerjee <[e-mail]>

Examples

# \donttest{
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))
# }