15
$\begingroup$

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

$\endgroup$
  • 1
    $\begingroup$ What's the inverse of the normal cdf? It can't be computed analytically, only if original CDF is approximated with piecewise linear function. $\endgroup$ – Artem Sobolev Jan 7 '15 at 14:13
  • $\begingroup$ Aren't the two are closely related? Box Muller, I believe, is a particular case for 2-variate generation. $\endgroup$ – ttnphns Jan 7 '15 at 14:13
  • $\begingroup$ 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. $\endgroup$ – user2350366 Jan 7 '15 at 14:41
  • 2
    $\begingroup$ There are versions of Box-Muller without sin and cosine. $\endgroup$ – Xi'an Jan 7 '15 at 14:45
  • 2
    $\begingroup$ @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. $\endgroup$ – whuber Jan 7 '15 at 19:58
15
$\begingroup$

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` `

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.
$\endgroup$
  • 3
    $\begingroup$ (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$. $\endgroup$ – whuber Jan 7 '15 at 19:51
  • 2
    $\begingroup$ 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. $\endgroup$ – whuber Jan 7 '15 at 20:16
  • 1
    $\begingroup$ I agree, qnorm(runif(N)) is even 20% faster than rnorm(N) $\endgroup$ – Xi'an Jan 7 '15 at 20:19
  • 3
    $\begingroup$ +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. $\endgroup$ – whuber Jan 7 '15 at 21:16
  • 1
    $\begingroup$ For comparison, using an i7-3740QM @ 2.7Ghz and R 3.12, for the following calls: RNGkind(kind = NULL, normal.kind = 'Inversion');At <- microbenchmark(A <- rnorm(1e5, 0, 1), times = 100L);RNGkind(kind = NULL, normal.kind = 'Box-Muller');Bt <- microbenchmark(B <- rnorm(1e5, 0, 1), times = 100L) I get mean 11.38363 median 11.18718 for inversion and mean 13.00401 median 12.48802 for Box-Muller $\endgroup$ – Avraham Jan 14 '15 at 18:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.