Please, anybody could explain the steps to compute the highest posterior density (HPD) interval, when the posterior distribution is known? For instance, when the posterior distribution is Beta distributed.
When the posteriori distribution is simulated, the Chen-Shao algorithm can be used to estimate the HPD interval.
An HPD region is defined as$$\mathfrak{h}_\tau \stackrel{\text{def}}{=} \{\theta;\ \pi(\theta|x)>\tau\}$$and it is an interval only when the parameter is unidimensional and the posterior is unimodal. Assuming this is the case and the posterior $\pi(\cdot|x)$ is available up to a multiplicative constant, finding an HPD interval consists in solving in $\theta$ the equation$$\pi(\theta|x)=\tau$$Since in most situations a coverage of $\alpha$ is requested, a second computational step consists in associating a coverage $\alpha(\tau)$ with the bound $\tau$, as in $$\int_{\{\theta;\ \pi(\theta|x)>\tau\}} \pi(\theta|x)\,\text{d}\theta = \alpha(\tau)a$$followed by the inversion of the function $\alpha(\tau)$ to find the value of $\tau$ guaranteeing the proper coverage.
In the case $\pi(\theta|x)$ is the density of a Beta $B(\delta,\beta)$distribution, the first step requires solving $$\theta^{\delta-1}(1-\theta)^{\beta-1}=\tau$$ which usually has no analytic solution [unless there exists $\gamma$ such that both $\gamma(\delta-1)$ and $\gamma(\beta-1)$ are integers. Hence a numerical resolution of the equation is required. For each $\tau$, the coverage $\alpha(\tau)$ can then be derived by calling the Beta cdf.
