Python SciPy qmc.Sobol用法及代碼示例

本文簡要介紹 python 語言中 scipy.stats.qmc.Sobol 的用法。

用法:

class  scipy.stats.qmc.Sobol(d, *, scramble=True, bits=None, seed=None, optimization=None)#

用於生成(加擾)Sobol 序列的引擎。

Sobol 的序列是 low-discrepancy、quasi-random 數字。可以使用兩種方法繪製點：

• random_base2 ：安全地繪製\(n=2^m\) 點。這種方法保證了序列的平衡特性。

• random ：從序列中繪製任意數量的點。請參閱下麵的警告。

參數 ：：

d： int

序列的維度。最大維數為 21201。

scramble： 布爾型，可選

如果為 True，則使用 LMS+shift 加擾。否則，不進行加擾。默認為真。

bits： 整數，可選

生成器的位數。控製可以生成的最大點數，即 2**bits 。最大值為 64。它與返回類型不對應，返回類型始終為np.float64，以防止點重複。默認值為 None，為了向後兼容，對應於 30。

optimization： {無，“random-cd”，“lloyd”}，可選

是否使用優化方案來提高采樣後的質量。請注意，這是一個後處理步驟，不能保證樣品的所有屬性都得到保留。默認為“無”。

• random-cd：坐標的隨機排列以降低中心差異。基於中心差異的最佳樣本不斷更新。與使用其他差異度量相比，基於中心差異的采樣對 2D 和 3D 子投影顯示出更好的 space-filling 魯棒性。

• lloyd：使用修改後的Lloyd-Max 算法擾動樣本。該過程收斂到等間隔的樣本。

seed： {無，整數， numpy.random.Generator }，可選

如果種子是 int 或 None，一個新的numpy.random.Generator是使用創建的np.random.default_rng(seed).如果種子已經是一個Generator實例，然後使用提供的實例。

注意：

Sobol 序列 [1] 在 \([0,1)^{d}\) 中提供了 \(n=2^m\) 低差異點。對它們進行置亂[3]使它們適用於奇異被積函數，提供了一種誤差估計方法，並且可以提高它們的收斂速度。所實施的加擾策略是(左)線性矩陣加擾(LMS)，然後是數字隨機移位(LMS+移位)[2]。

索博爾序列有許多版本，具體取決於它們的“方向號”。該代碼使用[4]中的方向號。因此，最大維度數為 21201。方向數已使用搜索標準 6 預先計算，並且可以在 https://web.maths.unsw.edu.au/~fkuo/sobol/ 處檢索。

警告

Sobol 序列是一個求積規則，如果使用的樣本量不是 2 的冪，或者跳過第一個點，或者對序列進行細化 [5]，它們就會失去平衡特性。

如果 \(n=2^m\) 點數不夠，則應將 \(2^M\) 點數用於 \(M>m\) 。加擾時，獨立重複的數量 R 不必是 2 的冪。

Sobol' 序列被生成到一定數量\(B\)位。後\(2^B\)點已經生成，序列將重複。因此，出現了錯誤。位數可以通過參數控製位.

參考：

[1]

I. M. Sobol，“立方體中點的分布和積分的準確評估。” ZH。維奇斯爾。填充。我在。物理學，7：784-802，1967。

[2]

J. Matousek，“關於錨定框的 L2 差異”。 J. of Complexity 14, 527-556, 1998。

[3]

Art B. Owen，“加擾 Sobol 和Niederreiter-Xing 點”。複雜性雜誌，14(4):466-489，1998 年 12 月。

[4]

S. Joe 和 F. Y. Kuo，“用更好的二維投影構建 sobol 序列”。 SIAM 科學計算雜誌，30(5):2635-2654，2008。

[5]

Art B. Owen，“關於放棄第一個 Sobol 點。” arXiv:2008.08051，2020 年。

例子：

從 Sobol 的低差異序列中生成樣本。

>>> from scipy.stats import qmc
>>> sampler = qmc.Sobol(d=2, scramble=False)
>>> sample = sampler.random_base2(m=3)
>>> sample
array([[0.   , 0.   ],
       [0.5  , 0.5  ],
       [0.75 , 0.25 ],
       [0.25 , 0.75 ],
       [0.375, 0.375],
       [0.875, 0.875],
       [0.625, 0.125],
       [0.125, 0.625]])

使用差異標準計算樣本的質量。

>>> qmc.discrepancy(sample)
0.013882107204860938

要繼續現有設計，可以通過再次調用 random_base2 來獲得加分。或者，您可以跳過一些要點，例如：

>>> _ = sampler.reset()
>>> _ = sampler.fast_forward(4)
>>> sample_continued = sampler.random_base2(m=2)
>>> sample_continued
array([[0.375, 0.375],
       [0.875, 0.875],
       [0.625, 0.125],
       [0.125, 0.625]])

最後，樣本可以縮放到邊界。

>>> l_bounds = [0, 2]
>>> u_bounds = [10, 5]
>>> qmc.scale(sample_continued, l_bounds, u_bounds)
array([[3.75 , 3.125],
       [8.75 , 4.625],
       [6.25 , 2.375],
       [1.25 , 3.875]])