Sensitivity Analysis

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]#

Compute derivatives and elasticities of a model statistic with respect to its parameters.

One backward pass of automatic differentiation through the model gives exact derivatives with respect to every selected parameter at once.

Gradients are reported with respect to the constrained parameter values, the ones the user reads. Each Param bijector acts elementwise, so the chain rule divides by the bijector’s own derivative per parameter. A bijector derivative of zero gives a nan gradient.

Only statistics with a pathwise derivative are available (see DIFFERENTIABLE_STATISTICS). A sample quantile or CVaR selects among draws and is piecewise constant in the parameters, so its automatic derivative would be zero almost everywhere.

Parameters:

  • compiled – The model.

  • attribute_name (str) – The attribute whose statistic is differentiated.

  • statistic (str) – "mean" (the expectation) or "cdf" (the probability P(attribute <= level)).

  • level (float | None) – Threshold for "cdf".

  • wrt (Iterable[str] | None) – Scalar parameter names to differentiate with respect to. Default every Param of the model.

  • params (Mapping[str, Any] | None) – Unconstrained parameter values to evaluate at, with missing names at their initial values. Pass a fit’s raw_params to get sensitivities at the fitted point.

Returns:

A LocalSensitivity report.

Raises:

  • KeyError – If a name in wrt is not a model parameter.

  • ValueError – If wrt is empty, statistic is not differentiable, level is missing for "cdf", or a selected parameter is not scalar.

Return type:

LocalSensitivity

Example

>>> sens = local_sensitivity(compiled, "energy", wrt=["hub_height", "rotor_radius"])
>>> sens.ranked()
class prophys.engine.LocalSensitivity(attribute, statistic, value, names, gradients, elasticities)[source]#

Derivatives and elasticities of one model statistic with respect to the parameters.

Returned by local_sensitivity(). ranked() orders the parameters by influence.

Parameters:
attribute: str#

The attribute the statistic is of.

statistic: str#

"mean" or "cdf".

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 in different units can be ranked against each other. nan where the parameter or the statistic is zero.

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]#

Compute variance based (Sobol) sensitivity indices of an attribute to its uncertain inputs.

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

Uses the A, B, AB design: two independent samples A and B from each input’s declared distribution, plus one hybrid per selected factor (A with that factor’s column taken from B). The model is evaluated (k + 2) * n_samples times, in k + 2 batched calls, for k factors.

Inputs are drawn independently of each other, as the variance decomposition assumes. A JointRandomVariable (a copula, a multivariate Gaussian) is one factor whose components share one draw, since resampling one component of a correlated pair on its own would evaluate the model off its input distribution.

Parameters:

  • compiled – The model.

  • attribute_name (str) – The attribute whose variance is decomposed.

  • factors (Sequence[str] | None) – Names of upstream random variables to analyse. Default all of them.

  • n_samples (int) – Rows in each of the two base samples. The estimator’s error falls as 1/sqrt(n_samples).

  • seed (int) – PRNG seed for the base samples.

  • params (Mapping[str, Any] | None) – Unconstrained parameter values, such as a fit’s raw_params. Missing names take their initial values.

Returns:

A SobolSensitivity report.

Raises:

  • ValueError – If the attribute has no upstream random variables or factors is empty.

  • KeyError – If a factor is not an upstream random variable of the attribute.

Return type:

SobolSensitivity

Example

>>> sobol = sobol_sensitivity(compiled, "energy", n_samples=8192)
>>> print(sobol)
class prophys.engine.SobolSensitivity(attribute, names, first_order, total_order, n_samples)[source]#

First order and total order Sobol indices of an attribute over its uncertain inputs.

Returned by sobol_sensitivity(). ranked() orders the inputs by total order index.

Parameters:
attribute: str#

The attribute whose variance is decomposed.

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 for k inputs.

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 no single input explains, carried by interactions alone. Clamped at zero, since a negative value only means the Monte Carlo noise of the estimator exceeds this share.

Type:

1 - sum(S_i)

ranked()[source]#

Inputs ordered by total order index, largest first.

Return type:

list[tuple[str, float]]

Shapley attribution#

Sobol indices split the variance by first and total order. Shapley effects split it into one share per input that sums to the variance, stays fair when inputs are dependent, and extends to tail statistics such as the conditional value at risk. See prophys.explain.shapley_effects() and Explainability. The epistemic random variable of a learned surrogate appears in both analyses, so the share of output uncertainty owed to a surrogate’s ignorance is reported next to the physical inputs.