contourlevel {cooltools} | R Documentation |
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 |
p |
vector of probabilities |
xmin , xmax |
(only used if |
neval |
(only used if |
napprox |
(only used if |
... |
(only used if |
Value
Returns a vector of levels l, which has the same length as p.
Author(s)
Danail Obreschkow
See Also
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)