inst/doc/gaussian-model-sensitivity.R

In gsaot: Compute Global Sensitivity Analysis Indices Using Optimal Transport

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## -----------------------------------------------------------------------------
# Define the input distribution parameters
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)

# Define the number of samples
N <- 1000

# Set the random number generator seed for reproducibility
set.seed(777)

# Sample from standard normals
x1 <- rnorm(N)
x2 <- rnorm(N)
x3 <- rnorm(N)

# Transform the standard normals into the required distribution
x <- cbind(x1, x2, x3)
x <- mx + x %*% chol(Sigmax)

# Define the model (matrix with coefficients)
A <- matrix(data = c(4, -2, 1, 2, 5, -1), nrow = 2, byrow = TRUE)

# Generate the output
y <- t(A %*% t(x))
colnames(y) <- c("y1", "y2")

## ----fig.width=4, fig.height=4------------------------------------------------
library(ggplot2)
library(patchwork)

# Set the theme
theme_set(theme_minimal())

# Prepare the marginal and the 2D density plots
p1 <- ggplot(as.data.frame(y), aes(x = y1, y = y2)) +
  geom_point(color = "black") +
  geom_density_2d()
p2 <- ggplot(as.data.frame(y), aes(x = y1)) +
  geom_density()
p3 <- ggplot(as.data.frame(y), aes(x = y2)) +
  geom_density() +
  coord_flip()

p2 + plot_spacer() + p1 + p3 +
  plot_layout(ncol = 2, nrow = 2, widths = c(4, 1), heights = c(1, 4), axes = "collect")

## -----------------------------------------------------------------------------
library(gsaot)

M <- 25
indices_wb <- ot_indices_wb(x, y, M)
indices_wb

## ----fig.width=6, fig.height=4------------------------------------------------
plot(indices_wb)

## ----fig.width=6, fig.height=4------------------------------------------------
# Enable bootstrap
boot <- TRUE
# Set the number of replicas
R <- 100
# Define the confidence level
conf <- 0.99
# Define the type of confidence interval
type <- "norm"

# Compute the indices
indices_wb <- ot_indices_wb(x, y, M, 
                            boot = boot, R = R, conf = conf, type = type)
plot(indices_wb)