1.反變換法

設需產(chǎn)生分布函數(shù)為F(x)的連續(xù)隨機數(shù)X。若已有[0,1]區(qū)間均勻分布隨機數(shù)R,則產(chǎn)生X的反變換公式為：

F(x)=r, 即x=F^-1(r)

反函數(shù)存在條件：如果函數(shù)y=f(x)是定義域D上的單調函數(shù),那么f(x)一定有反函數(shù)存在,且反函數(shù)一定是單調的。分布函數(shù)F(x)為是一個單調遞增函數(shù)，所以其反函數(shù)存在。從直觀意義上理解，因為r一一對應著x,而在[0,1]均勻分布隨機數(shù)R≤r的概率P(R≤r)=r。?因此,連續(xù)隨機數(shù)X≤x的概率P(X≤x)=P(R≤r)=r=F(x)

即X的分布函數(shù)為F(x)。

例子:下面的代碼使用反變換法在區(qū)間[0, 6]上生成隨機數(shù)，其概率密度近似為P(x)=e^-x??

1 import numpy as np 2 import matplotlib.pyplot as plt 3 4 # probability distribution we're trying to calculate 5 p = lambda x: np.exp(-x) 6 7 # CDF of p 8 CDF = lambda x: 1-np.exp(-x) 9 10 # invert the CDF 11 invCDF = lambda x: -np.log(1-x) 12 13 # domain limits 14 xmin = 0 # the lower limit of our domain 15 xmax = 6 # the upper limit of our domain 16 17 # range limits 18 rmin = CDF(xmin) 19 rmax = CDF(xmax) 20 21 N = 10000 # the total of samples we wish to generate 22 23 # generate uniform samples in our range then invert the CDF 24 # to get samples of our target distribution 25 R = np.random.uniform(rmin, rmax, N) 26 X = invCDF(R) 27 28 # get the histogram info 29 hinfo = np.histogram(X,100) 30 31 # plot the histogram 32 plt.hist(X,bins=100, label=u'Samples'); 33 34 # plot our (normalized) function 35 xvals=np.linspace(xmin, xmax, 1000) 36 plt.plot(xvals, hinfo[0][0]*p(xvals), 'r', label=u'p(x)') 37 38 # turn on the legend 39 plt.legend() 40 plt.show()

?一般來說，直方圖的外廓曲線接近于總體X的概率密度曲線。

?2.舍選抽樣法(Rejection Methold)

用反變換法生成隨機數(shù)時，如果求不出F^-1(x)的解析形式或者F(x)就沒有解析形式，則可以用F^-1(x)的近似公式代替。但是由于反函數(shù)計算量較大，有時也是很不適宜的。另一種方法是由Von Neumann提出的舍選抽樣法。下圖中曲線w(x)為概率密度函數(shù)，按該密度函數(shù)產(chǎn)生隨機數(shù)的方法如下：

基本的rejection methold步驟如下：

1. Draw x uniformly from [xmin ?xmax]

2.?Draw x uniformly from?[0, ymax]

3. if y < w(x)，accept the sample, otherwise reject it

4. repeat

即落在曲線w(x)和X軸所圍成區(qū)域內的點接受，落在該區(qū)域外的點舍棄。

例子:下面的代碼使用basic rejection sampling methold在區(qū)間[0, 10]上生成隨機數(shù)，其概率密度近似為P(x)=e^-x??

1 # -*- coding: utf-8 -*- 2 ''' 3 The following code produces samples that follow the distribution P(x)=e^?x 4 for x=[0, 10] and generates a histogram of the sampled distribution. 5 ''' 6 import numpy as np 7 import matplotlib.pyplot as plt 8 9 10 P = lambda x: np.exp(-x) 11 12 # domain limits 13 xmin = 0 # the lower limit of our domain 14 xmax = 10 # the upper limit of our domain 15 16 # range limit (supremum) for y 17 ymax = 1 18 19 N = 10000 # the total of samples we wish to generate 20 accepted = 0 # the number of accepted samples 21 samples = np.zeros(N) 22 count = 0 # the total count of proposals 23 24 # generation loop 25 while (accepted < N): 26 27 # pick a uniform number on [xmin, xmax) (e.g. 0...10) 28 x = np.random.uniform(xmin, xmax) 29 30 # pick a uniform number on [0, ymax) 31 y = np.random.uniform(0,ymax) 32 33 # Do the accept/reject comparison 34 if y < P(x): 35 samples[accepted] = x 36 accepted += 1 37 38 count +=1 39 40 print count, accepted 41 42 # get the histogram info 43 # If bins is an int, it defines the number of equal-width bins in the given range 44 (n, bins)= np.histogram(samples, bins=30) # Returns: n-The values of the histogram,n是直方圖中柱子的高度 45 46 # plot the histogram 47 plt.hist(samples,bins=30,label=u'Samples') # bins=30即直方圖中有30根柱子 48 49 # plot our (normalized) function 50 xvals=np.linspace(xmin, xmax, 1000) 51 plt.plot(xvals, n[0]*P(xvals), 'r', label=u'P(x)') 52 53 # turn on the legend 54 plt.legend() 55 plt.show()

>>>?

99552 10000

3.推廣的舍取抽樣法

從上圖中可以看出，基本的rejection methold法抽樣效率很低，因為隨機數(shù)x和y是在區(qū)間[xmin ?xmax]和區(qū)間[0 ?ymax]上均勻分布的，產(chǎn)生的大部分點不會落在w(x)曲線之下(曲線e^-x的形狀一邊高一邊低，其曲線下的面積占矩形面積的比例很小，則舍選抽樣效率很低)。為了改進簡單舍選抽樣法的效率，可以構造一個新的密度函數(shù)q(x)(called a proposal distribution from which we can readily draw samples)，使它的形狀接近p(x)，并選擇一個常數(shù)k使得kq(x)≥w(x)對于x定義域內的值都成立。對應下圖，首先從分布q(z)中生成隨機數(shù)z₀，然后按均勻分布從區(qū)間[0 ? kq(z₀)]生成一個隨機數(shù)u_0。?if?u_0?> p(z₀) then the sample is rejected,otherwise?u_0?is retained. ?即下圖中灰色區(qū)域內的點都要舍棄。可見，由于隨機點u₀只出現(xiàn)在曲線kq(x)之下，且在q(x)較大處出現(xiàn)次數(shù)較多，從而大大提高了采樣效率。顯然q(x)形狀越接近p(x)，則采樣效率越高。?

?

根據(jù)上述思想，也可以表達采樣規(guī)則如下：

1. Draw x from your proposal distribution q(x)

2. Draw y uniformly from [0 ?1]

3. if y < p(x)/kq(x) , accept the sample, otherwise reject it

4. repeat

下面例子中選擇函數(shù)p(x)=1/(x+1)作為proposal distribution，k=1。曲線1/(x+1)的形狀與e^-x相近。

1 import numpy as np 2 import matplotlib.pyplot as plt 3 4 p = lambda x: np.exp(-x) # our distribution 5 g = lambda x: 1/(x+1) # our proposal pdf (we're choosing k to be 1) 6 CDFg = lambda x: np.log(x +1) # generates our proposal using inverse sampling 7 8 # domain limits 9 xmin = 0 # the lower limit of our domain 10 xmax = 10 # the upper limit of our domain 11 12 # range limits for inverse sampling 13 umin = CDFg(xmin) 14 umax = CDFg(xmax) 15 16 N = 10000 # the total of samples we wish to generate 17 accepted = 0 # the number of accepted samples 18 samples = np.zeros(N) 19 count = 0 # the total count of proposals 20 21 # generation loop 22 while (accepted < N): 23 24 # Sample from g using inverse sampling 25 u = np.random.uniform(umin, umax) 26 xproposal = np.exp(u) - 1 27 28 # pick a uniform number on [0, 1) 29 y = np.random.uniform(0, 1) 30 31 # Do the accept/reject comparison 32 if y < p(xproposal)/g(xproposal): 33 samples[accepted] = xproposal 34 accepted += 1 35 36 count +=1 37 38 print count, accepted 39 40 # get the histogram info 41 hinfo = np.histogram(samples,50) 42 43 # plot the histogram 44 plt.hist(samples,bins=50, label=u'Samples'); 45 46 # plot our (normalized) function 47 xvals=np.linspace(xmin, xmax, 1000) 48 plt.plot(xvals, hinfo[0][0]*p(xvals), 'r', label=u'p(x)') 49 50 # turn on the legend 51 plt.legend() 52 plt.show()

>>>?

24051 10000

可以對比基本的舍取法和改進的舍取法的結果，前者產(chǎn)生符合要求分布的10000個隨機數(shù)運算了99552步，后者運算了24051步，可以看到效率明顯提高。

?

參考：

http://iacs-courses.seas.harvard.edu/courses/am207/blog/lecture-3.html

http://blog.csdn.net/xianlingmao/article/details/7768833

http://blog.sina.com.cn/s/blog_60b44d6a0101l45z.html

http://www.ruanyifeng.com/blog/2015/07/monte-carlo-method.html

轉載于:https://www.cnblogs.com/21207-iHome/p/5266399.html

總結

以上是生活随笔為你收集整理的蒙特卡洛法—非均匀随机数的产生的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內容還不錯，歡迎將生活随笔推薦給好友。