spwombling {nimblewomble}R Documentation

Posterior samples for wombling measures

Description

Posterior samples for wombling measures

Usage

spwombling(
  coords = NULL,
  curve = NULL,
  model = NULL,
  kernel = c("matern1", "matern2", "gaussian")
)

Arguments

coords

coordinates

curve

coordinates of the curve for wombling

model

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

kernel

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

Value

A list containing posterior samples of wombling measures and associated estimates and their 95

Author(s)

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

Examples

## Not run: 
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")

##########################
# 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")
############
# Wombling #
############
# Pick any curve (contour) of your choice
# curve = your contour
tvec = sapply(1:(nrow(curve) - 1), function(x){
sqrt(sum((curve[(x + 1),] - curve[x,])^2))})
umat = as.matrix(t(sapply(1:(nrow(curve) - 1), function(x){
 (curve[(x + 1),] - curve[x,])}))/tvec)

wm = spwombling(coords = coords,
                curve = curve,
                model = model,
                kernel = "matern2")

# Total wombling measure for gradient
colSums(wm$estimate.wm.1[,-4]); colSums(wm$estimate.wm.1[,-4])/sum(tvec)
# Total wombling measure for curvature
colSums(wm$estimate.wm.2[,-4]); colSums(wm$estimate.wm.2[,-4])/sum(tvec)


# Color code points based on significance
col.pts.1 = sapply(wm$estimate.wm.1$sig, function(x){
  if(x == 1) return("green")
  else if(x == -1) return("cyan")
  else return(NA)
  })

  col.pts.2 = sapply(wm$estimate.wm.2$sig, function(x){
  if(x == 1) return("green")
  else if(x == -1) return("cyan")
  else return(NA)
  })

p13 = sp_ggplot(data_frame = data.frame(coords, y))
p14 = p13 + geom_path(curve, mapping = aes(x, y), linewidth = 2)
p15 = p13 + geom_path(curve, mapping = aes(x, y), linewidth = 2) +
geom_path(curve, mapping = aes(x, y),
          colour = c(col.pts.1, NA), linewidth = 1, na.rm = TRUE)
p16 = p13 + geom_path(curve, mapping = aes(x, y), linewidth = 2) +
            geom_path(curve, mapping = aes(x, y),
            colour = c(col.pts.2, NA), linewidth = 1, na.rm = TRUE)


###############
# True Values #
###############
truth = matrix(0, nrow = nrow(curve) - 1, ncol = 2)
rule = seq(0, 1, by = 0.01)

for(i in 1:(nrow(curve) - 1)){
  u.perp = c(umat[i, 2], - umat[i, 1])
  s0 = curve[i,]

truth.lsegment = sapply(rule * tvec[i], function(x){
  s.t = s0 + x * umat[i,]
  true_sx = 20 * cos(sqrt(s.t[1]^2 + s.t[2]^2)) * s.t[1]/
            sqrt(s.t[1]^2 + s.t[2]^2)
  true_sy = 20 * cos(sqrt(s.t[1]^2 + s.t[2]^2)) * s.t[2]/
            sqrt(s.t[1]^2 + s.t[2]^2)
  true_sx * u.perp[1] + true_sy * u.perp[2]
  })

truth[i, 1] = sum(truth.lsegment * (tvec[i]/101))

truth.lsegment = sapply(rule * tvec[i], function(x){
  s.t = s0 + x * umat[i,]
  true_sxx = 20 * cos(sqrt(s.t[1]^2 + s.t[2]^2))/sqrt(s.t[1]^2 + s.t[2]^2) -
    20 * cos(sqrt(s.t[1]^2 + s.t[2]^2)) *
            s.t[1]^2/(s.t[1]^2 + s.t[2]^2)^(3/2) -
    20 * sin(sqrt(s.t[1]^2 + s.t[2]^2)) * s.t[1]^2/(s.t[1]^2 + s.t[2]^2)
  true_sxy = -20 * (cos(sqrt(s.t[1]^2 + s.t[2]^2)) -
                sin(sqrt(s.t[1]^2 + s.t[2]^2))) *
                      s.t[1] * s.t[2]/(s.t[1]^2 + s.t[2]^2)
  true_syy = 20 * cos(sqrt(s.t[1]^2 + s.t[2]^2))/sqrt(s.t[1]^2 + s.t[2]^2) -
    20 * cos(sqrt(s.t[1]^2 + s.t[2]^2)) *
               s.t[2]^2/(s.t[1]^2 + s.t[2]^2)^(3/2) -
    20 * sin(sqrt(s.t[1]^2 + s.t[2]^2)) * s.t[2]^2/(s.t[1]^2 + s.t[2]^2)
  true_sxx * u.perp[1]^2 + 2 * true_sxy * u.perp[1] * u.perp[2] +
                true_syy * u.perp[2]^2
  })
  truth[i, 2] = sum(truth.lsegment * (tvec[i]/101))
}
true.total = colSums(truth); true.total
true.avg.total = true.total/sum(tvec); true.avg.total

## End(Not run)

[Package nimblewomble version 0.1.0 Index]