统计学习方法第十九章作业：马尔可夫链蒙特卡罗法、吉布斯抽样算法(书上题目) 代码实现

马尔可夫链蒙特卡罗法

作业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()

本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！