2
$\begingroup$

I've implemented the Marsaglia polar method to generate random numbers that are normally distributed. Unfortunately the code as shown on this website returns numbers that are not within the range [0, 1).

Code:

     double x1, y1, w;

     do {
             x1 = 2.0 * ranf() - 1.0;
             y1 = 2.0 * ranf() - 1.0;
             w = x1 * x1 + y1 * y1;
     } while ( w >= 1.0 || w == 0 );

     w = sqrt( (-2.0 * ln( w ) ) / w );
     return x1 * w;

where ranf() returns values [0, 1).

For example, if (x1, y1) was (-0.43458512585358, -0.07521858582050478), this returns a value of -1.7830255550765148. Obviously not a value within my expected range.

I've seen some implementations that multiplies this by the standard deviation and adds the mean. But if I want to get numbers that range from [0, 1), what should I use as my standard deviation and mean? currently, I'm using a standard deviation of 1/8.125 ~= 0.121212... and 0.5 as the mean, but I've only stumbled on this by experimentation. What is the official way for me to get this appropriately ranged?

To answer the question in the comments below: I want to generate normal random variables within the range 0 and 1. Obviously I'd like the mean to be at 0.5 and values no be distributed around this mean.

More for my enlightenment, will the polar method really return numbers over the entire real line? If so, then it should be easy for me to scale the result down to [0, 1), assuming I don't hit floating point underflow due to the division to scale things down.

New edit: What I'm trying to do is take Aniko's advise in his answer to my other Power Factor Problem question. So my first thought was that simulate my distribution of velocities, was to implement my on GaussianRandom class by deriving from Random and overriding Random.Sample() and the 3 other methods as recommended by MSDN. Since Random.Sample() is suppose to return a double in the range [0, 1), then I've hit this stumbling block.

I've come to the same conclusion as mbq's answer after staring at the ln(x) graph where x is in (0, 1] and sqrt(x) where x is in (0, 2].

From a practical/programming standpoint, I can't override Random.Sample(), because later when I actually call GaussianRandom.NextDouble(), there is no way for me to pass in the actual standard deviations for my velocities. It looks like I'll have to use the Marsaglia polar method in the raw instead of wrapped within the Random class. My "cartridge" generator will have to return x1 * w * stdDev + mean.

$\endgroup$
  • $\begingroup$ Do you want to generate normal random variables truncated to range between 0 and 1 or uniform random variables between 0 and 1? The Marsagla polar method will return numbers over the entire real line as it generates standard normal random variables. $\endgroup$ – user28 Aug 25 '10 at 9:46
  • $\begingroup$ I don't know if this solves your problem, but you might wanna take a look at this question where some transformations that squash everything between 0 and 1 are discussed: stats.stackexchange.com/questions/1112/… $\endgroup$ – Henrik Aug 25 '10 at 10:40
  • 1
    $\begingroup$ @Ants: "Scaling down" from the real line to an interval will necessarily make the resulting generator non-Normal. $\endgroup$ – whuber Aug 25 '10 at 18:14
  • $\begingroup$ @Henrik: Thanks. I'll have to remember some of those tricks for future reference. $\endgroup$ – Ants Aug 26 '10 at 0:41
  • 1
    $\begingroup$ Your clarification completely confused me. You plan to generate a random normal variable x1 truncated to the 0-1 range, then plug it into the x1 * w * stdDev + mean formula (what is w?) to get a normal with mean mean and standard deviation stdDev? If so, that's totally wrong. Just take a standard normal z, and calculate mean + z*stdDev to get a non-standard normal variable. $\endgroup$ – Aniko Aug 26 '10 at 13:01
6
$\begingroup$

It is just impossible, since normal distribution has nonzero probability for the whole $R$. Of course because of the numerical issues it would be not the case for Marsaglia method (this obviously means that it is not creating normally distributed numbers), yet this numerical range is so huge that it is nonsense to scale it to [0,1].

EDIT: So, as long as I understand, you want this number to be in [0,1) to overload .net Random; you don't need to! Write your own class, say NormalRandom which will hold one Random object and will use its output to generate (using Marsaglia) normal random numbers from $-\infty$ to $\infty$ and just use them. Or rather use something already done; Google found http://www.codeproject.com/KB/recipes/Random.aspx .
On margin, remember that when it is said that something which cannot be negative is "normally distributed" (for instance mass) it does mean that it is not from normal but from "almost-normal" distribution which looks pretty normal, but still has one tail cut out (nevertheless simple rejection of negative outputs from generator is a fair way to deal with it).

$\endgroup$
  • $\begingroup$ Thanks for the pointer. LOL! My new implementation for NormalDistribution looks a lot like I'd taken a slice of that code project and focused on just normal random numbers. I've got the same abstract classes for RNG's and Distribution's. I've got just one concrete RNG that uses the .NET Random class, and a concrete NormalDistribution class. $\endgroup$ – Ants Aug 26 '10 at 0:48
  • $\begingroup$ But as I can see they are not overloading Microsoft's original, so they can break the [0,1) limit. $\endgroup$ – user88 Aug 26 '10 at 6:55
  • $\begingroup$ Yes, that is correct. $\endgroup$ – Ants Aug 26 '10 at 8:30
4
$\begingroup$

You also might want to consider that perhaps you don't want a normally distributed variable on [0,1) but rather a random variable with a bell-shaped distribution naturally restricted to the [0,1) interval. If so, consider the Beta distribution or the Kumaraswamy distribution.

$\endgroup$
  • $\begingroup$ I definitely want a normally distributed variable. I'll update the question with the specific problem I'm trying to solve. $\endgroup$ – Ants Aug 25 '10 at 20:08
2
$\begingroup$

Another possible interpretation is that you are confusing "normal" and "uniform". A uniform variate will have a mean of 0.5 and be evenly distributed on [0, 1). If this interpretation is correct, your code simplifies to

return ranf()
$\endgroup$
  • $\begingroup$ I want to get a "normal" distribution with the appropriate hump at 0.5. $\endgroup$ – Ants Aug 25 '10 at 21:03
  • 1
    $\begingroup$ Then you won't be able to restrict it to the interval [0, 1). As a practical matter, normal probabilities further out than about seven standard deviations are too small to be worth computing, anyway, so you could set the standard deviation to 1/2 / 7 = 1/14 = about 0.07. This isn't too far from the value you found empirically. With your value and enough draws, you are likely to get a number outside [0, 1); with 1/14, it's extremely unlikely. Thus, my recommendation is to take a value from the Marsaglia algorithm, divide by 14, and add 1/2. Adjust the 1/14 to suit your application. $\endgroup$ – whuber Aug 26 '10 at 1:08
1
$\begingroup$

Based on the update you want to generate truncated normal random variables. Scaling based on the maximum and minimum will 'destroy' the normal distribution properties of the draws. Instead, you should consider using one of the two methods described below:

Rejection Sampling

You draw from the standard normal and then accept the draw only if it is between 0 and 1.

Inverse Transform Sampling

Let:

$\phi(x)$ be the standard normal cdf.

Then a standard normal draw that is restricted to lie between 0 and 1 is given by:

$x = \phi^{-1}(\phi(0) + (\phi(1) - \phi(0)) U )$

where

$U$ is a uniform random draw between 0 and 1.

$\endgroup$
  • $\begingroup$ For rejection sampling, what will I use as my g(x) and M value to simulate the normal distribution? For inverse transform sampling, do I use D Ibbetson's Algorithm 209 (old.sigchi.org/~perlman/stat/doc/z.c) to compute the cdf? What do I use as the inverse? $\endgroup$ – Ants Aug 25 '10 at 20:21
  • $\begingroup$ @Ants If I am right, you don't need a cannon like ITS. $\endgroup$ – user88 Aug 25 '10 at 21:25
  • $\begingroup$ For rejection sampling set f(x) ~ standard normal pdf I(0<x<1), g(x) ~ standard normal and and M=1. Your condition will then become select x if u < I(0<x<1) which is basically what I stated in my one line sentence above under rejection sampling. $\endgroup$ – user28 Aug 25 '10 at 21:34
  • $\begingroup$ @mbq: I don't get the ITS reference. I've been missing out on a lot of pop culture and current events references the past few months. $\endgroup$ – Ants Aug 26 '10 at 0:37
  • $\begingroup$ @Ants ITS=Inverse Transform Sampling $\endgroup$ – user88 Aug 26 '10 at 6:53

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.