ABC_rejection {EasyABC}

Rejection sampling scheme for ABC

Description

This function launches a series of nb_simul model simulations with model parameters drawn in the prior distribution specified in prior_matrix.

Usage

  ABC_rejection(model, prior, nb_simul, prior_test=NULL, summary_stat_target=NULL,
  tol=NULL, use_seed=FALSE, seed_count=0, n_cluster=1, verbose=FALSE, progress_bar=FALSE)

Arguments

model

a R function implementing the model to be simulated. It must take as arguments a vector of model parameter values and it must return a vector of summary statistics. When using the option use_seed=TRUE, model must take as arguments a vector containing a seed value and the model parameter values. A tutorial is provided in the package's vignette to dynamically link a binary code to a R function. Users may alternatively wish to wrap their binary executables using the provided functions binary_model and binary_model_cluster. The use of these functions is associated with slightly different constraints on the design of the binary code (see binary_model and binary_model_cluster).

prior

a list of prior information. Each element of the list corresponds to a model parameter. The list element must be a vector whose first argument determines the type of prior distribution: possible values are "unif" for a uniform distribution on a segment, "normal" for a normal distribution, "lognormal" for a lognormal distribution or "exponential" for an exponential distribution. The following arguments of the list elements contain the characteritiscs of the prior distribution chosen: for "unif", two numbers must be given: the minimum and maximum values of the uniform distribution; for "normal", two numbers must be given: the mean and standard deviation of the normal distribution; for "lognormal", two numbers must be given: the mean and standard deviation on the log scale of the lognormal distribution; for "exponential", one number must be given: the rate of the exponential distribution. User-defined prior distributions can also be provided. See the vignette for additional information on this topic.

nb_simul

a positive integer equal to the desired number of simulations of the model.

prior_test

a string expressing the constraints between model parameters. This expression will be evaluated as a logical expression, you can use all the logical operators including "<", ">", ... Each parameter should be designated with "X1", "X2", ... in the same order as in the prior definition. If not provided, no constraint will be applied.

summary_stat_target

a vector containing the targeted (observed) summary statistics. If not provided, ABC_rejection only launches the simulations and outputs the simulation results.

tol

tolerance, a strictly positive number (between 0 and 1) indicating the proportion of simulations retained nearest the targeted summary statistics.

use_seed

logical. If FALSE (default), ABC_rejection provides as input to the function model a vector containing the model parameters used for the simulation. If TRUE, ABC_rejection provides as input to the function model a vector containing an integer seed value and the model parameters used for the simulation. In this last case, the seed value should be used by model to initialize its pseudo-random number generators (if model is stochastic).

seed_count

a positive integer, the initial seed value provided to the function model (if use_seed=TRUE). This value is incremented by 1 at each call of the function model.

n_cluster

a positive integer. If larger than 1 (the default value), ABC_rejection will launch model simulations in parallel on n_cluster cores of the computer.

verbose

logical. FALSE by default. If TRUE, ABC_rejection writes in the current directory intermediary results at the end of each step of the algorithm in the file "output". These outputs have a matrix format, in wich each raw is a different simulation, the first columns are the parameters used for this simulation, and the last columns are the summary statistics of this simulation.

progress_bar

logical, FALSE by default. If TRUE, ABC_rejection will output a bar of progression with the estimated remaining computing time. Option not available with multiple cores.

Value

The returned value is a list containing the following components:

param

The model parameters used in the model simulations.

stats

The summary statistics obtained at the end of the model simulations.

weights

The weights of the different model simulations. In the standard rejection scheme, all model simulations have the same weights.

stats_normalization

The standard deviation of the summary statistics across the model simulations.

nsim

The number of model simulations performed.

nrec

The number of retained simulations (if targeted summary statistics are provided).

computime

The computing time to perform the simulations.

Author(s)

Franck Jabot, Thierry Faure and Nicolas Dumoulin

References

Pritchard, J.K., and M.T. Seielstad and A. Perez-Lezaun and M.W. Feldman (1999) Population growth of human Y chromosomes: a study of Y chromosome microsatellites. Molecular Biology and Evolution, 16, 1791–1798.

See Also

binary_model, binary_model_cluster, ABC_sequential, ABC_mcmc

Examples

    ##### EXAMPLE 1 #####
    #####################
    set.seed(1)

    ## artificial example to show how to use the 'ABC_rejection' function.
    ## defining a simple toy model:
    toy_model<-function(x){ 2 * x + 5 + rnorm(1,0,0.1) }

    ## define prior information
    toy_prior=list(c("unif",0,1)) # a uniform prior distribution between 0 and 1

    ## only launching simulations with parameters drawn in the prior distributions
    set.seed(1)
    n=10
    ABC_sim<-ABC_rejection(model=toy_model, prior=toy_prior, nb_simul=n)
    ABC_sim

    ## launching simulations with parameters drawn in the prior distributions
    # and performing the rejection step
    sum_stat_obs=6.5
    tolerance=0.2
    ABC_rej<-ABC_rejection(model=toy_model, prior=toy_prior, nb_simul=n,
      summary_stat_target=sum_stat_obs, tol=tolerance)

    ## NB: see the package's vignette to see how to pipeline 'ABC_rejection' with the function
    # 'abc' of the package 'abc' to perform other rejection schemes.
 ## Not run:  

    ##### EXAMPLE 2 #####
    #####################

    ## this time, the model has two parameters and outputs two summary statistics.
    ## defining a simple toy model:
    toy_model2<-function(x){ c( x[1] + x[2] + rnorm(1,0,0.1) , x[1] * x[2] + rnorm(1,0,0.1) ) }

    ## define prior information
    toy_prior2=list(c("unif",0,1),c("normal",1,2))
    # a uniform prior distribution between 0 and 1 for parameter 1, and a normal distribution
    # of mean 1 and standard deviation of 2 for parameter 2.

    ## only launching simulations with parameters drawn in the prior distributions
    set.seed(1)
    n=10
    ABC_sim<-ABC_rejection(model=toy_model2, prior=toy_prior2, nb_simul=n)
    ABC_sim

    ## launching simulations with parameters drawn in the prior distributions
    # and performing the rejection step
    sum_stat_obs2=c(1.5,0.5)
    tolerance=0.2
    ABC_rej<-ABC_rejection(model=toy_model2, prior=toy_prior2, nb_simul=n,
      summary_stat_target=sum_stat_obs2, tol=tolerance)

    ## NB: see the package's vignette to see how to pipeline 'ABC_rejection' with the function
    # 'abc' of the package 'abc' to perform other rejection schemes.


    ##### EXAMPLE 3 #####
    #####################

    ## this time, the model is a C++ function packed into a R function -- this time, the option
    # 'use_seed' must be turned to TRUE.
    n=10
    trait_prior=list(c("unif",3,5),c("unif",-2.3,1.6),c("unif",-25,125),c("unif",-0.7,3.2))
    trait_prior

    ## only launching simulations with parameters drawn in the prior distributions
    ABC_sim<-ABC_rejection(model=trait_model, prior=trait_prior, nb_simul=n, use_seed=TRUE)
    ABC_sim

    ## launching simulations with parameters drawn in the prior distributions and performing
    # the rejection step
    sum_stat_obs=c(100,2.5,20,30000)
    tolerance=0.2
    ABC_rej<-ABC_rejection(model=trait_model, prior=trait_prior, nb_simul=n,
      summary_stat_target=sum_stat_obs, tol=tolerance, use_seed=TRUE)

    ## NB: see the package's vignette to see how to pipeline 'ABC_rejection' with the function
    # 'abc' of the package 'abc' to perform other rejection schemes.


    ##### EXAMPLE 4 - Parallel implementations #####
    ################################################

    ## NB: the option use_seed must be turned to TRUE.

    ## For models already running with the option use_seed=TRUE, simply change
    # the value of n_cluster:
    sum_stat_obs=c(100,2.5,20,30000)
    ABC_simb<-ABC_rejection(model=trait_model, prior=trait_prior, nb_simul=n,
      use_seed=TRUE, n_cluster=2)

    ## For other models, change the value of n_cluster and modify the model so that the first
    # parameter becomes a seed information value:
    toy_model_parallel<-function(x){ 
	set.seed(x[1])
	2 * x[2] + 5 + rnorm(1,0,0.1) }
    sum_stat_obs=6.5

    ABC_simb<-ABC_rejection(model=toy_model_parallel, prior=toy_prior, nb_simul=n,
      use_seed=TRUE, n_cluster=2)
 
## End(Not run)