时间戳

首页 > 解决方案 > MT19937 不会通过将种子值保持为常数来重现相同的伪随机序列

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

问题描述

我在 Fortran 90/95 中的 Monte Carlo 模拟中编写了一个检查点函数，我使用的编译器是 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()输出并不相同。

由于这个意想不到的结果，我在系统的任意状态下恢复模拟时遇到了困难，因为模拟没有再现我正在模拟的系统的最新配置状态。

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

标签： randomfortranmersenne-twister

解决方案

从您提供的源代码（参见 [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

第一个是“状态向量的数组”，第二个是一个标志，确保我们不会从未初始化的数组开始，第三个是一些位置标记，正如我从评论中所述的条件中猜测的那样。

看着subroutine init_genrand( s )，它设置mtinit标志，并从upto填充mt()数组。好吧。1n

看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 生成器将从不同的起点读取。我怀疑它是故意的，所以它可能是一个容易被忽视的小错误。

推荐阅读