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

最新推荐文章于 2022-12-25 07:00:00 发布

是算法不是法术

最新推荐文章于 2022-12-25 07:00:00 发布

阅读量807

点赞数

分类专栏：统计学习方法文章标签： python 算法机器学习

本文链接：https://blog.csdn.net/weixin_45839693/article/details/112077048

版权

统计学习方法专栏收录该内容

17 篇文章 2 订阅

订阅专栏

马尔可夫链蒙特卡罗法

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