统计学习方法第十九章作业:马尔可夫链蒙特卡罗法、吉布斯抽样算法(书上题目) 代码实现
作者:互联网
马尔可夫链蒙特卡罗法
作业19.7
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import beta
class MCMC:
def __init__(self,scale=0.5):
self.ta = np.random.random(1)
self.scale = 0.5
def 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_n
def 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 plt
class MCMC:
def __init__(self,p=None):
self.p = p
self.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()
标签:__,self,random,list,第十九章,x2,马尔可夫,蒙特卡罗,np 来源: https://blog.csdn.net/weixin_45839693/article/details/112077048