Papers

Species and feature sampling problems arise naturally whenever each observed unit is associated with one or more labels from a countable alphabet, and inference focuses on the unobserved portion of the distribution over such an alphabet. Within this framework, recent contributions have shifted attention from inference on the total probability mass over unseen labels to distribution-free confidence intervals for the largest unobserved label probability (i.e., the modal missing probability), thereby providing more refined information on whether the missing mass is concentrated on a few high-prevalence unseen labels or distributed across many negligible ones. Besides lacking a unified framework for species and features, these approaches employ a worst-case perspective that causes substantial information loss and produces overly-conservative intervals. We address these limitations through a unified model-based framework for inference on modal missing probabilities in both species and feature settings, which leverages a flexible Bayesian nonparametric formulation to localize uncertainty around models compatible with the observed data. This leads to sharper closed-form credible intervals that effectively exploit prior information and observed data, while preserving the theoretical frequentist properties and robustness of distribution-free intervals. Simulation studies confirm these improvements, while an organized crime application illustrates how our contribution has the potential to reshape law-enforcement decision-making in investigations.

We present tensor network representations for discrete maximum entropy distributions under expectation constraints. To this end, we introduce Computation-Activation Networks (CompActNets), a tensor network architecture that subsumes exponential families. By leveraging the geometry of the convex polytope of realizable expectation vectors, we represent any maximum entropy distribution in the same architecture. We exploit the fact that proper faces of this polytope correspond to the boundary closure of exponential families, which restricts the distribution's support. We then derive explicit representations for the support within the CompActNet architecture. The proposed framework suggests tensor network ranks as complexity measures for faces. Finally, a case study on Boolean statistics links the geometry of 0/1-polytopes directly to propositional formulas.

Aggregation of estimation procedures has found important applications in econometrics, statistics, and machine learning. Classical statistical aggregation theory has mainly focused on a pure aggregation setting, where the candidate estimators are either deterministic or constructed using a held-out sample independent of that used for aggregation. When the candidate estimators and the aggregation weights are estimated from the same dataset without sample splitting, Bellec [Ann. Statist. 46 (2018), 30–59] established the optimal rate for model-selection aggregation of affine estimators. Beyond this regime, however, fundamental questions regarding the optimal aggregation rates and the construction of aggregation rules attaining these rates remain largely unresolved. In this paper, we consider the problem of aggregating a finite collection of affine estimators to learn an optimal convex combination of them. This framework encompasses a rich class of estimators widely used in statistics and machine learning, including least squares estimators, kernel ridge estimators, random feature regression estimators, and among many others. We establish the minimax rate for convex aggregation of affine estimators. In particular, we show that estimating the weights by minimizing a Mallows' $C_p$ criterion attains the optimal rate. We further study the linear aggregation regime with unrestricted weights and establish matching minimax lower and upper bounds over suitable classes of affine estimators.

In clinical diagnosis, sequential screening, and cascaded prediction systems, information is often acquired in stages, with intermediate predictions guiding whether to act on the available evidence or defer action to gather further information. We study sequential predictive inference in these settings, aiming to construct a sequence of prediction sets that supports timely, resource-efficient decisions while remaining reliable throughout the acquisition process. Our objective is to maximize process-level utility subject to any-stage validity, which controls the probability of miscoverage anywhere in the prediction sequence. We develop COINS, a conformal framework for constructing such sequences, and establish a universality theorem identifying the common exclusion structure underlying all nested, calibration-symmetric, any-stage valid prediction sequences. Building on COINS, we develop Vopt-COINS for data-driven optimization of prediction sequences under task-specific utilities based on prediction-set size, stopping time and acquisition costs. Branchwise and localized extensions further support individualized acquisition and calibration. Numerical studies demonstrate the effectiveness and practical benefits of the proposed framework for reliable, resource-efficient sequential prediction.

We prove an impossibility theorem for uniformly consistent recovery of labeled edge probabilities from complete aggregated relational data, even when the population count law identifies every probability. For any known partition into two equal trait groups, we consider an independent-edge logistic network with a balanced rank-one signal and unknown activity effects with bounded total dyadic energy. The signal amplitude is known and fixed at a small positive value. Every node reports its counts to both groups. The minimax mean squared error for the probability matrix remains bounded away from zero as the network grows, whereas it tends to zero under full adjacency on the same parameter class. The lower bound accounts for the dependence across the entire count array through a two-node oracle from which all observed counts can be reconstructed. Conditional binomial smoothing bounds the information about a local sign orientation, and an anchored many-bit construction converts these ambiguities into a nonvanishing normalized matrix loss. A mixed-cumulant identity recovers every edge probability from the population degree law, locating the obstruction in estimation from a single aggregated network rather than in population identification.

Let $W=n^{-1/2}\sum_{i=1}^n X_i$, where the $X_i$ are independent centered random vectors in ${\mathbb R}^p$ with $|X_{ij}|\le B$ almost surely. Suppose that $\text{Cov}(W)$ has unit diagonal and smallest eigenvalue at least $b^2>0$. We prove that the distance between $W$ and a Gaussian vector with the same covariance, uniformly over axis-aligned rectangles, is at most $C\min\{1,b^{-2}Bn^{-1/2}\log^{3/2}(ep)\}$. For fixed $b$, the dependence on summand size and dimension matches known lower bounds in growing-dimensional regimes. The proof combines a concentration estimate near rectangle boundaries with a carefully chosen Gaussian comparison.

Estimating the first stage of an instrumental variables (IV) model with the least absolute shrinkage and selection operator (LASSO) requires choosing a dictionary of technical instruments and a penalty level. First-order asymptotic theory offers no guidance on these choices, as any consistent implementation yields a structural parameter estimator with the same limiting distribution. In finite samples, however, these choices can have a substantial impact on the resulting structural parameter estimate. Working in a model with a single endogenous regressor and homoskedastic Gaussian errors, we use first- and second-order Stein identities to derive the approximate mean squared error (AMSE) of the instrumental-variables LASSO (IV-LASSO) estimator, which can be consistently estimated and used to rank a prespecified list of dictionary-penalty candidates. The AMSE reveals a bias-variance trade-off: more complex first-stage fits better approximate the conditional mean of the endogenous variable but are also more correlated with the structural errors, with complexity measured by the degrees of freedom of the LASSO fit. The weight on this bias rises with the endogeneity of the regressor, a quantity that neither plug-in nor cross-validation penalty rules take into account. Despite the AMSE being derived in a Gaussian model, penalty selection by minimizing the feasible AMSE criterion delivers up to a one-third lower mean squared error compared to cross-validation and plug-in penalty rules in Gaussian and non-Gaussian simulation designs calibrated to the data of Gilchrist and Sands (2016).

Learning with group invariances is central to many scientific and geometric learning problems, yet its computational foundations remain poorly understood. Even for classical supervised regression settings, it has been unclear whether one can efficiently compute a regression function that is exactly invariant to a given group action. Recent work showed that exact invariance can be enforced in polynomial time when the underlying group is finite and known, but left open the cases of infinite groups and unknown symmetries. In this paper, we resolve both challenges. First, we present the first polynomial-time algorithm for learning with exact group invariances that applies uniformly to finite and infinite groups. The runtime is polynomial in the data dimension and sample size, and independent of the group, while achieving strong generalization guarantees. This provides a computational explanation for the empirical success of invariant and equivariant methods in geometric machine learning and partially answers a recent open question in the literature. Second, we study learning in the symmetry discovery setting, where the invariance group is unknown. Focusing on the subgroup lattice of a finite group, we show that exact symmetries can be identified from data and exploited for learning in polynomial time. For regression over finite-dimensional feature spaces, our algorithm provably recovers the underlying symmetry, matches the minimax-optimal sample complexity of the known-symmetry setting, and runs in time polynomial in the data dimension and sample size. Our analysis relies on tools from random Cayley graphs and expander theory, which may be of independent interest.

In this paper, the exact joint density of the top $k$ sample eigenvectors is derived for any $p\times p$ population covariance matrix $\varSigma$, sample size $n > p -1$, and $1 \le k \le p$. Using symmetric functions such as zonal polynomials and a dual summation identity by I. G. Macdonald, the result is expressed as a series in terms of determinants of differential operators. A matrix Kummer transformation converts the resulting local alternating expansion into a globally absolutely convergent series over the entire positive definite matrix space. For $k=2$, the ordered eigenvalue integrals reduce to the Gauss hypergeometric function ${}_2F_1$ with a simple closed form when $n=p+1$. Previously, explicit formulas were only available for a scalar matrix $\varSigma = \sigma^2 I_p$ or a general matrix $\varSigma$ with $p=2$ or $k=1$, obtained by T. W. Anderson and T. Sugiyama, respectively. The frame law also induces an exact Grassmann density that provides a finite-sample benchmark for principal subspace inference.

Pooling latent subpopulations can obscure relationships and yield misleading regression conclusions, including Simpson's paradox (SP). We propose a detector based on the steady-state dynamics of constant-step stochastic gradient descent (SGD). Unlike likelihood-based mixture tests and confounder-search methods, it requires neither a normal-mixture specification nor observed candidate confounders. Two linear pieces compete under a winner-take-all squared-error loss, and their normalized terminal separation forms the test statistic. Using diffusion approximations, we derive its asymptotic null distribution for general centered scalar covariates and for Gaussian multivariate covariates under symmetric noise. The distribution-dependent null center reduces to the dimension-free constant $4/\pi$ when both covariates and noise are Gaussian. We establish asymptotic size control and consistency against fixed mixture alternatives under regularity conditions. We extend the method to intercept and partial-mixture heterogeneity and study endogeneity, heteroskedasticity, and nonlinear misspecification. Simulations examine calibration, power, and robustness. Finally, a three-stage Detect–Screen–Verify toolkit separates evidence of heterogeneity from its substantive explanation. Across eight public datasets, it recovers four established SP benchmarks and identifies four cases that, to our knowledge, have not been documented previously. The detector requires neither latent-group labels nor a prespecified number of mixture components.

We prove a sharp lower bound for growing-degree polynomial estimation in the symmetric rank-one matrix estimation problem. The observation is \[ \boldsymbol{Y}=\frac{1}{\sqrt{n}} \boldsymbol{\theta} \boldsymbol{\theta}^{\top}+\boldsymbol{W}, \] where the coordinates of $\boldsymbol{\theta}$ are independent with a common, fixed sub-Gaussian distribution $\pi_\Theta$ satisfying $\mathbb E[\Theta]\neq0$, and $\boldsymbol{W}$ is symmetric with independent standard Gaussian upper-triangular entries. For every such prior, we give an explicit finite-dimensional bound implying that every sequence of polynomial estimators of degree $D(n)=o(n^{1/60})$ has normalized mean-square error with limit inferior at least $\mathbb E[\Theta^2]-q_{\mathsf{amp}}$, the limiting error of Bayes approximate message passing (AMP). This extends the constant-degree equivalence of Montanari and Wein~\cite{montanari2025equivalence} to growing degree under the same prior assumptions. Combined with their polynomial approximation of fixed-iteration AMP, the bound identifies the exact limiting low-degree MMSE whenever $D(n)\to\infty$ within this range. It thereby answers the growing-degree AMP-equivalence question discussed in~\cite{wein2025computational, maleki2026high} for nonzero-mean sub-Gaussian priors in the stated degree range. We also obtain a moment-dependent finite-dimensional lower bound for any nonzero-mean prior with all moments finite. The proof bounds low-degree MMSE using conditional joint cumulants and a refined application of Jensen's inequality. The conditional cumulants are defined using an auxiliary Gaussian channel $\boldsymbol{R}$ calibrated to the AMP fixed point. This retains signal dependence that is lost in unconditional cumulant bounds and produces the cancellations needed for quantitative control as the degree grows.

In the recent work of Barbier, Ko, and the second present author on sublinear-rank symmetric matrix factorization [Math. Stat. Learn. 9 (2026), 1-68], a key result is that, in the Bayes-optimal setting, the large-size limit of the free entropy of the finite-rank spiked Wigner model is the same as in the rank-one case when the signal has centered i.i.d. entries. In this paper, we show that this rank-one equivalence result extends to the case of finite-rank, even-order, symmetric tensor factorization. Moreover, we give a natural reformulation of a hypothesis that was stated in the aforementioned work to be necessary for this result. As in the matrix case, we use information-theoretic identities and replica symmetry to reduce a known multi-dimensional variational formula for the limiting free entropy to its one-dimensional analog. The novelty stems from the fact that said formula involves a replica symmetric potential containing Hadamard (entrywise) powers, rather than squares, of the matrix-valued variational parameter, so the eigenvalue-based approach used in the matrix case must be adjusted.

The Brier risk summarizes probability-prediction error but does not explain why two prediction rules differ. We develop Brier-UDAM, a framework for paired binary prediction-rule comparisons that separates the Brier-risk difference into mean-bias, dispersion, and Brier-scale alignment contributions. Although the observable contrast identity follows algebraically from Yates-type moment decompositions, the framework gives these terms a latent-risk interpretation and shows that paired contrasts can be identified from observed outcomes and paired predictions without estimating the conditional event-risk function. For two rules evaluated in the same target population, the exact identity $\Delta R = \Delta M + \Delta D + \Delta L$ holds, and the corresponding plug-in estimators reproduce the empirical paired Brier-risk difference exactly in finite samples. We derive influence-function-based inference for first-order regular contrasts and a projection-based confidence procedure for the mean-bias contrast when first-order inference degenerates. Simulations confirmed interpretable component behavior, recovery of the target contrasts, and generally satisfactory inference in regular settings; projection intervals retained coverage near degeneracy but were conservative. In a large public-data analysis of 30-day hospital readmission, models with similar held-out Brier risks showed substantially larger and opposing dispersion and alignment contributions, revealing cancellation that was hidden in the aggregate score. Brier-UDAM therefore provides a diagnostic account of paired Brier-risk differences, while calibration, discrimination, and information-based assessments continue to address their distinct performance questions.

We develop a new low-accuracy sampler, called smoothed Picard Hamiltonian Monte Carlo, which combines Gaussian smoothing, Picard iteration, and higher-order discretization. For a log-concave target $\pi \propto \exp(-V)$ in dimension $d$ satisfying $0 \prec \alpha I \preceq \nabla^2 V \preceq \beta I$, with condition number $\kappa := \beta/\alpha$, smoothed Picard HMC returns a sample with $\sqrt \alpha\,W_2(\cdot,\pi) \le \varepsilon$ using $\widetilde O(\kappa^2 + \kappa^{7/6} d^{1/6}/\varepsilon^{1/3})$ gradient queries. We also prove stronger $W_q$ bounds, and then develop an algorithmic framework, the recursive warm start generator, to upgrade these $W_q$ bounds to stronger divergence guarantees. This produces a warm start for the proximal bouncy particle sampler, introduced in a companion work, leading to a high-accuracy log-concave sampler with complexity $\widetilde O((\kappa^{7/6} d^{1/6} + \kappa^{1/2} d^{1/4})\mathrm{polylog}(1/\varepsilon))$.

We study the problem of sampling from $\mu(\mathrm{d}x)\propto e^{-V(x)}\,\mathrm{d}x$ on $\mathbb{R}^d$, where $V$ is $\alpha$-strongly convex and $\beta$-smooth, and write $\kappa:=\beta/\alpha$. We design and analyze the Proximal Bouncy Particle Sampler (Proximal BPS), a new sampler that combines ideas from the proximal sampler and the bouncy particle sampler. From a warm start initialization with $ O(1) $ Rényi divergence w.r.t. $\mu$, Proximal BPS returns a sample whose law is $\varepsilon$-close to $\mu$ in total variation distance using $\widetilde O(\sqrt\kappa\,d^{1/4} \,\mathrm{polylog}(1/\varepsilon))$ gradient queries in expectation.

Randomized dose-optimization trials may screen several candidate doses using binary efficacy and toxicity outcomes. A dose is inadmissible if efficacy is insufficient or toxicity is excessive, so each dose-specific null hypothesis is a union region and strong familywise error control must hold across arbitrary mixtures of inadmissible and promising doses. We study exact verification of multistage monitoring rules in which armwise decisions may be coupled through the number of active arms remaining at each interim analysis. For uncoupled rules, we derive an exact product representation and show that a complete-null boundary configuration is least favourable. Under active-count coupling, this factorization no longer holds because a promising dose can remain active and alter future boundaries for null doses. Assuming independent sampling across arms, prespecified per-dose analysis schedules, analysis timing not driven by the observed efficacy or toxicity outcomes, arm-specific monitoring statistics, and stagewise monotonicity, we establish an exact finite-boundary characterization without parametric restrictions on the within-arm joint efficacy-toxicity distribution. Each dose need only be evaluated at an efficacy-null boundary, a toxicity-null boundary, or a maximally favourable alternative. A deterministic finite-state recursion then verifies a fixed decision table without Monte Carlo error. Numerical studies confirmed the prespecified strong familywise error control and showed that active-count coupling can shift the least-favourable configuration away from the complete null while improving joint retention of multiple promising doses. A published randomized dose-ranging trial was used as a clinical illustration of how the framework could be prospectively implemented. The framework separates monitoring-rule construction from rigorous error verification.

The strictly standardized mean difference (SSMD) is an interpretable, unit-free effect size originally developed for quality control and hit selection in high-throughput screening and subsequently extended to effect-size assessment and error-rate control. Existing SSMD theory considers one variable at a time. Here, we extend SSMD to a linear combination of multiple variables and determine the coefficients that maximize the resulting SSMD. We show that this optimization can be formulated as a generalized Rayleigh quotient with a closed-form solution and identify conditions under which the optimal direction coincides with Fisher's linear discriminant. We derive method-of-moments, bias-corrected, and maximum-likelihood estimators of the maximized SSMD and construct confidence intervals based on Hotelling's T-squared statistic and the noncentral F distribution under different correlation and variance structures. We also characterize classification cutoffs for the resulting linear score and establish its relationship with the area under the receiver operating characteristic curve (AUROC). Simulations show that SSMD maximization coincides with Fisher's linear discriminant and logistic regression under equal covariance but diverges under unequal covariance and class imbalance. In paired designs, SSMD maximization uniquely incorporates cross-covariance information. Linear support vector machines are also evaluated as a modern comparator. This framework provides an interpretable approach to constructing multivariable biomarkers and signatures for high-throughput screening, cytokine profiling, metabolomics, salivary diagnostics, and continuous-monitoring applications.

We develop an expected shortfall factor model (ESFM) to estimate and price common variation in the severity of lower-tail losses in large panels of asset returns. Mean factor models describe common variation in average returns, while quantile factor models describe common movements in tail thresholds. ESFM instead captures common variation in the average severity of losses below those thresholds. The model combines observed risk exposures with latent common factors. We estimate ESFM using an orthogonalized two-step procedure under which first-stage quantile estimation error has no first-order effect on the ES coefficient estimates. We establish nonasymptotic error bounds for the ES coefficients, a finite-sample Gaussian approximation, and consistent selection of the number of latent factors. Applied to a large panel of equities, ESFM uncovers common factors that react sharply to market stress and contain information not captured by mean and quantile factors. Average returns increase across portfolios sorted on ESFM exposure; high-minus-low portfolios earn annualized returns of 8.0%–11.7% and Fama–French five-factor alphas of 10.3%–15.0%. These spreads remain positive and statistically significant after conditioning separately and jointly on mean- and quantile-factor exposures. Tail-by-tail spanning tests show that ESFM factors retain significant alphas after controlling for standard traded factors and the corresponding mean and quantile factors. Adding ESFM to these benchmark factor sets increases the maximum attainable Sharpe ratio. These findings identify common loss severity as a distinct and priced dimension of downside risk.

The Davis-Kahan theorem is a fundamental tool in spectral analysis, providing quantitative control over the distance between the eigenspaces of a symmetric matrix and its perturbation. However, when the matrix dimension is large, computing leading eigenvectors is computationally expensive, limiting the practical use of spectral methods in modern large-scale applications. This paper addresses this problem by proposing an independent Bernoulli sampling scheme and proves that the leading left singular vectors of the subsampled matrix faithfully approximate the target subspace of a low-rank symmetric matrix. Our main result is a subsampled Davis-Kahan bound that gives an explicit error bound depending directly on the sampling probability. The bound reveals the trade-off: the computational cost scales linearly with the sampling probability, while the statistical error scales as the inverse square root of the sampling probability. Our result thus extends the Davis-Kahan theorem to the subsampled setting, enabling scalable spectral analysis of large-scale symmetric matrices.

We study how a limited labeling budget should be allocated to minimize multiclass zero-one classification risk. We consider parametric classification problems in which features are observed for all sampling units while class labels can be acquired selectively. By combining the Fisher information supplied by an acquired label with the local geometry of multiclass excess risk, we derive an acquisition criterion that minimizes the leading asymptotic coefficient of expected multiclass excess risk. The resulting rule values a label according to how strongly its information is aligned with parameter directions that perturb the active Bayes decision boundary, rather than according to posterior uncertainty or global parameter information alone. We characterize the oracle acquisition design, establish its threshold structure, and derive face-specific and cost-sensitive extensions. An analytic example shows that posterior uncertainty and classification value can produce different, and even reversed, acquisition rankings. We further develop a two-stage adaptive procedure that attains the oracle leading-risk criterion under regularity conditions and provide explicit results for Gaussian discriminant analysis. Three-class QDA experiments illustrate the resulting acquisition geometry, while an application to the six-class Statlog Landsat Satellite data shows that classification-risk acquisition can differ materially from both uncertainty-based acquisition and the complete-classification-information comparator. The adaptive classification-risk design attains lower mean error than this Fisher comparator across the labeling budgets considered, although it does not uniformly outperform entropy or margin sampling and differences among the targeted strategies become small as the labeling budget increases.

We present a test statistic for the comparison of DNA sequences under some of the most popular evolutionary processes available in the literature. Theoretical properties for the test statistic as well as its empirical performance by stochastic simulations are presented. The proposed test statistic is a generalized $U$-statistics built for tests of distributional homogeneity under the null hypothesis. We show that a dicothomous situation exists here. Under the null hypothesis, the $U$-statistics kernel is first-order degenerated, this test statistic falls in the quasi $U$-statistics class and follows an asymptotic normal law, albeit of higher order than the standard case. Under heterogeneity, the asymptotic normality is attained on the more usual first-order asymptotics. Asymptotic normality is proven for the cases: high-dimension/large sample size, high-dimension/small sample size, low-dimension/large sample size. Moreover, the case of local alternatives is discussed, and the contiguity of the test statistic for them is established. Simulation studies are performed to assess some finite-dimensional properties of the test statistic, regarding issues such as balanced/unbalanced samples, dimension and sample size.

Cluster randomized crossover (CRXO) trials are often used when individual randomization is impractical and the number of available clusters is limited. However, statistical analysis of CRXO trials is complex because of the need to account for complex correlation structures over time. It becomes especially challenging when very few clusters are used because standard modeling assumptions may lead to unstable variance estimates, poor confidence interval coverage, and inflated type I error. This study evaluates individual-level mixed-effects and fixed-effects models with and without a cluster-period random effect, cluster-period summary analysis using normal- or \(t\)-based inference, and two-period crossover-difference estimators. Using extensive simulation studies under both nested exchangeable and discrete time decay correlation structures, we compare model performance in terms of bias, root mean squared error, coverage probability, type I error, and convergence. Across scenarios, all models produced approximately unbiased treatment effect estimates, but their inferential performance differed substantially. Models that explicitly accounted for cluster-period heterogeneity generally provided the most reliable control of coverage and type I error, whereas simpler exchangeable models performed adequately only when the true correlation structure closely matched their assumptions. Cluster-period level analysis performance improved with increasing numbers of periods but was unreliable in the sparsest designs. Overall, the findings suggest that in CRXO trials with very few clusters, accurate modeling of cluster-period correlation is more important than the choice between fixed and random cluster intercepts, and that results from extremely sparse designs should be interpreted with caution.

We propose the Hypercomplex Irregular Autoregressive (\HIAR) model, a quaternion extension of the \IAR/\CIAR/\BIAR{} family for four-component observations at irregular times. Temporal dependence is represented by a real power of a quaternion parameter and estimated in state-space form through the Gaussian innovation likelihood of the Kalman filter, conditional on the adopted covariance specifications. Across 12,000 Monte Carlo fits, mean absolute bias decreased with sample size and ranged from 0.0009 to 0.0023 at $N=300$. We applied the model to 30,824 Sentinel-2 pixel series covering the Mata de Santa Genebra ARIE from 2020 to 2023. The optimizer reported successful numerical termination in 98.73% of the fits; the median $\|{\Phihat}\|$ was 0.6248, 3,872 pixels met the operational high-persistence threshold ($\|{\Phihat}\}\geq0.95$), and the median residual RMSE for band B8 was 5.1551 percentage points. The results demonstrate the computational feasibility of \HIAR{} and its ability to generate descriptors of multispectral persistence, vector dynamics, predictive error, and spatial discontinuities potentially associated with forest edges.

We consider dimensionality reduction for high-dimensional observations accompanied by a supplied partition into two or more clusters. The objective is not to construct a low-rank projection, but to retain an interpretable subset of the original coordinates that preserves the distributional information distinguishing the clusters. For each coordinate, the proposed procedure compares the cluster-specific empirical distribution functions through a several-sample Kolmogorov-Smirnov separation statistic. We formalize the resulting marginal cluster support and establish simultaneous finite-sample concentration over all coordinates, explicit bounds for false inclusions and omissions, and exact support recovery when the minimum distributional separation dominates the high-dimensional stochastic error. We also quantify the dimension inflation induced by using an unadjusted testing level and give a familywise-error-controlled version. Under a conditional sufficiency condition, sure screening preserves the full-data posterior cluster probabilities, mutual information, and Bayes risk; an additional result characterizes robustness to imperfectly estimated cluster labels. The procedure is invariant to strictly increasing coordinate transformations and can retain low-variance cluster signals that principal components may discard. We further develop average dual information, a criterion combining partition agreement after transformation with structural coverage of cluster-relevant coordinates, and derive its basic properties and consistency. Simulations illustrate the theory, the interpretability of the selected coordinates, and the distinction between cluster-directed screening and variance-directed projection.

The mean squared displacement (MSD) is a cornerstone in the analysis of diffusion processes in complex media. When the system is heterogeneous and, in particular, when single-particle trajectories are short, it is essential to extract maximal information from each measured trajectory. This is typically done by time-averaging squared increments and examining the scaling of the time-averaged MSD in log-log space. However, classical regression methods perform poorly in this setting because time-averaging introduces correlations aggravated by those inherent to anomalous diffusion. We tackle these limitations by applying a generalized least squares framework, which substantially reduces variance and bias in diffusion parameter estimates, especially for short (around 100 points) and ultra-short (around 10 points) trajectories. The method is fully automated and requires no supervision. Furthermore, it enables prediction of estimation error probability density, which is asymptotically Gaussian, for both classical and enhanced approaches. Leveraging this prediction, we introduce a specialized deconvolution algorithm that reconstructs the underlying particle ensemble structure from experimental data.

Heterogeneous-agent New Keynesian (HANK) models characterize how entire cross-sectional distributions respond to structural shocks. Traditional representative-agent models are routinely disciplined by impulse responses from aggregate vector autoregressions (VARs). HANK models have no comparable established empirical benchmark because they make predictions not only about aggregates, but also about distributions of micro-level data. We propose a Bayesian benchmark that jointly models macroeconomic aggregates and several marginal distributions from repeated cross sections, including distributions observed in different surveys. Our approach can use both standard structural VAR identification approaches on macroeconomic aggregates and identification restrictions imposed on micro-level data. The model delivers a joint posterior of the distributional effects of shocks, without the need for household panel data or a separate first-stage density estimate.

Sports betting moves money continuously from a large public to a small number of firms. The most influential account of that flow, due to Levitt (2004), holds that books price away from the market-clearing point to exploit predictable public biases. It computes profit from two numbers (the probability that a side wins the proposition and the fraction of handle it attracts), treating the share on a side as independent of the outcome. Here we relax that assumption and derive a profit-bias identity: profit is affine and increasing in the expected share of handle on the losing side, with Levitt's expression as the special case of independence. The identity resolves the book's margin into exactly three channels: its hold, the product of its price shading and the public's lean, and the covariance between bet share and outcome. A lean is thus worthless without shading, and shading is worthless without a lean. Under a public-belief model, the profit driver is the public's Bayes error: the probability of a representative bettor selecting the losing side. We provide necessary and sufficient conditions for a "Goldilocks Zone": prices at which book and bettor both profit. Testing these predictions on 1,139 Major League Baseball games, we find that the apparent dependence between bet share and outcome is a Simpson's paradox: present when games are pooled, but absent once they are separated by which side the book favored. The public leans heavily toward favorites, but we detect no matching shading, and the realized margin is indistinguishable from the hold.

We study trend estimation in state-space models in which the trend has a fractional stochastic difference of order $d>0$ and the observation errors form a short-range-dependent stationary process. Using finite-sequence fractional summation and differencing operators, we analyze the penalized least-squares estimator obtained by shrinking the fractional differences of the trend. We derive asymptotic mean squared error (MSE) approximations for all $d>0$ and identify a sharp phase transition at $d=1/2$. When $d>1/2$, the estimator is consistent and its optimally balanced MSE has order $n^{-(2d-1)/(2d)}$. At the boundary $d=1/2$, we obtain a refined finite-sample approximation and show that the MSE decreases at the slower order $\log\log n/\log n$. When $0<d<1/2$, the MSE converges to an explicit positive limit, so consistent recovery of the trend is impossible under the considered scaling. We also describe a practical criterion for choosing the penalty parameter and differencing order, and numerical experiments illustrate the MSE approximations and the behavior of the selection procedure.

Bayes estimation has been extensively studied and widely used in statistics, decision theory, signal processing, machine learning, and system identification. Among its variants, empirical Bayes (EB) estimation has attracted considerable attention due to its favorable estimation performance and computational tractability. However, the direct plug-in dependence of an EB estimator on hyperparameters can make it locally sensitive to hyper-parameter perturbations. This paper considers the linear regression model and focuses on the EB estimator by employing the marginal maximum likelihood hyper-parameter estimator. For conciseness, this estimator is simply referred to as the EB estimator. Given a family of EB weighting functions, a generalized Bayes estimator is constructed with the same excess mean squared error (XMSE) as the corresponding EB estimator. Here, the XMSE is a second-order asymptotic measure of the mean squared error difference between the estimator of interest and the maximum likelihood estimator. Furthermore, the EB estimator is shown to be at most firstorder sensitive to hyper-parameter perturbations, whereas the constructed Bayes estimator is at most second-order sensitive, making it locally more robust. The computational complexities of these two estimators are also analyzed. In some cases, the constructed Bayes estimator can be computationally comparable to, or more efficient than, the EB estimator. These theoretical results are further supported by numerical simulations.

We derive an exact gradient-step representation of the RoPE-softmax forward pass. For every deterministic RoPE-softmax attention head with arbitrary affine projection weights, we construct a query-dependent effective matrix $\Delta M_i$ satisfying $y_i = \mu_i + u_i^\top \Delta M_i$, where $\mu_i$ is the uniform mean of the attended values and $u_i$ is the augmented query input. The construction applies the classical exponential divided difference $\rho = \phi_1$ to retain the softmax exactly. Its positive coefficients give a unit gradient-step representation on a query-conditioned quadratic objective. The same function connects the RoPE generator to exact positional finite differences. We derive a tokenwise formula for the error of reusing one query's matrix and prove that a nonconstant finite-cache head cannot admit a globally exact affine query readout. Reconstruction checks and frozen-reuse calibration on one pretrained Qwen2.5-0.5B layer verify the representation and quantify the correction required when one query's matrix is reused.