Pencil and paper example of fitting normal distribution to data with MCMC

Question

I've been trying to understand Markov Chain Monte Carlo methods for a while and even though I somewhat get the idea, when it comes to me applying MCMC, I'm not sure what I should do. Many times I've gotten the answer "use a package" from professors, but I don't want to use a package, I want to do it myself so I can understand! ;)

I've been trying to search internet for good examples but so far I have found none such example which explicitly shows step by step what is happening under the hood. Many examples which I've looked deal with kinda abstract examples which are difficult to fully grasp. In many of the examples what I'm left with is: "Okay very nice method, now I want to write a MCMC program myself which fits e.g. normal distribution to data". The problem is, I cannot program abstract concepts, I need a concrete example. Some examples I've found get lost in the details of the domain problem and forget to explain the MCMC method itself.

So my question is: Can you provide me with a very simple pencil and paper example? This could be e.g. a problem where you have three data points, say $(4,5), (3,3)$ and $(4,2)$ in the xy-plane and you need to fit this data to normal distribution. Lets also assume that for some reason we need to apply MCMC for this problem (Metropolis-Hastings would be nice). Show me detailed steps of how you do this, taking the derivatives, equating to zero, choosing the proposal distribution, taking a random sample from posterior (if you use Bayesian method) et cetera. Do this like one or two iterations, just to show everybody how the calculations proceed.

I think this example would provide valuable info for beginners in MCMC methods. Thank you!

"I hear and I forget. I see and I remember. I do and I understand."

Confucius

Answer 1

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

1. Specify the probability density function for the prior
2. Load the sample data
3. 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.