Logarithm of a random number

Stef1611 0
Viewed 81

I wrote a fortran program to code this algorithm (https://en.wikipedia.org/wiki/Reservoir_sampling#Algorithm_A-ExpJ).
It works on my computer. But after I asked these two questions (Intrinsic Rand, what is the interval [0,1] or ]0,1] or [0,1[ and How far can we trust calculus with infinity?), I think I could have a problem with log(random()) because call random_number(Xw); Xw = log(Xw) is used. Indeed, random_number(Xw) could return 0 and log(0)=-infinity.

Therefore, I plan to modify this line as follow call random_number(Xw); Xw = log(1-Xw) to change the random value interval from [0,1[ to ]0,1].

Is it a good idea or is there a best solution ?

infinity fortran random
2 Answers

Not a good idea. If you want your algorithm to be stable, you need to define bounds. your log function represents a priority, it can likely be as low as you want, but it must be a number. You can bind it to numerical precision:

program t
   use iso_fortran_env
   implicit none
     real(real64), parameter :: SAFE = exp(-0.5d0*huge(0.0_real64))

     print *, log(randoms_in_range(100,SAFE,1.0_real64))

   contains

      elemental real(real64) function in_range(f,low,hi) result(x)
         real(real64), intent(in) :: f ! in: [0:1]
         real(real64), intent(in) :: low,hi
         real(real64) :: frac
         frac = max(min(f,1.0_real64),0.0_real64)
         x = low+frac*(hi-low)
      end function in_range

      real(real64) function random_in_range(low,hi) result(x)
         real(real64), intent(in) :: low,hi
         call random_number(x)  ! [0,1]
         x = in_range(x,low,hi) ! [low,hi]
      end function random_in_range
      function randoms_in_range(n,low,hi) result(x)
         integer     , intent(in) :: n
         real(real64), intent(in) :: low,hi
         real(real64)             :: x(n)
         call random_number(x)  ! [0,1]
         x = in_range(x,low,hi) ! [low,hi]
      end function randoms_in_range
end program
Federico Perini0

While mathematically it is true that if X is uniformly distributed on (the real interval) [0,1) then 1-X is uniformly distributed on (0,1], this does not particularly help you.

As noted in the description of the algorithm to which you link, the underlying assumption is that the base uniform distribution is over the interval (0,1). This is not the same as (0,1].

You can use rejection sampling to generate X uniformly over (0,1) from random_number() (which is [0,1)): throw away all zero occurrences.

francescalus0
Related