LatinHypercube#
- class scipy.stats.qmc.LatinHypercube(d, *, scramble=True, strength=1, optimization=None, rng=None)[原始碼]#
拉丁超立方抽樣 (LHS)。
拉丁超立方樣本 [1] 在 \([0,1)^{d}\) 中生成 \(n\) 個點。每個單變量邊緣分佈都是分層的,在 \([j/n, (j+1)/n)\) 中恰好放置一個點,其中 \(j=0,1,...,n-1\)。當 \(n << d\) 時,它們仍然適用。
- 參數:
- dint
參數空間的維度。
- scramblebool,可選
當 False 時,將樣本置於多維網格的單元格中心。否則,樣本會隨機放置在網格的單元格內。
注意
設定
scramble=False
並不能確保確定性的輸出。為此,請使用 rng 參數。預設值為 True。
在版本 1.10.0 中新增。
- optimization{None, “random-cd”, “lloyd”}, 可選
是否使用最佳化方案來提高抽樣後的品質。請注意,這是一個後處理步驟,不保證樣本的所有屬性都將被保留。預設值為 None。
random-cd
:隨機排列座標以降低中心化差異。基於中心化差異的最佳樣本會不斷更新。與使用其他差異度量相比,基於中心化差異的抽樣在 2D 和 3D 子投影方面顯示出更好的空間填充穩健性。lloyd
:使用修改後的 Lloyd-Max 演算法擾動樣本。該過程會收斂到均勻間隔的樣本。
在版本 1.8.0 中新增。
在版本 1.10.0 中變更: 新增
lloyd
。- strength{1, 2}, 可選
LHS 的強度。
strength=1
產生純粹的 LHS,而strength=2
產生基於正交陣列的強度為 2 的 LHS [7], [8]。在這種情況下,只能採樣n=p**2
個點,其中p
是一個質數。它也限制d <= p + 1
。預設值為 1。在版本 1.8.0 中新增。
- rng
numpy.random.Generator
, 可選 偽隨機數生成器狀態。當 rng 為 None 時,會使用來自作業系統的熵建立新的
numpy.random.Generator
。numpy.random.Generator
以外的類型會傳遞給numpy.random.default_rng
以實例化Generator
。在版本 1.15.0 中變更: 作為從使用
numpy.random.RandomState
過渡到numpy.random.Generator
的 SPEC-007 轉換的一部分,此關鍵字已從 seed 變更為 rng。在過渡期間,這兩個關鍵字將繼續有效,但一次只能指定一個。在過渡期之後,使用 seed 關鍵字的函數呼叫將發出警告。經過棄用期後,seed 關鍵字將被移除。
另請參閱
註解
當 LHS 用於積分函數 \(f\) 超過 \(n\) 時,LHS 對於幾乎是可加性的被積函數非常有效 [2]。對於 \(n\) 個點的 LHS,積分的變異數始終低於 \(n-1\) 個點的純 MC [3]。關於積分的平均值和變異數,LHS 存在一個中心極限定理 [4],但對於最佳化的 LHS 則不一定,因為存在隨機化。
\(A\) 被稱為強度為 \(t\) 的正交陣列,如果在 \(A\) 的每個 n 行乘 t 列子矩陣中:所有 \(p^t\) 個可能的不同行出現的次數相同。\(A\) 的元素位於集合 \(\{0, 1, ..., p-1\}\) 中,也稱為符號。必須將 \(p\) 設定為質數的限制是為了允許模算術。增加強度會為樣本的子投影增加一些對稱性。強度為 2 時,樣本沿 2D 子投影的對角線對稱。這可能不受歡迎,但另一方面,樣本分散性得到了改善。
強度 1(純 LHS)比強度 0(MC)具有優勢,強度 2 是強度 1 的有用增量。強度 3 的增量較小,而諸如 Sobol'、Halton 等亂碼 QMC 效能更高 [7]。
為了建立強度為 2 的 LHS,正交陣列 \(A\) 通過應用符號集合到自身的一個隨機的、雙射的映射來隨機化。例如,在第 0 列中,所有 0 可能變成 2;在第 1 列中,所有 0 可能變成 1,等等。然後,對於每個列 \(i\) 和符號 \(j\),我們將大小為 \(p\) 的純粹一維 LHS 添加到 \(A^i = j\) 的子陣列中。最終得到的矩陣除以 \(p\)。
參考文獻
[1]Mckay 等人,「比較用於在電腦程式輸出分析中選擇輸入變數值的三種方法」。Technometrics,1979 年。
[2]M. Stein,「使用拉丁超立方抽樣的模擬的大樣本屬性」。Technometrics 29, no. 2: 143-151, 1987 年。
[3]A. B. Owen,「亂碼網格積分的蒙地卡羅變異數」。SIAM Journal on Numerical Analysis 34, no. 5: 1884-1910, 1997 年
[4]Loh, W.-L. 「關於拉丁超立方抽樣」。The annals of statistics 24, no. 5: 2058-2080, 1996 年。
[5]Fang 等人。「電腦實驗的設計和建模」。電腦科學與資料分析系列,2006 年。
[6]Damblin 等人,「空間填充設計的數值研究:拉丁超立方樣本的最佳化和子投影屬性」。Journal of Simulation, 2013 年。
[8]B. Tang,「基於正交陣列的拉丁超立方」。Journal of the American Statistical Association, 1993 年。
[9]Seaholm, Susan K. 等人 (1988)。拉丁超立方抽樣和蒙地卡羅流行病模型的靈敏度分析。Int J Biomed Comput, 23(1-2), 97-112。 DOI:10.1016/0020-7101(88)90067-0
範例
從拉丁超立方生成器產生樣本。
>>> from scipy.stats import qmc >>> sampler = qmc.LatinHypercube(d=2) >>> sample = sampler.random(n=5) >>> sample array([[0.1545328 , 0.53664833], # random [0.84052691, 0.06474907], [0.52177809, 0.93343721], [0.68033825, 0.36265316], [0.26544879, 0.61163943]])
使用差異準則計算樣本的品質。
>>> qmc.discrepancy(sample) 0.0196... # random
樣本可以縮放到邊界。
>>> l_bounds = [0, 2] >>> u_bounds = [10, 5] >>> qmc.scale(sample, l_bounds, u_bounds) array([[1.54532796, 3.609945 ], # random [8.40526909, 2.1942472 ], [5.2177809 , 4.80031164], [6.80338249, 3.08795949], [2.65448791, 3.83491828]])
以下是其他範例,展示了建構 LHS 的其他方法,即使空間覆蓋率更好。
使用基礎 LHS 作為基準線。
>>> sampler = qmc.LatinHypercube(d=2) >>> sample = sampler.random(n=5) >>> qmc.discrepancy(sample) 0.0196... # random
使用 optimization 關鍵字引數來產生具有較低差異的 LHS,但計算成本更高。
>>> sampler = qmc.LatinHypercube(d=2, optimization="random-cd") >>> sample = sampler.random(n=5) >>> qmc.discrepancy(sample) 0.0176... # random
使用 strength 關鍵字引數來產生基於正交陣列的強度為 2 的 LHS。在這種情況下,樣本點的數量必須是質數的平方。
>>> sampler = qmc.LatinHypercube(d=2, strength=2) >>> sample = sampler.random(n=9) >>> qmc.discrepancy(sample) 0.00526... # random
可以組合選項來產生最佳化的中心正交陣列 LHS。最佳化後,結果不能保證強度為 2。
真實世界範例
在 [9] 中,拉丁超立方抽樣 (LHS) 策略被用於對參數空間進行抽樣,以研究流行病模型中每個參數的重要性。這種分析也稱為靈敏度分析。
由於問題的維度很高 (6),因此覆蓋空間在計算上非常昂貴。當數值實驗成本很高時,QMC 能夠進行使用網格可能無法進行的分析。
該模型的六個參數代表疾病的機率、退出的機率和四個接觸機率。作者假設所有參數均為均勻分佈,並產生了 50 個樣本。
使用
scipy.stats.qmc.LatinHypercube
來複製協定,第一步是在單位超立方體中建立一個樣本>>> from scipy.stats import qmc >>> sampler = qmc.LatinHypercube(d=6) >>> sample = sampler.random(n=50)
然後可以將樣本縮放到適當的邊界
>>> l_bounds = [0.000125, 0.01, 0.0025, 0.05, 0.47, 0.7] >>> u_bounds = [0.000375, 0.03, 0.0075, 0.15, 0.87, 0.9] >>> sample_scaled = qmc.scale(sample, l_bounds, u_bounds)
這樣的樣本被用於運行模型 50 次,並建構了一個多項式響應面。這使作者能夠研究每個參數在所有其他參數的可能性範圍內的相對重要性。
在這個電腦實驗中,他們表明,與網格抽樣相比,將響應面上的誤差保持在 2% 以下所需的樣本數量減少了 14 倍。
方法
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 個。reset
()將引擎重設為基礎狀態。