Implementing Metropolis-Hastings algorithm

Question

I want to use this algorithm as a black-box, I'll be implementing it either in Python or R, but I don’t really understand it well to be able to turn it into a program.

How do we choose the initial values for parameters?

What exactly is q(x_t, ·) & q(y,x), I understand these are probabilities, but how do we calculate them? and what function is f(y)?

and finally, Do I use the average of the values in the output vectors for each parameter?

Thanks in advance!

Metropolis-Hasting Algorithm for parameter estimation

P – total no. of parameters

N – total no. of iterations

t – iteration variable

1. Initialize the parameter vector x0 = {x₀¹, x₀², …, X₀^P}

2. For iteration t = 1, 2, . . . N do

   • Generate a candidate ‘y’ for next sample i.e. draw from the distribution q(xt, ·)
   • Calculate acceptance ratio α(x_t, y) using α(x, y) = min⁡{1,(f(y) q(y,x))/(f(x) q(x,y) ) }
   • Draw u ~ U[0, 1] i.e. Uniform Distribution

   • Accept/Reject state depending on acceptance ratio

     If u < α(x_t, y)

     x_t+1 = y // accept the proposal

     else

     x_t+1 = x_t // reject the proposal

   • Save x_t¹, x_t²,…, x_t^P

3. Return {x₁¹, x₂¹, …., x_N¹} {x₁², x₂², …., x_N²} . . {x₁^P, x₂^P, …., x_N^P}

EDIT1: My question seeks a clear-cut translation of, let's say, α(x, y) = min⁡{1,(f(y) q(y,x))/(f(x) q(x,y) ) } How does this look like in a program? what exactly is the meaning of q(x, y)? how does the function look like? without any terminologies related to probability or statistics. Apologies if I was not clear before.

EDIT2: How should the functions 'proposal_func' & 'alpha_func' be implemented? (if rest of the program looks fine)

values <- rep(0, 100)
current <- 1

for (i in 1:length(values)) {
   proposal <- proposal_func(current)
   accept_ratio <- alpha_func(current, proposal)
   u <- runif(1,0,1)   
   if (u < accept_ratio) {
      values[i] <- proposal
      current <- proposal      
   } else {
      values[i] <- current
   }
}

alpha_func <- function(current, proposal) {

}

proposal_func <- function(current) {

}

Answer 1

Here is an illustrative implementation for simulating draws from a N(0,1)-distribution using a Metropolis-Hastings algorithm with a random walk kernel with Laplace distributed increments.

Hence, $P=1$ here, as we only simulate a univariate series. $f(y)$ is the standard normal density here, that is, the density from which you aim to generate draws. You take these as they are and do not average. $q(x,y)$ is the transition density, i.e., the probability of being in a new state $y$ given that the chain currently is in state $x$.

Of course, in practice you would never simulate from a standard normal in this fashion!

library(MASS)
library(VGAM)

A <- 1
N <- 50000
R <- 1000 # burn-in draws
x <- rep(NA,N+R)
x[1] <- 0 # initial value - as we have a burn-in, its choice usually is not important

for(i in 2:(N+R)) {
  u <- rlaplace(1, scale = A) # draws from Laplace distribution
  y <- x[i-1]+u # transition density q for generating y is random walk: previous value + Laplace increment
  x[i] <- ifelse(runif(1)<(dnorm(y)/dnorm(x[i-1])),y,x[i-1]) # runif(1) is uniform. Since it is l.th. 1 with prob. 1, we may omit the minimum qualifier
  # with a r.w. kernel and symmetric increments, q(x[i-1],y)=q(y,x[i-1]) so that part cancels out
  }

# Histogram of marginal distribution
x <- x[(R+1):(R+N)] # discard burn-in

truehist(x)
g <- seq(-4,4,length=500)
lines(g,dnorm(g))