Implementing a Metropolis Hastings Algorithm in R

Consider a univariate normal model with mean $µ$ and variance $τ$ . Suppose we use a Beta(2,2) prior for $µ$ (somehow we know µ is between zero and one) and a $log-normal(1,10)$ prior for $τ$ (recall that if a random variable $X$ is $log-normal(m,v)$ then $log X$ is $N(m,v))$. Assume a priori that $µ$ and $τ$ are independent. Implement a Metropolis-Hastings algorithm to evaluate the posterior distribution of $µ$ and $τ$ . Remember that you have to jointly accept or reject $µ$ and $τ$. Also compute the posterior probability that $µ$ is bigger than 0.5.

Here are the data:

2.3656491 2.4952035 1.0837817 0.7586751 0.8780483 1.2765341
1.4598699 0.1801679 -1.0093589 1.4870201 -0.1193149 0.2578262


Attempt: (Heres to enter my data)

mu <- rbeta(1,2,2)
tau <- rlnorm(1,1,10)
x <-
c(2.3656491, 2.4952035, 1.0837817, 0.7586751, 0.8780483, 1.2765341,
1.4598699, 0.1801679, -1.0093589, 1.4870201, -0.1193149, 0.2578262)


This is a MH algorithm found online. I am not sure how to apply the above data to it.

# starting point
a0 <- -5
b0 <- -10

# length of chain
nit <- 1000
a <- rep(0,nit)
b <- rep(0,nit)

# initialize
a[1] <- a0
b[1] <- b0

# tuning parameter
s0 <- 2.0 # maximum step size in random walk proposal
# function. try different s0, e.g., 0.1, 1.0, 2.0

# start chain
counter = 0 # monitor number of acceptances.
for( i in 1:nit) {
s <- s0*runif(1)
theta<-2*3.1415926*runif(1)
anew <- a[i] + s*cos(theta) # random walk
bnew <- b[i] + 2*s*sin(theta) # random walk
r <- pdf(anew,bnew)/pdf(a[i],b[i])# acceptance ratio
test <- runif(1)
if(test < r ) # accept proposed moved.
{
a[i+1] <- anew
b[i+1] <- bnew
counter = counter + 1;
}
else # reject proposed move, stay put.
{
a[i+1] <- a[i]
b[i+1] <- b[i]
}
}
# acceptance rate =
counter /nit

• Better to write your own if you want to learn... it's not that hard, just start with the definition of the acceptance ratio and work your way from the inside out, so to speak. – jbowman Oct 8 '13 at 18:36

It is indeed a very poor idea to start learning a topic just from an on-line code with no explanation. Better read a book (like our Introduction to Monte Carlo methods with R!) or an introductory paper and write your own code.

As written, this code proposes a random walk on the parameter $(a,b)$ which in your case could be $$(\text{logit}(\mu),\log(\tau))=\left(\log\left\{\dfrac{\mu}{1-\mu}\right\},\log(\tau)\right)$$ in order to account for the restricted supports of $\mu$ and $\tau$ that are incompatible with a random walk proposal. This random walk proposal is deduced from the lines

  theta<-2*3.1415926*runif(1)
anew <- a[i] + s*cos(theta) # random walk
bnew <- b[i] + 2*s*sin(theta) # random walk


which correspond to simulating a bivariate normal $\text{N}(0,1)$ (a direct call to rnorm(2) would work as well!) and scaling the first component by s and the second by 2s. (There is no reason for this difference, each component requiring its own scale based on acceptance rates.) The line

s <- s0*runif(1)


that changes the scale at each step is actually incorrect because it makes the usual Metropolis-Hasting ratio invalid: in this ratio

r <- pdf(anew,bnew)/pdf(a[i],b[i])# acceptance ratio


pdf means the density of the target distribution which, in your case, could be coded as

pdf <-function(a,b){

dbeta(1/(1+exp(-a)),2,2)*dnorm(b,1,10)*prod(
dnorm(x,mean=1/(1+exp(-a)),sd=exp(b)))*
exp(a)/(1+exp(a))^2}


which multiplies the prior by the likelihood by the Jacobian. This may sound complicated but the explanation is that

1. The random walk operates on the parametrisation $(\text{logit}(\mu),\log(\tau))$ and the Metropolis-Hastings ratio is the ratio of the targets for this parametrisation;
2. The prior is defined for the parametrisation $(\mu,\tau)$ or $(\mu,\log\,\tau)$ and must be transformed into a prior for the same parametrisation $(\text{logit}(\mu),\log(\tau))$, which means using the Normal model on $\log(\tau))$ and the logit transform of the Beta(2,2) prior on $\mu$, hence the Jacobian of the inverse logit transform:$$\dfrac{\exp(a)}{(1+\exp(a))^2}$$

Once you define this function, your program should run, at least from a theoretical perspective. If the acceptance rate is too high or too low you need to calibrate the scale s.