Smisc (version 0.3.9.1)

hpd: Calculate the highest posterior density credible interval for a unimodal density

Description

Calculate the highest posterior density credible interval for a unimodal density

Usage

hpd(pdf, support, prob = 0.95, cdf = NULL, njobs = 1, checkUnimodal = 0)
# S3 method for hpd
print(x, ...)
# S3 method for hpd
plot(x, ...)

Arguments

pdf

Function that takes a single numeric vector argument that returns a vector of probability density values

support

A numeric vector of length 2 giving the interval over which the random variable has support (i.e. for which the pdf is positive). For now, this must be a finite interval. Intervals for random variables within infinite support can still be calculated by setting the values of support to be suitably large and/or small. See examples.

prob

A numeric value in (0, 1] indicating the size of the desired probability for the credible interval.

cdf

A function that takes a single (not necessarily vector) argument and returns the cumulative probability. If NULL, the pdf is integrated as needed to calculate probabilities as needed. However, providing the cdf will speed up calculations.

njobs

The number of parallel jobs to spawn (where possible) using doCallParallel. This is helpful if pdf is expensive.

checkUnimodal

An integer that, when greater than 0, indicates the number of points in support for which pdf is evaluated to determine whether the function appears unimodal. This is done in parallel if njobs > 1. If checkUnimodal is not 0, it should be a large number (like 1000 or more).

x

object of class hpd, returned by hpd

…

For the plot method, these are additional arguments that may be passed to plotFun, plot.default, or abline

Value

A list of class hpd that contains the following elements:

lower

The lower endpoint of the highest posterior density interval

upper

The lower endpoint of the highest posterior density interval

prob

The acheived probability of the interval

cut

The horizontal cut point that gave rise to the interval

mode

The mode of the density

pdf

The probability density function

support

The support of the pdf

Methods (by generic)

• print: Prints the lower and upper limits of the credible interval, along with the achieved probabilty of that interval.

• plot: Plots the density, overlaying the lower and upper limits of the credible interval

Details

Parallel processing (via njobs > 1) may be advantageous if 1) pdf is a function that is computationally expensive, 2) the cdf is not provided, in which case pdf is integrated, and/or 3) when checkUnimodal is large.

Examples

# NOT RUN {
# A credible interval using the standard normal
int <- hpd(dnorm, c(-5,5), prob = 0.90, njobs = 2)
print(int)
plot(int)
# }
# NOT RUN {
# A credible interval with the gamma density
int <- hpd(function(x) dgamma(x, shape = 2, rate = 0.5), c(0, 20),
           cdf = function(x) pgamma(x, shape = 2, rate = 0.5), prob = 0.8)
print(int)
plot(int)

# }
# NOT RUN {
# A credible interval using the Beta density
dens <- function(x) dbeta(x, 7, 12)
dist <- function(x) pbeta(x, 7, 12)
int <- hpd(dens, c(0, 1), cdf = dist)
print(int)
plot(int)
# }