は、乱数生成器の上限

Question

これは実際に前のもののフォローアップの質問です：私は、前の質問の答えと私の問題の解決策だと思ったものを後 Rounding of double precision to single precision: Forcing an upper bound

、私は再び私のプログラムを実行してみました同じ問題があることがわかった。

私が使用しているMersenne Twisterの実装では、符号付き32ビットのランダムな整数が生成されます。

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 

そして、それはそう、私は次のように使用前の質問で提案を以下、完璧に動作：RNGを実装 男は[0,1）の範囲内のランダムな倍精度浮動小数点数を生成するには、この機能を作りましたランダム単精度浮動小数点数を生成する機能は、範囲内で私は[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の代わりに2.0 ** 32 + 1で数を分けることを考えましたが、理論的には正しいと思います。

Answer 1

私はこの質問を古い質問に投稿するのかこの質問を投稿するのかは不明です。いずれにしても、私は解決策を持つことができます（2番目のコードブロックにあります）。

私は約2年前から同じタスクに使用しているルーチンはこれです：

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になります。あなたはそれを試すことができます：もう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）が多い（私たちはすべての数百万個の数字を話している）、他の2つは、これまで持っていながら、決して失敗しなかった。 再帰バージョンは時間がかかりますが、可能な限り高い数が1に近いため、技術的に優れています。

また、スピードと "均一性"をテストし、[012]のサブルーチン[0,1]でも一様乱数を与える固有のrandom_numberサブルーチンと比較しました。 （慎重に、これは3×512 MBのファイルを作成します）

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 

結果はタイミングがdifferntているにもかかわらず、原則的には常に同じです：

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に依存するが、メルセンヌ・ツイスターは本来のものよりも遅くなる）。

私はまた、各メソッドは、（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() 

すべての3つのヒストグラムは、視覚的に非常に良いように見えますが、最高値はuniranの欠点を明らかに：

uniran: 
0.999999880790710 
resample uniran: 
0.999999940395355 
intrinsic: 
0.999999940395355 

私はこれを数回走らせました。結果は常に同じです。 resample_uniranであり、内在値は同じ最高値を持ちますが、uniranも常に同じですが低い値です。 Anderson-Darlingテストを試している間、KuiperのテストとKolmogorov-Smirnovのテストは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より小さくなる可能性があります。

32ビットを使用して問題を解決しようとすると、スケール係数を1.0以下に安全に引き上げるために微調整を行うと、依然として不均一な分布が問題になります。しかし、正確に表すことができるビット数（32ビット浮動小数点数の場合は24ビット）だけを使用して整数空間の範囲を修正すると、アンバランスな方法で値を切り上げたり縮小したりする心配はありません。