Papers

Reinforcement learning and contextual bandit algorithms have become increasingly common in sequential decision-making applications. When these methods are deployed in high-stakes domains, there is growing interest not only in learning effective policies, but also in conducting statistical inference for quantities learned under adaptive data collection. However, classical procedures applied naively in these settings can fail: even when estimators are unbiased, their variance becomes path-dependent and as a result may not be asymptotically normal. A growing literature has emerged to ameliorate this problem, but solutions tend to be problem specific and often rely on correct specification of a working model. In this work, we develop a unified framework for constructing asymptotically valid confidence intervals to cover smooth functionals of nonparametric M-estimands under adaptive sampling. Under Neyman orthogonality, we provide two novel methods for performing inference: (1) a self-normalized statistic based on the realized quadratic variation of the influence function and (2) a statistic using a plug-in estimate of the conditional variance based on reweighted influence function increments. Our results allow for flexible nonparametric estimation of nuisance parameters and remain valid under model misspecification. Our theory is supported by a simulation study for a dynamic pricing application which demonstrates that this method can produce asymptotically valid confidence intervals where standard methods fail.

Regression analysis for stochastic processes has been an important topic over the last decades. Typically, regressing a dependent process $Y$ on an explanatory $Z$ results in an integral relationship of the form $dY_t = \beta_t dZ_t + dX_t$ with some residual process $X$, and often with constant or piecewise constant $\beta$. In the case of jump processes such a relationship essentially boils down to $\Delta Y_t = \beta \Delta Z_t + \Delta X_t$ on the level of the jumps. Here the residual process $X$ is such that it never jumps together with $Z$, and so any jump in $Y$ is either exactly proportional to $Z$ or purely idiosyncratic. In this paper we discuss a related setup which appears more realistic and is closer to classical regression, namely $\Delta Y_t = (\beta + \eta_t) \Delta Z_t$ for i.i.d. $(\eta_t)_{t \ge 0}$. We propose an estimator for $\beta$ based on high-frequency observations of $(Y,Z)$ and discuss consistency and associated central limit theorems in two different asymptotic regimes.

For $\beta \in \mathbb{R}$, we define the partition function $$ Z_\beta(X) = \sum_{i=1}^N \exp(\beta X_i) $$ of a centered Gaussian random vector $X=(X_1, \ldots, X_N)$ with $\mathbb{E}[X_i^2]= 1$. For $q \in \mathbb{R}$, we prove a complete phase transition for the generalized $L_q$-means of $Z_\beta(X)$ at $q=1$. The centered Gaussian vector with covariance matrix $\Delta_N$ obtained from a regular simplex configuration of unit vectors maximizes these means for $q<1$, and minimizes them for $q>1$. The two regimes are governed by different principles. For $q > 1$ {and $\beta\ne0$}, the moment functional is globally strongly convex on the entire set of correlation matrices, with an explicit modulus of convexity and a quantitative centroid-shape stability estimate. For $q<1$, the strategy is different. We prove a universal comparison for log-concave profiles of reverse Brascamp–Lieb type. Specializing this result to the Gumbel profile yields a Laplace-transform comparison between the partition functions.

We study the problem of sampling from a continuous density $\pi\propto \exp(-V)$ on $\mathbb R^d$, where $V\in C^2(\mathbb R^d)$ has a $\beta$-Lipschitz gradient and $\pi$ satisfies a logarithmic Sobolev inequality with constant $\alpha^{-1}$, and write $\kappa = \beta/\alpha$. We introduce the Gaussian cloud sampler, which achieves total variation accuracy $\varepsilon$ using $\widetilde O(\kappa d^{1/5}\,\text{polylog}(1/\varepsilon))$ gradient queries in expectation. The algorithm uses first-order rejection sampling (FORS) to correct the law of smoothed Picard HMC trajectories. To do so, we represent the iterates of the ideal Picard iteration via Gaussian clouds, whose centers are never evaluated, and we develop a suite of likelihood correction gadgets which only use samples from this indirect cloud representation.

Imputation-based causal estimation is typically viewed as relying exclusively on an outcome model, in contrast to augmented inverse-probability weighting, whose consistency is protected by fitting two nuisance models. This paper argues that this view can be misleading by highlighting a hidden dual weighting structure in least-squares sieve regression imputation. Although only outcome regressions are explicitly fitted, the resulting imputation estimator admits an exact weighting representation whose induced weights balance every function in the sieve space and the corresponding population weighting functions are the $L^2$ projections of the inverse propensity scores onto the same sieve space. This projection structure yields an implicit form of double robustness and, under standard sieve approximation and growth conditions, asymptotic linearity with the efficient influence function. Thus, weighting, covariate balance, double robustness, and semiparametric efficiency can all emerge from imputation alone through the geometry of least-squares projection.

We develop a three-stage adaptive procedure for principal component analysis (PCA) when high-dimensional observations are collected sequentially and additional sampling incurs a cost. The procedure balances PCA compression loss against sampling cost while selecting the retained dimension through a prescribed explained-variance criterion. Starting from a pilot sample, an intermediate stage updates the PCA quantities before determining the final sample size, thereby avoiding reliance on unknown population eigenvalues. Under suitable regularity conditions, we establish both first- and second-order efficiency relative to the population oracle. Comparison with the corresponding two-stage rule shows that the additional recalibration yields sharper second-order control and reduces the influence of the pilot stage on the final sampling decision. The theory allows the ambient dimension to exceed the sample size under appropriate covariance and spectral conditions. Simulation studies demonstrate the strong finite-sample performance of the procedure across increasing dimensions and several dense covariance structures. As a real-data application, we conduct a retrospective study of gene-expression data from 32 cancer-type cohorts in The Cancer Genome Atlas, illustrating both cost-effective early stopping and settings in which additional observations are recommended.

We study multi-armed bandits (MAB) with multiple optimal arms, motivated by the fact that many practical decision making problems admit multiple correct answers. For $K$-armed bandits with $A$ optimal arms, we first provide a sharper analysis of previous sub-sampling algorithms (De Heide et al., 2021; Zhu and Nowak, 2020), establishing a $\tilde{O}\Big(\frac{K-A}{\sqrt{KA}}\sqrt{T} \Big)$ minimax regret, where $T$ is the total number of interactions and $\tilde O(\cdot)$ drops all constant and logarithmic factors, improving the previous $\tilde{O}(\sqrt{KT/A})$ regret. We then provide a matching lower bound up to logarithmic factors, indicating that our established rate is nearly minimax-optimal. We further show that the knowledge of $A$ up to $\tilde{O}(1)$ factors is necessary to achieve near-optimal regret, as near-optimal algorithms for one number of optimal arms must incur substantially larger regret than optimal regret for a smaller number. Overall, our results provide a comprehensive minimax characterization of $K$-armed bandits with $A$ over the entire range of $1 \leq A \leq K-1$.

High-frequency athlete-tracking data can be summarized as quantile cubes (Thomas and Hannig, Journal of Quantitative Analysis in Sports, 2026): counts of time spent in bins of velocity, acceleration and movement angle. In our application each athlete-match session has 100 such bins, with counts from under a hundred to several thousand. Whether a bin's negative binomial distribution is close to normal depends on its shape, the product of mean and dispersion, not on the count alone. We develop a hybrid generalized fiducial framework that uses exact negative binomial structural equations where the estimated shape is small and moment-matched normal approximations where it is large, providing uncertainty quantification for hundreds of regression and dispersion parameters without a prior. Only the normal components contribute to the fiducial Jacobian, which is block diagonal in closed form, and we prove a Bernstein-von Mises theorem for the correctly specified hybrid model. Simulations show near-nominal coverage, agreement with a flat-prior Bayesian analysis, and that the partition should be based on the estimated shape rather than on the mean or the observed responses. In 216 sessions from 17 professional women's soccer athletes, the estimated position and age effects concentrate at the velocity extremes and highest accelerations, but an athlete-level permutation check shows these patterns cannot be distinguished from label noise.

Grokking shows that fitting the training data and learning the underlying signal can occur at very different stages. However, existing theories offer limited quantitative insight into how this delayed generalization depends on inductive bias and signal structure. Our work addresses the gap by developing a statistical theory that characterizes how regularization geometry and signal sparsity govern generalization near interpolation. In particular, we focus on the prototypical setting of high-dimensional regression and identify regimes in which sparsity-promoting regularization makes exact interpolation much more accurate than approximate fitting. In strongly overparameterized noiseless problems, we prove a zero–one generalization law and construct a family of convex norms whose interpolators transition from the trivial risk of the all-zero predictor to exact recovery, while keeping the training error equal to $0$. Furthermore, when feature dimension and sample size are proportional, we provide a precise characterization of training and generalization errors along $\ell_r$-regularization paths. This in turn allows us to quantify the generalization gain that remains near interpolation: we show that this gain increases as the norm becomes more sparsity-promoting and as the target becomes sparser, with a sharp drop in generalization reached for noiseless data and $\ell_1$ regularization. Experiments on diagonal linear networks and transformers trained on modular arithmetic demonstrate the generality of our theoretical predictions. Finally, beyond grokking, our work reveals a statistical instability in minimum-norm interpolation: small perturbations in the regularization strength can lead to drastically different generalization, while preserving small training error.

We propose a novel local-polynomial estimator of the ratio $r=f/g$ of two $d$-dimensional densities $f$ and $g$, from which independent samples are available. The estimator is shown to achieve pointwise minimax optimal rates over Hölder classes of arbitrary smoothness index without additional logarithmic factors, and with smoothness being only assumed of $r$ but not of $f$ nor $g$. Our analysis remains valid for points on the boundary of the support of the distribution associated to $g$ under a mild geometric assumption on the (unknown) boundary. We also derive a concentration inequality for the estimator, which can be useful in applications to classification, and give a rate in the supremum norm, where the supremum is also taken over boundary support points. The smoothness class over which the rate is obtained is sufficiently large that the individual densities $f$ and $g$ cannot be consistently estimated uniformly over this class. Direct estimators of the partial derivatives of $r$ together with rates of convergence are also provided. We also obtain asymptotic normality of the estimators under quite general assumptions and again including boundary points, with a consistent variance estimator allowing for data-driven studentization. Moreover, we show how our estimator can be used to estimate the Kullback-Leibler information, in a construction which additionally uses debiasing. We provide the parametric rate together with asymptotic normality under sufficient smoothness of $r$ relative to the dimension $d$.

In the nearly unstable heavy-tailed regime of Jaisson and Rosenbaum, where the intensity of a linear Hawkes process converges to a rough square-root Volterra process, we ask what a record of length $T$ reveals about the relative distance to criticality $1-\rho$ and the exogenous rate $\mu$. For the first, the information is governed by $\mathcal I_T=T\mu(1-\rho)$, not by the number of events; $\mathcal I_T$ stays bounded although the event count diverges like $T^{2\gamma}$, so no estimator of the relative distance to criticality is uniformly locally consistent. For the second, the Fisher information diverges because the rough limit has an atom at zero: the asymptotic information geometry is singular in the exogenous direction and, unlike the light-tailed case, has no Feller-type threshold. We prove that the joint maximum likelihood estimator of $(\mu,\rho)$, whose likelihood is jointly concave, recovers $\mu$ consistently at the random rate $\{Q_T/\log\log Q_T\}^{1/2}$ given by the observed information $Q_T$, while its estimate of the relative margin is only tight: the exogenous rate is learned although the endogenous parameter it is coupled to is not. For empty-start records $Q_T$ diverges because the rough limit started at zero spends positive time at zero; on windows observed long after the start, still with the branching ratio unknown, learnability holds with probability tending to one as the window grows, through the ergodicity of the rough limit and the atom of its stationary law. We also prove finite-dimensional convergence of the stationary intensity to the stationary rough Volterra process.

The branching ratio $\rho$ of a self-exciting point process determines the branching margin $1-\rho$ to the critical boundary. We study inference on $\rho$ from a finite record for stationary Hawkes processes with completely monotone kernels whose slow relaxation mass obeys an envelope near zero rate. Slow excitation can then be exchanged for exogenous immigration; a spectral bound on the Kullback–Leibler rate and a coupling argument for the unobserved past give a minimax lower bound of order $T^{-\gamma/(2\gamma+1)}$ for a power-law envelope with tail exponent $\gamma$. A block-dispersion estimator whose block length grows as $T^{1/(2\gamma+1)}$ attains this rate. A quantitative central limit theorem, obtained from the cluster representation and a second-order Poincaré inequality on Poisson space, yields honest confidence intervals of nearly minimax expected length. For unknown $\gamma$, a Lepski-type rule adapts in probability at a logarithmic cost; free adaptation between memory classes is impossible, and honest intervals cannot adapt to lighter memory. Simulations are consistent with the rates.

We propose a general goodness-of-fit framework for distributions on separable metric spaces. Under suitable identifiability conditions, probability distributions are characterized by distance profiles, which motivates their use in goodness-of-fit testing, for simple and composite null hypotheses. For composite null hypotheses, parameter estimation is incorporated via a Bahadur-type expansion, and the asymptotic distribution of the empirical process for distance profiles is obtained under the null. We define test statistics based on this empirical process and derive their asymptotic null distributions. We further study the behavior of the proposed tests under fixed and local alternatives, establishing consistency results. Multiplier bootstrap procedures are developed, and their conditional asymptotic validity is established under both simple and composite null hypotheses. The methodology is illustrated with simulation studies and real-data applications for data on the sphere, hyperboloid, and simplex.

The classical information bottleneck (IB) measures the relevance of a representation $U$ of $X$ to a target $Y$ by $I(U;Y)$, which does not directly characterize the error of downstream decisions. For a binary hypothesis $Y$ inferred from many separately encoded observations, the optimal error exponent is the Chernoff information between the two conditional distributions of $U$ given $Y$. We study the mutual information constrained Chernoff bottleneck, which seeks an encoder that maximizes this Chernoff information subject to a rate constraint $I(U;X) \leq R$. We show that its optimal value $C(R)$ increases strictly up to $R = H(V)$, where $V$ merges the symbols of $X$ with equal likelihood ratio, remains at the uncompressed exponent beyond, and, unlike the IB curve, need not be concave. We further show that $k+1$ outputs suffice to attain $C(R)$, where $k$ is the cardinality of $V$. We propose an alternating algorithm that updates the encoder via a generalized Blahut–Arimoto algorithm and the Chernoff parameter $s$ via a nonlinear equation, and prove that its iterates remain feasible, with nondecreasing and convergent Chernoff information. Numerical experiments confirm the theory, and on real topic-detection data from the 20 Newsgroups corpus, compressing each word to only $17\%$ of its entropy retains $90\%$ of the error exponent and nearly the accuracy of the uncompressed classifier.

We introduce a novel proof strategy for semiparametric Bernstein-von Mises theorems based on Stein's method combined with an information-geometric framework. Rather than controlling the effect of the prior on the asymptotic marginal posterior distribution of the functional of interest through the stability of an integrated likelihood under a suitable perturbation, we characterise its influence through the prior-weighted divergence of suitable vector fields over the statistical model. Applying Stein's method in this setting instead of the usual Laplace-transform approach yields explicit non-asymptotic upper bounds on the bounded-Lipschitz distance between marginal posterior distributions and the corresponding Gaussian limits predicted by semiparametric efficiency theory. These bounds consist entirely of local scalar differential quantities associated with the functional, the prior and a chosen vector field-typically related to the efficient influence function-evaluated over posterior contraction sets. We apply the theory to quadratic functionals in Gaussian white-noise models and to linear, quadratic and general integral functionals in histogram density models, with both conjugate and non-conjugate priors. For linear and quadratic functionals, the same theory identifies the differential term responsible for posterior bias and allows us to remove it through a natural functional correction, yielding Bernstein-von Mises-type theorems in regimes where they may fail for the original uncorrected functional.

Projection pursuit searches for a direction along which the data look least Gaussian. When the observation space contains a large Gaussian complement, the empirical objective can be minimized by a direction that carries no signal, with empirical kurtosis as low as at the truth. Sample splitting exposes rather than repairs this failure. Appending coordinates independent of the latent regime degrades the search while leaving Bayes recoverability unchanged. Restricting the search to the column space of a known forward operator removes the failure exactly on the negative-kurtosis branch. Estimating a principal subspace from the data is the alternative. In a controlled two-component model, the leading sufficient scalings differ in the gain with which the operator transmits the discriminant: $\varsigma^{-4}$ for covariance-spike estimation and $\varsigma^{-8}$ for fourth-moment search. At fixed search dimension, the measured threshold ratio collapses onto $n/p^2$ with exponent $0.156$, close to the predicted $1/8$. This is an empirically supported scaling motivated by sufficient bounds, not a proved asymptotically tight law. When the search dimension is varied, the measured exponent is $0.325$, substantially larger than $1/8$, and the tested range does not identify its functional form. The crossing location also depends on calibration and model configuration. Under a downstream excess-error criterion, the scaling largely disappears.

Basis-restricted priors in Bayesian inverse problems can lose coverage when the truth has components outside the basis. We show that estimating the observation-noise variance can hide this loss. Under a linear forward model, when the in-span prior variance dominates the noise, the maximum-likelihood noise estimate absorbs the out-of-basis energy in the complement of the model range. Residual-magnitude and observation-coverage checks then stay near nominal while field coverage falls. We study repeated problems sharing one forward operator and one basis, fixed independently of the tested data. After projection onto the complement, and conditionally on the fitted noise scale, every exact test is a test of the scale-free direction of the residuals. We test the shape of their sample spectrum with John's sphericity statistic. Under Gaussian noise its null model is exact at finite sample size, and we derive its null mean and its power at proportional dimension. On synthetic problems and in a preregistered GEBCO topography study, the test detects structured out-of-basis variation that cross-validation and observation-coverage checks largely miss. It cannot detect Gaussian out-of-basis variation that is isotropic in the complement, since that is indistinguishable from a change of noise scale.

Goodness-of-fit testing is a basic tool for assessing whether a fitted procedure has captured the systematic information contained in the covariates. While traditional theory has largely focused on parametric regression models, modern data analysis increasingly relies on flexible black-box learners, whose predictive success alone is insufficient to assess model accuracy. In this paper, we propose SPARK, a general framework for goodness-of-fit testing that applies to traditional statistical models and general black-box learning procedures, continuous and binary responses, and low- and high-dimensional predictors. Based on a debiasing strategy, the residuals from an initial fit of a learning procedure are projected onto nearly orthogonal directions to extract any remaining signal. To capture information across all projection directions, we propose a kernel-based projection method and establish both its asymptotic properties and the consistency of a bootstrap procedure. Comprehensive simulations and real data analyses illustrate the effectiveness and flexibility of our proposed method.

In factorial randomized experiments with $K$ binary treatment components, we study which baseline covariates modify component main effects and interactions. We define effect-modifier cells $(j,S)$, $j\in[p]$ and $\emptyset\neq S\subseteq[K]$, through the coefficients $\theta^\ast_{jS}$ of population linear projections of conditional Walsh–Hadamard contrasts $\tau_{S}$ onto standardized covariates. Using baseline-robust direct scores, we develop frequentist inference for all $p(2^K-1)$ cells. First, we prove a debiased central limit theorem with global insensitivity: the residualization baseline enters the cross-fitted score through an exactly conditionally mean-zero perturbation. The influence-function remainder therefore requires only foldwise weighted $L_2$ consistency of the baseline estimator, with no prescribed polynomial rate or nuisance product-rate condition. A separate, explicit penalty-scale condition controls the baseline's effect on the score Lasso. Second, under explicitly stated uniform nuisance, baseline, studentization, and influence-array conditions derived from bounded sparse-regression primitives, we establish a high-dimensional Gaussian approximation over the full grid, calibrated by a Walsh-spectral multiplier bootstrap. This yields simultaneous confidence bands and Romano–Wolf step-down selection of modified cells with strong familywise error control, without effect-heredity assumptions. The factorial design is fixed in the asymptotic analysis, while the covariate dimension may grow; all $2^K-1$ contrasts are covered jointly.

We study multiple change-point detection in multivariate time series whose distributions change in a piecewise constant manner. Distributional changes can manifest across different moment orders, from shifts in the mean and covariance to changes in higher-order moments. Higher-order moments capture increasingly rich distributional features but become difficult to estimate in high dimensions. Our tensor representation unifies moments of different orders within a common linear algebraic framework, enabling a new method to detect changes in moments of all orders up to a prescribed fixed order $p$. The resulting procedure accommodates temporal dependence and allows the dimension of the time series to grow with the sample size. Under suitable regularity conditions, the proposed procedure achieves a localization error rate that matches a newly developed minimax lower bound. We further derive limiting distributions under both nonvanishing and vanishing moment jumps and construct asymptotically valid confidence intervals in the vanishing-jump regime. Numerical experiments and real-data analyses demonstrate the method's effectiveness in detecting moment changes and a range of distributional shifts.

Offline reinforcement learning seeks optimal decision rules from previously collected data. In some applications, a decision can be an entire function, such as a fluence map in radiation therapy or a smooth movement trajectory in robotics. In this paper, we study the finite-sample theory for fitted Q-iteration (FQI) with functional actions in a discounted infinite-horizon setting. Three major difficulties arise in this setting: first, the absence of a Lebesgue probability density for functional actions complicates coverage descriptions; second, conventional coverage requirements can be restrictive; and third, the large functional action space makes greedy optimization in FQI challenging. To address these difficulties, we study smoothness-regularized policy search under a critic-relative coverage condition. This condition measures how well logged data distinguish relevant action-value differences without requiring an action density. Our main theorem gives finite-sample guarantees for learned-policy regret relative to the best value within a fixed smooth class of functional-action policies. The results allow trajectory lengths to be either bounded or growing and Q-functions to be fitted by either functional-input kernel ridge regression or adaptive functional neural networks. For a few examples, we can obtain polynomially decaying regret bounds in the number of logged transitions, up to logarithmic factors, with logarithmically many FQI iterations. Numerical experiments show gains of learned functional-action policies over constant-action policies and support our adoption of a critic-relative coverage condition.

We study the computational complexity of hypothesis testing in the spiked Wigner model, a prototypical model for detecting low-rank structure in a large random matrix. Below the "BBP" eigenvalue transition, it is expected that strong detection — with both type I and II errors vanishing — requires exponential time. Assuming this as a conjecture, we determine the limits of polynomial-time weak detection, exactly characterizing the possible tradeoffs between type I and II errors. Specifically, the optimal tradeoff is achieved by a particular linear spectral statistic. Thus, the question of weak detection is entirely reduced to that of strong detection. The proof builds on ideas of Nagda-Raghavendra (2025) and Moitra-Wein (2025). The low-degree likelihood ratio (LDLR) plays a key role: any test that slightly beats the LDLR can be boosted to have an even higher success probability. This leads us to establish a computational analogue of the Neyman-Pearson lemma for a subclass of additive Gaussian models: for a given super-polynomial runtime, the best possible tradeoff between type I and II errors is either the one achieved by thresholding the LDLR, or the trivial tradeoff that results from strong detection.

Bell and Cover showed that an investor who multiplies the initial unit of capital by an independent uniform random variable on $(0,2)$, and then uses the log-optimal portfolio, wins a head-to-head wealth comparison with probability at least one half against every independently randomized competitor. We explain very simply how this result transfers to testing by betting: for any composite null $\mathcal P$ and simple alternative $Q$, denoting $E^*$ as the corresponding numeraire e-variable, we show that $UE^*$ exceeds any other e-variable $E$ with probability at least half. Interestingly, we show that this competitive optimality result is actually equivalent to the numeraire inequality $\mathbb E_Q[E/E^*]\leq1$, and in general randomization only helps the numeraire and fails to improve the competitive advantage of an arbitrary e-variable. Under optional stopping with or without knowledge of $U$, we emphasize a key distinction between e-process validity and competitive optimality. We also show that competitive optimality comes at the price of expected log wealth and power: thresholding $UE^*$ at $1/\alpha$ has sharp size at most $\alpha/2$, but the factor of two actually disappears under optional stopping. Even after correcting for this factor of two, the test is dominated in conditional rejection probability by randomizing the testing threshold (randomized Markov's inequality). Thus, Bell-Cover randomization is optimal for a specific competitive objective, at the cost of others.

For probability measures on $\mathbb{R}^d$ supported in fixed compact sets, we prove that quadratic optimal transport plans are quantitatively stable in Wasserstein distance under perturbation of both marginal measures, assuming that one of the initial measures satisfies an upper Ahlfors regularity condition with exponent strictly greater than $d-1$. Under the same assumption, we also prove novel quantitative stability results for optimal transport maps. Furthermore, we prove a coercivity theorem which states that any probability measure on the product space must be quantitatively close to the set of optimal plans, if it has similar marginals and a similar transport cost to the optimal value. Finally, we prove a quantitative uniqueness theorem which acts as a quantitative counterpart to Brenier's theorem. The Wasserstein diameter of the set of optimal plans is controlled by the distance of one marginal measure to a regular measure for which uniqueness holds. In this way, "almost uniqueness" of optimal plans is quantified by the source measure being "almost regular". Examples are provided which prove that the exponents of all estimates are sharp.

We study the supervised classification problem for interacting particle systems (IPS) from trajectory observations. The IPS belongs to one of K classes, each characterized by a distinct interaction drift. Given a learning dataset, the goal is then to predict the class label of a newly observed system based on the trajectory of a single particle. This setting raises several statistical challenges. First, the particles within each system interact and are therefore dependent. Second, although the full particle system is Markovian, the trajectory of a single particle is not. To address these difficulties, we exploit the McKean-Vlasov limit of the IPS, which describes the dynamics of a typical particle as the number of particles tends to infinity. We propose a plug-in classification procedure based on estimating the interaction drift associated with each class from discretely observed trajectory data. Hence, for each class, we provide a new nonparametric estimator of the drifts by minimizing a ridge-regularized least-squares contrast over a B-spline basis. In particular, our theoretical findings reveal that the convergence rate of the resulting classifier is of order N -1/6+$\epsilon$ for any $\epsilon$ > 0, where N denotes the number of particles in each system. Numerical experiments illustrate the performance of the proposed method and show that it outperforms an end-to-end neural network baseline that does not exploit the underlying particle-system structure.

Virtual simulators are widely used for studying complex physical phenomena, from particle collisions to rocket propulsion. Such "computer experiments" can be highly time-intensive, and a Bayesian surrogate model can be used for efficient emulation with reliable uncertainty quantification. To train accurate surrogates with a limited sample size $n$, recent work has explored the incorporation of monotonicity (or isotonicity) information, which can often be elicited from physical systems. In practical applications with many input variables, however, existing Bayesian isotonic models can face statistical and computational limitations, which may result in worse performance compared to models that do not incorporate isotonicity. We propose a new transformed additive isotonic model (TAIM), which aims to tame this "curse-of-dimensionality". TAIM makes use of a flexible transformed additive isotonic modeling framework, which leverages a data-estimated link transformation and a monotone basis model with spike-and-slab priors on basis weights. Prediction-wise, TAIM achieves (up to log factors) a posterior contraction rate of $O(n^{-1/3})$ when the true black-box function is in a transformed additive isotonic form with mild smoothness conditions. Such a rate does not depend on the input dimension $d$ for terms involving $n$, which softens the effect of dimensionality on posterior predictions. Computation-wise, TAIM allows for efficient posterior inference via a carefully designed Gibbs sampler, where each sampling iteration requires only linear work in $d$. We further present an extension of TAIM that can model potential deviations from transformed additivity. Numerical experiments and two applications show the effectiveness of TAIM for isotonic surrogate modeling with many input variables.

The Baringhaus–Henze–Epps–Pulley (BHEP) tests form a widely applicable class of affine invariant and consistent tests for multivariate normality. Under the null hypothesis, the test statistic converges to a weighted sum of independent chi-squared random variables. Although closed-form expressions for the first three cumulants of this limiting distribution are known for arbitrary dimension $d$ and smoothing parameter $\beta$, an explicit expression for the fourth cumulant has so far been available only in a special univariate case. We derive the fourth cumulant for arbitrary $d\geq1$ and $\beta>0$ and use the resulting first four moments to construct Johnson and Pearson approximations to the limiting null distribution. These yield essentially instantaneous analytic approximations to BHEP $p$-values, without numerical eigenvalue calculations. Using the recently obtained complete spectrum as a benchmark, we show that both four-moment approximations are highly accurate over a broad range of dimensions, smoothing parameters and upper-tail probabilities, and substantially improve on two- and three-parameter lognormal approximations. Monte Carlo results indicate good finite-sample calibration for many parameter combinations, although convergence to the limiting distribution may be slow in higher dimensions for extreme values of the smoothing parameter.

We study the maximum likelihood estimation of the coefficient {\beta} in well-specified Poisson regression. Using tools from empirical process theory and random conic geometry, we show that the probability of existence of the maximum likelihood estimator (MLE) exhibits a sharp phase transition at the threshold n > d. We then determine a minimum threshold exponential in the norm of {\beta} on the sample size n to guarantee with high probability an excess risk of the asymptotic order d/n. We reveal the existence of an intermediate regime in Poisson regression, when n is larger than d but smaller than this exponential threshold, where the MLE exists but does not achieve the optimal rate d/n. We close the gap between the two regimes up to a d^{1+\epsilon} term with {\epsilon} \in (0, 1) by providing an upper bound on the distance between the MLE and {\beta} whenever n > d^{1+\epsilon}. Along the way, we provide two generalizations of well-known PAC-Bayes inequalities regarding sub-Gamma random vectors and sub-Gamma random matrices that are of independent interest and that we use extensively to prove the main results of the present paper.

We propose a Bayesian nonparametric factor model that infers the number of factors, induces row-wise sparsity, and merges redundant dictionary elements via exact clustering. A Dirichlet process prior is placed on the columns of an overcomplete loading matrix and fully marginalized to an exact Pólya urn, avoiding stick-breaking and auxiliary variables. A spike-and-slab base measure allows entire factors to be exactly zero. The model uniquely combines exact zeros, exchangeability over columns, and exact merging within a single marginalized Dirichlet process, unlike CUSP, MGP, or the beta process. An exact Gibbs sampler with canonical relabeling and parallel C/MPI implementation is developed. We prove posterior contraction at rate $\sqrt{M s_0 \log n / n}$ for the covariance matrix, and in the fixed-dictionary setting obtain the minimax optimal rate $\sqrt{s_0 \log n / n}$ plus underfitting consistency; the overfitting direction is an open conjecture. The spike-and-slab is essential: without it the effective dimension scales as $pM$, yielding a slower rate. Simulations show the method is the only fully adaptive approach to recover the true rank, achieving the smallest covariance, loading, and signal-reconstruction errors, beating an oracle baseline. On van 't Veer breast cancer data ($n=97$, $p=1213$), the posterior concentrates on eight interpretable programmes; seven pass coherence and two pass Bonferroni-corrected Hallmark enrichment.

Classical statistical theory is insufficient to explain the successes of general-purpose AI models, because it depends on handcrafted inductive biases that it cannot justify. No Free Lunch (NFL) theorems force any learner that beats chance on some environments to underperform on others. We might hope that past experience informs which environments to expect, but NFL applies equally to meta-learning. Thus, any method that makes meaningful predictions necessarily begins with an inductive bias external to the data. Choosing to bias toward short programs yields Solomonoff induction (SI), whose performance is competitive against all computable learners - albeit up to "constants" that become large when comparing against specialized methods that exploit background information. We therefore relativize SI to an information vantage point, biasing toward short programs with access to all preexisting information. This reframes the inductive bias: instead of seeking some absolute notion of simplicity, we favor accessibility with respect to our vantage point. An algorithm can only outpredict the relativized SI to the extent that its code contains additional information about the data, and no algorithm can generate such information. While SI is incomputable and hence not a practical algorithm, it provides a formal optimum for inference in the limit of infinite compute, and there is evidence to suggest that frontier AI systems roughly approximate it. Thus, the only known answer to meta-NFL is rooted in algorithmic information theory, which we should expect to play a fundamental role in explaining the generalization behavior of modern (and future) AI systems.