**theoretical ecology » Submitted to R-bloggers**, and kindly contributed to R-bloggers)

Approximate Bayesian Computing and similar techniques, which are based on calculating approximate likelihood values based on samples from a stochastic simulation model, have attracted a lot of attention in the last years, owing to their promise to provide a general statistical technique for stochastic processes of any complexity, without the limitations that apply to “traditional” statistical models due to the problem of maintaining “tractable” likelihood functions. If you are unsure what all this means, I recommend you our recent review on statistical inference for stochastic simulation models, which aims at giving a pedagogical introduction to this exciting topic.

A colleague asked me now for a simple example of the Approximate Bayesian Computation MCMC (ABC-MCMC) algorithm that we discussed in our review. If you want to have more background on this algorithm, read the excellent paper by Marjoram et al. who proposed this algorithm for the first time. Below, I provide a minimal example, similar to my example for a simple Metropolis-Hastings MCMC in R, where the only main difference is that the Metropolis-Hastings acceptance has been changed for an ABC acceptance.

library(coda) # assuming the data are 10 samples of a normal distribution # with mean 5.3 and sd 2.7 data = rnorm(10, mean =5.3, sd = 2.7) # we want to use ABC to infer the parameters that were used. # we sample from the same model and use mean and variance # as summary statstitics. We return true for ABC acceptance when # the difference to the data is smaller than a certain threshold meandata <- mean(data) standarddeviationdata <- sd(data) ABC_acceptance <- function(par){ # prior to avoid negative standard deviation if (par[2] <= 0) return(F) # stochastic model generates a sample for given par samples <- rnorm(10, mean =par[1], sd = par[2]) # comparison with the observed summary statistics diffmean <- abs(mean(samples) - meandata) diffsd <- abs(sd(samples) - standarddeviationdata) if((diffmean < 0.1) & (diffsd < 0.2)) return(T) else return(F) } # we plug this in in a standard metropolis hastings MCMC, # with the metropolis acceptance exchanged for the ABC acceptance run_MCMC_ABC <- function(startvalue, iterations){ chain = array(dim = c(iterations+1,2)) chain[1,] = startvalue for (i in 1:iterations){ # proposalfunction proposal = rnorm(2,mean = chain[i,], sd= c(0.7,0.7)) if(ABC_acceptance(proposal)){ chain[i+1,] = proposal }else{ chain[i+1,] = chain[i,] } } return(mcmc(chain)) } posterior <- run_MCMC_ABC(c(4,2.3),300000) plot(posterior)

The result should look something like that:

* Figure*: Trace and marginal plots for the posterior sample. From the marginal plots to the right, you see that we are approximately retrieving the original parameter values, which were 5.3 and 2.7. If you are unsure how to read these plots, look at this older post.

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

**theoretical ecology » Submitted to R-bloggers**.

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