1.4. Sensitivity Analysis

1.4.1. Introduction

This example shows how to apply pSeven Core to a sensitivity analysis problem.

Consider “The jumping man” test case from Sensitivity analysis in practice: a guide to assessing scientific models (A. Saltelli, 2004, John Wiley & Sons Inc).

The model is represented by a simple function:

\[h_{min} = H - \frac{2 M g}{k \sigma},\]

which is the solution of the linear oscillator equation. \(H\) is the distance of the platform to the surface [\(m\)], \(M\) is our mass [\(kg\)], \(\sigma\) is the number of strands in the cord, \(g\) is the the standard gravity [\(m/s^2\)], \(k\) is the elastic constant of one strand [\(N/m\)], which is here considered fixed at 1.5.

\begin{eqnarray} &40\ m \leq H \leq 60\ m \\ &67\ kg \leq M \leq 74\ kg \\ &20\ strands \leq \sigma \leq 40\ strands \end{eqnarray}

We are interested in selecting for which input factor we would gain a better level of accuracy in order to have the highest reduction of the uncertainty of the risk.

Start by importing the Generic Tool for Sensitivity and Dependency Analysis (GTSDA) module:

from da.p7core import gtsda

1.4.2. Blackbox-Based Analysis

from da.p7core import blackbox

First, use blackbox-based technique. Derive from blackbox.Blackbox and describe the problem:

import numpy as np
def analytic_function(x):
  g = 9.8
  k = 1.5
  f = x[:, 0] - (2 * x[:, 1] * g) / (k * x[:, 2])
  return f

class ExampleBlackbox(blackbox.Blackbox):
  def prepare_blackbox(self):

    # add new variables to problem
    self.add_variable((40, 60), 'H')
    self.add_variable((67, 74), 'M')
    self.add_variable((20, 40), 'sigma')

    # add new response to problem
    self.add_response()

  def evaluate(self, queryx):
    result = analytic_function(queryx)
    return result

Create a gtsda.Analyzer instance:

analyzer = gtsda.Analyzer()

Set options and logger (see Options Interface, Loggers):

from da.p7core import loggers

analyzer.options.set("GTSDA/Seed", 100)
analyzer.set_logger(loggers.StreamLogger())

Create the problem:

bbox = ExampleBlackbox()

Perform sensitivity analysis:

result = analyzer.rank(blackbox=bbox, budget=500)

Print a readable result summary:

print(result)

Print sensitivity indices:

print("Results:")
for i, s in enumerate(result.scores):
  print('score for blackbox response[%d]: %s' % (i, s))

1.4.3. Sample-Based Analysis

GTSDA can solve the same problem using a sample-based technique.

Generate a sample:

import numpy as np
X = np.hstack((np.random.uniform(low=40.0, high=60.0, size=(10000, 1)),
               np.random.uniform(low=67.0, high=74.0, size=(10000, 1)),
               np.random.uniform(low=20.0, high=40.0, size=(10000, 1))))
Y = analytic_function(X)

Perform analysis:

result = analyzer.rank(x=X, y=Y)

Print the result:

print(result)
print("Results:")
for i, s in enumerate(result.scores):
  print('score for output[%d]: %s' % (i, s))

1.4.4. Full Example Code

import numpy as np

from da.p7core import gtsda
from da.p7core import blackbox
from da.p7core import loggers

def analytic_function(x):
  g = 9.8
  k = 1.5
  f = x[:, 0] - (2 * x[:, 1] * g) / (k * x[:, 2])
  return f

class ExampleBlackbox(blackbox.Blackbox):
  def prepare_blackbox(self):

    # add new variables to problem
    self.add_variable((40, 60), 'H')
    self.add_variable((67, 74), 'M')
    self.add_variable((20, 40), 'sigma')

    # add new response to problem
    self.add_response()

  def evaluate(self, queryx):
    result = analytic_function(queryx)
    return result

analyzer = gtsda.Analyzer()

analyzer.options.set("GTSDA/Seed", 100)
analyzer.set_logger(loggers.StreamLogger())

bbox = ExampleBlackbox()

result = analyzer.rank(blackbox=bbox, budget=500)

print(result)

print("Results:")
for i, s in enumerate(result.scores):
  print('score for blackbox response[%d]: %s' % (i, s))

X = np.hstack((np.random.uniform(low=40.0, high=60.0, size=(10000, 1)),
               np.random.uniform(low=67.0, high=74.0, size=(10000, 1)),
               np.random.uniform(low=20.0, high=40.0, size=(10000, 1))))
Y = analytic_function(X)

result = analyzer.rank(x=X, y=Y)

print(result)

print("Results:")
for i, s in enumerate(result.scores):
  print('score for output[%d]: %s' % (i, s))