SimInf (version 8.3.0)

abc: Approximate Bayesian computation

Description

Approximate Bayesian computation

Usage

abc(model, priors, ngen, npart, fn, ..., verbose = getOption("verbose", FALSE))
# S4 method for SimInf_model
abc(model, priors, ngen, npart, fn, ..., verbose = getOption("verbose", FALSE))

Arguments

model

The model to generate data from.

priors

The priors for the parameters to fit. Each prior is specified with a formula notation, for example, beta ~ uniform(0, 1) specifies that beta is uniformly distributed between 0 and 1. Use c() to provide more than one prior, for example, c(beta ~ uniform(0, 1), gamma ~ normal(10, 1). The following distributions are supported: gamma, normal and uniform. All parameters in priors must be either in gdata or ldata.

ngen

The number of generations of ABC-SMC to run.

npart

An integer specifying the number of particles.

fn

A function for calculating the summary statistics for a simulated trajectory. For each particle, the function must determine if the particle should be accepted (TRUE) or rejected (FALSE) and return that information using the abc_accept function. The first argument passed to the fn function is the result from a run of the model and it contains one trajectory. The second argument to fn is an integer with the generation of the particle(s). Depending on the underlying model structure, data for one or more particles have been generated in each call to fn. If the model only contains one node and all parameters to fit are in ldata, then that node will be replicated and each of the replicated nodes represent one particle in the trajectory (see ‘Examples’). On the other hand if the model contains multiple nodes or the parameters to fit are contained in gdata, then the trajectory in the result argument represents one particle.

...

Further arguments to be passed to fn.

verbose

prints diagnostic messages when TRUE. The default is to retrieve the global option verbose and use FALSE if it is not set.

Value

A SimInf_abc object.

References

T. Toni, D. Welch, N. Strelkowa, A. Ipsen, and M. P. H. Stumpf. Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems. Journal of the Royal Society Interface 6, 187--202, 2009. 10.1098/rsif.2008.0172

Examples

# NOT RUN {
## Let us consider an SIR model in a closed population with N = 100
## individuals of whom one is initially infectious and the rest are
## susceptible. First, generate one realisation (with a specified
## seed) from the model with known parameters \code{beta = 0.16} and
## \code{gamma = 0.077}. Then, use \code{abc} to infer the (known)
## parameters from the simulated data.
model <- SIR(u0 = data.frame(S = 99, I = 1, R = 0),
             tspan = 1:100,
             beta = 0.16,
             gamma = 0.077)

## Run the SIR model and plot the number of infectious.
set.seed(22)
infectious <- trajectory(run(model), "I")$I
plot(infectious, type = "s")

## The distance function to accept or reject a proposal. Each node
## in the simulated trajectory (contained in the 'result' object)
## represents one proposal. The 'generation' argument is the current
## generation of proposals.
acceptFun <- function(result, generation, tol, ptol, ...) {
    ## Determine the tolerance for this generation.
    tol <- tol * (ptol)^(generation - 1)

    ## Extract the time-series of infectious in each node as a
    ## data.frame.
    sim <- trajectory(result, "I")

    ## Split the 'sim' data.frame by node and calculate the sum of the
    ## squared distance at each time-point for each node.
    dist <- tapply(sim$I, sim$node, function(sim_infectious) {
        sum((infectious - sim_infectious)^2)
    })

    ## Return TRUE or FALSE for each node depending on if the distance
    ## is less than the tolerance.
    abc_accept(dist < tol, tol)
}

## Fit the model parameters using ABC-SMC. The priors for the
## paramters are specified in the second argument using a formula
## notation. Here we use a uniform distribtion for each parameter with
## lower bound = 0 and upper bound = 1. Note that we use a low number
## particles here to keep the run-time of the example short.  In
## practice you would want to use many more to ensure better
## approximations.
fit <- abc(model = model,
           priors = c(beta ~ uniform(0, 1), gamma ~ uniform(0, 1)),
           ngen = 4,
           npart = 100,
           fn = acceptFun,
           tol = 5000,
           ptol = 0.9)

## Print a brief summary.
fit

## Display the ABC posterior distribution.
plot(fit)

## Run one more generation.
fit <- continue(fit, tol = 5000, ptol = 0.9)

plot(fit)
# }