Sensitivity Analysis

Two questions, two functions. They answer different things and can disagree legitimately.

local_sensitivity()

sobol_sensitivity()

Question

How much does the answer move when this parameter moves, here?

How much of the answer’s variance does this uncertain input account for, across its whole range?

Method

Automatic differentiation of the model graph

Variance decomposition over Monte-Carlo draws

Cost

One backward pass, whatever the parameter count

(k + 2) * n_samples model evaluations

Valid where

Near the evaluation point

Over the entire input distribution

Reports

Derivatives and elasticities

First-order and total-order indices, hence interactions

A parameter with a large derivative that is nonetheless known precisely contributes little variance; an input whose effect is symmetric about its mean can have a first-order Sobol index of zero and a large local derivative.

Local sensitivity

import prophys as prp

report = prp.local_sensitivity(compiled, "power", params=result.raw_params)
print(report)
LocalSensitivity:
  mean[power] = 208.657
  d/d efficiency = +521.6  (elasticity +1.000)
  d/d noise = +0  (elasticity +0.000)

Two numbers per parameter:

Gradient\(\partial(\text{statistic}) / \partial(\text{parameter})\) in the parameter’s own units. Reported with respect to the constrained value, so the chain rule through Param.forward is already applied.

Elasticity\(\frac{\partial y}{\partial x}\cdot\frac{x}{y}\), the fractional change in the output per fractional change in the input. Being dimensionless, it ranks parameters measured in different units against each other; an elasticity of 1.0 means the output is proportional to that parameter. It is nan where either the parameter or the statistic is zero, and ranked() orders those last.

Only statistics with a pathwise derivative are available — the mean and the CDF. A sample quantile or CVaR is a selection among draws and is piecewise-constant in the parameters, so its automatic derivative would be zero almost everywhere; use Sobol indices for those.

Global (Sobol) sensitivity

Sobol indices decompose the variance of a model output over the joint distribution of its uncertain inputs — the attribute’s upstream RandomVariable leaves.

report = prp.sobol_sensitivity(compiled, "power", n_samples=20_000)
print(report)
SobolSensitivity:
  variance of power over 20000 samples
  wind: first 0.941, total 0.966 (interaction +0.025)
  direction: first 0.014, total 0.049 (interaction +0.035)
  unexplained by first-order effects: 0.045
  • First order \(S_i = \operatorname{Var}(E[Y \mid X_i]) / \operatorname{Var}(Y)\) — the share of output variance removed by learning \(X_i\) alone. Estimated by Saltelli (2010).

  • Total order \(S_{Ti} = E[\operatorname{Var}(Y \mid X_{\sim i})] / \operatorname{Var}(Y)\) — the share attributable to \(X_i\) through any route. Estimated by Jansen (1999).

  • Interaction \(S_{Ti} - S_i\) — the part of an input’s influence that acts only jointly with others. Near zero for a purely additive model.

  • Unexplained \(1 - \sum_i S_i\) — variance carried entirely by interactions.

The design is the standard three-matrix one: two independent samples A and B drawn from each input’s own declared distribution, plus one hybrid matrix per input. Inputs are sampled independently from their marginals, which is what the decomposition assumes; components of a JointRandomVariable (a copula, a multivariate Gaussian) count as one factor, since resampling one component of a correlated pair independently would evaluate the model off its own input distribution.

Scenario evaluation

The Sobol design is built on CompiledModel.conditional_expectation, which is directly useful on its own. Where expectation() marginalizes the uncertain inputs, this conditions on values you supply — a design scenario, a return-period event, a row of a sensitivity table:

import numpy as np

compiled.conditional_expectation("power", {"wind": np.array([4.0, 8.0, 12.0])})
# array([  32.,  256.,  864.], dtype=float32)

Every named leaf takes the values given for it, so passing arrays of length n evaluates n scenarios in one batched call. Every upstream random variable must be named; an unbound leaf raises rather than reverting to a sampled or default value.

prophys.engine.local_sensitivity(compiled, attribute_name, statistic='mean', *, level=None, wrt=None, params=None)[source]

Derivatives of a model statistic with respect to its parameters.

wrt selects parameters by name; the default is every Param in the model. params is the unconstrained parameter dict to evaluate at — pass a fit’s raw_params to get sensitivities at the fitted point rather than at the declared initial values.

Gradients are reported with respect to the constrained parameter values, the ones the user reads. Since each of Param’s bijectors acts elementwise, the chain rule there is a division by the bijector’s own derivative, applied per parameter.

Only statistics with a pathwise derivative are available (see DIFFERENTIABLE_STATISTICS): a sample quantile or CVaR is a selection among draws and is piecewise-constant in the parameters, so its automatic derivative would be zero almost everywhere rather than wrong-but-useful.

Parameters:
Return type:

LocalSensitivity

class prophys.engine.LocalSensitivity(attribute, statistic, value, names, gradients, elasticities)[source]

Derivatives of one model statistic with respect to the parameters.

Parameters:
value: float

The statistic itself, at the evaluation point.

names: list[str]

Parameters differentiated with respect to.

gradients: dict[str, float]

d(statistic) / d(parameter), in the parameter’s own units.

elasticities: dict[str, float]

the fractional change in the output per fractional change in the input. Dimensionless, so parameters measured in different units can be ranked against each other. nan where either the parameter or the statistic is zero, which makes the ratio undefined.

Type:

(d statistic / d parameter) * (parameter / statistic)

ranked()[source]

Parameters ordered by absolute elasticity, largest first. Parameters with an undefined elasticity are ordered last, by absolute gradient.

Return type:

list[tuple[str, float]]

prophys.engine.sobol_sensitivity(compiled, attribute_name, *, factors=None, n_samples=4096, seed=0, params=None)[source]

Variance-based sensitivity of an attribute to its uncertain inputs.

The uncertain inputs are the attribute’s upstream RandomVariable leaves; factors selects a subset by name. The output analysed is the attribute’s conditional expectation given those inputs — the part of the attribute that the inputs actually explain, with its own observation noise averaged out.

Uses the three-matrix design: two independent samples A and B drawn from each input’s own declared distribution, plus one hybrid matrix per input. The model is evaluated (k + 2) * n_samples times, batched into k + 2 calls.

Inputs are sampled independently from their marginals, which is what the variance decomposition assumes. For inputs coupled through a JointRandomVariable (a copula, a multivariate Gaussian), the joint draw is preserved — the components of one joint are one factor here, not several, since resampling one component of a correlated pair independently would evaluate the model off its own input distribution.

Parameters:
Return type:

SobolSensitivity

class prophys.engine.SobolSensitivity(attribute, names, first_order, total_order, n_samples)[source]

Variance-based sensitivity indices over a model’s uncertain inputs.

Parameters:
names: list[str]

The uncertain inputs, in index order.

first_order: dict[str, float]

the share of output variance removed by learning input i alone.

Type:

S_i

total_order: dict[str, float]

the share attributable to input i through any route, interactions included.

Type:

S_Ti

n_samples: int

Rows in each of the two base samples. The model was evaluated (k + 2) * n_samples times.

property interactions: dict[str, float]

the part of an input’s influence that acts only jointly with other inputs. Near zero for a purely additive model.

Type:

S_Ti - S_i

property unexplained: float

the share of variance not attributable to any single input, i.e. carried entirely by interactions. Negative values are clamped to zero; they mean the estimator’s Monte-Carlo noise exceeds the quantity being measured.

Type:

1 - sum(S_i)

ranked()[source]

Inputs ordered by total-order index, largest first.

Return type:

list[tuple[str, float]]