Estimating quantiles for posterior predictive distribution

Question

Let's say I have a parametrised distribution $X|\theta$ and a sample from the posterior distribution of $\theta$, $\hat{\theta}_1,...,\hat{\theta}_k$. If I was interested in the CDF of the predictive distribution evaluated at $x_0$, $F_X(x_0)$, I could use this sample to measure uncertainty around a central estimate given by

$$ \hat{F}_X(x_0) = \frac{1}{k} \sum_{j=1}^k F_{X|\theta}(x_0|\hat{\theta}_j) $$

and that would be an unbiased estimator since

$$ F_X(x_0) = \int_\Theta F_{X|\theta}(x_0|\theta)\cdot f(\theta) d\theta = \mathbb{E}_\theta[F_{X|\theta}(x_0|\theta)] $$

However if I was interested instead in the quantile of the predictive distribution at some $p \in (0,1)$, $F^{-1}_X(p)$, I don't think I can do the same because $F^{-1}_X(p) \neq \mathbb{E}_\theta[F^{-1}_{X|\theta}(p|\theta)]$. I can numerically find $\hat{q}$ such that

$$ \frac{1}{k} \sum_{j=1}^k F_{X|\theta}(\hat{q}|\hat{\theta}_j) = p $$

but then I wouldn't have any measure of variability around this estimate. What would be the proper way to do this?

Answer 1

To the Bayesian, there is no uncertainty around the predictive distribution (no hat). It is not something to be estimated. It has already incorporated all the uncertainty around the model parameters. Sampling from the posterior is simply a method of numerical integration, a convenience. The number of samples $k$ from the posterior would be very very large so that this numerical integration is close to exact. You can supply a percentage and identify the corresponding percentile from the predictive distribution. Likewise, you can supply a quantile and identify the corresponding percentage. They do not use hats on their samples from the posterior either.