library(sensitivity)
## Registered S3 method overwritten by 'sensitivity':
##   method    from 
##   print.src dplyr

The fixed parameter

E <- 0.2
S_h <- 0.003
S_V <- 3000

The testing parameter mean

C_m <- 1.5
a_m <- 0.1739
K_m <- 6.7208*10^5
F_z_m <- 9.0034*10^4

The testing parameter sd

C_sd <- 0.1
a_sd <- 0.036
K_sd <- 9.701*10^4
F_z_sd <- 5.5803*10^3

Sobol - generate parameter matrix

n <- 10000
C <- rnorm(n, C_m, C_sd)
a <- rnorm(n, a_m, a_sd)
K <- rnorm(n, K_m, K_sd)
F_z <- rnorm(n, F_z_m, F_z_sd)
X1 <- data.frame(matrix(c(C,a,K,F_z), nrow = n))
C <- rnorm(n, C_m, C_sd)
a <- rnorm(n, a_m, a_sd)
K <- rnorm(n, K_m, K_sd)
F_z <- rnorm(n, F_z_m, F_z_sd)
X2 <- data.frame(matrix(c(C,a,K,F_z), nrow = n))
x <- soboljansen(model = NULL, X1, X2, nboot = 100)

Claculate model output

for( i in seq(nrow(x$X))){
  D <- -3.4*10^-7*x$X[i,4]^2+0.74*x$X[i,4]
  B <- x$X[i,3]/D/x$X[i,1]
  F_y <- D * sin(x$X[i,1]*atan(B*(x$X[i,2]+S_h)-E*(B*(x$X[i,2]+S_h)-atan(B*x$X[i,2]+S_h))))+S_V
  if (i == 1) y <- F_y else y <- c(y, F_y)
}
hist(y)

Sobol - claculate Sobol indices (X1 = C; X2 = a; X3 = K; X4 = F_z)

tell(x,y)
print(x)
## 
## Call:
## soboljansen(model = NULL, X1 = X1, X2 = X2, nboot = 100)
## 
## Model runs: 60000 
## 
## First order indices:
##      original          bias std. error     min. c.i.  max. c.i.
## X1 0.03052948 -0.0025201811 0.01530504 -3.112346e-06 0.06720971
## X2 0.45468790 -0.0032829662 0.01272489  4.332534e-01 0.48427634
## X3 0.18339427 -0.0009877617 0.01644451  1.530832e-01 0.21648349
## X4 0.32118027 -0.0017089329 0.01231835  2.984512e-01 0.34563679
## 
## Total indices:
##      original          bias   std. error  min. c.i.  max. c.i.
## X1 0.01737156  3.419055e-05 0.0004438342 0.01641291 0.01818515
## X2 0.49307431 -1.390719e-04 0.0112479825 0.47084374 0.51534768
## X3 0.20720198  1.653553e-03 0.0050508063 0.19491080 0.21571724
## X4 0.32220133  7.532271e-04 0.0087575178 0.30496359 0.33812582
plot(x)