What is better in Monte Carlo integration: product of means or mean of products?

Question

Let $X$ and $Y$ be two independent continuous random variables with pdfs $f_X$ and $f_Y$, respectively. Let $\varphi_1$ and $\varphi_2$ be two continuous functions from ${\mathbb R}$ to ${\mathbb R}$. I want to calculate $E[\varphi_1(X)\varphi_2(Y)]< \infty$ numerically. This is,

$$I = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} \varphi_1(x)\varphi_2(y) f_X(x)f_Y(y)dxdy\\ =\int_{-\infty}^{\infty} \varphi_1(x)f_X(x)dx\int_{-\infty}^{\infty}\varphi_2(y)f_Y(y) dy< \infty.$$

Using Monte Carlo integration, I can either approximate $I$ using $$I \approx \frac{1}{n}\sum_{j=1}^{n} \varphi_1(x_j)\varphi_2(y_j),$$ where $(x_j,y_j)$ is an independent sample from the joint distribution of $(X,Y)$. Alternatively, I can use the approximation: $$I \approx \left[\frac{1}{n}\sum_{j=1}^{n} \varphi_1(x_j)\right]\left[\frac{1}{n}\sum_{j=1}^{n} \varphi_2(y_j)\right] .$$ In my case, and minding potential implementation errors, I am obtaining different results.

Question. Is there any reasons to prefer one approximation over the other?

Simulating from the pdfs is straightforward in my case.

Answer 1

Both approaches lead to unbiased estimators, hence comparing their variances amounts to comparing their second moments.

Without loss of generality, since the distributions of $X$ and $Y$ are arbitrary, we can take both $\varphi_1$ and $\varphi_2$ to be the identity transform. Then $$\mathbb E[(\bar X\bar Y)^2]=\mathbb E[\bar X^2]\mathbb E[\bar Y^2]=\frac{n\mathbb E[X]+\sigma^2_X}{n}\frac{n\mathbb E[Y]+\sigma^2_Y}{n}$$ while $$\mathbb E[\overline {XY}^2]=\frac{n\mathbb E[X]\mathbb E[X]+\sigma^2_{XY}}{n}$$ Therefore if both $\mathbb E[X]$ and $\mathbb E[Y]$ are close to zero $$\mathbb E[(\bar X\bar Y)^2]=\text O(n^{-2})\quad\text{and}\quad\mathbb E[\overline {XY}^2]=\text O(n^{-1})$$ while otherwise mileage may vary $$\mathbb E[(\bar X\bar Y)^2]=\text O(n^{-1})\quad\text{and}\quad\mathbb E[\overline {XY}^2]=\text O(n^{-1})$$

Using the same example as in the earlier answer shows that the variability is slightly lower for the product of the averages (left) than for the average of the products (right), obtained over 10³ replications of 10⁴ simulations (the 10³ curves are the cumulated means):

If we compare directly the variances estimated from the 10⁷ simulations in this example $$\text{var}(\bar X\bar Y)=\text{var}(\bar X)\mathbb E[\bar Y^2]+ \text{var}(\bar Y)\mathbb E[\bar X]^2\qquad\qquad\qquad\\=\frac{\text{var}(X)}{n}\frac{\mathbb E[Y^2]+(n-1)\mathbb E[Y]^2}{n}+\frac{\text{var}(Y)\mathbb E[X]^2}{n}$$ is estimated by

[1] 8.990575e-07

while the variance of $\overline{XY}$ is estimated by

[1] 1.097814e-06

If anything, I would favor the $\bar X\times\bar Y$ solution as it also writes as$$\frac{1}{n^2} \sum_{i=1}^n\sum_{j=1}^n X_iY_j$$giving the impression it exploits the independence between both samples in a more systematic manner.

Answer 2

Here is a trivial example in which both distributions are standard normal, the first function is $\phi_1(x) = x^\prime = x+2,$ and the second is $\phi_2(y) = y^\prime = y^2.$ So that $$E(X^\prime Y^\prime) = E(X^\prime)E(Y^\prime) = 2.$$

Letting $n = 10^7,$ simulations in R give essentially the same answer for both methods.

set.seed(1234)
x = rnorm(10^7);  y = rnorm(10^7)
x1 = x + 2;  y1 = y^2

mean(x1*y1)
[1] 1.998185

mean(x1); mean(y1); mean(x1)*mean(y1)
[1] 2.000025
[1] 0.9993708
[1] 1.998766