matlab和python中进行拉丁超立方抽样（逆变换法）
什么是拉丁超立方抽樣法？它和蒙特卡羅模擬有什么關系？
逆變換方法(The Inverse Transform Method)采樣
拉丁超立方抽樣
拉丁超立方抽樣(LHS)是一種從多維分布中生成參數值的近隨機樣本的統(tǒng)計方法。抽樣方法常用于構建計算機實驗或進行蒙特卡羅積分。在統(tǒng)計抽樣的上下文中，當(且僅當)每行和每列只有一個樣本時，包含樣本位置的方格為拉丁方格。拉丁超立方體是將這個概念推廣到任意數量的維數，因此每個樣本是包含它的每個軸向超平面上的唯一一個樣本。
在簡單隨機抽樣中，不考慮之前生成的樣本點，生成新的樣本點。我們并不需要事先知道需要多少采樣點。
在拉丁超立方抽樣中，首先必須決定使用多少采樣點，并且記住每個采樣點是在哪一行和哪一列中采樣的。請注意，這種配置類似于在不互相威脅的情況下，在棋盤上有N輛車。
拉丁超立方抽樣也是蒙特卡洛模擬法的一種，它改進了采樣策略能夠做到以較小的采樣規(guī)模獲得較高的采樣精度，在實際應用中具有高效性很受歡迎。核心步驟為分層抽樣、打亂排序。
在累積分布函數中，分層抽樣先將取值空間[0,1]做N等分，得到[0,1/N], [1/N,2/N], ..., [(N-1)/N,N]這N個子層，在每層中隨機選取采樣點，取反得到N個采樣值 ,..., 。排序會打亂采樣值的順序，避免出現第一層抽第一個樣本、第i層抽第f（x）個樣本這種尷尬無意有規(guī)律的場面，保持了樣本間的獨立性。
（注意CDF的值總是位于區(qū)間[0, 1]）
對于m個不同輸入隨機變量，分層采樣可以得到mxN個采樣值，樣本常以樣本矩陣的形式代入模型中進行進一步試驗。分層采樣的樣本值能夠覆蓋輸入隨機變量的整個分布區(qū)間，“抽樣不替換”，樣本重復少；相比蒙特卡羅模擬法的簡單隨機采樣，拉丁超立方抽樣法產生樣本的空間覆蓋率更高，其本質就是樣本的標準差較小。
概率密度函數PDF和累計概率密度函數CDF
拉丁超立方抽樣原理示意圖
蒙特卡羅抽樣技術完全是隨機的，即在輸入分布的范圍內，樣本可以落在任何位置。樣本更有可能從高發(fā)生概率的分布區(qū)域中抽取。在累積分布中，每個蒙特卡羅樣本使用一個0 和 1之間的新的隨機數。在足夠的迭代之后，蒙特卡羅抽樣通過抽樣“重建”輸入分布。但是，當執(zhí)行的迭代次數少的時候，會產生聚集的問題。
拉丁超立方體抽樣是抽樣技術的最新進展，和蒙特卡羅方法相比，它被設計成通過較少迭代次數的抽樣，準確地重建輸入分布。拉丁超立方體抽樣的關鍵是對輸入概率分布進行分層。分層在累積概率尺度（0，1.0）上把累積曲線分成相等的區(qū)間。然后，從輸入分布的每個區(qū)間或“分層”中隨機抽取樣本。抽樣被強制代表每個區(qū)間的值，于是，被強制重建輸入概率分布。
假設我們要在n維向量空間里抽取m個樣本。拉丁超立方體抽樣的步驟是：
（1）將每一維分成互不重迭的m個區(qū)間，使得每個區(qū)間有相同的概率
（通常考慮一個均勻分布，這樣區(qū)間的長度相同）。
（2）在每一維里的每一個區(qū)間中隨機的抽取一個點；
（3）再從每一維里隨機抽出（2）中選取的點，將它們組成向量。
在抽樣時，每個區(qū)間都抽取一個樣本。使用拉丁超立方體方法，樣本更加準確地反映輸入概率分布中值的分布。在拉丁超立方體抽樣中使用的技術是“抽樣不替換”，累積分布的分層數目等于執(zhí)行的迭代次數。
簡單總結一下: 由此可見這個逆采樣變換方法還是很牛逼的，只要你能求出隨機變量X的cdf，基本上就能用均勻分布生成它。缺點是要能求出隨機變量X的cdf。
matlab實現拉丁超立方抽樣（逆變換法）
1. 從離散均勻分布U[1, 1024]中抽樣10000個點
x = unidinv(linspace(0+eps, 1, 10000), 1024); idx = randperm(length(x)); y = x(idx); figure plot(y, '.') figure plot(unique(x), 'o') xlim([0, 1024])
抽樣結果
抽樣結果區(qū)間分布
直方圖（binwidth=3）
可以看到拉丁超立方的好處是覆蓋率高（因為這里即使僅抽取1024個點，也能覆蓋整個抽樣區(qū)間）。
2. 從高斯分布N(5, 3)中抽樣
x = norminv(linspace(0+eps, 1, 10000), 5, 3); idx = randperm(length(x)); y = x(idx); figure plot(y, '.') figure histogram(y, 'binwidth', 0.1)
抽樣結果
直方圖（binwidth=0.1）
3. matlab中常用分布的PDF、CDF和逆CDF
可以看到逆變換法需要用到CDF的逆，這里給出matlab中常見的概率分布。
來自：https://blog.csdn.net/leo2351960/article/details/39207299/
統(tǒng)計工具箱函數 Ⅰ-1 概率密度函數 函數名 對應分布的概率密度函數 betapdf 貝塔分布的概率密度函數 binopdf 二項分布的概率密度函數 chi2pdf 卡方分布的概率密度函數 exppdf 指數分布的概率密度函數 evpdf 最大值型的極值I型分布（Gumbel分布） fpdf f分布的概率密度函數 gampdf 伽瑪分布的概率密度函數 geopdf 幾何分布的概率密度函數 hygepdf 超幾何分布的概率密度函數 normpdf 正態(tài)（高斯）分布的概率密度函數 lognpdf 對數正態(tài)分布的概率密度函數 nbinpdf 負二項分布的概率密度函數 ncfpdf 非中心f分布的概率密度函數 nctpdf 非中心t分布的概率密度函數 ncx2pdf 非中心卡方分布的概率密度函數 poisspdf 泊松分布的概率密度函數 raylpdf 雷利分布的概率密度函數 tpdf 學生氏t分布的概率密度函數 unidpdf 離散均勻分布的概率密度函數 unifpdf 連續(xù)均勻分布的概率密度函數 weibpdf 威布爾分布的概率密度函數 Ⅰ-2 累加分布函數 函數名 對應分布的累加函數 betacdf 貝塔分布的累加函數 binocdf 二項分布的累加函數 chi2cdf 卡方分布的累加函數 expcdf 指數分布的累加函數 evcdf 最大值型的極值I型分布(Gumbel)的累加函數 fcdf f分布的累加函數 gamcdf 伽瑪分布的累加函數 geocdf 幾何分布的累加函數 hygecdf 超幾何分布的累加函數 logncdf 對數正態(tài)分布的累加函數 nbincdf 負二項分布的累加函數 ncfcdf 非中心f分布的累加函數 nctcdf 非中心t分布的累加函數 ncx2cdf 非中心卡方分布的累加函數 normcdf 正態(tài)（高斯）分布的累加函數 poisscdf 泊松分布的累加函數 raylcdf 雷利分布的累加函數 tcdf 學生氏t分布的累加函數 unidcdf 離散均勻分布的累加函數 unifcdf 連續(xù)均勻分布的累加函數 weibcdf 威布爾分布的累加函數 Ⅰ-3 累加分布函數的逆函數 函數名 對應分布的累加分布函數逆函數 betainv 貝塔分布的累加分布函數逆函數 binoinv 二項分布的累加分布函數逆函數 chi2inv 卡方分布的累加分布函數逆函數 expinv 指數分布的累加分布函數逆函數 evinv Gumbel分布的逆函數 finv f分布的累加分布函數逆函數 gaminv 伽瑪分布的累加分布函數逆函數 geoinv 幾何分布的累加分布函數逆函數 hygeinv 超幾何分布的累加分布函數逆函數 logninv 對數正態(tài)分布的累加分布函數逆函數 nbininv 負二項分布的累加分布函數逆函數 ncfinv 非中心f分布的累加分布函數逆函數 nctinv 非中心t分布的累加分布函數逆函數 ncx2inv 非中心卡方分布的累加分布函數逆函數 icdf norminv 正態(tài)（高斯）分布的累加分布函數逆函數 poissinv 泊松分布的累加分布函數逆函數 raylinv 雷利分布的累加分布函數逆函數 tinv 學生氏t分布的累加分布函數逆函數 unidinv 離散均勻分布的累加分布函數逆函數 unifinv 連續(xù)均勻分布的累加分布函數逆函數 weibinv 威布爾分布的累加分布函數逆函數
4. python庫pyDOE
pyDOE是python環(huán)境下進行試驗設計的包。
可以進行拉丁超立方抽樣
https://pythonhosted.org/pyDOE/randomized.html#latin-hypercube
https://github.com/tisimst/pyDOE
https://github.com/clicumu/pyDOE2
基本語法
from pyDOE import lhs lhs(n, [samples, criterion, iterations]) # n: 維度, integer # samples：采樣點數, integer # criterion：default "None" # 一個告訴 lhs 如何采樣點的字符串（默認值：無，它只是隨機化間隔內的點） # “center”或“c”：將采樣間隔內的點居中 # “maximin”或“m”：最大化點之間的最小距離，但將點放在其間隔內的隨機位置 # “centermaximin”或“cm”：與“maximin”相同，但在間隔內居中 # “correlation”或“corr”：最小化最大相關系數 # 輸出設計將所有變量范圍從零縮放到一，然后可以根據用戶的意愿進行轉換 #（例如使用 scipy.stats.distributions ppf（逆累積分布）函數轉換為特定的統(tǒng)計分布。
抽樣
例如，二維標準高斯分布抽樣：
>>> from scipy.stats import norm >>> lhd = lhs(2, samples=5) >>> lhd = norm(loc=0, scale=1).ppf(lhd) # this applies to both factors here
Out[28]:
array([[ 0.75863463, -0.3369805 ],
[-0.00529301, 0.99955206],
[-1.07653807, -0.1303521 ],
[-0.36624876, 0.3676775 ],
[ 0.8589196 , -0.85595459]])
拉丁超立方抽樣原理示意圖
計算過程也很直觀，首先基于cdf的逆進行分層區(qū)間抽樣，然后通過期望分布函數進行變換即可。
python環(huán)境下的統(tǒng)計分布從scipy.stats包導入。包括連續(xù)分布和離散分布。累計分布的逆采用ppf()得到。
https://docs.scipy.org/doc/scipy/tutorial/stats.html
|
pdf() |
Probability density function. |
|
cdf() |
Cumulative distribution function. |
|
ppf() |
Percent point function (inverse of |
從離散分布中抽樣也是一樣：
from scipy.stats import randint lhd = lhs(6, samples=10) lhd = randint(0, 5).ppf(lhd) # sample from U[0, 5)
Out[31]:
array([[1., 1., 2., 1., 3., 2.],
[4., 0., 1., 3., 1., 3.],
[2., 3., 2., 2., 2., 4.],
[2., 4., 0., 4., 4., 1.],
[3., 3., 3., 2., 0., 1.],
可以指定在每個區(qū)間采樣的準則：
lhs(4, criterion='center')
array([[ 0.875, 0.625, 0.875, 0.125],
[ 0.375, 0.125, 0.375, 0.375],
[ 0.625, 0.375, 0.125, 0.625],
[ 0.125, 0.875, 0.625, 0.875]])
還可以為每個維度指定不同的分布參數：
例如，現在，假設我們想將這些設計轉換為正態(tài)分布，均值=[1,2,3,4]，標準差= [0.1,0.5,1,0.25]:
>>> design = lhs(4, samples=10)
>>> from scipy.stats.distributions import norm
>>> means = [1, 2, 3, 4]
>>> stdvs = [0.1, 0.5, 1, 0.25]
>>> for i in xrange(4):
... design[:, i] = norm(loc=means[i], scale=stdvs[i]).ppf(design[:, i])
...
>>> design
array([[ 0.84947986, 2.16716215, 2.81669487, 3.96369414],
[ 1.15820413, 1.62692745, 2.28145071, 4.25062028],
[ 0.99159933, 2.6444164 , 2.14908071, 3.45706066],
[ 1.02627463, 1.8568382 , 3.8172492 , 4.16756309],
[ 1.07459909, 2.30561153, 4.09567327, 4.3881782 ],
[ 0.896079 , 2.0233295 , 1.54235909, 3.81888286],
[ 1.00415 , 2.4246118 , 3.3500082 , 4.07788558],
[ 0.91999246, 1.50179698, 2.70669743, 3.7826346 ],
[ 0.97030478, 1.99322045, 3.178122 , 4.04955409],
[ 1.12124679, 1.22454846, 4.52414072, 3.8707982 ]])
快去成為你想要的樣子！
總結
以上是生活随笔為你收集整理的matlab和python中进行拉丁超立方抽样（逆变换法）的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： Linux学习笔记（五）
- 下一篇： C++ 力扣剑指Offer16-数值的整