Why is pseudo-random sampling applicable for Monte Carlo integration, even though it does not satisfy the CLT requirements?

Question

Assume we have a function $f\left(x\right)$ defined on $\left[0, 1\right]$ that we want to integrate and estimate the error using Monte Carlo method. We generate realizations of uniformly distributed independent random variables $x_n$ on $\left[0, 1\right]$ and find $f_n = f\left(x_n\right)$ where $f_n$ happen to be independent and identically distributed random variables. Hence we may be able to apply the CLT to their sum or average. We can take $\overline f = \sum\limits_{n=1}^{N}f_n$ that may have an approximate normal distribution with expected value $\int\limits_0^1 f\left(x\right)dx$ and sample variance $\sigma^2 = \frac{1}{N-1} \sum\limits_{n=1}^{N} \left(f_n - \overline f\right)^2$. We can use the sample variance as an estimate of the integration error.

Now assume we do the same thing, but instead of random points $x_n$ we generate pseudo-random points $x_n$ (using Mersenne Twister for example). The sequence is determined now, and $f_n$ no longer satisfy the CLT requirements, but it seems like $\overline f$ and $\sigma$ still give good integral value and error estimations respectively.

So my question is how can one formally explain why this is happening. It is clear that pseudo-random points are "random enough" to be close to satisfying the CLT, but it is a heuristic explanation, and I am interested in an actual mathematical proof.

Answer 1

In this 2007 paper `Central limit behavior of deterministic dynamical systems by Tirnakli, Beck, and Tsallis, they state that

What is less known in the physics community but well-known in the mathematics community is the fact that there are also CLTs for the iterates of deterministic dynamical systems. The iterates of a deterministic dynamical system can never be completely independent, since they are generated by a deterministic algorithm. However, if the assumption of IID is replaced by the weaker property that the dynamical system is sufficiently strongly mixing, then various versions of CLTs can be proved for deterministic dynamical systems.

They provide as an example the logistic map $$x_{i+1}=T(x_i)=1-ax_i^2$$ When $a=2$, the sequence is strongly mixing, with invariant density $$\pi(x) = \frac{1}{\pi\sqrt{1−x^2}}$$ (that is $T\pi=\pi$). Then the limiting distribution of the average $$\frac{1}{\sqrt{N}}\sum_{i=1}^N x_i$$ is centred Gaussian, when the original value $x_1$ is a random variable with a absolutely continuous probability distribution. The (limiting) variance of this Gaussian distribution is $1/2$.

Note that the logistic map is chaotic for $a=2$ and hence can be used for a random generator.

Similar central limit theorems can be found in ergodic theory, for some dynamical systems with an infinite ergodic measure (see, e.g., these lecture notes, Theorem 76, or this paper).

Answer 2

There is no mathematical proof for the simple reason that a statement like the CLT (or even weaker conclusions such as "Arithmetic means will converge") is not true for pseudo random numbers. Nevertheless, there is theory available on how good you can approximate integrals by weighted sums of function evaluations. This is applicable to pseudo random numbers in particular but also to quasi random numbers or other designs such as grids.

In addition to integration, pseudo random numbers are used for other applications such as simulation. My answer does address only numerical integration not those other applications.

Counter example

The situation of the CLT is unique that it, among stronger conclusions, ensures the convergence of the estimate for the integral for every (square) integrable function. A similar statement is not possible for finite deterministic sequences (pseudo random or others). Here is a counter example:

• A function $f:[0,1]\rightarrow\mathbb{R}$ with $\int_0^1 f(x)dx=1$ but $\sum_i f(x_i)=0$ for all outputs from your random generator.

This is possible because your generator of pseudo random numbers will have a finite(!) number of possible outputs $x_1, \ldots, x_N$. This number is finite, because it is generated by a deterministic machine with finite memory. See this answer for further explanation. Now, given this finite set, it is a standard exercise in functional analysis to construct a function which is zero exactly on these points but sufficiently large everywhere else to make the integral one (see here). This function can actually be chosen to be infinitely differentiable, i.e. very smooth and "nice".

Error estimation

To rule out counterexamples as the above, error estimates have to take into account quantitative properties beyond smoothness. The best known example is the Koksma-Hlawka inequality:

$$ \left| \frac{1}{N}\sum_i f(x_i) - \int_0^1 f(x) dx \right|\leq V(f) \,D(x_1,\ldots,x_N).$$

$V(f)$ is the total variation of $f$ and $D(x_1,\ldots,x_N)$ the "star-discrepancy" of the set $\{x_1,\ldots,x_N\}$.

The star-discrepancy is a measure for the deviation of the points in the set from being uniform. It is defined as the largest difference between the relative number of points $x_i$ in any sub interval of $[0,1]$ and the Lebesgue measure of the interval, i.e. its length: $$ D(x_1,\ldots,x_N) = \sup_{I\subset [0,1]} \left| \frac{\text{# of $x_i$ in }I}{N} - \lambda(I)\right|.$$

What does it mean in practice?

The Koksma-Hlawka inequality provides a nice split of the integration error into two distinct contributions. One part related to the functions you would like to integrate and one part related to the design, i.e. the choice of points used for integration.

The inequality shows the problem with the counter example. Since one had to force the function down to zero and then up again for each and every $x_i$, the graph of the function will wiggle a lot. This creates a very large Total Variation and drives up the error, even though the choice of design points might be nearly uniform.

On the other hand the inequality explains why pseudo random numbers work (in most cases). Ultimately your pseudo random numbers will fill up the whole interval nearly uniform and your star discrepancy will become very, very small.

Furthermore, it explains why people did not stop at pseudo random numbers but went further to quasi Monte Carlo. Choosing points "at random" will lead to areas of higher and lower concentrations of points, smoothing this out will lower the star discrepancy and increase approximation quality.