Joint mean and variance: formula and validation ledger

This ledger covers pysht.mean_variance. All six functions assume independent observations from normal populations. SHT 0.1.9 was used only as a migration audit; published formulas and independently evaluated probability laws are the correctness oracles.

Formula-source ledger

pySHT function	Primary formula basis	Independent oracle	Status
as_1samp	Arnold and Shavelle (1998), pp. 133–140; normal likelihood ratio	Literal maximized log likelihood	Public; asymptotic calibration
pn_2samp	Pearson and Neyman (1930; reprinted 1967, pp. 99–115); exact first two moments of the likelihood ratio	Log-gamma moment calculation and exact ZXC comparison	Public; lower beta tail corrected
pl_2samp	Perng and Littell (1976), pp. 968–971	Independent SciPy pooled-t and F p-values followed by Fisher’s identity	Public
muirhead_2samp	Muirhead (1982), Theorem 10.8.4, p. 370, specialized to two univariate populations	Literal second-order survival expansion and null simulation	Public; approximation limits disclosed
zxc_2samp	Zhang, Xu, and Chen (2012), pp. 180–184; exact Dirichlet probability	High-accuracy conditional-beta quadrature and direct Dirichlet simulation	Public; exact
lrt_2samp	Wilks likelihood-ratio limit applied to the same normal likelihood	Literal maximized log likelihood	Public; asymptotic calibration

Shared likelihood ratio

For samples of sizes \(n\) and \(m\), define

\[ A=\sum_i(X_i-\bar X)^2,\qquad B=\sum_j(Y_j-\bar Y)^2, \]

and

\[ C=A+B+\frac{nm}{n+m}(\bar X-\bar Y)^2. \]

The likelihood ratio for equality of both normal-population parameters is

\[ \Lambda= \frac{(A/n)^{n/2}(B/m)^{m/2}} {(C/(n+m))^{(n+m)/2}}. \]

pySHT removes a shared location before numerical scaling. When both within- group sums of squares fit comfortably on that common scale, their ratio is formed directly. When their magnitudes differ by hundreds of orders, each group instead retains its own log sum of squares and the pooled null sum is formed with logsumexp. This makes every two-sample statistic invariant to group exchange, common translation, and a change of units without erasing the smaller group’s positive variance. A fixture with one group near \(10^{-87}\) and the other near \(10^{87}\) retains the finite value \(\log\Lambda=-799.9133180008080\) even though \(\Lambda\) itself underflows.

One-sample Arnold–Shavelle likelihood ratio

For null values \((\mu_0,\sigma_0^2)\) and the maximum-likelihood variance \(\widehat\sigma^2=n^{-1}\sum_i(X_i-\bar X)^2\),

\[ -2\log\Lambda=n\left[ \log\frac{\sigma_0^2}{\widehat\sigma^2} +\frac{\widehat\sigma^2}{\sigma_0^2} +\frac{(\bar X-\mu_0)^2}{\sigma_0^2}-1 \right]\ \xrightarrow{d}\ \chi^2_2. \]

The two R entries mvar1.1998AS and mvar1.LRT reduce algebraically to this same statistic and therefore map to one Python function. The fixed fixture is checked against the literal likelihood and after a joint location-scale transformation near \(10^{100}\). An additional null-scale fixture verifies that a nonconstant sample near \(10^{-300}\) against a variance of \(10^{300}\) retains the finite log-domain likelihood-ratio statistic 8284.413760573307 and zero p-value. The earlier standardize-then-square implementation returned infinity because both the mean and variance ratios underflowed before their logarithmic contribution was formed.

Pearson–Neyman moment match

Let \(a=E_0(\Lambda)\) and \(b=\operatorname{Var}_0(\Lambda)\). The beta shapes are

\[ p=-\frac{a}{b}(a^2-a+b),\qquad q=\frac{a-1}{b}(a^2-a+b). \]

The exact gamma-function expressions for \(E_0(\Lambda)\) and \(E_0(\Lambda^2)\) are evaluated with gammaln and expm1. Because small likelihood ratios contradict the null, the p-value is the lower tail \(I_{\lambda_{\rm obs}}(p,q)\). SHT 0.1.9 used the complementary tail and also mislabeled this routine as Muirhead in its printed method.

On the fixed ledger fixture, the corrected beta p-value is 0.6917901786656631; the exact ZXC probability at the same likelihood ratio is 0.6918190488406402. The 20,000-null release simulations were:

\((n,m)\)	seed	0.01 count/rate	0.05 count/rate	0.10 count/rate	gate
(3, 3)	20260821	154 (0.00770)	949 (0.04745)	1,903 (0.09515)	Pass
(3, 3)	20260822	170 (0.00850)	938 (0.04690)	1,907 (0.09535)	Pass
(7, 8)	20260821	182 (0.00910)	997 (0.04985)	1,951 (0.09755)	Pass
(7, 8)	20260822	190 (0.00950)	962 (0.04810)	1,955 (0.09775)	Pass
(20, 20)	20260821	201 (0.01005)	964 (0.04820)	1,994 (0.09970)	Pass
(20, 20)	20260822	185 (0.00925)	970 (0.04850)	1,902 (0.09510)	Pass

Perng–Littell combination

Under the joint normal null, the two-sided pooled-t p-value \(p_t\) is independent of the two-sided variance-ratio p-value \(p_F\). Therefore

\[ -2(\log p_t+\log p_F)\sim\chi^2_4. \]

The component p-values are recomputed independently in the tests and are reported as calibration diagnostics, not parameter estimates. Both component tails and Fisher’s combination are evaluated in the log domain. On the public heterogeneous-scale fixture \(X=\{-10^{-87},10^{-87}\}\) and \(Y=\{-10^{87},10^{87}\}\), the combined p-value remains the representable 5.110888469223042e-172 instead of collapsing to zero during an intermediate F tail calculation.

Muirhead second-order approximation

With \(N=n+m\), set

\[ \rho=1-\frac{22}{24N}\left(\frac{N}{n}+\frac{N}{m}-1\right) \]

and

\[ \gamma=\frac12\left[ \left(\frac{N}{n}\right)^2+\left(\frac{N}{m}\right)^2-1 \right] -\frac{121}{96}\left(\frac{N}{n}+\frac{N}{m}-1\right)^2. \]

For \(W=-2\rho\log\Lambda\), pySHT evaluates the survival approximation

\[ \Pr_0(W\geq w)\approx \bar G_2(w)+\frac{\gamma}{\rho^2N^2} \{\bar G_6(w)-\bar G_2(w)\}, \]

where \(\bar G_k\) is a chi-square survival function. SHT used the corresponding CDF, which increases toward one under strong alternatives. The finite expansion can leave \([0,1]\) in extreme tails, so pySHT truncates the result and reports \(\rho\) and the second-order coefficient as diagnostics.

This approximation is not advertised for very small samples. The complete passing \((n,m)=(20,20)\) and failing \((7,8)\) audits appear in the explicit tables below. Use zxc_2samp for exact small-sample inference.

Exact Zhang–Xu–Chen probability

Under the joint null, the normalized within-group and between-group sums of squares have a Dirichlet distribution with parameters

\[ \left(\frac{n-1}{2},\frac{m-1}{2},\frac12\right). \]

The exact p-value is the probability of a likelihood ratio no larger than the observed ratio. pySHT conditions on one Dirichlet coordinate, reducing the published two-dimensional region integral to beta tail masses plus one one-dimensional conditional-beta integral. Roots are solved in logit coordinates. Each half of the interior integrand is divided by its own log maximum before quadrature, and the four probability pieces are recombined by logsumexp. The smaller sample is assigned to the conditional coordinate, which makes group symmetry exact in floating point and avoids a thin unbalanced-sample boundary layer.

The fixed high-accuracy integral is 0.6918190488406402; the legacy 500-by-500 nested rule gives 0.6921061. A direct ten-million-draw Dirichlet calculation gave 0.6917072 with Monte Carlo standard error 0.000146, independently corroborating the stable integral.

Three independent far-tail fixtures exercise behavior that ordinary probability-space quadrature misses: \(\log\Lambda=-50\) with \((n,m)=(2,2)\) gives 1.2412994009417087e-10; the same log likelihood with \((n,m)=(20,200)\) gives 2.2176768591905327e-21; and \(\log\Lambda=-800\) with \((n,m)=(2,2)\) gives the still-representable 2.4572424773752458e-172. These values come from separate high-precision log-coordinate integrations.

Explicit asymptotic size gates

For every row below, 20,000 independent standard-normal null datasets were generated for each of two seeds. The one-sample null was \(N(0,1)\) with the default null mean and variance; each two-sample row used two independent \(N(0,1)\) samples of the displayed sizes. For each method, scenario, and seed, a fresh numpy.random.default_rng(seed) generated observations in sample and replication order, the public function was called once, and rejection meant result.pvalue < alpha. A row passes at all \(\alpha\in\{0.01,0.05,0.10\}\) when

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

Calibration	Scenario	Seed	0.01 count/rate	0.05 count/rate	0.10 count/rate	Gate
as_1samp, \(\chi^2_2\)	\(n=50\)	20260819	204 (0.01020)	1,029 (0.05145)	2,078 (0.10390)	Pass
as_1samp, \(\chi^2_2\)	\(n=50\)	20260820	221 (0.01105)	1,022 (0.05110)	2,036 (0.10180)	Pass
lrt_2samp, \(\chi^2_2\)	\(n=m=75\)	20260819	184 (0.00920)	1,016 (0.05080)	2,051 (0.10255)	Pass
lrt_2samp, \(\chi^2_2\)	\(n=m=75\)	20260820	204 (0.01020)	1,002 (0.05010)	2,082 (0.10410)	Pass
muirhead_2samp	\(n=m=20\)	20260819	210 (0.01050)	1,030 (0.05150)	2,065 (0.10325)	Pass
muirhead_2samp	\(n=m=20\)	20260820	206 (0.01030)	1,063 (0.05315)	2,050 (0.10250)	Pass

Smaller scenarios are not claimed as release-validated. These replayable rows show the failures rather than hiding them:

Calibration	Scenario	Seed	0.01 count/rate	0.05 count/rate	0.10 count/rate
as_1samp, \(\chi^2_2\)	\(n=20\)	20260819	251 (0.01255)	1,192 (0.05960)	2,268 (0.11340)
as_1samp, \(\chi^2_2\)	\(n=20\)	20260820	244 (0.01220)	1,121 (0.05605)	2,192 (0.10960)
lrt_2samp, \(\chi^2_2\)	\(n=m=20\)	20260819	263 (0.01315)	1,203 (0.06015)	2,358 (0.11790)
lrt_2samp, \(\chi^2_2\)	\(n=m=20\)	20260820	260 (0.01300)	1,241 (0.06205)	2,352 (0.11760)
muirhead_2samp	\((n,m)=(7,8)\)	20260819	474 (0.02370)	1,441 (0.07205)	2,556 (0.12780)
muirhead_2samp	\((n,m)=(7,8)\)	20260820	435 (0.02175)	1,374 (0.06870)	2,492 (0.12460)

The functions remain mathematically defined in these smaller designs, but their reported p-values must be interpreted as rough asymptotic approximations; zxc_2samp supplies the exact two-sample alternative.

Targeted alternative-power audit

At nominal 0.05, every public routine was run on 2,000 strong alternatives. For as_1samp, observations were \(N(0.75,1.5^2)\) and the tested null was the default \((0,1)\) at \(n=50\). Each two-sample row compared \(N(0,1)\) with \(N(1,2^2)\); the sample size stays in the passing calibration regime of the corresponding approximation. A fresh numpy.random.default_rng(20260829) was reset for every row, generated \(x\) then \(y\) in each replication, and the public function was called once.

Procedure	Design	Rejections/2,000	Rate
as_1samp	\(n=50\)	2,000	1.0000
pn_2samp	\(n=m=20\)	1,810	0.9050
pl_2samp	\(n=m=20\)	1,818	0.9090
muirhead_2samp	\(n=m=20\)	1,812	0.9060
zxc_2samp	\(n=m=20\)	1,810	0.9050
lrt_2samp	\(n=m=75\)	2,000	1.0000

The near agreement among PN, Muirhead, and exact ZXC in the common design is an additional directional check; no table entry asserts a uniform lower power bound over all joint alternatives.

Migration mapping

pySHT	SHT 0.1.9	Deliberate change
as_1samp	mvar1.1998AS, mvar1.LRT	Duplicate formulas share one function; null variance keyword is variance
pn_2samp	mvar2.1930PN	Lower beta tail, log-gamma moments, corrected method label
pl_2samp	mvar2.1976PL	Stable explicit component calculations
muirhead_2samp	mvar2.1982Muirhead	Survival tail and explicit truncation
zxc_2samp	mvar2.2012ZXC	Stable conditional-beta integral
lrt_2samp	mvar2.LRT	Log-domain likelihood ratio

Primary sources: Arnold and Shavelle (1998), Pearson and Neyman (1930), Perng and Littell (1976), Muirhead’s Aspects of Multivariate Statistical Theory (1982), and Zhang, Xu, and Chen (2012).