Metropolis Hastings for Poisson Distribution

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.

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)

  lam <- exp(rnorm(1,mean=log(MCMC[i]),sd=1))

• In your R code the proposal is a binomial lam <- rbinom(1,100,MCMC[i]/(100)) hence it should appear in the acceptance ratio. Commented Feb 23, 2021 at 21:31