Let
$$f(x,y) = \exp \left(- \frac{1}{2}a x^2 - \frac{1}{2}c y^2 + bxy \right)$$
where $a,b\ge 0$. I want to integrate numerically:
$$\frac{\int_{x_0}^{x_1}\mathrm{d}x \int_{y_0}^{y_1}\mathrm{d}y \, f(x,y) x^n}{\int_{x_0}^{x_1}\mathrm{d}x \int_{y_0}^{y_1}\mathrm{d}y \, f(x,y)}$$
where $-\infty < x_0 < x_1 < \infty$, $-\infty < y_0 < y_1 < \infty$, and $n\in\{1,2\}$.
A naive method has problems when the peak of $f(x,y)$ is far from the rectangle of integration.
Is there a method that I can use?