[2] Tests for Multivariate Mean

The functions in pysht.mean cover classical low-dimensional inference, high-dimensional trace and diagonal tests, unequal-covariance procedures, sparse alternatives, randomized dimension reduction, and multi-group designs. Observations are rows and variables are columns.

Choosing a procedure

Design and working assumptions	Function	R routine
One sample; Gaussian; invertible covariance	hotelling_1samp	mean1.1931Hotelling
Two samples; Gaussian; common invertible covariance	hotelling_2samp	mean2.1931Hotelling
One or two samples; Gaussian; covariance may be singular	dempster_1samp, dempster_2samp	mean1.1958Dempster, mean2.1958Dempster
One or two samples; dense high-dimensional shift; common covariance for two samples	bs_1samp, bs_2samp	mean1.1996BS, mean2.1996BS
One or two samples; diagonal standardization; common covariance for two samples	sd_1samp, sd_2samp	mean1.2008SD, mean2.2008SD
Two samples; dense high-dimensional shift; unequal covariance permitted	cq_2samp	–
One, two, or several samples; fixed small sample sizes and diverging dimension	li_1samp, li_2samp, li_ksamp	–
Two samples; unequal covariance; more rows than variables	yao_2samp, johansen_2samp, nvm_2samp, ky_2samp	mean2.1965Yao, mean2.1980Johansen, mean2.1986NVM, mean2.2004KY
Two samples; high dimension; one Gaussian projection	ljw_2samp	mean2.2011LJW
Two samples; high dimension; averaged random subspaces	thulin_2samp	mean2.2014Thulin
Two samples; sparse shift; Bayesian evidence	maximum_pairwise_bayes_factor_2samp	mean2.mxPBF
Several groups; common covariance	schott_ksamp	meank.2007Schott
Several groups; unequal covariance; Gaussian Scheffé transformation	zx_ksamp	meank.2009ZX
Several groups; unequal covariance; high-dimensional factor model	cph_ksamp	meank.2019CPH

An asymptotic p-value is not a guarantee that a test is suitable for a small dataset. Check the dimensional regime and covariance assumptions in the method ledger before interpreting the result. The reported alternative describes the scientific direction: small p-values are evidence that at least one mean vector differs. maximum_pairwise_bayes_factor_2samp instead returns maximum log Bayes-factor evidence and deliberately has no p-value.

For numerical stability, null vectors are removed before scaling in one-sample routines. Translation-invariant two- and multi-sample routines remove one deterministic feature-wise anchor, shared by every group, in the original input coordinates before choosing a global or per-feature scale. For procedures whose definitions do not split or pair rows, this also makes the deterministic arithmetic insensitive to row or group order and avoids letting a huge common location erase representable variation. The documented CPH split estimator and Zhang–Xu Scheffé transformation are exceptions because row splitting or pairing is part of their definitions. No numerical method can recover information already lost when values were rounded to float64.

Classical and unequal-covariance tests

pysht.mean.hotelling_1samp(x, *, popmean=None)

Perform the exact one-sample Hotelling \(T^2\) test.

Exact F calibration requires multivariate normal observations, n > p, and a positive-definite sample covariance matrix.

Parameters:

• x (ArrayLike)

• popmean (ArrayLike | None)

Return type:

HypothesisTestResult

pysht.mean.hotelling_2samp(x, y, *, paired=False)

Perform the exact equal-covariance two-sample Hotelling test.

Independent samples are assumed to share a positive-definite covariance matrix. With paired=True, the exact one-sample test is applied to the row-wise differences.

Parameters:

• x (ArrayLike)

• y (ArrayLike)

• paired (bool)

Return type:

HypothesisTestResult

pysht.mean.dempster_1samp(x, *, popmean=None)

Perform Dempster’s one-sample non-exact mean-vector test.

The test uses a Euclidean mean-square ratio and an approximate F law with trace-estimated effective degrees of freedom. It assumes multivariate normal observations but does not require an invertible covariance matrix.

Parameters:

• x (ArrayLike)

• popmean (ArrayLike | None)

Return type:

HypothesisTestResult

pysht.mean.dempster_2samp(x, y)

Perform Dempster’s equal-covariance two-sample non-exact test.

Unlike the legacy SHT routine, the positive mean-square ratio is calibrated with Dempster’s approximate F distribution, not a standard-normal tail.

Parameters:

• x (ArrayLike)

• y (ArrayLike)

Return type:

HypothesisTestResult

pysht.mean.yao_2samp(x, y)

Perform Yao’s approximate multivariate Behrens-Fisher test.

Both group covariance matrices may differ. The covariance of the mean difference must be positive definite.

Parameters:

• x (ArrayLike)

• y (ArrayLike)

Return type:

HypothesisTestResult

pysht.mean.johansen_2samp(x, y)

Perform Johansen’s Welch-James multivariate mean test.

Parameters:

• x (ArrayLike)

• y (ArrayLike)

Return type:

HypothesisTestResult

pysht.mean.nvm_2samp(x, y)

Perform the Nel-Van der Merwe multivariate Behrens-Fisher test.

Parameters:

• x (ArrayLike)

• y (ArrayLike)

Return type:

HypothesisTestResult

pysht.mean.ky_2samp(x, y)

Perform the Krishnamoorthy-Yu modified Nel-Van der Merwe test.

Parameters:

• x (ArrayLike)

• y (ArrayLike)

Return type:

HypothesisTestResult

Validation: classical mean, Dempster, and multivariate Behrens–Fisher.

High-dimensional trace and diagonal tests

pysht.mean.cq_2samp(x, y)

Perform the Chen–Qin dense high-dimensional two-sample mean test.

The two covariance matrices may differ. The Chen–Qin statistic is studentized using the translation-invariant unbiased trace estimators from Li–Chen (2012), Equations (2.1)–(2.2), rather than Chen–Qin’s different finite-sample leave-two-out estimate. At least four observations are required in each group. The asymptotic normal calibration requires the moment, factor, balance, and trace conditions in the validation ledger.

Parameters:

• x (ArrayLike)

• y (ArrayLike)

Return type:

HypothesisTestResult

pysht.mean.li_1samp(x, *, popmean=None)

Perform Li’s fixed-small-sample high-dimensional one-sample test.

The procedure treats the inner products from all unordered observation pairs as an asymptotically independent univariate sample as dimension diverges. Its Student calibration requires at least three observations and the factor-model and trace conditions in Li (2023).

Parameters:

• x (ArrayLike)

• popmean (ArrayLike | None)

Return type:

HypothesisTestResult

pysht.mean.li_2samp(x, y)

Perform Li’s fixed-small-sample high-dimensional two-sample test.

The defining Scheffe construction pairs the smaller sample with the first rows of the larger sample. Row order is consequently part of the realized statistic and must be unrelated to the measurements.

Parameters:

• x (ArrayLike)

• y (ArrayLike)

Return type:

HypothesisTestResult

pysht.mean.li_ksamp(*samples)

Perform Li’s fixed-small-sample high-dimensional ANOVA test.

The first smallest group in caller order is the paper’s reference population. Each other group is paired with that reference through the paper’s Scheffe construction, so row order is part of the realized statistic. Li’s theorem permits non-strictly ordered sample sizes. When several groups tie for the minimum, changing which tied group is supplied first can therefore change the realized statistic, while the theorem’s null calibration remains valid for every fixed reference choice.

Parameters:

samples (ArrayLike)

Return type:

HypothesisTestResult

Chen–Qin targets dense shifts without requiring equal covariances. The Li tests instead target genuinely fixed, small sample sizes while dimension diverges. Li’s two- and multi-sample Scheffé constructions pair the smaller sample with the first rows of larger samples, so row order must be arbitrary with respect to the measurements. The multi-sample method requires a uniquely smallest reference group; tied minima are rejected rather than selected by group input order.

pysht.mean.bs_1samp(x, *, popmean=None)

Perform the one-sample Bai-Saranadasa trace test.

The denominator uses the unbiased estimator of tr(Sigma @ Sigma); this corrects the missing n + 2 factor in legacy SHT’s one-sample code.

Parameters:

• x (ArrayLike)

• popmean (ArrayLike | None)

Return type:

HypothesisTestResult

pysht.mean.bs_2samp(x, y)

Perform the equal-covariance Bai-Saranadasa two-sample test.

Parameters:

• x (ArrayLike)

• y (ArrayLike)

Return type:

HypothesisTestResult

pysht.mean.sd_1samp(x, *, popmean=None)

Perform the Srivastava-Du one-sample diagonal mean test.

Parameters:

• x (ArrayLike)

• popmean (ArrayLike | None)

Return type:

HypothesisTestResult

pysht.mean.sd_2samp(x, y)

Perform the equal-covariance Srivastava-Du two-sample test.

Parameters:

• x (ArrayLike)

• y (ArrayLike)

Return type:

HypothesisTestResult

Validation: Chen–Qin and Li, Bai–Saranadasa, and Srivastava–Du.

Randomized and sparse tests

pysht.mean.ljw_2samp(x, y, *, calibration='asymptotic', n_resamples=999, rng=None)

Perform the Lopes-Jacob-Wainwright random-projection mean test.

A single Gaussian projection is drawn independently of the data, as in the defining procedure. Permutation calibration, when requested, holds that projection fixed for the observed and every permuted statistic. The pooled observations are canonicalized before seeded randomization, so row reordering or swapping sample labels cannot change a fixed-seed result.

The conditional F calibration assumes Gaussian sampling with a common covariance matrix. Unrestricted permutation calibration is exact under the stronger null that the pooled observations are exchangeable; equality of means alone is not sufficient when the distributions differ.

Parameters:

• x (ArrayLike)

• y (ArrayLike)

• calibration (Literal['asymptotic', 'monte-carlo'])

• n_resamples (int)

• rng (Generator | int | integer | None)

Return type:

HypothesisTestResult | ResamplingTestResult

calibration="asymptotic" uses the conditional F law for the single projection. calibration="monte-carlo" holds that same projection fixed for the observed statistic and every permutation. The F law assumes Gaussian samples with a common covariance. The permutation option requires pooled exchangeability; equal means alone do not make unrestricted relabeling exact. With an integer seed, canonical pooling makes the full Monte Carlo result invariant to row reordering and sample exchange.

pysht.mean.thulin_2samp(x, y, *, n_subspaces=100, n_resamples=999, rng=None)

Perform Thulin’s random-subspaces two-sample mean test.

The sampled feature subsets are drawn once and kept fixed throughout the permutation test. This both follows the conditional randomization design and removes the auxiliary-randomness bug in legacy SHT. The pooled rows and sample labels are canonicalized before seeded randomization.

Unrestricted permutation calibration is exact only when observations are exchangeable under the null. In particular, equality of means without a common distribution is not enough for exact finite-sample inference.

Parameters:

• x (ArrayLike)

• y (ArrayLike)

• n_subspaces (int)

• n_resamples (int)

• rng (Generator | int | integer | None)

Return type:

ResamplingTestResult

The subspace dimension is the paper’s floor((n_x + n_y - 2) / 2) and is reported as a diagnostic. The sampled subspaces are fixed across all permutations. The same pooled-exchangeability requirement applies. Canonical pooling makes fixed-integer-seed results invariant to row reordering and sample exchange.

pysht.mean.maximum_pairwise_bayes_factor_2samp(x, y, *, a0=0.0, b0=0.0, alpha=2.01, gamma=None)

Compute the Lee-You-Lin maximum pairwise Bayes factor for means.

The returned statistic is the maximum log Bayes factor. Component log Bayes factors are retained, and no frequentist p-value or automatic evidence threshold is manufactured. Residual sums of squares use the ordinary centered maximum-likelihood variances in Equation (4) of the primary paper. By default, gamma=max(n_x+n_y, p)**(-alpha) with the paper’s numerical choice alpha=2.01. A supplied gamma overrides that rate. Nonzero a0 or b0 selects pySHT’s optional shared inverse-gamma extension rather than the published Equation (4).

Parameters:

• x (ArrayLike)

• y (ArrayLike)

• a0 (float)

• b0 (float)

• alpha (float)

• gamma (float | None)

Return type:

BayesFactorTestResult

The primary statistic is the maximum component log Bayes factor. Positive values favor a coordinate-wise alternative over the null, but pySHT does not invent a universal evidence threshold. The published default is alpha=2.01 with gamma=max(n_x+n_y, p)**(-alpha); an explicit gamma overrides it. The paper’s Equation (4) corresponds to a0=b0=0. Nonzero a0 or b0 selects a clearly labeled pySHT inverse-gamma extension. The paper’s two-sample model uses a common covariance matrix.

Validation: Lopes–Jacob–Wainwright, Thulin, Lee–You–Lin.

Validation-blocked methods

SHT’s mean2.2014CLX (the Cai–Liu–Xia maximum mean test) has no public pySHT callable. Its practical-size Gumbel calibration failed the release gate, and neither estimated-precision branch has a passing advertised sparse-covariance scenario. The private implementation is retained only for scientific audit work; see the Cai–Liu–Xia validation ledger.

The Xue–Yao multiplier test is also withheld. Its formula implementation and batched fixed-multiplier oracle agree, but the requested default of 999 draws failed the project’s null-calibration gate at 0.01 and 0.05. There is no public xy_2samp callable; the exact pilot design and rejection counts are recorded in the Chen–Qin, Li, and Xue–Yao ledger.

Multi-group tests

pysht.mean.schott_ksamp(*samples)

Perform Schott’s high-dimensional one-way MANOVA mean test.

The error sum-of-products matrix uses (n_i - 1) S_i for every group. This is the defining MANOVA error matrix; legacy SHT instead used n_i S_i and thereby mismatched its claimed N - k error degrees of freedom.

Parameters:

samples (ArrayLike)

Return type:

HypothesisTestResult

pysht.mean.zx_ksamp(*samples, base_test='bai-saranadasa')

Perform the Zhang-Xu k-sample Behrens-Fisher mean test.

Scheffe’s transformation reduces the unequal-covariance k-sample problem to one sample. base_test="bai-saranadasa" applies the corrected Bai-Saranadasa calibration; base_test="hotelling" applies the exact Gaussian Hotelling calibration when the transformed dimension permits it.

The paper’s transformation uses the first n_min observations from larger groups. Their row order must therefore be arbitrary with respect to the measurements, as it is for an ordinary i.i.d. sample. If groups tie for the smallest size, their input order fixes the reference group and is part of the realized Scheffe transformation.

Parameters:

• samples (ArrayLike)

• base_test (Literal['bai-saranadasa', 'hotelling'])

Return type:

HypothesisTestResult

base_test="bai-saranadasa" is the high-dimensional default; base_test="hotelling" requires more transformed observations than transformed variables. The defining Scheffé construction uses the first n_min rows from larger groups. Consequently, changing their row pairing can change the finite-sample statistic even though its null law remains valid for i.i.d. rows. If groups tie for the smallest size, input order selects the reference group; changing the order of tied groups can also change the realized statistic.

pysht.mean.cph_ksamp(*samples, variance_estimator='original')

Perform the Cao-Park-He high-dimensional k-sample mean test.

variance_estimator selects the paper’s split-sample estimator or the Hu estimator. Both implementations use the current group in every loop; legacy SHT accidentally reused the last group in the split-sample branch.

Parameters:

• samples (ArrayLike)

• variance_estimator (Literal['original', 'hu'])

Return type:

HypothesisTestResult

variance_estimator="original" uses independent covariance estimates from the paper’s split; "hu" uses the corrected full-sample trace estimator. The split branch requires five rows per group and is intentionally sensitive to row order. Do not sort rows by the outcome before applying it.

Validation: Schott, Zhang–Xu, and Cao–Park–He.