Metropolis Hastings for Poisson Distribution

Question

Studying Bayesian Inference and Markov chain Monte Carlo (MCMC) algorithms, I am facing a self study question on a MCMC approach to a Poisson distribution with parameter $\lambda$.

Using R, my code is:

### Parameters
# observation
x <- 60 
# alpha parameter
a <- 10 
# beta parameter
b <- 15 

First step is to find the posterior distribution for Poisson with parameter $\lambda$ through a prior distribution and likelihood distribution. After some calculations, I got to find my posterior distribution as a Gamma with parameters alpha $(x+a)$ and beta $(b/(b+1))$.

ap <- x+a
bp <- b/(b+1)

Then, the last step is to perform a MCMC simulation. I am doing 5000 iterations with initial value 20.

inter = 5000
MCMC<- rep(0, inter+1)
MCMC[1] <- 20

for(i in 1:inter) {
  lam <- rbinom(1,100,MCMC[i]/(100)) 
  
  numerator <- dgamma(lam,ap,scale=bp)
  denominator <- dgamma(MCMC[i],ap,scale=bp)
    
  alfa <- min(1, (numerator/denominator))
  print(alfa)
  
  if (runif(1,0,1) <= alfa){
    (MCMC[i+1] <- lam)
  } else {
    (MCMC[i+1] = MCMC[i])
  }
}
plot(0:inter, MCMC, type="l", main="MCMC", ylim=c(0,100))

I am not sure if my approach is correct: I am trying to simulate a proposal for the acceptance through a binomial distribution. Then, I do the ratio between the first value of my MCMC and the proposal.

Answer 1

Guessing from the R code and the question it sounds like one observed $$X\sim\mathcal P(\lambda)$$ as $x=60$ and the prior distribution on $\lambda$ is a Gamma distribution $\mathcal Ga(a,b)$ [using the scale parameterisation of the Gamma distribution] meaning the posterior distribution is indeed $$\mathcal Ga(a+x,b/(b+1))$$ Running a Metropolis-Rosenbluth-Hastings algorithm means simulating a proposed value of $\lambda$ from a proposal $q(\lambda|\lambda^{(t)})$ and accepting with the Metropolis-Rosenbluth-Hastings ratio. However, the proposal used in the code $$\lambda|\lambda^{(t)}\sim\mathcal B(100,\lambda^{(t)}/100)$$ should not be used because it is supported on $\{0,1,\ldots,100\}$ rather than the positive real line. It has the wrong support as it only returns integer values. (Besides, the parameter $n=100$ in the Binomial sounds rather arbitrary for the scale of $\lambda$.) I would thus suggest using a continuous distribution $q(\lambda|\lambda^{(t)})$ on $\lambda$ as for instance a log-normal distribution. (Warning: there is a Jacobian appearing in the Metropolis-Rosenbluth-Hastings ratio.)

As a side remark the correct Metropolis-Rosenbluth-Hastings ratio associated with the Binomial proposal should be

numerador <- dgamma(lam,ap,scale=bp)*dbinom(MCMC[i],100,lam/100)
denominador<-dgamma(MCMC[i],ap,scale=bp)*dbinom(lam,100,MCMC[i]/100)

since the Binomial density is not symmetric and hence must appear in the ratio. With this correction, there is no visible difference in the path of the Markov chain when compared with a log-normal proposal, except that one chain is made of integers and the other one of real numbers:

The modification in the R code for the second (log-normal) case is

  lam <- exp(rnorm(1,mean=log(MCMC[i]),sd=1))
  numerador <- dgamma(lam,ap,scale=bp)*lam
  denominador <- dgamma(MCMC[i],ap,scale=bp)*MCMC[i]