I'm going to attempt an intuitive explanation in a similar style to the linked question:

---

>First of all why not simply exclude the uncertain data all together?

There might be several reasons for this. One might be interested in computing quantities that require that value be present, for example a regression or classification model. In these cases, "imputing" $z$ in this manner is more principled than picking an arbitrary value, particularly if the values of $x$ and $y$ might feasibly have something to do with why $z$ is missing.

This is (I believe) what the linked question in your linked question was discussing

---

>Now in order to use this data and make GP predictions on some new location (xnew,ynew) to obtain znew we must take into account all the uncertainty related to z∗, which now has pdf due to its uncertainty, and integrate it out!

Just to clarify a potential source of confusion here: the uncertainty that we "integrate out" is the *predictive* uncertainty associated with predicting $z_{\mathrm{new}}$; that is, for new points $(x_\mathrm{new}, y_\mathrm{new})$, we predict the unobserved $z_\mathrm{new}$, which *due to this prediction*, is uncertain.

Any uncertainty related to $z$ (that is; the observations we *do* have, *not* ones we might want to predict) should be included in the GP we use to model the *unobserved* $z_\mathrm{new}$.

To summarise:

 - Uncertainty in *measured* $z$, associated with noise or imprecise measurement, is included in the GP model (specifically, in the kernel).
 - Uncertainty in $z_\mathrm{new}$ is associated with predicting (inferring) an unobserved quantity. This prediction uncertainty takes into account any measurement uncertainty modelled by the GP; that is, measurement uncertainty is "included" in the model when formulating the predictive distribution.
 - The uncertainty we "integrate over" is the *predictive* uncertainty (which, given the previous point, "accounts for" any measurement uncertainty since we included it in the model used to compute the predictive distribution). You are correct that this is equivalent to taking the expectation over of $z_{\mathrm{new}}$.

---

>How does integrating or summing up all the values of uncertain $z^∗$'s  help us as described in the previous post?

Essentially what's happening here is we're taking into account the fact that $z^*$ has to be *predicted* and thus is fundamentally uncertain. This integration is an attempt at accounting for this prediction uncertainty in whatever it is that's being predicted. This is equivalent to taking the expectation over the predictive distribution of $z_{\mathrm{new}}$.

---

Your equations look fine to me at a glance, but someone else may like to confirm if it makes sense to write them like that.

Essentially to *marginalise* (that is, *integrate out*) uncertainty in $z^*_{n+1}, \dots, z^*_{n+k}$ you are computing a multi-dimensional integral. 

In practice, you might benefit from this by doing some kind of Monte-Carlo simulation: you could draw samples of $z^*_{n+1}, \dots, z^*_{n+k}$ from the predictive distribution of the GP (important, since the distributions of $z^*_{n+1}, \dots, z^*_{n+k}$ will be correlated!), and then compute an approximation to the expectation of $f$ using those samples. In other words, you need to sample all the $z^*$'s you wish to take the expectation over together to respect the fact that they are correlated (via the GP predictive distribution).