General Bayesian estimation using MHadaptive

February 6, 2012
By

(This article was first published on 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.


To 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...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Tags: ,

Comments are closed.