Sobol Sensitivity Analysis

This example demonstrates a Sobol sensitivity analysis using the saltelli sampler and sobol function from SALib (https://github.com/SALib/SALib). In this case, the sensitivity of the sum-of-squared errors (sse) to model parameters is evaluated. It is also possible to provide the name of an individual observation instead of the sse as an argument in the sobol method (i.e., argument “obsname”).

import sys,os
import matk
import numpy as np
from matplotlib import pyplot as plt

# define a simple decaying sinusoidal function
def sine_decay(params, x, data):
    """ model decaying sine wave, subtract data"""
    amp = params['amp']
    shift = params['shift']
    omega = params['omega']
    decay = params['decay']

    model = amp * np.sin(x * omega + shift) * np.exp(-x*x*decay)

    obsnames = ['obs'+str(i) for i in range(1,len(data)+1)]
    return dict(zip(obsnames,model))


# create noisy data
x = np.linspace(0, 15, 301)
np.random.seed(1000)
data = (5. * np.sin(2 * x - 0.1) * np.exp(-x*x*0.025) +
        np.random.normal(size=len(x), scale=0.2) )

# Create MATK object
p = matk.matk(model=sine_decay, model_args=(x,data,))

# Create parameters
p.add_par('amp', value=10, min=0., max=20.)
p.add_par('decay', value=0.1, min=0, max=0.2)
p.add_par('shift', value=0.0, min=-np.pi/2., max=np.pi/2.)
p.add_par('omega', value=3.0, min=0, max=6)

# Create observation names and set observation values to noisy data
for i in range(len(data)):
    p.add_obs('obs'+str(i+1), value=data[i])

# Create Saltelli sample where the argument indicates the number of
# samples per parameter. The actual number of samples will be N * (2*D + 2),
# where N is the specified number of samples per parameter and D is the number
# of parameters.
# The default is True for "calc_second_order", explicitly specified here for clarity.
ss1 = p.saltelli(1000, calc_second_order=True)

# Execute the model on the Saltelli sample parameter combinations.
ss1.run(verbose=False)

# Perform the Sobol analysis
SS = ss1.sobol(calc_second_order=True)
Parameter S1 S1_conf ST ST_conf
amp 0.434041 0.066804 0.669720 0.173768
decay 0.261124 0.090263 0.461670 0.129540
shift 0.011050 0.014759 0.047576 0.013243
omega 0.011556 0.018900 0.073801 0.029723

Parameter_1 Parameter_2 S2 S2_conf
amp decay 0.164074 0.217512
amp shift -0.023936 0.091280
amp omega -0.009315 0.092907
decay shift -0.012035 0.141097
decay omega 0.009343 0.158747
shift omega 0.035670 0.028652
# Manually print some results
print SS['S1']
[ 0.43404054  0.26112364  0.01105031  0.01155562]
print SS['S2']
[[        nan  0.16407351 -0.02393598 -0.00931453]
 [        nan         nan -0.01203482  0.00934286]
 [        nan         nan         nan  0.03567022]
 [        nan         nan         nan         nan]]

The results indicate that the model is most sensitive to “amp” followed by “decay”. The model is relatively insensitive to “shift” and “omega”. Considering parameter interactions, the model is most sensitive to interactions between “amp” and “decay”.