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

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.

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.
• 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}$.. Nov 18 '20 at 21:21
• 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? Nov 18 '20 at 22:36
• 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)$. Nov 19 '20 at 6:36
• 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. Nov 19 '20 at 14:46