The HSIC sensitivity indices: the Ishigami model

import openturns as ot
import openturns.viewer as otv
from openturns.usecases import ishigami_function

This example is a brief overview of the HSIC sensitivity indices classes and how to call them. HSIC estimators rely on a reproducing kernel of a Hilbert space. We can use them to compute sensitivity indices. We present the methods on the Ishigami function.

Definition of the model

We load the model from the usecases module.

im = ishigami_function.IshigamiModel()

We generate an input sample of size 100 (and dimension 3).

size = 100
X = im.distributionX.getSample(size)

We compute the output by applying the Ishigami model to the input sample.

Y = im.model(X)

Setting the covariance models

The HSIC algorithms use reproducing kernels defined on Hilbert spaces to estimate independence. For each input variable we choose a covariance kernel. Here we choose a SquaredExponential kernel for all input variables.

They are all stored in a list of covariance kernels where is the number of input variables. The remaining one is for the output variable.

covarianceModelCollection = []

for i in range(3):
    Xi = X.getMarginal(i)
    inputCovariance = ot.SquaredExponential(1)
    inputCovariance.setScale(Xi.computeStandardDeviation())
    covarianceModelCollection.append(inputCovariance)

Likewise we define a covariance kernel associated to the output variable.

outputCovariance = ot.SquaredExponential(1)
outputCovariance.setScale(Y.computeStandardDeviation())
covarianceModelCollection.append(outputCovariance)

The Global HSIC estimator

In this paragraph, we perform the analysis on the raw data: that is the global HSIC estimator.

Choosing an estimator

After having defined the covariance kernels one has to select an appropriate estimator for the computations.

Two estimators are proposed:

• an unbiased estimator through the HSICUStat class

• a biased, but asymptotically unbiased, estimator through the HSICVStat class

Beware that the conditional analysis used later cannot be performed with the unbiased estimator.

estimatorType = ot.HSICUStat()

We now build the HSIC estimator:

globHSIC = ot.HSICEstimatorGlobalSensitivity(
    covarianceModelCollection, X, Y, estimatorType
)

We get the R2-HSIC indices:

R2HSICIndices = globHSIC.getR2HSICIndices()
print("\n Global HSIC analysis")
print("R2-HSIC Indices: ", R2HSICIndices)

 Global HSIC analysis
R2-HSIC Indices:  [0.249305,-0.00429972,0.0437032]

and the HSIC indices:

HSICIndices = globHSIC.getHSICIndices()
print("HSIC Indices: ", HSICIndices)

HSIC Indices:  [0.0204961,-0.000366135,0.00366669]

The p-value by permutation.

pvperm = globHSIC.getPValuesPermutation()
print("p-value (permutation): ", pvperm)

p-value (permutation):  [0,0.50495,0.00990099]

We have an asymptotic estimate of the value for this estimator.

pvas = globHSIC.getPValuesAsymptotic()
print("p-value (asymptotic): ", pvas)

p-value (asymptotic):  [4.62161e-12,0.553716,0.0159241]

We vizualise the results.

graph1 = globHSIC.drawHSICIndices()
view1 = otv.View(graph1)

graph2 = globHSIC.drawPValuesAsymptotic()
view2 = otv.View(graph2)

graph3 = globHSIC.drawR2HSICIndices()
view3 = otv.View(graph3)

graph4 = globHSIC.drawPValuesPermutation()
view4 = otv.View(graph4)

• 
• 
• 
• 

The Target HSIC estimator

We now perform the target analysis which consists in using a filter function over the output.

Defining a filter function

We define a filter function on the output variable for the target analysis. In this example we use the function where is the distance to a well-chosen interval.

We first define a critical domain: in our case that is the interval.

criticalDomain = ot.Interval(5, float("inf"))

We have access to the distance to this domain thanks to the DistanceToDomainFunction class.

dist2criticalDomain = ot.DistanceToDomainFunction(criticalDomain)

We define the parameters in our function from the output sample

s = 0.1 * Y.computeStandardDeviation()[0]

We now define our filter function by composition of the parametrized function and the distance function.

f = ot.SymbolicFunction(["x", "s"], ["exp(-x/s)"])
phi = ot.ParametricFunction(f, [1], [s])
filterFunction = ot.ComposedFunction(phi, dist2criticalDomain)

We choose an unbiased estimator

estimatorType = ot.HSICUStat()

and build the HSIC estimator

targetHSIC = ot.HSICEstimatorTargetSensitivity(
    covarianceModelCollection, X, Y, estimatorType, filterFunction
)

We get the R2-HSIC indices:

R2HSICIndices = targetHSIC.getR2HSICIndices()
print("\n Target HSIC analysis")
print("R2-HSIC Indices: ", R2HSICIndices)

 Target HSIC analysis
R2-HSIC Indices:  [0.260374,0.00207302,-0.00658141]

and the HSIC indices:

HSICIndices = targetHSIC.getHSICIndices()
print("HSIC Indices: ", HSICIndices)

HSIC Indices:  [0.00108746,8.96771e-06,-2.80515e-05]

The p-value by permutation.

pvperm = targetHSIC.getPValuesPermutation()
print("p-value (permutation): ", pvperm)

p-value (permutation):  [0,0.237624,0.693069]

We have an asymptotic estimate of the value for this estimator.

pvas = targetHSIC.getPValuesAsymptotic()
print("p-value (asymptotic): ", pvas)

p-value (asymptotic):  [1.42697e-09,0.316344,0.59066]

We vizualise the results.

graph5 = targetHSIC.drawHSICIndices()
view5 = otv.View(graph5)

graph6 = targetHSIC.drawPValuesAsymptotic()
view6 = otv.View(graph6)

graph7 = targetHSIC.drawR2HSICIndices()
view7 = otv.View(graph7)

graph8 = targetHSIC.drawPValuesPermutation()
view8 = otv.View(graph8)

• 
• 
• 
• 

The Conditional HSIC estimator

In this last section we preprocess the input variables: that is the conditional analysis. To do so, one has to work with a weight function. Here the weight function is the filter function we used previously.

weightFunction = filterFunction

We have to select a biased -but asymptotically unbiased- estimator

estimatorType = ot.HSICVStat()

We build the conditional HSIC estimator

condHSIC = ot.HSICEstimatorConditionalSensitivity(
    covarianceModelCollection, X, Y, weightFunction
)

We get the R2-HSIC indices:

R2HSICIndices = condHSIC.getR2HSICIndices()
print("\n Conditional HSIC analysis")
print("R2-HSIC Indices: ", R2HSICIndices)

 Conditional HSIC analysis
R2-HSIC Indices:  [0.155438,0.017438,0.188395]

and the HSIC indices:

HSICIndices = condHSIC.getHSICIndices()
print("HSIC Indices: ", HSICIndices)

HSIC Indices:  [0.00492342,0.000779929,0.00858311]

For the conditional estimator we only have access to the p-value by permutation:

pvperm = condHSIC.getPValuesPermutation()
print("p-value (permutation): ", pvperm)
print("")

p-value (permutation):  [0.029703,0.712871,0]

We can vizualise the results thanks to the various drawing methods.

graph9 = condHSIC.drawHSICIndices()
view9 = otv.View(graph9)

graph10 = condHSIC.drawR2HSICIndices()
view10 = otv.View(graph10)

graph11 = condHSIC.drawPValuesPermutation()
view11 = otv.View(graph11)

• 
• 
• 

otv.View.ShowAll()