2
$\begingroup$

I attempt to sample the posterior distribution of parameters($\alpha,\alpha_1,\alpha_2,\lambda_1,\lambda_2$) using the Metropolis-Hastings within Gibbs.

Steps of Random Walk Metropolis-Hastings

  1. start with t=1 and the initial values {($\alpha=1,\alpha_1=1,\alpha_2=1,\lambda_1=1,\lambda_2$=1)}
  2. using the initial values generate candidate points from proposal densities

    alpha.can<-abs(rnorm(1,alpha[t-1],0.5)); alpha1.can<-abs(rnorm(1,alpha1[t-1],0.5)); alpha2.can<-abs(rnorm(1,alpha2[t-1],0.5)); lamda1.can<-abs(rnorm(1,lamda1[t-1],0.5)); lamda2.can<-abs(rnorm(1,lamda2[t-1],0.005))

  3. Generate u.$\alpha$~$\mathscr U$[0,1], u.$\alpha_1$~$\mathscr U$[0,1], u.$\alpha_2$~$\mathscr U$[0,1], u.$\lambda_1$~$\mathscr U$[0,1], u.$\lambda_2$~$\mathscr U$[0,1].

  4. Calculate the ratio of candidate point alpha.can and previous point alpha[t-1]

R1=$\frac{g(alpha.can \quad given \quad {\alpha_1[t-1],\alpha_2[t-1],\lambda_1[t-1],\lambda_2[t-1]})}{g(alpha[t-1] \quad given \quad {\alpha_1[t-1],\alpha_2[t-1],\lambda_1[t-1],\lambda_2[t-1]})}$ where g is the full conditional posterior of $\alpha$.

  1. If u.$\alpha \;$ <= min (R1,1), accept the candidate point, otherwise set alpha.can=alpha[t-1]

  2. Similarly from step 4, the ratio of candidate point alpha1.cand and previous point alpha1 [t-1] is

R2=$\frac{g(alpha1.can \quad given \quad {\alpha[t],\alpha_2[t-1],\lambda_1[t-1],\lambda_2[t-1]})}{g(alpha1[t-1] \quad given \quad {\alpha[t],\alpha_2[t-1],\lambda_1[t-1],\lambda_2[t-1]})}$ where g is the full conditional posterior of $\alpha_1$.

  1. If u.$\alpha_1 \;$ <= min (R2,1), accept the candidate point, otherwise set alpha1.can=alpha1[t-1]$\qquad$ and so on for $\alpha_2,\lambda_1,\lambda_2$.

  2. Repeat steps 2-7 (T=10000) times and obtain sample of {$\alpha,\alpha_1,\alpha_2,\lambda_1,\lambda_2 $}.

Thanks in Advance

$\endgroup$
  • $\begingroup$ @Xi'an, Thanks very much for your attention. Your correction give me results but Auto-correlation plot for[ four parameters (alpha,alpha1,alpha2 and lamda1)] decreases slowly. Also Histogram plot for alpha1 and alpha2 Concentrated around first bar in histogram and approximately zero for others. $\endgroup$ – Manal Sep 3 '17 at 16:19
  • $\begingroup$ I searched for decreases the Auto-correlation plot and I found that I need to Thinning the MCMC chain $\endgroup$ – Manal Sep 3 '17 at 16:25
  • $\begingroup$ In Acceptance probability, I used log posterior and multiplied Acceptance probability by the ratio suggested by you whereas for the parameters (parameter.can e.g. alpha.can I used alpha.can<-abs(rnorm(1,alphat,sd.Normal)) . when I try alpha.can<-dnorm((-1)^(0:1)*alphat,alphat,sd.Normal) give me this warning message In alpha1[t] <- alpha1.can number of items to replace is not a multiple of replacement length $\endgroup$ – Manal Sep 3 '17 at 16:40
3
$\begingroup$

The proposal of a absolute Gaussian/Normal variate $|X|$ when $X\sim\mathcal{N}(\mu,\sigma^2)$ has density $$\frac{1}{2}\varphi(x|\mu,\sigma^2)+\frac{1}{2}\varphi(-x|\mu,\sigma^2)$$ Hence, if you use the absolute value of a normal simulation for proposed move, your proposal density is no longer symmetric and it must appear in the Metropolis-Hastings acceptance ratio. Therefore, the Metropolis-Hastings ratios in your R code must be multiplied by the ratios $$\dfrac{\varphi(\theta_i^t|\theta_i,\sigma^2)+\varphi(-\theta_i^t|\theta_i,\sigma^2)}{\varphi(\theta_i|\theta_i^t,\sigma^2)+\varphi(-\theta_i|\theta_i^t,\sigma^2)}$$ For instance,

 R1=Post.alpha(Data,alpha.can,alpha1t,alpha2t,lamda1t,lamda2t)/
     Post.alpha(Data,alphat,alpha1t,alpha2t,lamda1t,lamda2t)*
      sum(dnorm((-1)^(0:1)*alphat,alpha.can,sd.Normal))/
       sum(dnorm((-1)^(0:1)*alpha.can,alphat,sd.Normal))

As for the appearance of NaN in the outcome of the Post.alpha functions, it is no surprise as they involve powers and exponentials. For instance,

> Post.alpha(Data,.0001,.0002,.0001,.0002,.0003)
[1] 8.512619e-223

You should work with log posteriors to reduce the possibility of underflows.

$\endgroup$

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.