betaHPD(alpha,beta,p=.95,plot=FALSE,xlim=NULL,debug=FALSE)
TRUE
then plot the density and
show the HDRNULL
, in which
case the function will confine plotting to where the density is
non-neglibleTRUE
produce messages to the
consoleNAs
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
.
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 $p1$>=1>
pbeta
, qbeta
,
dbeta
, uniroot
betaHPD(4,5)
betaHPD(2,120)
betaHPD(120,45,p=.75,xlim=c(0,1))
Run the code above in your browser using DataLab