Best method for transforming low discrepancy sequence into normal distribution?

Question

I've been using low discrepancy sequences for a while for Uniform Distributions, as I've found their properties useful (mainly in computer graphics for their random appearance and their ability to densely cover [0,1] in an incremental fashion).

For example, random values above, Halton sequence values below:

I was considering using them for some financial analysis planning, but I need different distributions than just uniform. I started off trying to generate a normal distribution from my uniform distributions via the Marsaglia polar algorithm, but the results don't seem as good as with the uniform distribution.

Another example, again random above, Halton below:

My question is: What is the best method for getting a normal distribution with the properties I get from a uniform low discrepancy sequence - coverage, incremental fill-in, non-correlation across multiple dimensions? Am I on the right track, or should I be taking a completely different approach?

(Python code for uniform and normal distributions I use above: Gist 2566569)

Answer 1

You can transform from $\mathcal U(0, 1)$ random variables to any other distribution using the inverse of the CDF, also called the percent point function or quantile function. It's implemented in scipy as scipy.stats.norm.ppf.

Answer 2

I recently stumbled upon this problem. Naively I thought that any transformation from uniform would work, so I plugged in a 1D Sobol (and Halton) sequence as if the sequence where a random number generator into an std::normal_distribution<> variate. To my surprise it didn't work, it obviously generated a non normal distribution.

Ok, then I took the Numerical Recipes Third Edition Chapter 7.3.9 Normal_dev function to generate normal numbers from the Sobol or Halton sequences by the method of "Ratio-of-Uniforms" and it failed in the same way. Then I though, ok, if you look at the code, it takes two uniform random numbers to generate two normally distributed random numbers. Perhaps if I used a Sobol (or Halton) 2D sequence it will work. Well, it failed again.

The I remembered about the "Box-Muller method" (mentioned in the comments) and since it has a more geometric interpretation then I though it could work. Well, it did work! I was very excited an starting doing other test, the distribution looks normal.

The problem I saw was that the distribution was no better than random, it terms of filling, so I was a bit disappointed, but ready to publish the result.

Then I did a deeper search (now that I knew what to look for), and it turn out that there is already a paper on this subject: http://www.sciencedirect.com/science/article/pii/S0895717710005935

In this paper it is actually claimed

Two well known methods used with pseudorandom numbers are the Box–Muller and the inverse transformation methods. Some researchers and financial engineers have claimed that it is incorrect to use the Box–Muller method with low-discrepancy sequences, and instead, the inverse transformation method should be used. In this paper we prove that the Box–Muller method can be used with low-discrepancy sequences, and discuss when its use could actually be advantageous.

So the overall conclusion is this:

1) You can use the Box-Muller on 2D low discrepancy sequences to obtain normally distributed sequences. But my few experiments seem to show that the low discrepancy/space, e.g. filling properties are lost in the normal-transformed sequence.

2) You can use the inverse method, presumably the low discrepancy/space filling properties will be preserved.

3) Ratio-of-Uniforms cannot be used.

EDIT: This https://mathoverflow.net/a/144234 points to the same conclusions.

I made an illustration (the first figure (Ratio-of-uniforms on Sobol) shows that the distribution obtained is not normal but the ohters (Box-Muller and random for comparison) are ):

EDIT2:

The main point is that, even if you find a method that can transform the "distribution" of a low discrepancy sequence, it is not obvious that you will preserve the good filling properties. So you are not better than with a truly random (standard) normal distribution. I have yet to find a method that is low discrepancy and yet it fills nicely with a non uniform distribution. I bet such method is very non-obvious and perhaps an open problem.

Answer 3

There are two good methods. First, as noted above, an accurate approximation to the inverse of the Gaussian distribution can be used. Then one can transform any low-discrepancy sequence into a Gaussian.

The second method is the Box-Muller. This method requires two in put numbers (R and A) and generates two outputs. A two-dimensional low-discrepancy sequence is needed. One takes (for example in the Halton Sequence), pairs of primes are used, one for the radial component (R) and one for the angular component (A). One gets Sqrt(-2*Log(R)) for the radial component and Sin(2*Pi*A) and Cos(2*Pi*A) for the angular components. Multiplying the radial by the two angular components (separately) gives two Gaussians. The efficiency is the same as above; two quasi-random inputs and two Gaussian outputs.

Any multi-dimensional low-discrepancy sequence can be used, depending on the dimensionality of the problem.

Answer 4

The most native method would be indeed using the inverse CDF to transform into normal Gaussian, but there are also a problems with this. If you have e.g. a LDS point set create by rank-1 lattices, then it would be that the starting point is always (0,0), so to transform it you need a little shift, best to have the same gap as for the corner (1,1).

So far no problem still, but for an ideal Gaussian distribution N(0,1)+N(0,1) should give the same distribution as the difference. However, this would be not the case using rank-1 lattice LDS and iCDF on each variable, because the starting point in each variable would give a certain iCDF, like $-3\sigma$ (depending on N), so the difference would be $-6\sigma$.

And that is a too extreme value, leading really to a systematic error (e.g. you will not get $+6\sigma$ at the other side). Best inspect your transformed LDS also for sum and differences, check for such extreme points and also for skew and kurtosis.