MT19937不通过保持种子值恒定来产生相同的伪随机序列。

Question

MT19937不通过保持种子值恒定来产生相同的伪随机序列。

5

我正在使用ifort 18.0.2编译器，在Fortran 90/95中编写Monte Carlo模拟程序的checkpoint功能。在详细说明之前，我想澄清一下我正在使用的伪随机生成器版本：

A C-program for MT19937, with initialization, improved 2002/1/26.
Coded by Takuji Nishimura and Makoto Matsumoto.

Code converted to Fortran 95 by Josi Rui Faustino de Sousa
Date: 2002-02-01

请参考mt19937获取源代码。

我的蒙特卡罗模拟代码的一般结构如下：

program montecarlo
 call read_iseed(...)
 call mc_subroutine(...)
end

在read_iseed函数中

subroutine read_iseed(...)
  use mt19937

    if (Restart == 'n') then

    call system('od -vAn -N4 -td4 < /dev/urandom > '//trim(IN_ISEED)
    open(unit=7,file=trim(IN_ISEED),status='old')
    read(7,*) i
    close(7)

    !This is only used to initialise the PRNG sequence
    iseed = abs(i)
    else if (Restart == 'y') then

    !Taking seed value from the latest iteration of previous simulation
    iseed = RestartSeed

    endif

    call init_genrand(iseed)
    print *, 'first pseudo-random value ',genrand_real3(), 'iseed ',iseed

    return
end subroutine

根据我的理解，如果种子值保持不变，伪随机数生成器应该能够每次重现相同的伪随机序列？

为了证明这一点，我使用相同的种子值运行了两个独立的模拟，它们能够重现完全相同的序列。到目前为止一切都很好！

基于之前的测试，我进一步假设无论在一个单独的模拟中调用多少次init_genrand()，PRNG也应该能够重现伪随机值序列？所以我对我的read_iseed()子程序进行了一些修改。

subroutine read_iseed(...)
  use mt19937

    if (Restart == 'n') then

    call system('od -vAn -N4 -td4 < /dev/urandom > '//trim(IN_ISEED)
    open(unit=7,file=trim(IN_ISEED),status='old')
    read(7,*) i
    close(7)

    !This is only used to initialise the PRNG sequence
    iseed = abs(i)
    else if (Restart == 'y') then

    !Taking seed value from the latest iteration of the previous simulation
    iseed = RestartSeed

    endif

    call init_genrand(iseed)
    print *, 'first time initialisation ',genrand_real3(), 'iseed ',iseed

    call init_genrand(iseed)
    print *, 'second time initialisation ',genrand_real3(), 'iseed ',iseed

    return
end subroutine

输出结果与我的预期不同，尽管在两次初始化之间iseed的输出是相同的，但genrand_real3()的输出却不同。由于这个意外的结果，我很难在系统的任意状态下恢复模拟，因为模拟无法再现我正在模拟的系统的最新配置状态。

我不确定是否提供了足够的信息，请让我知道如果这个问题的任何部分需要更具体的说明。

- Gvxfjørt

请在所有Fortran问题中使用标签[tagnfortran]。没有必要为Fortran90和Fortran95分别打标签，因为您的问题无论如何都不是版本特定的。 - Vladimir F Героям слава

init_genrand 如何填充完整的 PRNG 状态？ - francescalus

@francescalus 请查看[mt19937ar]{http://web.mst.edu/~vojtat/class_5403/mt19937/mt19937ar.f90}，我的第一次测试证明它能够通过提供一个恒定的种子值来重现序列，但在我的第二次测试中，当我使用相同的种子值重新初始化序列时，它似乎没有将序列状态重置回到起始位置。根据我阅读的子程序，函数`init_genrand`似乎是初始化序列开始的位置，给定一个种子值，而`genrand_real3()`提供(0,1)之间的值。您是否发现任何不同意见？ - Gvxfjørt

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- quetzalcoatl · Accepted Answer

从您提供的源代码中（请参见[mt19937]{http://web.mst.edu/~vojtat/class_5403/mt19937/mt19937ar.f90}获取源代码），init_genrand 不会清除 整个状态。

有3个关键状态变量：

integer( kind = wi )  :: mt(n)            ! the array for the state vector
logical( kind = wi )  :: mtinit = .false._wi   ! means mt[N] is not initialized
integer( kind = wi )  :: mti = n + 1_wi   ! mti==N+1 means mt[N] is not initialized

第一个是“状态向量数组”，第二个是确保我们不用未初始化的数组开始的标志，第三个是一些位置标记，我猜这是从注释中所述的条件。

看看子程序init_genrand(s)，它设置了mtinit标志，并从1到n填充了mt()数组。好的。

看看genrand_real3，它基于genrand_int32。

看看genrand_int32，它从头开始：

if ( mti > n ) then ! generate N words at one time
  ! if init_genrand() has not been called, a default initial seed is used
  if ( .not. mtinit ) call init_genrand( seed_d )

然后进行算术计算，然后开始得到结果：

y = mt(mti)
mti = mti + 1_wi

所以.. mti 是“状态数组”中的一个位置索引，每次从生成器读取整数后都会将其增加1。

回到init_genrand - 记得吗？它已经重置了数组mt()，但是它没有重置MTI回到其起始点mti = n + 1_wi。

我敢打赌这就是你观察到的现象的原因，因为在使用相同种子重新初始化后，数组将填充相同的值集合，但稍后int32生成器将从不同的起始点读取。我怀疑这并不是有意的，所以很可能是一个容易忽视的小错误。