contourlevel: Find contour levels of a d-dimensional density field

In obreschkow/cooltools: Practical Tools for Scientific Computations and Visualizations

Find contour levels of a d-dimensional density field

Description

Given a vector d-dimensional vector/array f or function f(x) of a d-element vector x, this routine evaluates the value l, such that the sum/integral of f over the domain f>l makes up a fraction p of the total sum/integral of f. The main purpose of this routine is to determine the iso-contour levels of a likelihood or density function, containing a cumulative probability-mass p.

Usage

contourlevel(
  f,
  p = c(0.6826895, 0.9544997),
  xmin = NULL,
  xmax = NULL,
  neval = 10000,
  napprox = 10,
  ...
)

Arguments

f	either a d-dimensional vector/array or a function f(x) of a d-element vector x. There is no need for f to be normalized.
p	vector of probabilities
xmin, xmax	(only used if f is a function) vectors with of lower and upper limits of x, defining the domain on which the function f(x) is evaluated. Outside this domain, f is assumed to vanish. These limits should be chosen as narrow as possible for the algorithm to converge quickly.
neval	(only used if f is a function) maximum number of function evaluations in numerical integrals
napprox	(only used if f is a function) number of points used interpolate the cumulative probability density. If set to 0, no interpolation is used.
...	(only used if f is a function) additional parameters to be passed to the function f.

Value

Returns a vector of levels l, which has the same length as p.

Author(s)

Danail Obreschkow

See Also

dpqr

Examples

## f(x) is a 1D PDF
# compute 1-sigma and 2-sigma contour levels of a normal distribution, i.e.
# the values l, such that int_{dnorm(x)>l} dnorm(x) dx = p (=68.3%, 95.4%).
l = contourlevel(dnorm,xmin=-10,xmax=10,napprox=0)
print(l)

# compare these values to the solutions dnorm(1), dnorm(2)
print(dnorm(c(1,2)))


## f(x) is a 2D likelihood function
# Produce 20%, 40%, 60% and 80% iso contours on the 2D likelihood function f(x)
f = function(x) cos(2*x[1]-x[2]-1)^2*exp(-x[1]^2-x[2]^2-x[1]*x[2])
p = c(0.2,0.4,0.6,0.8) # cumulative probability
l = contourlevel(f,p,c(-5,-5),c(5,5)) # values l, such that int_{f>l}(f(x)dx)=p int(f(x)*dx)

# Plot function and contours at the levels l
x = seq(-3,3,length=200)
m = pracma::meshgrid(x)
z = array(Vectorize(function(x,y) f(c(x,y)))(m$Y,m$X),dim(m$X))
image(x,x,z,col=terrain.colors(100))
contour(x,x,z,levels=l,add=TRUE,labels=sprintf('%.0f%%',p*100),labcex=0.7)


## f is a 20-by20 array representing a gridded pointset
# produce a set of 1000 points in 2D, drawn from a 2D normal distribution
set.seed(1)
x = MASS::mvrnorm(n=1000, mu=c(0,0), matrix(c(3,1,1,2),2,2))

# grid these points onto a regular 20-by-20 grid
g = griddata(x, min=-6, max=6)

# find 1-sigma and 2-sigma contour levels and draw contours at these levels
l = contourlevel(g$field)
plot(x, xlim=g$grid[[1]]$lim, ylim=g$grid[[2]]$lim, pch=20, cex=0.5)
contour(g$grid[[1]]$mid,g$grid[[2]]$mid,g$field,
        levels=l,add=TRUE,col='red',lwd=c(2,1),labels=NA)

obreschkow/cooltools documentation built on Nov. 16, 2024, 2:46 a.m.