How to compute integrals using any probability law with Monte Carlo? I am intrested in providing an estimation of :
$\iint C(x,y)dP_X(x)dP_Y(y)$
I am able to generate random numbers from the distribution of $P_Y$ and $P_X$. Therefore I generate a big number (n=10 000) of realizations from those distribution I obtain $x_1,...,x_n $~$P_X $ and $y_1,...,y_n$~$P_Y $ and compute 
$\frac{1}{n^2}\sum_{i,j} C(x_i,y_j) $
I wonder if that method is correct as I could not find precise examples on this topic.
Thank you for any help
 A: I'll derive a formula for the variance of your estimator and describe a couple of  special cases in which this variance is always less than that of the "standard" MC estimator. Of course this variance reduction is at the computational cost of doing $n^2$ (rather than only $n$) function evaluations.  
Your integral is evidently with respect to a product measure, so I assume $X$ and $Y$ are independent, and suppose $X,X',X_1,\ldots,X_n\sim P_X$ and $Y,Y',Y_1,\ldots,Y_n\sim P_Y$, and that all of these are mutually independent real random variables. For the expectation
$$\mu_f=\mathbb{E}[f(X,Y)]=\int\int f(x,y)\,dP_X(x)\,dP_Y(y)$$
the "standard" estimator is
$$\hat{\mu_1} = \frac{1}{n}\sum_{i}f(X_i,Y_i)
$$
whereas yours is
$$\hat{\mu_2} = \frac{1}{n^2}\sum_{(i,j)}f(X_i,Y_j).
$$
(All summation variables range over $1..n$.)
Using the properties of expectation and variance operators, it's straightforward to show the following:

$$\begin{align}
\mathbb{E}[\hat{\mu_1}]&=\mathbb{E}[f(X,Y)]\\
\mathbb{V}[\hat{\mu_1}]&=\frac{1}{n}\mathbb{V}[f(X,Y)]\\
\\
\mathbb{E}[\hat{\mu_2}]&=\mathbb{E}[f(X,Y)]\\
\mathbb{V}[\hat{\mu_2}]&= \frac{1}{n^2}\left(\mathbb{V}[f(X,Y)] + (n-1)(C_f + C_f')\right)
\end{align}
$$
  where 
  $$\begin{align}
C_f&=\text{cov}\left(f(X,Y),f(X,Y')\right)\\
C_f'&=\text{cov}\left(f(X,Y),f(X',Y)\right).
\end{align}
$$

Thus both $\hat{\mu_1}$ and $\hat{\mu_2}$ are unbiased consistent estimators of $\mu_f$.
Here's a derivation of the formula for $\mathbb{V}[\hat{\mu_2}]$ :
$$\begin{align}
\mathbb{V}[\hat{\mu_2}]&=\mathbb{V}\left[\frac{1}{n^2}\sum_{(i,j)}f(X_i,Y_j)\right]\\
&=\frac{1}{n^4}\left(\sum_{(i,j)}\mathbb{V}[f(X_i,Y_j)] + \sum_{\left((i,j),(i^{'},j^{'})\right):\atop (i,j)\ne(i^{'},j^{'})}\text{cov}\left( f(X_i,Y_j),\,f(X_{i^{'}},Y_{j^{'}})  \right) \right)\\
&=\frac{1}{n^4}\left(\sum_{(i,j)}\mathbb{V}[f(X_i,Y_j)] 
 + \sum_{\left((i,j),(i^{'},j^{'})\right):\atop i\ne i^{'}\land j\ne j^{'}}\text{cov}\left( f(X_i,Y_j),\,f(X_{i^{'}},Y_{j^{'}})  \right)
 + \sum_{\left((i,j),(i^{'},j^{'})\right):\atop i=i^{'}\land j\ne j^{'}}\text{cov}\left( f(X_i,Y_j),\,f(X_{i^{'}},Y_{j^{'}})  \right)
 + \sum_{\left((i,j),(i^{'},j^{'})\right):\atop i\ne i^{'}\land j=j^{'}}\text{cov}\left( f(X_i,Y_j),\,f(X_{i^{'}},Y_{j^{'}})  \right)  \right) \\
&=\frac{1}{n^4}\left(n^2\mathbb{V}[f(X,Y)] + 0 + n^2 (n-1) C_f + n^2 (n-1) C_f' \right)\\
&=\frac{1}{n^2}\left(\mathbb{V}[f(X,Y)] + (n-1)(C_f + C_f')\right)
\end{align}$$

Example 1: $\quad f(X,Y) = g(X)\,h(Y),\quad \mathbb{E}[g(X)]=0$.
In this case, we find
$\mathbb{V}[f(X,Y)] = \mathbb{E}[g(X)^2\,h(Y)^2]=\mathbb{V}[g(X)]\,\mathbb{E}[h(Y)^2]$
$C_f = \mathbb{E}[g(X)\,h(Y)\,g(X)\,h(Y')]=\mathbb{V}[g(X)]\,(\mathbb{E}[h(Y)])^2$
$C_f' = \mathbb{E}[g(X)\,h(Y)\,g(X')\,h(Y)]=0$
and the above formula yields 
$$\begin{align}\mathbb{V}[\hat{\mu_2}] &= \frac{1}{n^2}\left( \mathbb{V}[g(X)]\,\mathbb{E}[h(Y)^2] + (n-1)\mathbb{V}[g(X)](\mathbb{E}[h(Y)])^2 \right) \\
&= \frac{1}{n}\mathbb{V}[g(X)]\left( (\mathbb{E}[h(Y)])^2 + \frac{1}{n}\mathbb{V}[h(Y)] \right)\end{align}$$
giving 
$$\begin{align}\mathbb{V}[\hat{\mu_1}] - \mathbb{V}[\hat{\mu_2}] &= \frac{1}{n}\,(\mathbb{V}[g(X)]\,\mathbb{E}[h(Y)^2] - \frac{1}{n}\mathbb{V}[g(X)]\left( (\mathbb{E}[h(Y)])^2 + \frac{1}{n}\mathbb{V}[g(Y)]\right)\\
&=\left(\frac{1}{n}\right)\left(1-\frac{1}{n}\right)\,\mathbb{V}[g(X)]\,\mathbb{V}[h(Y)]\\
&>0\text{ for all }n>1.
\end{align}$$

Example 2: $\quad f(X,Y) = g(X)+h(Y),\quad \mathbb{E}[g(X)]=0$.
In this case, we find
$\mathbb{V}[f(X,Y)] = \mathbb{V}[g(X)]+\mathbb{V}[h(Y)]$
$C_f = \text{cov}(g(X)+h(Y),\,g(X)+h(Y')) =\mathbb{V}[g(X)]-(\mathbb{E}[h(Y)])^2$
$C_f' = \text{cov}(g(X)+h(Y),\,g(X')+h(Y))=\mathbb{V}[h(Y)]$
and the above formula yields 
$$V[\hat{\mu_2}]=\frac{1}{n}\left(\mathbb{V}[g(X)] + \mathbb{V}[h(Y)]\right)-\left(\frac{1}{n}\right)\left(1-\frac{1}{n}\right)(\mathbb{E}[h(Y)])^2$$
giving
$$\begin{align}V[\hat{\mu_1}]-V[\hat{\mu_2}] &= \left(\frac{1}{n}\right)\left(1-\frac{1}{n}\right)(\mathbb{E}[h(Y)])^2 \\
&> 0 \text{ for all }n>1.\end{align}$$

NB: MC simulations with various arbitrarily chosen functions $f$ have not yet revealed any exception to $\mathbb{V}[\hat{\mu_2}]\le\mathbb{V}[\hat{\mu_1}]$, but I doubt that this holds in general. (?)  
