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

• no they are not equivalent: in the first case you average independent terms, in the second case dependent ones, hence the variances will be different for different computing times of order $\text{O}(N)$ and $\text{O}(\sqrt{N})$. – Xi'an Feb 20 '18 at 10:55
• @Xi'an, this is a (true) qualitative statement, can you say something about the variance in the second case in terms of n, m? – Uri Cohen Feb 20 '18 at 15:26

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.

• Amazing answer! But why do you state "with the provision that I^1 requires $2N$ simulations while I^2 requires $2sqrt{N}$? I believe both require N evaluations of I(x,y). – Uri Cohen Feb 21 '18 at 6:58
• Evaluations of the function may take much less time than simulations, which is exactly what I stated. – Xi'an Feb 21 '18 at 7:10
• I'm sorry, can you elaborate what you mean by "simulations"? – Uri Cohen Feb 21 '18 at 7:12
• I assumed this was a Monte Carlo question. However, if design is considered (as suggested by your own answer where $m$ and $n$ can be chosen) then the sampling cost may impact the comparison as well. – Xi'an Feb 21 '18 at 8:51

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

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*}