根據(jù)概率密度函數(shù)生成隨機(jī)數(shù)的代碼

yidiantou 2016-12-14

展開(kāi)全文

我這里并不是要講“偽隨機(jī)”、“真隨機(jī)”這樣的問(wèn)題，而是關(guān)于如何生成服從某個(gè)概率分布的隨機(jī)數(shù)（或者說(shuō) sample）的問(wèn)題。比如，你想要從一個(gè)服從正態(tài)分布的隨機(jī)變量得到 100 個(gè)樣本，那么肯定抽到接近其均值的樣本的概率要大許多，從而導(dǎo)致抽到的樣本很多是集中在那附近的。當(dāng)然，要解決這個(gè)問(wèn)題，我們通常都假設(shè)我們已經(jīng)有了一個(gè) 生成 0 到 1 之間均勻分布的隨機(jī)數(shù)的工具，就好像給我們的結(jié)果那樣，事實(shí)上許多時(shí)候我們也并不太關(guān)心它們是真隨機(jī)數(shù)還是偽隨機(jī)數(shù)，看起來(lái)差不多就行了。

現(xiàn)在再回到我們的問(wèn)題，看起來(lái)似乎是很簡(jiǎn)單的，按照概率分布的話，只要在概率密度大的地方多抽一些樣本不就行了嗎？可是具體要怎么做呢？要真動(dòng)起手來(lái)，似乎有不是那么直觀了。實(shí)際上，這個(gè)問(wèn)題曾經(jīng)也是困擾了我很久，最近又被人問(wèn)起，那我們不妨在這里一起來(lái)總結(jié)一下。為了避免一下子就陷入抽象的公式推導(dǎo)，那就還是從一個(gè)簡(jiǎn)單的具體例子出發(fā)好了，假設(shè)我們要抽樣的概率分布其概率密度函數(shù)為 $p(x) = \frac{1}{9}x^2$ ，并且被限制在區(qū)間 [0, 3] 上，如右上圖所示。

好了，假設(shè)現(xiàn)在我們要抽 100 個(gè)服從這個(gè)分布的隨機(jī)數(shù)，直觀上來(lái)講，抽出來(lái)的接近 3 的數(shù)字肯定要比接近 0 的數(shù)字要多。那究竟要怎樣抽才能得到這樣的結(jié)果呢？由于我們實(shí)際上是不能控制最原始的隨機(jī)數(shù)生成過(guò)程的，我們只能得到一組均勻分布的隨機(jī)數(shù)，而這組隨機(jī)數(shù) 的生成過(guò)程對(duì)于我們完全是透明的，所以，我們能做的只有把這組均勻分布的隨機(jī)數(shù)做一些變換讓他符合我們的需求。找到下手的點(diǎn)了，可是究竟要怎樣變換呢？有一個(gè)變換相信大家都是很熟悉的，假設(shè)我們有一組 [0,1] 之間的均勻分布的隨機(jī)數(shù) X_0 ，那么令 X_1=3X_0 的話， X_1 就是一組在 [0,3] 之間均勻分布的隨機(jī)數(shù)了，不難想象， X_1 等于某個(gè)數(shù) x^* 的概率就是 X_0 等于 x^*/3 的概率（“等于某個(gè)數(shù)的概率”這種說(shuō)法對(duì)于連續(xù)型隨機(jī)變量來(lái)說(shuō)其實(shí)是不合適的，不過(guò)大概可以理解所表達(dá)的意思啦）。似乎有一種可以“逆轉(zhuǎn)回去”的感覺(jué)了。

于是讓我們來(lái)考慮更一般的變換。首先，我們知道 X_1 的概率密度函數(shù)是 $f(x) = 1/3, x\in[0,3]$ ，假設(shè)現(xiàn)在我們令 $Y = \phi (X_1)$ ，不妨先假定 $\phi(\cdot)$ 是嚴(yán)格單調(diào)遞增的函數(shù)，這樣我們可以求其逆函數(shù) $\phi^{-1}(\cdot)$ （也是嚴(yán)格單調(diào)遞增的）?，F(xiàn)在來(lái)看變換后的隨機(jī)變量會(huì)服從一個(gè)什么樣的分布呢？

這里需要小心，因?yàn)檫@里都是連續(xù)型的隨機(jī)變量，并不像離散型隨機(jī)變量那樣可以說(shuō)成“等于某個(gè)值的概率”，因此我們需要轉(zhuǎn)換為概率分布函數(shù)來(lái)處理，也就是求一個(gè)積分啦：

那么 X_1 的概率分布函數(shù)為 $F(x) = \frac{1}{3}x$ 。很顯然小于或等于某個(gè)特定的值 y^* 這件事情是等價(jià)于 $X_1=\phi^{-1}(Y)\leq\phi^{-1}(y^*)$ 這件事情的。換句話說(shuō)， $P(Y\leq y^*)$ 等于 $P(X_1 \leq \phi^{-1}(y^*))$ 。于是，的概率分布函數(shù)就可以得到了：

再求導(dǎo)我們就能得到的概率密度函數(shù)：

這樣一來(lái)，我們就得到了對(duì)于一個(gè)隨機(jī)變量進(jìn)行一個(gè)映射 $\phi(\cdot)$ 之后得到的隨即變量的分布，那么，回到我們剛才的問(wèn)題，我們想讓這個(gè)結(jié)果分布就是我們所求的，然后再反推得 $\phi(\cdot)$ 即可：

經(jīng)過(guò)簡(jiǎn)單的化簡(jiǎn)就可以得到 $\phi^{-1}(y) = \frac{1}{9} y^3$ ，亦即 $\phi(x) = (9x)^{1/3}$ 。也就是說(shuō)，把得到的隨機(jī)數(shù) X_1 帶入到到函數(shù) $\phi(\cdot)$ 中所得到的結(jié)果，就是符合我們預(yù)期要求的隨機(jī)數(shù)啦！讓我們來(lái)驗(yàn)證一下：

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

#!/usr/bin/python   import numpy as np import matplotlib.pyplot as plot   N = 10000 X0 = np.random.rand(N) X1 = 3*X0 Y = np.power(9*X1, 1.0/3)   t = np.arange(0.0, 3.0, 0.01) y = t*t/9  plot.plot(t, y, 'r-', linewidth=1) plot.hist(Y, bins=50, normed=1, facecolor='green', alpha=0.75)plot.show()

這就沒(méi)錯(cuò)啦，目的達(dá)成啦！讓我們來(lái)總結(jié)一下。問(wèn)題是這樣的，我們有一個(gè)服從均勻分布的隨機(jī)變量，它的概率密度函數(shù)為一個(gè)常數(shù) f(x)=C ，如果是 [0,1] 上的分布，那么常數(shù) 就直接等于 1 了?，F(xiàn)在我們要得到一個(gè)隨機(jī)變量使其概率密度函數(shù)為 g(y) ，做法就是構(gòu)造出一個(gè)函數(shù) $\phi(\cdot)$ 滿足（在這里加上了絕對(duì)值符號(hào)，這是因?yàn)?nbsp; $\phi(\cdot)$ 如果不是遞增而是遞減的話，推導(dǎo)的過(guò)程中有一處就需要反過(guò)來(lái)）

反推過(guò)來(lái)就是，對(duì)目標(biāo) 的概率密度函數(shù)求一個(gè)積分（其實(shí)就是得到它的概率分布函數(shù) CDF ，如果一開(kāi)始就拿到的是 CDF 當(dāng)然更好），然后求其反函數(shù)就可以得到需要的變換 $\phi(\cdot)$ 了。實(shí)際上，這種方法有一個(gè)聽(tīng)起來(lái)稍微專業(yè)一點(diǎn)的名字：Inverse Transform Sampling Method 。不過(guò)，雖然看起來(lái)很簡(jiǎn)單，但是實(shí)際操作起來(lái)卻比較困難，因?yàn)閷?duì)于許多函數(shù)來(lái)說(shuō)，求逆是比較困難的，求積分就更困難了，如果寫不出解析解，不得已只能用數(shù) 值方法來(lái)逼近的話，計(jì)算效率就很讓人擔(dān)心了?？墒聦?shí)上也是如此，就連我們最常見(jiàn)的一維標(biāo)準(zhǔn)正態(tài)分布，也很難用這樣的方法來(lái)抽樣，因?yàn)樗母怕拭芏群瘮?shù)

的不定積分沒(méi)有一個(gè)解析形式。這可真是一點(diǎn)也不好玩，費(fèi)了這么大勁，結(jié)果好像什么都干不了?？磥?lái)這個(gè)看似簡(jiǎn)單的問(wèn)題似乎還是比較復(fù)雜的，不過(guò)也不要灰心，至少對(duì)于高斯分布來(lái)說(shuō)，我們還有一個(gè)叫做 Box Muller 的方法可以專門來(lái)做這個(gè)事情。因?yàn)楦咚狗植急容^奇怪，雖然一維的時(shí)候概率分布函數(shù)無(wú)法寫出解析式，但是二維的情況卻可以通過(guò)一些技巧得出一個(gè)解析式來(lái)。

首先我們來(lái)考慮一個(gè)二維的且兩個(gè)維度相互獨(dú)立的高斯分布，它的概率密度函數(shù)為

這個(gè)分布是關(guān)于原點(diǎn)對(duì)稱的，如果考慮使用極坐標(biāo) $(\theta,r)$ （其中 $\theta\in[0,2\pi), r\in[0,\infty)$ ）的話，我們有 $x = r\cos\theta$ ， $y=r\sin\theta$ 這樣的變換。這樣，概率密度函數(shù)是寫成：

注意到在給定的情況下其概率密度是不依賴于 $\theta$ 的，也就是說(shuō)對(duì)于 $\theta$ 來(lái)說(shuō)是一個(gè)均勻分布，這和我們所了解的標(biāo)準(zhǔn)正態(tài)分布也是符合的：在一個(gè)圓上的點(diǎn)的概率是相等的。確定了 $\theta$ 的分布，讓我們?cè)賮?lái)看，用類似于前面的方法：

根據(jù)前面得出的結(jié)論，我現(xiàn)在得到了的概率分布函數(shù)，是不是只要求一下逆就可以得到一個(gè) $\phi(\cdot)$ 了？亦即 $\phi(t) = \sqrt{-2\log (1-t)}$ 。

現(xiàn)在只要把這一些線索串起來(lái)，假設(shè)我們有兩個(gè)相互獨(dú)立的平均分布在 [0,1] 上的隨機(jī)變量 T_1 和 T_2 ，那么 $2\pi T_1$ 就可以得到 $\theta$ 了，而 $\phi(T_2) = \sqrt{-2\log(1-T_2)}$ 就得到了（實(shí)際上，由于 T_2 和 1-T_2 實(shí)際上是相同的分布，所以通常直接寫為 $\sqrt{-2\log T_2}$ ）。再把極坐標(biāo)換回笛卡爾坐標(biāo)：

這樣我們就能得到一個(gè)二維的正態(tài)分布的抽樣了?？梢灾庇^地驗(yàn)證一下，二維不太好畫，就畫成 heatmap 了，看著比較熱的區(qū)域就是概率比較大的，程序如下：

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

#!/usr/bin/python   import numpy as np import matplotlib.pyplot as plot   N = 50000 T1 = np.random.rand(N) T2 = np.random.rand(N)   r = np.sqrt(-2*np.log(T2)) theta = 2*np.pi*T1 X = r*np.cos(theta) Y = r*np.sin(theta)   heatmap, xedges, yedges = np.histogram2d(X, Y, bins=80) extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]   plot.imshow(heatmap, extent=extent) plot.show()

畫出來(lái)的圖像這個(gè)樣子：

不太好看，但是大概的形狀是可以看出來(lái)的。其實(shí)有了二維的高斯分布，再注意到兩個(gè)維度在我們這里是相互獨(dú)立的，那么直接取其中任意一個(gè)維度，就是一個(gè)一維高斯分布了。如下：

如果 $X\sim N(0,1)$ 即服從標(biāo)準(zhǔn)正態(tài)分布的話，則有 $\sigma X+\mu \sim N(\mu, \sigma^2)$ ，也就是說(shuō)，有了標(biāo)準(zhǔn)正態(tài)分布，其他所有的正態(tài)分布的抽樣也都可以完成了。這下總算有點(diǎn)心滿意足了。不過(guò)別急，還有最后一個(gè)問(wèn)題：多元高斯分布。一般最常用不就是二元嗎？二元不是我們一開(kāi)始就推出來(lái)了嗎？推出來(lái)了確實(shí)沒(méi)錯(cuò)，不過(guò)我們考慮的是最簡(jiǎn)單的情形，當(dāng)然同樣可以通過(guò) $\sigma X+\mu$ 這樣的方式來(lái)處理每一個(gè)維度，不過(guò)高維的情形還有一個(gè)需要考慮的就是各個(gè)維度之間的相關(guān)性——我們之前處理的都是兩個(gè)維度相互獨(dú)立的情況。對(duì)于一般的多維正態(tài)分布 $X\sim N(\mathbf{\mu}, \Sigma)$ ，如果各個(gè)維度之間是相互獨(dú)立的，就對(duì)應(yīng)于協(xié)方差矩陣 $\Sigma$ 是一個(gè)對(duì)角陣，但是如果 $\Sigma$ 在非對(duì)角線的地方存在非零元素的話，就說(shuō)明對(duì)應(yīng)的兩個(gè)維度之間存在相關(guān)性。

這個(gè)問(wèn)題還是比較好解決的，高斯分布有這樣的性質(zhì)：類似于一維的情況，對(duì)于多維正態(tài)分布 $X\sim N(\mathbf{\mu}, \Sigma)$ ，那么新的隨機(jī)變量 $X_1=\mathbf{\mu}_1 + LX$ 將會(huì)滿足

所以，對(duì)于一個(gè)給定的高斯分布 $N(\mathbf{\mu}, \Sigma)$ 來(lái)說(shuō)，只要先生成一個(gè)對(duì)應(yīng)維度的標(biāo)準(zhǔn)正態(tài)分布 $X\sim N(0, I)$ ，然后令 $X_1 = \mu+LX$ 即可，其中是對(duì) $\Sigma$ 進(jìn)行 Cholesky Decomposition 的結(jié)果，即 $\Sigma = LL^T$ 。

結(jié)束之前讓我們來(lái)看看 matlab 畫個(gè) 3D 圖來(lái)改善一下心情：

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

N = 50000; T1 = rand(1, N); T2 = rand(1, N);   r = sqrt(-2*log(T2)); theta = 2*pi*T1; X =[r.*cos(theta); r.*sin(theta)]; mu = [1; 2]; Sigma = [5 2; 2 1]; L = chol(Sigma);   X1 = repmat(mu,1, N) + L*X;   nbin = 30;   hist3(X1', [nbin nbin]); set(gcf, 'renderer', 'opengl'); set(get(gca,'child'), 'FaceColor', 'interp', 'CDataMode', 'auto');   [z c] = hist3(X1', [nbin nbin]); [x y] =meshgrid(c{1}, c{2});   figure; surfc(x,y,-z);

下面兩幅圖，哪幅好看一些（注意坐標(biāo)比例不一樣，所以看不出形狀和旋轉(zhuǎn)了）？似乎都不太好看，不過(guò)感覺(jué)還是比前面的 heatmap 要好一點(diǎn)啦！

然后，到這里為止，我們算是把高斯分布弄清楚了，不過(guò)這只是給一個(gè)介紹性的東西，里面的數(shù)學(xué)推導(dǎo)也并不嚴(yán)格，而 Box Muller 也并不是最高效的高斯采樣的算法，不過(guò)，就算我們不打算再深入討論高斯采樣，采樣這個(gè)問(wèn)題本身也還有許多不盡人意的地方，我們推導(dǎo)出來(lái)的結(jié)論可以說(shuō)只能用于一小部分簡(jiǎn)單的分布，連高斯分布都要通過(guò) trick 來(lái)解決，另一些本身連概率密度函數(shù)都寫不出來(lái)或者有各種奇怪?jǐn)?shù)學(xué)特性的分布就更難處理了。所以本文的標(biāo)題里也說(shuō)了，這是上篇，如果什么時(shí)候有機(jī)會(huì)抽出時(shí)間來(lái)寫下篇的話，我將會(huì)介紹一些更加通用和強(qiáng)大的方法，諸如 Rejection Sampling 、Gibbs Sampling 以及 Markov Chain Monte Carlo (MCMC) 等方法。如果你比較感興趣，可以先自行 Google 一下解饞！

本站是提供個(gè)人知識(shí)管理的網(wǎng)絡(luò)存儲(chǔ)空間，所有內(nèi)容均由用戶發(fā)布，不代表本站觀點(diǎn)。請(qǐng)注意甄別內(nèi)容中的聯(lián)系方式、誘導(dǎo)購(gòu)買等信息，謹(jǐn)防詐騙。如發(fā)現(xiàn)有害或侵權(quán)內(nèi)容，請(qǐng)點(diǎn)擊一鍵舉報(bào)。

午夜视频在线网站,日韩视频精品在线,中文字幕精品一区二区三区在线,在线播放精品,1024你懂我懂的旧版人,欧美日韩一级黄色片,一区二区三区在线观看视频

根據(jù)概率密度函數(shù)生成隨機(jī)數(shù)的代碼