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

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!

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.

• Thanks for quick response, @Forgottenscience! I'm sure my probability density is normalised because when I calculate the expectation of function, i.e. the integral, I get the correct results within the error of the standard deviation. So, with the last comment if I'm to calculate the first moment for X I need to split f(x,y) and p(x,y) into two forms of f(x,y) $\propto$ f_x(x)f_y(y) and likewise for p(x,y) as I stated above? Then just do $\frac{1}{N} \sum_i (x_i f_x(x_i)/p_x(x_i))$? May 7 '20 at 23:32
• So, $p_x(x) \approx f_x(x) = e^{-0.5((x-3)^2)}$ and do $\frac{1}{N} \sum_{i}^{N} \frac{x_i f_x(x_i)}{p_x(x_i)}$ to calculate the 1st-order moment for x (likewise for y, and for higher order moments) May 8 '20 at 0:09
• The split was more to underscore the legality of your approach. but yes, in your particular case you don't even need to simulate $y$, per the fact that a marginal of a Gaussian is just a Gaussian again with the mean and covariance corresponding to the restricted mean and covariance. May 8 '20 at 19:17
• Right, thanks for the help! It turns out I can also just plug a x factor and integrate like before. So, $\langle x \rangle \approx \frac{1}{N} \sum_{i}^{N} \frac{x_{i}f(x_{i},y_{i})}{p(x_{i},y_{i})}$ and likewise for y and higher order moments for the distribution! May 11 '20 at 15:57