Covariance tests: formula and validation ledger#
This ledger covers the seven public functions in pysht.covariance and one
withheld Fisher implementation. The
pinned SHT 0.1.9 source was inspected for migration, but primary-paper
equations, intentionally literal implementations, algebraic invariants, and
null simulations are the correctness oracles.
Common numerical policy#
Except for the paper-specific Lee–You–Lin known-zero-mean model, each group has a feature-wise anchor removed in the original coordinates before centering and division by one common finite scale. An opposite-sign overflow fallback performs the subtraction in normalized coordinates. Covariance products are therefore evaluated without avoidable overflow, while standardized statistics remain unchanged. Tests cover separately shifted locations as large as \(10^{14}\) while preserving exactly representable within-group variation, common scales through \(10^{100}\), group exchange, group order, and feature permutation where the method implies those invariances. Undefined zero variance estimates and singular matrices are rejected before a p-value is formed.
Fisher one-sample test (withheld)#
fisher_1samp is not part of the 0.1.0 public API. The private
_fisher_1samp implementation is retained for continued validation. For \(N\)
observations, write \(n=N-1\), \(c=p/n\), and let \(S\) be the unbiased
sample covariance after centering and whitening by
\(\Sigma_0^{-1/2}\). Fisher’s unbiased spectral-moment estimators
\(a_1,\ldots,a_4\) are inserted into
and
Both use an upper standard-normal tail. A fixed nonidentity-null fixture is checked against an independent trace-power implementation for both variants. An invertible feature transformation applied simultaneously to the data and null covariance leaves the result unchanged. When the fourth-order polynomial lies outside float64, normalized trace estimators and a signed log-sum evaluation preserve its limiting sign instead of leaking an arithmetic overflow. The sign matters because an unbiased fourth-moment estimator need not be nonnegative in a very small finite sample, even though its population target is nonnegative.
Legacy correction. SHT calculated whitened observations but passed the unwhitened sample covariance to the moment estimator. Consequently its stated nonidentity null was not tested. pySHT computes every \(a_r\) from the whitened covariance.
The normal approximation is finite-sample sensitive. A historical optimized scatter simulation at \((N,p)=(201,800)\) was not retained as a reproducible public-call runner. Fresh public-path probes at the same aspect ratio but \((N,p)=(51,200)\) and \((101,400)\) failed calibration for at least one variant and level. Consequently neither null calibration nor targeted power has met the project’s reproducibility gate, and exposing the function would violate the 0.1.0 release policy. No historical rate is used as release evidence.
This caution also has primary-author support independent of the failed pySHT audit. Fisher’s dissertation, Section 2.3, derives the identity-covariance construction; on printed page 27 it says the fourth-moment estimator’s large variance can hinder the proposed statistic and that more work is necessary. The conclusion on printed page 28 says the construction could not then be justified. Those are PDF pages 36–37 in the repository copy. See Fisher (2009).
Primary source: Fisher [1].
Wu–Li random-projection tests#
For each Gaussian unit vector \(R_r\), the one-sample statistic uses
The two-sample statistic uses the variance-stabilized log ratio
Each group’s centered projections are normalized separately and their scalar units are restored after taking logarithms. Consequently, a finite log variance ratio remains evaluable even when forming \(s_{1r}^2/s_{2r}^2\) directly would overflow or one group’s squared working values would underflow.
Because covariance equality has an unrestricted alternative, pySHT reports \(M=\max_r|Z_r|\) and
The survival probability is evaluated from the small marginal tail
\(2\Phi(-M)\) with log1p/expm1. This avoids rounding
\(2\Phi(M)-1\) to one, which would otherwise report a zero p-value around
\(M=9\) even though the maximum-test tail remains representable.
This is a genuine two-sided maximum law under the paper’s asymptotic independence argument. The paper derives its procedure for known zero means; pySHT deliberately extends it to unknown means by centering and using \(N-1\) rather than \(N\) degrees of freedom. The two 20,000-run audits below validate that extension in the stated Gaussian regimes. Seed replay, global RNG isolation, group exchange, both low- and high-variance alternatives, and the literal maximum-normal probability are tested.
Legacy correction. SHT returned only max(Z_r) with a one-sided tail
while describing a two-sided equality test. That loses power when projected
variances are smaller, and the two-sample answer can change after exchanging
the groups. pySHT uses max(abs(Z_r)) and the corresponding two-sided tail.
Primary source: Wu and Li [2].
Fresh 20,000-replication public-call audits use the scenario keys and stream contract recorded below. The two-sample method failed at \((N_1,N_2,p,m)=(100,120,30,50)\) and is therefore advertised only at the larger validated sample-size regime.
scenario key |
seed |
rejection counts at 0.01, 0.05, 0.10 |
rates |
|---|---|---|---|
|
20260831 |
222, 987, 1963 |
0.01110, 0.04935, 0.09815 |
|
20260901 |
205, 990, 2019 |
0.01025, 0.04950, 0.10095 |
|
20260831 |
233, 1030, 2074 |
0.01165, 0.05150, 0.10370 |
|
20260901 |
214, 1042, 2033 |
0.01070, 0.05210, 0.10165 |
Li–Chen two-sample test#
The full path evaluates the mutually-distinct-index U-statistics
and Equation (2.2)’s four-term \(C_{n_1n_2}\). The unstandardized statistic is
Under the null, the paper’s standard-deviation estimator is
An independent four-loop fixture reproduces every ordered sum. The optimized implementation reduces the computation to Gram-matrix row sums without downgrading the U-statistic.
A former unbiased=False centered leading-term shortcut is deliberately not
part of the 0.1.0 signature. In two independent 2,000-replication Gaussian
null probes at \((n_1,n_2,p)=(30,30,50)\) (integer seeds 20260831 and 20260901,
one reset stream per seed), its nominal 0.05 rejection rates were 0.2185 and
0.2095 (437/2000 and 419/2000). That branch therefore failed before the
release-scale gate; pySHT exposes only the literal U-statistic path that passed
the 20,000-run audits.
Legacy correction and typesetting audit. SHT divided \(T\) by \(\sqrt{\widehat\sigma_0}\). Equation (2.7), physical units, scale invariance, and the authors’ published reference code all identify the displayed linear combination as the estimated standard deviation, despite an inconsistent square in nearby prose. pySHT divides by it directly. At the fixed ledger fixture, the corrected and square-rooted values differ materially.
A fresh public-call normal-null audit with \((n_1,n_2,p)=(30,30,50)\) produced
counts (190, 1031, 2060) and rates (0.00950, 0.05155, 0.10300) for seed
20260831, and counts (227, 1053, 2082) and rates
(0.01135, 0.05265, 0.10410) for seed 20260901. Both 20,000-run streams pass.
Primary source: Li and Chen [3].
Cai–Liu–Xia maximum test#
Let \(\widehat\sigma_{k,ij}\) be the maximum-likelihood sample covariance and
pySHT reports
with the published type-I extreme-value survival probability. The tail is
evaluated with expm1 and a log rate, avoiding cancellation near zero. A
literal double-loop fixture and all group, location, scale, and feature-order
invariances pass. Fresh 20,000-run public-call normal-null streams at
\((n_1,n_2,p)=(100,100,30)\) produced counts (163, 994, 2037) and rates
(0.00815, 0.04970, 0.10185) for seed 20260831, and counts
(176, 958, 2044) and rates (0.00880, 0.04790, 0.10220) for seed 20260901.
Both pass, while the output remains correctly labeled asymptotic.
Primary source: Cai et al. [4].
Lee–You–Lin maximum pairwise Bayes factor#
For every ordered pair \(i\ne j\), the method compares equality of the two conditional regressions under the paper’s known-zero-mean Gaussian model. Writing \(X_i\) and \(X_j\) for the raw variable columns, pySHT calculates
with ordinary least squares and no intercept, and inserts the two group
residuals and the pooled residual into Equations (12)–(14). No automatic
centering is performed. gammaln, logarithms, pairwise scaling, and logaddexp
are used throughout. Equation (15)’s maximum is retained in log form; the full
ordered matrix is immutable, its diagonal is -inf, and no frequentist p-value
or automatic threshold is invented.
An independent implementation of Equations (12)–(15) reproduces every
component, including the recommended a0=b0=0.01 and
gamma=max(n_1+n_2,p)^(-2.01) defaults. Tests cover group exchange, feature
permutation, explicit gamma overrides, extreme float64 scales, and
simultaneous conversion of b0 with squared data units. Translation invariance
is intentionally not claimed because it would contradict the known-zero-mean
model.
Legacy corrections. The old C++ kernel cast log arguments to float and
weakened the OLS projection by (1+gamma). The primary equations put gamma
only in the prior penalty. pySHT uses ordinary OLS residuals and float64
log-domain arithmetic.
Primary source: Lee et al. [5].
Schott multi-sample tests#
For the 2001 Wald test, let \(n_i=N_i-1\), \(S=\sum_i n_iS_i/\sum_i n_i\), and \(A_i=S^{-1/2}S_iS^{-1/2}\). The implemented nonnegative form is
It is algebraically identical to the paper’s double-trace expression and is checked against that literal form. A singular pooled covariance is outside the Wald statistic’s domain. SHT silently fell back to a pseudoinverse; pySHT rejects it.
The 2007 test uses the paper’s pairwise bias-corrected Frobenius statistic, pooled trace estimator, and normalizing \(\theta\). A separate transcription of every pairwise term agrees with the optimized implementation. Group order, separate locations, and common scale do not affect either Schott result.
The fresh 20,000-replication public-call Gaussian-null audit was:
scenario key |
seed |
rejection counts at 0.01, 0.05, 0.10 |
rates |
|---|---|---|---|
|
20260831 |
193, 1003, 2005 |
0.00965, 0.05015, 0.10025 |
|
20260901 |
183, 901, 1959 |
0.00915, 0.04505, 0.09795 |
|
20260831 |
243, 1100, 2110 |
0.01215, 0.05500, 0.10550 |
|
20260901 |
227, 1045, 2095 |
0.01135, 0.05225, 0.10475 |
Every entry satisfies the project’s release tolerance. A preliminary \(p=30\) audit missed the 0.05 boundary by two rejections for one seed, so that smaller-dimensional finite regime is not advertised as calibrated.
Reproducible stream and power evidence#
python -m tools.covariance_release_audits --all resets each scenario/seed
pair, constructs SeedSequence(seed), spawns the data stream first and the
auxiliary-randomness stream second, and advances both in replication order.
The runs above used Python 3.12.13, NumPy 2.5.1, SciPy 1.18.0, and pySHT
0.1.0. Every displayed count is from the public function named by its
scenario, not from a duplicated optimized formula.
Targeted alternatives used 2,000 replications, integer seed 20260902, and the
same reset/split policy with SeedSequence([seed, method_index]). The table
reports rejection counts and rates at level 0.05.
public function |
fully specified alternative |
count |
rate |
|---|---|---|---|
|
Gaussian, \(N=100,p=30,\Sigma=1.5I\) |
2000 |
1.0000 |
|
Gaussian, \(N_1=300,N_2=360,p=30,\Sigma_1=I,\Sigma_2=1.3I\) |
1997 |
0.9985 |
|
Gaussian, \(n_1=n_2=30,p=50,\Sigma_1=I,\Sigma_2=2I\) |
1799 |
0.8995 |
|
Gaussian, \(n_1=n_2=100,p=30\); second-group first SD 1.8 |
1817 |
0.9085 |
|
Gaussian, \(g=3,N_i=100,p=3\); third-group SD \((1.6,1,1)\) |
1947 |
0.9735 |
|
Gaussian, \(g=3,N_i=100,p=50\); third-group covariance \(1.5I\) |
1968 |
0.9840 |
For lyl_2samp, the same seed and 2,000 paired draws used
\(n_1=n_2=50,p=10\) and correlation 0.7 between the first two coordinates only
in the alternative second group. The median maximum log Bayes factor moved
from -5.38217 under the null comparator to 1.66247 under the alternative; the
alternative exceeded its paired null value in 1919/2000 draws (0.9595). This
is evidence direction, not a fabricated frequentist rejection rule.
Legacy mapping#
pySHT |
SHT 0.1.9 |
Deliberate change |
|---|---|---|
|
|
Corrected implementation retained privately; reproducible release gates incomplete |
|
|
Genuine two-sided maximum and local RNG |
|
|
Literal U-statistics; divide by estimated SD, not its square root |
|
|
Stable extreme-value tail and strict domains |
|
|
Genuine two-sided, group-symmetric maximum |
|
|
Published OLS residuals, log output, no fabricated p-value |
|
|
Reject singular pooled covariance; no pseudoinverse fallback |
|
|
Stable common scaling and literal correction factors |