Global Sensitiviy Analysis (GSA) test functions

How do I run this code?

The full code is found in the examples/Emulator/ directory of the github repository

In this example, we assess directly the performance of our machine learning emulators. The task is to learn a function for use in a global sensitivity analysis. In particular, we have two cases:

  • The Ishigami function

\[f(x; a, b) = (1 + bx_3^4)\sin(x_1) + a \sin(x_2), \forall x\in [-\pi,\pi]^3\]

with $a=7, b=0.1$.

  • The Sobol G-function

\[f(x; a) = \prod_{i=1}^d \frac{|4x_i - 2|+a_i}{1+x_i},\]

with the user-chosen input dimension $d$, and coefficients $a_i = \frac{i-1}{2} \geq 0$ for $i=1,\dots,d$, where small $i$ (and thus small $a_i$) imply larger first-order effects, and interactions are primarily present between these variables.

In this example, global sensitivity analysis refers to calculation of two Sobol indices. The first index collects proportions $V_i$ (a.k.a firstorder) of the variance of $f$ attributable to the input $x_i$, and the second index collects proportions $TV_i$ (a.k.a totalorder) of the residual variance having removed contributions attributable to inputs $x_j$ $\forall j\neq i$. The Ishigami function has an analytic formula for these Sobol indices, it is also known that one can obtain numerical approximation through quasi-Monto-Carlo methods by evaluating the Ishigami function on a special quasi-random Sobol sequence.

Ishigami

To emulate the Ishigami function, the data consists of 300 pairs $\{x,f(x)+\eta\}$ where $\eta \sim N(0,\Sigma)$ is additive noise, and the x are sampled from the Sobol sequence. The emulators are validated by evaluating the posterior mean function on the full 16000 points of the Sobol sequence and the Sobol indices are estimated. We rerun the experiment for many repeats of the random feature hyperparameter optimization and present the statistics of these indices, as well as plotting a realization of the emulator.

We use the package GlobalSensitivityAnalysis.jl for many of the GSA tools.

Walkthrough of the code

We first load some standard packages

using Distributions
using DataStructures # for `OrderedDict`
using Random
using LinearAlgebra
using CairoMakie, ColorSchemes

and the packages for providing the Ishigami function, Sobol sequence, and evaluation of the indices

using GlobalSensitivityAnalysis # for `SobolData`
const GSA = GlobalSensitivityAnalysis

then the CES packages for the emulators

using CalibrateEmulateSample.Emulators # for `SKLJL`, `GaussianProcess`, `SeparableKernel`, `LowRankFactor`, `OneDimFactor`, `ScalarRandomFeatureInterface`, `Emulator`
using CalibrateEmulateSample.DataContainers # for `PairedDataContainer`
using CalibrateEmulateSample.EnsembleKalmanProcesses # for `DataMisfitController`

We set up the sampling procedure, evaluate the true ishigami function on these points, and calculate the sobol indices

n_data_gen = 2000 

data = SobolData(
    params = OrderedDict(:x1 => Uniform(-π, π), :x2 => Uniform(-π, π), :x3 => Uniform(-π, π)),
    N = n_data_gen,
)

# To perform global analysis,
# one must generate samples using Sobol sequence (i.e. creates more than N points)
samples = GSA.sample(data)
n_data = size(samples, 1) # [n_samples x 3]
# run model (example)
y = GSA.ishigami(samples)
# perform Sobol Analysis
result = analyze(data, y)

Next we create the noisy training data from the sequence samples

n_train_pts = 300
ind = shuffle!(rng, Vector(1:n_data))[1:n_train_pts]
# now subsample the samples data
n_tp = length(ind)
input = zeros(3, n_tp)
output = zeros(1, n_tp)
Γ = 1e-2
noise = rand(rng, Normal(0, Γ), n_tp)
for i in 1:n_tp
    input[:, i] = samples[ind[i], :]
    output[i] = y[ind[i]] + noise[i]
end
iopairs = PairedDataContainer(input, output)

We have a few cases for the user to investigate

cases = [
    "Prior", # Scalar random feature emulator with no hyperparameter learning
    "GP", # Trained Sci-kit learn Gaussian process emulator
    "RF-scalar", # Trained scalar random feature emulator
]

Each case sets up a different machine learning configuration in the Emulator object.

For the random feature case, RF-scalar, we use a rank-3 kernel in the input space, and 500 features for prediction, though for efficiency we use only 200 when learning the hyperparameters. The relevant code snippets are

nugget = Float64(1e-12)
overrides = Dict(
    "scheduler" => DataMisfitController(terminate_at = 1e4),
    "n_features_opt" => 200,
)

kernel_structure = SeparableKernel(LowRankFactor(3, nugget), OneDimFactor())
n_features = 500
mlt = ScalarRandomFeatureInterface(
    n_features,
    3,
    rng = rng,
    kernel_structure = kernel_structure,
    optimizer_options = overrides,
)

For the gaussian process case GP we use the sci-kit learn package, a default squared exponential kernel with lengthscale learnt in each input dimensions. We do not learn an additional white kernel for the noise.

gppackage = Emulators.SKLJL()
pred_type = Emulators.YType()
mlt = GaussianProcess(
    gppackage;
    prediction_type = pred_type,
    noise_learn = false,
)

We finish by building the emulator object

emulator = Emulator(mlt, iopairs; obs_noise_cov = Γ * I, decorrelate = decorrelate)
optimize_hyperparameters!(emulator) # (only needed for some Emulator packages)

Results and plots

We validate the emulator by evaluating it on the entire Sobol sequence, and calculating the Sobol indices (presenting mean and std if using multiple repeats.

# predict on all Sobol points with emulator (example)    
y_pred, y_var = predict(emulator, samples', transform_to_real = true)

# obtain emulated Sobol indices
result_pred = analyze(data, y_pred')

Gaussian Process Emulator (sci-kit learn GP)

Here is the plot for one emulation

and the outputted table of Sobol indices

True Sobol Indices
******************
    firstorder: [0.31390519114781146, 0.44241114479004084, 0.0]
    totalorder: [0.5575888552099592, 0.44241114479004084, 0.24368366406214775]
 
Sampled truth Sobol Indices (# points 16000)
***************************
    firstorder: [0.31261591941512257, 0.441887746224135, -0.005810964687365922]
    totalorder: [0.5623611180844434, 0.44201284296404386, 0.24465876318633062]
 
Sampled Emulated Sobol Indices (# obs 300, noise var 0.01)
***************************************************************
    firstorder: [0.3094638183079643, 0.4518400892052567, -0.007351344957260407]
    totalorder: [0.5502469909342245, 0.4587734278791574, 0.23542404141319245]

Random feature emulator (Separable Low-rank kernel RF-scalar)

Here is the plot for one emulation

Table for 20 repeats of the algorithm

True Sobol Indices
******************
    firstorder: [0.31390519114781146, 0.44241114479004084, 0.0]
    totalorder: [0.5575888552099592, 0.44241114479004084, 0.24368366406214775]
 
Sampled truth Sobol Indices (# points 16000)
***************************
    firstorder: [0.31261591941512257, 0.441887746224135, -0.005810964687365922]
    totalorder: [0.5623611180844434, 0.44201284296404386, 0.24465876318633062]
 
Sampled Emulated Sobol Indices (# obs 300, noise var 0.01)
***************************************************************
(mean) firstorder: [0.33605548545490044, 0.41116050093679196, -0.0012213648484969539]
(std)  firstorder: [0.05909336956162543, 0.11484966121124164, 0.012908533302492602]
(mean) totalorder: [0.5670345355855254, 0.4716028261179354, 0.24108222433317]
(std)  totalorder: [0.10619345801872732, 0.1041023777237331, 0.07200225781785778]

Sobol G-function results

To emulate the Sobol function, a similar code script is used to set up the Ishigami emulation. The primary change is that the input dimension is now a user parameter n_dimension that can be adjusted, and some reasonable defaults are set within the script. As an output, plots are produced of the Sobol function values, and slices through the function.

For example, we repeat the scalar random feature emulation task 30 times over different training realizations.

For three input dimensions, one obtains the following plot of the analytic indices (X-true), qMC-approximated with true function (X-approx), and the 95% confidence interval of the qMC-approximated with emulator (X-emulate).

One also obtains the slices through the emulated G-function, with red being the training points and blue being the prediction

For ten input dimensions one obtains similar plots

Here we plot only slices through the three most sensitive dimensions