7.1 General Result

We state and prove a general result by abstracting a bit as follows. Note that the posterior density for a real valued parameter is a density on \(\mathbb{R}\). Also, note that the length of an arbitrary set \(C\) can be defined as \[ \int I_C(\theta) {\rm d}\theta, \] where \(I_C(\cdot)\) is the indicator function for the set \(C\), i.e. \(1\) on \(C\) and \(0\) otherwise.

Theorem Let \(f\) denote a density on \(\mathbb{R}\). Then for \(\alpha\in (0,1)\), a set \(C\) with the shortest length among those with \(f\) probability equal to \(1-\alpha\) equals \[ C_\alpha:=\left\{f>c_\alpha \right\} \cup E, \] from some \(c_\alpha>0\) and \(E\subseteq \left\{f=c_\alpha \right\}\). \(\blacksquare\)

We note that in applications, the set \(\left\{f=c \right\}\), for any \(c>0\), would contain atmost a finite many elements; in such cases, the set \(E\) can be taken to be the empty set, i.e. \[ C_\alpha:=\left\{f>c_\alpha \right\}, \] from some \(c_\alpha>0\). An extreme case where \(E\) is not empty for any \(\alpha\in(0,1)\) is when \(f\) is a uniform density - but a uniform density does not arise as a posterior density.

Proof of the theorem: We will break the proof into the following two steps. The second step is the more important one from a practical insight perspective; the first is for those motivated to strive towards a full understanding.
The existence of \(c_\alpha\): Consider the set \(D_c\) defined for \(c\geq 0\) as \[ D_c:=\left\{f>c \right\}. \] We note that \({\rm Pr}\left(D_c\right)\) under \(f\) is a right continuous function of \(c\) as it is the survival function of the density (as the random variable). The points of discontinuity are the level sets of \(f\) which have a positive probability under \(f\), i.e. values of \(c\) such that \[ \int I_{\left\{\eta:f(\eta)=c\right\}}(\theta) f(\theta) {\rm d}\theta >0. \] To see this let we note that since \(\psi\) defined by \(\psi(c):={\rm Pr}\left(D_c\right)\) is a decreasing function of \(c\) with \(\psi(0)=1\) \(\lim_{c\rightarrow\infty}\psi(c)=0\), there either exists a \(c_\alpha\) such that \(C_\alpha:=\left\{f>c_\alpha \right\}\) satisfies \[ \int I_{C_\alpha}(\theta) f(\theta) {\rm d}\theta =1-\alpha, \] or \[ \psi(c_\alpha-)\geq(1-\alpha)>\psi(c_\alpha). \] In the latter case, using a technical result from measure theory (that the Lebesgue measure is nonatomic; I teach this in my Fall PhD level course on probability theory), it can be shown that a subset \(E\) of \(\{\eta:f(\eta)=c_\alpha\}\) can be carved such that \[ C_\alpha:=\left\{f>c_\alpha \right\} \cup E, \] has \(f\) probability \(1-\alpha\).

The Optimality of \(C_\alpha\): Observe that since \(c_\alpha>0\), \[ \int I_{C_\alpha}(\theta) {\rm d}\theta\leq \int I_{C_\alpha}(\theta) \frac{f(\theta)}{c_\alpha}{\rm d}\theta\leq \frac{1-\alpha}{c_\alpha}<\infty. \] We begin by showing that for any set \(C\) with the same length as \(C_\alpha\) has \(f\) probability atmost \((1-\alpha)\), i.e. \({\rm Pr}\left(C_\alpha\right)\geq {\rm Pr}\left(C\right)\). Note that this is equivalent to showing that \({\rm Pr}\left(C_\alpha\cap {C}^{\mathsf{c}}\right)\geq {\rm Pr}\left(C\cap {C_\alpha}^{\mathsf{c}}\right)\). This is true since the lengths of \(C_\alpha\cap {C}^{\mathsf{c}}\) equals that of \(C\cap {C_\alpha}^{\mathsf{c}}\) and moreover, \[ C_\alpha\cap {C}^{\mathsf{c}}\subseteq \{\eta:f(\eta)\geq c_\alpha\} \quad \hbox{and} \quad C\cap {C_\alpha}^{\mathsf{c}}\subseteq \{\eta:f(\eta)\leq c_\alpha\}. \] More clearly, \[ \begin{aligned} {\rm Pr}\left(C_\alpha\cap {C}^{\mathsf{c}}\right)&=\int_{C_\alpha\cap {C}^{\mathsf{c}}} f(\theta){\rm d}\theta\\ &\geq c_\alpha \int_{C_\alpha\cap {C}^{\mathsf{c}}} {\rm d}\theta\\ &= c_\alpha \int_{C\cap {C_\alpha}^{\mathsf{c}}} {\rm d}\theta\\ &\geq \int_{C\cap {C_\alpha}^{\mathsf{c}}} f(\theta){\rm d}\theta={\rm Pr}\left(C\cap {C_\alpha}^{\mathsf{c}}\right).\\ \end{aligned} \] This implies that any other candidate \((1-\alpha)\) level credible set \(C\) has atleast the same length as \(C_\alpha\). This is so as other wise \(C\) can be augmented to increase its length to that of \(C_\alpha\) (using if need be non-atomicity of the Lebesgue measure), which will result in a post-augmentation probability under \(f\) of \(C\) being strictly greter than \((1-\alpha)\) - a contradiction to the above. \(\blacksquare\)

Note: It can be shown that any shortest length set \(C\) is essentially of the above form with a bit more effort, but since that is not practically relevant we will leave it for the motivated reader to either prove it themselves or discuss it during office hours.

7.2 Highest Posterior Density (HPD) Credible Set

Note that the above general result implies that the above \(C_\alpha\) is a shortest length \((1-\alpha)\) level credible set for a parameter \(\theta\), constructed with the posterior density of \(\theta\), denoted by \(\pi(\cdot|\tilde{x})\), replacing \(f\) above.

Now we turn to defining the above \(C_\alpha\) in practice via computing the \(c_\alpha\) - the cases we will look at will be such that the level sets of the posterior would be sets with finitely many elements, so the \(E\) above could be taken to be the empty set.

Note that mathematically the HPD set is easily defined if the posterior density is non-increasing from the left-end point of its support - in this case it will be an interval from this left-end point of its support to its \((1-\alpha)\)-th quantile. In the example above this observation can be used to construct the HPD set. But because of the Bayesian CLT (see Theorem 15.16 of the text), which is called the Bernstein von-Mises theorem, we see that but for these pathological cases the posterior will approach normal and hence its density is unlikely to be non-increasing. So this mathematically interesting case is not statistically or practically of much use.

On the other hand, asymptotic normality suggests that a good class of posteriors to look at will be those which are unimodal. But before we do this, we comment that if the posterior is a symmetric unimodal then the HPD will be constructed by solving for \(y\) such that \[ (1-\alpha)=F(\hbox{posterior mode}+y)-F(\hbox{posterior mode}-y)=2*F(mode+y)-1, \] where \(F\) denotes the posterior. Note that suc a \(y\) is given by \[ y=F^{-1}(1-\alpha/2)-\hbox{posterior mode}, \] which can be computed if \(F^{-1}\) is available (like qnorm) or by using a root finding function like uniroot.