Lee–You–Lin maximum pairwise Bayes factor for means

Status: Equation (4), log-domain, translation, scaling, group-exchange, large-magnitude, and degenerate-feature gates pass.

For each feature, lyl_2samp uses the ordinary centered maximum-likelihood variances from Equation (4). The published model is Gaussian with a common covariance matrix under the two samples. At a0=b0=0, its component log Bayes factor is

0.5 log(gamma / (1 + gamma)) + n / 2 * log(RSS_pooled / RSS_within).

The maximum component log Bayes factor is returned without exponentiation and without a fabricated p-value or universal cutoff. All component values remain available on the immutable result.

By default, alpha=2.01 and gamma=max(n_x+n_y, p)**(-alpha), matching the paper’s numerical choice and dimension-dependent rate. Supplying gamma explicitly overrides that rule and is recorded in diagnostics.

Audit finding

The legacy C++ helper weakened centering by dividing the squared total by 1 + gamma. That made the result change after adding the same constant to both groups and does not match Equation (4). pySHT uses exact centered residual sums; gamma appears only in the prior penalty. A shared inverse-gamma a0,b0 extension is retained for compatibility and reduces algebraically to Equation (4) at zero. It is explicitly labeled as a pySHT extension in the result diagnostics and documentation, not attributed to Equation (4). Scaling b0 by the square of the data scale preserves the result.

The two samples use a shared deterministic feature-wise anchor before working-scale selection. This preserves group exchange and translation invariance at large common locations without changing the log Bayes factors.

Targeted alternative-evidence gate

Because lyl_2samp returns a Bayes factor rather than a p-value, this gate does not invent a rejection cutoff. Seed 2026090322 initializes a SeedSequence; child 0 drives one persistent PCG64 data stream and child 1 is reserved as an independent auxiliary stream. In each of 1,000 replications, child 0 draws x and then y as independent 30-by-60 standard-normal matrices. The matched null value is mean.lyl_2samp(x, y).statistic; the alternative value is mean.lyl_2samp(x, y + shift).statistic, where shift equals 1.0 in the first three coordinates and zero elsewhere. The same realized (x,y) is used in each matched pair, and the default published gamma rate is retained.

Replications	Null median maximum log BF	Alternative median maximum log BF	Alternative exceeds matched null
1,000	-0.817807	5.961739	992 / 1,000 = 0.992

The paired ordering and median displacement establish directional evidence under this sparse strong alternative; neither number is proposed as a universal evidence threshold. The audit is reproduced by python -m tools.mean_power_audits mean.lyl_2samp under Python 3.12.13, NumPy 2.5.1, and SciPy 1.18.0.

Primary reference: K. Lee, K. You, and L. Lin, Bayesian Optimal Two-Sample Tests for High-Dimensional Gaussian Populations, Bayesian Analysis 19 (2024), 869–893, doi:10.1214/23-BA1373.