# Implementing Metropolis-Hastings algorithm [duplicate]

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(xt, ·) & 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?

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 = {x01, x02, …, X0P}
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 α(xt, 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 < α(xt, y)

xt+1 = y // accept the proposal

else

xt+1 = xt // reject the proposal

• Save xt1, xt2,…, xtP

3. Return {x11, x21, …., xN1} {x12, x22, …., xN2} . . {x1P, x2P, …., xNP}

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) {

}


## marked as duplicate by John, Juho Kokkala, Michael Chernick, Peter Flom♦Jun 27 '17 at 21:14

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))