MT19937不会通过保持种子值为常数来重现相同的伪随机序列

2023年4月8日 375次阅读

我在Fortran 90/95的蒙特卡罗模拟中编写了一个检查点函数,我正在使用的编译器是ifort 18.0.2,在详细介绍我正在使用的伪随机生成器的版本之前：

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

根据我的理解,如果种子值保持不变,PRNG应该能够每次都重现伪随机序列吗？

为了证明这种情况,我使用相同的种子值运行了两次单独的模拟,它们能够重现精确的序列.到现在为止还挺好！

基于之前的测试,我进一步假设无论在一个单独的模拟中调用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()输出不相同.

由于这种意想不到的结果,我努力在系统的任意状态下恢复模拟,因为模拟没有再现我正在模拟的系统的最新配置状态.

我不确定我是否提供了足够的信息,请告诉我这个问题的任何部分是否需要更具体？

最佳答案 从您提供的源代码(参见[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标志,并将mt()数组从1加到n.好的.

看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是’state array’中的位置索引,并且在从生成器读取的每个整数之后递增1.

回到init_genrand – 还记得吗？它一直在重置数组mt(),但它没有将MTI重置回其起始位置mti = n 1_wi.

我敢打赌这是你观察到的现象的原因,因为在用相同的种子重新初始化之后,数组将被填充相同的值集,但是后来int32生成器将从不同的起始点读取.我怀疑它是有意的,所以它可能是一个容易被忽视的小虫子.

点赞