Calculate continuous domain excursion and credible contour sets

Calculates continuous domain excursion and credible contour sets

Usage

continuous(
  ex,
  geometry,
  alpha,
  method = c("log", "linear", "step"),
  output = c("sp", "fm", "inla"),
  subdivisions = 1,
  calc.credible = TRUE
)

Arguments

ex: An excurobj object generated by a call to excursions() or contourmap().
geometry: Specification of the lattice or triangulation geometry of the input. One of list(x, y), list(loc, dims), fm_lattice_2d, inla.mesh.lattice, fm_mesh_2d, or inla.mesh, where x and y are vectors, loc is a two-column matrix of coordinates, and dims is the lattice size vector. The first three versions are all treated topologically as lattices, and the lattice boxes are assumed convex.
alpha: The target error probability. A warning is given if it is detected that the information ex isn't sufficient for the given alpha. Defaults to the value used when calculating ex.
method: The spatial probability interpolation transformation method to use. One of log, linear, or step. For log, the probabilities are interpolated linearly in the transformed scale. For step, a conservative step function is used.
output: Specifies what type of object should be generated. sp gives a SpatialPolygons object, and fm or inla gives a fm_segm object.
subdivisions: The number of mesh triangle subdivisions to perform for the interpolation of the excursions or contour function. 0 is no subdivision. The setting has a small effect on the evaluation of P0 for the log method (higher values giving higher accuracy) but the main effect is on the visual appearance of the interpolation. Default=1.
calc.credible: Logical, if TRUE (default), calculate credible contour region objects in addition to avoidance sets.

Value

A list with elements

M: SpatialPolygons or fm_segm object. The subsets are tagged, so that credible regions are tagged "-1", and regions between levels are tagged as.character(0:nlevels).
F: Interpolated F function.
G: Contour and inter-level set indices for the interpolation.
F.geometry: Mesh geometry for the interpolation.
P0: P0 measure based on interpolated F function (only for contourmap input).

References

Bolin, D. and Lindgren, F. (2017) Quantifying the uncertainty of contour maps, Journal of Computational and Graphical Statistics, vol 26, no 3, pp 513-524.

Bolin, D. and Lindgren, F. (2018), Calculating Probabilistic Excursion Sets and Related Quantities Using excursions, Journal of Statistical Software, vol 86, no 1, pp 1-20.

Author

Finn Lindgren finn.lindgren@gmail.com

Examples

if (FALSE) { # \dontrun{
if (require("fmesher") &&
  require("sp")) {
  # Generate mesh and SPDE model
  n.lattice <- 10 # Increase for more interesting, but slower, examples
  x <- seq(from = 0, to = 10, length.out = n.lattice)
  lattice <- fm_lattice_2d(x = x, y = x)
  mesh <- fm_rcdt_2d_inla(lattice = lattice, extend = FALSE, refine = FALSE)

  # Generate an artificial sample
  sigma2.e <- 0.1
  n.obs <- 100
  obs.loc <- cbind(
    runif(n.obs) * diff(range(x)) + min(x),
    runif(n.obs) * diff(range(x)) + min(x)
  )
  Q <- fm_matern_precision(mesh, alpha = 2, rho = 3, sigma = 1)
  x <- fm_sample(n = 1, Q = Q)
  A <- fm_basis(mesh, loc = obs.loc)
  Y <- as.vector(A %*% x + rnorm(n.obs) * sqrt(sigma2.e))

  ## Calculate posterior
  Q.post <- (Q + (t(A) %*% A) / sigma2.e)
  mu.post <- as.vector(solve(Q.post, (t(A) %*% Y) / sigma2.e))
  vars.post <- excursions.variances(chol(Q.post), max.threads = 1)

  ## Calculate contour map with two levels
  map <- contourmap(
    n.levels = 2, mu = mu.post, Q = Q.post,
    alpha = 0.1, F.limit = 0.1, max.threads = 1
  )

  ## Calculate the continuous representation
  sets <- continuous(map, mesh, alpha = 0.1)

  ## Plot the results
  reo <- mesh$idx$lattice
  cols <- contourmap.colors(map,
    col = heat.colors(100, 1, rev = TRUE),
    credible.col = grey(0.5, 1)
  )
  names(cols) <- as.character(-1:2)

  par(mfrow = c(2, 2))
  image(matrix(mu.post[reo], n.lattice, n.lattice),
    main = "mean", axes = FALSE, asp = 1
  )
  image(matrix(sqrt(vars.post[reo]), n.lattice, n.lattice),
    main = "sd", axes = FALSE, asp = 1
  )
  image(matrix(map$M[reo], n.lattice, n.lattice),
    col = cols, axes = FALSE, asp = 1
  )
  idx.M <- setdiff(names(sets$M), "-1")
  plot(sets$M[idx.M], col = cols[idx.M])
}
} # }