This question already has an answer here:

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?

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


      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

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.


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!


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

g <- seq(-4,4,length=500)

enter image description here


Not the answer you're looking for? Browse other questions tagged or ask your own question.