Gaussian random generator

Question

I have quite a straight-forward question. What I aim for is the generation of a certain set of random numbers with a normal distribution (mu = 0, sigma = 1).

Now, the best way to approach the above mentioned gaussian bell is to generate quite a large number of samples.

Unfortunately, in my case, I can only generate a retained subset of samples: usually around 2048 samples, which lead to a poor fit.

I am actually wondering if there is any way to improve the fitness, by first generating a larger set of samples or more subsets of the same size (as 10 times 2048 samples) and then pick those 2048 which approximate the best for desired mean and standard deviation (and the gaussian of course).

What would you suggest or how would you proceed?

Answer 1

To elaborate on Dilip's answer (which is perfectly correct, though in practice the Ziggurat method is much more computationally efficient than Box-Mueller):

One key ingredient missing from your reasoning is whether you want your samples to be independent. It is not clear from your question, but this is the most common situation...

If you want your samples to be independent, then you will have to accept that their empirical mean and variance will not exactly be 0 and 1 - what is known is that the larger the sample, the closer the empirical values will be to 0, 1 - this is absolutely fine - this is a "feature" of the normal distribution rather than a "bug". It is a bad idea to try to select the samples that look more "random" or "well-behaved"; because this very selection process is making them less random! For example, if you come up with a process that always results in 2048-sample blocks with exactly $\mu = 0$ and $\sigma = 1$, it would mean that there are only 2046 possible degrees of freedom in your generation process, and thus that your samples are not independent (given the 2046 first samples, you can guess the values of the 2047th and 2048th using the $\mu = 0$ and $\sigma = 1$ equation... so these last two values are not random!).

Answer 2

Dilip and pichenettes have already pointed out two methods for generating Gaussian random variables (the Box-Muller transform and the Ziggurat algorithm). For completeness, I will point out another: inverse transform sampling. I recently had the need to create a maximum-throughput software Gaussian random number generator, and after evaluating all of the methods that I could find, I settled on this one, as it proved to be the fastest on the target system.

The idea of inverse transform sampling is that you start with a base random generator that yields uniform values over some interval, typically $[0, 1)$ (or some variation of whether the interval is open or closed on either end). You then apply those uniform values to a function that is chosen such that the resulting outputs will have the desired distribution. The setup for this technique goes like this:

• Select the distribution that you want to generate random numbers from. Determine its cumulative distribution function (cdf) $F_x(x)$. For the Gaussian distribution, this function is: $$ F_x(x) = \frac{1}{2}\left(1 + \operatorname{erf}\left(\frac{x-\mu}{\sqrt{2\sigma^2}}\right)\right) $$ where $\operatorname{erf}(z)$ is the error function. For the standard Gaussian distribution, $\mu = 0$ and $\sigma = 1$.

• Invert the distribution's CDF to yield its inverse CDF $F^{-1}_x(x)$, sometimes called the quantile function. This is often difficult or impossible to do in a compact analytical form, so you might need to rely upon numerical approximations. For the standard Gaussian distribution, this is the probit function, which can be obtained from the above expression for the cdf: $$ F^{-1}_x(x) = \sqrt{2}\operatorname{erf}^{-1}(2x - 1) $$

  If your computing environment has an implementation of the inverse error function, then this should be easy to evaluate. Otherwise, you'll need to rely upon some numerical approximation of the function.

Armed with the above information, you can now generate values from the desired distribution. The simple procedure is as follows:

1. Generate a uniform value $u$ on the interval $[0, 1)$.

2. Calculate $y = F^{-1}_x(u)$.

3. Output $y$.

Assuming that you have a good base uniform generator (and a good implementation of the inverse CDF function), the resulting samples of $y$ should be distributed very closely to the desired distribution. There are a couple nice things about this technique:

• It is very general; while I talked about Gaussian random number generation above, you can use this successfully with many distributions by just changing the function you use for the inverse CDF.

• If you have a simple implementation of the inverse CDF function, it can be very fast (although the speed of the uniform generator is important). I was able to squeeze more speed out of this method than the Ziggurat algorithm by using an inverse CDF approximation that didn't use any logical checks or branching.

Answer 3

Box Müller is mostly good but not necessarily superior to inverse CDF, as it suffers from the Neave Effect, see H. R. Neave, On using the Box-Muller transformation with multiplicative congruential pseudorandom number generators Applied Statistics, 22, 92-97, 1973, or Tezuka https://www.researchgate.net/publication/3528180_NEAVE_EFFECT_ALSO_OCCURS_WITH_TAUSWORTHE_SEQUENCES

This stack overflow post is related https://stackoverflow.com/questions/2325472/generate-random-numbers-following-a-normal-distribution-in-c-c