Papers

In medical research, patients often have multiple interdependent outcomes, such as posttraumatic stress disorder (PTSD), depression, and pain among trauma survivors. Most existing research uses multinomial regression to analyze these interdependent outcomes separately, which ignores correlations between concurrent conditions. This omission may lead to loss of information and reduced predictive accuracy. Accounting for correlations between multiple categorical outcomes requires a high-dimensional parameter space, making model estimation challenging. In this paper, we propose a multivariate multinomial logit model that captures outcome correlations and uses the ANOVA decomposition of the parameter space to reduce the number of parameters. The ANOVA decomposition enables explicit conditional model formulations, which allow for a computationally much simpler composite likelihood for model estimation. We develop an efficient Minorization-Maximization (MM) algorithm to maximize the composite likelihood, which also incorporates variable selection via a bridge penalty. Simulation studies are conducted to evaluate our method, demonstrating its accuracy in parameter estimation and variable selection. We further illustrate our method using data from the AURORA study.

We analyze approximate message passing (AMP) with an independent Gaussian initialization for noiseless phase retrieval in the proportional asymptotic regime. A random initialization has overlap of order $d^{-1/2}$ with the signal, and AMP requires a growing number of iterations to attain non-vanishing overlap. Thus, its precise behavior cannot be characterized by classical fixed-time state evolution. We prove a Gaussian decomposition of the AMP trajectory and control its error over the horizons required for recovery. The resulting analysis shows that random initialization attains the weak-recovery threshold $\delta_{\rm weak}=1/2$. For $\delta\in(\delta_{\rm weak},\delta_{\rm str})$, where $\delta_{\rm str}\approx1.13$, the signal strength follows state evolution and approaches its stable finite fixed point uniformly for \(n^{1/3}/\operatorname{polylog}(n)\) iterations. For $\delta>\delta_{\rm str}$, AMP reaches any prescribed fixed recovery accuracy within $O_{\delta,\varepsilon}(\log n)$ iterations. The majority of our analysis applies more generally to generalized AMP for single-index models.

Datasets in many settings naturally partition into clusters arising from sub-populations, batch effects, or aggregation across multiple sources. A common response to such heterogeneity is to ensemble learners trained on each cluster rather than fit a single model to the pooled data. Prior work motivating such approaches has typically considered settings in which both the covariate distribution and the conditional outcome model differ across clusters; the role of cluster-aware partitioning and ensembling based solely on the covariate distribution remains to be explored. We address this case for ridge-regularized least-squares regression under a linear outcome model and consider all ridge penalty values $\lambda \geq 0$, including the special case of the ridgeless predictor at $\lambda = 0$. By considering both fixed-effects and random-effects models, we argue that under random effects, an optimally tuned pooled ridge predictor always outperforms ensembles of individually optimally tuned predictors. For fixed effects, we derive a general formula for the pooled and ensembled predictors to characterize the role of both regression coefficients as well as the predictor distribution shifts. Together, these results generalize prior risk analyses of bagging and random-partition estimation using ridge and ridgeless regression predictors from the i.i.d. setting to encompass covariate shift and heterogeneity-aware partition structure.

Target trial emulation structures observational causal analyses around the protocol of an ideal randomized trial. By aligning eligibility, treatment assignment, time zero, and follow-up, it can reduce avoidable biases, but implementation still requires coordinated decisions about data construction, inverse probability weighting, diagnostics, outcome models, standardization, competing risks, and uncertainty estimation. This article provides a self-contained methodological guide and practical tutorial for TTE, an R package for target trial emulations with longitudinal observational data. We describe target trial protocols, intention-to-treat and per-protocol estimands, identification assumptions, baseline and person-period data structures, temporal ordering for longitudinal weights, stabilized treatment and censoring weights, weight truncation, balance and effective-sample-size diagnostics, weighted pooled discrete-time survival models, model-based standardization, competing-risk analysis, weighted Kaplan-Meier and Aalen-Johansen estimation, and cluster bootstrap at the original-individual level. Two fully synthetic examples illustrate end-to-end workflows: sodium-glucose cotransporter 2 inhibitor versus dipeptidyl peptidase-4 inhibitor initiation with all-cause death, and sequentially nested angiotensin receptor blocker versus calcium channel blocker trials with heart-failure hospitalization and competing death. The examples show how to obtain and diagnose estimates in R and interpret relative hazards, absolute risks, cumulative incidence, and differences between intention-to-treat and per-protocol effects.

We propose a robust deep mixture model based on a pathway-wise shared scale-mixture construction. Layer-specific component indicators are independently distributed according to their corresponding mixing proportions and jointly define a complete pathway through the latent hierarchy. Conditional on the selected pathway, a single gamma-distributed latent precision variable is shared across the deepest latent distribution, every intermediate latent transition, and the observation model. Integrating out this shared precision yields an exact multivariate Student-$t$ distribution for each complete pathway, allowing robustness to propagate coherently throughout the entire latent hierarchy rather than being introduced separately within individual latent layers. Model parameters are estimated using a stochastic expectation–maximisation algorithm. Complete-pathway responsibilities are evaluated analytically, whereas the shared latent precision variables and latent Gaussian variables are generated from their conditional distributions before updating the model parameters. The pathway-specific degrees-of-freedom parameters are estimated by one-dimensional numerical optimisation. Simulation studies demonstrate accurate recovery of the pathway-specific degrees-of-freedom parameters together with consistently improved clustering performance relative to the deep Gaussian mixture model under heavy-tailed and contaminated settings. Real-data applications further illustrate the ability of the proposed model to identify heterogeneous latent structures while reducing the influence of atypical observations. The proposed framework retains the hierarchical representation and parsimonious parameter-sharing structure of the deep Gaussian mixture model while providing coherent pathway-wise robustness.

We propose method-of-moments estimators for the eigenvalues of variance component covariance matrices in multivariate mixed effects models. Assuming a parametric form for the eigenvalue distribution, we focus on the high-dimensional regime where the number of predictors is large and comparable to the number of realizations of each random effect. In this setting, we show that the empirical moments of sum-of-squares matrices (e.g., MANOVA estimators of the covariance matrices) can be closely approximated by deterministic functions of the underlying parameters. This relationship enables the construction of consistent and asymptotically normal estimators via moment matching. Our approach is motivated by applications in quantitative genetics, where estimating genetic covariance components across multiple phenotypic traits is of central interest. We implement our method in a new python package mlmm-mom, and demonstrate how our method adapts to several common experimental designs in this domain.

Machine-learning benchmarks often pair a label that aggregates a long temporal horizon with input observed through one or a few short windows. Their apparent performance ceiling may therefore be an acquisition-protocol ceiling rather than a model-capacity ceiling. We study labels of the form $\Theta_{g,T}=T^{-1}\int_0^T g\{Z(t)\}\,\mathrm{d}t$ when the latent Gaussian process contains both a stable individual trait and a correlated within-individual state. An exact protocol-conditioned Bayes-risk identity provides a common tool. First, we decompose label variance into an $O(1)$ trait component and an $O(T^{-1})$ state component, explaining why a snapshot can retain cross-sectional predictability while poorly tracking within-person change. Second, we derive task-dependent effective temporal spans: mean labels depend on the ordinary correlation time, whereas occupation-time labels depend on an entire spectrum of higher-order correlation times. Third, state-driven occupation-label variance is maximal when the stable trait lies at the threshold; window efficiency decays much more slowly away from that boundary. Under an equal segment budget, exact risks and Monte Carlo experiments show that repeated segments at one time rapidly saturate, whereas temporally dispersed observations continue to increase state explainability. The trait ceiling uses quantities available from ordinary test-retest data; only the state ceiling requires short-lag temporal calibration. The results distinguish architectural limits from protocol limits and show that the label, rather than duration or segment count alone, defines the relevant timescale.

In recent years, numerous open-source software libraries have been developed for computing sets of features from univariate time series. The type and number of features vary across these feature sets, which have been constructed with varying disciplinary perspectives on quantifying structure in time-series data. To date, the relative strengths and weaknesses of these feature sets on time-series classification problems remains largely unexplored. Here we aimed to understand the relative performance of six open-source feature sets and three baseline feature sets (based on distributional and/or basic spectral structure) across 124 univariate time-series classification problems using a normalization-based approach to problem-level benchmarking that better indexes the relative strengths and weaknesses of different algorithms compared to prior rank-based approaches. Despite their dramatic differences in size, composition, and computation time, we found that feature sets performed relatively similarly overall (85.3% of pairwise comparisons resulted in ties), with the largest feature set, tsfresh, exhibiting the strongest overall performance (29.03% wins across all pairwise comparisons against other feature sets). We also highlighted specific problems on which the specific composition of a given feature set gave it a substantial performance advantage or disadvantage, and problems where simple baselines comprised of Fourier coefficients and quantiles were sufficient to achieve strong performance. Our results demonstrate the need to consider problem-level performance when benchmarking time-series feature sets, and highlight the importance of feature make-up in driving relative classification performance.

Drifting objectives compare a target and model distribution through a vector field observed noisily at finitely many locations. We ask what distributional conclusion such a frozen measurement system warrants. For integrable antisymmetric interactions and absolutely continuous laws in a declared finite density basis, the unnormalized sampled numerator satisfies $\operatorname{vec}(V_X)=Mc$, where $c$ is an antisymmetric mismatch and $M$ is probe-dependent. This identity yields an a posteriori total-variation (TV) upper confidence bound accounting for held-out field noise, estimated-operator error, and externally validated $L^1$ residual radii around normalized density approximants in the span; a nonpositive observability margin returns the trivial TV bound and abstains. The audit recomputes this numerator from held-out samples; a normalized drift statistic requires a separate joint numerator–denominator analysis. For Gaussian-RBF interactions, a global envelope supports distribution-free and empirical-Bernstein radii without truncation, with companion bounds for the Laplace similarity in the original drifting objective. We characterize random-probe observability by a population Gram matrix, identify rank and symmetry degeneracies, and prove large-bandwidth collapse toward mean matching. Synthetic studies exercise Gaussian and Laplace numerators, separately prespecified bounded-vector and variance-adaptive radii, Monte Carlo-calibrated operators, nonzero residual radii around normalized finite-basis approximants, outward-rounded observability bounds, and designed abstention. A joint basis-size/dimension stress path extends evaluation through $m=8$. The result is a conditional diagnostic for a finite density class, or for normalized finite-basis density approximants with external residual radii, not a universal guarantee from small training drift.

Background. Aircraft-wildlife collisions are a documented aviation hazard. The FAA National Wildlife Strike Database is the main public source of U.S. reports, but raw extracts are large and rarely released as versioned aggregates. Methods. We convert the public FAA Access export (STRIKE_REPORTS) to SQLite, define a report as damaging when DAMAGE_LEVEL is one of M, M?, S, D, and fold rows into open aggregates (year, airport, species, aircraft, animal family) via hierarchical SPECIES_ID prefixes. We analyze all civil reports and a commercial subset (AC_MASS 3-5). Rates are per reported strike, not per flight hour, so every rate is conditional on a report having been filed. Trend statistics use complete incident years 1990-2025; cross-sectional facets cover the full seed, incident years 1990-2026. Results. Numbers are locked to the 2026-07-04 freeze. Reported strikes in 2025 reach a series high (24,458 total; 901 damaging). Species that dominate strike counts are not those with the highest damage rates: perching songbirds account for 41.1% of identified strikes but damage aircraft in 1.5% of reports (95% CI 1.4-1.6), whereas deer account for 0.7% but damage aircraft in 81.8% of reports (95% CI 79.8-83.7). Airports ranked by total strikes (led by Denver Intl) diverge from rankings by damaging strikes (led by Sacramento Intl); across the 452 airports with at least 50 reports, volume and conditional damage rate are negatively rank-correlated (Spearman rho = -0.44, p < 0.001); large stations plausibly capture minor events more completely, so this is not a hazard ranking. The share of reports labeled damaging declines over decades, consistent with improved reporting of minor events, while absolute damaging counts peak in 2025. Conclusions. Versioned open aggregates and a documented damage-rate analysis enable independent reuse; released under CC BY 4.0.

We study the problem of identifying the dominant arm in multi-armed bandits, where the objective is to find the action with the highest probability of exceeding the realized rewards of all other actions. Conventional mean-based and pairwise comparison-based algorithms often fail to identify the arm with the highest realized reward. To address this challenge, we introduce a novel dominant arm criterion and an efficient estimator with theoretical guarantees. Our approach relies on two key technical innovations: (i) a dominance score criterion that an arm beats the locally dominant over the partitioned reward space and (ii) a joint mixing and recycling mechanism coupled with a doubly robust estimator that guarantees simultaneous convergence of the empirical distribution functions for all arms. These key innovations pave a way to efficient computation of global arm dominance. Our proposed elimination algorithm identifies the best dominant arm with nearly optimal rate of sample complexity. Numerical experiments demonstrate that our algorithm consistently achieves exact recovery of the true dominant arm, outperforming existing baselines.

Methods for estimating causal effects in longitudinal, quasi-experimental settings are widely used in economics, public health, political science, and other fields. However, studies evaluating effects of health policies often rely on limited sample sizes, both in terms of study units and time periods analyzed. Synthetic difference-in-differences (SDiD) and the augmented synthetic control method (ASCM) are recently developed methods for evaluating the effects of policy interventions. Although SDiD and ASCM generally rely on weaker assumptions than both DiD and SCM, they lack theoretical performance guarantees with respect to bias or coverage in relevant settings with small sample sizes or short panel lengths. To evaluate the performance of SDiD and ASCM in realistic, small-sample settings, we employ a calibrated simulation strategy that allows the injection of a known treatment effect into existing data in a setting of interest. Drawing on findings from these empirical investigations, we offer practical guidance for researchers and policymakers on when SDiD and ASCM are likely to yield reliable estimates and inferences under realistic scenarios.

Covariance matrices play a central role across diverse domains such as finance, machine learning, neuroscience, and genetics, where they are used for tasks including dimensionality reduction, connectivity inference, and risk estimation. Many of these applications require testing whether a covariance matrix follows a specific structure. While several R packages provide partial coverage of such tests, Python, an important language in machine learning, lacks a comprehensive, well-tested implementation. To address this gap, we introduce scikit-covtest, a Python package implementing a variety of hypothesis tests for covariance matrices spanning four categories: identity, sphericity, proportionality, and two-sample equality. The package provides a consistent SciPy-style API, extensive documentation, and supporting functionality for multiple testing correction, synthetic data generation, and diagnostic evaluation. scikit-covtest is open source, available through PyPI, and lowers the barrier to applying modern covariance-testing methods in scientific applications.

Sequential subset search – forward selection with floating backtracking and its descendants – remains the quality reference in feature selection, but every member of the family sweeps the full pool of remaining candidate features at each step, which excludes it from very-high-dimensional problems; there, only individual-feature ranking remains practical, and it models feature interplay weakly or not at all. We introduce a budgeted sampled step operator pair that replaces the full sweeps by a fixed number of candidate evaluations per step. Candidates are drawn by temperature-controlled softmax sampling from dependency-aware per-feature statistics learned online from every criterion evaluation the search performs, guarded by a uniform exploration floor; per-step cost becomes independent of dimensionality. Substituting the operators turns any sequential method into its stochastic counterpart, defining the Stochastic Sequential Search (SSS) family; we study the stochastic counterpart of floating search, sSFFS. On 500-dimensional madelon, sSFFS retains at least 97% of the full-SFFS criterion value at every subset size at about a quarter of its evaluations, while uniform sampling at the same budget collapses on madelon's synergistic features. On 5,000-dimensional gisette, far beyond full-SFFS reach, sSFFS exceeds the saturated criterion level of DAF and BIF ranking at matched budgets; holdout validation shows that at 500 training samples the binding constraint beyond the sequential frontier becomes the criterion, not the search. On 10,105-dimensional reuters, under a trustworthy multinomial filter criterion, sSFFS dominates BIF and DAF on the search objective and on holdout accuracy at every subset size, in about two minutes of single-core evaluation work. A verified standalone implementation accompanies the paper.

In this paper, we investigate a computational class of random inverse problems that incorporates model uncertainties through random variable parameters nonlinearly in the forward model as well as additive observational uncertainty. Random inverse problems with nonlinear parameter dependencies may arise in engineering, geophysics, image processing or uncertainty quantification. We study a perspective on structural ambiguities due to the non-injectivity of the forward model together with probabilistic ambiguities by assigning mixture model densities with separate components to the parameters, which leads to an in general complex forward model, observation model and posterior. As a result, the mixture-model parameters in the forward model can be viewed as simultaneously describing aspects of both the nonlinear ambiguity and the uncertainty of the inverse problem. The presented solution algorithms are based on Bayesian inversion and Monte Carlo computation for three observation scenarios leading to posterior densities for the input for given output samples or an observed output density. We apply the derived algorithms to analytic 1D and 2D quadratic models, to an epidemiological inverse parameter estimation and to a heat equation for inverse source localization. We numerically demonstrate in which observation scenarios the proposed algorithm can resolve probabilistic ambiguities in the solution and in which they cannot and show the interaction patterns between these two types of ambiguities. The results suggest that the proposed perspective is useful in making residual structural ambiguities visible in highly ambiguous inverse problems, including cases with a finite and an infinite number of solutions.

Normally the statistical mechanics of learning treats constraints on weight distributions as restrictions that shrink the space of possible solutions. Therefore, it reduces model capacity. In this paper we would like to take a contrary approach, which, however, is based on the earlier work on distribution-constrained perceptrons. Rather than treating a prescribed weight distribution as a mere restriction, we propose that it defines the intrinsic geometry upon which learning naturally unfolds. We formulate both deep neural networks and variational quantum circuits as gradient flows on a product of Wasserstein manifolds – one classical Wasserstein space for each layer and one quantum Wasserstein space for the circuit parameters. Within this geometry, the capacity reduction, which was previously associated with distributional constraints, appears as the metric structure of the constraint manifold itself. We develop a hierarchical mean-field description for deep networks, extend the framework to the quantum setting using the quantum Wasserstein distance of order 1, and introduce two such practical algorithms, Hierarchical DisCo-SGD and Quantum DisCo, that follow approximate geodesics on the manifold of the product itself. Experiments on teacher-student problems, standard image classification tasks, and small variational quantum classifiers show that respecting these distributional geometries improves generalization, stabilizes training, and reduces the severity of barren plateaus compared with unconstrained and purely norm-based baselines. This approach firstly reframes structural constraints as geometric priors and suggests a route for incorporating biological, spectral, or hardware-derived distributional information into both learning systems, viz., classical and quantum learning.

Sequential trial emulation (STE) estimates the effect of a sustained treatment by stacking nested emulated trials with inverse-probability weighting. Additive-hazards STE estimators of the marginal risk difference recommend the nonparametric bootstrap without evaluating its coverage. Using a correctly-specifiable mechanism (up to a small, disclosed residual), we compare analytic and bootstrap standard errors for both estimands. Exploiting the closed-form linearity of the additive-hazards estimating equation, we derive the influence functions and prove the default row-level robust variance inconsistent for the marginal risk-difference curve, and show the same failure empirically for the constant hazard difference: the row-level variance omits a within-subject cross-trial covariance that is positive under a sign condition we verify across our mechanisms, and is anticonservative at every horizon except the first - only there, where the covariance is zero, is the row-level standard error unimpaired. The subject-clustered variance and multiplier bootstrap are consistent for the fixed-weight linearisation and support simultaneous confidence bands, whose measured coverage is 0.88. Across our simulations the model-based and row-level robust intervals are anticonservative and worsen with sample size, coverage falling to 0.71 at n=5000; clustering leaves the constant-hazard-difference coverage near 0.86 at n=5000, and the multiplier bootstrap leaves the risk-difference-curve coverage near 0.90 (0.86 at the longest horizon). The STE constant hazard difference is a design-weighted summary of a time-varying effect, dependent on the trial structure. We illustrate on the Stanford heart transplant data and provide them in the steCI R package.

Distributed statistical analyses often aggregate local posterior or predictive densities by repeated pairwise fusion. When the binary fusion rule is nonassociative, changing the ordered aggregation tree can change the final density and the reported posterior summaries. We study this aggregation-order variation for smooth $f$-divergence balancing in a local density chart around a common reference density. Under square-root-transformed supplied weights and first-order nondegeneracy at each generated merge, every tree in a fixed finite family shares the same first-order coefficient, whereas the first probable tree-dependent term appears at second order. We derive the explicit second-order coefficient, propagate it through every tree in the family by a recursion for the tree-indexed second-order density coefficient, and show that a finite three-state posterior contrast detects the resulting discrepancy. The same coefficient determines density-level and quantity-of-interest diameters, rotation bounds, lower bounds, a limitation of scalar calibration, and a local corrected chart representation. The resulting expansions are uniform over each fixed finite tree family satisfying the stated local conditions.

We systematically evaluate a suite of off-the-shelf variational inference (VI) software packages from the perspective of a standard practitioner. Using simple analytic benchmark models, we assess the accuracy and stability of the default VI settings in PyMC, NumPyro, and TensorFlow Probability. Unlike previous research focusing on methodological advances, our evaluation emphasizes software implementation and the default configurations that typical users encounter. Our results show that default settings can yield biased approximations of posterior summaries even for simple one-dimensional conjugate models, controls of initialization and transformations differ between software implementations, and relying on defaults may yield silent failures or poor approximations.

Clerical review of candidate record pairs remains the de facto gold standard for evaluating record linkage, but it is resource-intensive and often designed informally. We propose a design-based framework that treats clerical review as finite-population sampling over fine-grained strata defined by match weight, comparison pattern, record-level ambiguity, and demographic group. Match-probability bands are constructed from model-based score deciles. Within bands, strata combine comparison patterns with an ambiguity factor derived from matchability and conditional candidate perplexity. A band-specific margin-of-error profile encodes substantive priorities, such as tighter precision in high-score bands, while a single scaling parameter enforces the overall clerical budget. We evaluate the framework using a labelled dataset deduplicated in Splink, comprising 50,000 records and approximately 478,000 candidate pairs. We compare a baseline design reviewing about 23% of pairs with a budget-constrained design reviewing about 7%. The baseline accurately estimates global and band-specific match rates, while sampled distributions of comparison patterns, gender, and ambiguity broadly track the population. Under the budget design, global error approximately doubles, with the largest band-level errors in middle-score bands where matches, non-matches, and ambiguous cases are intermixed. Accuracy in the highest-score bands and the band-wise ambiguity profile are largely preserved, although representativeness by comparison pattern and gender declines. The framework generalises to other clerical-review objectives and can incorporate gold-standard data as prior information for design and calibration. It makes explicit the negotiable trade-offs between workload, precision, representativeness, and coverage of linkage uncertainty.

We consider mixed-integer convex optimization problems in which binary indicators control continuous variables. We introduce the Coordinate Optimality Reformulation (CORe) framework, which augments standard indicator formulations by incorporating coordinate-wise optimality information. The resulting reformulations preserve global optimality while substantially improving branch-and-bound performance, particularly in sparse and structured settings where the coordinate-wise optimality conditions expose exploitable problem structure. We first develop the main components of CORe, including coordinate-wise optimality conditions, closed-form characterizations, and disjunctive reformulations. We then demonstrate the framework across multiple problem families, including quadratic problems and robust single-index models. Computational experiments show that CORe can substantially improve solver performance compared with standard big-$M$ formulations.

Masked prediction learns to infer missing variables from visible context selected by a mask schedule. When does small excess risk guarantee recovery of the true joint distribution? We study this question on finite product spaces under masked-block log loss, with conditionals induced by a single joint distribution. To quantify recovery, we introduce an $\varepsilon$-identifiability modulus measuring the worst-case error among joint distributions with population excess risk at most $\varepsilon$. For data with separated modes pinned down by sufficiently large visible contexts, schedules that always retain such contexts can permit nonvanishing mode-weight errors while incurring exponentially small excess risk. An information decomposition explains why: the loss captures mode-weight mismatch only where the visible context leaves the mode uncertain. We prove two-sided bounds showing that, over a fixed range of mode weights, sensitivity to mode reweighting is governed by residual mode uncertainty averaged over the mask schedule. Assigning schedule mass to low-visibility masks that retain this uncertainty yields recovery bounds within the reweighting family. Beyond this family, positive full-mask probability characterizes uniform control of joint KL divergence by excess risk. Exact calculations and controlled gradient optimization validate these predictions.

The unbiased sample versions of squared distance covariance and the Hilbert-Schmidt independence criterion (HSIC) are fourth-order U-statistics, yet U-centering evaluates them from pairwise arrays in $O(n^2)$ operations. We show that U-centering is exactly the least-squares residual obtained after fitting additive endpoint effects to a symmetric hollow array. This interpretation explains the zero row sums and the denominator $n(n-3)$ through the residual degrees of freedom. The same pairwise residualization also gives useful regression identities. After endpoint effects are removed from both arrays, the U-centered dependence $t$-statistic is the ordinary slope $t$-statistic obtained by regressing one adjusted array on the other. In the two-sample problem, pooling the observations and using the between-group pair indicator as the predictor shows that the generalized-energy statistic is twice the fitted slope. The common-endpoint and fully interacted regressions give the same slope but use different residual standard errors. For $n\ge2r$, we extend the construction to arrays indexed by $r$-subsets. Higher-order U-centering removes all effects involving fewer than $r$ sample labels, leaves zero $(r-1)$-way margins, and projects onto a residual space of dimension $\binom nr-\binom n{r-1}$. For two symmetric kernels with $r$ arguments, the normalized inner product of the centered arrays is unbiased for the cross-moment of their $r$th Hoeffding components. A direct estimator can involve products spanning as many as $2r$ observations, but subset-margin inversion or higher-order U-centering evaluates the same quantity in $O(n^r)$ operations for fixed $r$. When both arrays are formed from the same kernel and sample, this becomes a nonnegative unbiased estimator of the variance of the highest-order Hoeffding component.

Predictive models employing artificial intelligence (AI) and machine learning (ML) are increasingly being used for decision support in healthcare settings. These models may exhibit differential performance across population subgroups defined by race, age, sex, and other factors and cause disparate clinical impacts, leading to intensive recent study of what has been termed "model fairness". While many methods have been proposed to assess risk prediction model fairness, these techniques generally require that the end user pre-specify the groups across which fairness is to be evaluated. In real-world settings, however, important model performance disparities may arise in unknown subgroups defined by multiple intersecting characteristics. To address this problem, we propose the unfairness tree (utree), a data-driven recursive partitioning framework for identifying subgroups with differential model performance. In simulations, the utree exhibits nominal empirical type I error rates and good ability to detect, quantify, and characterize performance discrepancies defined by higher-order variable interactions. In six mortality risk models fit to the GUSTO-I acute myocardial infarction trial dataset, utrees identified subgroup-specific performance patterns, with age, sex, blood pressure, and Killip class consistently associated with differential model performance.

We study honest adaptive confidence sets for the regression function in random-design binary regression under $L^2(dx)$ loss. Assuming only known bounds $0<c\leq g\leq C<\infty$ on the unknown design density, we construct asymptotically honest, rate-adaptive confidence sets without requiring $g$ to be smooth. Full adaptation is possible when the range of regression-function smoothness spans at most a factor of two. Over wider smoothness ranges, adaptation is achieved on the usual separated classes at the corresponding testing rates $n^{-2s/(4s+d)}$. A lower bound under the uniform design shows that these separation rates are rate-optimal. This answers a question raised by Mukherjee and Sen (2018).

Multimodal domain adaptation uses labeled source datasets to predict outcomes in an unlabeled target population. Here, different sources may observe different subsets of modalities and have distribution shifts from the target. In this scenario, the cross-modal patterns used for alignment depend on the outcome distribution, so aligning unadjusted source data can misrepresent the target. Nevertheless, the distribution shift of the outcome cannot be measured directly when labels are unavailable in the target sample and the missing modality blocks prevent simple pooling. We propose a reference-anchored domain adaptation method. A reference modality observed in every source and the target is used to estimate the target outcome distribution and reweight source observations before alignment. The auxiliary modalities are then mapped to a common, target-defined representation obtained by canonical correlation analysis (CCA) in the target and reproduced in each source by ridge-regression maps. Outcome information is transferred through density-ratio models on the aligned representation. When gold-standard source labels are sparse, a surrogate-label-assisted approach is developed to enable robust domain adaptation. We establish a unified-rotation match-up result for the target CCA and source ridge maps and consistency of the target conditional outcome distribution. Simulations and a renal cell carcinoma (RCC) application show improved calibration and stable prediction under distributional shift and blockwise missing modalities.

Detecting that a stream of high-dimensional embeddings has changed is usually framed as a choice of statistic. We give a scale law that constrains any moment-based choice and test it against topological alternatives. The law: certifying a feature of spatial scale eps carrying mass fraction f requires polynomial tests of degree N* >= log(1/f)/(2 eps), proved via the Chebyshev extremal problem; a Gauss-quadrature construction gives N* >= 4b-1 for a b-scale topology, so cost is set by feature fineness, not feature count. The law is one-sided: we exhibit an annulus whose mean, covariance and all fourth-order moments equal those of a filled disk, yet H_1 is nonzero. Its practical content is a calibration rule. The upper bound is attained by Gaussian test functions, the RKHS witness of an RBF kernel, so the law predicts which bandwidth an MMD test should use: the feature scale. On real embedding streams we measure sigma*/eps with median 1.12 (IQR 1.01-1.52, n=26) over three settings and three scales, and a data-driven bandwidth reaches AUC >= 0.95. Against an adversary optimised against the defender's statistics (mean, covariance, k-NN, kurtosis), only a bandwidth-matched kernel test still detects. For persistent homology the verdict is mixed and depends on choices usually left implicit. The summary matters more than the filtration: total persistence attains recall 0.75 at FPR 1% where the first persistence landscape attains 0.00. What survives is a cost gap, not a power gap: where persistence works it costs 116x kurtosis, which works at least as well. We conclude not that topological summaries are useless, but that on this task a kernel test whose bandwidth the law sets dominates them.

For structural components whose failure is governed by the weakest-link theory, existing reliability models typically either assume that the underlying mechanical model is known or neglect to exploit the spatial information contained in observed failure locations. In practice, however, idealized mechanical models may systematically deviate from the actual stress due to simplifying or incorrect assumptions. To address this limitation, we propose a hierarchical Bayesian weakest-link model that jointly estimates the latent material strength and stress profile from paired failure load and failure zone observations. In our formulation, the stress profile is estimated via a B-spline basis expansion, and the non-differentiable weakest-link mechanism is approximated by a differentiable Softmin function to account for unobserved material flaws and facilitate Bayesian inference. Simulation studies demonstrate that the proposed framework provides accurate and robust estimation under various experiment configurations. Applied to a real-data analysis of Douglas-fir crossarms, the proposed model identifies systematic deviations from idealized beam theory that cannot be captured by deterministic stress derivations.

Outlier detection is a fundamental task for ensuring reliable statistical modeling and inference. Interval-valued functional data (IVFD), in which each observation is represented by an interval-valued curve that preserves the variability and uncertainty within the observation, have attracted increasing attention in statistics and related applications. Developing effective outlier detection procedures for IVFD is therefore an important methodological problem. To address this issue, we develop a robust projection-based outlier detection framework. We first represent each interval-valued functional observation through its center and log-radius functions and apply interval-valued functional principal component analysis (IFPCA) to obtain a joint low-dimensional representation. We then introduce the interval-valued least trimmed functional scores (ILTFS) method, which identifies a robust reference subset by minimizing a trimmed aggregate of standardized IFPCA score distances. Finally, we proposed the ILTFS-FDR outlier detection procedure by converting the resulting projection distances into empirical $p$-values and adjusting using the Benjamini–Hochberg procedure at a prespecified target false discovery rate level. Theoretically, we derive the finite-sample breakdown point of the ILTFS mean estimator and establish the descent property of the concentration-step algorithm. Simulation studies and an empirical application to high-frequency ETF data demonstrate the effectiveness and robustness of the proposed ILTFS-FDR procedure in detecting abnormal interval-valued functional observations.

We study the Wasserstein gradient flow of the squared Maximum Mean Discrepancy (MMD) generated by the nonsmooth energy kernels $K(z)=-|z|^q$, $0<q<2$. In dimensions $d\ge2$, the corresponding energies are not displacement semiconvex, so standard Wasserstein-gradient-flow theory does not apply. When $d+q-2>0$, we prove global well-posedness on $\mathbb{R}^d$ for probability densities in subcritical $L^p$ spaces, with targets in the same integrability class and with finite moments. We also include the one-dimensional Coulomb endpoint $d=q=1$. For the associated $N$-particle system, we prove global noncollision and fixed-$N$ convergence to the collision-free critical set, a particle-to-continuum criticality principle, and a modulated-energy mean-field estimate that yields convergence of the particle dynamics to the continuum flow as $N\to\infty$ on every finite time interval. We also construct collision-free saddle equilibria, showing that deterministic particle trajectories need not approach global empirical minimizers. For $1\le q<2$, every continuum solution in our class has a narrowly relatively compact orbit, every $\omega$-limit point is Lagrangian critical, and the orbit approaches the Lagrangian critical set. For $0<q<1$, the same conclusions hold under uniform-in-time moment and subcritical $L^p$ bounds. We prove that an absolutely continuous Lagrangian critical point equals the target when the source and target have finite moments of order $q$, except when $0<q<1$ and $d\in\{1,3\}$. Under the preceding uniform bounds, rigidity gives convergence of the continuum flow to the target throughout the rigid part of the well-posedness range. Finally, we show that no initial-data-independent multiplicative MMD decay modulus exists on $\mathbb{R}^d$, and that global Polyak–Łojasiewicz inequalities fail in several whole-space and periodic Riesz/Coulomb regimes.