fips.Posterior.chi2_tests#

Posterior.chi2_tests()[source]#

Test whether the innovation, the residual and the increment are as large as the covariances say: Rodgers’ four chi-squared tests.

Each test takes a vector that is \(N(0, C)\) when the prior error and model-data mismatch covariances are right, and its chi-squared \(v^T C^{+} v\) (\(^{+}\) the pseudo-inverse), which then follows a chi-squared distribution whose degrees of freedom are the rank of \(C\): it averages them. With \(G = H S_0 H^T + S_z / \gamma\) (innovation_cov) and \(r\) the rank of \(H S_0\), the number of directions of the state the observations constrain at all:

innovation: the observations against what the prior models, \(d = z - H x_0 - c\), with \(C = G\): \(m\) degrees of freedom. Its chi-squared is chi2: this row is chi2 and reduced_chi2 with their p-value.
residual: the observations against what the estimate models, \(z - H \hat{x} - c\), with \(C = (S_z/\gamma) G^{-1} (S_z/\gamma)\): \(m\) degrees of freedom.
modeled_increment: how far the estimate moves the modeled observations from the prior’s, \(H (\hat{x} - x_0)\), with \(C = H S_0 H^T G^{-1} H S_0 H^T\): \(r\) degrees of freedom.
increment: the estimate against the prior, \(\hat{x} - x_0\), with \(C = S_0 - \hat{S} = S_0 H^T G^{-1} H S_0\): \(r\) degrees of freedom.

Rodgers (2000, chapter 12) names them the measurement against the prior, the retrieval against the measurement, and the retrieval against the prior in the measurement and in the state space, as pyOptimalEstimation computes them. The first two test both covariances at once; the last two whether the observations move the estimate as far from the prior as the prior error covariance allows.

In a linear problem the four are two pairs of the same numbers, so a table with duplicate rows is expected: the residual is \((S_z/\gamma) G^{-1} d\), so its chi-squared is \(d^T G^{-1} d\), the innovation’s; and the two increments are the same components of the innovation, scaled, so their chi-squared is one sum (below). When the observations constrain as many directions as there are observations (\(r = m\), usual when \(n \ge m\)), all four are the same number. They differ when the estimate is not linear in the observations (a nonlinear method, whose \(H\) is the forward model linearized at the estimate): the table has four rows so that every method gives it the same shape.

The degrees of freedom are the ranks of the four \(C\). \(G\) and \(S_z\) are positive definite, so the first two have rank \(m\). For a positive-definite \(G^{-1}\), \(B^T G^{-1} B\) has the rank of \(B\): \(B = H S_0 H^T\) for the third and \(B = H S_0\) for the fourth, both of rank \(r\) (with a positive-definite \(S_0\), \(r\) is the rank of \(H\)). So \(r \le \min(m, n)\): with fewer state elements than observations the increments have at most \(n\), and observations that see nothing (\(H S_0 = 0\)) give them none.

The analytic method computes the increments from the eigendecomposition of its factors: with \(s_k^2\) its eigenvalues times \(\gamma\) and \(e'_k\) the innovation whitened by \(S_z / \gamma\) and rotated onto its eigenvectors (variance \(1 + s_k^2\)), the increment has the components \(s_k e'_k / (1 + s_k^2)\) in the state whitened by \(S_0\) and the modeled increment \(s_k^2 e'_k / (1 + s_k^2)\) in the observations whitened by \(S_z / \gamma\); both chi-squared are \(\sum_k e_k'^2 / (1 + s_k^2)\) over the \(r\) eigenvalues above zero. After fips.solve() this eigendecomposes once, as the first with_gamma() does, and keeps it; a Monte Carlo over with_data() should start from post.with_gamma(gamma), whose eigendecomposition every trial shares. An eigenvalue counts when it is above roundoff relative to the largest, so the units of the problem do not change the degrees of freedom.

Returns:: One row per test (innovation, residual, modeled_increment, increment), indexed by test: chi2, dof (the degrees of freedom), reduced_chi2 (chi2 / dof, near 1 when the covariances are right) and p_value, the probability of a chi-squared at least this large when they are right. A small p-value says the vector is larger than the covariances allow (they are too small); one near 1 that it is smaller (they are too large). With no degrees of freedom (the observations constrain nothing) the last two are NaN.
Return type:: DataFrame

References

Rodgers, C. D. (2000). Inverse Methods for Atmospheric Sounding: Theory and Practice, chapter 12. World Scientific.

See also

chi2: The chi-squared of the innovation.
desroziers: The observation error variances the residuals diagnose.

Examples

Three observations of two state elements: the innovation and the residual have three degrees of freedom, the two increments two. The problem is linear, so the rows come in pairs of the same numbers.

>>> import numpy as np
>>> import pandas as pd
>>> import fips
>>> cells = pd.Index(["a", "b"], name="cell")
>>> obs_id = pd.Index(["o1", "o2", "o3"], name="obs_id")
>>> prior = fips.Vector.from_series(pd.Series(0.0, index=cells), block="flux")
>>> obs = fips.Vector.from_series(
...     pd.Series([1.0, 2.0, 4.0], index=obs_id), block="co2"
... )
>>> problem = fips.Problem(
...     prior=prior,
...     prior_cov=fips.Independent().with_std(1.0),
...     obs=obs,
...     obs_cov=fips.Independent().with_std(np.sqrt(0.1)),
...     forward=fips.Operator(
...         np.array([[1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]),
...         obs.space,
...         prior.space,
...     ),
... )
>>> fips.solve(problem).chi2_tests().round(3).reset_index()
                test    chi2  dof  reduced_chi2  p_value
0         innovation  10.293    3         3.431    0.016
1           residual  10.293    3         3.431    0.016
2  modeled_increment   6.960    2         3.480    0.031
3          increment   6.960    2         3.480    0.031