馬爾可夫鏈蒙特卡羅法

作業19.7

import numpy as np import matplotlib.pyplot as plt from scipy.stats import betaclass MCMC:def __init__(self,scale=0.5):self.ta = np.random.random(1)self.scale = 0.5def update_ta(self):ta_n = np.random.normal(loc=self.ta, scale=self.scale, size=1)[0]a = min(1,beta.pdf(ta_n,1,1)/beta.pdf(self.ta,1,1))u = np.random.random(1)if u <= a :self.ta = ta_ndef fit(self,n,m):self.sample_list = []for i in range(m):self.update_ta()if i > n :self.sample_list.append(self.ta)def plot(self):plt.hist(self.sample_list,bins=50,alpha=0.3)plt.show()def main():mc = MCMC(0.2)mc.fit(5000,10000)print(np.mean([beta.pdf(x,1,1) for x in mc.sample_list])*0.4)mc.plot()if __name__ == '__main__':main()

例題19.10

import numpy as np import matplotlib.pyplot as pltclass MCMC:def __init__(self,p=None):self.p = pself.x1 = np.random.random(1)[0]self.x2 = np.random.random(1)[0]def update_x1(self):self.x1 = np.random.normal(loc=self.p*self.x2, scale=np.sqrt(1-self.p**2), size=1)[0]def update_x2(self):self.x2 = np.random.normal(loc=self.p*self.x1, scale=np.sqrt(1-self.p**2), size=1)[0]def fit(self,n,m):self.sample_list = []self.x1_list = []self.x2_list = []for i in range(m):self.update_x1()self.update_x2()if i > n :self.sample_list.append((self.x1,self.x2))self.x1_list.append(self.x1)self.x2_list.append(self.x2)def plot(self):plt.hist(self.x1_list,bins=50,alpha=0.3)plt.hist(self.x2_list,bins=50,alpha=0.3)plt.hist(np.random.normal(loc=0, scale=np.sqrt(1-self.p**2), size=5000),bins=50,alpha=.3)plt.show()def main():mc = MCMC(0.5)mc.fit(5000,10000)print(mc.sample_list)mc.plot()if __name__ == '__main__':main()

總結

以上是生活随笔為你收集整理的统计学习方法第十九章作业：马尔可夫链蒙特卡罗法、吉布斯抽样算法(书上题目) 代码实现的全部內容，希望文章能夠幫你解決所遇到的問題。

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