Markov Chain Monte Carlo

6.1 Metropolis-Hastings Algorithm

python
Run Inline
import numpy as np
# Sample from N(5, 1) using Metropolis-Hastings
target_mean, target_std = 5, 1
n_samples = 10000
samples = np.zeros(n_samples)
current = 0.0
for i in range(n_samples):
    proposal = current + np.random.randn() * 2
    # acceptance probability
    log_ratio = ( -0.5*(proposal-target_mean)**2/target_std**2
                  + 0.5*(current-target_mean)**2/target_std**2 )
    if np.log(np.random.uniform()) < log_ratio:
        current = proposal
    samples[i] = current
print(f"Mean: {samples.mean():.2f}, Std: {samples.std():.2f}")

6.2 Gibbs Sampling