What is the correct procedure for posterior inference about the cumulative distribution $F(y|\theta)$? 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

*

*Sample $\theta_k \sim  p(\theta|D) $

*Sample $\tilde{x}_k \sim p(\tilde{x})$

*Evaluate $F(\tilde{y}|\tilde{x}_k,\theta_k)$ to get $\hat{F}_K$

*Repeat many times. Then $\hat{F}_1,...,\hat{F}_K$ are the posterior samples of $F$ at $\tilde{y}$ integrated across $\theta, x$

*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

*

*Sample $\theta_k \sim  p(\theta|D) $

*Sample $\tilde{x}_k \sim p(\tilde{x})$

*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}$

*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

*

*Sample $\theta_k \sim  p(\theta|D) $

*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})$.

*Sample $\tilde{y}_j \sim f(\tilde{y}|\theta_k,\tilde{x}_j), j=1,...,J $.

*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.
 A: 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.
