Automated Calibration for Complex Models

Description

Automated Calibration for Complex Models

Details

calibrar package: Automated Calibration for Complex Models

This package allows the parameter estimation (i.e. calibration) of complex models, including stochastic ones. It implements generic functions that can be used for fitting any type of models, especially those with non-differentiable objective functions, with the same syntax as base::optim. It supports multiple phases estimation (sequential parameter masking), constrained optimization (bounding box restrictions) and automatic parallel computation of numerical gradients. Some common maximum likelihood estimation methods and automated construction of the objective function from simulated model outputs is provided. See <https://roliveros-ramos.github.io/calibrar/> for more details.

Author(s)

Ricardo Oliveros-Ramos Maintainer: Ricardo Oliveros-Ramos <ricardo.oliveros@gmail.com>

References

calibrar: an R package for the calibration of ecological models (Oliveros-Ramos and Shin 2014)

Examples

## Not run: 
require(calibrar)
set.seed(880820)
path = NULL # NULL to use the current directory
# create the demonstration files
demo = calibrar_demo(model="PoissonMixedModel", L=5, T=100) 
# get calibration information
calibrationInfo = calibration_setup(file=demo$path)
# get observed data
observed = calibration_data(setup=calibrationInfo, path=demo$path)
# read forcings for the model
forcing = read.csv(file.path(demo$path, "master", "environment.csv"), row.names=1)
# Defining 'runModel' function
runModel = function(par, forcing) {
output = calibrar:::.PoissonMixedModel(par=par, forcing=forcing)
# adding gamma parameters for penalties
output = c(output, list(gammas=par$gamma)) 
return(output)
}
# real parameters
cat("Real parameters used to simulate data\n")
print(demo$par)
# objective functions
obj  = calibration_objFn(model=runModel, setup=calibrationInfo, 
                               observed=observed, forcing=forcing)
cat("Starting calibration...\n")
control = list(weights=calibrationInfo$weights, maxit=3.6e5) # control parameters
cat("Running optimization algorithms\n", "\t", date(), "\n")
cat("Running optim AHR-ES\n")
ahr = calibrate(par=demo$guess, fn=obj, lower=demo$lower, upper=demo$upper, control=control)
summary(ahr)

## End(Not run) 

Get an specific argument from the command line

Description

Get an specific argument from the command line

Usage

.get_command_argument(
  x,
  argument,
  prefix = "--",
  default = FALSE,
  verbose = FALSE
)

Arguments

Value

The value of the argument, assumed to be followed after '=' or, TRUE if nothing but the argument was found. If the argument is not found, FALSE is returned.

Examples

.get_command_argument(commandArgs(), "interactive")
.get_command_argument(commandArgs(), "RStudio")
.get_command_argument(commandArgs(), "RStudio", prefix="")
.get_command_argument(commandArgs(), "vanilla")
.get_command_argument("--control.file=baz.txt", "control.file")

Read a configuration file.

Description

File is expected to have lines of the form 'key SEP value' where key is the name of the parameter, SEP a separator (can be '=' ',', ';') and value the value of the parameter itself. The SEP for each line is determined and parameters values are returned as a list.

Usage

.read_configuration(
  file,
  recursive = TRUE,
  keep.names = TRUE,
  conf.key = NULL,
  ...
)

Arguments

Adaptative Hierarchical Recombination Evolutionary Strategy (AHR-ES) for derivative-free and black-box optimization

Description

This function performs the optimization of a function using the Adaptative Hierarchical Recombination Evolutionary Strategy (AHR-ES, Oliveros & Shin, 2015).

Usage

ahres(
  par,
  fn,
  gr = NULL,
  ...,
  lower = -Inf,
  upper = +Inf,
  active = NULL,
  control = list(),
  hessian = FALSE,
  parallel = FALSE
)

Arguments

Value

A list with components:

par

The best set of parameters found.

value

The value of fn corresponding to par.

counts

A two-element integer vector giving the number of calls to fn and gr respectively. This excludes those calls needed to compute the Hessian, if requested, and any calls to fn to compute a finite-difference approximation to the gradient.

convergence

An integer code. 0 indicates successful completion.

message

A character string giving any additional information returned by the optimizer, or NULL.

hessian

Only if argument hessian is true. A symmetric matrix giving an estimate of the Hessian at the solution found. Note that this is the Hessian of the unconstrained problem even if the box constraints are active.

Author(s)

Ricardo Oliveros-Ramos

See Also

Other optimisers: calibrate(), optim2(), optimh()

Examples

## Not run: ahres(par=rep(1, 5), fn=sphereN)

Defunct functions in package calibrar.

Description

The functions listed below are defunct. When possible, alternative functions with similar functionality are also mentioned. Help pages for defunct functions are available at help("<function>-defunct").

Usage

getObservedData(info, path, data.folder = "data", ...)

getCalibrationInfo(
  path,
  file = "calibrationInfo.csv",
  stringsAsFactors = FALSE,
  ...
)

createObjectiveFunction(
  runModel,
  info,
  observed,
  aggFn = .weighted.sum,
  aggregate = FALSE,
  ...
)

getObservedData

Deprecated in v0.3, defunct in v0.9. For getObservedData, use calibration_data.

getCalibrationInfo

Deprecated in v0.3, defunct in v0.9. For getCalibrationInfo, use calibration_setup.

createObjectiveFunction

Deprecated in v0.3, defunct in v0.9. For createObjectiveFunction, use calibration_objFn.

Demos for the calibrar package

Description

Creates demo files able to be processed for a full calibration using the calibrar package

Usage

calibrar_demo(path = NULL, model = NULL, ...)

Arguments

Details

Current implemented models are:

PoissonMixedModel

Poisson Autoregressive Mixed model for the dynamics of a population in different sites:

log(\mu_{i, t+1}) = log(\mu_{i, t}) + \alpha + \beta X_{i, t} + \gamma_t

where \mu_{i, t} is the size of the population in site i at year t, X_{i, t} is the value of an environmental variable in site i at year t. The parameters to estimate were \alpha, \beta, and \gamma_t, the random effects for each year, \gamma_t \sim N(0,\sigma^2), and the initial population at each site \mu_{i, 0}. We assumed that the observations N_{i,t} follow a Poisson distribution with mean \mu_{i, t}.

PredatorPrey

Lotka Volterra Predator-Prey model. The model is defined by a system of ordinary differential equations for the abundance of prey $N$ and predator $P$:

\frac{dN}{dt} = rN(1-N/K)-\alpha NP

\frac{dP}{dt} = -lP + \gamma\alpha NP

The parameters to estimate are the prey’s growth rate r, the predator’s mortality rate l, the carrying capacity of the prey K and \alpha and \gamma for the predation interaction. Uses deSolve package for numerical solution of the ODE system.

SIR

Susceptible-Infected-Recovered epidemiological model. The model is defined by a system of ordinary differential equations for the number of susceptible $S$, infected $I$ and recovered $R$ individuals:

\frac{dS}{dt} = -\beta S I/N

\frac{dI}{dt} = \beta S I/N -\gamma I

\frac{dR}{dt} = \gamma I

The parameters to estimate are the average number of contacts per person per time \beta and the instant probability of an infectious individual recovering \gamma. Uses deSolve package for numerical solution of the ODE system.

IBMLotkaVolterra

Stochastic Individual Based Model for Lotka-Volterra model. Uses ibm package for the simulation.

Value

A list with the following elements:

Author(s)

Ricardo Oliveros–Ramos

References

Oliveros-Ramos and Shin (2014)

Examples

## Not run: 

summary(ahr)
set.seed(880820)
path = NULL # NULL to use the current directory
# create the demonstration files
demo = calibrar_demo(path=path, model="PredatorPrey", T=100) 
# get calibration information
calibration_settings = calibration_setup(file = demo$setup)
# get observed data
observed = calibration_data(setup = calibration_settings, path=demo$path)
# Defining 'run_model' function
run_model = calibrar:::.PredatorPreyModel
# real parameters
cat("Real parameters used to simulate data\n")
print(unlist(demo$par)) # parameters are in a list
# objective functions
obj  = calibration_objFn(model=run_model, setup=calibration_settings, observed=observed, T=demo$T)
obj2 = calibration_objFn(model=run_model, setup=calibration_settings, observed=observed, 
T=demo$T, aggregate=TRUE)
cat("Starting calibration...\n")
cat("Running optimization algorithms\n", "\t")
cat("Running optim AHR-ES\n")
ahr = calibrate(par=demo$guess, fn=obj, lower=demo$lower, upper=demo$upper, phases=demo$phase)
summary(ahr)

## End(Not run) 

Sequential parameter estimation for the calibration of complex models

Description

This function performs the optimization of a function, possibly in sequential phases of increasing complexity, and it is designed for the calibration of a model, by minimizing the error function fn associated to it.

Usage

calibrate(
  par,
  fn,
  gr,
  ...,
  method,
  lower,
  upper,
  phases,
  control,
  hessian,
  replicates,
  parallel
)

## Default S3 method:
calibrate(
  par,
  fn,
  gr = NULL,
  ...,
  method = NULL,
  lower = NULL,
  upper = NULL,
  phases = NULL,
  control = list(),
  hessian = FALSE,
  replicates = 1,
  parallel = FALSE
)

Arguments

Details

In the control list, aggFn is a function to aggregate fn to a scalar value if the returned value is a vector. Some optimization algorithm can exploite the additional information provided by a vectorial output from fn.

Author(s)

Ricardo Oliveros-Ramos

See Also

Other optimisers: ahres(), optim2(), optimh()

Examples

calibrate(par=rep(NA, 5), fn=sphereN)
## Not run: 
calibrate(par=rep(NA, 5), fn=sphereN, replicates=3)
calibrate(par=rep(0.5, 5), fn=sphereN, replicates=3, lower=-5, upper=5)
calibrate(par=rep(0.5, 5), fn=sphereN, replicates=3, lower=-5, upper=5, phases=c(1,1,1,2,3))
calibrate(par=rep(0.5, 5), fn=sphereN, replicates=c(1,1,4), lower=-5, upper=5, phases=c(1,1,1,2,3))

## End(Not run)

Get observed data for the calibration of a model

Description

Create a list with the observed data with the information provided by its main argument.

Usage

calibration_data(setup, path = ".", file = NULL, verbose = TRUE, ...)

Arguments

Value

A list with the observed data needed for a calibration, to be used in combination with the calibration_objFn.

Author(s)

Ricardo Oliveros-Ramos

See Also

calibration_objFn, calibration_setup.

Create an objective function to be used with optimization routines

Description

Create a new function, to be used as the objective function in the calibration, given a function to run the model within R, observed data and information about the comparison with data.

Usage

calibration_objFn(model, setup, observed, aggFn = NULL, aggregate = FALSE, ...)

Arguments

Value

A function, integrating the simulation of the model and the comparison with observed data.

Author(s)

Ricardo Oliveros-Ramos

See Also

calibration_data, calibration_setup.

Get information to run a calibration using the calibrar package.

Description

A wrapper for read.csv checking column names and data types for the table with the calibration information.

Usage

calibration_setup(file, control = list(), ...)

Arguments

Value

A data.frame with the information for the calibration of a model, to be used with the calibration_objFn and calibration_data.

Author(s)

Ricardo Oliveros-Ramos

See Also

calibration_objFn, calibration_data.

Create an objective function to be used with optimization routines

Description

Create a new function, to be used as the objective function in the calibration, given a function to run the model within R, observed data and information about the comparison with data.

Arguments

Value

A function, integrating the simulation of the model and the comparison with observed data.

Author(s)

Ricardo Oliveros-Ramos

See Also

calibrar-defunct

Calculate a discretization of the 2D Gaussian Kernel

Description

Calculate a discretization of the 2D Gaussian Kernel

Usage

gaussian_kernel(par, lower, upper, n = 10, checkSymmetry = TRUE, ...)

Arguments

Value

A list, with 'x', 'y' and 'z' components.

Get information to run a calibration using the calibrar package.

Description

A wrapper for read.csv checking column names and data types for the table with the calibration information.

Arguments

Value

A data.frame with the information for the calibration of a model, to be used with the createObjectiveFunction and getObservedData.

Author(s)

Ricardo Oliveros-Ramos

See Also

calibrar-defunct

Get observed data for the calibration of a model

Description

Create a list with the observed data with the information provided by its main argument.

Arguments

Value

A list with the observed data needed for a calibration, to be used in combination with the createObjectiveFunction.

Author(s)

Ricardo Oliveros-Ramos

See Also

calibrar-defunct

Numerical computation of the gradient, with parallel capabilities

Description

This function calculates the gradient of a function, numerically, including the possibility of doing it in parallel.

Usage

gradient(fn, x, method, control, parallel, ...)

Arguments

Value

The gradient of fn at x.

Examples

gradient(fn=function(x) sum(x^3), x=0)

Calcuted error measure between observed and simulated data

Description

Calcuted error measure between observed and simulated data

Usage

objFn(obs, sim, FUN, ...)

fitness(obs, sim, FUN, ...)

Arguments

Value

the value of FUN(obs, sim, ...)

General-purpose optimization with parallel numerical gradient computation

Description

General-purpose optimization with parallel numerical gradient computation

Usage

optim2(
  par,
  fn,
  gr = NULL,
  ...,
  method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN", "Brent", "nlm", "nlminb",
    "Rcgmin", "Rvmmin", "hjn", "spg", "LBFGSB3", "AHR-ES"),
  lower = -Inf,
  upper = +Inf,
  active = NULL,
  control = list(),
  hessian = FALSE,
  parallel = FALSE
)

Arguments

Value

A list with components:

par

The best set of parameters found.

value

The value of fn corresponding to par.

counts

A two-element integer vector giving the number of calls to fn and gr respectively. This excludes those calls needed to compute the Hessian, if requested, and any calls to fn to compute a finite-difference approximation to the gradient.

convergence

An integer code. 0 indicates successful completion.

message

A character string giving any additional information returned by the optimizer, or NULL.

hessian

Only if argument hessian is true. A symmetric matrix giving an estimate of the Hessian at the solution found. Note that this is the Hessian of the unconstrained problem even if the box constraints are active.

Author(s)

Ricardo Oliveros-Ramos

See Also

Other optimisers: ahres(), calibrate(), optimh()

Examples

optim2(par=rep(NA, 5), fn=sphereN)

General-purpose optimization using heuristic algorithms

Description

General-purpose optimization using heuristic algorithms

Usage

optimh(
  par,
  fn,
  gr = NULL,
  ...,
  method = c("AHR-ES", "Nelder-Mead", "SANN", "hjn", "bobyqa", "CMA-ES", "genSA", "DE",
    "soma", "genoud", "PSO", "hybridPSO", "mads", "hjk", "hjkb", "nmk", "nmkb"),
  lower = -Inf,
  upper = +Inf,
  active = NULL,
  control = list(),
  hessian = FALSE,
  parallel = FALSE
)

Arguments

Value

A list with components:

par

The best set of parameters found.

value

The value of fn corresponding to par.

counts

A two-element integer vector giving the number of calls to fn and gr respectively. This excludes those calls needed to compute the Hessian, if requested, and any calls to fn to compute a finite-difference approximation to the gradient.

convergence

An integer code. 0 indicates successful completion.

message

A character string giving any additional information returned by the optimizer, or NULL.

hessian

Only if argument hessian is true. A symmetric matrix giving an estimate of the Hessian at the solution found. Note that this is the Hessian of the unconstrained problem even if the box constraints are active.

Author(s)

Ricardo Oliveros-Ramos

See Also

Other optimisers: ahres(), calibrate(), optim2()

Examples

optim2(par=rep(NA, 5), fn=sphereN)

Sphere function with random noise

Description

This function calculates the Euclidian distance from a point to the origin after a random displacement of its position.

Usage

sphereN(x, sd = 0.1, aggregate = TRUE)

Arguments

Value

The distance from the point x to the origin after a random displacement.

Author(s)

Ricardo Oliveros–Ramos

Examples

sphereN(rep(0, 10))

Predict time-varying parameters using splines.

Description

Predict time-varying parameters using splines.

Usage

spline_par(par, n, knots = NULL, periodic = FALSE, period = NULL)

Arguments

Value

A list with the interpolates values as 'x' and 'time'.

Summary for calibration results object

Description

Summary for calibration results object

Usage

## S3 method for class 'calibrar.results'
summary(object, ..., show_par = NULL, par.only = FALSE)

Arguments