Верхняя граница генератор случайных чисел

Question

Это фактически прослеживание вопрос о предыдущей: Rounding of double precision to single precision: Forcing an upper bound

После того, что я думал, было решение моих проблем с ответом предыдущий вопрос, я попытался запустить свою программу еще раз и обнаружил, что у меня такая же проблема.

Реализация Mersenne Twister, которую я использую, генерирует подписанное 32-битное случайное целое число. Парень, который реализовал ГСЧ сделал эту функцию, чтобы генерировать случайные двойной точности с плавающей точкой в ​​диапазоне [0,1):

function genrand_real2() 
    double precision genrand_real2,r 
    integer genrand_int32 
    r=dble(genrand_int32()) 
    if(r.lt.0.d0)r=r+2.d0**32 
    genrand_real2=r/4294967296.d0 
    return 
    end 

И это работает безупречно, так что следующее предложение в предыдущем вопросе я использовал следующее функция для генерации случайных одинарной точности с плавающей точкой, в диапазоне я подумал бы [0,1):

function genrand_real() 
    real genrand_real, r 
    integer genrand_int32 
    r = real(genrand_int32()) 
    if (r .lt. 0.0) r = r + 2.0**32 
    genrand_real = r/4294967296.0 
    return 
    end 

Однако я получил ту же ошибку я получил перед тем, вызванным 1,0 числа. Поэтому я написал небольшую программу, чтобы показать, что мой genrand_real на самом деле генерирует 1.0 и обнаружил, что я прав, а 1.0 создан. Это приводит к тому, что я использую для генерации целого числа в диапазоне [1, MAX] (в этом примере [1,5]), чтобы сгенерировать значение MAX + 1, среди других неудобств в коде, над которым я работаю.

i = 0 
    do while (.true.) 
    r = genrand_real() 
    if (r .gt. 0.99999) then 
     i = i + 1 
     print *, 'number is:', r 
     print *, 'conversion is: ', int(5*r)+1 
    endif 
    if (i .gt. tot_large) exit 
    enddo 

Вопрос в том, почему он работает для двойной точности, но не для одиночного прецизионного поплавка? Я не вижу причины, по которой он терпит неудачу, поскольку 2 ** 32 помещается в один прецизионный поплавок. Кроме того, что я должен сделать, чтобы исправить это? Я думал о делении числа на 2.0 ** 32 + 1 вместо 2.0 ** 32, но я не уверен, что это теоретически правильно и что числа будут одинаковыми.

Answer 1

Я не уверен, следует ли публиковать этот ответ по старому вопросу или здесь. В любом случае у меня может быть решение (во втором блоке кода).

Рутина, что я использовал для одной и той же задачи, так как около двух лет назад это:

function uniran() 
    implicit none 
    integer, parameter :: dp = selected_real_kind(15, 307) 
    real(dp) :: tmp 
    real :: uniran 
    tmp = 0.5_dp + 0.2328306e-9_dp * genrand_int32() 
    uniran = real(tmp) 
end function uniran 

Я забыл, где код от и всегда, хотя это просто, но есть тонкий трюк для него , который я только сейчас понял. Очевидное различие заключается в умножении вместо деления, но это только потому, что быстрее умножаться с фиксированным числом, чем разделить (0.2328306e-9 = 1/4294967296).
Трюк: это не так. 1/4294967296 = 0.23283064365386962890625e-9, поэтому программа использует менее значимые цифры, чем может удерживать двойная точность (15, а всего 7). Если вы увеличиваете количество цифр, итоговое число приближается к 1 и становится точно одним во время последующего преобразования. Вы можете попробовать: если вы используете только одну цифру, она начинает сбой (= 1.0). Видимо, это решение несколько рубить, поэтому я также попробовал другой подход, передискретизацию, если результат точно 1:

recursive function resample_uniran() result(res) 
    implicit none 
    integer, parameter :: dp = selected_real_kind(15, 307) 
    real(dp) :: tmp 
    real :: res 
    tmp = 0.5_dp + 0.23283064365386962890625e-9_dp * genrand_int32() 
    res = real(tmp) 
    if (res == 1.0) then 
     res = resample_uniran() 
    end if 
end function resample_uniran 

Я написал программу, которая проверяет функцию (модуль, который содержит функции и подпрограммы в конце поста, это относительно длинный):

program prng_fail 
use mod_prngtest 
implicit none 
integer(kind=16) :: i, j, k 

! loop counters 
i = 0 
j = 0 
k = 0 

call init_genrand_int32() 

do 
    i = i + 1 
    j = j + 1 
    k = k + 1 
    if (genrand_real() == 1.0) then 
     print*, 'genrand_real fails after ', i, ' iterations' 
     i = 0 
    end if 
    if (uniran() == 1.0) then 
     print*, 'uniran fails after ', j, ' iterations' 
     j = 0 
    end if 
    if (resample_uniran() == 1.0) then 
     print*, 'resample_uniran fails after ', k, ' iterations' 
     k = 0 
    end if 
end do 

end program prng_fail 

с результатом, что genrand_real терпит неудачу (= 1,0) часто (мы говорим каждые несколько миллионов номеров), в то время как две другие имеют до сих пор никогда не проваливался. Рекурсия-версия стоит вам времени, но технически лучше, потому что максимально возможное число ближе к 1.

Я также тестировал скорость и «единообразие» и сравнивал с внутренней функцией random_number, что также дает равномерные случайные числа в [0,1]. (Осторожно, это создает 3 х 512 МБ файлов)

program prng_uniformity 
use mod_prngtest 
implicit none 
integer, parameter :: n = 2**27 
real, dimension(n) :: uniran_array, resamp_array, intrin_array 
integer :: array_recl, i 
real :: start_time, end_time 

call init_genrand_int32() 
call init_random_seed() 

! first check how long they take to produce PRNs 
call cpu_time(start_time) 
do i=1,n 
    uniran_array(i) = uniran() 
end do 
call cpu_time(end_time) 
print*, 'uniran took ', end_time - start_time, ' s to produce ', n, ' PRNs' 

call cpu_time(start_time) 
do i=1,n 
    resamp_array(i) = resample_uniran() 
end do 
call cpu_time(end_time) 
print*, 'resamp took ', end_time - start_time, ' s to produce ', n, ' PRNs' 

call cpu_time(start_time) 
do i=1,n 
    call random_number(resamp_array(i)) 
end do 
call cpu_time(end_time) 
print*, 'intrin took ', end_time - start_time, ' s to produce ', n, ' PRNs' 

! then save PRNs into files. Use both() to have the same random 
! underlying integers, reducing the difference purely to 
! the scaling into the interval [0,1) 
inquire(iolength=array_recl) uniran_array 
open(11, file='uniran.out', status='replace', access='direct', action='write', recl=array_recl) 
open(12, file='resamp.out', status='replace', access='direct', action='write', recl=array_recl) 
open(13, file='intrin.out', status='replace', access='direct', action='write', recl=array_recl) 
do i=1,n 
    call both(uniran_array(i), resamp_array(i)) 
    call random_number(intrin_array(i)) 
end do 
write(11, rec=1) uniran_array 
write(12, rec=1) resamp_array 
write(13, rec=1) intrin_array 

end program prng_uniformity 

Результаты всегда одинаковы в принципе, даже если тайминги различно:

uniran took 0.700139999  s to produce 134217728 PRNs 
resamp took 0.737253010  s to produce 134217728 PRNs 
intrin took 0.773686171  s to produce 134217728 PRNs 

uniran быстрее, чем resample_uniran, который быстрее, чем внутренняя (хотя это в значительной степени зависит от PRNG, Mersenne twister будет медленнее, чем внутренняя).

Я также посмотрел на выходе каждый метод обеспечивает (с Python):

import numpy as np 
import matplotlib.pyplot as plt 

def read1dbinary(fname, xdim): 
    with open(fname, 'rb') as fid: 
     data = np.fromfile(file=fid, dtype=np.single) 
    return data 

if __name__ == '__main__': 
    n = 2**27 
    data_uniran = read1dbinary('uniran.out', n) 
    print('uniran:') 
    print('{0:.15f}'.format(max(data_uniran))) 
    plt.hist(data_uniran, bins=1000) 
    plt.show() 

    data_resamp = read1dbinary('resamp.out', n) 
    print('resample uniran:') 
    print('{0:.15f}'.format(max(data_resamp))) 
    plt.hist(data_resamp, bins=1000) 
    plt.show() 

    data_intrin = read1dbinary('intrin.out', n) 
    print('intrinsic:') 
    print('{0:.15f}'.format(max(data_intrin))) 
    plt.hist(data_intrin, bins=1000) 
    plt.show() 

Все три гистограммы выглядят очень хорошо визуально, но наибольшее значение показывает недостатки uniran:

uniran: 
0.999999880790710 
resample uniran: 
0.999999940395355 
intrinsic: 
0.999999940395355 

Я провел это несколько раз, и результат всегда идентичен. resample_uniran и внутренние значения имеют одинаковую высоту, тогда как uniran также всегда одинаковы, но ниже. Я хотел бы иметь некоторый надежный статистический тест, который показывает, насколько однородным является выход на самом деле, но, пытаясь проверить Андерсон-Дарлинг, тест Койпера и тест Колмогорова-Смирнова, я столкнулся с this problem. По сути, чем больше образцов у вас есть, тем выше вероятность того, что тесты найдут что-то неправильно с выходом. Может быть, нужно что-то сделать, например, this, но пока я еще не добрался до этого.

Для полноты module:

module mod_prngtest 
implicit none 
integer :: iseed_i, iseed_j, iseed_k, iseed_n 
integer, dimension(4) :: seed 

contains 

    function uniran() 
    ! Generate uniformly distributed random numbers in [0, 1) from genrand_int32 
    ! New version 
     integer, parameter :: dp = selected_real_kind(15, 307) 
     real(dp) :: tmp 
     real :: uniran 
     tmp = 0.5_dp + 0.2328306e-9_dp * genrand_int32() 
     uniran = real(tmp) 
    end function uniran 

    recursive function resample_uniran() result(res) 
    ! Generate uniformly distributed random numbers in [0, 1) from genrand_int32 
    ! New version, now recursive 
     integer, parameter :: dp = selected_real_kind(15, 307) 
     real(dp) :: tmp 
     real :: res 
     tmp = 0.5_dp + 0.23283064365386962890625e-9_dp * genrand_int32() 
     res = real(tmp) 
     if (res == 1.0) then 
      res = resample_uniran() 
     end if 
    end function resample_uniran 

    recursive subroutine both(uniran, resamp) 
     integer, parameter :: dp = selected_real_kind(15, 307) 
     real(dp) :: tmp1, tmp2 
     integer :: prn 
     real :: uniran, resamp 

     prn = genrand_int32() 

     tmp1 = 0.5_dp + 0.2328306e-9_dp * prn 
     uniran = real(tmp1) 

     tmp2 = 0.5_dp + 0.23283064365386962890625e-9_dp * prn 
     resamp = real(tmp2) 
     if (resamp == 1.0) then 
      call both(uniran, resamp) 
     end if 
    end subroutine both 

    function genrand_real() 
    ! Generate uniformly distributed random numbers in [0, 1) from genrand_int32 
    ! Your version, modified by me earlier 
     real genrand_real, r 
     r = real(genrand_int32()) 
     if (r .lt. 0.0) r = r + 2.0**32 
     genrand_real = r/4294967296.0 
     return 
    end 

    subroutine init_genrand_int32() 
    ! seed the PRNG, if you don't have /dev/urandom comment out this block ... 
     open(11, file='/dev/urandom', form='unformatted', access='stream') 
     read(11) seed 
     iseed_i=1+abs(seed(1)) 
     iseed_j=1+abs(seed(2)) 
     iseed_k=1+abs(seed(3)) 
     iseed_n=1+abs(seed(4)) 

    ! ... and use this block instead (any integer > 0) 
     !iseed_i = 1253795357 
     !iseed_j = 520466003 
     !iseed_k = 68202083 
     !iseed_n = 1964789093 
    end subroutine init_genrand_int32 

    function genrand_int32() 
    ! From Marsaglia 1994, return pseudorandom integer over the 
    ! whole range. Fortran doesn't have a function like that intrinsically. 
    ! Replace this with your Mersegne twister PRNG 
     implicit none 
     integer :: genrand_int32 
     genrand_int32=iseed_i-iseed_k 
     if(genrand_int32.lt.0)genrand_int32=genrand_int32+2147483579 
     iseed_i=iseed_j 
     iseed_j=iseed_k 
     iseed_k=genrand_int32 
     iseed_n=69069*iseed_n+1013904243 
     genrand_int32=genrand_int32+iseed_n 
    end function genrand_int32 

    subroutine init_random_seed() 
     use iso_fortran_env, only: int64 
     implicit none 
     integer, allocatable :: seed(:) 
     integer :: i, n, un, istat, dt(8), pid 
     integer(int64) :: t 

     call random_seed(size = n) 
     allocate(seed(n)) 
     ! First try if the OS provides a random number generator 
     open(newunit=un, file="/dev/urandom", access="stream", & 
      form="unformatted", action="read", status="old", iostat=istat) 
     if (istat == 0) then 
      read(un) seed 
      close(un) 
     else 
      ! Fallback to XOR:ing the current time and pid. The PID is 
      ! useful in case one launches multiple instances of the same 
      ! program in parallel. 
      call system_clock(t) 
      if (t == 0) then 
       call date_and_time(values=dt) 
       t = (dt(1) - 1970) * 365_int64 * 24 * 60 * 60 * 1000 & 
        + dt(2) * 31_int64 * 24 * 60 * 60 * 1000 & 
        + dt(3) * 24_int64 * 60 * 60 * 1000 & 
        + dt(5) * 60 * 60 * 1000 & 
        + dt(6) * 60 * 1000 + dt(7) * 1000 & 
        + dt(8) 
      end if 
      pid = getpid() 
      t = ieor(t, int(pid, kind(t))) 
      do i = 1, n 
       seed(i) = lcg(t) 
      end do 
     end if 
     call random_seed(put=seed) 
    contains 
     ! This simple PRNG might not be good enough for real work, but is 
     ! sufficient for seeding a better PRNG. 
     function lcg(s) 
      integer :: lcg 
      integer(int64) :: s 
      if (s == 0) then 
       s = 104729 
      else 
       s = mod(s, 4294967296_int64) 
      end if 
      s = mod(s * 279470273_int64, 4294967291_int64) 
      lcg = int(mod(s, int(huge(0), int64)), kind(0)) 
     end function lcg 
     end subroutine init_random_seed 
end module mod_prngtest 

Answer 2

Я не знаю, Fortran на всех, но попробовать что-то вроде этого:

function genrand_real() 
    real genrand_real, r 
    integer genrand_int32 
    r = real(IAND(genrand_int32(), 16777215)) 
    genrand_real = r/16777216.0 
    return 
end 

Я рискую искажая тонкости плавающего Точечное округление на языке, который я не знаю, но я все равно попробую ...

Ваша проблема в том, что вы пытаетесь сжать слишком много бит в мантиссу 32-битное значение с плавающей запятой. Это вызывает проблемы округления, которые могут подтолкнуть значение, близкое к 1.0, равное 1.0. В то же время это может привести к округлению округленных значений от 0.0, и потому что нет ничего ниже 0, чтобы округлить до 0, это оставляет вам меньше шансов получить 0.0.

Если вы попытаетесь устранить проблему, используя 32 бита и настроив масштабный коэффициент, чтобы безопасно довести его до уровня ниже 1.0, вы все равно столкнетесь с проблемой неравномерного распределения. Но если вы исправите диапазон в целочисленном пространстве, используя только столько битов, сколько вы можете точно представлять (24 бита для 32-битного поплавка), вам не нужно беспокоиться о том, что значения округляются вверх или вниз несбалансированным образом ,