Talk:Box–Muller transform

the C++ code

I don't understand why epsilon need be brought in. The test for u1 being less than or equal to epsilon could be replaced by a test for u1 being greater than 0. Keep it simple! Jswd (talk) 05:53, 10 December 2023 (UTC)

On my x86-64 system, epsilon is 2.22045e-16 which you think of as the granularity of 64-bit doubles around 1. The smallest non-zero number usefully available is MUCH smaller - about 1e-308. Attempting to calculate the logarithm of zero will fail. However, 1e-308 has a decimal logarithm of -308 and a natural logarithm of about -710 which can be handled without any error. Jswd (talk) 19:44, 20 December 2023 (UTC)

Nice article

Nice article and clearly written. — Preceding unsigned comment added by 2001:4643:E6E3:0:65F3:23E8:7C14:9334 (talk) 13:28, 4 October 2019 (UTC)

Untitled

Michael, thanks for catching my mistake regarding phi & varphi. Strange that the html markup & phi should be \varphi in tex. Any ideas on correctly numbering the two methods? --snoyes 22:28 Feb 20, 2003 (UTC)

I wonder wether those two numbers will be linked or independent from each other? May anyone tell? -- Feb 23, 2004

they are linked; they both transform a uniform distribution over the unit circle to a gaussian distribution; the difference is that the first is given in polar coordinates, while the second is in cartesian coordinates

JeffBobFrank 04:29, 4 Mar 2004 (UTC)

---

The second method is typically faster because it uses only one transcendental function instead of three
I can only see the need for two transcendental functions in the first method, since $\sin ^{2}(x)+\cos ^{2}(x)=1$. I assume ×π doesn't count. --Henrygb 14:36, 20 Jul 2004 (UTC)

---

The discussion after the formulas says we get out uniformly distributed values, but I thought we wanted normally distributed ones. Also, I still don't quite understand the characteristics of the two formulas. The first seems to go from polar coordinates to... cartesian? More polar coords? The second seems to go from cart -> cart, but what range can I expect? The discussion suggests output values are in a unit circle, but that's not what I get (fortunately, since with a stddev of 1 it'd cut off the tails). Lunkwill 02:18, 16 Oct 2004 (UTC)

getting the uniformly distributed numbers is one step in both methods. they're then multiplied by the thing in the square root, which gives them a normal distribution. and the first one goes from polar to cart, the second from cart to cart. Frencheigh 02:46, 16 Oct 2004 (UTC)

https://archive.org/details/bitsavers_mitreESDTe69266ANewMethodofGeneratingGaussianRando_2706065 — Preceding unsigned comment added by 123.16.159.239 (talk) 01:12, 17 March 2016 (UTC)

Be Careful with the Second Method ("Polar form")

Warning: Be very careful with the second Box-Muller method given here ("Polar form"). It may not be correct as written. It would be correct if the quantities labelled as "r" under the square root signs were replaced by r-squared, because r-squared is uniformly distributed in (0,1], x/r is a cosine, and x/r is a sine. TopKatz 22:48, 9 October 2005 (UTC)

r is defined here as $x^{2}+y^{2}$. You were probably expecting $r$ to be the norm of $(x,y)$. Deco 16:53, 26 May 2006 (UTC)

The Second Paragraph is in Error

The basic method takes points uniformly distributed in the unit square (not circle) and transforms them as given by the formulas listed. The polar method takes points uniformly distributed in the unit circle (note that you throw away all points with a radius or radius squared greater than 1) and transforms them. As noted above, R is the radius-squared, which is uniform (0,1]. The polar coordinates, (r, theta) are not uniformly distributed. Theta is uniformly distributed, but r is not (r-squared is uniform). In my previous sentence, I used the standard notation r for radius, not r for a uniform variable as is done in the article.

It might be helpful the write the polar form in a way that corresponds to the basic form, i.e, x/sqrt (R) is the cosine of a uniform angle ( R is the radius-squared again) and sqrt (-2 ln R) is equivalent to sqrt (-2 ln r), since both R, the radius squared, and r, which has nothing to do with a radius, are uniform (0,1]. The formulas as given should remain; I am suggesting something of the form

z = ( x/ sqrt (b)) * sqrt (c) = a sqrt (b/c)

(I have never edited before and don't feel comfortable editing the article itself)

72.134.75.20 00:49, 21 September 2006 (UTC) Al

correct myself

z = x/ sqrt (b)) * sqrt (c) = x sqrt (c/b)

72.134.75.20 00:57, 21 September 2006 (UTC) Al

Rejection probabiility

Hello,

I screwed up, so reverted, just to clear it uo. The Box-Muller transform has, as the output from 2 uniform deviates (UDs), 2 random numbers (Knuth). The acceptance rate for input UD pairs being pi*r^2 / (2r)^2 = pi/4 therefore the rejection rate is (1-pi/4) = 0.214. The acceptance being pi/4, the numbers generated per output are 4/pi. reverting my bad... User A1 04:53, 19 April 2007 (UTC)

Numerically more stable?

The article, in the section on comparing the two forms states: "it is typically 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"

I do not believe that Knop's polar method is more numerically stable than the basic transform. The cited reference does not convince me; it is not from a numerical analyst but from someone writing for a general audience of FORTH programmers. It doesn't give a proof or any reference. All it does is just state that the algorithm is more robust.

A hand-wavy argument about why it may be less stable is as follows:

In both forms of the algorithm, the argument of ln is uniformly distributed. Thus it is small with equal probability.

In the polar form, sqrt((-ln s)/s) is multiplied by a number in the interval 0..1. In the original form sqrt(-ln s) is multiplied by a number in the interval 0..1. The only real difference in these two is the difference between ln s/s and ln s. (ln s)/s becomes very negative much faster than ln s as s goes to 0. And it is always true that ln s/s < ln s on the interval (0,1) Thus if a certain value s is too close to 0 in ln s. Then ln s/s will be worse. Further ln s/s has a larger slope meaning that little errors in the input will be more exaggerated in the output.

Neither function should have any problem near 1 since the slope of both ln s/s and ln s is pretty horizontal there, meaning that small errors in the input yield even smaller errors in the output.

Thus, I believe that contrary to the current form of the article that the original box muller is better in terms of stability. The degree of this increase and whether this outweighs (possible) speed trade-offs is an issue the implementer must deal with. But I don't think the article should state "more robust" without some better proof.

BrotherE (talk) 21:58, 2 May 2011 (UTC)

Mapping diagram

Hi all, I removed the mapping diagram I made from the intro after receiving the following e-mail about it:

If I'm right you are the author of this (File:Box_Muller.svg) picture on Wikipedia. I'm using Box-Muller Transform in my student's project, and at the moment I'm trying to understand how this circles are created. On the first graph you have circles of uniform distribution, so:

r=sqrt(u^2+v^2)

and r is from 0.1 to 1.0 with step 0.1.

s=r^2

When we look to the second graph you use this function:

R=sqrt(-2*ln(r)/r).

Could you tell me why it is so simple? When we look at the equations:

z0 = u * sqrt(-2*ln(r)/r)

z1 = v * sqrt(-2*ln(r)/r)

Because of this, I thought that R should be calculated as:

R = sqrt(z0^2 + z1^2) = sqrt(u^2+v^2)*sqrt(-2ln(r)/r))

Why we can remove this "first" element = sqrt(u^2+v^2)?

I think they're right, my diagram appears to be incorrect. Unfortunately I'm not entirely sure how to correct it. I'd like to somehow visualize the two-dimensional mapping properly. Any suggestions would be welcome. Dcoetzee 05:12, 27 October 2011 (UTC)

Incorrect implementation

I don't know who did it but the implementation seems incorrect (unused variables, return with uninitialized variables). I'm not a mathematician so I cannot correct the code myself. — Preceding unsigned comment added by 82.225.160.248 (talk) 16:21, 26 June 2016 (UTC)

The thread_locals behave like statics, but are thread local, so nothing's uninitialised, your worries are unfounded. The u1 rejection loop is a bit dumb, one could simply use (rand()+1)*U1_SCALE_FACTOR instead, where U1_SCALE_FACTOR would be 1./(1.+RAND_MAX), another static const. — Preceding unsigned comment added by 2001:7D0:8310:5B80:DACB:8AFF:FEA7:FCDD (talk) 07:48, 8 September 2017 (UTC)

Value of sigma?

The Basic Method section gives a useful formula, but I'm not clear what the standard deviation of the resulting distribution will be. The link to the standard normal distribution article suggests the standard deviation is 1.0, but I'd rather see the article say this explicitly instead of letting me figure it out. It's been a while since I've looked at these distributions, and my math is very rusty. —MiguelMunoz (talk) 06:16, 5 May 2021 (UTC)