How to calculate individual moments of a 2-dimensional distribution via Monte Carlo Integration

Question

I've recently been using Monte Carlo integration to calculate a particular integrals which I can do fine but I've hit a problem where I can't calculate the individual moments of my distribution for each dimension, and I'm not entirely sure what to do.

For example, let's say I have a unnormalised 2-dimensional Gaussian centred at (3,3) with a standard deviation of (1,1), i.e. the covariance is the identity matrix. This can be described by,

$$ f(x,y) = e^{-0.5((x-3)^2+(y-3)^2)}$$

and to calculate the integral, $I$, via Monte Carlo can be done by,

$$I = \int_{-\infty}^{\infty} f(x,y) \ dx \ dy = \int_{-\infty}^{\infty} \frac{f(x,y)}{p(x,y)}p(x,y) \ dx \ dy = \langle f(x,y) \rangle _{x,y \sim p(x,y)} \approx \frac{1}{N} \sum_{i}^{N} \frac{f(x_{i},y_{i})}{p(x_{i},y_{i})}$$

and by the law of large numbers this will equal the integral, $I$, as $N \rightarrow \infty$.

Now, in order to calculate the moments of my distribution I essentially repeat the same process but with a factor of $x^{a}$ or $y^{a}$ where $a$ is the order of the moment of the distribution. In the case of $a = 1$ for x would be the centre of x, and would be calculated by,

$$I = \int_{-\infty}^{\infty} x f(x,y) \ dx \ dy = \int_{-\infty}^{\infty} \frac{x f(x,y)}{p(x,y)}p(x,y) \ dx \ dy = \langle x f(x,y) \rangle _{x,y \sim p(x,y)} \approx \frac{1}{N} \sum_{i}^{N} \frac{x_{i} f(x_{i},y_{i})}{p(x_{i},y_{i})}$$.

But the thing is, when I calculate this I get the answer to be 9, instead of 3 which is the centre for x. Interestingly enough, this method seems to give the product of the centres. As in, I changed the centre to be (3,2) and the answer is 6 in this case. For reference, the samples of x and y are drawn from a model probability distribution which approximates f(x,y). This is because the probability distribution with minimal variance is proportional to the function you're integrating.

Has anyone calculated moments of distributions via Monte Carlo Integration? If so, how can one integrate just purely along $x$ or $y$ when the function returns just $f(x,y)$. Do I need to times the expression by p(y) and assume p(x,y) can be approximated by $p(x,y) \propto p_x(x) p_y(y)$. If that's the case I'm not sure how to do that as the function for p(x,y) returns just a scalar value for the joint probability at that point in the xy-plane.

Many thanks in advance for the help!

Answer 1

Monte Carlo is a very broad set of techniques, what you are doing is importance sampling. And unless you make adjustments, importance sampling requires access to probability densities and not just regular unnormalized functions. Are you making this adjustment? If not, your true density is

$$\hat f(x,y) = \frac{f(x,y)}{Z}, $$

where $Z$ is the normalizing constant. If you don't have access to the normalized density, you should adjust for it via

$$\sum_{i}^N \frac{x_i \frac{f(x_i, y_i)}{p(x_i, y_i)}}{\sum_j^N \frac{f(x_j, y_j)}{p(x_j, y_j)}}.$$

In this way, the numerator and denominator both contain the ratio of normalizing constants, and they subsequently cancel.

But if you are doing things correctly, it should work. You are estimating the mean of the function $g(x,y): \mathbb R^2 \rightarrow \mathbb R$ given by $g(x,y) = x$ under your distribution. By Fubini's theorem, you can therefore just switch the order of the integrals and

$$\int g(x,y) \hat f(x,y)dx dy = \int x \left (\int \hat f(x,y) dy \right )dx = \int x \hat f(x), $$

where I by $\hat f$ denote the normalized density (not that it matters for Fubini's theorem to be applicable here.) Thus it is legal to estimate the mean of $g(x,y)$ just by looking at the marginal. In practice, we can just discard samples from the dimensions we do not care about and just directly use the remaining.