CoordTransform Cadastral & High-Precision Python

Chi-Square Testing of Transformation Residuals

A residual set can look acceptable and still be telling you the stochastic model is wrong, and the chi-square test on the variance factor is what turns that suspicion into a decision. This guide, part of error distribution modeling in Python, applies the global model test to transformation residuals: what statistic to compute, what its distribution should be, and what each side of a rejection actually means for a cadastral deliverable.

The Statistic and Its Distribution

After a least-squares fit with nn observations and uu parameters, the a posteriori variance factor is

The two tails of the global model test The scaled weighted sum of squares is compared against two critical values. Above the upper value the residuals are larger than the weights predicted, meaning either optimistic weights or unmodelled systematic error, and the reported uncertainties are too small — the dangerous direction. Below the lower value the residuals are smaller than predicted, meaning pessimistic weights and uncertainties that are too large. Between them the stochastic model matches the data. (n-u) x sigma0^2 weighted sum of squares Below the lower value weights pessimistic Between model matches data Above the upper value optimistic or biased Investigate systematics BEFORE rescaling weights

Figure — both tails are informative, and they mean opposite things.

σ^02=vTWvn−u\hat{\sigma}_0^2 = \frac{\mathbf{v}^{\mathsf{T}} \mathbf{W} \mathbf{v}}{n - u}

where v\mathbf{v} is the residual vector and W\mathbf{W} the weight matrix. If the stochastic model is correct — if the weights really do describe the observation uncertainties — then σ^02\hat{\sigma}_0^2 should be near one, and the scaled quantity (n−u) σ^02(n-u)\,\hat{\sigma}_0^2 follows a chi-square distribution with n−un-u degrees of freedom.

The test compares that quantity against the two-sided critical values. Both tails matter, and they mean opposite things. A value above the upper critical value means the residuals are larger than the weights predicted: either the weights are optimistic, or there is systematic error the functional model does not describe. A value below the lower critical value means the residuals are smaller than predicted, which is not good news — it means the weights were pessimistic, and every uncertainty derived from them is too large.

Complete Runnable Implementation 

from __future__ import annotations

from dataclasses import dataclass

import numpy as np
from scipy import stats


@dataclass(frozen=True)
class GlobalTest:
    """Result of the chi-square test on the a posteriori variance factor."""

    variance_factor: float
    dof: int
    statistic: float
    lower_critical: float
    upper_critical: float
    verdict: str                # "accept" | "weights_optimistic" | "weights_pessimistic"

    def describe(self) -> str:
        return (f"sigma0^2 = {self.variance_factor:.3f} on {self.dof} dof; "
                f"chi2 = {self.statistic:.1f} against "
                f"[{self.lower_critical:.1f}, {self.upper_critical:.1f}] -> "
                f"{self.verdict}")


def global_model_test(residuals: np.ndarray, weights: np.ndarray,
                      n_parameters: int, alpha: float = 0.05) -> GlobalTest:
    """Two-sided chi-square test on the variance factor.

    `weights` is the diagonal of the weight matrix (inverse variances). Both tails
    are tested because both directions are informative, and reporting only the
    upper tail hides the case where the survey was assigned uncertainties it
    comfortably beat.
    """
    residuals = np.asarray(residuals, dtype=np.float64).ravel()
    weights = np.asarray(weights, dtype=np.float64).ravel()
    if residuals.shape != weights.shape:
        raise ValueError("residuals and weights must have the same shape")
    if np.any(weights <= 0.0):
        raise ValueError("weights must be positive; a zero weight excludes the row")
    dof = residuals.size - n_parameters
    if dof <= 0:
        raise ValueError(
            f"{residuals.size} observations and {n_parameters} parameters leave "
            f"no redundancy; there is nothing to test"
        )
    vtwv = float(residuals @ (weights * residuals))
    sigma0_sq = vtwv / dof
    lo = float(stats.chi2.ppf(alpha / 2.0, dof))
    hi = float(stats.chi2.ppf(1.0 - alpha / 2.0, dof))
    if vtwv > hi:
        verdict = "weights_optimistic"
    elif vtwv < lo:
        verdict = "weights_pessimistic"
    else:
        verdict = "accept"
    return GlobalTest(sigma0_sq, dof, vtwv, lo, hi, verdict)

Parameter Reference

Name Type Units Note
residuals np.ndarray m Post-fit; one entry per observation component
weights np.ndarray m⁻² Inverse variances, not standard deviations
n_parameters int — 7 for a Helmert fit, 6 for an affine, etc.
alpha float — 0.05 gives a two-sided 95 per cent test
verdict str — Three outcomes, not two

Worked Example 

import numpy as np

rng = np.random.default_rng(11)
sigma = 0.015                                  # 15 mm, as assigned
residuals = rng.normal(0.0, sigma, size=40)
weights = np.full(40, 1.0 / sigma ** 2)

print(global_model_test(residuals, weights, n_parameters=7).describe())
# sigma0^2 = 0.907 on 33 dof; chi2 = 29.9 against [19.0, 50.7] -> accept

# The same survey with uncertainties assigned twice as optimistically:
print(global_model_test(residuals, weights * 4.0, n_parameters=7).describe())
# sigma0^2 = 3.626 on 33 dof; chi2 = 119.7 against [19.0, 50.7] -> weights_optimistic

The second call did not change the survey at all — only the claim about how good it was. That is exactly what the test detects, and it is why a failing global test is a question about the stochastic model before it is a question about the field work.

Acceptance window for the variance factor against degrees of freedom Two traces bounding the acceptance region for the a posteriori variance factor at a two-sided 95 per cent level, against degrees of freedom from 3 to 200. At 3 degrees of freedom the window runs from about 0.07 to 3.12, which accepts almost anything. By 30 it has narrowed to 0.56 to 1.57, and by 200 to 0.81 to 1.21. The test only becomes informative with real redundancy. 0.00 1.00 2.00 3.00 3 10 30 60 100 200 degrees of freedom sigma0^2 upper bound lower bound

Figure — the acceptance window for the variance factor narrows with redundancy.

Validation Check 

def assert_test_is_calibrated(n: int = 40, u: int = 7, trials: int = 2000) -> None:
    """Under a correct model, the test should reject about alpha of the time."""
    rng = np.random.default_rng(3)
    sigma = 0.015
    rejects = 0
    for _ in range(trials):
        r = rng.normal(0.0, sigma, size=n)
        w = np.full(n, 1.0 / sigma ** 2)
        if global_model_test(r, w, u).verdict != "accept":
            rejects += 1
    rate = rejects / trials
    assert 0.03 <= rate <= 0.07, f"rejection rate {rate:.3f} is not near alpha=0.05"

Common Mistakes

Rescaling the weights to force the factor to one. It always works and it destroys the diagnostic: the model now “fits” by construction, and any systematic error that caused the original failure is absorbed into an inflated uncertainty rather than investigated. Rescaling is legitimate only after the systematic explanation has been ruled out, and it should be recorded when done.

Responding to a rejected global model test Four steps in order. Screen for blunders first, because one bad observation inflates the sum of squares enough to reject a sound model. Then plot the residual field: a uniform offset, a rotation or a radial pattern each name a missing term in the functional model. Then check the units and the assignment of the weights, which is where a factor-of-a-thousand error hides. Only when all three are exhausted is rescaling the stochastic model the answer, and it is recorded when done. 1. screen for blunders first one blunder rejects a sound model 2. plot the residual field second structure names the missing term 3. check weight units third sigma vs variance is a x1000 error 4. rescale, and record it last never the first response

Figure — what to do about a rejection, in the order that finds the cause.

Testing only the upper tail. A variance factor of 0.2 passes a one-sided test comfortably and means the reported uncertainties are more than twice as large as the data supports. That is conservative in the deliverable and misleading in the report, and it usually points at observation uncertainties assigned by class rather than measured.

Applying the test with almost no redundancy. With three degrees of freedom the chi-square distribution is so wide that almost nothing is rejected, so a passing test says very little. The test earns its keep from about ten degrees of freedom upward, which is another argument for the redundancy least squares adjustment for control networks recommends.

Reporting the Test in a Deliverable

The global test earns its place in a report only if what is reported is enough to re-derive it. Four numbers and one sentence do that: the variance factor, the degrees of freedom, the two critical values, and the verdict — followed by a sentence on what was concluded and what, if anything, was changed as a result.

The last part is what a reviewer actually reads. “Variance factor 1.83 on 41 degrees of freedom, above the upper critical value; residual field inspected and found structureless; observation uncertainties revised from the assumed class values to the session estimates and the adjustment re-run” is a complete account of a decision. “Variance factor 1.02” with no context is a number that could have been produced by rescaling.

Where a test result was overridden — where the model was accepted despite a rejection — the override belongs in the report in the same words. That is not a confession; it is the normal outcome when a known systematic effect is small enough not to matter for the tolerance in play, and it is only indefensible when it is invisible.

Frequently Asked Questions

What do I do when the test rejects?

Look for systematic error before touching the weights. Plot the residual field: a uniform offset, a rotation or a radial pattern each point at a missing term in the functional model, as debugging polynomial shift residuals in GIS describes. Only when the field looks structureless is an optimistic stochastic model the likely explanation.

Should the test be run before or after outlier screening?

After, and then again if anything was removed. One blunder inflates the weighted sum of squares enough to reject a perfectly good model, so the first run mostly detects blunders. The second run, on the screened set, is the one that tests the model.

Is a variance factor of exactly one expected?

No — it is a random variable, and with thirty degrees of freedom values between about 0.6 and 1.5 are entirely ordinary. That spread is what the critical values encode, and quoting the factor without its degrees of freedom leaves a reader unable to judge it.

Does this test replace checking residuals against tolerance?

No. The chi-square test asks whether the model and the data agree with each other; the tolerance check asks whether the result meets the specification. A survey can pass the global test with a well-described but insufficient precision, and fail the tolerance — those are different conclusions with different remedies.

Can the test be applied to a transformation with no adjustment behind it?

Only if you have a stochastic model to test against, and applying a published operation accuracy to a set of control residuals is exactly that. Treat the operation’s declared accuracy as the assumed observation sigma, compute the weighted sum of squares over the control residuals, and test it on the number of control points rather than on degrees of freedom from a fit. A rejection then says the operation is not performing to its declared accuracy in your area — which is a legitimate and useful finding, and one that belongs in the deliverable rather than in a private note, since it bears directly on the uncertainty budget.