[This article was first published on R snippets, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Bootstrapping is a very popular statistical technique. However, its Bayesian analogue proposed by Rubin (1981) is not very common. I was looking for an example of its implementation in GNU R and could not find one so I decided to write a snippet presenting it.Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
In standard bootstrapping observations are sampled with replacement. This implies that observation weights follow multinomial distribution. In Bayesian bootstrap multinomial distribution is replaced by Dirichlet distribution.
This observation leads to very simple implementation of Bayesian bootstrap using gtools package. Here is the code presenting it with a simple application giving frequentist and Bayesian 95% confidence interval for mean in fbq and bbq variables:
library(gtools)< o:p>
# Bayesian bootstrap
mean.bb <- function(x, n) {< o:p>
apply(rdirichlet(n, rep(1, length(x))), 1, weighted.mean, x = x)< o:p>
}< o:p>
# standard bootstrap
mean.fb <- function(x, n) {< o:p>
mean.fb <- function(x, n) {< o:p>
replicate(n, mean(sample(x, length(x), TRUE)))< o:p>
}< o:p>
set.seed(1)< o:p>
reps <- 100000< o:p>
x <- cars$dist< o:p>
system.time(fbq <- quantile((mean.fb(x, reps)), c(0.025, 0.075)))< o:p>
system.time(bbq <- quantile((mean.bb(x, reps)), c(0.025, 0.075)))As it can be seen implementation of Bayesian bootstrap is fairly simple.
On my computer Bayesian bootstrap is approximately 80% slower than standard bootstrap, but its performance probably could be improved.
To leave a comment for the author, please follow the link and comment on their blog: R snippets.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.