mcs: Monte-Carlo simulation
In dkneis/mcu: Model calibration utilities

Description Usage Arguments Value Note Author(s) Examples

Performs a Monte-Carlo simulation with parameters sampled from a uniform distribution using a latin hypercube method.

1	mcs(fn, p, nRuns = 10, silent = TRUE, parallel = FALSE, ...)

`fn`	Function of interest. It must accept as its first argument a named numeric vector of parameters. It must return a named numeric vector (which can be of length 1).
`p`	Data frame specifying ranges and defaults for the varied parameters. There must be four columns: 'name', 'default', 'min', and 'max'.
`nRuns`	Desired number of parameter samples. The total number of evaluations of `fn` is `nRuns` + 1.
`silent`	If `TRUE`, diagnostic messages are suppressed.
`parallel`	If `TRUE`, parameter sets are processed in parallel if supported by the machine. The number of threads is controlled automatically (by calling `registerDoParallel()` with no arguments).
`...`	Additional arguments passed to function `fn`.

A list of 3 elements. Elements p and out are matrices with nRuns + 1 rows. p holds the tested parameter values with column names taken from the 'name' field of ranges. out holds the return values of fn. Each row in out corresponds to the same row of p. In the common case where fn returns a scalar result, out contains just a single column. The first row in both data frames corresponds to the default parameter set. The third list element, cpu is a vector of length nRuns + 1, holding the times spent on the evaluation of fn.

If fn generated an error (or a warning) when called with a parameter set, the corresponding row in the result matrix out. is set to NA. The same is true for the corresponding element of the returned vector cpu.

David Kneis david.kneis@tu-dresden.de

# Analysis of the residuals' sum of squares for a linear model
obs= data.frame(x=c(1,2), y=c(1,2))
model= function(p, x) { p["slope"] * x + p["intercept"] }
objfun= function(p, obs) { c(sse= sum((obs$y - model(p, obs$x))^2)) }
p= data.frame(
  name=c("slope","intercept"),
  default= c(1, 0),
  min= c(0.5, -1),
  max= c(2, 1)
)
x= mcs(fn=objfun, p=p, obs=obs)
layout(matrix(1:2, ncol=2))
plot(x$p[,"slope"], x$out[,"sse"], xlab="slope", ylab="SSE")
plot(x$p[,"intercept"], x$out[,"sse"], xlab="intercept", ylab="SSE")
layout(matrix(1))