# \donttest{
library(parallel)
library(doParallel)
library(foreach)
library(gtools)
library(boot)
library(RANN)
###########################################################
# Linear Model with Gaussian correlated inputs
library(mvtnorm)
set.seed(1234)
n <- 1000
beta<-c(1,-1,0.5)
sigma<-matrix(c(1,0,0,
0,1,-0.8,
0,-0.8,1),
nrow=3,
ncol=3)
X <-rmvnorm(n, rep(0,3), sigma)
colnames(X)<-c("X1","X2", "X3")
y <- X%*%beta + rnorm(n,0,2)
# Without Bootstrap confidence intervals
x<-shapleysobol_knn(model=NULL, X=X,
n.knn=3,
noise=TRUE)
tell(x,y)
print(x)
plot(x)
#Using the extract method to get first-order and total Sobol' indices
extract(x)
# With Boostrap confidence intervals
x<-shapleysobol_knn(model=NULL, X=X,
nboot=10,
n.knn=3,
noise=TRUE,
boot.level=0.7,
conf=0.95)
tell(x,y)
print(x)
plot(x)
#####################
# Extracting Sobol' indices with Bootstrap confidence intervals
nboot <- 10 # put nboot=50 for consistency
#Total Sobol' indices
x<-shapleysobol_knn(model=NULL, X=X,
nboot=nboot,
n.knn=3,
U=0,
noise=TRUE,
boot.level=0.7,
conf=0.95)
tell(x,y)
print(x)
plot(x)
#First-order Sobol' indices
x<-shapleysobol_knn(model=NULL, X=X,
nboot=nboot,
n.knn=3,
U=1,
noise=TRUE,
boot.level=0.7,
conf=0.95)
tell(x,y)
print(x)
plot(x)
#Closed Sobol' indices for specific subsets (list)
x<-shapleysobol_knn(model=NULL, X=X,
nboot=nboot,
n.knn=3,
U=list(c(1,2), c(1,2,3), 2),
noise=TRUE,
boot.level=0.7,
conf=0.95)
tell(x,y)
print(x)
plot(x)
#####################################################
# Test case: the non-monotonic Sobol g-function
# Example with a call to a numerical model
# First compute first-order indices with ranking
n <- 1000
X <- data.frame(matrix(runif(8 * n), nrow = n))
x <- shapleysobol_knn(model = sobol.fun, X = X, U = 1, method = "rank")
print(x)
plot(x)
library(ggplot2) ; ggplot(x)
# We can use the output sample generated for this estimation to compute total indices
# without additional calls to the model
x2 <- shapleysobol_knn(model = NULL, X = X, U = 0, method = "knn", n.knn = 5)
tell(x2,x$y)
plot(x2)
ggplot(x2)
#####################################################
# Test case: the Ishigami function
# Example with given data and the use of approximate nearest neighbour search
n <- 5000
X <- data.frame(matrix(-pi+2*pi*runif(3 * n), nrow = n))
Y <- ishigami.fun(X)
x <- shapleysobol_knn(model = NULL, X = X, U = NULL, method = "knn", n.knn = 5,
n.limit = 2000)
tell(x,Y)
plot(x)
library(ggplot2) ; ggplot(x)
# Extract first-order and total Sobol indices
x1 <- extract(x) ; print(x1)
######################################################
# Test case : Linear model (3 Gaussian inputs including 2 dependent) with scaling
# See Iooss and Prieur (2019)
library(mvtnorm) # Multivariate Gaussian variables
library(whitening) # For scaling
modlin <- function(X) apply(X,1,sum)
d <- 3
n <- 10000
mu <- rep(0,d)
sig <- c(1,1,2)
ro <- 0.9
Cormat <- matrix(c(1,0,0,0,1,ro,0,ro,1),d,d)
Covmat <- ( sig %*% t(sig) ) * Cormat
Xall <- function(n) mvtnorm::rmvnorm(n,mu,Covmat)
X <- Xall(n)
x <- shapleysobol_knn(model = modlin, X = X, U = NULL, method = "knn", n.knn = 5,
rescale = TRUE, n.limit = 2000)
print(x)
plot(x)
# }
Run the code above in your browser using DataLab