Package 'HDElliptical'

Description

Evaluates the angular central Gaussian log-likelihood for nonzero row vectors. Row lengths need not equal one because they cancel from the quadratic contribution.

Usage

acg_loglik(x, shape, include_constant = FALSE)

Arguments

Value

A scalar log-likelihood.

References

Tyler, D. E. (1987). A distribution-free M-estimator of multivariate scatter. Annals of Statistics, 15, 234–251. doi:10.1214/aos/1176350263.

Examples

directions <- rbind(c(1, 0), c(0, 1), c(1, 1))
acg_loglik(directions, diag(2))

Description

Tests equality of two multivariate means under a common covariance matrix using the original Bai–Saranadasa statistic. This is deliberately a separate API from covariance-unrestricted Chen–Qin tests: its numerator is

$\xi = \|\bar X_1 - \bar X_2\|^2 - \frac{N}{n_1 n_2}\operatorname{tr}(S_p),$

where $S_p$ is the unbiased pooled covariance matrix and $N=n_1+n_2$.

Usage

bai_saranadasa_two_sample_test(x, y)

Arguments

Details

The traces of $S_p$ and $S_p^2$ are computed using the smaller of a primal $p \times p$ matrix and a dual $N \times N$ Gram matrix. Both groups are divided by the same global numerical scale before any moments are formed. The standardised statistic is computed on that safe scale, while reported means, traces, numerator, and variance are converted back according to their physical dimensions. Consequently a reported fourth-order component can be infinite when its value exceeds double range, or zero when it is below the representable range, without making the standardised statistic or p-value non-finite. The finite internal normalising variance is retained in diagnostics. The estimated null variance must be strictly positive; no ridge, absolute-value repair, or variance flooring is applied.

Value

An object of class c("hd_location_test", "htest"). Its components include the original numerator xi, pooled covariance traces, the bias-corrected estimate of $\operatorname{tr}(\Sigma^2)$, and the estimated null variance.

References

Bai, Z. and Saranadasa, H. (1996). Effect of high dimension: by an example of a two sample problem. Statistica Sinica, 6, 311–329.

Examples

set.seed(22)
x <- matrix(rnorm(24), nrow = 8, ncol = 3)
y <- matrix(rnorm(27, mean = 0.2), nrow = 9, ncol = 3)
bai_saranadasa_two_sample_test(x, y)

Description

Applies the real-valued BASIC approximate inverse map independently to the eigenvalues of the fitted SSCM. For a shape eigenvalue $\lambda$, the map used by Raninen and Ollila is

$\widetilde\delta(\lambda)=\frac12\int_0^1 \frac{\lambda}{1-t+t\lambda}(1-t)^{p/2-1}\,dt.$

The implementation evaluates the equivalent smooth integral after $1-t=v^2$ with paired Gauss–Legendre rules, brackets each inverse root, and uses bisection. It is independent of the official lookup table: no interpolation, spline extrapolation, eigenvalue floor, or out-of-range repair is used. Raw inverted eigenvalues are finally normalized to sum p.

Usage

basic_shape(
  x,
  center = c("spatial", "mean", "none"),
  tol = 1e-08,
  max_iter = 1000L,
  strict = TRUE,
  keep_signs = FALSE,
  quadrature_order = 64L,
  inversion_tol = 1e-10,
  integration_tol = 1e-10,
  inversion_max_iter = 100L,
  max_bracket_iter = 100L
)

Arguments

Details

Boundary SSCM eigenvalues equal to zero map exactly to zero and are reported. For p = 1, shape is identically one and no numerical inversion is needed.

Value

An hd_shape_estimator with the BASIC shape, SSCM eigensystem, raw and normalized inverse eigenvalues, mapped values, integration and inversion errors, brackets, iterations, and boundary diagnostics.

References

Raninen, E. and Ollila, E. (2022). Bias adjusted sign covariance matrix. IEEE Signal Processing Letters, 29, 339–343. doi:10.1109/LSP.2021.3134940.

Examples

x <- matrix(c(-2, 0, -1, 2, 0, -2, 1, 1, 2, -1, 3, 2),
            ncol = 2, byrow = TRUE)
basic_shape(x, quadrature_order = 32)

Description

Computes the RSSCM data weight, forms $\widehat\Lambda_{RSSCM}=\widehat\alpha pS_{sgn}+ (1-\widehat\alpha)I_p$, and applies the same self-contained BASIC inverse map as basic_shape() to the eigenvalues of $\widehat\Lambda_{RSSCM}/p$. This follows the BASICS construction while avoiding the official code's spline extrapolation. The RSSCM and BASIC numerical diagnostics are both retained.

Usage

basics_shape(
  x,
  center = c("spatial", "mean", "none"),
  tol = 1e-08,
  max_iter = 1000L,
  strict = TRUE,
  keep_signs = FALSE,
  quadrature_order = 64L,
  inversion_tol = 1e-10,
  integration_tol = 1e-10,
  inversion_max_iter = 100L,
  max_bracket_iter = 100L
)

Arguments

Value

An hd_shape_estimator containing the BASICS shape, RSSCM, its raw/projected shrinkage weight, inverse-map details, and full diagnostics.

References

Raninen, E. and Ollila, E. (2022). Bias adjusted sign covariance matrix. IEEE Signal Processing Letters, 29, 339–343. doi:10.1109/LSP.2021.3134940.

Examples

x <- matrix(c(-2, 0, -1, 2, 0, -2, 1, 1, 2, -1, 3, 2),
            ncol = 2, byrow = TRUE)
basics_shape(x, quadrature_order = 32)

Description

Forms the centered sample covariance with divisor $n$, as in Bickel and Levina (2008), and applies

$s_{ij}\,1\{|s_{ij}|>\lambda\}.$

When threshold = NULL, $\lambda=C\sqrt{\log(p)/n}$. The book draft first defines a divisor-$n-1$ covariance and then cites the primary divisor-$n$ formula. divisor makes this finite-sample distinction explicit; the default is the primary-paper convention. Thresholding need not preserve positive definiteness and no eigenvalue repair is performed.

Usage

bickel_levina_covariance_threshold(
  x,
  threshold = NULL,
  constant = 1,
  center = TRUE,
  divisor = c("n", "n-1"),
  threshold_diagonal = TRUE
)

Arguments

Value

An hd_covariance_estimator list whose estimate field is the thresholded covariance matrix.

References

Bickel, P. J. and Levina, E. (2008). Annals of Statistics, 36, 2577–2604.

Examples

set.seed(31)
x <- matrix(rnorm(40), nrow = 10, ncol = 4)
bickel_levina_covariance_threshold(x, threshold = 0.2)

Description

Implements the formula-complete Gaussian benchmark printed in Chapter 4, kept distinct from gaussian_alpha_combination_test(), which implements the primary Feng–Lan–Liu–Ma Bonferroni rule. With unrestricted OLS residual variances $\hat\sigma_i^2$, residual correlation $\hat R$, and intercept estimates $\hat\alpha$, this function uses

$T_{sum} = \frac{T\sum_i\hat\alpha_i^2/\hat\sigma_i^2-N} {\{2\mathrm{tr}(\hat R^2)\}^{1/2}}$

and the maximum squared OLS intercept t statistic. Their upper-tail normal and extreme-value p-values are combined by the equal-weight Cauchy rule.

Usage

book_gaussian_alpha_cauchy_test(returns, trace_R2, factors = NULL)

Arguments

Details

The manuscript denotes the trace term abstractly and does not prescribe a unique finite-sample estimator for this Gaussian benchmark. Consequently, trace_R2 is required rather than silently replacing it by a residual plug-in or a bias-corrected alternative.

Value

An upper-tail alpha-test object explicitly labelled as a book benchmark rather than a primary named procedure.

References

High-Dimensional Data Analysis for Elliptically Symmetric Distributions, Chapter 4, equations for the Gaussian alpha benchmark.

Examples

f <- cbind(seq(-1, 1, length.out = 12))
y <- outer(seq_len(12), 1:3, function(i, j) cos(i / 2 + j))
book_gaussian_alpha_cauchy_test(y, trace_R2 = 3, factors = f)

Description

The primary-paper pilot covariance and variability estimate are

$\hat\sigma_{ij}=n^{-1}\sum_k Z_{ki}Z_{kj},\qquad \hat\theta_{ij}=n^{-1}\sum_k(Z_{ki}Z_{kj}-\hat\sigma_{ij})^2,$

with entrywise threshold

$\lambda_{ij}=\delta \sqrt{\hat\theta_{ij}\log(p)/n}.$

Both quantities therefore use the same divisor-$n$ covariance by default. Selecting divisor = "n-1" is explicit and is reported in the diagnostics; it is not the primary finite-sample formula.

Usage

cai_liu_adaptive_covariance_threshold(
  x,
  delta = 2,
  rule = c("hard", "soft", "scad", "adaptive_lasso"),
  scad_a = 3.7,
  adaptive_eta = 1,
  center = TRUE,
  divisor = c("n", "n-1"),
  threshold_diagonal = TRUE
)

Arguments

Value

An hd_covariance_estimator object retaining $\hat\theta$ and every entrywise threshold.

References

Cai, T. and Liu, W. (2011). Journal of the American Statistical Association, 106, 672–684. arXiv:1102.2237.

Examples

set.seed(33)
x <- matrix(rnorm(40), nrow = 10, ncol = 4)
cai_liu_adaptive_covariance_threshold(x, delta = 2)

Description

Tests equality of two high-dimensional mean vectors using the precision-adjusted maximum statistic of Cai, Liu, and Xia (2014). The calibration is the type-I extreme-value limit

$F(g)=\exp\{-\pi^{-1/2}\exp(-g/2)\},$

for

$G=M-2\log(p)+\log\{\log(p)\}.$

Usage

cai_liu_xia_two_sample_test(
  x,
  y,
  precision = NULL,
  precision_source = c("adaptive", "oracle", "feasible"),
  threshold_delta = 2,
  eigen_floor = sqrt(.Machine$double.eps)
)

Arguments

Details

precision_source deliberately separates three statistically different paths. With "oracle", precision is a known population precision matrix and the denominator is its diagonal. With "feasible", precision is a user-supplied estimate and the denominator is formed from empirical within-group variances of the transformed observations, each with divisor $n_k$, as in equations (6)–(7) of the paper. With "adaptive", the function estimates a precision matrix by hard-thresholding the pooled covariance and then uses that same feasible denominator. It never uses diag(precision) to standardise an estimated-precision statistic.

For the adaptive backend, the pooled covariance and $\widehat\theta_{ij}$ both use divisor $N=n_1+n_2$, and

$\lambda_{ij}=\delta \sqrt{\widehat\theta_{ij}\log(p)/N}.$

If the thresholded covariance is not numerically positive definite, its eigenvalues are raised to eigen_floor times a data-scale reference before inversion. The floor, the unadjusted and adjusted extreme eigenvalues, and the number of modified eigenvalues are all returned in diagnostics. Both samples are centered relative to observed anchors before precision transformation, avoiding cancellation under a large common translation. The adaptive covariance and fourth-order threshold moments are also formed on an internal common scale and then mapped back to the original units.

Value

An object of class c("hd_location_test", "htest"). The reported statistic G has the extreme-value calibration, while raw.statistic contains M. Coordinatewise transformed scores, their denominators, the maximizing coordinate index and name, and the precision matrix used by the test are retained in components. diagnostics identifies the precision path, denominator construction, and every adaptive-threshold eigenvalue adjustment.

References

Cai, T. T., Liu, W., and Xia, Y. (2014). Two-sample test of high dimensional means under dependence. Journal of the Royal Statistical Society: Series B, 76, 349–372.

Cai, T. T. and Liu, W. (2011). Adaptive thresholding for sparse covariance matrix estimation. Journal of the American Statistical Association, 106, 672–684.

Examples

set.seed(29)
x <- matrix(rnorm(80), 20, 4)
y <- matrix(rnorm(96, 0.2), 24, 4)
cai_liu_xia_two_sample_test(x, y)

omega <- diag(4)
cai_liu_xia_two_sample_test(
  x, y, precision = omega, precision_source = "oracle"
)

Description

Tests the tail null that canonical correlations after null_rank are zero. For null_rank = k, the statistic uses

$\Lambda_k=\prod_{j=k+1}^{m}(1-\widehat\rho_j^2),\qquad -\{n-1-(p_x+p_y+1)/2\}\log\Lambda_k,$

with chi-squared degrees of freedom (p_x-k)(p_y-k). The special case k = 0 is the full independence test printed in the book.

Usage

cca_bartlett_test(object, null_rank = 0L)

Arguments

Value

An object of class htest with Wilks' Lambda and the Bartlett approximation diagnostics.

References

Anderson, T. W. (2003). An Introduction to Multivariate Statistical Analysis, 3rd ed. Wiley.

Examples

basis <- stats::poly(seq_len(12), degree = 4)
x <- basis[, 1:2, drop = FALSE]
y <- cbind(0.8 * basis[, 1] + 0.6 * basis[, 3],
           0.3 * basis[, 2] + sqrt(0.91) * basis[, 4])
colnames(x) <- c("x1", "x2")
colnames(y) <- c("y1", "y2")
fit <- classical_cca(x, y)
cca_bartlett_test(fit, null_rank = 0)

Description

Tests equality of two mean vectors with the covariance-unrestricted U-statistic of Chen and Qin (2010). Its numerator is

$T_{CQ}=\|\bar X-\bar Y\|^2- \operatorname{tr}(S_X)/n_1-\operatorname{tr}(S_Y)/n_2.$

The variance estimator is the original paper's leave-out estimator of the two within-sample covariance traces and their cross trace. This is distinct from the common-covariance Bai–Saranadasa test.

Usage

chen_qin_two_sample_test(x, y)

Arguments

Details

The numerator is exactly invariant to a common translation. The original finite-sample leave-out variance estimator, however, is not translation invariant; this known property is retained rather than silently replacing the published estimator. Both samples are internally divided by one common global scale. The returned components include both the original-unit quantities and their finite internally scaled counterparts. An original-unit fourth-order quantity can legitimately overflow to Inf or underflow to zero when the input itself spans the limits of double precision; the reported statistic is reconstructed from the scaled quantities and remains well defined. No ridge, absolute-value repair, or variance flooring is applied.

Value

An object of class c("hd_location_test", "htest"). The components field includes the U-statistic numerator, leave-out trace estimates A1, A2, and A12, the estimated null variance, and their .scaled counterparts used for stable calibration.

References

Chen, S. X. and Qin, Y.-L. (2010). A two-sample test for high-dimensional data with applications to gene-set testing. Annals of Statistics, 38, 808–835.

Examples

set.seed(24)
x <- matrix(rnorm(21), nrow = 7, ncol = 3)
y <- matrix(rnorm(24, mean = 0.2), nrow = 8, ncol = 3)
chen_qin_two_sample_test(x, y)

Description

Implements the two rank procedures for which the primary paper supplies direct, scalable statistics and complete extreme-value calibrations: Spearman's rho and Kendall's tau. The paper does not establish rank-based sum or adaptive Fisher tests; those constructions in the book draft are not used here.

Usage

chen_song_feng_rank_white_noise_test(
  x,
  lag = 1L,
  measure = c("spearman", "kendall"),
  keep_lag = FALSE
)

Arguments

Value

An htest object. Exact ties are rejected because the published null variance and distribution-free calibration assume continuous margins.

References

Chen, D., Song, F. and Feng, L. Rank Based Tests for High Dimensional White Noise. Statistica Sinica 35, 1323–1347. doi:10.5705/ss.202022.0382.

Examples

t <- seq_len(18)
x <- cbind(sin(t), cos(t / 2), sin(t / 3 + 0.2))
chen_song_feng_rank_white_noise_test(x, lag = 2, measure = "spearman")

Description

Constructs the location-invariant unbiased estimators

$T_{1n}=Y_{1n}-Y_{3n},\qquad T_{2n}=Y_{2n}-2Y_{4n}+Y_{5n}$

of $\operatorname{tr}(\Sigma)$ and $\operatorname{tr}(\Sigma^2)$. For sphericity it uses $U_n=pT_{2n}/T_{1n}^2-1$; for identity it uses $V_n=T_{2n}/p-2T_{1n}/p+1$. In either case the primary null calibration is $nU_n/2$ or $nV_n/2$, respectively, against an upper standard-normal tail.

Usage

chen_zhang_zhong_covariance_test(x, null = c("sphericity", "identity"))

Arguments

Value

An hd_covariance_test object retaining all five U-statistic building blocks.

References

Chen, S. X., Zhang, L.-X. and Zhong, P.-S. (2010). Journal of the American Statistical Association, 105, 810–819. doi:10.1198/jasa.2010.tm09560.

Examples

set.seed(42)
x <- matrix(rnorm(48), nrow = 12, ncol = 4)
chen_zhang_zhong_covariance_test(x, null = "sphericity")

Description

Tests equality of two population spatial-sign covariance matrices (SSCMs), equivalently proportionality of the two scatter matrices under elliptical symmetry. With separately estimated spatial medians and centered signs, the feasible statistic is

$T=p(A+B-2C),$

where $A$ and $B$ average squared inner products over ordered within-group pairs and $C$ averages them over all cross pairs.

Usage

cheng_sscm_equality_test(
  x,
  y,
  alpha = 0.05,
  tol = 1e-08,
  max_iter = 500L,
  zero_tol = 0,
  strict = TRUE
)

Arguments

Details

The implementation uses the complete finite-sample bias expression in Cheng et al. (2019), Equation (3.3). For each group it estimates the inverse radial moments through order three, evaluates the $n^{-2}$ and $n^{-3}$ terms separately, and subtracts $p\widehat\delta$. Under the null, Remark 3.1 estimates

$q=p^{-1}\mathrm{tr}(\Lambda^2) =\frac{p}{n_1+n_2}\{n_1(A-\delta_1)+n_2(B-\delta_2)\}.$

Consequently the standard-error denominator and its square are

$\widehat v_0= \frac{2(n_1^{-1}+n_2^{-1})}{p+2}(p\widehat q),\qquad \widehat v_0^2=\{\widehat v_0\}^2.$

The reported statistic is $Z=(T-p\widehat\delta)/\widehat v_0$ and the paper's rejection rule is the upper standard-normal tail.

A non-positive trace estimate, a zero fitted radius, or an unconverged spatial median invalidates calibration. No absolute value, variance floor, radius perturbation, or replacement center is used. With strict = FALSE the raw components are returned, but statistic, p.value, and the rejection decision are NA and diagnostics$calibrated is FALSE. The plug-in bias calibration has the additional regime stated in Remark 3.2 of the primary paper: under bounded eigenvalues it requires $p=o(n_l^{3/2})$ for both groups. This is stronger than merely using the theorem's $p=O(n_l^2)$ expansion condition.

Value

An object of class c("hd_proportionality_test", "htest"). components contains all ordered-pair components, bias terms, scaled inverse radial moments, the trace estimate, the standard-error denominator and its squared variance, fitted centers, signs, and SSCMs.

References

Cheng, G., Liu, B., Peng, L., Zhang, B., and Zheng, S. (2019). Testing the equality of two high-dimensional spatial sign covariance matrices. Scandinavian Journal of Statistics, 46, 257–271. doi:10.1111/sjos.12350.

Examples

x <- matrix(c(-2, 0, 1, 2, 1, -1, 0, 2, -1, -2, 2, 0), 6, 2)
y <- matrix(c(-1, 1, 2, -2, 0, 3, 1, -1, 2, 0, -2, 1), 6, 2)
cheng_sscm_equality_test(x, y)

Description

Fits the CHIME iteration of Cai, Ma, and Zhang for a two-component Gaussian mixture with a common covariance. At stage $t$, the sparse discriminant vector is the certified solution of

$\min_\beta\;\tfrac12\beta^T\widehat\Sigma_t\beta- \beta^T(\widehat\mu_{1t}-\widehat\mu_{2t})+ \lambda_t\|\beta\|_1.$

The covariance is never inverted. Coordinate descent uses the actual diagonal of the quadratic metric and must satisfy the exact lasso KKT equations before an iterate is accepted.

Usage

chime_clustering(
  x,
  lambda,
  initial,
  K = 2L,
  max_iter = length(lambda) - 1L,
  lasso_tol = 1e-08,
  lasso_max_iter = 10000L,
  support_tol = 0,
  psd_tol = sqrt(.Machine$double.eps)
)

Arguments

Details

This function deliberately differs from the authors' numerical script in three ways required for a reusable unsupervised method: it applies no hidden covariance ridge, it never uses true labels to choose a penalty, and it requires the complete penalty path (including the initial beta fit) from the caller. The method is only defined here for K = 2, matching the stated model and the available primary implementation.

Value

A chime_fit object containing labels, responsibilities, all parameter stages, the explicit lambda path, covariance matrices, and per-stage objective, convergence, and KKT certificates. Score ties use the printed weak inequality and are assigned to component 1.

References

Cai, T. T., Ma, J., and Zhang, L. (2019). CHIME: Clustering of high-dimensional Gaussian mixtures with EM algorithm and its optimality. Annals of Statistics, 47, 1234–1267. doi:10.1214/18-AOS1711.

Examples

x <- rbind(
  c(-2.2, -1.0), c(-1.8, -1.2), c(-2.0, -0.8), c(-1.7, -1.1),
  c( 1.8,  1.0), c( 2.2,  1.1), c( 2.0,  0.9), c( 1.7,  1.2)
)
initial <- list(
  omega = 0.5, mu1 = c(-1.5, -0.8), mu2 = c(1.5, 0.8),
  covariance = diag(2)
)
chime_clustering(x, lambda = c(0.2, 0.15, 0.1), initial = initial,
                 max_iter = 2)

Description

Forms all three covariance blocks with one explicit divisor, constructs the symmetric inverse square roots of the two within-block covariance matrices, and applies an SVD to the whitened cross-covariance matrix. Both within-block matrices must be strictly positive definite and full rank at rank_tol. No ridge, pseudoinverse, or eigenvalue floor is used.

Usage

classical_cca(
  x,
  y,
  components = NULL,
  center = c("mean", "none"),
  covariance_divisor = c("n-1", "n"),
  rank_tol = sqrt(.Machine$double.eps)
)

Arguments

Value

A classical_cca_fit object with canonical coefficients, scores, structure loadings, correlations, covariance blocks, and whitened singular-vector projectors.

References

Hotelling, H. (1936). Relations between two sets of variates. Biometrika, 28, 321–377.

Examples

basis <- stats::poly(seq_len(8), degree = 4)
x <- basis[, 1:2, drop = FALSE]
y <- cbind(0.8 * basis[, 1] + 0.6 * basis[, 3],
           0.3 * basis[, 2] + sqrt(0.91) * basis[, 4])
colnames(x) <- c("x1", "x2")
colnames(y) <- c("y1", "y2")
fit <- classical_cca(x, y, components = 1)
fit$canonical.correlations

Description

Computes the univariate CUSUM

$C_n(k)=n^{-1/2}\{S_k-kS_n/n\}$

and rejects for a large trimmed maximum of $\{(k/n)(1-k/n)\}^{-\gamma}|C_n(k)|/\widehat\sigma$. For gamma = 0, calibration = "asymptotic" uses the classical two-sided Brownian-bridge (Kolmogorov) limit. calibration = "grid-gaussian" simulates the method's finite scan-grid Gaussian bridge; this is null calibration, not replication of a paper simulation study.

Usage

classical_cusum_test(
  x,
  gamma = 0,
  trim = 1L,
  sigma = NULL,
  variance = c("difference", "sample"),
  calibration = c("asymptotic", "grid-gaussian"),
  calibration_draws = 4999L,
  alpha = 0.05,
  seed = NULL,
  strict = TRUE
)

Arguments

Details

If sigma is not supplied, the scale is either the sample standard deviation or Rice's adjacent-difference estimate. This software choice is returned explicitly because the book's classical display assumes a population $\sigma$ and does not prescribe its estimator.

Value

An htest-compatible object with the complete CUSUM path.

References

Page (1954); Csorgo and Horvath (1997).

Examples

x <- c(-1, 0, 1, -1, 0, 3, 4, 3)
classical_cusum_test(x)

Description

Estimates the common covariance with the pooled unbiased estimator and inserts it into the Gaussian LDA log-likelihood ratio. The covariance must be strictly positive definite. No generalized inverse or ridge is used.

Usage

classical_lda_classifier(x, y, prior = "equal", strict = TRUE)

Arguments

Value

An hd_classifier_fit.

References

Anderson, T. W. (2003). An Introduction to Multivariate Statistical Analysis, 3rd ed. Wiley.

Examples

x <- rbind(c(2, 0), c(1, 1), c(1, -1),
           c(-2, 0), c(-1, 1), c(-1, -1))
classical_lda_classifier(x, rep(c("A", "B"), each = 3))

Description

Computes PCA from the centered second-moment matrix with a user-visible divisor. The book's displayed estimator uses covariance_divisor = "n"; "n-1" gives the usual unbiased sample covariance. No variance floor or rank repair is applied. components = NULL returns all numerically usable positive-variance components.

Usage

classical_pca(
  x,
  components = NULL,
  center = c("mean", "none"),
  covariance_divisor = c("n-1", "n"),
  eigen_tol = sqrt(.Machine$double.eps)
)

Arguments

Details

Eigenvector signs are anchored by making the first coordinate attaining the largest absolute loading positive. Individual eigenvectors inside a repeated-eigenvalue block are not identified; the returned projector is the invariant object for such a block.

Value

A classical_pca_fit and hd_pca_fit object containing loadings, scores, covariance operator, eigenvalues, projector, reconstruction, and complete centering/divisor diagnostics.

References

Anderson, T. W. (2003). An Introduction to Multivariate Statistical Analysis, 3rd ed. Wiley.

Examples

x <- rbind(c(2, 0), c(1, 1), c(-1, -1), c(-2, 0))
fit <- classical_pca(x, components = 1, covariance_divisor = "n")
fit$loadings

Description

Uses classwise unbiased covariance matrices with divisors n1 - 1 and n2 - 1, then forms the canonical Gaussian QDA log-likelihood ratio. Both covariance matrices must be strictly positive definite. Ordinary positive determinants are used; no pseudo-determinant, absolute value, ridge, or generalized inverse is substituted.

Usage

classical_qda_classifier(x, y, prior = "equal", strict = TRUE)

Arguments

Value

An hd_classifier_fit.

References

Anderson, T. W. (2003). An Introduction to Multivariate Statistical Analysis, 3rd ed. Wiley.

Examples

x <- rbind(c(2, 0), c(1, 1), c(1, -1), c(2, 2),
           c(-2, 0), c(-1, 2), c(-1, -2), c(-2, -1))
classical_qda_classifier(x, rep(c("A", "B"), each = 4))

Description

Implements the fixed-dimensional spatial-sign test, the one-sample spatial signed-rank test, and the pooled two-sample spatial-rank test. All three use the zero-direction convention $U(0)=0$ and a Cholesky solve; a singular studentizing matrix is reported rather than replaced by a generalized inverse.

Usage

spatial_sign_test(x, mu = NULL, tol = sqrt(.Machine$double.eps), zero_tol = 0)

spatial_signed_rank_test(
  x,
  mu = NULL,
  tol = sqrt(.Machine$double.eps),
  zero_tol = 0
)

spatial_rank_test(x, y, tol = sqrt(.Machine$double.eps), zero_tol = 0)

Arguments

Details

For centered observations $Z_i=X_i-\mu_0$, the sign statistic uses

$\bar U=n^{-1}\sum_i U(Z_i),\quad B_U=n^{-1}\sum_i U(Z_i)U(Z_i)^T,\quad Q_{sign}=n\bar U^T B_U^{-1}\bar U.$

The signed ranks are

$R_i=n^{-1}\sum_j U(Z_i+Z_j).$

Since the first-order projection of their average has covariance $4B_R$, where $B_R=n^{-1}\sum_i R_iR_i^T$, its statistic is

$Q_{SR}=\frac{n}{4}\bar R^T B_R^{-1}\bar R.$

The factor 1/4 is essential.

For two samples, every observation is ranked against the pooled sample:

$R(Y_i)=N^{-1}\sum_j U(Y_i-Y_j).$

If $C=(N-1)^{-1}\sum_i R(Y_i)R(Y_i)^T$, the sample covariance of the pooled ranks (whose average is exactly zero), the statistic is

$Q_{2SR}=\frac{n_1n_2}{N}(\bar R_1-\bar R_2)^T C^{-1}(\bar R_1-\bar R_2).$

These calibrations are asymptotic chi-squared laws for fixed dimension, not high-dimensional approximations. Euclidean spatial signs make the tests invariant to translations, common positive rescaling, and orthogonal transformations, but not to arbitrary nonspherical affine transformations.

Value

An object of classes hd_location_test and htest. statistic is the chi-squared statistic, parameter is its asymptotic degrees of freedom, variance is the studentizing second-moment matrix (or ⁠4 B_R⁠ for the signed-rank test), and components contains the directional mean and unscaled moment matrices. Numerical rank, reciprocal condition number, solver, zero counts, and calibration type are in diagnostics.

References

Mottonen, J. and Oja, H. (1995). Multivariate spatial sign and rank methods. Journal of Nonparametric Statistics, 5, 201–213.

Oja, H. (2010). Multivariate Nonparametric Methods with R. Springer.

Examples

set.seed(12)
x <- matrix(rnorm(80), 20, 4)
spatial_sign_test(x)
spatial_signed_rank_test(x)

y <- matrix(rnorm(96, 0.25), 24, 4)
spatial_rank_test(x, y)

Description

Solves the Cai–Liu–Luo CLIME program column by column,

$\min_b\|b\|_1\quad\text{subject to}\quad \|S_n b-e_j\|_\infty\leq\lambda,$

using a primal–dual algorithm. Every raw column must pass primal feasibility, dual feasibility, lasso stationarity, and relative duality-gap certificates. The final symmetric estimator retains the entry of smaller absolute value from each transposed pair; exact ties retain the row-column entry before mirroring, making the convention deterministic.

Usage

clime_precision(
  x,
  lambda,
  center = TRUE,
  divisor = c("n", "n-1"),
  solver_tol = 1e-07,
  solver_max_iter = 100000L,
  strict = TRUE
)

Arguments

Details

Primary CLIME does not guarantee that this post-LP symmetrization remains feasible, nor that the final matrix is positive definite. Both facts are reported and neither is repaired. Consequently valid certifies the raw column programs, not an invented SPD condition. SCIO and scaled-lasso precision estimation remain review-only because their additional programs are not specified by the short book display.

Value

A gaussian_precision_fit list. estimate is the primary smaller-absolute-value symmetrization of the certified raw columns.

References

Cai, T. T., Liu, W. and Luo, X. (2011). A constrained $\ell_1$ minimization approach to sparse precision matrix estimation. Journal of the American Statistical Association, 106, 594–607. doi:10.1198/jasa.2011.tm10155.

Examples

x <- rbind(c(-2, 0), c(-1, -1), c(1, 1), c(2, 0))
clime_precision(x, lambda = 0)

Description

Implements the two-sample Composite $T^2$ test of Feng, Zou, Wang, and Zhu (2017). Observations are rows and variables are columns. Despite the wording at ch2_location.tex:778--784 in the current book draft, this is not a Behrens–Fisher test: the paper assumes that the two groups have the same unknown covariance matrix. Nor does the method combine existing test p-values or estimate correlations among component tests. Instead, it builds a block-diagonal approximation to the common pooled covariance and combines the resulting block Hotelling quadratic forms.

Usage

composite_t2_two_sample_test(x, y, block_size = 2L, selection = "paper_greedy")

Arguments

Details

For a fixed leave-out covariance, the paper's practical block rule starts with the pair having the largest absolute sample correlation and repeatedly adds the variable with the largest sum of absolute correlations to the current block. Variables already assigned to a block are removed and the procedure is repeated. selection = "paper_greedy" implements that rule, with lexicographic tie-breaking. The final block contains all remaining variables when fewer than block_size remain. The paper's notation and simulations use equal-size blocks; this deterministic remainder convention extends the displayed rule to dimensions not divisible by the block size.

Let $\widehat{\Sigma}_{O_K,i_1,i_2,j_1,j_2}^{-1}$ be the block-diagonal inverse obtained after deleting observations $i_1,i_2$ from group 1 and $j_1,j_2$ from group 2, recomputing the unbiased pooled covariance and its correlation matrix, and rebuilding the blocks. The published statistic is

$Q_n={1\over n_1n_2(n_1-1)(n_2-1)} \sum_{i_1\ne i_2}\sum_{j_1\ne j_2} (X_{1i_1}-X_{2j_1})^T \widehat{\Sigma}_{O_K,i_1,i_2,j_1,j_2}^{-1} (X_{1i_2}-X_{2j_2}).$

The upper-tail normal calibration is

$Z={Q_n\over {2(n_1^{-1}+n_2^{-1})^2 \widehat{\mathrm{tr}(\Lambda_K^2)}}^{1/2}}.$

Following Theorem 2 and the paragraph immediately after it, the feasible trace estimate is the paper's one-sample leave-four-out estimator computed from group 1 only:

$\widehat{\mathrm{tr}(\Lambda_K^2)}= {1\over 2P_{n_1}^4}\sum^* (X_{i_1}-X_{i_2})^T\widehat{\Sigma}_{O_K,\setminus4}^{-1} (X_{i_3}-X_{i_4}) (X_{i_1}-X_{i_4})^T\widehat{\Sigma}_{O_K,\setminus4}^{-1} (X_{i_3}-X_{i_2}).$

Thus $Q_n$ and its population scaling are symmetric in the two samples, but the paper's feasible finite-sample standardisation is label-asymmetric. Exchanging the samples can change the reported $Z$ and p-value.

The first group must satisfy nrow(x) >= block_size + 5: after leaving out four rows, a block of size block_size then has enough residual degrees of freedom to be nonsingular. The second group must contain at least three rows for the pooled 2+2 leave-out covariance. These count conditions are necessary, not sufficient. Every marginal variance used for selection and every selected full or leave-out covariance block must be finite and strictly positive definite. Failure produces an error. No ridge, generalized inverse, absolute-value repair, or variance floor is applied.

The C++ kernel uses sufficient statistics and sums unordered pairs and quadruples as exact symmetry reductions of the displayed ordered formulas. Both samples are internally translated by a common anchor and divided by a common positive scale in each column. This is algebraically neutral and protects the scale-invariant method under very large or small units.

Value

An object of class c("hd_location_test", "htest"). The statistic is the standardised $Z$ and the p-value is its upper standard-normal tail probability. components retains $Q_n$, the group-1 leave-four-out trace estimate, every calibration coefficient and ordered denominator, combination counts, and a full-sample diagnostic partition. diagnostics records the common-covariance assumption, dynamic block selection, first-sample calibration, conditioning summaries, and the absence of numerical repair.

References

Feng, L., Zou, C., Wang, Z., and Zhu, L. (2017). Composite T-squared test for high-dimensional data. Statistica Sinica, 27, 1419–1436. doi:10.5705/ss.202015.0199.

Examples

set.seed(2717)
x <- matrix(rnorm(21), 7, 3)
y <- matrix(rnorm(15, 0.2), 5, 3)
composite_t2_two_sample_test(x, y)

Description

This helper performs only the algebra specified by the conditional-alpha papers. Rows are time observations. If center_alpha = TRUE, its first block is the column-centered supplied basis; otherwise it is the original basis. The remaining blocks are, in factor order, the rowwise products $f_{jt}B(t/T)$. It does not choose spline order, knots, basis dimension, or a BIC rule.

Usage

conditional_alpha_sieve_design(
  basis,
  factors = NULL,
  center_alpha = TRUE,
  alpha_contrast = NULL
)

Arguments

Details

A normalized B-spline basis usually contains the constant function, so all of its centered columns are linearly dependent. The papers write an ordinary inverse despite this identity. This function never silently drops a column: pass an explicit full-column-rank alpha_contrast, or supply a reduced basis. The contrast changes only the coordinates, not the spanned centered-alpha space, when it has the intended range.

Value

A full-column-rank design matrix carrying its basis/factor metadata.

References

Ma, S., Lan, W., Su, L. and Tsai, C.-L. (2020). Testing alphas in conditional time-varying factor models with high-dimensional assets. Journal of Business & Economic Statistics, 38, 214–227. doi:10.1080/07350015.2018.1482758.

Examples

tt <- seq(0, 1, length.out = 12)
basis <- cbind(tt, tt^2)
factors <- cbind(market = sin(seq_len(12)))
conditional_alpha_sieve_design(basis, factors)

Description

Fits, independently for every asset, the no-intercept null regression on a supplied full-rank sieve design. The residualized intercept $h=M_Z1_T$ is retained because it enters every feasible calibration. The function uses one common response scale and independent design-column scales for numerical stability; both leave the fitted residual subspace and every reported test unchanged. It never uses a generalized inverse or ridge.

Usage

conditional_alpha_sieve_fit(returns, design)

Arguments

Value

A reusable conditional_alpha_sieve_fit object containing the restricted residuals, h, design metadata, and projection diagnostics.

References

Ma, S., Lan, W., Su, L. and Tsai, C.-L. (2020). Testing alphas in conditional time-varying factor models with high-dimensional assets. Journal of Business & Economic Statistics, 38, 214–227. doi:10.1080/07350015.2018.1482758.

Examples

tt <- seq(0, 1, length.out = 14)
z <- conditional_alpha_sieve_design(cbind(tt, tt^2))
y <- cbind(sin(1:14), cos(1:14), sin(1:14 / 2))
conditional_alpha_sieve_fit(y, z)

Description

Evaluates the Chapter 4 supplied-estimate benchmark

$W=T\hat\delta^T\hat\Omega_\delta^{-1}\hat\delta,$

with a chi-squared reference distribution having length(delta) degrees of freedom. This function does not estimate spline or kernel nuisances.

Usage

conditional_factor_wald_test(delta, covariance, sample_size)

Arguments

Value

A strict Cholesky-based upper-tail Wald htest object.

References

Li, D. and Yang, L. (2011). Nonparametric tests of conditional factor models. Ang, A. and Kristensen, D. (2012). Testing conditional factor models.

Examples

conditional_factor_wald_test(c(0.1, -0.2), diag(c(1, 2)), 40)

Description

Fits the Mai–Zou–Yuan DSDA lasso with the exact response coding $-n/n_1$ for class 1 and $+n/n_2$ for class 2 and objective

$n^{-1}\sum_i(y_i-a-x_i^\top\beta)^2+ \lambda\|\beta\|_1.$

The regression coefficient points toward class 2. The package reverses it only when constructing its documented score, so positive always means class 1. Following the primary equal-prior rule, classification uses the midpoint intercept rather than the penalized regression intercept. The coordinate-descent solution must pass both KKT and Fenchel dual-gap checks.

Usage

dsda_classifier(
  x,
  y,
  lambda,
  prior = c(0.5, 0.5),
  solver_tol = 1e-07,
  solver_max_iter = 100000L,
  strict = TRUE,
  tie = c("class1", "class2")
)

Arguments

Value

An object of class hd_classifier_fit.

References

Mai, Q., Zou, H. and Yuan, M. (2012). A direct approach to sparse discriminant analysis in ultra-high dimensions. Biometrika, 99, 29–42. doi:10.1093/biomet/asr066.

Examples

x <- rbind(
  c(2, 1), c(1, 2), c(2, 2), c(3, 1),
  c(-2, -1), c(-1, -2), c(-2, -2), c(-3, -1)
)
y <- factor(rep(c("first", "second"), each = 4))
dsda_classifier(x, y, lambda = 0.2)

Description

Implements the convex off-diagonal-lasso branch of Liu, Wang and Zhao's EC2 estimator. The primary method first forms the empirical correlation matrix $R$, solves

$\min_{C:\,\operatorname{diag}(C)=1,\,\lambda_{\min}(C)\geq\tau} \frac12\|R-C\|_F^2+\lambda\sum_{i\ne j}|c_{ij}|,$

and restores marginal empirical standard deviations. This differs from the book's shortened direct-covariance display. The default covariance divisor is $n$; selecting "n-1" is explicit and changes only the marginal covariance scale, not the sample correlation.

Usage

ec2_covariance(
  x,
  lambda,
  tau,
  penalty = "l1",
  center = TRUE,
  divisor = c("n", "n-1"),
  rho = 1,
  solver_tol = 1e-07,
  solver_max_iter = 50000L,
  strict = TRUE
)

Arguments

Details

The ISP/ADMM solver uses the paper's soft-threshold and spectral-projection updates. A fit is valid only when the equality, fixed-point, off-diagonal subgradient, spectral-dual, complementarity, diagonal, and eigenvalue certificates all pass. The method-defining tau projection is not a numerical repair. Adaptive EC2 and MC+ EC2 require additional weight or shape choices and remain review-only; this function never guesses them.

Value

An ec2_covariance_fit list. estimate is the covariance matrix only for a certified fit; correlation is its EC2 correlation estimate.

References

Liu, H., Wang, L. and Zhao, T. (2014). Sparse covariance matrix estimation with eigenvalue constraints. Journal of Computational and Graphical Statistics, 23, 439–459. doi:10.1080/10618600.2013.782818.

Examples

x <- rbind(c(-2, -1), c(-1, 0), c(1, 0), c(2, 1))
ec2_covariance(x, lambda = 0.2, tau = 0.1)

Description

Computes the trace-$p$ spatial-sign pilot centered at the sample spatial median and applies the eigenvalue-ratio (ER) or growth-ratio (GR) selector of Xu, Ma, Wang and Feng. For ordered eigenvalues $\lambda_j$,

$\widehat m_{ER}=\arg\max_{1\le j\le M} \lambda_j/\lambda_{j+1}.$

The GR criterion uses $V_j=\sum_{l=j+1}^{\min(n,p)-1}\lambda_l$ exactly as in the primary paper. max_factors is deliberately supplied: the paper treats $M$ as a predetermined upper bound and its numerical choice is not a universal tuning rule.

Usage

elliptical_factor_number(
  x,
  max_factors,
  method = c("er", "gr"),
  tol = 1e-08,
  max_iter = 500L,
  zero_tol = 0,
  strict = TRUE
)

Arguments

Value

An elliptical_factor_number object containing the selected count, all pilot eigenvalues, criterion values, and spatial-median diagnostics.

References

Xu, X., Ma, H., Wang, H. and Feng, L. (2025). High dimensional matrix estimation through elliptical factor models. arXiv:2512.19325. https://arxiv.org/abs/2512.19325.

Examples

x <- rbind(c(-3, -2, 0), c(-2, -1, 1), c(-1, 0, -1),
           c(1, 0, 1), c(2, 1, -1), c(3, 2, 0))
elliptical_factor_number(x, max_factors = 1)

Description

Applies the primary CLIME or GLASSO problem to the raw idiosyncratic complement of a spatial_sign_poet() or poet_tme() fit, and reconstructs the full inverse scatter through the Woodbury identity,

$V=V_u-V_u\Gamma(\Lambda^{-1}+\Gamma^T V_u\Gamma)^{-1} \Gamma^T V_u.$

CLIME must pass per-column primal, dual, stationarity, gap, and final symmetrized-feasibility certificates. GLASSO must remain positive definite and satisfy its full elementwise-$\ell_1$ KKT equation. The tuning parameter is explicit because the theoretical constant multiplying the rate is not an observable default.

Usage

elliptical_factor_precision(
  fit,
  lambda,
  method = c("sclime", "sglasso"),
  solver_tol = 1e-07,
  solver_max_iter = 100000L,
  initial_step = 1,
  max_backtracking = 100L,
  strict = TRUE
)

Arguments

Value

An elliptical_factor_precision_fit containing the idiosyncratic and full precision estimates and complete solver certificates.

References

Xu, X., Ma, H., Wang, H. and Feng, L. (2025). arXiv:2512.19325.

Examples

x <- rbind(c(-3, -2), c(-2, -1), c(-1, 1),
           c(1, -1), c(2, 1), c(3, 2))
f <- spatial_sign_poet(x, factors = 0, threshold = 0.2)
elliptical_factor_precision(f, lambda = 0.5, method = "sglasso")

Description

Computes the exact generator-aware score

$\log\pi_1-\tfrac12\log|\Lambda_1|+\log g_1(d_1) -\log\pi_2+\tfrac12\log|\Lambda_2|-\log g_2(d_2),$

where $d_k=(z-\mu_k)'\Lambda_k^{-1}(z-\mu_k)$. The log-generator functions must be vectorized. They may return -Inf for a zero density, but not NA, NaN, or Inf.

Usage

elliptical_oracle_classifier(
  location1,
  location2,
  shape1,
  shape2,
  log_generator1,
  log_generator2 = log_generator1,
  prior = "equal",
  levels = c("class1", "class2"),
  feature_names = NULL
)

Arguments

Details

A common generator and common shape do not generally reduce this rule to a midpoint linear rule under unequal priors. That shortcut is exact for equal priors, and for the exponential radial generator underlying Gaussian LDA.

Value

A generator-aware hd_classifier_fit.

References

Fang, K.-T. and Anderson, T. W. (1990). Statistical Inference in Elliptically Contoured and Related Distributions. Allerton Press. Wakaki, H. (1994). Discriminant analysis under elliptical populations. Hiroshima Mathematical Journal, 24, 257–298.

Examples

normal_log_generator <- function(d) -d / 2
fit <- elliptical_oracle_classifier(
  c(1, 0), c(-1, 0), diag(2), diag(2), normal_log_generator
)
predict(fit, rbind(c(1, 0), c(-1, 0)))

Description

Computes the feasible ERHT statistic over a fixed deterministic ridge grid and combines its upper-tail normal p-values with the analytic Cauchy rule of Feng, Zhou, and Wang (2026). Equal weights are used by default. The default grid $\{0.1,0.2,\ldots,1\}$ is the grid used in the paper's numerical implementation; this function does not reproduce any simulation design from that paper.

Usage

elliptical_regularized_hotelling_cauchy_test(
  x,
  mu = NULL,
  rho = seq(0.1, 1, by = 0.1),
  weights = NULL,
  alpha = 0.05,
  tol = 1e-08,
  max_iter = 1000L,
  strict = TRUE,
  keep_companion = FALSE
)

Arguments

Details

For positive weights $\varpi_k$ summing to one,

$T_{CC}=\sum_k\varpi_k\tan[\pi\{1/2-p_k\}],\qquad p_{CC}=1/2-\pi^{-1}\arctan(T_{CC}).$

The computation uses signed-log cotangents and an angle representation to avoid overflow and cancellation in extreme tails. The article proves dependence-robust small-tail validity for this analytic p-value; except in special dependence structures it does not claim exact fixed-level calibration at an ordinary level such as 0.05.

Value

An object of class c("hd_location_test", "htest"). Its p-value is the analytic Cauchy combination. The marginal table contains every raw ERHT statistic, feasible center and variance, standardized statistic, p-value, log-p-value, and fixed-ridge decision.

References

Feng, L., Zhou, L., and Wang, X. (2026). Elliptical regularized Hotelling testing for high dimensional data. arXiv:2606.25942. doi:10.48550/arXiv.2606.25942.

Examples

set.seed(2)
x <- matrix(rnorm(160), 40, 4)
elliptical_regularized_hotelling_cauchy_test(
  x, rho = c(0.2, 0.5, 1)
)

Description

Tests a high-dimensional location vector under elliptical symmetry and pervasive dependence using the fixed-ridge ERHT statistic of Feng, Zhou, and Wang (2026). Observations are rows and variables are columns.

Usage

elliptical_regularized_hotelling_test(
  x,
  mu = NULL,
  rho = 0.5,
  alpha = 0.05,
  tol = 1e-08,
  max_iter = 1000L,
  strict = TRUE,
  keep_companion = FALSE
)

Arguments

Details

Let $\widehat\theta$ be the sample spatial median, $\widehat Y_i=\sqrt p\,U(X_i-\widehat\theta)$, and

$\widehat R_n=n^{-1}\sum_i\widehat Y_i\widehat Y_i^{\mathsf T}.$

For $\rho>0$, the raw statistic is

$T_n(\rho)=n(\widehat\theta-\theta_0)^{\mathsf T} (\widehat R_n+\rho I)^{-1}(\widehat\theta-\theta_0).$

The function implements the paper's direct feasible companion-matrix centering $n\widehat\mu_n(\rho)$ and variance $n\widehat\sigma_{D,n}^2(\rho)$, returning the upper-tail normal calibration of

$Z_n(\rho)=\{T_n(\rho)-n\widehat\mu_n(\rho)\}/ \{n\widehat\sigma_{D,n}^2(\rho)\}^{1/2}.$

The C++ kernel uses one economy SVD of the sign matrix. It evaluates the ridge quadratic in row and null spaces and switches between the companion matrix $A$ and its complement $I-A$ for the feasible calibration. Thus it never forms a persistent $p\times p$ inverse or subtracts nearly equal Woodbury, diagonal-weight, or companion terms. A zero fitted residual makes the required inverse distance undefined and is reported as an error; no observation is omitted or perturbed and no additional ridge, pseudoinverse, absolute-value repair, or variance floor is used.

The paper's numerical section mentions a separate Bartlett center correction with $a=8n/p$, but neither the article nor its public arXiv source specifies the formula that maps $a$ into the feasible center. This function therefore implements the fully stated equations only and deliberately does not guess that simulation-program modification.

Value

An object of class c("hd_location_test", "htest"). The reported statistic is $Z_n(\rho)$; raw.statistic is $T_n(\rho)$. Components include all feasible functionals, raw and internally scaled values, inverse distances, spatial-median diagnostics, and companion eigenvalues. Raw dimensional fields can underflow to zero or overflow to infinity when their physical units exceed the double-precision range; the separately normalized fields in components$scaled remain the numerical-audit representation used to compute the reported statistic.

References

Feng, L., Zhou, L., and Wang, X. (2026). Elliptical regularized Hotelling testing for high dimensional data. arXiv:2606.25942. doi:10.48550/arXiv.2606.25942.

Examples

set.seed(1)
x <- matrix(rnorm(120), 30, 4)
elliptical_regularized_hotelling_test(x, rho = 0.5)

Description

Implements the single-change path or the discretized adjacent-triple scan of Song, Wen, and Feng (2026). For adjacent segments with sizes n1, n2, the raw statistic is

$V_\rho^{raw}=N\widehat\Delta^T (\widehat R+\rho I)^{-1}\widehat\Delta, \qquad N=n_1n_2/(n_1+n_2).$

The common pool is the full sample. Its centered spatial signs are columns of $Y$; with $A=Y^TQY/n$, the primary studentization is

$\widehat\kappa=\sum_i\beta_i^2A_{ii},\qquad \widehat\sigma^2=2n\sum_{i\ne j}\beta_i^2\beta_j^2A_{ij}^2,$

$Z_\rho=(V_\rho^{raw}-n\widehat\kappa)/ \{n\widehat\sigma^2\}^{1/2}.$

The ordered off-diagonal variance is accumulated directly. It is never obtained by flooring a cancellation-prone difference.

Usage

erht_change_point_test(
  x,
  ridge,
  scan = c("single", "multiple"),
  epsilon = 0.1,
  calibration = c("gaussian-supremum", "time-permutation", "none"),
  calibration_draws = 4999L,
  alpha = 0.05,
  seed = NULL,
  tol = 1e-08,
  max_iter = 1000L,
  zero_tol = 0,
  cauchy_weights = NULL,
  keep_calibration = FALSE,
  strict = TRUE
)

Arguments

Details

ridge contains the actual positive $\rho$ values in the primary formula. The public ERHTCP reproduction code instead accepts empirical ratios rho / (p/n); this function does not silently make that conversion. Users who want that convention can supply ridge_ratio * p / n explicitly, and all resolved ridge values are returned.

calibration = "gaussian-supremum" simulates the primary marginal Gaussian-process limit on the exact candidate grid. "time-permutation" is the paper's practical exchangeability calibration and preserves each multivariate row. It is invalid under unaddressed serial dependence. With multiple ridge values, the returned Cauchy transform combines the marginal p-values, but calibration.exact is FALSE: the paper explicitly distinguishes this analytic rule from exact calibration by a joint-limit quantile involving the unknown cross-ridge correlation r_E.

The primary inverse-distance average is exactly mean(sqrt(p) / distance), with a coincident observation assigned weight zero. The optional (p - 1) / p factor used by the public reproduction repository is not part of the displayed primary statistic and is not used.

Value

An htest object with every local raw/center/variance/Z component.

References

Song, Wen and Feng (2026), arXiv:2607.22162.

Examples

x <- rbind(
  c(-2, 0), c(2, 0), c(0, -2), c(0, 2),
  c(-1, -1), c(1, 1), c(-1, 1), c(1, -1),
  c(-2, 1), c(2, -1), c(-1, 2), c(1, -2)
)
erht_change_point_test(
  x, ridge = 0.5, epsilon = 0.25,
  calibration_draws = 99, seed = 1
)

Description

Implements the primary WBS–ERHT algorithm of Song, Wen, and Feng (2026). On an interval I, every (or, when triple_step > 1, an explicitly thinned) adjacent triple compares two segments of length at least ceiling(epsilon * length(I)); the common scatter pool is the whole interval. Scores are maximized over the supplied actual ridge grid.

Usage

erht_wbs(
  x,
  ridge,
  threshold,
  min_interval,
  refinement_radius,
  deletion_radius,
  intervals = NULL,
  M = NULL,
  epsilon = 0.1,
  triple_step = 1L,
  max_changes = Inf,
  seed = NULL,
  tol = 1e-08,
  max_iter = 1000L,
  zero_tol = 0,
  strict = TRUE
)

Arguments

Details

The recursion uses the paper's narrowest-over-threshold rule: among sampled intervals contained in the current segment, plus the current segment itself, retain those with score strictly greater than threshold; choose the narrowest, break a length tie by the larger score, refine within the requested radius, delete the requested neighborhood, and recurse.

The primary theory does not choose the ridge grid, WBS threshold, number of intervals, minimum length, refinement radius, or deletion radius. Hence this API requires those quantities (or explicit intervals) rather than inventing automatic tuning constants. Random WBS intervals are generated independently of the observations from uniformly sampled unordered endpoint pairs. Scanning/WBS is intrinsic to the estimator and is not replication of a paper simulation study.

triple_step = 1 is the literal exhaustive adjacent-triple collection. Larger values are an explicit computational approximation and are recorded. Any invalid spatial median or local variance invalidates the fit; failed triples are never silently removed.

Value

A change-point fit with sorted estimates and selection history.

References

Song, Wen and Feng (2026), arXiv:2607.22162, Algorithm 1.

Examples

x <- rbind(
  c(-2, 0), c(2, 0), c(0, -2), c(0, 2),
  c(-1, -1), c(1, 1), c(-1, 1), c(1, -1),
  c(-2, 1), c(2, -1), c(-1, 2), c(1, -2)
)
erht_wbs(
  x, ridge = 0.5, threshold = 10,
  min_interval = 8, refinement_radius = 3, deletion_radius = 1,
  intervals = matrix(c(1, 12), ncol = 2), epsilon = 0.25
)