# R/base.R In sensitivity: Global Sensitivity Analysis of Model Outputs

# Base functions of the sensitivity package
# Gilles Pujol, 2006 - 2008
# Modified by Bertrand Iooss.
# Modified by Frank Weber (2016): support model functions
# returning a matrix or a 3-dimensional array.

tell <- function(x, y = NULL, ...)
UseMethod("tell")

extract <- function(x, ...) # fct added in 2020 for sobolshap_knn()
UseMethod("extract")

response <- function(x, loop = FALSE, other_types_allowed = FALSE, ...) {
id <- deparse(substitute(x))

if (class(x\$model) == "function") {
if (loop) {
n <- nrow(x\$X)
y <- sapply(1:n, function(i){
x\$model(x\$X[i,], ...)
}, simplify = "array")
if(is.matrix(y)){
if(typeof(y) == "list"){
stop("The model function returns a list when applied to a row of x\$X")
} else{
# This means the model function returned vectors, so y is a
# (m times n)-matrix (with unknown m). For better consistency with the
# default case (that y is a vector of length n), y is transposed to a
# (n times m)-matrix.
y <- t(y)
}
} else if(is.array(y) && length(dim(y)) == 3){
# This means the model function returned matrices of the same sizes.
# Let m be the number of rows and z be the number of columns.
# For better consistency with the other cases (y is a vector or a
# matrix), y is transformed to an array of dimensions c(n, m, z).

# Change the order of the dimensions:
y <- aperm(y, perm = c(3, 1, 2))
} else if(!is.numeric(y)){
stop("x\$model returned an object that can't be handled")
}
} else {
y <- x\$model(x\$X, ...)
}
} else if (TRUE %in% (paste("predict.", class(x\$model), sep="") %in% methods(predict))) {
y <- predict(x\$model, x\$X, ...)
} else {
stop("The model isn't a function or does not have a predict method")
}

if(other_types_allowed){
#    if (!class(y)[1] %in% c("numeric", "matrix", "array") ||
if ((!inherits(y, "numeric") && !inherits(y, "matrix") && !inherits(y, "array")) ||
(is.array(y) && typeof(y) == "list")) {
y <- as.numeric(y)
warning("Conversion of the response to numeric")
} else if(inherits(y, "array") && length(dim(y)) > 3){
stop("If the model returns an array, it must not have more ",
"than 3 dimensions")
}
} else{
if(!inherits(y, "numeric")) {
y <- as.numeric(y)
warning("Conversion of the response to numeric")
}
}

# Assign column names resp. dimnames if not existing:
if(inherits(y, "matrix") &&
(is.null(colnames(y)) || any(nchar(colnames(y)) == 0))){
colnames(y) <- paste0("ycol", 1:ncol(y))
} else if(inherits(y, "array")){
if(is.null(dimnames(y))){
dimnames(y) <- list(NULL, paste0("ycol", 1:dim(y)[2]),
paste0("ydim3_", 1:dim(y)[3]))
} else{
if(is.null(dimnames(y)[[2]]) || any(nchar(dimnames(y)[[2]]) == 0)){
dimnames(y)[[2]] <- paste0("ycol", 1:dim(y)[2])
}
if(is.null(dimnames(y)[[3]]) || any(nchar(dimnames(y)[[3]]) == 0)){
dimnames(y)[[3]] <- paste0("ydim3_", 1:dim(y)[3])
}
}
}

x\$y <- y
assign(id, x, parent.frame())
}

########################################
# function added in 2014 (for sobolTIIlo() and sobolTIIpf())

plotFG <- function(x)
UseMethod("plotFG")

########################################
# function added in 2016 (tests of PoincareConstant())

# *********************************************************************
# Functions for the Truncated Normal Distribution

# density
dnorm.trunc <- function(x, mean = 0, sd = 1, min = -1e6, max = 1e6){
out = dnorm(x, mean, sd) / (pnorm(max, mean, sd) - pnorm(min, mean, sd))
out[(x < min) | (x > max)] = 0
return(out)
}

# distribution function
pnorm.trunc <- function(q, mean = 0, sd = 1, min = -1e6, max = 1e6){
out = (pnorm(q, mean, sd) - pnorm(min, mean, sd)) /
(pnorm(max, mean, sd) - pnorm(min, mean, sd))
out[q < min] = 0
out[q > max] = 1
return(out)
}

# quantile function
qnorm.trunc <- function(p, mean = 0, sd = 1, min = -1e6, max = 1e6){
return(qnorm((1 - p) * pnorm(min, mean, sd) + p * pnorm(max, mean, sd),
mean, sd))
}

# random generation
rnorm.trunc <- function(n, mean = 0, sd = 1, min = -1e6, max = 1e6){
return(qnorm.trunc(runif(n), mean, sd, min, max))
}

# ***************************************
# Functions for the Truncated Gumbel Distribution

# density
dgumbel.trunc <- function(x, loc = 0, scale = 1, min = -1e6, max = 1e6){
out = evd::dgumbel(x, loc, scale) / (evd::pgumbel(max, loc, scale) - evd::pgumbel(min, loc, scale))
out[(x < min) | (x > max)] = 0
return(out)
}

# distribution function
pgumbel.trunc <- function(q, loc = 0, scale = 1, min = -1e6, max = 1e6){
out = (evd::pgumbel(q, loc, scale) - evd::pgumbel(min, loc, scale)) /
(evd::pgumbel(max, loc, scale) - evd::pgumbel(min, loc, scale))
out[q < min] = 0
out[q > max] = 1
return(out)
}

# quantile function
qgumbel.trunc <- function(p, loc = 0, scale = 1, min = -1e6, max = 1e6){
return(evd::qgumbel((1 - p) * evd::pgumbel(min, loc, scale) + p * evd::pgumbel(max, loc, scale),
loc, scale))
}

# random generation
rgumbel.trunc <- function(n, loc = 0, scale = 1, min = -1e6, max = 1e6){
return(qgumbel.trunc(runif(n), loc, scale, min, max))
}

########################################
# function added in 2017 (for support())

scatterplot <- function(x, i, xprob = FALSE, p = NULL, p.arg = NULL, cex = 1, cex.lab = 1, ...)
UseMethod("scatterplot")

########################################
# function added in 2020 (for sobolMultOut())

plotMultOut <- function(x, ylim = c(0, 1), ...)
UseMethod("plotMultOut")

## Try the sensitivity package in your browser

Any scripts or data that you put into this service are public.

sensitivity documentation built on Sept. 5, 2021, 5:13 p.m.