Am 25.08.2010 06:13, schrieb Anthony D. Baye:

> Second my random number function. (I've attached the source file, but I'll
> reiterate the function here) I'm using the rand() function from cstdlib:
>
> double Random(double r)
> {
>       double rNum;
>
>       rNum = r - (2*r)/(rand()%65535);
>
>       return rNum;
> }

What type of random distribution do you want to achieve? If you're 
attempting to get a uniform distribution in the interval -r to +r (i.e. 
each value in this interval is equally likely), then you're pretty much 
off the track.

First, note that rand() /is/ uniformly distributed, in the interval 0 to 
RAND_MAX (which, BTW, is typically 32767, but you shouldn't rely on that).

Now, taking this number modulo 65535 doesn't change a thing, provided 
that RAND_MAX < 65535. Otherwise, it will give you a uniformly 
distributed value in the range 0 to 65534 (sic!).

You take the inverse of this value, which gives you a /non-uniformly/ 
distributed number. Not sure from the top of my hat /what/ distribution 
this gives you, but I think it's an exponential one, with the value x 
having a probability of c*e^-x. At any rate, the range is also modified. 
In a nutshell, at this point you get random values in the range 
1/RAND_MAX (about 0.00003) to 1/0 (1.#INF, positive infinity), and while 
each individual value in practice still has roughly the same 
probability, the values clump together near 0.0 (for instance, the 
highest values are ..., 0.20, 0.25, 0.333, 0.50, 1.00, 1.#INF).

The final operations - multiplying by 2*r and subtracting the result 
from r - just scale and shift this non-uniform distribution.


What you probably really want is

     rNum = r - (2*r) * ((double)rand() / (double)RAND_MAX);

This will give you a uniformly distributed value in the range -r to +r 
(excluding +r). Or, to get rid of the small bias towards -r:

     rNum = r - (2*r) * ((rand() + 0.5) / (double)RAND_MAX);


OTOH, if the exponential distribution is intentional, then what you want is:

     rNum = r / (1 - 2 * ((rand() + 0.5) / (double)RAND_MAX));

This will give you a uniformly distributed value in the range -1.#INF to 
+1.#INF (actually -r*RAND_MAX to +r*RAND_MAX if I'm not mistaken), 
clustering around 0.0.


>       srand(65535);

Are you sure you know what seed() does? The parameter, though not being 
illegal, looks suspicious for a random seed.

Post a reply to this message