马尔可夫链蒙特卡罗积分和无限循环

Question

我正在使用 Metropolis 和 Barker 的 α 实现马尔可夫链蒙特卡罗，用于数值积分。我创建了一个名为MCMCIntegrator() 的类。在__init__() 方法下方，是我正在集成的函数的PDF g(x) 和alpha 方法，实现Metropolis 和Barker α's。

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

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

Answer 1

虽然我之前没有接触过 Python，也没有直接解释无限循环，但代码中存在一些问题：

1.  while i != len(sample):
   

   只有当 Uniform 变量低于接受概率 if al > u: 时，循环才会增加值 i 这不是 Metropolis-Hastings 的操作方式。如果 Uniform 变量高于接受概率，则链的当前值被复制。但这并不能解释无限循环，因为建议的值最终应该被接受。

2. 如果目标密度是

   st.gamma.pdf(x, 1, scale=1.378008857)*np.abs(np.cos(1.10257704))  
   

   那么 (i) 第二个常数项 np.abs(np.cos(1.10257704)) 的意义何在以及 (ii) 哪里需要这么多数字？

3. 提案分布

      new = np.random.normal(loc=old, scale=self.std)
      new = abs(new)
   

   是折叠法线，其密度不是对称的。因此它应该出现在 Metropolis-Hastings 概率中，但由于规模较小，可能影响不大。

这是我对 Python 代码的 R 渲染（已编辑和更正）

self.size = 1e5
self.std = 0.6
self.real_int = 0.06496359

g <- function(x){dgamma(x, shape=1, scale=1.378)}

alpha <- function(a, b, method=1)ifelse(method,
            min(1, r <- g(b) / g(a)), 1 / (1 + 1 / r))

old = 0 
smple = rep(0,self.size * 1.3)
for (i in 1:length(smple)){
            new = abs(old+self.std*rnorm(1))
            al = alpha(old, new, 0)
            old=smple[i]=ifelse(al > runif(1), new, old)
             }

ind = trunc(length(smple)*0.3)
smple = sample[ind:length(smple)]

hist(smple,pro=TRUE,nclass=10*log2(self.size),col="wheat")
curve(g(x),add=TRUE,lwd=2,col="sienna")

清晰再现 Gamma 目标：

没有对非对称提案进行修正。更正将是

q <- function(a, b)
        dnorm(b-a,sd=self.std)+dnorm(-b-a,sd=self.std)

alpha <- function(a, b, method=1){
        return(ifelse(method,
            min(1, r <- g(b) * q(b,a) / g(a) / q(a,b)),
            1 / (1 + 1/r)))}

old = 0 
smple = rep(0,self.size * 1.3)
for (i in 1:length(smple)){
            new = abs(old+self.std*rnorm(1))
            al = alpha(old, new, 3)
            old=smple[i]=ifelse(al > runif(1), new, old)
             }

和上图没有区别。 （Metropolis ratio 的录取率为 85%，而 Baker's 的录取率为 48%。）