python - 如何使用 Metropolis Hastings 从 Gamma(26,6) 分布中采样

Question

我正在研究 MCMC，并决定测试我对共轭贝叶斯问题的 Metropolis Hastings 算法的理解。给定泊松似然和 Gamma 先验，我试图绘制后验分布。给出数据如下，y = (4, 4, 5, 8, 3)，gamma先验参数为a = b = 1。

因为这是一个共轭模型，所以我想出了真正的后验概率，Ga(a + n*y_bar, b + n) = Ga(25, 6)。这是它的代码-

##### The real distribution ####
import numpy as np
import scipy.stats as st
import matplotlib.pyplot as plt
data = np.array([4, 4, 5, 8, 3])
n = len(data)
#Gamma(shape = 25, rate = 6)
lik = st.gamma.pdf(a = 25, scale = 1/6, x = np.linspace(0, 10, num=100))
plt.plot(lik)
plt.show()

这是我的 Metropolis Hastings 代码。但是，我无法使用它获得真正的后验。我从这里借鉴了很多想法 - https://people.duke.edu/~ccc14/sta-663/MCMC.html

def target(lambda_val, log_val = True):
    if log_val == False:
        return ((lambda_val)**(n*np.mean(data)+a+1))*np.exp(-lambda_val*(b+n))
    else:
        return np.exp((n*np.mean(data)+a+1)*np.log(lambda_val) - lambda_val*(n+b))

### Using MCMC to sample from it #####
niters = 1000
theta = 40

a = 1
b = 1
samples = np.zeros(niters+1)
samples[0] = 5
naccept = 0
for i in range(niters):
    theta_p = theta + st.norm(0,3).rvs()
    rho= min(1, target(theta_p)/target(theta))
    u = np.random.uniform()
    if u<rho:
        naccept += 1
        theta = theta_p
    samples[i+1] = theta
nmcmc = len(samples)//2
print("Efficiency = ", naccept/niters)

plt.hist(samples[nmcmc:], 40)
plt.show()

如果我遗漏了任何细节，请告诉我。

Answer 1

可能性没有正确定义。具体来说，需要通过返回 0(或 -Inf for logspace)来明确排除支持(非正)之外的提案。

考虑到实际位置 (~4)，步长 (sd=3) 也相当大。这导致经常提出负数。

最后，为什么初始化是theta=40？这可能是由于误解了绘图，它没有使用正确的 x 轴，而只使用了一个索引:

使用有效的 x 轴绘图:

xs = np.linspace(0, 10, 100)
lik = st.gamma.pdf(a=25, scale=1/6, x=xs)
plt.plot(xs, lik)
plt.show()

有人认为更合适的初始化可能是 4。