如何使用scipy.stats.qmc实现多次随机化拟蒙特卡洛采样

问题原因说明 #

你遇到的报错是因为单个Sobol实例是有状态的，每次调用random_base2会沿着序列顺序向后生成新的点，累计生成点数不是2的幂时就会触发该错误。

最优实现方案（复用基础序列，仅执行加扰） #

完全符合你「避免每次采样重新生成Sobol序列、仅执行加扰」需求的实现如下，只需初始化一次实例，复用基础序列参数，每次仅重置状态并重做加扰即可：

import numpy as np
from scipy.stats import qmc

# 配置参数
dim = 2
m = 10  # 每次生成2^10=1024个点
n_randomizations = 1000  # 随机化次数
rng = np.random.default_rng(seed=42)  # 固定种子可复现结果

# 仅初始化一次Sobol实例，基础序列参数全局复用
sobol_engine = qmc.Sobol(d=dim, scramble=True, seed=rng)

# 存储所有随机化结果
scrambled_results = np.empty((n_randomizations, 2**m, dim))

for i in range(n_randomizations):
    # 重新生成加扰参数
    sobol_engine.scramble()
    # 重置序列生成指针，回到起始位置
    sobol_engine.reset()
    # 生成当前加扰版本的序列
    scrambled_results[i] = sobol_engine.random_base2(m=m)

替代方案：手动随机移位（更轻量） #

如果追求极致性能，也可以预先生成未加扰的基础Sobol序列，通过随机移位模1的标准随机化QMC方法生成多次实现，该方法完全保留Sobol序列的低偏差性质：

# 预先生成固定的未加扰基础序列
sobol_engine = qmc.Sobol(d=dim, scramble=False)
base_sobol = sobol_engine.random_base2(m=m)

scrambled_results = np.empty((n_randomizations, 2**m, dim))
for i in range(n_randomizations):
    # 生成随机移位向量
    shift = rng.uniform(size=dim)
    # 基础序列加移位后模1得到随机化序列
    scrambled_results[i] = (base_sobol + shift) % 1

低开销替代方案：每次新建实例 #

Sobol实例的初始化开销极低，就算创建1000次也不会有明显性能损耗，也可以通过每次创建独立实例实现：

scrambled_results = []
for _ in range(n_randomizations):
    engine = qmc.Sobol(d=dim, scramble=True, seed=rng.integers(0, 10**6))
    scrambled_results.append(engine.random_base2(m=m))
scrambled_results = np.array(scrambled_results)

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