3
votes

What's the best way to generate pseudo-random numbers in the closed interval [0,1] instead of the usual [0,1)? One idea I've came up with is to reject values in (1/2,1), then double the number. I wonder if there is a better method.

real x
do
   call random_number(x)
   if (x <= 0.5) exit
end do
x = 2*x
print *, x
end

The most important requirement is that the algorithm should not make a worse distribution (in terms of uniformity and correlation) than that generated by random_number(). Also I'd favour simplicity. A wrapper around random_number() would be perfectly good, I'm not looking to implement a whole new generator.


As @francescalus points out in the comments, with the algorithm above lots of numbers in [0,1] will have zero probability of appearing. The following code implements a slightly different approach: the interval is enlarged a bit, then values in excess of 1 are cut out. It should behave better in that aspect.

real x
do
   call random_number(x)
   x = x*(1 + 1e-6)
   if (x <= 1.) exit
end do
print *, x
end
1
"Best" will be too broad, so perhaps you can explain what properties you require in this new distribution. For example, in the case of the question, lots of numbers in the interval [0,1] now have zero probability of appearing. - francescalus
@francescalus The most important requirement is that it should not "ruin" the distribution generated by random_number(). I'm not looking to implement a whole new generator, just something to wrap random_number() in. I'd like it to be concise. Judging by the first requirement, the second algorithm I've added should be better, if I understood the problem with the first one correctly. - Arch Stanton
@ArchStanton The wiki page mentions random numbers generated according to U(0,1), this does not mean that 1 has to be included. I also think the page is a bit messy with the limits. If you look in numerical recipes it is nicely explained. In short, you can use random_number for the Box-Muller transform without any problem. You will have a very nice Gaussian distribution. - kvantour
The cumulative distribution functions of the 4 uniform distributions U((0,1)), U([0,1)), U((0,1]) and U([0,1]) are identical, so you are worried about something which is virtually nothing. - John Coleman
@kvantour It wouldn't hurt to include 1 (at the cost of excluding some other representable real in [0,1) ), but any code that made explicit use of its inclusion (e.g. by explicitly testing for equality with 1) is probably bad code (IMHO). - John Coleman

1 Answers

4
votes

What about swapping x and 1-x? Sorry, my Fortran is rusty

real function RNG()
real ::    x
logical, save :: swap = .TRUE.

call random_number(x)
if (swap .EQV. .TRUE.)
    RNG = x
    swap = .FALSE.
else
    RNG = 1.0 - x
    swap = .TRUE.
end if

end

And if you want to use Box-Muller, use 1-U everywhere and it should work

z0 = sqrt(-2.0*log(1.0-U1))*sin(TWOPI*U2)
z1 = sqrt(-2.0*log(1.0-U1))*cos(TWOPI*U2)

same for rejection version of Box-Muller