Joint mean and covariance: formula and validation ledger

This ledger covers pysht.mean_covariance. Each alternative is the complement of the joint null: a departure in the mean, covariance, or both is evidence against the null. SHT 0.1.9 was used only as migration evidence; published formulas, literal independent implementations, and null simulations are the scientific gates.

Formula-source ledger

pySHT function	Primary formula basis	Independent oracle	Status
llzs_1samp	Liu, Liu, Zheng, and Shi (2017), Section 2, pp. 84–87	Literal raw-second-moment formula and high-dimensional null simulations	Public; legacy centering and tail corrected
lrt_1samp	Multivariate normal likelihood and Wilks’ fixed-dimension limit	Original-coordinate solve and log-determinant identity	Public; fixed-dimensional asymptotic calibration
hn_2samp	Hyodo and Nishiyama (2018), Sections 2.1–2.2; equations (2.1)–(2.2) in the open technical report	Literal trace-estimator transcription and high-dimensional null simulations	Public

Stable null whitening

For the one-sample methods, define

\[ Z_i=\Sigma_0^{-1/2}(X_i-\mu_0). \]

pySHT obtains \(Z_i\) with a Cholesky solve, never by forming an inverse or matrix square root. Tests generate \(X_i=\mu_0+LZ_i\) with \(\Sigma_0=LL^\mathsf{T}\) and recover the identity-null statistic and p-value to floating-point tolerance. A supplied null covariance must be finite, symmetric, and positive definite.

LLZS one-sample high-dimensional test

For \(n\) whitened observations in \(p\) dimensions, let

\[ \bar Z=n^{-1}\sum_i Z_i,\qquad M_n=n^{-1}\sum_i Z_iZ_i^\mathsf{T}, \]

and define

\[ T_L=\|\bar Z\|^2+\operatorname{tr}(M_n-I_p)^2. \]

With \(y=p/n\) and marginal excess-kurtosis parameter \(\beta\), Liu, Liu, Zheng, and Shi derive

\[ \mu_0=y(p+\beta+2),\qquad \sigma_0^2=4y^2\{y(2+\beta)+1\}, \]

and \((T_L-\mu_0)/\sigma_0\Rightarrow N(0,1)\). pySHT uses the paper’s moment estimator

\[ \widehat\beta=(np)^{-1}\sum_{i,j}Z_{ij}^4-3 \]

and the upper normal tail. The population targets behind both squared departures are nonnegative, so alternatives move the statistic to the right. The realized aspect ratio and estimated kurtosis are calibration diagnostics, not scientific parameter estimates returned by the procedure.

When \(p>n\), the implementation evaluates \(\operatorname{tr}(M_n^2)\) through the \(n\)-by-\(n\) observation Gram matrix. It therefore does not allocate a \(p\)-by-\(p\) matrix in the regime for which the method was designed. A separate fixture compares this route with the literal feature-space formula.

SHT 0.1.9 replaced \(M_n\) with the covariance matrix centered at the sample mean while retaining the null center derived for \(M_n\). When \(p/n\to y>0\), that changes the statistic by an order-one amount. It also used a two-sided normal p-value. The fixed regression fixture distinguishes raw and centered second moments and checks the corrected upper tail.

The advertised regime has independent rows, the paper’s coordinate and spectral regularity, a finite fourth moment, and \(p/n\) bounded away from zero and infinity. The explicit release gates below cover Gaussian, bounded non-Gaussian, and finite-fourth-moment heavy-tailed coordinates.

Classical one-sample likelihood-ratio test

Let \(\bar Z\) be the whitened sample mean and

\[ S_n=n^{-1}\sum_i(Z_i-\bar Z)(Z_i-\bar Z)^\mathsf{T}. \]

The normal likelihood gives

\[ -2\log\Lambda=n\left\{ \operatorname{tr}(S_n)-\log|S_n|-p+\|\bar Z\|^2 \right\} \ \xrightarrow{d}\ \chi^2_{p(p+3)/2}. \]

The implementation requires \(n>p\) and a positive-definite fitted covariance; it never substitutes a pseudoinverse. Eigenvalue contributions are evaluated as \(\lambda-1-\log1p(\lambda-1)\) to preserve accuracy near the null. The fixed-data oracle computes the same likelihood directly in the original coordinates with solves and determinant ratios.

This is a multivariate-normal, fixed-dimension limit. It is not a fallback for a singular or proportional high-dimensional design.

Hyodo–Nishiyama two-sample test

For group \(g\in\{1,2\}\), let \(S_g\) be the unbiased sample covariance and

\[ K_g=(n_g-1)^{-1}\sum_i\|X_{gi}-\bar X_g\|^4. \]

The unbiased trace-square estimator is

\[ A_g= \frac{n_g-1}{n_g(n_g-2)(n_g-3)} \left[ (n_g-1)(n_g-2)\operatorname{tr}(S_g^2) +\operatorname{tr}(S_g)^2-n_gK_g \right]. \]

With \(C=\operatorname{tr}(S_1S_2)\), the squared-distance estimators are

\[ \widehat d^2= \|\bar X_1-\bar X_2\|^2 -\frac{\operatorname{tr}S_1}{n_1} -\frac{\operatorname{tr}S_2}{n_2}, \]

\[ \widehat D^2=A_1+A_2-2C. \]

Under the joint null, define

\[ V_1= \frac{2A_1}{n_1^2}+\frac{2A_2}{n_2^2} +\frac{4C}{n_1n_2}, \]

\[ V_2= \frac{4A_1^2}{n_1^2}+\frac{4A_2^2}{n_2^2} +\frac{8C^2}{n_1n_2}. \]

Hyodo and Nishiyama’s statistic is

\[ T=\frac{\widehat d^2}{\sqrt{V_1}} +\frac{\widehat D^2}{\sqrt{V_2}}, \qquad \frac{T}{\sqrt2}\Rightarrow N(0,1), \]

with an upper-tail p-value. At least four rows per group are needed to define \(A_g\), and nonpositive finite-sample variance estimates are rejected.

The statistic is evaluated after a common translation and scaling to protect fourth-order arithmetic. The two distance estimates are converted back to the original squared and fourth-power measurement units. If an original-unit estimate lies outside finite float64 range, that estimate is omitted rather than mislabeled in computational units.

Trace products are evaluated adaptively in feature space when \(p\) is small and through observation Gram matrices when \(p\) is large. Both paths are checked against literal covariance-matrix calculations.

Tests independently transcribe the paper’s equations and cover group exchange, common translation and scaling, orthogonal feature transformations, and alternatives that change only the mean or only the covariance. The advertised null limit requires uniformly bounded eighth moments, the paper’s mixed-moment factorization for distinct latent coordinates, both group sizes and dimension increasing, and assumptions A1–A2 on sample and trace growth.

Explicit asymptotic size gates

Every release row uses 20,000 null datasets and is assessed simultaneously at \(\alpha\in\{0.01,0.05,0.10\}\) against

\[ |\widehat\alpha-\alpha| \leq \max\left(0.005,\, 4\sqrt{\frac{\alpha(1-\alpha)}{20000}}\right). \]

LLZS proportional-growth regime

scenario key	seed	counts at 0.01, 0.05, 0.10	rates	gate
mean-covariance.llzs.normal-n100-p50	20260831	235, 1047, 1964	0.01175, 0.05235, 0.09820	Pass
mean-covariance.llzs.normal-n100-p50	20260901	255, 1044, 1972	0.01275, 0.05220, 0.09860	Pass
mean-covariance.llzs.normal-n100-p100	20260831	198, 1007, 2009	0.00990, 0.05035, 0.10045	Pass
mean-covariance.llzs.normal-n100-p100	20260901	229, 989, 1917	0.01145, 0.04945, 0.09585	Pass
mean-covariance.llzs.normal-n100-p200	20260831	183, 963, 1961	0.00915, 0.04815, 0.09805	Pass
mean-covariance.llzs.normal-n100-p200	20260901	195, 958, 1939	0.00975, 0.04790, 0.09695	Pass
mean-covariance.llzs.t8-n100-p100	20260831	193, 938, 1871	0.00965, 0.04690, 0.09355	Pass
mean-covariance.llzs.t8-n100-p100	20260901	196, 987, 1915	0.00980, 0.04935, 0.09575	Pass
mean-covariance.llzs.uniform-n100-p100	20260831	221, 995, 1974	0.01105, 0.04975, 0.09870	Pass
mean-covariance.llzs.uniform-n100-p100	20260901	222, 1030, 2018	0.01110, 0.05150, 0.10090	Pass

The smaller \((n,p)=(50,25)\) Gaussian scenario failed one seed at level 0.01 with rejection rate 0.01515; pySHT therefore does not claim broad small-sample calibration from the asymptotic LLZS law.

Fixed-dimensional likelihood-ratio regime

scenario key	seed	counts at 0.01, 0.05, 0.10	rates	gate
mean-covariance.lrt.normal-n200-p2	20260831	200, 1021, 2023	0.01000, 0.05105, 0.10115	Pass
mean-covariance.lrt.normal-n200-p2	20260901	222, 1042, 2056	0.01110, 0.05210, 0.10280	Pass
mean-covariance.lrt.normal-n300-p3	20260831	201, 997, 2021	0.01005, 0.04985, 0.10105	Pass
mean-covariance.lrt.normal-n300-p3	20260901	213, 1035, 2041	0.01065, 0.05175, 0.10205	Pass
mean-covariance.lrt.normal-n500-p5	20260831	194, 1041, 2084	0.00970, 0.05205, 0.10420	Pass
mean-covariance.lrt.normal-n500-p5	20260901	201, 1030, 2041	0.01005, 0.05150, 0.10205	Pass

At \((n,p)=(50,5)\), the uncorrected Wilks calibration yielded approximately \((0.015,0.071,0.131)\) and failed. The documentation therefore treats lrt_1samp as genuinely fixed-dimensional and large-sample.

HN high-dimensional regime

scenario key	seed	counts at 0.01, 0.05, 0.10	rates	gate
mean-covariance.hn.normal-n100-n100-p200	20260831	215, 1083, 2055	0.01075, 0.05415, 0.10275	Pass
mean-covariance.hn.normal-n100-n100-p200	20260901	236, 1030, 1992	0.01180, 0.05150, 0.09960	Pass
mean-covariance.hn.normal-n150-n180-p300	20260831	261, 1074, 2084	0.01305, 0.05370, 0.10420	Pass
mean-covariance.hn.normal-n150-n180-p300	20260901	240, 1008, 2016	0.01200, 0.05040, 0.10080	Pass
mean-covariance.hn.uniform-n150-n180-p300	20260831	251, 1054, 2050	0.01255, 0.05270, 0.10250	Pass
mean-covariance.hn.uniform-n150-n180-p300	20260901	222, 1047, 2034	0.01110, 0.05235, 0.10170	Pass

The smaller unequal scenario \((100,120,200)\) narrowly failed one of two seeds at level 0.05 (0.05655 versus a tolerance endpoint of about 0.05616). This is evidence of finite-sample convergence error, not a formula substitution. The public method is documented for the larger high-dimensional regime rather than being described as a generic procedure for all arrays with four rows.

Stream contract and targeted power

All fresh tables above are produced by python -m tools.covariance_release_audits --all. Each scenario/seed pair resets SeedSequence(seed), spawns the data stream first and auxiliary stream second, and advances in replication order. The runs used Python 3.12.13, NumPy 2.5.1, SciPy 1.18.0, and pySHT 0.1.0. Every p-value comes from the public function itself.

Targeted alternatives used 2,000 replications and integer seed 20260902, with SeedSequence([seed, method_index]) reset for each function. At level 0.05:

public function	fully specified alternative	count	rate
llzs_1samp	Gaussian, \(n=100,p=100\); every mean coordinate shifted by 0.2	2000	1.0000
lrt_1samp	Gaussian, \(n=300,p=3\); mean shift \((0.3,0,0)\)	1945	0.9725
hn_2samp	Gaussian, \(n_1=150,n_2=180,p=300\); every second-group mean coordinate shifted by 0.1	2000	1.0000

Run python -m tools.covariance_power_audits to reproduce these counts. This is targeted power evidence for the advertised regimes, not a claim that the same power holds for smaller samples or arbitrary alternatives.

Migration mapping

pySHT	SHT 0.1.9	Deliberate change
llzs_1samp	sim1.2017Liu	Raw second moment, upper tail, correct logical alternative
lrt_1samp	sim1.LRT	Stable whitening, strict SPD/rank domain, correct logical alternative
hn_2samp	sim2.2018HN	Stable scaling, original-unit estimates, correct logical alternative

Primary sources are Liu, Liu, Zheng, and Shi (2017) and Hyodo and Nishiyama (2018). The open Hyodo–Nishiyama technical report provides the equation-numbered oracle used in the independent tests.