Goodness of Fit Tests Based on Joint Densities of Multiple Sample Statistics

Published 2 Jul 2026 in stat.ME | (2607.02285v1)

Abstract: We propose goodness-of-fit tests based on simulated confidence sets for joint distributions of multiple sample statistics, focusing on absolutely continuous null distributions with known parameters. One class of tests uses hyperrectangular confidence sets for principal components of order statistics and related statistic vectors. Extending earlier work on horizontal and vertical confidence bands for cumulative distribution functions, these tests are compared with some classical, Zhang, and related graphical tests. Simulations show that the proposed procedures are competitive with, and often more powerful than, existing methods. We also study the geometry of principal-component-based statistics; under a normal null distribution, the first principal component corresponds to the sample mean, while the second is related to a linear analogue of variance. A second class of tests uses confidence sets of arbitrary shape constructed through highest density regions. Unlike earlier kernel-density-based approaches, we use a k-nearest-neighbor method for detecting highest density regions, which is better suited to higher-dimensional statistic vectors. We study tests based on order statistics, empirical distribution function values, moments, and combinations of classical goodness-of-fit statistics. The resulting procedures are powerful against a wide range of alternatives. We also outline a two-sample extension via permutation tests based on joint distributions of several statistics and compare moment-based versions with energy-distance permutation tests. Finally, we discuss transformations other than the probability integral transform, showing that mapping data to another target distribution, such as the standard normal, can be advantageous when powerful tests are available for that distribution.

Abstract PDF Upgrade to Chat

Authors (1)

Roman Guchenko

Summary

The paper develops a comprehensive GoF testing framework by estimating the joint distribution of multiple sample statistics through simulated confidence regions.
It employs PCA-based hyperrectangle confidence sets with variance weighting, enhancing sensitivity compared to classical univariate methods.
The study demonstrates robust performance using k-nearest neighbor HDR estimation and offers extensions for two-sample testing in high dimensions.

Goodness of Fit Tests Based on Joint Densities of Multiple Sample Statistics

The paper "Goodness of Fit Tests Based on Joint Densities of Multiple Sample Statistics" (2607.02285) develops a broad computational framework for constructing and evaluating goodness-of-fit (GoF) tests grounded in simulated confidence regions for the joint distribution of multiple sample statistics. Two major methodologies are developed: one based on high-dimensional hyperrectangular confidence sets, potentially in the space of principal components of statistics (notably order statistics), and another based on the direct estimation of highest density regions (HDRs) of these joint statistics, using $k$ -nearest neighbor (kNN) procedures adapted for high dimensionality.

Framework and Problem Formulation

The work considers the universal GoF setup: given an IID sample $x_1, ..., x_n \sim P$ , with $P$ unknown, and a fixed absolutely continuous null distribution $P_0$ , construct tests of $H_0: P=P_0$ . Rather than limit attention to classical univariate GoF test statistics or marginal statistics, the authors focus on the entire joint distribution $\mathcal{P}_{T_1,...,T_k}(P_0)$ of $k$ symmetric sample statistics $T_j$ , for $j=1,...,k$ , computed from random samples of size $n$ drawn from $x_1, ..., x_n \sim P$ 0. The hypothesis is rejected if the vector of statistics for the observed data is untypical with respect to this $x_1, ..., x_n \sim P$ 1-dimensional empirical null distribution.

A core computational insight is to approximate $x_1, ..., x_n \sim P$ 2 via a large sample-of-samples: simulate $x_1, ..., x_n \sim P$ 3 datasets of size $x_1, ..., x_n \sim P$ 4 from $x_1, ..., x_n \sim P$ 5, compute the $x_1, ..., x_n \sim P$ 6 statistics for each, and use the resulting empirical distribution as a surrogate for the true $x_1, ..., x_n \sim P$ 7-variate law (Figure 1).

Figure 1: Histograms for empirical distributions of order statistics for samples from standard normal distribution of size 10.

Construction of High-Dimensional Confidence Sets

Hyperrectangular Confidence Sets: Order Statistics and Principal Component Rotation

Historically, works such as Aldor-Noiman et al. (2013) and Sailynoja et al. (2022) computed marginal quantile bounds for order statistics or eCDF values and defined the acceptance region as the product (hyperrectangle) of these individual bounds, possibly with equal weights per coordinate. This paper generalizes the construction in two crucial directions:

Optimization of Hyperrectangle Orientation via PCA: Rather than default to the canonical coordinates (e.g., order statistics), the paper proposes rotating the data into the principal component (PC) basis of the simulated null sample-of-samples, and constructing hyperrectangular bounds along these axes. For the case of order statistics under normality, the first principal component corresponds to the sample mean, the second captures "linear analog of variance," and higher components reflect finer features (skewness, kurtosis, etc.) (Figure 2).
Figure 2: Histograms for empirical distributions of principal components of order statistics for samples from standard normal distribution of size 10.
Weighting of Margins According to Variance ("Importance"): The width of the hyperrectangle along each PC direction is set in proportion to the corresponding eigenvalue (i.e., the variability of the statistic under the null), making the test more sensitive to deviations detected by high-variance components.

Numerical studies demonstrate that, for normal and uniform nulls, the PC-based hyperrectangle test outperforms both unrotated hyperrectangles and classical GoF tests (e.g., KS, CvM, AD) against a breadth of alternatives, including shift, scale, skew, and heavy-tail distributions (Figures 3-7).

Figure 3: Pairwise contour plots for 2d projections of empirical joint density of order statistics for standard normal with pairwise correlations of order statistics.

Figure 4: Pairwise contour plots for 2d projections of empirical joint density of principal components of standard normal order statistics. Confidence rectangles are constructed using equal weights for principal components. Principal components are uncorrelated by construction, which we see in lower triangle.

Figure 5: Pairwise hexbin plots for 2d projections of empirical joint density of principal components of order statistics for standard normal samples of size 10. Confidence rectangles are constructed using equal weights for principal components.

Direct Estimation of Highest Density Regions

For alternatives where the convexity constraints or the marginal unimodality assumptions do not hold, the authors propose a nonparametric approach: directly estimate the highest density regions (HDRs) of the $x_1, ..., x_n \sim P$ 8-dimensional null statistic distribution. This is accomplished using the $x_1, ..., x_n \sim P$ 9-nearest neighbor density estimate at the observed statistic value; if the observed vector lies in a region with density below the ( $P$ 0)-quantile of the null’s density values, $P$ 1 is rejected. This method inherits robustness to the "curse of dimensionality" relative to kernel density methods and enables potent composite tests by combining disparate statistics (e.g., mean, variance, skewness, kurtosis) in higher dimensions.

Optimality and Power Analysis

Extensive simulation studies establish the following empirical findings:

PC-based hyperrectangle tests (with variance-based weights) are markedly more powerful than marginal (coordinate-aligned) hyperrectangles for normal and uniform nulls across location, scale, skewed, and heavy-tailed alternatives.
For arbitrary alternatives, PC-based tests are generally on par with the best likelihood (e.g., Zhang’s $P$ 2, $P$ 3) and metric-based (KS, AD, CvM) tests, but can dramatically outperform classical tests in situations where deviations impact higher-order components or the global geometry of the distribution.
Type I error is properly controlled at the nominal level, provided sufficiently large sample-of-samples and correct Monte Carlo implementation of empirical quantiles.
For nulls that are not normal or uniform, the probability integral transform followed by application of a powerful test for a convenient reference null (e.g., normal) is shown to be surprisingly effective. The invariance properties of order-statistics-based hyperrectangles under such transformations are formally demonstrated and illustrated in power tables and figures.
Principal components of order statistics for uniform laws correspond to sine and cosine-type eigenfunctions, in line with Karhunen-Loève expansions of the Brownian bridge—the geometric structure is substantiated both numerically and analytically (see Figure 2, Figure 4).
Figure 6: Order statistics from N(0,1) (red hexes) and N(1,1) (blue hexes), hyperrectangle bounds for N(0,1) (red rectangles).

Figure 7: Principal components of order statistics from N(0,1) (red hexes) and N(1,1) (blue hexes), hyperrectangle bounds for N(0,1) (red rectangles).

Theoretical and Computational Implications

This joint-density simulation framework represents a natural extension of "graphical" GoF methods (confidence bands for CDFs) and can, in principle, exploit any informative set of statistics jointly. The results substantiate that the geometry of the statistic space (as revealed by PCA or similar rotations) is critical: naive coordinate-wise bounds can admit implausible statistic combinations, lowering power, while PC-rotation "boxes" better match the actual high-density region of the null.

The methodology is computationally intensive, but the paper includes extensive R/C++ code as a reproducible backbone for further extensions—practically, the storage of large sample-of-samples and the explicit calculation of multivariate quantiles are feasible for moderate $P$ 4 and $P$ 5. For large-scale testing, algorithmic acceleration (e.g., through fast approximate nearest neighbors for HDRs) remains future work.

Extensions: Arbitrary Transforms, Two-Sample Testing, and Permutation-Invariant Procedures

The probability integral transform and its inverse enable recycling of optimal tests developed for a reference null into arbitrary nulls, leveraging the same joint-density simulation procedures. This flexibility is both a practical and theoretical advance, facilitating the use of distribution-tailored tests universally.

Additionally, the paper briefly sketches the extension of the framework to two-sample problems using permutation tests based on the joint distribution of multidimensional sample statistics, showing parity with energy-based distances in several empirical settings.

Implications and Future Directions

The construction of GoF tests based on joint density regions (hyperrectangular or HDR) of multiple sample statistics offers significant gains in sensitivity and interpretability, especially in small to moderate sample sizes and for high-dimensional tests involving composite features (mean, variance, higher moments, order statistics). The geometric perspective—aligning confidence regions with the actual data structure—represents an advance over previous univariate and marginal approaches. Furthermore, the empirical finding that powerful GoF tests for a canonical null distribution (e.g., normal) can be deployed against arbitrary nulls via probability integral transforms should stimulate future investigations into optimal universal GoF methodologies.

Potential future work involves:

Algorithmic scaling for large $P$ 6 and $P$ 7 via approximate HDR computation.
Extension to discrete and mixed-type data using suitable statistics and sampling schemes.
Characterization of minimax or local asymptotic optimality for various statistics and test constructions.
Systematic empirical evaluation against other modern GoF methods (e.g., kernel, graph, or deep learning-based tests) for complex and structured data.

Conclusion

The proposed methodology for GoF testing—via simulated confidence sets for the joint distribution of high-dimensional sample statistics, with rotation and weighting optimized by PCA, and robust HDR detection—substantially advances the flexibility and power of goodness-of-fit inference. The numerical results demonstrate clear quantitative advantages over coordinatewise and classical tests. The geometric and algorithmic perspectives developed provide a unifying foundation for high-dimensional, interpretable, and transform-invariant GoF testing, paving the way for robust universal hypothesis tests in modern data analysis.

Markdown Report Issue