gsaot

The package gsaot provides a set of tools to compute and plot Optimal Transport (OT) based sensitivity indices. The core functions of the package are:

• ot_indices(): compute OT indices for multivariate outputs using different solvers for OT (network simplex, Sinkhorn, and so on).

• ot_indices_wb(): compute OT indices for univariate or multivariate outputs using the Wasserstein-Bures semi-metric.

• ot_indices_1d(): compute OT indices for univariate outputs using OT solution in one dimension. The package provides also functions to plot the resulting indices and the inner statistics.

Installation

You can install the development version of gsaot from GitHub with:

# install.packages("devtools")
devtools::install_github("pietrocipolla/gsaot")

❗ ❗ Installation note

The sinkhorn and sinkhorn_log solvers in gsaot greatly benefit from optimization in compilation. To add this option (before package installation), edit your .R/Makevars file with the desired flags. Even though different compilers have different options, a common flag to enable a safe level of optimization is

CXXFLAGS+=-O2

More detailed information on how to customize the R packages compilation can be found in the R guide.

Example

A basic application of the functions implemented in the package:

library(gsaot)

N <- 1000

# Define the inputs distribution
mx <- c(1, 1, 1)
Sigmax <- matrix(data = c(1, 0.5, 0.5, 0.5, 1, 0.5, 0.5, 0.5, 1), nrow = 3)

# Sample the inputs
x1 <- rnorm(N)
x2 <- rnorm(N)
x3 <- rnorm(N)

x <- cbind(x1, x2, x3)
x <- mx + x %*% chol(Sigmax)

# Define the (linear)  model and the output
A <- matrix(data = c(4, -2, 1, 2, 5, -1), nrow = 2, byrow = TRUE)
y <- t(A %*% t(x))

x <- data.frame(x)

M <- 25

# Compute the sensitivity indices using Sinkhorn's solver and default parameters
sensitivity_indices <- ot_indices(x, y, M)
#> Using default values for solver sinkhorn
sensitivity_indices
#> Method: sinkhorn 
#> 
#> Indices:
#>        X1        X2        X3 
#> 0.6024491 0.5979907 0.2734895 
#> 
#> Upper bound: 93.27764

# Compute the sensitivity indices using the Network Simplex solver and default parameters
sensitivity_indices <- ot_indices(x, y, M, solver = "transport")
#> Using default values for solver transport
sensitivity_indices
#> Method: transport 
#> 
#> Indices:
#>        X1        X2        X3 
#> 0.5231955 0.5104445 0.1825898 
#> 
#> Upper bound: 93.27764

# Compute the Wasserstein-Bures indices
sensitivity_indices <- ot_indices_wb(x, y, M)
sensitivity_indices
#> Method: wasserstein-bures 
#> 
#> Indices:
#>        X1        X2        X3 
#> 0.5040783 0.4852300 0.1397265 
#> 
#> Advective component:
#>        X1        X2        X3 
#> 0.3050624 0.3105509 0.1221156 
#> 
#> Diffusive component:
#>         X1         X2         X3 
#> 0.19901591 0.17467902 0.01761093 
#> 
#> Upper bound: 93.27764

# Compute the sensitivity map using 1-dimensional solver
sensitivity_indices <- ot_indices_smap(x, y, M)
sensitivity_indices
#>             X1         X2        X3
#> [1,] 0.6090166 0.04683973 0.1790092
#> [2,] 0.3234579 0.70364026 0.1424398