smfsb (version 1.5)

abcSmc: Run an ABC-SMC algorithm for infering the parameters of a forward model

Description

Run an ABC-SMC algorithm for infering the parameters of a forward model. This sequential Monte Carlo algorithm often performs better than simple rejection-ABC in practice.

Usage

abcSmc(N, rprior, dprior, rdist, rperturb, dperturb, factor=10,
                   steps=15, verb=FALSE)

Value

A matrix (or vector) with rows (or elements) representing samples from the approximate posterior distribution.

Arguments

N

An integer representing the number of simulations to pass on at each stage of the SMC algorithm. Note that the TOTAL number of forward simulations required by the algorithm will be (roughly) 'N*steps*factor'.

rprior

A function without arguments generating a single parameter (vector) from prior distribution.

dprior

A function with required argument a model parameter (such as generated by 'rprior') and optional parameter 'log' returing the (log) density of the parameter under the prior distribution.

rdist

A function taking a parameter (vector) as argument and returning a scalar "distance" representing a measure of how good the chosen parameter is. This will typically be computed by first using the parameter to run a forward model, then computing required summary statistics, then computing a distance. See the example for details.

rperturb

A function which takes a parameter as its argument and returns a perturbed parameter from an appropriate kernel.

dperturb

A function which takes a pair of parameters as its first two arguments (new first and old second), and has an optional argument 'log' for whether to return the log of the density associated with this perturbation kernel.

factor

At each step of the algorithm, 'N*factor' proposals are generated and the best 'N' of these are weighted and passed on to the next stage. Note that the effective sample size of the parameters passed on to the next step may be (much) smaller than 'N', since some of the particles may be assigned small (or zero) weight.

steps

The number of steps of the ABC-SMC algorithm. Typically, somewhere between 5 and 100 steps seems to be used in practice.

verb

Boolean indicating whether some progress should be printed to the console (the number of steps remaining).

See Also

pfMLLik, StepGillespie, abcRun, simTs, stepLVc

Examples

# \donttest{
data(LVdata)
rprior <- function() { c(runif(1, -3, 3), runif(1, -8, -2), runif(1, -4, 2)) }
dprior <- function(x, ...) { dunif(x[1], -3, 3, ...) + 
                dunif(x[2], -8, -2, ...) + dunif(x[3], -4, 2, ...) }
rmodel <- function(th) { simTs(c(50,100), 0, 30, 2, stepLVc, exp(th)) }
rperturb <- function(th){th + rnorm(3, 0, 0.5)}
dperturb <- function(thNew, thOld, ...){sum(dnorm(thNew, thOld, 0.5, ...))}
sumStats <- identity
ssd = sumStats(LVperfect)
distance <- function(s) {
    diff = s - ssd
    sqrt(sum(diff*diff))
}
rdist <- function(th) { distance(sumStats(rmodel(th))) }
out = abcSmc(5000, rprior, dprior, rdist, rperturb,
             dperturb, verb=TRUE, steps=6, factor=5)
print(summary(out))
# }