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

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

用法:

class  scipy.stats.qmc.LatinHypercube(d, *, scramble=True, strength=1, optimization=None, seed=None)#

拉丁超立方抽樣 (LHS)。

拉丁超立方體樣本 [1] 在 \([0,1)^{d}\) 中生成 \(n\) 點。每個單變量邊際分布都是分層的，在 \([j/n, (j+1)/n)\) 中為 \(j=0,1,...,n-1\) 放置一個點。它們在 \(n << d\) 時仍然適用。

參數 ：：

d： int

參數空間的維度。

scramble： 布爾型，可選

當為 False 時，將樣本集中在多維網格的單元格內。否則，樣本會隨機放置在網格的單元格內。

注意

環境scramble=False不確保確定性輸出。為此，請使用種子範圍。

默認為 True。

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

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

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

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

strength： {1, 2}，可選

LHS 的強度。 strength=1 產生一個普通的 LHS，而strength=2 產生一個基於正交陣列的強度為 2 [7]、[8] 的 LHS。在這種情況下，隻能對n=p**2 點進行采樣，而p 是一個素數。它還約束 d <= p + 1 。默認值為 1。

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

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

注意：

當 LHS 用於對函數 \(f\) 對 \(n\) 進行積分時，LHS 對於幾乎可加的被積數非常有效 [2]。對於 \(n\) 點的 LHS，積分方差始終低於 \(n-1\) 點上的普通 MC [3]。 LHS 在積分的均值和方差上存在中心極限定理 [4]，但由於隨機化，優化的 LHS 不一定有中心極限定理。

如果在 \(A\) 的每個 n-row-by-t 列子矩陣中，\(A\) 被稱為強度 \(t\) 的正交數組：所有 \(p^t\) 可能的不同行出現相同的次數。 \(A\) 的元素位於集合 \(\{0, 1, ..., p-1\}\) 中，也稱為符號。 \(p\) 必須是素數的約束是為了允許模運算。增加強度會增加樣本 sub-projections 的一些對稱性。強度為 2 時，樣本沿 2D sub-projections 的對角線對稱。這可能是不希望的，但另一方麵，樣品分散得到改善。

強度 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 等人，“在計算機代碼輸出分析中選擇輸入變量值的三種方法的比較”。技術計量學，1979 年。

[2]

M. Stein，“使用拉丁超立方抽樣的模擬的大樣本屬性”。技術計量學 29，沒有。 2：143-151，1987。

[3]

A. B. Owen，“加擾網絡正交的蒙特卡洛方差”。 SIAM 數值分析雜誌 34，沒有。 5: 1884-1910, 1997

[4]

羅，W.-L。 “關於拉丁超立方抽樣。”統計年鑒24，沒有。 5：2058-2080，1996。

[5]

方等人。 “計算機實驗的設計和建模”。計算機科學與數據分析係列，2006 年。

[6]

Damblin 等人，“空間填充設計的數值研究：拉丁超立方體樣本和子投影屬性的優化”。模擬雜誌，2013 年。

[7] (1,2)

A. B. Owen，“用於計算機實驗、集成和可視化的正交陣列”。中央統計局，1992 年。

[8]

B. Tang，“基於正交陣列的拉丁超立方體”。美國統計協會雜誌，1993 年。

[9]

蘇珊·K·西霍姆等人。 “拉丁超立方抽樣和蒙特卡羅流行病模型的敏感性分析”。 Int 生物醫學計算雜誌，23(1-2), 97-112，DOI:10.1016/0020-7101(88)90067-0，1988。

例子：

在[9]中，使用拉丁超立方采樣策略對參數空間進行采樣，以研究流行病模型每個參數的重要性。這種分析也稱為敏感性分析。

由於問題的維數很高 (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 倍。

以下是其他示例，展示了構建具有更好空間覆蓋範圍的 LHS 的替代方法。

使用基本 LHS 作為基線。

>>> sampler = qmc.LatinHypercube(d=2)
>>> sample = sampler.random(n=5)
>>> qmc.discrepancy(sample)
0.0196...  # random

使用優化關鍵字參數以較高的計算成本產生具有較小差異的 LHS。

>>> sampler = qmc.LatinHypercube(d=2, optimization="random-cd")
>>> sample = sampler.random(n=5)
>>> qmc.discrepancy(sample)
0.0176...  # random

使用強度關鍵字參數來生成強度為 2 的基於正交數組的 LHS。在這種情況下，樣本點的數量必須是素數的平方。

>>> sampler = qmc.LatinHypercube(d=2, strength=2)
>>> sample = sampler.random(n=9)
>>> qmc.discrepancy(sample)
0.00526...  # random

可以組合選項以產生優化的基於中心正交陣列的 LHS。優化後的結果不能保證強度為2。