史上最全采樣方法詳細解讀與代碼實現
bitcarmanlee 2018-09-20 23:21:08 13503 收藏 40
分類專欄： ml algorithm
版權
1.什么是采樣
在信號系統、數字信號處理中，采樣是每隔一定的時間測量一次聲音信號的幅值，把時間連續的，模擬信號轉換成時間離散、幅值連續的采樣信號。如果采樣的時間間隔相等，這種采樣稱為均勻采樣。
在計算機系統中，有一個重要的問題就是給定一個概率分布p(x) ， 我們如何在計算機中生成它的樣本。平時我們接觸比較多的場景是，給定一堆樣本數據，求出這堆樣本的概率分布p(x)。而采樣剛好是個逆命題：給定一個概率分布p(x)，如何生成滿足條件的樣本？
2.均勻分布采樣
均勻分布式是一種最簡單的分布。在計算機中生成[0,1]之間的偽隨機數序列，就可以看做是一種均勻分布。隨機數生成的方法有很多，比較簡單的一種方式比如：
xn+1=(axn+c) mod mx_{n+1} = (ax_n + c) \ mod \ m
x
n+1
=(ax
n
+c) mod m
當然計算機中產生的隨機數一般都是偽隨機數，不過在絕大部分場景下也夠用了。
3.離散分布采樣
離散分布是比均勻分布稍微復雜一點的情況。例如令p(x)=[0.1,0.2,0.3,0.4]p(x) = [0.1, 0.2, 0.3, 0.4]p(x)=[0.1,0.2,0.3,0.4]，那么我們可以把概率分布的向量看做一個區間段，然后從x∼U(0,1)x \sim U_{(0,1)}x∼U
(0,1)
分布中采樣，判斷x落在哪個區間內。區間長度與概率成正比，這樣當采樣的次數越多，采樣就越符合原來的分布。
舉個例子，假設p(x)=[0.1,0.2,0.3,0.4]p(x) = [0.1, 0.2, 0.3, 0.4]p(x)=[0.1,0.2,0.3,0.4]，分別為"hello", “java”, “python”, "scala"四個詞出現的概率。我們按照上面的算法實現看看最終的結果
import numpy as np
from collections import defaultdict
dic = defaultdict(int)
def sample():
u = np.random.rand()
if u <= 0.1:
dic["hello"] += 1
elif u <= 0.3:
dic["java"] += 1
elif u <= 0.6:
dic["python"] += 1
else:
dic["scala"] += 1
def sampleNtimes():
for i in range(10000):
sample()
for k,v in dic.items():
print k,v
sampleNtimes()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
代碼輸出結果：
python 3028
java 2006
hello 971
scala 3995
1
2
3
4
上面的采樣方法，基本滿足我們的要求。
4. Box-Muller算法
如果概率密度分布函數p(x)p(x)p(x)是連續分布，如果這個分布可以計算累積分布函數(Cumulative Distribution Function, CDF)，可以通過計算CDF的反函數，獲得采樣。
如果隨機變量U1U_1U
1
,U2U_2U
2
是IID獨立同分布，且U1,U2∼Uniform[0,1]U_1, U_2 \sim Uniform[0, 1]U
1
,U
2
∼Uniform[0,1]，那么有:
Z0=−2lnU1−−−−−−√cos(2πU2)Z_0 = \sqrt{-2ln U_1} cos(2 \pi U_2)
Z
0
=
−2lnU
1
cos(2πU
2
)
Z1=−2lnU1−−−−−−√sin(2πU2)Z_1 = \sqrt{-2ln U_1} sin(2 \pi U_2)
Z
1
=
−2lnU
1
sin(2πU
2
)
Z0,Z1Z_0, Z_1Z
0
,Z
1
獨立且服從標准正態分布。
具體證明過程可以見參考文獻2。
Box–Muller變換最初由 George E. P. Box 與 Mervin E. Muller 在1958年提出。George E. P. Box 是統計學的一代大師，統計學中的很多名詞術語都以他的名字命名。Box 之於統計學的家學淵源相當深厚，他的導師是 統計學開山鼻祖 皮爾遜的兒子，英國統計學家Egon Pearson，同時Box還是統計學的另外一位巨擘級奠基人 費希爾 的女婿。統計學中的名言“all models are wrong, but some are useful”（所有模型都是錯的，但其中一些是有用的）也出自Box之口
利用Box-Muller變換生成高斯分布隨機數的方法可以總結為以下步驟：
1.生成兩個隨機數 U1,U2∼U[0,1]U_1, U_2 \sim U[0,1]U
1
,U
2
∼U[0,1]
2.令 R=−2ln(u1)−−−−−−−√R = \sqrt{-2ln(u_1)}R=
−2ln(u
1
)
, θ=2πU2\theta = 2 \pi U_2θ=2πU
2
3.z0=Rcosθ,z1=Rsinθz_0 = Rcos \theta, z_1 = Rsin \thetaz
0
=Rcosθ,z
1
=Rsinθ
按上面的步驟生成正態分布的抽樣：
import numpy as np
from pylab import *
def sample():
x = np.random.rand() # 兩個均勻分布分別為x, y
y = np.random.rand()
R = np.sqrt(-2 * np.log(x))
theta = 2 * np.pi * y
z0 = R * np.cos(theta)
#z1 = R * np.sin(theta)
return z0
def sampleNtimes():
list = []
n = 100000
for i in range(n):
x = sample()
list.append(x)
y = np.reshape(list, n, 1)
hist(y, normed=1, fc='c') # 直方圖
x = arange(-4, 4, 0.1)
plot(x, 1 / np.sqrt(2 * np.pi) * np.exp(-0.5 * x ** 2), 'g', lw=6) # 標准正態分布
xlabel('x', fontsize = 24)
ylabel('p(x)', fontsize = 24)
show()
sampleNtimes()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
最后代碼輸出圖像：
5.拒絕采樣(Rejection Sampling)
假設我們已經可以抽樣高斯分布q(x)（如Box–Muller_transform 算法），我們按照一定的方法拒絕某些樣本，達到接近p(x)分布的目的:
具體的步驟如下：
1.首先，確定常量k，使得p(x)總在kq(x)的下方。
2.x軸的方向：從q(x)分布抽樣取得a。但是a不一定留下，會有一定的幾率被拒絕。
3.y軸的方向：從均勻分布(0, kq(a))中抽樣得到u。如果u>p(a)，也就是落到了灰色的區域中，拒絕，否則接受這次抽樣。
代碼如下：
import numpy as np
import matplotlib.pyplot as plt
def qsample():
"""使用均勻分布作為q(x)，返回采樣"""
return np.random.rand() * 4.
def p(x):
"""目標分布"""
return 0.3 * np.exp(-(x - 0.3) ** 2) + 0.7 * np.exp(-(x - 2.) ** 2 / 0.3)
def rejection(nsamples):
M = 0.72 # 0.8 k值
samples = np.zeros(nsamples, dtype=float)
count = 0
for i in range(nsamples):
accept = False
while not accept:
x = qsample()
u = np.random.rand() * M
if u < p(x):
accept = True
samples[i] = x
else:
count += 1
print "reject count: ", count
return samples
x = np.arange(0, 4, 0.01)
x2 = np.arange(-0.5, 4.5, 0.1)
realdata = 0.3 * np.exp(-(x - 0.3) ** 2) + 0.7 * np.exp(-(x - 2.) ** 2 / 0.3)
box = np.ones(len(x2)) * 0.75 # 0.8
box[:5] = 0
box[-5:] = 0
plt.plot(x, realdata, 'g', lw=3)
plt.plot(x2, box, 'r--', lw=3)
import time
t0 = time.time()
samples = rejection(10000)
t1 = time.time()
print "Time ", t1 - t0
plt.hist(samples, 15, normed=1, fc='c')
plt.xlabel('x', fontsize=24)
plt.ylabel('p(x)', fontsize=24)
plt.axis([-0.5, 4.5, 0, 1])
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
最終結果：
在高維的情況下，Rejection Sampling有兩個問題：
1.合適的q分布很難找
2.很難確定一個合理的k值
導致拒絕率很高。
6.MCMC采樣(Markov Chain Monte Carlo)
MCMC采樣方法中的兩個MC，分別指的是馬爾可夫鏈Markov Chain與蒙特卡洛算法。關於蒙特卡洛算法，前面寫的蒙特卡洛算法有專門分析過，見參考文獻3。關於馬爾可夫鏈見參考文獻4。
之前提到的抽樣方法都是可以並行的。換句話說，后面的樣本跟前面的樣本沒關系。而MCMC是一個鏈式抽樣的過程，每一個抽樣的樣本跟且只跟前面的一個樣本有關系。
由平穩馬爾可夫過程的結論可知：
1.不論初始狀態為啥，只要該馬爾科夫鏈收斂，第n次的抽樣概率p(xn)p(x^n)p(x
n
)一定會收斂到我們預期的分布p(x)p(x)p(x)
2.如果非周期馬爾科夫鏈的狀態轉移矩陣P和概率分布π(x)對於所有的i,j滿足：
π(i)P(i,j)=π(j)P(j,i)\pi(i) P(i, j) = \pi(j) P(j, i)
π(i)P(i,j)=π(j)P(j,i)
則稱概率分布π(x)是狀態轉移矩陣P的平穩分布。
由於馬氏鏈能收斂到平穩分布， 於是一個很的漂亮想法是：如果我們能構造一個轉移矩陣為PPP的馬氏鏈，使得該馬氏鏈的平穩分布恰好是p(x),那么我們從任何一個初始狀態x0x_0x
0
出發沿着馬氏鏈轉移, 得到一個轉移序列Unexpected text node: ' 'Unexpected text node: ' 'x
0
,x
1
,⋯,x
n
,x
n+1
,⋯。如果馬氏鏈在第n步已經收斂，那么xn+1x_{n+1}x
n+1
以后的樣本必然都滿足p(x)分布，就達到了我們的抽樣樣本滿足概率為p(x)分布的條件。
這個絕妙的想法在1953年被 Metropolis想到了，為了研究粒子系統的平穩性質， Metropolis 考慮了物理學中常見的波爾茲曼分布的采樣問題，首次提出了基於馬氏鏈的蒙特卡羅方法，即Metropolis算法，並在最早的計算機上編程實現。Metropolis 算法是首個普適的采樣方法，並啟發了一系列 MCMC方法，所以人們把它視為隨機模擬技術騰飛的起點。 Metropolis的這篇論文被收錄在《統計學中的重大突破》中， Metropolis算法也被遴選為二十世紀的十個最重要的算法之一。
(以上兩小段來自參考文獻5)
由於一般情況下，目標平穩分布π(x)\pi(x)π(x)和某一個馬爾科夫鏈狀態轉移矩陣Q不滿足細致平穩條件
π(i)Q(i,j)≠π(j)Q(j,i)\pi(i) Q(i, j) \ne \pi(j) Q(j,i)
π(i)Q(i,j)
̸
=π(j)Q(j,i)
為了使細致平穩條件成立，可以這樣做：
π(i)Q(i,j)α(i,j)=π(j)Q(j,i)α(j,i)\pi(i) Q(i, j) \alpha(i,j) = \pi(j) Q(j,i) \alpha(j,i)
π(i)Q(i,j)α(i,j)=π(j)Q(j,i)α(j,i)
為了兩邊相等，可以取
α(i,j)=π(j)Q(j,i)\alpha(i, j) = \pi(j) Q(j,i)
α(i,j)=π(j)Q(j,i)
α(j,i)=π(i)Q(i,j)\alpha(j, i) = \pi(i) Q(i, j)
α(j,i)=π(i)Q(i,j)
這樣，轉移矩陣PPP最后就變成為:
P(i,j)=Q(i,j)α(i,j)P(i,j) = Q(i, j) \alpha(i, j)
P(i,j)=Q(i,j)α(i,j)
也就是說，目標矩陣PPP可以通過任意一個馬爾科夫鏈狀態轉移矩陣Q乘以α(i,j)\alpha(i,j)α(i,j)得到，α(i,j)\alpha(i,j)α(i,j)一般稱為接受率，取值在[0,1]之間，即目標矩陣P
可以通過任意一個馬爾科夫鏈狀態轉移矩陣Q以一定的接受率獲得。與前面的拒絕采樣比較，拒絕采樣是以一個常用的分布通過一定的拒絕率獲得一個非常見的分布。而MCMC則是以一個常見的馬爾科夫鏈狀態轉移矩陣Q通過一定的拒絕率得到目標轉移矩陣P，兩者思路類似。
總結一下上面的過程，MCMC采樣算法的過程為：
1.初始化馬爾科夫鏈初始狀態為X0=x0X_0 = x_0X
0
=x
0
2.循環采樣過程如下：
2.1 第t個時刻馬爾科夫鏈狀態為Xt=xtX_t = x_tX
t
=x
t
，由轉移矩陣Q(x)Q(x)Q(x)中采樣得到y∼q(x∣xt)y \sim q(x | x_t)y∼q(x∣x
t
)
2.2 從均勻分布采樣u∼Uniform[0,1]u \sim Uniform[0,1]u∼Uniform[0,1]
2.3 如果u<α(xt,y)=p(y)q(y,xt)u <\alpha(x_t, y) = p(y) q(y, x_t)u<α(x
t
,y)=p(y)q(y,x
t
)，則接受轉移xt→yx_t \to yx
t
→y, Xt+1=yX_{t+1} = yX
t+1
=y
2.4 否則不接受轉移，即Xt+1=xtX_{t+1} = x_tX
t+1
=x
t
以上的 MCMC 采樣算法已經能很漂亮的工作了，不過它有一個小的問題:馬氏鏈Q在轉移的過程中的接受率 α(i,j) 可能偏小，這樣采樣過程中馬氏鏈容易原地踏步，拒絕大量的跳轉，這使得馬氏鏈遍歷所有的狀態空間要花費太長的時間，收斂到平穩分布p(x)的速度太慢。有沒有辦法提升一些接受率呢?
假設 α(i,j)=0.1,α(j,i)=0.2, 此時滿足細致平穩條件，於是
p(i)q(i,j)×0.1=p(j)q(j,i)×0.2
上式兩邊擴大5倍，我們改寫為
p(i)q(i,j)×0.5=p(j)q(j,i)×1
看，我們提高了接受率，而細致平穩條件並沒有打破！這啟發我們可以把細致平穩條件中的α(i,j),α(j,i) 同比例放大，使得兩數中最大的一個放大到1，這樣我們就提高了采樣中的跳轉接受率。所以我們可以取
α(i,j)=min{p(j)q(j,i)p(i)q(i,j),1}\alpha(i, j) = min\{\frac{p(j)q(j, i)}{p(i)q(i,j)}, 1\}
α(i,j)=min{
p(i)q(i,j)
p(j)q(j,i)
,1}
經過以上細微的變化，我們就得到了如下教科書中最常見的 Metropolis-Hastings 算法。
1.初始化馬爾科夫鏈初始狀態為X0=x0X_0 = x_0X
0
=x
0
2.循環采樣過程如下：
2.1 第t個時刻馬爾科夫鏈狀態為Xt=xtX_t = x_tX
t
=x
t
，由轉移矩陣Q(x)Q(x)Q(x)中采樣得到y∼q(x∣xt)y \sim q(x | x_t)y∼q(x∣x
t
)
2.2 從均勻分布采樣u∼Uniform[0,1]u \sim Uniform[0,1]u∼Uniform[0,1]
2.3 如果u<α(xt,y)=min{p(j)q(j,i)p(i)q(i,j),1}u <\alpha(x_t, y) = min\{\frac{p(j)q(j, i)}{p(i)q(i,j)}, 1\}u<α(x
t
,y)=min{
p(i)q(i,j)
p(j)q(j,i)
,1}，則接受轉移xt→yx_t \to yx
t
→y, Xt+1=yX_{t+1} = yX
t+1
=y
2.4 否則不接受轉移，即Xt+1=xtX_{t+1} = x_tX
t+1
=x
t
(以上部分內容來自參考文獻5)
說了這么多的理論，直接看代碼
from __future__ import division
import numpy as np
import matplotlib.pylab as plt
mu = 3
sigma = 10
# 轉移矩陣Q,因為是模擬數字，只有一維，所以Q是個數字(1*1)
def q(x):
return np.exp(-(x-mu)**2/(sigma**2))
# 按照轉移矩陣Q生成樣本
def qsample():
return np.random.normal(mu, sigma)
# 目標分布函數p(x)
def p(x):
return 0.3*np.exp(-(x-0.3)**2) + 0.7* np.exp(-(x-2.)**2/0.3)
def mcmcsample(n = 20000):
sample = np.zeros(n)
sample[0] = 0.5 # 初始化
for i in range(n-1):
qs = qsample() # 從轉移矩陣Q(x)得到樣本xt
u = np.random.rand() # 均勻分布
alpha_i_j = (p(qs) * q(sample[i])) / (p(sample[i]) * qs) # alpha(i, j)表達式
if u < min(alpha_i_j, 1):
sample[i + 1] = qs # 接受
else:
sample[i + 1] = sample[i] # 拒絕
return sample
x = np.arange(0, 4, 0.1)
realdata = p(x)
sampledata = mcmcsample()
plt.plot(x, realdata, 'g', lw = 3) # 理想數據
plt.plot(x,q(x),'r') # Q(x)轉移矩陣的數據
plt.hist(sampledata,bins=x,normed=1,fc='c') # 采樣生成的數據
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
最后結果：
7.吉布斯采樣(Gibbs)
吉布斯采樣（Gibbs sampling）是統計學中用於馬爾科夫蒙特卡洛（MCMC）的一種算法，用於在難以直接采樣時從某一多變量概率分布中近似抽取樣本序列。該序列可用於近似聯合分布、部分變量的邊緣分布或計算積分（如某一變量的期望值）。某些變量可能為已知變量，故對這些變量並不需要采樣。
在高維的情況下，由於接受率α\alphaα的存在，Metropolis-Hastings算法效率不夠高。如果能找到一個轉移矩陣Q使得α=1\alpha = 1α=1，那么抽樣的效率會大大提高。
先看看二維的情況。假設有一個概率分布p(x,y)p(x,y)p(x,y)，x坐標相同的兩個點A(x1,y1),B(x1,y2)A(x_1, y_1), B(x_1, y_2)A(x
1
,y
1
),B(x
1
,y
2
)，有:
p(x1,y1)p(y2∣x1)=p(x1)p(y1∣x1)p(y2∣x1)p(x1,y2)p(y1∣x1)=p(x1)p(y2∣x1)p(y1∣x1)p(x_1, y_1) p(y_2 | x_1) = p(x_1)p(y_1|x_1)p(y_2|x_1) \\p(x_1, y_2) p(y_1 | x_1) = p(x_1)p(y_2|x_1)p(y_1|x_1)
p(x
1
,y
1
)p(y
2
∣x
1
)=p(x
1
)p(y
1
∣x
1
)p(y
2
∣x
1
)
p(x
1
,y
2
)p(y
1
∣x
1
)=p(x
1
)p(y
2
∣x
1
)p(y
1
∣x
1
)
很容易得到
p(x1,y1)p(y2∣x1)=p(x1,y2)p(y1∣x1)p(x_1, y_1) p(y_2 | x_1) = p(x_1, y_2) p(y_1 | x_1)
p(x
1
,y
1
)p(y
2
∣x
1
)=p(x
1
,y
2
)p(y
1
∣x
1
)
即
p(A)p(y2∣x1)=p(B)p(y1∣x1)p(A) p(y_2 | x_1) = p(B) p(y_1 | x_1)
p(A)p(y
2
∣x
1
)=p(B)p(y
1
∣x
1
)
通過上面的式子，我們發現在x=x1x = x_1x=x
1
這條平行於y軸的直線上，如果使用條件分布p(y∣x1)p(y|x_1)p(y∣x
1
)作為任意兩個點之間的轉移概率，那么任意兩個點之間的轉移滿足細致平穩條件。同理，如果在y=y1y = y_1y=y
1
這條直線上任意去兩個點A(x1,y1),C(x2,y1)A(x_1,y_1), C(x_2, y_1)A(x
1
,y
1
),C(x
2
,y
1
)，有如下等式
p(A)p(x2∣y1)=p(C)p(x1∣y1)p(A) p(x_2 | y_1) = p(C) p(x_1 | y_1)
p(A)p(x
2
∣y
1
)=p(C)p(x
1
∣y
1
)
我們可以如下構造平面上任意兩點之間的轉移概率矩陣Q
Q(A→B)=p(yB∣x1)ifxA=xB=x1Q(A→C)=p(xC∣y1)ifyA=yC=y1Q(A→D)=0othersQ(A \to B) = p(y_B|x_1) \qquad if x_A = x_B = x_1 \\Q(A \to C) = p(x_C|y_1) \qquad if y_A = y_C = y_1 \\Q(A \to D) = 0 \qquad others
Q(A→B)=p(y
B
∣x
1
)ifx
A
=x
B
=x
1
Q(A→C)=p(x
C
∣y
1
)ify
A
=y
C
=y
1
Q(A→D)=0others
於是這個二維空間上的馬氏鏈將收斂到平穩分布 p(x,y)。而這個算法就稱為 Gibbs Sampling 算法,是 Stuart Geman 和Donald Geman 這兩兄弟於1984年提出來的，之所以叫做Gibbs Sampling 是因為他們研究了Gibbs random field, 這個算法在現代貝葉斯分析中占據重要位置。
則二維Gibbs Sampling的算法為：
1隨機初始化X0=x0,Y0=y0X_0 = x_0, Y_0 = y_0X
0
=x
0
,Y
0
=y
0
2.循環采樣
2.1 yt+1∼p(y∣xt)y_{t+1} \sim p(y|x_t)y
t+1
∼p(y∣x
t
)
2.2 xt+1∼p(x∣yt+1)x_{t+1} \sim p(x|y_{t+1})x
t+1
∼p(x∣y
t+1
)
總結起來，Gibbs sampling是 MCMC方法的一種特例。當我們采樣幾個變量，但是不知道其聯合分布，但是知道幾個變量相互之間所有的條件分布時，我們可以用其來獲取一些近似觀察樣本。
直接看代碼
from pylab import *
from numpy import *
def pXgivenY(y, m1, m2, s1, s2):
return random.normal(m1 + (y - m2) / s2, s1)
def pYgivenX(x, m1, m2, s1, s2):
return random.normal(m2 + (x - m1) / s1, s2)
def gibbs(N=5000):
k = 20
x0 = zeros(N, dtype=float)
m1 = 10
m2 = 20
s1 = 2
s2 = 3
for i in range(N):
y = random.rand(1)
# 每次采樣需要迭代 k 次
for j in range(k):
x = pXgivenY(y, m1, m2, s1, s2)
y = pYgivenX(x, m1, m2, s1, s2)
x0[i] = x
return x0
def f(x):
return exp(-(x - 10) ** 2 / 10)
# 畫圖
N = 10000
s = gibbs(N)
x1 = arange(0, 17, 1)
hist(s, bins=x1, fc='b')
x1 = arange(0, 17, 0.1)
px1 = zeros(len(x1))
for i in range(len(x1)):
px1[i] = f(x1[i])
plot(x1, px1 * N * 10 / sum(px1), color='g', linewidth=3)
show()
效果如下
參考文獻：
1.https://applenob.github.io/1_MCMC.html#MCMC
2.https://blog.csdn.net/baimafujinji/article/details/6492982
3.https://blog.csdn.net/bitcarmanlee/article/details/82716641 小白都能看懂的蒙特卡洛方法以及python實現
4.https://blog.csdn.net/bitcarmanlee/article/details/82819860 小白都能看懂的馬爾可夫鏈詳解
5.http://www.cnblogs.com/ywl925/archive/2013/06/05/3118875.html
6.https://zh.wikipedia.org/wiki/吉布斯采樣
————————————————
版權聲明：本文為CSDN博主「bitcarmanlee」的原創文章，遵循CC 4.0 BY-SA版權協議，轉載請附上原文出處鏈接及本聲明。
原文鏈接：https://blog.csdn.net/bitcarmanlee/java/article/details/82795137