returns 更新种子值的随机数生成器 (RNG/PRNG)
Random number generator (RNG/PRNG) that returns updated value of seed
我正在尝试编写一个 RNG，它也 returns 更新种子的值。这样做的一个明显原因可能是新的随机变量可以在以后添加到程序中，而无需更改现有 RNG 的值。有关此问题的 python/numpy 版本，请参阅示例：
这是一个使用（暂定）建议解决方案的示例：
program main
! note that I am assuming size=12 for the random
! seed but this is implementation specific and
! could be different in your implementation
! (you can query this with random_seed(size) btw)
integer :: iseed1(12) = 123456789
integer :: iseed2(12) = 987654321
do i = 1, 2
x = ran(iseed1)
y = ran(iseed2)
end do
end program main
function ran( seed )
integer, intent(inout) :: seed(12)
real :: ran
call random_seed( put=seed )
call random_number( ran )
call random_seed( get=seed )
end function ran
请注意，此解决方案（以及任何其他解决方案）的一个重要方面是 如果我们在上面添加 seed3 和 x3，则不会有任何变化x1 和 x2 的实现值。同样，可以从代码中删除 x1 或 x2，而不会影响另一个的值。
如果有帮助，这就是 ran()（扩展）函数在 Compaq/HP 和英特尔编译器上的工作方式，我基本上是在尝试在 gfortran 上实现相同的 API。但是请注意，该函数的种子是标量，而它是使用 Fortran 90 子例程 random_seed 的一维数组（数组的 size/length 是特定于实现的）。
我提供我当前的解决方案作为基准错误，但希望其他人可以批评该答案或提供更好的答案。
如果能根据标准 PRNG 测试对基准测试结果进行任何分析，尤其是对种子设置方式的分析，我将不胜感激。在基准测试中，我只是使用标量广播到一个非常简单明了的数组，并希望避免显式提供整个整数数组。
所以我想要一个稍微严格的确认这很好，或者一个更简洁的方法来设置种子（以可重复的方式！）而不是写出整个数组，比如 12 或 33 个整数。例如。也许有一些简单明了的方法可以从一个整数中产生一个很好的 12 个伪随机整数流（因为数组长度很短，这可能很容易？）。
编辑添加： 关于在 fortran 中设置种子的后续问题：
你提出的解决方案在我看来应该可行 - 你正在记录生成器的整个状态（通过 get），并在必要时在流之间交换（通过 put）。不过请注意，我还没有实际测试过您的代码。
之所以出现这个答案，是因为之前的一个答案（现已删除）只保存了状态数组的第一个元素，并用它来设置整个状态数组。它看起来像：
integer :: seed_scalar, state_array(12)
...
! -- Generate a random number from a thread that is stored as seed_scalar
state_array(:) = seed_scalar
call random_seed( put=state_array )
call random_number( ran )
call random_seed( get=state_array )
seed_scalar = state_array(1)
! -- discard state_array, store seed_scalar
此答案旨在证明，对于某些编译器（gfortran 4.8 和 8.1，pgfortran 15.10），这种仅通过标量维护单独线程的方法会导致错误的 RNG 行为。
考虑以下代码来初步测试随机数生成器。生成了许多随机数——本例中为 100M——并以两种方式监控性能。首先，记录随机数增加或减少的次数。其次，生成 bin 宽度为 0.01 的直方图。 （这显然是一个原始的测试用例，但事实证明足以证明这一点。）。最后，还生成了所有概率的估计 one-sigma 标准差。这使我们能够确定变化何时是随机的或具有统计意义的。
program main
implicit none
integer, parameter :: wp = selected_real_kind(15,307)
integer :: N_iterations, state_size, N_histogram
integer :: count_incr, count_decr, i, loc, fid
integer, allocatable, dimension(:) :: state1, histogram
real(wp) :: ran, oldran, p, dp, re, rt
! -- Some parameters
N_iterations = 100000000
N_histogram = 100
call random_seed(size = state_size)
allocate(state1(state_size), histogram(N_histogram))
write(*,'(a,i3)') '-- State size is: ', state_size
! -- Initialize RNG, taken as today's date (also tested with other initial seeds)
state1(:) = 20180815
call random_seed(put=state1)
count_incr = 0
count_decr = 0
histogram = 0
ran = 0.5_wp ! -- To yield proprer oldran
! -- Loop to grab a batch of random numbers
do i=1,N_iterations
oldran = ran
! -- This preprocessor block modifies the RNG state in a manner
! -- consistent with storing only the first value of it
#ifdef ALTER_STATE
call random_seed(get=state1)
state1(:) = state1(1)
call random_seed(put=state1)
#endif
! -- Get the random number
call random_number(ran)
! -- Process Histogram
loc = CEILING(ran*N_histogram)
histogram(loc) = histogram(loc) + 1
if (oldran<ran) then
count_incr = count_incr + 1
elseif (oldran>ran) then
count_decr = count_decr + 1
else
! -- This should be statistically impossible
write(*,*) '** Error, same value?', ran, oldran
stop 1
endif
enddo
! -- For this processing, we have:
! re Number of times this event occurred
! rt Total number
! -- Probability of event is just re/rt
! -- One-sigma standard deviation is sqrt( re*(rt-re)/rt**3 )
! -- Write histogram
! -- For each bin, compute probability and uncertainty in that probability
open(newunit=fid, file='histogram.dat', action='write', status='replace')
write(fid,'(a)') '# i, P(i), dP(i)'
rt = real(N_iterations,wp)
do i=1,N_histogram
re = real(histogram(i),wp)
p = re/rt
dp = sqrt(re*(rt-1)/rt**3)
write(fid,'(i4,2es26.18)') i, p, dp
enddo
close(fid)
! -- Probability of increase and decrease
re = real(count_incr,wp)
p = re/rt
dp = sqrt(re*(rt-1)/rt**3)
write(*,'(a,f20.15)') 'Probability of increasing: ', p
write(*,'(a,f20.15)') ' one-sigma deviation: ', dp
re = real(count_decr,wp)
p = re/rt
dp = sqrt(re*(rt-1)/rt**3)
write(*,'(a,f20.15)') 'Probability of decreasing: ', p
write(*,'(a,f20.15)') ' one-sigma deviation: ', dp
write(*,'(a)') 'complete'
end program main
无修改
没有预处理器指令 ALTER_STATE，我们按预期使用 gfortran 的内置 PRNG，结果符合预期：
enet-mach5% gfortran --version
GNU Fortran (SUSE Linux) 4.8.3 20140627 [gcc-4_8-branch revision 212064]
Copyright (C) 2013 Free Software Foundation, Inc.
GNU Fortran comes with NO WARRANTY, to the extent permitted by law.
You may redistribute copies of GNU Fortran
under the terms of the GNU General Public License.
For more information about these matters, see the file named COPYING
enet-mach5% gfortran -cpp -fcheck=all main.f90 && time ./a.out
-- State size is: 12
Probability of increasing: 0.499970710000000
one-sigma deviation: 0.000070708606619
Probability of decreasing: 0.500029290000000
one-sigma deviation: 0.000070712748851
complete
real 0m2.414s
user 0m2.408s
sys 0m0.004s
increasing/decreasing 的预期概率为 0.5，两者都具有估计的不确定性（0.49997 小于 0.00007，从 0.5 开始）。带有误差条的直方图是
对于每个 bin，预期概率 (0.01) 的变化很小，通常在估计的不确定性范围内。因为我们生成了很多数字，所以所有变化都很小（0.1% 量级）。本次测试基本上没有发现任何可疑行为。
有修改
如果我们启用 ALTER_STATE 中的块，我们将在每次生成数字时修改随机数生成器的内部状态。这是为了模仿现在已删除的解决方案，该解决方案仅存储状态的第一个值。结果是：
enet-mach5% gfortran -cpp -DALTER_STATE -fcheck=all main.f90 && time ./a.out
-- State size is: 12
Probability of increasing: 0.501831930000000
one-sigma deviation: 0.000070840096343
Probability of decreasing: 0.498168070000000
one-sigma deviation: 0.000070581021884
complete
real 0m16.489s
user 0m16.492s
sys 0m0.000s
观察到的增加概率 远 超出预期变化（26 西格玛！）。这已经表明出了问题。直方图是：
请注意，y 的范围发生了很大变化。在这里，我们的变化比之前的情况大大约两个数量级，远远超出预期的变化。在这里很难看到误差线，因为 y 范围要大得多。如果我的随机数生成器是这样执行的，我会觉得用它做任何事情都不舒服，即使是抛硬币也不行。
关闭
random_seed 的 put 和 get 选项访问随机数生成器的处理器相关内部状态。这通常比单个数字具有更多的熵，并且其形式取决于处理器。不能保证第一个数字代表整个州。
如果您要初始化一个随机种子一次，然后生成多次，那么使用单个标量就可以了。但是，如果您打算使用该状态生成每个数字。
，则显然有必要存储多个数字。
坦率地说，我有点惊讶这个原始测试能够证明不良行为。 RNG 的有效性是一个复杂的课题，我绝不是专家。结果也依赖于编译器：
- 显示的结果和直方图适用于状态大小为 12 的 gfortran 4.8。
- Intel 16.0 使用状态大小 2，此测试未显示任何不良行为。
- gfortran 8.1.0 的状态大小为 33，PGI 15.10 的状态大小为 34。它们的行为是一样的，而且是最糟糕的。当将整个状态设置为单个值时，后续随机数 总是 相同。当从单个种子初始化时，这些生成器需要 'priming' 大约 30 代才能使数字开始看起来合理。当在状态下只保存一个种子时，这种启动永远不会发生。
当仅保存单个值时，较大的状态大小可能会导致更多的熵损失，因此它与较差的行为相关。这当然符合我的观察。但是，如果不了解每个生成器的内部结构，就无法分辨。
我正在尝试编写一个 RNG，它也 returns 更新种子的值。这样做的一个明显原因可能是新的随机变量可以在以后添加到程序中，而无需更改现有 RNG 的值。有关此问题的 python/numpy 版本，请参阅示例：
这是一个使用（暂定）建议解决方案的示例：
program main
! note that I am assuming size=12 for the random
! seed but this is implementation specific and
! could be different in your implementation
! (you can query this with random_seed(size) btw)
integer :: iseed1(12) = 123456789
integer :: iseed2(12) = 987654321
do i = 1, 2
x = ran(iseed1)
y = ran(iseed2)
end do
end program main
function ran( seed )
integer, intent(inout) :: seed(12)
real :: ran
call random_seed( put=seed )
call random_number( ran )
call random_seed( get=seed )
end function ran
请注意，此解决方案（以及任何其他解决方案）的一个重要方面是 如果我们在上面添加 seed3 和 x3，则不会有任何变化x1 和 x2 的实现值。同样，可以从代码中删除 x1 或 x2，而不会影响另一个的值。
如果有帮助，这就是 ran()（扩展）函数在 Compaq/HP 和英特尔编译器上的工作方式，我基本上是在尝试在 gfortran 上实现相同的 API。但是请注意，该函数的种子是标量，而它是使用 Fortran 90 子例程 random_seed 的一维数组（数组的 size/length 是特定于实现的）。
我提供我当前的解决方案作为基准错误，但希望其他人可以批评该答案或提供更好的答案。
如果能根据标准 PRNG 测试对基准测试结果进行任何分析，尤其是对种子设置方式的分析，我将不胜感激。在基准测试中，我只是使用标量广播到一个非常简单明了的数组，并希望避免显式提供整个整数数组。
所以我想要一个稍微严格的确认这很好，或者一个更简洁的方法来设置种子（以可重复的方式！）而不是写出整个数组，比如 12 或 33 个整数。例如。也许有一些简单明了的方法可以从一个整数中产生一个很好的 12 个伪随机整数流（因为数组长度很短，这可能很容易？）。
编辑添加： 关于在 fortran 中设置种子的后续问题：
你提出的解决方案在我看来应该可行 - 你正在记录生成器的整个状态（通过 get），并在必要时在流之间交换（通过 put）。不过请注意，我还没有实际测试过您的代码。
之所以出现这个答案，是因为之前的一个答案（现已删除）只保存了状态数组的第一个元素，并用它来设置整个状态数组。它看起来像：
integer :: seed_scalar, state_array(12)
...
! -- Generate a random number from a thread that is stored as seed_scalar
state_array(:) = seed_scalar
call random_seed( put=state_array )
call random_number( ran )
call random_seed( get=state_array )
seed_scalar = state_array(1)
! -- discard state_array, store seed_scalar
此答案旨在证明，对于某些编译器（gfortran 4.8 和 8.1，pgfortran 15.10），这种仅通过标量维护单独线程的方法会导致错误的 RNG 行为。
考虑以下代码来初步测试随机数生成器。生成了许多随机数——本例中为 100M——并以两种方式监控性能。首先，记录随机数增加或减少的次数。其次，生成 bin 宽度为 0.01 的直方图。 （这显然是一个原始的测试用例，但事实证明足以证明这一点。）。最后，还生成了所有概率的估计 one-sigma 标准差。这使我们能够确定变化何时是随机的或具有统计意义的。
program main
implicit none
integer, parameter :: wp = selected_real_kind(15,307)
integer :: N_iterations, state_size, N_histogram
integer :: count_incr, count_decr, i, loc, fid
integer, allocatable, dimension(:) :: state1, histogram
real(wp) :: ran, oldran, p, dp, re, rt
! -- Some parameters
N_iterations = 100000000
N_histogram = 100
call random_seed(size = state_size)
allocate(state1(state_size), histogram(N_histogram))
write(*,'(a,i3)') '-- State size is: ', state_size
! -- Initialize RNG, taken as today's date (also tested with other initial seeds)
state1(:) = 20180815
call random_seed(put=state1)
count_incr = 0
count_decr = 0
histogram = 0
ran = 0.5_wp ! -- To yield proprer oldran
! -- Loop to grab a batch of random numbers
do i=1,N_iterations
oldran = ran
! -- This preprocessor block modifies the RNG state in a manner
! -- consistent with storing only the first value of it
#ifdef ALTER_STATE
call random_seed(get=state1)
state1(:) = state1(1)
call random_seed(put=state1)
#endif
! -- Get the random number
call random_number(ran)
! -- Process Histogram
loc = CEILING(ran*N_histogram)
histogram(loc) = histogram(loc) + 1
if (oldran<ran) then
count_incr = count_incr + 1
elseif (oldran>ran) then
count_decr = count_decr + 1
else
! -- This should be statistically impossible
write(*,*) '** Error, same value?', ran, oldran
stop 1
endif
enddo
! -- For this processing, we have:
! re Number of times this event occurred
! rt Total number
! -- Probability of event is just re/rt
! -- One-sigma standard deviation is sqrt( re*(rt-re)/rt**3 )
! -- Write histogram
! -- For each bin, compute probability and uncertainty in that probability
open(newunit=fid, file='histogram.dat', action='write', status='replace')
write(fid,'(a)') '# i, P(i), dP(i)'
rt = real(N_iterations,wp)
do i=1,N_histogram
re = real(histogram(i),wp)
p = re/rt
dp = sqrt(re*(rt-1)/rt**3)
write(fid,'(i4,2es26.18)') i, p, dp
enddo
close(fid)
! -- Probability of increase and decrease
re = real(count_incr,wp)
p = re/rt
dp = sqrt(re*(rt-1)/rt**3)
write(*,'(a,f20.15)') 'Probability of increasing: ', p
write(*,'(a,f20.15)') ' one-sigma deviation: ', dp
re = real(count_decr,wp)
p = re/rt
dp = sqrt(re*(rt-1)/rt**3)
write(*,'(a,f20.15)') 'Probability of decreasing: ', p
write(*,'(a,f20.15)') ' one-sigma deviation: ', dp
write(*,'(a)') 'complete'
end program main
无修改
没有预处理器指令 ALTER_STATE，我们按预期使用 gfortran 的内置 PRNG，结果符合预期：
enet-mach5% gfortran --version
GNU Fortran (SUSE Linux) 4.8.3 20140627 [gcc-4_8-branch revision 212064]
Copyright (C) 2013 Free Software Foundation, Inc.
GNU Fortran comes with NO WARRANTY, to the extent permitted by law.
You may redistribute copies of GNU Fortran
under the terms of the GNU General Public License.
For more information about these matters, see the file named COPYING
enet-mach5% gfortran -cpp -fcheck=all main.f90 && time ./a.out
-- State size is: 12
Probability of increasing: 0.499970710000000
one-sigma deviation: 0.000070708606619
Probability of decreasing: 0.500029290000000
one-sigma deviation: 0.000070712748851
complete
real 0m2.414s
user 0m2.408s
sys 0m0.004s
increasing/decreasing 的预期概率为 0.5，两者都具有估计的不确定性（0.49997 小于 0.00007，从 0.5 开始）。带有误差条的直方图是
对于每个 bin，预期概率 (0.01) 的变化很小，通常在估计的不确定性范围内。因为我们生成了很多数字，所以所有变化都很小（0.1% 量级）。本次测试基本上没有发现任何可疑行为。
有修改
如果我们启用 ALTER_STATE 中的块，我们将在每次生成数字时修改随机数生成器的内部状态。这是为了模仿现在已删除的解决方案，该解决方案仅存储状态的第一个值。结果是：
enet-mach5% gfortran -cpp -DALTER_STATE -fcheck=all main.f90 && time ./a.out
-- State size is: 12
Probability of increasing: 0.501831930000000
one-sigma deviation: 0.000070840096343
Probability of decreasing: 0.498168070000000
one-sigma deviation: 0.000070581021884
complete
real 0m16.489s
user 0m16.492s
sys 0m0.000s
观察到的增加概率 远 超出预期变化（26 西格玛！）。这已经表明出了问题。直方图是：
请注意，y 的范围发生了很大变化。在这里，我们的变化比之前的情况大大约两个数量级，远远超出预期的变化。在这里很难看到误差线，因为 y 范围要大得多。如果我的随机数生成器是这样执行的，我会觉得用它做任何事情都不舒服，即使是抛硬币也不行。
关闭
random_seed 的 put 和 get 选项访问随机数生成器的处理器相关内部状态。这通常比单个数字具有更多的熵，并且其形式取决于处理器。不能保证第一个数字代表整个州。
如果您要初始化一个随机种子一次，然后生成多次，那么使用单个标量就可以了。但是，如果您打算使用该状态生成每个数字。
，则显然有必要存储多个数字。坦率地说，我有点惊讶这个原始测试能够证明不良行为。 RNG 的有效性是一个复杂的课题，我绝不是专家。结果也依赖于编译器：
- 显示的结果和直方图适用于状态大小为 12 的 gfortran 4.8。
- Intel 16.0 使用状态大小 2，此测试未显示任何不良行为。
- gfortran 8.1.0 的状态大小为 33，PGI 15.10 的状态大小为 34。它们的行为是一样的，而且是最糟糕的。当将整个状态设置为单个值时，后续随机数 总是 相同。当从单个种子初始化时，这些生成器需要 'priming' 大约 30 代才能使数字开始看起来合理。当在状态下只保存一个种子时，这种启动永远不会发生。
当仅保存单个值时，较大的状态大小可能会导致更多的熵损失，因此它与较差的行为相关。这当然符合我的观察。但是，如果不了解每个生成器的内部结构，就无法分辨。