Introduction

This vignette presents a simple stochastic process model in which the quantity of interest is the time needed to reach a given number of successes. The model is defined at the level of individual units and is driven by exponential and Bernoulli random variables.

The goal is to compute Sobol indices for this model when the distributional parameters are uncertain.

Process model

The elementary unit is implemented by sobol4r_one_unit. The individual level model sobol4r_process_indiv aggregates units up to a target number of successes.

Sobol4R:::one_unit(
  lambda1 = 1 / 60,
  lambda2 = 1 / 1012,
  lambda3 = 1 / 72,
  p1      = 0.18,
  p2      = 0.5
)

process_fun_indiv(
  X_indiv = c(
    lambda1 = 1 / 60,
    lambda2 = 1 / 1012,
    lambda3 = 1 / 72,
    p1      = 0.18,
    p2      = 0.5
  ),
  M = 50
)

Design for distributional parameters

We build two designs X1 and X2 for the uncertain distributional parameters, which are interpreted row wise by the process model.

n <- 200

draw_params <- function(n) {
  data.frame(t(replicate(
    n,
    c(
      1 / runif(1, min = 20,  max = 100),
      1 / runif(1, min = 24,  max = 2000),
      1 / runif(1, min = 24,  max = 120),
      runif(1,  min = 0.05, max = 0.3),
      runif(1,  min = 0.3,  max = 0.7)
    )
  )))
}

X1 <- draw_params(n)
X2 <- draw_params(n)

gensol_proc <- sobol4r_design(
  X1    = X1,
  X2    = X2,
  order = 2,
  nboot = 100
)
gensol_proc_s2007 <- sensitivity::sobol2007(
  model=NULL,
  X1    = X1,
  X2    = X2,
  nboot = 100
)

Sobol indices based on a single trajectory

Yproc_1 <- process_fun_row_wise(gensol_proc$X, M = 50)

Yproc_s2007 <- process_fun_row_wise(gensol_proc_s2007$X, M = 50)

xproc_1 <- tell(gensol_proc, Yproc_1)
xproc_s2007 <- tell(gensol_proc_s2007, Yproc_s2007)

print(xproc_1)
Sobol4R::autoplot(xproc_1, ncol = 1)

print(xproc_s2007)
Sobol4R::autoplot(xproc_s2007)

Sobol indices based on a quantity of interest

Instead of relying on a single trajectory, we can define a quantity of interest for each parameter vector, for instance the mean time to reach the target number of successes. This is implemented by sobol4r_process_qoi.

Yproc_2 <- process_fun_mean_to_M(gensol_proc$X, M = 50)

xproc_2 <- tell(gensol_proc, Yproc_2)

print(xproc_2)
Sobol4R::autoplot(xproc_2, ncol = 1)

Qoi Mean

res_sobol_mean <- sobol4r_qoi_indices(
  model = process_fun_row_wise,
  X1      = X1,
  X2      = X2,
  qoi_fun = base::mean,
  nrep    = 1000, 
  order   = 2,
  nboot   = 20,
  M       = 50,
  type  = "sobol"
)
print(res_sobol_mean)
Sobol4R::autoplot(res_sobol_mean, ncol = 1)

res_sobol2007_mean <- sobol4r_qoi_indices(
  model = process_fun_row_wise,
  X1      = X1,
  X2      = X2,
  qoi_fun = base::mean,
  nrep    = 1000, 
  nboot   = 20,
  M       = 50,
  type  = "sobol2007"
)
print(res_sobol2007_mean)
Sobol4R::autoplot(res_sobol2007_mean)

Qoi Median

res_sobol_median <- sobol4r_qoi_indices(
  model = process_fun_row_wise,
  X1      = X1,
  X2      = X2,
  qoi_fun = stats::median,
  nrep    = 1000, 
  order   = 2,
  nboot   = 20,
  M       = 50,
  type    = "sobol"
)
print(res_sobol_median)
Sobol4R::autoplot(res_sobol_median, ncol = 1)

res_sobol2007_median <- sobol4r_qoi_indices(
  model = process_fun_row_wise,
  X1      = X1,
  X2      = X2,
  qoi_fun = stats::median,
  nrep    = 1000, 
  nboot   = 20,
  M       = 50,
  type    = "sobol2007"
)
print(res_sobol2007_median)
Sobol4R::autoplot(res_sobol2007_median)

Summary

This vignette illustrates how Sobol4R can be used to compute Sobol indices for models in which both the mechanisms and their distributional parameters are stochastic. The same workflow can be applied to more complex simulators and to various quantities of interest (QoI).

Sobol4R: Sobol indices for a stochastic process model