# Gaussian random generator

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?

• First, read about the Box-Muller method of generating a pair of independent Gaussian random variables with mean zero and variance $1$. Then, use it $1024$ times to get $2048$ random numbers which will, with high probability, have a sample mean of $0$ and sample variance of $1$. Generating a larger sample and then cherry-picking $2048$ "best samples" is likely to end in disaster. – Dilip Sarwate Apr 24 '13 at 21:27

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!).

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.

• The method described by Jason R is a standard method for generating random samples with a specified CDF. The difficulty arises when the CDF is not invertible analytically which the Gaussian CDF (and the $\text{erf}$ function) are not. At best, one can have approximations to the inverse CDF, and how difficult it is to compute the approximation is an issue. The Box-Muller method avoids this difficulty for Gaussian random variables by generating independent Rayleigh and $U[0.2\pi)$ random variables and then taking $R\cos\theta$ and $R\sin\theta$ as independent Gaussian samples. – Dilip Sarwate Apr 25 '13 at 11:36
• @DilipSarwate: Does it really matter that much that ${\rm erf}$ is not "analytically invertible"? I know of no general purpose CPU that has built-in functions for $\sin$, $\cos$, $\exp$, etc. All of them are "approximations" when it gets to implementing a psuedo-random number generator. Sure, you have to be careful, especially with the "tails", but this is a valid approach. – Peter K. Apr 25 '13 at 12:15
• "Does it matter..."? Well, that depends. The Box-Muller approach was developed because some statisticians felt that the inverse CDF approach via approximations was not providing satisfactory results: the distribution of the values generated differed enough from the Gaussian for their purposes. It is indeed the tails where the problems usually occur, and unfortunately, it is accuracy in the tails that is most needed in most dsp and comm applications. For the kinds of uses needed in STAT 101-type applications, the sum of $12$ independent $U[0,1]$ minus $6$ is a good enough $N(0,1)$ sample. – Dilip Sarwate Apr 25 '13 at 13:21
• @DilipSarwate: Understood! The sum-of-uniforms approach is clearly not good for the "tails" of the distribution. Agreed also that getting the inverse CDF approximation right (i.e. accurate in the tails) might be more trouble that it's worth compared with doing Box-Muller for Gaussians. – Peter K. Apr 25 '13 at 16:43

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