betaHPD function

compute and optionally plot beta HDRs

Compute and optionally plot highest density regions for the Beta distribution.

betaHPD(alpha,beta,p=.95,plot=FALSE,xlim=NULL,debug=FALSE)

Arguments

• alpha: scalar, first shape parameter of the Beta density. Must be greater than 1, see details
• beta: scalar, second shape parameter of the Beta density. Must be greater than 1, see details
• p: scalar, content of HPD, must lie between 0 and 1
• plot: logical flag, if TRUE then plot the density and show the HDR
• xlim: numeric vector of length 2, the limits of the density's support to show when plotting; the default is NULL, in which case the function will confine plotting to where the density is non-negligible
• debug: logical flag, if TRUE produce messages to the console

Details

The Beta density arises frequently in Bayesian models of binary events, rates, and proportions, which take on values in the open unit interval. For instance, the Beta density is a conjugate prior for the unknown success probability in binomial trials. With shape parameters $\alpha > 1$ and $\beta > 1$, the Beta density is unimodal.

In general, suppose $\theta \in \Theta \subseteq R^k$

is a random variable with density $f(\theta)$. A highest density region (HDR) of $f(\theta)$ with content c("$p \\in\n$", "$(0,1]$") is a set $\mathcal{Q} \subseteq \Theta$ with the following properties:

$$\int_\mathcal{Q} f(\theta) d\theta = p$$

and

$$f(\theta) > f(\theta^*) \, \forall\\theta \in \mathcal{Q},\theta^* \not\in \mathcal{Q}.$$

For a unimodal Beta density (the class of Beta densities handled by this function), a HDR of content $0 < p < 1$ is simply an interval $\mathcal{Q} \in (0,1)$.

This function uses numerical methods to solve for the end points of a HDR for a Beta density with user-specified shape parameters, via repeated calls to the functions dbeta, pbeta and qbeta. The function optimize is used to find points $v$ and $w$

such that

$$f(v) = f(w)$$

subject to the constraint

$$\int_v^w f(\theta; \alpha, \beta) d\theta = p,$$

where $f(\theta; \alpha, \beta)$ is a Beta density with shape parameters $\alpha$ and $\beta$.

In the special case of $\alpha = \beta > 1$, the end points of a HDR with content $p$ are given by the $(1 \pm p)/2$

quantiles of the Beta density, and are computed with the qbeta function.

Again note that the function will only compute a HDR for a unimodal Beta density, and exit with an error if alpha<=1 | beta <=1. Note that the uniform density results with $\alpha = \beta = 1$, which does not have a unique HDR with content c("$0 < p <\n$", "$1$"). With shape parameters $\alpha<1$ and $\beta>1$ (or vice-versa, respectively), the Beta density is infinite at 0 (or 1, respectively), but still integrates to one, and so a HDR is still well-defined (but not implemented here, at least not yet). Similarly, with $0 < \alpha, \beta < 1$ the Beta density is infinite at both 0 and 1, but integrates to one, and again a HDR of content $p<1$ is well-defined in this case, but will be a set of two disjoint intervals (again, at present, this function does not cover this case).

Returns

If the numerical optimization is successful an vector of length 2, containing $v$ and $w$, defined above. If the optimization fails for whatever reason, a vector of NAs is returned.

The function will also produce a plot of the density with area under the density supported by the HDR shaded, if the user calls the function with plot=TRUE; the plot will appear on the current graphics device.

Debugging messages are printed to the console if the debug

logical flag is set to TRUE.

Author(s)

Simon Jackman simon.jackman@sydney.edu.au . Thanks to John Bullock who discovered a bug in an earlier version.

See Also

pbeta, qbeta, dbeta, uniroot

Examples

betaHPD(4,5)
betaHPD(2,120)
betaHPD(120,45,p=.75,xlim=c(0,1))