Implement Box-Mueller random number generator in C #

Question

From this question: a random number generator that gravitates numbers to any given number in the range? I did some research, since I came across such a random number generator before, All I remember was the name "Muller", so I think I found it here:

• Box-Mueller Conversion

I can find many of its implementations in other languages, but I cannot correctly implement it in C #.

This page, for example, the Box-Muller Method for generating Gaussian random numbers, says the code should look like this (this is not C #):

#include <stdio.h> #include <stdlib.h> #include <math.h> #include <time.h> double gaussian(void) { static double v, fac; static int phase = 0; double S, Z, U1, U2, u; if (phase) Z = v * fac; else { do { U1 = (double)rand() / RAND_MAX; U2 = (double)rand() / RAND_MAX; u = 2. * U1 - 1.; v = 2. * U2 - 1.; S = u * u + v * v; } while(S >= 1); fac = sqrt (-2. * log(S) / S); Z = u * fac; } phase = 1 - phase; return Z; } 

Now, here is my implementation above in C #. Note that the conversion gives 2 numbers, hence the trick with the “phase” above. I just discard the second value and return the first.

 public static double NextGaussianDouble(this Random r) { double u, v, S; do { u = 2.0 * r.NextDouble() - 1.0; v = 2.0 * r.NextDouble() - 1.0; S = u * u + v * v; } while (S >= 1.0); double fac = Math.Sqrt(-2.0 * Math.Log(S) / S); return u * fac; } 

My question is in the following specific scenario, where my code does not return a value in the range from 0 to 1, and I cannot figure out what the source code might look like.

• u = 0.5, v = 0.1
• S becomes 0.5*0.5 + 0.1*0.1 = 0.26
• fac becomes ~ 3.22
• Thus, the return value is ~ 0.5 * 3.22 or ~ 1.6

It is not within 0 .. 1 .

What am I doing wrong / don't understand?

If I change my code so that instead of multiplying fac by u , I multiply by S , I get a value that ranges from 0 to 1 but has an incorrect distribution (it seems to have a maximum distribution around 0.7-0.8, and then tapers in both directions.)

Answer 1

Your code is ok. Your mistake is that it should return values ​​exclusively within [0, 1] . A normal distribution (standard) is a distribution with a non-zero weight over the entire real line. That is, values ​​outside of [0, 1] possible. In fact, the values ​​inside [-1, 0] just as probable as the values ​​within [0, 1] , and, in addition, the addition [0, 1] has about 66% of the weight of the normal distribution. Therefore, 66% of the time we expect a value outside of [0, 1] .

In addition, I think this is not a Box-Mueller transformation, but it is actually a polar Marsaglia method.

Answer 2

I am not a mathematician or statistician, but if I think about it, I would not expect the Gaussian distribution to return numbers in the exact range. Given your implementation, the average is 0 and the standard deviation is 1, so I expect the values ​​distributed along the call curve from 0 in the center and then decrease as numbers deviate from 0 on each side. Thus, the sequence will definitely cover as +/- numbers.

Then, since it is statistical, why is it difficult to limit it to -1.1 just because std.dev is 1? Statistically there can be some kind of game on both sides and still fulfill the statistical requirements.

Answer 3

The uniform random variable is indeed within 0..1, but the Gaussian random variable (which the Box-Muller algorithm generates) can be anywhere on the real line. See wiki / NormalDistribution for more details .

Answer 4

I think the function returns polar coordinates. Therefore, you need both values ​​to get the correct results.

In addition, the Gaussian distribution is not between 0 .. 1 . This can easily end as 1000, but the likelihood of such an occurrence is extremely low.

Answer 5

This is the monte carlo method, so you cannot pinch the result, but what you can do is ignore the samples.

 // return random value in the range [0,1]. double gaussian_random() { double sigma = 1.0/8.0; // or whatever works. while ( 1 ) { double z = gaussian() * sigma + 0.5; if (z >= 0.0 && z <= 1.0) return z; } }