ABC: Perform an ABC analysis
In dennisprangle/rainfallABC: ABC analysis of rainfall data

Description Usage Arguments Value Examples

Perform an ABC analysis, including simulation of data

1
2
3

ABC(nsims, prior_function, sim_function, sumstat_functions, obs_data,
  dist = "weightedEuclidean", nacc = 200, output = "accepted",
  nCores = parallel::detectCores(), cutoff = 0)

`nsims`	How many simulations to perform.
`prior_function`	A function with no arguments returing a vector of parameters drawn from the prior. Ideally this vector should be named.
`sim_function`	A function taking a vector of parameters as an argument and returning a model data object.
`sumstat_functions`	A list of functions. Each should take a model data object as argument and return a summary statistic (or a vector, see later), which ideally is named. Sometimes it is convenient to use the same calculations for several summaries. In this case a function can return a vector of summaries. Again these should ideally be named.
`obs_data`	Observed model data object.
`dist`	What distance function to use. Currently choose either "weightedEuclidean" or "rank".
`nacc`	How many simulations to accept.
`output`	Specify what output to return. Currently choose either "accepted" or "all".
`nCores`	Number of parallel cores to use.
`cutoff`	If cutoff>0, rain day in real data defined as day where rain>=cutoff, otherwise rain day in real data defined as day where rain>0.

A data frame containing parameters and summaries. If output=="accepted" then only the accepted simulations are returned. Otherwise all simulations are returns and an output=="accepted" column shows if they were accepted.

library(ismev)
data(rain)
doSim = function(pars) {
    rainfallABC::sim_simple_data(pars[1], pars[2], pars[3], 17531, 1) ##TO DO: un-hardcode dur=17531
}
rprior = function() {
    c(lambda = runif(1, 0, 2),
      mu_x = runif(1, 0, 50),
      eta = runif(1, 0, 10))
}
abcout = ABC(200, rprior, doSim, list(mean_overall, prop_dry, cor_overall, max_rain), rain, nacc=50)
par(mfrow=c(1,3))
for (i in 1:3) {
    pname = names(abcout)[i]
    hist(abcout[,i], main=pname, xlab=pname, freq=FALSE)
}