In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
%matplotlib inline
plt.rcParams['figure.figsize'] = (5.0, 5.0)
# the model parameters
a = np.pi
b = 1.6818
# my arbitrary constants
mu_x = np.exp(1.0) # see definitions above
tau_x = 1.0
s = 1.0
N = 50 # number of data points
# get some x's and y's
x = mu_x + tau_x*np.random.randn(N)
y = a + b*x + s*np.random.randn(N)
plt.plot(x, y, 'o');
In [ ]:
def lnPost(params, x, y):
# This is written for clarity rather than numerical efficiency. Feel free to tweak it.
a = params[0]
b = params[1]
lnp = 0.0
# Using informative priors to achieve faster convergence is cheating in this exercise!
# But this is where you would add them.
lnp += -0.5*np.sum((a+b*x - y)**2)
return lnp
In [ ]:
class ExactPosterior:
def __init__(self, x, y, a0, b0):
X = np.matrix(np.vstack([np.ones(len(x)), x]).T)
Y = np.matrix(y).T
self.invcov = X.T * X
self.covariance = np.linalg.inv(self.invcov)
self.mean = self.covariance * X.T * Y
self.a_array = np.arange(0.0, 6.0, 0.02)
self.b_array = np.arange(0.0, 3.25, 0.02)
self.P_of_a = np.array([self.marg_a(a) for a in self.a_array])
self.P_of_b = np.array([self.marg_b(b) for b in self.b_array])
self.P_of_ab = np.array([[self.lnpost(a,b) for a in self.a_array] for b in self.b_array])
self.P_of_ab = np.exp(self.P_of_ab)
self.renorm = 1.0/np.sum(self.P_of_ab)
self.P_of_ab = self.P_of_ab * self.renorm
self.levels = scipy.stats.chi2.cdf(np.arange(1,4)**2, 1) # confidence levels corresponding to contours below
self.contourLevels = self.renorm*np.exp(self.lnpost(a0,b0)-0.5*scipy.stats.chi2.ppf(self.levels, 2))
def lnpost(self, a, b): # the 2D posterior
z = self.mean - np.matrix([[a],[b]])
return -0.5 * (z.T * self.invcov * z)[0,0]
def marg_a(self, a): # marginal posterior of a
return scipy.stats.norm.pdf(a, self.mean[0,0], np.sqrt(self.covariance[0,0]))
def marg_b(self, b): # marginal posterior of b
return scipy.stats.norm.pdf(b, self.mean[1,0], np.sqrt(self.covariance[1,1]))
exact = ExactPosterior(x, y, a, b)
In [ ]:
plt.plot(exact.a_array, exact.P_of_a);
In [ ]:
plt.plot(exact.b_array, exact.P_of_b);
In [ ]:
plt.contour(exact.a_array, exact.b_array, exact.P_of_ab, colors='blue', levels=exact.contourLevels);
plt.plot(a, b, 'o', color='red');
Normally, you should always use quantitative tests of convergence in addition to visual inspection, as you saw on Tuesday. For this class (only), let's save some time by relying only on visual impressions and comparison to the exact posterior.
(see the snippets farther down)
In [ ]:
Nsamples = # fill in a number
samples = np.zeros((Nsamples, 2))
# put any more global definitions here
for i in range(Nsamples):
a_try, b_try = proposal() # propose new parameter value(s)
lnp_try = lnPost([a_try,b_try], x, y) # calculate posterior density for the proposal
if we_accept_this_proposal(lnp_try, lnp_current):
# do something
else:
# do something else
In [ ]:
plt.rcParams['figure.figsize'] = (12.0, 3.0)
plt.plot(samples[:,0]);
plt.plot(samples[:,1]);
In [ ]:
plt.rcParams['figure.figsize'] = (5.0, 5.0)
plt.plot(samples[:,0], samples[:,1]);
In [ ]:
plt.rcParams['figure.figsize'] = (5.0, 5.0)
plt.hist(samples[:,0], 20, normed=True, color='cyan');
plt.plot(exact.a_array, exact.P_of_a, color='red');
In [ ]:
plt.rcParams['figure.figsize'] = (5.0, 5.0)
plt.hist(samples[:,1], 20, normed=True, color='cyan');
plt.plot(exact.b_array, exact.P_of_b, color='red');
In [ ]:
# If you know how to easily overlay the 2D sample and theoretical confidence regions, by all means do so.