data = rng.normal(1.1, 1.0, 40)
def logpost(theta):
return stats.norm.logpdf(theta, 0, 3) + stats.norm.logpdf(data, theta, 1).sum()
logpost(0.0)np.float64(-76.53621846778138)Sampling-Based Approaches to Calculating Marginal Densities
Alan Gelfand and Adrian Smith (1990)
Core Contribution
Bayesian inference often needs marginal posterior quantities
\[ p(\theta_j\mid y)=\int p(\theta\mid y)\,d\theta_{-j}, \]
which can be analytically intractable. Gelfand and Smith helped make simulation-based posterior calculation routine: construct a Markov chain with stationary density \(p(\theta\mid y)\) and estimate posterior expectations by
\[ E[g(\theta)\mid y]\approx \frac{1}{M}\sum_{m=1}^M g(\theta^{(m)}). \]
Minimal Implementation
Define the posterior density \(p(\theta\mid y)\) up to an additive constant.
Generate a Markov chain \(\theta^{(m)}\) with a random-walk Metropolis step.
draws = np.zeros(5000)
theta = 0.0; accept = 0
for i in range(len(draws)):
prop = theta + rng.normal(0, 0.35)
if np.log(rng.uniform()) < logpost(prop) - logpost(theta):
theta = prop; accept += 1
draws[i] = theta
post = draws[800:]
post.mean(), post.std()(np.float64(0.9789782951079137), np.float64(0.1566021053267428))Plot the chain and the posterior marginal estimated from simulation.
fig, ax = plt.subplots(1, 2, figsize=(8, 3.3))
ax[0].plot(draws[:900], lw=0.8)
ax[0].set(title=f"trace, accept={accept/len(draws):.2f}")
ax[1].hist(post, bins=40, density=True, alpha=0.65)
grid = np.linspace(post.min() - 0.3, post.max() + 0.3, 200)
ax[1].plot(grid, stats.gaussian_kde(post)(grid), color="#ffcc66", lw=2.3)
ax[1].set(title="posterior marginal")
plt.show()