針對離散數(shù)據(jù)概率分布的MCMC算法python實現(xiàn)

---

對于mcmc算法，如何選擇狀態(tài)轉(zhuǎn)移矩陣對實驗結(jié)果是否有影響，設(shè)計下面幾組實驗

輸入為p=[0.9,0.05,0.05]，取q為[1/3,1/3,1/3],輸出10000000個樣本，觀察樣本的概率分布，代碼如下：

#coding=utf-8 ''' Created on 2016年9月14日 基本MCMC算法以及M-H算法的python實現(xiàn) @author: whz p:輸入的概率分布，離散情況采用元素為概率值的數(shù)組表示 N:認(rèn)為迭代N次馬爾可夫鏈?zhǔn)諗?Nlmax:馬爾可夫鏈?zhǔn)諗亢笥秩〉姆膒分布的樣本數(shù) isMH:是否采用MH算法，默認(rèn)為True ''' from __future__ import division import matplotlib.pyplot as plt import matplotlib as mpl import numpy as np from array import array def mcmc(p,N=10000,Nlmax=10000,isMH=True):A = np.array([[1.0 / len(p) for x in range(len(p))] for y in range(len(p))], dtype=np.float64) X0 = np.random.randint(len(p))count = 0samplecount = 0L = array("d",[X0])l = array("d")while True:X = L[samplecount]cur = np.argmax(np.random.multinomial(1,A[X]))count += 1if isMH:a = (p[cur]*A[cur][X])/(p[X]*A[X][cur])alpha = min(a,1)else:alpha = p[cur]*A[cur][X]u = np.random.uniform(0,1) if u<alpha:samplecount += 1L.append(cur)if count>N:l.append(cur)if len(l)>=Nlmax:breakelse:continueLa = np.frombuffer(L)la = np.frombuffer(l)return La,la def count(q,n):L = array("d")l1 = array("d")l2 = array("d")for e in q:L.append(e)for e in range(n):l1.append(L.count(e))for e in l1:l2.append(e/sum(l1))return l1,l2p = np.array([0.9,0.05,0.05]) a = mcmc(p,Nlmax=10000000)[1] l1 = ['state%d'% x for x in range(len(p))] plt.pie(count(a,len(p))[0],labels=l1,labeldistance=0.3,autopct='%1.2f%%') plt.title("sampling") plt.show()

圖1可見該樣本收斂于0.7691,0.1154,0.1155，與我們的輸入的差距相差甚遠(yuǎn)

2. 把輸入更改成p=[0.8,0.1,0.1]，以同樣的迭代次數(shù)，輸出同樣數(shù)量樣本可得到
圖2也是同樣有很大差距，可以用方差來衡量差距

3. 把輸入更改成p=[0.6,0.2,0.2]，以同樣的迭代次數(shù)，輸出同樣數(shù)量樣本可得到
圖3同樣還是有著很大的差距

4. 最后把輸入更改成p=[1/3,1/3,1/3]，以同樣的迭代次數(shù)，輸出同樣數(shù)量樣本可得到
圖4這次準(zhǔn)確率卻很高

由以上四組實驗，猜測使構(gòu)造的轉(zhuǎn)移矩陣的每一行與輸入的分布p(x)相等時，可以達(dá)到較高的精度，所以，如下更改代碼：

#coding=utf-8 ''' Created on 2016年9月15日 基本MCMC算法以及M-H算法的python實現(xiàn)(離散數(shù)據(jù)的情況） @author: whz p:輸入的概率分布，離散情況采用元素為概率值的數(shù)組表示 構(gòu)造的轉(zhuǎn)移矩陣每一行都與輸入的p相同 N:認(rèn)為迭代N次馬爾可夫鏈?zhǔn)諗?Nlmax:馬爾可夫鏈?zhǔn)諗亢笥秩〉姆膒分布的樣本數(shù) isMH:是否采用MH算法，默認(rèn)為True ''' from __future__ import division import matplotlib.pyplot as plt import numpy as np from array import array def mcmc(p,N=10000,Nlmax=10000,isMH=True):A = np.array([p for y in range(len(p))], dtype=np.float64) X0 = np.random.randint(len(p))count = 0samplecount = 0L = array("d",[X0])l = array("d")while True:X = L[samplecount]cur = np.argmax(np.random.multinomial(1,A[X]))count += 1if isMH:a = (p[cur]*A[cur][X])/(p[X]*A[X][cur])alpha = min(a,1)else:alpha = p[cur]*A[cur][X]u = np.random.uniform(0,1) if u<alpha:samplecount += 1L.append(cur)if count>N:l.append(cur)if len(l)>=Nlmax:breakelse:continueLa = np.frombuffer(L)la = np.frombuffer(l)return La,la def count(q,n):L = array("d")l1 = array("d")l2 = array("d")for e in q:L.append(e)for e in range(n):l1.append(L.count(e))for e in l1:l2.append(e/sum(l1))return l1,l2if __name__ == '__main__': p = np.array([0.9,0.05,0.05])a = mcmc(p)[1]l1 = ['state%d'% x for x in range(len(p))]plt.pie(count(a,len(p))[0],labels=l1,labeldistance=0.3,autopct='%1.2f%%')plt.title("sampling")plt.show()

得到如下：

與圖1相比，很明顯，圖5準(zhǔn)確率提高很多，依照同樣程序，分別更改p，同樣也可以得到對應(yīng)于圖2、3、4的準(zhǔn)確率高的結(jié)果。

結(jié)論

要將狀態(tài)轉(zhuǎn)移矩陣的每一行都設(shè)計成與輸入概率分布相一致，這樣該算法才能得到更好的結(jié)果。
該算法該實現(xiàn)的適用范圍：輸入概率分布為離散的，有限的。
思考：雖然無限的離散型的應(yīng)該也能按照上述的模式計算，但是對于無窮矩陣的計算，本小白還沒研究，等著日后思考；而對于連續(xù)分布，則沒有所謂的概率分布，而替代以概率密度，而且，由于在實現(xiàn)過程中使用了隨機數(shù)生成，所以可以輕易通過變換隨機數(shù)的方式生成的連續(xù)性分布的樣本如果使用mcmc算法未免有點大炮轟蚊子，但若是考慮那些無法顯示表出概率密度的連續(xù)型分布，我實在是還沒想到好的解決辦法，待后續(xù)解決。

總結(jié)

以上是生活随笔為你收集整理的MCMC算法—MH算法的Python实现（一）的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

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