How to use Monte Carlo simulation to get the conditional mean Given the following assumptions:

*

*$Z,Z'\in\mathbb{R}^4$ where $(Z,Z')\sim N(0,\Sigma)$, for some known $\Sigma\in\mathbb{R}^{8\times 8}$.

*$Y=f(Z,u,\epsilon)=Z_1\boldsymbol{1}\Big[u<\frac{\exp(Z_3)}{\exp(Z_3)+\exp(Z_4)}\Big]+Z_2\boldsymbol{1}\Big[\frac{\exp(Z_3)}{\exp(Z_3)+\exp(Z_4)}<u\leq 1\Big]+\epsilon$, where $Y\in\mathbb{R}$, $Z=(Z_1,Z_2,Z_3,Z_4)$, $u\sim\text{Uniform}[0,1]$, $\epsilon\sim N(0,\sigma^2)$ for some known $\sigma$, and $\boldsymbol{1}[\cdot]$ is the indicator function. We can think of $Y$ as just some complicated function that depends on $Z$.

Question: Given that we have access to $Z^k$, $Y$, $\Sigma$, and $\sigma$, can we use Monte Carlo to get an approximation for $\mathbb{E}[Z|Z',Y]$? If so, then how can we do it? In other words, given $Z'$ and $Y$, how can we sample $Z$? I've tried thinking about this but wasn't able to come up with a solution.
Some thoughts: If $Y$ is a Gaussian (e.g., $Y=Z_1+Z_2+\epsilon$) then we can get $\mathbb{E}[Z|Z',Y]$ explicitly. But the problem now is that $Y$ is not a Gaussian, requiring me to resort to some form of Monte Carlo simulation. Essentially, I just want to obtain some form of approximation of $\mathbb{E}[Z|Z',Y]$ so any suggestion (even if it has nothing to do with Monte Carlo) is welcomed.
Idea: The best I can think of is this very naive method:

*

*Sample $z$ from the distribution of $Z$, sample $u$ from $\text{Uniform}[0,1]$, and sample $\epsilon$ from $N(0,\sigma^2)$.

*Set $\zeta=0.1$ and select $z$ iff $|f(Z,u,\epsilon)-Y|<\zeta$.

*Take the average of all the selected $z$'s and return it as an approximate of $\mathbb{E}[Z|Z',Y]$.

Can we do any better?

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