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

Question

我正在用大都会和巴克的α来实现马尔可夫链蒙特卡罗，以进行数值积分。我创建了一个名为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='
')
    mcmc = MCMCIntegrator()
    print(f"	{mcmc.integrate()}", end='
')
    print(f"	{np.abs(mcmc.integrate() - mcmc.real_int) / mcmc.real_int}")


if __name__ == "__main__":
    main()

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

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