准蒙特卡洛法模拟器 — Sobol序列 返回
数值分析

准蒙特卡洛法模拟器 — Sobol序列

使用Sobol序列、Halton序列等"低偏差序列"进行数值积分的工具。改变样本点数、序列类型和被积函数时,决定性点列可均匀填充单位正方形,误差比伪随机数收敛得更快,可实时观察这一过程。

参数设置
样本点数 N
在单位正方形中放置的点数
序列类型
点如何排列
被积函数
在单位正方形[0,1]²上积分的函数
计算结果
积分估计值
真实值
绝对误差
序列类型
相对伪随机的精度提升
样本点数
单位正方形上的点配置

单位正方形[0,1]²中的采样点。伪随机数会产生聚集和空隙,而Sobol序列和Halton序列均匀分散。四分圆情况下,圆内的点用不同颜色区分。

误差收敛(对数-对数)
推定值收敛
理论和主要公式

$$I=\int_{[0,1]^2}f(\mathbf{x})\,d\mathbf{x}\;\approx\;\frac{1}{N}\sum_{i=1}^{N}f(\mathbf{x}_i)$$

数值积分用在样本点x_i处的被积函数f的平均值来近似。点序列的选择决定了收敛速度。

$$\text{误差}_{MC}\sim\frac{1}{\sqrt N}\qquad\text{vs}\qquad \text{误差}_{QMC}\sim\frac{(\log N)^d}{N}$$

伪随机蒙特卡洛法仅以1/√N收敛,但低偏差序列均匀填充区域,准蒙特卡洛法误差收敛速度快得多(d为维数)。

$$\left|I-\frac{1}{N}\sum f(\mathbf{x}_i)\right|\le V(f)\,D_N^*$$

Koksma–Hlawka不等式。积分误差由函数的总变分V(f)与点序列的星形偏差D*_N的乘积从上界限制。

准蒙特卡洛法是什么

🙋
\"准蒙特卡洛法\"和普通蒙特卡洛法有什么区别?名称相似让我有点困惑。
🎓
问得很好。普通蒙特卡洛积分像掷骰子一样用\"随机数\"散布点,然后取那些点处函数值的平均来估算积分。很方便,但随机数其实有偏差。把左边的序列切换为\"伪随机数\",看画布就会发现,点在某些地方聚集得像丸子一样,而其他地方留下空隙。准蒙特卡洛法就是把这些随机数换成故意\"均匀散开的决定性点序列\"。
🙋
看起来随机数应该\"更公平\"地散开,却反而聚集了。当我切换到Sobol序列时,点确实像格子一样整齐排列。
🎓
正是这里最违反直觉的地方。随机数\"下一个点打在哪\"和过去的点无关,所以无法避免碰巧打在附近。因此会聚集。Sobol序列和Halton序列叫\"低偏差序列\",会故意避开已打过的点,在各个尺度上都没有缝隙地填满。积分的关键是\"用多少均匀的代表点采样区域\",所以这种均匀性直接变成精度。
🙋
\"偏差量\"这个词出现了,它衡量的是什么?
🎓
简单说就是\"点有多偏斜\"。在区域里随便取一个大小的框,框内的点比例和框的面积完全一致最理想。偏差量就是对所有框衡量这个偏离的最坏情况。伪随机数因为有聚集就偏差大,Sobol序列被故意设计得很小。Koksma–Hlawka不等式告诉我们\"积分误差≤函数变化×偏差量\"。所以减小偏差量就能减小误差上界。
🙋
下面\"误差收敛\"图中,Sobol序列的线下降得更陡。这就是效果吗?
🎓
完全正确。在两对数图上,斜率表示收敛速度。伪随机数的误差是1/√N,点数翻倍误差只少一半。准蒙特卡洛接近(log N)ᵈ/N,基本上就是1/N,点数翻倍误差就少一半。同样精度的话需要样本点数差别巨大。这就是\"计算金融的期权定价和渲染中Sobol序列是标准\"的原因。
🙋
那是不是就该永远用Sobol序列,没有缺点?
🎓
不是万能的。误差上界里有(log N)ᵈ这一项,当d特别大时,对数的幂就会起作用,优势会衰退。还有如果被积函数不连续或总变分是无穷大,Koksma–Hlawka就保证不了。实务中常见的做法是\"以Sobol序列为基础,加入微量随机位移来做乱置\",既保留了决定性的好处,又便于误差评估。先从这个工具开始,切换序列比较收敛差异,体会一下吧。

常见问题

普通蒙特卡洛积分使用伪随机数散布点,导致点在某些地方聚集,其他地方留下空隙。误差仅以1/√N的速率减小。准蒙特卡洛法用Sobol序列或Halton序列等"低偏差序列"(low-discrepancy sequence)代替随机数。这些序列在各个尺度上都被设计为均匀填充区域,误差以接近(log N)ᵈ/N的速率减小。换言之,用相同的精度需要远少得多的样本数。
偏差量(discrepancy)是衡量点序列"均匀"填充空间程度的指标。在区域内任取任意小的框,框内点数占总点数的比例与框体积(面积)的接近程度,对所有框求最坏情况。伪随机数由于聚集和空隙,偏差量较大;Sobol序列和Halton序列故意设计成偏差量很小。根据Koksma–Hlawka不等式,积分误差受被积函数的变分与点序列偏差量乘积的上界。
Sobol序列对每一维预先准备"方向数"(direction numbers)整数序列,用格雷码表示点的索引,按位进行XOR累积生成。第一维的方向数与基数2的基本反演(van der Corput序列)一致,第二维及以后由原始多项式确定的方向数推导。仅用整数运算的XOR就能快速计算,且逐点添加下一点也不会破坏平衡——这称为"网络性"。本工具使用标准格雷码构造生成二维Sobol序列。
在需要高精度、低成本进行高维积分的场景广泛应用。代表例子是计算金融中的期权定价评估(数十至数百维路径积分)、物理基础渲染(光线路径采样)和工程不确定性量化(多个输入变化对响应统计)。这些场景中,纯随机数收敛缓慢计算耗时,但用Sobol序列替代可用远少的样本获得相同精度。不过次元非常高或被积函数变动很大时,优势会减弱。

实世界应用

计算金融(期权定价评估):亚式期权和篮子期权的价格表示为数十至数百个时间步长的概率路径期望值,即高维积分。伪随机蒙特卡洛法收敛缓慢(1/√N),精度提升一位数需要100倍的路径。用Sobol序列替换能用少得多的路径获得相同精度,大幅缩短衍生品评估时间。这是金融界广泛采用准蒙特卡洛法的代表原因。

物理基础渲染(CG):电影和游戏渲染中,每个像素的颜色计算为\"沿光线路径的积分\"。用伪随机采样会留下明显的颗粒噪声,去除需要海量样本。用Sobol序列等低偏差序列,同样样本数下噪声明显减少,缩短渲染时间。许多光线追踪器内部采用Sobol序列。

工程不确定性量化(UQ):结构分析和流体分析中,材料属性、尺寸公差、载荷等输入变化如何传递到结果。这是对多个不确定输入的响应期望值和方差的高维积分。用Sobol序列采样,用少的有限元试算就能精确推算置信区间和敏感指标(Sobol灵敏度指标的名字就源于此),节省昂贵模拟的计算成本。

全局优化和试验设计:设计变量众多的优化中,初始探索点需均匀散布在空间。网格当维数增大会爆炸,随机数有偏差。Sobol序列和Halton序列在任意维数和任意点数下提供均匀覆盖,广泛用于机器学习超参数搜索和代理模型学习点配置(空间填充设计)的初始采样。

常见误解和注意事项

最常见误解是\"准蒙特卡洛法在任何维数都一直很快\"。误差上界(log N)ᵈ/N含有维数d的幂次。超过数百维时,除非N非常大,否则(log N)ᵈ项占主导,看起来和1/√N差不多。现代理解是:如果被积函数\"实际上低维\"(少数变量有效),准蒙特卡洛在高维也有效。维数本身不是关键,函数结构才是。

其次\"Sobol序列决定性所以完全没法评估误差\"这个想法是错的。固定点序列确实每次结果一样,不能像随机数那样算标准差。实务中用\"随机化准蒙特卡洛法\":给整个Sobol序列加均匀随机位移,或对比特做乱置,保持均匀覆盖的好处同时注入随机性。多次独立运行就能用散布估计误差。本工具展示纯序列行为,所以用了无位移的原始点序列。

最后\"序列既然均匀就不用关心顺序\"是误解。Sobol序列和Halton序列在\"2的幂\"点数(Sobol)或基数幂组合(Halton)这种"舒适"点数特别均匀。点数中途停下来最后几个点可能有偏差。而且低质实现的高维Halton列有臭名昭著的相关性(条纹图案)。实用中要用验证过的库实现,点数也选舒适值,这样才安全。

使用指南

  1. 在256~65536范围内设置nNum(样本点数),生成Sobol序列的准蒙特卡洛采样点
  2. 在0~π区间指定nRange(积分范围),执行sin(x)、cos(x)等被积函数的数值积分
  3. 点击\"开始计算\"按钮,同时模拟Sobol序列和伪随机数两种方法,对比观察收敛曲线

具体计算例

在[0, π]区间积分sin(x),解析解为2.0。用1024个点的Sobol序列绝对误差被抑制在0.0018,同样样本数的伪随机数(Mersenne Twister)产生约0.0045的误差。16384个点时Sobol序列误差降至0.00008,比伪随机数的0.0012好2个数量级。偏差以O(log(N)/N)低速收敛,而Sobol序列实现O(1/N)高速收敛。

实务注意