Double integral, monte carlo estimation Suppose I have pairs of random variables where $X_i$~$U[0,1]$ and $Y_i$~$U[0,1]$ and I want to estimate it
$$\theta=\int_{0.5}^{1}\int_0^{0.5}e^{xy}xydxdy$$
but $\theta$ needs to have variance less than $0.001$, and from what I understand through a book of examples, I would need a variance type $\frac{c}{n}$ where $c<0.001$ I do not know if that's possible, but exercise does not specify how large n can be.
 A: From the way the question is presently phrased, it appears the objective is to estimate
$$\theta = \int_{1/2}^1 \int_0^{1/2} e^{xy} xy dx dy$$
using nothing more than the ability to generate independent uniform variates and some arithmetic.  If we let $X$ have a Uniform$(1/2,1)$ distribution and $Y$ independently have a Uniform$(0,1/2)$ distribution, then this integral is seen to equal
$$\theta = \frac{1}{4}\mathbb{E}(XYe^{XY}).$$
Inspired by your question at Control Variates, Monte Carlo integration, we might look for a simplification of the function $u(X,Y) = XY e^{XY}$ whose expectation is easy to compute exactly.  Eliminating the exponential for the time being, consider
$$\tau = \frac{1}{4}\mathbb{E}(XY) = \frac{1}{4}\mathbb{E}(X)\mathbb{E}(Y).$$
That's dead easy to compute, since (by symmetry) $\mathbb{E}(X) = 3/4$ and $\mathbb{E}(Y)=1/4$, whence
$$\tau = \frac{1}{4}\left(\frac{3}{4}\right)\left(\frac{1}{4}\right) = \frac{3}{64}.$$
Consider, then, $n$ independent realizations of $(X,Y)$ indexed by $i=1,2,\ldots,n$  For each of them we may (easily) compute both $T_i = X_i Y_i$ and $U_i = u(X_i,Y_i)$.  Look at the statistics
$$\bar U = \frac{1}{n}\sum_i U_i,\ \bar T = \frac{1}{n}\sum_i T_i.$$
We have already determined that $\mathbb{E}(\bar T) = 4\tau = 3/16$.  Then, according to the argument in the Wikipedia article on Control Variates, for any real constant $c$
$$\theta = \frac{1}{4}\mathbb{E}(u(X,Y)) = \frac{1}{4}\mathbb{E}(
\bar U) = \frac{1}{4}\mathbb{E}\left(\bar U - c\left(\bar T - \frac{3}{16}\right)\right).$$
The value of $c$ that minimizes the variance of the right hand side is
$$c^{*} = \frac{\text{Cov}(\bar U, \bar T)}{\text{Var}(\bar T)} = \frac{\text{Cov}(U,  T)}{\text{Var}( T)}.$$
The numerator and denominator can both be estimated from the sample to approximate $c^{*}$.  Although this value might not be optimal, for all but extremely small $n$ it should be close to optimal and thereby allow a (substantial) variance reduction.
Finally, how large must $n$ be to achieve a standard error (square root of the estimation variance) that is equal to or better some target value?  We expect the variance of the estimator to be proportional to $1/n$ in the worst case, .  Thus, we could perform estimates for small values of $n$, fit a linear function of these values against $1/n$, and predict the value of $n$ that will result in a standard error equal to the target.  We then use that value of $n$ and see whether it works.  If not, we improve our guess by including this new datum in the fit and keep trying until we get it.
The example in the question requires only two samples to achieve a standard error of $\sqrt{0.001}$.  Here, for instance, is a histogram of $10,000$ independent estimates using just two $(X,Y)$ samples each.  Its variance is $0.00002$, far better than the target of $0.001$.  (It actually is biased a little low because $n=2$ is so small: the true value of the integral is shown by the vertical red line.)

A: $\int_{0.5}^1\int_0^{0.5} xy\,dx\,dy\approx \frac{1}{n}\sum_{i=1}^n x_i y_iI(x_i<0.5)\cdot I(y_i>0.5)$
where $I$ is the indicator function. That is,
$\phi(x, y) = x yI(x<0.5)\cdot I(y>0.5)$. In the denominator there is an $n$ since you have $n$ observations of the variable $\phi(X, Y)$.
If the variables are pairwise independent then
take $Z_i = X_i/2$ and $W_i = Y_i/2 + 0.5$. Then
$\mathbb{E}[Z_iW_i] = \mathbb{E}[Z_i]\mathbb{E}[W_i]
= 4\int_0^{0.5} z\,dz\int_{0.5}^1 w\,dw = 4\int_{0.5}^1\int_0^{0.5} zw\,dz\,dw$ so
$\int_{0.5}^1\int_0^{0.5} xy\,dx\,dy = \frac{1}{4}\mathbb{E}[Z_i]\mathbb{E}[W_i]$.
Edit:
Since my post the question has changed!
