**bayesianbiologist » Rstats**, and kindly contributed to R-bloggers)

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.

(log)Likelihood:

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

Prior:

prior_reg<-function(pars) { 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. prior_epsilon<-dgamma(epsilon,1,1/100,log=TRUE) return(prior_a + prior_b + prior_epsilon) }

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

x<-runif(30,5,15) y<-x+rnorm(30,0,5) ##Slope=1, intercept=0, epsilon=5 d<-cbind(x,y) plot(x,y)

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.

mcmc_r<-Metro_Hastings(li_func=li_reg,pars=c(0,1,1), par_names=c('a','b','epsilon'),data=d)

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

plotMH(mcmc_r)

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

4) Print posterior Credible Intervals.

BCI(mcmc_r) # 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.

**leave a comment**for the author, please follow the link and comment on his blog:

**bayesianbiologist » Rstats**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...