hpd: Calculate the highest posterior density credible interval for...
In Smisc: Sego Miscellaneous

Description Usage Arguments Details Value Methods (by generic) Author(s) Examples

Calculate the highest posterior density credible interval for a unimodal density

hpd(pdf, support, prob = 0.95, cdf = NULL, njobs = 1, checkUnimodal = 0)

## S3 method for class 'hpd'
print(x, ...)

## S3 method for class 'hpd'
plot(x, ...)

`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`

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.

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

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

Landon Sego

# A credible interval using the standard normal
int <- hpd(dnorm, c(-5,5), prob = 0.90, njobs = 2)
print(int)
plot(int)


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


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