
Last chance! 50% off unlimited learning
Sale ends in
This function launches a series of nb_simul
model simulations with model parameters drawn in the prior distribution specified in prior_matrix
.
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)
The returned value is a list containing the following components:
The model parameters used in the model
simulations.
The summary statistics obtained at the end of the model
simulations.
The weights of the different model
simulations. In the standard rejection scheme, all model
simulations have the same weights.
The standard deviation of the summary statistics across the model
simulations.
The number of model
simulations performed.
The number of retained simulations (if targeted summary statistics are provided).
The computing time to perform the simulations.
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
).
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.
a positive integer equal to the desired number of simulations of the model.
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.
a vector containing the targeted (observed) summary statistics.
If not provided, ABC_rejection
only launches the simulations and outputs the simulation results.
tolerance, a strictly positive number (between 0 and 1) indicating the proportion of simulations retained nearest the targeted summary statistics.
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).
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
.
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.
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.
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.
Franck Jabot, Thierry Faure and Nicolas Dumoulin
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.
binary_model
, binary_model_cluster
, ABC_sequential
, ABC_mcmc
##### 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.
if (FALSE) {
##### 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)
}
Run the code above in your browser using DataLab