FME (version 1.3.6.3)

modCRL: Monte Carlo Analysis

Description

Given a model consisting of differential equations, estimates the global effect of certain (sensitivity) parameters on selected sensitivity variables.

This is done by drawing parameter values according to some predefined distribution, running the model with each of these parameter combinations, and calculating the values of the selected output variables at each output interval.

This function is useful for ``what-if'' scenarios.

If the output variables consist of a time-series or spatially dependent, use sensRange instead.

Usage

modCRL(func, parms = NULL, sensvar = NULL, dist = "unif",
       parInput = NULL, parRange = NULL, parMean = NULL, parCovar = NULL,
       num = 100, ...)
# S3 method for modCRL
summary(object,  ...)
# S3 method for modCRL
plot(x, which = NULL, trace = FALSE, ask = NULL, ...)
# S3 method for modCRL
pairs(x, which = 1:ncol(x), nsample = NULL, ...)
# S3 method for modCRL
hist(x, which = 1:ncol(x), ask = NULL, ...)

Value

a data.frame of type modCRL containing the parameter(s) and the corresponding values of the sensitivity output variables.

The list returned by modCRL has a method for the generic functions

summary and plot -- see note.

Arguments

func

an R-function that has as first argument parms and that returns a vector with variables whose sensitivity should be estimated.

parms

parameters passed to func; should be either a vector, or a list with named elements. If NULL, then the first element of parInput is taken.

sensvar

the output variables for which the sensitivity needs to be estimated. Either NULL, the default=all output variables, or a vector with output variable names (which should be present in the vector returned by func), or a vector with indices to output variables as present in the output vector returned by func.

dist

the distribution according to which the parameters should be generated, one of "unif" (uniformly random samples), "norm", (normally distributed random samples), "latin" (latin hypercube distribution), "grid" (parameters arranged on a grid).

The input parameters for the distribution are specified by parRange (min,max), except for the normally distributed parameters, in which case the distribution is specified by the parameter means parMean and the variance-covariance matrix, parCovar. Note that, if the distribution is "norm" and parRange is given, then a truncated distribution will be generated. (This is useful to prevent for instance that certain parameters become negative). Ignored if parInput is specified.

parRange

the range (min, max) of the sensitivity parameters, a matrix or (preferred) a data.frame with one row for each parameter, and two columns with the minimum (1st) and maximum (2nd) value. The rownames of parRange should be parameter names that are known in argument parms. Ignored if parInput is specified.

parInput

a matrix with dimension (*, npar) with the values of the sensitivity parameters.

parMean

only when dist is "norm": the mean value of each parameter. Ignored if parInput is specified.

parCovar

only when dist is "norm": the parameter's variance-covariance matrix.

num

the number of times the model has to be run. Set large enough. If parInput is specified, then num parameters are selected randomly (from the rows of parInput.

object

an object of class modCRL.

x

an object of class modCRL.

which

the name or the index to the variables and parameters that should be plotted. Default = all variables and parameters.

nsample

the number of xy pairs to be plotted on the upper panel in the pairs plot. When NULL all xy pairs plotted. Set to a lower number in case the graph becomes too dense (and the exported picture too large). This does not affect the histograms on the diagonal plot (which are estimated using all values).

trace

if TRUE, adds smoothed line to the plot.

ask

logical; if TRUE, the user is asked before each plot, if NULL the user is only asked if more than one page of plots is necessary and the current graphics device is set interactive, see par(ask=.) and dev.interactive.

...

additional arguments passed to function func or to the methods.

Author

Karline Soetaert <[e-mail]>.

References

Soetaert, K. and Petzoldt, T. 2010. Inverse Modelling, Sensitivity and Monte Carlo Analysis in R Using Package FME. Journal of Statistical Software 33(3) 1--28. tools:::Rd_expr_doi("10.18637/jss.v033.i03")

Examples

## =======================================================================
## Bacterial growth model as in Soetaert and Herman, 2009
## =======================================================================

pars <- list(gmax = 0.5,eff = 0.5,
              ks = 0.5, rB = 0.01, dB = 0.01)

solveBact <- function(pars) {
  derivs <- function(t,state,pars) {    # returns rate of change
    with (as.list(c(state,pars)), {
      dBact <-  gmax*eff * Sub/(Sub + ks)*Bact - dB*Bact - rB*Bact
      dSub  <- -gmax     * Sub/(Sub + ks)*Bact + dB*Bact
      return(list(c(dBact, dSub)))
    })
  }

 state <- c(Bact = 0.1, Sub = 100)
 tout  <- seq(0, 50, by = 0.5)
 ## ode solves the model by integration...
 return(as.data.frame(ode(y = state, times = tout, func = derivs,
                          parms = pars)))
}

out <- solveBact(pars)

plot(out$time, out$Bact, main = "Bacteria",
  xlab = "time, hour", ylab = "molC/m3", type = "l", lwd = 2)

## Function that returns the last value of the simulation
SF <- function (p) {
  pars["eff"] <- p
  out <- solveBact(pars)
  return(out[nrow(out), 2:3])
}

parRange <- matrix(nr = 1, nc = 2, c(0.2, 0.8),
  dimnames = list("eff", c("min", "max")))
parRange

CRL <- modCRL(func = SF, parRange = parRange)

plot(CRL)  # plots both variables
plot(CRL, which = c("eff", "Bact"), trace = FALSE) #selects one