I will try to phrase the question in a general way, then give my specific case as an example.
Suppose I want to evaluate $Q = \mathbb E \left[ f\left(X, Y \right) \right]$ where $X$ and $Y$ are independent random variables. Usually we draw $N$ samples $\left\{\left(x_i,y_i\right)\right\}_{i=1}^{i=N}\,\,\,x_i,y_i\sim X,Y$, and estimate: $$ \hat {Q}_N =\frac {1}{N}\sum _{{i=1}}^{N}f(x_i,y_i) $$ In this case the error bars would be the standard error of the mean (see for instance: Monte Carlo integration).
Suppose that instead we look at the sample $\left\{\left(x_i^{(j)},y_i\right)\right\}_{i=1,j=1}^{i=N,j=M}$, that is, for each sample from $Y$ we have $M$ samples from $X$ (so the samples are now correlated). How should I calculate the error bars if I estimate: $$ \hat {Q}_{N,M} =\frac {1}{N}\sum _{{i=1}}^{N}\frac {1}{M}\sum _{{j=1}}^{M}f(x_i^{(j)},y_i) $$
Is that estimation even valid?
The reason I need this is because calculating $f$ for arbitrary inputs is expansive, but for a given $y$ I can easily calculate $f(x,y)$ for many $x$ values.
