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

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.

• A variant of the second approximation is better, replacing the first factor by $$\frac1m \sum_{i=1}^m \varphi_1(x_i).$$ Then you can get an answer accurate to within (e.g.) 1% by choosing $m$ and $n$ separately so each factor is accurate to within 0.5%. Commented Nov 11, 2019 at 22:20
• You must have particular distributions and functions in mind because you have done the simulations. Depending on the choices there might be various difficulties with simulation. Please show your work and results so far so someone can answer based on context. Commented Nov 11, 2019 at 22:39
• @BruceET The model I am working on is relatively complex. I have tried to summarize my main question as other aspects of the model are irrelevant. Simulating for the pdfs is not a problem at all, in fact, it is straightforward. Commented Nov 11, 2019 at 22:42
• If you're in doubt, maybe you should let us decide what may be surprisingly problematic and what may seem irrelevant and actually be important. As it stands, I'll stick with my Answer. Commented Nov 11, 2019 at 23:27

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.

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