Understanding different Monte Carlo approximation notations

Question

Currently working on a project involving Monte Carlo integrals. I haven't had any prior studies of this method, so hence the following question.

Consider the following expectation:

$$E[f(X)]=\int_A f(x)g(x)dx.$$

Let $X$ be a random variable taking values in $A\subseteq\mathbb{R}^n$. Let $g:A\to\mathbb{R}_+$ be the probability density of $X$, and $f:A\to\mathbb{R}$ a function such that the expectation above is finite.

If $X_1,X_2,...X_N$ be independent random variables with probability density $g$, then by the law of large numbers,

$$\frac{1}{N} \sum_{i=1}^N f(X_i) \to E[f(X)] \quad \text{as N} \to \infty.$$

As far I understand, the sum above is a general Monte Carlo approximation of the integral.

Does the above approximation make any assumption on the pdf, i.e. uniformity and normalisation? If it is a general approximation, it should hold for any pdf, but I have seen different approximations such as $V\frac{1}{N}\sum_{i=1}^N f(X_i)$ and$\frac{1}{N}\sum_{i=1}^N \frac{f(X_i)}{g(X_i)}$, where in the former $V$ denotes the definite integral over the pdf. How are these related and derived?

Answer 1

Yes, the formula you provide should converge to a true answer for arbitrary probability distribution $g(x)$ given infinite sample points. The problem is that you don't want to wait infinitely long. So instead a more interesting question is whether it is likely to converge to a value close to the true value given a finite number of samples. And here the answer depends on the distribution of $f(x)$ in space. For distributions of $f(x)$ that are more or less uniform on the domain of interest basic MC sampling works very well. However, if majority of the stuff in $f(x)$ is concentrated in a small region, especially in higher dimensions, basic MC is completely unfeasible. This problem is actually relatively frequent in real life, where $f(x)$ is a narrow multidimensional gaussian. MC sampling over a cube containing that gaussian is a very bad idea in high dimensions.

In order to solve this problem, people have designed many methods to "sample where it matters". The simplest of those is so called importance sampling. The idea is that you have prior knowledge on how $f(x)$ might be distributed, and sample using some compromise between $g(x)$ and that prior distribution, but then you also have to correct the resulting answer to adjust for the fact that your were not sampling exactly from $g(x)$. That is the last formula you have provided. The middle formula I have not seen before.

Finally, importance sampling depends on the prior. Even in the absence of prior it is possible to do better than basic MC by adaptively finding the prior distribution. However, this is an actively researched open problem.

So, to summarize, there are multiple formulas for MC that all work for arbitrary $f(x)$ and $g(x)$ but have different convergence speeds and are thus better or worse in specific scenarios

Answer 2

In probabilistic terms, the Monte Carlo method (or its justification) is called the Law of Large Numbers. The convergence $$\frac{1}{N} \sum_{i=1}^N f(X_i) \stackrel{\text{a.s.}}{\to} \mathbb E_g[f(X)]\tag{1}$$ does not assume anything but iid-ness of the $X_i$'s and the existence of the expectation.

A more precise characterisation of the convergence requires further properties of the pair $(f,g)$. For instance, the variance of the l.h.s. in (1) goes to zero with $N$ provided the variance$$\text{var}_g(f(X))$$exists (in dimension one). The speed it goes to zero is precisely $\text{O}(\sqrt{N})$ no matter what the dimension of $X$ and no matter what the Monte Carlo method.

The second part of the question alludes to other forms of Monte Carlo approximations. They are a consequence of the non-identifiability of the pair $(f,g)$ in the integral$$\mathfrak I=\int_A f(x)g(x)\text{d}x$$ which can be equally written as$$\mathfrak I=\int_A \frac{f(x)g(x)}{h(x)} h(x)\text{d}x$$for an arbitrary density $h$ with support including $A$ (i.e. positive over $A$). Due to this lack of identifiability, the choice of $h$ is mostly free and the optimal choice of $h$ is $$h^\star(x) = \dfrac{|f(x)|g(x)}{\int_A |f(x)|g(x)\text{d}x}$$ as it achieves the minimal variance. This variance is zero when $f$ is non-negative (or non-positive) over the whole set $A$. Obviously, in practice, this choice of $h$ is unavailable but it explains why simulating from $g$ is rarely the optimal choice.