From c56a2af0e363f1cafb694ede4e23ba6b3eb6ea18 Mon Sep 17 00:00:00 2001 From: Rasilgon Date: Thu, 30 Apr 2015 16:59:40 +0200 Subject: [PATCH 1/5] Update adaptativeMetropolis.Rmd --- CommentedCode/02-Samplers/MCMC/adaptativeMetropolis.Rmd | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/CommentedCode/02-Samplers/MCMC/adaptativeMetropolis.Rmd b/CommentedCode/02-Samplers/MCMC/adaptativeMetropolis.Rmd index 83c82ba..e4357b7 100644 --- a/CommentedCode/02-Samplers/MCMC/adaptativeMetropolis.Rmd +++ b/CommentedCode/02-Samplers/MCMC/adaptativeMetropolis.Rmd @@ -79,7 +79,7 @@ prior <- function(param){ ## The MCMC -Now comes the mcmc to sample from the posterior. We are using a Metropolis-Hastings MCMC. I you are uncertain how this works, first read http://theoreticalecology.wordpress.com/. +Now comes the mcmc to sample from the posterior. We are using a Metropolis-Hastings MCMC. In this previous version, the proposal function was simlply creating independent normal draws $proposalfunction <- function(param){return(rnorm(3,mean = param, sd= c(0.1,0.5,0.3))) }$ . @@ -116,8 +116,7 @@ run_metropolis_MCMC <- function(iterations){ } ``` -But first of all, we will just run the MCMC like before and check the convergence as explained in http://theoreticalecology.wordpress.com/2011/12/09/mcmc-chain-analysis-and-convergence-diagnostics-with-coda-in-r/ - +But first of all, we will just run the MCMC like before and check the convergence. ```{r convergence} @@ -142,7 +141,7 @@ gelman.plot(combinedchains)
#### Adaptation steps -We use the samples obtained already to adjust the proposal for exaplanations why the scaling factor of 2.38^2/d is optional (see references) +We use the samples obtained already to adjust the proposal for exaplanations why the scaling factor of $2.38^2/d$ is optional (see references) From 2eed11df3bb635f2df46a6d152ce60d54481ed2f Mon Sep 17 00:00:00 2001 From: Rasilgon Date: Thu, 30 Apr 2015 17:03:21 +0200 Subject: [PATCH 2/5] Create easyABC.Rmd --- .../08-ApproximateBayesian/easyABC.Rmd | 65 +++++++++++++++++++ 1 file changed, 65 insertions(+) create mode 100644 CommentedCode/08-ApproximateBayesian/easyABC.Rmd diff --git a/CommentedCode/08-ApproximateBayesian/easyABC.Rmd b/CommentedCode/08-ApproximateBayesian/easyABC.Rmd new file mode 100644 index 0000000..16a51ca --- /dev/null +++ b/CommentedCode/08-ApproximateBayesian/easyABC.Rmd @@ -0,0 +1,65 @@ +--- +output: + html_document: + fig_caption: yes + keep_md: yes +--- +The EasyABC package for Approximate Bayesian Computation in R +==== + +```{r global_options, include=FALSE} +knitr::opts_chunk$set(fig.width=12, fig.height=8) + +``` + +```{r, echo = F} +set.seed(123) +``` + + +Approximate Bayesian Computation (ABC) is a relatively new method that allows treating any stochastic model (IBM, stochastic population model, …) in a statistical framework by generating “approximate” likelihood values through simulating from the model. We provide a gentle introduction to ABC and some alternative approaches in our recent Ecology Letters review on “statisitical inference for stochastic simulation models”. + +ABC has a huge potential as a solution for many typical ecological problems, but to make this more widely known is currently hindered by the fact that you have to code everything by hand, which excludes a large number of users. + +The EasyABC package, available from CRAN (developed by Franck Jabot, Thierry Faure, Nicolas Dumoulin and maintained by Nicolas Dumoulin.), implements a number of algorithms for the three main sampling strategies used in ABC, namely Rejection Sampling, Sequential Monte Carlo (SMC) and Markov Chain Monte Carlo (MCMC). All those are also discussed in our review. The use of the package is relatively straightforward. + + +```{r main} +library(EasyABC) + +# 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 for the model and the data. + +# observed summary statistics +summarydata = c(mean(data), sd(data)) + +# stochastic model generates a sample for given par and returns summary statistics +model <- function(par){ + samples <- rnorm(10, mean =par[1], sd = par[2]) + return(c(mean(samples), sd(samples))) +} + +# call to EasyABC with the ABC-MCMC algorithm Marjoram, P.; Molitor, J.; Plagnol, V. & +# Tavare, S. (2003) Markov chain Monte Carlo without likelihoods. Proc. Natl. Acad. Sci. USA, 100, 15324-15328. +# with some automatic adjustment options +ABC_Marjoram_original<-ABC_mcmc(method="Marjoram", model=model, + prior=list(c("unif",0,10),c("unif",1,5)), + summary_stat_target=summarydata, n_rec = 10000) + + +str(ABC_Marjoram_original) +par(mfrow=c(2,1)) +hist(ABC_Marjoram_original$param[5000:10000,1], main = "Posterior for mean") +hist(ABC_Marjoram_original$param[5000:10000,2], main = "Posterior for standard deviation") +``` +
+
+ +## References for further reading + +Hartig, F., Calabrese, J. M., Reineking, B., Wiegand, T. and Huth, A. (2011), Statistical inference for stochastic simulation models – theory and application. Ecology Letters, 14: 816–827. From 934d0d933ca3c8ecf399549b0ab604e53cbeadf8 Mon Sep 17 00:00:00 2001 From: Rasilgon Date: Thu, 30 Apr 2015 17:11:43 +0200 Subject: [PATCH 3/5] Update easyABC.Rmd --- .../08-ApproximateBayesian/easyABC.Rmd | 30 +++++++++++-------- 1 file changed, 17 insertions(+), 13 deletions(-) diff --git a/CommentedCode/08-ApproximateBayesian/easyABC.Rmd b/CommentedCode/08-ApproximateBayesian/easyABC.Rmd index 16a51ca..968f684 100644 --- a/CommentedCode/08-ApproximateBayesian/easyABC.Rmd +++ b/CommentedCode/08-ApproximateBayesian/easyABC.Rmd @@ -24,17 +24,19 @@ ABC has a huge potential as a solution for many typical ecological problems, but The EasyABC package, available from CRAN (developed by Franck Jabot, Thierry Faure, Nicolas Dumoulin and maintained by Nicolas Dumoulin.), implements a number of algorithms for the three main sampling strategies used in ABC, namely Rejection Sampling, Sequential Monte Carlo (SMC) and Markov Chain Monte Carlo (MCMC). All those are also discussed in our review. The use of the package is relatively straightforward. -```{r main} +Let's create some training data, assuming that there are 10 samples of a normal distribution with mean 5.3 and sd 2.7 +```{r first} +# Load library library(EasyABC) - -# assuming the data are 10 samples of a normal distribution -# with mean 5.3 and sd 2.7 + +#create data 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 for the model and the data. -# 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 for the model and the data. - +```{r second} + # observed summary statistics summarydata = c(mean(data), sd(data)) @@ -43,15 +45,17 @@ model <- function(par){ samples <- rnorm(10, mean =par[1], sd = par[2]) return(c(mean(samples), sd(samples))) } - -# call to EasyABC with the ABC-MCMC algorithm Marjoram, P.; Molitor, J.; Plagnol, V. & -# Tavare, S. (2003) Markov chain Monte Carlo without likelihoods. Proc. Natl. Acad. Sci. USA, 100, 15324-15328. -# with some automatic adjustment options +``` + +
+Now we call to EasyABC with the ABC-MCMC algorithm with some automatic adjustment options + + +```{r ABC} ABC_Marjoram_original<-ABC_mcmc(method="Marjoram", model=model, prior=list(c("unif",0,10),c("unif",1,5)), summary_stat_target=summarydata, n_rec = 10000) - str(ABC_Marjoram_original) par(mfrow=c(2,1)) hist(ABC_Marjoram_original$param[5000:10000,1], main = "Posterior for mean") From 9c6375276ca602a3ecd43de000ac33ec0a86c5a8 Mon Sep 17 00:00:00 2001 From: Rasilgon Date: Thu, 30 Apr 2015 17:13:53 +0200 Subject: [PATCH 4/5] Update easyABC.Rmd --- CommentedCode/08-ApproximateBayesian/easyABC.Rmd | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/CommentedCode/08-ApproximateBayesian/easyABC.Rmd b/CommentedCode/08-ApproximateBayesian/easyABC.Rmd index 968f684..06a8a02 100644 --- a/CommentedCode/08-ApproximateBayesian/easyABC.Rmd +++ b/CommentedCode/08-ApproximateBayesian/easyABC.Rmd @@ -48,7 +48,7 @@ model <- function(par){ ```
-Now we call to EasyABC with the ABC-MCMC algorithm with some automatic adjustment options +Now we call to EasyABC with the ABC-MCMC algorithm (Marjoram et al., 2003) with some automatic adjustment options ```{r ABC} @@ -67,3 +67,5 @@ hist(ABC_Marjoram_original$param[5000:10000,2], main = "Posterior for standard d ## References for further reading Hartig, F., Calabrese, J. M., Reineking, B., Wiegand, T. and Huth, A. (2011), Statistical inference for stochastic simulation models – theory and application. Ecology Letters, 14: 816–827. + +Marjoram, P.; Molitor, J.; Plagnol, V. & Tavaré, S. (2003) Markov chain Monte Carlo without likelihoods From 8a27529a8ef2744ef5d7fedbcfce8be2ed6cb322 Mon Sep 17 00:00:00 2001 From: Rasilgon Date: Thu, 30 Apr 2015 17:18:14 +0200 Subject: [PATCH 5/5] Create abcMCMC.Rmd --- .../08-ApproximateBayesian/abcMCMC.Rmd | 92 +++++++++++++++++++ 1 file changed, 92 insertions(+) create mode 100644 CommentedCode/08-ApproximateBayesian/abcMCMC.Rmd diff --git a/CommentedCode/08-ApproximateBayesian/abcMCMC.Rmd b/CommentedCode/08-ApproximateBayesian/abcMCMC.Rmd new file mode 100644 index 0000000..33478e9 --- /dev/null +++ b/CommentedCode/08-ApproximateBayesian/abcMCMC.Rmd @@ -0,0 +1,92 @@ +--- +output: + html_document: + fig_caption: yes + keep_md: yes +--- +A simple Approximate Bayesian Computation MCMC (ABC-MCMC) in R +==== + +```{r global_options, include=FALSE} +knitr::opts_chunk$set(fig.width=12, fig.height=8) + +``` + +```{r, echo = F} +set.seed(123) +``` + +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 want to have more background on this algorithm, read the excellent paper by Marjoram et al. (2003) who proposed this algorithm for the first time. + + +Let's create the training data, assuming the data are 10 samples of a normal distribution with mean 5.3 and sd 2.7. + +```{r data} +library(coda) + +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. + + +```{r ABC} +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) +} + +``` + +
+Now we plug this in in a standard metropolis hastings MCMC, with the metropolis acceptance exchanged for the ABC acceptance. + +```{r MCMC_ABC} +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) + +``` +
+The result should look something like that: +```{r plot_posterior} +plot(posterior) +``` +
+
+ +## References for further reading + +Marjoram, P.; Molitor, J.; Plagnol, V. & Tavaré, S. (2003) Markov chain Monte Carlo without likelihoods