为了举例,我们假设我们不知道如何计算后验,因此我们将 Metropolis-Hastings 算法实现到 Python 代码中以获得近似答案。 我们将借助 SciPy 统计函数来完成此操作:
def post(θ, Y, α=1, β=1):
if 0 <= θ <= 1:
prior = stats.beta(α, β).pdf(θ)
like = stats.bernoulli(θ).pmf(Y).prod()
prob = like * prior
else:
prob = -np.inf
return prob
Y = stats.bernoulli(0.7).rvs(20)
n_iters = 1000
can_sd = 0.05
α = β = 1
θ = 0.5
trace = {"θ":np.zeros(n_iters)}
p2 = post(θ, Y, α, β)
for iter in range(n_iters):
θ_can = stats.norm(θ, can_sd).rvs(1)
p1 = post(θ_can, Y, α, β)
pa = p1 / p2
if pa > stats.uniform(0, 1).rvs(1):
θ = θ_can
p2 = p1
trace["θ"][iter] = θ