Computing Normalizing Constants with Monte Carlo Integration

Question

I am trying to verify a Monte Carlo integrator but I'm having some trouble.

Given the distribution $p(x)$, I can compute $\int p(x) dx=1$ using existing method to confirm that $p(x)$ is implemented correctly.

I then sample $x_i \sim p$ and compute $\frac{1}{N}\sum_{i=1}^N p(x_i)$ and expect that as $N \rightarrow \infty$, that this quantity is $1.$

However, I'm finding that my code is giving $\approx 0.28$. Does anyone see what I'm doing wrong?

To clarify a bit---my ultimate goal is to compute $\int f(x) \pi(x) dx$ with importance sampling (i.g., $\int f \frac{\pi}{p} p dx$). As a sanity check, the test described above is a test case with $f=1$ and $p=\pi$.

Answer 1

You should not expect that the average of a collection of of density values at values sampled from the same density would be close to 1.

You're in effect computing a Monte Carlo approximation for the integrated squared density. This would not usually be 1. It might be greater than 1 or less than 1.

For example, if you do this for a standard normal, you get $\frac{1}{2\sqrt{\pi}}\approx 0.28$ (perhaps that's what you tried):

\begin{eqnarray} \int_{-\infty}^{\infty}\left(\frac{1}{\sqrt{2\pi}}e^{-\frac{_1}{^2}x^{2}}\right)^{2}dx & =&\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{-x^{2}}dx\\ & =&\frac{\sqrt{\pi}}{2\pi}\\ & =&\frac{1}{2\sqrt{\pi}} \end{eqnarray}

Simulations do indeed give values close to the theoretical calculation. In R:

 mean(dnorm(rnorm(10000)))
 [1] 0.281636

If you do it for a standard exponential you get $\frac12$. If you do it for a standard uniform, you get $1$.

Answer 2

Just to clarify ... If the aim is to approximate: $$ \hat{\phi} = \int_E \phi(x)f(x)dx $$

for some density $f$ on space $E$ and some $\phi : E → \mathbb{R}$. Since you are not taking the integral of $f(x)$ as clearly shown above then you shouldn't expect to get a value of $1$.

Then you take appropriate $i.i.d$ samples $X_i$ from $f$ and approximate $\phi$ with

$$ \hat{\phi} = \frac1{n}\sum_{i=1}^n \phi(X_i)$$

However, typically the variance is $\mathcal{O}(1/n)$ so there will be an error of $\mathcal{O}_\mathbb{P}(1/ \sqrt{n})$.