逆采樣(Inverse Sampling)和拒絕采樣(Reject Sampling)原理詳解

在蒙特卡羅方法中，有一個關鍵的問題需要解決，即如何基于機率密度函數去采的 n n n個 x x x的樣本集。

逆采樣(Inverse Sampling)和拒絕采樣(Reject Sampling)就是用于解決這個問題的。

其中，對于常見的分布，如均勻分布，高斯分布，指數分布，t分布，F分布，Beta分布，Gamma分布等，可以采用逆采樣的方法進行采樣；不過很多時候，我們的 x x x的機率分布不是常見的分布，這些分布的機率分布函數CDF 不可逆，是以沒有辦法用逆采樣來采樣，這意味着我們沒法友善的得到這些非常見的機率分布的樣本集。拒絕采樣就是用來解決這個問題。

一、逆采樣(Inverse Sampling)

我們知道，對于常見的均勻分布 u n i f o r m ( 0 , 1 ) uniform(0,1) uniform(0,1)是非常容易采樣樣本的，一般通過線性同餘發生器可以很友善的生成(0,1)之間的僞随機數樣本。而其他常見的機率分布，無論是離散的分布還是連續的分布，它們的樣本都可以通過 u n i f o r m ( 0 , 1 ) uniform(0,1) uniform(0,1)的樣本轉換而得。那麼應該如何得到呢？這就是逆采樣。

下面我們以指數分布為例，說明如何通過均勻分布來采樣服從指數分布的樣本集。指數分布的機率密度函數PDF為：

p e x p ( x ) = { λ e x p ( − λ x ) , x ≥ 0 0 , x &lt; 0 (1) p_{exp}(x) = \begin{cases} \lambda exp(-\lambda x) &amp;,x\ge0\\ 0&amp;, x\lt0 \end{cases}\qquad \text{(1)} pexp​(x)={λexp(−λx)0​,x≥0,x<0​(1)

那麼它的機率分布函數CDF為：

F ( x ) = ∫ − ∞ x p ( x ) d x (2) F(x) = \int_{-\infty}^xp(x)dx\qquad \text{(2)} F(x)=∫−∞x​p(x)dx(2)

下圖為指數分布和均勻分布的CDF圖。從左圖上看，在 x ≥ 0 x\ge0 x≥0的部分是一個單調遞增的函數(在定義域上單調非減)，定義域和值域是 [ 0 , + ∞ ) → [ 0 , 1 ) [0,+\infty) \to [0,1) [0,+∞)→[0,1)，在 p ( x ) p(x) p(x)大的地方它增長快，反之亦然。

因為它是唯一映射的（在>0的部分，接下來我們隻考慮這一部分），是以它的反函數可以表示為 F − 1 ( a ) ， a ∈ [ 0 , 1 ) , 值 域 為 [ 0 , + ∞ ) F^{-1}(a)，a\in [0,1), 值域為[0,+\infty) F−1(a)，a∈[0,1),值域為[0,+∞)。(上圖的左圖就是 F ( x ) F(x) F(x)的圖)

因為 F F F單調遞增，是以 F − 1 F^{-1} F−1也是單調遞增的：

x ≤ y ⇒ F ( x ) ≤ F ( y ) a ≤ b ⇒ F − 1 ( a ) ≤ F − 1 ( b ) (3) \begin{aligned} x\le y &amp;\Rightarrow F(x) \le F(y)\\ a\le b &amp;\Rightarrow F^{-1}(a) \le F^{-1}(b)\qquad \text{(3)} \end{aligned} x≤ya≤b​⇒F(x)≤F(y)⇒F−1(a)≤F−1(b)(3)​

利用反函數的定義，我們有：

F − 1 ( a ) &lt; x , i f f    a &lt; F ( x ) (4) F^{-1}(a)&lt;x,iff\ \ a&lt;F(x)\qquad \text{(4)} F−1(a)<x,iff  a<F(x)(4)

接下來，我們定義一下 [ 0 , 1 ] [0,1] [0,1]均勻分布的CDF,這個很好了解：

P ( a ≤ x ) = H ( x ) = { 1 , x ≥ 1 x , 0 ≤ x ≤ 1 0 , x &lt; 0 (5) P(a\le x) = H(x) = \begin{cases} 1 &amp;,x \ge 1\\ x &amp;,0\le x\le1\\ 0&amp;, x\lt0 \end{cases}\qquad \text{(5)} P(a≤x)=H(x)=⎩⎪⎨⎪⎧​1x0​,x≥1,0≤x≤1,x<0​(5)

根據公式（4）和公式（5）,有：

P ( F − 1 ( a ) ≤ x ) = P ( a ≤ F ( x ) ) = H ( F ( x ) ) (6) P(F^{-1}(a) \le x) = P(a \le F(x)) = H(F(x))\qquad \text{(6)} P(F−1(a)≤x)=P(a≤F(x))=H(F(x))(6)

因為 F ( x ) F(x) F(x)的值域 [ 0 , 1 ) [0,1) [0,1)，根據公式（5），公式（6）可以改寫為：

P ( F − 1 ( a ) ≤ x ) = H ( F ( x ) ) = F ( x ) (7) P(F^{-1}(a) \le x) = H(F(x)) = F(x)\qquad \text{(7)} P(F−1(a)≤x)=H(F(x))=F(x)(7)

據 F ( x ) F(x) F(x)的定義，它是exp分布的CDF，是以公式（7）的意思是 F − 1 ( a ) F^{-1}(a) F−1(a)符合exp分布，我們通過 F F F的反函數将一個0到1均勻分布的随機數轉換成了符合exp分布的随機數，注意，以上推導對于CDF可逆的分布都是一樣的。對于exp來說，它的反函數的形式是：

F e x p − 1 ( a ) = − 1 λ ∗ l o g ( 1 − a ) (8) F^{-1}_{exp}(a) = -\frac{1}{\lambda}*log(1-a)\qquad \text{(8)} Fexp−1​(a)=−λ1​∗log(1−a)(8)

具體的映射關系可以看下圖(a)，我們從 y 軸 0-1 的均勻分布樣本（綠色）映射得到了服從指數分布的樣本（紅色）。

我們寫一點代碼來看看效果:

def sampleExp(Lambda = 2,maxCnt = 50000):
    ys = []
    standardXaxis = []
    standardExp = []
    for i in range(maxCnt):
        u = np.random.random()
        y = -1/Lambda*np.log(1-u) #F-1(X)
        ys.append(y)
    for i in range(1000):
        t = Lambda * np.exp(-Lambda*i/100)
        standardXaxis.append(i/100)
        standardExp.append(t)
    plt.plot(standardXaxis,standardExp,'r')
    plt.hist(ys,1000,normed=True)
    plt.show()

sampleExp()
           

最後繪制出來的直方圖可以看出來就是 exp 分布圖，見上圖(b)。可以看到随着采樣數量的變多，機率直方圖和真實的 CDF 就越接近。

以上就是逆采樣的過程。我們的結論是：因為CDF是單調函數（累積的機率隻能越來越大，直到為1），是以，隻要某分布的CDF可逆，那麼就可以通過均勻分布來采樣服從該分布的樣本集。

二、拒絕采樣(Reject Sampling)

2.1 原理

對于那些複雜的不可求逆的分布，拒絕采樣就是針對此類複雜問題的一種随機采樣方法。

我們以求圓周率 π \pi π的例子入手，講解拒絕采樣的思想。通過采樣的方法來計算 π \pi π值，也就是在一個 1 × 1 1×1 1×1的範圍内随機采樣一個點，如果它到原點的距離小于1,則說明它在1/4圓内，則接受它，最後通過接受的占比來計算1/4圓形的面積，進而根據公式反算出預估 π ^ \hat{\pi} π^的值，随着采樣點的增多，最後的結果 π ^ \hat{\pi} π^會越精準。

上面這個例子裡說明一個問題，我們想求一個空間裡均勻分布的集合面積，可以嘗試在更大範圍内按照均勻分布随機采樣，如果采樣點在集合中，則接受，否則拒絕。最後的接受機率就是集合在”更大範圍“的面積占比。

接下來，我們來形式化地說明拒絕采樣。

給定一個機率分布 p ( z ) = 1 Z p p ~ ( z ) p(z)=\frac{1}{Z_p} \tilde{p}(z) p(z)=Zp​1​p~​(z)，其中， p ~ ( z ) \tilde{p}(z) p~​(z)已知， Z p Z_p Zp​為歸一化常數，未知。要對該分布 p ( z ) p(z) p(z)進行拒絕采樣，首先需要借用一個簡單的參考分布(proposal distribution)，記為 q ( x ) q(x) q(x)，該分布的采樣易于實作，如均勻分布、高斯分布。然後引入常數 k k k，使得對所有的的 z z z，滿足 k q ( z ) ≥ p ~ ( z ) kq(z)\geq\tilde{p}(z) kq(z)≥p~​(z)，如下圖所示，紅色的曲線為 p ~ ( z ) \tilde{p}(z) p~​(z)，藍色的曲線為 k q ( z ) kq(z) kq(z)。

在每次采樣中，首先從 q ( z ) q(z) q(z)采樣一個數值 z 0 z_0 z0​，然後在區間 [ 0 , k q ( z 0 ) ] [0,kq(z_0)] [0,kq(z0​)]進行均勻采樣，得到 u 0 u_0 u0​。如果 u 0 &lt; p ~ ( z 0 ) u_0&lt;\tilde{p}(z_0) u0​<p~​(z0​)，則保留該采樣值，否則舍棄該采樣值。最後得到的資料就是對 p ( z ) {p}(z) p(z)分布的一個近似采樣。 結合圖，直覺來了解上述的過程：在 x = z 0 x=z_0 x=z0​這條線上，從 [ 0 , k q ( z 0 ) ] [0,kq(z_0)] [0,kq(z0​)]均勻采樣中一個值，如果這個值小于 p ~ ( z 0 ) \tilde{p}(z_0) p~​(z0​)，即這個均勻采樣的這個值落在了 p ~ ( z 0 ) \tilde{p}(z_0) p~​(z0​)的下方，我們就接受 z 0 z_0 z0​這個采樣值。

我們知道，每次采樣的接受機率計算如下：

p ( a c c e p t ) = ∫ p ~ ( z ) k q ( z ) q ( z ) d z = 1 k ∫ p ~ ( z ) d z p(accept)=\displaystyle\int{\dfrac{\tilde{p}(z)}{kq(z)}}q(z)dz = \dfrac{1}{k} \int{\tilde{p}(z)}dz p(accept)=∫kq(z)p~​(z)​q(z)dz=k1​∫p~​(z)dz

是以，為了提高接受機率，防止舍棄過多的采樣值而導緻采樣效率低下， k k k的選取應該在滿足 k q ( z ) ≥ p ~ ( z ) kq(z)\geq\tilde{p}(z) kq(z)≥p~​(z)的基礎上盡可能小。

拒絕采樣問題可以這樣了解， p ~ ( z ) \tilde{p}(z) p~​(z)與 x x x軸之間的區域為要估計的問題，類似于上面提到的圓形區域， k q ( z ) kq(z) kq(z)與 x x x軸之間的區域為參考區域，類似于上面提到的正方形。由于 k q ( z ) kq(z) kq(z)與 x x x軸之間的區域面積為 k k k(原來的機率密度函數的面積為1，現在擴大了 k k k倍，故 k × 1 k\times1 k×1)，是以， p ~ ( z ) \tilde{p}(z) p~​(z)與 x x x軸之間的區域面積除以 k k k即為對 p ( z ) p(z) p(z)的估計。在每一個采樣點，以 [ 0 , k q ( z 0 ) ] [0,kq(z_0)] [0,kq(z0​)]為界限，落在 p ~ ( z ) \tilde{p}(z) p~​(z)曲線以下的點就是服從 p ( z ) p(z) p(z)分布的點。

2.2 實驗

我們寫一段代碼來看看拒絕采樣是如何對複雜分布進行采樣的。假設，我們要采樣的分布是：

p ~ ( z ) = 0.3 e x p ( − ( z − 0.3 ) 2 ) + 0.7 e x p ( − ( z − 2 ) 2 / 0.3 ) \tilde{p}(z)=0.3exp(-(z-0.3)^2)+0.7exp(-(z-2)^2/0.3) p~​(z)=0.3exp(−(z−0.3)2)+0.7exp(−(z−2)2/0.3)

其歸一化常數 Z p ≈ 1.2113 Z_p \approx 1.2113 Zp​≈1.2113,參考分布為高斯分布 q ( z ) = G a s s i a n ( 1.4 , 1.2 ) q(z)=Gassian(1.4,1.2) q(z)=Gassian(1.4,1.2)，其中均值和方差是經過計算和嘗試得到的，以滿足 k q ( z ) ≥ p ~ ( z ) kq(z)\geq\tilde{p}(z) kq(z)≥p~​(z)。python代碼如下：

import numpy as np
import matplotlib.pyplot as plt

def f1(x):
    return (0.3*np.exp(-(x-0.3)**2) + 0.7* np.exp(-(x-2.)**2/0.3))/1.2113
x = np.arange(-4.,6.,0.01)
plt.plot(x,f1(x),color = "red")

size = int(1e+07)
sigma = 1.2
z = np.random.normal(loc = 1.4,scale = sigma, size = size)
qz = 1/(np.sqrt(2*np.pi)*sigma)*np.exp(-0.5*(z-1.4)**2/sigma**2)
k = 2.5
#z = np.random.uniform(low = -4, high = 6, size = size)
#qz = 0.1
#k = 10
u = np.random.uniform(low = 0, high = k*qz, size = size)

pz =  0.3*np.exp(-(z-0.3)**2) + 0.7* np.exp(-(z-2.)**2/0.3)
sample = z[pz >= u]
plt.hist(sample,bins=150, normed=True, edgecolor='black')
plt.show()
           

得到的采樣分布如下圖：

可以看到，采樣結果完全符合原分布。另外如果我們把參考分布換為均勻分布(代碼中z,q,k換為注釋部分)，仍然得到較好的采樣結果，如下圖:

可以看到，采樣結果也是完全符合原分布。

以上的實驗說明了拒絕采樣的有效性。

2.3 注意事項

通過在2.1節的分析，我們知道：必須知道複雜分布的機率密度函數PDF，才可以進行拒絕采樣。然而，在現實情況中：

1）對于一些二維分布 p ( x , y ) p(x,y) p(x,y)，有時候我們隻能得到條件分布 p ( x ∣ y ) p(x|y) p(x∣y)和 p ( y ∣ x ) p(y|x) p(y∣x)，卻很難得到二維分布的機率密度函數 p ( x , y ) p(x,y) p(x,y)的一般形式，這時我們無法用拒絕采樣得到其樣本集。

2）對于一些高維的複雜非常見分布 p ( x 1 , x 2 , . . . , x n ) p(x_1,x_2,...,x_n) p(x1​,x2​,...,xn​)，我們要找到一個合适的 q ( x ) q(x) q(x)和 k k k非常困難。

是以，實際上，我們仍然要找到一種方法可以解決如何友善得到各種複雜機率分布的對應的采樣樣本集的問題。馬爾科夫鍊有能力幫助找到這些複雜機率分布的對應的采樣樣本集。而這也是MCMC采樣的基礎，我們回頭會講解。

參考文獻

【1】拒絕采樣(reject sampling)原理詳解

【2】采樣方法（一）

【3】蒙特卡洛采樣之拒絕采樣（Reject Sampling）

【4】一文了解采樣方法