如何初始化PRNG使Julia与Python生成相同数据集？

问题原因 #

Julia的Random模块与NumPy默认的随机数生成器实现不同（Julia用MersenneTwister，NumPy默认用PCG64），即使设置相同种子，生成的随机数序列也不一致，导致最终u_meas结果不同。

解决方案1：硬编码Julia生成的随机数（结果完全一致） #

直接使用Julia代码生成的随机数值，在Python中复现计算过程：

import numpy as np

def f(x, t):
    return 3.2 * (x + 0.2 * t)

n_meas = 20
# Julia生成的X第一列（对应x）
x = np.array([0.5331830160153941, 0.4540291384680807, 0.01768684265536442,
              0.16008711641342648, 0.2723255096534325, 0.2986142831666028,
              0.23881082143275582, 0.8709608367886603, 0.737396059590325,
              0.6523466502079163, 0.2505121861650398, 0.2877759568670005,
              0.4190074290995679, 0.2960415950133171, 0.7112745832391288,
              0.5662251655185582, 0.04269782521750889, 0.12630699927892007,
              0.6005517413497284, 0.1264084457427827])
# Julia生成的X第二列（对应t，是2*rand(n_meas)的结果）
t = np.array([1.2276999094664066, 0.35717915096776506, 1.8667119700408316,
              1.088767709338047, 1.150521635131777, 0.6502591841527363,
              0.04067050120077841, 1.9608977412308806, 0.7975869978170758,
              1.4796665140226275, 1.9767697119308938, 1.0354908624277027,
              0.4821772050977875, 0.5322919267987108, 1.4005138566262444,
              0.8021581088737053, 0.1410205100771572, 0.7432092960456113,
              0.2987941936785101, 1.564067929519556])

# 计算基础值
u_meas = f(x, t)
# Julia生成的标准正态噪声（randn(20)结果）
noise = np.array([0.364645607605806, -0.601706612818598, 0.428660341168872,
                  -0.728435726781429, 0.400420242259798, -0.089902896165849,
                  -0.722160652320023, 0.237600666313441, 1.94392293868275,
                  -0.130704628331112, -0.124135348080726, 0.773250830640705,
                  0.675385666202709, -0.303724282114892, -0.294794760301971,
                  -1.34272227025067, 0.98877389352771, 0.287461605383667,
                  -0.575742840714578, 0.410984335301736])
# 添加噪声
u_meas += 0.05 * np.mean(u_meas) * noise

print(u_meas)

解决方案2：使用Python的MersenneTwister匹配Julia随机数生成 #

Julia的Random.seed!(42)使用MersenneTwister 19937生成器，Python的numpy.random可以切换到相同生成器，生成结果会非常接近Julia代码（因实现细节差异可能存在微小浮点误差）：

import numpy as np

# 切换到MersenneTwister生成器，匹配Julia的随机数生成器
rng = np.random.MT19937(seed=42)
rng = np.random.Generator(rng)

def f(x, t):
    return 3.2 * (x + 0.2 * t)

n_meas = 20
X = np.zeros((n_meas, 2))
# Julia的rand()生成[0,1)均匀分布，与numpy的random一致
X[:, 0] = rng.random(n_meas)
X[:, 1] = 2 * rng.random(n_meas)

u_meas = f(X[:, 0], X[:, 1])
# Julia的randn()生成标准正态分布，对应numpy的standard_normal
u_meas += 0.05 * np.mean(u_meas) * rng.standard_normal(n_meas)

print(u_meas)

内容的提问来源于stack exchange，提问作者NGA