Algorithms for the Simulated Annealing approach to Approximate Bayesian Computation (SABC).
SABC(r.model, r.prior, d.prior, n.sample, eps.init, iter.max,
v=ifelse(method=="informative",0.4,1.2), beta=0.8,
delta=0.1, resample=5*n.sample, verbose=n.sample,
method="noninformative", adaptjump=TRUE,
summarystats=FALSE, y=NULL, f.summarystats=NULL, ...)Returns a list with the following components:
Matrix with ensemble of samples.
Matrix with prior ensemble of samples.
Value of tolerance (temperature) at final iteration.
Effective sample size, due to final bias correction (informative algorithm only).
Function that returns either a random sample from the likelihood or a scalar distance between such a sample and the data. The first argument must be the parameter vector.
Function that returns a random sample from the prior.
Function that returns the density of the prior distribution.
Size of the ensemble.
Initial tolerance or temperature.
Total number of simulations from the likelihood.
Tuning parameter that governs the annealing speed. Defaults to 1.2, for the noninformative algorithm and 0.4, for the informative algorithm.
Tuning parameter that governs the mixing in parameter space. Defaults to 0.8.
Tuning parameter for the resampling steps. Defaults to 0.1.
Number of accepted particle updates after which a resampling step is performed. Defaults to 5*n.sample.
Shows the iteration progress each verbose simulations from the likelihood. NULL for no output. Defaults to verbose = n.sample.
Whether to adapt covariance of jump distribution. Default is TRUE.
Argument to select algorithm. Accepts noninformative or informative.
Whether summary statistics shall be calculated (semi-) automatically. Defaults to FALSE.
Data vector. Needs to be provided if either summarystats = TRUE or if r.model returns a sample from the likelihood.
If summarystats = TRUE this function is needed for the calculation of the summary statistics. Defaults to f.summarystats(x)=(x,x^2,x^3), where the powers are to be understood element-wise.
further arguments passed to r.model
Carlo Albert <carlo.albert@eawag.ch>, Andreas Scheidegger, Tobia Fasciati. Package initially compiled by Lukas M. Weber.
SABC defines a class of algorithms for particle ABC that are inspired by Simulated Annealing. Unlike other algorithms, this class is not based on importance sampling, and hence does not suffer from a loss of effective sample size due to re-sampling. The approach is presented in detail in Albert, Kuensch, and Scheidegger (2014; see references).
This package implements two versions of SABC algorithms, for the cases of a non-informative or an informative prior. These are described in detail in the paper. The algorithms can be selected using the method argument: method=noninformative or method=informative.
In the informative case, the algorithm corrects for the bias caused by an over- or under-representation of the prior.
The argument adaptjump allows a choice of whether to adapt the covariance of the jump distribution. Default is TRUE.
Furthermore, the package allows for three different ways of using the data.
If y is not provided, the algorithm expects r.model to return a scalar measuring the distance between a random sample from the likelihood and the data.
If y is provided and summarystats = FALSE, the algorithm expects r.model to return a random sample from the likelihood and uses the relative sum of squares to measure the distances between y and random likelihood samples.
If summarystats = TRUE the algorithm calculates summary statistics semi-automatically, as described in detail in the paper by Fearnhead et al. (2012; see references).
The summary statistics are calculated by means of a linear regression applied to a sample from the prior and the image of f.summarystats of an associated sample from the likelihood.
C. Albert, H. R. Kuensch and A. Scheidegger, Statistics and Computing 0960-3174 (2014), arXiv:1208.2157, A Simulated Annealing Approach to Approximate Bayes Computations.
P. Fearnhead and D. Prangle, J. R. Statist. Soc. B 74 (2012), Constructing summary statistics for approximate Bayesian computation: semi-automatic approximate Bayesian computation.
if (FALSE) {
## Example for "noninformative" case
# Prior is uniform on [-10,10]
d.prior <- function(par)
dunif(par,-10,10)
r.prior <- function()
runif(1,-10,10)
# Model is the sum of two normal distributions. Return distance to observation 0:
f.dist <- function(par)
return( abs(rnorm( 1 , par , ifelse(runif(1)<0.5,1,0.1 ) )))
# Run algorithm ("noninformative" case)
res <- SABC(f.dist,r.prior,d.prior,n.sample=500,eps.init=2,iter.max=50000)
}
if (FALSE)
# Histogram of results
hist(res$E[,1],breaks=200)
Run the code above in your browser using DataLab