Approximate marginal posterior distribution via sampling

Question

I'm a beginner in Bayesian inference and I have some confusion about posterior distributions, in particular sampling from it vs. approximating its value at some given point.

Suppose I have a model with parameters $\theta = (\theta_1, \theta_2)$, and observed data $y$. Suppose I can sample from the posterior distribution $p(\theta| y)$, and let $\mathcal S = \{\theta^{(s)} = (\theta_1^{(s)}, \theta_2^{(s)}):s=1,\ldots, S\}$ be samples.

Then, regarding the posterior predictive distribution $$p(\widetilde y| y) = \int p(\widetilde y|\theta)p(\theta|y)d\theta$$ I can (taking from Hoff's Bayesian Statistical Method Section 4.3):

1. Sample from $p(\widetilde y|y)$: starting from samples $\mathcal S$ I draw $\widetilde y^{(s)} \sim p(\widetilde y|\theta^{(s)})$ obtaining samples $\{(\widetilde y^{(s)}, \theta^{(s)})\}$ from the joint posterior $p(\widetilde y, \theta|y)$ and then I just ignore the $\theta$'s obtaining a sample $\{\widetilde y^{(s)}\}$ from the posterior predictive distribution.
2. Approximate the value of $p(\widetilde y|y)$ for any given $\widetilde y$ by Monte Carlo approximation: since $$p(\widetilde y|y) = \int p(\widetilde y|\theta)p(\theta|y)d\theta = \mathrm E_{p(\theta|y)}[p(\widetilde y|\theta)] \approx \frac{1}{S}\sum_{s=1}^Sp(\widetilde y| \theta^{(s)})$$

My question is if a similar argument allows to approximate the marginal posterior distribution \begin{equation} p(\theta_1|y) = \int p(\theta_1|\theta_2, y)p(\theta_2|y)d\theta_2\tag{1} \end{equation} How do I approximate the value of $p(\theta_1|y)$ for a given value of $\theta_1$?

Reasoning as above:

• Sampling from $p(\theta_2|y)$: take $\mathcal S$ and just ignore the $\theta_1$'s, obtaining a sample $\{\theta_2^{(s)}\}$ from the posterior marginal $p(\theta_2|y)$
• Approximating $p(\theta_1|y)$ for a given $\theta_1$: my guess is that I can consider $p(\theta_1|y) = \mathrm E_{p(\theta_2|y)}[p(\theta_1|\theta_2,y)]$ so that again by Monte Carlo simulation I can approximate $p(\theta_1|y) \approx \frac{1}{S}\sum_{s=1}^Sp(\theta_1|\theta_2^{(s)}, y)$. Is this argument correct?

It seems to me that this case is analogous to the above for $p(\widetilde y|y)$ however I could not find this clearly stated in books so I wonder if I am missing something basic. For example in BDA3, Section 3.1 they say:

We rarely evaluate integral (1) explicitly [...]. Posterior distributions [i.e., $p(\theta_1, \theta_2|y)$] can be computed by marginal and conditional simulation, first drawing $\theta_2$ from its marginal posterior distribution [i.e., $p(\theta_2| y)$] and then $\theta_1$ from its conditional posterior distribution [i.e., $p(\theta_1|\theta_2, y)$] given the drawn value of $\theta_2$. In this way the integration embodied in (1) is performed indirectly

I don't understand in which sense this amounts to computing the integral (1) for some given value of $\theta_1$, it looks like a way to sample from the distribution it defines. I would greatly appreciate any reference (I guess it could appear as well when integrating nuisance parameters).

Answer 1

Using $$p(\theta_1|y) \approx \frac{1}{S}\sum_{s=1}^Sp(\theta_1|\theta_2^{(s)}, y)$$ as an approximation to the marginal posterior density is called a Rao-Blackwellisation in the MCMC literature, started in Gelfand & Smith (1990) foundational¹ paper, appropriately entitled "Sampling-Based Approaches to Calculating Marginal Densities". This density estimator is converging at a parametric speed, $\text{O}(1/\sqrt S)$.

When proposing the Gibbs sampling algorithm as a way to simulating from marginal densities and (hence) posterior distributions, Gelfand and Smith (1990) explicitly relate to the Rao-Blackwell theorem, as shown by the following quote²

...we consider the problem of calculating a final form of marginal density from the final sample produced by either the substitution or Gibbs sampling algorithms. Since for any estimated marginal the corresponding full conditional has been assumed available, efficient inference about the marginal should clearly be based on using this full conditional distribution. In the simplest case of two variables, this implies that $[X \mid Y]$ and the $y^{(i)}_j$'s $(j = 1, \ldots, m)$ should be used to make inferences about $[X]$, rather than imputing $X^{(i)}_j$ $(j = 1, \ldots, nm)$ and basing inference on these $X^{(i)}_j$'s. Intuitively, this follows, because to estimate $[X]$ using the $x^{(i)}_j$'s requires a kernel density estimate. Such an estimate ignores the known form $[X \mid Y]$ that is mixed to obtain $[X]$. The formal argument is essentially based on the Rao--Blackwell theorem. We sketch a proof in the context of the density estimator itself. If $X$ is a continuous p-dimensional random variable, consider any kernel density estimator of $[X]$ based on the $X^{(i)}_j$'s (e.g., see Devroye and Györfi, 1985) evaluated at $x_0$:$$\Delta^{(i)}_{x_0} = (1/h_m^p) \sum_{j=1}^m K[(X_0 - X^{(i)}_j)/h_m]$$ say, where $K$ is a bounded density on $\mathbb R^p$ and the sequence $\{h_m\}$ is such that as $m \rightarrow \infty$, $h_m \rightarrow 0$, whereas $mh_m \rightarrow \infty$. To simplify notation, set$$Q_{m,x_0}(X) = (1/h_m^p) K[(X - X^{(i)}_j)/h_m]$$ so that$$\Delta^{(i)}_{x_0} = (1/m) \sum_{j=1}^m Q_{m,x_0}(X_j^{(i)})$$ Define$$\gamma_{x_0}^i = (1/m) \sum_{j=1}^m \mathbb E[Q_{m,x_0}(X)\mid Y^{(i)}_j]$$ By our earlier theory, both $\Delta^{(i)}_{x_0}$ and $\gamma_{x_0}^i$ have the same expectation. By the Rao--Blackwell theorem,$$\text{var}\, \mathbb E[Q_{m,x_0}(X)] \mid Y) \le \text{var}\, Q_{m,x_0}(X)$$ and hence$$\text{MSE}(\gamma_{x_0}^i)\le\text{MSE}(\Delta^{(i)}_{x_0})$$ where MSE denotes the mean squared error of the estimate of $[X]$.

---

¹ "foundational" as it launched the MCMC revolution.

² The text has been retyped by me and may hence contains typos. The notations are those introduced by Gelfand and Smith (1990) and used for a while in the literature with $[X \mid Y]$ denoting the conditional density of $X$ given $Y$. The double indexation of the sequence is explained in the quote.