Adding markdowns for ABC

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)
 <br />
 
 #### 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)
 
 
 

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)
+``` 
+<br /> 
+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)
+}
+
+```
+
+<br />
+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)
+
+```
+<br />
+The result should look something like that:
+```{r plot_posterior}
+plot(posterior)
+```
+<br />
+<br />
+
+## References for further reading
+
+Marjoram, P.; Molitor, J.; Plagnol, V. & Tavaré, S. (2003) Markov chain Monte Carlo without likelihoods

CommentedCode/08-ApproximateBayesian/easyABC.Rmd

@@ -0,0 +1,71 @@
+---
+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. 
+
+
+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)
+
+#create data
+data =  rnorm(10, mean =5.3, sd = 2.7)
+```
+<br />
+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))
+
+# 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)))
+}
+```
+
+<br />
+Now we call to EasyABC with the ABC-MCMC algorithm (Marjoram et al., 2003) 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")
+hist(ABC_Marjoram_original$param[5000:10000,2], main = "Posterior for standard deviation")
+```
+<br />
+<br />
+
+## 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