马尔可夫链 Monte Carlo 积分和无限循环
Markov Chain Monte Carlo integration and infinite while loop
我正在使用 Metropolis 和 Barker 的 α 实现马尔可夫链 Monte Carlo 以进行数值积分。我创建了一个名为 MCMCIntegrator()
的 class。 __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
是 class 必须生成的样本大小,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()
我陷入了无限循环,不知道为什么会这样。
虽然我之前没有接触过 Python 也没有对无限循环的直接解释,但代码中存在一些问题:
-
while i != len(sample):
仅当统一变量低于接受概率时,循环才增加值i
if al > u:
这不是 Metropolis-Hastings 的运作方式。 如果统一变量高于接受概率,链的当前值是重复的。但这并不能解释无限循环,因为建议的值最终应该被接受。
如果目标密度是
st.gamma.pdf(x, 1, scale=1.378008857)*np.abs(np.cos(1.10257704))
那么 (i) 第二个常数项 np.abs(np.cos(1.10257704))
的意义何在 (ii) 哪里需要这么多数字?
提案分布
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%。)
我正在使用 Metropolis 和 Barker 的 α 实现马尔可夫链 Monte Carlo 以进行数值积分。我创建了一个名为 MCMCIntegrator()
的 class。 __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
是 class 必须生成的样本大小,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()
我陷入了无限循环,不知道为什么会这样。
虽然我之前没有接触过 Python 也没有对无限循环的直接解释,但代码中存在一些问题:
-
while i != len(sample):
仅当统一变量低于接受概率时,循环才增加值
i
if al > u:
这不是 Metropolis-Hastings 的运作方式。 如果统一变量高于接受概率,链的当前值是重复的。但这并不能解释无限循环,因为建议的值最终应该被接受。 如果目标密度是
st.gamma.pdf(x, 1, scale=1.378008857)*np.abs(np.cos(1.10257704))
那么 (i) 第二个常数项
np.abs(np.cos(1.10257704))
的意义何在 (ii) 哪里需要这么多数字?提案分布
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%。)