我正在用大都会和巴克的α来实现马尔可夫链蒙特卡罗,以进行数值积分。我创建了一个名为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,也没有对无限循环的直接解释,但是代码中存在一些问题: