一種通用隨機(jī)數(shù)生成方式 -- 逆變換采樣

2021-08-09 16:57 作者:nyasyamorina 0人讀過 | 我要投稿

前言

這里討論的隨機(jī)數(shù)(random number)為隨機(jī)實(shí)數(shù),? 對(duì)于隨機(jī)整數(shù)(random integer)會(huì)叫做...隨機(jī)整數(shù),? 嗯.

實(shí)際上自己并沒有學(xué)過概率論什么的,? 自己本質(zhì)上就是一個(gè)小會(huì)高數(shù)的高中生罷了.? 以下內(nèi)容賭的成分偏多 (

在實(shí)際應(yīng)用中,? 生成任意隨機(jī)數(shù)都可以分為三步:? 1)? 生成隨機(jī)整數(shù).? 2)? 從隨機(jī)整數(shù)得出均勻分布隨機(jī)數(shù).? 3)? 從均勻分布隨機(jī)數(shù)里得出其他分布隨機(jī)數(shù).

隨機(jī)整數(shù)

生成隨機(jī)整數(shù)的方法有很多,? 有簡單的線性同余發(fā)生器(Linear congruential generator, LCG),? 有復(fù)雜的梅森旋轉(zhuǎn)器(Mersenne twister),? 甚至可以使用現(xiàn)實(shí)電路產(chǎn)生真隨機(jī)數(shù).? 這里并不會(huì)詳細(xì)介紹產(chǎn)生隨機(jī)整數(shù)的方法,? 有興趣的可以去翻翻wiki(見底部).

大部分隨機(jī)整數(shù)生成器會(huì)指定整數(shù)分布的區(qū)間?[MIN, MAX],? 給出的隨機(jī)整數(shù)將會(huì)均勻分布在區(qū)間里.? 常用的隨機(jī)整數(shù)生成器區(qū)間為 $%5B0%2C%202%5E%7B32%7D-1%5D$ ,? 或是 $%5B0%2C%202%5E%7B64%7D-1%5D$ ,? 剛好是uint32和uint64可以表示的全部范圍.? 確保隨機(jī)整數(shù)分布的均勻性是非常重要的,? 因?yàn)檫@是為了方便后續(xù)步驟變換隨機(jī)數(shù)的分布方式.

均勻分布隨機(jī)數(shù)

常用的均勻分布隨機(jī)數(shù)的分布區(qū)間為 [0, 1),? 需要注意的是半開區(qū)間,? 也就是說實(shí)數(shù)1是不能取到的.??

把均勻分布的隨機(jī)整數(shù)轉(zhuǎn)化為均勻分布的隨機(jī)數(shù)是非常簡單的,? 把隨機(jī)數(shù)減去MIN再除以MAX-MIN+1就可以得到分布在[0,1)里的隨機(jī)數(shù)了.? 因?yàn)橐婚_始生成的隨機(jī)整數(shù)不是連續(xù)分布的,? 所以使用這種方法轉(zhuǎn)化的隨機(jī)數(shù)也不是連續(xù)分布的,? 但是對(duì)于實(shí)際使用來說已經(jīng)足夠好了.

了解了如何生成均勻分布的隨機(jī)數(shù),? 接下來就是生成其他分布形式的隨機(jī)數(shù)了.? 因?yàn)檫@里是說通用的生成方法,? 所以有必要了解什么是隨機(jī)數(shù)的分布.

概率密度函數(shù)

概率密度函數(shù)(probability density function, pdf)是一個(gè)用來描述隨機(jī)數(shù)分布在給定區(qū)間的可能性的函數(shù).??對(duì)指定區(qū)間下的pdf進(jìn)行積分即為隨機(jī)數(shù)分布在這個(gè)區(qū)間里的概率. 即

$Prob(a%5Cleq%20x%3Cb)%3D%5Cint_a%5Eb%5Crho(x)dx$

例子: 給出一個(gè)隨機(jī)數(shù)生成方式:? 2 * random(),? 對(duì)應(yīng)的pdf為

$%5Crho(x)%3D%5Cleft%5C%7B%5Cbegin%7Bmatrix%7D0.5%2Cx%5Cin%5B0%2C2)%5C%5C0%2Cotherwise%5Cend%7Bmatrix%7D%5Cright.$

則生成隨機(jī)數(shù)會(huì)以50%概率分布在0到1里,? 以100%概率分布在-2到2里,? 以此類推.

另外一個(gè)著名的例子就是正態(tài)分布,? 下面給出標(biāo)準(zhǔn)正態(tài)分布的pdf:

$%5Crho(x)%3D%5Cfrac%7B1%7D%7B%5Csqrt%7B2%5Cpi%7D%7Dexp(-%5Cfrac%7Bx%5E2%7D%7B2%7D)$

符合標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)數(shù)會(huì)以68.27%概率分布在-1到1之間,? 95.45%分布在-2到2之間,? 并且以50%分布在0到+∞之間.

累積分布函數(shù)

累積分布函數(shù)(Cumulative?Distribution?Function, cdf)定義為隨機(jī)分布小于給定數(shù)字的概率,? 根據(jù)pdf的定義,? cdf為

$F(x)%3D%5Cint_%7B-%5Cinfty%7D%5Ex%5Crho(t)dt$

根據(jù)cdf,? 隨機(jī)數(shù)分布在特定區(qū)間下的概率可以表示為

$Prob(a%5Cleq%20x%3Cb)%3DF(b)-F(a)$

例子: 以2*random()生成的隨機(jī)數(shù)的cdf為

$F(x)%3D%5Cleft%5C%7B%5Cbegin%7Bmatrix%7D0%2Cx%3C0%5C%5C0.5x%2Cx%5Cin%5B0%2C2)%5C%5C1%2Cx%5Cgeq2%5Cend%7Bmatrix%7D%5Cright.$

而標(biāo)準(zhǔn)正態(tài)分布的cdf為

$F(x)%3D%5Cfrac%7B1%7D%7B%5Csqrt%7B2%5Cpi%7D%7D%5Cint_%7B-%5Cinfty%7D%5Exexp(-%5Cfrac%7Bt%5E2%7D%7B2%7D)dt$

逆累積分布函數(shù)

因?yàn)閜df為隨機(jī)數(shù)的概率密度函數(shù),? 所以pdf會(huì)符合關(guān)于"概率"的規(guī)則:? pdf必定大于等于0 (不存在負(fù)數(shù)概率),? pdf在全體實(shí)數(shù)上積分為1(概率總和為1).? 因?yàn)閏df是pdf的積分,? pdf的第一條特性定義了cdf必定是單調(diào)遞增的,? pdf的第二條特性定義了cdf有下界0和上界1.

因?yàn)閏df是單調(diào)遞增的,? 所以cdf存在逆函數(shù),? 并且cdf的逆也是單調(diào)遞增的.? cdf有下界0和上界1說明cdf的逆定義域?yàn)閇0, 1]

從均勻隨機(jī)數(shù)里生成任意分布的隨機(jī)數(shù)

暴論:? 已知某種特定分布的逆累積分布函數(shù)inv-cdf,? 則把均勻分布的隨機(jī)數(shù)經(jīng)過inv-cdf處理后得到的隨機(jī)數(shù)符合這種分布,? 即

$Y%3DF%5E%7B-1%7D(X)$

其中X是均勻分布隨機(jī)數(shù),? 則Y符合特定分布.

這是十分容易證明的:?? $Prob(F%5E%7B-1%7D(X)%5Cleq%20Y)%3DProb(X%5Cleq%20F(Y))$ ,? 因?yàn)閄是均勻分布的,? 所以有 $Prob(X%5Cleq%20F(Y))%3DF(Y)$ ,? 如此可知Y滿足特定分布.

例子:? 三角形分布

三角形分布因?yàn)樗膒df形如三角形而得名.? 三角形分布需要給出下限a, 上限b, 和眾數(shù)c.??

三角形分布的pdf定義為

$%5Crho(x)%3D%5Cleft%5C%7B%5Cbegin%7Bmatrix%7D%5Cfrac%7B2(x-a)%7D%7B(b-a)(c-a)%7D%2Cx%5Cin%5Ba%2Cc)%5C%5C%5Cfrac%7B2(b-x)%7D%7B(b-a)(b-c)%7D%2Cx%5Cin%5Bc%2Cb)%5C%5C0%2Cotherwise%5Cend%7Bmatrix%7D%5Cright.$

對(duì)pdf積分求逆得到逆累積分布函數(shù):

$F%5E%7B-1%7D(x)%3D%5Cleft%5C%7B%5Cbegin%7Bmatrix%7Da%2B%5Csqrt%7B(b-a)(c-a)x%7D%2Cx%3C%5Cfrac%7Bc-a%7D%7Bb-a%7D%5C%5Cb-%5Csqrt%7B(b-a)(b-c)(1-x)%7D%2Cotherwise%5Cend%7Bmatrix%7D%5Cright.$

如果X是均勻分布的隨機(jī)數(shù), 那么 $F%5E%7B-1%7D(X)$ 是三角形分布,? 可以簡單地寫一個(gè)jio本進(jìn)行驗(yàn)證

統(tǒng)計(jì)區(qū)間寬度0.01, 采樣次數(shù)1000000

例子2: 正態(tài)分布

正態(tài)分布的pdf為

$%5Crho(x%3B%5Csigma%2C%5Cmu)%3D%5Cfrac%7B1%7D%7B%5Csigma%5Csqrt%7B2%5Cpi%7D%7D%5Cexp(-%5Cfrac%7B(x-%5Cmu)%5E2%7D%7B2%5Csigma%5E2%7D)$

那么cdf為

$F(x%3B%5Csigma%2C%5Cmu)%3D%5Cfrac%7B1%7D%7B2%7D%5Cleft(1%2Berf(%5Cfrac%7Bx-%5Cmu%7D%7B%5Csigma%5Csqrt2%7D)%5Cright)$

其中? $erf(x)%3D%5Cfrac%7B2%7D%7B%5Csqrt%5Cpi%7D%5Cint_0%5Ex%5Cexp(-t%5E2)dt$ ?稱為誤差函數(shù)(error function),? 則cdf的逆為

$F%5E%7B-1%7D(x%3B%5Csigma%2C%5Cmu)%3D%5Csigma%5Csqrt2erf%5E%7B-1%7D(2x-1)%2B%5Cmu$

需要注意到的是,? erf并不是初等函數(shù),? erf的逆也不能直接求出.? 但幸運(yùn)的是,? 目前erf的逆有許多近似的方法.??使用一種快速但低準(zhǔn)確度的逆erf進(jìn)行驗(yàn)證隨機(jī)數(shù)的正確性:

關(guān)于逆變換采樣的討論

1) 逆變換采樣具有通用性,? 所有隨機(jī)分布都有相應(yīng)的分布函數(shù)pdf,? 由pdf推導(dǎo)出cdf的逆即可生成任意分布的隨機(jī)數(shù).? 但是需要手動(dòng)計(jì)算積分和求逆,? 很有可能不存在初等函數(shù)可以表示的解.? 2) 大部分隨機(jī)過程都是由隨機(jī)整數(shù)生成器為基礎(chǔ),? 而隨機(jī)整數(shù)生成器存在周期性,? 比如說LCG只有 $2%5E%7B32%7D$ 長度的周期,? 而梅森旋轉(zhuǎn)器的周期長達(dá) $2%5E%7B19937%7D-1$ .? 使用其他產(chǎn)生特定分布的方法可能會(huì)使得周期產(chǎn)生變化,? 比如使用LCG與C++TR1里的方法產(chǎn)生正態(tài)分布會(huì)使得周期縮短至 $2%5E%7B16%7D$ ,? 這是因?yàn)樵贑++TR1里的生成正態(tài)分布的方法需要取兩次隨機(jī)數(shù)才可以生成一個(gè)正態(tài)分布隨機(jī)數(shù).? 而對(duì)于逆變換采樣,? 不管什么分布都只需要取一次隨機(jī)數(shù).? 3) 逆變換采樣具有對(duì)比拒絕采樣(rejection sampling)更快的運(yùn)算速度.? 對(duì)于拒絕采樣以及其衍生算法,? 需要采用大量的判斷語句,? 而分支對(duì)運(yùn)行速度是不好的 (特別是GPU運(yùn)算).? 因?yàn)槟孀儞Q采樣需要用到大量初等函數(shù),? 所以面對(duì)其他的優(yōu)化算法就顯得特別累贅了.

示例代碼在下面貼出,? 運(yùn)行需要Python3 + numpy + matplotlib.

END.

吐槽

為什么會(huì)有一篇這樣的東西出來呢.? 其實(shí)是最近自己在用C艸寫光追,? 打算自己造輪子的時(shí)候遇上了"如何生成正態(tài)分布隨機(jī)數(shù)"這個(gè)難題(baidu-less).? 然后靈機(jī)一動(dòng)自己發(fā)現(xiàn)了這種計(jì)算方法,? 但是直到專欄寫到一半才看到了原來這種方法叫逆變換采樣.? 嗨,? 還以為發(fā)現(xiàn)了什么新大陸.

關(guān)于上面的示例代碼,? 原本是打算用C艸寫的,? 順便做一波性能分析,? 但是突然發(fā)現(xiàn)C艸畫圖實(shí)在太麻煩的(相對(duì)于matplotlib),? 于是還是滾回了心愛的python下.? 原本代碼是打算發(fā)到gayhub上的,? 但是感覺為了一個(gè)100多行的文件去開一個(gè)倉庫也太小題大做了,? 而且發(fā)百度云有點(diǎn)low,? 也不知道有沒有其他分享代碼的方法,? 就直接在文章結(jié)尾貼出了.

更多資訊:

梅森旋轉(zhuǎn)器:?en.wikipedia.org/wiki/Mersenne_Twister
線性同余發(fā)生器:?en.wikipedia.org/wiki/Linear_congruential_generator
概率密度函數(shù):?en.wikipedia.org/wiki/Probability_density_function
累積分布函數(shù):?en.wikipedia.org/wiki/Cumulative_distribution_function
逆變換采樣:?en.wikipedia.org/wiki/Inverse_transform_sampling
如何產(chǎn)生正態(tài)分布的隨機(jī)數(shù)？zhihu.com/question/29971598
逆誤差函數(shù)(1)?Anthony J. Strecok -?On the Calculation of the Inverse of the Error Function
逆誤差函數(shù)(2) Mathematics of Computation -?Rational Chebyshev Approximations for the Inverse of the Error Function
逆誤差函數(shù)(3)?Mike Giles - Approximating the erfinv function

標(biāo)簽：

最美情侣中文字幕电影,在线麻豆精品传媒,在线网站高清黄,久久黄色视频