How to create the sampling matrixes for Sobol sensitivity analysis in R (package "sensitivity")

Viewed 141

I would like to perform a Sobol sensitivity analysis in R

The package "sensitivity" should allow me to do so, but I don't understand how to generate the sampling matrixes (X1, X2). I have a model that runs outside of R. I have 6 parameters with uniform distribution.

In my text book: N = (2k+2)*M ; M = 2^b ; b=[8,12] (New sampling method : Wu et al. 2012)

I had the feeling that I should create two sampling matrix and feed the two to the sobol function X1_{M,k} X2_{M,k}.

The dimension of final sampling matrix x$X is then (k+2)*M. because:

  X <- rbind(X1, X2)
  for (i in 1:k) {
    Xb <- X1
    Xb[, i] <- X2[, i]
    X <- rbind(X, Xb)
  }

How should I conduct my sampling to get the right number of runs as (2*k+2)*M ?

This script is for the old method but does someone know if the new method is already implemented yet in the sensitivity package? Feel free to comment this procedure

name = c("a" , "b" , "c" , "d" , "e", "f")
vals <- list(list(var="a",dist="unif",params=list(min=0.1,max=1.5)),
             list(var="b",dist="unif",params=list(min=-0.3,max=0.4)),
             list(var="c",dist="unif",params=list(min=-0.3,max=0.3)),
             list(var="d",dist="unif",params=list(min=0,max=0.5)),
             list(var="e",dist="unif",params=list(min=2.4E-5,max=2.4E-3)),
             list(var="f",dist="unif",params=list(min=3E-5,max=3E-3)))
k = 6
b = 8
M = 2^b
n <- 2*M
X1 <- makeMCSample(n,vals, p = 1)
X2 <- makeMCSample(n,vals, p = 2)

x <- sobol2007(model = NULL, X1, X2, nboot = 200)

if I understand correctly, I should provide a y for each x$X sampling combination

then I can use the function "tell" which will generate the Sobol' first-order indices as well as the total indices

tell(x,y)
ggplot(x)

Supplemental R function SobolR

makeMCSample <- function(n, vals) {
  # Packages to generate quasi-random sequences
  # and rearrange the data
  require(randtoolbox)
  require(plyr)
  
  # Generate a Sobol' sequence
  if (p == 2){ sob <- sobol(n, length(vals), seed = 4321, scrambling = 1)
  }else{sob <- sobol(n, length(vals), seed = 1234, scrambling = 1)}
  
  # Fill a matrix with the values
  # inverted from uniform values to
  # distributions of choice
  samp <- matrix(rep(0,n*(length(vals)+1)), nrow=n)
  samp[,1] <- 1:n
  for (i in 1:length(vals)) {
    # i=1
    l <- vals[[i]]
    dist <- l$dist
    params <- l$params
    fname <- paste("q",dist,sep="")
    samp[,i+1] <- do.call(fname,c(list(p=sob[,i]),params))
  }
  
  # Convert matrix to data frame and add labels
  samp <- as.data.frame(samp)
  names(samp) <- c("n",laply(vals, function(l) l$var))
  return(samp)
}

ref: Qiong-Li Wu, Paul-Henry Cournède, Amélie Mathieu, 2012, Efficient computational method for global sensitivity analysis and its application to tree growth modelling

0 Answers
Related