hdrlevels: Highest Density Region (HDR) Levels

Description Usage Arguments Details Value Author(s) References See Also Examples

View source: R/util.R

Description

Compute the levels of Highest Density Regions (HDRs) for any density and probability levels.

Usage

1

Arguments

density

A vector of density values computed on a set of (observed) evaluation points.

prob

A vector of probability levels in the range [0,1].

Details

From Hyndman (1996), let f(x) be the density function of a random variable X. Then the 100(1-α)\% HDR is the subset R(f_α) of the sample space of X such that

R(f_α) = {x : f(x) ≥ f_α }

where f_α is the largest constant such that Pr( X \in R(f_α)) ≥ 1-α

Value

The function returns a vector of density values corresponding to HDRs at given probability levels.

Author(s)

L. Scrucca

References

Rob J. Hyndman (1996) Computing and Graphing Highest Density Regions. The American Statistician, 50(2):120-126.

See Also

plot.densityMclust

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
# Example: univariate Gaussian
x <- rnorm(1000)
f <- dnorm(x)
a <- c(0.5, 0.25, 0.1)
(f_a <- hdrlevels(f, prob = 1-a))

plot(x, f)
abline(h = f_a, lty = 2)
text(max(x), f_a, labels = paste0("f_", a), pos = 3)

mean(f > f_a[1])
range(x[which(f > f_a[1])])
qnorm(1-a[1]/2)

mean(f > f_a[2])
range(x[which(f > f_a[2])])
qnorm(1-a[2]/2)

mean(f > f_a[3])
range(x[which(f > f_a[3])])
qnorm(1-a[3]/2)

# Example 2: univariate Gaussian mixture
set.seed(1)
cl <- sample(1:2, size = 1000, prob = c(0.7, 0.3), replace = TRUE)
x <- ifelse(cl == 1, 
            rnorm(1000, mean = 0, sd = 1),
            rnorm(1000, mean = 4, sd = 1))
f <- 0.7*dnorm(x, mean = 0, sd = 1) + 0.3*dnorm(x, mean = 4, sd = 1)

a <- 0.25
(f_a <- hdrlevels(f, prob = 1-a))

plot(x, f)
abline(h = f_a, lty = 2)
text(max(x), f_a, labels = paste0("f_", a), pos = 3)

mean(f > f_a)

# find the regions of HDR
ord <- order(x)
f <- f[ord]
x <- x[ord]
x_a <- x[f > f_a]
j <- which.max(diff(x_a))
region1 <- x_a[c(1,j)]
region2 <- x_a[c(j+1,length(x_a))]
plot(x, f, type = "l")
abline(h = f_a, lty = 2)
abline(v = region1, lty = 3, col = 2)
abline(v = region2, lty = 3, col = 3)

Example output

Package 'mclust' version 5.4.3
Type 'citation("mclust")' for citing this R package in publications.
       50%        75%        90% 
0.31878769 0.19683029 0.08894239 
[1] 0.5
[1] -0.6671619  0.6681078
[1] 0.6744898
[1] 0.75
[1] -1.184999  1.187529
[1] 1.150349
[1] 0.9
[1] -1.726935  1.731267
[1] 1.644854
       75% 
0.09291342 
[1] 0.75

mclust documentation built on Nov. 20, 2020, 5:09 p.m.