Simulating a Probit model using Metropolis-Hastings Algorithm (MCMC)

Question

I am working in R and here is my code:

n <- 1000
trueB0 <- 0
trueB1 <- 0.1
prior.mean0 <- 0
prior.sd0 <- 10
prior.mean1 <- 0
prior.sd1 <- 10

x <- rnorm( n, 0, 1 )
eta <- trueB0 + trueB1 * x
P <- pnorm(eta)
Y <- rbinom(n, 1, P)

prior <- function(param){
  b0=param[1]
  b1=param[2]
  B0prior <- dnorm(b0, mean = prior.mean0, sd = prior.sd0, log = T)
  B1prior <- dnorm(b1, mean = prior.mean1, sd = prior.sd1, log = T)
  return(B0prior + B1prior)
}

likelihood <- function(param){
  b0=param[1]
  b1=param[2]
  eta.l <- b0 + b1 * x
  P.l <- pnorm(eta.l)
  singlelikelihoods = dbinom(Y, 1, P.l, log = T)
  sumall=sum(singlelikelihoods)
  return(sumall)
}

posterior <- function(param){
  return(prior(param)+likelihood(param))
}

rho <- 0
sd0 <- 1
sd1 <- 1
sigma <- matrix(c(sd0^2, sd0*sd1*rho, sd0*sd1*rho, sd1^2), ncol = 2)

proposalfunction <- function(param){
  return(rmvnorm(1, mean = param, sigma))
}

q <- function(m1, m2){
  return(dmvnorm(m1, mean = m2, sigma, log = T))
}

run_metropolis_MCMC <- function(initialvalue, iterations){
  chain = array(dim = c(iterations+1,2))
  chain[1,] = initialvalue
  for (i in 1:iterations){
    proposal = proposalfunction(chain[i, ])
    #bug here in the q evaluatoin part
    alpha = exp(posterior(proposal) - posterior(chain[i, ]) + q(chain[i, ], proposal) - q(proposal, chain[i, ]))
    if (runif(1) < min(1,alpha)) 
    {
      chain[i+1,] = proposal
    }
    else
    {
      chain[i+1,] = chain[i,]
    }
  }  
  return(chain)
}

initialvalue <- c(0 , 0.1)
iterations <- 10000
MH = run_metropolis_MCMC(initialvalue, iterations)

I guess the issue is that when computing the alpha in the MH part, I should be getting q(beta|beta*)/q(beta*|beta) (q is the proposal function, beta is the current chain value and beta* is the value generated by q using beta) . However I cannot really figure out how can I get it in R code. (I know q should be a bivariate normal distribution density but I dont know how to implement beta and beta* into the function and get the evaluation of q)

Answer 1

I run your code and if I add library(mvtnorm) I do not get any bug in the q evaluation. I am actually confused what your question is because your q function gives exactly that: q(chain[i, ], proposal) = q(beta|beta*) for a bivariate normal (well it gives the log density for that).

A couple general comments for your code:

• runif(1) < min(1,alpha) would be the same with runif(1) < alpha

• It's good if your code counts the number of accepted proposals, so you know the percentage of that. If that percentage is small then your proposal sd is large and vice versa. Add n.accept <-0 before your for loop and n.accept <- n.accept + 1 at the first if statement. Then add before return

  print(c("percentage accepted draws for betas:", 100*n.accept/iterations))

Please mark this answer as accepted if it helped you.