While not strictly "pen and paper" I imagine you are planning to code this up in some language, so I am going to write out an example in python code, but hopefully written in such a way that I sacrifice efficiency and cleanliness of the code in exchange for something that is easy to follow and translate into a language of your choosing - this is not good code but I hope it is readable code!
So first basic imports, specifying the prior distributions and their hyperparameters as well as loading the data and specifying the log-likelihood. In this code block we
- Specify the probability density function for the prior
- Load the sample data
- Specify the model loglikelihood
import numpy as np
from scipy.stats import multivariate_normal, wishart
#######
# Change this section to use your own priors
# with suitable hyperparameters
def mu_prior_logpdf(m):
m0 = np.zeros(2)
C0 = np.diag(np.ones(2))
return multivariate_normal.logpdf(m, mean=m0, cov=C0)
def S_prior_logpdf(S):
df0 = 2
V0 = np.diag(np.ones(2))
return wishart.logpdf(S, df=df0, scale=V0)
#######
# Your sample data
nData = 3
X = np.array([[4, 5], [3, 3], [4, 2] ])
###
# Your specified model
#
# X_i ~ N(m,S), with X_i i.i.d
def loglik(X, m, S):
ll = 0.
for n in range(nData):
ll += multivariate_normal.logpdf( X[n, ], mean=m, cov=S )
return ll
Step 2: Specify the proposal distribution
There are multiple options at this point with regards to what MCMC algorithm to choose, in particular there are good reasons to use Gibbs sampling - however because you want a fairly general method we will demonstrate a Metropolis-Hastings MCMC algorithm which requires no special knowledge of the conditional distributions.
Recall for MH-MCMC we need to specify a proposal distribution $q(\theta^* | \theta^{(n)} )$ which gives the conditional distribution of a new sample $\theta^*$ conditional on the current value of $\theta^{(n)}$.
Now some pseudo-code for the proposal distributions, again should be clear how to modify all of this, in this particular example the proposal $q(\mu^* | \mu^{(n)})$ is given by a random walk centered on the current value. While the proposal for the covariance is just a basic independent sampler $q(\Sigma^* | \Sigma^{(n)}) = q(\Sigma^*)$, this is likely to display poor convergence and sampling properties but again this is just to demonstrate the basic set up
# Proposal distributions:
# - return a random variable
# - able to evaluate the pdf of q(z' | z)
# mean proposal is simple random walk proposal
# centered at the current value so need to specify the
# std. dev of the step size
scale = 0.1
mpropCov = np.array([[ scale*scale, 0.], [0., scale*scale]])
def mu_proposal_rvs(mcur, size=1):
d = np.random.normal(size=2, scale=scale) # pair of ind. mean 0 steps
return mcur + d
def mu_proposal_logpdf(mnew, mcur):
return multivariate_normal.logpdf(mnew, mean=mcur, cov=mpropCov)
# proposal for the covariance, just going to put a simple independent
# sampler here, so change this for something more suitable
def S_proposal_rvs(scur, size=1):
return wishart.rvs(df=2, scale=np.diag(np.ones(2)))
def S_proposal_logpdf(Snew, Scur):
return wishart.logpdf(Snew, df=2, scale=np.diag(np.ones(2)))
Step 3: Single MCMC step
And now a single MH step for your data can proceed as follows; draw a random variable from the proposal and then accept with probability
$$
A(\theta^*, \theta^{(n)}) = \min \left(1, \frac{p(\mathbf{X}|\theta^*)\pi(\theta^{*}) q(\theta^{(n)}|\theta^*)}{p(\mathbf{X}|\theta^{(n)})\pi(\theta^{(n)})q(\theta^* | \theta^{(n)})} \right)
$$
def MHstep(mcur, Scur):
# propose a new value of the mean parameter
mnew = mu_proposal_rvs(mcur)
# evaluate the acceptance ratio
# note that since we are not updating S at this point the prior
# for S and the proposal of S will cancel in the numerator and
# denominator of A so we don't need to evaluate it here
logAnum = loglik(X, mnew, Scur) + mu_prior_logpdf(mnew) + mu_proposal_logpdf(mcur, mnew)
logAden = loglik(X, mcur, Scur) + mu_prior_logpdf(mcur) + mu_proposal_logpdf(mnew, mcur)
A = np.exp( logAnum - logAden )
# Now we accept the proposed value of the new mean
# with probability A
U = np.random.uniform()
if U <= A:
# Accept the proposed update
mcur = mnew
else:
# do nothing
pass
# Now repeat the steps to update the covariance term with the
# obvious changes
#
# ...
# TO DO!
# ...
return mcur, Scur
And repeat...
So to get a sample of size $N$ from the posterior all you now need to do is supply some initial conditions and start sampling, you can also modify in the obvious way to add a burn in and keep only a subset of the posterior samples etc.
# Initialise some parameters for the chain
mcur = np.array([0., 0.])
Scur = np.array([[1., 0.], [0., 1.]])
meanRVs = []
covRVs = []
nt = 0; N = 10
while nt < N:
mcur, Scur = MHstep(mcur, Scur)
meanRVs.append( mcur )
covRVs.append( Scur )
nt += 1
You can now modify the code above with some print ... statements, or even add user input/output to slow each step down and see exactly what is happening.
