Markov Chain蒙特卡洛积分和无限while循环

问题描述 投票:1回答:1

我正在用大都会和巴克的α来实现马尔可夫链蒙特卡罗,以进行数值积分。我创建了一个名为MCMCIntegrator()的类。在__init__()方法的下面,是g(x)我要集成的函数的PDF和alpha方法,用于实现Metropolis和Barkerα。

import numpy as np
import scipy.stats as st

class MCMCIntegrator:

    def __init__(self):
        self.size = 1000
        self.std = 0.6
        self.real_int = 0.06496359
        self.sample = None

    @staticmethod
    def g(x):
        return st.gamma.pdf(x, 1, scale=1.378008857)*np.abs(np.cos(1.10257704))

    def alpha(self, a, b, method):

        if method:
            return min(1, self.g(b) / self.g(a))

        else:
            return self.g(b) / (self.g(a) + self.g(b))

size是该类必须生成的样本的大小,std是正常内核的标准偏差,您将在几秒钟内看到它。 real_int是我们要积分的函数的1到2之间的积分值。我已经用R脚本生成了它。现在,解决问题。

     def _chain(self, method):

        """
            Markov chain heat-up with burn-in

        :param method: Metropolis or Barker alpha
        :return: np.array containing the sample
        """
        old = 0
        sample = np.zeros(int(self.size * 1.3))
        i = 0

        while i != len(sample):
            new = np.random.normal(loc=old, scale=self.std)
            new = abs(new)
            al = self.alpha(old, new, method=method)
            u = np.random.uniform()

            if al > u:
                sample[i] = new
                i += 1
                old = new

        return np.array(sample)

此方法下面是一种integrate()方法,用于计算[1、2]区间中数字的比例:

    def integrate(self, method=None):
        """
            Integration step

        """

        sample = self._chain(method=method)

        # discarding 30% of the sample for the burn-in
        ind = int(len(sample)*0.3)
        sample = sample[ind:]
        setattr(self, "sample", sample)

        sample = [1 if 1 < v < 2 else 0 for v in sample]
        return np.mean(sample)

这是主要功能:

def main():

    print("-- RESULTS --".center(20), end='\n')
    mcmc = MCMCIntegrator()
    print(f"\t{mcmc.integrate()}", end='\n')
    print(f"\t{np.abs(mcmc.integrate() - mcmc.real_int) / mcmc.real_int}")


if __name__ == "__main__":
    main()

我陷入了无限的while循环,不知道为什么会这样。

我正在用大都会和巴克的α来实现马尔可夫链蒙特卡罗,以进行数值积分。我创建了一个名为MCMCIntegrator()的类。 __init __()方法下面是g(x)...

python mcmc
1个回答
2
投票

虽然我以前没有接触过Python,也没有对无限循环的直接解释,但是代码中存在一些问题:

© www.soinside.com 2019 - 2024. All rights reserved.