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

1 |

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

`xlim` |
numeric vector of length 2, the limits of the density's
support to show when plotting; the default is |

`debug` |
logical flag, if |

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 *α > 1* and *β > 1*, the Beta density is
unimodal.

In general, suppose *θ \in Θ \subseteq R^k*
is a random variable with density *f(θ)*. A highest
density region (HDR) of *f(θ)* with content *p \in
(0,1]* is a set *\mathcal{Q} \subseteq Θ* with the
following properties:

*\int_\mathcal{Q} f(θ) dθ = p*

and

*f(θ) > f(θ^*) \, \forall\
θ \in \mathcal{Q},
θ^* \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(θ; α, β) dθ = p,*

where *f(θ; α, β)* is a Beta density with shape
parameters *α* and *β*.

In the special case of *α = β > 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 *α = β = 1*,
which does not have a unique HDR with content *0 < p <
1*. With shape parameters *α<1* and *β>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 < α, β < 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).

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`

.

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

`pbeta`

, `qbeta`

,
`dbeta`

, `uniroot`

1 2 3 |

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.

All documentation is copyright its authors; we didn't write any of that.