马尔科夫链的蒙太卡洛积分和无限while循环

Question

世界！

我正在用metropolis和barkes alphas来实现一个Markov Chain Montecarlo的数值积分。我创建了一个叫做 MCMCIntegrator(). 我已经用一些属性加载了它，其中一个是我们试图整合的函数（一个lambda）的pdf，叫做 g.

import numpy as np
import scipy.stats as st


class MCMCIntegrator:

    def __init__(self):

        self.g = lambda x: st.gamma.pdf(x, 0, 1, scale=1 / 1.23452676)*np.abs(np.cos(1.123454156))
        self.size = 10000
        self.std = 0.6
        self.real_int = 0.06496359

在这个类中还有其他方法，即 size 是该类必须生成的样本大小。std 是正态核的标准差，你将在几秒钟内看到。你将在几秒钟内看到。real_int 是我们要积分的函数从1到2的积分值。我已经用一个R脚本生成了它。现在，说说问题。

 def _chain(self, method=None):

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

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

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

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

        while i != len(sample):
            new = st.norm(loc=old, scale=self.std).rvs()
            new = abs(new)
            al = alpha(old, new)
            u = st.uniform.rvs()

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

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

我卡在了一个无穷大的循环里 我不知道为什么会这样 谁能帮帮我？

Answer 1

有几件事... 你在链式方法中被挂住了，因为α计算正在返回NaN，因为g()正在返回NaN。 看看我在你的代码中插入的打印语句，然后运行它... ...

提示

1. 把g()变成一个类函数，就像 chain.

2. 在一些测试值上测试g()，明显有问题。

3. 不要有条件地定义一个函数，如 alpha. Wayyy混淆和错误pronetough故障排除。 只要把它需要的东西传给alpha，你也可以把它变成一个类函数。alpha(a, b, method=None)
4. 看看你在"_chain "函数中更新i的地方...... 你意识到你正在冒着一个漫长的循环过程的风险，因为你只是有条件地更新i。
5. 你在使用你的"_chain "函数的时候，你意识到你正在冒着漫长的循环过程的风险，因为你只是有条件地更新了i！你正在为你的 numpy 数组。 你可能会在你的实际数据后面有一堆尾部的零，因为你过度写入了大的零列表。 你不需要 numpy 在这里，只需使用一个 python 空列表，然后根据任何情况向其添加新的值，可以是 0 或 1...。

当你进行故障排除（或单元测试你的函数）时，添加几个打印语句。 试试我在下面的函数中添加的内容......我就是用它来找出问题所在的。

def _chain(self, method=None, verbose=True):

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

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

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

    else:
        def alpha(a, b): 
            if verbose: print(f'g(a): {self.g(a)}, g(b): {self.g(b)}')
            return self.g(b) / (self.g(a) + self.g(b))

    while i != len(sample):
        new = st.norm(loc=old, scale=self.std).rvs()
        new = abs(new)
        al = alpha(old, new)
        u = st.uniform.rvs()
        if verbose: print(f'old: {old:.3f} new: {new:.3f} alpha: {al:.3f} u: {u:.3f}')
        if al > u:
            sample[i] = new
            old = new
            i += 1              # do you really want to conditionally update i?
        sys.exit(-1)            # to breakout of infinite loop...


    return np.array(sample)