How Can I Model Poisson(2) with a Metropolis Hastings Algorithm in R?

[chewingum]

Hi everyone I'm in my final year of university and doing a project on MCMC mainly for applications to bayesian statistics, I think I understand the concept of it so far however I'm struggling to actually make examples.

I'm trying to model Poisson(2) with a random walk metropolis hastings algorithm in R but really don't know where to even start as I've never used R and just confusing myself.

theta <- 0 
set.seed(1)
k<-1
lambda<-2
lik<-function(lambda)(k^(lambda)*exp(-lambda))/factorial(k)
alpha<-function(theta,phi) min(lik(phi)/lik(theta),1)
THETA<-NULL
b<-runif(1, min=0, max=1)
for(i in 1: 100)
{

    if( alpha(theta,phi) < b) 
    {theta.star<-theta
   }
    else{theta.star<-theta.star} 
    theta = c(theta, theta.star) 
}    
plot(density(theta))

I know this is wrong but I don't think I'm too far off, could anyone point me in the right direction? :)

Thanks

 

[Valkarie]

a lot!It looks like you're on the right track. However, there are some mistakes in your code that need to be corrected. Firstly, you should create a function that generates the Metropolis-Hastings algorithm to sample from the posterior distribution. This would take the current position (theta), the prior distribution (lambda) and the likelihood function (lik) as inputs, and return the next position (theta.star). A good starting point is the following:mh_sample <- function(theta, lambda, lik){ alpha <- min(lik(theta)/lik(lambda), 1) b <- runif(1, min = 0, max = 1) if (alpha < b) { theta.star <- theta } else { theta.star <- rnorm(1, mean = theta, sd = 1) } return(theta.star)}Then, you can use this to generate a sequence of samples by looping through the Metropolis-Hastings algorithm. Here's an example of how you could do this:set.seed(1)k <- 1lambda <- 2lik <- function(lambda)(k^(lambda)*exp(-lambda))/factorial(k)theta <- 0THETA <- NULLfor (i in 1:100){ theta.star <- mh_sample(theta, lambda, lik) theta <- theta.star THETA <- c(THETA, theta.star)}plot(density(THETA))Hopefully this helps to clarify how to use the Metropolis-Hastings algorithm for Bayesian statistics.