fips.Posterior.recovery#

Posterior.recovery(world=None, *, truth=None, W=None, level=0.95)[source]#

Return how well the estimate recovers the truth when the data come from world: exactly, without sampling.

The exact answer of the Monte Carlo loop over world.simulate(rng) and with_data(). The estimate is linear in the observations, \(\hat{x} = x_0 + K (z - H x_0 - c)\), with the gain \(K\) of this posterior (its problem’s \(H\), \(c\), \(S_0\) and \(S_z / \gamma\)). When the world makes the data, with a truth \(x \sim N(\mu, P)\) (its prior, or a given truth with \(P = 0\)) and observations \(z = H_w x + c_w + e\), \(e \sim N(0, R_w)\) (its forward operator, offset and model-data mismatch, not divided by \(\gamma\)), the error \(\hat{x} - x\) is Gaussian, with

\[b = (x_0 - \mu) + K (H_w \mu + c_w - H x_0 - c), \qquad V = (K H_w - I) P (K H_w - I)^T + K R_w K^T.\]

For an aggregation \(W\), \(W b\) and \(W V W^T\), from \(W K\) (k x m) and \(W K H_w - W\) (k x n), with \(P\) and \(R_w\) applied to k columns in their own forms (through the factors of a Kronecker covariance): nothing \(n \times n\) or \(m \times m\) beyond what the posterior holds. With the world the problem itself and \(\gamma = 1\), the bias is zero and \(V\) is \(\hat{S}\) (in Joseph form), so std equals posterior_std.

Parameters:

  • world (Problem | None) – The problem that makes the data, with exactly this posterior’s state and observation labels (in any order); its own observations are not used. By default the posterior’s problem. A world on other state labels (the whole state where the solve was restricted) is not taken: use world.simulate and with_data() in a loop.

  • truth (Vector | None) – A fixed truth, on the state’s labels in any order, instead of the world’s prior (another inventory, 1.5 times the prior). Only the observation error is then random.

  • W (Operator | None) – An aggregation, as cov() takes it: one row per total. By default, one row per state element.

  • level (float) – The probability of the posterior’s interval that coverage is computed for, \(\hat{x} \pm z_{level} \hat{\sigma}\).

Returns:

One row per aggregate (or state element), on the aggregation’s rows (or the state’s labels): prior (\(W x_0\)), truth (\(W \mu\)), estimate (the expected \(W \hat{x}\)), bias (estimate - truth), std (the standard deviation of the error, \(\sqrt{\operatorname{diag} W V W^T}\)), posterior_std (what the posterior claims, \(\sqrt{\operatorname{diag} W \hat{S} W^T}\)) and coverage (the probability that the level interval of the posterior holds the truth, for an error \(N(\text{bias}, \text{std}^2)\)). A variance below zero by roundoff gives a standard deviation of 0; one further below zero, which only a covariance that is not positive semidefinite gives, gives NaN (and a NaN coverage), with a warning logged. Roundoff is sqrt(eps) (1.5e-8) times the largest variance the row could have under covariances with the same variances, \((\sum_i |x_i| \sigma_i)^2\).

Return type:

DataFrame

Raises:

  • TypeError – If the method’s estimate is not linear in the observations (its solution does not set estimate_is_linear), or world is not a Problem, or truth not a Vector.

  • ValueError – If the world’s state or observations, the truth or the columns of W do not have exactly this posterior’s labels (the message names those on one side only), or level is not between 0 and 1.

See also

fips.Problem.simulate

One draw of a truth and its observations.

with_data

The posterior of one draw’s observations.

cov

The posterior covariance, of the state or of aggregates.

Examples

Paper 1’s experiment: a truth of 1.5 times the prior, observed without noise, and the fraction of the true change of each monthly total that the inversion recovers. The expected estimate is the one from noise-free observations, since the noise has mean zero; std is the spread the noise adds.

>>> import numpy as np
>>> import pandas as pd
>>> import fips
>>> state = pd.MultiIndex.from_product(
...     [pd.Index([1, 2], name="month"), pd.Index(["a", "b"], name="cell")]
... )
>>> prior = fips.Vector.from_series(pd.Series(1.0, index=state), block="flux")
>>> obs = fips.Vector.from_series(
...     pd.Series(0.0, index=pd.RangeIndex(6, name="obs_id")), block="ch4"
... )
>>> rng = np.random.default_rng(0)
>>> problem = fips.Problem(
...     prior=prior,
...     prior_cov=fips.Independent().with_std(0.5),
...     obs=obs,
...     obs_cov=fips.Independent().with_std(1.0),
...     forward=fips.Operator(
...         rng.uniform(0.0, 1.0, (6, 4)), obs.space, prior.space
...     ),
... )
>>> post = fips.solve(problem)
>>> monthly = fips.aggregate.sum(prior.space, by="month")
>>> truth = fips.Vector(prior.values * 1.5, prior.space)
>>> table = post.recovery(truth=truth, W=monthly)
>>> table = table.droplevel("block")
>>> table.round(3).reset_index()
   month  prior  truth  estimate   bias    std  posterior_std  coverage
0      1    2.0    3.0     2.667 -0.333  0.309          0.528     0.988
1      2    2.0    3.0     2.538 -0.462  0.285          0.586     0.992
>>> moved = table["estimate"] - table["prior"]
>>> needed = table["truth"] - table["prior"]
>>> (moved / needed).round(3).tolist()  # the recovered fraction
[0.667, 0.538]

In its own world (the default), an inversion at \(\gamma = 1\) has no bias, and its intervals hold the truth as often as they claim:

>>> honest = post.recovery(W=monthly).droplevel("block")
>>> honest[["bias", "std", "posterior_std", "coverage"]].round(3).reset_index()
   month  bias    std  posterior_std  coverage
0      1   0.0  0.528          0.528      0.95
1      2   0.0  0.586          0.586      0.95