Over the past few days I have been trying to understand how Markov Chain Monte Carlo (MCMC) works. In particular I have been trying to understand and implement the Metropolis-Hastings algorithm. So far I think I have an overall understanding of the algorithm but there are a couple of things that are not clear to me yet. I want to use MCMC to fit some models to data. Because of this I will describe my understanding of Metropolis-Hastings algorithm for fitting a straight line $f(x)=ax$ to some observed data $D$:
1) Make an initial guess for $a$. Set this $a$ as our current $a$ ($a_0$). Also add $a$ at the end of the Markov Chain ($C$).
2) Repeat the steps bellow a number of times.
3) Evaluate current likelihood (${\cal L_0}$) given $a_0$ and $D$.
4) Propose a new $a$ ($a_1$) by sampling from a normal distribution with $\mu=a_0$ and $\sigma=stepsize$. For now, $stepsize$ is constant.
5) Evaluate new likelihood (${\cal L_1}$) given $a_1$ and $D$.
6) If ${\cal L_1}$ is bigger than ${\cal L_0}$, accept $a_1$ as the new $a_0$, append it at the end of $C$ and go to step 2.
7) If ${\cal L_1}$ is smaller than ${\cal L_0}$ generate a number ($U$) in range [0,1] from a uniform distribution
8) If $U$ is smaller than the difference between the two likelihoods (${\cal L_1}$ - ${\cal L_0}$), accept $a_1$ as the new $a_0$, append it at the end of $C$ and go to step 2.
9) If $U$ is bigger than the difference between the two likelihoods (${\cal L_1}$ - ${\cal L_0}$), append the $a_0$ at the end of $C$, keep using the same $a_0$, go to step 2.
10) End of Repeat.
11) Remove some elements from the start of $C$ (burn-in phase).
12) Now take the average of the values in $C$. This average is the estimated $a$.
Now I have some questions regarding the above steps:
- How do I construct the likelihood function for $f(x)=ax$ but also for any arbitrary function?
- Is this a correct implementation of Metropolis-Hastings algorithm?
- How the selection of the random number generation method at Step 7 can change the results?
- How is this algorithm going to change if I have multiple model parameters? For example, if I had the model $f(x)=ax+b$.
Notes/Credits: The main structure of the algorithm described above is based on the code from an MPIA Python Workshop.