excursions (version 2.5.8)

continuous: Calculate continuous domain excursion and credible contour sets

Description

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
)

Value

A list with elements

M

SpatialPolygons or inla.mesh.segment 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).

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.

Author

Finn Lindgren

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.

Examples

if (FALSE) {
if (require("fmesher")) {
  # 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))

  ## 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])
}
}