17

In order to simulate a normal distribution from a set of uniform variables, there are several techniques:

  1. The Box-Muller algorithm, in which one samples two independent uniform variates on $(0,1)$ and transforms them into two independent standard normal distributions via: $$ Z_0 = \sqrt{-2\text{ln}U_1}\text{cos}(2\pi U_0)\\ Z_1 = \sqrt{-2\text{ln}U_1}\text{sin}(2\pi U_0) $$

  2. the CDF method, where one can equate the normal cdf $(F(Z))$ to a Uniform variate: $$ F(Z) = U $$ and derive $$Z = F^{-1}(U)$$

My question is: which is computationally more efficient? I would think its the latter method - but most of the papers I read use Box-Muller - why?

Additional information:

The inverse of the normal CDF is know and given by:

$$F^{-1}(Z)\; =\; \sqrt2\;\operatorname{erf}^{-1}(2Z - 1), \quad Z\in(0,1).$$

Hence: $$ Z = F^{-1}(U)\; =\; \sqrt2\;\operatorname{erf}^{-1}(2U - 1), \quad U\in(0,1).$$

Xi'an
  • 90,397
  • 9
  • 157
  • 575
user2350366
  • 674
  • 5
  • 12
  • Aren't the two are closely related? Box Muller, I believe, is a particular case for 2-variate generation. – ttnphns Jan 07 '15 at 14:13
  • 2
    What's the inverse of the normal cdf? It can't be computed analytically, only if original CDF is approximated with piecewise linear function. – Artem Sobolev Jan 07 '15 at 14:13
  • Hi Barmaley, I've added some more information above. The Inverse CDF has an expression - however the $\text{erf}^{-1}$ must be calculated computationally - so that might be why box Muller is preferred? I assumed $\text{erf}^1$ would be calculated in lookup tables, much like values of $\text{sin}$ and $\text{cosine}$ are? Hence not that much more computationally expensive? I may be wrong however. – user2350366 Jan 07 '15 at 14:41
  • 3
    There are versions of Box-Muller without sin and cosine. – Xi'an Jan 07 '15 at 14:45
  • _Are_ the values of the sine and cosine functions found by consulting look-up tables? But even if they are, and you want to do something similar and just look up the value of $\operatorname{erf}^{-1}(x)$ in a look-up table, consider that the function is unbounded, and table look-up will have huge gaps. – Dilip Sarwate Jan 07 '15 at 14:50
  • 2
    @Dilip For very low precision applications, such as computer graphics, sine and cosine can indeed be optimized by use of suitable lookup tables. For statistical applications, though, such optimization is never used. Ultimately it's not really any harder to compute $\text{erf}^{-1}$ than $\log$ or $\text{sqrt}$, but on modern computing systems elementary functions related to $\exp$--including the trig functions--tend to be optimized ($\cos$ and $\log$ were basic instructions way back on the Intel 8087 chip!), whereas erf is either unavailable or has been coded at a higher (=slower) level. – whuber Jan 07 '15 at 19:58
  • My first normal random generator was like this: norm=runif;for(i=1:n-1)norm=nrom+runif();norm/sqrt(n) – Aksakal Jan 07 '15 at 21:07

1 Answers1

18

From a purely probabilistic perspective, both approaches are correct and hence equivalent. From an algorithmic perspective, the comparison must consider both the precision and the computing cost.

Box-Muller relies on a uniform generator and costs about the same as this uniform generator. As mentioned in my comment, you can get away without sine or cosine calls, if not without the logarithm:

  • generate $$U_1,U_2\stackrel{\text{iid}}{\sim}\mathcal{U}(-1,1)$$ until $S=U_1^2+U_2^2\le 1$
  • take $Z=\sqrt{-2\log(S)/S}$ and define $$X_1=ZU_1\,,\ X_2=Z U_2$$

The generic inversion algorithm requires the call to the inverse normal cdf, for instance qnorm(runif(N)) in R, which may be more costly than the above and more importantly may fail in the tails in terms of precision, unless the quantile function is well-coded.

To follow on comments made by whuber, the comparison of rnorm(N)and qnorm(runif(N))is at the advantage of the inverse cdf, both in execution time:

> system.time(qnorm(runif(10^8)))
sutilisateur     système      écoulé
 10.137           0.120      10.251 
> system.time(rnorm(10^8))
utilisateur     système      écoulé
 13.417           0.060      13.472` `

using the more accurate R benchmark

        test replications elapsed relative user.self sys.self
3 box-muller          100   0.103    1.839     0.103        0
2    inverse          100   0.056    1.000     0.056        0
1    R rnorm          100   0.064    1.143     0.064        0

and in terms of fit in the tail: enter image description here

Following a comment of Radford Neal on my blog, I want to point out that the default rnorm in R makes use of the inversion method, hence that the above comparison is reflecting on the interface and not on the simulation method itself! To quote the R documentation on RNG:

‘normal.kind’ can be ‘"Kinderman-Ramage"’, ‘"Buggy
 Kinderman-Ramage"’ (not for ‘set.seed’), ‘"Ahrens-Dieter"’,
 ‘"Box-Muller"’, ‘"Inversion"’ (the default), or ‘"user-supplied"’.
 (For inversion, see the reference in ‘qnorm’.)  The
 Kinderman-Ramage generator used in versions prior to 1.7.1 (now
 called ‘"Buggy"’) had several approximation errors and should only
 be used for reproduction of old results.  The ‘"Box-Muller"’
 generator is stateful as pairs of normals are generated and
 returned sequentially.  The state is reset whenever it is selected
 (even if it is the current normal generator) and when ‘kind’ is
 changed.
Xi'an
  • 90,397
  • 9
  • 157
  • 575
  • 3
    (1) How are $\log$ and $\sqrt{}$ implemented and how do they compare to the implementation of $\Phi^{-1}$? It is plausible that the latter could have comparable speed. (A detailed analysis of how $\Phi^{-1}$ has to be implemented shows otherwise; on most platforms there is an order of magnitude difference in timing.) (2) Box-Mueller is *far* slower than generating uniform variates on many systems. (3) Your method generates only nonnegative values of $X_1$ and $X_2$. You probably intended for the $U_i$ to be uniform between $-1$ and $1$ rather than $0$ and $1$. – whuber Jan 07 '15 at 19:51
  • @whuber: thank you for (3), I corrected the mistake. – Xi'an Jan 07 '15 at 20:02
  • 3
    On my system, running `R 3.0.2`, the Box-Mueller method (without sine and cosine and coded using `rowSums` to compute $S$ very quickly) is nevertheless 75% slower than `qnorm(runif(N))`. When running *Mathematica* 9, Box-Mueller is ten times faster than the direct `InverseCDF[NormalDistribution[], #] &`. The point is that *the answer depends on the capabilities of the computational platform* (as well as one's coding skills). This accords with what you say in the first paragraph but might cause readers to re-interpret the rest of your answer. – whuber Jan 07 '15 at 20:16
  • 2
    I agree, `qnorm(runif(N))` is even 20% faster than `rnorm(N)` – Xi'an Jan 07 '15 at 20:19
  • 4
    +1 for the new information. (Actually, your original post was worth upvoting in its own right.) I would like to re-emphasize, though, that the conclusions can change on a different platform. For instance, I would definitely use Box-Mueller if coding in C, Fortran, or Assembly language, where I would have access to extremely fast implementations of algebraic and elementary transcendental functions which would be much faster than *any* implementation of $\Phi^{-1}$, even an approximate one. Since $\sin$ and $\cos$ would be so efficient, I would use them rather than rejection sampling, too. – whuber Jan 07 '15 at 21:16
  • 1
    For comparison, using an i7-3740QM @ 2.7Ghz and R 3.12, for the following calls: `RNGkind(kind = NULL, normal.kind = 'Inversion');At – Avraham Jan 14 '15 at 18:04