Since $$\mathbb E[Z|Z'=z',Y=y]=\int_{\mathbb R^4} z f(z|z',y)\,\text dz$$ a Monte Carlo approach requires simulating from the conditional distribution of $Z$ given $Z',Y$. Since $$f(z|z',y)\propto f(z,z',y)=\underbrace{f(z,z'|\Sigma)}_\text{Gaussian} \times f(y|z,\sigma)\tag{1}$$ simulation is feasible (via, e.g., MCMC) if the rhs of (1) is available in closed form. Now, $$y|z,\sigma\sim p(z)\mathcal N(z_1,\sigma^2)+(1-p(z))\mathcal N(z_2,\sigma^2)$$ with $$p(z)=\dfrac{e^{z_3}}{e^{z_3}+e^{z_4}}$$ therefore $f(y|z,\sigma)$ is clearly available.

Since $$\mathbb E[Z|Z'=z',Y=y]=\int_{\mathbb R^4} z f(z|z',y)\,\text dz$$ a Monte Carlo approach requires simulating from the conditional distribution of $Z$ given $Z',Y$. Since $$f(z|z',y)\propto f(z,z',y)=\underbrace{f(z,z'|\Sigma)}_\text{Gaussian} \times f(y|z,\sigma)\tag{1}$$ simulation is feasible (via, e.g., MCMC) if the rhs of (1) is available in closed form. Now, $$y|z,\sigma\sim p(z)\mathcal N(z_1,\sigma^2)+(1-p(z))\mathcal N(z_2,\sigma^2)$$ with $$p(z)=\dfrac{e^{z_3}}{e^{z_3}+e^{z_4}}$$ therefore $f(y|z,\sigma)$ is clearly available.

This development leads to $$f(z|z',y)\propto p(z)f(z|z',\Sigma)\varphi(\{y-z_1\}\sigma^{-1})+(1-p(z)) f(z|z',\Sigma)\varphi(\{y-z_2\}\sigma^{-1})$$