Box–Muller transform

The Box–Muller transform is commonly expressed in two forms. The basic form as given by Box and Muller takes two samples from the uniform distribution on the interval $(0, 1)$ and maps them to two standard, normally distributed samples. The polar form takes two samples from a different interval, $[- 1, 1]$ , and maps them to two normally distributed samples without the use of sine or cosine functions.

The Box–Muller transform was developed as a more computationally efficient alternative to the inverse transform sampling method. The ziggurat algorithm gives a more efficient method for scalar processors (e.g. old CPUs), while the Box–Muller transform is superior for processors with vector units (e.g. GPUs or modern CPUs).

Basic form

Suppose $U_{1}$ and $U_{2}$ are independent samples chosen from the uniform distribution on the unit interval $(0, 1)$ . Set

Z_{0} = \sqrt{- 2 \ln U_{1}} \cos (2 π U_{2}) = R \cos (Θ)

and

Z_{1} = \sqrt{- 2 \ln U_{1}} \sin (2 π U_{2}) = R \sin (Θ) .

Then $Z_{0}$ and $Z_{1}$ are independent random variables with a standard normal distribution.

The derivation is based on a property of a two-dimensional Cartesian system, where $X$ and $Y$ coordinates are described by two independent and normally distributed random variables, the random variables for R² and Θ (shown above) in the corresponding polar coordinates are also independent and can be expressed as $R^{2} = - 2 \cdot \ln U_{1}$ and $Θ = 2 π U_{2} .$

Because R² is the square of the norm of the standard bivariate normal variable (X, Y), it has the chi-squared distribution with two degrees of freedom. In the special case of two degrees of freedom, the chi-squared distribution coincides with the exponential distribution, and the equation for R² above is a simple way of generating the required exponential variate.

4 more chapters are free with a Lably account.

Polar form

The polar form was first proposed by James Bell and then modified by R. Knop. While several different versions of the polar method have been described, the version of R. Knop will be described here because it is the most widely used, in part due to its inclusion in Numerical Recipes. A slightly different form is described as "Algorithm P" by D. Knuth in The Art of Computer Programming.

Given u and v, independent and uniformly distributed in the closed interval [−1, +1], set s = R² = u² + v². If s = 0 or s ≥ 1, discard u and v, and try another pair (u, v). Because u and v are uniformly distributed and because only points within the unit circle have been admitted, the values of s will be uniformly distributed in the open interval (0, 1), too. The latter can be seen by calculating the cumulative distribution function for s in the interval (0, 1). This is the area of a circle with radius $\sqrt{s}$ , divided by $π$ . From this we find the probability density function to have the constant value 1 on the interval (0, 1). Equally so, the angle θ divided by $2 π$ is uniformly distributed in the interval [0, 1) and independent of s.

We now identify the value of s with that of U₁ and $θ / (2 π)$ with that of U₂ in the basic form. As shown in the figure, the values of $\cos θ = \cos 2 π U_{2}$ and $\sin θ = \sin 2 π U_{2}$ in the basic form can be replaced with the ratios $\cos θ = u / R = u / \sqrt{s}$ and $\sin θ = v / R = v / \sqrt{s}$ , respectively. The advantage is that calculating the trigonometric functions directly can be avoided. This is helpful when trigonometric functions are more expensive to compute than the single division that replaces each one.

Just as the basic form produces two standard normal deviates, so does this alternate calculation. $z_{0} = \sqrt{- 2 \ln U_{1}} \cos (2 π U_{2}) = \sqrt{- 2 \ln s} (\frac{u}{\sqrt{s}}) = u \cdot \sqrt{\frac{- 2 \ln s}{s}}$ and $z_{1} = \sqrt{- 2 \ln U_{1}} \sin (2 π U_{2}) = \sqrt{- 2 \ln s} (\frac{v}{\sqrt{s}}) = v \cdot \sqrt{\frac{- 2 \ln s}{s}} .$

Contrasting the two forms

The polar method differs from the basic method in that it is a type of rejection sampling. It discards some generated random numbers, but can be faster than the basic method because it is simpler to compute (provided that the random number generator is relatively fast) and is more numerically robust. Avoiding the use of expensive trigonometric functions improves speed over the basic form. It discards 1 − π/4 ≈ 21.46% of the total input uniformly distributed random number pairs generated, i.e. discards 4/π − 1 ≈ 27.32% uniformly distributed random number pairs per Gaussian random number pair generated, requiring 4/π ≈ 1.2732 input random numbers per output random number.

The basic form requires two multiplications, 1/2 logarithm, 1/2 square root, and one trigonometric function for each normal variate. On some processors, the cosine and sine of the same argument can be calculated in parallel using a single instruction. Notably for Intel-based machines, one can use the fsincos assembler instruction or the expi instruction (usually available from C as an intrinsic function), to calculate complex $\exp (i z) = e^{i z} = \cos z + i \sin z,$ and just separate the real and imaginary parts.

Note: To explicitly calculate the complex-polar form use the following substitutions in the general form,

Let $r = \sqrt{- \ln (u_{1})}$ and $z = 2 π u_{2} .$ Then $r e^{i z} = \sqrt{- \ln (u_{1})} e^{i 2 π u_{2}} = \sqrt{- 2 \ln (u_{1})} [\cos (2 π u_{2}) + i \sin (2 π u_{2})] .$

The polar form requires 3/2 multiplications, 1/2 logarithm, 1/2 square root, and 1/2 division for each normal variate. The effect is to replace one multiplication and one trigonometric function with a single division and a conditional loop.

Tails truncation

When a computer is used to produce a uniform random variable it will inevitably have some inaccuracies because there is a lower bound on how close numbers can be to 0. If the generator uses 32 bits per output value, the smallest non-zero number that can be generated is $2^{- 32}$ . When $U_{1}$ and $U_{2}$ are equal to this the Box–Muller transform produces a normal random deviate equal to $δ = \sqrt{- 2 \ln (2^{- 32})} \cos (2 π 2^{- 32}) \approx 6.660$ . This means that the algorithm will not produce random variables more than 6.660 standard deviations from the mean. This corresponds to a proportion of $2 (1 - Φ (δ)) ≃ 2.738 \times 10^{- 11}$ lost due to the truncation, where $Φ (δ)$ is the standard cumulative normal distribution. With 64 bits the limit is pushed to $δ = 9.419$ standard deviations, for which $2 (1 - Φ (δ)) < 5 \times 10^{- 21}$ .

The smallest positive number representable in IEEE-standard single- and double-precision floating point numbers are subnormal and are $2^{- 149}$ and $2^{- 1074}$ respectively. The smallest normal numbers are $(2^{- 126} - 2^{- 149})$ and $(2^{- 1022} - 2^{- 1074})$ . Most random number generators that output a uniformly distributed floating point number do not, however, use this precision because they work by simply converting a random L-bit integer to floating point, then dividing it by 2^L or 2^L - 1. It is possible to write an algorithm that uses this extra precision, though to access this precision more random bits will need to be consumed.

Implementation

C++

The standard Box–Muller transform generates values from the standard normal distribution (i.e. standard normal deviates) with mean 0 and standard deviation 1. The implementation below in standard C++ generates values from any normal distribution with mean $μ$ and variance $σ^{2}$ . If $Z$ is a standard normal deviate, then $X = Z σ + μ$ will have a normal distribution with mean $μ$ and standard deviation $σ$ . The random number generator has been seeded to ensure that new, pseudo-random values will be returned from sequential calls to the generateGaussianNoise function.

<span id="mw7g"></span><span class="cp" id="mw7w">#include</span><span class="w" id="mw8A"> </span><span class="cpf" id="mw8Q"><cmath></span> <span class="cp" id="mw8g">#include</span><span class="w" id="mw8w"> </span><span class="cpf" id="mw9A"><limits></span> <span class="cp" id="mw9Q">#include</span><span class="w" id="mw9g"> </span><span class="cpf" id="mw9w"><random></span> <span class="cp" id="mw-A">#include</span><span class="w" id="mw-Q"> </span><span class="cpf" id="mw-g"><utility></span> <span class="c1" id="mw-w">//"mu" is the mean of the distribution, and "sigma" is the standard deviation.</span> <span class="n" id="mw_A">std</span><span class="o" id="mw_Q">::</span><span class="n" id="mw_g">pair</span><span class="o" id="mw_w"><</span><span class="kt" id="mwAQA">double</span><span class="p" id="mwAQE">,</span><span class="w" id="mwAQI"> </span><span class="kt" id="mwAQM">double</span><span class="o" id="mwAQQ">></span><span class="w" id="mwAQU"> </span><span class="n" id="mwAQY">generateGaussianNoise</span><span class="p" id="mwAQc">(</span><span class="kt" id="mwAQg">double</span><span class="w" id="mwAQk"> </span><span class="n" id="mwAQo">mu</span><span class="p" id="mwAQs">,</span><span class="w" id="mwAQw"> </span><span class="kt" id="mwAQ0">double</span><span class="w" id="mwAQ4"> </span><span class="n" id="mwAQ8">sigma</span><span class="p" id="mwARA">)</span> <span class="p" id="mwARE">{</span> <span class="w" id="mwARI"> </span><span class="k" id="mwARM">constexpr</span><span class="w" id="mwARQ"> </span><span class="kt" id="mwARU">double</span><span class="w" id="mwARY"> </span><span class="n" id="mwARc">two_pi</span><span class="w" id="mwARg"> </span><span class="o" id="mwARk">=</span><span class="w" id="mwARo"> </span><span class="mf" id="mwARs">2.0</span><span class="w" id="mwARw"> </span><span class="o" id="mwAR0">*</span><span class="w" id="mwAR4"> </span><span class="n" id="mwAR8">M_PI</span><span class="p" id="mwASA">;</span> <span class="w" id="mwASE"> </span><span class="c1" id="mwASI">//initialize the random uniform number generator (runif) in a range 0 to 1</span> <span class="w" id="mwASM"> </span><span class="k" id="mwASQ">static</span><span class="w" id="mwASU"> </span><span class="n" id="mwASY">std</span><span class="o" id="mwASc">::</span><span class="n" id="mwASg">mt19937</span><span class="w" id="mwASk"> </span><span class="nf" id="mwASo">rng</span><span class="p" id="mwASs">(</span><span class="n" id="mwASw">std</span><span class="o" id="mwAS0">::</span><span class="n" id="mwAS4">random_device</span><span class="p" id="mwAS8">{}());</span><span class="w" id="mwATA"> </span><span class="c1" id="mwATE">// Standard mersenne_twister_engine seeded with rd()</span> <span class="w" id="mwATI"> </span><span class="k" id="mwATM">static</span><span class="w" id="mwATQ"> </span><span class="n" id="mwATU">std</span><span class="o" id="mwATY">::</span><span class="n" id="mwATc">uniform_real_distribution</span><span class="o" id="mwATg"><></span><span class="w" id="mwATk"> </span><span class="n" id="mwATo">runif</span><span class="p" id="mwATs">(</span><span class="mf" id="mwATw">0.0</span><span class="p" id="mwAT0">,</span><span class="w" id="mwAT4"> </span><span class="mf" id="mwAT8">1.0</span><span class="p" id="mwAUA">);</span> <span class="w" id="mwAUE"> </span><span class="c1" id="mwAUI">//create two random numbers, make sure u1 is greater than zero</span> <span class="w" id="mwAUM"> </span><span class="kt" id="mwAUQ">double</span><span class="w" id="mwAUU"> </span><span class="n" id="mwAUY">u1</span><span class="p" id="mwAUc">,</span><span class="w" id="mwAUg"> </span><span class="n" id="mwAUk">u2</span><span class="p" id="mwAUo">;</span> <span class="w" id="mwAUs"> </span><span class="k" id="mwAUw">do</span> <span class="w" id="mwAU0"> </span><span class="p" id="mwAU4">{</span> <span class="w" id="mwAU8"> </span><span class="n" id="mwAVA">u1</span><span class="w" id="mwAVE"> </span><span class="o" id="mwAVI">=</span><span class="w" id="mwAVM"> </span><span class="n" id="mwAVQ">runif</span><span class="p" id="mwAVU">(</span><span class="n" id="mwAVY">rng</span><span class="p" id="mwAVc">);</span> <span class="w" id="mwAVg"> </span><span class="p" id="mwAVk">}</span> <span class="w" id="mwAVo"> </span><span class="k" id="mwAVs">while</span><span class="w" id="mwAVw"> </span><span class="p" id="mwAV0">(</span><span class="n" id="mwAV4">u1</span><span class="w" id="mwAV8"> </span><span class="o" id="mwAWA">==</span><span class="w" id="mwAWE"> </span><span class="mi" id="mwAWI">0</span><span class="p" id="mwAWM">);</span> <span class="w" id="mwAWQ"> </span><span class="n" id="mwAWU">u2</span><span class="w" id="mwAWY"> </span><span class="o" id="mwAWc">=</span><span class="w" id="mwAWg"> </span><span class="n" id="mwAWk">runif</span><span class="p" id="mwAWo">(</span><span class="n" id="mwAWs">rng</span><span class="p" id="mwAWw">);</span> <span class="w" id="mwAW0"> </span><span class="c1" id="mwAW4">//compute z0 and z1</span> <span class="w" id="mwAW8"> </span><span class="k" id="mwAXA">auto</span><span class="w" id="mwAXE"> </span><span class="n" id="mwAXI">mag</span><span class="w" id="mwAXM"> </span><span class="o" id="mwAXQ">=</span><span class="w" id="mwAXU"> </span><span class="n" id="mwAXY">sigma</span><span class="w" id="mwAXc"> </span><span class="o" id="mwAXg">*</span><span class="w" id="mwAXk"> </span><span class="n" id="mwAXo">sqrt</span><span class="p" id="mwAXs">(</span><span class="mf" id="mwAXw">-2.0</span><span class="w" id="mwAX0"> </span><span class="o" id="mwAX4">*</span><span class="w" id="mwAX8"> </span><span class="n" id="mwAYA">log</span><span class="p" id="mwAYE">(</span><span class="n" id="mwAYI">u1</span><span class="p" id="mwAYM">));</span> <span class="w" id="mwAYQ"> </span><span class="k" id="mwAYU">auto</span><span class="w" id="mwAYY"> </span><span class="n" id="mwAYc">z0</span><span class="w" id="mwAYg"> </span><span class="o" id="mwAYk">=</span><span class="w" id="mwAYo"> </span><span class="n" id="mwAYs">mag</span><span class="w" id="mwAYw"> </span><span class="o" id="mwAY0">*</span><span class="w" id="mwAY4"> </span><span class="n" id="mwAY8">cos</span><span class="p" id="mwAZA">(</span><span class="n" id="mwAZE">two_pi</span><span class="w" id="mwAZI"> </span><span class="o" id="mwAZM">*</span><span class="w" id="mwAZQ"> </span><span class="n" id="mwAZU">u2</span><span class="p" id="mwAZY">)</span><span class="w" id="mwAZc"> </span><span class="o" id="mwAZg">+</span><span class="w" id="mwAZk"> </span><span class="n" id="mwAZo">mu</span><span class="p" id="mwAZs">;</span> <span class="w" id="mwAZw"> </span><span class="k" id="mwAZ0">auto</span><span class="w" id="mwAZ4"> </span><span class="n" id="mwAZ8">z1</span><span class="w" id="mwAaA"> </span><span class="o" id="mwAaE">=</span><span class="w" id="mwAaI"> </span><span class="n" id="mwAaM">mag</span><span class="w" id="mwAaQ"> </span><span class="o" id="mwAaU">*</span><span class="w" id="mwAaY"> </span><span class="n" id="mwAac">sin</span><span class="p" id="mwAag">(</span><span class="n" id="mwAak">two_pi</span><span class="w" id="mwAao"> </span><span class="o" id="mwAas">*</span><span class="w" id="mwAaw"> </span><span class="n" id="mwAa0">u2</span><span class="p" id="mwAa4">)</span><span class="w" id="mwAa8"> </span><span class="o" id="mwAbA">+</span><span class="w" id="mwAbE"> </span><span class="n" id="mwAbI">mu</span><span class="p" id="mwAbM">;</span> <span class="w" id="mwAbQ"> </span><span class="k" id="mwAbU">return</span><span class="w" id="mwAbY"> </span><span class="n" id="mwAbc">std</span><span class="o" id="mwAbg">::</span><span class="n" id="mwAbk">make_pair</span><span class="p" id="mwAbo">(</span><span class="n" id="mwAbs">z0</span><span class="p" id="mwAbw">,</span><span class="w" id="mwAb0"> </span><span class="n" id="mwAb4">z1</span><span class="p" id="mwAb8">);</span> <span class="p" id="mwAcA">}</span>

JavaScript

<span id="mwAcQ"></span><span class="cm" id="mwAcU">/* Syntax:</span> <span class="cm" id="mwAcY"> *</span> <span class="cm" id="mwAcc"> * [ x, y ] = rand_normal();</span> <span class="cm" id="mwAcg"> * x = rand_normal()[0];</span> <span class="cm" id="mwAck"> * y = rand_normal()[1];</span> <span class="cm" id="mwAco"> */</span> <span class="kd" id="mwAcs">function</span><span class="w" id="mwAcw"> </span><span class="nx" id="mwAc0">rand_normal</span><span class="p" id="mwAc4">()</span><span class="w" id="mwAc8"> </span><span class="p" id="mwAdA">{</span> <span class="w" id="mwAdE"> </span><span class="kd" id="mwAdI">let</span><span class="w" id="mwAdM"> </span><span class="nx" id="mwAdQ">theta</span><span class="w" id="mwAdU"> </span><span class="o" id="mwAdY">=</span><span class="w" id="mwAdc"> </span><span class="mf" id="mwAdg">2</span><span class="w" id="mwAdk"> </span><span class="o" id="mwAdo">*</span><span class="w" id="mwAds"> </span><span class="nb" id="mwAdw">Math</span><span class="p" id="mwAd0">.</span><span class="nx" id="mwAd4">PI</span><span class="w" id="mwAd8"> </span><span class="o" id="mwAeA">*</span><span class="w" id="mwAeE"> </span><span class="nb" id="mwAeI">Math</span><span class="p" id="mwAeM">.</span><span class="nx" id="mwAeQ">random</span><span class="p" id="mwAeU">();</span> <span class="w" id="mwAeY"> </span><span class="kd" id="mwAec">let</span><span class="w" id="mwAeg"> </span><span class="nx" id="mwAek">R</span><span class="w" id="mwAeo"> </span><span class="o" id="mwAes">=</span><span class="w" id="mwAew"> </span><span class="nb" id="mwAe0">Math</span><span class="p" id="mwAe4">.</span><span class="nx" id="mwAe8">sqrt</span><span class="p" id="mwAfA">(</span><span class="o" id="mwAfE">-</span><span class="mf" id="mwAfI">2</span><span class="w" id="mwAfM"> </span><span class="o" id="mwAfQ">*</span><span class="w" id="mwAfU"> </span><span class="nb" id="mwAfY">Math</span><span class="p" id="mwAfc">.</span><span class="nx" id="mwAfg">log</span><span class="p" id="mwAfk">(</span><span class="nb" id="mwAfo">Math</span><span class="p" id="mwAfs">.</span><span class="nx" id="mwAfw">random</span><span class="p" id="mwAf0">()));</span> <span class="w" id="mwAf4"> </span><span class="kd" id="mwAf8">let</span><span class="w" id="mwAgA"> </span><span class="nx" id="mwAgE">x</span><span class="w" id="mwAgI"> </span><span class="o" id="mwAgM">=</span><span class="w" id="mwAgQ"> </span><span class="nx" id="mwAgU">R</span><span class="w" id="mwAgY"> </span><span class="o" id="mwAgc">*</span><span class="w" id="mwAgg"> </span><span class="nb" id="mwAgk">Math</span><span class="p" id="mwAgo">.</span><span class="nx" id="mwAgs">cos</span><span class="p" id="mwAgw">(</span><span class="nx" id="mwAg0">theta</span><span class="p" id="mwAg4">);</span> <span class="w" id="mwAg8"> </span><span class="kd" id="mwAhA">let</span><span class="w" id="mwAhE"> </span><span class="nx" id="mwAhI">y</span><span class="w" id="mwAhM"> </span><span class="o" id="mwAhQ">=</span><span class="w" id="mwAhU"> </span><span class="nx" id="mwAhY">R</span><span class="w" id="mwAhc"> </span><span class="o" id="mwAhg">*</span><span class="w" id="mwAhk"> </span><span class="nb" id="mwAho">Math</span><span class="p" id="mwAhs">.</span><span class="nx" id="mwAhw">sin</span><span class="p" id="mwAh0">(</span><span class="nx" id="mwAh4">theta</span><span class="p" id="mwAh8">);</span> <span class="w" id="mwAiA"> </span><span class="k" id="mwAiE">return</span><span class="w" id="mwAiI"> </span><span class="p" id="mwAiM">[</span><span class="w" id="mwAiQ"> </span><span class="nx" id="mwAiU">x</span><span class="p" id="mwAiY">,</span><span class="w" id="mwAic"> </span><span class="nx" id="mwAig">y</span><span class="w" id="mwAik"> </span><span class="p" id="mwAio">];</span> <span class="p" id="mwAis">}</span>

Julia

<span id="mwAi8"></span><span class="s" id="mwAjA">"""</span> <span class="s" id="mwAjE"> boxmullersample(N)</span> <span class="s" id="mwAjI">Generate `2N` samples from the standard normal distribution using the Box-Muller method.</span> <span class="s" id="mwAjM">"""</span> <span class="k" id="mwAjQ">function</span><span class="w" id="mwAjU"> </span><span class="n" id="mwAjY">boxmullersample</span><span class="p" id="mwAjc">(</span><span class="n" id="mwAjg">N</span><span class="p" id="mwAjk">)</span> <span class="w" id="mwAjo"> </span><span class="n" id="mwAjs">z</span><span class="w" id="mwAjw"> </span><span class="o" id="mwAj0">=</span><span class="w" id="mwAj4"> </span><span class="kt" id="mwAj8">Array</span><span class="p" id="mwAkA">{</span><span class="kt" id="mwAkE">Float64</span><span class="p" id="mwAkI">}(</span><span class="nb" id="mwAkM">undef</span><span class="p" id="mwAkQ">,</span><span class="n" id="mwAkU">N</span><span class="p" id="mwAkY">,</span><span class="mi" id="mwAkc">2</span><span class="p" id="mwAkg">);</span> <span class="w" id="mwAkk"> </span><span class="k" id="mwAko">for</span><span class="w" id="mwAks"> </span><span class="n" id="mwAkw">i</span><span class="w" id="mwAk0"> </span><span class="k" id="mwAk4">in</span><span class="w" id="mwAk8"> </span><span class="n" id="mwAlA">axes</span><span class="p" id="mwAlE">(</span><span class="n" id="mwAlI">z</span><span class="p" id="mwAlM">,</span><span class="mi" id="mwAlQ">1</span><span class="p" id="mwAlU">)</span> <span class="w" id="mwAlY"> </span><span class="n" id="mwAlc">z</span><span class="p" id="mwAlg">[</span><span class="n" id="mwAlk">i</span><span class="p" id="mwAlo">,</span><span class="o" id="mwAls">:</span><span class="p" id="mwAlw">]</span><span class="w" id="mwAl0"> </span><span class="o" id="mwAl4">.=</span><span class="w" id="mwAl8"> </span><span class="n" id="mwAmA">sincospi</span><span class="p" id="mwAmE">(</span><span class="mi" id="mwAmI">2</span><span class="w" id="mwAmM"> </span><span class="o" id="mwAmQ">*</span><span class="w" id="mwAmU"> </span><span class="n" id="mwAmY">rand</span><span class="p" id="mwAmc">());</span> <span class="w" id="mwAmg"> </span><span class="n" id="mwAmk">z</span><span class="p" id="mwAmo">[</span><span class="n" id="mwAms">i</span><span class="p" id="mwAmw">,</span><span class="o" id="mwAm0">:</span><span class="p" id="mwAm4">]</span><span class="w" id="mwAm8"> </span><span class="o" id="mwAnA">.*=</span><span class="w" id="mwAnE"> </span><span class="n" id="mwAnI">sqrt</span><span class="p" id="mwAnM">(</span><span class="o" id="mwAnQ">-</span><span class="mi" id="mwAnU">2</span><span class="w" id="mwAnY"> </span><span class="o" id="mwAnc">*</span><span class="w" id="mwAng"> </span><span class="n" id="mwAnk">log</span><span class="p" id="mwAno">(</span><span class="n" id="mwAns">rand</span><span class="p" id="mwAnw">()));</span> <span class="w" id="mwAn0"> </span><span class="k" id="mwAn4">end</span> <span class="w" id="mwAn8"> </span><span class="n" id="mwAoA">vec</span><span class="p" id="mwAoE">(</span><span class="n" id="mwAoI">z</span><span class="p" id="mwAoM">)</span> <span class="k" id="mwAoQ">end</span> <span class="s" id="mwAoU">"""</span> <span class="s" id="mwAoY"> boxmullersample(n,μ,σ)</span> <span class="s" id="mwAoc">Generate `n` samples from the normal distribution with mean `μ` and standard deviation `σ` using the Box-Muller method.</span> <span class="s" id="mwAog">"""</span> <span class="k" id="mwAok">function</span><span class="w" id="mwAoo"> </span><span class="n" id="mwAos">boxmullersample</span><span class="p" id="mwAow">(</span><span class="n" id="mwAo0">n</span><span class="p" id="mwAo4">,</span><span class="n" id="mwAo8">μ</span><span class="p" id="mwApA">,</span><span class="n" id="mwApE">σ</span><span class="p" id="mwApI">)</span> <span class="w" id="mwApM"> </span><span class="n" id="mwApQ">μ</span><span class="w" id="mwApU"> </span><span class="o" id="mwApY">.+</span><span class="w" id="mwApc"> </span><span class="n" id="mwApg">σ</span><span class="o" id="mwApk">*</span><span class="n" id="mwApo">boxmullersample</span><span class="p" id="mwAps">(</span><span class="n" id="mwApw">cld</span><span class="p" id="mwAp0">(</span><span class="n" id="mwAp4">n</span><span class="p" id="mwAp8">,</span><span class="mi" id="mwAqA">2</span><span class="p" id="mwAqE">))[</span><span class="mi" id="mwAqI">1</span><span class="o" id="mwAqM">:</span><span class="n" id="mwAqQ">n</span><span class="p" id="mwAqU">];</span> <span class="k" id="mwAqY">end</span>