一、直接采樣

直接采樣的思想是，通過對(duì)均勻分布采樣，實(shí)現(xiàn)對(duì)任意分布的采樣。因?yàn)榫鶆蚍植疾蓸雍貌?#xff0c;我們想要的分布采樣不好采，那就采取一定的策略通過簡(jiǎn)單采取求復(fù)雜采樣。
假設(shè)y服從某項(xiàng)分布p(y)，其累積分布函數(shù)CDF為h(y)，有樣本z~Uniform(0,1)，我們令 z = h(y)，即 y = h(z)^(-1)，結(jié)果y即為對(duì)分布p(y)的采樣。

直接采樣的核心思想在與CDF以及逆變換的應(yīng)用。在原分布p(y)中，如果某個(gè)區(qū)域[a, b]的分布較多，然后對(duì)應(yīng)在CDF曲線中，[h(a), h(b)]的曲線斜率便較大。那么，經(jīng)過逆變換之后，對(duì)y軸（z）進(jìn)行均勻分布采樣時(shí)，分布多的部分（占據(jù)y軸部分多）對(duì)應(yīng)抽樣得到的概率便更大，

局限性

實(shí)際中，所有采樣的分布都較為復(fù)雜，CDF的求解和反函數(shù)的求解都不太可行。

二、拒絕采樣

拒絕采樣是由一個(gè)易于采樣的分布出發(fā)，為一個(gè)難以直接采樣的分布產(chǎn)生采樣樣本的通用算法。既然 p(x) 太復(fù)雜在程序中沒法直接采樣，那么便一個(gè)可抽樣的分布 q(x) 比如高斯分布，然后按照一定的方法拒絕某些樣本，達(dá)到接近 p(x) 分布的目的。

計(jì)算步驟

設(shè)定一個(gè)方便抽樣的函數(shù) q(x)，以及一個(gè)常量 k，使得 p(x) 總在 k*q(x) 的下方。（參考上圖）

• x 軸方向：從 q(x) 分布抽樣得到 a；
• y 軸方向：從均勻分布（0, k*q(a)) 中抽樣得到 u；
• 如果剛好落到灰色區(qū)域： u > p(a), 拒絕， 否則接受這次抽樣；
• 重復(fù)以上過程。

計(jì)算步驟（BN）

• 根據(jù)網(wǎng)絡(luò)指定的先驗(yàn)概率分布生成采樣樣本；
• 拒絕所有與證據(jù)不匹配的樣本；
• 在剩余樣本中對(duì)事件X=x的出現(xiàn)頻繁程度計(jì)數(shù)從而得到估計(jì)概率、

局限性

• 拒絕了太多的樣本！隨著證據(jù)變量個(gè)數(shù)的增多，與證據(jù)e相一致的樣本在所有樣本中所占的比例呈指數(shù)級(jí)下降，所以對(duì)于復(fù)雜問題這種方法是完全不可用的。
• 難以找到合適的k*q(a)，接受概率可能會(huì)很低。

三、重要性采樣（似然加權(quán)）

Likelihood Weighting

重要性采樣主要是用于求一個(gè)復(fù)雜分布p(x)的均值，最后并沒有得到樣本。

重要性采樣的思想是借助一個(gè)易于采樣的簡(jiǎn)單分布q(x)，對(duì)這個(gè)簡(jiǎn)單分布q(x)所得到的樣本全部接受。但是以此得到的樣本肯定不滿足分布p(x)，因此需要對(duì)每一個(gè)樣本附加相應(yīng)的重要性權(quán)重。在重要性采樣中，以p(x0)/q(x0)的值作為每個(gè)樣本的權(quán)重。這樣，當(dāng)樣本和分布p(x)相近時(shí)，對(duì)應(yīng)的權(quán)重大；與分布p(x)相差過多時(shí)，對(duì)應(yīng)的權(quán)重小。這個(gè)方法采樣得到的是帶有重要性權(quán)重的服從q(z)分布的樣本，這個(gè)權(quán)重乘以樣本之后的結(jié)果其實(shí)就是服從p(z)分布的。

通過上述公式，我們可以知道重要性采樣可以用于近似復(fù)雜分布的均值。

四、吉布斯采樣

假設(shè)有一個(gè)例子：E：吃飯、學(xué)習(xí)、打球；時(shí)間T：上午、下午、晚上；天氣W：晴朗、刮風(fēng)、下雨。樣本（E，T，W）滿足一定的概率分布。現(xiàn)要對(duì)其進(jìn)行采樣，如：打球+下午+晴朗。

問題是我們不知道p(E,T,W)，或者說，不知道三件事的聯(lián)合分布。當(dāng)然，如果知道的話，就沒有必要用吉布斯采樣了。但是，我們知道三件事的條件分布。也就是說，p(E|T,W), p(T|E,W), p(W|E,T)。現(xiàn)在要做的就是通過這三個(gè)已知的條件分布，再用Gibbs sampling的方法，得到聯(lián)合分布。
具體方法：首先隨便初始化一個(gè)組合,i.e. 學(xué)習(xí)+晚上+刮風(fēng)，然后依條件概率改變其中的一個(gè)變量。具體說，假設(shè)我們知道晚上+刮風(fēng)，我們給E生成一個(gè)變量，比如，學(xué)習(xí)→吃飯。我們?cè)僖罈l件概率改下一個(gè)變量，根據(jù)學(xué)習(xí)+刮風(fēng)，把晚上變成上午。類似地，把刮風(fēng)變成刮風(fēng)（當(dāng)然可以變成相同的變量）。這樣學(xué)習(xí)+晚上+刮風(fēng)→吃飯+上午+刮風(fēng)。同樣的方法，得到一個(gè)序列，每個(gè)單元包含三個(gè)變量，也就是一個(gè)馬爾可夫鏈。然后跳過初始的一定數(shù)量的單元（比如100個(gè)），然后隔一定的數(shù)量取一個(gè)單元（比如隔20個(gè)取1個(gè)）。這樣sample到的單元，是逼近聯(lián)合分布的。

五、蓄水池采樣

蓄水池抽樣（Reservoir Sampling ），即能夠在o（n）時(shí)間內(nèi)對(duì)n個(gè)數(shù)據(jù)進(jìn)行等概率隨機(jī)抽取，例如：從1000個(gè)數(shù)據(jù)中等概率隨機(jī)抽取出100個(gè)。另外，如果數(shù)據(jù)集合的量特別大或者還在增長(zhǎng)（相當(dāng)于未知數(shù)據(jù)集合總量），該算法依然可以等概率抽樣。

算法步驟：

• 先選取數(shù)據(jù)流中的前k個(gè)元素，保存在集合A中；
• 從第j（k + 1 <= j <= n）個(gè)元素開始，每次先以概率p = k/j選擇是否讓第j個(gè)元素留下。若j被選中，則從A中隨機(jī)選擇一個(gè)元素并用該元素j替換它；否則直接淘汰該元素；
• 重復(fù)步驟2直到結(jié)束，最后集合A中剩下的就是保證隨機(jī)抽取的k個(gè)元素。

六、MCMC算法

隨機(jī)采樣和隨機(jī)模擬：吉布斯采樣Gibbs Sampling

MCMC算法學(xué)習(xí)總結(jié)

【重點(diǎn)】采樣方法（二）MCMC相關(guān)算法介紹及代碼實(shí)現(xiàn)

\(\lim_{p\to\infty}P_{ij}^n\)存在且與i無(wú)關(guān)，記\(\lim_{p\to\infty}P_{ij}^n = \pi(j)\)，我們有：

其中\(\pi = [\pi(1), \pi(2), ... , \pi(j), ...], \sum_{i=0}^{\infty}\pi_i = 1, \pi\)稱為馬氏鏈的平穩(wěn)分布。

所有的 MCMC(Markov Chain Monte Carlo) 方法都是以這個(gè)定理作為理論基礎(chǔ)的。

說明：

該定理中馬氏鏈的狀態(tài)不要求有限，可以是有無(wú)窮多個(gè)的；

定理中的“非周期“這個(gè)概念不解釋，因?yàn)槲覀冇龅降慕^大多數(shù)馬氏鏈都是非周期的；

細(xì)致平穩(wěn)條件

針對(duì)一個(gè)新的分布，如何構(gòu)造對(duì)應(yīng)的轉(zhuǎn)移矩陣？

對(duì)于一個(gè)分布\(\pi(x)\)，根據(jù)細(xì)致平穩(wěn)條件，如果構(gòu)造的轉(zhuǎn)移矩陣P滿足\(\pi(i)P_{ij} = \pi(j)P_{ji}\)，那么\(\pi(x)\)即為該馬氏鏈的平穩(wěn)分布，因此可以根據(jù)這個(gè)條件去構(gòu)造轉(zhuǎn)移矩陣。

通常情況下，初始的轉(zhuǎn)移矩陣\(P\)一般不滿足細(xì)致平穩(wěn)條件，因此我們通過引入接受率構(gòu)造出新的轉(zhuǎn)移矩陣\(P'\)，使其和\(\pi(x)\)滿足細(xì)致平穩(wěn)條件。由此，我們便可以用任何轉(zhuǎn)移概率矩陣（均勻分布、高斯分布）作為狀態(tài)間的轉(zhuǎn)移概率。

如果我們假設(shè)狀態(tài)之間的轉(zhuǎn)移概率是相同的，那么在算法實(shí)現(xiàn)時(shí)，接收率可以簡(jiǎn)單得用\(\pi(j)/\pi(i)\)表示。

Metropolis-Hastings采樣

對(duì)于給定的概率分布p(x),我們希望能有便捷的方式生成它對(duì)應(yīng)的樣本。由于馬氏鏈能收斂到平穩(wěn)分布， 于是一個(gè)很的漂亮想法是：如果我們能構(gòu)造一個(gè)轉(zhuǎn)移矩陣為P的馬氏鏈，使得該馬氏鏈的平穩(wěn)分布恰好是p(x), 那么我們從任何一個(gè)初始狀態(tài)x0出發(fā)沿著馬氏鏈轉(zhuǎn)移, 得到一個(gè)轉(zhuǎn)移序列 x0,x1,x2,?xn,xn+1?,， 如果馬氏鏈在第n步已經(jīng)收斂了，于是我們就得到了 π(x) 的樣本xn,xn+1?。

在馬爾科夫狀態(tài)鏈中，每一個(gè)狀態(tài)代表一個(gè)樣本\(x_n\)，即所有變量的賦值情況。

通過分析MCMC源碼，可以知道：假設(shè)狀態(tài)間的轉(zhuǎn)移概率相同，那么下一個(gè)樣本的采樣會(huì)依賴于上一個(gè)樣本。假設(shè)上一個(gè)樣本所對(duì)應(yīng)的原始分布概率\(\pi(x)\)很小，那么下一個(gè)樣本的接受率很大概率為1；反之如果上一個(gè)樣本的原始分布概率\(\pi(x)\)很大，那么下一個(gè)樣本會(huì)有挺大概率被拒絕。這樣的機(jī)制保證了生成的樣本服從分布\(\pi(x)\)。

從上述分析可以看出，假如初始狀態(tài)的樣本對(duì)應(yīng)的分布概率很小，那么在算法剛開始運(yùn)行時(shí)所產(chǎn)生的樣本（即使是分布概率很小的樣本）很大可能都會(huì)被接收，從而使得算法剛開始運(yùn)行時(shí)采樣的樣本不滿足原始分布\(\pi(x)\)。只要算法采樣到分布概率大的樣本（此時(shí)即為收斂！），那么之后所采樣得到的樣本就會(huì)基本服從原始分布。當(dāng)然，從初始狀態(tài)遍歷到分布概率大的狀態(tài)時(shí)需要運(yùn)行一段時(shí)間，這段過程即為收斂的過程。MCMC算法在收斂之后，保證了在分布概率\(\pi(x)\)大的地方產(chǎn)生更多的樣本，在分布概率\(\pi(x)\)小的地方產(chǎn)生較少的樣本。

一個(gè)馬爾可夫鏈需要經(jīng)過多次的狀態(tài)轉(zhuǎn)移過程采用達(dá)到一個(gè)穩(wěn)定狀態(tài)，這時(shí)候采樣才比較接近真實(shí)的分布。此過程稱為burn in。一般可通過丟棄前面的N個(gè)采樣結(jié)果來(lái)達(dá)到burn in。

疑問

• MCMC的收斂是什么意思？這個(gè)過程中是什么參數(shù)會(huì)更新導(dǎo)致收斂？如何確定何時(shí)收斂？

  收斂過程沒有參數(shù)會(huì)更新，收斂的思想類似于大數(shù)定理。應(yīng)用MCMC算法采樣時(shí)，初始的樣本量少，服從的分布可能和復(fù)雜分布\(\pi(x)\)相差挺遠(yuǎn)，但隨著狀態(tài)轉(zhuǎn)移數(shù)的增加（轉(zhuǎn)移矩陣P的應(yīng)用），根據(jù)上述定理的證明，最終的樣本分布會(huì)逐漸服從復(fù)雜分布\(\pi(x)\)。

• \(\pi\)是每個(gè)狀態(tài)所對(duì)應(yīng)的概率分布嗎？如果是的話，初始選定一個(gè)狀態(tài)后，這個(gè)\(\pi\)如何設(shè)定？或則在MCMC證明過程中，初始\(\pi\)的概率分布如何設(shè)置？

  在MCMC的證明過程中，\(\pi\)是每個(gè)狀態(tài)所對(duì)應(yīng)的概率分布。證明中所給定的初始\(\pi\)應(yīng)該只是為了證明無(wú)論初始樣本符合什么分布，在經(jīng)過一定數(shù)量的轉(zhuǎn)移之后，得到的樣本會(huì)服從復(fù)雜分布\(\pi (x)\)，在實(shí)際代碼實(shí)現(xiàn)中，不用對(duì)這個(gè)\(\pi\)進(jìn)行設(shè)定。

七、代碼

import numpy as np import random import matplotlib.pyplot as plt import pandas as pdf

Rejection Sampling

def f(x):if 0 <= x and x <= 0.25:y = 8 * xelif 0.25 < x and x <= 1:y = (1 - x) * 8/3else:y = 0return y def g(x):if 0 <= x and x <= 1:y = 1else:y = 0return y def plot(fun):X = np.arange(0, 1.0, 0.01)Y = []for x in X:Y.append(fun(x))plt.plot(X, Y)plt.xlabel("x")plt.ylabel("y")plt.show() plot(f) plot(g)

def rejection_sampling(N=10000):M = 3cnt = 0samples = {}while cnt < N:x = random.random()acc_rate = f(x) / (M * g(x))u = random.random()if acc_rate >= u:if samples.get(x) == None:samples[x] = 1else:samples[x] = samples[x] + 1cnt = cnt + 1return samples s = rejection_sampling(100000) X = [] Y = [] for k, v in s.items():X.append(k)Y.append(v) plt.hist(X, bins=100, edgecolor='None')

MCMC Sampling

Metropolis-Hastings Algorithm

參考：MCMC相關(guān)算法介紹及代碼實(shí)現(xiàn)

PI = 3.1415926def get_p(x):# 模擬pi函數(shù)return 1/(2*PI)*np.exp(- x[0]**2 - x[1]**2)def get_tilde_p(x):# 模擬不知道怎么計(jì)算Z的PI，20這個(gè)值對(duì)于外部采樣算法來(lái)說是未知的，對(duì)外只暴露這個(gè)函數(shù)結(jié)果return get_p(x) def domain_random(): #計(jì)算定義域一個(gè)隨機(jī)值return np.random.random()*3.8-1.9 def metropolis(x):new_x = (domain_random(),domain_random()) #新狀態(tài)#計(jì)算接收概率acc = min(1,get_tilde_p((new_x[0],new_x[1]))/get_tilde_p((x[0],x[1])))#使用一個(gè)隨機(jī)數(shù)判斷是否接受u = np.random.random()if u<acc:return new_xreturn x def testMetropolis(counts = 100,drawPath = False):plt.figure()#主要邏輯x = (domain_random(),domain_random()) #x0xs = [x] #采樣狀態(tài)序列for i in range(counts):xs.append(x)x = metropolis(x) #采樣并判斷是否接受#在各個(gè)狀態(tài)之間繪制跳轉(zhuǎn)的線條幫助可視化X1 = [x[0] for x in xs]X2 = [x[1] for x in xs]if drawPath: plt.plot(X1, X2, 'k-',linewidth=0.5)##繪制采樣的點(diǎn)plt.scatter(X1, X2, c = 'g',marker='.')plt.show() testMetropolis(5000)

def metropolis(x):new_x = domain_random()#計(jì)算接收概率acc = min(1,f(new_x)/f(x))#使用一個(gè)隨機(jī)數(shù)判斷是否接受u = np.random.random()if u<acc:return new_xreturn x def testMetropolis(counts = 100,drawPath = False):plt.figure()#主要邏輯x = domain_random()xs = [x] #采樣狀態(tài)序列for i in range(counts):xs.append(x)x = metropolis(x) #采樣并判斷是否接受#在各個(gè)狀態(tài)之間繪制跳轉(zhuǎn)的線條幫助可視化plt.hist(xs, bins=100, edgecolor='None') # plt.plot(xs)plt.show() testMetropolis(100000)

Gibbs Sampling

def partialSampler(x,dim):xes = []for t in range(10): #隨機(jī)選擇10個(gè)點(diǎn)xes.append(domain_random())tilde_ps = []for t in range(10): #計(jì)算這10個(gè)點(diǎn)的未歸一化的概率密度值tmpx = x[:]tmpx[dim] = xes[t]tilde_ps.append(get_tilde_p(tmpx))#在這10個(gè)點(diǎn)上進(jìn)行歸一化操作，然后按照概率進(jìn)行選擇。norm_tilde_ps = np.asarray(tilde_ps)/sum(tilde_ps)u = np.random.random()sums = 0.0for t in range(10):sums += norm_tilde_ps[t]if sums>=u:return xes[t] def gibbs(x):rst = np.asarray(x)[:]path = [(x[0],x[1])]for dim in range(2): #維度輪詢，這里加入隨機(jī)也是可以的。new_value = partialSampler(rst,dim)rst[dim] = new_valuepath.append([rst[0],rst[1]])#這里最終只畫出一輪輪詢之后的點(diǎn)，但會(huì)把路徑都畫出來(lái)return rst,path def testGibbs(counts = 100,drawPath = False):plt.figure()x = (domain_random(),domain_random())xs = [x]paths = [x]for i in range(counts):xs.append([x[0],x[1]])x,path = gibbs(x)paths.extend(path) #存儲(chǔ)路徑p1 = [x[0] for x in paths]p2 = [x[1] for x in paths]xs1 = [x[0] for x in xs]xs2 = [x[1] for x in xs]if drawPath: plt.plot(p1, p2, 'k-',linewidth=0.5)##繪制采樣的點(diǎn)plt.scatter(xs1, xs2, c = 'g',marker='.')plt.show() testGibbs(5000)

轉(zhuǎn)載于:https://www.cnblogs.com/CSLaker/p/9962912.html

總結(jié)

以上是生活随笔為你收集整理的MCMC等采样算法的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

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