continuous: Calculate continuous domain excursion and credible contour...

In finnlindgren/excursions: Excursion Sets and Contour Credibility Regions for Random Fields

Description

Calculates continuous domain excursion and credible contour sets

Usage

continuous(
  ex,
  geometry,
  alpha,
  method = c("log", "linear", "step"),
  output = c("sp", "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), inla.mesh.lattice, 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 inla gives a inla.mesh.segment 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:

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).

Author(s)

Finn Lindgren [e-mail]

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

## Not run: 
if (require.nowarnings("INLA")) {
#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=inla.mesh.lattice(x=x,y=x)
mesh=inla.mesh.create(lattice=lattice, extend=FALSE, refine=FALSE)
spde <- inla.spde2.matern(mesh, alpha=2)

#Generate an artificial sample
sigma2.e = 0.01
n.obs=100
obs.loc = cbind(runif(n.obs)*diff(range(x))+min(x),
runif(n.obs)*diff(range(x))+min(x))
Q = inla.spde2.precision(spde, theta=c(log(sqrt(0.5)), log(sqrt(1))))
x = inla.qsample(Q=Q)
A = inla.spde.make.A(mesh=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),
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)
image(matrix(sqrt(vars.post[reo]),n.lattice,n.lattice),
main="sd", axes = FALSE)
image(matrix(map$M[reo],n.lattice,n.lattice),col=cols,axes=FALSE)
idx.M = setdiff(names(sets$M), "-1")
plot(sets$M[idx.M], col=cols[idx.M])
}
## End(Not run)

finnlindgren/excursions documentation built on Jan. 17, 2021, 12:20 a.m.