The MH code is basically directly copied from code made available by Chris Fonnesbeck. I'm just using it to learn. I adapated the example he used to work with my own model to infer the mean and precision of some data.
Resources used:
I am going to do this first with a very simple situation, estimating the mean and variance of some data.
Priors
$\mu \sim Normal(0,\sigma=1000)$
$\tau \sim Gamma(0.01, 0.01)$
Likelihood
$data \sim Normal(\mu, tau)$
import numpy as np
true_mean = 2
true_sigma = 1.5
true_tau = 1 / true_sigma**2
print(true_tau)
N = 15
data = np.random.normal(true_mean, true_sigma, N)
0.444444444444
print(data)
[ 0.84416686 1.20806974 3.95816697 2.66058809 3.59166034 1.87905828 3.04081119 2.70481818 1.77847004 1.5228081 4.13710834 1.68940245 -1.31600811 2.31972634 1.79109786]
%matplotlib inline
import matplotlib.pyplot as plt
plt.rc('font', family='serif')
plt.rc('font', serif='Times New Roman')
plt.hist(data);
import scipy.stats as sp
def calc_posterior(mu, tau, x=data):
'''
This is basically the model. This function calculates the log
posterior of a given set of parameters.
'''
# prior for mu
logp = sp.distributions.norm.logpdf(mu, 0, 1000)
# prior for tau
logp += sp.distributions.gamma.logpdf(tau, 0.01, 0.01)
# calculate sigma from tau, needed for likelihood
sigma = (1.0/tau)**0.5
# data likeligood
logp += sum(sp.distributions.norm.logpdf(x, mu, sigma))
return logp
For fun, to test this, calculate some log posteriors for some paramer combinations.
print(calc_posterior(mu=-2, tau=1./1**2, x=data))
print(calc_posterior(mu=-1, tau=1./0.5**2, x=data))
print(calc_posterior(mu=1, tau=1./0.1**2)) # don't HAVE to supply data as it's a KWARG
print(calc_posterior(mu=2, tau=1./2**2, x=data))
-167.652197482 -365.780793581 -2349.42645371 -38.7401449306
So now we have some working code that calculates the log posterior probability of some parameter combination with the data
rnorm = np.random.normal
runif = np.random.rand
def metropolis(n_iterations, initial_values, prop_var=1):
n_params = len(initial_values)
# Initial proposal standard deviations
prop_sd = [prop_var]*n_params
# Initialize trace for parameters
trace = np.empty((n_iterations+1, n_params))
# Set initial values
trace[0] = initial_values
# Calculate joint posterior for initial values
current_log_prob = calc_posterior(*trace[0], x=data)
# Initialize acceptance counts
accepted = [0]*n_params
for i in range(n_iterations):
if not i%1000: print 'Iteration', i
# Grab current parameter values
current_params = trace[i]
for j in range(n_params):
# Get current value for parameter j
p = trace[i].copy()
# Propose new value
if j==2:
# Ensure tau is positive
theta = np.exp(rnorm(np.log(current_params[j]), prop_sd[j]))
else:
theta = rnorm(current_params[j], prop_sd[j])
# Insert new value
p[j] = theta
# Calculate log posterior with proposed value
proposed_log_prob = calc_posterior(*p, x=data)
# Log-acceptance rate
alpha = proposed_log_prob - current_log_prob
# Sample a uniform random variate
u = runif()
# Test proposed value
if np.log(u) < alpha:
# Accept
trace[i+1,j] = theta
current_log_prob = proposed_log_prob
accepted[j] += 1
else:
# Reject
trace[i+1,j] = trace[i,j]
return trace, accepted
n_iter = 5000
initial_params = (-3,0.2)
trace, acc = metropolis(n_iter, initial_params, prop_var=0.1)
Iteration 0 Iteration 1000 Iteration 2000 Iteration 3000 Iteration 4000
plt.plot(trace)
plt.axhline(y=true_tau)
plt.axhline(y=true_mean)
<matplotlib.lines.Line2D at 0x10ba63450>
plt.figure(figsize=(8,8))
plt.plot(trace[:,0], trace[:,1])
plt.xlabel('$\mu$')
plt.ylabel('tau')
<matplotlib.text.Text at 0x10ba805d0>
plt.figure(1, figsize=(12,6))
x=np.linspace(-15,15,500)
for i in range(3000):
mu, tau = trace[i,:]
sigma = (1.0/tau)**0.5
plt.plot(x, sp.norm.pdf(x,mu,sigma), alpha=0.01, color='r')
plt.plot(data, np.zeros(data.shape)
,'+'
, c='k'
, markersize=20)
#plt.savefig('MH_example.png', bbox_inches='tight')
[<matplotlib.lines.Line2D at 0x10ba9c650>]
trace.shape
(5001, 2)
trace[:,0]
array([-3. , -3.02442006, -2.99259942, ..., 1.57084789, 1.57084789, 1.57084789])
mu, tau = trace[1,:]
mu
-2.9287384801694802
tau
0.20000000000000001
Create a range of initial starting values
muv, tauv = np.meshgrid(np.linspace(-20, 20, 4), np.linspace(-20, 20, 4))
def run_chains(muv, tauv, n_iter):
# number of chains decided by number of initial parameters provided
n_chains = muv.size
fig = plt.figure(figsize=(8,8))
for chain in range(n_chains):
print 'chain: ' , chain
intial_params=(muv.flat[chain], tauv.flat[chain])
print intial_params
trace, acc = metropolis(n_iter, initial_params)
plt.plot(trace[:,0], trace[:,1], alpha=0.5)
plt.xlabel('$\mu$')
plt.ylabel('tau')
plt.axhline(y=true_tau, linewidth=3, c='r', ls='--')
plt.axvline(x=true_mean, linewidth=3, c='r', ls='--')
run_chains(muv, tauv, n_iter=500)
chain: 0 (-20.0, -20.0) Iteration 0 chain: 1 (-6.6666666666666661, -20.0) Iteration 0 chain: 2 (6.6666666666666679, -20.0) Iteration 0 chain: 3 (20.0, -20.0) Iteration 0 chain: 4 (-20.0, -6.6666666666666661) Iteration 0 chain: 5 (-6.6666666666666661, -6.6666666666666661) Iteration 0 chain: 6 (6.6666666666666679, -6.6666666666666661) Iteration 0 chain: 7 (20.0, -6.6666666666666661) Iteration 0 chain: 8 (-20.0, 6.6666666666666679) Iteration 0 chain: 9 (-6.6666666666666661, 6.6666666666666679) Iteration 0 chain: 10 (6.6666666666666679, 6.6666666666666679) Iteration 0 chain: 11 (20.0, 6.6666666666666679) Iteration 0 chain: 12 (-20.0, 20.0) Iteration 0 chain: 13 (-6.6666666666666661, 20.0) Iteration 0 chain: 14 (6.6666666666666679, 20.0) Iteration 0 chain: 15 (20.0, 20.0) Iteration 0