Metropolis-Hastings 接受-拒绝实施
Metropolis-Hastings accept-reject implementation
我一直在阅读 Metropolis-Hastings (MH) 算法。从理论上讲,我了解该算法的工作原理。现在,我正在尝试使用 python.
实现 MH 算法
我遇到了以下 notebook。它完全适合我的问题,因为我想通过一条直线来拟合我的数据,同时考虑到我的数据的测量误差。我将粘贴我发现难以理解的代码:
# initial m, b
m,b = 2, 0
# step sizes
mstep, bstep = 0.1, 10.
# how many steps?
nsteps = 10000
chain = []
probs = []
naccept = 0
print 'Running MH for', nsteps, 'steps'
# First point:
L_old = straight_line_log_likelihood(x, y, sigmay, m, b)
p_old = straight_line_log_prior(m, b)
prob_old = np.exp(L_old + p_old)
for i in range(nsteps):
# step
mnew = m + np.random.normal() * mstep
bnew = b + np.random.normal() * bstep
# evaluate probabilities
# prob_new = straight_line_posterior(x, y, sigmay, mnew, bnew)
L_new = straight_line_log_likelihood(x, y, sigmay, mnew, bnew)
p_new = straight_line_log_prior(mnew, bnew)
prob_new = np.exp(L_new + p_new)
if (prob_new / prob_old > np.random.uniform()):
# accept
m = mnew
b = bnew
L_old = L_new
p_old = p_new
prob_old = prob_new
naccept += 1
else:
# Stay where we are; m,b stay the same, and we append them
# to the chain below.
pass
chain.append((b,m))
probs.append((L_old,p_old))
print 'Acceptance fraction:', naccept/float(nsteps)
代码简单易懂,但我对MH的实现方式有一定的理解。
我的问题在chain.append
(倒数第三行)。作者附上 m
和 b
是否被接受或拒绝。为什么?他不应该 append
只有接受的分数吗?
快速阅读算法说明:当一个候选人被拒绝时,它仍然算作一个步骤,但值与旧步骤相同。 IE。 b, m
以任何一种方式附加,但只有在候选人被接受的情况下,它们才会 更新 (到 bnew, mnew
)。
以下 R 代码演示了为什么捕获被拒绝的案例很重要:
# 20 samples from 0 or 1. 1 has an 80% probability of being chosen.
the.population <- sample(c(0,1), 20, replace = TRUE, prob=c(0.2, 0.8))
# Create a new sample that only catches changes
the.sample <- c(the.population[1])
# Loop though the.population,
# but only copy the.population to the.sample if the value changes
for( i in 2:length(the.population))
{
if(the.population[i] != the.population[i-1])
the.sample <- append(the.sample, the.population[i])
}
这段代码运行时,the.population
得到20个值,例如:
0 1 1 1 1 1 1 1 1 0 1 1 1 1 0 0 1 1 1 1
此总体中 1 的概率为 16/20 或 0.8。正是我们预期的概率...
另一方面,仅记录更改的样本如下所示:
0 1 0 1 0 1
样本中出现 1 的概率为 3/6 或 0.5。
我们正在尝试构建分布,拒绝新值意味着旧值比新值 更 的可能性。这需要被捕获,以便我们的分布是正确的。
我一直在阅读 Metropolis-Hastings (MH) 算法。从理论上讲,我了解该算法的工作原理。现在,我正在尝试使用 python.
实现 MH 算法我遇到了以下 notebook。它完全适合我的问题,因为我想通过一条直线来拟合我的数据,同时考虑到我的数据的测量误差。我将粘贴我发现难以理解的代码:
# initial m, b
m,b = 2, 0
# step sizes
mstep, bstep = 0.1, 10.
# how many steps?
nsteps = 10000
chain = []
probs = []
naccept = 0
print 'Running MH for', nsteps, 'steps'
# First point:
L_old = straight_line_log_likelihood(x, y, sigmay, m, b)
p_old = straight_line_log_prior(m, b)
prob_old = np.exp(L_old + p_old)
for i in range(nsteps):
# step
mnew = m + np.random.normal() * mstep
bnew = b + np.random.normal() * bstep
# evaluate probabilities
# prob_new = straight_line_posterior(x, y, sigmay, mnew, bnew)
L_new = straight_line_log_likelihood(x, y, sigmay, mnew, bnew)
p_new = straight_line_log_prior(mnew, bnew)
prob_new = np.exp(L_new + p_new)
if (prob_new / prob_old > np.random.uniform()):
# accept
m = mnew
b = bnew
L_old = L_new
p_old = p_new
prob_old = prob_new
naccept += 1
else:
# Stay where we are; m,b stay the same, and we append them
# to the chain below.
pass
chain.append((b,m))
probs.append((L_old,p_old))
print 'Acceptance fraction:', naccept/float(nsteps)
代码简单易懂,但我对MH的实现方式有一定的理解。
我的问题在chain.append
(倒数第三行)。作者附上 m
和 b
是否被接受或拒绝。为什么?他不应该 append
只有接受的分数吗?
快速阅读算法说明:当一个候选人被拒绝时,它仍然算作一个步骤,但值与旧步骤相同。 IE。 b, m
以任何一种方式附加,但只有在候选人被接受的情况下,它们才会 更新 (到 bnew, mnew
)。
以下 R 代码演示了为什么捕获被拒绝的案例很重要:
# 20 samples from 0 or 1. 1 has an 80% probability of being chosen.
the.population <- sample(c(0,1), 20, replace = TRUE, prob=c(0.2, 0.8))
# Create a new sample that only catches changes
the.sample <- c(the.population[1])
# Loop though the.population,
# but only copy the.population to the.sample if the value changes
for( i in 2:length(the.population))
{
if(the.population[i] != the.population[i-1])
the.sample <- append(the.sample, the.population[i])
}
这段代码运行时,the.population
得到20个值,例如:
0 1 1 1 1 1 1 1 1 0 1 1 1 1 0 0 1 1 1 1
此总体中 1 的概率为 16/20 或 0.8。正是我们预期的概率...
另一方面,仅记录更改的样本如下所示:
0 1 0 1 0 1
样本中出现 1 的概率为 3/6 或 0.5。
我们正在尝试构建分布,拒绝新值意味着旧值比新值 更 的可能性。这需要被捕获,以便我们的分布是正确的。