**Wiekvoet**, and kindly contributed to R-bloggers)

There are many ways to run general Bayesian calculations in or from R. The best known are JAGS, OpenBUGS and STAN. Then some time ago Rasmus Bååth had a post Three ways to run Bayesian models in R in which he mentioned LaplacesDemon (not on CRAN) on top of those. A check of the Bayes task view gives ‘MCMCpack (…) contains a generic Metropolis sampler that can be used to fit arbitrary models’, ‘The mcmc package consists of an R function for a random-walk Metropolis algorithm for a continuous random vector’ and ‘Note that rcppbugs is a package that attempts to provide a pure R alternative to using OpenBUGS/WinBUGS/JAGS for MCMC’.

Clearly there is a wealth of approaches and each will have its strengths and weaknesses. To get an idea on their comparison I decided to run a number of calculations through all of them. The source of these calculations will be the fine SAS manual, where PROC MCMC has 21 examples. The first example, which I will try this week, is drawing samples from a number of distributions, from which I selected the mixture of three normal distributions.

### Task

Draw samples from a mixture of three normal distributions. Let me start by stating that I would not use anything as complex as an MCMC sampler for this task, my personal SAS solution would be:

data p1;

do i=1 to 1000;

s=rand(‘table’, .3, .4);

select (s);

when (1) r=rand(‘normal’, -3, 2);

when (2) r=rand(‘normal’, 2, 1);

when (3) r=rand(‘normal’, 10, 4);

end;

output;

end;

keep r;

run;

proc sgplot data=p1;

density r / type=kernel (c=.5);

run;

The equivalent R code:

library(plyr)

library(dplyr)

nu <- c(-3,2,10)

sd <- c(2,1,4)

tau <- sd^-2

wgt <- c(.3,.4,.3)

sampls <- sample(1:3,1000,replace=TRUE,wgt)

sampls <- rnorm(1000,nu[sampls],sd[sampls])

density(sampls) %>% plot

#### MCMCpack

Just having a density was enough for MCMCpack. The samples have quite a long autocorrelation but do follow the distribution(not shown).

library(MCMCpack)

mydens <- function(x,…) {

step1 <- sum(wgt*dnorm(x,nu,sd))

log(step1)

}

mcmcpack <- MCMCmetrop1R(mydens,theta.init=1)

acf(mcmcpack)

#### mcmc

mcmc has the property that its result is not the samples, but rather summary statistics of the samples. Hence it does not give samples from the desired distribution.

#### JAGS

The JAGS samples had no autocorrelation and followed the distribution.

library(R2jags)

jmodel <- function() {

i ~ dcat(wgt)

nunu <- (i==1)*nu[1]+

(i==2)*nu[2]+

(i==3)*nu[3]

tautau <- (i==1)*tau[1]+

(i==2)*tau[2]+

(i==3)*tau[3]

p ~ dnorm(nunu,tautau)

}

datain <- list(nu=nu,tau=tau,wgt=wgt)

parameters <- c(‘i’,’p’)

jj <- jags(model.file=jmodel,data=datain,parameters=parameters,DIC=FALSE)

#### STAN

I did not manage to get STAN to give me the samples.

#### rccpbugs

I did not manage to get rccpbugs to give me the samples.

#### LaplacesDemon

LaplacesDemon feels like the most complex of direct R solutions.

library(‘LaplacesDemon’)

mon.names <- “LP”

parm.names <- ‘x’

MyData <- list(mon.names=mon.names, parm.names=parm.names)

N <-1

Model <- function(parm, Data)

{

x <- parm[1]

LL <- sum(wgt*dnorm(x,nu,sd)) %>% log

Modelout <- list(LP=LL, Dev=-2*LL, Monitor=LL,

yhat=x, parm=parm)

return(Modelout)

}

Initial.Values <- 1

Fit <- LaplacesDemon(Model, Data=MyData, Initial.Values)

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

**Wiekvoet**.

R-bloggers.com offers

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