To best get in touch with the main ideas behind MCMC, the thing to do is to code the father of all of these algorithms, Metropolis-Hastings, on a very simple case.
So here is a set of (independent) measured values ("data") assuming to originate from a Gaussian distribution (the "model") of null mean value but unknown standard deviation 
("the parameter") that we want to infer. Copy the file and then:
- write the likelihood function L or rather chi2=2*log(L) : what is its arguments? draw the function. Following Main.Lik#toc5 give the 68% confidence-level interval.
- can take any value? how can it be descrived in terms of prior knowledge?
- code the Metropolis-Hastings algorithm described in the introduction: for the proposal use a simple Gaussian shot. It is centered on 0 but what should be its standard deviation (that we will call "scaling" in the following to avoid confusion) ? There is no clear answer. Some "rules of thumbs" exist. There are based on studying the acceptance rate (AR) which is the mean number of times the proposed move was accepted. So compute it also in the code. The rule of thumb for gaussians is the following: try to reach an acceptance rate of about 50% in 1D (so here), and about 25% for several dimensions by playing with the scaling .
- the concatenated samples forms the "chain" (do not reject repeated samples which are absolutely necessary to sample really the distribution). To visually check the quality of sampling simply plot the evolution of the samples and possibly their likelihood values.
When you are satisfied just do the histogram of your samples which represents the posterior distribution of the parameter , or rather an estimate of it.
How do you obtain a 68% credible interval? Compare to the chi2 function
back
