Compute and optionally plot highest density regions for the Beta distribution.
betaHPD(alpha,beta,p=.95,plot=FALSE,xlim=NULL,debug=FALSE)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.
scalar, first shape parameter of the Beta density. Must be greater than 1, see details
scalar, second shape parameter of the Beta density. Must be greater than 1, see details
scalar, content of HPD, must lie between 0 and 1
logical flag, if TRUE then plot the density and
show the HDR
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
logical flag, if TRUE produce messages to the
console
Simon Jackman simon.jackman@sydney.edu.au. Thanks to John Bullock who discovered a bug in an earlier version.
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).
betaHPD(4,5)
betaHPD(2,120)
betaHPD(120,45,p=.75,xlim=c(0,1))
Run the code above in your browser using DataLab