---
title: "Sobol4R: Sobol indices for a stochastic process model"
shorttitle: "Sobol4R: Sobol indices for a stochastic process model"
author:
- name: "Frédéric Bertrand"
  affiliation:
  - Cedric, Cnam, Paris
  email: frederic.bertrand@lecnam.net
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{Sobol4R: Sobol indices for a stochastic process model}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "figures/sobol-process-",
  fig.width = 7,
  fig.height = 5,
  dpi = 150,
  message = FALSE,
  warning = FALSE,
  eval=FALSE
)

LOCAL <- identical(Sys.getenv("LOCAL"), "TRUE")

library(sensitivity)
library(Sobol4R)
set.seed(4669)
```

# 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.

```{r process-demo, cache=TRUE, eval=LOCAL}
Sobol4R:::one_unit(
  lambda1 = 1 / 60,
  lambda2 = 1 / 1012,
  lambda3 = 1 / 72,
  p1      = 0.18,
  p2      = 0.5
)
```

```{r process-indiv-demo, cache=TRUE, eval=LOCAL}
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.

```{r process-design, cache=TRUE, eval=LOCAL}
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)
```

```{r, cache=TRUE, eval=LOCAL}
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

```{r process-sobol-single, cache=TRUE, eval=LOCAL}
Yproc_1 <- process_fun_row_wise(gensol_proc$X, M = 50)
```


```{r process-sobol-single2, cache=TRUE, eval=LOCAL}
Yproc_s2007 <- process_fun_row_wise(gensol_proc_s2007$X, M = 50)
```

```{r process-sobol-tell, cache=TRUE, eval=LOCAL}
xproc_1 <- tell(gensol_proc, Yproc_1)
xproc_s2007 <- tell(gensol_proc_s2007, Yproc_s2007)
```

```{r process-sobol-single-plot, cache=TRUE, eval=LOCAL}
print(xproc_1)
Sobol4R::autoplot(xproc_1, ncol = 1)
```

```{r process-sobol-single-plot2, cache=TRUE, eval=LOCAL}
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`.

```{r process-sobol-qoi, eval=LOCAL, cache=TRUE}
Yproc_2 <- process_fun_mean_to_M(gensol_proc$X, M = 50)
```

```{r process-sobol-qoi2, cache=TRUE, eval=LOCAL}
xproc_2 <- tell(gensol_proc, Yproc_2)
```

```{r process-sobol-qoi-plot, cache=TRUE, eval=LOCAL}
print(xproc_2)
Sobol4R::autoplot(xproc_2, ncol = 1)
```

## Qoi Mean

```{r, cache=TRUE, eval=LOCAL}
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)
```

```{r, cache=TRUE, eval=LOCAL}
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

```{r, cache=TRUE, eval=LOCAL}
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)
```

```{r, cache=TRUE, eval=LOCAL}
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).