Nothing
## -----------------------------------------------------------------------------
## Monte Carlo runs
## -----------------------------------------------------------------------------
modCRL <- function(func, parms = NULL, sensvar = NULL, dist = "unif",
parInput = NULL, parRange = NULL, parMean = NULL,
parCovar = NULL, num = 100, ...) {
if (is.null(parms) & ! is.null(parInput))
parms <- parInput[1,]
if (is.null(parms))
parms <- parMean
if (is.null(parms)) {
if (is.vector(parRange))
parms <- mean(parRange)
else parms <- rowMeans(parRange)
}
if (is.null(parms))
stop ("'parms' not known")
if (is.matrix(parms) && nrow(parms)>1)
stop ("'parms' should be a vector")
Solve <- function(parms) func(parms, ...)
if (! is.null(parInput)) {
dist <- "input"
nr <- nrow(parInput)
num <- min(num,nr)
if (num == nr)
ii <- 1:nr
else ii <- sample(1:nr, size = num, replace = FALSE)
parset <- as.matrix(parInput[ii,])
if(is.null(parms)) parms <- parInput[1,]
}
Parms <- parms
## reference run
yRef <- func(parms, ...)
if (is.matrix(yRef) | is.data.frame(yRef))
if(nrow(yRef)>1)
stop("func should return a vector or a matrix/data.frame with one row")
sens <- sensRange(func, parms, sensvar, dist, parInput, parRange, parMean,
parCovar, map = NULL, num = num, ...)
class(sens) <- c("modCRL", "data.frame")
return (sens)
}
## -----------------------------------------------------------------------------
## S3 methods of modCRL
## -----------------------------------------------------------------------------
summary.modCRL <- function(object, ...) {
SumSens <- summary.sensRange(object, ...)
class(SumSens) <- c("summary.modCRL", "data.frame")
return(SumSens)
}
## -----------------------------------------------------------------------------
hist.modCRL <- function(x, which = 1:ncol(x), ask = NULL, ...) {
hh <- list()
hh$pars <- x
hist.modMCMC(hh, which = which, Full = FALSE, ask = ask, ...)
}
## -----------------------------------------------------------------------------
pairs.modCRL <- function(x, which = 1:ncol(x), nsample = NULL, ...) {
panel.main <- function(x, y, ...)
points(x[ii], y[ii], ...)
var <- colnames(x)
which <- selectvar(which, var, "x", Nall = TRUE)
X <- x[,which]
X <- as.matrix(X)
if (is.null(nsample))
ii <- 1:nrow(X) else
ii <- sample((1:nrow(X)), nsample)
labels <- colnames(X)
dots <- list(...)
dots$diag.panel <- if(is.null(dots$diag.panel)) panel.hist else dots$diag.panel
dots$lower.panel <- if(is.null(dots$lower.panel)) panel.cor else dots$lower.panel
dots$upper.panel <- if(is.null(dots$upper.panel)) panel.main else dots$upper.panel
dots$gap <- if(is.null(dots$gap)) 0 else dots$gap
dots$labels <- if(is.null(dots$labels)) labels else dots$labels
do.call("pairs", c(alist(X), dots))
}
## -----------------------------------------------------------------------------
plot.modCRL<-function(x, which = NULL, trace = FALSE, ask = NULL, ...) {
vars <- attr(x,"var")
var <- colnames(x)
NP <- attr(x,"npar")
parnames <- var[1:NP]
varnames <- var[-(1:NP)]
## which selection has to stay the way it is: selecting from TWO sets...
if (!is.null(which)) {
if (! is.numeric(which)) {
ln <- length(which)
Select <- NP + which(varnames %in% which)
SelPar <- which(parnames %in% which)
if(length(Select) + length(SelPar) != ln)
stop("not all variables or parameters in 'which' are in 'x'")
} else {
Select <- which[which > NP]
SelPar <- which[which <= NP]
if (max(which) > length(var))
stop("index in 'which' too large")
}
} else {
Select <- NP+(1:length(varnames))
SelPar <- 1:NP
}
nvar <- length(Select)
np <- length(SelPar)
if (np == 0)
stop("cannot plot Monte Carlo results: select at least one parameter")
if (nvar == 0)
stop("cannot plot Monte Carlo results: select at least one variable")
if (np > 1 & nvar > 1)
stop("cannot plot Monte Carlo results for more than 1 parameter and more than one variable")
pars <- attr(x, "names")
dots <- list(...)
nmdots <- names(dots)
nv <- max(nvar, np)
## Set par mfrow and ask.
ask <- setplotpar (nmdots, dots, nv, ask)
## interactively wait if there are remaining figures
if (ask) {
oask <- devAskNewPage(TRUE)
on.exit(devAskNewPage(oask))
}
Main <- is.null(dots$main)
Xlab <- is.null(dots$xlab)
dots$ylab <- if(is.null(dots$ylab)) "" else dots$ylab
if (np == 1) {
ip <- SelPar[1]
dots$xlab <- if(is.null(dots$xlab)) pars[ip] else dots$xlab
for(i in Select) {
if (Main) dots$main <- var[i]
do.call("plot", c(alist(x[,ip], x[,i]), dots))
if (trace) lines(lowess(x[,ip], x[,i]), lwd = 2)
}
} else { # nvar==1
iv <- Select[1]
for(i in SelPar) {
if (Xlab) dots$xlab <- pars[i]
if (Main) dots$main <- var[iv]
do.call("plot", c(alist(x[,i], x[,iv]), dots))
if (trace) lines(lowess(x[,i], x[,iv]), lwd = 2)
}
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.