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