# compute and optionally plot beta HDRs Compute and optionally plot highest density regions for the Beta distribution. ```r 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 ```r betaHPD(4,5) betaHPD(2,120) betaHPD(120,45,p=.75,xlim=c(0,1)) ```

betaHPD function