1

Suppose I perform an experiment where the data have a $Poisson(\lambda)$ sampling density. My uncertainty about $\lambda$ using a Gamma prior density with parameters $\alpha$ and $\beta$. I also describe our uncertainty about $\alpha$ and $\beta$ using independent Gamma prior densities.

I'm trying to implement an MCMC algorithm to calculate posterior densities for $\lambda$, $\alpha$, and $\beta$. What is wrong with my code below? I'm getting an acceptance ratio of 0.

obs <- rpois(50, 5)

m <- 100000
x <- matrix(NA, nrow=m, ncol=3)
x[1,] <- c(5, 1, 1)
accept <- 0
sd1 <- 2
sd2 <- 2
sd3 <- 2

post <- function(lambda, alpha, beta) {
  prod(dpois(x, lambda))*dgamma(lambda, 
 alpha,beta)*dgamma(alpha,1,1)*dgamma(beta,1,1)} 

for(i in 2:m){
  xc <- c(rnorm(2, x[i-1,1]-sd1,x[i-1,1]+sd1),runif(1,x[i-1,2]-sd2,x[i-
1,2]+sd2), runif(1,x[i-1,3]-sd3, x[i-1,3]+sd3))
  if (post(xc[1], xc[2], xc[3])/post(x[i-1,1],x[i-1,2], x[i-1,3]) > runif(1)){
    x[i,] <- xc; accept <- accept + 1
  } else{
    x[i,] <- x[i-1,]
 }
}
accept/m
Amanda R.
  • 287
  • 1
  • 2
  • 17
  • For starters.. `xc[1]` is drawn from a normal distribution (hence, can be negative) and is being used as `lambda`, which is invalid. – Julius Vainora Mar 01 '18 at 23:07
  • @Julius so do I want to draw from unif(0,1) such as: `xc <- c(runif(1, x[i-1,1]-sd1,x[i-1,1]+sd1),runif(1,x[i-1,2]-sd2,x[i- 1,2]+sd2), runif(1,x[i-1,3]-sd3, x[i-1,3]+sd3)) – Amanda R. Mar 01 '18 at 23:16
  • Also `prod(dpois(x, lambda))` returns NA, so whatever `post` returns is going to be NA. Unless I'm misunderstanding something it's never going to be greater than `runif(1)`. – be_green Mar 01 '18 at 23:19
  • @be_green So should it be x[i,1] instead? – Amanda R. Mar 01 '18 at 23:25
  • You don't have any data for the Poisson likelihood – erocoar Mar 01 '18 at 23:52
  • @erocoar oops I forgot to include it, I've added it to the code above – Amanda R. Mar 02 '18 at 00:05

0 Answers0