Strategies for sampling values from two-dimensional probability space Consider $I:X\times Y\to R$ a function of $x\in X$ and $y\in Y$. Assuming $x,y$ are statistically independent, we seek to estimate $E_{X,Y}[I(x,y)]$ using $N$ sampled pairs $(x,y)$. 
Are the two following strategies equivalent?
(1) Sample $N$ pairs $(x_i,y_i)$ and estimate $\hat{I}_N=\frac{1}{N}\sum_i I(x_i,y_i)$
(2) Choose $n$ values of $x_i$ and $m$ values $y_j$ with $nm=N$ and estimate $\hat{I}_N=\frac{1}{N}\sum_i \sum_j I(x_i,y_j)$
 A: Comparison between
$$\hat{I}^1_N=\frac{1}{N}\sum_i I(x_i,y_i)\quad\text{and}\quad\hat{I}^2_N=\frac{1}{N}\sum_i \sum_j I(x_i,y_j)$$is possible since
$$\text{var} \hat{I}^1_N=\frac{1}{N}\,\text{var}\{I(X,Y)\}$$while
\begin{align*}\text{var} \hat{I}^2_N 
&=\frac{1}{N^2}\,\text{var}\left\{\sum_{i,j}I(X_i,Y_j)\right\}\\
&=\frac{1}{N^2}\,\mathbb{E}\left[\text{var}\left\{\sum_{i,j}I(X_i,Y_j)\Big|X_1,\ldots,X_n\right\}\right]\qquad\text{[decomposition of variance]}\\
&\qquad+\frac{1}{N^2}\,\text{var}\left[\mathbb{E}\left\{\sum_{i,j}I(X_i,Y_j)\Big|X_1,\ldots,X_n\right\}\right]\\
&=\frac{1}{N^2}\,\mathbb{E}\left[m\,\text{var}^Y\left\{\sum_{i}I(X_i,Y)\right\}\right]\qquad\qquad\text{[the $Y_j$'s are i.i.d.]}\\
&\qquad+\frac{1}{N^2}\,\text{var}\left[m\,\mathbb{E}^Y\left\{\sum_{i}I(X_i,Y)\right\}\right]\\
&=\frac{m}{N^2}\,\mathbb{E}\left[\text{var}^Y\left\{\sum_{i}I(X_i,Y)\right\}\right]\\&\qquad +\frac{m^2n}{N^2}\,\text{var}^X\left[\mathbb{E}^Y\left\{I(X,Y)\right\}\right]\qquad\qquad\text{[the $X_i$'s are i.i.d.]}\\
&=\frac{mn}{N^2}\,\mathbb{E}^X\left[\text{var}^Y\left\{I(X,Y)\right\}\right]\qquad\quad\qquad\text{[variance of sum decomposition]}\\ 
&\qquad+\frac{mn(n-1)}{N^2}\,\mathbb{E}^{X_1,X_2}\left[\text{cov}^Y\left\{I(X_1,Y),I(X_2,Y)\right\}\right]\\&\qquad +\frac{m^2n}{N^2}\,\text{var}^X\left[\mathbb{E}^Y\left\{I(X,Y)\right\}\right]\\
&=\frac{mn}{N^2}\,\mathbb{E}^X\left[\text{var}^Y\left\{I(X,Y)\right\}\right]+\frac{mn}{N^2}\,\text{var}^X\left[\mathbb{E}^Y\left\{I(X,Y)\right\}\right]\qquad\qquad\text{[add &]}\\ 
&\qquad+\frac{mn(n-1)}{N^2}\,\mathbb{E}^{X_1,X_2}\left[\text{cov}^Y\left\{I(X_1,Y),I(X_2,Y)\right\}\right]\\
&\qquad +\frac{m^2n}{N^2}\,\text{var}^X\left[\mathbb{E}^Y\left\{I(X,Y)\right\}\right]-\frac{mn}{N^2}\,\text{var}^X\left[\mathbb{E}^Y\left\{I(X,Y)\right\}\right]\qquad\text{[substract]}\\
&=\underbrace{\frac{1}{N}\,\text{var}\left\{I(X,Y)\right\}}_{\text{var} (\hat{I}^1_N)}\qquad\qquad\qquad\qquad\qquad\text{[recomposition of variance]}\\
&\qquad\qquad+\frac{mn(n-1)}{N^2}\,\mathbb{E}^{X_1,X_2}\left[\text{cov}^Y\left\{I(X_1,Y),I(X_2,Y)\right\}\right]\\ 
&\qquad\qquad +\frac{m(m-1)n}{N^2}\,\text{var}^X\left[\mathbb{E}^Y\left\{I(X,Y)\right\}\right]\\
&=\text{var} (\hat{I}^1_N)+\frac{mn(n-1)}{N^2}\,\mathbb{E}^{X_1,X_2}\left[\text{cov}^Y\left\{I(X_1,Y),I(X_2,Y)\right\}\right]\\ 
&\qquad\qquad +\frac{m(m-1)n}{N^2}\,\text{var}^X\left[\mathbb{E}^Y\left\{I(X,Y)\right\}\right]\\
\end{align*}
And
\begin{align*}
\mathbb{E}^{X_1,X_2}\left[\text{cov}^Y\left\{I(X_1,Y),I(X_2,Y)\right\}\right]&=\mathbb{E}^{X_1,X_2}\left[\mathbb{E}^{Y}\left\{I(X_1,Y)I(X_2,Y)\right\}\right]\\
&\qquad\qquad-\mathbb{E}^{X_1,X_2}\left[\mathbb{E}^{Y}\left\{I(X_1,Y)\right\}\mathbb{E}^{Y}\left\{I(X_2,Y)\right\}\right]\\
&=\mathbb{E}^Y\left[\mathbb{E}^{X_1,X_2}\left\{I(X_1,Y)I(X_2,Y)\right\}\right]\\
&\qquad\qquad-\mathbb{E}\left\{I(X_1,Y)\right\}\mathbb{E}\left\{I(X_1,Y)\right\}\\
&=\mathbb{E}^Y\left[\mathbb{E}^{X_1}\left\{I(X_1,Y)\right\}\mathbb{E}^{X_2}\left\{I(X_2,Y)\right\}\right]-\mathbb{E}\left\{I(X,Y)\right\}^2\\
&=\mathbb{E}^Y\left[\mathbb{E}^{X_1}\left\{I(X_1,Y)\right\}^2\right]-\mathbb{E}\left\{I(X,Y)\right\}^2\\
&=\text{var}^Y\left\{ \mathbb{E}^{X}\left\{I(X,Y)\right\}\right\}
\end{align*}
Ergo,
$$\text{var} \hat{I}^2_N=\text{var} (\hat{I}^1_N)+\frac{mn(n-1)}{N^2}\,\text{var}^Y\left\{ \mathbb{E}^{X}\left\{I(X,Y)\right\}\right\}+\frac{m(m-1)n}{N^2}\,\text{var}^X\left[\mathbb{E}^Y\left\{I(X,Y)\right\}\right]$$(which makes the formula beautifully symmetric in $X$ and $Y$) and thus
$$\text{var} (\hat{I}^1_N) \le \text{var} (\hat{I}^2_N)$$
with the provision that [for Monte Carlo purposes] $\hat{I}^1_N$ requires $2N$ simulations while $\hat{I}^2_N$ requires $2\sqrt{N}$ simulations. An interesting feature of the variance $\text{var} \hat{I}^2_N$ is that it is approximately
$$\text{var} \hat{I}^2_N\approx\text{var} (\hat{I}^1_N)+\frac{1}{m}\,\text{var}^Y\left\{ \mathbb{E}^{X}\left\{I(X,Y)\right\}\right\}+\frac{1}{n}\,\text{var}^X\left[\mathbb{E}^Y\left\{I(X,Y)\right\}\right]$$and thus $\hat{I}^2_N$ has these two extra variation factors due to the recycling of the $x_i$'s and $y_j$'s. 
A: Note that @Xi'an answer above also provide us with the optimal choice of $n,m$ when using the second estimator. Denoting:
$$a=Var^Y [E^X [I(x,y)]]$$
$$b=Var^X [E^Y [I(x,y)]]$$
by multiplying the expression with $N=mn$ we seek to minimize
$$\min_{nm=N} na+mb$$
which is solved using Lagrange multipliers and yields:
$$n=\sqrt{Nb/a}$$
$$m=\sqrt{Na/b}$$
A: Here is a different way of deriving the same result as @Xi'an, simply using the formula for the variance of a sum:
\begin{align*}
Var(\hat I^2_N)&=\frac{1}{N^2}Var\left(\sum_i\sum_j I(x_i, y_j)\right)\\
&=\frac{1}{N^2}N\cdot Var(I(X, Y) + \frac{1}{N^2} \left(\sum_i\sum_{j\neq j'} Cov\left\{I(x_i, y_j), I(x_i, y_{j'})\right\}\right)
+ \frac{1}{N^2} \left(\sum_{i\neq i'}\sum_{j} Cov\left\{I(x_i, y_j), I(x_{i'}, y_{j})\right\}\right) \qquad\text{[the other covariances are 0 by independence]}\\
&=\frac 1N Var(I(X, Y)) + \frac{nm(m-1)}{N^2} Cov\left\{I(X, Y), I(X, Y')\right\} + \frac{n(n-1)m}{N^2}Cov\left\{I(X, Y), I(X',Y)\right\}\\
&=Var(\hat I^1_N)+ \frac{m-1}{N} Cov\left\{I(X, Y), I(X, Y')\right\} + \frac{n-1}{N}Cov\left\{I(X, Y), I(X',Y)\right\}\\
\end{align*}
