General Bayesian estimation using MHadaptive

If you can write the likelihood function for your model, MHadaptive will take care of the rest (ie. all that MCMC business). I wrote this R package to simplify the estimation of posterior distributions of arbitrary models. Here’s how it works:

1) Define your model (ie the likelihood * prior). In this example, lets build a simple linear regression model.


a<-pars[1]      #intercept
b<-pars[2]      #slope
sd_e<-pars[3]   #error (residuals)
pred <- a + b * data[,1]
log_likelihood<-sum( dnorm(data[,2],pred,sd_e, log=TRUE) )
prior<- prior_reg(pars)
return(log_likelihood + prior)


a<-pars[1]          #intercept
b<-pars[2]          #slope
epsilon<-pars[3]    #error

prior_a<-dnorm(a,0,100,log=TRUE)     ## non-informative (flat) priors on all
prior_b<-dnorm(b,0,100,log=TRUE)     ## parameters.

return(prior_a + prior_b + prior_epsilon)

Now lets just simulate some data to give it a test run:

y<-x+rnorm(30,0,5) ##Slope=1, intercept=0, epsilon=5

2) The function Metro_Hastings() does all the work. Just pass your model to this function, sit back and let MC Emmcee take it away.


3) You can view the posteriors of all model parameters using plotMH()


By default, this will also plot the pair-wise correlation between all parameters.

4) Print posterior Credible Intervals.

#              0.025    0.975
# a       -5.3345970 6.841016
# b        0.4216079 1.690075
# epsilon  3.8863393 6.660037

Tadda! Easy-Peasy.

For a more introductory look at Bayesian modeling, check out the slides and script file from an R workshop I gave at McGill on the subject. There are some additional optional arguments you can pass to Metro_Hastings() which you can read all about in the package manual. You can get the package source from cran, or from my github page.

*notes: MHadaptive is only recommended for relatively low dimensional models. Mixing can be quite slow for higher order models.


9 thoughts on “General Bayesian estimation using MHadaptive

  1. This is great, and looks very simple to use! To address the slow-mixing issue, we’ve implemented something very similar, whereby the user inputs their density function and 2 algorithms (parallel adaptive MH and a parallel adaptive Wang-Landau algorithm) conduct automatic estimation of the function. You can find the R package here:

    We’ve also written a bit about why automatic Monte Carlo methods are so crucial:

    • Very cool. I’m checking out PAWL right now. Yours is the first method I have seen which can handle multi-modality in a reliable way. MHadaptive assumes that the target distribution can be adequately approximated by a multivariate normal. Thanks for the note!

  2. Pingback: Visualising the Metropolis-Hastings algorithm « bayesianbiologist

    • There are two ways I can think of doing this. I’m not 100% sure of either method as in each case it is possible to get R^2 < 0. Both approaches treat R^2 as 1-(residual var of y given model)/(total var of y).

      ## Using the posterior estimated residual variance (sigma^2):

      ## Using the actual residual variance over the posterior distribution of slope and intercept:

      • Thanks very much for responding so quickly. The second method does not give r-squared values anything like what I’m getting from lm()… Your first method seems to be right on, thanks again for that.

  3. There are several adaptive MH, then can you give me a reference that explain this method with more details?

  4. I try the code by with mi function and always this error occur: Error in optim(pars, li_func, control = list(fnscale = -1), hessian = TRUE, :
    function cannot be evaluated at initial parameters
    I really don’t know what should I do?please help me.

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s