Uncertainty

Propagating the uncertainty of inputs to results: the standard deviation of a total of uncertain pixel values, of a product of two uncertain layers, of a result with respect to uncertain parameters (the delta method, checked by Monte Carlo), and confidence intervals. The ls.uncertainty module; the forms are those plangea uses in its restoration plans, with their errors corrected (each function's notes say which).

import landschaft as ls

ls.uncertainty.sum_sd(carbon_sd, weights=area * restored)          # a total's deviation
ls.uncertainty.product_sd(x, sd_x, v, sd_v)                         # each product's
d = ls.uncertainty.delta_method(f, params, sd)                      # d.value, d.sd
low, high = ls.uncertainty.interval(d.value, d.sd, level=0.95)

Totals and products

sum_sd

sum_sd(
    sd, weights=None, *, method: str = "quadratic", correlation: float | None = None
) -> float

Standard deviation of a total of uncertain values, sum of w_i v_i.

Each value v_i (a pixel's carbon, say) has a standard deviation sigma_i, and the total weighs it by w_i (the pixel's area times the share restored). The total's standard deviation is

• "quadratic": errors independent, sqrt(sum (sigma_i w_i)^2);
• "linear": errors that add up, sum |sigma_i w_i| (an upper bound, reached when they all go the same way);
• with correlation rho (one value for every pair): sqrt((1 - rho) sum s_i^2 + rho (sum s_i)^2), s_i = sigma_i w_i, which is the quadratic form for 0 and |sum s_i| for 1.

Parameters:

Name	Type	Description	Default
sd		The standard deviations (an array or raster; NaN left out).	required
weights		The weights (same shape, or a number; default 1).	None
method	str	"quadratic" or "linear" (ignored with correlation).	'quadratic'
correlation	float | None	The correlation of every pair of errors, -1 to 1.	None

Returns:

Type	Description
float

Notes

plangea's prop_quad and prop_lin are the two methods.

Examples:

>>> import landschaft as ls
>>> ls.uncertainty.sum_sd([3.0, 4.0])
5.0
>>> ls.uncertainty.sum_sd([3.0, 4.0], method="linear")
7.0
>>> ls.uncertainty.sum_sd([3.0, 4.0], correlation=1.0)
7.0

product_sd

product_sd(x, sd_x, v, sd_v, *, total: bool = False, method: str = 'quadratic')

Standard deviation of the product of two uncertain values, x v.

With x and v independent, the exact variance of the product (Goodman 1960) is x^2 sigma_v^2 + v^2 sigma_x^2 + sigma_x^2 sigma_v^2.

Parameters:

Name	Type	Description	Default
x		The values and their standard deviations (arrays or rasters of one shape, or numbers).	required
sd_x		The values and their standard deviations (arrays or rasters of one shape, or numbers).	required
v		The values and their standard deviations (arrays or rasters of one shape, or numbers).	required
sd_v		The values and their standard deviations (arrays or rasters of one shape, or numbers).	required
total	bool	Return the standard deviation of the sum of the products instead of each product's: quadratic, sqrt(sum of the variances); linear, the sum of |x| sigma_v + |v| sigma_x + sigma_x sigma_v.	False
method	str	For total: "quadratic" (independent products) or "linear".	'quadratic'

Returns:

Type	Description
ndarray or DataArray or float

Notes

plangea's aux_calc_var_sd adds the cross term sigma_x sigma_v without squaring it in the quadratic form, and by default takes 10 % of the deviation of x as the deviation of v; here the deviations are always given and the variance is Goodman's.

Examples:

>>> import landschaft as ls
>>> round(float(ls.uncertainty.product_sd(2.0, 0.1, 3.0, 0.2)), 4)
0.5004

Parameters

delta_method

delta_method(
    f: Callable,
    params,
    sd,
    *,
    covariance=None,
    method: str = "quadratic",
    step=None,
    scheme: str = "central",
) -> DeltaResult

Uncertainty of a result from the uncertainty of its parameters (delta method).

The result f(params) (a number, an array, a raster) changes with each parameter p_k by df/dp_k, found here by finite differences; its standard deviation is

• "quadratic" (independent parameters): sqrt(sum_k (df/dp_k sigma_k)^2), or sqrt(J C J^T) with a covariance matrix C;
• "linear": sum_k |df/dp_k sigma_k|.

It is the first-order propagation of the GUM; monte_carlo() checks it when f is far from linear over the parameters' spread.

Parameters:

Name	Type	Description	Default
f	Callable	The function: f(params) with the parameters as a 1-D array.	required
params		The parameters' values.	required
sd		Their standard deviations (ignored with covariance).	required
covariance		The parameters' covariance matrix.	None
method	str	"quadratic" or "linear".	'quadratic'
step		The finite-difference step of every parameter (default: the cube root (central) or square root (forward) of the machine epsilon, times max(|p|, 1)).	None
scheme	str	"central" (f at p + h and p - h) or "forward" (p + h and p), plangea's.	'central'

Returns:

Type	Description
DeltaResult	The value, the standard deviations and the Jacobian.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> r = ls.uncertainty.delta_method(lambda p: p[0] * p[1], [2.0, 3.0], [0.1, 0.2])
>>> round(float(r.sd), 4)
0.5

DeltaResult dataclass

DeltaResult(value: object, sd: object, jacobian: ndarray)

Result of delta_method().

Attributes:

Name	Type	Description
value	object	The function at the parameters.
sd	object	The standard deviation of every output.
jacobian	ndarray	The derivatives of the outputs (rows, flattened) with respect to the parameters (columns).

monte_carlo

monte_carlo(
    f: Callable,
    params,
    sd,
    *,
    n: int = 1000,
    covariance=None,
    seed: int = 0,
    quantiles=(0.025, 0.975),
    keep: bool = False,
    on=None,
    chunk: int | None = None,
) -> MonteCarloResult

Uncertainty of a result by drawing its parameters (Monte Carlo).

The parameters are drawn n times from normal distributions (their sd, or the covariance matrix) with a fixed seed, and the result is computed for every draw: its mean, standard deviation and quantiles describe its uncertainty without assuming f linear.

Parameters:

Name	Type	Description	Default
f	Callable	The function: f(params) with the parameters as a 1-D array.	required
params		The parameters' means and standard deviations.	required
sd		The parameters' means and standard deviations.	required
n	int	The number of draws.	1000
covariance		The parameters' covariance matrix (instead of sd).	None
seed	int	The seed of the draws (NumPy's PCG64).	0
quantiles		The quantiles to report.	(0.025, 0.975)
keep	bool	Keep every draw's outputs in samples.	False
on		ls.helpers(): compute the draws here and on the helpers of this computer's group (parallel), with the same results; f is then a tool of landschaft with its arguments, call(), the drawn parameters marked by param().	None
chunk	int | None	Draws per task with on (default: about a hundredth of them).	None

Returns:

Type	Description
MonteCarloResult

Examples:

>>> import landschaft as ls
>>> r = ls.uncertainty.monte_carlo(lambda p: p[0] * p[1], [2.0, 3.0], [0.1, 0.2], n=20000, seed=1)
>>> round(float(r.sd), 2)
0.5

MonteCarloResult dataclass

MonteCarloResult(
    mean: object, sd: object, quantiles: ndarray, samples: ndarray | None = None
)

Result of monte_carlo().

Attributes:

Name	Type	Description
mean, sd		The mean and the standard deviation of every output over the draws.
quantiles	ndarray	The quantiles asked for, one row per quantile.
samples	ndarray | None	The outputs of every draw (draws x outputs), if kept.

Intervals

interval

interval(value, sd, level: float = 0.95)

The normal confidence interval of a value: value -/+ z sd, z the normal quantile of the level.

plangea fixes z at 1.96 (the 95 % level).

Parameters:

Name	Type	Description	Default
value		The value and its standard deviation (numbers, arrays or rasters).	required
sd		The value and its standard deviation (numbers, arrays or rasters).	required
level	float	The confidence level, 0 to 1.	0.95

Returns:

Type	Description
tuple	The lower and the upper bound.

Examples:

>>> import landschaft as ls
>>> low, high = ls.uncertainty.interval(10.0, 1.0)
>>> round(low, 3), round(high, 3)
(8.04, 11.96)