| betaHPD {pscl} | R Documentation |
compute and optionally plot beta HDRs
Description
Compute and optionally plot highest density regions for the Beta distribution.
Usage
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 |
xlim |
numeric vector of length 2, the limits of the density's
support to show when plotting; the default is |
debug |
logical flag, if |
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 p \in
(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 0 < p <
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).
Value
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
Examples
betaHPD(4,5)
betaHPD(2,120)
betaHPD(120,45,p=.75,xlim=c(0,1))