Let $y$ follow $f(y|\theta,x)$ with $\theta$ parameters with prior $\pi(\theta)$ and $x$ covariates with distribution $p(x)$. Let $p(\theta| D)$ be the posterior distribution of the parameters, where data $D=(y,x)$. For a new value $\tilde{x}$, the posterior predictive distribution is $$f(\tilde{y}|\tilde{x},D) = \int_\theta f(\tilde{y}|\theta, \tilde{x}) p(\theta| D) d \theta$$

This is the setting of the ppd in regression analysis, for example. Now what I am interested in is posterior inference about the cumulative distribution of $\tilde{y}$: $$F(\tilde{y} | \theta ) = \int_x F(\tilde{y}|\theta,\tilde{x}) p(\tilde{x}) d\tilde{x}. $$

One way to sample from the posterior predictive distribution of $F$ evaluated at $y$ seems to be $$ \int_x \int_\theta F(\tilde{y}|\tilde{x},\theta) p(\theta|D) p(\tilde{x}) d\theta d\tilde{x}.$$ So if we can sample from $p(\theta|D)$, and we know / can calculate $F$ we would

  1. Sample $\theta_k \sim p(\theta|D) $
  2. Sample $\tilde{x}_k \sim p(\tilde{x})$
  3. Evaluate $F(\tilde{y}|\tilde{x}_k,\theta_k)$ to get $\hat{F}_K$
  4. Repeat many times. Then $\hat{F}_1,...,\hat{F}_K$ are the posterior samples of $F$ at $\tilde{y}$ integrated across $\theta, x$
  5. To get the posterior mean estimate of $F(\tilde{y}|\theta)$, for example, we could take the average of $F_k$ across $k$

Although this procedure seems valid, an alternative seems to be based on

$$f(\tilde{y}|D) = \int_x \int_\theta f(\tilde{y}|\theta, \tilde{x}) p(\theta| D) p(\tilde{x}) d \theta d\tilde{x},$$ yielding procedure

  1. Sample $\theta_k \sim p(\theta|D) $
  2. Sample $\tilde{x}_k \sim p(\tilde{x})$
  3. Sample $\tilde{y}_k \sim f(\tilde{y}|\theta_k,\tilde{x}_k) $. These are the samples from the posterior predictive distribution integrated across $\theta, \tilde{x}$
  4. Evaluate the empirical cumulative distribution of the samples $y_1,...,y_K$

At first I thought the two procedures are equivalent because they both seem to yield $F(\tilde{y}|D)$. However they yield different results. In particular the second approach does not seem to yield valid inference on $F(\tilde{y}|\theta)$. Am I right, and if so, why?

Edit After the reply by Xi'an and some more own thought I have come up with a third, and potentially superior approach which is a modification of the non-parametric sampling algorithm 2 above. Would welcome thoughts.

Let $q_\alpha(\tilde{y}| \theta)$ denote the $\alpha$ quantile of $\tilde{y}$ given $\theta$. Then

  1. Sample $\theta_k \sim p(\theta|D) $
  2. Now we estimate $q_\alpha(\tilde{y}| \theta_k)$ by monte carlo. For this we sample J times, j=1,...,J, from $\tilde{x}_j \sim p(\tilde{x})$.
  3. Sample $\tilde{y}_j \sim f(\tilde{y}|\theta_k,\tilde{x}_j), j=1,...,J $.
  4. Estimate $q_\alpha(\tilde{y}| \theta_k)$ as empirical quantile of $\tilde{y}_1,...,\tilde{y}_J$

We now have the posterior distribution of $q_\alpha(\tilde{y}| \theta_1),...,q_\alpha(\tilde{y}| \theta_K)$. Let the posterior median be $\bar{q}$, then $\alpha$ is a posterior estimate for $F(\tilde{y}|\theta)$ at $\tilde{y}=\bar{q}$. A credible interval can also be obtained.


1 Answer 1


The two quantities $$ F(\tilde y|D) = \int_\Theta \int_{\mathcal X} F(\tilde{y}|\tilde{x},\theta) p(\theta|D) p(\tilde{x})\, d\tilde{x}\, d \theta \tag{1} $$ and $$ f(\tilde{y}|D) = \int_\Theta \int_{\mathcal X} f(\tilde{y}|\theta, \tilde{x}) p(\theta| D) p(\tilde{x}) \, d\tilde{x}\, d \theta \tag{2} $$ are identical in that one is the cdf and the other one the pdf. The difference in the simulation schemes is that (1) produces a (parametric) Rao-Blackwellised approximation to the marginal cdf, $$ \frac{1}{N}\sum_{i=1}^n F(\tilde y|\theta_i,\tilde x_i) $$ while (2) leads to a (non-parametric) empirical approximation, hence less efficient.

Note that both aim at approximating the marginal predictive distribution on $\tilde Y$, which does not bring information about $F(\tilde y|\theta)$ since $\theta$ is integrated out. To find the posterior of the quantity $F(\tilde y|\theta)$ (for a fixed $\tilde y$) simply requires to find the push-forward transform of $p(\theta|D)$ by $F(\tilde y|\cdot)$. In that sense, version 1 of the algorithm does produce a Monte Carlo sample from that distribution.

  • $\begingroup$ It seems that with algorithm 2 (push forward) I can also get a credible interval by evaluating the percentiles of $F_1,...,F_k$ at $\tilde{y}$, whereas with algorithm 1 I only obtain something like an approximation of the cumulative distribution of $\tilde{Y}$.. $\endgroup$
    – tomka
    Nov 18, 2020 at 21:21
  • $\begingroup$ And on the same note, I could perhaps also use the median of $F_1,...,F_k$ at $\tilde{y}$ instead of the mean; in other words I can draw posterior inference on $F$ as usual. Is that correct? $\endgroup$
    – tomka
    Nov 18, 2020 at 22:36
  • 1
    $\begingroup$ The sample of $F_i$'s produced by the first algorithm in your question is an iid sample from the posterior distribution of $F(\tilde y|\theta)$, hence can be used in any way you wish. The second algorithm in your question does not produce information about $F(\tilde y|\theta)$. $\endgroup$
    – Xi'an
    Nov 19, 2020 at 6:36
  • $\begingroup$ Thanks. I have come up with another method to estimate the marginal cdf, please see above. Would welcome comments. This also works when $F$ does not have a known parametrical form and gives a posterior interval at each quantile. $\endgroup$
    – tomka
    Nov 19, 2020 at 14:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.