1
$\begingroup$

I am trying to set up a Metropolis-Hastings algorithm in Matlab in order to estimate the parameters ${\theta}$ (it is a vector of 5 elements) to fit a curve to a set of data $D={X_i,Y_i,\delta_i}$. $X$ and $Y$ are the data to relate, $\delta$ indicates if the data is censored or not.

I can estimate the log-likelihood function $log(L(\theta|D))=\delta_i log[f_y(y_i,x_i,\theta)]+(1-\delta_i) log[1-F_w(y_i,x_i,\theta)] $

Moreover, I don't have prior information on the distribution of the parameters. The only things that I know are:

$0<\theta_1<100$; $-5<\theta_2<0$; $0<\theta_3<1.5$; $1<\theta_4<6$; $0<\theta_5<1.5$.

For these reasons I pick a random estimation of $\theta$ in this boundaries, let's call it $\theta_0$. Moreover, I can define the candidate distribution $q_i$ (for each parameter) as a uniform distribution with lower and upper bound as specified. Supposing that there is no statistical correlation, the $q$ is given by the product of $q_i$.

Starting from $\theta_0$, i calculate the log-likelihood $log(L(\theta_0|D))$ and $q_0$, then I perturbate the $\theta_{0_{i}}$ using a normal distribution $N(\theta_{0_{i}};j_i)$ where $j_i=\theta_{0_{i}}*0.1$; what I find is a new vector $\theta_1$ for which I calculate the log-likelihood and the candidate distribution $q$ as specified before.

$\theta_1$ will be accepted as it is if the ratio $r=[log(L(\theta_1|D))*q_1]/[log(L(\theta_0|D))*q_0]$ is bigger than unity or if it bigger than $U(0,1)$ otherwise it will be set equal to $\theta_0$ and a new perturbation will be done to find $\theta_2$ ... and so on.

If this algorithm is right, there is a problem. my log-likelihood is negative then if $\theta_1$ has a higher likelihood, the value of $r$ will be smaller than unity. if I invert the ratio, since there is the risk that $q_1=0$, the ratio can even be not defined.

The approach that I am trying to replicate is from an article in which they define the algorithm as follows:

1: set

  1. set an initial value for the chain: $\theta_c=\theta_0$ and choose $j$
  2. compute a=loglikelihood($\theta_c$)+logprior($\theta_c$)
  3. draw $\theta_p$ from $N(\theta_c;j_i)$
  4. compute b==loglikelihood($\theta_p$)+logprior($\theta_p$)
  5. let H=min(1,exp(b-a)) and draw r from U(0,1)
  6. if H>r then
  7. $\theta_c=\theta_p$
  8. a=b
  9. repeat steps from 3 to 8 until 100000 posteriors are sampled

Is it the same algorithm? Why did they use the exp(b-a) instead then the ratio that is usually used in books?

$\endgroup$
  • $\begingroup$ Because they used log's: mathwords.com/l/logarithm_rules.htm $\endgroup$ – Tim Jul 12 '16 at 13:48
  • $\begingroup$ thank you! So, if I am not wrong, the logprior that they define is the log(q), right? $\endgroup$ – D.Leo Jul 12 '16 at 14:04
  • 1
    $\begingroup$ I don't know what do you mean... But candidate distribution is not the same as prior distribution. $\endgroup$ – Tim Jul 12 '16 at 14:08
  • $\begingroup$ OK, thank you again. the prior distribution is the Uniform distribution that I set for each parameter of $\theta$, right? Did the authors of the algorithm that I reported use any candidate distribution? $\endgroup$ – D.Leo Jul 12 '16 at 14:13
2
$\begingroup$

Since a=loglikelihood($\theta_c$)+logprior($\theta_c$) and b==loglikelihood($\theta_p$)+logprior($\theta_p$), it follows that

$$\exp(b-a) = \dfrac{likelihood(\theta_p)\times prior(\theta_p)}{likelihood(\theta_c)\times prior(\theta_c)},$$

which is the "usual ratio" you refer to. The instrumental distribution used by the authors is normal (this is where you sample the $\theta$s from).

Moreover, you don't need to come up with intricate instrumental distributions, just define your log posterior as:

if($0<\theta_1<100$; $-5<\theta_2<0$; $0<\theta_3<1.5$; $1<\theta_4<6$; $0<\theta_5<1.5$.) return log posterior

else return $-\infty$

This will do the truncation job for you.

I also recommend the following R packages for an easier implementation:

  1. mcmc. Look at the command metrop().
  2. Rtwalk. A self adaptive MCMC sampler based on Metropolis steps.
$\endgroup$
  • $\begingroup$ Thank you. Is logprior($theta_p$)=log((1/100)*(1/5)*(1/1.5)*(1/5)*(1/1.5)) if all the parameters are in the ranges? otherwise it is 0, right? $\endgroup$ – D.Leo Jul 12 '16 at 14:23
  • $\begingroup$ @D.Leo No, it is not. $\theta_p$ is a random number generated from the instrumental distribution. I recommend you to read carefully how the Metropolis Hastings works: en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm $\endgroup$ – Tube Jul 12 '16 at 14:25
  • $\begingroup$ Yes, $\theta_p$ results from the perturbation of $\theta_c$: for each $\theta_{c_{i}}$ I define a normal distribution with mean $\theta_{c_{i}}$ and standard deviation $\theta_{c_{i}}*j_i$ and sample the value of each $\theta_{p_{i}}$. The problem is that I don't understand how to calculate $logprior(\theta_c)$ and $logprior(\theta_p)$ $\endgroup$ – D.Leo Jul 12 '16 at 14:35
  • $\begingroup$ @D.Leo You need to come up with a prior distribution first! You can start with a uniform prior (uniform distribution on the corresponding support), but this is something the user has to specify: en.wikipedia.org/wiki/Prior_probability $\endgroup$ – Tube Jul 12 '16 at 14:40
  • $\begingroup$ I would define a prior Uniform distribution that for each element of theta has the lower bound (ub) and upper bound (ub) specified before. And this is straightforward because the value of the probability density is 1/(ub-lb) if the variable is between lb and ub. Then the value of the prior distribution is given by the product of all the probability density obtained is there is no correlation. The thing that I don't understand is how to formulate the logarithmic prior distribution for $\theta$, as they did. $\endgroup$ – D.Leo Jul 12 '16 at 14:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.