scipy.stats.qmc.

Sobol#

類別 scipy.stats.qmc.Sobol(d, *, scramble=True, bits=None, rng=None, optimization=None)[原始碼]#

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

Sobol' 序列是低差異、準隨機數字。可以使用兩種方法繪製點

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

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

參數:
dint

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

scramblebool, 選項性

如果為 True,則使用 LMS+shift 加擾。否則,不進行加擾。預設為 True。

bitsint, 選項性

產生器的位元數。控制可以產生的最大點數,即 2**bits。最大值為 64。它不對應於傳回類型,傳回類型始終為 np.float64,以防止點重複自身。預設值為 None,為了向後相容性,對應於 30。

在版本 1.9.0 中新增。

optimization{None, “random-cd”, “lloyd”}, 選項性

是否使用最佳化方案來提高取樣後的品質。請注意,這是一個後處理步驟,不保證樣本的所有屬性都將被保留。預設值為 None。

  • random-cd: 座標的隨機排列,以降低中心差異。基於中心差異的最佳樣本會不斷更新。與使用其他差異度量相比,基於中心差異的取樣在 2D 和 3D 子投影方面顯示出更好的空間填充穩健性。

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

在版本 1.10.0 中新增。

rngnumpy.random.Generator, 選項性

偽隨機數產生器狀態。當 rng 為 None 時,會使用作業系統的熵建立新的 numpy.random.Generator。除了 numpy.random.Generator 之外的類型會傳遞給 numpy.random.default_rng 以實例化 Generator

在版本 1.15.0 中變更:作為從使用 numpy.random.RandomState 過渡到 numpy.random.GeneratorSPEC-007 一部分,此關鍵字已從 seed 變更為 rng。在過渡期間,這兩個關鍵字將繼續運作,但一次只能指定一個。在過渡期之後,使用 seed 關鍵字的功能呼叫將發出警告。在棄用期之後,將移除 seed 關鍵字。

筆記

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

Sobol' 序列有很多版本,具體取決於它們的「方向數」。此程式碼使用來自 [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\) 個點後,序列將重複。因此,會引發錯誤。位元數可以使用參數 bits 控制。

參考文獻

[1]

I. M. Sobol’, “The distribution of points in a cube and the accurate evaluation of integrals.” Zh. Vychisl. Mat. i Mat. Phys., 7:784-802, 1967.

[2]

J. Matousek, “On the L2-discrepancy for anchored boxes.” J. of Complexity 14, 527-556, 1998.

[3]

Art B. Owen, “Scrambling Sobol and Niederreiter-Xing points.” Journal of Complexity, 14(4):466-489, December 1998.

[4]

S. Joe and F. Y. Kuo, “Constructing sobol sequences with better two-dimensional projections.” SIAM Journal on Scientific Computing, 30(5):2635-2654, 2008.

[5]

Art B. Owen, “On dropping the first Sobol’ point.” 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]])

方法

fast_forward(n)

將序列快速前進 n 個位置。

integers(l_bounds, *[, u_bounds, n, ...])

l_bounds (包含) 到 u_bounds (不包含) 繪製 n 個整數,或者如果 endpoint=True,則從 l_bounds (包含) 到 u_bounds (包含)。

random([n, workers])

在半開區間 [0, 1) 中繪製 n

random_base2(m)

從 Sobol' 序列繪製點。

重設()

將引擎重設為基本狀態。