matlab和python中进行拉丁超立方抽样（逆变换法）
什么是拉丁超立方抽樣法？它和蒙特卡羅模擬有什么關(guān)系？
逆變換方法(The Inverse Transform Method)采樣
拉丁超立方抽樣
拉丁超立方抽樣(LHS)是一種從多維分布中生成參數(shù)值的近隨機(jī)樣本的統(tǒng)計(jì)方法。抽樣方法常用于構(gòu)建計(jì)算機(jī)實(shí)驗(yàn)或進(jìn)行蒙特卡羅積分。在統(tǒng)計(jì)抽樣的上下文中，當(dāng)(且僅當(dāng))每行和每列只有一個(gè)樣本時(shí)，包含樣本位置的方格為拉丁方格。拉丁超立方體是將這個(gè)概念推廣到任意數(shù)量的維數(shù)，因此每個(gè)樣本是包含它的每個(gè)軸向超平面上的唯一一個(gè)樣本。
在簡(jiǎn)單隨機(jī)抽樣中，不考慮之前生成的樣本點(diǎn)，生成新的樣本點(diǎn)。我們并不需要事先知道需要多少采樣點(diǎn)。
在拉丁超立方抽樣中，首先必須決定使用多少采樣點(diǎn)，并且記住每個(gè)采樣點(diǎn)是在哪一行和哪一列中采樣的。請(qǐng)注意，這種配置類似于在不互相威脅的情況下，在棋盤上有N輛車。
拉丁超立方抽樣也是蒙特卡洛模擬法的一種，它改進(jìn)了采樣策略能夠做到以較小的采樣規(guī)模獲得較高的采樣精度，在實(shí)際應(yīng)用中具有高效性很受歡迎。核心步驟為分層抽樣、打亂排序。
在累積分布函數(shù)中，分層抽樣先將取值空間[0,1]做N等分，得到[0,1/N], [1/N,2/N], ..., [(N-1)/N,N]這N個(gè)子層，在每層中隨機(jī)選取采樣點(diǎn)，取反得到N個(gè)采樣值 ,..., 。排序會(huì)打亂采樣值的順序，避免出現(xiàn)第一層抽第一個(gè)樣本、第i層抽第f（x）個(gè)樣本這種尷尬無(wú)意有規(guī)律的場(chǎng)面，保持了樣本間的獨(dú)立性。
（注意CDF的值總是位于區(qū)間[0, 1]）
對(duì)于m個(gè)不同輸入隨機(jī)變量，分層采樣可以得到mxN個(gè)采樣值，樣本常以樣本矩陣的形式代入模型中進(jìn)行進(jìn)一步試驗(yàn)。分層采樣的樣本值能夠覆蓋輸入隨機(jī)變量的整個(gè)分布區(qū)間，“抽樣不替換”，樣本重復(fù)少；相比蒙特卡羅模擬法的簡(jiǎn)單隨機(jī)采樣，拉丁超立方抽樣法產(chǎn)生樣本的空間覆蓋率更高，其本質(zhì)就是樣本的標(biāo)準(zhǔn)差較小。
概率密度函數(shù)PDF和累計(jì)概率密度函數(shù)CDF
拉丁超立方抽樣原理示意圖
蒙特卡羅抽樣技術(shù)完全是隨機(jī)的，即在輸入分布的范圍內(nèi)，樣本可以落在任何位置。樣本更有可能從高發(fā)生概率的分布區(qū)域中抽取。在累積分布中，每個(gè)蒙特卡羅樣本使用一個(gè)0 和 1之間的新的隨機(jī)數(shù)。在足夠的迭代之后，蒙特卡羅抽樣通過(guò)抽樣“重建”輸入分布。但是，當(dāng)執(zhí)行的迭代次數(shù)少的時(shí)候，會(huì)產(chǎn)生聚集的問(wèn)題。
拉丁超立方體抽樣是抽樣技術(shù)的最新進(jìn)展，和蒙特卡羅方法相比，它被設(shè)計(jì)成通過(guò)較少迭代次數(shù)的抽樣，準(zhǔn)確地重建輸入分布。拉丁超立方體抽樣的關(guān)鍵是對(duì)輸入概率分布進(jìn)行分層。分層在累積概率尺度（0，1.0）上把累積曲線分成相等的區(qū)間。然后，從輸入分布的每個(gè)區(qū)間或“分層”中隨機(jī)抽取樣本。抽樣被強(qiáng)制代表每個(gè)區(qū)間的值，于是，被強(qiáng)制重建輸入概率分布。
假設(shè)我們要在n維向量空間里抽取m個(gè)樣本。拉丁超立方體抽樣的步驟是：
（1）將每一維分成互不重迭的m個(gè)區(qū)間，使得每個(gè)區(qū)間有相同的概率
（通常考慮一個(gè)均勻分布，這樣區(qū)間的長(zhǎng)度相同）。
（2）在每一維里的每一個(gè)區(qū)間中隨機(jī)的抽取一個(gè)點(diǎn)；
（3）再?gòu)拿恳痪S里隨機(jī)抽出（2）中選取的點(diǎn)，將它們組成向量。
在抽樣時(shí)，每個(gè)區(qū)間都抽取一個(gè)樣本。使用拉丁超立方體方法，樣本更加準(zhǔn)確地反映輸入概率分布中值的分布。在拉丁超立方體抽樣中使用的技術(shù)是“抽樣不替換”，累積分布的分層數(shù)目等于執(zhí)行的迭代次數(shù)。
簡(jiǎn)單總結(jié)一下: 由此可見(jiàn)這個(gè)逆采樣變換方法還是很牛逼的，只要你能求出隨機(jī)變量X的cdf，基本上就能用均勻分布生成它。缺點(diǎn)是要能求出隨機(jī)變量X的cdf。
matlab實(shí)現(xiàn)拉丁超立方抽樣（逆變換法）
1. 從離散均勻分布U[1, 1024]中抽樣10000個(gè)點(diǎn)
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])
抽樣結(jié)果
抽樣結(jié)果區(qū)間分布
直方圖（binwidth=3）
可以看到拉丁超立方的好處是覆蓋率高（因?yàn)檫@里即使僅抽取1024個(gè)點(diǎn)，也能覆蓋整個(gè)抽樣區(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)
抽樣結(jié)果
直方圖（binwidth=0.1）
3. matlab中常用分布的PDF、CDF和逆CDF
可以看到逆變換法需要用到CDF的逆，這里給出matlab中常見(jiàn)的概率分布。
來(lái)自：https://blog.csdn.net/leo2351960/article/details/39207299/
統(tǒng)計(jì)工具箱函數(shù) Ⅰ-1 概率密度函數(shù) 函數(shù)名 對(duì)應(yīng)分布的概率密度函數(shù) betapdf 貝塔分布的概率密度函數(shù) binopdf 二項(xiàng)分布的概率密度函數(shù) chi2pdf 卡方分布的概率密度函數(shù) exppdf 指數(shù)分布的概率密度函數(shù) evpdf 最大值型的極值I型分布（Gumbel分布） fpdf f分布的概率密度函數(shù) gampdf 伽瑪分布的概率密度函數(shù) geopdf 幾何分布的概率密度函數(shù) hygepdf 超幾何分布的概率密度函數(shù) normpdf 正態(tài)（高斯）分布的概率密度函數(shù) lognpdf 對(duì)數(shù)正態(tài)分布的概率密度函數(shù) nbinpdf 負(fù)二項(xiàng)分布的概率密度函數(shù) ncfpdf 非中心f分布的概率密度函數(shù) nctpdf 非中心t分布的概率密度函數(shù) ncx2pdf 非中心卡方分布的概率密度函數(shù) poisspdf 泊松分布的概率密度函數(shù) raylpdf 雷利分布的概率密度函數(shù) tpdf 學(xué)生氏t分布的概率密度函數(shù) unidpdf 離散均勻分布的概率密度函數(shù) unifpdf 連續(xù)均勻分布的概率密度函數(shù) weibpdf 威布爾分布的概率密度函數(shù) Ⅰ-2 累加分布函數(shù) 函數(shù)名 對(duì)應(yīng)分布的累加函數(shù) betacdf 貝塔分布的累加函數(shù) binocdf 二項(xiàng)分布的累加函數(shù) chi2cdf 卡方分布的累加函數(shù) expcdf 指數(shù)分布的累加函數(shù) evcdf 最大值型的極值I型分布(Gumbel)的累加函數(shù) fcdf f分布的累加函數(shù) gamcdf 伽瑪分布的累加函數(shù) geocdf 幾何分布的累加函數(shù) hygecdf 超幾何分布的累加函數(shù) logncdf 對(duì)數(shù)正態(tài)分布的累加函數(shù) nbincdf 負(fù)二項(xiàng)分布的累加函數(shù) ncfcdf 非中心f分布的累加函數(shù) nctcdf 非中心t分布的累加函數(shù) ncx2cdf 非中心卡方分布的累加函數(shù) normcdf 正態(tài)（高斯）分布的累加函數(shù) poisscdf 泊松分布的累加函數(shù) raylcdf 雷利分布的累加函數(shù) tcdf 學(xué)生氏t分布的累加函數(shù) unidcdf 離散均勻分布的累加函數(shù) unifcdf 連續(xù)均勻分布的累加函數(shù) weibcdf 威布爾分布的累加函數(shù) Ⅰ-3 累加分布函數(shù)的逆函數(shù) 函數(shù)名 對(duì)應(yīng)分布的累加分布函數(shù)逆函數(shù) betainv 貝塔分布的累加分布函數(shù)逆函數(shù) binoinv 二項(xiàng)分布的累加分布函數(shù)逆函數(shù) chi2inv 卡方分布的累加分布函數(shù)逆函數(shù) expinv 指數(shù)分布的累加分布函數(shù)逆函數(shù) evinv Gumbel分布的逆函數(shù) finv f分布的累加分布函數(shù)逆函數(shù) gaminv 伽瑪分布的累加分布函數(shù)逆函數(shù) geoinv 幾何分布的累加分布函數(shù)逆函數(shù) hygeinv 超幾何分布的累加分布函數(shù)逆函數(shù) logninv 對(duì)數(shù)正態(tài)分布的累加分布函數(shù)逆函數(shù) nbininv 負(fù)二項(xiàng)分布的累加分布函數(shù)逆函數(shù) ncfinv 非中心f分布的累加分布函數(shù)逆函數(shù) nctinv 非中心t分布的累加分布函數(shù)逆函數(shù) ncx2inv 非中心卡方分布的累加分布函數(shù)逆函數(shù) icdf norminv 正態(tài)（高斯）分布的累加分布函數(shù)逆函數(shù) poissinv 泊松分布的累加分布函數(shù)逆函數(shù) raylinv 雷利分布的累加分布函數(shù)逆函數(shù) tinv 學(xué)生氏t分布的累加分布函數(shù)逆函數(shù) unidinv 離散均勻分布的累加分布函數(shù)逆函數(shù) unifinv 連續(xù)均勻分布的累加分布函數(shù)逆函數(shù) weibinv 威布爾分布的累加分布函數(shù)逆函數(shù)
4. python庫(kù)pyDOE
pyDOE是python環(huán)境下進(jìn)行試驗(yàn)設(shè)計(jì)的包。
可以進(jìn)行拉丁超立方抽樣
https://pythonhosted.org/pyDOE/randomized.html#latin-hypercube
https://github.com/tisimst/pyDOE
https://github.com/clicumu/pyDOE2
基本語(yǔ)法
from pyDOE import lhs lhs(n, [samples, criterion, iterations]) # n: 維度, integer # samples：采樣點(diǎn)數(shù), integer # criterion：default "None" # 一個(gè)告訴 lhs 如何采樣點(diǎn)的字符串（默認(rèn)值：無(wú)，它只是隨機(jī)化間隔內(nèi)的點(diǎn)） # “center”或“c”：將采樣間隔內(nèi)的點(diǎn)居中 # “maximin”或“m”：最大化點(diǎn)之間的最小距離，但將點(diǎn)放在其間隔內(nèi)的隨機(jī)位置 # “centermaximin”或“cm”：與“maximin”相同，但在間隔內(nèi)居中 # “correlation”或“corr”：最小化最大相關(guān)系數(shù) # 輸出設(shè)計(jì)將所有變量范圍從零縮放到一，然后可以根據(jù)用戶的意愿進(jìn)行轉(zhuǎn)換 #（例如使用 scipy.stats.distributions ppf（逆累積分布）函數(shù)轉(zhuǎn)換為特定的統(tǒng)計(jì)分布。
抽樣
例如，二維標(biāo)準(zhǔn)高斯分布抽樣：
>>> 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]])
拉丁超立方抽樣原理示意圖
計(jì)算過(guò)程也很直觀，首先基于cdf的逆進(jìn)行分層區(qū)間抽樣，然后通過(guò)期望分布函數(shù)進(jìn)行變換即可。
python環(huán)境下的統(tǒng)計(jì)分布從scipy.stats包導(dǎo)入。包括連續(xù)分布和離散分布。累計(jì)分布的逆采用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.],
可以指定在每個(gè)區(qū)間采樣的準(zhǔn)則：
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]])
還可以為每個(gè)維度指定不同的分布參數(shù)：
例如，現(xiàn)在，假設(shè)我們想將這些設(shè)計(jì)轉(zhuǎn)換為正態(tài)分布，均值=[1,2,3,4]，標(biāo)準(zhǔn)差= [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 ]])
快去成為你想要的樣子！
總結(jié)
以上是生活随笔為你收集整理的matlab和python中进行拉丁超立方抽样（逆变换法）的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。
- 上一篇： Linux学习笔记（五）
- 下一篇： 算法设计与分析——动态规划之矩阵连乘