马尔可夫链蒙特卡洛积分和无限循环
Markov Chain montecarlo integration and infinite while loop
我正在实施一个马尔可夫链蒙特卡洛与大都会和 barkes alphas 用于数值积分。我创建了一个名为 MCMCIntegrator()
的 class。我已经为它加载了一些属性,其中之一是我们试图集成的函数(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
这个class还有其他方法,size
是class必须生成的样本大小,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()
我陷入了无限 while 循环,我不知道为什么会这样。
两件事...您挂在链方法中,因为 alpha 计算返回 NaN,因为 g() 返回 NaN。看一下我插入到您的代码中的打印语句,运行 它...
小贴士:
使 g() 成为 class 函数,就像 chain
.
在一些测试值上测试 g(),显然有些不对劲
- 不要像
alpha
这样有条件地定义函数。 Wayyy 混淆和错误 prone/tough 进行故障排除。只需传递它需要的 alpha 即可,您还可以将其设为 class 函数 alpha(a, b, method=None)
- 看看您在 `_chain' 函数中更新 i 的位置....您意识到您正在冒着长时间循环过程的风险,因为您只是有条件地更新 i!
- 使用
numpy
阵列,您将面临灾难。您可能在实际数据之后有一堆尾随零,因为您正在覆盖大零列表。您在这里不需要 numpy
数组,只需使用 python 空列表并向其附加新值,零或一个...基于任何内容。
在进行故障排除(或对函数进行单元测试)时添加几个打印语句。试试我在下面添加到你的函数中……这是我用来弄清楚发生了什么的方法
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)
我正在实施一个马尔可夫链蒙特卡洛与大都会和 barkes alphas 用于数值积分。我创建了一个名为 MCMCIntegrator()
的 class。我已经为它加载了一些属性,其中之一是我们试图集成的函数(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
这个class还有其他方法,size
是class必须生成的样本大小,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()
我陷入了无限 while 循环,我不知道为什么会这样。
两件事...您挂在链方法中,因为 alpha 计算返回 NaN,因为 g() 返回 NaN。看一下我插入到您的代码中的打印语句,运行 它...
小贴士:
使 g() 成为 class 函数,就像
chain
.在一些测试值上测试 g(),显然有些不对劲
- 不要像
alpha
这样有条件地定义函数。 Wayyy 混淆和错误 prone/tough 进行故障排除。只需传递它需要的 alpha 即可,您还可以将其设为 class 函数alpha(a, b, method=None)
- 看看您在 `_chain' 函数中更新 i 的位置....您意识到您正在冒着长时间循环过程的风险,因为您只是有条件地更新 i!
- 使用
numpy
阵列,您将面临灾难。您可能在实际数据之后有一堆尾随零,因为您正在覆盖大零列表。您在这里不需要numpy
数组,只需使用 python 空列表并向其附加新值,零或一个...基于任何内容。
在进行故障排除(或对函数进行单元测试)时添加几个打印语句。试试我在下面添加到你的函数中……这是我用来弄清楚发生了什么的方法
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)