# How to optimally spread draws when calculating multiple expectations

Suppose we want to calculate some expectation:

$$E_YE_{X|Y}[f(X,Y)]$$

Suppose we want to approximate this using Monte Carlo simulation.

$$E_YE_{X|Y}[f(X,Y)] \approx \frac1{RS}\sum_{r=1}^R\sum_{s=1}^Sf(x^{r,s},y^r)$$

BUT suppose it is costly to draw samples from both distributions, so that we only can afford to draw a fixed number $K$.

How should we allocate $K$? Examples include $K/2$ draws to each distribution, or in the extreme, one draw in the outer and $K-1$ draws in the inner, vice versa etc.....

My intuition tells me that it will have to do with the variance/entropy of the distributions relative to each other. Suppose the outer one is a mass point, then the division of $K$ that minimizes MC error would be draw 1 of the $Y$ and draw $K-1$ of the $X|Y$.

Hopefully this was clear.

• Fixed it for you Aug 29 '16 at 13:04
• The "vice versa" and your comment to @Xi'ans answer seem to indicate that you consider it possible to draw the outer variable more times than the inner variable, but how could that make sense -- aren't all outers for which $0$ inners are drawn wasted? Aug 29 '16 at 15:58
• Fair enough, minimum one draw per outer I guess. Or you could think of programming it to save the draw I suppose Aug 29 '16 at 16:18
• @robertevansanders Please confirm whether the interpretation of your question in the first two sentences of Xi'ans answer is correct Aug 29 '16 at 19:38
• Like you said, yeah, but switch y and x Aug 29 '16 at 19:46

This is a very interesting question with little documentation in the Monte Carlo literature, except in connection with stratification and Rao-Blackwellisation. This is possibly due to the fact that the computations of the expected conditional variance and of the variance of the conditional expectation are rarely feasible.

First, let us assume you run $R$ simulations from $\pi_X$, $x_1,\ldots,x_R$ and for each simulated $x_r$, you run $S$ simulations from $\pi_{Y|X=x_r}$, $y_{1r},\ldots,y_{sr}$. Your Monte Carlo estimate is then $$\delta(R,S)=\frac{1}{RS}\sum_{r=1}^R\sum_{s=1}^S f(x_r,y_{rs})$$ The variance of this estimate is decomposed as follows \begin{align*} \text{var} \{\delta(R,S)\} &= \frac{1}{R^2S^2} R\text{var} \left\{\sum_{s=1}^S f(x_r,y_{rs})\right\}\\ &= \frac{1}{RS^2} \text{var}_X\mathbb{E}_{Y|X}\left\{\sum_{s=1}^S f(x_r,y_{rs})\big|x_r\right\}+\frac{1}{RS^2}\mathbb{E}_{X}\text{var}_{Y|X} \left\{\sum_{s=1}^S f(x_r,y_{rs})\big|x_r\right\}\\ &=\frac{1}{RS^2} \text{var}_X\{ S \mathbb{E}_{Y|X}[f(x_r,Y)|x_r]\}+ \frac{1}{RS^2} \mathbb{E}_{X}[S\text{var}_{Y|X}\{f(x_r,Y)|x_r\}]\\ &=\frac{1}{R} \text{var}_X\{\mathbb{E}_{Y|X}[f(x_r,Y)|x_r]\}+ \frac{1}{RS} \mathbb{E}_{X}[\text{var}_{Y|X}\{f(x_r,Y)|x_r\}]\\ &\stackrel{K=RS}{=}\frac{1}{R}\text{var}_X\{\mathbb{E}_{Y|X}[f(x_r,Y)|x_r]\}+ \frac{1}{K} \mathbb{E}_{X}[\text{var}_{Y|X}\{f(x_r,Y)|x_r\}] \end{align*} Therefore if one wants to minimise this variance the optimal choice is $R=K$. Implying that $S=1$. Except when the first variance term is null, in which case it does not matter. However, as discussed in the comments, the assumption $K=RS$ is unrealistic as it does not account for the production of one $x_r$ [or assumes this comes for free].

Now let us assume different simulation costs and the budget constraint $R+aRS=b$, meaning that the $y_{rs}$'s cost $a$ times more to simulate than the $x_r$'s. The above decomposition of the variance is then $$\frac{1}{R}\text{var}_X\{\mathbb{E}_{Y|X}[f(x_r,Y)|x_r]\}+ \frac{1}{R(b-R)/aR} \mathbb{E}_{X}[\text{var}_{Y|X}\{f(x_r,Y)|x_r\}]$$ which can be minimised in $R$ as $$R^*=b\big/1+\{a\mathbb{E}_{X}[\text{var}_{Y|X}\{f(x_r,Y)|x_r\}/\text{var}_X\{\mathbb{E}_{Y|X}[f(x_r,Y)|x_r]\}\}^{1/2}$$ [the closest integer under the constraints $R\ge 1$ and $S\ge 1$], except when the first variance is equal to zero, in which case $R=1$. When $\mathbb{E}_{X}[\text{var}_{Y|X}\{f(x_r,Y)|x_r\}]=0$, the minimum variance corresponds to a maximum $R$, which leads to $S=1$ in the current formalism.

Note also that this solution should be compared with the symmetric solution when the inner integral is in $X$ given $Y$ and the outer integral is against the marginal in $Y$ (assuming the simulations are also feasible in this order).

An interesting extension to the question would be to consider a different number of simulations $S(x_r)$ for each simulated $x_r$, depending on the value $\text{var}_{Y|X}\{f(x_r,Y)|x_r\}$.

• In the final conclusion, you seem to be assuming $K=RS$ but in the setting of the question $K=RS+R$ since draws of the outer variable are supposed to be counted, too. The result here says that if sampling the outer variable was free, of course one should sample a new outer for every inner. (Also, the role of $x$ and $y$ are switched here compared to the question but that of course does not matter). Aug 29 '16 at 8:14
• Yes but we can decide the value of $R$... Consider the degenerate setting where the outer var $X$ is a.s. constant. It's better to sample the constant once and $Y$ $K-1$ times, rather than the constant $K/2$ times and $Y$ $K/2$ times (which is what $S=1$ would imply)? Or am I completely misunderstanding the question? (I only now read the second sentence of your comment -- isn't the assumption stated in the question that they have the same cost) Aug 29 '16 at 8:23
• @Xi'an yes Kolkata is correct, your solution cannot generally hold. Suppose now that the inner variable has a degenerate distribution and the outer has meaningful variance, then you'd want to sample as few inner draws as possible Aug 29 '16 at 11:20
• I think your answer can't be right. Suppose the inside distribution is degenerate and the outside is large variance, how can S be 1 Aug 29 '16 at 21:22
• @robertevansanders: if the inner distribution is degenerate, $\text{var}_{Y|X}\{f(x_r,Y)|x_r\}=0$, hence $R^*=b$ and we pick the closest integer $R$ under the constraints $S\ge 1$ and $R(1+aS)\le b$, which means taking $S=1$ to make $R$ as close as possible to $b$. Aug 29 '16 at 21:53