Papers

In many applications, covariates are naturally organized into groups, yet conventional group-selection methods often retain all variables within selected groups, failing to achieve within-group sparsity. To address this limitation, we propose an adaptive two-stage sparse-group selection procedure for high-dimensional quantile regression. The proposed method integrates convolution-smoothed quantile loss with $\ell_{2,0}$ and $\ell_0$ penalties to identify important groups and relevant variables within them, achieving sparsity at both levels. We establish theoretical guarantees for sure screening and exact support recovery at both stages, as well as consistent model-size selection using information criteria, with explicit bounds on failure probabilities. The proposed procedure is computationally efficient by restricting optimization to selected groups and active variables. Simulation studies and a breast cancer gene-expression analysis demonstrate its competitive performance in variable selection, estimation accuracy, and computational efficiency, while identifying biologically relevant genes and delivering accurate quantile predictions.

In the literature on multiple testing methods that control the familywise error rate, a fundamental concept is "coherence". A procedure is coherent if it is "logically consistent" (Hochberg and Tamhane, 1987) in the sense that if a hypothesis is rejected, then all hypotheses that imply it are also rejected. If a method is not coherent, then it can be made coherent without increasing the familywise error rate. Recently, the concept of coherence was generalized to a much broader class of multiple testing procedures, namely all methods that provide simultaneous confidence bounds for false discovery proportions. Such procedures were defined to be coherent if they cannot be further improved through a process called "interpolation". However, we show that such a coherent procedure can sometimes be further improved by exploiting remaining logical inconsistencies. This leads to what we call a "perfectly coherent" method. Thus, interpolation does not necessarily lead to perfect coherence. However, we provide conditions under which it does. Fortunately, many methods based on interpolation turn out to be perfectly coherent.

This article proposes rational quadratic kriging, a kriging estimator that addresses the limitations of existing approaches: the mean reversion issue of ordinary kriging and the diverging variance issue of rational kriging. Through the rational quadratic kernel, we establish a connection between rational quadratic kriging and kernel regression. The connection allows us to prove desirable theoretical properties of the new estimator and further propose a novel estimation procedure that utilizes factor weights from kernel regression. We also derive expressions to enable efficient leave-one-out estimation and conformal prediction. Extensive simulations across both deterministic functions and noisy datasets show significant improvements in both prediction and uncertainty quantification. A fast implementation is available in the open-source R package RQK.

In this paper, we study the estimation of the signal subspace dimension $k$ in a $p$-variate signal-plus-noise complex Gaussian model, when the noise covariance matrix is assumed to be diagonal or Toeplitz. We consider the high-dimensional asymptotic regime in which both the data dimension $p$ and the available sample size $n$ diverge to infinity at the same rate, while we also allow $k$ to diverge at a slower speed. In this context, the classical estimates of the literature, consistent in the high-dimensional regime, often rely on the a priori knowledge of an upper bound $k_*$ of $k$, which can be misspecified in practice. Our proposed approach relies on building estimates of the tracial moments of the noise covariance matrix using the $p-k_*$ smallest eigenvalues of the sample covariance matrix (SCM), which are in turn used to build a threshold above which the $k$ largest eigenvalues are isolated. The resulting estimate of $k$ is proved to be consistent in the high-dimensional regime, even for the misspecified situations where $k_* < k$. Numerical experiments in different challenging scenarios illustrate the improvement of the proposed approach over the aforementioned methods.

Nested dichotomies (NDs) reduce multiclass classification to a tree of binary classifiers. Tree structure can change accuracy, bias, and variance, but the full space of structures has not been mapped in aggregate. We introduce ndscape, a two-dimensional mapping of ND structures based on depth-weighted structural tokens. We use this embedding to conduct what is, to our knowledge, the first bias-variance analysis across the ND structure space. In simulations with 6 to 10 classes, NDs with similar tree structures exhibit similar squared bias, variance, and squared probability loss, regardless of whether the binary classifier is logistic regression or a decision tree. Variance, however, follows mean root-to-leaf depth when logistic regression is used at each node, so flatter, more balanced NDs are more stable, whereas with decision-tree node models, variance is not associated with depth and differs little across ND structures. On six datasets with 6 to 19 classes, held-out accuracy, log loss, and prediction variance remain spatially organised on the same topology-only map. None of the six standard construction heuristics tested selects a tree from the best-performing cluster on any dataset, although some select high-performing individual trees. The ndscape embedding turns ND selection from one opaque construction step into a readable comparison of class groupings, exposing the upper-level splits that distinguish high-performing and stable regions.

Coordinate descent solvers for regularized generalized linear models (GLMs) differ in how they update the intercept. This seemingly routine choice can substantially affect convergence when the response is imbalanced. For direct coordinate descent on the original GLM loss, we compare a conservative step based on worst-case curvature, a safeguarded Newton step based on current curvature, and repeated Newton steps that solve the intercept subproblem. The conservative update can leave much of the interaction between the intercept and the coefficients unresolved, producing a slowdown when that interaction controls convergence. A single Newton step avoids this attenuation and, in our experiments, follows nearly the same outer trajectory as exact minimization while requiring fewer inner steps. Simulated and real-data experiments, including comparisons across six production solvers and a controlled intervention within skglm, support this explanation. We also show how the same mechanism appears in quadratic-approximation methods and extends to multinomial models with rare classes. We recommend a safeguarded Newton update for direct coordinate descent and full intercept optimization within frozen local quadratic surrogates.

In modern sequential Monte Carlo (SMC) applications, such as diffusion model samplers, propagating a particle is expensive, so only a few particles can be afforded. In this regime, each resampling step decides which regions of the state space survive, so a better resampling scheme is worth extra computation. However, classical resampling methods only replicate existing particles and suffer from sample impoverishment, especially when the particle size is small relative to the dimension. To create new particles, we relax the unbiasedness of resampling and control the mean squared error via a Wasserstein criterion. Specifically, we propose Wasserstein-Dirichlet resampling (WDR), which constructs an equal-weight measure that optimally approximates the weighted particle distribution by solving a free-support Wasserstein barycenter problem. We then randomize the resulting optimal coupling by Dirichlet mixing. Theoretically, we provide non-asymptotic mean squared error bounds that quantify the bias-variance trade-off and show that the rate depends on the intrinsic dimension. A truncated variant preserves high-weight particles while applying WDR to the remaining mass, targeting mild weight degeneracy. Experiments on state-space models and conditional image generation demonstrate that WDR consistently improves estimation accuracy, stability, and sample diversity. The improvements persist at matched wall-clock cost when the propagation cost dominates.

Relating observations of different categories across space requires methods that accommodate distinct sampling locations and spatial dependence. Spatial kernel restricted canonical correlation analysis (skrCCA) uses spatial weights to identify associated feature combinations across data categories. Here, we develop a statistical theory of its covariance form after feature reduction. Conditional on the observed locations and under a mean-zero model, the pairwise association matrices target kernel-weighted averages of cross-covariances. For prespecified directions, we derive the exact null variance under independence between blocks. We establish recovery and stability bounds in two- and multiblock settings under explicit population models, and show that exact cyclic updates have convergent objective values. The variance characterization shows how spatial weighting and within-block dependence jointly determine the effective information available for estimating association. The recovery bounds separate statistical and optimization error, relating stability to population separation and the strength of connections between blocks. These results clarify what fitted spatial association axes estimate and provide conditions under which they can be recovered reliably.

Least squares summarizes the joint distribution of a response and regressor through second-order information. When an observed regressor is composed of latent components carrying different relationships with the response, this summary can average over distinct linear effects. We study the additional linear information contained in higher-order mixed cumulants and formalize it through the latent slope spectrum. Ordinary least squares is a variance-weighted barycenter of this spectrum, whereas higher-order cumulants can separate its components when the spectrum is sufficiently parsimonious. In the scalar one-confounder model, fourth-order cumulants identify the unordered pair $\{\beta,\beta+\gamma\}$ through a $2\times2$ linear system and a quadratic equation. We develop estimation and root-$n$ inference under regular separation and show that detection and estimation diverge near degeneracy: slope multiplicity remains detectable on an $n^{-1/4}$ local scale even when the individual slopes are not uniformly regularly estimable. In multiple regression, a quadratic tensor recurrence yields a test of one versus multiple effective coefficient vectors that does not require the vectors themselves to be identified; multiplicity can therefore be consistently detectable even when the recurrence operator is noninjective. The same quadratic diagnostic can retain signal in directions where a first-order higher-cumulant restriction cancels exactly. The results distinguish representation, identification, detection, and estimation of latent linear heterogeneity, while making explicit the limits of structural and causal interpretation from observational data alone.

Expected Shortfall (ES) measures the average size of tail losses. Its sequential backtests typically also require a Value-at-Risk (VaR) forecast. We develop two e-backtesting methods using only observed losses and ES forecasts: one for iid losses only and the other for time series under an AR-GARCH conditional VaR structure. Under the stated assumptions, both control false rejections at any stopping time. For the time-series method, validity requires a fixed parameter–quantile pair representing the conditional VaR path and a positive report margin. The time-series statistic takes an infimum over possible parameter and VaR values. We choose the betting process $(\lambda_t)$ by maximizing an empirical growth rate averaged across candidate parameters. We compute the statistic by continuous VaR profiling on a dense parameter grid. For grid approximations, we prove convergence to the continuous infimum, together with convergence of the probability of false rejection, as the grid is refined over a fixed monitoring period. Lower and upper bounds for the VaR infimum on the parameter grid control the accuracy of the numerical profiling. Simulations, including ES reports that change from over-reporting to under-reporting, and an application to daily NASDAQ losses provide numerical comparisons with VaR–ES e-processes and fixed-sample ES backtests, assessing detection frequencies and the timing of rejections.

Watermarking embeds an imperceptible signal in text from large language models that is statistically detectable with a secret key. The recent textGrain watermark (Li et. al. (2026b)) calibrates this signal against the model's sampling entropy through an explicit entropy budget, combining the stochasticity of Red-Green and the unbiasedness of Gumbel-max watermarks, and its accompanying report proposes a detector that sums pivotal statistics. We study when, and how efficiently, textGrain can be detected after human edits, modeled as a sparse mixture in which a fraction $\varepsilon_n\asymp n^{-p}$ of positions retain the watermark and token uncertainty decays as $n^{-q}$. First, entropy calibration comes at a cost: the minimax detection boundary of textGrain is $p+q=1/2$, so its detectable region is strictly smaller than that of Gumbel-max and maximal-coupling watermarks, whose boundary is $2p+q=1$. This cost stems from textGrain's construction rather than from entropy calibration itself: an unbiased scheme with the same entropy loss remains detectable where textGrain is not. Second, perhaps surprisingly, this loss makes the sum detector optimal at the boundary when the number of transport columns is fixed, unlike for Gumbel-max. When the number of columns grows and the entropy budget does not vanish, the detectable region expands while that of the sum detector does not. Third, for a fixed edit fraction, we characterize the optimal exponential rate of type II error decay and show that the sum detector falls strictly short of it at every fixed alternative. Finally, we propose a weighted multiscale binomial scan, valid in finite samples and needing no knowledge of the next-token distributions, entropy budget, or edit fraction; it attains every detection boundary above and the optimal exponent.

The Hybrid forecast covariance model is an operationally successful way to build a background error covariance matrix. It is constructed as a weighted sum, with scalar weights, of the climatological covariance matrix and a localized ensemble covariance matrix, with a localization matrix typically depending on a localization length. For best performance, the Hybrid needs to be optimized by tuning those weights and, possibly, other parameters, such as the localization length. Here, a systematic and automated method for this tuning operation is introduced. This technique is referred to as hypervariational method, and it is formulated in terms of a cost function, built using a time series of innovations as an essential ingredient. This cost function, which depends on the parameters, is minimized when they have optimal values. Equations for its minimization, based on Newton's method, are derived. Besides the aforementioned, commonly used Hybrid, the hypervariational method is applicable to any known Hybrid (and to any set of parameters on which it depends analytically), including recently developed Hybrids, not necessarily constructed as weighted sums. Here, the hypervariational method is applied to tune the ensemble covariance matrix weight in the seminal-formulation Hybrid (which does not feature localization). This is done in a synthetic experimental setting where a theoretically-derived expression for the optimal value for that weight is known. Moreover, in that experimental setting, a further method, based on a linear fit procedure, to determine that optimal value is possible. Both the hypervariational method and that linear fit procedure are shown to produce an optimal weight value matching (up to a tolerance of $0.1$) the optimal value known from the theory. Furthermore, Newton's minimization in the hypervariational method is shown to converge immediately (at the first recursive step).

The Moran effect traditionally refers to the synchronization of population dynamics driven by shared environmental forcing rather than direct interactions among populations. Although extensively studied in ecology, the Moran effect has received limited attention in geography, and its applicability has yet to be fully explored. In this study, we build upon the classical Moran effect in population ecology as a conceptual foundation, developing a more comprehensive conceptual framework in a geographic context. We introduce the new concept Geographic Moran Effect (GME), defined as temporally coherent and synchronized variation across geographic entities induced by common external forcing, with the magnitude and potentially the direction of local responses varying according to a spatially varying sensitivity field. A two-step spatiotemporal factor framework is developed to identify and interpret the GME. The framework first tests for the possible presence of GME-type synchrony in spatiotemporal data and then examines the GME and its spatiotemporal evolution by investigating its structure through spatial and temporal dimensions separately. A synthetic experiment validates the framework's reliability, while an application to U.S. county-level presidential election data (2008-2024) demonstrates its ability to identify the GME and its spatial and temporal components. By redefining the Moran effect as a geographic theory and a distinct object of inference, this work provides a mechanism-oriented perspective on spatiotemporal analysis with broad applicability in geography, ecology, and public health, where explicitly identifying and interpreting GME-type synchrony is essential for understanding system-level dynamics.

Unsupervised clustering is widely used to discover underlying group structure, after which researchers often seek to validate these groups via post-hoc hypothesis testing on the population means. However, using the same data both to form clusters and to test for a difference in means inflates the Type I error rate. Recent work in selective inference has proposed solutions when the underlying data are drawn from a multivariate Gaussian distribution with a known covariance matrix or an unknown isotropic covariance matrix. We propose a new selective test for the difference in means between two clusters that provably controls the Type I error rate under an arbitrary unknown covariance matrix. Our test statistic does not require estimating the feature covariance matrix; instead, it leverages a Monte Carlo sampling procedure over the set of orthogonal matrices. We show that this new selective test can be computed efficiently, and we derive closed-form solutions when clusters are formed using k-means or hierarchical clustering. In simulations, we demonstrate that our proposed p-value controls the selective Type I error rate and achieves substantial power, and we illustrate its application on gene expression and drug review data sets.

Bridge sampling estimates of marginal likelihoods are routinely reported with a delta-method Monte Carlo standard error (MCSE). Existing theory is asymptotic and cannot say whether a given run is long enough to trust the estimate or its MCSE. We give finite-sample answers. The optimal bridge estimator fits a logistic regression that separates posterior from proposal draws, and its accuracy is governed by their overlap, within a factor of two of the classification's Bayes error. Because the estimating function is a sum of bounded terms, no moment conditions are needed, and all bounds hold uniformly over targets, proposals and dimensions. A number of draws inversely proportional to the overlap suffices and is necessary in the worst case, paralleling the exponential-in-Kullback–Leibler requirement of importance sampling. A Berry–Esseen bound with an explicit constant holds at the optimal rate, and the MCSE is accurate once the effective overlap is large; it is estimated at no cost. For Markov chain draws, the required sample size grows in proportion to the sampler's relaxation time, unavoidably in the worst case. The estimated bridge, the fixed-point iteration and the fitted proposal matter only at second order, and we give free diagnostics for sampler failure and misestimated autocorrelation.

Although the frequentist versus Bayesian debate has a long-standing history, and both are well-established approaches for inference, they are often considered separately in pharmacometric literature, with their theoretical and practical differences receiving limited direct comparison. We address this gap by providing a conceptual background and comparative evaluation of frequentist and Bayesian methods for parameter inference in pharmacometrics. First, we highlight the theoretical and algorithmic differences between approaches within a general and pharmacometric context. We then apply both to classical compartmental pharmacokinetic models using simulated and real data; initially through evaluating four inference methods (two frequentist, two Bayesian) across a range of simulated data conditions, followed by investigating each inference approach in the practical setting. The simulation studies facilitate evaluation of the sensitivity of the methods to sample size, number of observations, residual noise and inter-individual variability, comparing the performance with respect to relative error, confidence/credible interval coverage and bias in the estimated variability. The real-world data analysis further highlights how each approach quantifies and presents uncertainty estimates by examining this for both the parameter estimates and system dynamics. By providing a theoretical and practical comparison that highlights the similarities and differences between approaches, this work aims to lower the barrier to understanding and implementing Bayesian and frequentist inference in pharmacometrics.

The Kullback–Leibler (KL) divergence is the expectation value of the log-likelihood ratio between two specified models, evaluated with respect to one of the two models. It provides an average measure of the distinguishability of the models. We examine its connection to discovery potential and the accuracy of parameter estimates when one model represents a null hypothesis, and the alternative is a family of models characterized by one or more continuous parameters. For data drawn from a multivariate Gaussian distribution, we derive the log-likelihood-ratio distribution and its cumulants under two such hypotheses. We then search over the alternative model's parameter space and use peak prominence, the height of a likelihood-profile peak above its surrounding baseline, as an operational detection statistic. Repeating this search on simulated data drawn from the null model allows us to calibrate the prominence threshold for a specified false-positive rate. We study a toy Gaussian model as well as three cosmological examples involving cosmic topology, the tensor-to-scalar ratio, and dark energy. In all these cases, detectability of a particular model depends not only on KL divergence, but also on the size of the model parameter space, the likelihood-ratio distributions, and the distance of the model from any parameter space boundaries. The relationship between KL divergence, discovery potential, and parameter estimation must therefore be calibrated for each family of models and detection statistic under consideration.

Variational inference is a flexible tool for modeling experimental data and has proven useful in crystallography. Pathwise gradient estimators provide low-variance gradients for reparameterized variational inference, but are not readily available for some distributions central to crystallography. This work derives and implements implicit and explicit reparameterized gradient estimators for the folded Normal and Rice distributions in \prog{PyTorch}, and compares their variance analytically and empirically. In each case the implicit estimator achieves lower variance for the location derivative, particularly at small location-to-scale ratios, while the scale-derivative variances are nearly equal. For both distributions, implicit estimator is proved to be the conditional expectation of the explicit estimator given the sampled magnitude, so its variance is never larger. As a practical application, French-Wilson scaling is formulated as a variational inference problem, using the folded Normal and Rice distributions as variational families for centric and acentric reflections, respectively. Both estimators recover the exact posterior to better than $1\%$ on the amplitude, and converge at the same rate. The implementations are available in the open-source package \prog{rs-distributions}.

Shapley-based measures are widely used for measuring variable importance in statistics and machine learning. However, these measures have a serious issue–a nonactive variable can get a positive importance measure when it is correlated with an active variable. Herin et al. recently proposed a measure called proportional marginal effects (PME) that mitigates this spurious effect. However, PME is developed for deterministic models and has high computational cost. In this article, we propose a new measure called proportional inclusive marginal effects (PRIME) that overcomes these two limitations of PME. The key idea in PRIME is to use a novel model-free variable selection technique to identify the active variables and include only them in the computation of the importance measure. Analysis of test functions and real datasets show that PRIME can accurately and reliably quantify variable importance from noisy input-output data.

This paper develops a mathematical framework connecting structural metadata, logical admissibility, and correction in statistical data editing. A normalized routing graph generates a canonical indexed representation and an associated Horn theory, incorporating explicit semantic edits through deterministically evaluated auxiliary variables. Under specified record semantics, the theory characterizes admissibility, distinguishes structural minimality from logical irredundancy, and permits validation without materializing the deductive closure. Canonical indices preserve the provenance and dependency scopes of violations. With complete diagnostics and strictly positive weights, a nonnegative violation Hamiltonian vanishes exactly on admissible corrections within the permitted candidate domain. Verified zero-energy outputs therefore certify feasibility and global minimization of violation energy. Minimum modification cost is formulated separately, with sufficient conditions for an equivalent penalized objective and optimality certificates obtained by matching attained costs with valid lower bounds. A household survey application verifies admissible corrections and certifies minimum cost for a subset of returned records through a retrospective audit. The framework connects structural requirements to correction certificates independently of the numerical search procedure, while statistical accuracy of replacement values requires separate assessment.

We study a stochastic optimization problem in which the random vector follows a semiparametric elliptical distribution with unknown location and shape parameters and an unspecified radial distribution. Assuming that the decision maker only has access to samples of the random vector, we develop a framework to jointly address parameter uncertainty and distributional ambiguity arising from this composite structure. To update the posterior distribution of the parameters, a key challenge is that the density of an elliptical distribution depends on the unknown radial distribution, and hence conventional likelihood-based Bayesian inference is infeasible. To overcome this challenge, we propose a novel estimating-equation-based Gibbs inference framework that updates posterior distributions of model parameters solely through estimating equations. For elliptical distributions, we construct a pair of estimating equations for the location and shape parameters that are invariant to the radial component, which enables posterior learning without specifying the density. For the distributional uncertainty, we construct an ambiguity set for the radial distribution by transforming observed samples into pseudo-samples of the radial variable. Combining the above components, we formulate a composite Bayesian risk–distributionally robust optimization model, in which the inner problem captures worst-case radial distributions within the ambiguity set, and the outer problem evaluates the decision risk induced by the parameter uncertainty through a user-specified risk measure. We further develop a cutting-plane algorithm to solve the resulting optimization problem. Numerical experiments demonstrate that the proposed approach provides substantial performance improvements over benchmark methods, including sample average approximation and distributionally robust optimization.

Selecting weights and projection directions from the observations used for testing can invalidate a reference distribution. We develop a framework for exact conditional inference under rotation symmetry. A maximal invariant permits adaptation while leaving a uniform common phase for calibration. This yields explicit Cauchy references for invariantly selected weights and projections. Under separate block rotations, the same principle gives exact normal-score combinations with coefficients selected from all blockwise invariants. The calibration results require neither densities nor moments. To illustrate the framework, we study direction learning and power in spherical Gaussian location models, identifying a regime in which consistent detection is possible without consistent direction recovery. For heterogeneous functional data, a matched mean–residual construction produces the required null symmetry under a conditional Gaussian model. Simulations compare Radius and Blend, two representative procedures, with sample splitting, and an electrocardiogram analysis illustrates their application.

Multilevel regression with post-stratification (MrP) is the dominant method for estimating public opinion in subnational units. The method requires the full joint distribution of demographic variables in a subnational unit, which is often unavailable from public census tabulations. Existing approaches fill this gap by imposing assumptions to combine marginal distributions of census variables into a synthetic joint distribution. We propose sharp MrP bounds instead: the smallest interval containing all estimates of subnational public opinion consistent with the marginal distributions alone. The bounds solve a pair of linear programs and are computationally feasible for typical applications in political science. Applied to two prominent examples in the literature, the bounds are tight enough to be informative and often more accurate than MrP restricted to a coarser, observed joint distribution of census variables. The bounds serve as margin-consistent interval estimates and a diagnostic when researchers use a synthetic joint distribution to construct MrP estimates.

We consider quickest change-point detection in generalized linear models when the post-change parameter is unknown. Such problems arise, for example, when monitoring whether the relationship between a learned representation and an observed response changes over time. Since detection must occur soon after a change, the post-change parameter must be estimated from a limited number of observations, making finite-sample estimation accuracy particularly important. We propose a window-limited CUSUM procedure that estimates the post-change parameter using a variational-inequality (VI) formulation. For noncanonical links, VI estimation can offer favorable finite-sample accuracy over short windows despite the asymptotic efficiency advantage of maximum likelihood estimation. We characterize how estimation error reduces the effective post-change information rate and affects detection delay, and derive finite-window guarantees accounting for both statistical error and incomplete numerical optimization. The resulting procedure retains false-alarm control and achieves first-order asymptotic optimality under suitable conditions. Numerical experiments compare detection delay and update cost for fixed-budget estimation and benchmark VI-WLCUSUM against alternative detectors.

Exact reversible locally adaptive step sizes are introduced by Bou-Rabee et al. (2026) in the context of the No-U-Turn Hamiltonian Monte Carlo (HMC) sampler for the purpose of high performance HMC sampling in challenging target distributions. This paper carefully explains and illustrates the locally adaptive step size procedure under more general assumptions. Alternative error measures are introduced with the aim of further enhancing MCMC performance by reducing the prevalence of reversibility failure rejections. In addition, the paper shows how to introduce the locally adaptive step sizes also within standard HMC and Metropolis-corrected kinetic Langevin/generalized HMC (GHMC). Finally, a numerical study of the relative merits of the proposed error measures in the context of different Hamiltonian flow-based MCMC methods is provided.

We propose a heteroskedasticity-robust and dimention-adaptive test for multiple linear restric- tions in sparse high-dimensional linear regression models. It combines desparsified LASSO estimation with a ridge-regularized quadratic form statistic and critical values that adapt to the restriction dimension. Under stated regularity conditions, our approach is asymptotically valid under a fixed or diverging number of restrictions without specifying the asymptotic regime in advance. Our random-design theory allows independent non-identically distributed covari- ate rows with marginal sub-Weibull tails and heteroskedastic errors satisfying finite polynomial moments. Monte Carlo simulations demonstrate improved finite-sample size control relative to conventional Wald tests, particularly as the number of restrictions increases.

We study online individualized treatment selection when only the outcome under the assigned treatment is observed. Our objective is to preferentially treat individuals whose potential outcomes under treatment exceed a covariate-specific control reference, while limiting false discoveries among treated individuals. We propose Partial-Feedback Adaptive Individualized Treatment Selection (PAITS), which initializes learning from historical data and combines sequential estimation of local null risks, randomized exploration, and feedback-driven updates of the treatment selection threshold. Under regularity conditions, we establish asymptotic false discovery rate control and sublinear gap-weighted assignment regret relative to a population oracle. We further establish a minimax lower bound on assignment regret under prescribed uniform nuisance-estimation requirements and give conditions under which PAITS attains the corresponding constrained minimax rate up to logarithmic factors. These results quantify the tradeoff between accurate treatment assignment and nuisance learning under partial feedback. Simulation studies and a real-data analysis using the Jester recommendation dataset illustrate the finite-sample performance of PAITS.

Bagging, which refers to the repeated subsampling of data and averaging the resulting outputs, is a common approach to reduce sampling variability of unstable learning and data-driven optimization procedures. However, uncertainty quantification for these bagged estimators, namely to generate confidence intervals for their target quantities, faces challenges in the simultaneous attainment of coverage validity while being computationally affordable. The latter refers to the practical consideration of overall resampling effort which, when using conventional bootstraps for instance, would lead to nested resampling procedures and thus can be multiplicatively large. In this paper, we develop an interval construction approach that is both asymptotically coverage-valid and computationally light in that it does not use much more than the resampling effort in the bagging scheme itself. Our approach builds on the recent so-called cheap bootstrap idea, a method designed to substantially reduce the number of resamples based on an integrative treatment between data noises and Monte Carlo resampling noises, and an in-depth derivation of the algorithmic configurations under which the Monte Carlo error from incomplete bagging is negligible. Our theory shows that our required computation budget is significantly smaller than existing methods in relation to the bag size and number of Monte Carlo subsamples. Numerical experiments in random forests, as well as an optimality bound estimation in stochastic optimization where bagging is useful, show that our proposed intervals achieve reliable coverage with low computation cost.

Variants of targeted maximum likelihood estimation (TMLE) and augmented inverse-probability weighting (AIPW) with machine-learning nuisance models differ in which observations train them, how bias is corrected, and whether splits are repeated. In known-truth simulations of a binary-outcome average treatment effect, eight headline estimators shared each setting's library, inner folds and truncation bound; coverage differences were judged against a prespecified two-percentage-point margin. With a low-complexity library and good overlap (n = 3000), all but a half-sample baseline covered equivalently (0.928-0.939). With one random forest on the primary data-generating mechanism, estimators fitting the propensity in-sample undercovered severely (0.12-0.48), whereas cross-fitting both nuisances raised coverage to 0.91-1.00, two overcovering. Under near-positivity, every AIPW estimator covered better than every TMLE estimator with the unweighted fluctuation (0.907-0.917 versus 0.656-0.806 at the primary setting, all undercovering); looser truncation widened the gap. At the primary settings, pooled-versus-fold-specific and single-versus-double cross-fit constructions (involving four additional estimators) were usually equivalent; under extreme weights, pooled-targeting estimates were 1.7-1.8 times as dispersed as fold-specific ones. Repetition reduced partition dependence but made several intervals more conservative. Two real-data analyses illustrate these sensitivities. Results support cross-fitting both nuisances when empirical-process conditions cannot be justified, and assessing truncation and AIPW-versus-TMLE sensitivity under poor overlap.

This paper develops a unified framework for distributional regression with arbitrary response variables. Rather than relying on scalar thresholds or Euclidean structure, we model conditional event probabilities indexed by measurable sets, which extends distributional regression to general response spaces while retaining a direct connection to classical event-based formulations. Within this framework, we study two complementary estimators: a preliminary estimator based on separate regression of event indicators, and a constrained estimator that enforces inclusion-consistent predictions through a directed acyclic graph (DAG) over a finite set family. We establish general error bounds for both estimators and further specialize the results to dense ReLU networks, yielding explicit convergence rates overcoming the curse of dimensionality. Simulation studies and real-data experiments further demonstrate the practical effectiveness of the proposed framework.